# Diff of /trunk/doc/user/Models.tex

revision 1709 by lgraham, Fri Aug 15 01:09:45 2008 UTC revision 1966 by jfenwick, Wed Nov 5 04:59:22 2008 UTC
# Line 1  Line 1
1
2    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3  %  %
4  % $Id: Models.tex 1316 2007-09-25 03:18:30Z ksteube$  % Copyright (c) 2003-2008 by University of Queensland
5  %  % Earth Systems Science Computational Center (ESSCC)
6  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  % http://www.uq.edu.au/esscc
%
%           Copyright 2003-2007 by ACceSS MNRF
%       Copyright 2007 by University of Queensland
%
%                http://esscc.uq.edu.au
7  %  %
8  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  % Primary Business: Queensland, Australia
11  %  %
12    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14
15  \chapter{Models}  \chapter{Models}
16
17  The following sections give a breif overview of the model classes and their corresponding methods.  The following sections give a breif overview of the model classes and their corresponding methods.
18
19  \section{Stokes Cartesian (Saddle Point Problem)}  \section{Stokes Problem}
20    The velocity field $v$ and pressure $p$ of an incompressible fluid is given as the solution of their
21  \subsection{Description}  Stokes problem
22    \label{Stokes 1}
23    -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}
24
25    where $\eta$ is the viscosity and $F_i$ defines an internal force. We assume an incompressible media:
26    \label{Stokes 2}
27    -v\hackscore{i,i}=0
28
29    Natural boundary conditions are taken in the form
30    \label{Stokes Boundary}
31    \left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=f
32
33    which can be overwritten by a constraint
34    \label{Stokes Boundary0}
35    v\hackscore{i}(x)=0
36
37    where the index $i$ may depend on the location $x$ on the bondary.
38
Saddle point type problems emerge in a number of applications throughout physics and engineering. Finite element discretisation of the Navier-Stokes (momentum) equations for incompressible flow leads to equations of a saddle point type, which can be formulated as a solution of the following operator problem for $u \in V$ and $p \in Q$ with suitable Hilbert spaces $V$ and $Q$:
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
43  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
44  A     & B \\  A     & B^{*} \\
45  b^{*} & 0 \\  B & 0 \\
46  \end{array} \right]  \end{array} \right]
47  \left[ \begin{array}{c}  \left[ \begin{array}{c}
48  u \\  u \\
49  p \\  p \\
50  \end{array} \right]  \end{array} \right]
51  =\left[ \begin{array}{c}  =\left[ \begin{array}{c}
52  f \\  F \\
53  g \\  0 \\
54  \end{array} \right]  \end{array} \right]
56
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
60    \left[ \begin{array}{cc}
61    A^{-1}     & 0 \\
62    S^{-1} B A^{-1}  & -S^{-1} \\
63    \end{array} \right]
65
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
69    \begin{array}{rcl}
70    A dv & = & F-Av-B^{*}p \\
71    S dp & = & B(v+dv) \\
72    \end{array}
74
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
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}
84
85    To accelerate the convergence we are using the restarted $GMRES$ method using the norm
86
87    \|(v,p)\|^2= \int\hackscore{\Omega} \eta^2 v\hackscore{i,j}^2 + p^2 \; dx
89
90    or alternatively the $PCG$ method on the pressure only using the norm
91
92    \|p\|^2= \int\hackscore{\Omega} \frac{1}{\eta^2} p^2 \; dx
94
95
96
97
where $A$ is coercive, self-adjoint linear operator in $V$, $B$ is a linear operator from $Q$ into $V$ and $B^{*}$ is the adjoint operator of $B$. $f$ and $g$ are given elements from $V$ and $Q$ respectivitly. For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.
98
The Uzawa scheme scheme is used to solve the momentum equation with the secondary condition of incompressibility \cite{GROSS2006,AAMIRBERKYAN2008}.
99
100  \begin{classdesc}{StokesProblemCartesian}{domain,debug}  \begin{classdesc}{StokesProblemCartesian}{domain,debug}
101  opens the stokes equations on the \Domain domain. Setting debug=True switches the debug mode to on.  opens the stokes equations on the \Domain domain. Setting debug=True switches the debug mode to on.
# Line 55  solution.setTolerance(TOL) \\ Line 108  solution.setTolerance(TOL) \\
109  velocity,pressure=solution.solve(velocity,pressure,max\_iter=max\_iter,solver=solver) \\  velocity,pressure=solution.solve(velocity,pressure,max\_iter=max\_iter,solver=solver) \\
110
111  \subsection{Benchmark Problem}  % \subsection{Benchmark Problem}
112    %
113  Convection problem  % Convection problem
114
115
116  \section{Temperature Cartesian}  \section{Temperature Cartesian}
# Line 107  Update level set function; advection and Line 160  Update level set function; advection and
160  Rayleigh-Taylor instability problem  Rayleigh-Taylor instability problem
161
162
163  \section{Drucker Prager Model}  % \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    \label{IKM-EQU-2}
169    D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
170
171    with the elastic strain given as
172    \label{IKM-EQU-3}
173    D\hackscore{ij}'^{el}=\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
174
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    \label{IKM-EQU-4}
178    D\hackscore{ij}'^{vp}=\sum\hackscore{q} D\hackscore{ij}'^{q}
179
180    where $D\hackscore{ij}^{q}$ is the strain in material $q$ given as
181    \label{IKM-EQU-5}
182    D\hackscore{ij}'^{q}=\frac{1}{2 \eta^{q}} \sigma'\hackscore{ij}
183
184    where $\eta^{q}$ is the viscosity of material $q$. We assume the following
185    betwee the the strain in material $q$
186    \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
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    \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
196    With
197    \label{IKM-EQU-8}
198    \dot{\gamma}=\sqrt{2 D\hackscore{ij} D\hackscore{ij}}
199
200    one gets
201    \label{IKM-EQU-8b}
202    \tau = \eta^{vp} \dot{\gamma}^{vp} \;.
203
204    With the Drucker-Prager cohesion factor $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$ we want to achieve
205    \label{IKM-EQU-8c}
206    \tau \le \tau\hackscore{Y} + \beta \; p
207
208    which leads to the condition
209    \label{IKM-EQU-8d}
210    \eta^{vp} \le \frac{\tau\hackscore{Y} + \beta \; p}{ \dot{\gamma}^{vp}} \; .
211
212    Therefore we modify the definition of $\eta^{vp}$ to the form
213    \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
216    The deviatoric stress needs to fullfill the equilibrion equation
217    \label{IKM-EQU-1}
218    -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
219
220    where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:
221    \label{IKM-EQU-2}
222    -v\hackscore{i,i}=0
223
224    Natural boundary conditions are taken in the form
225    \label{IKM-EQU-Boundary}
226    \sigma'\hackscore{ij}n\hackscore{j}-n\hackscore{i}p=f
227
228    which can be overwritten by a constraint
229    \label{IKM-EQU-Boundary0}
230    v\hackscore{i}(x)=0
231
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    \label{IKM-EQU-3b}
237    D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt } \left( \sigma\hackscore{ij}' - \sigma\hackscore{ij}^{'-} \right)
238
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    \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
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    \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
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
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
263    where $v^+=v+dv$. As this problem is non-linear the Jacobi-free Newton-GMRES method is used with the norm
264
265    \|(\sigma, p)\|^2= \int\hackscore{\Omega} \sigma\hackscore{i,j}^2 + p^2 \; dx
266    \label{IKM iteration 3}
267
268
269
\section{Plate Mantel}

Legend:
 Removed from v.1709 changed lines Added in v.1966