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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3288 - (show annotations)
Tue Oct 19 01:08:31 2010 UTC (8 years, 11 months ago) by caltinay
File MIME type: application/x-tex
File size: 9184 byte(s)
Corrections to user guide sections 1.5 and 1.6.
Converted backgrounds of the very complex stokes flow eps files by raster
images leaving the vectors untouched so the change is barely visible but
makes loading the PDF page some orders of magnitude faster.

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

  ViewVC Help
Powered by ViewVC 1.1.26