Solution scheme in DES3D

Equation of motion

DES3D solves the equation of motion, i.e., linear momentum balance equation, in the following fully dynamic form:

$$\rho \dot{\mathbf{u}}=\nabla \cdot \mathbf{\sigma }+\rho \mathbf{g}$$

where $\rho$ is the material density, $\mathbf{u}$ is the velocity vector, $\mathbf{\sigma }$ is the total (Cauchy) stress tensor, and $\mathbf{g}$ is the acceleration of gravity. The dot above $\mathbf{u}$ indicates total time derivative while bold face indicates a vector or tensor. The spatial gradient is denoted by $\nabla$, the inner product between vectors is denoted by $\cdot$, while $\nabla \cdot$ represents the divergence operator.

This equation must be complemented with appropriate initial and boundary conditions. The motion is described using a Lagrangian formulation.

Discretization

The momentum equation is discretized using a two- or three-dimensional (2D or 3D), unstructured mesh. The displacement $\mathbf{x}$, velocity $\mathbf{u}$, acceleration $\mathbf{a}$, force $\mathbf{f}$, and temperature $T$ are nodal values on linear (P1) elements while other physical quantities (e.g., stress $\mathbf{\sigma }$ and strain $\mathbf{ϵ}$) and material properties (e.g., density $\rho$ and viscosity $\eta$) are evaluated on the one-point quadrature in the elements; and thus are piecewise constant (P0).

A weak form of the momentum balance equation is constructed by the standard finite element method. The velocity function $\mathbf{u}$ is expanded with piecewise linear basis functions:

$$\mathbf{u}(\mathbf{x})=\sum _a{\mathbf{u}}_a{\mathbb{N}}_a(\mathbf{x}),$$

where ${\mathbb{N}}_a(\mathbf{x})$ is the basis function associated with global node number $a$ in the mesh with the property, ${\mathbb{N}}_i({\mathbf{x}}_j)={\delta }_{ij}$ for the coordinates of $j$-th node, ${\mathbf{x}}_j$.

${\mathbb{N}}_a(\mathbf{x})$'s are also used as weighting functions and multipled to both sides of the equation, and the product is integrated over the domain:

$${\int }_{\Omega }{\mathbb{N}}_a\rho \dot{\mathbf{u}}d\Omega ={\int }_{\Omega }{\mathbb{N}}_a\left(\nabla \cdot \mathbf{\sigma }+\rho \mathbf{g}\right)d\Omega ,a=0,\dots ,N,$$

where $N+1$ is the total number of nodes in a mesh.

Plugging in the finite dimensional representation of $\mathbf{u}$ and carrying out the element-to-global assemablage, we obtain the following equation for the acceleration ${\mathbf{a}}_a$ of every node $a$:

$$m_a{\mathbf{a}}_a={\mathbf{f}}_a={\mathbf{f}}_a^{int}+{\mathbf{f}}_a^{bc}+{\mathbf{f}}_a^{ext}.$$

$m_a$ is the nodal mass given by

$$m_a=\sum _e\left({\int }_{{\Omega }_e}{N}_e^a{\rho }_{f}d\Omega \right)=\sum _e\left({\rho }_f{\int }_{{\Omega }_e}{N}_e^{a}d\Omega \right)=\sum _e^{a\in e}\left(\frac{1}{M}{\rho }_f{\Omega }_e\right),$$

where ${\Omega }_e$ is the area (volume in 3D) of the element $e$, $N_a^e$ is $N_a$ confined to element $e$, and $M$ is the number of apexes of an element (i.e., $M=3$ for 2D triangles and $M=4$ for 3D tetrahedra.) The summation should be understood as done for all the elements having node $a$ as an apex.

☝ Two features to note
Fictitious density, ${\rho }_f$, instead of the true density, $\rho$, is used in the definition of $m_a$.
Row-sum mass lumping is applied to obtain a diagonal mass matrix.

The total force ${\mathbf{f}}_a$ is composed of three parts: the internal, boundary, and external forces. The internal force, ${\mathbf{f}}_a^{int}$, is defined as:

$$\begin{array}{rl}{\mathbf{f}}_a^{int}& =\sum _e\left({\int }_{{\Omega }_e}\frac{\partial N_a^e}{\partial \mathbf{x}}\cdot \mathbf{\sigma }d\Omega \right)=\sum _e\left(\frac{\partial N_a^e}{\partial \mathbf{x}}\cdot \mathbf{\sigma }{\int }_{{\Omega }_e}d\Omega \right)\\ & =\sum _e^{a\in e}\left(\frac{\partial N_a^e}{\partial \mathbf{x}}\cdot \mathbf{\sigma }{\Omega }_e\right).\end{array}$$

Neumann boundary conditions are tractions prescribed on the surface of the body. These tractions yield a boundary force denoted ${\mathbf{f}}_a^{bc}$:

$${\mathbf{f}}_a^{bc}=\sum _e\left({\oint }_{\partial {\Omega }_e}-N_a^e\mathbf{\sigma }\cdot \mathbf{n}dL\right)=-\sum _e^{a\in s,s\in \partial {\Omega }_e}\left(\frac{1}{M-1}{\mathbf{\sigma }}_s\cdot {\mathbf{n}}_s{L}_s\right).$$

The summation is over the boundary segment $s$, which has a length $L_s$ (surface area in 3D), the outward, unit normal vector $\mathbf{n}$, and a prescribed (constant) stress ${\mathbf{\sigma }}_s$ on the Neumann boundary. The external force, ${\mathbf{f}}_a^{ext}$, is given by:

$${\mathbf{f}}_a^{ext}=\sum _e\left({\int }_{{\Omega }_e}{N}_a^e\rho \mathbf{g}d\Omega \right)=\sum _e\left(\rho \mathbf{g}{\int }_{{\Omega }_e}{N}_a^{e}d\Omega \right)=\sum _e^{a\in e}\left(\frac{1}{M}\rho \mathbf{g}{\Omega }_e\right).$$

When deriving the equations above, we utilize the fact that ${\rho }_f$, $\rho$, $\frac{\partial N_a^e}{\partial \mathbf{x}}$, $\mathbf{\sigma }$, and $\mathbf{g}$ are constants on each element, and these identities:

$${\int }_{{\Omega }_e}d\Omega ={\Omega }_e,{\int }_{{\Omega }_e}{N}_a^{e}d\Omega =\frac{1}{M}{\Omega }_e,{\oint }_{\partial {\Omega }_e}{N}_a^{e}dL=\frac{1}{M-1}{L}_s.$$

The Rest of the solution scheme

We are interested in tectonic deformation, which can be properly simulated in a quasi-static fashion. Thus, we apply a technique called dynamic relaxation, which enables us to achieve a static equilibrium from the dynamic momentum equation by damping out the intertial force. Additionally, using mass scaling, we substitute the true density by a fictitious scaled density that allows us to increase the size of admissible stable time steps in the explicit time integration scheme. That is, using the resulting scaled acceleration and velocity, we compute an instantaneous velocity and position of each node in the mesh, which updates the model geometry at each time step.