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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3071 - (hide annotations)
Wed Jul 21 05:37:30 2010 UTC (11 years, 3 months ago) by gross
File MIME type: application/x-tex
File size: 12306 byte(s)
new darcy flux
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 ksteube 1811
15 lgraham 1702 \chapter{Models}
16 lgraham 2128 \label{MODELS CHAPTER}
17 lgraham 1702
18     The following sections give a breif overview of the model classes and their corresponding methods.
19    
20 gross 2719 \input{stokessolver}
21 gross 3071 \input{darcyfluxNew}
22 lgraham 1702
23 gross 2432 %\input{levelsetmodel}
24 gross 2100
25 lgraham 1702
26    
27 gross 1841 \section{Isotropic Kelvin Material \label{IKM}}
28 gross 2415 As proposed by Kelvin~\cite{Muhlhaus2005} material strain $D\hackscore{ij}=\frac{1}{2}(v\hackscore{i,j}+v\hackscore{j,i})$ can be decomposed into
29 gross 1859 an elastic part $D\hackscore{ij}^{el}$ and visco-plastic part $D\hackscore{ij}^{vp}$:
30 gross 1841 \begin{equation}\label{IKM-EQU-2}
31 gross 1859 D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
32 gross 1841 \end{equation}
33 gross 1859 with the elastic strain given as
34 gross 1841 \begin{equation}\label{IKM-EQU-3}
35 gross 2252 D\hackscore{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
36 gross 1841 \end{equation}
37 gross 1859 where $\sigma'\hackscore{ij}$ is the deviatoric stress (Notice that $\sigma'\hackscore{ii}=0$).
38     If the material is composed by materials $q$ the visco-plastic strain can be decomposed as
39 gross 1841 \begin{equation}\label{IKM-EQU-4}
40 gross 2252 D\hackscore{ij}^{vp'}=\sum\hackscore{q} D\hackscore{ij}^{q'}
41 gross 1841 \end{equation}
42 gross 1859 where $D\hackscore{ij}^{q}$ is the strain in material $q$ given as
43 gross 1841 \begin{equation}\label{IKM-EQU-5}
44 gross 2252 D\hackscore{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'\hackscore{ij}
45 gross 1841 \end{equation}
46 gross 1859 where $\eta^{q}$ is the viscosity of material $q$. We assume the following
47     betwee the the strain in material $q$
48     \begin{equation}\label{IKM-EQU-5b}
49 gross 2438 \eta^{q}=\eta^{q}\hackscore{N} \left(\frac{\tau}{\tau\hackscore{t}^q}\right)^{1-n^{q}}
50 gross 1859 \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'\hackscore{ij} \sigma'\hackscore{ij}}
51     \end{equation}
52 gross 2371 for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau\hackscore{t}^q$, see~\cite{Muhlhaus2005}.
53 gross 1859 Notice that $n^{q}=1$ gives a constant viscosity.
54 gross 1841 After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
55 gross 1859 \begin{equation}\label{IKM-EQU-6}
56 gross 2100 D\hackscore{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum\hackscore{q} \frac{1}{\eta^{q}} \;.
57 gross 1841 \end{equation}
58 gross 2300 and finally with~\ref{IKM-EQU-2}
59 gross 2371 \begin{equation}\label{IKM-EQU-2bb}
60 gross 2300 D\hackscore{ij}'=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij}+\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
61 gross 1859 \end{equation}
62 gross 2300 The total stress $\tau$ needs to fullfill the yield condition \index{yield condition}
63 gross 1859 \begin{equation}\label{IKM-EQU-8c}
64     \tau \le \tau\hackscore{Y} + \beta \; p
65     \end{equation}
66 gross 2300 with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
67 gross 1859 The deviatoric stress needs to fullfill the equilibrion equation
68 gross 1841 \begin{equation}\label{IKM-EQU-1}
69 gross 1859 -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
70 gross 1841 \end{equation}
71 gross 1859 where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:
72 gross 2371 \begin{equation}\label{IKM-EQU-2bbb}
73 gross 1859 -v\hackscore{i,i}=0
74 gross 1841 \end{equation}
75 gross 1878 Natural boundary conditions are taken in the form
76     \begin{equation}\label{IKM-EQU-Boundary}
77     \sigma'\hackscore{ij}n\hackscore{j}-n\hackscore{i}p=f
78     \end{equation}
79     which can be overwritten by a constraint
80     \begin{equation}\label{IKM-EQU-Boundary0}
81     v\hackscore{i}(x)=0
82     \end{equation}
83     where the index $i$ may depend on the location $x$ on the bondary.
84    
85     \subsection{Solution Method \label{IKM-SOLVE}}
86     By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form
87     \begin{equation}\label{IKM-EQU-3b}
88 gross 2300 \dot{\sigma}\hackscore{ij}=\frac{1}{dt } \left( \sigma\hackscore{ij} - \sigma\hackscore{ij}^{-} \right)
89 gross 1878 \end{equation}
90 gross 2300 and
91     \begin{equation}\label{IKM-EQU-2b}
92 gross 2415 D\hackscore{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma\hackscore{ij}'-\frac{1}{2 \mu dt } \sigma\hackscore{ij}^{-'}
93 gross 2300 \end{equation}
94     where $\sigma\hackscore{ij}^{-}$ is the stress at the precious time step. With
95     \begin{equation}\label{IKM-EQU-2c}
96 gross 2415 \dot{\gamma} = \sqrt{ 2 \left( D\hackscore{ij}' +
97     \frac{1}{ 2 \mu \; dt} \sigma\hackscore{ij}^{-'}\right)^2}
98 gross 2300 \end{equation}
99     we have
100     \begin{equation}
101     \tau = \eta\hackscore{eff} \cdot \dot{\gamma}
102     \end{equation}
103     where
104     \begin{equation}
105     \eta\hackscore{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
106     , \eta\hackscore{max}) \mbox{ with }
107     \eta\hackscore{max} = \left\{
108     \begin{array}{rcl}
109     \frac{\tau\hackscore{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
110     &\mbox{ if } \\
111     \infty & & \mbox{otherwise}
112     \end{array}
113     \right.
114     \end{equation}
115     The upper bound $\eta\hackscore{max}$ makes sure that yield condtion~\ref{IKM-EQU-8c} holds. With this setting the eqaution \ref{IKM-EQU-2b} takes the form
116 gross 1878 \begin{equation}\label{IKM-EQU-10}
117 gross 2415 \sigma\hackscore{ij}' = 2 \eta\hackscore{eff} \left( D\hackscore{ij}' +
118 gross 2300 \frac{1}{ 2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)
119 gross 1878 \end{equation}
120 gross 1859 After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
121     \begin{equation}\label{IKM-EQU-1ib}
122 gross 2100 -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
123 gross 2415 \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
124     \left(\frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij}^{'-} \right)\hackscore{,j}
125 gross 1859 \end{equation}
126 gross 2252 Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a
127     Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.
128 gross 2300
129     If we set
130     \begin{equation}\label{IKM-EQU-44}
131     \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
132     \end{equation}
133     we need to solve the nonlinear problem
134 gross 1878 \begin{equation}
135 gross 2300 \eta\hackscore{eff} - min(\eta( \dot{\gamma} \cdot \eta\hackscore{eff})
136     , \eta\hackscore{max}) =0
137     \end{equation}
138     We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
139     \begin{equation}\label{IKM-EQU-45}
140     \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
141     \eta\hackscore{eff}^{(n)} -
142     \frac{\eta\hackscore{eff}^{(n)} - \eta( \tau^{(n)}) }
143     {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
144     =\min(\eta\hackscore{max},
145     \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot \eta'( \tau^{(n)} ) }
146     {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
147 gross 2100 \end{equation}
148 gross 2300 where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
149     and $\tau^{(n)} = \dot{\gamma} \cdot \eta\hackscore{eff}^{(n)}$.
150    
151     Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
152     the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
153     In fact we have
154 gross 2100 \begin{equation}
155 gross 2300 \eta' = - \frac{\Theta'}{\Theta^2}
156 gross 2100 \mbox{ with }
157 gross 2300 \Theta' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'
158 gross 2100 \label{IKM iteration 7}
159     \end{equation}
160 gross 2300 As
161     \begin{equation}\label{IKM-EQU-47}
162 gross 2100 \left(\frac{1}{\eta^{q}} \right)'
163 gross 2438 = \frac{n^{q}-1}{\eta^{q}\hackscore{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau\hackscore{t}^q)^{n^{q}-1}}
164     = \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau}
165 gross 2100 \end{equation}
166 gross 2300 we have
167     \begin{equation}\label{IKM-EQU-48}
168 gross 2438 \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum\hackscore{q}\frac{n^{q}-1}{\eta^{q}}
169 gross 2300 \end{equation}
170     which leads to
171 gross 2371 \begin{equation}\label{IKM-EQU-49}
172 gross 2300 \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
173     \eta\hackscore{eff}^{(n)}
174     \frac{\Theta^{(n)} + \omega^{(n)} }
175     {\eta\hackscore{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
176 gross 2432 \end{equation}
177    
178    
179     \subsection{Functions}
180    
181     \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
182     domain
183     \optional{, stress=0
184     \optional{, v=0
185     \optional{, p=0
186     \optional{, t=0
187     \optional{, numMaterials=1
188 gross 2474 \optional{, verbose=True
189     \optional{, adaptSubTolerance=True
190     }}}}}}}}
191 gross 2433 opens an incompressible, isotropic flow problem in Cartesian cooridninates
192     on the domain \var{domain}.
193 gross 2432 \var{stress},
194     \var{v},
195     \var{p}, and
196     \var{t} set the initial deviatoric stress, velocity, pressure and time.
197     \var{numMaterials} specifies the number of materials used in the power law
198     model. Some progress information are printed if \var{verbose} is set to
199 gross 2474 \True. If \var{adaptSubTolerance} is equal to True the tolerances for subproblems are set automatically.
200 gross 2793
201     The domain
202     needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.
203     For instance one can use second order polynomials for velocity and
204     first order polynomials for the pressure on the same element. Alternativly, one can use
205     macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided
206     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
207     \begin{python}
208     velocity=Vector(0.0, Solution(mesh))
209     pressure=Scalar(0.0, ReducedSolution(mesh))
210     \end{python}
211 gross 2432 \end{classdesc}
212    
213 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{}
214     returns the domain.
215     \end{methoddesc}
216 gross 2432
217 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
218     Returns current time.
219 gross 2432 \end{methoddesc}
220    
221 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
222     Returns current stress.
223 gross 2432 \end{methoddesc}
224    
225 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
226     Returns current deviatoric stress.
227     \end{methoddesc}
228 gross 2432
229 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
230     Returns current pressure.
231 gross 2432 \end{methoddesc}
232 gross 2433
233     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
234     Returns current velocity.
235 gross 2432 \end{methoddesc}
236 gross 2433
237     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
238     Returns deviatoric strain of current velocity
239 gross 2432 \end{methoddesc}
240 gross 2433
241     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
242     Returns current second invariant of deviatoric stress
243 gross 2432 \end{methoddesc}
244 gross 2433
245     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{}
246     Returns current second invariant of deviatoric strain
247 gross 2432 \end{methoddesc}
248 gross 2433
249    
250     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
251     Sets the tolerance used to terminate the iteration on a time step.
252 gross 2432 \end{methoddesc}
253    
254 gross 2433 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}
255     Sets the relative tolerance for the incompressible solver, see \class{StokesProblemCartesian} for details.
256     \end{methoddesc}
257 gross 2432
258 gross 2474 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
259 gross 2433 Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
260     \end{methoddesc}
261    
262     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}
263     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}.
264     \end{methoddesc}
265    
266     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}
267     {\optional{tau_Y=None, \optional{friction=None}}}
268     Sets the parameters $\tau\hackscore{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
269     condition is not applied.
270     \end{methoddesc}
271    
272     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
273     Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
274     \end{methoddesc}
275    
276    
277     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}
278     Sets the parameters of the power-law for all materials as defined in
279     equation~\ref{IKM-EQU-5b}.
280     \var{eta_N} is the list of viscosities $\eta^{q}\hackscore{N}$,
281     \var{tau_t} is the list of reference stresses $\tau\hackscore{t}^q$,
282     and \var{power} is the list of power law coefficients $n^{q}$.
283     \end{methoddesc}
284    
285    
286     \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt
287     \optional{, iter_max=100
288     \optional{, inner_iter_max=20
289     }}}
290     Updates stress, velocity and pressure for time increment \var{dt}.
291     where \var{iter_max} is the maximum number of iteration steps on a time step to
292     update the effective viscosity and \var{inner_iter_max} is the maximum
293     number of itertion steps in the incompressible solver.
294     \end{methoddesc}
295    
296     \subsection{Example}
297     later
298    
299 gross 2647 \input{faultsystem}
300    
301 gross 2432 % \section{Drucker Prager Model}

  ViewVC Help
Powered by ViewVC 1.1.26