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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years ago) by jfenwick
File MIME type: application/x-tex
File size: 11677 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26