# Diff of /trunk/doc/user/Models.tex

revision 3328 by caltinay, Fri Oct 29 00:30:51 2010 UTC revision 3329 by caltinay, Fri Oct 29 05:55:55 2010 UTC
# Line 23  corresponding methods. Line 23  corresponding methods.
23  %\input{levelsetmodel}  %\input{levelsetmodel}
24
25  \section{Isotropic Kelvin Material \label{IKM}}  \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  As proposed by Kelvin~\cite{Muhlhaus2005} material strain
27  an elastic part $D_{ij}^{el}$ and visco-plastic part $D_{ij}^{vp}$:  $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  \label{IKM-EQU-2}  \label{IKM-EQU-2}
30  D_{ij}=D_{ij}^{el}+D_{ij}^{vp}  D_{ij}=D_{ij}^{el}+D_{ij}^{vp}
31
32  with the elastic strain given as  with the elastic strain given as
33  \label{IKM-EQU-3}  \label{IKM-EQU-3}
34  D_{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}_{ij}'  D_{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}_{ij}'
35
36  where $\sigma'_{ij}$ is the deviatoric stress (Notice that $\sigma'_{ii}=0$).  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 decomposed as  If the material is composed by materials $q$ the visco-plastic strain can be
38    decomposed as
39  \label{IKM-EQU-4}  \label{IKM-EQU-4}
40  D_{ij}^{vp'}=\sum_{q} D_{ij}^{q'}  D_{ij}^{vp'}=\sum_{q} D_{ij}^{q'}
41
42  where $D_{ij}^{q}$ is the strain in material $q$ given as  where $D_{ij}^{q}$ is the strain in material $q$ given as
43  \label{IKM-EQU-5}  \label{IKM-EQU-5}
44  D_{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'_{ij}  D_{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'_{ij}
45
46  where $\eta^{q}$ is the viscosity of material $q$. We assume the following  and $\eta^{q}$ is the viscosity of material $q$.
47  betwee the the strain in material $q$  We assume the following between the strain in material $q$
48  \label{IKM-EQU-5b}  \label{IKM-EQU-5b}
49  \eta^{q}=\eta^{q}_{N} \left(\frac{\tau}{\tau_{t}^q}\right)^{1-n^{q}}  \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}}  \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'_{ij} \sigma'_{ij}}
51
52  for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau_{t}^q$, see~\cite{Muhlhaus2005}.  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.  Notice that $n^{q}=1$ gives a constant viscosity.
55  After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:  After inserting \eqn{IKM-EQU-5} into \eqn{IKM-EQU-4} one gets:
56  \label{IKM-EQU-6}  \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}} \;.  D_{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'_{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum_{q} \frac{1}{\eta^{q}} \;.
58
59  and finally with~\ref{IKM-EQU-2}  and finally with~\ref{IKM-EQU-2}
60  \label{IKM-EQU-2bb}  \label{IKM-EQU-2bb}
61  D_{ij}'=\frac{1}{2 \eta^{vp}} \sigma'_{ij}+\frac{1}{2 \mu} \dot{\sigma}_{ij}'  D_{ij}'=\frac{1}{2 \eta^{vp}} \sigma'_{ij}+\frac{1}{2 \mu} \dot{\sigma}_{ij}'
62
63  The total stress $\tau$ needs to fullfill the yield condition \index{yield condition}  The total stress $\tau$ needs to fulfill the yield condition\index{yield condition}
64  \label{IKM-EQU-8c}  \label{IKM-EQU-8c}
65  \tau \le \tau_{Y} + \beta \; p  \tau \le \tau_{Y} + \beta \; p
66
67  with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau_{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.  with the Drucker-Prager\index{Druck-Prager} cohesion factor\index{cohesion factor}
68  The deviatoric stress needs to fullfill the equilibrion equation  $\tau_{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
69    The deviatoric stress needs to fulfill the equilibrium equation
70  \label{IKM-EQU-1}  \label{IKM-EQU-1}
71  -\sigma'_{ij,j}+p_{,i}=F_{i}  -\sigma'_{ij,j}+p_{,i}=F_{i}
72
73  where $F_{j}$ is a given external fource. We assume an incompressible media:  where $F_{j}$ is a given external force. We assume an incompressible medium:
74  \label{IKM-EQU-2bbb}  \label{IKM-EQU-2bbb}
75  -v_{i,i}=0  -v_{i,i}=0
76
77  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
78  \label{IKM-EQU-Boundary}  \label{IKM-EQU-Boundary}
79  \sigma'_{ij}n_{j}-n_{i}p=f  \sigma'_{ij}n_{j}-n_{i}p=f
80
81  which can be overwritten by a constraint  which can be overwritten by a constraint
82  \label{IKM-EQU-Boundary0}  \label{IKM-EQU-Boundary0}
83  v_{i}(x)=0  v_{i}(x)=0
84
85  where the index $i$ may depend on the location $x$ on the bondary.  where the index $i$ may depend on the location $x$ on the boundary.
86
87  \subsection{Solution Method \label{IKM-SOLVE}}  \subsection{Solution Method \label{IKM-SOLVE}}
88  By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form  By using a first order finite difference approximation with step size
89    $dt>0$ \eqn{IKM-EQU-3} is transformed to
90  \label{IKM-EQU-3b}  \label{IKM-EQU-3b}
91  \dot{\sigma}_{ij}=\frac{1}{dt } \left( \sigma_{ij} - \sigma_{ij}^{-} \right)  \dot{\sigma}_{ij}=\frac{1}{dt } \left( \sigma_{ij} - \sigma_{ij}^{-} \right)
92
93  and  and
94  \label{IKM-EQU-2b}  \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}^{-'}  D_{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma_{ij}'-\frac{1}{2 \mu dt } \sigma_{ij}^{-'}
96
97  where $\sigma_{ij}^{-}$ is the stress at the precious time step. With  where $\sigma_{ij}^{-}$ is the stress at the previous time step. With
98  \label{IKM-EQU-2c}  \label{IKM-EQU-2c}
99  \dot{\gamma} = \sqrt{ 2 \left( D_{ij}' +  \dot{\gamma} = \sqrt{ 2 \left( D_{ij}' +
100  \frac{1}{  2 \mu \; dt} \sigma_{ij}^{-'}\right)^2}  \frac{1}{  2 \mu \; dt} \sigma_{ij}^{-'}\right)^2}
# Line 110  where Line 115  where
115  \end{array}  \end{array}
116  \right.  \right.
117
118  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  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  \label{IKM-EQU-10}  \label{IKM-EQU-10}
121  \sigma_{ij}' =  2 \eta_{eff}  \left( D_{ij}' +  \sigma_{ij}' =  2 \eta_{eff}  \left( D_{ij}' +
122  \frac{1}{  2 \mu \; dt} \sigma_{ij}^{'-}\right)    \frac{1}{  2 \mu \; dt} \sigma_{ij}^{'-}\right)
# Line 121  After inserting~\ref{IKM-EQU-10} into~\r Line 127  After inserting~\ref{IKM-EQU-10} into~\r
127  \right)_{,j}+p_{,i}=F_{i}+  \right)_{,j}+p_{,i}=F_{i}+
128   \left(\frac{\eta_{eff}}{\mu dt } \sigma_{ij}^{'-} \right)_{,j}   \left(\frac{\eta_{eff}}{\mu dt } \sigma_{ij}^{'-} \right)_{,j}
129
130  Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a  Combining this with the incompressibility condition~\ref{IKM-EQU-2} we need to
131  Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.  solve a Stokes problem as discussed in \Sec{STOKES SOLVE} in each time step.
132
133  If we set  If we set
134  \label{IKM-EQU-44}  \label{IKM-EQU-44}
# Line 133  we need to solve the nonlinear problem Line 139  we need to solve the nonlinear problem
139  \eta_{eff} -  min(\eta( \dot{\gamma} \cdot \eta_{eff})  \eta_{eff} -  min(\eta( \dot{\gamma} \cdot \eta_{eff})
140  , \eta_{max}) =0  , \eta_{max}) =0
141
142  We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem  We use the Newton-Raphson scheme\index{Newton-Raphson scheme} to solve this
143    problem:
144  \label{IKM-EQU-45}  \label{IKM-EQU-45}
145  \eta_{eff}^{(n+1)} = \min(\eta_{max},  \eta_{eff}^{(n+1)} = \min(\eta_{max},
146  \eta_{eff}^{(n)} -  \eta_{eff}^{(n)} -
# Line 143  We use the Newton-Raphson Scheme \index{ Line 150  We use the Newton-Raphson Scheme \index{
150  \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot  \eta'( \tau^{(n)} )  }  \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot  \eta'( \tau^{(n)} )  }
151  {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )  {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
152
153  where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$  where $\eta'$ denotes the derivative of $\eta$ with respect to $\tau$
154  and $\tau^{(n)} = \dot{\gamma} \cdot \eta_{eff}^{(n)}$.  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  Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate  formulate the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
157  the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.  In fact we have
In fact we have
158
159  \eta' = - \frac{\Theta'}{\Theta^2}  \eta' = - \frac{\Theta'}{\Theta^2}
160  \mbox{ with }  \mbox{ with }
# Line 165  we have Line 171  we have
171  \label{IKM-EQU-48}  \label{IKM-EQU-48}
172  \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum_{q}\frac{n^{q}-1}{\eta^{q}}  \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum_{q}\frac{n^{q}-1}{\eta^{q}}
173
175  \label{IKM-EQU-49}  \label{IKM-EQU-49}
176  \eta_{eff}^{(n+1)} = \min(\eta_{max},  \eta_{eff}^{(n+1)} = \min(\eta_{max},
177  \eta_{eff}^{(n)}  \eta_{eff}^{(n)}
179  {\eta_{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })  {\eta_{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
180
181

182  \subsection{Functions}  \subsection{Functions}
183
184  \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{  \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
# Line 186  domain Line 191  domain
191  \optional{, verbose=True  \optional{, verbose=True
193  }}}}}}}}  }}}}}}}}
194  opens an incompressible, isotropic flow problem in Cartesian cooridninates  opens an incompressible, isotropic flow problem in Cartesian coordinates on
195  on the domain \var{domain}.  the domain \var{domain}.
196  \var{stress},  \var{stress}, \var{v}, \var{p}, and \var{t} set the initial deviatoric stress,
197  \var{v},  velocity, pressure and time.
\var{p}, and
\var{t} set the initial deviatoric stress, velocity, pressure and time.
198  \var{numMaterials} specifies the number of materials used in the power law  \var{numMaterials} specifies the number of materials used in the power law
199  model. Some progress information are printed if \var{verbose} is set to  model. Some progress information is printed if \var{verbose} is set to \True.
200  \True. If \var{adaptSubTolerance} is equal to True the tolerances for subproblems are set automatically.  If \var{adaptSubTolerance} is equal to \True the tolerances for subproblems
201    are set automatically.
202  The domain
203  needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.  The domain needs to support LBB compliant elements for the Stokes problem,
204  For instance one can use second order polynomials for velocity and  see~\cite{LBB} for details\index{LBB condition}.
205  first order polynomials for the pressure on the same element. Alternativly, one can use  For instance one can use second order polynomials for velocity and first order
206  macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided  polynomials for the pressure on the same element. Alternatively, one can use
207  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  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}  \begin{python}
213  velocity=Vector(0.0, Solution(mesh))    velocity=Vector(0.0, Solution(mesh))
214  pressure=Scalar(0.0, ReducedSolution(mesh))    pressure=Scalar(0.0, ReducedSolution(mesh))
215  \end{python}  \end{python}
216  \end{classdesc}  \end{classdesc}
217
# Line 213  returns the domain. Line 220  returns the domain.
220  \end{methoddesc}  \end{methoddesc}
221
222  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
223  Returns current time.  returns current time.
224  \end{methoddesc}  \end{methoddesc}
225
226  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
227  Returns current stress.  returns current stress.
228  \end{methoddesc}  \end{methoddesc}
229
230  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
231  Returns current deviatoric stress.  returns current deviatoric stress.
232  \end{methoddesc}  \end{methoddesc}
233
234  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
235  Returns current pressure.  returns current pressure.
236  \end{methoddesc}  \end{methoddesc}
237
238  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
239  Returns current velocity.  returns current velocity.
240  \end{methoddesc}  \end{methoddesc}
241
242  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
243  Returns deviatoric strain of current velocity  returns deviatoric strain of current velocity
244  \end{methoddesc}  \end{methoddesc}
245
246  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
247  Returns current second invariant of deviatoric stress  returns current second invariant of deviatoric stress
248  \end{methoddesc}  \end{methoddesc}
249
251  Returns current second invariant of deviatoric strain  returns current second invariant of deviatoric strain
252  \end{methoddesc}  \end{methoddesc}
253

254  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
255  Sets the tolerance used to terminate the iteration on a time step.  sets the tolerance used to terminate the iteration on a time step.
256  \end{methoddesc}  \end{methoddesc}
257
258  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}
259  Sets the relative tolerance for the incompressible solver, see \class{StokesProblemCartesian} for details.  sets the relative tolerance for the incompressible solver, see
260    \class{StokesProblemCartesian} for details.
261  \end{methoddesc}  \end{methoddesc}
262
263  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
264  Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.  sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default)
265    elasticity is not applied.
266  \end{methoddesc}  \end{methoddesc}
267
268  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}
269  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}.  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}  \end{methoddesc}
273
274  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}
275  {\optional{tau_Y=None, \optional{friction=None}}}  {\optional{tau_Y=None, \optional{friction=None}}}
276  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  sets the parameters $\tau_{Y}$ and $\beta$ for the Drucker-Prager model in
277  condition is not applied.  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}  \end{methoddesc}
280
281  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
282  Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.  sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default)
283    elasticity is not applied.
284  \end{methoddesc}  \end{methoddesc}
285
286
287  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}
288  Sets the parameters of the power-law for all materials as defined in  sets the parameters of the power-law for all materials as defined in \eqn{IKM-EQU-5b}.
equation~\ref{IKM-EQU-5b}.
289  \var{eta_N} is the list of viscosities $\eta^{q}_{N}$,  \var{eta_N} is the list of viscosities $\eta^{q}_{N}$,
290  \var{tau_t} is the list of reference stresses  $\tau_{t}^q$,  \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}$.  and \var{power} is the list of power law coefficients $n^{q}$.
292  \end{methoddesc}  \end{methoddesc}
293

294  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt  \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt
295  \optional{, iter_max=100  \optional{, iter_max=100
296  \optional{, inner_iter_max=20  \optional{, inner_iter_max=20
297  }}}  }}}
298  Updates stress, velocity and pressure for time increment \var{dt}.  updates stress, velocity and pressure for time increment \var{dt}, where
299  where \var{iter_max} is the maximum number of iteration steps on a time step to  \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  update the effective viscosity and \var{inner_iter_max} is the maximum
301  number of itertion steps in the incompressible solver.  number of iteration steps in the incompressible solver.
302  \end{methoddesc}  \end{methoddesc}
303
304  \subsection{Example}  %\subsection{Example}
305  later  %later
306
307  \input{faultsystem}  \input{faultsystem}
308
309  % \section{Drucker Prager Model}  %\section{Drucker Prager Model}
310

Legend:
 Removed from v.3328 changed lines Added in v.3329