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

Revision 3325 - (show annotations)
Fri Oct 29 00:30:51 2010 UTC (11 years, 8 months ago) by caltinay
File MIME type: application/x-tex
File size: 11530 byte(s)
Section 6.1.


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