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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1966 - (show annotations)
Wed Nov 5 04:59:22 2008 UTC (11 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 10552 byte(s)
Some more fixes to user guide.
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2008 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
17 The following sections give a breif overview of the model classes and their corresponding methods.
18
19 \section{Stokes Problem}
20 The velocity field $v$ and pressure $p$ of an incompressible fluid is given as the solution of their
21 Stokes problem
22 \begin{equation}\label{Stokes 1}
23 -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}
24 \end{equation}
25 where $\eta$ is the viscosity and $F_i$ defines an internal force. We assume an incompressible media:
26 \begin{equation}\label{Stokes 2}
27 -v\hackscore{i,i}=0
28 \end{equation}
29 Natural boundary conditions are taken in the form
30 \begin{equation}\label{Stokes Boundary}
31 \left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=f
32 \end{equation}
33 which can be overwritten by a constraint
34 \begin{equation}\label{Stokes Boundary0}
35 v\hackscore{i}(x)=0
36 \end{equation}
37 where the index $i$ may depend on the location $x$ on the bondary.
38
39
40 \subsection{Solution Method \label{STOKES SOLVE}}
41 In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem
42 \begin{equation}
43 \left[ \begin{array}{cc}
44 A & B^{*} \\
45 B & 0 \\
46 \end{array} \right]
47 \left[ \begin{array}{c}
48 u \\
49 p \\
50 \end{array} \right]
51 =\left[ \begin{array}{c}
52 F \\
53 0 \\
54 \end{array} \right]
55 \label{SADDLEPOINT}
56 \end{equation}
57 where $A$ is coercive, self-adjoint linear operator in $V$, $B$ is a divergence operator and $B^{*}$ is it adjoint operator (=gradient operator)). For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.
58 We apply the preconditioner matrix
59 \begin{equation}
60 \left[ \begin{array}{cc}
61 A^{-1} & 0 \\
62 S^{-1} B A^{-1} & -S^{-1} \\
63 \end{array} \right]
64 \label{SADDLEPOINT PRECODITIONER}
65 \end{equation}
66 with the Schur complement $S=BA^{-1}B^{*}$ to solve the problem iteratively. The updates $dv$ and $dp$ for
67 given velocity $v$ and pressure $p$ is given as
68 \begin{equation}
69 \begin{array}{rcl}
70 A dv & = & F-Av-B^{*}p \\
71 S dp & = & B(v+dv) \\
72 \end{array}
73 \label{SADDLEPOINT iteration}
74 \end{equation}
75 This scheme is called the Uzawa scheme.
76
77 In the case of the Stokes problem it can be shown that $S$ can be approximated by $\frac{1}{\eta}$. With this the iteration scheme can be implemented as
78 \begin{equation}
79 \begin{array}{rcl}
80 -\left(\eta(dv\hackscore{i,j}+ dv\hackscore{i,j})\right)\hackscore{,j} & = & F\hackscore{i}+\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}-p\hackscore{,i} \\
81 \frac{1}{\eta} dp & = & - \left(v\hackscore{,i}+dv\hackscore{,i} \right)\hackscore{,i} \\
82 \end{array}
83 \label{SADDLEPOINT iteration 2}
84 \end{equation}
85 To accelerate the convergence we are using the restarted $GMRES$ method using the norm
86 \begin{equation}
87 \|(v,p)\|^2= \int\hackscore{\Omega} \eta^2 v\hackscore{i,j}^2 + p^2 \; dx
88 \label{SADDLEPOINT iteration 3}
89 \end{equation}
90 or alternatively the $PCG$ method on the pressure only using the norm
91 \begin{equation}
92 \|p\|^2= \int\hackscore{\Omega} \frac{1}{\eta^2} p^2 \; dx
93 \label{SADDLEPOINT iteration 4}
94 \end{equation}
95
96
97
98
99
100 \begin{classdesc}{StokesProblemCartesian}{domain,debug}
101 opens the stokes equations on the \Domain domain. Setting debug=True switches the debug mode to on.
102 \end{classdesc}
103
104 example usage:
105
106 solution=StokesProblemCartesian(mesh) \\
107 solution.setTolerance(TOL) \\
108 solution.initialize(fixed\_u\_mask=b\_c,eta=eta,f=Y) \\
109 velocity,pressure=solution.solve(velocity,pressure,max\_iter=max\_iter,solver=solver) \\
110
111 % \subsection{Benchmark Problem}
112 %
113 % Convection problem
114
115
116 \section{Temperature Cartesian}
117
118 \begin{equation}
119 \rho c\hackscore{p} \left (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right ) = k \nabla^{2}T
120 \label{HEAT EQUATION}
121 \end{equation}
122
123 where $\vec{v}$ is the velocity vector, $T$ is the temperature, $\rho$ is the density, $\eta$ is the viscosity, $c\hackscore{p}$ is the specific heat at constant pressure and $k$ is the thermal conductivity.
124
125 \subsection{Description}
126
127 \subsection{Method}
128
129 \begin{classdesc}{TemperatureCartesian}{dom,theta=THETA,useSUPG=SUPG}
130 \end{classdesc}
131
132 \subsection{Benchmark Problem}
133
134
135 \section{Level Set Method}
136
137 \subsection{Description}
138
139 \subsection{Method}
140
141 Advection and Reinitialisation
142
143 \begin{classdesc}{LevelSet}{mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth}
144 \end{classdesc}
145
146 %example usage:
147
148 %levelset = LevelSet(mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth)
149
150 \begin{methoddesc}[LevelSet]{update\_parameter}{parameter}
151 Update the parameter.
152 \end{methoddesc}
153
154 \begin{methoddesc}[LevelSet]{update\_phi}{paramter}{velocity}{dt}{t\_step}
155 Update level set function; advection and reinitialization
156 \end{methoddesc}
157
158 \subsection{Benchmark Problem}
159
160 Rayleigh-Taylor instability problem
161
162
163 % \section{Drucker Prager Model}
164
165 \section{Isotropic Kelvin Material \label{IKM}}
166 As proposed by Kelvin~\ref{KELVN} material strain $D\hackscore{ij}=v\hackscore{i,j}+v\hackscore{j,i}$ can be decomposed into
167 an elastic part $D\hackscore{ij}^{el}$ and visco-plastic part $D\hackscore{ij}^{vp}$:
168 \begin{equation}\label{IKM-EQU-2}
169 D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
170 \end{equation}
171 with the elastic strain given as
172 \begin{equation}\label{IKM-EQU-3}
173 D\hackscore{ij}'^{el}=\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
174 \end{equation}
175 where $\sigma'\hackscore{ij}$ is the deviatoric stress (Notice that $\sigma'\hackscore{ii}=0$).
176 If the material is composed by materials $q$ the visco-plastic strain can be decomposed as
177 \begin{equation}\label{IKM-EQU-4}
178 D\hackscore{ij}'^{vp}=\sum\hackscore{q} D\hackscore{ij}'^{q}
179 \end{equation}
180 where $D\hackscore{ij}^{q}$ is the strain in material $q$ given as
181 \begin{equation}\label{IKM-EQU-5}
182 D\hackscore{ij}'^{q}=\frac{1}{2 \eta^{q}} \sigma'\hackscore{ij}
183 \end{equation}
184 where $\eta^{q}$ is the viscosity of material $q$. We assume the following
185 betwee the the strain in material $q$
186 \begin{equation}\label{IKM-EQU-5b}
187 \eta^{q}=\eta^{q}\hackscore{N} \left(\frac{\tau}{\tau\hackscore{t}^q}\right)^{\frac{1}{n^{q}}-1}
188 \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'\hackscore{ij} \sigma'\hackscore{ij}}
189 \end{equation}
190 for a given power law coefficients $n^{q}$ and transition stresses $\tau\hackscore{t}^q$, see~\ref{POERLAW}.
191 Notice that $n^{q}=1$ gives a constant viscosity.
192 After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
193 \begin{equation}\label{IKM-EQU-6}
194 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}} \;.
195 \end{equation}
196 With
197 \begin{equation}\label{IKM-EQU-8}
198 \dot{\gamma}=\sqrt{2 D\hackscore{ij} D\hackscore{ij}}
199 \end{equation}
200 one gets
201 \begin{equation}\label{IKM-EQU-8b}
202 \tau = \eta^{vp} \dot{\gamma}^{vp} \;.
203 \end{equation}
204 With the Drucker-Prager cohesion factor $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$ we want to achieve
205 \begin{equation}\label{IKM-EQU-8c}
206 \tau \le \tau\hackscore{Y} + \beta \; p
207 \end{equation}
208 which leads to the condition
209 \begin{equation}\label{IKM-EQU-8d}
210 \eta^{vp} \le \frac{\tau\hackscore{Y} + \beta \; p}{ \dot{\gamma}^{vp}} \; .
211 \end{equation}
212 Therefore we modify the definition of $\eta^{vp}$ to the form
213 \begin{equation}\label{IKM-EQU-6b}
214 \frac{1}{\eta^{vp}}=\max(\sum\hackscore{q} \frac{1}{\eta^{q}}, \frac{\dot{\gamma}^{vp}} {\tau\hackscore{Y} + \beta \; p})
215 \end{equation}
216 The deviatoric stress needs to fullfill the equilibrion equation
217 \begin{equation}\label{IKM-EQU-1}
218 -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
219 \end{equation}
220 where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:
221 \begin{equation}\label{IKM-EQU-2}
222 -v\hackscore{i,i}=0
223 \end{equation}
224 Natural boundary conditions are taken in the form
225 \begin{equation}\label{IKM-EQU-Boundary}
226 \sigma'\hackscore{ij}n\hackscore{j}-n\hackscore{i}p=f
227 \end{equation}
228 which can be overwritten by a constraint
229 \begin{equation}\label{IKM-EQU-Boundary0}
230 v\hackscore{i}(x)=0
231 \end{equation}
232 where the index $i$ may depend on the location $x$ on the bondary.
233
234 \subsection{Solution Method \label{IKM-SOLVE}}
235 By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form
236 \begin{equation}\label{IKM-EQU-3b}
237 D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt } \left( \sigma\hackscore{ij}' - \sigma\hackscore{ij}^{'-} \right)
238 \end{equation}
239 where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.
240 Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get
241 \begin{equation}\label{IKM-EQU-10}
242 \sigma\hackscore{ij}' = 2 \eta\hackscore{eff} D\hackscore{ij}' +
243 \frac{\eta\hackscore{eff}}{\mu \; dt}
244 \sigma\hackscore{ij}^{'-} \mbox{ with }
245 \frac{1}{\eta\hackscore{eff}}=\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
246 \end{equation}
247 Notice that $\eta\hackscore{eff}$ is a function of stress.
248 After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
249 \begin{equation}\label{IKM-EQU-1ib}
250 -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
251 \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}
252 \end{equation}
253 Together with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a problem with a form almost identical
254 to the Stokes problem discussed in section~\ref{STOKES SOLVE} but with the difference that $\eta\hackscore{eff}$ is depending on the solution. Analog to the iteration scheme~\ref{SADDLEPOINT iteration 2} we can run
255 \begin{equation}
256 \begin{array}{rcl}
257 -\left(\eta\hackscore{eff}(dv\hackscore{i,j}+ dv\hackscore{i,j})\right)\hackscore{,j} & = & F\hackscore{i}+ \sigma\hackscore{ij,j}' \\
258 \frac{1}{\eta\hackscore{eff}} dp & = & - v\hackscore{i,i}^+ \\
259 d\sigma\hackscore{ij}' & = & 2 \eta\hackscore{eff} D\hackscore{ij}^{+'} + \frac{\eta\hackscore{eff}}{\mu \; dt} \sigma\hackscore{ij}' - \sigma\hackscore{ij}\\
260 \end{array}
261 \label{IKM iteration 2}
262 \end{equation}
263 where $v^+=v+dv$. As this problem is non-linear the Jacobi-free Newton-GMRES method is used with the norm
264 \begin{equation}
265 \|(\sigma, p)\|^2= \int\hackscore{\Omega} \sigma\hackscore{i,j}^2 + p^2 \; dx
266 \label{IKM iteration 3}
267 \end{equation}
268
269

  ViewVC Help
Powered by ViewVC 1.1.26