1 
lgraham 
2191 

2 
jfenwick 
3989 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
jfenwick 
6651 
% Copyright (c) 20032018 by The University of Queensland 
4 
jfenwick 
3989 
% http://www.uq.edu.au 
5 
caltinay 
3325 
% 
6 


% Primary Business: Queensland, Australia 
7 
jfenwick 
6112 
% Licensed under the Apache License, version 2.0 
8 


% http://www.apache.org/licenses/LICENSE2.0 
9 
caltinay 
3325 
% 
10 
jfenwick 
3989 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
jfenwick 
4657 
% Development 20122013 by School of Earth Sciences 
12 


% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
jfenwick 
3989 
% 
14 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 
caltinay 
3325 

16 
lgraham 
2191 
\section{Stokes Flow} 
17 


\label{STOKES FLOW CHAP} 
18 
acodd 
6928 
In this section we look at Computational Fluid Dynamics (CFD) to simulate 
19 
caltinay 
3288 
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 
lgraham 
2191 
% 
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 
caltinay 
3288 
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 
lgraham 
2191 
% 
44 


\begin{equation} 
45 
acodd 
6928 
(\eta(v_{i,j} + v_{j,i})),_{j} + p,_{i} = f_{i}, 
46 
lgraham 
2191 
\label{GENERAL NAVIER STOKES COM} 
47 


\end{equation} 
48 


% 
49 


with the incompressibility condition 
50 


% 
51 


\begin{equation} 
52 
jfenwick 
3295 
v_{i,i} = 0. 
53 
lgraham 
2191 
\label{INCOMPRESSIBILITY COM} 
54 


\end{equation} 
55 


% 
56 
jfenwick 
3295 
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 
lgraham 
2191 
% 
59 


%\begin{equation} 
60 
jfenwick 
3295 
%\sigma^{'}_{ij} = 2\eta D^{'}_{ij}, 
61 
lgraham 
2191 
%\label{STRESS} 
62 


%\end{equation} 
63 


% 
64 
jfenwick 
3295 
%where the deviatoric stretching $D^{'}_{ij}$ is defined as 
65 
lgraham 
2191 
% 
66 


%\begin{equation} 
67 
jfenwick 
3295 
%D^{'}_{ij} = D^{'}_{ij}  \frac{1}{3}D_{kk}\delta_{ij}. 
68 
lgraham 
2191 
%\label{DEVIATORIC STRETCHING} 
69 


%\end{equation} 
70 


% 
71 
sshaw 
4552 
%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 
caltinay 
3288 
The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in 
73 
jfenwick 
3295 
the $x_{3}$ direction and is given as $f=g\rho\delta_{i3}$. 
74 
caltinay 
3288 
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 
caltinay 
5296 
In order to keep numerical stability and satisfy the CourantFriedrichsLewy condition (CFL condition)\index{Courant number}\index{CFL condition}, the 
79 
caltinay 
3288 
timestep size needs to be kept below a certain value. 
80 
gross 
3379 
The Courant number \index{Courant number} is defined as: 
81 
lgraham 
2191 
% 
82 


\begin{equation} 
83 
caltinay 
3291 
C = \frac{v \delta t}{h} 
84 
lgraham 
2191 
\label{COURANT} 
85 


\end{equation} 
86 


% 
87 
caltinay 
3288 
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 
lgraham 
2191 

92 
caltinay 
3288 
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 
caltinay 
5296 
solver used is PCG (Preconditioned Conjugate Gradients). 
101 


The mesh is defined as a rectangle to represent the body of fluid. 
102 
caltinay 
3288 
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 
caltinay 
3291 
on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}. 
109 
caltinay 
3288 
The fact that pressure and velocity are represented in different ways is 
110 


expressed by 
111 
gross 
2793 
\begin{python} 
112 
caltinay 
3288 
velocity=Vector(0., Solution(mesh)) 
113 


pressure=Scalar(0., ReducedSolution(mesh)) 
114 
gross 
2793 
\end{python} 
115 
caltinay 
3288 
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 
jfenwick 
3295 
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 
caltinay 
3288 
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 
sshaw 
4552 
The timestep size is calculated based on the CourantFriedrichsLewy condition 
128 


(CFL condition)\index{Courant number}\index{CFL condition}, to ensure stable solutions. 
129 
caltinay 
3288 
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 
lgraham 
2191 
% 
134 


\begin{python} 
135 
caltinay 
3288 
from esys.escript import * 
136 


import esys.finley 
137 


from esys.escript.linearPDEs import LinearPDE 
138 


from esys.escript.models import StokesProblemCartesian 
139 
caltinay 
3348 
from esys.weipa import saveVTK 
140 
lgraham 
2191 

141 
caltinay 
3288 
# physical constants 
142 


eta=1. 
143 


rho=100. 
144 


g=10. 
145 
lgraham 
2191 

146 
caltinay 
3288 
# 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 
lgraham 
2191 

154 
caltinay 
3288 
# 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 
lgraham 
2191 

161 
caltinay 
3288 
# gravitational force 
162 


Y=Vector(0., Function(mesh)) 
163 


Y[1] = rho*g 
164 
lgraham 
2191 

165 
caltinay 
3288 
# element spacing 
166 


h = Lsup(mesh.getSize()) 
167 
lgraham 
2191 

168 
caltinay 
3288 
# boundary conditions for slip at base 
169 
jfenwick 
6678 
boundary_cond=whereZero(coordinates[1])*[0.0,1.0] 
170 


+whereZero(coordinates[0])*[1.0,0.0] 
171 
lgraham 
2191 

172 
caltinay 
3288 
# velocity and pressure vectors 
173 


velocity=Vector(0., Solution(mesh)) 
174 


pressure=Scalar(0., ReducedSolution(mesh)) 
175 
lgraham 
2191 

176 
caltinay 
3288 
# Stokes Cartesian 
177 


solution=StokesProblemCartesian(mesh) 
178 


solution.setTolerance(tolerance) 
179 
lgraham 
2191 

180 
caltinay 
3288 
while t <= t_end: 
181 


print("  Time step = %s "%t) 
182 


print("Time = %s seconds"%time) 
183 
lgraham 
2191 

184 
caltinay 
3288 
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 
lgraham 
2191 

188 
caltinay 
3288 
print("Max velocity =", Lsup(velocity), "m/s") 
189 
lgraham 
2191 

190 
gross 
3379 
# CFL condition 
191 
caltinay 
3288 
dt=0.4*h/(Lsup(velocity)) 
192 


print("dt =", dt) 
193 



194 


# displace the mesh 
195 


displacement = velocity * dt 
196 


coordinates = mesh.getX() 
197 
jfenwick 
3914 
newx=interpolate(coordinates + displacement, ContinuousFunction(mesh)) 
198 


mesh.setX(newx) 
199 
caltinay 
3288 

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 
lgraham 
2191 
\end{python} 
208 
lgraham 
2193 
% 
209 
caltinay 
3288 
The results from the simulation can be viewed with \mayavi, by executing the 
210 


following command: 
211 


\begin{verbatim} 
212 
caltinay 
5296 
mayavi2 d vel.00.vtu m Surface 
213 
caltinay 
3288 
\end{verbatim} 
214 
lgraham 
2191 
% 
215 
caltinay 
3288 
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 
lgraham 
2193 
% 
222 
gross 
2654 
\begin{figure}[ht] 
223 
lgraham 
2191 
\center 
224 
caltinay 
3288 
\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 
caltinay 
3279 
%\includegraphics[scale=0.25]{stokesfluidcolorbar} 
235 
caltinay 
3288 
\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 
lgraham 
2191 
\end{figure} 
242 
lgraham 
2193 
% 
243 
caltinay 
3288 
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 
gross 
2371 
%The Level Set Method is discussed in Section \ref{LEVELSET CHAP}. 
249 
caltinay 
3288 
