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

Revision 3329 - (show annotations)
Fri Oct 29 05:55:55 2010 UTC (8 years, 10 months ago) by caltinay
File MIME type: application/x-tex
File size: 11498 byte(s)
Faultsystems2d. not finished yet.


 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2010 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 \chapter{Models} 15 \label{MODELS CHAPTER} 16 17 The following sections give a brief overview of the model classes and their 18 corresponding methods. 19 20 \input{stokessolver} 21 %\input{darcyfluxNew} 22 \input{darcyflux} 23 %\input{levelsetmodel} 24 25 \section{Isotropic Kelvin Material \label{IKM}} 26 As proposed by Kelvin~\cite{Muhlhaus2005} material strain 27 $D_{ij}=\frac{1}{2}(v_{i,j}+v_{j,i})$ can be decomposed into an elastic part 28 $D_{ij}^{el}$ and a visco-plastic part $D_{ij}^{vp}$: 29 \begin{equation}\label{IKM-EQU-2} 30 D_{ij}=D_{ij}^{el}+D_{ij}^{vp} 31 \end{equation} 32 with the elastic strain given as 33 \begin{equation}\label{IKM-EQU-3} 34 D_{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}_{ij}' 35 \end{equation} 36 where $\sigma'_{ij}$ is the deviatoric stress (notice that $\sigma'_{ii}=0$). 37 If the material is composed by materials $q$ the visco-plastic strain can be 38 decomposed as 39 \begin{equation}\label{IKM-EQU-4} 40 D_{ij}^{vp'}=\sum_{q} D_{ij}^{q'} 41 \end{equation} 42 where $D_{ij}^{q}$ is the strain in material $q$ given as 43 \begin{equation}\label{IKM-EQU-5} 44 D_{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'_{ij} 45 \end{equation} 46 and $\eta^{q}$ is the viscosity of material $q$. 47 We assume the following between the strain in material $q$ 48 \begin{equation}\label{IKM-EQU-5b} 49 \eta^{q}=\eta^{q}_{N} \left(\frac{\tau}{\tau_{t}^q}\right)^{1-n^{q}} 50 \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'_{ij} \sigma'_{ij}} 51 \end{equation} 52 for given power law coefficients $n^{q}\ge1$ and transition stresses 53 $\tau_{t}^q$, see~\cite{Muhlhaus2005}. 54 Notice that $n^{q}=1$ gives a constant viscosity. 55 After inserting \eqn{IKM-EQU-5} into \eqn{IKM-EQU-4} one gets: 56 \begin{equation}\label{IKM-EQU-6} 57 D_{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'_{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum_{q} \frac{1}{\eta^{q}} \;. 58 \end{equation} 59 and finally with~\ref{IKM-EQU-2} 60 \begin{equation}\label{IKM-EQU-2bb} 61 D_{ij}'=\frac{1}{2 \eta^{vp}} \sigma'_{ij}+\frac{1}{2 \mu} \dot{\sigma}_{ij}' 62 \end{equation} 63 The total stress $\tau$ needs to fulfill the yield condition\index{yield condition} 64 \begin{equation}\label{IKM-EQU-8c} 65 \tau \le \tau_{Y} + \beta \; p 66 \end{equation} 67 with the Drucker-Prager\index{Druck-Prager} cohesion factor\index{cohesion factor} 68 $\tau_{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$. 69 The deviatoric stress needs to fulfill the equilibrium equation 70 \begin{equation}\label{IKM-EQU-1} 71 -\sigma'_{ij,j}+p_{,i}=F_{i} 72 \end{equation} 73 where $F_{j}$ is a given external force. We assume an incompressible medium: 74 \begin{equation}\label{IKM-EQU-2bbb} 75 -v_{i,i}=0 76 \end{equation} 77 Natural boundary conditions are taken in the form 78 \begin{equation}\label{IKM-EQU-Boundary} 79 \sigma'_{ij}n_{j}-n_{i}p=f 80 \end{equation} 81 which can be overwritten by a constraint 82 \begin{equation}\label{IKM-EQU-Boundary0} 83 v_{i}(x)=0 84 \end{equation} 85 where the index $i$ may depend on the location $x$ on the boundary. 86 87 \subsection{Solution Method \label{IKM-SOLVE}} 88 By using a first order finite difference approximation with step size 89 $dt>0$ \eqn{IKM-EQU-3} is transformed to 90 \begin{equation}\label{IKM-EQU-3b} 91 \dot{\sigma}_{ij}=\frac{1}{dt } \left( \sigma_{ij} - \sigma_{ij}^{-} \right) 92 \end{equation} 93 and 94 \begin{equation}\label{IKM-EQU-2b} 95 D_{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma_{ij}'-\frac{1}{2 \mu dt } \sigma_{ij}^{-'} 96 \end{equation} 97 where $\sigma_{ij}^{-}$ is the stress at the previous time step. With 98 \begin{equation}\label{IKM-EQU-2c} 99 \dot{\gamma} = \sqrt{ 2 \left( D_{ij}' + 100 \frac{1}{ 2 \mu \; dt} \sigma_{ij}^{-'}\right)^2} 101 \end{equation} 102 we have 103 \begin{equation} 104 \tau = \eta_{eff} \cdot \dot{\gamma} 105 \end{equation} 106 where 107 \begin{equation} 108 \eta_{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1} 109 , \eta_{max}) \mbox{ with } 110 \eta_{max} = \left\{ 111 \begin{array}{rcl} 112 \frac{\tau_{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\ 113 &\mbox{ if } \\ 114 \infty & & \mbox{otherwise} 115 \end{array} 116 \right. 117 \end{equation} 118 The upper bound $\eta_{max}$ makes sure that yield condition~\ref{IKM-EQU-8c} 119 holds. With this setting the equation \ref{IKM-EQU-2b} takes the form 120 \begin{equation}\label{IKM-EQU-10} 121 \sigma_{ij}' = 2 \eta_{eff} \left( D_{ij}' + 122 \frac{1}{ 2 \mu \; dt} \sigma_{ij}^{'-}\right) 123 \end{equation} 124 After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get 125 \begin{equation}\label{IKM-EQU-1ib} 126 -\left(\eta_{eff} (v_{i,j}+ v_{i,j}) 127 \right)_{,j}+p_{,i}=F_{i}+ 128 \left(\frac{\eta_{eff}}{\mu dt } \sigma_{ij}^{'-} \right)_{,j} 129 \end{equation} 130 Combining this with the incompressibility condition~\ref{IKM-EQU-2} we need to 131 solve a Stokes problem as discussed in \Sec{STOKES SOLVE} in each time step. 132 133 If we set 134 \begin{equation}\label{IKM-EQU-44} 135 \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}} 136 \end{equation} 137 we need to solve the nonlinear problem 138 \begin{equation} 139 \eta_{eff} - min(\eta( \dot{\gamma} \cdot \eta_{eff}) 140 , \eta_{max}) =0 141 \end{equation} 142 We use the Newton-Raphson scheme\index{Newton-Raphson scheme} to solve this 143 problem: 144 \begin{equation}\label{IKM-EQU-45} 145 \eta_{eff}^{(n+1)} = \min(\eta_{max}, 146 \eta_{eff}^{(n)} - 147 \frac{\eta_{eff}^{(n)} - \eta( \tau^{(n)}) } 148 {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} ) 149 =\min(\eta_{max}, 150 \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot \eta'( \tau^{(n)} ) } 151 {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} ) 152 \end{equation} 153 where $\eta'$ denotes the derivative of $\eta$ with respect to $\tau$ 154 and $\tau^{(n)} = \dot{\gamma} \cdot \eta_{eff}^{(n)}$. 155 Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense to 156 formulate the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$. 157 In fact we have 158 \begin{equation} 159 \eta' = - \frac{\Theta'}{\Theta^2} 160 \mbox{ with } 161 \Theta' = \sum_{q} \left(\frac{1}{\eta^{q}} \right)' 162 \label{IKM iteration 7} 163 \end{equation} 164 As 165 \begin{equation}\label{IKM-EQU-47} 166 \left(\frac{1}{\eta^{q}} \right)' 167 = \frac{n^{q}-1}{\eta^{q}_{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau_{t}^q)^{n^{q}-1}} 168 = \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau} 169 \end{equation} 170 we have 171 \begin{equation}\label{IKM-EQU-48} 172 \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum_{q}\frac{n^{q}-1}{\eta^{q}} 173 \end{equation} 174 which leads to 175 \begin{equation}\label{IKM-EQU-49} 176 \eta_{eff}^{(n+1)} = \min(\eta_{max}, 177 \eta_{eff}^{(n)} 178 \frac{\Theta^{(n)} + \omega^{(n)} } 179 {\eta_{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} }) 180 \end{equation} 181 182 \subsection{Functions} 183 184 \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{ 185 domain 186 \optional{, stress=0 187 \optional{, v=0 188 \optional{, p=0 189 \optional{, t=0 190 \optional{, numMaterials=1 191 \optional{, verbose=True 192 \optional{, adaptSubTolerance=True 193 }}}}}}}} 194 opens an incompressible, isotropic flow problem in Cartesian coordinates on 195 the domain \var{domain}. 196 \var{stress}, \var{v}, \var{p}, and \var{t} set the initial deviatoric stress, 197 velocity, pressure and time. 198 \var{numMaterials} specifies the number of materials used in the power law 199 model. Some progress information is printed if \var{verbose} is set to \True. 200 If \var{adaptSubTolerance} is equal to \True the tolerances for subproblems 201 are set automatically. 202 203 The domain needs to support LBB compliant elements for the Stokes problem, 204 see~\cite{LBB} for details\index{LBB condition}. 205 For instance one can use second order polynomials for velocity and first order 206 polynomials for the pressure on the same element. Alternatively, one can use 207 macro elements\index{macro elements} using linear polynomials for both 208 pressure and velocity but with a subdivided element for the velocity. 209 Typically, the macro element method is more cost effective. 210 The fact that pressure and velocity are represented in different ways is 211 expressed by 212 \begin{python} 213 velocity=Vector(0.0, Solution(mesh)) 214 pressure=Scalar(0.0, ReducedSolution(mesh)) 215 \end{python} 216 \end{classdesc} 217 218 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{} 219 returns the domain. 220 \end{methoddesc} 221 222 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{} 223 returns current time. 224 \end{methoddesc} 225 226 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{} 227 returns current stress. 228 \end{methoddesc} 229 230 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{} 231 returns current deviatoric stress. 232 \end{methoddesc} 233 234 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{} 235 returns current pressure. 236 \end{methoddesc} 237 238 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{} 239 returns current velocity. 240 \end{methoddesc} 241 242 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{} 243 returns deviatoric strain of current velocity 244 \end{methoddesc} 245 246 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{} 247 returns current second invariant of deviatoric stress 248 \end{methoddesc} 249 250 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{} 251 returns current second invariant of deviatoric strain 252 \end{methoddesc} 253 254 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4} 255 sets the tolerance used to terminate the iteration on a time step. 256 \end{methoddesc} 257 258 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4} 259 sets the relative tolerance for the incompressible solver, see 260 \class{StokesProblemCartesian} for details. 261 \end{methoddesc} 262 263 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None} 264 sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) 265 elasticity is not applied. 266 \end{methoddesc} 267 268 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8} 269 sets the relative tolerance for the effective viscosity. Iteration on a time 270 step is completed if the relative of the effective viscosity is less than 271 \var{rtol}. 272 \end{methoddesc} 273 274 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw} 275 {\optional{tau_Y=None, \optional{friction=None}}} 276 sets the parameters $\tau_{Y}$ and $\beta$ for the Drucker-Prager model in 277 condition~\ref{IKM-EQU-8c}. If \var{tau_Y} is set to None (default) then the 278 Drucker-Prager condition is not applied. 279 \end{methoddesc} 280 281 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None} 282 sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) 283 elasticity is not applied. 284 \end{methoddesc} 285 286 287 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power} 288 sets the parameters of the power-law for all materials as defined in \eqn{IKM-EQU-5b}. 289 \var{eta_N} is the list of viscosities $\eta^{q}_{N}$, 290 \var{tau_t} is the list of reference stresses $\tau_{t}^q$, 291 and \var{power} is the list of power law coefficients $n^{q}$. 292 \end{methoddesc} 293 294 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt 295 \optional{, iter_max=100 296 \optional{, inner_iter_max=20 297 }}} 298 updates stress, velocity and pressure for time increment \var{dt}, where 299 \var{iter_max} is the maximum number of iteration steps on a time step to 300 update the effective viscosity and \var{inner_iter_max} is the maximum 301 number of iteration steps in the incompressible solver. 302 \end{methoddesc} 303 304 %\subsection{Example} 305 %later 306 307 \input{faultsystem} 308 309 %\section{Drucker Prager Model} 310