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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3989 - (hide annotations)
Tue Sep 25 02:21:54 2012 UTC (6 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 9745 byte(s)
More copyright fixes.
pyvisi traces removed.
Some install doco
1 lgraham 2191
2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 3911 % Copyright (c) 2003-2012 by University of Queensland
4 jfenwick 3989 % http://www.uq.edu.au
5 caltinay 3325 %
6     % Primary Business: Queensland, Australia
7     % Licensed under the Open Software License version 3.0
8     % http://www.opensource.org/licenses/osl-3.0.php
9     %
10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11     % Development since 2012 by School of Earth Sciences
12     %
13     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14 caltinay 3325
15 lgraham 2191 \section{Stokes Flow}
16     \label{STOKES FLOW CHAP}
17 caltinay 3288 In this section we will look at Computational Fluid Dynamics (CFD) to simulate
18     the flow of fluid under the influence of gravity.
19     The \class{StokesProblemCartesian} class will be used to calculate the velocity
20     and pressure of the fluid.
21     The fluid dynamics is governed by the Stokes equation. In geophysical problems
22     the velocity of fluids is low; that is, the inertial forces are small compared
23     with the viscous forces, therefore the inertial terms in the Navier-Stokes
24     equations can be ignored.
25     For a body force $f$, the governing equations are given by:
26 lgraham 2191 %
27     \begin{equation}
28     \nabla \cdot (\eta(\nabla \vec{v} + \nabla^{T} \vec{v})) - \nabla p = -f,
29     \label{GENERAL NAVIER STOKES}
30     \end{equation}
31     %
32     with the incompressibility condition
33     %
34     \begin{equation}
35     \nabla \cdot \vec{v} = 0.
36     \label{INCOMPRESSIBILITY}
37     \end{equation}
38     %
39 caltinay 3288 where $p$, $\eta$ and $f$ are the pressure, viscosity and body forces, respectively.
40     Alternatively, the Stokes equations can be represented in Einstein summation
41     tensor notation (compact notation):
42 lgraham 2191 %
43     \begin{equation}
44 jfenwick 3295 -(\eta(v_{i,j} + v_{j,i})),_{j} - p,_{i} = f_{i},
45 lgraham 2191 \label{GENERAL NAVIER STOKES COM}
46     \end{equation}
47     %
48     with the incompressibility condition
49     %
50     \begin{equation}
51 jfenwick 3295 -v_{i,i} = 0.
52 lgraham 2191 \label{INCOMPRESSIBILITY COM}
53     \end{equation}
54     %
55 jfenwick 3295 The subscript comma $i$ denotes the derivative of the function with respect to $x_{i}$.
56     %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}:
57 lgraham 2191 %
58     %\begin{equation}
59 jfenwick 3295 %\sigma^{'}_{ij} = 2\eta D^{'}_{ij},
60 lgraham 2191 %\label{STRESS}
61     %\end{equation}
62     %
63 jfenwick 3295 %where the deviatoric stretching $D^{'}_{ij}$ is defined as
64 lgraham 2191 %
65     %\begin{equation}
66 jfenwick 3295 %D^{'}_{ij} = D^{'}_{ij} - \frac{1}{3}D_{kk}\delta_{ij}.
67 lgraham 2191 %\label{DEVIATORIC STRETCHING}
68     %\end{equation}
69     %
70 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$).
71 caltinay 3288 The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in
72 jfenwick 3295 the $x_{3}$ direction and is given as $f=-g\rho\delta_{i3}$.
73 caltinay 3288 The Stokes equation is a saddle point problem, and can be solved using a Uzawa scheme.
74     A class called \class{StokesProblemCartesian} in \escript can be used to solve
75     for velocity and pressure. A more detailed discussion of the class can be
76     found in Chapter \ref{MODELS CHAPTER}.
77 gross 3379 In order to keep numerical stability and satisfy the Courant–Friedrichs–Lewy condition (CFL condition) \index{Courant number}\index{CFL condition}, the
78 caltinay 3288 time-step size needs to be kept below a certain value.
79 gross 3379 The Courant number \index{Courant number} is defined as:
80 lgraham 2191 %
81     \begin{equation}
82 caltinay 3291 C = \frac{v \delta t}{h}
83 lgraham 2191 \label{COURANT}
84     \end{equation}
85     %
86 caltinay 3288 where $\delta t$, $v$, and $h$ are the time-step, velocity, and the width of
87     an element in the mesh, respectively. The velocity $v$ may be chosen as the
88     maximum velocity in the domain. In this problem the time-step size was
89     calculated for a Courant number of $0.4$.
90 lgraham 2191
91 caltinay 3288 The following \PYTHON script is the setup for the Stokes flow simulation, and
92     is available in the \ExampleDirectory as \file{fluid.py}.
93     It starts off by importing the classes, such as the \class{StokesProblemCartesian}
94     class, for solving the Stokes equation and the incompressibility condition for
95     velocity and pressure.
96     Physical constants are defined for the viscosity and density of the fluid,
97     along with the acceleration due to gravity.
98     Solver settings are set for the maximum iterations and tolerance; the default
99     solver used is PCG.
100     The mesh is defined as a rectangle, to represent the body of fluid.
101     We are using $20 \times 20$ elements with piecewise linear elements for the
102     pressure and for velocity but the elements are subdivided for the velocity.
103     This approach is called \textit{macro elements}\index{macro elements} and
104     needs to be applied to make sure that the discretized problem has a unique
105     solution, see~\cite{LBB} for details\footnote{Alternatively, one can use
106     second order elements for the velocity and first order elements for pressure
107 caltinay 3291 on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}.
108 caltinay 3288 The fact that pressure and velocity are represented in different ways is
109     expressed by
110 gross 2793 \begin{python}
111 caltinay 3288 velocity=Vector(0., Solution(mesh))
112     pressure=Scalar(0., ReducedSolution(mesh))
113 gross 2793 \end{python}
114 caltinay 3288 The gravitational force is calculated based on the fluid density and the
115     acceleration due to gravity.
116     The boundary conditions are set for a slip condition at the base and the left
117 jfenwick 3295 face of the domain. At the base fluid movement in the $x_{0}$-direction
118     is free, but fixed in the $x_{1}$-direction, and similarly at the left
119     face fluid movement in the $x_{1}$-direction is free but fixed in
120     the $x_{0}$-direction.
121 caltinay 3288 An instance of the \class{StokesProblemCartesian} class is defined for the
122     given computational mesh, and the solver tolerance set.
123     Inside the \code{while} loop, the boundary conditions, viscosity and body
124     force are initialized.
125     The Stokes equation is then solved for velocity and pressure.
126 gross 3379 The time-step size is calculated based on the Courant–Friedrichs–Lewy condition (CFL condition) \index{Courant number}\index{CFL condition}, to ensure stable solutions.
127 caltinay 3288 The nodes in the mesh are then displaced based on the current velocity and
128     time-step size, to move the body of fluid.
129     The output for the simulation of velocity and pressure is then saved to a file
130     for visualization.
131 lgraham 2191 %
132     \begin{python}
133 caltinay 3288 from esys.escript import *
134     import esys.finley
135     from esys.escript.linearPDEs import LinearPDE
136     from esys.escript.models import StokesProblemCartesian
137 caltinay 3348 from esys.weipa import saveVTK
138 lgraham 2191
139 caltinay 3288 # physical constants
140     eta=1.
141     rho=100.
142     g=10.
143 lgraham 2191
144 caltinay 3288 # solver settings
145     tolerance=1.0e-4
146     max_iter=200
147     t_end=50
148     t=0.0
149     time=0
150     verbose=True
151 lgraham 2191
152 caltinay 3288 # define mesh
153     H=2.
154     L=1.
155     W=1.
156     mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20)
157     coordinates = mesh.getX()
158 lgraham 2191
159 caltinay 3288 # gravitational force
160     Y=Vector(0., Function(mesh))
161     Y[1] = -rho*g
162 lgraham 2191
163 caltinay 3288 # element spacing
164     h = Lsup(mesh.getSize())
165 lgraham 2191
166 caltinay 3288 # boundary conditions for slip at base
167     boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
168 lgraham 2191
169 caltinay 3288 # velocity and pressure vectors
170     velocity=Vector(0., Solution(mesh))
171     pressure=Scalar(0., ReducedSolution(mesh))
172 lgraham 2191
173 caltinay 3288 # Stokes Cartesian
174     solution=StokesProblemCartesian(mesh)
175     solution.setTolerance(tolerance)
176 lgraham 2191
177 caltinay 3288 while t <= t_end:
178     print(" ----- Time step = %s -----"%t)
179     print("Time = %s seconds"%time)
180 lgraham 2191
181 caltinay 3288 solution.initialize(fixed_u_mask=boundary_cond, eta=eta, f=Y)
182     velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter, \
183     verbose=verbose)
184 lgraham 2191
185 caltinay 3288 print("Max velocity =", Lsup(velocity), "m/s")
186 lgraham 2191
187 gross 3379 # CFL condition
188 caltinay 3288 dt=0.4*h/(Lsup(velocity))
189     print("dt =", dt)
190    
191     # displace the mesh
192     displacement = velocity * dt
193     coordinates = mesh.getX()
194 jfenwick 3914 newx=interpolate(coordinates + displacement, ContinuousFunction(mesh))
195     mesh.setX(newx)
196 caltinay 3288
197     time += dt
198    
199     vel_mag = length(velocity)
200    
201     #save velocity and pressure output
202     saveVTK("vel.%2.2i.vtu"%t, vel=vel_mag, vec=velocity, pressure=pressure)
203     t = t+1.
204 lgraham 2191 \end{python}
205 lgraham 2193 %
206 caltinay 3288 The results from the simulation can be viewed with \mayavi, by executing the
207     following command:
208 lgraham 2191 %
209 caltinay 3288 \begin{verbatim}
210     mayavi2 -d vel.00.vtu -m SurfaceMap
211     \end{verbatim}
212 lgraham 2191 %
213 caltinay 3288 Colour-coded scalar maps and velocity flow fields can be viewed by selecting
214     them in the menu.
215     The time-steps can be swept through to view a movie of the simulation.
216     \fig{FLUID OUTPUT} shows the simulation output.
217     Velocity vectors and a colour map for pressure are shown.
218     As the time progresses the body of fluid falls under the influence of gravity.
219 lgraham 2193 %
220 gross 2654 \begin{figure}[ht]
221 lgraham 2191 \center
222 caltinay 3288 \subfigure[t=1]{\label{FLOW OUTPUT 01}\includegraphics[height=5cm]{stokes-fluid-t01}}
223     \hspace{1.6cm}
224     \subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokes-fluid-t10}}
225     \hspace{1.6cm}
226     \subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokes-fluid-t20}}\\
227     \subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokes-fluid-t30}}
228     \hspace{1cm}
229     \subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokes-fluid-t40}}
230     \hspace{1cm}
231     \subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokes-fluid-t50}}
232 caltinay 3279 %\includegraphics[scale=0.25]{stokes-fluid-colorbar}
233 caltinay 3288 \caption{Simulation output for Stokes flow. Fluid body starts off as a
234     rectangular shape, then progresses downwards under the influence of gravity.
235     Colour coded distribution represents the scalar values for pressure.
236     Velocity vectors are displayed at each node in the mesh to show the flow field.
237     Computational mesh used was 20$\times$20 elements.}
238     \label{FLUID OUTPUT}
239 lgraham 2191 \end{figure}
240 lgraham 2193 %
241 caltinay 3288 The view used here to track the fluid is the Lagrangian view, since the mesh
242     moves with the fluid. One of the disadvantages of using the Lagrangian view is
243     that the elements in the mesh become severely distorted after a period of time
244     and introduce solver errors. To get around this limitation the Level Set
245     Method can be used, with the Eulerian point of view for a fixed mesh.
246 gross 2371 %The Level Set Method is discussed in Section \ref{LEVELSET CHAP}.
247 caltinay 3288

  ViewVC Help
Powered by ViewVC 1.1.26