Arc-length solvers (pyfe3d.solver.arc_length)

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

pyfe3d.solver.arc_length.DESIRED_NUM_ITER = 4

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

pyfe3d.solver.arc_length.MAX_GROWTH_FACTOR = 1.5

Maximum factor to increase the arc length after a converged step

pyfe3d.solver.arc_length.MAX_NUM_STEPS = 1000

Maximum number of converged steps of an arc-length analysis

pyfe3d.solver.arc_length.METHODS = ('riks', 'crisfield')

Names of the arc-length methods, see _solver_arc_length()

pyfe3d.solver.arc_length.MIN_GROWTH_FACTOR = 0.5

Minimum factor to reduce the arc length after a converged step

pyfe3d.solver.arc_length.PSI = 1.0

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

pyfe3d.solver.arc_length._solver_arc_length(an, pb, method, silent=False)

Arc-length solver

Traces the equilibrium path of the model 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. Used by Analysis.static() when NL_method is 'arc_length_riks' or '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 check_convergence() and newton_raphson._check_divergence(); 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 1.0, the last solution being computed at a load factor of exactly 1.0;

• 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.

Parameters:

anAnalysis

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

pbobject

Object with the callables, an itself or a FreeDofs.

methodstr

'riks' or 'crisfield'.

silentbool, optional

Whether the log messages are suppressed. Default False.

Raises:

ValueError

When method is not one of METHODS.

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 a uniform change of the units of the displacements; it is not independent of a change of the units of the rotations only. 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 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\).

Predictor

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 c_q, 1), (\delta c_q, 1) \rangle}\) at the origin, corresponds to a load factor increment of an.initialInc along the initial tangent.

Corrections

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\).

The step fails when \(\delta\lambda\) cannot be computed or is not finite, and when the converged increment points against the tangent 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.

The convergence criterion is the one of the Newton-Raphson solver, see check_convergence(), at the current load factor \(\lambda + \Delta\lambda\), the arc-length method changing only how the load factor is updated and not the residual. Both criteria also consider the external force vector of the increment, \(\Delta\lambda \boldsymbol{F}_{ext}\), as a reference.

Adaptation of the arc length

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.

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.

Configuration- and load-factor-dependent loads

When 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 \(\boldsymbol{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 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.

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.