Resolving mixed loading for time independent plasticity and damage

Mostly, the so-called function resolve_flux_grad(E, deel, deto, devi) is used in Zmat routines to compute the increment of elastic strain eel knowing the Hooke’s tensor E, the increment of the total strain eto and the increment of the plastic strain evi. However, this function is not suitable in the case of time independent plasticity and it is even less if damage is introduced into the model.

The problem to be solved is the determination of unknown components of the total strain rate based on the Hooke’s law

\[\dot\sigma=(1-D)E\left(\dot\varepsilon-\dot{\ten{\varepsilon}}^p\right) -E\left(\ten{\varepsilon}-\ten{\varepsilon}^p\right)\dot{D}\]

and the consistency condition

\[(1-D)nE\dot{\varepsilon}=\left(H+nEn\right)\dot{p}\]

A linear system is built and solved on the light of the two previous equations. The unkowns are the strain components that are not involved in the loading.

CODE	DESCRIPTION
resolve_flux_grad_consistency(E,norm,H,deto,devcum,eel,D,dD dp)

Numerical backgrounds

In the following, subscripts “k” and “u” refer to as “known” and “unknown” total strain respectively. Thus, vectors containing total strain components and stress components together with the Hooke’s tensor write:

\[\begin{split}\dot\varepsilon= \left(\begin{array}{l} \dot{\varepsilon}^k: \textrm{known}\\ \dot{\varepsilon}^u: \textrm{unknown} \end{array} \right) \qquad \dot\sigma= \left(\begin{array}{l} \dot\sigma^k: \textrm{unknown}\\ \dot\sigma^u: \textrm{known} \end{array} \right) \qquad E= \left(\begin{array}{ll} E_{kk} & E_{ku}\\ E_{uk} & E_{uu} \end{array} \right)\end{split}\]

The linear system allowing to compute the strain rates \(\dot{\varepsilon}^u\) and the plastic multiplier \(\dot p\) is presented in the case of time independent model without damage and then including damage:

• time independent model without damage

  The Hooke’s law writes:

  \[\dot\sigma=E\left(\dot\varepsilon-\dot{p}n\right)\]
  \[\begin{split}\left(\begin{array}{l} \dot{\sigma}^k\\ \dot{\sigma}^u \end{array} \right) = \left(\begin{array}{ll} E_{kk} & E_{ku}\\ E_{uk} & E_{uu} \end{array} \right) \left(\begin{array}{l} \dot{\varepsilon}^k- \dot{p}n^k\\ \dot{\varepsilon}^u- \dot{p}n^u \end{array} \right) \qquad\end{split}\]
  \[\dot{\sigma}^u=E_{uk}\dot{\varepsilon}^k+E_{uu}\dot{\varepsilon}^u -\dot{p}\left(E_{uk}n^k+E_{uu}n^u\right)\]

  “u” équations are then be obtained:

  \[\boxed{ E_{uu}\dot{\varepsilon}^u-\left(En\right)^u\dot{p}=\dot{\sigma}^u-E_{uk}\dot{\varepsilon}^k }\]

  The consistency condition \(\dot f = 0\) gives:

  \[n\dot{\sigma}=H\dot{p} \qquad n\left(E\dot{\varepsilon}\right)=\left(H+nEn\right)\dot{p}\]
  \[\begin{split}\left(\begin{array}{ll} n^{k}\\ n^{u} \end{array} \right) \left(\begin{array}{ll} E_{kk} & E_{ku}\\ E_{uk} & E_{uu} \end{array} \right) \left(\begin{array}{l} \dot{\varepsilon}^k\\ \dot{\varepsilon}^u \end{array} \right) = \left(H+nEn\right)\dot{p}\end{split}\]
  \[\begin{split}\left(\begin{array}{ll} n^{k}\\ n^{u} \end{array} \right) \left(\begin{array}{ll} E_{kk}\dot{\varepsilon}^k+ E_{ku}\dot{\varepsilon}^u\\ E_{uk}\dot{\varepsilon}^k+ E_{uu}\dot{\varepsilon}^u \end{array} \right) = \left(H+nEn\right)\dot{p}\end{split}\]

  and leads to the scalar equation:

  \[\boxed{ \left(n^{k}E_{ku}+n^{u}E_{uu}\right)\dot{\varepsilon}^u-\left(H+nEn\right)\dot{p}= -\left(n^{k}E_{kk}+n^{u}E_{uk}\right)\dot{\varepsilon}^k }\]

  The system to be solved in the case of time independent plasticity is as follows:

  \[\begin{split}\left(\begin{array}{ll} H_{cc} & H_{cp}\\ H_{pc} & H_{pp} \end{array} \right) \left(\begin{array}{l} \dot{\varepsilon}^u\\ \dot{p} \end{array} \right) = \left(\begin{array}{l} b_c\\ b_p \end{array} \right)\end{split}\]
  \[\begin{split}\begin{aligned} H_{cc}&=E_{uu}\\ H_{cp}&=-\left(En\right)^u\\ H_{pc}&=\left(n^{k}E_{ku}+n^{u}E_{uu}\right)\\ H_{pp}&=-\left(H+nEn\right)\\ b_c&=\dot{\sigma}^u-E_{uk}\dot{\varepsilon}^k\\ b_p&=-\left(n^{k}E_{kk}+n^{u}E_{uk}\right)\dot{\varepsilon}^k \end{aligned}\end{split}\]

• time independent model with damage

  The Hooke’s law is now such that:

  \[\begin{split}\begin{aligned} \sigma&=(1-D)E\left(\varepsilon-\varepsilon^p\right)\\ \dot\sigma&=(1-D)E\left(\dot\varepsilon-\dot{p}n/(1-D)\right) -E\left(\varepsilon-\varepsilon^p\right)\dot{D} \end{aligned}\end{split}\]
  \[\begin{split}\left(\begin{array}{l} \dot{\sigma}^k\\ \dot{\sigma}^u \end{array} \right) =(1-D) \left(\begin{array}{ll} E_{kk} & E_{ku}\\ E_{uk} & E_{uu} \end{array} \right) \left(\begin{array}{l} \dot{\varepsilon}^k- \dot{p}n^k/(1-D)\\ \dot{\varepsilon}^u- \dot{p}n^u/(1-D) \end{array} \right) - \left(\begin{array}{ll} E_{kk} & E_{ku}\\ E_{uk} & E_{uu} \end{array} \right) \left(\begin{array}{l} \varepsilon^k- \varepsilon_p^k\\ \varepsilon^u- \varepsilon_p^u \end{array} \right) \dot{D} \qquad\end{split}\]
  \[\dot{\sigma}^u=(1-D)\left(E_{uk}\dot{\varepsilon}^k+E_{uu}\dot{\varepsilon}^u\right) -\dot{p}\left(E_{uk}n^k+E_{uu}n^u\right) -\dot{D}\Big( E_{uk}\left(\varepsilon^k-\varepsilon_p^k\right)+ E_{uu}\left(\varepsilon^u-\varepsilon_p^u\right) \Big)\]

  The “u” equations depend on the damage law. If, for example, the damage evolution rule \(\dot{D}=(Y/S)^s\dot{p}\) is used:

  \[\boxed{ (1-D)E_{uu}\dot{\varepsilon}^u -\Big[ \left(En\right)^u + (Y/S)^s \left(E\varepsilon^e\right)^u \Big] \dot{p} =\dot{\sigma}^u-(1-D)E_{uk}\dot{\varepsilon}^k }\]

  The damage variable leads to the modification of the consistency condition as:

  \[(1-D)n\left(E\dot{\varepsilon}\right)=\left(H+nEn\right)\dot{p}\]
  \[\begin{split}(1-D)\left(\begin{array}{ll} n^{k}\\ n^{u} \end{array} \right) \left(\begin{array}{ll} E_{kk} & E_{ku}\\ E_{uk} & E_{uu} \end{array} \right) \left(\begin{array}{l} \dot{\varepsilon}^k\\ \dot{\varepsilon}^u \end{array} \right) = \left(H+nEn\right)\dot{p}\end{split}\]
  \[\begin{split}(1-D)\left(\begin{array}{ll} n^{k}\\ n^{u} \end{array} \right) \left(\begin{array}{ll} E_{kk}\dot{\varepsilon}^k+ E_{ku}\dot{\varepsilon}^u\\ E_{uk}\dot{\varepsilon}^k+ E_{uu}\dot{\varepsilon}^u \end{array} \right) = \left(H+nEn\right)\dot{p}\end{split}\]
  \[\boxed{ (1-D)\left(n^{k}E_{ku}+n^{u}E_{uu}\right)\dot{\varepsilon}^u-\left(H+nEn\right)\dot{p}= -(1-D)\left(n^{k}E_{kk}+n^{u}E_{uk}\right)\dot{\varepsilon}^k }\]

  Finally, the system to be solved is of the form:

  \[\begin{split}\left(\begin{array}{ll} H_{cc} & H_{cp}\\ H_{pc} & H_{pp} \end{array} \right) \left(\begin{array}{l} \dot{\varepsilon}^u\\ \dot{p} \end{array} \right) = \left(\begin{array}{l} b_c\\ b_p \end{array} \right)\end{split}\]

  where:

  \[\begin{split}\begin{aligned} H_{cc}&=(1-D)E_{uu}\\ H_{pc}&=(1-D)\left(n^{k}E_{ku}+n^{u}E_{uu}\right)\\ H_{pp}&=-\left(H+nEn\right)\\ b_p&=-(1-D)\left(n^{k}E_{kk}+n^{u}E_{uk}\right)\dot{\varepsilon}^k \end{aligned}\end{split}\]

  If a damage evolution law of the form \(\dot{D}=(Y/S)^s\dot{p}\) is chosen:

  \[\begin{split}\begin{aligned} H_{cp}&= -\Big[ \left(En\right)^u + (Y/S)^s \left(E\varepsilon^e\right)^u \Big]\\ b_c&=\dot{\sigma}^u-(1-D)E_{uk}\dot{\varepsilon}^k \end{aligned}\end{split}\]

Example

Here’s an example of a time independent model involving damage

\[f = J(\ten\sigma- \ten X)-R-R_0 = 0\]

\[\ten\sigma=\tenf \Lambda:\ten{\varepsilon}^e \qquad \ten X=\frac{\displaystyle 2}{\displaystyle 3}C \ten{\alpha}\qquad R=bQr \qquad Y=\frac{\displaystyle 1}{\displaystyle 2}\ten{\varepsilon}^e:\tenf \Lambda:\ten{\varepsilon}^e\]

\[\frac{\partial f}{\partial \ten\sigma}:\overset{.}{\ten\sigma}+ \frac{\partial f}{\partial \ten X}:\overset{.}{\ten X}+ \frac{\partial f}{\partial R}:\dot R =0\]

\[\textrm{Notation} \qquad \ten n= \frac{\partial f}{\partial \ten\sigma}= \frac{\displaystyle 3}{\displaystyle 2}\ \frac{\left(\displaystyle\ten{s}-\ten X\right)}{\displaystyle J(\ten\sigma- \ten X)}\]

\[\frac{\partial f}{\partial \ten\sigma}=\ten n \qquad \frac{\partial f}{\partial \ten X}=-\ten n \qquad \frac{\partial f}{\partial R}=-1\]

\[\overset{.}{\ten\sigma}=\tenf \Lambda:\overset{.}{\overline{\varepsilon}}^e \qquad \overset{.}{\ten X}=(2/3)C\dot{\ten\alpha}\qquad \dot R=bQ\dot r \qquad \dot D=\left(\frac{Y}{S}\right)^s{\dot\lambda}_p\]

\[\dot{\ten{\varepsilon}}^p= \frac{1}{1-D}\ten n {\dot\lambda}_p \qquad \dot{\ten{\alpha}}= \frac{1}{1-D}\left(\ten n-\gamma\ten{\alpha}\right) \dot\lambda_p \qquad \dot r = \frac{1-br}{1-D}{\dot\lambda}_p\]

\[{\dot\lambda}_p=(1-D)\frac{\ten n:\dot{\ten\sigma}}{H} \qquad {\dot\lambda}_p=(1-D)\frac{\ten n:\tenf{\Lambda}:\dot{\ten{\varepsilon}}}{H+\ten n:\tenf{\Lambda}:\ten n} \qquad\]

\[H =(2/3)C\left(\ten n -\gamma\ten{\alpha}\right)\ten n +bQ(1-br)\]

@Class  TEST_CDM : BASIC_NL_BEHAVIOR, BASIC_SIMULATOR {
   @Name     test_cdm;
   @SubClass ELASTICITY elasticity;
   @tVarInt  evi,eel,alpha;
   @sVarInt  evcum,D,r;
   @Coefs    C,B,b,Q,R0;
   @Coefs    S,s;
};
@StrainPart {
  double umD = 1.0;
  if(D>0.0) {
     umD = 1.0 - Zmin(0.99,D);
  }
 sig = umD*(*elasticity*eel);
 m_tg_matrix=*elasticity;
}
@Derivative {

   double    J,Reff,Y,crit,umD,H,miso;
   TENSOR2   m,norm,sigeff,sprime,Xeff;
   ELASTICITY& E=*elasticity;
   umD = 1.0;
   if(D>0.0) {
     umD = 1.0 - Zmin(0.99,D);
   }
   sig          = umD*(E*eel);
   sigeff       = E*eel;
   Xeff         = (2./3.)*C*alpha;
   Reff         = b*Q*r;
   sprime       = deviator(sigeff-Xeff);
   J            = sqrt(1.5*(sprime|sprime));
   crit         = J - Reff - R0 ;
   if(crit<=0.){
      devcum=0.; dalpha=0.; dr=0.0; dD=0.0;
      resolve_flux_grad(*elasticity*umD, deel, deto);
}
   else {
      norm   = 1.5*sprime/J;
      m      = norm - B*alpha;
      miso   = 1. - b*r;
      H      = norm|((C/1.5)*m);
      H     += b*Q*miso;
      Y      = (1./2.)*(eel|sigeff);
      double dD_dp = pow(Y/S,s);
      resolve_flux_grad_consistency(E,norm,H,deto,devcum,eel,D,dD_dp);
      devi      = devcum*norm/umD;
      deel      = deto - devi;
      dalpha    = devcum*(norm-B*alpha)/umD;
      dr        = devcum*(1.-b*r)/umD;
      dD        = devcum*dD_dp;
  }
}