Elasto-plastic numerical models

In this analysis some computational models are introduced for the finite element analysis of isotropic and kinematic hardening behaviors. The response to a cyclic load is analyzed for some stainless steels. These simulations exploit a linear quadratic finite element with elastic-plastic behavior with isotropic and kinematic hardening. These models are fundamental for analyzing the material response following a repeated cyclic load.

Background

The best known material behavior is the elastic behavior in which the observed phenomenology is that if we apply a deformation and observe a certain tension and once removed we return to the initial state.

Elasto-plasticity

An extension is given by the possibility of also taking into account a plasticity model for which, once a certain tension σ_y is exceeded, a plasticity mechanism is activated. That is, the tension will also be governed by the plastic deformation:

\sigma =\mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)

This requires, in order to be solved, the knowledge of the evolution of the plasticity mechanisms and therefore of the variation of ε_p. That is, it requires the introduction of an evolutionary law that allows to express the phenomena of plasticity and at the same time respect the thermodynamic laws. to do this, a yield function is defined, expressed in the space of tensions, inside which we will be in a regime of elasticity while outside plastic mechanisms will be triggered.

A choice of strain energy that satisfies both the material behavior and the laws of thermodynamics will be:

\psi\left(\underline{\varepsilon}, \underline{\varepsilon}_{p}\right)=\frac{1}{2}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)^{t} \mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)

Then the importance lies in the introduction of an evolutionary law that satisfies the thermodynamic requirements and the typical choice consists in an associative law. That is, a law such as to be associated with the yield function through an unknown multiplier, where f is chosen as a convex function.

 \dot{\varepsilon}_{p}=\dot{\lambda} \frac{\partial f}{\partial \sigma}

This explain as the plastic deformation is a vector with a direction normal to the yield surface. Following the evolution of the tension, once the yield value has been reached, this must remain at this tension level, so once the plasticity is activated it will remain on the surface as it moves I will move on the edge.

Then this, from the computational point of view, introduces a new function (beyond the elastic equilibrium) which will have to be solved by the same stress state through an iterative procedure.

Hardening

An even more sophisticated model includes the possibility of hardening, i.e. the plastic deformation will no longer be linear but will have a certain slope, expressed precisely by the hardening module.

Therefore a first model of hardening is the isotropic one where, once the load is removed, one goes back along the same elastic branch. Some plastic deformation will therefore remain. Continuing the loading process in the opposite direction will activate the opposite plasticity. The work hardening will have determined a greater resistance of the material in all directions.

On the other hand, in kinematic hardening, by varying the direction, the new level of yield stress at which the plasticity is activated can vary.

In the ideal case, the stress state is measured by verifying that it is within the yield function

f(\underline{\sigma})=\|\underline{\sigma}\|-\sigma_{y_{0}}

While in isotropic hardening, the yield function will vary:

f(\underline{\sigma}, q)=\|\underline{\sigma}\|-\left(\sigma_{y_{0}}-q\right)

It is a free energy function of a scalar field α which takes into account the Hardening module, H, that is:

\psi\left(\underline{\varepsilon}, \underline{\varepsilon}_{p}, \alpha\right)=\frac{1}{2}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)^{T} \mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)+\frac{1}{2} H \alpha^{2}

Instead in the kinematic hardening instead it will simply be translated:

f(\underline{\sigma}, \underline{q})=\|\underline{\sigma}-\underline{q}\|-\sigma_{y_{o}}

It is a free energy similar to the previous one but which takes into account a tensor quantity:

\psi\left(\underline{\varepsilon}, \underline{\varepsilon}_{p}, \underline{\alpha}\right)=\frac{1}{2}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)^{T} \mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)+\frac{1}{2} \underline{\alpha}^T \mathbb{H} \underline{\alpha}

This leads to consider the strain of hardening that takes the name of backstress to which the associative law is referred:

\underline{q}=-\frac{\partial \psi}{\partial \underline{\alpha}}=-\mathbb H\underline{\alpha}
\\\:\\
\dot{\alpha}=\dot{\lambda} \frac{\partial f}{\partial q}

FEM formulations

This hardening model must merge with the finite element formulation to solve the problem of elastic equilibrium. Starting from the weak form we go back to the relative geometry and introduce the interpolation on the nodal coordinates.

\int_{V} \underline{\sigma} \cdot \delta \underline{\varepsilon} d V-\int_{V} \underline{f} \cdot \delta \underline{u} d V-\int_{A} t \cdot \delta u d A=0

That is, it goes back to the i-th component of the residue

[R]_i=\sum_{g} B_{I}^{T}(\xi, \eta) \mathbb{C}\left(\sum_{k} B_{k}\left(\xi_{g}, \eta_{g}\right) \underline{u}_{k}-\left.\underline{\varepsilon}\right|_{g}\right) J_{e, g} d V

This requires the nodal coordinates and the plastic deformation regulated by the flow rule, variable over time, for which a discretization in time steps is used.

\underline{\varepsilon}_{p}=\underline{\varepsilon}_{p_{n}}+\Delta t \cdot \underline{\dot{\varepsilon}}_{p_{n}}

That is:

\dot{\varepsilon}_{p}=\frac{\dot{\varepsilon}_{p}-\underline{\varepsilon}_{p_{n}}}{\Delta t}=\dot{\lambda} \cdot \underline{n}_{\sigma}=\frac{\lambda-\lambda_{n}}{\Delta t} \underline{n}_{\sigma}

And from this relation it is possible to consider the evolutionary laws representing an independent time plasticity that is not linked to the speed with which the load is applied:

\left\{\begin{array}{l}\underline{\varepsilon}_{p}=\underline{\varepsilon}_{p_{n}}+\left(\lambda-\lambda_{n}\right) \underline{n}_{\sigma}(\underline{\sigma}, q) \\ \underline{\alpha}=\underline{\alpha}+\left(\lambda-\lambda_{n}\right) \underline{n}_{q}(\underline{\sigma}, q) \\ f(\underline{\sigma} \cdot \underline{q})=0\end{array}\right.

Therefore these relations constitute a second algebraic system that will have to be solved together with the one representing the elastic equilibrium.

R(\underline{u}, \underline{h})=0
\\\:\\
\underline{Q}(\underline{u}, \underline{h})=0

In solving the equilibrium a predictor-corrector algorithm is exploited which calculates the stresses by initially assuming an elastic deformation and verifying that the stress state is within the limits established by the yield function. If the condition is not satisfied, the prediction is corrected by finding the solution from the system Q which reports on the edge of the yield function. Then a new stress state is found which will also have to satisfy the equilibrium on the R system. We proceed with iterative steps satisfying both algebraic systems.

Simulations

The simulation testing is conducted in AceFEM. It focuses on two different structures. First of all, a patch test is considered, to verify the correctness of the code and to have a confirmation of how the variation of the hardening characteristics influences the structural response and the trend of the internal stresses. The analysis is then moved to a plate (L x H) with a slot of variable size. The analysis starts using a finite element Q1, quadrilateral with linear shape functions, with plastic elastic behavior and isotropic hardening. This will then be transformed into an element capable of reproducing an elasto-platic behavior with kinematic hardening

Q1EPS

To implement a kinematic hardening behavior it is necessary to consider that the variable α will become a vector variable and the evolutionary law will vary taking into account the translation of the yield function in the stress space.

\underline q=-{d\psi \over d\underline \alpha}
\\\:\\
f(\underline \sigma,\underline q)=\sqrt{{3\over 2}\|\underline \sigma_{\operatorname{dev}}-\underline q}-\sigma_{yo}

Then the system for the corrector predictor algorithm will also be calculated differently:

\bold Q_{\alpha}=\underline \alpha-\underline \alpha_n-(\lambda-\lambda_n){\partial f\over \partial \underline q}

This is reflected in the introduction of new labor variables, updating the definition of the deformation energy by taking into account the internal product and in redefining the evolutionary law.

SMSFreeze[\[Alpha], {{\[DoubleStruckH]g[[4]], \[DoubleStruckH]g[[
    6]]}, {\[DoubleStruckH]g[[6]], \[DoubleStruckH]g[[5]]}}, 
 "Symmetric" -> True]
 \[Lambda] = \[DoubleStruckH]g[[7]];

\[CapitalPsi]e \[DoubleRightTee] \[Lambda]e/ 2 (Tr[\[DoubleStruckCapitalD]e])^2 + \[Mu]e Tr[\\[DoubleStruckCapitalD]e . \[DoubleStruckCapitalD]e] +  1/2 H Tr[Transpose[\[Alpha]] . \[Alpha]];

fg \[DoubleRightTee] SMSSqrt[3/2 Tr[Transpose[\[Sigma]dev - q] . (\[Sigma]dev - q)]] - (\[Sigma]YO);

Furthermore, the system for the corrector predictor algorithm is also updated taking into account the new values:

\[Alpha]n = {{\[DoubleStruckH]gnIO[[4]], \[DoubleStruckH]gnIO[[
   5]]}, {\[DoubleStruckH]gnIO[[5]], \[DoubleStruckH]gnIO[[6]]}};
\[Lambda]n = \[DoubleStruckH]gnIO[[7]];

SMSFreeze[q, -SMSD[\[CapitalPsi]e, \[Alpha], "Symmetric" -> True], 
  "Symmetric" -> True];

\[DoubleStruckCapitalQ] = {\[DoubleStruckCapitalQ]\[Epsilon][[1, 
    1]], \[DoubleStruckCapitalQ]\[Epsilon][[2, 
    2]], \[DoubleStruckCapitalQ]\[Epsilon][[1, 2]], Qq[[1, 1]], 
   Qq[[2, 2]], Qq[[1, 2]], Q\[Lambda]};

This allow to generate a new finite element like:

• Q1EPS, quadrilateral of order 1 with elastic-plastic behavior and isotropic hardening
• Q1EPS2, quadrilateral of order 1 with elastic-plastic behavior and kinematic hardening

Hardening features

It is possible to consider a stainless steel given the excellent mechanical properties and the wide field of application in medical devices, surgical instruments, prostheses and much more. It is precisely their chemical composition with chromium present at at least 10.5% and a reduced carbon content below 1.2% that makes them highly resistant to corrosion and leads to the generation of microstructures that are typically not present in steels and led to these mechanical properties. Clearly, as the composition of the steel changes, different mechanical properties will vary, including the structural response. We can therefore consider some average parameters and undertake a parametric study. It is useful to remember that we can define the hardening modulus starting from the yield strength and the critical load as:

H={\sigma_u-\sigma_y\over {\varepsilon_{u,\%}\over 100} - \varepsilon _y}

Stainless steels are typically more ductile and have the ability to accommodate very large plastic deformations. These plastic properties are reflected in hardening and an increase in yield stress, ultimate stress and hardness. A variable hardening modulus between 600 and 1600 MPa and a yield stress between 175 and 450 MPa is then chosen.

Structural response and isotropic hardening

Patch test

The patch test is a test that is carried out on a very simple domain that is discretized in an irregular way but with a very loose mesh. This allows, by exploiting boundary conditions of which the analytical solution is known, to test new algorithms. Furthermore, it also allows us to perceive how the variations of some material parameters influence the structural behavior by acting on material responses that are simpler in themselves.

In this case a square domain discretized in 5 elements is analyzed. A uniaxial traction is required for which we will obtain a constant tension. We will then apply a time-varying load condition that alternates compression and traction through the multiplier λ.

It is then possible to analyze the trend of stresses over time by observing exactly a behavior with isotropic hardening.

This behavior shows us how the stress σyy grows linearly both with the increase of σy and of the hardening modulus.

Plate

Considering instead a plate L x H with a half-height slot of length equal to L / 4 it is possible to observe how its structural response varies. By applying a vertical tensile load, a crack widening is observed with a concentration of stresses right next to the crack apex. That is, there will be a separation.

Therefore it is possible to analyze the deformation as the hardening characteristics variability.

Yielding stress

A reduction in crack enlargement is observed as the yield stress increases. A material with a higher yield stress will exhibit less displacement. The results in the figure consider an average Hardening coefficient.

Hardening coefficient

As the hardening coefficient increases, a reduction in vertical displacement is observed. This is due to the fact that considering the same load conditions there will be a stiffer material that will respond with less plastic deformation.

The overall effect of both a variation in the yield stress and the hardening coefficient can be observed with a contour plot. The following figure shows how, as the yield stress or the hardening modulus increases, there is an increasingly smaller residual plastic deformation. In particular, the figure plots the vertical displacement of the upper edge of the crack as the hardening parameters vary.

Variations

Considering the case with the hardening parameters in average conditions, it is possible to see how the load multiplier influences the displacement as time varies.

Another important variation is given by the tensions σyy by reading them in the middle of the plate vertically with a coordinate that varies from the end of the crack to the right edge.

Effort intensification factor

Considering the hardening parameters set at an average value, it is possible to observe what happens when the crack length changes. A very interesting parameter to observe is the effort intensification factor. With the presence of the crack, a distribution of strees inside the plate is created that is different from that which would be considered a homogeneous plate. Stress is observed to increase as we approach the crack.

Different cases are analyzed in which the crack has length:

• L/8
• L/4
• L/2
• 3L/4

Kinematic hardening

The results obtained so far consider the finished element Q1EPS, that is an elastic plastic behavior with isotropic hardening. Using the Q1EPS2 element instead, it is possible to analyze a kinematic hardening. The behavior is immediately evident by carrying out the simulation on a patch test.

This behavior more closely mirrors what actually happens in metals where the tensile yield strength increases at the expense of the compressive yield strength. This is known as the Bauschinger effect. That is, the compressive yield stress of a material brought to yield by traction is lower than that of the original material. The same effect is achieved with a reverse load history. This effect is witnessed by the asymmetry of the hysteresis cycle of a ductile material.

Macroscopic response

Evidence of a kinematic hardening also appears by carrying out the simulation with a time-varying load according to a triangular law on the plate with the crack or on a Cook membrane.

The variation of the compression yield stress is reflected in a greater yield effect in the compression phase.

In isotropic hardening, the material remembers the previous loading phase and the subsequent yield is represented by a point with a higher yield stress but increased homothetically. On the other hand, in isotropic hardening, there is a translation of the boundary surface. This backstress is reflected in a reduction of the yield stress in a different direction than that of the first load cycle.

You can also export the animation to see how the displacement and deformation of the Cook membrane varies along the load history (time steps).

Code

The code is available at the GitHub repository: https://github.com/mastroalex/elasto-plastic

In particular:

• Variations of the hardening properties ../final_test_ripetuti_11_passi/
• Variations of the crack dimensions ../final_test_vari/
• Kinematic hardening ../final_kinematic/

References

• von Mises Material Properties – Autodesk Knowledge
• Mechanical properties of structural stainless steels – Tylek, I and Kuchta, K
• Materiali per l’ingegneria civile ed industriale – Callisters, W and Rethwisc, D
• Scienza e tecnologia dei materiali – Smith, W and Hashemi, J
• Stainless Steels for Medical and Surgical Applications – Nutt, M