Benchmark solution
=====================


Calculation method used for the reference solution
--------------------------------------------------------

Analytical solution was determined by F. Voldoire (EDF R&D/AMA):

**Axisymmetric case (** :math:`\mathrm{2D}` **)**

Move fields :math:`u={U}_{r}(r){e}_{r}` (blocked in :math:`z`)

Deformation field: :math:`\varepsilon (u)=(\begin{array}{ccc}{u}_{r}\text{'}& 0& 0\\ 0& 0& 0\\ 0& 0& \frac{{u}_{r}}{r}\end{array})` according to :math:`(\begin{array}{c}r\\ z\\ \theta \end{array})`

Constraint fields: :math:`\sigma (u)={\sigma }_{l}(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array})(\mathrm{cf.}\mathrm{conditions}\mathrm{aux}\mathrm{limites})` according to :math:`(\begin{array}{c}r\\ z\\ \theta \end{array})`

**Parallelepipedic case**

Move fields :math:`u={U}_{x}(x){e}_{x}+{U}_{y}(y){e}_{y}` (blocked in :math:`z`)

Deformation field: :math:`\varepsilon (u)=(\begin{array}{ccc}{u}_{x}\text{'}& 0& 0\\ 0& 0& 0\\ 0& 0& {u}_{y}\text{'}\end{array})` according to :math:`(\begin{array}{c}x\\ z\\ y\end{array})`

Constraint fields: :math:`\sigma (u)=\sigma (\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array})(\mathrm{cf.}\mathrm{conditions}\mathrm{aux}\mathrm{limites})` according to :math:`(\begin{array}{c}x\\ z\\ y\end{array})`

The case can be studied under plane constraints and in :math:`\mathrm{3D}`.


.. image:: images/100002000000018A0000012A11AB208A2C4730E1.png
    :width: 2.6138in
    :height: 1.9689in


.. _RefImage_100002000000018A0000012A11AB208A2C4730E1.png:

:math:`2\nu =\frac{E}{(1+\nu )}` :math:`K=\frac{E}{(1-2\nu )}`

**The law of behavior** can be written (scalar internal variable :math:`p`):

:math:`\varepsilon =\frac{1}{\mathrm{9K}}\mathrm{tr}\sigma \mathrm{Id}+\frac{1}{2\mu }{\sigma }^{D}+{\varepsilon }^{P}+\alpha (T–{T}^{0})\mathrm{Id}`

with: :math:`{\sigma }^{D}=\sigma –\frac{1}{3}\mathrm{tr}(\sigma )\mathrm{Id}` (constraint deviator)

:math:`{\dot{\sigma }}^{P}=\frac{3}{2}\dot{p}\frac{{\sigma }^{D}}{∥{\sigma }_{\mathrm{éq}}∥}\text{et}∥{\sigma }_{\mathrm{éq}}∥=\sqrt{\frac{3}{2}{\sigma }^{D}\mathrm{.}{\sigma }^{D}}`

:math:`\dot{p}=0` if :math:`f(\sigma ,p)=∥{\sigma }_{\mathrm{éq}}∥–R(p)<0`

:math:`\dot{p}\ge 0` if :math:`f(\sigma ,p)=0`

:math:`R(p)` refers to the work hardening function: :math:`R(p)={\sigma }_{y}+\frac{{\mathrm{EE}}_{T}}{E-{E}_{T}}p`

The :math:`\dot{p}` rate can be expressed, when :math:`f(\sigma ,p)=0`. Indeed, from :math:`\dot{p}f` that is identical to zero, we draw: :math:`\dot{p}\dot{f}+\ddot{p}f=0`

So, when we are on the criterion, necessarily :math:`\dot{f}=0`. That is to say:

 :math:`\frac{3}{2}\frac{{\sigma }^{D}\mathrm{.}{\dot{\sigma }}^{D}}{∥{\sigma }_{\mathit{éq}}∥}–\frac{\mathrm{\partial }R}{\mathrm{\partial }T}\mathrm{.}\dot{T}–\frac{\mathrm{\partial }R}{\mathrm{\partial }p}\dot{p}\mathrm{=}0`

 :math:`\frac{3}{2}\frac{{\sigma }^{D}\mathrm{.}{\dot{\sigma }}^{D}}{∥{\sigma }_{\mathit{éq}}∥}+{\sigma }_{y}^{o}s\dot{T}\mathrm{-}\frac{{\mathit{EE}}_{T}}{E\mathrm{-}{E}_{T}}\dot{p}\mathrm{=}0`

From where:

 :math:`\dot{p}=\frac{E-{E}_{T}}{{\mathrm{EE}}_{T}}(\frac{3}{2}\frac{{\sigma }^{D}{\dot{\sigma }}^{D}}{∥{\sigma }_{\mathrm{éq}}∥}+{\sigma }_{y}^{o}\mathrm{.}s\mathrm{.}\dot{T})` if :math:`\dot{p}\ge 0`, for :math:`∥{\sigma }_{\mathrm{éq}}∥=R(p)`

Since the stress field is uniaxial, we have:

 :math:`{\sigma }^{D}=\frac{{\sigma }_{L}}{3}(\begin{array}{ccc}-1& 0& 0\\ 0& 2& 0\\ 0& 0& -1\end{array})`

So:

 :math:`∥{\sigma }_{\mathrm{éq}}∥=∣{\sigma }_{L}∣`

and:

 :math:`{\dot{\varepsilon }}^{P}=\frac{\dot{p}}{2}\mathrm{.}\mathrm{sgn}({\sigma }_{L})\mathrm{.}(\begin{array}{ccc}-1& 0& 0\\ 0& 2& 0\\ 0& 0& -1\end{array})`

The behavioral relationship leads to:

 :math:`\{\begin{array}{c}\dot{{\varepsilon }_{\mathrm{rr}}}=\dot{{\varepsilon }_{\theta \theta }}=\frac{-\nu }{E}\dot{{\sigma }_{L}}-\frac{\dot{p}}{2}\mathrm{sgn}({\sigma }_{L})+\alpha \dot{T}(={\dot{\varepsilon }}_{\mathrm{xx}}={\dot{\varepsilon }}_{\mathrm{yy}}\text{pour le cas du parallélépipède})\\ \dot{{\varepsilon }_{\mathrm{rr}}}=\dot{{\varepsilon }_{\theta \theta }}=\frac{-\nu }{E}\dot{{\sigma }_{L}}-\frac{\dot{p}}{2}\mathrm{sgn}({\sigma }_{L})+\alpha \dot{T}\end{array}`

From where:



 :math:`\{\begin{array}{c}\dot{{\varepsilon }_{\mathrm{rr}}}=\dot{{\varepsilon }_{\theta \theta }}=\frac{-3}{2}\alpha \dot{T}+\frac{1-2\nu }{2E}\dot{{\sigma }_{L}}\\ \dot{p}=\mathrm{sgn}({\sigma }_{L})(-\alpha \dot{T}-\frac{\dot{{\sigma }_{L}}}{E})=0\text{si}∣{\sigma }_{L}∣<R(p)\\ \dot{p}=\mathrm{Max}[0,\frac{E-{E}_{T}}{{\mathrm{EE}}_{T}}(\frac{3}{2}\frac{{\sigma }^{D}\dot{{\sigma }^{D}}}{∥{\sigma }_{\mathrm{éq}}∥}+{\sigma }_{y}^{o}\mathrm{.}s\mathrm{.}\dot{T})]\text{sinon}\end{array}`

That is to say, in case :math:`∣{\sigma }_{L}∣=R(p)` (criterion met):

:math:`\dot{p}=\mathrm{max}\left[0,\frac{E-{E}_{T}}{{\mathrm{EE}}_{T}}(\frac{3}{2}\frac{{\sigma }^{D}\mathrm{.}{\dot{\sigma }}^{D}}{∥{\sigma }_{\mathrm{éq}}∥}+{\sigma }_{y}^{o}\mathrm{.}s\mathrm{.}\dot{T})\right]`

Elastic phase
~~~~~~~~~~~~~~~~

At the start of thermal loading, :math:`∣{\sigma }_{L}∣` being less than :math:`{\sigma }_{y}`, :math:`\dot{p}` is zero.

From where: :math:`\dot{{\sigma }_{L}}=-E\alpha \dot{T};\dot{{\varepsilon }_{\mathrm{rr}}}=\dot{{\varepsilon }_{\theta \theta }}=\alpha \dot{T}(1+\nu )`

So: :math:`\{\begin{array}{c}{\sigma }_{L}=-E\alpha \theta t\\ {\varepsilon }_{\mathrm{rr}}={\varepsilon }_{\theta \theta }=\alpha \theta (1+\nu )t\end{array}` (compression :math:`{\sigma }_{L}<0`)

     This corresponds to the reference solution for the orthotropic elasticity test

**Validity of the elastic solution**

The criterion is: :math:`∣{\sigma }_{L}(t)∣-{\sigma }_{y}(t)=E=\theta t-{\sigma }_{y}^{o}(1-s\theta t)\le 0`

The criterion is not met for :math:`t=[\mathrm{0,}{t}_{y}]`, with: :math:`{t}_{y}=\frac{{\sigma }_{y}^{o}}{\theta (E\alpha +{\sigma }_{y}^{o}s)}`


.. image:: images/10000200000001E60000012AD278620AEBD0A733.png
    :width: 2.5752in
    :height: 1.5752in


.. _RefImage_10000200000001E60000012AD278620AEBD0A733.png:

At the moment :math:`{t}_{y}`: :math:`{\sigma }_{L}({t}_{y})=\frac{-E\alpha {\sigma }_{y}^{o}}{E\alpha +{\sigma }_{y}^{o}s}`

.. code-block:: text

    The deformation energy density is: :math:`w({t}_{y})=\frac{-1}{2}E{(\alpha \theta )}^{2}`

    In the parallelepiped case we have: :math:`w({t}_{y})=\frac{-1}{2}E{(\alpha \theta )}^{2.}({x}_{B}-{x}_{A})H`

    In the axisymmetric case we have: :math:`w({t}_{y})=\frac{-1}{2}E{(\alpha \theta )}^{2}\mathrm{.}\frac{({r}_{B}^{2}-{r}_{A}^{2})}{2}H` (for 1 radian)


Elastoplastic phase
~~~~~~~~~~~~~~~~~~~~~~~~

:math:`t\ge {t}_{y}`. We are on the criteria. So:

 :math:`\dot{p}=\mathrm{max}\left[0,\frac{E-{E}_{T}}{{\mathrm{EE}}_{T}}(\dot{{\sigma }_{L}}\mathrm{sgn}({\sigma }_{L})+{\sigma }_{y}^{o}\mathrm{.}s\mathrm{.}\dot{T})\right]`

Assuming that you are "in charge" :math:`\dot{p}\ge 0`, then you eliminate :math:`\dot{p}` to have:

 :math:`\dot{{\sigma }_{L}}=-{E}_{T}\mathrm{.}\dot{T}(\alpha +\mathrm{sgn}({\sigma }_{L})\mathrm{.}\frac{E-{E}_{T}}{{\mathrm{EE}}_{T}}{\sigma }_{y}^{o}\mathrm{.}s)`

then:

 :math:`\dot{p}=\frac{E-{E}_{T}}{E}\dot{T}(-\alpha \mathrm{sgn}({\sigma }_{L})+\frac{{\sigma }_{y}^{o}\mathrm{.}s}{E})`

:math:`t={t}_{y}`, :math:`{\sigma }_{L}=-E\alpha \theta {t}_{y}<0`; we then integrate these expressions for :math:`t\ge {t}_{y}(\dot{T}=\theta )`:


*
    * :math:`\{\begin{array}{c}{\sigma }_{L}(t)=-{E}_{T}\mathrm{.}\theta (t-{t}_{y})[\alpha -\frac{E-{E}_{T}}{{\mathrm{EE}}_{T}}s{\sigma }_{y}^{o}]-{\sigma }_{L}(t-y)\\ p(t)=\frac{{\sigma }_{y}^{o}(E-{E}_{T})}{{E}^{2}}(\frac{t}{{t}_{y}}-1)\end{array}`

Or, after rearrangement, :math:`t\ge {t}_{y}`:

 :math:`\mathrm{\{}\begin{array}{c}{\sigma }_{L}(t)\mathrm{=}{\sigma }_{y}^{o}(s\theta t\mathrm{-}1+\frac{{E}_{T}}{E}(1\mathrm{-}\frac{t}{{t}_{y}}))\\ p(t)\mathrm{=}\frac{{\sigma }_{y}^{o}(E\mathrm{-}{E}_{T})}{{E}^{2}}(\frac{t}{{t}_{y}}\mathrm{-}1)\end{array}`

**Validity of this elastoplastic solution**

You have to make sure that :math:`{\sigma }_{L}(t)` stays negative. Knowing that :math:`s\theta t<1`, and :math:`t>{t}_{y}`, the previous result confirms that :math:`{\sigma }_{L}(t)<0`. Finally, we note that:

 :math:`\mathrm{sgn}({\sigma }_{L})\frac{1-2\nu }{2}\dot{p}+\dot{{\varepsilon }_{\mathrm{rr}}}=\alpha (1+\nu )\dot{T}`

from where, since :math:`{\sigma }_{L}(t)<0`:

 :math:`{\varepsilon }_{\mathrm{rr}}(t)={\varepsilon }_{\theta \theta }(t)=\alpha \theta (1+\nu )t+\frac{1-2\nu }{2}p(t),\forall t\in \left[{t}_{y},{t}_{\mathrm{fin}}\right]`


Special case of H and I models
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The thermal deformation is replaced by a given anelastic deformation. Like :math:`s=0`,

:math:`{\sigma }_{L}=-E\alpha \theta` :math:`{\varepsilon }_{\mathrm{xx}}={\varepsilon }_{\mathrm{zz}}=\alpha \theta (1+\nu )t`

The solution remains elastic as long as :math:`t<{t}_{y}=\frac{{\sigma }_{0}}{\theta E\alpha }=200s`


Benchmark results
----------------------

:math:`{\varepsilon }_{\mathrm{rr}}` Or :math:`{\varepsilon }_{\mathrm{xx}}`, :math:`{\sigma }_{\mathrm{zz}}`, and :math:`p` in :math:`{t}_{y}` and beyond:

**Elastic phase**: for :math:`t<{t}_{y}`:

:math:`{\sigma }_{L}=-E\alpha \theta t` :math:`{\varepsilon }_{\mathrm{rr}}={\varepsilon }_{\theta \theta }=\alpha \theta (1+\nu )t` in axisymmetric

     :math:`{\varepsilon }_{\mathrm{xx}}=\alpha \theta (1+\nu )t` in plane constraints

The elastic limit is reached in :math:`{t}_{y}=\frac{{\sigma }_{0}}{\theta (E\alpha +{\sigma }_{0}s)}=66.666s` from where :math:`{\sigma }_{L}({t}_{y})=\frac{\sigma }{(1+{\sigma }_{0}\frac{s}{E\alpha })}`

**Elastoplastic phase**: for :math:`t\ge {t}_{y}`

 :math:`{\sigma }_{L}(t)={\sigma }_{0}(s\theta t-1+\frac{{E}_{t}}{E}(1-\frac{t}{{t}_{y}}))`

 :math:`p(t)=\frac{{\sigma }_{0}(E-{E}_{T})}{{E}^{2}}(\frac{t}{{t}_{y}}-1)`

 :math:`{\varepsilon }_{\mathrm{rr}}={\varepsilon }_{\theta \theta }=\alpha \theta (1+\nu )t+\frac{1-2\nu }{2}p(t)` in axisymmetric

or :math:`{\varepsilon }_{\mathrm{xx}}={\varepsilon }_{\theta \theta }=\alpha \theta (1+\nu )t+\frac{1-2\nu }{2}p(t)` in plane constraints

From where:

**elastic phase**

:math:`\begin{array}{c}\{\begin{array}{c}{t}_{y}=66.666s\\ {\sigma }_{L}({t}_{y})=-133.333\mathrm{MPa}\end{array}\\ {\varepsilon }_{\mathrm{rr}}({t}_{y})={\varepsilon }_{\theta \theta }({t}_{y})=0.86666{E}^{-3}\end{array}\}\text{phase élastique}`

:math:`W=4.44410-3\mathrm{MPa}`

:math:`W=0.17778{\mathrm{MPa.mm}}^{2}` (PLAN or 3D)

:math:`W=0.26666{\mathrm{Mpa.mm}}^{3}\mathrm{.}\mathrm{rad}-1` (axi)

**Then elastoplastic phase**:

:math:`àt=80s`: :math:`\sigma (80)=-\mathrm{100.0MPa}`

     :math:`p(80)=0.300{E}^{-3}`

     :math:`{\varepsilon }_{\mathrm{rr}}(80)={\varepsilon }_{\theta \theta }(80)=1.1{E}^{-3}`

to :math:`t\mathrm{=}90s`: :math:`\sigma (90)=-\mathrm{75.0MPa}`

     :math:`p(90)=0.525{E}^{-3}`

     :math:`{\varepsilon }_{\mathrm{rr}}(90)={\varepsilon }_{\theta \theta }(90)=1.275{E}^{-3}`


Special case of H and I models
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

To :math:`t\mathrm{=}66.67s` :math:`{\sigma }_{L}(66.67)=-133.33\mathrm{MPa}`

     :math:`{\varepsilon }_{\mathrm{xx}}(66.67)={\varepsilon }_{\mathrm{zz}}(66.67)=8.667{E}^{-4}`

To :math:`t\mathrm{=}80s` :math:`{\sigma }_{L}(80)=-160\mathrm{MPa}`

     :math:`{\varepsilon }_{\mathrm{xx}}(80)={\varepsilon }_{\mathrm{zz}}(80)=8.667{E}^{-4}`

To :math:`t\mathrm{=}90s` :math:`{\sigma }_{L}(90)=-180\mathrm{MPa}`

     :math:`{\varepsilon }_{\mathrm{xx}}(90)={\varepsilon }_{\mathrm{zz}}(90)=8.667{E}^{-4}`


Uncertainty about the solution
---------------------------

Analytical solution.