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 |
|