/[escript]/trunk/doc/user/Models.tex
ViewVC logotype

Contents of /trunk/doc/user/Models.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3989 - (show annotations)
Tue Sep 25 02:21:54 2012 UTC (6 years, 10 months ago) by jfenwick
File MIME type: application/x-tex
File size: 11615 byte(s)
More copyright fixes.
pyvisi traces removed.
Some install doco
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2012 by University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Open Software License version 3.0
8 % http://www.opensource.org/licenses/osl-3.0.php
9 %
10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 % Development since 2012 by School of Earth Sciences
12 %
13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14
15 \chapter{Models}
16 \label{MODELS CHAPTER}
17
18 The following sections give a brief overview of the model classes and their
19 corresponding methods.
20
21 \input{stokessolver}
22 %\input{darcyfluxNew}
23 \input{darcyflux}
24 %\input{levelsetmodel}
25
26 \section{Isotropic Kelvin Material \label{IKM}}
27 As proposed by Kelvin~\cite{Muhlhaus2005} material strain
28 $D_{ij}=\frac{1}{2}(v_{i,j}+v_{j,i})$ can be decomposed into an elastic part
29 $D_{ij}^{el}$ and a visco-plastic part $D_{ij}^{vp}$:
30 \begin{equation}\label{IKM-EQU-2}
31 D_{ij}=D_{ij}^{el}+D_{ij}^{vp}
32 \end{equation}
33 with the elastic strain given as
34 \begin{equation}\label{IKM-EQU-3}
35 D_{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}_{ij}'
36 \end{equation}
37 where $\sigma'_{ij}$ is the deviatoric stress (notice that $\sigma'_{ii}=0$).
38 If the material is composed by materials $q$ the visco-plastic strain can be
39 decomposed as
40 \begin{equation}\label{IKM-EQU-4}
41 D_{ij}^{vp'}=\sum_{q} D_{ij}^{q'}
42 \end{equation}
43 where $D_{ij}^{q}$ is the strain in material $q$ given as
44 \begin{equation}\label{IKM-EQU-5}
45 D_{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'_{ij}
46 \end{equation}
47 and $\eta^{q}$ is the viscosity of material $q$.
48 We assume the following between the strain in material $q$
49 \begin{equation}\label{IKM-EQU-5b}
50 \eta^{q}=\eta^{q}_{N} \left(\frac{\tau}{\tau_{t}^q}\right)^{1-n^{q}}
51 \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'_{ij} \sigma'_{ij}}
52 \end{equation}
53 for given power law coefficients $n^{q}\ge1$ and transition stresses
54 $\tau_{t}^q$, see~\cite{Muhlhaus2005}.
55 Notice that $n^{q}=1$ gives a constant viscosity.
56 After inserting \eqn{IKM-EQU-5} into \eqn{IKM-EQU-4} one gets:
57 \begin{equation}\label{IKM-EQU-6}
58 D_{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'_{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum_{q} \frac{1}{\eta^{q}} \;.
59 \end{equation}
60 and finally with~\ref{IKM-EQU-2}
61 \begin{equation}\label{IKM-EQU-2bb}
62 D_{ij}'=\frac{1}{2 \eta^{vp}} \sigma'_{ij}+\frac{1}{2 \mu} \dot{\sigma}_{ij}'
63 \end{equation}
64 The total stress $\tau$ needs to fulfill the yield condition\index{yield condition}
65 \begin{equation}\label{IKM-EQU-8c}
66 \tau \le \tau_{Y} + \beta \; p
67 \end{equation}
68 with the Drucker-Prager\index{Druck-Prager} cohesion factor\index{cohesion factor}
69 $\tau_{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
70 The deviatoric stress needs to fulfill the equilibrium equation
71 \begin{equation}\label{IKM-EQU-1}
72 -\sigma'_{ij,j}+p_{,i}=F_{i}
73 \end{equation}
74 where $F_{j}$ is a given external force. We assume an incompressible medium:
75 \begin{equation}\label{IKM-EQU-2bbb}
76 -v_{i,i}=0
77 \end{equation}
78 Natural boundary conditions are taken in the form
79 \begin{equation}\label{IKM-EQU-Boundary}
80 \sigma'_{ij}n_{j}-n_{i}p=f
81 \end{equation}
82 which can be overwritten by a constraint
83 \begin{equation}\label{IKM-EQU-Boundary0}
84 v_{i}(x)=0
85 \end{equation}
86 where the index $i$ may depend on the location $x$ on the boundary.
87
88 \subsection{Solution Method \label{IKM-SOLVE}}
89 By using a first order finite difference approximation with step size
90 $dt>0$ \eqn{IKM-EQU-3} is transformed to
91 \begin{equation}\label{IKM-EQU-3b}
92 \dot{\sigma}_{ij}=\frac{1}{dt } \left( \sigma_{ij} - \sigma_{ij}^{-} \right)
93 \end{equation}
94 and
95 \begin{equation}\label{IKM-EQU-2b}
96 D_{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma_{ij}'-\frac{1}{2 \mu dt } \sigma_{ij}^{-'}
97 \end{equation}
98 where $\sigma_{ij}^{-}$ is the stress at the previous time step. With
99 \begin{equation}\label{IKM-EQU-2c}
100 \dot{\gamma} = \sqrt{ 2 \left( D_{ij}' +
101 \frac{1}{ 2 \mu \; dt} \sigma_{ij}^{-'}\right)^2}
102 \end{equation}
103 we have
104 \begin{equation}
105 \tau = \eta_{eff} \cdot \dot{\gamma}
106 \end{equation}
107 where
108 \begin{equation}
109 \eta_{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
110 , \eta_{max}) \mbox{ with }
111 \eta_{max} = \left\{
112 \begin{array}{rcl}
113 \frac{\tau_{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
114 &\mbox{ if } \\
115 \infty & & \mbox{otherwise}
116 \end{array}
117 \right.
118 \end{equation}
119 The upper bound $\eta_{max}$ makes sure that yield condition~\ref{IKM-EQU-8c}
120 holds. With this setting the equation \ref{IKM-EQU-2b} takes the form
121 \begin{equation}\label{IKM-EQU-10}
122 \sigma_{ij}' = 2 \eta_{eff} \left( D_{ij}' +
123 \frac{1}{ 2 \mu \; dt} \sigma_{ij}^{'-}\right)
124 \end{equation}
125 After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
126 \begin{equation}\label{IKM-EQU-1ib}
127 -\left(\eta_{eff} (v_{i,j}+ v_{i,j})
128 \right)_{,j}+p_{,i}=F_{i}+
129 \left(\frac{\eta_{eff}}{\mu dt } \sigma_{ij}^{'-} \right)_{,j}
130 \end{equation}
131 Combining this with the incompressibility condition~\ref{IKM-EQU-2} we need to
132 solve a Stokes problem as discussed in \Sec{STOKES SOLVE} in each time step.
133
134 If we set
135 \begin{equation}\label{IKM-EQU-44}
136 \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
137 \end{equation}
138 we need to solve the nonlinear problem
139 \begin{equation}
140 \eta_{eff} - min(\eta( \dot{\gamma} \cdot \eta_{eff})
141 , \eta_{max}) =0
142 \end{equation}
143 We use the Newton-Raphson scheme\index{Newton-Raphson scheme} to solve this
144 problem:
145 \begin{equation}\label{IKM-EQU-45}
146 \eta_{eff}^{(n+1)} = \min(\eta_{max},
147 \eta_{eff}^{(n)} -
148 \frac{\eta_{eff}^{(n)} - \eta( \tau^{(n)}) }
149 {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
150 =\min(\eta_{max},
151 \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot \eta'( \tau^{(n)} ) }
152 {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
153 \end{equation}
154 where $\eta'$ denotes the derivative of $\eta$ with respect to $\tau$
155 and $\tau^{(n)} = \dot{\gamma} \cdot \eta_{eff}^{(n)}$.
156 Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense to
157 formulate the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
158 In fact we have
159 \begin{equation}
160 \eta' = - \frac{\Theta'}{\Theta^2}
161 \mbox{ with }
162 \Theta' = \sum_{q} \left(\frac{1}{\eta^{q}} \right)'
163 \label{IKM iteration 7}
164 \end{equation}
165 As
166 \begin{equation}\label{IKM-EQU-47}
167 \left(\frac{1}{\eta^{q}} \right)'
168 = \frac{n^{q}-1}{\eta^{q}_{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau_{t}^q)^{n^{q}-1}}
169 = \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau}
170 \end{equation}
171 we have
172 \begin{equation}\label{IKM-EQU-48}
173 \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum_{q}\frac{n^{q}-1}{\eta^{q}}
174 \end{equation}
175 which leads to
176 \begin{equation}\label{IKM-EQU-49}
177 \eta_{eff}^{(n+1)} = \min(\eta_{max},
178 \eta_{eff}^{(n)}
179 \frac{\Theta^{(n)} + \omega^{(n)} }
180 {\eta_{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
181 \end{equation}
182
183 \subsection{Functions}
184
185 \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
186 domain
187 \optional{, stress=0
188 \optional{, v=0
189 \optional{, p=0
190 \optional{, t=0
191 \optional{, numMaterials=1
192 \optional{, verbose=True
193 \optional{, adaptSubTolerance=True
194 }}}}}}}}
195 opens an incompressible, isotropic flow problem in Cartesian coordinates on
196 the domain \var{domain}.
197 \var{stress}, \var{v}, \var{p}, and \var{t} set the initial deviatoric stress,
198 velocity, pressure and time.
199 \var{numMaterials} specifies the number of materials used in the power law
200 model. Some progress information is printed if \var{verbose} is set to \True.
201 If \var{adaptSubTolerance} is equal to \True the tolerances for subproblems
202 are set automatically.
203
204 The domain needs to support LBB compliant elements for the Stokes problem,
205 see~\cite{LBB} for details\index{LBB condition}.
206 For instance one can use second order polynomials for velocity and first order
207 polynomials for the pressure on the same element. Alternatively, one can use
208 macro elements\index{macro elements} using linear polynomials for both
209 pressure and velocity but with a subdivided element for the velocity.
210 Typically, the macro element method is more cost effective.
211 The fact that pressure and velocity are represented in different ways is
212 expressed by
213 \begin{python}
214 velocity=Vector(0.0, Solution(mesh))
215 pressure=Scalar(0.0, ReducedSolution(mesh))
216 \end{python}
217 \end{classdesc}
218
219 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{}
220 returns the domain.
221 \end{methoddesc}
222
223 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
224 returns current time.
225 \end{methoddesc}
226
227 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
228 returns current stress.
229 \end{methoddesc}
230
231 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
232 returns current deviatoric stress.
233 \end{methoddesc}
234
235 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
236 returns current pressure.
237 \end{methoddesc}
238
239 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
240 returns current velocity.
241 \end{methoddesc}
242
243 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
244 returns deviatoric strain of current velocity
245 \end{methoddesc}
246
247 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
248 returns current second invariant of deviatoric stress
249 \end{methoddesc}
250
251 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{}
252 returns current second invariant of deviatoric strain
253 \end{methoddesc}
254
255 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
256 sets the tolerance used to terminate the iteration on a time step.
257 \end{methoddesc}
258
259 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}
260 sets the relative tolerance for the incompressible solver, see
261 \class{StokesProblemCartesian} for details.
262 \end{methoddesc}
263
264 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
265 sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default)
266 elasticity is not applied.
267 \end{methoddesc}
268
269 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}
270 sets the relative tolerance for the effective viscosity. Iteration on a time
271 step is completed if the relative of the effective viscosity is less than
272 \var{rtol}.
273 \end{methoddesc}
274
275 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}
276 {\optional{tau_Y=None, \optional{friction=None}}}
277 sets the parameters $\tau_{Y}$ and $\beta$ for the Drucker-Prager model in
278 condition~\ref{IKM-EQU-8c}. If \var{tau_Y} is set to None (default) then the
279 Drucker-Prager condition is not applied.
280 \end{methoddesc}
281
282 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
283 sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default)
284 elasticity is not applied.
285 \end{methoddesc}
286
287
288 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}
289 sets the parameters of the power-law for all materials as defined in \eqn{IKM-EQU-5b}.
290 \var{eta_N} is the list of viscosities $\eta^{q}_{N}$,
291 \var{tau_t} is the list of reference stresses $\tau_{t}^q$,
292 and \var{power} is the list of power law coefficients $n^{q}$.
293 \end{methoddesc}
294
295 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt
296 \optional{, iter_max=100
297 \optional{, inner_iter_max=20
298 }}}
299 updates stress, velocity and pressure for time increment \var{dt}, where
300 \var{iter_max} is the maximum number of iteration steps on a time step to
301 update the effective viscosity and \var{inner_iter_max} is the maximum
302 number of iteration steps in the incompressible solver.
303 \end{methoddesc}
304
305 %\subsection{Example}
306 %later
307
308 \input{faultsystem}
309
310 %\section{Drucker Prager Model}
311

  ViewVC Help
Powered by ViewVC 1.1.26