Arc-length solvers (structsolve.arc_length)

Arc-length solvers for non-linear static analyses

The Riks and Crisfield methods share the implementation of _solver_arc_length(), and are selected with Analysis.NL_method = 'arc_length_riks' or Analysis.NL_method = 'arc_length_crisfield', as does the indirect displacement control, Analysis.NL_method = 'displacement_control'. The module-level constants control the arc-length scaling and adaptation.

structsolve.arc_length._solver_arc_length(an, method, silent=False)

Arc-length solver

Traces the equilibrium path of the model of an with an arc-length method, in which the load factor is an unknown of each step besides the solution vector, such that limit points can be passed. method selects how the load factor of each iteration is computed, see arc_length_riks._solver_arc_length_riks() and arc_length_crisfield._solver_arc_length_crisfield().

The steps are controlled by the attributes initialInc, minInc, maxInc and maxArcLength of an, and by the module-level constants; the iterations of each step by convergence, relTOL, absTOL, maxNumIter and too_slow_TOL, with the checks of newton_raphson._check_convergence() and newton_raphson._check_divergence(), where the convergence criterion also considers the external force vector of the increment of the load factor of the current step; and an.modified_NR, an.compute_every_n and an.kT_initial_state control how often the tangent stiffness matrix is updated. A step that fails is repeated with the arc length multiplied by INC_CUT_FACTOR. The analysis stops when:

• the load factor reaches an.final_load_factor, 1.0 by default, the last solution being computed at exactly this load factor; never with an.final_load_factor = None

• an.stop_criterion(an) returns True after a converged step, see stopping

• the cumulative arc length reaches an.maxArcLength

• the arc length becomes smaller than an.minInc

• MAX_NUM_STEPS steps converged

• the tangent stiffness matrix at the last converged state is singular

The callables of an that declare the optional keyword arguments inc, c or lbd receive the current load factor and configuration, see callbacks.

Parameters:

anAnalysis

Analysis object, the converged load factors and solutions are appended to an.increments and an.cs.

methodstr

'riks', 'crisfield', 'displacement_control' or 'generalized_displacement_control'.

silentbool, optional

A boolean to tell whether the log messages should be printed.

Notes

The equilibrium path \(\boldsymbol{F}_{int}(\boldsymbol{c}) = \lambda \boldsymbol{F}_{ext}\) is traced with the arc-length constraint

\[\frac{\Delta \boldsymbol{c}^\top \Delta \boldsymbol{c}}{\boldsymbol{u}_{ref}^\top \boldsymbol{u}_{ref}} + \psi^2 \Delta\lambda^2 = \Delta s^2\]

where \(\Delta \boldsymbol{c}\) and \(\Delta\lambda\) are the increments of the current step and \(\boldsymbol{u}_{ref}\) is the linear solution for \(\lambda=1\), which makes the arc length \(\Delta s\) independent of the units of the model. In the linear regime \(\Delta s \approx \sqrt{1 + \psi^2} \Delta\lambda\). It is the spherical arc-length constraint of [4] (Section 9.3.2, Eq. (9.19), p. 271), where the displacement term is not scaled and the load term is \(\psi^2 \Delta\lambda^2 \boldsymbol{F}_{ext}^\top \boldsymbol{F}_{ext}\). \(\psi\) is the constant PSI, and \(\boldsymbol{u}_{ref}^\top \boldsymbol{u}_{ref}\) is replaced by 1 when it is null.

With an.arc_length_weights \(= \boldsymbol{w}\), the products of displacements in the constraint, and in the inner product below, are weighted, \(\boldsymbol{a}^\top \boldsymbol{W} \boldsymbol{b}\) with \(\boldsymbol{W} = \mathrm{diag}(\boldsymbol{w})\), also in the scale \(\boldsymbol{u}_{ref}^\top \boldsymbol{W} \boldsymbol{u}_{ref}\). When the unknowns mix translations with their derivatives or with rotations, of other units and magnitudes, the weights 1 for the translations and 0 for the others measure the arc length on the translations only. The arc length is then independent of the units of the degrees of freedom of null weight, as is Crisfield’s convergence criterion, see convergence.

With the inner product

\[\langle (a, \alpha), (b, \beta) \rangle = \frac{\boldsymbol{a}^\top \boldsymbol{b}}{\boldsymbol{u}_{ref}^\top \boldsymbol{u}_{ref}} + \psi^2 \alpha \beta\]

the constraint reads \(g = \frac{1}{2}(\langle (\Delta c, \Delta\lambda), (\Delta c, \Delta\lambda) \rangle - \Delta s^2) = 0\).

Each step starts from a tangent predictor, whose direction follows the previous increment, allowing limit points to be passed. With \(\delta \boldsymbol{c}_q = \boldsymbol{K}_T^{-1} \boldsymbol{F}_{ext}\) at the last converged state, the predictor

\[\Delta\lambda = \pm \frac{\Delta s}{\sqrt{\langle (\delta c_q, 1), (\delta c_q, 1) \rangle}}, \qquad \Delta \boldsymbol{c} = \Delta\lambda \delta \boldsymbol{c}_q\]

satisfies the constraint, with the sign of \(\langle (\delta c_q, 1), (\Delta c_{prev}, \Delta\lambda_{prev}) \rangle\), where \((\Delta c_{prev}, \Delta\lambda_{prev})\) is the increment of the previous converged step, and a positive sign at the first step. The initial arc length, \(\Delta s_0 = initialInc \sqrt{\langle (\delta \boldsymbol{c}_q, 1), (\delta \boldsymbol{c}_q, 1) \rangle}\) at the origin, corresponds to a load factor increment of an.initialInc along the initial tangent.

At each iteration, the residual is \(\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{c} + \Delta \boldsymbol{c}) - (\lambda + \Delta\lambda) \boldsymbol{F}_{ext}\), with \(\partial \boldsymbol{R}/\partial \lambda = -\boldsymbol{F}_{ext}\), and the bordered system \(\boldsymbol{K}_T \delta \boldsymbol{c} - \delta\lambda \, \boldsymbol{F}_{ext} = -\boldsymbol{R}\) is solved with two solutions of the tangent stiffness matrix, \(\delta \boldsymbol{c}_R = -\boldsymbol{K}_T^{-1} \boldsymbol{R}\) and \(\delta \boldsymbol{c}_q = \boldsymbol{K}_T^{-1} \boldsymbol{F}_{ext}\); the increment is updated to \((\Delta \boldsymbol{c} + \delta \boldsymbol{c}_R + \delta\lambda \, \delta \boldsymbol{c}_q, \Delta\lambda + \delta\lambda)\), where the correction of the load factor \(\delta\lambda\) comes from the constraint:

• method='riks': the constraint is linearized at the current increment (updated normal plane), which is Newton-Raphson applied to the augmented system of the equilibrium equations and the constraint ([4], Section 9.3.2, p. 272), with quadratic convergence when \(\boldsymbol{K}_T\) is updated at every iteration and is the exact Jacobian. It follows Riks [1], who applies Newton’s method to the equilibrium equations augmented with an equation that controls the arc length, and states its quadratic convergence (Sections 3.2 and 3.3, Eqs. (3.7), (3.11) and (3.17), pp. 537-539). The equation of Riks is linear (Section 3.1, Eqs. (3.4) and (3.9), pp. 537-538): the plane normal to the unit tangent of the path at the last converged point, at the distance \(\Delta s\) from it, in the unscaled Euclidean space of the load parameter and the generalized coordinates; the plane is fixed during the iterations of a step and redefined at each converged point. Without the mismatch \(g\) of the constraint, it is the updated normal plane of Ramm [3] (p. 75); [4] (Section 9.3.2.2, Eq. (9.29), pp. 274-275) attributes the mismatch term to Schweizerhof and Wriggers and to Riks:

  \[\delta\lambda = -\frac{g + \langle (\Delta c, \Delta\lambda), (\delta c_R, 0) \rangle}{\langle (\Delta c, \Delta\lambda), (\delta c_q, 1) \rangle}\]

• method='crisfield': the quadratic (spherical) constraint is solved exactly, choosing the root with the smallest angle to the current increment. When the roots are complex, the linearized constraint is used, see [2]. With \(\boldsymbol{u} = \Delta \boldsymbol{c} + \delta \boldsymbol{c}_R\), the constraint \(\langle (u + \delta\lambda \, \delta c_q, \Delta\lambda + \delta\lambda), (u + \delta\lambda \, \delta c_q, \Delta\lambda + \delta\lambda) \rangle = \Delta s^2\) is the quadratic equation \(a \, \delta\lambda^2 + b \, \delta\lambda + c = 0\) with \(a = \langle (\delta c_q, 1), (\delta c_q, 1) \rangle\), \(b = 2 \langle (u, \Delta\lambda), (\delta c_q, 1) \rangle\) and \(c = \langle (u, \Delta\lambda), (u, \Delta\lambda) \rangle - \Delta s^2\). Its roots are computed as \(q/a\) and \(c/q\), with \(q = -\frac{1}{2}(b + \mathrm{sign}(b) \sqrt{b^2 - 4 a c})\), and the root that maximizes the inner product of the new increment with the current one \((\Delta c, \Delta\lambda)\) is chosen, i.e. the smallest angle, since both have the length \(\Delta s\).

• method='displacement_control': indirect displacement control, the arc-length constraint being replaced by the linear constraint

  \[\frac{\boldsymbol{\psi}^\top \Delta \boldsymbol{c}}{s_{ref}} = \Delta s\]

  with \(\boldsymbol{\psi}\) = an.control_vector and \(s_{ref}\) = an.control_scale, by default \(|\boldsymbol{\psi}^\top \boldsymbol{u}_{ref}|\), such that in the linear regime \(\Delta s = \Delta\lambda\), which must then not be orthogonal to \(\boldsymbol{u}_{ref}\) (cosine above CONTROL_MIN_COS); a unit vector \(\boldsymbol{\psi}\) controls the displacement of one degree of freedom, e.g. the end shortening, and a buckling mode controls its amplitude. The correction of the load factor satisfies the constraint exactly:

  \[\delta\lambda = \frac{s_{ref} \Delta s - \boldsymbol{\psi}^\top (\Delta \boldsymbol{c} + \delta \boldsymbol{c}_R)}{\boldsymbol{\psi}^\top \delta \boldsymbol{c}_q}\]

  and the predictors are scaled to satisfy it, the tangent predictor \(\Delta\lambda = s_{ref} \Delta s/(\boldsymbol{\psi}^\top \delta \boldsymbol{c}_q)\), without the sign rule, the sign of the load factor increment following from the one of \(\boldsymbol{\psi}^\top \delta \boldsymbol{c}_q\). The step length \(\Delta s\) plays the role of the arc length, with an.initialInc the first one, \(\Delta s\) keeping the sign for which the controlled quantity increases with the load factor at the origin. The path can be traced through the limit points of the load and of the other degrees of freedom, e.g. snap-backs, as long as the controlled quantity increases monotonically along it; at a limit point of the controlled quantity \(\boldsymbol{\psi}^\top \delta \boldsymbol{c}_q\) vanishes and the iterations fail or converge to another part of the path. For the indirect displacement control, see de Borst [5]. The check against the doubling back below does not apply, since the previous converged point does not satisfy the constraint.

• method='generalized_displacement_control': the generalized displacement control of Yang and Shieh [6]. With \(\delta \boldsymbol{c}_q^{(k)}\) the tangent \(\boldsymbol{K}_T^{-1} \boldsymbol{F}_{ext}\) at the start of step \(k\) and \(\langle a, b \rangle_W = \boldsymbol{a}^\top \boldsymbol{W} \boldsymbol{b}\), with the weights of an.arc_length_weights, the predictor of step \(k\) is \(\Delta\lambda = s_k \Delta s \sqrt{|GSP_k|}\), \(\Delta \boldsymbol{c} = \Delta\lambda \, \delta \boldsymbol{c}_q^{(k)}\), with the generalized stiffness parameter

  \[GSP_k = \frac{\langle \delta c_q^{(1)}, \delta c_q^{(1)} \rangle_W}{\langle \delta c_q^{(k-1)}, \delta c_q^{(k)} \rangle_W}, \qquad GSP_1 = 1\]

  which decreases with the stiffness along the path and changes sign where the tangent turns back, at the limit points of the load; the sign \(s_k = s_{k-1} \, \mathrm{sign}(GSP_k)\), \(s_1 = 1\), then reverses the load increment. The corrections satisfy

  \[\langle \delta c_q^{(k-1)}, \delta \boldsymbol{c}_R + \delta\lambda \, \delta \boldsymbol{c}_q \rangle_W = 0, \qquad \delta\lambda = -\frac{\langle \delta c_q^{(k-1)}, \delta c_R \rangle_W}{\langle \delta c_q^{(k-1)}, \delta c_q \rangle_W}\]

  with \(\delta c_q^{(0)} = \delta c_q^{(1)}\), i.e. they are orthogonal to the tangent of the previous step, which does not depend on the units of the load nor on a scale of the displacements. The step length \(\Delta s\) is a load factor, an.initialInc the first one and the largest, cut and adapted as the arc length, but never increased beyond an.initialInc, since growing it as the arc length would defeat the scaling by \(GSP_k\), which shortens the steps where the structure softens: with longer steps the iterations may converge beyond a limit point on another part of the path. an.maxArcLength bounds the sum of the \(\Delta s\). With an.arc_length_weights = None, \(GSP_k\), the predictor and \(\delta\lambda\) are Eqs. (33), (34) and (30) of [6], with \(\delta c_q^{(0)} = \delta c_q^{(1)}\) as there. [6] multiplies the load increment by \(-1\) when \(GSP_k < 0\), its step 2d, which keeps the displacements going forward only if the reversal is relative to the previous step, \(s_k = s_{k-1} \, \mathrm{sign}(GSP_k)\) as here, since \(GSP_k\) is negative only at the step after each limit point. The weights and the adaptation of \(\Delta s\) are additions of this implementation: in [6] \(\Delta s\) is the first load increment throughout. The check against the doubling back below does not apply.

With an.arc_length_line_search = True, the whole correction \((\delta \boldsymbol{c}_R + \delta\lambda \, \delta \boldsymbol{c}_q, \delta\lambda)\) is multiplied by the first factor \(\eta\) of LINE_SEARCH_ETAS, \(1, 1/2, 1/4\), for which the norm of the residual of the convergence criterion, see newton_raphson._check_convergence(), decreases, or by the factor of the smallest norm when none decreases it. The internal force vector of the accepted point is reused in the residual of the next iteration, such that the full correction costs no additional evaluation. With \(\eta < 1\) the iterate leaves the constraint surface, the next correction returning to it; the converged point may thus be closer than \(\Delta s\) to the last one when the last correction was reduced. Contrary to the line search of newton_raphson._NR_iterations(), there is no sufficient-decrease test: the residual of the equilibrium equations alone is not the merit function of the augmented system of the equilibrium equations and the constraint, along whose Newton direction it need not decrease.

The step fails when \(\delta\lambda\) cannot be computed or is not finite, when the iterations are stopped by newton_raphson._check_divergence(), earlier with an.early_divergence = True, and when the converged increment points against the predictor \((\Delta c_{pred}, \Delta\lambda_{pred})\) of the step, \(\langle (\Delta c, \Delta\lambda), (\Delta c_{pred}, \Delta\lambda_{pred}) \rangle \le 0\). The previous converged point is also at the distance \(\Delta s\) from the last one, and with long steps or near a limit point the iterations may converge to it and then trace the path backwards, the “doubling back” that [2] (p. 56) avoids with the choice of the root. The accepted points are thus on the side of the tangent predictor of the hyperplane normal to it, in the inner product above, through the last converged point, as the points of the plane of Riks [1], normal to the tangent at the distance \(\Delta s\) from the last converged point. With full Newton-Raphson, \(\boldsymbol{K}_T\) is updated at every iteration at the current iterate \((c + \Delta c, \lambda + \Delta\lambda)\); with an.modified_NR = True the factorization at the last converged state is used for the first an.compute_every_n iterations and the tangent is then updated every an.compute_every_n iterations; at the first step with an.kT_initial_state = False the linear constitutive stiffness matrix replaces the tangent of the initial state.

After each converged step, the arc length is multiplied by sqrt(DESIRED_NUM_ITER/num_iter), bounded by MIN_GROWTH_FACTOR and MAX_GROWTH_FACTOR and limited by an.maxInc, where num_iter is the number of corrections of the step. The factor is at most 1 after a step that was repeated, and the limit is in fact max(an.maxInc, initial arc length). When a step fails, the arc length is multiplied by INC_CUT_FACTOR and the step is repeated from the same tangent predictor direction. The square root of the ratio is the rule recommended by Ramm [3] (p. 76), who found that the ratio itself, as used by Crisfield [2] (Eq. 15), makes the number of iterations oscillate from step to step.

The tangent stiffness matrix is factorized with an.linear_solver, see static._factorize(). With an.inertia = True, it is factorized at every converged state, the factorization being that of the predictor of the next step, and its number of negative eigenvalues and the sign of its determinant are recorded in an.num_negative_eigvals and an.det_signs, see newton_raphson._converged_state(): on a stable branch without prestress-induced instabilities the count is zero, it increases by one at each limit point of the load beyond which the branch is unstable, and at each bifurcation point crossed without branching. A perturbation whose shape keeps the path on a branch of a given symmetry, e.g. axisymmetric, may thus leave many negative eigenvalues unnoticed without this count.

With an.max_inertia_jump \(= j\), which requires an.inertia and an.linear_solver = 'pardiso', a converged step whose number of negative eigenvalues differs by more than \(j\) from that of the last converged state is repeated with a shorter arc length, as a step that failed, unless the shorter arc length would be smaller than an.minInc. With \(j = 1\) the eigenvalues cross zero one at a time, and a step over several critical points, e.g. beyond the limit load of a cylinder under axial compression, is cut back to them.

With an.branch_switching \(= n > 0\), which requires an.inertia, the path is switched to the bifurcated branch at the first \(n\) bifurcation points crossed, with every method. A critical point is crossed when the number of negative eigenvalues changes between two converged states, increases or decreases, or, when it is not known (an.linear_solver = 'splu'), when the sign of the determinant changes. It is a bifurcation point when the load factor varies in the same direction along the tangent predictors of the inner-product rule at the two converged states, before and after the critical point, a limit point when the direction is reversed. At a bifurcation, the eigenvector \(\boldsymbol{\phi}\) of the eigenvalue that became negative, or positive when the count decreased, is computed with _critical_eigenvector(), with the sign for which \(\boldsymbol{\phi}^\top \boldsymbol{c} \ge 0\), i.e. the direction of the deformation already present, e.g. of an imperfection, and the next step starts from the predictor \((\Delta s \hat{\phi}, 0)\) along it, \(\hat{\phi}\) normalized in the inner product above, with \(\Delta s\) the larger of the current arc length and that of the last step, which crossed the critical point. With method='displacement_control' and 'generalized_displacement_control', whose step lengths are not arc lengths, \(\Delta s\) is the length of the last increment in the inner product above, multiplied by the growth of the step length since then when it grew, as for the arc length. The step is corrected with the spherical constraint of method='crisfield' of radius \(\Delta s\), whose root closest to the current increment is on the bifurcated branch, which passes through the critical point, at a distance smaller than \(\Delta s\), and is tangent to \(\boldsymbol{\phi}\) at it, whereas the linearized constraint of method='riks' is singular there for a symmetric bifurcation. The step is accepted when the cosine between its increment and the predictor is at least BRANCH_SWITCH_MIN_COS, otherwise it is repeated with a shorter \(\Delta s\), at most BRANCH_SWITCH_MAX_TRIES times, after which the switch is abandoned and the path continues on the branch already traced. The index in an.increments of the state from which the path was switched is appended to an.branch_points. The step after the switch starts along the tangent in the direction of the switching increment: with the inner-product rule for method='riks' and 'crisfield'; with 'displacement_control', the sign of \(\Delta s\) for which the predictor makes a positive inner product with the switching increment; with 'generalized_displacement_control', the same sign of the load increment, \(GSP_k = \langle \delta c_q^{(1)}, \delta c_q^{(1)} \rangle_W / \langle \delta c_q^{(k)}, \delta c_q^{(k)} \rangle_W\), i.e. \(\delta c_q^{(k-1)}\) replaced by \(\delta c_q^{(k)}\), and the corrections orthogonal to \(\delta c_q^{(k)}\), as at the first step, since the tangent before the switch is that of the other branch. The detection by the count of negative pivots and the switch along the critical eigenvector follow Wagner and Wriggers [7], whose algorithm (p. 106, step (c)) detects a stability point when the number of negative diagonal elements of the factorization changes, and whose classification by the current stiffness parameter, zero at the limit points and not at the bifurcations, their Eq. (22), is replaced here by the direction of the load along the tangents. In [7] the critical point is first located by bisection, the eigenvector is added to the displacements of the last converged state, \(\bar{\boldsymbol{v}} + \xi \boldsymbol{\phi}/\|\boldsymbol{\phi}\|\) with \(|\xi| = \|\bar{\boldsymbol{v}}\|/\tau\) and \(\tau\) of the order of 100, Eqs. (24) and (25), and an equilibrium state is then found by the arc-length method, repeated with another \(\tau\), up to about three times, when it falls back on the primary branch. Here the critical point is not located, the perturbation is the predictor of an arc-length step from the first state beyond the critical point, with the length of the arc-length steps instead of \(|\xi|\), and the shorter steps of the repeated attempts play the role of the other values of \(\tau\).

With an.restart_state, the dictionary an.path_state saved after a converged step of an analysis with the same method, the analysis continues from that step, with its arc length and the directions of its predictors, as the analysis that saved it would have, appending the next converged increments to the empty an.increments and an.cs. an.path_state holds arrays and numbers, e.g. to be saved with pickle and the analysis restarted with other options or after the end of a job.

When the load factor 1.0 is exceeded within a converged step, from \(\lambda < 1\) to \(\lambda + \Delta\lambda > 1\), the step is replaced by Newton-Raphson iterations at exactly \(\lambda = 1\), see newton_raphson._NR_iterations(), from the predictor \(t \Delta \boldsymbol{c}\) with \(t = (1 - \lambda)/\Delta\lambda\). When these iterations do not converge, the step is repeated with a shorter arc length.

When the load depends on the configuration, i.e. calc_fext accepts the keyword argument c, \(\boldsymbol{F}_{ext}\) in \(\delta \boldsymbol{c}_q\) and in the tangent predictor is the load vector of the current configuration calc_fext(inc=1., c=c), minus the derivative \(\partial \boldsymbol{R}/\partial \lambda\) of the residual, whereas \(\boldsymbol{u}_{ref}\) remains the linear solution of the reference vector. The callables that accept inc receive the current load factor, see callbacks.

When calc_fext also accepts the keyword argument lbd, e.g. for prescribed displacements imposed exactly, it receives the load factor together with the configuration, \(\boldsymbol{F}_{ext}\) above being then \(-\partial \boldsymbol{R}/\partial \lambda\) at the current \(c\) and \(\lambda\); \(\boldsymbol{u}_{ref}\) is the linear solution of \(-\partial \boldsymbol{R}/\partial \lambda\) at \(c = 0\), \(\lambda = 0\), since the reference vector may be null, e.g. under a prescribed displacement alone; and \(-\lambda \, \partial \boldsymbol{R}/\partial \lambda\) is an additional reference of the relative convergence criterion. \(-\partial \boldsymbol{R}/\partial \lambda\) is evaluated together with the tangent stiffness matrix, at the same iterate, i.e. once per update of the tangent stiffness matrix, and used in the bordered solution; the last one evaluated is the reference of the next convergence checks. At the converged state it is kept with the tangent stiffness matrix when the step is repeated with a shorter arc length.

With an.constraints and an.constraint_values, the constraints \(\boldsymbol{C}^\top \boldsymbol{c} = \lambda \boldsymbol{d}\), e.g. a displacement prescribed in proportion to the load factor, linear in \((\boldsymbol{c}, \lambda)\): the load solutions \(\delta \boldsymbol{c}_q\) (and \(\boldsymbol{u}_{ref}\)) are those of the system bordered by \(\boldsymbol{d}\), \(\boldsymbol{C}^\top \delta \boldsymbol{c}_q = \boldsymbol{d}\), and the residual solutions \(\delta \boldsymbol{c}_R\) those of the null border, see static._factorize(), such that every predictor, tangent, of the displacement controls or the jump to \(\lambda = 1\), and every correction \(\delta \boldsymbol{c}_R + \delta\lambda \, \delta \boldsymbol{c}_q\) satisfies \(\boldsymbol{C}^\top \Delta \boldsymbol{c} = \Delta\lambda \boldsymbol{d}\); the eigenvector predictor of the branch switching, at a constant load, has the null border. Under a prescribed displacement alone, \(\boldsymbol{F}_{ext} = \boldsymbol{0}\), \(\boldsymbol{u}_{ref}\) is the displacement of the linear solution, not null, and the reactions of the constraints, in the internal force vector, are the reference of the convergence criterion.

References

[1] (1,2)

Riks, E., An incremental approach to the solution of snapping and buckling problems, International Journal of Solids and Structures, 15(7), 529-551, 1979. https://doi.org/10.1016/0020-7683(79)90081-7

[2] (1,2,3)

Crisfield, M.A., A fast incremental/iterative solution procedure that handles “snap-through”, Computers & Structures, 13(1-3), 55-62, 1981. https://doi.org/10.1016/0045-7949(81)90108-5

[3] (1,2)

Ramm, E., Strategies for tracing the nonlinear response near limit points, in: Wunderlich, W., Stein, E., Bathe, K.J. (eds.), Nonlinear Finite Element Analysis in Structural Mechanics, Springer, 63-89, 1981. https://doi.org/10.1007/978-3-642-81589-8_5

[4] (1,2,3)

Crisfield, M.A., Non-linear Finite Element Analysis of Solids and Structures, Volume 1: Essentials, John Wiley & Sons, Chichester, 1991.

[5]

de Borst, R., Computation of post-bifurcation and post-failure behavior of strain-softening solids, Computers & Structures, 25(2), 211-224, 1987. https://doi.org/10.1016/0045-7949(87)90144-1

[6] (1,2,3,4)

Yang, Y.-B., Shieh, M.-S., Solution method for nonlinear problems with multiple critical points, AIAA Journal, 28(12), 2110-2116, 1990. https://doi.org/10.2514/3.10529

[7] (1,2)

Wagner, W., Wriggers, P., A simple method for the calculation of postcritical branches, Engineering Computations, 5(2), 103-109, 1988. https://doi.org/10.1108/eb023727

Riks and Crisfield methods

structsolve.arc_length_riks._solver_arc_length_riks(an, silent=False)

Arc-Length solver using the Riks method

Runs arc_length._solver_arc_length() with method='riks'.

Notes

The arc-length constraint is linearized at every iteration (updated normal plane), see arc_length._solver_arc_length().

structsolve.arc_length_crisfield._solver_arc_length_crisfield(an, silent=False)

Arc-Length solver using Crisfield`s method

Runs arc_length._solver_arc_length() with method='crisfield'.

Notes

The quadratic (spherical) arc-length constraint is solved exactly at every iteration, see arc_length._solver_arc_length().

Critical points

structsolve.arc_length._critical_eigenvector(solve_kT, size, num_neg, num_new)

Eigenvector of the tangent stiffness matrix whose eigenvalue just became negative

Subspace iteration with the inverse of the factorization solve_kT, see static._factorize(), and the Rayleigh-Ritz projection of the inverse on the orthonormalized subspace, which gives the eigenvalues of \(\boldsymbol{K}_T\) closest to zero and their eigenvectors in the Euclidean inner product, see _eigs_near_zero(). The subspace has num_neg + EIGVEC_NUM_GUARD vectors, num_neg being the number of negative eigenvalues of \(\boldsymbol{K}_T\), and is iterated until the Ritz values change by less than EIGVEC_TOL, at most EIGVEC_MAX_ITER times. When fewer than num_neg of the Ritz values are negative, i.e. positive eigenvalues are closer to zero than some negative ones, e.g. of weakly stiff terms of a Ritz model, the size of the subspace is doubled and the iteration repeated, up to EIGVEC_MAX_SIZE vectors. Without a count of negative eigenvalues ('splu'), num_neg is 1 when the determinant is negative and 0 otherwise. With an.constraints the factorization is the bordered one, whose solutions satisfy the constraints, and the eigenpairs are those of \(\boldsymbol{K}_T\) restricted to them.

Of the negative eigenvalues found, the num_new closest to zero are taken as those that became negative in the last step, and the most negative of them is the critical one, the first that became negative. Without negative eigenvalue, the one closest to zero is taken. With num_new < 0, -num_new eigenvalues became positive in the last step, the subspace has -num_new more vectors, and the critical one is the largest of the -num_new smallest positive eigenvalues, the first that became positive.

Returns:

phiarray-like

The eigenvector, of unit Euclidean norm.

mufloat

Its eigenvalue.

structsolve.arc_length._eigs_near_zero(solve_kT, size, p, start=None)

The p eigenvalues of \(\boldsymbol{K}_T\) closest to zero and their eigenvectors, by the subspace iteration of _critical_eigenvector(), with at most EIGVEC_MAX_ITER iterations, from the vectors start, columns, completed by random vectors. With 'pardiso' the solutions of the factorization are used without the check of their residual and GMRES, see static._factorize(): the tangent stiffness matrix near a critical point is almost singular, the GMRES of its solutions took about one second each for the DOE09 cylinder (57000 degrees of freedom), against 0.03 s for those of the factorization alone, and the subspace iteration converges to the eigenvectors closest to zero without them.

Returns:

musarray-like

The eigenvalues, not sorted.

vecsarray-like

The eigenvectors, the columns, orthonormal.

Constants

structsolve.arc_length.PSI = 1.0

Scaling factor of the load-factor term in the arc-length constraint

structsolve.arc_length.DESIRED_NUM_ITER = 4

Desired number of corrections per step, used to adapt the arc length

structsolve.arc_length.MAX_GROWTH_FACTOR = 1.5

Maximum factor to increase the arc length after a converged step

structsolve.arc_length.MIN_GROWTH_FACTOR = 0.5

Minimum factor to reduce the arc length after a converged step

structsolve.arc_length.MAX_NUM_STEPS = 1000

Maximum number of converged steps of an arc-length analysis

structsolve.arc_length.BRANCH_SWITCH_MIN_COS = 0.5

Minimum cosine between the converged increment of a branch-switching step and its predictor, the critical eigenvector

structsolve.arc_length.BRANCH_SWITCH_MAX_TRIES = 3

Number of attempts of a branch-switching step, each with the arc length multiplied by INC_CUT_FACTOR, before the switch is abandoned

structsolve.arc_length.EIGVEC_NUM_GUARD = 2

Number of eigenpairs computed besides the negative eigenvalues by the subspace iteration of the critical eigenvector

structsolve.arc_length.EIGVEC_MAX_ITER = 100

Maximum number of iterations of the subspace iteration

structsolve.arc_length.EIGVEC_TOL = 1e-10

Relative tolerance of the eigenvalues of the subspace iteration

structsolve.arc_length.EIGVEC_MAX_SIZE = 64

Largest subspace of _critical_eigenvector(), whose size is doubled until it holds the negative eigenvalues