Static analyses (structsolve.analysis)#
- class structsolve.Analysis(calc_fext=None, calc_fint=None, calc_kC=None, calc_kG=None)#
Class that embodies all data required for linear/non-linear analysis
The structural model is defined by the callables passed to the constructor, which return the force vectors and stiffness matrices. The static analysis, linear or non-linear, is run with
static(), whose solution is stored in the attributesincrementsandcs.The non-linear analyses find the solution vectors
cfor which the internal force vectorcalc_fint(c)balances the external force vectorcalc_fext(inc=1.)scaled by the load factor, using the tangent stiffness matrixcalc_kC(c, NLgeom=True) + calc_kG(c, NLgeom=True).The analysis parameters are attributes of this class, described in the following tables together with their default values:
Non-Linear Algorithm
Description
NL_methodstr,'NR'(default) for the Newton-Raphson method, seenewton_raphson._solver_NR(), or'arc_length_riks'and'arc_length_crisfield'for the arc-length methods, or'displacement_control'for the indirect displacement control, or'generalized_displacement_control', seearc_length._solver_arc_length()control_vectorarray, the vector \(\boldsymbol{\psi}\) of the controlled quantity \(\boldsymbol{\psi}^\top \boldsymbol{c}\) ofNL_method='displacement_control', e.g. a unit vector that picks one degree of freedom, or a buckling mode. Default isNonecontrol_scalefloat, the reference value of the controlled quantity, which makes the step dimensionless. Default isNone, i.e. \(|\boldsymbol{\psi}^\top \boldsymbol{u}_{ref}|\), its value in the linear solution for a load factor of 1line_searchbool, activates a safeguarding line-search, for the Newton-Raphson method only, seenewton_raphson._NR_iterations(). Default isFalsemax_iter_line_searchint, maximum number of iterations of the line-search. Default is20arc_length_line_searchbool, activates a backtracking line search on the corrections of the arc-length methods, seearc_length._solver_arc_length(). Default isFalsebranch_switchingint, number of bifurcation points at which the arc-length and the displacement control methods switch to the bifurcated branch, along the critical eigenvector, where the number of negative eigenvalues increases or decreases, which requiresinertia, seearc_length._solver_arc_length(). Default is0arc_length_weightsarray, non-negative weights of the degrees of freedom in the norm of the arc-length constraint, e.g. 1 for the translations and 0 for their derivatives or the rotations, seearc_length._solver_arc_length(). Default isNone, all equal to 1early_divergencebool, stops a diverging step earlier, such that its increment is cut before the iterations are spent, seenewton_raphson._check_divergence(). Default isFalsemodified_NRbool, activates the modified Newton-Raphson method, where the tangent stiffness matrix is not updated at every iteration. Default isFalse, i.e. full Newton-Raphsoncompute_every_nint, ifmodified_NR=True, the tangent stiffness matrix is updated at every \(n\) iterations. Default is6kT_initial_statebool, ifmodified_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 isTruelinear_solverstr, factorization of the tangent stiffness matrix of the non-linear solvers,'splu'(default), the sparse LU of SuperLU, or'pardiso', the symmetric indefinite factorization of PARDISO, which requires oneMKL, seepardiso, and gives the number of negative eigenvalues, seestatic._factorize()inertiabool, the tangent stiffness matrix is factorized at every converged increment of the non-linear analyses, and its number of negative eigenvalues ('pardiso') and the sign of its determinant are recorded, such that limit points, bifurcations crossed and unstable branches are visible. Default isFalseinertia_num_eigvalsint\(k\), withinertia: when the count of PARDISO is not available, e.g. for an unsymmetric tangent stiffness matrix, the number of its negative real eigenvalues among the \(k\) of smallest magnitude, scaled by the diagonal of the linear stiffness matrix, by shift-invert with ARPACK and the LU factorization of the converged state, seenewton_raphson._eigvals_inertia(). Default isNone, not computedmax_inertia_jumpint, a converged increment whose number of negative eigenvalues changed by more than this value is repeated with a smaller increment, which requiresinertiaandlinear_solver='pardiso'orinertia_num_eigvals. Default isNoneconstraintsarrayor sparse matrix \(\boldsymbol{C}\), of shape(n, m), a 1-D array for one constraint: the linear constraints \(\boldsymbol{C}^\top \boldsymbol{c} = \boldsymbol{0}\) of the solution, e.g. an inertia relief or multipoint constraints, imposed by bordering every factorization of the stiffness matrices, seestatic._factorize(), in the linear and the non-linear analyses. Default isNoneconstraint_valuesarray\(\boldsymbol{d}\), of shape(m,), the values of the constraints, \(\boldsymbol{C}^\top \boldsymbol{c} = \lambda \boldsymbol{d}\) along the non-linear path and \(\boldsymbol{C}^\top \boldsymbol{c} = \boldsymbol{d}\) in the linear analysis, e.g. a displacement or a rotation prescribed in proportion to the load factor, with or without external load, which requiresconstraints. Default isNone, all nullstop_criterioncallable,stop_criterion(an)is called after every converged increment of the non-linear analyses, which stop when it returnsTrue, e.g. a number of increments beyond the first limit point or a displacement measure, seestopping. Default isNonefinal_load_factorfloat, the load factor at which the non-linear analyses finish, exactly, or, for the arc-length and displacement control methods,None: the path is traced beyond any load factor until another stop, e.g.stop_criterionormaxArcLength. Default is1.restart_statedict, thepath_stateof a previous analysis with the sameNL_method, from whose last converged increment the non-linear analysis continues. Default isNoneIncrementation
Description
initialIncinitial load increment. In the arc-length methods it defines the initial arc-length increment, corresponding to a load factor increment of
initialIncalong the initial tangent. Default is0.1minIncminimum increment; the analysis stops when the load increment (Newton-Raphson) or the arc-length increment (arc-length methods) becomes smaller than
minInc. Default is1.e-4maxIncmaximum load increment, or maximum arc-length increment for the arc-length methods. Default is
1.maxArcLengthmaximum cumulative arc length covered by the arc-length methods, dimensionless. Default is
18Convergence Criteria
Description
convergencestr,'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, seeconvergencerelTOLtolerance \(\epsilon\) of the relative force imbalance: Crisfield’s \(\zeta \le\)
relTOL, or, with'relative', the norm of the residual force vector smaller thanrelTOLtimes the largest norm between the external and internal force vectors. Not used ifNone. Default is1.e-9, the tight tolerance of Crisfield’s criterion with a consistent tangent; until version 0.8.0 the relative criterion was the only one, with the default1.e-6absTOLthe 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 isNonemaxNumItermaximum number of iterations (corrections) of a step; if achieved the increment is reduced. Default is
30too_slow_TOLa 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 withinc=1.and scale the returned vector by the load factor. It may declare the optional keyword argumentsc, the current configuration, andlbd, the current load factor, and must then return the load vector of that configuration and load factor for a unit load factor,inc=1., seecallbacks. The linear static analysis never passeslbd.- calc_fintcallable, optional
calc_fint(c, silent=False), must return a 1-D array with the internal force vector for the solution vectorc. Required for non-linear analyses. It may declare the optional keyword argumentinc, which then receives the current load factor, seecallbacks.- calc_kCcallable, optional
calc_kC(c=None, NLgeom=False, silent=False), must return a sparse matrix with the constitutive stiffness matrix. Withc=NoneandNLgeom=Falseit must return the linear stiffness matrix, and withNLgeom=Truethe constitutive part of the tangent stiffness matrix atc. Required for linear and non-linear static analyses. It may declare the optional keyword argumentinc, which then receives the current load factor in the calls withNLgeom=True, seecallbacks.- calc_kGcallable, optional
calc_kG(c=None, NLgeom=False, silent=False), must return a sparse matrix with the geometric stiffness matrix atc. It is called withNLgeom=True. Required for non-linear analyses. It may declare the optional keyword argumentinc, which then receives the current load factor, seecallbacks. The sum ofcalc_kCandcalc_kGmay be unsymmetric.
- Attributes:
- incrementslist
Load factors of the converged increments, filled by
static().- cslist
Solution vectors of the converged increments, filled by
static().- num_negative_eigvalslist
With
inertia=True, for each converged increment of a non-linear analysis, the number of negative eigenvalues of the tangent stiffness matrix on the degrees of freedom whose row or column is not null,Nonewithlinear_solver='splu', when PARDISO fell back to SuperLU for an unsymmetric matrix, or for a singular matrix, unlessinertia_num_eigvalsis given: the number of negative real eigenvalues among theinertia_num_eigvalsof smallest magnitude is then recorded for the SuperLU factorizations. Empty otherwise.- det_signslist
With
inertia=True, for each converged increment of a non-linear analysis, the sign of the determinant of the same matrix,1or-1,0when it is singular. Empty otherwise.- num_perturbed_pivotslist
With
inertia=Trueandlinear_solver='pardiso', the number of pivots perturbed by PARDISO at each converged increment, which makes the count of negative eigenvalues uncertain when not zero;Noneotherwise.- branch_pointslist
With
branch_switching, the indices inincrementsof the converged increments from which the path was switched to a bifurcated branch.- path_statedict or None
After each converged increment of a non-linear analysis, the state from which it can be continued, see
restart_state: the solution, the load factor, the next increment and the directions of the predictors.- last_analysisstr
Type of the last analysis run,
'static'afterstatic().
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_kCandcalc_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. The convergence is checked by default with Crisfield’s criterion, the residual and the reference force vectors scaled with the diagonal of the linear stiffness matrixcalc_kC()over the free degrees of freedom, such that it does not depend on the units of each degree of freedom, e.g. translations and rotations, seeconvergence;convergence = 'relative'keeps the Euclidean norm of the versions up to 0.8.0.The line-search tries the full step first and only reduces it when it fails a sufficient-decrease test on the norm of the residual of the convergence criterion, 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, seearc_length._solver_arc_length().Configuration-dependent loads, e.g. follower pressures, 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)andcalc_kG(c, NLgeom=True, inc=lbd)receive the current load factor, such thatcalc_fintmay contain the configuration-dependent part of the load, scaled byinc, andcalc_kC + calc_kGits derivative.calc_fext(inc=1., c=c)must return the load vector of the configurationcfor 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 declareslbd, 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 vectorcalc_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 criterion, seecallbacks. The linear static analysis never passeslbd.
The tangent stiffness matrix may then be unsymmetric, which the solvers support, since they use LU factorizations.
- static(NLgeom=False, silent=False)#
General solver for static analyses
The linear analysis solves the system of the linear stiffness matrix
calc_kC()and the external force vectorcalc_fext()usingsolve(). When the load depends on the configuration, i.e.calc_fextdeclares the keyword argumentc, seecallbacks, 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 criteriarelTOLandabsTOLare satisfied, seenewton_raphson._check_convergence(), ormaxNumIteriterations 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 theNL_methodattribute.- Parameters:
- NLgeombool, optional
Flag to indicate whether a linear or a non-linear analysis is to be performed.
- silentbool, optional
A boolean to tell whether the log messages should be printed.
- Returns:
- incrementslist
Load factors of the converged increments. A linear analysis returns
[1.]. Non-linear analyses finish at a load factor of exactly1.unless they are stopped earlier, see the documentation of each solver.- cslist
Solution vectors of the converged increments.
- Raises:
- ValueError
When
NLgeom=Trueand theNL_methodattribute is not'NR','arc_length_riks','arc_length_crisfield','displacement_control'or'generalized_displacement_control'.
Notes
The linear analysis solves \(\boldsymbol{K}_C \boldsymbol{c} = \boldsymbol{F}_{ext}\). When the load depends on the configuration, the structure is still geometrically linear but the load is evaluated in the current configuration, \(\boldsymbol{K}_C \boldsymbol{c} = \boldsymbol{F}_{ext}(\boldsymbol{c})\), which is solved by the load correction iterations
\[\boldsymbol{c}_{k+1} = \boldsymbol{c}_k + \boldsymbol{K}_C^{-1}(\boldsymbol{F}_{ext}(\boldsymbol{c}_k) - \boldsymbol{K}_C \boldsymbol{c}_k)\]starting from \(\boldsymbol{c}_0 = \boldsymbol{K}_C^{-1} \boldsymbol{F}_{ext}\), the solution of the reference load
calc_fext(), with \(\boldsymbol{F}_{ext}(\boldsymbol{c}_k)\) given bycalc_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 \(c\), \(\boldsymbol{F}(\boldsymbol{c}) = \boldsymbol{F}_0 - \boldsymbol{K}_f \boldsymbol{c}\), this is the solution of \((\boldsymbol{K}_C + \boldsymbol{K}_f) \boldsymbol{c} = \boldsymbol{F}_0\), with the convergence rate of the spectral radius of \(\boldsymbol{K}_C^{-1} \boldsymbol{K}_f\), i.e. of the ratio of the load to the lowest critical load of the linearized problem; beyond it the iterations diverge. The error is multiplied by \(-\boldsymbol{K}_C^{-1} \boldsymbol{K}_f\) at each iteration, whose spectral radius is \(1/\min|\lambda|\) over the load multipliers \(\lambda\) of \((\boldsymbol{K}_C + \lambda \boldsymbol{K}_f) \boldsymbol{u} = \boldsymbol{0}\), including the negative and complex ones.