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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 ksteube 1811
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 lgraham 1702 %
4 jfenwick 2881 % Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 % Earth Systems Science Computational Center (ESSCC)
6     % http://www.uq.edu.au/esscc
7 lgraham 1702 %
8 ksteube 1811 % 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 lgraham 1702 %
12 ksteube 1811 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13 lgraham 1702
14     \chapter{Models}
15 lgraham 2128 \label{MODELS CHAPTER}
16 lgraham 1702
17 caltinay 3325 The following sections give a brief overview of the model classes and their
18     corresponding methods.
19 lgraham 1702
20 gross 2719 \input{stokessolver}
21 gross 3072 %\input{darcyfluxNew}
22     \input{darcyflux}
23 gross 2432 %\input{levelsetmodel}
24 gross 2100
25 gross 1841 \section{Isotropic Kelvin Material \label{IKM}}
26 jfenwick 3295 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 gross 1841 \begin{equation}\label{IKM-EQU-2}
29 jfenwick 3295 D_{ij}=D_{ij}^{el}+D_{ij}^{vp}
30 gross 1841 \end{equation}
31 gross 1859 with the elastic strain given as
32 gross 1841 \begin{equation}\label{IKM-EQU-3}
33 jfenwick 3295 D_{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}_{ij}'
34 gross 1841 \end{equation}
35 jfenwick 3295 where $\sigma'_{ij}$ is the deviatoric stress (Notice that $\sigma'_{ii}=0$).
36 gross 1859 If the material is composed by materials $q$ the visco-plastic strain can be decomposed as
37 gross 1841 \begin{equation}\label{IKM-EQU-4}
38 jfenwick 3295 D_{ij}^{vp'}=\sum_{q} D_{ij}^{q'}
39 gross 1841 \end{equation}
40 jfenwick 3295 where $D_{ij}^{q}$ is the strain in material $q$ given as
41 gross 1841 \begin{equation}\label{IKM-EQU-5}
42 jfenwick 3295 D_{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'_{ij}
43 gross 1841 \end{equation}
44 gross 1859 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 jfenwick 3295 \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 gross 1859 \end{equation}
50 jfenwick 3295 for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau_{t}^q$, see~\cite{Muhlhaus2005}.
51 gross 1859 Notice that $n^{q}=1$ gives a constant viscosity.
52 gross 1841 After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
53 gross 1859 \begin{equation}\label{IKM-EQU-6}
54 jfenwick 3295 D_{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'_{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum_{q} \frac{1}{\eta^{q}} \;.
55 gross 1841 \end{equation}
56 gross 2300 and finally with~\ref{IKM-EQU-2}
57 gross 2371 \begin{equation}\label{IKM-EQU-2bb}
58 jfenwick 3295 D_{ij}'=\frac{1}{2 \eta^{vp}} \sigma'_{ij}+\frac{1}{2 \mu} \dot{\sigma}_{ij}'
59 gross 1859 \end{equation}
60 gross 2300 The total stress $\tau$ needs to fullfill the yield condition \index{yield condition}
61 gross 1859 \begin{equation}\label{IKM-EQU-8c}
62 jfenwick 3295 \tau \le \tau_{Y} + \beta \; p
63 gross 1859 \end{equation}
64 jfenwick 3295 with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau_{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
65 gross 1859 The deviatoric stress needs to fullfill the equilibrion equation
66 gross 1841 \begin{equation}\label{IKM-EQU-1}
67 jfenwick 3295 -\sigma'_{ij,j}+p_{,i}=F_{i}
68 gross 1841 \end{equation}
69 jfenwick 3295 where $F_{j}$ is a given external fource. We assume an incompressible media:
70 gross 2371 \begin{equation}\label{IKM-EQU-2bbb}
71 jfenwick 3295 -v_{i,i}=0
72 gross 1841 \end{equation}
73 gross 1878 Natural boundary conditions are taken in the form
74     \begin{equation}\label{IKM-EQU-Boundary}
75 jfenwick 3295 \sigma'_{ij}n_{j}-n_{i}p=f
76 gross 1878 \end{equation}
77     which can be overwritten by a constraint
78     \begin{equation}\label{IKM-EQU-Boundary0}
79 jfenwick 3295 v_{i}(x)=0
80 gross 1878 \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 jfenwick 3295 \dot{\sigma}_{ij}=\frac{1}{dt } \left( \sigma_{ij} - \sigma_{ij}^{-} \right)
87 gross 1878 \end{equation}
88 gross 2300 and
89     \begin{equation}\label{IKM-EQU-2b}
90 jfenwick 3295 D_{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma_{ij}'-\frac{1}{2 \mu dt } \sigma_{ij}^{-'}
91 gross 2300 \end{equation}
92 jfenwick 3295 where $\sigma_{ij}^{-}$ is the stress at the precious time step. With
93 gross 2300 \begin{equation}\label{IKM-EQU-2c}
94 jfenwick 3295 \dot{\gamma} = \sqrt{ 2 \left( D_{ij}' +
95     \frac{1}{ 2 \mu \; dt} \sigma_{ij}^{-'}\right)^2}
96 gross 2300 \end{equation}
97     we have
98     \begin{equation}
99 jfenwick 3295 \tau = \eta_{eff} \cdot \dot{\gamma}
100 gross 2300 \end{equation}
101     where
102     \begin{equation}
103 jfenwick 3295 \eta_{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
104     , \eta_{max}) \mbox{ with }
105     \eta_{max} = \left\{
106 gross 2300 \begin{array}{rcl}
107 jfenwick 3295 \frac{\tau_{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
108 gross 2300 &\mbox{ if } \\
109     \infty & & \mbox{otherwise}
110     \end{array}
111     \right.
112     \end{equation}
113 jfenwick 3295 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 gross 1878 \begin{equation}\label{IKM-EQU-10}
115 jfenwick 3295 \sigma_{ij}' = 2 \eta_{eff} \left( D_{ij}' +
116     \frac{1}{ 2 \mu \; dt} \sigma_{ij}^{'-}\right)
117 gross 1878 \end{equation}
118 gross 1859 After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
119     \begin{equation}\label{IKM-EQU-1ib}
120 jfenwick 3295 -\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 gross 1859 \end{equation}
124 gross 2252 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 gross 2300
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 gross 1878 \begin{equation}
133 jfenwick 3295 \eta_{eff} - min(\eta( \dot{\gamma} \cdot \eta_{eff})
134     , \eta_{max}) =0
135 gross 2300 \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 jfenwick 3295 \eta_{eff}^{(n+1)} = \min(\eta_{max},
139     \eta_{eff}^{(n)} -
140     \frac{\eta_{eff}^{(n)} - \eta( \tau^{(n)}) }
141 gross 2300 {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
142 jfenwick 3295 =\min(\eta_{max},
143 gross 2300 \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot \eta'( \tau^{(n)} ) }
144     {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
145 gross 2100 \end{equation}
146 gross 2300 where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
147 jfenwick 3295 and $\tau^{(n)} = \dot{\gamma} \cdot \eta_{eff}^{(n)}$.
148 gross 2300
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 gross 2100 \begin{equation}
153 gross 2300 \eta' = - \frac{\Theta'}{\Theta^2}
154 gross 2100 \mbox{ with }
155 jfenwick 3295 \Theta' = \sum_{q} \left(\frac{1}{\eta^{q}} \right)'
156 gross 2100 \label{IKM iteration 7}
157     \end{equation}
158 gross 2300 As
159     \begin{equation}\label{IKM-EQU-47}
160 gross 2100 \left(\frac{1}{\eta^{q}} \right)'
161 jfenwick 3295 = \frac{n^{q}-1}{\eta^{q}_{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau_{t}^q)^{n^{q}-1}}
162 gross 2438 = \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau}
163 gross 2100 \end{equation}
164 gross 2300 we have
165     \begin{equation}\label{IKM-EQU-48}
166 jfenwick 3295 \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum_{q}\frac{n^{q}-1}{\eta^{q}}
167 gross 2300 \end{equation}
168     which leads to
169 gross 2371 \begin{equation}\label{IKM-EQU-49}
170 jfenwick 3295 \eta_{eff}^{(n+1)} = \min(\eta_{max},
171     \eta_{eff}^{(n)}
172 gross 2300 \frac{\Theta^{(n)} + \omega^{(n)} }
173 jfenwick 3295 {\eta_{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
174 gross 2432 \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 gross 2474 \optional{, verbose=True
187     \optional{, adaptSubTolerance=True
188     }}}}}}}}
189 gross 2433 opens an incompressible, isotropic flow problem in Cartesian cooridninates
190     on the domain \var{domain}.
191 gross 2432 \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 gross 2474 \True. If \var{adaptSubTolerance} is equal to True the tolerances for subproblems are set automatically.
198 gross 2793
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 gross 2432 \end{classdesc}
210    
211 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{}
212     returns the domain.
213     \end{methoddesc}
214 gross 2432
215 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
216     Returns current time.
217 gross 2432 \end{methoddesc}
218    
219 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
220     Returns current stress.
221 gross 2432 \end{methoddesc}
222    
223 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
224     Returns current deviatoric stress.
225     \end{methoddesc}
226 gross 2432
227 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
228     Returns current pressure.
229 gross 2432 \end{methoddesc}
230 gross 2433
231     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
232     Returns current velocity.
233 gross 2432 \end{methoddesc}
234 gross 2433
235     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
236     Returns deviatoric strain of current velocity
237 gross 2432 \end{methoddesc}
238 gross 2433
239     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
240     Returns current second invariant of deviatoric stress
241 gross 2432 \end{methoddesc}
242 gross 2433
243     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{}
244     Returns current second invariant of deviatoric strain
245 gross 2432 \end{methoddesc}
246 gross 2433
247    
248     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
249     Sets the tolerance used to terminate the iteration on a time step.
250 gross 2432 \end{methoddesc}
251    
252 gross 2433 \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 gross 2432
256 gross 2474 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
257 gross 2433 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 jfenwick 3295 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 gross 2433 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 jfenwick 3295 \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 gross 2433 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 gross 2647 \input{faultsystem}
298    
299 gross 3072 % \section{Drucker Prager Model}

  ViewVC Help
Powered by ViewVC 1.1.26