Nodal Mixed Discretization

The linear triangular elements used in DynEarthSol2D are known to suffer volumetric locking when subject to incompressible deformations. Since incompressible plastic or viscous flow are often needed in LTM, we adopt an anti-volumetric locking correction based on the nodal mixed discretization (NMD) methodology (e.g., De Micheli and Mocellin, 2009).

The strain rate of element $e$, ${\dot{\mathbf{ϵ}}}_e$, is computed from the velocity:

$${\dot{ϵ}}_{e,ij}^{t+\Delta t}=\frac{1}{2}\sum _{a\in e}\left(\frac{\partial N_a^e}{\partial x_i}{u}_{a,j}^{t+\Delta t}+\frac{\partial N_a^e}{\partial x_j}{u}_{a,i}^{t+\Delta t}\right),$$

where $i$, $j$ are spatial indices. The strain rate tensor can be decomposed into the deviatoric and the isotropic parts:

$${\dot{\mathbf{ϵ}}}_e=\text{dev}({\dot{\mathbf{ϵ}}}_e)+\frac{1}{D}\text{tr}({\dot{\mathbf{ϵ}}}_e)\mathbf{I},$$

where $\text{dev}(\cdot )$ represents an operator returning the deviatoric tensor, $\text{tr}(\cdot )$ is an operator returning the trace of the tensor, $D$ is the number of diagonal terms of the tensor (2 for 2D case and 3 for 3D or plain strain cases), and $\mathbf{I}$ is an appropriate identity tensor. (When plane strain description is used, that is, ${ϵ}_{yy}=0$ and ${\dot{ϵ}}_{yy}=0$, but ${\sigma }_{yy}$ can be non-zero and must be included in the calculation.)

The basic idea is to average volumetric strain rate over a group of neighboring elements and then replace each element's volumetric strain rate with the averaged one. The NMD method first assigns an area (volume in 3D) average of the trace of ${\dot{\mathbf{ϵ}}}_e$ to each node $a$:

$${\dot{\epsilon }}_a=\frac{\sum _e^{a\in e}\text{tr}({\dot{\mathbf{ϵ}}}_e){\Omega }_e}{\sum _e^{a\in e}{\Omega }_e}.$$

Then the nodal field ${\dot{\epsilon }}_a$ is interpolated back to the element to retrieve an averaged volumetric strain rate for an element $e$:

$${\overline{\dot{ϵ}}}_e=\sum _{a\in e}\frac{1}{M}{\dot{\epsilon }}_a,$$

where, as before, $M$ is the number of apexes in an element.

Finally, the averaged volumetric strain rate of an element is used to modify the original strain rate tensor. The anti-locking modification replaces the isotropic part with ${\overline{\dot{ϵ}}}_e$:

$${\dot{\mathbf{ϵ}}}_e^{\prime }=\text{dev}({\dot{\mathbf{ϵ}}}_e)+\frac{1}{D}{\overline{\dot{ϵ}}}_e\mathbf{I}$$

This modified strain rate tensor substitutes the original strain rate tensor when updating strain tensor and in defining constitutive update. For the sake of brevity, we drop the prime and use $\dot{\mathbf{ϵ}}$ to refer the modified strain rate tensor from now on.

The strain tensor $\mathbf{ϵ}$ is accumulated:

$${\mathbf{ϵ}}_e^{t+\Delta t}={\mathbf{ϵ}}_e^t+\Delta t{\dot{\mathbf{ϵ}}}_e^{t+\Delta t}$$