2. Quasi-static nonlinear study of the structure#
This step is justified if the structure has strong nonlinearities, which Euler’s analysis cannot take into account. The operator for solving nonlinear problems in quasistatics is called STAT_NON_LINE [bib7].
These nonlinearities can be linked to the material, which may have elastoplastic behavior [bib8], as in the example that follows. Taking contact, or even friction, into account is another source of non-linearities. We can also cite the case of follower loadings, such as pressure ([bib1] and [bib2] for solid shell elements), which require a non-linear approach in Code_Aster.
Two types of non-linear stability analyses can be distinguished, which can be combined.
On the one hand, the generalization of the Euler analysis presented above: we will monitor critical loads and associated modes during the incremental nonlinear calculation. This results in an Euler-type analysis on the updated stiffness matrices. This type of analysis is often done on a structure with no initial defect.
On the other hand, we can take into account defects introduced on the perfect model, in order to « force » the bifurcation into solution and to do branch monitoring to analyze the post-critical response.
Obviously, this post-critical solution monitoring can be initiated by the eigenvalue analysis on updated stiffness matrices, in particular in order to properly detect the bifurcation and to define the defect then introduced based on the observed buckling mode.
2.1. Stability analysis on updated stiffness matrices#
Whether in quasistatics (operator STAT_NON_LINE) or in dynamics (operator DYNA_NON_LINE), Code_Aster allows you to conduct incremental stability analyses in the sense of buckling on current stiffness matrices. These calculation steps are managed by the keyword factor CRIT_STAB with the option TYPE = “FLAMBEMENT” (cf. [U4.51.03] and [U4.53.01]).
Notes for transient analysis
Option CRIT_STAB acts as for the quasistatic case: we always conduct a buckling analysis (therefore based only on the study of stiffness matrices), not a dynamic instability analysis (damping becoming negative).
With coupled fluid-structure modeling \((u,p,\phi )\) [R4.02.02], you have to modify the assembled stiffness matrix (as well as the geometric stiffness when it is used). To do this, you must enter the following keywords, sub CRIT_STAB :
MODI_RIGI = “OUI”,
DDL_EXCLUS =( “PHI”, “PRES”, “DH”,).
The list of excluded degrees of freedom must include all types of degrees of freedom related to the fluid model: in the example of test case FDNV100, we therefore have the potential PHI , the pressure PRES and the vertical displacement at the level of the free surface DH. If we don’t do this treatment, then the call to CRIT_STAB will crash because of a singular matrix and no lag strategy can overcome this.
In quasistatics, this problem does not arise because coupled fluid-structure modeling then makes no sense.
This keyword makes it possible to trigger the calculation, at the end of each time increment, of a stability criterion. This criterion is useful for detecting, during loading, the point at which stability is lost (for example by buckling).
This criterion is calculated as follows: at the end of a time step, in small disturbances, we solve \(\mathit{det}\left({K}^{T}+\mathrm{\mu }\mathrm{.}{K}^{g}\right)=0\). \({K}^{T}\) is the coherent tangent matrix at this point in time. \({K}^{g}\) is the geometric rigidity matrix, calculated from the stress field at that moment and \(\mathrm{\mu }\) is the critical load coefficient.
In practice, the load is unstable if:math: |mathrm {mu} |<1 (actually:math: 0<mathrm {mu} <1). The eigenvalues are calculated using the Sorensen method (cf. CALC_MODES [U4.52.02]). This can be quite expensive for larger problems.
The CHAR_CRIT keyword saves time by only doing a Sturm test in the frequency band provided. If we find at least one frequency, then we actually calculate the values of the critical loads in this range.
For large displacements and large deformations, we solve \(\mathit{det}\left({K}^{T}+\mathrm{\mu }\mathrm{.}I\right)=0\) because \({K}^{T}\) then contains \({K}^{g}\).
The criterion is then an instability criterion: when \(\mathrm{\mu }\) changes sign (therefore passes through 0) the load is unstable.
Notes on the geometric stiffness matrix
The stability analysis will be more accurate if we have the geometric stiffness matrix: in fact, in its absence, certain instabilities cannot be detected as, for example, in purely elastic cases. Between sufficiently realistic models, it is therefore necessary to favor the use, for precise stability analyses, of models allowing the calculation of this geometric stiffness.
The keyword NB_MODE refers to the number of critical loads to be calculated. Often the first is enough but there can be multiple modes.
The eigenmode corresponding to the lowest critical load (in absolute value) is stored in the result object, under the name MODE_FLAMB. This clean mode can be extracted and visualized (like a travel field or a classic clean mode). It is normalized to 1 on the largest displacement component. All calculated critical loads are shown in the.mess file.
In practice, in order to limit the additional calculation cost, it is recommended to optimize the calls to CRIT_STAB. To do this, you can use the keywords LIST_INST/INST/PAS_CALCUL in the keyword factor CRIT_STAB. It is thus possible to specify at which time steps the buckling modes will be calculated.
In addition, it is a good idea to use CRIT_STAB only on time intervals where there is a suspicion of the possibility of instabilities.
Finally, if we want a very good assessment of critical loads, it is necessary to refine the time step when approaching this zone. This advice is all the more relevant in semi-static cases because the user then often resorts to longer time steps than in dynamics.
It is possible to properly stop (reusable base in progress) a calculation with STAT_NON_LINE or DYNA_NON_LINE when an instability is detected. It is not the step by default where the code will try to continue to solve the problem, and if we have convergence, then it means that we have succeeded in following one of the branches of solution.
To manage this stop (cf. test case ERREU10), you must have previously used DEFI_LIST_INST with the following arguments:
ECHEC =_F (EVENEMENT =” INSTABILITE “, ACTION =” ARRET”,),
This means that in the event of an instability-type event, the action triggered will be the stop.
Under CRIT_STAB, the following optional options allow you to control this stopping criterion:
PREC_INSTAB to define the (dimensionless) precision of the stopping criterion,
SIGNE to specify the critical values to consider.
The second keyword is only used when the geometric stiffness matrix is used. By default, the solution is considered to be unstable if the critical load becomes between 1 and -1, but if necessary, only the positive or negative part of this interval can be taken into account.
Without a geometric stiffness matrix, instability will be detected when an eigenvalue of the overall stiffness matrix assembled, i.e.:
will tend to 0 (with a relative precision given by PREC_INSTAB),
will change sign.
For example, in the case of a reservoir filled with water under an earthquake, we can start the incremental or transitory calculation with a critical load equal to 0.8 (analysis with geometric stiffness): which means that the reservoir would explode if we imposed a depression equal to 0.8 times the imposed hydrostatic pressure (the positive value of the critical load corresponds to a reversal of the direction of the load in question). So if we don’t specify anything, the calculation would be considered unstable and would stop. Since, in this case, we assume that there will be no depressurization (for example by emptying), then we clear the \(\mathrm{[}\mathrm{0,1}\mathrm{]}\) interval in the stability analysis. So the problem will become unstable if the critical load reaches interval \(\mathrm{[}\mathrm{-}\mathrm{1,0}\mathrm{]}\).
For analyses in monotonic evolution, this type of reasoning is easily conceivable, which is the case for the static part of the load, but for the dynamic part, the load being cyclical, and unless specific information is available, it is safer and more conservative to keep the option by default and therefore to consider the structure to be unstable if the critical load becomes less than 1 in absolute value.
During a stop due to instability, the calculation will stop by properly closing the database: the user will be able to use it continuously.
2.2. Unstable solution tracking#
For the study of a structure that is potentially unstable or likely to experience an end point, which is therefore likely to encounter a bifurcation in solution during the evolution of the load, it is often useful to be able to choose a particular branch of solution (often the physical solution when it is defined a priori without ambiguity). To do this, the user may have to introduce an initial defect that will « force » the structure to branch off to the particular solution branch.
Several methods exist to define this defect.
One of the most suitable is to pre-deform the structure slightly according to the shape of the buckling mode corresponding to the branch you want to follow. The amplitude of this pre-deformation must be low, for example less than 1/10th of the thickness for a thin structure. The ideal is to find the minimum defect that is compatible with a satisfactory performance of the residue algorithm in balance. In fact, a defect that is too weak can cause a difficulty in convergence of the residue, mainly in the case of force control. This initial defect can very wisely be constructed from the instability mode calculated with option CRIT_STAB of STAT_NON_LINE/DYNA_NON_LINE (*cf. previous paragraph). This mode then takes into account all the nonlinearities introduced in the complete model. The more economical alternative is to use the Euler buckling mode, but which corresponds to the linear case.
The geometric defect can also be defined by experimental measurements of the real part whose geometry cannot be perfect.
The defect can also take the form of a disturbance in the loading (misalignment, addition of a localized load,…) or in the mechanical characteristics of the material (local weakening of the Young’s modulus, for example). However, it may then be more difficult to adapt the defect to the desired buckling mode, especially if the structure has relatively similar modes.
Note
In some cases, even on the undisturbed problem, the load is such that it causes the desired bifurcation.
One of the other particular points, linked to instability, is the choice of the technique for controlling algorithm STAT_NON_LINE. In fact, conventional effort control is no longer suitable because it cannot capture an unstable branch of solution. Likewise, when approaching a limit point, convergence with force control will become more and more difficult, the tangent stiffness matrix becoming singular. It is then necessary to reduce the load increment and to increase the maximum number of iterations in order to continue the calculation.
It is also possible to use the possibility of stopping properly in the event of instability (cf. previous paragraph) to continuously manage the bifurcation on the chosen branch of solution by initiating the rest of the calculation by a disturbance according to this mode of instability.
There are control techniques [bib9] to circumvent these numerical difficulties. Among the methods proposed by Code_Aster, the one called by arc length [bib12] (option TYPE =” LONG_ARC “from the key word PILOTAGE in STAT_NON_LINE), which is the most suitable for buckling instabilities, in the case of possible « soft » snap-backs [bib13]. For more brutal snap-backs, Crisfield offers a variant [bib13], not available in Code_Aster.
Other methods exist, such as Riks’s [bib14] (not available either), which also treats the dynamic case.
If you only want to obtain the limit point, even with good precision, load control may suffice, provided you manage the parameters of no load increment (remember to use the DEFI_LIST_INST command) and the maximum authorized number of iterations (ITER_GLOB_MAXI of CONVERGENCE). It may also be useful, when approaching the limit point, to no longer use the updated tangent matrix for the solver, since it is almost singular. We can then content ourselves with not updating this matrix at each calculation (parameters REAC_INCR and REAC_ITER) or, in the worst case, adopt the basic elastic matrix (PREDICTION =” ELASTIQUE “and MATRICE =” ELASTIQUE” and =” “of the keyword NEWTON).
Here is an example of the use of STAT_NON_LINE for an elastoplastic calculation in large displacements ([bib4] for the elements used, which are of the solid shell type), with force control:
RESU = STAT_NON_LINE (…
EXCIT = (_F (CHARGE = CONDLIM,
TYPE_CHARGE = “FIXE_CSTE”,),
_F (CHARGE = PESA,
TYPE_CHARGE = “FIXE_CSTE”,),
_F (CHARGE = PRESPH,
FONC_MULT = FONCMUL2,
TYPE_CHARGE = “SUIV”,),
_F (CHARGE = PRESPS1,
FONC_MULT = FONCMUL,
TYPE_CHARGE = “SUIV”,),),
…)
Remarks
The updated tangent matrix is used for each calculation, allowing the sub-division of the load step.
The pressures imposed are follow-up efforts (TYPE_CHARGE =” SUIV “) .
In the case of modeling with massive elements, the recommended deformation tensor in large displacements is” SIMO_MIEHE “ .
If you want to replace force control by an arc-length method, all you have to do is add:
RESU = STAT_NON_LINE (…
PILOTAGE = _F (GROUP_NO = “G”,
TYPE = “LONG_ARC”,
NOM_CMP = (“DY”,),
COEF_MULT = 7.),
…)
Remarks
In Code_Aster, you cannot control follower forces.
For arc-length control, it is, in general, recommended that GROUP_NOcontienne the entire structure.
Finally, let’s cite two articles by Crisfield that give a good general overview of the problems and methods related to nonlinear calculations that can present various types of instabilities ([bib15] and [bib11]).
The [U2.06.11] documentation shows an example of using CRIT_STAB to study the behavior of a metal tank.
Some Code_Aster test cases dealing with buckling:
Euler modes:
sdls504
sdls505
ssll103
ssll105
ssll403
ssll404
ssls110
Euler modes and nonlinear calculation:
ssnl123
Nonlinear modes (CRIT_STAB):
sdnv106 (also present MODE_VIBR for the calculation of vibration modes on updated stiffness)
ssll105
ssnl126
SSNP306
Nonlinear calculation:
ssnl502
ssnp305: calculation up to a snap-through