[cig-commits] r16574 - short/3D/PyLith/trunk/doc/userguide/governingeqns
brad at geodynamics.org
brad at geodynamics.org
Wed Apr 21 17:13:07 PDT 2010
Author: brad
Date: 2010-04-21 17:13:06 -0700 (Wed, 21 Apr 2010)
New Revision: 16574
Modified:
short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx
Log:
More work on governing equations.
Modified: short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx 2010-04-22 00:07:59 UTC (rev 16573)
+++ short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx 2010-04-22 00:13:06 UTC (rev 16574)
@@ -765,6 +765,28 @@
Finite-Element Formulation of Elasticity Equation
\end_layout
+\begin_layout Standard
+We formulate a set of algebraic equations using Galerkin's method.
+ We consider a trial solution,
+\begin_inset Formula $\vec{u}$
+\end_inset
+
+, that is a piecewise differntiable vector field and satisfies the Dirichlet
+ boundary conditions on
+\begin_inset Formula $S_{u}$
+\end_inset
+
+, and a weighting function,
+\begin_inset Formula $\vec{\phi}$
+\end_inset
+
+, that is a piecewise differentialble vector field and is zero on
+\begin_inset Formula $S_{u}$
+\end_inset
+
+.
+\end_layout
+
\begin_layout Subsection
Index Notation
\end_layout
@@ -784,20 +806,7 @@
\end_inset
We construct the weak form by computing the dot product of the wave equation
- and a trial function and setting the integral over the domain to zero.
- The trial function is a piecewise differentiable vector field,
-\begin_inset Formula $\phi_{i}$
-\end_inset
-
-, where
-\begin_inset Formula $\phi_{i}=0$
-\end_inset
-
- on
-\begin_inset Formula $S_{u}.$
-\end_inset
-
- Hence our weak form is
+ and weighting function and setting the integral over the domain to zero:
\begin_inset Note Greyedout
status open
@@ -810,12 +819,12 @@
\begin_inset Formula \begin{gather}
\int_{V}\left(\sigma_{ij,j}+f_{i}-\rho\ddot{u}_{i}\right)\phi_{i}\, dV=0\text{, or }\\
-\int_{V}\sigma_{ij,j}\phi_{i}\: dV+\int_{V}f_{i}\phi_{i}\: dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\: dV=\vec{0}.\end{gather}
+\int_{V}\sigma_{ij,j}\phi_{i}\: dV+\int_{V}f_{i}\phi_{i}\: dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\: dV=0.\end{gather}
\end_inset
Consider the divergence theorem applied to the dot product of the stress
- tensor and the trial function,
+ tensor and the weighting function,
\begin_inset Formula $\sigma_{ij}\phi_{i}$
\end_inset
@@ -848,39 +857,43 @@
\end_inset
-Now,
-\begin_inset Formula $\sigma_{ij}\phi_{i,j}$
+Turning our attention to the second term, we separate the integration over
+
+\begin_inset Formula $S$
\end_inset
- is a scalar, so it is symmetric,
-\begin_inset Formula \begin{equation}
-\sigma_{ij}\phi_{i,j}=\sigma_{ji}\phi_{j,i},\end{equation}
-
+ into integration over
+\begin_inset Formula $S_{T}$
\end_inset
-and we know that
-\begin_inset Formula $\sigma_{ij}$
+ and
+\begin_inset Formula $S_{u}$
\end_inset
- is symmetric, so
-\begin_inset Formula \begin{equation}
-\sigma_{ij}\phi_{i,j}=\sigma_{ij}\phi_{j,i},\end{equation}
+,
+\begin_inset Note Greyedout
+status open
+\begin_layout Plain Layout
+Add fault constraint
+\end_layout
+
\end_inset
-which means
+
\begin_inset Formula \begin{equation}
-\phi_{i,j}=\phi_{j,i},\end{equation}
+-\int_{V}\sigma_{ij}\phi_{i,j}\, dV+\int_{S_{T}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{S_{u}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0,\end{equation}
\end_inset
-which we can write as
-\begin_inset Formula \begin{equation}
-\phi_{i,j}=\frac{1}{2}(\phi_{i,j}+\phi_{j,i}).\end{equation}
+and recognize that
+\begin_inset Formula \begin{gather}
+\sigma_{ij}n_{i}=T_{i}\text{ on }S_{T}\text{ and}\\
+\phi_{i}=0\text{ on }S_{u},\end{gather}
\end_inset
-Substituting into the first term gives
+so that the equation reduces to
\begin_inset Note Greyedout
status open
@@ -892,47 +905,101 @@
\begin_inset Formula \begin{equation}
--\int_{V}\frac{1}{2}\sigma_{ij}\left(\phi_{i,j}+\phi_{j,i}\right)\, dV+\int_{S}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.\end{equation}
+-\int_{V}\sigma_{ij}\phi_{i,j}\: dV+\int_{S_{T}}T_{i}\phi_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.\label{eq:elasticity:integral}\end{equation}
\end_inset
-Turning our attention to the second term, we separate the integration over
-
-\begin_inset Formula $S$
+We express the trial solution and weighting function as linear combinations
+ of basis functions,
+\begin_inset Formula \begin{gather*}
+u=\sum_{m}a_{i}^{m}N^{m},\\
+\phi=\sum_{n}c_{i}^{n}N^{n}.\end{gather*}
+
\end_inset
- into integration over
-\begin_inset Formula $S_{T}$
+Note that because the trial solution satisfies the Dirichlet boundary condition,
+ the number of basis functions for
+\begin_inset Formula $u$
\end_inset
- and
-\begin_inset Formula $S_{u}$
+ is generally greater than the number of basis functions for
+\begin_inset Formula $\phi$
\end_inset
-,
-\begin_inset Note Greyedout
-status open
+, i.e.,
+\begin_inset Formula $m>n$
+\end_inset
-\begin_layout Plain Layout
-Add fault constraint
+.
+ Substituting in the expressions for the trial solution and weighting function
+ yields
+\begin_inset Formula \begin{gather*}
+-\int_{V}\sigma_{ij}\sum_{n}c_{i}^{n}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}\sum_{n}c_{i}^{n}N^{n}\, dS+\int_{V}f_{i}\sum_{n}c_{i}^{n}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}\sum_{n}c_{i}^{n}N^{n}\ dV=0,\text{ or}\\
+\sum_{n}c_{i}^{n}(-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}N^{n}\ dV)=0.\end{gather*}
+
+\end_inset
+
+ Because the weighting function is aribtrary, this equation must hold for
+ all
+\begin_inset Formula $c^{n}$
+\end_inset
+
+, so that the quantity in parathesis is zero for each
+\begin_inset Formula $c^{n}$
+\end_inset
+
+
+\begin_inset Formula \begin{equation}
+-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}N^{n}\ dV=\vec{0}.\label{eq:elasticity:integral-1}\end{equation}
+
+\end_inset
+
+We want to solve this equation for the unknown coefficients
+\begin_inset Formula $a_{i}^{m}$
+\end_inset
+
+.
+
\end_layout
+\begin_layout Subsection
+Old
+\end_layout
+
+\begin_layout Standard
+Now,
+\begin_inset Formula $\sigma_{ij}\phi_{i,j}$
\end_inset
+ is a scalar, so it is symmetric,
+\begin_inset Formula \begin{equation}
+\sigma_{ij}\phi_{i,j}=\sigma_{ji}\phi_{j,i},\end{equation}
+\end_inset
+
+and we know that
+\begin_inset Formula $\sigma_{ij}$
+\end_inset
+
+ is symmetric, so
\begin_inset Formula \begin{equation}
--\int_{V}\frac{1}{2}\sigma_{ij}(\phi_{i,j}+\phi_{j,i})\, dV+\int_{S_{T}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{S_{u}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0,\end{equation}
+\sigma_{ij}\phi_{i,j}=\sigma_{ij}\phi_{j,i},\end{equation}
\end_inset
-and recognize that
-\begin_inset Formula \begin{gather}
-\sigma_{ij}n_{i}=T_{i}\text{ on }S_{T}\text{ and}\\
-\phi_{i}=0\text{ on }S_{u},\end{gather}
+which means
+\begin_inset Formula \begin{equation}
+\phi_{i,j}=\phi_{j,i},\end{equation}
\end_inset
-so that the equation reduces to
+which we can write as
+\begin_inset Formula \begin{equation}
+\phi_{i,j}=\frac{1}{2}(\phi_{i,j}+\phi_{j,i}).\end{equation}
+
+\end_inset
+
+Substituting into the first term gives
\begin_inset Note Greyedout
status open
@@ -944,12 +1011,11 @@
\begin_inset Formula \begin{equation}
--\int_{V}\frac{1}{2}\sigma_{ij}(\phi_{i,j}+\phi_{j,i})\: dV+\int_{S_{T}}T_{i}\phi_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.\label{eq:elasticity:integral}\end{equation}
+-\int_{V}\frac{1}{2}\sigma_{ij}\left(\phi_{i,j}+\phi_{j,i}\right)\, dV+\int_{S}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.\end{equation}
\end_inset
-This is the equation we want to solve.
-
+
\end_layout
\begin_layout Section
@@ -960,13 +1026,13 @@
In quasi-static problems we neglect the inertial terms, so equation
\begin_inset CommandInset ref
LatexCommand eqref
-reference "eq:elasticity:integral"
+reference "eq:elasticity:integral-1"
\end_inset
reduces to
\begin_inset Formula \[
--\int_{V}\frac{1}{2}\sigma_{ij}(\phi_{i,j}+\phi_{j,i})\: dV+\int_{S_{T}}T_{i}\phi_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV=0.\]
+-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV=\vec{0}.\]
\end_inset
@@ -978,7 +1044,7 @@
,
\begin_inset Formula \begin{equation}
--\int_{V}\frac{1}{2}\sigma_{ij}(t+\Delta t)(\phi_{i,j}+\phi_{j,i})\: dV+\int_{S_{T}}T_{i}(t+\Delta t)\phi_{i}\, dS+\int_{V}f_{i}(t+\Delta t)\phi_{i}\, dV=0.\label{eq:elasticity:integral:t+dt}\end{equation}
+-\int_{V}\sigma_{ij}(t+\Delta t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV=\vec{0}.\label{eq:elasticity:integral:t+dt}\end{equation}
\end_inset
@@ -986,26 +1052,14 @@
of equations, involving the residual and Jacobian.
The residual is simply
\begin_inset Formula \[
-R=-\int_{V}\frac{1}{2}\sigma_{ij}(t+\Delta t)(\phi_{i,j}+\phi_{j,i})\: dV+\int_{S_{T}}T_{i}(t+\Delta t)\phi_{i}\, dS+\int_{V}f_{i}(t+\Delta t)\phi_{i}\, dV.\]
+\vec{R}=-\int_{V}\sigma_{ij}(t+\Delta t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV.\]
\end_inset
-In the finite-element discretization, the trial functions are scalar quantities,
- i.e., the same trial function is used for each component, and we employ quadratur
-e to integrate the quantities over the domain of each finite-element cell.
- We formulate the residual using
-\begin_inset Note Greyedout
-status open
-
-\begin_layout Plain Layout
-Need to refine the finite-element implementation stuff.
-\end_layout
-
-\end_inset
-
-
+We employ numerical quadrature in the finite-element discretization and
+ replace the integrals with sums over the cells and quadrature points,
\begin_inset Formula \[
-R=-\sum_{\text{vol cells}}\sum_{\text{quad pt}}\frac{1}{2}\sigma_{ij}(x_{q},t+\Delta t)(\phi_{,j}(x_{q})+\phi_{,i}(x_{q}))\: w_{q}|J_{m}(x_{q})|+\sum_{\text{vol cells}}\sum_{\text{quad pt}}f_{i}(x_{q},t+\Delta t)\phi(x_{q})\, w_{q}|J_{m}(x_{q})|+\sum_{\text{tract cells}}\sum_{\text{quad pt}}T_{i}(x_{q},t+\Delta t)\phi_{i}(x_{q})\, w_{q}|J_{m}(x_{q})|,\]
+\vec{R}=-\sum_{\text{vol cells}}\sum_{\text{quad pts}}\sigma_{ij}(x_{q},t+\Delta t)N^{n}(x_{q})\: w_{q}|J_{cell}(x_{q})|+\sum_{\text{vol cells}}\sum_{\text{quad pt}s}f_{i}(x_{q},t+\Delta t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|+\sum_{\text{tract cells}}\sum_{\text{quad pts}}T_{i}(x_{q},t+\Delta t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|,\]
\end_inset
@@ -1018,7 +1072,7 @@
\end_inset
are the weights of the quadrature points, and
-\begin_inset Formula $|J_{m}(x_{q})|$
+\begin_inset Formula $|J_{cell}(x_{q})|$
\end_inset
is the determinant of the Jacobian matrix evaluated at the quadrature points
@@ -1037,7 +1091,7 @@
Substituting into equation
\begin_inset CommandInset ref
LatexCommand eqref
-reference "eq:elasticity:integral:t+dt"
+reference "eq:elasticity:integral:dynamic:t"
\end_inset
@@ -1125,7 +1179,7 @@
\begin_inset Formula \[
-J=\sum_{\text{vol cells }}\sum_{\text{quad pt}}\frac{1}{4}C_{ijkl}(\phi_{,l}(x_{q})+\phi_{,k}(x_{q}))(\phi_{,j}(x_{q})+\phi_{,i}(x_{q}))\ w_{q}|J_{m}(x_{q})|\]
+J=\sum_{\text{vol cells }}\sum_{\text{quad pts}}\frac{1}{4}C_{ijkl}(\phi_{,l}(x_{q})+\phi_{,k}(x_{q}))(\phi_{,j}(x_{q})+\phi_{,i}(x_{q}))\ w_{q}|J_{m}(x_{q})|\]
\end_inset
@@ -1420,6 +1474,61 @@
\end_layout
\begin_layout Standard
+Time-dependence enters through the constitutive relationships, loading condition
+s, and the intertial terms.
+ We consider the deformation at time
+\begin_inset Formula $t$
+\end_inset
+
+,
+\begin_inset Formula \begin{equation}
+-\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t)N^{n}\, dS+\int_{V}f_{i}(t)N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ dV=\vec{0}.\label{eq:elasticity:integral:dynamic:t}\end{equation}
+
+\end_inset
+
+We solve this equation through formulation of a linear algebraic system
+ of equations, involving the residual and Jacobian.
+ The residual is simply
+\begin_inset Formula \[
+\vec{R}=-\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t)N^{n}\, dS+\int_{V}f_{i}(t)N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ dV.\]
+
+\end_inset
+
+We employ numerical quadrature in the finite-element discretization and
+ replace the integrals with sums over the cells and quadrature points,
+\begin_inset Formula \[
+\vec{R}=-\sum_{\text{vol cells}}\sum_{\text{quad pts}}\sigma_{ij}(x_{q},t+\Delta t)N^{n}(x_{q})\: w_{q}|J_{cell}(x_{q})|+\sum_{\text{vol cells}}\sum_{\text{quad pt}s}f_{i}(x_{q},t+\Delta t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|+\sum_{\text{tract cells}}\sum_{\text{quad pts}}T_{i}(x_{q},t+\Delta t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|,\]
+
+\end_inset
+
+where
+\begin_inset Formula $x_{q}$
+\end_inset
+
+ are the coordinates of the quadrature points,
+\begin_inset Formula $w_{q}$
+\end_inset
+
+ are the weights of the quadrature points, and
+\begin_inset Formula $|J_{cell}(x_{q})|$
+\end_inset
+
+ is the determinant of the Jacobian matrix evaluated at the quadrature points
+ associated with mapping the reference cell to the actual cell.
+ The quadrature scheme for the integral over the tractions is one dimension
+ lower than the one used in integrating the terms for the volume cells.
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}N^{n}\ dV=\vec{0}.\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
Using the central difference method to approximate the acceleration (and
velocity),
\begin_inset Formula \begin{gather}
More information about the CIG-COMMITS
mailing list