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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3348 - (show annotations)
Fri Nov 12 05:55:50 2010 UTC (8 years, 10 months ago) by caltinay
File MIME type: application/x-tex
File size: 9394 byte(s)
Added weipa import statements to the relevant parts in user guide and cookbook.
Documented the deprecation of esys.escript.saveVTK in the changes chapter.

1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 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/osl-3.0.php
11 %
12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14 \section{Stokes Flow}
15 \label{STOKES FLOW CHAP}
16 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 Navier-Stokes
23 equations can be ignored.
24 For a body force $f$, the governing equations are given by:
25 %
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 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 %
42 \begin{equation}
43 -(\eta(v_{i,j} + v_{j,i})),_{j} - p,_{i} = f_{i},
44 \label{GENERAL NAVIER STOKES COM}
45 \end{equation}
46 %
47 with the incompressibility condition
48 %
49 \begin{equation}
50 -v_{i,i} = 0.
51 \label{INCOMPRESSIBILITY COM}
52 \end{equation}
53 %
54 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 %
57 %\begin{equation}
58 %\sigma^{'}_{ij} = 2\eta D^{'}_{ij},
59 %\label{STRESS}
60 %\end{equation}
61 %
62 %where the deviatoric stretching $D^{'}_{ij}$ is defined as
63 %
64 %\begin{equation}
65 %D^{'}_{ij} = D^{'}_{ij} - \frac{1}{3}D_{kk}\delta_{ij}.
66 %\label{DEVIATORIC STRETCHING}
67 %\end{equation}
68 %
69 %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 The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in
71 the $x_{3}$ direction and is given as $f=-g\rho\delta_{i3}$.
72 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 time-step size needs to be kept below a certain value.
78 The Courant number is defined as:
79 %
80 \begin{equation}
81 C = \frac{v \delta t}{h}
82 \label{COURANT}
83 \end{equation}
84 %
85 where $\delta t$, $v$, and $h$ are the time-step, 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 time-step size was
88 calculated for a Courant number of $0.4$.
89
90 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 on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}.
107 The fact that pressure and velocity are represented in different ways is
108 expressed by
109 \begin{python}
110 velocity=Vector(0., Solution(mesh))
111 pressure=Scalar(0., ReducedSolution(mesh))
112 \end{python}
113 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 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 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 time-step 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 time-step 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 %
131 \begin{python}
132 from esys.escript import *
133 import esys.finley
134 from esys.escript.linearPDEs import LinearPDE
135 from esys.escript.models import StokesProblemCartesian
136 from esys.weipa import saveVTK
137
138 # physical constants
139 eta=1.
140 rho=100.
141 g=10.
142
143 # solver settings
144 tolerance=1.0e-4
145 max_iter=200
146 t_end=50
147 t=0.0
148 time=0
149 verbose=True
150
151 # 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
158 # gravitational force
159 Y=Vector(0., Function(mesh))
160 Y[1] = -rho*g
161
162 # element spacing
163 h = Lsup(mesh.getSize())
164
165 # boundary conditions for slip at base
166 boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
167
168 # velocity and pressure vectors
169 velocity=Vector(0., Solution(mesh))
170 pressure=Scalar(0., ReducedSolution(mesh))
171
172 # Stokes Cartesian
173 solution=StokesProblemCartesian(mesh)
174 solution.setTolerance(tolerance)
175
176 while t <= t_end:
177 print(" ----- Time step = %s -----"%t)
178 print("Time = %s seconds"%time)
179
180 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
184 print("Max velocity =", Lsup(velocity), "m/s")
185
186 # 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 \end{python}
203 %
204 The results from the simulation can be viewed with \mayavi, by executing the
205 following command:
206 %
207 \begin{verbatim}
208 mayavi2 -d vel.00.vtu -m SurfaceMap
209 \end{verbatim}
210 %
211 Colour-coded scalar maps and velocity flow fields can be viewed by selecting
212 them in the menu.
213 The time-steps 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 %
218 \begin{figure}[ht]
219 \center
220 \subfigure[t=1]{\label{FLOW OUTPUT 01}\includegraphics[height=5cm]{stokes-fluid-t01}}
221 \hspace{1.6cm}
222 \subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokes-fluid-t10}}
223 \hspace{1.6cm}
224 \subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokes-fluid-t20}}\\
225 \subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokes-fluid-t30}}
226 \hspace{1cm}
227 \subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokes-fluid-t40}}
228 \hspace{1cm}
229 \subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokes-fluid-t50}}
230 %\includegraphics[scale=0.25]{stokes-fluid-colorbar}
231 \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 \end{figure}
238 %
239 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 %The Level Set Method is discussed in Section \ref{LEVELSET CHAP}.
245

  ViewVC Help
Powered by ViewVC 1.1.26