| 1 |
|
| 2 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 3 |
% Copyright (c) 2003-2018 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/LICENSE-2.0 |
| 9 |
% |
| 10 |
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) |
| 11 |
% Development 2012-2013 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 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 Navier-Stokes |
| 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 wi-th 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 Courant-Friedrichs-Lewy condition (CFL condition)\index{Courant number}\index{CFL condition}, the |
| 79 |
time-step 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 time-step, 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 time-step 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 time-step size is calculated based on the Courant-Friedrichs-Lewy 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 |
time-step 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.0e-4 |
| 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 |
Colour-coded scalar maps and velocity flow fields can be viewed by selecting |
| 216 |
them in the menu. |
| 217 |
The time-steps 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]{stokes-fluid-t01}} |
| 225 |
\hspace{1.6cm} |
| 226 |
\subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokes-fluid-t10}} |
| 227 |
\hspace{1.6cm} |
| 228 |
\subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokes-fluid-t20}}\\ |
| 229 |
\subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokes-fluid-t30}} |
| 230 |
\hspace{1cm} |
| 231 |
\subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokes-fluid-t40}} |
| 232 |
\hspace{1cm} |
| 233 |
\subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokes-fluid-t50}} |
| 234 |
%\includegraphics[scale=0.25]{stokes-fluid-colorbar} |
| 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 |
|
Parent Directory
| 
Revision Log