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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 1 month ago) by jfenwick
File MIME type: application/x-tex
File size: 11677 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26