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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26