1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032018 by The University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Apache License, version 2.0 
8 
% http://www.apache.org/licenses/LICENSE2.0 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development 20122013 by School of Earth Sciences 
12 
% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

16 
\section{Stokes Flow} 
17 
\label{STOKES FLOW CHAP} 
18 
In this section we will look at Computational Fluid Dynamics (CFD) to simulate 
19 
the flow of fluid under the influence of gravity. 
20 
The \class{StokesProblemCartesian} class will be used to calculate the velocity 
21 
and pressure of the fluid. 
22 
The fluid dynamics is governed by the Stokes equation. In geophysical problems 
23 
the velocity of fluids is low; that is, the inertial forces are small compared 
24 
with the viscous forces, therefore the inertial terms in the NavierStokes 
25 
equations can be ignored. 
26 
For a body force $f$, the governing equations are given by: 
27 
% 
28 
\begin{equation} 
29 
\nabla \cdot (\eta(\nabla \vec{v} + \nabla^{T} \vec{v}))  \nabla p = f, 
30 
\label{GENERAL NAVIER STOKES} 
31 
\end{equation} 
32 
% 
33 
with the incompressibility condition 
34 
% 
35 
\begin{equation} 
36 
\nabla \cdot \vec{v} = 0. 
37 
\label{INCOMPRESSIBILITY} 
38 
\end{equation} 
39 
% 
40 
where $p$, $\eta$ and $f$ are the pressure, viscosity and body forces, respectively. 
41 
Alternatively, the Stokes equations can be represented in Einstein summation 
42 
tensor notation (compact notation): 
43 
% 
44 
\begin{equation} 
45 
(\eta(v_{i,j} + v_{j,i})),_{j}  p,_{i} = f_{i}, 
46 
\label{GENERAL NAVIER STOKES COM} 
47 
\end{equation} 
48 
% 
49 
with the incompressibility condition 
50 
% 
51 
\begin{equation} 
52 
v_{i,i} = 0. 
53 
\label{INCOMPRESSIBILITY COM} 
54 
\end{equation} 
55 
% 
56 
The subscript comma $i$ denotes the derivative of the function with respect to $x_{i}$. 
57 
%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}: 
58 
% 
59 
%\begin{equation} 
60 
%\sigma^{'}_{ij} = 2\eta D^{'}_{ij}, 
61 
%\label{STRESS} 
62 
%\end{equation} 
63 
% 
64 
%where the deviatoric stretching $D^{'}_{ij}$ is defined as 
65 
% 
66 
%\begin{equation} 
67 
%D^{'}_{ij} = D^{'}_{ij}  \frac{1}{3}D_{kk}\delta_{ij}. 
68 
%\label{DEVIATORIC STRETCHING} 
69 
%\end{equation} 
70 
% 
71 
%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$). 
72 
The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in 
73 
the $x_{3}$ direction and is given as $f=g\rho\delta_{i3}$. 
74 
The Stokes equation is a saddle point problem, and can be solved using a Uzawa scheme. 
75 
A class called \class{StokesProblemCartesian} in \escript can be used to solve 
76 
for velocity and pressure. A more detailed discussion of the class can be 
77 
found in Chapter \ref{MODELS CHAPTER}. 
78 
In order to keep numerical stability and satisfy the CourantFriedrichsLewy condition (CFL condition)\index{Courant number}\index{CFL condition}, the 
79 
timestep size needs to be kept below a certain value. 
80 
The Courant number \index{Courant number} is defined as: 
81 
% 
82 
\begin{equation} 
83 
C = \frac{v \delta t}{h} 
84 
\label{COURANT} 
85 
\end{equation} 
86 
% 
87 
where $\delta t$, $v$, and $h$ are the timestep, velocity, and the width of 
88 
an element in the mesh, respectively. The velocity $v$ may be chosen as the 
89 
maximum velocity in the domain. In this problem the timestep size was 
90 
calculated for a Courant number of $0.4$. 
91 

92 
The following \PYTHON script is the setup for the Stokes flow simulation, and 
93 
is available in the \ExampleDirectory as \file{fluid.py}. 
94 
It starts off by importing the classes, such as the \class{StokesProblemCartesian} 
95 
class, for solving the Stokes equation and the incompressibility condition for 
96 
velocity and pressure. 
97 
Physical constants are defined for the viscosity and density of the fluid, 
98 
along with the acceleration due to gravity. 
99 
Solver settings are set for the maximum iterations and tolerance; the default 
100 
solver used is PCG (Preconditioned Conjugate Gradients). 
101 
The mesh is defined as a rectangle to represent the body of fluid. 
102 
We are using $20 \times 20$ elements with piecewise linear elements for the 
103 
pressure and for velocity but the elements are subdivided for the velocity. 
104 
This approach is called \textit{macro elements}\index{macro elements} and 
105 
needs to be applied to make sure that the discretized problem has a unique 
106 
solution, see~\cite{LBB} for details\footnote{Alternatively, one can use 
107 
second order elements for the velocity and first order elements for pressure 
108 
on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}. 
109 
The fact that pressure and velocity are represented in different ways is 
110 
expressed by 
111 
\begin{python} 
112 
velocity=Vector(0., Solution(mesh)) 
113 
pressure=Scalar(0., ReducedSolution(mesh)) 
114 
\end{python} 
115 
The gravitational force is calculated based on the fluid density and the 
116 
acceleration due to gravity. 
117 
The boundary conditions are set for a slip condition at the base and the left 
118 
face of the domain. At the base fluid movement in the $x_{0}$direction 
119 
is free, but fixed in the $x_{1}$direction, and similarly at the left 
120 
face fluid movement in the $x_{1}$direction is free but fixed in 
121 
the $x_{0}$direction. 
122 
An instance of the \class{StokesProblemCartesian} class is defined for the 
123 
given computational mesh, and the solver tolerance set. 
124 
Inside the \code{while} loop, the boundary conditions, viscosity and body 
125 
force are initialized. 
126 
The Stokes equation is then solved for velocity and pressure. 
127 
The timestep size is calculated based on the CourantFriedrichsLewy condition 
128 
(CFL condition)\index{Courant number}\index{CFL condition}, to ensure stable solutions. 
129 
The nodes in the mesh are then displaced based on the current velocity and 
130 
timestep size, to move the body of fluid. 
131 
The output for the simulation of velocity and pressure is then saved to a file 
132 
for visualization. 
133 
% 
134 
\begin{python} 
135 
from esys.escript import * 
136 
import esys.finley 
137 
from esys.escript.linearPDEs import LinearPDE 
138 
from esys.escript.models import StokesProblemCartesian 
139 
from esys.weipa import saveVTK 
140 

141 
# physical constants 
142 
eta=1. 
143 
rho=100. 
144 
g=10. 
145 

146 
# solver settings 
147 
tolerance=1.0e4 
148 
max_iter=200 
149 
t_end=50 
150 
t=0.0 
151 
time=0 
152 
verbose=True 
153 

154 
# define mesh 
155 
H=2. 
156 
L=1. 
157 
W=1. 
158 
mesh = esys.finley.Rectangle(l0=L, l1=H, order=1, n0=20, n1=20) 
159 
coordinates = mesh.getX() 
160 

161 
# gravitational force 
162 
Y=Vector(0., Function(mesh)) 
163 
Y[1] = rho*g 
164 

165 
# element spacing 
166 
h = Lsup(mesh.getSize()) 
167 

168 
# boundary conditions for slip at base 
169 
boundary_cond=whereZero(coordinates[1])*[0.0,1.0] 
170 
+whereZero(coordinates[0])*[1.0,0.0] 
171 

172 
# velocity and pressure vectors 
173 
velocity=Vector(0., Solution(mesh)) 
174 
pressure=Scalar(0., ReducedSolution(mesh)) 
175 

176 
# Stokes Cartesian 
177 
solution=StokesProblemCartesian(mesh) 
178 
solution.setTolerance(tolerance) 
179 

180 
while t <= t_end: 
181 
print("  Time step = %s "%t) 
182 
print("Time = %s seconds"%time) 
183 

184 
solution.initialize(fixed_u_mask=boundary_cond, eta=eta, f=Y) 
185 
velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter, \ 
186 
verbose=verbose) 
187 

188 
print("Max velocity =", Lsup(velocity), "m/s") 
189 

190 
# CFL condition 
191 
dt=0.4*h/(Lsup(velocity)) 
192 
print("dt =", dt) 
193 

194 
# displace the mesh 
195 
displacement = velocity * dt 
196 
coordinates = mesh.getX() 
197 
newx=interpolate(coordinates + displacement, ContinuousFunction(mesh)) 
198 
mesh.setX(newx) 
199 

200 
time += dt 
201 

202 
vel_mag = length(velocity) 
203 

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