1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 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/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 
\section{Stokes Flow} 
15 
\label{STOKES FLOW CHAP} 
16 
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 NavierStokes 
23 
equations can be ignored. 
24 
For a body force $f$, the governing equations are given by: 
25 
% 
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 
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 
% 
42 
\begin{equation} 
43 
(\eta(v_{i,j} + v_{j,i})),_{j}  p,_{i} = f_{i}, 
44 
\label{GENERAL NAVIER STOKES COM} 
45 
\end{equation} 
46 
% 
47 
with the incompressibility condition 
48 
% 
49 
\begin{equation} 
50 
v_{i,i} = 0. 
51 
\label{INCOMPRESSIBILITY COM} 
52 
\end{equation} 
53 
% 
54 
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 
% 
57 
%\begin{equation} 
58 
%\sigma^{'}_{ij} = 2\eta D^{'}_{ij}, 
59 
%\label{STRESS} 
60 
%\end{equation} 
61 
% 
62 
%where the deviatoric stretching $D^{'}_{ij}$ is defined as 
63 
% 
64 
%\begin{equation} 
65 
%D^{'}_{ij} = D^{'}_{ij}  \frac{1}{3}D_{kk}\delta_{ij}. 
66 
%\label{DEVIATORIC STRETCHING} 
67 
%\end{equation} 
68 
% 
69 
%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 
The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in 
71 
the $x_{3}$ direction and is given as $f=g\rho\delta_{i3}$. 
72 
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 
timestep size needs to be kept below a certain value. 
78 
The Courant number is defined as: 
79 
% 
80 
\begin{equation} 
81 
C = \frac{v \delta t}{h} 
82 
\label{COURANT} 
83 
\end{equation} 
84 
% 
85 
where $\delta t$, $v$, and $h$ are the timestep, 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 timestep size was 
88 
calculated for a Courant number of $0.4$. 
89 

90 
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 
on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}. 
107 
The fact that pressure and velocity are represented in different ways is 
108 
expressed by 
109 
\begin{python} 
110 
velocity=Vector(0., Solution(mesh)) 
111 
pressure=Scalar(0., ReducedSolution(mesh)) 
112 
\end{python} 
113 
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 
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 
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 timestep 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 
timestep 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 
% 
131 
\begin{python} 
132 
from esys.escript import * 
133 
import esys.finley 
134 
from esys.escript.linearPDEs import LinearPDE 
135 
from esys.escript.models import StokesProblemCartesian 
136 

137 
# physical constants 
138 
eta=1. 
139 
rho=100. 
140 
g=10. 
141 

142 
# solver settings 
143 
tolerance=1.0e4 
144 
max_iter=200 
145 
t_end=50 
146 
t=0.0 
147 
time=0 
148 
verbose=True 
149 

150 
# define mesh 
151 
H=2. 
152 
L=1. 
153 
W=1. 
154 
mesh = esys.finley.Rectangle(l0=L, l1=H, order=1, n0=20, n1=20) 
155 
coordinates = mesh.getX() 
156 

157 
# gravitational force 
158 
Y=Vector(0., Function(mesh)) 
159 
Y[1] = rho*g 
160 

161 
# element spacing 
162 
h = Lsup(mesh.getSize()) 
163 

164 
# boundary conditions for slip at base 
165 
boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0] 
166 

167 
# velocity and pressure vectors 
168 
velocity=Vector(0., Solution(mesh)) 
169 
pressure=Scalar(0., ReducedSolution(mesh)) 
170 

171 
# Stokes Cartesian 
172 
solution=StokesProblemCartesian(mesh) 
173 
solution.setTolerance(tolerance) 
174 

175 
while t <= t_end: 
176 
print("  Time step = %s "%t) 
177 
print("Time = %s seconds"%time) 
178 

179 
solution.initialize(fixed_u_mask=boundary_cond, eta=eta, f=Y) 
180 
velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter, \ 
181 
verbose=verbose) 
182 

183 
print("Max velocity =", Lsup(velocity), "m/s") 
184 

185 
# Courant condition 
186 
dt=0.4*h/(Lsup(velocity)) 
187 
print("dt =", dt) 
188 

189 
# displace the mesh 
190 
displacement = velocity * dt 
191 
coordinates = mesh.getX() 
192 
mesh.setX(coordinates + displacement) 
193 

194 
time += dt 
195 

196 
vel_mag = length(velocity) 
197 

198 
#save velocity and pressure output 
199 
saveVTK("vel.%2.2i.vtu"%t, vel=vel_mag, vec=velocity, pressure=pressure) 
200 
t = t+1. 
201 
\end{python} 
202 
% 
203 
The results from the simulation can be viewed with \mayavi, by executing the 
204 
following command: 
205 
% 
206 
\begin{verbatim} 
207 
mayavi2 d vel.00.vtu m SurfaceMap 
208 
\end{verbatim} 
209 
% 
210 
Colourcoded scalar maps and velocity flow fields can be viewed by selecting 
211 
them in the menu. 
212 
The timesteps can be swept through to view a movie of the simulation. 
213 
\fig{FLUID OUTPUT} shows the simulation output. 
214 
Velocity vectors and a colour map for pressure are shown. 
215 
As the time progresses the body of fluid falls under the influence of gravity. 
216 
% 
217 
\begin{figure}[ht] 
218 
\center 
219 
\subfigure[t=1]{\label{FLOW OUTPUT 01}\includegraphics[height=5cm]{stokesfluidt01}} 
220 
\hspace{1.6cm} 
221 
\subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokesfluidt10}} 
222 
\hspace{1.6cm} 
223 
\subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokesfluidt20}}\\ 
224 
\subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokesfluidt30}} 
225 
\hspace{1cm} 
226 
\subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokesfluidt40}} 
227 
\hspace{1cm} 
228 
\subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokesfluidt50}} 
229 
%\includegraphics[scale=0.25]{stokesfluidcolorbar} 
230 
\caption{Simulation output for Stokes flow. Fluid body starts off as a 
231 
rectangular shape, then progresses downwards under the influence of gravity. 
232 
Colour coded distribution represents the scalar values for pressure. 
233 
Velocity vectors are displayed at each node in the mesh to show the flow field. 
234 
Computational mesh used was 20$\times$20 elements.} 
235 
\label{FLUID OUTPUT} 
236 
\end{figure} 
237 
% 
238 
The view used here to track the fluid is the Lagrangian view, since the mesh 
239 
moves with the fluid. One of the disadvantages of using the Lagrangian view is 
240 
that the elements in the mesh become severely distorted after a period of time 
241 
and introduce solver errors. To get around this limitation the Level Set 
242 
Method can be used, with the Eulerian point of view for a fixed mesh. 
243 
%The Level Set Method is discussed in Section \ref{LEVELSET CHAP}. 
244 
