[cigcommits] r13139  in doc/snac: . figures
echoi at geodynamics.org
echoi at geodynamics.org
Tue Oct 28 10:38:45 PDT 2008
Author: echoi
Date: 20081028 10:38:45 0700 (Tue, 28 Oct 2008)
New Revision: 13139
Added:
doc/snac/figures/snacmesh.pdf
Modified:
doc/snac/snac.lyx
Log:
Updated the section 1.1 SNAC Implementation.
Added: doc/snac/figures/snacmesh.pdf
===================================================================
(Binary files differ)
Property changes on: doc/snac/figures/snacmesh.pdf
___________________________________________________________________
Name: svn:mimetype
+ application/octetstream
Modified: doc/snac/snac.lyx
===================================================================
 doc/snac/snac.lyx 20081028 16:14:29 UTC (rev 13138)
+++ doc/snac/snac.lyx 20081028 17:38:45 UTC (rev 13139)
@@ 1,4 +1,4 @@
#LyX 1.5.1 created this file. For more info see http://www.lyx.org/
+#LyX 1.5.3 created this file. For more info see http://www.lyx.org/
\lyxformat 276
\begin_document
\begin_header
@@ 259,13 +259,17 @@
\begin_layout Itemize
E.
 Choi, Thoutireddy, P., Lavier, L., Quenette, S., Tan, E., Gurnis, M., Aivazis
 M., and Appelbe, B., GeoFramework Part II: Coupling models of crustal deformation
 and mantle convection with a computational framework, to be submitted to

+ Choi, Lavier, L., Gurnis, M., Thermomechanics of midocean ridge segmentation,
+ accepted by Physics of the Earth and Planetary Interiors, 2008,
\emph on
Geochem., Geophys., Geosys.
+\begin_inset LatexCommand htmlurl
+name "http://dx.doi.org/10.1016/j.pepi.2008.08.010"
+target "http://dx.doi.org/10.1016/j.pepi.2008.08.010"
+
+\end_inset
+
+.
\end_layout
\begin_layout Standard
@@ 319,7 +323,7 @@
\begin_layout Standard
SNAC (StGermaiN Analysis of Continua) is an updated Lagrangian explicit
 finite element code for modeling a finitely deforming viscoplastic solid
+ finite difference code for modeling a finitely deforming viscoplastic solid
in 3D, released under the GNU General Public License (see Appendix
\begin_inset LatexCommand vref
reference "cha:License"
@@ 329,6 +333,9 @@
).
In this code, nodal velocities satisfying a weakform of the momentum balance
are obtained as the nodal solution.
+ SNAC shares a mathematical foundation, and thus major advantages, with
+ a standard finite element method (FEM).
+ However, it departs from the FEM by not making explicit use of shape functions.
A Cartesian mesh consisting of 4node linear or constantstrain tetrahedral
elements is used to represent a discretized domain, although a spherical
domain can also be used.
@@ 348,7 +355,7 @@
\end_layout
\begin_layout Subsection
Weak Formulation
+Governing Equations
\end_layout
\begin_layout Standard
@@ 359,16 +366,12 @@
\begin_layout Standard
\begin_inset Formula \begin{equation}
\begin{array}{c}
\frac{\partial\sigma_{ij}}{\partial x_{j}}+\rho g_{i}=\frac{Dv_{i}}{Dt}\,\, in\,\,\Omega\\
v_{i}=g_{i}\,\, on\,\,\Gamma_{g}\\
t_{i}=h_{i}\,\, on\,\,\Gamma_{h}\end{array}\label{eq:Momentum balance for continuum}\end{equation}
+\frac{\partial\sigma_{ij}}{\partial x_{j}}+\rho g_{i}=\frac{Dv_{i}}{Dt}\,\, in\,\,\Omega,\\
+v_{i}=g_{i}\,\, on\,\,\Gamma_{g},\\
+t_{i}=h_{i}\,\, on\,\,\Gamma_{h},\end{array}\label{eq:Momentum balance for continuum}\end{equation}
\end_inset

\end_layout

\begin_layout Standard
where
\begin_inset Formula $\partial\Omega$
\end_inset
@@ 401,329 +404,734 @@
is the material or total derivative, and equal to partial derivative with
respect to time in SNAC since the updated Lagrangian viewpoint is taken.
 The corresponding weak form of the momentum balance is
+
\end_layout
+\begin_layout Subsection
+Spatial Descritization
+\end_layout
+
\begin_layout Standard
\begin_inset Formula \begin{equation}
\int\delta\nu_{i,j}\sigma_{ij}d\Omega\int\delta\nu_{i}\rho g_{i}d\Omega+\int\delta\nu_{i}\rho\frac{Dv_{i}}{Dt}d\Omega\int_{\Gamma_{h}}\delta\nu_{i}\left(t_{i}h_{i}\right)d\Gamma=0\label{eq:weak form of momentum balance}\end{equation}
+A 3D domain is discretized into hexahedral elements, each of which is filled
+ with two sets of 5 tetrahedra (Fig.
+
+\begin_inset LatexCommand ref
+reference "fig:snac_mesh"
\end_inset
+a).
+ In this mesh hierarchy, called the mixed discretization
+\begin_inset ERT
+status collapsed
+\begin_layout Standard
+
+
+\backslash
+citep{MartCund1982}
\end_layout
+\end_inset
+
+, hexahedral elements are used only as an averaging 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 different
+ from those of the other set because of the differently orientated faces
+ of tetrahedra in each set
+\begin_inset ERT
+status collapsed
+
\begin_layout Standard
where
\begin_inset Formula $\delta\nu_{i}$
+
+
+\backslash
+citep[e.g.,]{Zien_etal1995}
+\end_layout
+
\end_inset
 are components of virtual velocity and a comma stands for partial differentiati
on.
 This weak form is also known as principle of virtual power.
 Introducing the isoparametric discretization for the velocity in each element
+.
\end_layout
\begin_layout Standard
\begin_inset Formula \begin{equation}
v_{i}=v_{i}^{a}N_{a}\left(X\right)\label{eq:isoparametric discretization for the velocity}\end{equation}
+\begin_inset Float figure
+wide false
+sideways false
+status open
+\begin_layout Standard
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+centering
+\end_layout
+
\end_inset
+
+\begin_inset Graphics
+ filename figures/snacmesh.pdf
+ scale 60
+\end_inset
+
+
\end_layout
\begin_layout Standard
where
\begin_inset Formula $v_{i}^{a}$
+\begin_inset Caption
+
+\begin_layout Standard
+\begin_inset OptArg
+status open
+
+\begin_layout Standard
+Configurations of tetrahedra and conventions for the notation
+\end_layout
+
\end_inset
 is the
\emph on
i
\emph default
th component velocity of node
\begin_inset Formula $a$
+(a) Two configurations of five tetrahedra in a hexahedral element used in
+ the mixed discretization.
+ Numbers next to apexes indicate the local node numbering.
+ (b) Conventions for the notation.
+
+\begin_inset Formula $A_{l}$
\end_inset
,
\begin_inset Formula $N_{a}$
+ and
+\begin_inset Formula $n_{l}$
\end_inset
 is the shape function associated with node
\begin_inset Formula $a$
+ denote the face and the unit normal vector, respectively, associated with
+ a local node
+\begin_inset Formula $l$
\end_inset
, and
\begin_inset Formula $X$
+.
+\end_layout
+
\end_inset
 material coordinate.
 Since the discretization is isoparametric, the same discretization holds
 for virtual velocities.
 This discretization then provides a discrete version of the principle of
 virtual power:
+
\end_layout
\begin_layout Standard
\begin_inset Formula \begin{equation}
\delta\nu_{i}^{a}\left(\int\sigma_{ij}N_{a,j}d\Omega\int N_{a}\rho g_{i}d\Omega+\int N_{a}\rho\frac{Dv_{i}}{Dt}d\Omega\int_{\Gamma_{h}}N_{a}\left(t_{i}h_{i}\right)d\Gamma\right)=0\label{eq:discrete version of virtual power}\end{equation}
+\begin_inset LatexCommand label
+name "fig:snac_mesh"
\end_inset
+
+\end_layout
+\end_inset
+
+
\end_layout
+\begin_layout Subsection
+Approximation of partial derivatives
+\end_layout
+
\begin_layout Standard
Since the virtual velocities,
\begin_inset Formula $\delta\nu_{i}$
\end_inset
+The approximation of partial derivatives with respect to spatial variables
+ follows the integral definitions
+\begin_inset ERT
+status collapsed
, are arbitrary everywhere except on Dirichlet boundary
+\begin_layout Standard
+
+
+\backslash
+citep[e.g.,]{Wilkins1964}
\end_layout
\begin_layout Standard
+\end_inset
+
+:
\begin_inset Formula \begin{equation}
r_{i}^{a}=\int N_{a}\rho g_{i}d\Omega\int\sigma_{ij}N_{a,j}d\Omega\int N_{a}\rho\frac{Dv_{i}}{Dt}d\Omega\int_{\Gamma_{h}}N_{a}\left(t_{i}h_{i}\right)d\Gamma=0\label{eq:5}\end{equation}
+\int_{\Omega}f_{,i}dV=\int_{\partial\Omega}fn_{i}d\Gamma,\label{eq:partialapprox}\end{equation}
\end_inset
+where
+\shape italic
\end_layout
+\begin_inset Formula $\Omega$
+\end_inset
\begin_layout Standard
+
+\shape default
+ represent a tetrahedron as an integration domain,
+\shape italic
+
+\begin_inset Formula $\partial\Omega$
+\end_inset
+
+
+\shape default
+ is the boundary surfaces of the tetrahedron,
+\begin_inset Formula $f_{,i}$
+\end_inset
+
+ is the partial derivative of a variable
+\begin_inset Formula $f$
+\end_inset
+
+ with respect to
+\begin_inset Formula $i$
+\end_inset
+
+th spatial coordinate,
+\begin_inset Formula $n_{i}$
+\end_inset
+
+ is the
+\begin_inset Formula $i$
+\end_inset
+
+th component of the unit normal vector of the surface.
+ If the partial derivative is constant within a tetrahedron, it is evaluated
+ as
+\begin_inset Formula \begin{equation}
+f_{,i}=\frac{1}{V}\int_{\partial\Omega}fn_{i}d\Gamma,\label{eq:def_fi}\end{equation}
+
+\end_inset
+
where
\begin_inset Formula $r_{i}^{a}$
+\begin_inset Formula $V$
\end_inset
 is the residual force vector.
 The above set of equations corresponds to discrete momentum balance and
 we solve them for the nodal velocities,
\begin_inset Formula $v_{i}^{a}$
+ is the volume of the tetrahedron.
+ By further substituting an algebraic expression for the surface integral,
+ reordering terms, and using
+\begin_inset Formula $\int_{\partial\Omega}n_{i}d\Gamma=0$
\end_inset
.
 Because SNAC seeks a quasistatic solution, the inertia term in Eq.

+ (when
+\begin_inset Formula $f=1$
+\end_inset
+
+ in (
\begin_inset LatexCommand ref
reference "eq:5"
+reference "eq:def_fi"
\end_inset
 can be dropped.
 This weak form of momentum balance equation, Eq.
+)),
+\begin_inset Formula \begin{equation}
+\begin{split}f_{,i} & =\frac{1}{V}\sum_{l=1}^{4}\bar{f}^{l}n_{i}^{l}A^{l}=\frac{1}{V}\sum_{l=1}^{4}\frac{1}{3}\sum_{m=1,\neq l}^{4}f^{m}n_{i}^{l}A^{l}\\
+ & =\frac{1}{3V}\sum_{m=1}^{4}f^{m}\sum_{l=1,\neq m}^{4}n_{i}^{l}A^{l}\\
+ & =\frac{1}{3V}\sum_{m=1}^{4}f^{m}n_{i}^{m}A^{m},\end{split}
+\label{eq:formula_fi}\end{equation}
+
+\end_inset
+
+ where
+\begin_inset Formula $l$
+\end_inset
+
+ is the local node index varying from 1 to 4,
+\begin_inset Formula $A^{l}$
+\end_inset
+
+ and
+\begin_inset Formula $n^{l}$
+\end_inset
+
+ are the area and the unit normal vector of the triangular surface not having
+ the node
+\begin_inset Formula $l$
+\end_inset
+
+ as one of its apexes (Fig.
\begin_inset LatexCommand ref
reference "eq:5"
+reference "fig:snac_mesh"
\end_inset
, is equivalent to that acquired from the original Finite Difference formulation
 of FLAC [Cundall, 1989].
 In the finite difference formulation, all the partial derivatives are replaced
 by the integral definitions [Wilkins, 1964], so that the partial derivatives
 of shape functions are not involved explicitly.
+b).
+ Hereafter, we call such a face a
+\emph on
+corresponding
+\emph default
+ face to node
+\begin_inset Formula $l$
+\end_inset
+
+.
+\begin_inset Formula $\bar{f}^{l}$
+\end_inset
+
+ is the averaged
+\begin_inset Formula $f$
+\end_inset
+
+ on the surface
+\begin_inset Formula $l$
+\end_inset
+
+.
\end_layout
\begin_layout Subsection
Solution Scheme
+Nodal Assemblage
\end_layout
\begin_layout Standard
In SNAC, the solution to velocity is obtained by an explicit version of
 dynamic relaxation.
 Consider the autonomous set of equations
\end_layout
+We can convert the differential equation for momentum balance (
+\begin_inset LatexCommand ref
+reference "eq:Momentum balance for continuum"
\begin_layout Standard
+\end_inset
+
+) (the following description is applied to the heat equation in the same
+ fashion) to a principle of minimum work rate as in the standard finite
+ element formulation:
\begin_inset Formula \begin{equation}
M\frac{D^{2}\nu}{Dt^{2}}+C\frac{Dv}{Dt}+K\nu=f\label{eq:autonomous set of equations}\end{equation}
+\int_{\Omega}\delta v_{i}\rho\frac{Dv_{i}}{Dt}dV=\int_{\Omega}\delta v_{i}\rho g_{i}dV+\int_{\Omega}\delta\xi_{ij}\sigma_{ij}dV,\label{eq:momentumweak}\end{equation}
\end_inset
+where
+\shape italic
\end_layout
+\begin_inset Formula $\xi_{ij}$
+\end_inset
\begin_layout Standard
where
\series bold
\emph on
M
\series default
\emph default
,
\series bold
\emph on
C
\series default
\emph default
, and
\series bold
\emph on
K
\series default
\emph default
 are mass, damping and stiffness matrix respectively and
\series bold
\emph on
\begin_inset Formula $f$
+\shape default
+ are components of the strain rate tensor,
+\shape italic
+
+\begin_inset Formula $\delta v_{i}$
\end_inset
\series default
\emph default
 is independent of time.
 For
\series bold
\emph on
C
\series default
\emph default
, positive definite, the transient velocity solution converges to that of
 equilibrium velocity solution
\begin_inset Formula $K\nu=f$
+\shape default
+ and
+\shape italic
+
+\begin_inset Formula $\delta\xi_{ij}$
\end_inset
 as
\begin_inset Formula $t\rightarrow\infty$
+
+\shape default
+ represent variations of velocity and strain rate, and
+\shape italic
+
+\begin_inset Formula $\Omega$
\end_inset
.
 Critical damping is desired for the fastest convergence to equilibrium,
 but the accurate evaluation of critical damping involves the expensive
 computation of eigenvalues.
 So, we adopt a local nonviscous damping [Cundall, 1987].
 In this approach, the magnitude of the damping force on a node is proportional
 to the magnitude of the unbalanced force, and the direction of damping
 force is such that energy is always dissipated.
 The velocity update can then be written as follows:
\end_layout
\begin_layout Standard
+\shape default
+ here corresponds to the whole domain.
+ The local contribution to nodes corresponding to each term can be computed
+ by following the standard finite element procedure for linear tetrahedral
+ elements.
+ However, our method does not need to construct coefficient matrices such
+ as mass and stiffness matrices since it adopts an explicit time discretization.
+ The resultant momentum equation is
\begin_inset Formula \begin{equation}
\nu_{ia}^{t+\frac{\Delta t}{2}}=\nu_{ia}^{t\frac{\Delta t}{2}}+\left(r_{ia}f_{ia}^{d}\right)\frac{\Delta t}{M_{a}}\label{eq:velocity update}\end{equation}
+M^{n}\frac{Dv_{i}^{n}}{Dt}=\frac{1}{3}T_{i}^{[n]}+\frac{1}{4}\rho^{[n]}g_{i}V^{[n]},\label{eq:momentumdiscrete}\end{equation}
\end_inset
+where the superscript
+\begin_inset Formula $n$
+\end_inset
\end_layout
+ represents values evaluated at the global node
+\begin_inset Formula $n$
+\end_inset
\begin_layout Standard
where
\begin_inset Formula $M_{a}$
+, the superscript
+\begin_inset Formula $[n]$
\end_inset
 is the fictitious nodal mass for node
\begin_inset Formula $a$
+ means the sum of contributions from all the tetrahedra having the global
+ node n as an apex,
+\begin_inset Formula $T_{i}$
\end_inset
, and
\begin_inset Formula $f_{ia}^{d}$
+ is the traction that is defined as
+\shape italic
+
+\begin_inset Formula $\sigma_{ij}n_{j}$
\end_inset
 is the
\emph on
i
\emph default
th component damping force on node
\begin_inset Formula $a$
+
+\shape default
+ and evaluated on a face of one of the contributing tetrahedra.
+ Any applied traction boundary conditions should be explicitly computed
+ and added on the right hand side of (
+\begin_inset LatexCommand ref
+reference "eq:momentumdiscrete"
+
\end_inset
+).
+ The nodal mass
+\begin_inset Formula $M^{n}$
+\end_inset
+
+ is not the actual inertial mass but an adjusted one to satisfy a local
+ stability criterion discussed in the section
+\begin_inset LatexCommand ref
+reference "sec:appendix_massscaling"
+
+\end_inset
+
.
 The damping force,
\begin_inset Formula $f_{ia}^{d}$
+ The correspondence between an apex and a face for the traction calculation
+ is determined as in the derivation of the expression, (
+\begin_inset LatexCommand ref
+reference "eq:formula_fi"
+
\end_inset
, is given by
+).
+ Note that the factor of 1/3 in the traction term is inherited from (
+\begin_inset LatexCommand ref
+reference "eq:formula_fi"
+
+\end_inset
+
+) and the factor of 1/4 in the body force term implies that the nodal contributi
+on takes one quarter of a tetrahedron's volumedependent quantity.
\end_layout
\begin_layout Standard
+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 flux) at each node is used to update velocities
+ and node coordinates (or temperature).
+\end_layout
+
+\begin_layout Subsection
+Solution Scheme
+\end_layout
+
+\begin_layout Standard
+We seek static or quasistatic solutions through a dynamic relaxation method.
+ Instead of adding a usual velocitydependent friction term, we adopt a
+ local nonviscous damping scheme
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep{Cundall1987}
+\end_layout
+
+\end_inset
+
+:
\begin_inset Formula \begin{equation}
f_{ia}^{d}=\alphar_{ia}\mathrm{sgn}\left(\nu_{ia}^{t\frac{\Delta t}{2}}\right)\label{eq:damping forcec}\end{equation}
+F_{i}^{damped}=F_{i}\alpha\text{ sgn}(v_{i})F_{i},\label{eq:forcedamped}\end{equation}
\end_inset
+ where
+\begin_inset Formula $F_{i}$
+\end_inset
\end_layout
+ is the
+\begin_inset Formula $i$
+\end_inset
\begin_layout Standard
where
+th component of the residual force vector (the right hand side of (
+\begin_inset LatexCommand ref
+reference "eq:momentumdiscrete"
+
+\end_inset
+
+)),
+\shape italic
+
\begin_inset Formula $\alpha$
\end_inset
 is the proportionality constant for damping and ??? (ends in G3 paper)
 TODO: what need to happen here?
+
+\shape default
+ is a positive coefficient less than 1, sgn
+\begin_inset Formula $(v_{i})$
+\end_inset
+
+ returns the sign of the
+\begin_inset Formula $i$
+\end_inset
+
+th component of velocity,
+\begin_inset Formula $v_{i}$
+\end_inset
+
+.
+ Once net forces are assembled and damped, velocity at that node is updated
+ using a forward Euler method:
+\begin_inset Formula \begin{equation}
+v(t+\frac{\Delta t}{2})=v(t\frac{\Delta t}{2})+\Delta t\frac{F_{i}^{damped}}{M}\label{eq:velupdate}\end{equation}
+
+\end_inset
+
+
+\begin_inset Formula \begin{equation}
+x(t+\Delta t)=x(t)+\Delta tv(t+\frac{\Delta t}{2}).\label{eq:posupdate}\end{equation}
+
+\end_inset
+
+ Damping is irrelevant to the update of temperature field, but the same
+ forward Euler method is used.
\end_layout
\begin_layout Standard
Due to the explicit nature of the solution scheme, the solution is conditionally
 stable when the numerical information speed of the mesh is more than the
 physical information speed.
 This condition can be expressed as
+\begin_layout Subsection
+Mass scaling for numerical stability
\end_layout
\begin_layout Standard
\begin_inset Formula \[
\Delta t\leq\frac{\Delta x}{C_{p}}\]
+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 userspecified time step can be locally satisfied.
+ The stability condition to be satisfied, however, is not the same as in
+ the CFL condition, i.e.,
+\shape italic
+\begin_inset Formula $\Delta t$
\end_inset
+\shape default
+
+\begin_inset Formula $\leq$
+\end_inset
+
+ (
+\begin_inset Formula $l_{min}/v_{p}$
+\end_inset
+
+), where
+\shape italic
+
+\begin_inset Formula $\Delta t$
+\end_inset
+
+
+\shape default
+ is the time step,
+\begin_inset Formula $l_{min}$
+\end_inset
+
+ is the minimum element size, and
+\begin_inset Formula $v_{p}$
+\end_inset
+
+ is the P wave velocity.
+ Instead, through an analogy of continuum to an infinite massspring system,
+ we use a criterion that does not explicitly include length scale and P
+ wave velocity
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep[see Ch.
+ 9
\end_layout
\begin_layout Standard
+
+in]{Bathe1996}
+\end_layout
+
+\end_inset
+
+:
+\begin_inset Formula \begin{equation}
+\Delta t\leq\frac{T}{\pi},\label{eq:dtcrit}\end{equation}
+
+\end_inset
+
+ where
+\begin_inset Formula $T$
+\end_inset
+
+ is the period of system,
+\begin_inset Formula $2\pi(m/k)^{1/2}$
+\end_inset
+
+,
+\begin_inset Formula $m$
+\end_inset
+
+ is a point mass, and
+\begin_inset Formula $k$
+\end_inset
+
+ is the stiffness of the spring attached to the point mass.
+ Now, reducing the infinite series of mass and springs in one dimension
+ to a single massspring system, the stiffness of that single system becomes
+
+\begin_inset Formula $4k$
+\end_inset
+
+, leading to an expression for the mass scaling:
+\begin_inset Formula \begin{equation}
+m\geq k(\Delta t)^{2}.\label{eq:critmass}\end{equation}
+
+\end_inset
+
+ For a given size of
+\begin_inset Formula $\Delta t$
+\end_inset
+
+, the nodal mass is adjusted according to (
+\begin_inset LatexCommand ref
+reference "eq:critmass"
+
+\end_inset
+
+) to automatically satisfy the stability critetion, (
+\begin_inset LatexCommand ref
+reference "eq:dtcrit"
+
+\end_inset
+
+).
+ The value of
+\begin_inset Formula $k$
+\end_inset
+
+ is computed by equating internal force contribution at a node with
+\begin_inset Formula $ku_{i}$
+\end_inset
+
+:
+\begin_inset Formula \begin{equation}
+\begin{split}\frac{1}{3}T_{i} & =ku_{i}\Rightarrow\\
+\frac{1}{3}(\lambda+2\mu)(\dot{\epsilon}_{ii}dt)n_{i}S & =k(v_{i}dt)\text{ (no sum)},\end{split}
+\label{eq:eqfork}\end{equation}
+
+\end_inset
+
+where only the volumetric contribution from internal forces is taken into
+ account.
+ By substituting the approximation for the partial derivative (
+\begin_inset LatexCommand ref
+reference "eq:formula_fi"
+
+\end_inset
+
+) into the above equation and dividing both sides by
+\begin_inset Formula $v_{i}dt$
+\end_inset
+
+, we obtain
+\begin_inset Formula \begin{equation}
+k_{i}^{l}=\frac{1}{9V}(\lambda+2\mu)(n_{i}^{l}S^{l})^{2},\label{eq:kdef}\end{equation}
+
+\end_inset
+
where
\begin_inset Formula $C_{p}=\sqrt{\frac{\kappa+4G/3}{\rho}}$
+\begin_inset Formula $l$
\end_inset
 is the Pwave speed and
\begin_inset Formula $\Delta x$
+ is the local index for apexes of a tetrahedron and the surfacerelated
+ quantities are computed on the corresponding face of the tetrahedron.
+ Finally, a tetrahedron's contribution to the scaled mass is given as
+\begin_inset Formula \begin{equation}
+m^{l}=\frac{\lambda+2\mu}{9V}\max[(n_{i}^{l}S^{l})^{2},i=1,\dots,3].\label{eq:scaledmass}\end{equation}
+
\end_inset
 is the element size.
 Since the mass matrix is fictitious, we can choose a mass matrix in such
 a way that convergence is optimal.
 The optimal convergence is obtained when local values of critical time
 steps for each element are equal.
 The condition for this can be obtained as
+
\end_layout
\begin_layout Standard
\begin_inset Formula \begin{equation}
M_{ia}=\sum\frac{\left(\kappa+4G/3\right)\Delta x_{\mathrm{max}}^{2}}{6V}\label{eq:condition for optimal convergence}\end{equation}
+As in the standard FEM, appropriate mappings between local and global indices
+ are required.
+\end_layout
\end_inset
+\begin_layout Subsection
+Constitutive update
+\end_layout
+\begin_layout Standard
+SNAC uses a general elastoviscoplastic rheological model to update the
+ Cauchy stress tensor
+\begin_inset ERT
+status collapsed
+\begin_layout Standard
+
+
+\backslash
+citep[e.g.,]{Albe_etal2000}
\end_layout
\begin_layout Standard
where summation is over all the tetrahedral elements sharing the node a,
 and the
\begin_inset Formula $\Delta x_{\mathrm{max}}$
\end_inset
 is the maximum propagation distance.
+.
+ First, the initial guess of stress is acquired by the Maxwell viscoelastic
+ constitutive law
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep{Poli_etal1993a}
\end_layout
\begin_layout Subsection
Constitutive Update
+\end_inset
+
+.
+ If this initial guess exceeds a given yield stress, it is projected onto
+ the yield surface using a return mapping method
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep{SimoHugh2004}
\end_layout
\begin_layout Standard
SNAC uses a general elastoviscoplastic (or viscoplastic) rheological model
 to update Cauchy stress based on strainrate.
 This viscoplastic model can deal with various constitutive laws that are
 typically used for the Earth's crustal and mantle as its limiting cases.
+\end_inset
+
+; otherwise, the viscoelastic stress update is retained.
+ This elastoviscoplastic model can deal with various constitutive laws
+ that are typically used for the Earth's crustal and mantle material as
+ its limiting cases.
For example, elastic, viscoelastic and elastoplastic material are realized
in the following cases:
\end_layout
\begin_layout Enumerate
Elastic material corresponds to the limit of infinite viscosity and infinite
 yield stress.
+Elastic material corresponds to the limit of infinite viscosity and yield
+ stress.
\end_layout
@@ 733,76 +1141,67 @@
\end_layout
\begin_layout Enumerate
Elasoplastic material corresponds to the infinite viscosity.
+Elasoplastic material corresponds to the inifinte viscosity.
\end_layout
\begin_layout Standard
Using the viscoplastic rheology is physically more realistic than using
one of the limiting cases listed above since all materials have dissipative
 mechanisms, however small they might be, and hence have viscosity.
+ mechanisms and hence viscosity.
This viscosity also provides a length scale for the problem of localization,
which in turn enables physically meaningful mesh independent solution when
 the mesh size is smaller than this length scale [e.g., Needleman, 1988].

+ the mesh size is smaller than this length scale.
\end_layout
\begin_layout Standard
Since the nodal variables are velocities and whose spatial gradients are
 deformation rates, we formulate constitutive update in terms of strain
+ deformation rates, we formulate the constitutive update in terms of strain
rate.
 The objective stress rate of choice is Jaumann or the corotational stress
 rate (
\begin_inset Formula $\sigma^{\Delta J}$
+ The objective stress rate of choice is the Jaumann or the corotational
+ stress rate (
+\begin_inset Formula $\Delta\sigma^{\Delta J}$
\end_inset
 ),
\begin_inset Formula $\sigma^{\Delta J}=C^{\Delta J}D$
\end_inset
+)
+\begin_inset ERT
+status collapsed

\end_layout

\begin_layout Standard
where
\end_layout
\begin_layout Standard
\begin_inset Formula \begin{equation}
D\,(D_{i,j}=\left(\frac{\partial\nu_{i}}{\partial x_{j}}+\frac{\partial\nu_{j}}{\partial x_{i}}\right)),\,\, W\,(W_{i,j}=\frac{1}{2}\left(\frac{\partial\nu_{i}}{\partial x_{j}}\frac{\partial\nu_{j}}{\partial x_{i}}\right))\label{eq:deformation and spin tensors}\end{equation}
\end_inset
+\backslash
+citep{RudnRice1975}
+\end_layout
 are the deformation and spin tensors and
\begin_inset Formula $C^{\Delta J}$
\end_inset
 is Jaumann tangent modulus.
 For the isotropic material Jaumann tangent modulus is the same as the elastic
 tangent modulus, in which case correction to the stresses due to rotation
 can be given as
\end_layout

\begin_layout Standard
+
\begin_inset Formula \begin{equation}
\Delta\sigma=\Delta\sigma^{C}+\Delta\sigma^{R}=\Delta\sigma^{C}+\left(W\sigma\sigma W\right)\Delta t,\label{eq:lots of deltas}\end{equation}
+\Delta\sigma^{\Delta J}=\frac{\partial(\Delta\sigma)}{\partial t}W\cdot\Delta\sigma\Delta\sigma\cdot W^{T},\label{eq:corotstress}\end{equation}
\end_inset
+ where
+\begin_inset Formula $W_{ij}=(1/2)(\partial v_{i}/\partial x_{j}$
+\end_inset
\end_layout
+
+\begin_inset Formula $\partial v_{j}/\partial x_{i})$
+\end_inset
\begin_layout Standard
where
\begin_inset Formula $\Delta\sigma^{C}$
+ are the components of spin tensor and
+\begin_inset Formula $\Delta\sigma$
\end_inset
 denotes the stress update by constitutive law and
\begin_inset Formula $\Delta\sigma^{R}$
+ is the increment of stress tensor.
+ Correction to the stresses due to rotation can be given as
+\begin_inset Formula \begin{equation}
+\sigma^{t+\Delta t}=\sigma^{t}+\Delta\sigma^{\Delta J}.\Delta t\label{eq:stressupdate}\end{equation}
+
\end_inset
 stands for the correction due to finite rotation, of the stress update.

+
\end_layout
\begin_layout Standard
@@ 897,8 +1296,8 @@
\begin_inset Formula $\sigma^{n+1}$
\end_inset
 onto the yield surface using a returnmapping algorithm [TODO need citation:
 Ortiz and Simo, 1986].
+ onto the yield surface using a returnmapping algorithm [Simo and Hughes,
+ 1999].
\end_layout
@@ 979,8 +1378,7 @@
\begin_layout Standard
In general, flow rule for frictional materials is nonassociative, i.e., flow
direction differs from the normal of the yield surface normal.
 The plastic flow potential for plastic flow [TODO  is there some duplication
 in the foregoing phrase?] can be given as
+ The plastic flow potential can be given as
\end_layout
\begin_layout Standard
@@ 1314,8 +1712,7 @@
\begin_layout Standard
SNAC is built up from the StGermain framework, an environment for the developmen
t of physically based scientific codes [Quenette, et al.
 TODOneed year].
+t of physically based scientific codes [Quenette, et al., 2005].
As frameworks, StGermain and Pyre differ in that StGermain provides the
infrastructural components needed to build complete codes while Pyre provides
the superstructure to couple codes together; for this reason the former
More information about the CIGCOMMITS
mailing list