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

Revision 1878 - (show annotations)
Tue Oct 14 03:39:13 2008 UTC (13 years, 8 months ago) by gross
File MIME type: application/x-tex
File size: 10554 byte(s)
new version of JacobiFree Newton GMRES + test added.

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2008 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 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~\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 269