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} |