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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26