| 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 |
|
Parent Directory
| 
Revision Log