| 1 |
|
| 2 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 3 |
% |
| 4 |
% Copyright (c) 2003-2010 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 |
\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} |
Parent Directory
| 
Revision Log