BASIC_SIMULATOR

Description

The BASIC_SIMULATOR ZebFront class type is used for FEM behaviors which are meant to be run in the simulation mode as well. Classes of this type must be valid BASIC_NL_BEHAVIOR classes, and require special code to calculate the mixed-mode loading case, but also allow definition of yield surface functions.

Syntax

The class declaration is the same as that for BASIC_NL_BEHAVIOR with some extensions. No new class methods (main functions) are provided for the simulation mode.

The class definition for BASIC_SIMULATOR is the following:

@class NAME_OF_CLASS : BASIC_NL_BEHAVIOR, BASIC_SIMULATOR { \(~\,~\,~\,\) standard behavior options from page \(~\,~\,~\,\) @Criterion list; \(~\,~\,~\,\) additional C++ code };

Resolving mixed flux-grad loading

One of the main benefits of the simulator is the ability to solve mixed loadings exactly, while displacement based FE solutions are approximate and require iterations in the non-linear case (and therefore desiring a good tangent matrix calculation). The disadvantage is a formulation must be made to solve the mixed rate loading. To address the later, special methods are given in the BASIC_SIMULATOR base class for resolving such a difficulty (this is in contrast to the solution in the class SIMUL_MODEL). Using the notation \({\bf f}\) for the flux and \(\bf g\) for the gradient, the following forms are allowed. See the developers manual for a description of how these methods are implemented.

Code	Description
resolve_flux_grad(E, de, dg, de2)	\(\dot{{\bf f}}=\bf E(\dot{\bf g}-\dot{\bf e}-\dot{\bf e_2})\)
resolve_flux_grad(E, de, dg)	\(\dot{{\bf f}}=\bf E(\dot{\bf g}-\dot{\bf e})\)

Example

Here’s an example of a combined FEM-simulator model with a criterion object.

@Class FEM_SIM_BEHAVIOR : BASIC_NL_BEHAVIOR, BASIC_SIMULATOR {
   @Name example;
   @Coefs    E, poisson;
   @Coefs    R0, H, Q, b;
   @Coefs    alpha, beta, A, k, r;
   @Coefs    dmax;
   @tVarInt  eel;
   @sVarInt  evcum, D;
   @sVarAux  R,j0,j1,j2,chi;
   @tVarAux  evi;
   @Criterion yield, damage;
};

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;
  }
}