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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3325 - (show annotations)
Fri Oct 29 00:30:51 2010 UTC (7 years, 1 month 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 % 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}

  ViewVC Help
Powered by ViewVC 1.1.26