Reference solutions
======================

Calculation method used for reference solutions
----------------------------------------------------------

Linear workpieces
~~~~~~~~~~~~~~~~~~~~~~~~

**Isotropic work hardening**

For uniaxial traction, the plasticity criterion is written as:

:math:`∣{\sigma }_{L}∣-R(p)\le 0`

where :math:`p` is the cumulative plastic deformation

:math:`R(p)=R\text{'}p+{\sigma }^{y}` and :math:`R\text{'}=\frac{E{E}_{1}}{E-{E}_{1}}`

The criterion is then written:

:math:`∣{\sigma }_{L}∣-R\text{'}p-{\sigma }^{y}\le 0`

The stress tensor is obtained by:

:math:`\sigma =A\mathrm{.}(\varepsilon (u)-{\varepsilon }^{p})-\mathrm{3K}\alpha (T-{T}^{\mathrm{ref}})\mathrm{Id}`

So we deduce the expression for :math:`{\sigma }_{L}`

:math:`{\sigma }_{L}=E(\varepsilon -\alpha T)-E{\varepsilon }^{p}` :math:`({T}^{\mathrm{ref}}=0)`

In our case, :math:`\varepsilon =0` so:

:math:`{\sigma }_{L}=E{\varepsilon }_{L}-E{\varepsilon }^{p}` with :math:`{\varepsilon }_{L}=-\alpha T`

So:

• If :math:`∣{\sigma }_{L}∣-R(p)<0`:

*p=0* and :math:`{\sigma }_{L}=E{\varepsilon }_{L}`

• If :math:`∣{\sigma }_{L}∣-R(p)=0`:

:math:`p=(\frac{∣\sigma ∣-{\sigma }^{y}}{R\text{'}})`

:math:`{\sigma }_{L}=E{\varepsilon }_{L}-E{\varepsilon }^{p}`

**Application to the loading path**

**Instant 1:**

:math:`\sigma =E\varepsilon =\mathrm{200MPa}` and :math:`R(p)=R\text{'}p+{\sigma }^{y}=\mathrm{100MPa}` because*p=0*.

We have :math:`{\sigma }_{L}-R(p)\le 0`.

The criterion has not been met, the evolution is elastic: :math:`{\sigma }_{L}=\mathrm{100MPa}` and :math:`N=\mathrm{100kN}`

**Instant 2:**

The criterion is met:

:math:`\begin{array}{}{\sigma }_{L}=\frac{E}{E+R\text{'}}(R\text{'}{\varepsilon }_{L}+{\sigma }^{y})=\frac{2\cdot {10}^{11}}{2\cdot {10}^{11}+2.02\cdot \mathrm{10⁹}}(2.02\cdot \mathrm{10⁹}\times 3.5\cdot {10}^{-3}+2\cdot \mathrm{10⁸})=\mathrm{205MPa}\\ N=\mathrm{102.5kN}\end{array}`

and :math:`p=2.475\cdot {10}^{-3}`

**Instant 3:**

We discharge elastically:

:math:`\begin{array}{}{\sigma }_{L}=E{\varepsilon }_{L}-E{\varepsilon }^{p}=2\cdot {10}^{11}(1.5\cdot {10}^{-3}-2.475\cdot {10}^{-3})=-195\mathrm{MPa}\\ N=-97.5\mathrm{kN}\end{array}`

**Instant 4:**

We laminate again:

Criteria is written as: :math:`\mid \sigma \mid -R\text{'}p-{\sigma }^{y}=0` with :math:`p={p}_{1}+{p}_{2}` where :math:`{p}_{1}=2.475\cdot {10}^{-3}`

So we get:

:math:`{p}_{2}=\frac{\mid \sigma \mid -{\sigma }^{y}}{R\text{'}}-{p}_{1}`

:math:`\begin{array}{}\sigma =-E{\varepsilon }^{p}=-E({p}_{1}-{p}_{2})\\ \sigma =(\frac{R\text{'}}{R\text{'}+E})(-2E{p}_{1}-\frac{E{\sigma }^{y}}{R\text{'}})=-\mathrm{207.9MPa}\end{array}`

And so :math:`N=-\mathrm{103.95kN}`

**Instant 5:**

We discharge elastically:

:math:`\begin{array}{}{\sigma }_{L}=E{\varepsilon }_{L}-E{\varepsilon }^{p}=2\cdot {10}^{11}(2\cdot {10}^{-3}-1.0395\cdot {10}^{-3})=192.1\mathrm{MPa}\\ N=96.05\mathrm{kN}\end{array}`

**Moments 6 and 7:**

The reasoning is the same.

We find:

:math:`\begin{array}{}{N}_{(\mathrm{inst.6})}=\mathrm{105.87kN}\\ {N}_{(\mathrm{inst.7})}=-\mathrm{44.13kN}\end{array}`

**Kinematic work hardening**

The calculation method is the same, but in this case, the plasticity criterion is written as:

:math:`\sigma -{X}_{\mathrm{eq}}-{\sigma }^{y}\le 0` with :math:`{X}_{\mathrm{eq}}=C{({\varepsilon }^{p})}_{\mathrm{eq}}=\frac{3}{2}C{\varepsilon }^{p}=\frac{E{E}_{t}}{E-{E}_{t}}{\varepsilon }^{p}`

With the previous notations, the criterion is written:

:math:`∣{\sigma }_{L}-R\text{'}{\varepsilon }^{p}∣-{\sigma }^{y}\le 0` from where :math:`{\sigma }_{L}=R\text{'}{\varepsilon }^{p}\pm {\sigma }^{y}` (depending on the direction of flow).

**Application to the loading path**

**Instant 1:**

The criterion has not been met, the evolution is elastic: :math:`{\sigma }_{L}=\mathrm{100MPa}` and :math:`N=\mathrm{100kN}`

**Instant 2:**

The criterion is met: :math:`\mid {\sigma }_{L}-R\text{'}{\varepsilon }^{p}\mid -{\sigma }^{y}=0`

:math:`{\sigma }_{L}=R\text{'}{\varepsilon }^{p}+{\sigma }^{y}=2.02\cdot \mathrm{10⁹}\times 2.475\cdot {10}^{-3}+2\cdot \mathrm{10⁸}=\mathrm{205MPa}`

**Instant 3:**

We discharge elastically:

:math:`\begin{array}{}{\sigma }_{L}=E{\varepsilon }_{L}-E{\varepsilon }^{p}=2\cdot {10}^{11}(1.5\cdot {10}^{-3}-2.475\cdot {10}^{-3})=-\mathrm{195MPa}\\ N=-\mathrm{97.5kN}\end{array}`

**Instant 4:**

:math:`\mid \sigma -R\text{'}{\varepsilon }^{p}\mid -{\sigma }^{y}=0` with :math:`{p}_{1}=2.475\cdot {10}^{-3}`

:math:`\begin{array}{}{\varepsilon }^{p}={p}_{1}-{p}_{2}\\ {p}_{2}={p}_{1}-\frac{\mid \sigma +{\sigma }^{y}\mid }{R\text{'}}\\ \sigma =-E{\varepsilon }^{p}=-E({p}_{1}-{p}_{2})\\ \sigma =-E(\frac{\mid \sigma +{\sigma }^{y}\mid }{R\text{'}})=-\mathrm{198MPa}\\ N=-\mathrm{99kN}\end{array}`

**Instant 5:**

We discharge elastically:

:math:`\begin{array}{}{\sigma }_{L}=E{\varepsilon }_{L}-E{\varepsilon }^{p}=2\cdot {10}^{11}(2\cdot {10}^{-3}-9.9\cdot {10}^{-4})=202\mathrm{MPa}\\ N=\mathrm{101kN}\end{array}`

**Moments 6 and 7:**

The reasoning is the same. We find:

:math:`\begin{array}{}{N}_{(\mathrm{inst.6})}=\mathrm{103kN}\\ {N}_{(\mathrm{inst.7})}=-\mathrm{47kN}\end{array}`


Pinto-Menegotto model
~~~~~~~~~~~~~~~~~~~~~~~~~~~

This model is described in the *Code_Aster* Reference Manual [:ref:`R5.03.09 <R5.03.09>`] [:ref:`bib1 <bib1>`]. The constitutive law of steels is composed of two distinct parts: monotonic loading composed of three successive zones (linear elasticity, plastic level and work hardening) and cyclic loading where the path between two points of reversal (half-cycle) is described by an analytical expression curve of the type :math:`\sigma =f(\varepsilon )`.

As before, the deformations imposed are thermal deformations: :math:`\varepsilon =-\alpha T`


Case without buckling
^^^^^^^^^^^^^^^^^^

**First load**

• Linear elasticity: :math:`\sigma =E\varepsilon`

**Instant 1:**

:math:`N=E\varepsilon S=2\cdot {10}^{11}\times 1\cdot {10}^{-3}\times 5\cdot {10}^{-4}=\mathrm{100kN}`

• Plastic bearing: :math:`\sigma ={\sigma }^{y}`

• Degree 4 polynomial: :math:`\sigma ={\sigma }_{u}-({\sigma }_{\mathrm{su}}-{\sigma }_{y}^{0}){(\frac{{\varepsilon }_{u}-\varepsilon }{{\varepsilon }_{u}-{\varepsilon }_{h}})}^{4}`

**Instant 2:**

:math:`\varepsilon =3.5\cdot {10}^{-3}>{\varepsilon }_{h}=2.3\cdot {10}^{-3}`, we use the degree 4 polynomial:

:math:`\sigma =\mathrm{209.416MPa}` and :math:`N=\mathrm{104.708kN}`

**Cycles**

**Half cycle 1:**

We determine :math:`{\zeta }_{p}^{0}`:

:math:`{\zeta }_{p}^{0}={\varepsilon }_{r}^{0}-{\varepsilon }_{y}^{0}=3.5\cdot {10}^{-3}-{1.10}^{-3}=2.5\cdot {10}^{-3}` because :math:`{\varepsilon }_{r}^{0}={\varepsilon }_{(\mathrm{inst.2})}`

Then :math:`\Delta {\sigma }^{0}`:

:math:`\Delta {\sigma }^{0}={E}_{h}{\zeta }_{p}^{0}=2\cdot {10}^{9}\times 2.5\cdot {10}^{-3}=\mathrm{5MPa}`

Where :math:`{\sigma }_{y}^{1}={\sigma }_{y}^{0}\cdot \mathrm{sign}(-{\zeta }_{p}^{0})+\Delta {\sigma }^{0}=-200+5=-\mathrm{195MPa}`

We then calculate :math:`{\varepsilon }_{y}^{1}`:

:math:`{\varepsilon }_{y}^{1}={\varepsilon }_{r}^{0}+\frac{{\sigma }_{y}^{1}-{\sigma }_{r}^{0}}{E}=3.5\cdot {10}^{-3}+\frac{(-195-209.416)\mathrm{10⁶}}{2.0\cdot {10}^{11}}=1.477\cdot {10}^{-3}`

We thus determine :math:`\sigma *=f(\varepsilon*)`, defined by:

:math:`{\sigma }^{}=b{\varepsilon }^{}+(\frac{1-b}{{(1+{({\varepsilon }^{})}^{R})}^{(1/R)}}){\varepsilon }^{}`, with :math:`b=\frac{{E}_{h}}{E}`

:math:`\begin{array}{}{\varepsilon }^{}=\frac{\varepsilon -{\varepsilon }_{r}^{0}}{{\varepsilon }_{y}^{1}-{\varepsilon }_{r}^{0}}\\ \\ {\sigma }^{}=\frac{\sigma -{\sigma }_{r}^{0}}{{\sigma }_{y}^{1}-{\sigma }_{r}^{0}}\end{array}`

:math:`{\xi }_{p}^{0}=\frac{{\zeta }_{p}^{0}}{\mid {\varepsilon }_{y}^{1}-{\varepsilon }_{r}^{0}\mid }` and :math:`\mathrm{R¹}={R}_{0}-\frac{{A}_{1}\cdot {\xi }_{p}^{0}}{{A}_{2}+{\xi }_{p}^{0}}`

We get :math:`{\xi }_{p}^{0}=-1.23` and :math:`\mathrm{R¹}=3.51`

We can then calculate the value of :math:`\sigma` at times 3 and 4:

**Instant 3:**

:math:`{\varepsilon }^{}=\frac{{\varepsilon }_{(\mathrm{inst.3})}-{\varepsilon }_{r}^{0}}{{\varepsilon }_{y}^{1}-{\varepsilon }_{r}^{0}}=\frac{1.5\cdot {10}^{-3}-3.5\cdot {10}^{-3}}{1.477\cdot {10}^{-3}-3.5\cdot {10}^{-3}}=0.988`

:math:`{\sigma }^{}=b{\varepsilon }^{}+(\frac{1-b}{{(1+{({\varepsilon }^{})}^{R})}^{(1/R)}}){\varepsilon }^{}=0.01\times 0.988+(\frac{1-0.01}{{(1+{(0.988)}^{3.51})}^{(1/3.51)}})=0.82`

and :math:`\sigma ={\sigma }^{}({\sigma }_{y}^{1}-{\sigma }_{r}^{0})+{\sigma }_{r}^{0}=0.82\times (-195-209.416)+209.416=-122\mathrm{MPa}`

Hence :math:`N=-\mathrm{61kN}`

**Instant 4:**

We use the same method, with :math:`\varepsilon =0`.

:math:`\begin{array}{}{\varepsilon }^{}=1.73\\ {\sigma }^{}=0.56\\ \sigma =-20\mathrm{MPa}\\ N=-\mathrm{10kN}\end{array}`

**Half cycle 2: Instant 5 and 6:**

The calculation method is identical, we determine:

:math:`{\zeta }_{p}^{1}`, :math:`{\sigma }_{y}^{2}`, :math:`{\varepsilon }_{y}^{2}`, :math:`{\xi }_{p}^{1}`, :math:`\mathrm{R²}`, then :math:`{\sigma }_{(\mathrm{inst.5})}^{}=f({\varepsilon }_{(\mathrm{inst.5})}^{})` and :math:`{\sigma }_{(\mathrm{inst.6})}^{}=f({\varepsilon }_{(\mathrm{inst.6})}^{})`

and finally :math:`{\sigma }_{(\mathrm{inst.5})}` and :math:`{\sigma }_{(\mathrm{inst.6})}`.

**Half-cycle 3: Instant 7**: Same


Case with buckling
^^^^^^^^^^^^^^^^^^

**First load**

Same as the previous case.

**Cycles**

**Half cycle 1 (compression) :**

The calculation method is the same, but the value of the asymptote slope is changed:

A new coefficient :math:`{b}_{c}` is calculated:

:math:`{b}_{c}=a(5.0-L/D)e(b\xi \text{'}\frac{E}{{\sigma }^{y}-{\sigma }^{\infty }})=0.006\times (5.0-5.9){e}^{(0.01\times 1.477\cdot {10}^{-3}\frac{2\cdot {10}^{11}}{2\cdot {10}^{8}-1.36\cdot {10}^{8}})}=-0.0057`

Next, as in the model without buckling, we must determine :math:`{\sigma }_{y}^{n}`. The reasoning is the same, but an additional constraint :math:`{\sigma }_{s}^{}` is added in order to correctly position the curve in relation to the asymptote.

:math:`{\sigma }_{s}^{}={\gamma }_{s}bE\frac{b-{b}_{c}}{1-{b}_{c}}=0.028\times 0.01\times 2\cdot {10}^{11}\times \frac{0.01+0.0057}{1+0.0057}=0.87\mathrm{MPa}`

where :math:`{\gamma }_{s}` is given by: :math:`{\gamma }_{s}=\frac{11.0-L/D}{10({e}^{c(L/D)}-1.0)}=0.028`

**Half cycle 2 (traction) :**

• In traction, a reduced Young's modulus is adopted:

:math:`{E}_{r}=E({a}_{5}+(1.0-{a}_{5}){e}^{-{a}_{6}\xi {"}^{2}})=2\cdot {10}^{11}\times (\mathrm{0,88}+(1-0.88){e}^{(-620\times 1.473\cdot {10}^{-6})})=1.99\cdot {10}^{11}\mathrm{MPa}`

with :math:`{a}_{5}=1.0+(5.0-L/D)/7.5=0.88`

The rest of the method is the same.


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

Constant :math:`N` normal effort on the bar


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

None, the solution is analytical


Bibliographical references
---------------------------

[:ref:`1 <1>`] *Code_Aster* reference manual [:ref:`R5.03.09 <R5.03.09>`].

[:ref:`2 <2>`] S. ANDRIEUX: TD 1 Three perfect Von Mises thermoelastoplastic bars. In "Introduction to thermoplasticity in *Code_Aster*", HI-74/96/013 November 1996 (course reference manual).