1. Presentation#
1.1. Generalities#
STAT_NON_LINE is the Code_Aster operator for performing non-linear mechanical calculations when inertia effects are neglected (if inertia effects are to be taken into account, use DYNA_NON_LINE, see [R5.05.05]).
In principle, the calculation only concerns mechanical variables (displacements, constraints, internal variables), excluding any coupling with other physical phenomena (thermal,…). Consequently, the associated fields influencing mechanical behavior (thermal, hydraulic, metallurgical fields) are calculated beforehand by other operators (THER_LINEAIRE [U4.33.01], THER_NON_LINE [U4.33.02]), or even by other codes (for example CODE_SATURNE for fluid mechanics,…).
There is an exception with regard to thermo-hydro-mechanical modeling (modeling called “THM”) for which STAT_NON_LINE treats the entire problem coupled with the thermal diffusion equations, the pressure of the fluid (s) and the mechanical balance [R7.01.10].
It should be noted that when we talk about the moment of calculation in this document, we almost always refer to a pseudotime, which has no physical meaning and which only serves to configure the incremental algorithm. However, the moment retains physical significance in visco-plasticity and when the control variables depend on it.
1.2. Types of nonlinearities#
1.2.1. Nonlinear behaviors#
Nonlinear behavioral relationships are described in documents [R5.03….], for general behaviors, and [R7.01….] for geo-materials. The choice of the type of behavior is made by the keyword factor COMPORTEMENT. Through the equilibrium equation, this leads to a non-linear system that can be of two forms:
The system depends explicitly on the displacement field \(\mathrm{u}\) in relation to the reference configuration, and parameterized by the moment of calculation (through, among other things, thermal evolution).
The system depends implicitly on the displacement field :math:`mathrm{u}` by an**implicit differential equation* (for example elasto-plasticity, visco‑plasticity, hypo-elasticity, etc.). In this case, the behavior relationship is integrated as presented for example in [R5.03.02]: by relating a displacement increment \(\Delta \mathrm{u}\) calculated from a given mechanical state (the mechanical state being represented by a displacement field \(\mathrm{u}\), a stress field \(\sigma\) and a field of internal variables and a field of internal variables \(\alpha\)) to the stress field at the time \(t\) of calculation. The equilibrium equation therefore leads to a non-linear system in \(\Delta \mathrm{u}\), but which is also parameterized by the moment of calculation through the data of the problem (mechanical load variation and thermal evolution for example).
1.2.2. Major transformations#
In cases where the hypothesis of small disturbances (moderate displacements and deformations) is not verified, the method of solving the problem must then integrate the evolution of the geometry of the problem, manipulate a particular kinematics and use an adequate formulation of the law of behavior.
In practice, the hypothesis of small deformations can be applied as long as the square of the deformation modulus remains less than the precision of the calculations envisaged. Likewise, the hypothesis of small rotations can be applied as long as the product between the square of the angle of rotation and the deformation modulus remains less than the precision of the calculations envisaged.
Various alternatives exist within Code_Aster; our aim is not to present them in detail here and we refer to the various documents dealing specifically with each problem. A distinction should be made between formulations that operate on massive isoparametric elements (2D or 3D) and formulations used for structural elements (beams, plates and shells). For the cases of massive isoparametric elements, there are three main types of formulation of kinematics for the case of large deformations:
The kinematics DEFORMATION =” PETIT_REAC “makes it possible to treat any law of behavior in large deformations. The law is written in small deformations and large deformations are taken into account only by updating the geometry. This formulation is not incrementally objective. Moreover, it can only be used if the behavior is isotropic, if the elastic deformations are low compared to the plastic deformations, if the rotations remain low (less than \(10°\)) and if a sufficiently fine discretization is adopted in time. In addition, the absence of the geometric contribution in the tangent matrix can sometimes make convergence difficult (see [R5.03.21 for more details).
The kinematics DEFORMATION =” SIMO_MIEHE “makes it possible to treat a law of elastoplastic behavior with isotropic work hardening in large deformations, the law of ductile rupture called « Rousselier » or elasto (visco) plasticity with phase change for metallurgy (see [R5.03.21], [R5.03.06] and [R4.04.03]). This formulation is incrementally objective, with no limitation on the level of transformations applied, but it is only available for the three laws of behavior mentioned and only makes it possible to treat cases where the behavior is isotropic.
The kinematics DEFORMATION =” GDEF_LOG “also makes it possible to treat any law of hypo-elastoplastic behavior. It is incrementally objective, without limitation on the level of transformations applied and makes it possible to deal with the case of anisotropic behaviors (see [R5.03.24]).
To deal with large elastic deformations, another formalism must be used, called via DEFORMATION = “GROT_GDEP”, which can be used for non-linear hyper-elastic behavior relationships in large displacements (see [R5.03.20] and [R5.03.22] for the case of small deformations) or for the law of hyperelastic behavior (see [R5.03.19]).
Finally, for structural elements (beams, plates or shells), there are specific formulations. We can cite:
Beams in large movements (see [R5.03.40]) or multifibre beams in large movements (see [R3.08.09]). Keyword DEFORMATION =” GROT_GDEP “.
Solid shell elements in geometric nonlinear (see [R3.07.05]).
There are no structural elements (beam, plate or shell) that can be used in large deformations in Code_Aster.
1.2.3. Unilateral contact and friction#
For contact and friction, reference will be made to three documents: [R5.03.51] on discrete contact with friction, [R5.03.52] for the hybrid formulation using contact/friction elements and [R5.03.53] on contact in large sliding movements with method XFEM.
1.2.4. Limit conditions and loads#
It is possible to define non-linear boundary conditions, i.e. dependent on the movement of the structure. There are two main categories of non-linear boundary conditions:
For Neumann loads, it is essentially pressure: this pressure must remain normal at the surface, the fact that the structure undergoes large transformations implies that the load*depends on displacements. In this regard, we will consult the documents [R3.03.07], [R3.03.04] and § 3.2;
For Dirichlet loads, Code_Aster currently only allows loading of type LIAISON_SOLIDE (see [U4.44.01]) which exactly rigidifies part of the structure. We will consult [R3.02.02] on this subject;
In both cases, it is necessary to specify explicitly that this load is non-linear (we talk about follower loading) using the keyword TYPE_CHARGE in the EXCIT keyword factor defining the loads (see [U4.01.03]).
Note: in general, contact and friction can be assimilated to a mixed limit condition (in movement and in pressure) which is non-linear.
1.3. Position of the nonlinear quasistatic problem#
1.3.1. General problem#
As a result of paragraph 1.1, we can see that it is legitimate to consider that the non-linear problem has a displacement as unknown and that it is parameterized by time. Let us therefore consider the nonlinear, quasi-static problem using the classical expression of the principle of virtual work.
where:
\(t\) represents the instant variable;
\(u\) is the displacement field taken from a reference configuration;
\(v\) is the kinematically admissible virtual displacement field;
\({L}^{\text{ext}}\) is the external mechanical load to which the structure is subjected (pressure, imposed force, … );
\({L}^{\text{int}}\) represents the internal forces of the nonlinear quasistatic mechanics problem. In the linear case, we have \({L}^{\text{int}}(u,t)=K\text{.}u\), where \(K\) is the stiffness matrix of the structure.
For the time being, the Dirichlet boundary conditions are not taken into account. In fact, more precisely, \({L}^{\text{int}}(u,t)\) is linked to the stress field \(\sigma\) by the virtual deformation work operator \({Q}^{T}\) according to the following relationship:
For small trips, \({Q}^{T}\) is independent of trips; for large trips, this is no longer the case. We give ourselves a discretization of the time interval to be calculated:
The constraint field \({\sigma }_{i}\) at time \({t}_{i}\) is written \(\sigma ({u}_{i},{\beta }_{i},{t}_{i},{H}_{i-1})\), if we write the command variable fields \({\beta }_{i}\) and \({H}_{i-1}\) the past history of the structure. For elastic behaviors, history does not intervene: set \({H}_{i-1}\) is therefore empty. For incremental behaviors, the story is the set of states (travel fields, constraints and internal variables) at the previous moment: \({H}_{i-1}=\left\{{u}_{i-1},{\sigma }_{i-1},{\alpha }_{i-1},{t}_{i-1}\right\}\).
In the general case, the dependence of the operator \({L}^{\text{int}}\) is, as we saw in [§ 1.1], implicit in relation to time: it results from the integration of the behavioral relationship over time (for elasto-plasticity problems for example). Explicit dependence on time appears in particular in the case of behavioral relationships taking into account a phenomenon of hardening by time (time-hardening) or in the case of aging.
The formulation of the quasi-static problem consists in expressing the balance of the structure (the internal forces are equal to the external forces) for a series of calculation moments \({\left\{{t}_{i}\right\}}_{1\le i\le I}\) which parameterize the loading, we will note the quantities at time \({t}_{i}\) by the index \(i\) (for example \({\mathrm{L}}^{\text{int}}({\mathrm{u}}_{\mathrm{i}},{t}_{i})={\mathrm{L}}_{i}^{\text{int}}\)):
: label: eq-4
{mathrm {L}} _ {i} ^ {text {int}} = {mathrm {L}} _ {i} ^ {text {ext}} ^ {text {ext}}
In \(({\mathrm{u}}_{i},{t}_{i})\), this amounts to cancelling the \({\mathrm{R}}_{i}({\mathrm{u}}_{i},{t}_{i})\) vector, known as the residual equilibrium vector, defined by:
: label: eq-5
{mathrm {R}} _ {i} ({mathrm {u}}} _ {i}, {t} _ {i}) = {mathrm {L}} _ {i} _ {i} ^ {text {int}} ^ {text {int}} ^ {text {int}} ^ {text {int}}} - {mathrm {L}} _ {i}} _ {text {ext}}
The state of the structure in \({t}_{0}\) is assumed to be known. We perform \(I\) load increments (or no steps). The unknowns are calculated incrementally by the global resolution algorithm (even for elastic behaviors). Starting with \({\mathrm{u}}_{i-1}\), a solution satisfying the balance in \({t}_{i-1}\), we determine \(\Delta {\mathrm{u}}_{i}\) which will make it possible to obtain the solution in \({t}_{i}\):
Increment \(\Delta {u}_{i}\) is first estimated by linearizing the problem with respect to time around \(({\mathrm{u}}_{i-1},{t}_{i-1})\) (so-called prediction or Euler phase). Then we use a Newton method or one of its variants to solve equation () iteratively (we calculate a sequence \(\Delta {\mathrm{u}}_{i}^{n}\) where the exponent \(n\) is the iteration number). In addition to these variables, for incremental behavioral relationships, we need to know in \({t}_{i-1}\) the constraint field \({\sigma }_{i\mathrm{-}1}\) and the internal variables field \({\alpha }_{i-1}\) (see [R5.03.02] for an example).
1.3.2. Dirichlet boundary conditions#
We have imposed that the virtual travel field be cinematically eligible, that is to say that it meets the travel limit conditions (or Dirichlet conditions). In general, linear Dirichlet conditions are written as:
With \({u}_{j}\) the list of \(r\) degrees of freedom, \({\alpha }_{j}\) the coefficients and \(\beta\) the second member that can depend on time.
In Code_Aster, there are two ways to take these conditions into account:
By a method called « kinematic » (definition by AFFE_CHAR_CINE) which uses a method of pseudo-elimination in the resulting linear system (see [U4.44.03]). In this case, one can only define simple boundary conditions such as \(r=1\);
By a dualization method (definition by AFFE_CHAR_MECA) that uses the Lagrange double technique (see [R3.01.01]).
It is possible to eliminate dualized boundary conditions by using a special method described in [R3.03.05], but only by using the PETSc solver.
We now describe how to take dualised boundary conditions into account when solving a non-linear problem. First, we will assume that the boundary conditions () can be written in the following discrete form:
\(\mathrm{B}\) is a linear operator for the space of displacement fields on a space of functions defined on a part of the edge of the structure, \({u}^{d}\) is a function given on this part. The dualization of the conditions at the Dirichlet \(B\mathrm{.}u={u}^{d}(t)\) limits leads to modification () in the following way:
The unknowns are now, at all times \(t\) * , the \((u,\lambda )\) couple, where \(\lambda\) represents the Lagrange multipliers of the Dirichlet boundary conditions. The vector \({B}^{T}\text{.}\lambda\) is then interpreted as the opposite of the support reactions at the corresponding nodes.
By noting the quantities at time \({t}_{i}\) by the index \(i\):
The equilibrium residue vector \({R}_{i}({u}_{i},{\lambda }_{i},{t}_{i})\) is now equal to:
The unknowns are always calculated incrementally by the global resolution algorithm. Starting with \(({u}_{i-1},{\lambda }_{i-1})\), a solution satisfying the balance in \({t}_{i-1}\), we determine \(\Delta {u}_{i}\) and \(\Delta {\lambda }_{i}\) which will make it possible to obtain the solution in \({t}_{i}\):