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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC revision 1966 by jfenwick, Wed Nov 5 04:59:22 2008 UTC
# Line 16  Line 16 
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    \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    
 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  \begin{equation}  \begin{equation}
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]
55  \label{SADDLEPOINT}  \label{SADDLEPOINT}
56  \end{equation}  \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    
 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 53  solution.setTolerance(TOL) \\ Line 108  solution.setTolerance(TOL) \\
108  solution.initialize(fixed\_u\_mask=b\_c,eta=eta,f=Y) \\  solution.initialize(fixed\_u\_mask=b\_c,eta=eta,f=Y) \\
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 105  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    \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    
 \section{Plate Mantel}  

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

  ViewVC Help
Powered by ViewVC 1.1.26