1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032012 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/osl3.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 NavierStokes 
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 
timestep 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 timestep, 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 timestep 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 timestep 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 
timestep 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.0e4 
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 
Colourcoded scalar maps and velocity flow fields can be viewed by selecting 
214 
them in the menu. 
215 
The timesteps 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]{stokesfluidt01}} 
223 
\hspace{1.6cm} 
224 
\subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokesfluidt10}} 
225 
\hspace{1.6cm} 
226 
\subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokesfluidt20}}\\ 
227 
\subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokesfluidt30}} 
228 
\hspace{1cm} 
229 
\subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokesfluidt40}} 
230 
\hspace{1cm} 
231 
\subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokesfluidt50}} 
232 
%\includegraphics[scale=0.25]{stokesfluidcolorbar} 
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 
