Constitutive update

The stress tensor is updated using the strain rate and strain tensors according to an appropriate constitutive relationship. Since the stress update calculations are performed at the element level, we drop the subscript $e$ to simplify notation. The EVP material model is approximated by a composite rheology which uses visco-elastic and elasto-plastic sub-models. With the bulk modulus $K_s$, shear modulus $G$, viscosity $\eta$, cohesion $C$, and internal friction angle $\varphi$, we calculate the visco-elastic stress ${\mathbf{\sigma }}_{ve}$ and the elasto-plastic stress ${\mathbf{\sigma }}_{ep}$.

The visco-elastic stress increment $\Delta {\mathbf{\sigma }}_{ve}$ is calculated assuming a linear Maxwell material, where a total deviatoric strain increment $\Delta \mathbf{ϵ}$ is composed of the elastic and the viscous components while the deviatoric stress increment is identical for each component:

$$\text{dev}(\Delta \mathbf{ϵ})=\frac{\text{dev}(\Delta {\mathbf{\sigma }}_{ve})}{2G}+\frac{\text{dev}({\mathbf{\sigma }}_{ve})\Delta t}{2\eta }$$

Substituting $\Delta \mathbf{ϵ}$ with ${\mathbf{ϵ}}^{t+\Delta t}-{\mathbf{ϵ}}^t,$ $\Delta {\mathbf{\sigma }}_{ve}$ with ${\mathbf{\sigma }}_{ve}^{t+\Delta t}-{\mathbf{\sigma }}^t$, and ${\mathbf{\sigma }}_{ve}$ with $({\mathbf{\sigma }}_{ve}^{t+\Delta t}+{\mathbf{\sigma }}^t)/2$, the equation above is reduced to:

$$\text{dev}({\mathbf{\sigma }}_{ve}^{t+\Delta t})=\frac{\left(1-\frac{G\Delta t}{2\eta }\right)\text{dev}({\mathbf{\sigma }}^t)+2G\text{dev}({\mathbf{ϵ}}^{t+\Delta t}-{\mathbf{ϵ}}^t)}{1+\frac{G\Delta t}{2\eta }}$$

The isotropic stress components are updated based on the volume change. As a result, the visco-elastic stress is:

$${\mathbf{\sigma }}_{ve}^{t+\Delta t}=\text{dev}({\mathbf{\sigma }}_{ve}^{t+\Delta t})+\Delta t{K}_s\text{tr}({\dot{\mathbf{ϵ}}}^{t+\Delta t})\mathbf{I}.$$

The elasto-plastic stress ${\mathbf{\sigma }}_{ep}$ is computed using linear elasticity and the Mohr-Coulomb (MC) failure criterion with a general (associative or non-associative) flow rule. Following a standard operator-splitting scheme [e.g., @Lubliner1990; @SimoHugh2004; @Wilkins1964a], an elastic trial stress ${\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}$ is first calculated as

$${\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}={\mathbf{\sigma }}^t+(K_s-\frac{2}{3}G)\text{tr}({\dot{\mathbf{ϵ}}}^{t+\Delta t})\mathbf{I}\Delta t+2G{\dot{\mathbf{ϵ}}}^{t+\Delta t}\Delta t.$$

If the elastic trial stress, ${\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}$, is on or within a yield surface, that is, $f\left({\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}\right)\ge 0,$ where $f$ is the yield function, then the stress does not need a plastic correction. So, ${\mathbf{\sigma }}_{ep}^{t+\Delta t}$ is set to be equal to ${\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}$. However, if ${\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}$ is outside the yield surface, we project it onto the yield surface using a return-mapping algorithm [@SimoHugh2004].

In the case of a Mohr-Coulomb material, it is convenient to express the yield function for shear failure in terms of principal stresses:

$$f_s({\sigma }_1,{\sigma }_3)={\sigma }_1-N_{\varphi }{\sigma }_3+2C\sqrt{N_{\varphi }},$$

where ${\sigma }_1$ and ${\sigma }_3$ are the maximal and minimal compressive principal stresses with the sign convention that tension is positive (i.e., ${\sigma }_1\le {\sigma }_2\le {\sigma }_3$), $C$ is the material's cohesion, $N_{\varphi }=\frac{1+\sin \varphi }{1-\sin \varphi }$, $\sqrt{N_{\varphi }}=\frac{\cos \varphi }{1-\sin \varphi }$, and $\varphi$ is an internal friction angle ($<90^{\circ }$). The yield function for tensile failure is defined as

$$f_t({\sigma }_3)={\sigma }_3-{\sigma }_t,$$

where ${\sigma }_t$ is the tension cut-off. If a value for the tension cut-off is given as a parameter, the smallest value between the theoretical limit ($C/\tan \varphi$) and the given value is assigned to ${\sigma }_t$. This comparison is required because the theoretical limit is not constant in the strain weakening case, where the material cohesion, $C$, and the friction angle $\varphi$ may change.

To guarantee a unique decision on the mode of yielding (shear versus tensile), we define an additional function, $f_h({\sigma }_1,{\sigma }_3)$, which bisects the obtuse angle made by two yield functions on the ${\sigma }_1$-${\sigma }_3$ plane, as

$$\begin{array}{rl}{f}_h({\sigma }_1,{\sigma }_3)& ={\sigma }_3-{\sigma }_t+\left(\sqrt{N_{\varphi }^2+1}+N_{\varphi }\right)\left({\sigma }_1-N_{\varphi }{\sigma }_t+2C\sqrt{N_{\varphi }}\right).\end{array}$$

Once yielding is identified, that is, $f_s({\sigma }_{el,1},{\sigma }_{el,3})<0$ or $f_t({\sigma }_{el,3})>0$, the mode of failure (shear or tensile) is decided based on the value of $f_h$, in other words, shear failure occurs if $f_h({\sigma }_{el,1},{\sigma }_{el,3})<0$, tensile failure occurs otherwise.

The flow rule for frictional materials is in general non-associative, that is, the direction of plastic flow in the principal stress space during plastic flow is not the same as the direction of the vector normal to the yield surface. As in the definitions of yield functions, the plastic flow potential for shear failure in the Mohr-Coulomb model can be defined as

$$g_s\left({\sigma }_1,{\sigma }_3\right)={\sigma }_1-\frac{1+\sin \psi }{1-\sin \psi }{\sigma }_3,$$

where $\psi$ is the dilation angle. Likewise, the tensile flow potential is given as

$$g_t\left({\sigma }_3\right)={\sigma }_3-{\sigma }_t.$$

In the presence of plasticity, the total strain increment $\Delta \mathbf{ϵ}$ is given by

$$\Delta \mathbf{ϵ}=\Delta {\mathbf{ϵ}}_{\text{el}}+\Delta {\mathbf{ϵ}}_{\text{pl}},$$

where $\Delta {\mathbf{ϵ}}_{\text{el}}$ and $\Delta {\mathbf{ϵ}}_{\text{pl}}$ are the elastic and plastic strain increments, respectively. The plastic strain increment is normal to the flow potential surface and can be written as

$$\Delta {\mathbf{ϵ}}_{\text{pl}}=\beta \frac{\partial g}{\partial \mathbf{\sigma }},$$

where $\beta$ is the plastic flow magnitude. $\beta$ is computed by requiring that the updated stress state lies on the yield surface,

$$f\left({\mathbf{\sigma }}_{ep}^{t+\Delta t}\right)=f\left({\mathbf{\sigma }}^t+\Delta {\mathbf{\sigma }}_{ep}\right)=0.$$

In the principal component representation, ${\sigma }_A=E_{AB}{ϵ}_B$ where ${\sigma }_A$ and ${ϵ}_A$ are the principal stress and strain, respectively, and $\mathbf{E}$ is a corresponding elastic moduli matrix with components:

$$\begin{array}{rlrl}{E}_{AB}& =\left(K_s-\frac{2}{3}G\right)& & \text{if }A\ne B,\\ E_{AB}& =\left(K_s+\frac{4}{3}G\right)& & \text{otherwise.}\end{array}$$

By applying the consistency condition and using ${\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}={\mathbf{\sigma }}^t+\mathbf{E}\cdot \Delta \mathbf{ϵ}$ (in the principal component representation), we obtain the following formulae for $\beta$

$$\beta =\frac{{\sigma }_{\text{el},1}^{t+\Delta t}-N_{\varphi }{\sigma }_{\text{el},3}^{t+\Delta t}+2C\sqrt{N_{\varphi }}}{{\sum }_B\left(E_{1B}\frac{\partial g_s}{\partial {\sigma }_B}-N_{\varphi }{E}_{3B}\frac{\partial g_s}{\partial {\sigma }_B}\right)}\text{for shear failure,}$$

and

$$\beta =\frac{{\sigma }_{\text{el},3}^{t+\Delta t}-{\sigma }_t}{\frac{\partial g_t}{\partial {\sigma }_3}}\text{for tensile failure.}$$

Likewise, $\partial g/\partial \mathbf{\sigma }$ takes different forms according to the failure mode:

$$\begin{array}{rl}\partial g/\partial {\sigma }_1& =1\\ \partial g/\partial {\sigma }_2& =0\\ \partial g/\partial {\sigma }_3& =-\frac{1+\sin \psi }{1-\sin \psi }\end{array}\text{(for shear failure,)}$$

and

$$\begin{array}{rl}\partial g/\partial {\sigma }_1& =0\\ \partial g/\partial {\sigma }_2& =0\\ \partial g/\partial {\sigma }_3& =1\end{array}\text{(for tensile failure.)}$$

Once $\Delta {\mathbf{ϵ}}_{pl}$ is computed, ${\mathbf{\sigma }}_{ep}$ is updated as

$${\mathbf{\sigma }}_{ep}={\mathbf{\sigma }}_{\text{el}}^{t+\Delta t}-\mathbf{E}\cdot \Delta {\mathbf{ϵ}}_{\text{pl}}.$$

in the principal component representation and transformed back to the orignal coordinate system.

After the visco-elastic stress ${\mathbf{\sigma }}_{ve}$ and elasto-plastic stress ${\mathbf{\sigma }}_{ep}$ are evaluated, we compute the second invariant of the deviatoric components of each. If the visco-elastic stress has a smaller second invariant ($J_2$), ${\mathbf{\sigma }}_{ve}$ is be used as the updated stress; otherwise, ${\mathbf{\sigma }}_{ep}$ is used.

The fundamental deformation measures in DynEarthSol are strain rates. Thus, the stress update by rate-independent constitutive models like elasto-plastic stresses need to be considered as the time-integration of the rate form of the corresponding stresses. Since a stress rate is not frame-indifferent in general, an objective (or co-rotational) stress rate needs to be constructed and integrated instead. The Jaumann stress rate is our choice for DynEarthSol among the possible objective rates because of its simplicity.

The Jaumann stress rate ($\check{\dot{\sigma }}$) is defined as

$$\check{\dot{\mathbf{\sigma }}}=\dot{\mathbf{\sigma }}-\mathbf{\omega }\cdot \mathbf{\sigma }+\mathbf{\sigma }\cdot \mathbf{\omega },$$

where $\mathbf{\omega }$ is the spin tensor, which is defined as,

$${\omega }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right).$$

Based on this definition, the new objective stress (${\check{\mathbf{\sigma }}}^{t+\Delta t}$) is,

$${\check{\mathbf{\sigma }}}^{t+\Delta t}={\mathbf{\sigma }}^{t+\Delta t}+\Delta t({\mathbf{\sigma }}^{t+\Delta t}\cdot {\mathbf{w}}^{t+\Delta t}-{\mathbf{w}}^{t+\Delta t}\cdot {\mathbf{\sigma }}^{t+\Delta t}),$$

where ${\mathbf{\sigma }}^{t+\Delta t}$ is the updated stress equal to either ${\mathbf{\sigma }}_{ve}$ or ${\mathbf{\sigma }}_{ep}$, depending on which has a lower value of $J_2$.