This article shows the implementation of an isotropic linear viscous model based on Skorohold-Olevsky Viscous Sintering (SOVS) model. This page is inspired by the paper of Lester [1]. For a more comprenhensive reading about the model, refer to the original paper of Olevsky [2].

This behaviour will be integrated using an implicit scheme.Such an example illustrates:

• the usage of StandardElasticity brick (see this page).
• how to use an auxiliary state variable and an integration variable to reduce the size of the implicit system.

The full implementation is available here

About the Skorohold-Olevsky Viscous Sintering (SOVS) model

SOVS model predicts sintering shrinkage and densification of ceramics and composites. It serves to enhance our ability to understand, predict, and control sintering. Additionally, it provides the potencial for more cost-effective net or near-net-shape manufacturing [3].

SOVS model has been succesfully used to predict the sintering of ZnO [4] and porous bilayered bodies [5, 6]. Regarding the ceramic tiles, viscosity models have been developed to predict the shrinkage and the residual stresses that appear when materials of the porcelain tile type are cooled [7, 8]. SOVS model can extend the functionality of those viscous models to predict the sintering, whose investigation is in process.

Description

The behaviour is described by a standard decomposition of the strain \(\underline{\varepsilon}^{\mathrm{to}}\) in an elastic and a viscous component, respectively denoted \(\underline{\varepsilon}^{\mathrm{el}}\) and \(\underline{\varepsilon}^{\mathrm{in}}\):

\[ \underline{\varepsilon}^{\mathrm{to}}= \underline{\varepsilon}^{\mathrm{el}}+ \underline{\varepsilon}^{\mathrm{in}} \qquad(1)\]

This equation can be written in rate form as follows:

\[ \underline{\dot{\varepsilon}}^{\mathrm{to}}= \underline{\dot{\varepsilon}}^{\mathrm{el}}+ \underline{\dot{\varepsilon}}^{\mathrm{in}} \qquad(2)\]

Elasticity

The stress \(\underline{\sigma}\) is related to the the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) by the standard Hooke law expressed using the Lamé coefficients:

\[ \underline{\sigma}= \lambda\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\, \underline{I} + 2\,\mu\,\underline{\varepsilon}^{\mathrm{el}} \]

Viscosity

Assuming linear viscosity, the inelastic strain rate, \(\underline{\dot{\varepsilon}}^{\mathrm{in}}\), may be written as:

\[ \underline{\dot{\varepsilon}}^{\mathrm{in}}= \frac{\underline{s}}{2\,\eta_{0}(T)\,\phi(\rho)} + \frac{{\mathrm{tr}{\left(\underline{\sigma}\right)}}-3\,\sigma_{s}(\rho)}{18\,\eta_{0}(T)\,\psi(\rho)}\underline{I} \qquad(3)\]

where:

• \(\eta_{0}(T)\) is the shear viscosity of the fully dense skeleton phase.
• \(\phi(\rho)\) is the normalized shear viscosity.
• \(\psi(\rho)\) is the normalized bulk viscosity.
• \(\sigma_s\) is the sintering stress.
  
• \(T\) is the temperature
• \(\rho\) is the relative density
• \(\underline{s}\) is the deviatoric stress tensor.

The dependency of the viscosities and the sintering stress used are the following:

\[ \phi(\rho)=a_{1}\rho^{b_{1}}, \qquad \psi(\rho)=a_{2}\frac{\rho^{b_{2}}}{(1-\rho)^{c_{2}}} \]

\[ \sigma_{s}(\rho)=\sigma_{s0}\overline{\sigma}_{s}(\rho), \qquad \overline{\sigma}_{s}(\rho)=a_{3}\rho^{b_{3}} \]

\[ \eta_{0}(T)=a_{4}e^{b_{4}/T} \]

where:

• \(a_{i}\), \(b_{i}\) and \(c_{i}\) are fitting coefficients.
• \(\sigma_{s0}\) is the local sintering stress given by \(3\alpha/r_{0}\).
• \(\alpha\) is the surface tension.
• \(r_{0}\) the average grain size.

By projecting Equation (3) on the hydrostatic and deviatoric spaces, one gets:

\[ \begin{aligned} {\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}}&=\frac{{\mathrm{tr}{\left(\underline{\sigma}\right)}}-3\,\sigma_{s}(\rho)}{6\,\eta_{0}(T)\,\psi(\rho)}\\ \underline{s}_{\underline{\dot{\varepsilon}}^{\mathrm{in}}}&=\frac{\underline{s}}{2\,\eta_{0}(T)\,\phi(\rho)}\\ \end{aligned} \qquad(4)\]

where \(\underline{s}_{\underline{\dot{\varepsilon}}^{\mathrm{in}}}\) stands for the deviatoric part of \(\underline{\dot{\varepsilon}}^{\mathrm{in}}\).

Evolution of the relative density

Finally, the evolution of the relative density is:

\[ \dot{\rho}=-\rho \, \, {\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}} \qquad(5)\]

Implicit scheme

The state variables are:

• the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\).
• the inelastic strain \(\underline{\varepsilon}^{\mathrm{in}}\).
• the relative density \(\rho\).

Their evolutions is governed by a a system of ordinary differential equations given by Equations (2), (3) and (5).

In its general form, a system of ordinary differential equations may be written as:

\[ \vec{\dot{X}}=G{\left(\vec{X},t\right)} \qquad(6)\]

where \(\vec{X}\) is the vector of the unknowns. Here \(\vec{X}\) is decomposed as follows:

\[ \vec{X}= \left( \begin{aligned} \underline{\varepsilon}^{\mathrm{el}}\\ \underline{\varepsilon}^{\mathrm{in}}\\ \rho \\ \end{aligned} \right) \]

In Equation (6), the time appears as a placeholder for external state variables that varies independently of the mechanics. Here the only external state behaviour is the temperature.

An implicit \(\theta\)-scheme turns this system of ordinary differential equations into a system of of non-linear equations by introducing the increment \(\Delta\,\vec{X}\) over a finite time step increment \(\Delta\,t\) as follows:

\[ f{\left(\Delta\,X\right)}=\Delta\,\vec{X}-\Delta\,t\,G{\left(\vec{X}+\theta\,\Delta\vec{X},t+\theta\,\Delta\,t\right)} \] where \(\theta\) is a numerical parameter (\(0<\theta\leq 1\)).

Again \(\Delta\,X\) and \(f{\left(\Delta\,X\right)}\) can be decomposed as follows:

\[ \Delta\,\vec{X}= \left( \begin{aligned} \Delta\,\underline{\varepsilon}^{\mathrm{el}}\\ \Delta\,\underline{\varepsilon}^{\mathrm{in}}\\ \Delta\,\rho \\ \end{aligned} \right) \quad\text{and}\quad f{\left(\Delta\,\vec{X}\right)}= \left( \begin{aligned} f_{\underline{\varepsilon}^{\mathrm{el}}} \\ f_{\underline{\varepsilon}^{\mathrm{in}}} \\ f_{\rho} \\ \end{aligned} \right) \]

This equation is solved by a Newton-Raphson algorithm.

In the following, the increment of the state varaibles, respectively denoted \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\), \(\Delta\,\underline{\varepsilon}^{\mathrm{in}}\) and \(\Delta\,\rho\) and he equations associated with those increment are respectively denoted \(f_{\underline{\varepsilon}^{\mathrm{el}}}\), \(f_{\underline{\varepsilon}^{\mathrm{in}}}\) and \(f_{\rho}\).

We now start by making an explicit derivation of \(f_{\underline{\varepsilon}^{\mathrm{el}}}\), \(f_{\underline{\varepsilon}^{\mathrm{in}}}\) and \(f_{\rho}\), and we then discuss how the system of implicit equations can be reduced.

Implicit equation associated with the elastic strain

The implicit equation associated with the the elastic strain derives from the split of strain (Equation (2)) as follows:

\[ f_{\underline{\varepsilon}^{\mathrm{el}}}=\Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,\underline{\varepsilon}^{\mathrm{in}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}} \qquad(7)\]

Implicit equation associated with the trace of the inelastic strain

Using Equation (3) can be discretised as follows:

\[ f_{\underline{\varepsilon}^{\mathrm{in}}}= \Delta\,\underline{\varepsilon}^{\mathrm{in}}- \Delta\,t\,\frac{{\left.\underline{s}\right|_{t+\theta\,\Delta\,t}}}{2\,\eta_{0}({\left.T\right|_{t+\theta\,\Delta\,t}})\,\phi{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}\right)}} - \Delta\,t\,\frac{{\left.{\mathrm{tr}{\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}}\right|_{t+\theta\,\Delta\,t}}-3\,\sigma_{s}{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}\right)}}{18\,\eta_{0}({\left.T\right|_{t+\theta\,\Delta\,t}})\,\psi{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}\right)}}\,\underline{I} \qquad(8)\]

where:

• \({\left.\underline{s}\right|_{t+\theta\,\Delta\,t}}=2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,\underline{s}_{{\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}= 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}{\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}\right)}}\,\underline{I}\right)}\)
• \({\left.\rho\right|_{t+\theta\,\Delta\,t}}={\left.\rho\right|_{t}}+\theta\,\Delta\,\rho\)

Evolution of the relative density

Equation (5) can be discretized as the other equations as follows:

\[ f_{\rho}=\Delta\,\rho+{\left.\rho\right|_{t+\theta\,\Delta\,t}}\,{\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{in}}\right)}} \qquad(9)\]

However, using this equation has a two major disadvantages:

1. The relative density \(\rho\) is bounded (i.e. positive) but it fairly complex to take this into account during the Newton resolution or event to enforce for large time increment.
2. It is trivial to eliminate this equation from the implicit system because the previous equation allows to express \(\Delta\,\rho\) as an explicit function of \(\Delta\,\underline{\varepsilon}^{\mathrm{in}}\). So keeping, Equation (9) leads to a less efficient integration scheme.

We now propose an alternative integration scheme for Equation (5) that remove the limitations associated with the first point and still allows to eliminate \(\Delta\,\rho\) from the implicit system.

Indeed, Equation (5) can be integrated exactly over the time step as follows:

\[ \begin{aligned} \dot{\rho}&=-\rho \, {\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}} \\ {{\displaystyle \frac{\displaystyle \dot{\rho}}{\displaystyle \rho}}}&=-{\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}} \\ \int_{t}^{t+dt}{{\displaystyle \frac{\displaystyle \dot{\rho}}{\displaystyle \rho}}}\,dt&=-\int_{t}^{t+dt}{\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}}\,dt \\ \log{\left({{\displaystyle \frac{\displaystyle {\left.\rho\right|_{t+\Delta\,t}}}{\displaystyle {\left.\rho\right|_{t}}}}}\right)}&=-{\mathrm{tr}{\left(\Delta\,\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}} \\ \end{aligned} \]

This leads to an explicit expression of \({\left.\rho\right|_{t+\Delta\,t}}\) as a function of \({\left.\rho\right|_{t}}\) and the trace of the inelastic strain increment \({\mathrm{tr}{\left(\Delta\,\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}}\):

\[ {\left.\rho\right|_{t+\Delta\,t}}={\left.\rho\right|_{t}}\,\exp{\left(-{\mathrm{tr}{\left(\Delta\,\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}}\right)} \qquad(10)\]

Similarly, the value of the relative density at \(t+\theta\,\Delta\,t\) is equal to:

\[ {\left.\rho\right|_{t+\theta\,\Delta\,t}}={\left.\rho\right|_{t}}\,\exp{\left(-{\mathrm{tr}{\left(\theta\,\Delta\,\underline{\dot{\varepsilon}}^{\mathrm{in}}\right)}}\right)} \qquad(11)\]

Equation (11) will be used to eliminate the relative density from the implicit system. Equation (10) will be used to compute the final value of the relative density.

About auxiliary state variables

Such a variable, whose value is kept from one time to another but is not part of the implicit system, is called an auxiliary state variable in MFront. A dedicated code block named @UpdateAuxiliaryStateVariables is meant to update those variables after the behaviour integration. At this stage, the state variables have been updated to their final values and the stress at the end of the time step has been computed.

Reduction of the implicit scheme

A close look at Equations (7), (8), and (11) shows that it can be efficient to eliminate part of the unknowns and keep the two following variables \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\) and \({\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{in}}\right)}}\).

In the following, \(e\) denotes \({\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{in}}\right)}}\): \[ e={\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{in}}\right)}} \]

About integration variables

Here \(e\) is not an internal state variable, because there is no reason to keep its value from one step to the other. Such a variable is called in MFront an integration variable.

In pratice, for this integration variable \(e\), MFront automatically defines:

• a variable e meant to store the value of \(e\) at the beginning of the time step if required (this is not the case here). This initial value is meant to be computed in the @InitLocalVariables block from the values of the state variables or auxiliary state variables.
• its increment de
• an associated equation in the implicit system called fe.
• the jacobian blocks associated with fe (here dfe_dde and dfe_ddeel).

Combining Equations (7), (8), and (11) leads to the new system of equations:

\[ \left\{ \begin{aligned} f_{\underline{\varepsilon}^{\mathrm{el}}} &= \Delta\,\underline{\varepsilon}^{\mathrm{el}}+ \Delta\,t\,\frac{{\left.\underline{s}\right|_{t+\theta\,\Delta\,t}}}{2\,\eta_{0}({\left.T\right|_{t+\theta\,\Delta\,t}})\,\phi{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}} + {{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,\Delta\,e\,\underline{I} -\Delta\,\underline{\varepsilon}^{\mathrm{to}}\\ f_{e}&=\Delta\,e- 3\,\Delta\,t\,\frac{{\left.{\mathrm{tr}{\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}}\right|_{t+\theta\,\Delta\,t}}-3\,\sigma_{s}{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}}{18\,\eta_{0}({\left.T\right|_{t+\theta\,\Delta\,t}})\,\psi{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}} \end{aligned} \right. \qquad(12)\]

where \({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\) is given by Equation (11).

For the sake of clarity, the previous equation can be rewritten as follows: \[ \left\{ \begin{aligned} f_{\underline{\varepsilon}^{\mathrm{el}}} &= \Delta\,\underline{\varepsilon}^{\mathrm{el}}+ \Delta\,t\,A{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}\,{\left.\underline{s}\right|_{t+\theta\,\Delta\,t}}+ {{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,\Delta\,e\,\underline{I} -\Delta\,\underline{\varepsilon}^{\mathrm{to}}\\ f_{e}&=\Delta\,e- 3\,\Delta\,t\,{\left(B{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}\,{\mathrm{tr}{\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}}-C{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}\right)} \end{aligned} \right. \qquad(13)\]

where:

\[ \left\{ \begin{aligned} A{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}&={{\displaystyle \frac{\displaystyle 1}{\displaystyle {2\,\eta_{0}({\left.T\right|_{t+\theta\,\Delta\,t}})\,\phi{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}}}}}\\ B{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}&={{\displaystyle \frac{\displaystyle 1}{\displaystyle 18\,\eta_{0}({\left.T\right|_{t+\theta\,\Delta\,t}})\,\psi{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}}}}\\ C{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}&=3\,B{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}\,\sigma_{s}{\left({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\right)}\\ \end{aligned} \right. \]

Since the deviatoric part of the inelastic strain is linear with the stress, the System (12) may be further reduced to a scalar equation with \(\Delta\,e\) as the only unknowns. However, this considerably complicates the implementation of the behaviour, and in particular, the derivation of the consistent tangent operator. Hence, for the sake simplicity, this reduction to a scalar equation is not considered here. However, the reader may refer to [9] for an example of such a reduction to a scalar equation may be done in MFront.

To solve System (13), one must compute its jacobian as follows:

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}&=\underline{\underline{\mathbf{I}}}+2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,\theta\,A\,{\left(\underline{\underline{\mathbf{I}}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\underline{I}\,\otimes\,\underline{I}\right)}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,e}}&={{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,\underline{I}+\Delta\,t\,{\displaystyle \frac{\displaystyle \partial A}{\displaystyle \partial {\left.\rho\right|_{t+\theta\,\Delta\,t}}}}\,{\displaystyle \frac{\displaystyle \partial {\left.\rho\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,e}}\,{\left.\underline{s}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{e}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}&=-3\,\theta\,\Delta\,t\,B\,{\left(3\,{\left.\lambda\right|_{t+\theta\,\Delta\,t}}+2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\right)}\,\underline{I}\\ {\displaystyle \frac{\displaystyle \partial f_{e}}{\displaystyle \partial \Delta\,e}}&=1-3\,\Delta\,t{\left({\displaystyle \frac{\displaystyle \partial B}{\displaystyle \partial {\left.\rho\right|_{t+\theta\,\Delta\,t}}}}\,{\mathrm{tr}{\left({\left.\sigma\right|_{t+\theta\,\Delta\,t}}\right)}}-{\displaystyle \frac{\displaystyle \partial C}{\displaystyle \partial {\left.\rho\right|_{t+\theta\,\Delta\,t}}}}\right)}\,{\displaystyle \frac{\displaystyle \partial {\left.\rho\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,e}}\\ \end{aligned} \right. \]

where we omitted the dependency of \(A\), \(B\) and \(C\) to \({\left.\rho\right|_{t+\theta\,\Delta\,t}}{\left(\Delta\,e\right)}\).

The formal derivatives of \(A\), \(B\) and \(C\) are a bit tiedous and won’t be detailled here.

Implementation

Choice of domain specific language

While not mandatory (the @DSL keyword can be place anywhere in the file), its is convenient to start the implementation by declaring the domain specific language to be used. For an integration by a \(\theta\)-scheme, the Implicit domain specific language is choosen:

@DSL Implicit;

Name of the behaviour

The @Behaviour keyword is used to give the name of the behaviour.

@Behaviour sovs;

Metadata

The following instructions give some information about the author, the date and some description of the behaviour.

@Author Joan Balaguer Osuna;
@Date 04/2019;
@Description{
  "Verification of the Skorohod - Olevsky Viscous Sintering(SOVS) Model."
  "Brian Lester"
  "Sandia National Laboritories"
  "November 16, 2017"
}

Numerical parameters

The following instruction changes the default value of the stopping criterion used by the Newton-Raphson method:

@Epsilon 1.e-16;

This can be changed at runtime by modifying the epsilon parameter.

The default value of the \(\theta\) parameter is left unchanged (its default value is one half). This can be changed at runtime by modifying the theta parameter.

Usage of the StandardElasticity brick

The implicit scheme used satisfies the requirements of the StandardElasticity brick as described here.

The StandardElasticity brick which provides:

• Automatic computation of the stress tensor at various stages of the behaviour integration.
• Automatic computation of the consistent tangent operator.
• Automatic support for plane stress and generalized plane stress modelling hypotheses (The axial strain is defined as an additional state variable and the associated equation in the implicit system is added to enforce the plane stess condition).
• Automatic addition of the standard terms associated with the elastic strain state variable. The usage of the StandardElasticity is introduced as follows:

@Brick StandardElasticity;

Supported modelling hypothesis

Thanks to the StandardElasticity brick, all the modelling hypotheses can be supported. The following statement, starting with the @ModellingHypotheses, enables all the modelling hypotheses:

@ModellingHypotheses {".+"};

Elastic properties

@ElasticMaterialProperties {123.7e9,0.356};

Material parameters

@Parameter a1 = 1;
@Parameter b1 = 2;

@Parameter a2 = 0.666666666666666;
@Parameter b2 = 3;
@Parameter c2 = 1;

@Parameter a3 = 1;
@Parameter b3 = 2;
@Parameter sig_s0 = 3.81e6;

@Parameter eta0_A = 0.0000166;
@Parameter eta0_B = 48259.583;

The relative density as an auxiliary state variable

@AuxiliaryStateVariable real r;
r.setEntryName("RelativeDensity");

The trace of the viscoplastic strain as an integration variable

//! trace of the viscoplastic strain
@IntegrationVariable strain e;

Local variable

In MFront, an integration variable is defined to store a variable and use it in various code block.

Here a local variable called eta0 is declared to later computes \(\eta_{0}{\left(T+\theta\,\Delta\,T\right)}\) once for all.

@LocalVariable stress eta0;

Initialisation a the local variable

The @InitLocalVariables code block is called before the behaviour integration. It is used here to compute \(\eta_{0}{\left(T+\theta\,\Delta\,T\right)}\) once for all:

@InitLocalVariables {
  eta0 = eta0_A*exp(eta0_B/(T + theta * dT));
}

Implicit system implementation

The implementation of the implicit system and its derivative is done in the @Integrator code block:

@Integrator{

First, we define the identity for symmetric second order tensors, as follows:

  constexpr const auto id = StrainStensor::Id();

The following lambda function (this function is local to the @Integration code block) will be useful later to computed some of the derivatives efficiently:

  const auto is_zero = [=](const real v) {
    constexpr const real eps0 = 1e-8;
    return (abs(v) - eps0 < 0) ? true : false;
  };

An estimation of the relative density \({\left.\rho\right|_{t+\theta\,\Delta\,t}}\) is computed from the current estimation of \(\Delta\,e\):

  const auto r_ = r*exp(-theta * de);

\(\phi\), \(\psi\) and \(sigma_s\) can now be computed:

  const auto phi = a1 * pow(r_, b1);
  const auto psi = a2 * pow(r_, b2) / pow(1 - r_, c2);
  const auto sig_s = sig_s0 * a3 * pow(r_, b3);

\(A\), \(B\) and \(C\) can now be computed:

  const auto A = 1 / (2 * eta0 * phi);
  const auto B = 1 / (18 * eta0 * psi);
  const auto C = 3 * sig_s * B;

To express the viscoplastic flow, one need the deviatoric part of the stress:

  const auto s = deviator(sig);

Computation of the stress

sig stands for the current estimation of the stress at \(t+\theta\,\Delta\,t\) using \({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}\). It is automatically computed by the StandardElasticity brick.

The implicit system is then readily implemented:

  feel += de/3 * id + A * s;
  fe   -= 3 * dt * (B * trace(sig) - C);

Implicit initialisation of the implicit equations

This implementation takes into account the fact that the Implicit DSL automatically initializes feel to the current estimation of \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\) and fe to the current estimation of \(\Delta\,e\). Moreover, the StandardElasticity brick automatically subtracts \(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\) to feel.

The previous lines are thus equivalent to

  feel = deel - deto + de/3 * id + A * s;
  fe   = de + 3 * dt * (B * trace(sig) - C);

We now computes the derivatives of \(A\), \(B\) and \(C\) with respect to \(\Delta\,e\):

  const auto dr__dde = -theta * r_;
  const auto dphi_dr_ = is_zero(r_) ? real(0) : b1 * phi / r_;
  const auto dpsi_dr_ =
      is_zero(r_) ? real(0) : psi * ((c2 - b2) * r_ + b2) / (r_ * (1 - r_));
  const auto dpsi_dde = dpsi_dr_ * dr__dde;
  const auto dphi_dde = dphi_dr_ * dr__dde;
  const auto dA_dphi = -1 / (2 * eta0 * phi * phi);
  const auto dB_dpsi = -1 / (18 * eta0 * psi * psi);
  const auto dA_dde = dA_dphi * dphi_dde;
  const auto dB_dde = dB_dpsi * dpsi_dde;
  const auto dsig_s_dr_ = is_zero(r_) ? real(0) : b3 * sig_s / r_;
  const auto dsig_s_dde = dsig_s_dr_ * dr__dde;
  const auto dC_dde = 3 * sig_s * dB_dde + 3 * B * dsig_s_dde;

The termes of the jacobian can now be readily implemented:

  dfeel_ddeel += 2 * mu * theta * A * Stensor4::K();
  dfeel_dde    = id / 3 + dA_dde * s;
  dfe_dde     -= 3 * dt * (dB_dde * trace(sig) + dC_dde);
  dfe_ddeel    =  -3 * (3 * lambda + 2 * mu) * theta * B * dt * id;

Implicit initialization of the jacobian blocks

The Implicit DSL automatically initializes the jacobian to the identity. The previous lines are thus equivalent to:

  dfeel_ddeel = Stensor4::Id() + 2 * mu * theta * A * Stensor4::K();
  dfeel_dde   = id / 3 + dA_dde * s;
  dfe_dde     = 1 -3 * dt * (dB_dde * trace(sig) + dC_dde);
  dfe_ddeel   =  -3 * (3 * lambda + 2 * mu) * theta * B * dt * id;

The following curly bracket closes the @Integrator block:

}

The relative density is updated after the behaviour integration as follows:

@UpdateAuxiliaryStateVariables {
  r *= exp(-de);
}

References