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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3348 - (hide annotations)
Fri Nov 12 05:55:50 2010 UTC (8 years, 10 months ago) by caltinay
File MIME type: application/x-tex
File size: 9394 byte(s)
Added weipa import statements to the relevant parts in user guide and cookbook.
Documented the deprecation of esys.escript.saveVTK in the changes chapter.

1 lgraham 2191
2 caltinay 3325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 lgraham 2191 \section{Stokes Flow}
15     \label{STOKES FLOW CHAP}
16 caltinay 3288 In this section we will look at Computational Fluid Dynamics (CFD) to simulate
17     the flow of fluid under the influence of gravity.
18     The \class{StokesProblemCartesian} class will be used to calculate the velocity
19     and pressure of the fluid.
20     The fluid dynamics is governed by the Stokes equation. In geophysical problems
21     the velocity of fluids is low; that is, the inertial forces are small compared
22     with the viscous forces, therefore the inertial terms in the Navier-Stokes
23     equations can be ignored.
24     For a body force $f$, the governing equations are given by:
25 lgraham 2191 %
26     \begin{equation}
27     \nabla \cdot (\eta(\nabla \vec{v} + \nabla^{T} \vec{v})) - \nabla p = -f,
28     \label{GENERAL NAVIER STOKES}
29     \end{equation}
30     %
31     with the incompressibility condition
32     %
33     \begin{equation}
34     \nabla \cdot \vec{v} = 0.
35     \label{INCOMPRESSIBILITY}
36     \end{equation}
37     %
38 caltinay 3288 where $p$, $\eta$ and $f$ are the pressure, viscosity and body forces, respectively.
39     Alternatively, the Stokes equations can be represented in Einstein summation
40     tensor notation (compact notation):
41 lgraham 2191 %
42     \begin{equation}
43 jfenwick 3295 -(\eta(v_{i,j} + v_{j,i})),_{j} - p,_{i} = f_{i},
44 lgraham 2191 \label{GENERAL NAVIER STOKES COM}
45     \end{equation}
46     %
47     with the incompressibility condition
48     %
49     \begin{equation}
50 jfenwick 3295 -v_{i,i} = 0.
51 lgraham 2191 \label{INCOMPRESSIBILITY COM}
52     \end{equation}
53     %
54 jfenwick 3295 The subscript comma $i$ denotes the derivative of the function with respect to $x_{i}$.
55     %A linear relationship between the deviatoric stress $\sigma^{'}_{ij}$ and the stretching $D_{ij} = \frac{1}{2}(v_{i,j} + v_{j,i})$ is defined as \cite{GROSS2006}:
56 lgraham 2191 %
57     %\begin{equation}
58 jfenwick 3295 %\sigma^{'}_{ij} = 2\eta D^{'}_{ij},
59 lgraham 2191 %\label{STRESS}
60     %\end{equation}
61     %
62 jfenwick 3295 %where the deviatoric stretching $D^{'}_{ij}$ is defined as
63 lgraham 2191 %
64     %\begin{equation}
65 jfenwick 3295 %D^{'}_{ij} = D^{'}_{ij} - \frac{1}{3}D_{kk}\delta_{ij}.
66 lgraham 2191 %\label{DEVIATORIC STRETCHING}
67     %\end{equation}
68     %
69 jfenwick 3295 %where $\delta_{ij}$ is the Kronecker $\delta$-symbol, which is a matrix with ones for its diagonal entries ($i = j$) and zeros for the remaining entries ($i \neq j$).
70 caltinay 3288 The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in
71 jfenwick 3295 the $x_{3}$ direction and is given as $f=-g\rho\delta_{i3}$.
72 caltinay 3288 The Stokes equation is a saddle point problem, and can be solved using a Uzawa scheme.
73     A class called \class{StokesProblemCartesian} in \escript can be used to solve
74     for velocity and pressure. A more detailed discussion of the class can be
75     found in Chapter \ref{MODELS CHAPTER}.
76     In order to keep numerical stability and satisfy the Courant condition, the
77     time-step size needs to be kept below a certain value.
78     The Courant number is defined as:
79 lgraham 2191 %
80     \begin{equation}
81 caltinay 3291 C = \frac{v \delta t}{h}
82 lgraham 2191 \label{COURANT}
83     \end{equation}
84     %
85 caltinay 3288 where $\delta t$, $v$, and $h$ are the time-step, velocity, and the width of
86     an element in the mesh, respectively. The velocity $v$ may be chosen as the
87     maximum velocity in the domain. In this problem the time-step size was
88     calculated for a Courant number of $0.4$.
89 lgraham 2191
90 caltinay 3288 The following \PYTHON script is the setup for the Stokes flow simulation, and
91     is available in the \ExampleDirectory as \file{fluid.py}.
92     It starts off by importing the classes, such as the \class{StokesProblemCartesian}
93     class, for solving the Stokes equation and the incompressibility condition for
94     velocity and pressure.
95     Physical constants are defined for the viscosity and density of the fluid,
96     along with the acceleration due to gravity.
97     Solver settings are set for the maximum iterations and tolerance; the default
98     solver used is PCG.
99     The mesh is defined as a rectangle, to represent the body of fluid.
100     We are using $20 \times 20$ elements with piecewise linear elements for the
101     pressure and for velocity but the elements are subdivided for the velocity.
102     This approach is called \textit{macro elements}\index{macro elements} and
103     needs to be applied to make sure that the discretized problem has a unique
104     solution, see~\cite{LBB} for details\footnote{Alternatively, one can use
105     second order elements for the velocity and first order elements for pressure
106 caltinay 3291 on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}.
107 caltinay 3288 The fact that pressure and velocity are represented in different ways is
108     expressed by
109 gross 2793 \begin{python}
110 caltinay 3288 velocity=Vector(0., Solution(mesh))
111     pressure=Scalar(0., ReducedSolution(mesh))
112 gross 2793 \end{python}
113 caltinay 3288 The gravitational force is calculated based on the fluid density and the
114     acceleration due to gravity.
115     The boundary conditions are set for a slip condition at the base and the left
116 jfenwick 3295 face of the domain. At the base fluid movement in the $x_{0}$-direction
117     is free, but fixed in the $x_{1}$-direction, and similarly at the left
118     face fluid movement in the $x_{1}$-direction is free but fixed in
119     the $x_{0}$-direction.
120 caltinay 3288 An instance of the \class{StokesProblemCartesian} class is defined for the
121     given computational mesh, and the solver tolerance set.
122     Inside the \code{while} loop, the boundary conditions, viscosity and body
123     force are initialized.
124     The Stokes equation is then solved for velocity and pressure.
125     The time-step size is calculated based on the Courant condition, to ensure stable solutions.
126     The nodes in the mesh are then displaced based on the current velocity and
127     time-step size, to move the body of fluid.
128     The output for the simulation of velocity and pressure is then saved to a file
129     for visualization.
130 lgraham 2191 %
131     \begin{python}
132 caltinay 3288 from esys.escript import *
133     import esys.finley
134     from esys.escript.linearPDEs import LinearPDE
135     from esys.escript.models import StokesProblemCartesian
136 caltinay 3348 from esys.weipa import saveVTK
137 lgraham 2191
138 caltinay 3288 # physical constants
139     eta=1.
140     rho=100.
141     g=10.
142 lgraham 2191
143 caltinay 3288 # solver settings
144     tolerance=1.0e-4
145     max_iter=200
146     t_end=50
147     t=0.0
148     time=0
149     verbose=True
150 lgraham 2191
151 caltinay 3288 # define mesh
152     H=2.
153     L=1.
154     W=1.
155     mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20)
156     coordinates = mesh.getX()
157 lgraham 2191
158 caltinay 3288 # gravitational force
159     Y=Vector(0., Function(mesh))
160     Y[1] = -rho*g
161 lgraham 2191
162 caltinay 3288 # element spacing
163     h = Lsup(mesh.getSize())
164 lgraham 2191
165 caltinay 3288 # boundary conditions for slip at base
166     boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
167 lgraham 2191
168 caltinay 3288 # velocity and pressure vectors
169     velocity=Vector(0., Solution(mesh))
170     pressure=Scalar(0., ReducedSolution(mesh))
171 lgraham 2191
172 caltinay 3288 # Stokes Cartesian
173     solution=StokesProblemCartesian(mesh)
174     solution.setTolerance(tolerance)
175 lgraham 2191
176 caltinay 3288 while t <= t_end:
177     print(" ----- Time step = %s -----"%t)
178     print("Time = %s seconds"%time)
179 lgraham 2191
180 caltinay 3288 solution.initialize(fixed_u_mask=boundary_cond, eta=eta, f=Y)
181     velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter, \
182     verbose=verbose)
183 lgraham 2191
184 caltinay 3288 print("Max velocity =", Lsup(velocity), "m/s")
185 lgraham 2191
186 caltinay 3288 # Courant condition
187     dt=0.4*h/(Lsup(velocity))
188     print("dt =", dt)
189    
190     # displace the mesh
191     displacement = velocity * dt
192     coordinates = mesh.getX()
193     mesh.setX(coordinates + displacement)
194    
195     time += dt
196    
197     vel_mag = length(velocity)
198    
199     #save velocity and pressure output
200     saveVTK("vel.%2.2i.vtu"%t, vel=vel_mag, vec=velocity, pressure=pressure)
201     t = t+1.
202 lgraham 2191 \end{python}
203 lgraham 2193 %
204 caltinay 3288 The results from the simulation can be viewed with \mayavi, by executing the
205     following command:
206 lgraham 2191 %
207 caltinay 3288 \begin{verbatim}
208     mayavi2 -d vel.00.vtu -m SurfaceMap
209     \end{verbatim}
210 lgraham 2191 %
211 caltinay 3288 Colour-coded scalar maps and velocity flow fields can be viewed by selecting
212     them in the menu.
213     The time-steps can be swept through to view a movie of the simulation.
214     \fig{FLUID OUTPUT} shows the simulation output.
215     Velocity vectors and a colour map for pressure are shown.
216     As the time progresses the body of fluid falls under the influence of gravity.
217 lgraham 2193 %
218 gross 2654 \begin{figure}[ht]
219 lgraham 2191 \center
220 caltinay 3288 \subfigure[t=1]{\label{FLOW OUTPUT 01}\includegraphics[height=5cm]{stokes-fluid-t01}}
221     \hspace{1.6cm}
222     \subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokes-fluid-t10}}
223     \hspace{1.6cm}
224     \subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokes-fluid-t20}}\\
225     \subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokes-fluid-t30}}
226     \hspace{1cm}
227     \subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokes-fluid-t40}}
228     \hspace{1cm}
229     \subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokes-fluid-t50}}
230 caltinay 3279 %\includegraphics[scale=0.25]{stokes-fluid-colorbar}
231 caltinay 3288 \caption{Simulation output for Stokes flow. Fluid body starts off as a
232     rectangular shape, then progresses downwards under the influence of gravity.
233     Colour coded distribution represents the scalar values for pressure.
234     Velocity vectors are displayed at each node in the mesh to show the flow field.
235     Computational mesh used was 20$\times$20 elements.}
236     \label{FLUID OUTPUT}
237 lgraham 2191 \end{figure}
238 lgraham 2193 %
239 caltinay 3288 The view used here to track the fluid is the Lagrangian view, since the mesh
240     moves with the fluid. One of the disadvantages of using the Lagrangian view is
241     that the elements in the mesh become severely distorted after a period of time
242     and introduce solver errors. To get around this limitation the Level Set
243     Method can be used, with the Eulerian point of view for a fixed mesh.
244 gross 2371 %The Level Set Method is discussed in Section \ref{LEVELSET CHAP}.
245 caltinay 3288

  ViewVC Help
Powered by ViewVC 1.1.26