1 
lgraham 
2191 

2 
caltinay 
3325 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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 
lgraham 
2191 
\section{Stokes Flow} 
15 


\label{STOKES FLOW CHAP} 
16 
caltinay 
3288 
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 
lgraham 
2191 
% 
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 
caltinay 
3288 
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 
lgraham 
2191 
% 
42 


\begin{equation} 
43 
jfenwick 
3295 
(\eta(v_{i,j} + v_{j,i})),_{j}  p,_{i} = f_{i}, 
44 
lgraham 
2191 
\label{GENERAL NAVIER STOKES COM} 
45 


\end{equation} 
46 


% 
47 


with the incompressibility condition 
48 


% 
49 


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


\end{equation} 
53 


% 
54 
jfenwick 
3295 
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 
lgraham 
2191 
% 
57 


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


%\end{equation} 
61 


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


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


%\end{equation} 
68 


% 
69 
jfenwick 
3295 
%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 
caltinay 
3288 
The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in 
71 
jfenwick 
3295 
the $x_{3}$ direction and is given as $f=g\rho\delta_{i3}$. 
72 
caltinay 
3288 
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 
lgraham 
2191 
% 
80 


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


\end{equation} 
84 


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

90 
caltinay 
3288 
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 
caltinay 
3291 
on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}. 
107 
caltinay 
3288 
The fact that pressure and velocity are represented in different ways is 
108 


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


pressure=Scalar(0., ReducedSolution(mesh)) 
112 
gross 
2793 
\end{python} 
113 
caltinay 
3288 
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 
jfenwick 
3295 
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 
caltinay 
3288 
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 
lgraham 
2191 
% 
131 


\begin{python} 
132 
caltinay 
3288 
from esys.escript import * 
133 


import esys.finley 
134 


from esys.escript.linearPDEs import LinearPDE 
135 


from esys.escript.models import StokesProblemCartesian 
136 
caltinay 
3348 
from esys.weipa import saveVTK 
137 
lgraham 
2191 

138 
caltinay 
3288 
# physical constants 
139 


eta=1. 
140 


rho=100. 
141 


g=10. 
142 
lgraham 
2191 

143 
caltinay 
3288 
# solver settings 
144 


tolerance=1.0e4 
145 


max_iter=200 
146 


t_end=50 
147 


t=0.0 
148 


time=0 
149 


verbose=True 
150 
lgraham 
2191 

151 
caltinay 
3288 
# define mesh 
152 


H=2. 
153 


L=1. 
154 


W=1. 
155 


mesh = esys.finley.Rectangle(l0=L, l1=H, order=1, n0=20, n1=20) 
156 


coordinates = mesh.getX() 
157 
lgraham 
2191 

158 
caltinay 
3288 
# gravitational force 
159 


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


Y[1] = rho*g 
161 
lgraham 
2191 

162 
caltinay 
3288 
# element spacing 
163 


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

165 
caltinay 
3288 
# boundary conditions for slip at base 
166 


boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0] 
167 
lgraham 
2191 

168 
caltinay 
3288 
# velocity and pressure vectors 
169 


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


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

172 
caltinay 
3288 
# Stokes Cartesian 
173 


solution=StokesProblemCartesian(mesh) 
174 


solution.setTolerance(tolerance) 
175 
lgraham 
2191 

176 
caltinay 
3288 
while t <= t_end: 
177 


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


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

180 
caltinay 
3288 
solution.initialize(fixed_u_mask=boundary_cond, eta=eta, f=Y) 
181 


velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter, \ 
182 


verbose=verbose) 
183 
lgraham 
2191 

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

186 
caltinay 
3288 
# Courant condition 
187 


dt=0.4*h/(Lsup(velocity)) 
188 


print("dt =", dt) 
189 



190 


# displace the mesh 
191 


displacement = velocity * dt 
192 


coordinates = mesh.getX() 
193 


mesh.setX(coordinates + displacement) 
194 



195 


time += dt 
196 



197 


vel_mag = length(velocity) 
198 



199 


#save velocity and pressure output 
200 


saveVTK("vel.%2.2i.vtu"%t, vel=vel_mag, vec=velocity, pressure=pressure) 
201 


t = t+1. 
202 
lgraham 
2191 
\end{python} 
203 
lgraham 
2193 
% 
204 
caltinay 
3288 
The results from the simulation can be viewed with \mayavi, by executing the 
205 


following command: 
206 
lgraham 
2191 
% 
207 
caltinay 
3288 
\begin{verbatim} 
208 


mayavi2 d vel.00.vtu m SurfaceMap 
209 


\end{verbatim} 
210 
lgraham 
2191 
% 
211 
caltinay 
3288 
Colourcoded scalar maps and velocity flow fields can be viewed by selecting 
212 


them in the menu. 
213 


The timesteps can be swept through to view a movie of the simulation. 
214 


\fig{FLUID OUTPUT} shows the simulation output. 
215 


Velocity vectors and a colour map for pressure are shown. 
216 


As the time progresses the body of fluid falls under the influence of gravity. 
217 
lgraham 
2193 
% 
218 
gross 
2654 
\begin{figure}[ht] 
219 
lgraham 
2191 
\center 
220 
caltinay 
3288 
\subfigure[t=1]{\label{FLOW OUTPUT 01}\includegraphics[height=5cm]{stokesfluidt01}} 
221 


\hspace{1.6cm} 
222 


\subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokesfluidt10}} 
223 


\hspace{1.6cm} 
224 


\subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokesfluidt20}}\\ 
225 


\subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokesfluidt30}} 
226 


\hspace{1cm} 
227 


\subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokesfluidt40}} 
228 


\hspace{1cm} 
229 


\subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokesfluidt50}} 
230 
caltinay 
3279 
%\includegraphics[scale=0.25]{stokesfluidcolorbar} 
231 
caltinay 
3288 
\caption{Simulation output for Stokes flow. Fluid body starts off as a 
232 


rectangular shape, then progresses downwards under the influence of gravity. 
233 


Colour coded distribution represents the scalar values for pressure. 
234 


Velocity vectors are displayed at each node in the mesh to show the flow field. 
235 


Computational mesh used was 20$\times$20 elements.} 
236 


\label{FLUID OUTPUT} 
237 
lgraham 
2191 
\end{figure} 
238 
lgraham 
2193 
% 
239 
caltinay 
3288 
The view used here to track the fluid is the Lagrangian view, since the mesh 
240 


moves with the fluid. One of the disadvantages of using the Lagrangian view is 
241 


that the elements in the mesh become severely distorted after a period of time 
242 


and introduce solver errors. To get around this limitation the Level Set 
243 


Method can be used, with the Eulerian point of view for a fixed mesh. 
244 
gross 
2371 
%The Level Set Method is discussed in Section \ref{LEVELSET CHAP}. 
245 
caltinay 
3288 
