Plane stress in explicit schemes

[RonanRiverie]

Hello Thomas,

I encountered several convergence issues when running simulations with the EKMC model, integrated with an implicit scheme. Therefore, I switched to an explicit one and tested the different algorithms proposed by Mfront. However, when I use the RK54 or RKCastem methods for uniaxial tension in plane stress, I got this error in Castem (even when using a large number or time increments):

When I use the RK4 algorithm, I don’t have this error.

I read the documentations about the Umat interface and how to consider plane stress hypothesis. If I understand well, since the constitutive law does not treat the plane stress case, Castem runs its own loop to impose the hypothesis. Nevertheless, I don’t know how to manage this inner Castem loop.

The Mfront code:

@DSL RungeKutta;
@Behaviour ekmc;

@Algorithm rk54;

@RequireStiffnessTensor;

@MaterialProperty stress young ;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu ;
nu.setGlossaryName("PoissonRatio");
@MaterialProperty strainrate de0;
@MaterialProperty real Ea;
@MaterialProperty real Va;
@MaterialProperty stress P1;
@MaterialProperty real P2;
@MaterialProperty real omega;
@MaterialProperty stress R0;
@MaterialProperty stress Q;
@MaterialProperty real b;
@MaterialProperty stress H;
@MaterialProperty real nta;
@MaterialProperty real palpha;

@StateVariable strain p;
p.setGlossaryName("EquivalentViscoplasticStrain");

@StateVariable time ta;
ta.setEntryName("StrainAgingTime");

@LocalVariable strainrate de;

@LocalVariable real sinh0;

@InitLocalVariables {
  constexpr auto kb = PhysicalConstants::kb;
  de = de0 * exp(-Ea / (kb * T));
  sinh0 = kb * T / Va;
}

@ComputeStress {
  sig = D*eel;
}

@Derivative {
  const auto seq = sigmaeq(sig);
  Stensor n = Stensor(0);
  if (seq > 1.e-12) {
    n = 1.5 * deviator(sig) / seq;
  }
  p = max(p, strain(0));
  ta = max(ta, time(0));
  const auto Cs = 1 - exp(-P2 * pow(p, palpha) * pow(ta, nta));
  const auto Ra = P1 * Cs;
  const auto R = R0 + Q * (1 - exp(-b * p)) + H * p + Ra;
  const auto f = max(seq - R, stress{0});
  deel = deto - dp * n;
  dp = de * sinh(f / sinh0);
  dta = 1 - (ta/(omega/dp));
}

@TangentOperator {
  Dt = D;
}

Thank you for your help

[thelfer]

Hi @RonanRiverie,

The plane stress hypothesis is supported using an external loop around the behaviour integration which is not optimal and maybe not robust enough for your needs: the behaviour is integrated under the generalized plane stress assumption and the axial strain is determined iteratively.

It is much better to take the plane stress hypothesis inside the behaviour in an implicit scheme.

[thelfer]

If the behaviour do not work even in MTest, there are little chance to make it workin Cast3M. Would you share your implementation and test ?

[thelfer]

Yes, it would be much better in implicit but I don't manage to make the computations converge in implicit

Let us try to make the implicit work :)

[RonanRiverie]

Yes sure.
Here is the Mfront code:

@DSL Implicit;
@Behaviour ekmc;

@Algorithm NewtonRaphson_NumericalJacobian;
@PerturbationValueForNumericalJacobianComputation 1e-8;
@Epsilon 1e-14;
@Theta 1;

@ModellingHypotheses {".+"};

@Brick StandardElasticity;

@MaterialProperty strainrate de0;
@MaterialProperty real Ea;
@MaterialProperty real Va;
@MaterialProperty stress P1;
@MaterialProperty real P2;
@MaterialProperty real omega;
@MaterialProperty stress R0;
@MaterialProperty stress Q;
@MaterialProperty real b;
@MaterialProperty stress H;
@MaterialProperty real nta;
@MaterialProperty real palpha;

@StateVariable strain p;
p.setGlossaryName("EquivalentViscoplasticStrain");

@AuxiliaryStateVariable time ta;
ta.setEntryName("StrainAgingTime");

@LocalVariable time dta;

@LocalVariable strainrate de;

@LocalVariable real sinh0;

@InitLocalVariables {
  constexpr auto kb = PhysicalConstants::kb;
  const auto T_theta = T + theta * dT;
  de = de0 * exp(-Ea / (kb * (T_theta)));
  sinh0 = kb * T_theta / Va;
}

@Integrator {
  const auto seq = sigmaeq(sig);
  Stensor n = Stensor(0);
  if (seq > 1.e-12) {
    n = 1.5 * deviator(sig) / seq;
  }
  const auto p_theta = max(p + theta * dp, 0.);
  dta = (omega * dt - ta * dp) / (omega + theta * dp);
  const auto ta_theta = max(ta + theta * dta, time(0));
  const auto Cs = 1 - exp(-P2 * pow(p_theta, palpha) * pow(ta_theta, nta));
  const auto Ra = P1 * Cs;
  const auto R = R0 + Q * (1 - exp(-b * p_theta)) + H * p_theta + Ra;
  const auto f = max(seq - R, stress{0});
  feel += dp * n;
  fp -= de * sinh(f / sinh0) * dt;
}

@UpdateAuxiliaryStateVariables {
  ta += dta;
}

And on of the mtest test case:

@Behaviour<umat> 'src\libUmatBehaviour.dll' 'umatekmc' ;

@AccelerationAlgorithm 'Cast3M' ;

@MaterialProperty<constant> 'YoungModulus' 209.e9 ;
@MaterialProperty<constant> 'PoissonRatio' 0.3 ;
@MaterialProperty<constant> 'MassDensity' 0.0 ;
@MaterialProperty<constant> 'ThermalExpansion' 0.0 ;
@MaterialProperty<constant> 'ReferenceTemperatureForThermalExpansionCoefficient' 0.0 ;
@MaterialProperty<constant> 'ReferenceTemperatureForThermalExpansion' 0.0 ;

@MaterialProperty<constant> 'R0' 200.e6 ;
@MaterialProperty<constant> 'Q' 270.e6 ;
@MaterialProperty<constant> 'b' 23.3 ;
@MaterialProperty<constant> 'H' 440.e6 ;
@MaterialProperty<constant> 'de0' 1.0e-07 ;
@MaterialProperty<constant> 'Ea' 3.52e-22 ;
@MaterialProperty<constant> 'Va' 9.48e-28 ;
@MaterialProperty<constant> 'nta' 0.66 ;
@MaterialProperty<constant> 'P1' 120.e6 ;
@MaterialProperty<constant> 'omega' 0.0015 ;
@MaterialProperty<constant> 'P2' 0.0003 ;
@MaterialProperty<constant> 'palpha' 0.1 ;
@InternalStateVariable 'StrainAgingTime' 5000000. ;
@ExternalStateVariable 'Temperature' 293.15 ;

@Real 'edot' 1.e-2 ;
@Real 'EXXmax' 0.1 ;
@Real 'tmax' 'EXXmax/edot';	
@ImposedStrain 'EXX' { 0. :  0., 'tmax' : 'EXXmax'};
@Times {0., 'tmax' in 50000};

I tried to change, the epsilon parameter or the maximum number of iterations but it does not help. Actually, for some other values of the parameter, it converges but they are not those from which I am trying to compare in order to check my implementation with a paper. So, it seems that the values of the parameters do not help.

This constitutive law is supposed to reproduce Lüders or Portevin- Le Chatelier instabilities on structures, so it is highly non-linear.

[thelfer]

In your code, Cs's derivative is singular for $p=0$. This fixes the issue:

@Integrator {
  constexpr auto peps = strain{1e-8};
  ....
  const auto Cs =
      1 - exp(-P2 * pow(p_theta + peps, palpha) * pow(ta_theta, nta));
  ...
}

With this, the MTest test converges smoothly in 50 time steps

[RonanRiverie]

I have just tested to add these lines and both the mtest cases and the simulation on the specimen are running with few iterations. I will check if it also works for orther simulations with other values of the parameters.

For the singulartity, I understand that the derivative with respect to $p$ can indeed gives a division by $0$ when $p=0$ for the selected value of $\alpha$ and this derivative is taken into account in the jacobian ?

[thelfer]

I have just tested to add these lines and both the mtest cases and the simulation on the specimen are running with few iterations.

Great

this derivative is taken into account in the jacobian ?

Not really, the numerical jacobian is however very bad in this case, which leads to oscillations in the Newton algorithm

[RonanRiverie]

Not really, the numerical jacobian is however very bad in this case, which leads to oscillations in the Newton algorithm

So, I don't really understand where this derivative takes place

[thelfer]

You are using a Newton algorithm with a numerical jacobian.
MFront thus evaluates the jacobian using a finite centered difference scheme (since you have 7 internal state variables in 3D, it calls the residual 1 + 7 * 2 = 15 times per iterations). Since the real jacobian is not defined (infinite), the numerical jacobian is very bad and leads to oscillations of the algorithm

[RonanRiverie]

I am trying to implement the analytical jacobian. However, I have an issue with $\frac{\partial f_p}{\partial \Delta p}$ since it requires to compute ${\frac{\partial t_a\left(p\right)}{\partial p}}$ but $t_a$ is itself a function of $p$ through the differential equation. For the other terms, I followed this page:
https://thelfer.github.io/tfel/web/isotropicplasticityamstrongfrederickinematichardening.html

The equations are:

$$ \begin{cases} f_{el} : \Delta \underline{\varepsilon}^{\mathrm{el}} - \Delta \underline{\varepsilon}^{\mathrm{tot}} - \Delta p \underline{n} = 0 \\ f_{p} : \Delta p - \underbrace{\varepsilon_0 \exp{\left(-\frac{E_a}{k_B T}\right)}\sinh{\left(\frac{V_a \max(f(\sigma), 0)}{k_B T}\right)}}_{g\left(f\right)} \Delta t = 0 \end{cases} $$

$$ f\left(\underline{\sigma}\right) =\sigma_{eq} - R\left(p\right) - R_a\left(p, t_a\right) $$

$$ \dot{t}_a = 1 - \frac{\dot{p}}{\omega}t_a $$

$$ \frac{\partial f_p}{\partial \Delta p} = 1 - \Delta t \frac{\partial g}{\partial \Delta p} $$

$$ \frac{\partial g}{\partial \Delta p} = \frac{\partial g}{\partial f} \frac{\partial f}{\partial \Delta p} = g'\left(f\right)\frac{\partial f}{\partial \Delta p} $$

$$ \frac{\partial f}{\partial \Delta p} = \frac{\partial}{\partial \Delta p}\left(-R\left(p\right) - R_a\left(p\right)\right) = \frac{\partial p}{\partial \Delta p} \frac{\partial}{\partial p}\left(-R\left(p\right) - R_a\left(p\right)\right) =- \theta \left(\frac{\partial R}{\partial p} +\frac{\partial R_a}{\partial p}\right) $$

$$ \frac{\partial R_a}{\partial p} = -P_1 P_2 \exp\left(-P_2 p^\alpha t_a\left(p\right)^n\right) \left[ \alpha p^{\alpha - 1} t_a\left(p\right)^n + n p^\alpha t_a\left(p\right)^{n-1} \frac{\partial t_a\left(p\right)}{\partial p} \right] $$

So for the moment, I have just replaced $\frac{\partial t_a\left(p\right)}{\partial p}$ by $\frac{\Delta t_a}{\Delta p}$ but it does not seem to be the best idea. The results seem ok but the convergence is degraded since mtest is sometimes diving the time step by two.

I have checked with the @CompareToNumericalJacobian tool and there are some differences on the derivatives, especially on the one I am concerned with, i.e. $\frac{\partial f_p}{\partial \Delta p}$.

Here is the mfront code:

@DSL Implicit;
@Behaviour ekmc;

@Algorithm NewtonRaphson;
// @CompareToNumericalJacobian true;
// @JacobianComparisonCriterion 1.e-6;
// @Algorithm NewtonRaphson_NumericalJacobian;
// @PerturbationValueForNumericalJacobianComputation 1e-8;
@Epsilon 1e-14;
@Theta 1;

@ModellingHypotheses {".+"};

@Brick StandardElasticity;

@MaterialProperty stress young; 
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;   
nu.setGlossaryName("PoissonRatio");
@MaterialProperty strainrate de0;
@MaterialProperty real Ea;
@MaterialProperty real Va;
@MaterialProperty stress P1;
@MaterialProperty real P2;
@MaterialProperty real omega;
@MaterialProperty stress R0;
@MaterialProperty stress Q;
@MaterialProperty real b;
@MaterialProperty stress H;
@MaterialProperty real nta;
@MaterialProperty real palpha;

@StateVariable strain p;
p.setGlossaryName("EquivalentViscoplasticStrain");

@AuxiliaryStateVariable time ta;
ta.setEntryName("StrainAgingTime");

@LocalVariable time dta;

@LocalVariable strainrate de;

@LocalVariable real sinh0;

@LocalVariable stress mu;

@InitLocalVariables {
  constexpr auto kb = PhysicalConstants::kb;
  const auto T_theta = T + theta * dT;
  de = de0 * exp(-Ea / (kb * (T_theta)));
  sinh0 = kb * T_theta / Va;

  mu = computeMu(young,nu);
}

@Integrator {
  const auto seq = sigmaeq(sig);
  constexpr const real eps = real(1.e-12);
  const auto iseq = 1/(max(seq, eps*young));
  constexpr auto peps = strain{1e-8};
  const auto n = 1.5 * deviator(sig) * iseq;
  const auto p_theta = max(p + theta * dp, peps);
  dta = (omega * dt - ta * dp) / (omega + theta * dp);
  const auto ta_theta = max(ta + theta * dta, time(0));
  const auto Cs_exp =  exp(-P2 * pow(p_theta, palpha) * pow(ta_theta, nta));
  const auto Cs = 1 - Cs_exp;
  const auto Ra = P1 * Cs;
  const auto R = R0 + Q * (1 - exp(-b * p_theta)) + H * p_theta + Ra;
  const auto f = max(seq - R, stress{0});
  feel += dp * n;
  fp -= de * sinh(f / sinh0) * dt;
  // Jacobian matrix
  const auto Jmn = eval(Stensor4::M()-(n^n));
  const auto dg_df = de * (1/sinh0) * cosh(f /sinh0);
  const auto dR = Q * b * exp(-b * p_theta) + H;
  const auto dCs_dp = palpha * pow(p_theta, (palpha - 1)) * pow(ta_theta, nta);
  const auto dta_dp = dta/(max(dp, peps));
  const auto dCs_dta = nta * pow(p_theta, palpha) * pow(ta_theta, (nta - 1)) * dta_dp;
  const auto dRa = P1 * (-P2*Cs_exp * (dCs_dp + dCs_dta));
  dfeel_ddeel += 2*mu*theta*dp*Jmn*iseq ;
  dfeel_ddp = n;
  dfp_ddeel = - dg_df*dt * 2*mu*theta*n;
  dfp_ddp += theta*dg_df * (dR + dRa) * dt;

  // std::cout << std::endl;
  // std::cout << "dg_df: " << dg_df << std::endl;
  // std::cout << "dR: " << dR << std::endl;
  // std::cout << "dRa: " << dRa << std::endl;
  // std::cout << "dCs_dp: " << dCs_dp << std::endl;
  // std::cout << "dCs_dta: " << dCs_dta << std::endl;
  // std::cout << "p_theta: " << p_theta << std::endl;
  // std::cout << "dta_dp: " << dta_dp << std::endl;
  // std::cout << "pow(p_theta, palpha): " << pow(p_theta, palpha) << std::endl;
  // std::cout << "pow(ta_theta, (nta - 1)): " << pow(ta_theta, (nta - 1)) << std::endl;
  // std::cout << "dta: " << dta << std::endl;
  // std::cout << "dp: " << dp << std::endl;
  // std::cout << std::endl;
}

@UpdateAuxiliaryStateVariables {
  ta += dta;
}

[RonanRiverie]

Actually it was when imposing 500 time steps, mtest was using 512 with 8 substeps while before it was not dividing the time step but by imposing 50 time steps mtest keeps 50 periods with the analytical jacobian, so everything seems alright.

If you want the mtest file:

@Behaviour<umat> 'src\libUmatBehaviour.dll' 'umatekmc' ;
@AccelerationAlgorithm 'Cast3M' ;

@MaterialProperty<constant> 'YoungModulus' 202.e9 ;
@MaterialProperty<constant> 'PoissonRatio' 0.3 ;
@MaterialProperty<constant> 'MassDensity' 0.0 ;
@MaterialProperty<constant> 'ThermalExpansion' 0.0 ;
@MaterialProperty<constant> 'ReferenceTemperatureForThermalExpansionCoefficient' 0.0 ;
@MaterialProperty<constant> 'ReferenceTemperatureForThermalExpansion' 0.0 ;

@MaterialProperty<constant> 'R0' 184.e6 ;
@MaterialProperty<constant> 'Q' 212.e6 ;
@MaterialProperty<constant> 'b' 24.8 ;
@MaterialProperty<constant> 'H' 440.e6 ;
@MaterialProperty<constant> 'de0' 1.0e-07 ;
@MaterialProperty<constant> 'Ea' 3.52e-22 ;
@MaterialProperty<constant> 'Va' 1.38e-27 ;
@MaterialProperty<constant> 'nta' 0.66 ;
@MaterialProperty<constant> 'P1' 110.e6 ;
@MaterialProperty<constant> 'omega' 0.0015 ;
@MaterialProperty<constant> 'P2' 0.17 ;
@MaterialProperty<constant> 'palpha' 0.17 ;
@InternalStateVariable 'StrainAgingTime' 300. ;
@ExternalStateVariable 'Temperature' 423.15 ;

@Real 'edot' 1.e-4 ;
@Real 'EXXmax' 0.1 ;
@Real 'tmax' 'EXXmax/edot';	
@ImposedStrain 'EXX' { 0. :  0., 'tmax' : 'EXXmax'};
@Times {0., 'tmax' in 500};

[thelfer]

Normally, everything is better when using small time steps !
I'll have a look. Thanks

[RonanRiverie]

By the way, in the current version, there is not any distinction of cases with $f>0$ because I was simply thinking that when we are in elasticity we get $\dot{p}=0$ but it might be more optimal with the distinction. So, I tried that with the example above and now I don't get the extra substeps when imposing 500 time steps by having 500 periods.
Maybe some iterations are launched while not required in elasticity and so with this correction the jacobian is computed only in the plastic domain.

[thelfer]

I changed p0 from 1e-8 to 1e-6 and made some optimisations. It seems more robust:

ekmc.zip

[RonanRiverie]

Ok, I see your changes, I just don't get the point with the if since with the definition of $f$ using $\max{}$, it is always positive but I might be wrong.

Finally, I just have a last question, is it correct to do like that to be able to save $\dot{p}$ as a variable, so thus I can get its value either in mtest or in CAST3M:
So, first, I just add a new auxiliary state variable,

@AuxiliaryStateVariable strain dp_dt;
dp_dt.setEntryName("EquivalentViscoplasticStrainRate");

an then, I compute its value in the @integrator block dp_dt = dp / dt;