[cig-commits] commit: Worked on governing equations section.

[Mercurial]

changeset:   53:924b986f337e
tag:         tip
user:        Brad Aagaard <baagaard at usgs.gov>
date:        Mon Aug 29 17:54:05 2011 -0700
files:       faultRup.tex
description:
Worked on governing equations section.


diff -r b11cd13e7c17 -r 924b986f337e faultRup.tex
--- a/faultRup.tex	Fri Aug 26 17:31:57 2011 -0700
+++ b/faultRup.tex	Mon Aug 29 17:54:05 2011 -0700
@@ -11,8 +11,9 @@
 \newcommand\matt[1]{{\color{blue}\bf [MATT: #1]}}
 \newcommand\charles[1]{{\color{green}\bf [CHARLES: #1]}}
 
-\newcommand{\refeqn}[1]{~(\ref{#1})}
 \newcommand{\tensor}[1]{\underline{#1}}
+\newcommand{\tpdt}{t+\Delta t}
+\newcommand{\tmdt}{t-\Delta t}
 
 \widowpenalty=999999
 \clubpenalty=999999
@@ -165,60 +166,213 @@ and verification of its implementation u
 % ------------------------------------------------------------------
 \section{Numerical Model of Fault Slip}\brad{Rough draft}
 
+In this section we summarize the formulation of the governing
+equations using the finite-element method. We augment the conventional
+finite-element formulation for elasticity with a domain decomposition
+approach \cite{Zienkiewicz:Taylor:2005} to implement the fault slip.
+\citeN{PyLith:manual:1.6.2} provides a step-by-step description of the
+formulation.
+
 We solve the elasticity equation including inertial terms,
 \begin{gather}
   \rho \frac{\partial^2\vec{u}}{\partial t^2} - \vec{f} 
   - \nabla \cdot \tensor{\sigma} = \vec{0} \text{ in }V, \\
   \tensor{\sigma} \cdot \vec{n} = \vec{T} \text{ on }S_T, \\
   \vec{u} = \vec{u}_0 \text{ on }S_u, \\
-  \tensor{R} \cdot \left( \vec{u}^{+} - \vec{u}^{-}\right) = \vec{d}
+  \tensor{R} \left( \vec{u}^{+} - \vec{u}^{-}\right) = \vec{d}
   \text{ on }S_f,
 \end{gather}
 where $\vec{u}$ is the displacement vector, $\rho$ is the mass
 density, $f$ is the body force vector, $\tensor{\sigma}$ is the Cauchy
 stress tensor, and $t$ is time. We specify tractions $\vec{T}$ on
-surface $S_f$, displacements $\vec{u_0}$ on surface $S_u$, and slip
-$\vec{d}$ on fault surface $S_f$ (we also consider the case of a fault
-constitutive model later in this section). The rotation matrix
-$\tensor{R}$ transforms vectors from the global coordinate system to
-the fault coordinate system. Because both $\vec{T}$ and $\vec{u}$ are
-vector quantities, there can be some spatial overlap of the surfaces
-$S_T$ and $S_u$; however, a degree of freedom cannot simultaneously be
-associated with both types of boundary conditions.
+surface $S_T$, displacements $\vec{u_0}$ on surface $S_u$, and slip
+$\vec{d}$ on fault surface $S_f$. The rotation matrix $\tensor{R}$
+transforms vectors from the global coordinate system to the fault
+coordinate system. Because both $\vec{T}$ and $\vec{u}$ are vector
+quantities, there can be some spatial overlap of the surfaces $S_T$
+and $S_u$; however, a degree of freedom cannot be associated with both
+types of boundary conditions simultaneously.
 
-Construct the weak form ...
+Following a conventional finite-element formulation (ignoring the
+fault surface for a moment), we construct the
+weak form by taking the dot product with a weighting function and
+setting the integral over the domain equal to zero,
+\begin{equation}
+  \int_{V} \left( \nabla \cdot \tensor{\sigma} + \vec{f} -
+    \rho\frac{\partial^{2}\vec{u}}{\partial t^{2}} \right) 
+  \cdot \vec{\phi}\, dV=0.
+\end{equation}
+The weighting function $\vec{\phi}$ is a piecewise differential vector
+field with $\vec{\phi} = \vec{0}$ on $S_u$. After some algebra and
+use of the boundary conditions, we
+have
+\begin{equation}
+- \int_{V} \tensor{\sigma} : \nabla\vec{\phi} \, dV
++ \int_{S_T} \vec{T} \cdot \vec{\phi}\, dS
++ \int_{V} \vec{f} \cdot \vec{\phi} \, dV 
+- \int_{V} \rho \frac{\partial^{2}\vec{u}}{\partial t^{2}} \cdot
+\vec{\phi}\, dV
+=0.
+\end{equation}
 
-We express the trial solution and weighting function as linear
-combinations of basis functions ...
+Using a domain decomposition approach, we consider the fault surface
+as an interior boundary between two domains. We assign a fault normal
+direction to this interior boundary and ``positive'' and ``negative''
+labels to the two sides of the fault, such that the fault normal is
+the vector from the negative side of the fault to the positive side of
+the fault. Slip on the fault is the motion on the positive side
+relative to the motion on the negative side. Slip on the fault also
+corresponds to equal and opposite tractions on the positive and
+negative sides of the fault, which we impose using Lagrange
+multipliers with $\vec{l}^{+} - \vec{l}^{-} = 0$.
 
-We want to solve this equation for the unknown coefficients ...
-subject to ...
+Recognizing that the tractions on the fault surface are analogous to
+the boundary tractions, we add in the contributions from integrating
+the Lagrange multipliers (fault tractions) over the fault surface,
+\begin{equation}
+- \int_{V} \tensor{\sigma} : \nabla\vec{\phi} \, dV
++ \int_{S_T} \vec{T} \cdot \vec{\phi}\, dS
++ \int_{S_f^{+}} \vec{l} \cdot \vec{\phi}\, dS
+- \int_{S_f^{-}} \vec{l} \cdot \vec{\phi}\, dS
++ \int_{V} \vec{f} \cdot \vec{\phi} \, dV 
+- \int_{V} \rho \frac{\partial^{2}\vec{u}}{\partial t^{2}} \cdot
+\vec{\phi}\, dV
+=0.
+\end{equation}
+We also construct the weak form for the constraint associated with
+slip on the fault by taking the dot product of the constraint equation
+with the weighting function and setting the integral over the fault
+surface to zero,
+\begin{equation}
+  \int_{S_f} \left( \vec{u}^{+} - \vec{u}^{-} -
+    \tensor{R}_T \vec{d} \right) \cdot \vec{\phi} \, dS = 0.
+\end{equation}
+The rotation matrix $\tensor{R}_T$ is the inverse of
+$\tensor{R}$ and is composed of blocks that are the transpose of the
+corresponding blocks in $\tensor{R}$.
 
-Quasi-static
+We express the trial solution $\vec{u}$, weighting function
+$\vec{\phi}$, and Lagrange multipliers $\vec{l}$ as linear combinations
+of basis functions,
+\begin{gather}
+\vec{u} = \sum_{m} \vec{u^{m}} N^{m}, \\
+\vec{\phi} = \sum_{n} \vec{a^{n}} N^{n}, \\
+\vec{l} = \sum_{p} \vec{l^{p}} N^{p}.
+\end{gather}
+Because the weighting function is zero on $S_u$, the number of basis
+functions for the trial solution $\vec{u}$ is generally greater than
+the number of basis functions for the weighting function $\vec{\phi}$,
+i.e., $m > n$. The basis functions for the Lagrange multipliers are
+associated with the fault surface, which is a lower dimension than the
+domain, so $p \ll n$ in most cases. The weighting function is
+arbitrary, so the integrands must be zero for all $\vec{a^n}$, which
+leads to
+\begin{gather}
+- \int_{V} \tensor{\sigma} : \nabla N^{n} \, dV
++ \int_{S_T} \vec{T} N^n \, dS
++ \int_{S_f^{+}} \sum_{p} \vec{l^p} N^p N^n \, dS
+- \int_{S_f^{-}} \sum_{p} \vec{l^p} N^p N^n \, dS
++ \int_{V} \vec{f} N^n \, dV
+- \int_{V} \rho\sum_{m} \frac{\partial^2 \vec{u^{m}}}{\partial
+  t^{2}}N^{m}N^{n}\, dV
+=\vec{0}, \\
+%
+  \int_{S_f} \left( \sum_m \vec{u^{m+}} N^m - \sum_m \vec{u^{m-}} N^m -
+    \tensor{R}_T \vec{d} \right) N^p \, dS = \vec{0}.
+\end{gather}
+We want to solve these equations for the coefficients $\vec{u^m}$ and
+$\vec{l^p}$ subject to $\vec{u} = \vec{u}_0 \text{ on }S_u$.
 
-Dynamic
+For nonlinear bulk rheologies it is convenient to work with the
+increment in stress and strain, so we formulate the solution of the
+equations in terms of the increment in the solution from time $t$ to
+$\tpdt$ rather than the solution at time $\tpdt$.
+Consequently, rather than constructing a system with the form
+$\tensor{A} \vec{u}(\tpdt) = \vec{b}(\tpdt)$, we construct a system
+with the form $\tensor{A} \vec{du} = \vec{b}(\tpdt) - \tensor{A}
+\vec{u}(t)$. 
 
-Impose fault slip using Lagrange multipliers ...
+\subsection{Quasi-static Simulations}
 
+For quasi-static simulations we ignore the inertial term and
+time-dependence only enters through the constitutive models
+and the loading conditions. Considering the deformation at time
+$\tpdt$,
+\begin{gather}
+\label{eqn:quasi-static:residual:domain}
+- \int_{V} \tensor{\sigma}(\tpdt) : \nabla N^{n} \, dV
++ \int_{S_T} \vec{T}(\tpdt) N^n \, dS
++ \int_{S_f^{+}} \sum_{p} \vec{l^p}(\tpdt) N^p N^n \, dS
+- \int_{S_f^{-}} \sum_{p} \vec{l^p}(\tpdt) N^p N^n \, dS
++ \int_{V} \vec{f}(\tpdt) N^n \, dV
+=\vec{0}, \\
+%
+\label{eqn:quasi-static:residual:fault}
+  \int_{S_f} \left( \sum_m \vec{u^{m+}}(\tpdt) N^m - \sum_m \vec{u^{m-}}(\tpdt) N^m -
+    \tensor{R}_T \vec{d}(\tpdt) \right) N^p \, dS = \vec{0}.
+\end{gather}
+We solve these equations using the Portable, Extensible Toolkit for
+Scientific Computation (PETSc), which provides a suite of tools for
+solving linear systems of algebraic equstions with parallel processing
+\cite{PETSc}. In solving the system, we compute the residual (i.e.,
+$\vec{r} = \vec{b} - \tensor{A} \vec{u}$ and the Jacobian of the
+system ($A$). The residual is simply the left hand sides of
+equations~(\ref{eqn:quasi-static:residual:domain})
+and~(\ref{eqn:quasi-static:residual:fault}).
 
-Kinematic rupture model
+The Jacobian of the system is the action (operations) on the increment
+in the solution.  To find the portion of the Jacobian associated with
+equation~(\ref{eqn:quasi-static:residual:domain}), we let
+$\tensor{\sigma}(\tpdt) = \tensor{\sigma}(t) +
+\tensor{d\sigma}(t)$. The action on the increment of the solution is
+associated with the term $\tensor{d\sigma}(t)$. We approximate the
+increment in the stress tensor using linear elasticity and
+infinitesimal strains,
+\begin{equation}
+  \tensor{d\sigma}(t) = \frac{1}{2} \tensor{C}(t) (\nabla + \nabla^T)
+  \vec{u}(t),
+\end{equation}
+where $\tensor{C}$ is the fourth order tensor of elastic
+constants. Substituting into the first term in
+equation~(\ref{eqn:quasi-static:residual:domain}), after some algebra,
+we find this portion of the Jacobian is
+\begin{equation}
+  \tensor{K} = \frac{1}{2} \int_v \tensor{C} (\nabla +
+  \nabla^T) : \nabla N^n  \, dV.
+\end{equation}
+The portion of the Jacobian associated with
+equation~(\ref{eqn:quasi-static:residual:fault}) is
+\begin{equation}
+  \tensor{L} = \int_{S_f} (N^{m+}-N^{m-})N^p \, dS.
+\end{equation}
+Thus, the Jacobian of the entire system has the form,
+\begin{equation}
+  \tensor{A} = 
+  \left( \begin{array}{cc}
+      \tensor{K} & \tensor{L}^T \\ \tensor{L} & \tensor{0} 
+    \end{array} \right).
+\end{equation}
+
+\subsection{Dynamic Simulations}
+
+\subsection{Prescribed Fault Rupture}
 
 Slip-time functions
 
-Dynamic (spontaneous) rupture model
+\subsection{Spontaneous Fault Rupture}
 
 Fault constitutive models
 
 % ------------------------------------------------------------------
-\subsection{Finite-Element Mesh Processing}\matt{Rough draft}
+\section{Finite-Element Mesh Processing}\matt{Rough draft}
   adjusting topology
 
 % ------------------------------------------------------------------
-\subsection{Fault Discretization}\brad{Rough draft}
+\section{Fault Discretization}\brad{Rough draft}
   cohesive cells
 
 % ------------------------------------------------------------------
-\subsection{Solver Customization}\matt{Rough draft}
+\section{Solver Customization}\matt{Rough draft}
   \begin{itemize}
     \item saddle point
     \item custom preconditioner