4. Sensitivity calculations#

The sensitivity analysis focuses only on the associated version of the formulation described in chapter 2.2.

4.1. Sensitivity to material data#

4.1.1. The direct problem#

We place ourselves in this part in the framework of the resolution of non-linear calculations.

The nonlinear static calculation is solved incrementally. It therefore requires the resolution of the non-linear equation system at each load step \(i\in \{1,I\}\):

\(\mathrm{\{}\begin{array}{c}R({u}_{i},{t}_{i})+{B}^{t}{\lambda }_{i}\mathrm{=}{L}_{i}\\ {\text{Bu}}_{i}\mathrm{=}{u}_{i}^{d}\end{array}\) eq 4.1.1-1

with

\((\mathrm{R}({\mathrm{u}}_{\mathrm{i}},{t}_{i}){)}_{k}\mathrm{=}{\mathrm{\int }}_{\Omega }\sigma ({\mathrm{u}}_{\mathrm{i}})\mathrm{:}\varepsilon ({\mathrm{w}}_{\mathrm{k}})d\Omega\) eq 4.1.1-2

\({\mathrm{w}}_{\mathrm{k}}\) is the shape function of the \(k\) th degree of freedom of the modelled structure,
\((\mathrm{R}({\mathrm{u}}_{\mathrm{i}},{t}_{i}))\) is the vector of nodal forces.

The resolution of this system is done by the Newton-Raphson method:

\(\{\begin{array}{ccc}{K}_{i}^{n}{\mathrm{\delta u}}_{i}^{n+1}+{B}^{t}{\mathrm{\delta \lambda }}_{i}^{n+1}& & {L}_{i}-R({u}_{i}^{n},{t}_{i})+{B}^{t}{\lambda }_{i}^{n}\\ {\mathrm{B\delta u}}_{i}^{n+1}& & -{\mathrm{Bu}}_{i-1}^{n}\end{array}\) eq 4.1.1-3

where \({K}_{i}^{n}=\frac{\partial R}{\partial u}{\mid }_{({u}_{i}^{n},{t}_{i})}\) is the matrix tangent to the load step \(i\) and to the Newton iteration \(n\).

The solution is therefore given by:

\(\{\begin{array}{ccc}{u}_{i}& & {u}_{i-1}+\sum _{n=0}^{N}{\mathrm{\delta u}}_{i}^{n}\\ {\lambda }_{i}& & {\lambda }_{i-1}+\sum _{n=0}^{N}{\mathrm{\delta \lambda }}_{i}^{n}\end{array}\) eq 4.1.1-4

with \(N\), the number of Newton iterations that were required to reach convergence.

4.1.2. Derivative calculus#

4.1.2.1. Foreplay#

As part of the sensitivity calculation, it is necessary to emphasize the dependencies of one quantity in relation to others. We will thus explain that the results of the previous calculation depend on a given parameter \(\mathrm{\Phi }\) (Young’s modulus, elastic limit, density,…) in the following way:

\({u}_{i}\mathrm{=}{u}_{i}(\Phi )\), \({\lambda }_{i}\mathrm{=}{\lambda }_{i}(\Phi )\).

But that is not enough. So we place ourselves in the framework of an incremental calculation with a Drucker-Prager law of behavior. If we consider the inter-dependencies of the parameters at an algorithmic level, we can write:

\(R=R({\sigma }_{i-1}(\Phi ),{p}_{i-1}(\Phi ),\Delta u(\Phi ))\)

\({\sigma }_{i}={\sigma }_{i-1}(\Phi )+\Delta \sigma ({\sigma }_{i-1}(\Phi ),{p}_{i-1}(\Phi ),\Delta u(\Phi ),\Phi )\)

\({p}_{i}={p}_{i-1}(\Phi )+\mathrm{\Delta p}({\sigma }_{i-1}(\Phi ),{p}_{i-1}(\Phi ),\Delta u(\Phi ),\Phi )\)

Where \(\Delta u\) is the displacement increment at convergence at load step \(i\).

Let’s specify the meaning of the notations we’ll use for derivatives:

:math:`frac{partial X}{partial Y}` refers to the**explicit* partial derivative of \(X\) with respect to \(Y\),

:math:`X{,}_{Y}` refers to the**total* variation of \(X\) compared to \(Y\).

4.1.2.2. Equilibrium derivation#

Given the previous remarks, let’s express the total variation of [éq 2.1-1] compared to \(\Phi\):

\(\{\begin{array}{ccc}\frac{\partial R}{\partial \Phi }+\frac{\partial R}{\partial \mathrm{\Delta u}}\cdot \mathrm{\Delta u}{,}_{\Phi }+\frac{\partial R}{\partial {\sigma }_{i-1}}\cdot {\sigma }_{i-1}{,}_{\Phi }+\frac{\partial R}{\partial {p}_{i-1}}\cdot {p}_{i-1}{,}_{\Phi }+{B}^{t}{\lambda }_{i}{,}_{\Phi }& & 0\\ \mathrm{B\Delta u}{,}_{\Phi }& & -{\text{Bu}}_{i-\mathrm{1,}\Phi }\end{array}\) eq 4.1.2.2-1

Note that here \(\frac{\partial R}{\partial \Phi }=0\): \(R\) does not depend explicitly on \(\Phi\) but implicitly as we will see in detail later.

That is:

\(\{\begin{array}{ccc}{K}_{i}^{N}\mathrm{\Delta u}{,}_{\Phi }+{B}^{t}{\lambda }_{i}{,}_{\Phi }& & -R{,}_{\Phi }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}\\ \mathrm{B\Delta u}{,}_{\Phi }& & -{\text{Bu}}_{i-\mathrm{1,}\Phi }\end{array}\) eq 4.1.2.2-2

Where

\({K}_{i}^{N}\) is the last tangent matrix used to reach convergence in Newton iterations,
\(R{,}_{\Phi }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}\) is the total change in \(R\), not taking into account the dependence of \(\mathrm{\Delta u}\) on \(\Phi\).

The problem now lies in calculating \(R{,}_{\Phi }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}\).

Note:

In [éq 4.1.2.2-2], we used the fact that \({K}_{i}^{N}=\frac{\partial R({u}_{i},{t}_{i})}{\partial \Delta u}\) while in [éq 4.1.1-3] we defined it by \({K}_{i}^{N}=\frac{\partial R({u}_{i},{t}_{i})}{\partial {u}_{i}^{N}}\). We have true equivalence of these two definitions inas insofar as \({u}_{i}={u}_{i-1}+\mathrm{\Delta u}\) and \(R\) actually depends on \(\Delta u\) (and also of course on \({\sigma }_{i-1}\) and \({p}_{i-1}\) ) .

Note:

If we derive with \(\Phi\) directly [éq 4.1.1-3], we find \({K}^{n}=\frac{\partial {u}^{n+1}}{\partial \Phi }+{B}^{t}\lambda ,\Phi =-R{,}_{\Phi /\mathrm{\Delta u}\ne \mathrm{\Delta u}/\Phi }-{K}^{n}{,}_{\Phi }{\mathrm{\delta u}}^{n+1}\). Which is the same **converging and makes it appear that the error on* \(\frac{\partial u}{\partial \Phi }\) depends on \({K}^{-1}K{,}_{\Phi }\) .

4.1.2.3. Calculation of the derivative of the law of behavior#

In the future, for the sake of clarity, we’ll drop the \(i-1\) hints.

According to [éq 4.1.1-2], \(R{,}_{\Phi }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}\) can be rewritten as:

math:: `mathrm {R} {,} {,} _ {Phi} {|} |} _ {\ mathit {\ Delta u}\ne\ mathit {\ Delta u} (\ Phi)} = {\ int} _ {\ int} _ _ {\ int} _ _ {\ Phi} _ {\ Phi} _ {\ Phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {\ phi} _ {} {,} _ {\ Phi} {|} |} _ {Deltamathrm {u}neDeltamathrm {u} (Phi)}right)mathrm {:}mathrm {epsilon}} _ {mathrm {w}}emathrm {omega}} eq 4.1.2.epsilon} ({epsilon}} ({mathrm {w}}}) _ {mathrm {k}}) dmathrm {Omega}} eq 4.1.2.epsilon} ({epsilon}} _ {mathrm {w}}} 3-1

So we need to calculate \(\Delta \sigma {,}_{\Phi }{\mid }_{\Delta u\ne \Delta u(\Phi )}\). To do this, we will use the expressions that are involved in the numerical integration of the law of behavior.

4.1.2.4. The case of linear elasticity#

In the context of linear elasticity, the law of behavior is expressed by:

\(\{\begin{array}{}\Delta \tilde{\sigma }=\mathrm{2\mu }\text{.}\tilde{\epsilon }(\Delta u)\\ \text{Tr}(\Delta \sigma )=\mathrm{3K}\text{.}\text{Tr}(\epsilon (\Delta u))\end{array}\)

or else:

\(\Delta \sigma =\mathrm{2\mu }\text{.}\tilde{\epsilon }(\Delta u)+K\text{.}\text{Tr}(\epsilon (\Delta u))\text{.}\text{Id}\) eq 4.1.2.4-1

where \(\text{Id}\) is the 2nd order identity tensor.

So, calculating the total variation of [éq 4.1.2.4-1] with respect to \(\Phi\) we get:

\(\Delta \sigma {,}_{\Phi }=\mathrm{2\mu }{,}_{\Phi }\text{.}\tilde{\epsilon }(\Delta u)+K{,}_{\Phi }\text{.}\text{Tr}(\epsilon (\Delta u))\text{.}\text{Id}+\mathrm{2\mu }\text{.}\tilde{\epsilon }(\Delta u{,}_{\Phi })+K\text{.}\text{Tr}(\epsilon (\Delta u{,}_{\Phi }))\text{.}\text{Id}\) eq 4.1.2.4-2

That is:

\(\Delta \sigma {,}_{\Phi }{\mid }_{\Delta u\ne \Delta u(\Phi )}=\mathrm{2\mu }{,}_{\Phi }\text{.}\tilde{\epsilon }(\Delta u)+K{,}_{\Phi }\text{.}\text{Tr}(\epsilon (\Delta u))\text{.}\text{Id}\) eq 4.1.2.4-3

4.1.2.5. Case of elastoplasticity of the Drucker-Prager type#

The Drucker-Prager law of behavior is written as:

\(\{\begin{array}{ccc}\epsilon (\mathrm{\Delta u})-S:\sigma & & \frac{3}{2}\cdot \mathrm{\Delta p}\cdot \frac{\tilde{\sigma }+\Delta \tilde{\sigma }}{(\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}}\\ (\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}+A\cdot \text{Tr}(\sigma +\mathrm{\Delta \sigma })& & R(p+\mathrm{\Delta p})\end{array}\) eq 4.1.2.5-1

where \(S\) is the elastic flexibility tensor and \(R\) is the plasticity criterion defined by:

in the case of linear work hardening:

\(\begin{array}{}R(p)=h\cdot p+{\sigma }^{y}\text{pour}0\le p\le {p}_{\text{ultm}}\\ R(p)=h\cdot {p}_{\text{ultm}}\text{pour}p\ge {p}_{\text{ultm}}\end{array}\)

in the case of parabolic work hardening:

\(\begin{array}{}R(p)={\sigma }^{y}\cdot (\text{1-}(\text{1-}\sqrt{\frac{{\sigma }_{\text{ultm}}^{y}}{{\sigma }^{y}}})\cdot \frac{p}{{p}_{\text{ultm}}}{)}^{2}\text{pour}0\le p\le {p}_{\text{ultm}}\\ R(p)={\sigma }_{\text{ultm}}^{y}\text{pour}p\ge {p}_{\text{ultm}}\end{array}\)

In numerical terms, this law of behavior is integrated using a radial feedback algorithm: an elastic prediction is made (noted \({\sigma }^{e}\)) that is corrected if the threshold is violated. So we write:

\(\{\begin{array}{ccc}\Delta \tilde{\sigma }& & \mathrm{2\mu }\text{.}\tilde{\epsilon }(\mathrm{\Delta u})-\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}}{{\sigma }_{\text{eq}}^{e}}\cdot {\tilde{\sigma }}^{e}\\ \text{Tr}(\mathrm{\Delta \sigma })& & \mathrm{3K}\cdot \text{Tr}(\epsilon (\mathrm{\Delta u}))-\mathrm{9K}\cdot A\cdot \mathrm{\Delta p}\\ \mathrm{\Delta p}& & \text{solution de}{\sigma }_{\text{eq}}^{e}-(\mathrm{3\mu }+\mathrm{9K}\cdot {A}^{2})\cdot \mathrm{\Delta p}+A\cdot \text{Tr}({\sigma }^{e})-R({p}^{-}+\mathrm{\Delta p})=0\end{array}\) eq 4.1.2.5-2

We are going to distinguish between two cases.

1st case: \(\mathrm{\Delta p}=0\)

This means that during the present no load, the Gauss point in question did not see an increase in its plasticization. We then find ourselves in the case of linear elasticity:

2nd case: \(\mathrm{\Delta p}>0\)

Given the dependencies between variables in [éq 4.1.2.5-1], we can write:

\(\{\begin{array}{ccc}\mathrm{\Delta \sigma }{,}_{\Phi }& & \frac{\partial \mathrm{\Delta \sigma }}{\partial \Phi }+\frac{\partial \mathrm{\Delta \sigma }}{\partial \sigma }\cdot \sigma {,}_{\Phi }+\frac{\partial \mathrm{\Delta \sigma }}{\partial p}\cdot p{,}_{\Phi }+\frac{\partial \mathrm{\Delta \sigma }}{\partial \epsilon (\mathrm{\Delta u})}\cdot \epsilon (\mathrm{\Delta u}{,}_{\Phi })\\ \mathrm{\Delta p}{,}_{\Phi }& & \frac{\partial \mathrm{\Delta p}}{\partial \Phi }+\frac{\partial \mathrm{\Delta p}}{\partial \sigma }\cdot \sigma {,}_{\Phi }+\frac{\partial \mathrm{\Delta p}}{\partial p}\cdot p{,}_{\Phi }+\frac{\partial \mathrm{\Delta p}}{\partial \epsilon (\mathrm{\Delta u})}\cdot \epsilon (\mathrm{\Delta u}){,}_{\Phi }\end{array}\) eq 4.1.2.5-4

In addition, in accordance with the algorithmic integration of the law, we will separate deviatoric and hydrostatic parts.

\(\{\begin{array}{ccc}\mathrm{\Delta \sigma }{,}_{\Phi }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}& & \frac{\partial \Delta \tilde{\sigma }}{\partial \Phi }+\frac{1}{3}\cdot \frac{\partial \text{Tr}(\mathrm{\Delta \sigma })}{\partial \Phi }\cdot \text{Id}\\ & +& \frac{\partial \Delta \tilde{\sigma }}{\partial \sigma }\cdot \sigma {,}_{\Phi }+\frac{1}{3}\cdot \frac{\partial \text{Tr}(\mathrm{\Delta \sigma })}{\partial \sigma }\cdot \text{Id}\cdot \sigma {,}_{\Phi }\\ & +& \frac{\partial \Delta \tilde{\sigma }}{\partial p}\cdot p{,}_{\Phi }+\frac{1}{3}\cdot \frac{\partial \text{Tr}(\mathrm{\Delta \sigma })}{\partial p}\cdot \text{Id}\cdot p{,}_{\Phi }\\ \mathrm{\Delta p}{,}_{\Phi }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}& & \frac{\partial \mathrm{\Delta p}}{\partial \Phi }+\frac{\partial \mathrm{\Delta p}}{\partial \sigma }\cdot \sigma {,}_{\Phi }+\frac{\partial \mathrm{\Delta p}}{\partial p}\cdot p{,}_{\Phi }\end{array}\) eq 4.1.2.5-5

And so, we calculate:

\(\frac{\partial \Delta \sigma }{\partial \Phi }\)

\(\frac{\partial \Delta \tilde{\sigma }}{\partial \Phi }=\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon }(\mathrm{\Delta u})-\frac{\partial \mathrm{3\mu }}{\partial \Phi }\cdot \frac{\mathrm{\Delta p}}{{\sigma }_{\text{eq}}^{e}}\cdot {\tilde{\sigma }}^{e}-\mathrm{3\mu }\cdot \frac{\frac{\partial \mathrm{\Delta p}}{\partial \Phi }}{{\sigma }_{\text{eq}}^{e}}\cdot {\tilde{\sigma }}^{e}+\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}\cdot \frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial \Phi }}{{\sigma }_{\text{eq}}^{{e}^{2}}}\cdot {\tilde{\sigma }}^{e}-\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}}{{\sigma }_{\text{eq}}^{e}}\cdot \frac{\partial {\tilde{\sigma }}^{e}}{\partial \Phi }\)

\(\frac{\partial \text{Tr}(\Delta \sigma )}{\partial \Phi }=\frac{\partial \mathrm{3K}}{\partial \Phi }\cdot \text{Tr}(\epsilon (\Delta u))-\frac{\partial \mathrm{9K}}{\partial \Phi }\cdot {\rm A}\cdot \Delta p-\mathrm{9K}\cdot \frac{\partial {\rm A}}{\partial \Phi }\cdot \mathrm{\Delta p}-\mathrm{9K}\cdot {\rm A}\cdot \frac{\partial \mathrm{\Delta p}}{\partial \Phi }\)

\(\frac{\partial \Delta \sigma }{\partial \sigma }\)

\(\frac{\partial \Delta \tilde{\sigma }}{\partial \sigma }=-\mathrm{3\mu }\cdot \frac{\frac{\partial \mathrm{\Delta p}}{\partial \sigma }}{{\sigma }_{\text{eq}}^{e}}\otimes {\tilde{\sigma }}^{e}+\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}}{{\sigma }_{\text{eq}}^{{e}^{2}}}\cdot \frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial \sigma }\otimes {\tilde{\sigma }}^{e}-\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}}{{\sigma }_{\text{eq}}^{e}}\cdot J\)

where \(J\) is the deviatory operator defined by: \(J:\sigma =\tilde{\sigma }\)

\(\frac{\partial \text{Tr}(\Delta \sigma )}{\partial \sigma }=-\mathrm{9K}\cdot {\rm A}\cdot \frac{\partial \mathrm{\Delta p}}{\partial \sigma }\)

\(\frac{\partial \mathrm{\Delta \sigma }}{\partial p}\)

\(\frac{\partial \Delta \tilde{\sigma }}{\partial p}=-\frac{\mathrm{3\mu }}{{\sigma }_{\text{eq}}^{e}}\cdot \frac{\partial \mathrm{\Delta p}}{\partial p}\cdot {\tilde{\sigma }}^{e}\)

\(\frac{\partial \text{Tr}(\mathrm{\Delta \sigma })}{\partial p}=-9\cdot K\cdot {\rm A}\cdot \frac{\partial \mathrm{\Delta p}}{\partial p}\)

\(\mathrm{\Delta p}{,}_{\Phi }\)

We use the fact that: \((\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}=(\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}^{e}-3\mu \cdot \mathit{\Delta p}\)

\({\mathrm{\Delta p}}_{,\Phi }=\frac{1}{\mathrm{3\mu }}\cdot ((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\Phi }}-(\sigma +\mathrm{\Delta \sigma }{)}_{{\text{eq}}_{,\Phi }}-\frac{\partial \mathrm{3\mu }}{\partial \Phi }\cdot \mathrm{\Delta p})\)

Note:

In these calculations the following results were or should be used:
\(\frac{\partial {\tilde{\sigma }}^{e}}{\partial \Phi }=\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon }(\mathrm{\Delta u})\)	\(\frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial \Phi }=\frac{3}{2}\cdot \frac{(\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon }(\mathrm{\Delta u})):(\tilde{\sigma }+\mathrm{2\mu }\text{.}\tilde{\epsilon }(\mathrm{\Delta u}))}{{\sigma }_{\text{eq}}^{e}}\)
Order tensor 2	Scalar
\(\frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial \sigma }=\frac{3}{2}\cdot \frac{{\tilde{\sigma }}^{e}}{{\sigma }_{\text{eq}}^{e}}\)	\(\frac{\partial {\tilde{\sigma }}^{e}}{\partial \sigma }=J\)
Order tensor 2	4 order tensor
\(\frac{\partial \text{Tr}({\sigma }^{e})}{\partial \sigma }=\text{Id}\)	\(\frac{\partial \text{Tr}({\sigma }^{e})}{\partial \Phi }=\frac{\partial \mathrm{3K}}{\partial \Phi }\cdot \text{Tr}(\epsilon (\mathrm{\Delta u}))\)
Order tensor 2	Scalar

It is also necessary to calculate the partial derivatives of the cumulative plastic deformation increment with respect to the material parameters, the stresses and the cumulative plastic deformation (see Appendix)

These are obtained by deriving the equation solved to calculate the cumulative plastic deformation increment during the direct calculation.

4.1.2.6. Calculating the derivative of displacement#

Once: math: Deltamathrm {mathrm {sigma} {sigma} {,} {sigma} {,} {,} _ {deltamathrm {u} (Phi)}} has been calculated, we can create the second member:math: R {,} _ {Phi} {mid} _ {Delta unePhiePhi)} `calculated, we can create the second member:math: `R {,}} _ {Phi} {mid} _ {Delta uneDelta unePhi)} calculated, we can create the second member:math: R {,}} _ {Phi} {mid} _ {Delta unePhi)} `calculated} using [éq 4.1.2.3-1]. We then solve the [éq 4.1.2.2-2] system and we get the displacement increment derived with respect to:math: Phi.

4.1.2.7. Calculation of the derivative of other quantities#

Now that we have \(\Delta u{,}_{\Phi }\), we need to calculate the derivative of the other quantities. Two more cases are separated:

linear elasticity

According to [éq 4.1.2.5-1], the derivative of the stress increment is calculated as follows:

\(\Delta \sigma {,}_{\Phi }=\Delta \sigma {,}_{\Phi }{\mid }_{\Delta u\ne \Delta u(\Phi )}+\mathrm{2\mu }\text{.}\tilde{\epsilon }(\Delta u{,}_{\Phi })+K\text{.}\text{Tr}(\epsilon (\Delta u{,}_{\Phi }))\text{.}\text{Id}\)

The cumulative plastic deformation increment, for its part, does not see any evolution:

\(\mathrm{\Delta p}{,}_{\Phi }=0\)

Drucker-Prager elastoplasticity

If \(\mathrm{\Delta p}=0\), we find the previous case.

Otherwise, we get from [éq 4.1.2.5-2]:

\(\Delta \sigma {,}_{\Phi }=\Delta \sigma {,}_{\Phi }{\mid }_{\Delta u\ne \Delta u(\Phi )}+\frac{\partial \Delta \sigma }{\partial \epsilon (\Delta u)}:\epsilon (\Delta u{,}_{\Phi })\)

And for the cumulative plastic deformation, the following relationship is used:

\((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}=(\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{e}-\mathrm{3\mu }\cdot \mathrm{\Delta p}\)

This allows us to write that:

The equivalent sensitive stresses are calculated as follows:

\((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\Phi }}=\frac{3}{2(\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{e}}\cdot ({\tilde{\sigma }}_{,\Phi }+\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon }(\mathrm{\Delta u})+\mathrm{2\mu }\cdot \tilde{\epsilon }({\mathrm{\Delta u}}_{,\Phi })):(\tilde{\sigma }+\mathrm{2\mu }\cdot \tilde{\epsilon }(\mathrm{\Delta u}))\)

\((\sigma +\mathrm{\Delta \sigma }{)}_{{\text{eq}}_{,\Phi }}=\frac{3}{2(\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}}\cdot ({\tilde{\sigma }}_{,\Phi }+\Delta {\tilde{\sigma }}_{,\Phi }):(\tilde{\sigma }+\Delta \tilde{\sigma })\)

Once all these calculations are complete, we update all the derived quantities and we move on to the next load step.

4.1.2.8. Synthesis#

To summarize the preceding paragraphs, the various stages of the calculation are represented by the following diagram:

4.2. Load sensitivity#

The approach here is quite similar to that in the previous paragraph. However, we develop it entirely for the sake of clarity, so that this paragraph can be read independently.

4.2.1. The direct problem: expression of loading#

So far we have expressed the direct problem in the form of:

\(\{\begin{array}{ccc}R({u}_{i},{t}_{i})+{B}^{t}{\lambda }_{i}& & {L}_{i}\\ {\text{Bu}}_{i}& & {u}_{i}^{d}\end{array}\) eq 4.2.1-1

The loads are collected at the second member and include imposed forces \({L}_{i}\) and enforced displacements \({u}_{i}^{d}\).

Assume that the forced load \({L}_{i}\) depends on a scalar parameter \(\alpha\) in the following way:

\({L}_{i}(\alpha )={L}_{i}^{1}+{L}_{i}^{2}(\alpha )\) eq 4.2.1-2

Where

\({L}_{i}^{1}\) is a vector independent of \(\alpha\),
\({L}_{i}^{2}\) depends linearly on \(\alpha\).

It is desired to calculate the sensitivity of the results of the direct calculation to a variation of the parameter \(\alpha\).

4.2.2. The derived problem#

4.2.2.1. Equilibrium derivation#

As in the previous chapter, taking into account the dependencies between the different fields, we derive the balance [éq 4.2.1-1] with respect to \(\alpha\):

\(\{\begin{array}{ccc}\frac{\partial R}{\partial \alpha }+\frac{\partial R}{\partial \mathrm{\Delta u}}\cdot \mathrm{\Delta u}{,}_{\alpha }+\frac{\partial R}{\partial {\sigma }_{i-1}}\cdot {\sigma }_{i-1}{,}_{\alpha }+\frac{\partial R}{\partial {p}_{i-1}}\cdot {p}_{i-1}{,}_{\alpha }+{B}^{t}{\lambda }_{i}{,}_{\alpha }& & {L}_{i}^{2}(1)\\ \mathrm{B\Delta u}{,}_{\alpha }& & -{\text{Bu}}_{\iota -\mathrm{1,}\alpha }\end{array}\) eq 4.2.2.1-1

We used the fact that \({L}_{i}^{2}\) depends linearly on \(\alpha\).

That is:

\(\{\begin{array}{ccc}{K}_{i}^{N}\mathrm{\Delta u}{,}_{\alpha }+{B}^{t}{\lambda }_{i}{,}_{\alpha }& & {L}_{i}^{2}(1)-R{,}_{\alpha }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\alpha )}\\ \mathrm{B\Delta u}{,}_{\alpha }& & -{\text{Bu}}_{\iota -\mathrm{1,}\alpha }\end{array}\) eq 4.2.2.1-2

Where

\({K}_{i}^{N}\) is the last tangent matrix used to reach convergence in Newton iterations,
\(R{,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}\) is the total change in \(R\), not taking into account the dependence of \(\Delta u\) on \(\alpha\).

As before, the problem lies in calculating \(R{,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}\).

4.2.2.2. Calculation of the derivative of the law of behavior#

According to [éq 4.1.1-2], \(R{,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}\) can be rewritten as:

\(R{,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}={\int }_{\Omega }^{}(\sigma {,}_{\alpha }+\Delta \sigma {,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}):\epsilon ({w}_{k})d\Omega\) eq 4.2.2.2-1

To do this, we will use the expressions that are involved in the numerical integration of the law of behavior to calculate \(\Delta \sigma {,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}\).

4.2.2.3. The case of linear elasticity#

In the context of linear elasticity, the law of behavior is expressed by:

\(\Delta \sigma =\mathrm{2\mu }\text{.}\tilde{\epsilon }(\Delta u)+K\text{.}\text{Tr}(\epsilon (\Delta u))\text{.}\text{Id}\) eq 4.2.2.3-1

where \(\text{Id}\) is the 2nd order identity tensor.

So, calculating the total variation of [éq 4.2.2.3-1] with respect to \(\alpha\) we get:

\(\begin{array}{c}\mathit{\Delta \sigma }{,}_{\alpha }\mathrm{=}\mathrm{2\mu }{,}_{\alpha }\text{.}\tilde{\epsilon }(\mathit{\Delta u})+K{,}_{\alpha }\text{.}\text{Tr}(\epsilon (\mathit{\Delta u}))\text{.}\text{Id}+\mathrm{2\mu }\text{.}\tilde{\epsilon }(\mathit{\Delta u}{,}_{\alpha })+K\text{.}\text{Tr}(\epsilon (\mathit{\Delta u}{,}_{\alpha }))\text{.}\text{Id}\\ \mathrm{=}0\text{.}+0\text{.}+\mathrm{2\mu }\text{.}\tilde{\epsilon }(\mathit{\Delta u}{,}_{\alpha })+K\text{.}\text{Tr}(\epsilon (\mathit{\Delta u}{,}_{\alpha }))\text{.}\text{Id}\end{array}\) eq 4.2.2.3-2

That is:

\(\Delta \sigma {,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}=0\text{.}\)

4.2.2.4. Case of elastoplasticity of the Drucker-Prager type#

As before, we will distinguish between two cases.

1st case: \(\mathrm{\Delta p}=0\)

This means that during the present no load, the Gauss point in question did not see an increase in its plasticization. We then find ourselves in the case of linear elasticity:

\(\Delta \sigma {,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}=0\text{.}\)

2nd case: \(\mathrm{\Delta p}>0\)

Taking into account the dependencies between variables, we can write:

\(\mathrm{\{}\begin{array}{ccc}\mathit{\Delta \sigma }{,}_{\alpha }& \mathrm{=}& \frac{\mathrm{\partial }\mathit{\Delta \sigma }}{\mathrm{\partial }\alpha }+\frac{\mathrm{\partial }\mathit{\Delta \sigma }}{\mathrm{\partial }\sigma }\mathrm{\cdot }\sigma {,}_{\alpha }+\frac{\mathrm{\partial }\mathit{\Delta \sigma }}{\mathrm{\partial }p}\mathrm{\cdot }p{,}_{\alpha }+\frac{\mathrm{\partial }\mathit{\Delta \sigma }}{\mathrm{\partial }\epsilon (\mathit{\Delta u})}\mathrm{\cdot }\epsilon (\mathit{\Delta u}{,}_{\alpha })\\ \mathit{\Delta p}{,}_{\alpha }& \mathrm{=}& \frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }\alpha }+\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }\sigma }\mathrm{\cdot }\sigma {,}_{\alpha }+\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }p}\mathrm{\cdot }p{,}_{\alpha }+\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }\epsilon (\mathit{\Delta u})}\mathrm{\cdot }\epsilon (\mathit{\Delta u}){,}_{\alpha }\end{array}\)

In addition, in accordance with the algorithmic integration of the law, we will separate deviatoric and hydrostatic parts.

\(\mathrm{\{}\begin{array}{ccc}\mathit{\Delta \sigma }{,}_{\alpha }{\mathrm{\mid }}_{\mathit{\Delta u}\mathrm{\ne }\mathit{\Delta u}(\alpha )}& \mathrm{=}& \frac{\mathrm{\partial }\Delta \tilde{\sigma }}{\mathrm{\partial }\alpha }+\frac{1}{3}\mathrm{\cdot }\frac{\mathrm{\partial }\text{Tr}(\mathit{\Delta \sigma })}{\mathrm{\partial }\alpha }\mathrm{\cdot }\text{Id}\\ & +& \frac{\mathrm{\partial }\Delta \tilde{\sigma }}{\mathrm{\partial }\sigma }\mathrm{\cdot }\sigma {,}_{\alpha }+\frac{1}{3}\mathrm{\cdot }\frac{\mathrm{\partial }\text{Tr}(\mathit{\Delta \sigma })}{\mathrm{\partial }\sigma }\mathrm{\cdot }\text{Id}\mathrm{\cdot }\sigma {,}_{\alpha }\\ & +& \frac{\mathrm{\partial }\Delta \tilde{\sigma }}{\mathrm{\partial }p}\mathrm{\cdot }p{,}_{\alpha }+\frac{1}{3}\mathrm{\cdot }\frac{\mathrm{\partial }\text{Tr}(\mathit{\Delta \sigma })}{\mathrm{\partial }p}\mathrm{\cdot }\text{Id}\mathrm{\cdot }p{,}_{\alpha }\\ \mathit{\Delta p}{,}_{\alpha }{\mathrm{\mid }}_{\mathit{\Delta u}\mathrm{\ne }\mathit{\Delta u}(\alpha )}& \mathrm{=}& \frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }\alpha }+\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }\sigma }\mathrm{\cdot }\sigma {,}_{\alpha }+\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }p}\mathrm{\cdot }p{,}_{\alpha }\end{array}\)

And so, we calculate:

\(\frac{\mathrm{\partial }\Delta \mathrm{\sigma }}{\mathrm{\partial }\alpha }\)

Since we don’t have an explicit dependency of \(\Delta \mathrm{\sigma }\) on \(\alpha\), we get:

\(\frac{\mathrm{\partial }\Delta \tilde{\mathrm{\sigma }}}{\mathrm{\partial }\alpha }\mathrm{=}0\text{.}\)

\(\frac{\mathrm{\partial }\text{Tr}(\Delta \mathrm{\sigma })}{\mathrm{\partial }\alpha }\mathrm{=}0\text{.}\)

\(\frac{\mathrm{\partial }\Delta \mathrm{\sigma }}{\mathrm{\partial }\mathrm{\sigma }}\)

\(\frac{\mathrm{\partial }\Delta \tilde{\sigma }}{\mathrm{\partial }\sigma }\mathrm{=}\mathrm{-}\mathrm{3\mu }\mathrm{\cdot }\frac{\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }\sigma }}{{\sigma }_{\text{eq}}^{e}}\otimes {\tilde{\sigma }}^{e}+\mathrm{3\mu }\mathrm{\cdot }\frac{\mathit{\Delta p}}{{\sigma }_{\text{eq}}^{{e}^{2}}}\mathrm{\cdot }\frac{\mathrm{\partial }{\sigma }_{\text{eq}}^{e}}{\mathrm{\partial }\sigma }\otimes {\tilde{\sigma }}^{e}\mathrm{-}\mathrm{3\mu }\mathrm{\cdot }\frac{\mathit{\Delta p}}{{\sigma }_{\text{eq}}^{e}}\mathrm{\cdot }J\)

where \(\mathrm{J}\) is the deviatory operator.

\(\frac{\mathrm{\partial }\text{Tr}(\mathit{\Delta \sigma })}{\mathrm{\partial }\sigma }\mathrm{=}\mathrm{-}\mathrm{9K}\mathrm{\cdot }{\rm A}\mathrm{\cdot }\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }\sigma }\)

\(\frac{\mathrm{\partial }\Delta \mathrm{\sigma }}{\mathrm{\partial }p}\)

\(\frac{\mathrm{\partial }\Delta \tilde{\sigma }}{\mathrm{\partial }p}\mathrm{=}\mathrm{-}\frac{\mathrm{3\mu }}{{\sigma }_{\text{eq}}^{e}}\mathrm{\cdot }\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }p}\mathrm{\cdot }{\tilde{\sigma }}^{e}\)

\(\frac{\mathrm{\partial }\text{Tr}(\mathit{\Delta \sigma })}{\mathrm{\partial }p}\mathrm{=}\mathrm{-}9\mathrm{\cdot }K\mathrm{\cdot }{\rm A}\mathrm{\cdot }\frac{\mathrm{\partial }\mathit{\Delta p}}{\mathrm{\partial }p}\)

\(\mathit{\Delta p}{,}_{\alpha }\)

We use the fact that: \((\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}\mathrm{=}(\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}^{e}\mathrm{-}\mathrm{3\mu }\mathrm{\cdot }\mathit{\Delta p}\)

\({\mathit{\Delta p}}_{,\alpha }\mathrm{=}\frac{1}{\mathrm{3\mu }}\mathrm{\cdot }((\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\alpha }}\mathrm{-}(\sigma +\mathit{\Delta \sigma }{)}_{{\text{eq}}_{,\alpha }}\mathrm{-}\frac{\mathrm{\partial }\mathrm{3\mu }}{\mathrm{\partial }\alpha }\mathrm{\cdot }\mathit{\Delta p})\)

We will refer again to the remark at the end of [§4.1.2.5] for quantities whose calculation has not been detailed here.

4.2.2.5. Calculating the derivative of displacement#

Once \(\Delta \mathrm{\sigma }{,}_{\alpha }{\mathrm{\mid }}_{\Delta \mathrm{u}\mathrm{\ne }\Delta \mathrm{u}(\alpha )}\) has been calculated, we can form the second member \(R{,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}\). The [éq 4.2.2.1-1] system is then solved and the displacement increment derived with respect to \(\alpha\) is obtained.

4.2.2.6. Calculation of the derivative of other quantities#

Now that we have \(\Delta u{,}_{\alpha }\), we need to calculate the derivative of the other quantities. Two more cases are separated:

linear elasticity

According to [éq 4.2.2.3-1], the derivative of the stress increment is calculated as follows:

\(\Delta \sigma {,}_{\alpha }=0\text{.}+\mathrm{2\mu }\text{.}\tilde{\epsilon }(\Delta u{,}_{\alpha })+K\text{.}\text{Tr}(\epsilon (\Delta u{,}_{\alpha }))\text{.}\text{Id}\)

The cumulative plastic deformation increment, for its part, does not see any evolution:

\(\mathrm{\Delta p}{,}_{\alpha }=0\)

Drucker Prager type elastoplasticity

If \(\mathrm{\Delta p}=0\), we find the previous case.

If not, we get:

\(\mathrm{\Delta \sigma }{,}_{\alpha }=\mathrm{\Delta \sigma }{,}_{\alpha }{\mid }_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\alpha )}+\frac{\partial \mathrm{\Delta \sigma }}{\partial \epsilon (\mathrm{\Delta u})}:\epsilon (\mathrm{\Delta u}{,}_{\alpha })\)

And for the cumulative plastic deformation:

\({\mathrm{\Delta p}}_{,\alpha }=\frac{1}{\mathrm{3\mu }}\cdot ((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\alpha }}-(\sigma +\mathrm{\Delta \sigma }{)}_{{\text{eq}}_{,\alpha }}-\frac{\partial \mathrm{3\mu }}{\partial \alpha }\cdot \mathrm{\Delta p})\)

Once all these calculations are complete, we update all the derived quantities and we move on to the next load step.

4.2.2.7. Synthesis#

To summarize the preceding paragraphs, the various stages of the calculation are represented by the following diagram:

Assembly of \(R{,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}\)

Calculating terms \(\Delta \sigma {,}_{\alpha }{\mid }_{\Delta u\ne \Delta u(\alpha )}\)

Convergence of the load step \(n\) of the direct calculation

Passage to the step

Charging \(n+1\)

Calculation of

\(\Delta \sigma {,}_{\alpha }\) and \(\mathrm{\Delta p}{,}_{\alpha }\)

System resolution [éq4.2.2.1-1]

\(\Rightarrow \Delta u{,}_{\alpha }\)