/[escript]/trunk/doc/user/Models.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26