This article shows how to implement an orthotropic plastic behaviour with isotropic linear hardening in MFront. Such an example illustrates:
StandardElasticity brick (see this page).@ComputeStiffnessTensor keyword to define the orthotropic stiffness tensor.@HillTensor keyword to define the Hill tensor.The whole implementation is available here.
The behaviour is described by a standard decomposition of the strain \(\underline{\varepsilon}^{\mathrm{to}}\) in an elastic and a plastic component, respectively denoted \(\underline{\varepsilon}^{\mathrm{el}}\) and \(\underline{\varepsilon}^{\mathrm{p}}\):
\[ \underline{\varepsilon}^{\mathrm{to}}=\underline{\varepsilon}^{\mathrm{el}}+\underline{\varepsilon}^{\mathrm{p}} \]
The stress \(\underline{\sigma}\) is related to the the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) by a the orthotropic elastic stiffness \(\underline{\underline{\mathbf{D}}}\):
\[ \underline{\sigma}= \underline{\underline{\mathbf{D}}}\,\colon\,\underline{\varepsilon}^{\mathrm{el}} \]
The plastic part of the behaviour is described by the following yield surface: \[ f{\left(\sigma_{H},p\right)} = \sigma_{H}-\sigma_{0}-R\,p \]
where \(\sigma_{H}\) is the Hill stress defined below, \(p\) is the equivalent plastic strain, \(\sigma_{0}\) is the yield stress and \(R\) is the hardening slope.
The Hill stress \(\sigma_{H}\) is defined using the fourth order Hill tensor \(H\): \[ \sigma_{H}=\sqrt{\underline{\sigma}\,\colon\,\underline{\underline{\mathbf{H}}}\colon\,\underline{\sigma}} \]
The plastic flow is assumed to be associated, so the flow direction \(\underline{n}\) is given by \({\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{\sigma}}}\):
\[ \underline{n} = {\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{\sigma}}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle \sigma_{H}}}}\,\underline{\underline{\mathbf{H}}}\,\colon\,\underline{\sigma} \]
The previous constitutive equations will be integrated using a standard implicit scheme.
Assuming a plastic loading, the system of equations to be solved is: \[ \left\{ \begin{aligned} \Delta\,\underline{\varepsilon}^{\mathrm{el}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}+\Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} &= 0 \\ f{\left({\left.\sigma_{H}\right|_{t+\theta\,\Delta\,t}},{\left.p\right|_{t+\theta\,\Delta\,t}}\right)} &= 0 \\ \end{aligned} \right. \]
where \({\left.X\right|_{t+\theta\,\Delta\,t}}\) is the value of \(X\) at \(t+\theta\,\Delta\,t\), \(\theta\) being a numerical parameter. In the following, the first (tensorial) equation is noted \(f_{\underline{\varepsilon}^{\mathrm{el}}}\) and the second (scalar) equation is noted \(f_{p}\).
In practice, it is physically sound to make satisfy exactly the yield condition at the end of the time step (otherwise, stress extrapolation can lead to stress state outside the yield surface and spurious oscillations can also be observed). This leads to the choice \(\theta=1\).
The jacobian \(J\) of the implicit system can be decomposed by blocks: \[ J= \begin{pmatrix} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p}} & \\\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,p}} \\ \end{pmatrix} \]
The expression of the previous terms is given by:
\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= \underline{I} + {{\displaystyle \frac{\displaystyle \theta\,dp}{\displaystyle \sigma_{H}}}}\,{\left(\underline{\underline{\mathbf{H}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}\,\underline{\underline{\mathbf{D}}} \\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p}} &= {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= -\theta\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\colon\,\underline{\underline{\mathbf{D}}}\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,p}} &= -R\,\theta \end{aligned} \right. \]
Assuming an elastic loading, the system of equations to be solved is trivially: \[ \left\{ \begin{aligned} \Delta\,\underline{\varepsilon}^{\mathrm{el}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}&= 0 \\ \Delta\,p &= 0 \\ \end{aligned} \right. \]
The jacobian associated with this system is the identity matrix.
In the plastic loading case, the system of equations to be solved is non-linear. We choose the Newton-Raphson algorithm with an analytical jacobian. Compared to other algorithm available in MFront (Runge-Kutta, Broyden, Newton-Raphson with numerical jacobian, etc..), this algorithm is generally the most efficient in pratice.
The implementation begins with the choice of the Implicit domain specific language (dsl):
@DSL Implicit;Note that this dsl automatically declares the elastic strain eel as a state variable.
As discussed before, we explicit state that a fully implicit integration will be used by default:
@Theta 1;This value can be changed at runtime by modifying the theta parameter.
The stopping criterion is chosen low, to ensure the quality of the consistent tangent operator (the default value, \(10^{-8}\) is enough to ensure a good estimation of the state variable evolution, but is not enough to have a proper estimation of the consistent tangent operator):
@Epsilon 1e-14;We then declare the behaviour to be orthotropic. We choose the pipe orthotropic convention to have a consistent definition of the elastic stiffness tensor and hill tensor for all modelling hypotheses (axes are automatically permuted for plane modelling hypotheses):
@OrthotropicBehaviour<Pipe>;We then declare that we want to support all the modelling hypotheses:
@ModellingHypotheses {".+"};The support of plane stress modelling hypotheses are handled by the StandardElasticity brick and will not be discussed here.
StandardElasticity brickTo implement this behaviour, we will use the StandardElasticity brick which provides:
This behaviour brick is fully described here.
The usage of the StandardElasticity is declared as follows:
@Brick StandardElasticity;The elastic stiffness tensor \(D\) is defined using @ComputeStiffnessTensor keyword by giving the elastic material properties as constants:
@ComputeStiffnessTensor<UnAltered> {
// YoungModulus1 YoungModulus2 YoungModulus3
7.8e+10,2.64233e+11,3.32e+11,
// PoissonRatio12 PoissonRatio23 PoissonRatio13
0.13,0.24,0.18,
// ShearModulus12 ShearModulus23 ShearModulus13
4.8e+10,1.16418e+11,7.8e+10
};This computed stiffness is stored in a variable D. A second variable D_tdt is also introduced. As the material properties are constants, D_tdt is an alias to D.
The elastic material properties can be changed at runtime time by modifying the following parameters: YoungModulus1, YoungModulus2,YoungModulus3, PoissonRatio12, PoissonRatio23, PoissonRatio13, ShearModulus12, ShearModulus23 and ShearModulus13.
Rather than constants, one can also use correlations implemented in seperate MFront files. This allows to take into account the dependency of the material properties with the temperature for example. In this case, the variable D contains the stiffness tensor at \(t+\theta\,\Delta\,t\) and the variable D_tdt contains the stiffness tensor at \(t+\Delta\,t\).
Another possibility is to use the @RequireStiffnessTensor keyword. In this case, the elastic material properties must be computed by the calling solver at the end of the time step (and furnished to the mechanical behaviours through hidden material properties).
For the computation of the Hill tensor, we make use of the @HillTensor keyword:
@HillTensor H {0.371,0.629,4.052,1.5,1.5,1.5};As stated earlier, the state variable eel is automatically declared by the Implicit dsl.
The equivalent plastic strain state variable p is declared as:
@StateVariable real p;We then associate the appropriate glossary name to this variable:
p.setGlossaryName("EquivalentPlasticStrain");The definition of yield surface introduces two material coefficients \(\sigma_{0}\) and \(R\) that we declare as parameters:
@Parameter s0 = 150e6;
s0.setGlossaryName("YieldStress");
@Parameter R = 150e9;
R.setEntryName("HardeningSlope");The YieldStress is an entry of the glossary (see here). The HardeningSlope name is not declared in the glossary name (yet) and is then associated to the \(R\) variable with the setEntryName method.
To select the implicit system associated either with the elastic or plastic loading case, we introduce a boolean local variable b.
@LocalVariable bool b; // if true, plastic loadingBefore solving the implicit system, the code block introduced by the @InitLocalVariables keyword is executed. For this behaviour, this block will select either the elastic or plastic loading case.
We first make an elastic prediction of the stress at the end of the time step. We use the computeElasticPrediction method introduced by the StandardElasticity brick. This method takes into account the modelling hypothesis, which is mandatory for plane stress modelling hypotheses. We then make an elastic prediction of the Hill equivalent stress and check whether or not this elastic prediction is inside the elastic domain. The latter information is stored in the boolean value b which will be false (no plastic loading) if the loading is elastic.
@InitLocalVariables{
const auto s = computeElasticPrediction();
const auto seq = sqrt(s|H*s);
b = seq-s0-R*p > 0;
}Finally, we describe how the implicit system and the computation of the jacobian is written in a code block introduced by the @Integrator keyword.
We use the following facts:
feel is initialized to deel).dfeel_ddeel is initialized to the identity tensor).feel by the StandardElasticity brick.Apart from those facts, the code is an almost direct translation of the mathematical expression described in previous sections and boils down to the following lines of code:
@Integrator{
if(b){
const auto seq = sqrt(sig|H*sig);
const auto iseq = 1/(max(seq,real(1.e-10*D(0,0))));
const auto n = iseq*H*sig;
// elasticity
feel += dp*n;
dfeel_ddeel += theta*dp*iseq*(H-(n^n))*D;
dfeel_ddp = n;
// plasticity
fp = (seq-s0-R*(p+theta*dp))/D(0,0);
dfp_ddp = -theta*(R/D(0,0));
dfp_ddeel = theta*(n|D)/D(0,0);
}
}