/[escript]/trunk/doc/user/stokesflow.tex
ViewVC logotype

Contents of /trunk/doc/user/stokesflow.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6678 - (show annotations)
Mon Jun 11 04:05:30 2018 UTC (15 months, 1 week ago) by jfenwick
File MIME type: application/x-tex
File size: 9834 byte(s)
Fixing overfull hboxes

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

  ViewVC Help
Powered by ViewVC 1.1.26