Static analyses of a model given by callables (pyfe3d.solver.analysis)

The class Analysis receives the callables that compute the force vectors and the stiffness matrices of a model and runs a linear static analysis, or a non-linear static analysis with the Newton-Raphson solver of newton_raphson or the arc-length solvers of arc_length, with Crisfield’s convergence criterion by default, see convergence. The optional keyword arguments of the callables are described in callbacks. With bu, the callables work with the vectors and matrices of all the degrees-of-freedom, as assembled by the elements of this package, and FreeDofs restricts them to the free ones and imposes the prescribed displacements uk.

class pyfe3d.solver.analysis.Analysis(calc_fext=None, calc_fint=None, calc_kC=None, calc_kG=None, bu=None, uk=None)

Linear and non-linear static analyses of a model given by callables

The structural model is defined by the callables passed to the constructor, which return the force vectors and the stiffness matrices. The static analysis, linear or non-linear, is run with static(), whose solution is stored in the attributes increments and cs.

The non-linear analyses find the solution vectors c for which the internal force vector calc_fint(c) balances the external force vector calc_fext(inc=1.) scaled by the load factor, using the tangent stiffness matrix calc_kC(c, NLgeom=True) + calc_kG(c, NLgeom=True). With the elements of this package, calc_kC normally returns KC0 + KCNL and calc_kG returns KG, both updated at c, and calc_fint the vector of update_fint(..., nonlinear=1).

The analysis parameters are attributes of this class, described in the following tables together with their default values:

Non-Linear Algorithm	Description
NL_method	str, 'NR' (default) for the Newton-Raphson method, see newton_raphson._solver_NR(), or 'arc_length_riks' and 'arc_length_crisfield' for the arc-length methods, see arc_length._solver_arc_length()
line_search	bool, activates a safeguarding line-search, for the Newton-Raphson method only, see newton_raphson._NR_iterations(). Default is False
max_iter_line_search	int, maximum number of iterations of the line-search. Default is 20
modified_NR	bool, activates the modified Newton-Raphson method, where the tangent stiffness matrix is not updated at every iteration. Default is False, i.e. full Newton-Raphson
compute_every_n	int, if modified_NR=True, the tangent stiffness matrix is updated every \(n\) iterations. Default is 6
kT_initial_state	bool, if modified_NR=True, tells if the tangent stiffness matrix should be calculated already at the first iteration of the analysis, which is required for example when initial imperfections take place. Otherwise the linear constitutive stiffness matrix is used. Default is True

Incrementation	Description
initialInc	initial load increment. In the arc-length methods it defines the initial arc-length increment, corresponding to a load factor increment of initialInc along the initial tangent. Default is 0.1
minInc	minimum increment; the analysis stops when the load increment (Newton-Raphson) or the arc-length increment (arc-length methods) becomes smaller than minInc. Default is 1.e-4
maxInc	maximum load increment, or maximum arc-length increment for the arc-length methods. Default is 1.
maxArcLength	maximum cumulative arc length covered by the arc-length methods, dimensionless. Default is 18

Convergence Criteria	Description
convergence	str, 'crisfield' (default) for Crisfield’s criterion, the residual scaled with the diagonal of the linear stiffness matrix over the free degrees-of-freedom, or 'relative' for the Euclidean norm of the residual, see convergence
relTOL	tolerance \(\epsilon\) of the relative force imbalance: Crisfield’s \(\zeta \le\) relTOL, or, with 'relative', the norm of the residual force vector smaller than relTOL times the largest norm between the external and internal force vectors. Not used if None. Default is 1.e-9, the tight tolerance of Crisfield’s criterion with a consistent tangent, see convergence; the same default applies with 'relative'
absTOL	the convergence is also achieved when the maximum absolute residual force is smaller than this value, which depends on the units of the model. Not used if None. Default is None
maxNumIter	maximum number of iterations (corrections) of a step; if achieved the increment is reduced. Default is 30
too_slow_TOL	a step is considered too slow when the smallest residual norm is not reduced by this fraction over the last iterations; the increment is then reduced. Default is 0.005

Parameters:

calc_fextcallable, optional

calc_fext(inc=1., silent=False), must return a 1-D array with the external force vector. Required for linear and non-linear static analyses. The non-linear solvers call it with inc=1. and scale the returned vector by the load factor. It may declare the optional keyword arguments c, the current configuration, and lbd, the current load factor, and must then return the load vector of that configuration and load factor for a unit load factor, inc=1., see callbacks. The linear static analysis never passes lbd.

calc_fintcallable, optional

calc_fint(c, silent=False), must return a 1-D array with the internal force vector for the solution vector c. Required for non-linear analyses. It may declare the optional keyword argument inc, which then receives the current load factor, see callbacks.

calc_kCcallable, optional

calc_kC(c=None, NLgeom=False, silent=False), must return a sparse matrix with the constitutive stiffness matrix. With c=None and NLgeom=False it must return the linear stiffness matrix, which also gives the scaling of Crisfield’s criterion, and with NLgeom=True the constitutive part of the tangent stiffness matrix at c. Required for linear and non-linear static analyses. It may declare the optional keyword argument inc, which then receives the current load factor in the calls with NLgeom=True, see callbacks.

calc_kGcallable, optional

calc_kG(c=None, NLgeom=False, silent=False), must return a sparse matrix with the geometric stiffness matrix at c. It is called with NLgeom=True. Required for non-linear analyses. It may declare the optional keyword argument inc, which then receives the current load factor, see callbacks. The sum of calc_kC and calc_kG may be unsymmetric.

buarray_like or None, optional

Boolean mask of the free degrees-of-freedom, 1-D array of size N. With bu, the callables receive and return the vectors and matrices of all the N degrees-of-freedom, e.g. assembled with the elements of this package, and the solvers work with the free ones only, see FreeDofs. With None, the callables work with the vectors and matrices of the solver, i.e. of the free degrees-of-freedom only, the constrained ones being partitioned out: a null row or column of the stiffness matrices raises an error, see check_null_rows_cols(). Default None.

ukarray_like or None, optional

Prescribed displacements for a unit load factor, 1-D array of size N, used with bu: the displacements of the constrained degrees-of-freedom are lbd*uk[~bu], with lbd the load factor. None for null prescribed displacements. Default None.

Attributes:

incrementslist

Load factors of the converged increments, filled by static().

cslist

Solution vectors of the converged increments, filled by static(), of size N with bu.

last_analysisstr

Type of the last analysis run, 'static' after static().

Methods

static([NLgeom, silent])

Linear or non-linear static analysis

Notes

For non-linear analyses, the internal force vector \(\boldsymbol{F}_{int}(\boldsymbol{c})\) must be in equilibrium with the external force vector scaled by the load factor \(\lambda\):

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

and the tangent stiffness matrix \(\boldsymbol{K}_T = \boldsymbol{K}_C + \boldsymbol{K}_G\), given by the callables calc_kC and calc_kG, should be the exact derivative of \(\boldsymbol{F}_{int}\) with respect to \(\boldsymbol{c}\), \(\boldsymbol{K}_T = \partial \boldsymbol{R}/\partial \boldsymbol{c}\), such that the Newton-Raphson iterations, \(\boldsymbol{K}_T \delta \boldsymbol{c} = -\boldsymbol{R}\), converge quadratically ([1], Section 1.2.2, p. 9). The convergence is checked by default with Crisfield’s criterion, see convergence.

The line-search tries the full step first and only reduces it when it fails a sufficient-decrease test on the residual norm, see newton_raphson._NR_iterations().

The arc length is dimensionless, with displacements scaled by the linear solution for a load factor of 1, such that in the linear regime an arc length of about sqrt(2) corresponds to a load factor increment of 1, see arc_length._solver_arc_length().

Configuration-dependent loads

Configuration-dependent loads, e.g. follower pressures [2], are supported through optional keyword arguments, passed only to the callables that declare them explicitly, see callbacks:

• calc_fint(c, inc=lbd), calc_kC(c, NLgeom=True, inc=lbd) and calc_kG(c, NLgeom=True, inc=lbd) receive the current load factor, such that calc_fint may contain the configuration-dependent part of the load, scaled by inc, and calc_kC + calc_kG its derivative, the load stiffness matrix.

• calc_fext(inc=1., c=c) must return the load vector of the configuration c for a unit load factor, minus the derivative of the residual with respect to the load factor, \(-\partial \boldsymbol{R}/\partial \lambda\), used by the arc-length solvers.

• calc_fext(inc=1., c=c, lbd=lbd) also receives the current load factor when it declares lbd, for a derivative \(\partial \boldsymbol{R}/\partial \lambda\) that depends on the load factor, e.g. prescribed displacements imposed exactly as linear constraints, where the reference vector calc_fext(inc=1.) may be null. The solvers then use \(-\partial \boldsymbol{R}/\partial \lambda\) for the displacement scale of the arc-length methods, for the predictor of the first Newton-Raphson increment and as a reference of the convergence criteria. The linear static analysis never passes lbd.

The tangent stiffness matrix may then be unsymmetric, which the solvers support, since they use LU factorizations. With bu and uk, the prescribed displacements are imposed with lbd, see FreeDofs.

References

[1]

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

[2]

Schweizerhof, K. and Ramm, E., Displacement dependent pressure loads in nonlinear finite element analyses, Computers & Structures, 18(6), 1099-1114, 1984.

static(NLgeom=False, silent=False)

Linear or non-linear static analysis

The linear analysis solves the system of the linear stiffness matrix calc_kC() and the external force vector calc_fext() using solve(). When the load depends on the configuration, i.e. calc_fext declares the keyword argument c, see callbacks, the load is evaluated in the current configuration and corrected by iterations that reuse one factorization of the linear stiffness matrix, starting from the solution of the reference load, until the convergence criterion is satisfied, see check_convergence(), or maxNumIter iterations are done; a warning is printed when they do not converge and the last iterate is returned. The non-linear analysis uses the solver selected by the NL_method attribute.

Parameters:

NLgeombool, optional

Whether a non-linear analysis is performed. Default False.

silentbool, optional

Whether the log messages are suppressed. Default False.

Returns:

incrementslist

Load factors of the converged increments, dimensionless. A linear analysis returns [1.]. Non-linear analyses finish at a load factor of exactly 1. unless they are stopped earlier, see the documentation of each solver.

cslist

Solution vectors of the converged increments, 1-D arrays, of size N with bu, including the prescribed displacements.

Raises:

ValueError

When NLgeom=True and the NL_method attribute is not 'NR', 'arc_length_riks' or 'arc_length_crisfield', or when uk is given without bu.

numpy.linalg.LinAlgError

When the linear stiffness matrix of the solver calc_kC(), with bu that of the free degrees-of-freedom, has a null row or column, see check_null_rows_cols(), or when it is exactly singular. A singular tangent stiffness matrix during the non-linear iterations does not raise, it stops the iterations, see the documentation of each solver.

Notes

The linear analysis solves \(\boldsymbol{K}_{C0} \boldsymbol{c} = \boldsymbol{F}_{ext}\); with bu and uk, \(\boldsymbol{K}_{{C0}_{ff}} \boldsymbol{c} = \boldsymbol{F}_{{ext}_f} - \boldsymbol{K}_{{C0}_{fk}} \bar{\boldsymbol{u}}_k\). When the load depends on the configuration, the structure is still geometrically linear but the load is evaluated in the current configuration, \(\boldsymbol{K}_{C0} \boldsymbol{c} = \boldsymbol{F}_{ext}(\boldsymbol{c})\), which is solved by the load correction iterations

\[\boldsymbol{c}_{k+1} = \boldsymbol{c}_k + \boldsymbol{K}_{C0}^{-1}(\boldsymbol{F}_{ext}(\boldsymbol{c}_k) - \boldsymbol{K}_{C0} \boldsymbol{c}_k)\]

starting from \(\boldsymbol{c}_0 = \boldsymbol{K}_{C0}^{-1} \boldsymbol{F}_{ext}\), the solution of the reference load calc_fext(), with \(\boldsymbol{F}_{ext}(\boldsymbol{c}_k)\) given by calc_fext(inc=1., c=c_k). This is a fixed-point iteration without the load stiffness in the iteration matrix. For a follower load whose force vector is affine in \(\boldsymbol{c}\), \(\boldsymbol{F}(\boldsymbol{c}) = \boldsymbol{F}_0 - \boldsymbol{K}_F \boldsymbol{c}\), this is the solution of \((\boldsymbol{K}_{C0} + \boldsymbol{K}_F) \boldsymbol{c} = \boldsymbol{F}_0\), the error being multiplied by \(-\boldsymbol{K}_{C0}^{-1} \boldsymbol{K}_F\) at each iteration, whose spectral radius is \(1/\min|\lambda|\) over the load multipliers \(\lambda\) of \((\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_F) \boldsymbol{u} = \boldsymbol{0}\), including the negative and complex ones; beyond the load of the smallest \(|\lambda|\) the iterations diverge.

class pyfe3d.solver.analysis.FreeDofs(calc_fext, calc_fint, calc_kC, calc_kG, bu, uk=None)

Callables of the free degrees-of-freedom, with prescribed displacements imposed exactly

Wraps the callables of the full model, which receive and return the vectors and matrices of all the N degrees-of-freedom, such that the solvers see the free degrees-of-freedom bu only. The displacements of the constrained degrees-of-freedom are lbd*uk[~bu], prescribed and proportional to the load factor lbd, or null when uk is None. Used by Analysis when bu is given.

The wrapped calc_fint, calc_kC and calc_kG always declare inc, which they need to build the full configuration, and pass it to the callables of the full model only when they declare it, see callbacks. The wrapped calc_fext declares c when the one of the full model declares it or when displacements are prescribed, and lbd when the one of the full model declares it or when displacements are prescribed.

Parameters:

calc_fext, calc_fint, calc_kC, calc_kGcallable

The callables of the full model, see Analysis.

buarray_like

Boolean mask of the free degrees-of-freedom, 1-D array of size N.

ukarray_like or None, optional

Prescribed displacements for a unit load factor, 1-D array of size N whose entries at the free degrees-of-freedom are ignored, of the dimension of the displacements and rotations. None for null prescribed displacements. Default None.

Attributes:

free, knownnp.ndarray

Indices of the free and of the constrained degrees-of-freedom.

prescribedbool

Whether a prescribed displacement is not null.

Methods

calc_fint(c[, silent, inc])	Internal force vector of the free degrees-of-freedom
calc_kC([c, NLgeom, silent, inc])	Constitutive stiffness matrix of the free degrees-of-freedom
calc_kG([c, NLgeom, silent, inc])	Geometric stiffness matrix of the free degrees-of-freedom
full(c, lbd)	Configuration of the full model
linear_fext([c])	Load vector of the linear static analysis, at a unit load factor

Raises:

ValueError

When bu is not a boolean array, or when uk does not have the size of bu.

Notes

With \(f\) the free and \(k\) the constrained degrees-of-freedom, the configuration of the full model is

\[\begin{split}\boldsymbol{u} = \boldsymbol{Z} \boldsymbol{c} + \lambda \boldsymbol{u}_k, \qquad \boldsymbol{Z} = \begin{bmatrix} \boldsymbol{I}_{ff} \\ \boldsymbol{0}_{kf} \end{bmatrix}\end{split}\]

the null-space reduction of the constraints \(\boldsymbol{u}_k = \lambda \bar{\boldsymbol{u}}_k\) imposed exactly, see callbacks. The solvers see the reduced residual \(\boldsymbol{R}_f = \boldsymbol{F}_{{int}_f}(\boldsymbol{u}) - \lambda \boldsymbol{F}_{{ext}_f}\), its Jacobian \(\boldsymbol{K}_{T_{ff}}\) and

\[-\frac{\partial \boldsymbol{R}_f}{\partial \lambda} = \boldsymbol{F}_{{ext}_f} - \boldsymbol{K}_{T_{fk}} \bar{\boldsymbol{u}}_k\]

with \(\boldsymbol{K}_T\) at the current configuration and load factor, reused from the last evaluation of the tangent stiffness matrix when it was done at the same configuration and load factor. The norms of the convergence criteria are thus over the free degrees-of-freedom only. The reaction forces are \(\boldsymbol{F}_{{int}_k}(\boldsymbol{u}) - \lambda \boldsymbol{F}_{{ext}_k}\), given by the callables of the full model at a converged configuration.

calc_fint(c, silent=True, inc=1.0)

Internal force vector of the free degrees-of-freedom

Parameters:

cnp.ndarray

Configuration of the free degrees-of-freedom.

silentbool, optional

Passed to the callable of the full model. Default True.

incfloat, optional

Load factor, dimensionless. Default 1..

Returns:

fintnp.ndarray

calc_fint of the full model at full(), restricted to the free degrees-of-freedom.

calc_kC(c=None, NLgeom=False, silent=True, inc=1.0)

Constitutive stiffness matrix of the free degrees-of-freedom

Parameters:

cnp.ndarray or None, optional

Configuration of the free degrees-of-freedom. None gives the linear stiffness matrix. Default None.

NLgeombool, optional

Passed to the callable of the full model; False gives the linear stiffness matrix. Default False.

silentbool, optional

Passed to the callable of the full model. Default True.

incfloat, optional

Load factor, dimensionless. Default 1..

Returns:

kCsparse_matrix

Free-free block of calc_kC of the full model.

calc_kG(c=None, NLgeom=False, silent=True, inc=1.0)

Geometric stiffness matrix of the free degrees-of-freedom

Parameters:

cnp.ndarray or None, optional

Configuration of the free degrees-of-freedom. Default None.

NLgeombool, optional

Passed to the callable of the full model. Default False.

silentbool, optional

Passed to the callable of the full model. Default True.

incfloat, optional

Load factor, dimensionless. Default 1..

Returns:

kGsparse_matrix

Free-free block of calc_kG of the full model.

config_dependent

whether the external load of the full model depends on the configuration

full(c, lbd)

Configuration of the full model

Parameters:

cnp.ndarray

Configuration of the free degrees-of-freedom.

lbdfloat

Load factor, dimensionless.

Returns:

unp.ndarray

Configuration of all the degrees-of-freedom, 1-D array of size N, with lbd*uk at the constrained ones.

linear_fext(c=None)

Load vector of the linear static analysis, at a unit load factor

Parameters:

cnp.ndarray or None, optional

Configuration of the free degrees-of-freedom, passed to the calc_fext of the full model when it declares c. Default None.

Returns:

fnp.ndarray

\(\boldsymbol{F}_{{ext}_f} - \boldsymbol{K}_{{C0}_{fk}} \bar{\boldsymbol{u}}_k\), with \(\boldsymbol{K}_{C0}\) the linear stiffness matrix of the full model.