2. Problem position#

We place ourselves in the general framework of solving a typical problem of non-linear mechanics discretized in space \((K\in \Omega )\) and time \((t\in [\mathrm{0,}T])\), and written in motion. Its resolution, at moment \({t}_{n\text{+}1}\), consists in determining the displacements \({U}_{n\text{+}1}\) (therefore the deformations \(\varepsilon ({U}_{n\text{+}1})\)), the internal variables \({\alpha }_{n\text{+}1}\) and the constraints \({\sigma }_{n\text{+}1}\) verifying:

  1. At the global level, the balance of forces:

\({F}_{\text{ext}}({t}_{n\text{+}1})\mathrm{-}{F}_{\text{int}}({\sigma }_{n\text{+}1}({U}_{n\text{+}1}),{t}_{n\text{+}1})\) (1)

  1. At the local level, the equations constituting the law of behavior considered for the material:

State law: \({\sigma }_{n\text{+}1}\mathrm{=}\mathrm{\sum }(\varepsilon ({U}_{n\text{+}1}),{\alpha }_{n\text{+}1})\)

Law of evolution: \(\mathrm{\{}\begin{array}{cc}f({\alpha }_{n\text{+}1},{\sigma }_{n\text{+}1})\le 0& \text{Convexe de réversibilité}\\ {\dot{\alpha }}_{n\text{+}1}\mathrm{=}{\dot{\lambda }}_{n\text{+}1}\ge 0& \text{Evolution des variables internes}\\ {\dot{\lambda }}_{n\text{+}1}f({\alpha }_{n\text{+}1},{\sigma }_{n\text{+}1})\mathrm{=}0& \text{Condition de Kuhn-Tucker}\end{array}\) (2)

with \({F}_{\text{ext}}\) and \({F}_{\text{int}}\) the external and internal forces, and \(\dot{{\lambda }_{n\text{+}1}}\) the multiplier (of damage or plasticity). In the following, a behavior independent of physical time is considered (no viscosity or dynamic effect); the « pseudo-instant » \({t}_{n\text{+}1}\) therefore represents the load factor applied.

In the vast majority of cases, each of the equations (1) and (2) is relatively easy to solve.

At the global level, Newton’s method leads to a succession of linear problems of the type:

\({K}_{n\text{+}1}^{i}\delta {U}_{n\text{+}1}^{i\text{+}1}={F}_{\text{ext}}({t}_{n\text{+}1})-{F}_{\text{int}}({\sigma }_{n\text{+}1}^{i}({U}_{n\text{+}1}^{i}))\) (3)

where the high index represents the Newton iteration under consideration, \(\delta {U}_{n\text{+}1}^{i\text{+}1}\) is the displacement increment between two successive Newton iterations and \({K}_{n\text{+}1}^{i}\) is the global tangent matrix. Note that this global tangent matrix can be written (cf. [2])

\(\begin{array}{ccc}{K}_{n\text{+}1}& \text{=}& \frac{\mathrm{\partial }{F}_{\text{int}}({\sigma }_{n\text{+}1}({U}_{n\text{+}1}))}{\mathrm{\partial }{U}_{n\text{+}1}}\\ & \text{=}& {\text{A}}_{e\text{=}1}^{p}({\mathrm{\int }}_{{\Omega }_{\varepsilon }}\mathrm{\nabla }{N}_{e}\mathrm{.}{C}_{{e}_{n\text{+}1}}\mathrm{.}\mathrm{\nabla }{N}_{e}d\Omega )\end{array}\) (4)

where \(\text{A}\) is the assembly operator, \(p\) the number of elements in the mesh, \({C}_{{e}_{n\text{+}1}}\) the local tangent operator (resulting from the behavior) and \({N}_{e}\) the shape functions (for the element in question, noted \(e\)). Quantity \({F}_{\text{ext}}({t}_{n\text{+}1})-{F}_{\text{int}}({\sigma }_{n\text{+}1}^{i}({U}_{n\text{+}1}^{i}))\) represents the equilibrium residue, which must be cancelled from a favourable point of view. Determining \(\delta {U}_{n\text{+}1}^{i\text{+}1}\) is therefore limited to the reversal of the \({K}_{n\text{+}1}^{i}\) matrix. If therefore the local elements resulting from the behavior (constraints, tangent matrix) are known and fixed, and generate a well-conditioned \(K\) matrix, the resolution is easy.

Moreover locally, if the displacements \({U}_{n\text{+}1}\), and therefore the local deformations \(\varepsilon ({U}_{n\text{+}1})\), are known, the resolution of the majority of the laws of behavior, that is to say the determination of the constraints \({\sigma }_{n\text{+}1}\), of the internal variables \({\alpha }_{n\text{+}1}\) and of the local tangent matrix \({C}_{{e}_{n\text{+}1}}\), is relatively easy.

It is therefore in the simultaneous nature of these two laws that the difficulty of implicit resolution lies. The most classical method of resolution is Newton’s iterative algorithm (cf. [2]). This iterative pattern does not always converge. Although this remains subject to discussion (there are other possible factors for loss of convergence, such as leaving Newton’s basin of attraction or the presence of bifurcations in the solution), the authors mainly attribute this loss of robustness to the singular nature of the matrices, in particular for the laws of damage, in particular for the laws of damage, in particular for the laws of damage, whose internal variables impact the matrices of rigidity of behavior: in cases of significant damage, close to or reaching rupture, softening leads to locally very tangent matrices » weak », which may in the extreme no longer be defined as positive. Through the assembly process (4), as damage progresses in the structure, the overall matrix \(K\) is deteriorated: it becomes too « soft », and its minimum eigenvalue tends to zero. It then becomes singular and the algorithm diverges. The robustness of the calculation is then no longer ensured.

To increase the robustness of the calculation in these situations, Oliver*et al.* [1] propose a particular method of resolution, called IMPLEX. The global balance is then verified approximately using a local tangent matrix (secant in the case of damage laws) extrapolated explicitly and assembled, and the laws of behavior are resolved implicitly from the approximate displacement field. The method is shown below.