A 117
Appendix Algorithm of SNAC
A.1 Governing equations
The software package SNAC solves the momentum and the heat energy balance equations in the following dierential form: Dvi ij , + gi = xj Dt qi DT . + r = Cp xi Dt (A.1)
(A.2)
In the momentum balance equation, is the mass density, vi is velocity, ij is the Cauchy stress tensor, and gi is the gravitational acceleration. T is temperature, Cp is the specic heat at constant pressure, qi is the heat ux vector, and r is the volumetric heat source. D/Dt represents the material time derivative. In this study, no heat sources are considered, including shear heating. Viscosity is temperature- and/or stress-dependent. The elastic component of stress has an extra contribution from thermal stress.
A.2
Spatial discretization
A 3-D domain is discretized into hexahedral elements, each of which is lled with two sets of 5 tetrahedra (Fig. A.1a). In this mesh hierarchy, called the mixed discretization (Marti and Cundall, 1982), hexahedral elements are used only as an averaging
118
Figure A.1: (a) Two congurations of ve tetrahedra in a hexahedral element used in the mixed discretization. Numbers next to apexes indicate the local node numbering. (b) Conventions for the notation. Al and nl denote the face and the unit normal vector, respectively, associated with a local node l. unit for volumetric strain. The averaging is enforced at all times, for incompressible viscoelastic or plastic constitutive laws. The use of two equivalent sets of tetrahedra is required to ensure a symmetric response. For a given loading, responses of one set of tetrahedra can be dierent from those of the other set because of the dierently orientated faces of tetrahedra in each set (e.g., Zienkiewicz et al., 1995). The approximation of partial derivatives with respect to spatial variables follows the integral denitions (e.g., Wilkins, 1964): f,i dV =
f ni d,
(A.3)
where represent a tetrahedron as an integration domain, is the boundary surfaces of the tetrahedron, f,i is the partial derivative of a variable f with respect to i-th spatial coordinate, ni is the i-th component of the unit normal vector of the surface.
119 If the partial derivative is constant within a tetrahedron, it is evaluated as f,i = 1 V f ni d,
(A.4)
where V is the volume of the tetrahedron. By further substituting an algebraic expression for the surface integral, reordering terms, and using f = 1 in (A.4)), 1 f,i = V = 1 3V
4
ni d = 0 (when
l=1 4
1 f l nl Al = i V
4
4
l=1
1 f m nl Al i 3 m=1,=l (A.5)
4
fm
m=1 4 l=1,=m
nl Al i
=
1 3V
f m nm Am , i
m=1
where l is the local node index varying from 1 to 4, Al and nl are the area and the unit normal vector of the triangular surface not having the node l as one of its apexes (Fig. A.1b). Hereafter, we call such a face a corresponding face to node l. f l is the averaged f on the surface l.
A.3
Nodal assemblage
We can convert the dierential equation for momentum balance (A.1) (the following description is applied to the heat equation in the same fashion) to a principle of minimum work rate as in the standard nite element formulation: vi
Dvi dV = Dt
vi gi dV +
ij ij dV,
(A.6)
where ij are components of the strain rate tensor, vi and ij represent variations of velocity and strain rate, and here corresponds to the whole domain. The local contribution to nodes corresponding to each term can be computed by following the standard nite element procedure for linear tetrahedral elements. However, our
120 method does not need to construct coecient matrices such as mass and stiness matrices since it adopts an explicit time discretization. The resultant momentum equation is Mn
n 1 [n] 1 Dvi = Ti + [n] gi V [n] , Dt 3 4
(A.7)
where the superscript n represents values evaluated at the global node n, the superscript [n] means the sum of contributions from all the tetrahedra having the global node n as an apex, Ti is the traction that is dened as ij nj and evaluated on a face of one of the contributing tetrahedra. The nodal mass M n is not the actual inertial mass but an adjusted one to satisfy a local stability criterion discussed in the section A.5. The correspondence between an apex and a face for the traction calculation is determined as in the derivation of the expression, (A.5). Note that the factor of 1/3 in the traction term is inherited from (A.5) and the factor of 1/4 in the body force term implies that the nodal contribution takes one quarter of a tetrahedrons volume-dependent quantity. While looping over the entire set of nodes, mass and nodal forces are assembled by adding up the contributions from boundary conditions and all the tetrahedra sharing that node as one of their apexes. The structured mesh of SNAC renders the assemblage step conveniently static. The acquired net force (or heat ux) at each node is used to update velocities and node coordinates (or temperature).
A.4
Damping and explicit time marching
We seek static or quasi-static solutions through a dynamic relaxation method. Instead of adding a usual velocity-dependent friction term, we adopt a local non-viscous damping scheme (Cundall, 1987): Fidamped = Fi sgn(vi )|Fi |, (A.8)
where Fi is the i-th component of the residual force vector, is a positive coecient less than 1, sgn(vi ) returns the sign of the i-th component of velocity, vi . Once net
121 forces are assembled and damped, velocity at that node is updated using a forward Euler method: t t Fidamped v(t + ) = v(t ) + t 2 2 M x(t + t) = x(t) + tv(t + t ). (A.10)
Damping 2
(A.9) is irrelevant to the update of temperature eld, but the same forward Euler method is used.
A.5
Mass scaling for numerical stability
The conventional CourantFriedrichsLewy (CFL) condition imposes a stringent upper limit for the time step size such that dynamic relaxation takes long time to get quasistatic solution over a geological time scale. To overcome this limit, a mass scaling technique is applied. This technique adjusts each nodal mass such that the stability condition for a user-specied time step can be locally satised. The stability condition to be satised, however, is not the same as in the CFL condition, i.e., t (lmin /vp ), where t is the time step, lmin is the minimum element size, and vp is the P wave velocity. Instead, through an analogy of continuum to an innite mass-spring system, we use a criterion that does not explicitly include length scale and P wave velocity (see Ch. 9 in Bathe, 1996): t T , (A.11)
where T is the period of system, 2(m/k)1/2 , m is a point mass, and k is the stiness of the spring attached to the point mass. Now, reducing the innite series of mass and springs in one dimension to a single mass-spring system, the stiness of that single system becomes 4k, leading to an expression for the mass scaling: m k(t)2 . (A.12)
For a given size of t, the nodal mass is adjusted according to (A.12) to automatically satisfy the stability critetion, (A.11). The value of k is computed by equating internal
122 force contribution at a node with kui : 1 Ti = kui 3
1 ( + 2)(ii dt)ni S = k(vi dt) (no sum), 3
(A.13)
where only the volumetric contribution from internal forces is taken into account. By substituting the approximation for the partial derivative (A.5) into the above equation and dividing both sides by vi dt, we obtain
l ki =
1 ( + 2)(nl S l )2 , i 9V
(A.14)
where l is the local index for apexes of a tetrahedron and the surface-related quantities are computed on the corresponding face of the tetrahedron. Finally, a tetrahedrons contribution to the scaled mass is given as ml = + 2 max[(nl S l )2 , i = 1, . . . , 3]. i 9V (A.15)
As in the standard FEM, appropriate mappings between local and global indices are required.
A.6
Constitutive update
SNAC uses a general elasto-visco-plastic rheological model to update the Cauchy stress tensor (e.g., Albert et al., 2000). First, the initial guess of stress is acquired by the Maxwell viscoelastic constitutive law (Poliakov et al., 1993). If this initial guess exceeds a given yield stress, it is projected onto the yield surface using a return mapping method (Simo and Hughes, 2004); otherwise, the viscoelastic stress update is retained. This elasto-visco-plastic model can deal with various constitutive laws that are typically used for the Earths crustal and mantle material as its limiting cases. For example, elastic, viscoelastic and elastoplastic material are realized in the following cases:
123 1. Elastic material corresponds to the limit of innite viscosity and yield stress. 2. Viscoelastic material corresponds to the limit of innite yield strength. 3. Elasoplastic material corresponds to the ininte viscosity. Using the viscoplastic rheology is physically more realistic than using one of the limiting cases listed above since all materials have dissipative mechani...