Non-linear static analyses#

The geometrically non-linear static analyses are run with Analysis, which receives four callables of the model and finds the configurations \(\boldsymbol{u}\) in equilibrium with the external load scaled by the load factor \(\lambda\):

\[\boldsymbol{R}(\boldsymbol{u}, \lambda) = \boldsymbol{F}_{int}(\boldsymbol{u}) - \lambda \boldsymbol{F}_{ext} = \boldsymbol{0}\]

The elements give the terms of the residual and of its Jacobian, the solvers of pyfe3d.solver the increments and the iterations.

From the elements to the callables#

Callable	Returns	With the elements
`calc_fext(inc=1., silent=True)`	\(\boldsymbol{F}_{ext}\), scaled by `inc`	the nodal loads of the model
`calc_fint(c, silent=True)`	\(\boldsymbol{F}_{int}(\boldsymbol{c})\)	`update_fint(fint, prop, nonlinear=1)` of every element, after `update_probe_xe(x)` and `update_probe_ue(c)`
`calc_kC(c=None, NLgeom=False, silent=True)`	\(\boldsymbol{K}_{C0}\), or \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}(\boldsymbol{c}) + \boldsymbol{K}_{GNL}(\boldsymbol{c})\) with `NLgeom=True`	`update_KC0` once, `update_KCNL` at `c`
`calc_kG(c=None, NLgeom=False, silent=True)`	\(\boldsymbol{K}_G(\boldsymbol{c})\)	`update_KG` at `c`

calc_kC() + calc_kG() is then the exact Jacobian of calc_fint, \(\boldsymbol{K}_T = \partial \boldsymbol{R}/\partial \boldsymbol{u}\), and the full Newton-Raphson method converges quadratically. The class Model of tests/test_solver_nonlinear_static.py, see Non-linear static analysis with pyfe3d.solver, wraps a list of elements of one type into these four callables. calc_kC(), the linear stiffness matrix, also gives the scaling of the convergence criterion.

Boundary conditions and prescribed displacements#

With Analysis(..., bu=bu), the callables receive and return the vectors and matrices of all the N dofs, as assembled by the elements, and FreeDofs restricts them to the free dofs bu. The solution vectors an.cs have size N, with zeros at the constrained dofs.

With uk, the constrained dofs have the prescribed displacements \(\lambda \boldsymbol{u}_k\), imposed exactly, proportional to the load factor: uk[~bu] are the values at \(\lambda = 1\), and a load-free analysis, calc_fext returning zeros, is displacement controlled. The reaction at a prescribed dof is calc_fint(u) at that dof minus the load applied there, see test_prescribed_displacement_and_reaction in Non-linear static analysis with pyfe3d.solver.

Without bu, the callables return the vectors and matrices of the free dofs only. The solvers do not remove null rows and columns: a null row or column of the linear stiffness matrix, a dof without stiffness, raises a numpy.linalg.LinAlgError naming the dofs, see check_null_rows_cols().

Solvers#

an.NL_method selects the solver, and Analysis.static() with NLgeom=True runs it, returning the converged load factors an.increments and solutions an.cs:

'NR', the default: load control with the Newton-Raphson method, newton_raphson._solver_NR(), full or modified (an.modified_NR), see Crisfield (1991), Sections 1.2.2 and 1.2.3, with an optional line search (an.line_search). The first increment is an.initialInc; a step that fails is repeated with half the increment, and the increment grows after a step that converged without being cut, up to an.maxInc. Beyond a limit point of the load it may converge to another branch of the equilibrium path.
'arc_length_riks' and 'arc_length_crisfield': the arc-length methods following Riks (1979) and Crisfield (1981), arc_length._solver_arc_length(), in which the load factor is an unknown of each step, such that limit points, snap-through and snap-back are traced. The first arc length corresponds to the load factor increment an.initialInc along the initial tangent.

Both stop at the load factor 1.0, or earlier when the increment falls below an.minInc or, for the arc-length methods, when the cumulative arc length reaches an.maxArcLength. A run that did not finish is detected with:

increments, cs = an.static(NLgeom=True)
if not increments or increments[-1] != 1.:
    ...  # stopped at increments[-1], see the log messages

The linear static analysis, an.static(NLgeom=False), solves \(\boldsymbol{K}_{C0} \boldsymbol{u} = \boldsymbol{F}_{ext}\) with solve() and returns a single increment.

Convergence criterion#

The iterations of each increment stop when the residual is small according to an.convergence, see convergence:

'crisfield', the default: the residual scaled with the diagonal of the linear stiffness matrix of the free dofs, \(d_j = K_{{C0}_{jj}}\), the entries of the diagonal scaling vector \(\boldsymbol{d}\),

\[\zeta = \frac{\|\boldsymbol{R}\|_{\boldsymbol{d}}}{\max( \|\lambda \boldsymbol{F}_{ext}\|_{\boldsymbol{d}}, \|\boldsymbol{F}_{int}\|_{\boldsymbol{d}}, \ldots)} \le \text{relTOL}, \qquad \|v\|_{\boldsymbol{d}} = \sqrt{\sum_j v_j^2/|d_j|}\]

with an.relTOL = 1e-9 by default, see crisfield_zeta() for the references in the denominator. It does not depend on the units of each dof, e.g. on whether the rotations are in radians or in milliradians, and it ignores the support reactions;
'relative': the Euclidean norm of the residual over the largest Euclidean norm of the references, which adds the squares of forces and of moments.

an.absTOL also accepts the iterations whose largest absolute residual is smaller, in the units of force of the model. With both criteria, an iteration has also converged at the round-off floor of the residual,

\[\|\boldsymbol{R}\|_{\boldsymbol{d}} \le 3 \varepsilon \, \| |\boldsymbol{K}_{C0}| \, |\boldsymbol{u}| \|_{\boldsymbol{d}}\]

with \(\varepsilon\) the machine precision, \(|\boldsymbol{K}_{C0}|\) and \(|\boldsymbol{u}|\) the elementwise absolute values of the linear stiffness matrix and of the displacements of the free dofs, and the Euclidean norm with 'relative', see roundoff_floor(): the residual cannot be computed more accurately, and, under large self-equilibrated axial forces, e.g. of an arch near a null load, \(\zeta\) would otherwise stagnate above relTOL and the analysis stop, see test_fung_kaplan_beam_arch in Non-linear static analysis with pyfe3d.solver. The log messages mark the iterations converged at the round-off floor. With a consistent tangent stiffness matrix and full Newton-Raphson, \(\zeta\) falls quadratically, as verified in Non-linear static analysis with pyfe3d.solver; a slow decrease of \(\zeta\) reveals an inconsistent tangent stiffness matrix, e.g. a missing calc_kG, which changes the iterations but not the converged solution.

Configuration-dependent loads#

A load that depends on the configuration, e.g. a follower pressure or the aerodynamic load of the piston theory, is given through optional keyword arguments, passed only to the callables that declare them, see callbacks:

calc_fint(c, inc=...) receives the load factor and includes the configuration-dependent part of the load, \(-\lambda \boldsymbol{F}(\boldsymbol{c})\), and calc_kC(c, NLgeom=True, inc=...) its derivative, the load stiffness matrix \(\boldsymbol{K}_F = \lambda \boldsymbol{K}_{F_{ref}}\), with \(\boldsymbol{K}_{F_{ref}} = -\partial \boldsymbol{F}/\partial \boldsymbol{c}\), which is unsymmetric in general, see below; the solvers use LU factorizations;
calc_fext(inc=1., c=...) returns the load vector of the configuration c for a unit load factor, \(-\partial \boldsymbol{R}/\partial \lambda\), used by the arc-length methods.

The aerodynamic load of the piston theory, \(\lambda \beta \boldsymbol{K}_{A_\beta} \boldsymbol{u}\), with update_KA_beta of Quad4R or Quad4, is such a load, see test_piston_theory_load_unsymmetric_tangent in Non-linear static analysis with pyfe3d.solver. The keyword argument lbd of calc_fext gives a load vector that depends on the load factor, which is how FreeDofs imposes the prescribed displacements.

\(\boldsymbol{K}_F\) is symmetric, and the load conservative, only in particular cases. For a pressure normal to the deformed surface, Hibbitt (1979, p. 938) splits the load stiffness of a uniform pressure into a symmetric surface integral and an antisymmetric integral along the boundary of the loaded surface, which vanishes where the displacements are prescribed: a uniform pressure on a fully enclosed volume, or on a surface whose boundaries have prescribed displacements, gives a symmetric \(\boldsymbol{K}_F\), whereas free edges, where the pressure changes, and a non-uniform pressure give unsymmetric terms (p. 939). In Table 2 of Schweizerhof and Ramm (1984), the domain terms are symmetric for a uniform pressure, and for a non-uniform one only when the pressure field is fixed in space, and the boundary terms vanish under one of the conditions of their Sec. 3.2.1, Eqs. (32)-(35): the pressure is zero on the whole boundary; at least two displacement components are prescribed at the edge; the edge is restrained normal to a continuous surface, e.g. the undeformed surface, \(w = 0\) for a plate on the \(XY\) plane; or the boundary terms add up to zero. pyfe3d does not symmetrise \(\boldsymbol{K}_F\): the tangent stiffness matrix, with the matrix returned by calc_kC, is factorized with scipy.sparse.linalg.splu(), as in solve(), such that a consistent unsymmetric \(\boldsymbol{K}_F\) keeps the quadratic convergence of the Newton-Raphson method, as verified by test_newton_raphson_quadratic_with_follower_load of tests/test_solver_follower_loads.py.

Log messages#

The solvers print their progress, the increments and the \(\zeta\) of each iteration, with pyfe3d.solver.logger, unless silent=True is passed to Analysis.static(). The callables must accept silent.

References#

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
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
Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 1: Essentials, John Wiley & Sons, Chichester, 1991.
Hibbitt, H. D., Some follower forces and load stiffness, International Journal for Numerical Methods in Engineering, 14(6), 937-941, 1979. https://doi.org/10.1002/nme.1620140613
Schweizerhof, K. and Ramm, E., Displacement dependent pressure loads in nonlinear finite element analyses, Computers & Structures, 18(6), 1099-1114, 1984. https://doi.org/10.1016/0045-7949(84)90154-8