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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1878 - (hide annotations)
Tue Oct 14 03:39:13 2008 UTC (11 years, 4 months ago) by gross
File MIME type: application/x-tex
File size: 10554 byte(s)
new version of JacobiFree Newton GMRES + test added.
1 ksteube 1811
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 lgraham 1702 %
4 ksteube 1811 % Copyright (c) 2003-2008 by University of Queensland
5     % Earth Systems Science Computational Center (ESSCC)
6     % http://www.uq.edu.au/esscc
7 lgraham 1702 %
8 ksteube 1811 % 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 lgraham 1702 %
12 ksteube 1811 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13 lgraham 1702
14 ksteube 1811
15 lgraham 1702 \chapter{Models}
16    
17     The following sections give a breif overview of the model classes and their corresponding methods.
18    
19 gross 1878 \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 lgraham 1702
39    
40 gross 1878 \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 lgraham 1702 \begin{equation}
43     \left[ \begin{array}{cc}
44 gross 1878 A & B^{*} \\
45     B & 0 \\
46 lgraham 1702 \end{array} \right]
47     \left[ \begin{array}{c}
48     u \\
49     p \\
50     \end{array} \right]
51     =\left[ \begin{array}{c}
52 gross 1878 F \\
53     0 \\
54 lgraham 1702 \end{array} \right]
55     \label{SADDLEPOINT}
56     \end{equation}
57 gross 1878 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 lgraham 1702
77 gross 1878 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 lgraham 1702
96    
97 gross 1878
98    
99    
100 lgraham 1702 \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 gross 1859 % \subsection{Benchmark Problem}
112     %
113     % Convection problem
114 lgraham 1702
115    
116     \section{Temperature Cartesian}
117    
118     \begin{equation}
119 lgraham 1709 \rho c\hackscore{p} \left (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right ) = k \nabla^{2}T
120 lgraham 1702 \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 gross 1859 % \section{Drucker Prager Model}
164 lgraham 1702
165 gross 1841 \section{Isotropic Kelvin Material \label{IKM}}
166 gross 1859 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 gross 1841 \begin{equation}\label{IKM-EQU-2}
169 gross 1859 D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
170 gross 1841 \end{equation}
171 gross 1859 with the elastic strain given as
172 gross 1841 \begin{equation}\label{IKM-EQU-3}
173 gross 1878 D\hackscore{ij}'^{el}=\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
174 gross 1841 \end{equation}
175 gross 1859 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 gross 1841 \begin{equation}\label{IKM-EQU-4}
178 gross 1859 D\hackscore{ij}'^{vp}=\sum\hackscore{q} D\hackscore{ij}'^{q}
179 gross 1841 \end{equation}
180 gross 1859 where $D\hackscore{ij}^{q}$ is the strain in material $q$ given as
181 gross 1841 \begin{equation}\label{IKM-EQU-5}
182 gross 1859 D\hackscore{ij}'^{q}=\frac{1}{2 \eta^{q}} \sigma'\hackscore{ij}
183 gross 1841 \end{equation}
184 gross 1859 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 gross 1841 After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
193 gross 1859 \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 gross 1841 \end{equation}
196 gross 1859 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 gross 1841 \begin{equation}\label{IKM-EQU-1}
218 gross 1859 -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
219 gross 1841 \end{equation}
220 gross 1859 where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:
221 gross 1841 \begin{equation}\label{IKM-EQU-2}
222 gross 1859 -v\hackscore{i,i}=0
223 gross 1841 \end{equation}
224 gross 1878 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 gross 1859 After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
249     \begin{equation}\label{IKM-EQU-1ib}
250 gross 1878 -\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 gross 1859 \end{equation}
253 gross 1878 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~\label{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 gross 1841
269 gross 1878

  ViewVC Help
Powered by ViewVC 1.1.26