1. Theoretical formulation of limit analysis#

1.1. Definition of limit load#

We consider a solid occupying a bounded domain \(\Omega\) subjected to surface loads \(\lambda F+{F}_{0}\) on the edge \({\Gamma }_{f}\) and volume loads \(\lambda f+{f}_{0}\) on \(\Omega\). A distinction is made between loading \((F,f)\), set by the positive real number \(\lambda\), and loading « permanent » \(({F}_{0},{f}_{0})\). The homogeneous Dirichlet conditions that are applied to the complementary edge \({\Gamma }_{u}\) of \(\partial \Omega\) (an imposed displacement or an initial anelastic deformation —thermal, plastic… — have no effect on the range of permissible loads). Several other useful properties can be found in [bib5].

The material constituting the solid has a strength criterion expressed by a scalar stress function, which is negative for the admissible stresses. The criterion used for a material of the perfect elastoplastic type with a VonMises threshold and retained here is:

\(g(\sigma )=J(\sigma )-{\sigma }_{y}=\sqrt{}\text{.}\sqrt{{\sigma }^{D}\text{.}{\sigma }^{D}}-{\sigma }_{y}=\frac{\sqrt{2}}{2}\text{.}\sqrt{{({\sigma }_{1}-{\sigma }_{2})}^{2}+{({\sigma }_{2}-{\sigma }_{3})}^{2}+{({\sigma }_{1}-{\sigma }_{3})}^{2}}-{\sigma }_{y}\)

\({\sigma }^{D}\)

is the stress tensor deviator,

\({\sigma }_{y}\)

is the simple tensile strength threshold (such as an elastic limit), which may vary depending on the areas of the solid in question.

\({\sigma }_{i}\)

being the main constraints of \(\sigma\).

Given this resistance criterion, we seek to calculate the limit value of \(\lambda\), called the limit load \({\lambda }_{\text{lim}}\), for which the structure can support the loads \({\lambda }_{\text{lim}}F+{F}_{0}\) and \({\lambda }_{\text{lim}}f+{f}_{0}\).

Strictly speaking, the value \({\lambda }_{\text{lim}}\) refers to the limit of bearable loads, but for materials that obey the Principle of Maximum Plastic Work, this value is the limit of loads supported.

1.2. Calculation of the limit load using a kinematic approach#

In fracture calculation, two approaches are possible: the static approach (in stress variables) and the kinematic approach (in speed variables). These approaches provide limits of the load limit: a minor for the static approach and an increase for the kinematic approach. When both provide the same result, the limit load obtained is accurate.

The kinematic approach is the one used in*Code_Aster* using finite elements in displacement. For the given load

_images/Object_21.svg

, the space of kinematically admissible and normalized velocities is defined by:

\({V}_{a}^{1}=\left\{v\text{admissible,}v=0\text{sur}{\Gamma }_{u}\text{,}L(v)={\int }_{\Omega }f\text{.}vd\Omega +{\int }_{{\Gamma }_{f}}F\text{.}v\text{ds}=1\right\}\)

This standardization forces the work of loading.

_images/Object_23.svg

to be unitary. The power of « permanent » charging \(({F}_{0},{f}_{0})\) is noted: \({L}_{0}(v)\).

Based on the stress resistance criterion \(g(\sigma )\), we define:

  • the set of constraints allowed by: \({G}_{(x)}=\left\{\sigma (x)\text{,}g(\sigma (x))\le 0\text{}\right\}\)

(\({G}_{(x)}\) is convex for criterion \(g\))

  • the indicator function: \({\Psi }_{G}(\sigma (x))=\{\begin{array}{c}0,\text{si}\sigma (x)\in {G}_{(x)}\\ +\infty ,\text{si}\sigma (x)\notin {G}_{(x)}\end{array}\)

  • the support function: \(\pi (\varepsilon )=\underset{\sigma \in {\mathrm{IR}}^{6}}{\text{Sup}}\left[\sigma \text{.}\varepsilon -{\Psi }_{G}(\sigma )\right]\)

The Sup in \(\pi (\varepsilon )\) can only be achieved if \(\sigma\) is chosen in \({G}_{(x)}\), such as: \(\sigma =\lambda {\varepsilon }^{D}+\mu \text{Id}\) (which ensures \(\sigma \text{//}{\varepsilon }^{D}\)). The optimum corresponds to \(g(\stackrel{ˉ}{\sigma })=0\Rightarrow \stackrel{ˉ}{\lambda }={\sigma }_{y}\sqrt{}\text{.}{({\varepsilon }^{D}\text{.}{\varepsilon }^{D})}^{\text{-}1/2}\)

_images/Object_38.svg

Figure 1.2-a: Optimum \(\stackrel{ˉ}{\sigma }\) and graph of the \(\pi (\varepsilon )\) n function in 1D

Hence the support function:

\(\pi (\varepsilon (v))={\sigma }_{y}\text{.}\sqrt{\frac{2}{3}}\sqrt{\epsilon (v)\text{.}\epsilon (v)}+\underset{\mu \in \text{IR}}{\text{Sup}}(\mu \text{.}\text{div}v)={\pi }_{R}(\varepsilon (v))+\underset{\mu \in \text{IR}}{\text{Sup}}(\mu \text{.}\text{div}v)\)

We observe that the function \(\pi (\varepsilon )\), which is interpreted as the power density that can be dissipated at the hardware point, is not differentiable in \(0\).

To date, Code_Aster treats the possible internal surfaces of discontinuity within solid \(\Omega\) [bib 4].

The kinematic approach is defined using the convex functional \({S}_{e}(v)\), positively homogeneous of degree one, for \(v\in {V}_{a}^{1}\) defined over the entire domain:

\(\begin{array}{cc}{S}_{e}(v)& ={\int }_{\Omega }\pi (\varepsilon (v))d\Omega \end{array}-{L}_{0}(v)\)

This functional is the integral over the domain of the support function \(\pi\) of the convex \({G}_{(x)}\), calculated in \(\varepsilon (v)\) and is interpreted as the maximum resistive power in the speed field

_images/Object_51.svg

(the interface resistance contribution on discontinuous surfaces is assumed to be zero). The support function \(\pi\) is positively homogeneous by degree 1, and therefore the functional \({S}_{e}(v)\) also as a result.

With the Von Mises criterion the power functional \({S}_{e}(v)\) is:

\({S}_{e}(v)={\int }_{\Omega }\left[{\sigma }_{y}\text{.}\sqrt{\frac{2}{3}}\sqrt{\varepsilon (v)\text{.}\varepsilon (v)}+\underset{q\in \text{IR}}{\text{Sup}}(q\text{.}\text{div}v)\right]d\Omega -{\text{L}}_{0}(v)\) eq 1.2-1

where we see that only fields \(v\) belonging to \(C=\left\{v\in {V}_{a}^{1},\text{div}v=0\text{dans}\Omega \right\}\) provide finite values. Fields \(v\) must therefore verify the so-called incompressibility condition \(\text{div}v=\text{tr}\varepsilon (v)=0\). This is why it is necessary to use quasiincompressible elements for a limit load calculation with the Von Mises criterion [R3.06.08].

The load limit \({\lambda }_{\text{lim}}\) given by the kinematic approach is:

\({\mathrm{\lambda }}_{\text{lim}}=\underset{\mathrm{v}\in {V}_{a}^{1}}{\text{Inf}}{S}_{e}\left(\mathrm{v}\right)=\underset{\begin{array}{c}\mathrm{v}\in {V}_{a}^{1}\\ L\left(v\right)>0\end{array}}{\text{Inf}}\frac{{S}_{e}\left(v\right)}{L\left(v\right)}=\underset{l>0}{\text{Sup}}\underset{\mathrm{v}\in {V}_{a}}{\text{Inf}}\left({S}_{e}\left(\mathrm{v}\right)-l\left(L\left(\mathrm{v}\right)-1\right)\right)\)

where \(l\) is the Lagrange multiplier.

At the optimum, a solution \(u\) and the limit charge \({\lambda }_{\text{lim}}\) are obtained (not the uniqueness of \(u\) but the uniqueness of \({\lambda }_{\text{lim}}\)). So any \({L}_{0}(v)+\lambda L(v)\) load with \(0\le \lambda \le {\lambda }_{\text{lim}}\) is bearable. Beyond \({\lambda }_{\text{lim}}\), the balance problem has no solution.

Note:

There are situations where, even if \({L}_{0}(v)\) is not bearable alone, the combination \({L}_{0}(v)+\lambda L(v)\) , for \({\lambda }_{1}\le \lambda \le {\lambda }_{2}\) , becomes bearable over a certain interval, and not only for two collinear loads.

Note:

The limit load calculated for a two-dimensional problem, in plane deformations, is necessarily greater than that obtained for this problem modelled in plane stresses. This result therefore provides a benefit. If one wishes to deal with a problem in plane constraints, it is then necessary to do the kinematic approach on a three-dimensional model.

1.3. Regularization of the kinematic approach by the Norton-Hoff-Friaâ method#

The numerical implementation of the kinematic approach requires the minimization of the non-differentiable functional \({S}_{e}(v)\). Numerous regularization techniques exist [bib4]. The Norton-Hoff-Friâa method is used here [bib2], [bib7]. It is based on groundbreaking work by Casciaro in 1971. It consists in replacing the support function \(\pi (\varepsilon )\) by the regularized and differentiable support function \({\pi }^{\mathrm{NH}}(\varepsilon )\). It is adjustable by a regularization parameter \(m\) (\(1\le m\le 2\)), whose limit value \(m\to {1}^{+}\) leads to convergence towards the support function \(\pi (\varepsilon )\):

\({\pi }^{\text{NH}}(\varepsilon )=\frac{{k}^{1\text{-}m}}{m}{(\pi (\varepsilon ))}^{m}\) eq 1.3-1

The scalar \(k\) in [éq 1.3-1] is homogeneous to a constraint. We note the space of the admissible velocities adapted to the viscous flow problem for the Norton-Hoff law of order \(m\):

\({V}_{a}^{\mathrm{m1}}=\left\{v\in {L}^{m}(\Omega ),\text{et}\varepsilon (v)\in {L}^{m}(\Omega ),v=0\text{sur}{\Gamma }_{u},L(v)=1\right\}\)

In this space, we define the regulated functional \({S}_{e}^{m}(v)\):

\({S}_{e}^{m}(v)={\int }_{\Omega }\frac{{k}^{1\text{-}m}}{m}\pi {(\varepsilon (v))}^{m}d\Omega -{L}_{0}(v)\)

The minimization problem \(\underset{v\in {V}_{a}^{\mathrm{m1}}}{\text{Inf}}\left[{S}_{e}^{m}(v)\right]\) is well posed thanks to the properties of \({L}^{m}(O)\) spaces and admits a unique solution \({u}_{m}\), for which the value reached by the

_images/Object_88.svg

is precisely \({\lambda }_{m}\). We denote the space of the incompressible fields of \({V}_{a}^{m1}\) by:

\(\tilde{{V}_{a}^{m1}}=\left\{v\in {V}_{a}^{m1}\text{tel}\text{que}\text{div}v=0\right\}\)

We then show that this problem can also be written in the form of looking for saddle point \(({\lambda }_{m},{u}_{m})\) of the following Lagrangian:

\(\underset{l\in \text{IR}}{\text{Max}}\left[\underset{\mathrm{v}\in \stackrel{~}{{V}_{a}^{m}}}{\text{Inf}}\phantom{\rule{0.5em}{0ex}}\left\{{\int }_{\mathrm{\Omega }}\frac{A\left(m\right)}{m}\text{.}{\left(\sqrt{\mathrm{\epsilon }(\mathrm{v})\text{.}\mathrm{\epsilon }(\mathrm{v})}\right)}^{m}d\mathrm{\Omega }-{L}_{0}(\mathrm{v})-l\left(L(\mathrm{v})-1\right)\right\}\right]\) eq 1.3-2

with: \(A(m)={k}^{1\text{-}m}{(\frac{2}{3})}^{m/2}\text{.}{\sigma }_{y}^{m}={\sigma }_{y}^{2\text{-}m}\text{.}{(3\mu )}^{m\text{-}1}\text{.}{(\frac{2}{3})}^{m/2}\).

Note that \(A(m)\) is increasing with \(m\) (if \(E\ge {\sigma }_{y}\), which is the case in practice) and homogeneous to a constraint, and remains bounded when \(m\to {1}^{+}\). If we choose \(k={\sigma }_{y}\), then \(A(m)={\sigma }_{y}{(\frac{2}{3})}^{m/2}\) and we find the incompressible elastic problem when \(m=2\), if we choose a Young’s modulus \(E={\sigma }_{y}\).

We can therefore see that this potential [éq 1.3-2] defines a law of behavior giving the stress tensor \(\sigma (u)\) by the Norton-Hoff behavior relationship, see [§2.1].

A decreasing sequence of \({\lambda }_{m}\) is thus constructed and the limit load \({\lambda }_{\text{lim}}\) is the limit of this sequence when \(m\to {1}^{+}\) (i.e. \(n\to \text{+}\infty\)):

\({\lambda }_{\text{lim}}=\underset{m\to 1}{\text{lim}}(\underset{v\in {V}_{a}^{\mathrm{m1}}}{\text{Inf}}\left[{S}_{e}^{m}(v)\right])=\underset{m\to 1}{\text{lim}}({S}_{e}^{m}({u}_{m}))\) eq 1.3-3

For the demonstration we will refer to [bib4] and [bib7].

Note:

If we amplify the intensity of the loading \(L\to \beta L\) (while we do not consider a permanent load \({L}_{0}=0\) ), the solutions depend on the factor \(\beta\) according to the following relationships:

\({u}_{m}(\beta )={\beta }^{\text{-}1}{u}_{m}(1);{\sigma }^{D}({u}_{m}(\beta ))={\beta }^{1\text{-}m}{\sigma }^{D}({u}_{m}(1));{S}_{e}^{m}({u}_{m}(\beta ))={\beta }^{\text{-}m}{S}_{e}^{m}({u}_{m}(1))\) .

At the convergence for \(m\to {1}^{+}\) , the conclusion given by the solution \({u}_{m}(\beta )\) is indeed the same as that given by \({u}_{m}(1)\), since \({\lambda }_{\text{lim}}(\beta )={\lambda }_{\text{lim}}(1)/\beta\) .