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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3295 - (hide annotations)
Fri Oct 22 01:56:02 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 8958 byte(s)


1 lgraham 2191
2     \section{Stokes Flow}
3     \label{STOKES FLOW CHAP}
4 caltinay 3288 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 lgraham 2191 %
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 caltinay 3288 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 lgraham 2191 %
30     \begin{equation}
31 jfenwick 3295 -(\eta(v_{i,j} + v_{j,i})),_{j} - p,_{i} = f_{i},
32 lgraham 2191 \label{GENERAL NAVIER STOKES COM}
33     \end{equation}
34     %
35     with the incompressibility condition
36     %
37     \begin{equation}
38 jfenwick 3295 -v_{i,i} = 0.
39 lgraham 2191 \label{INCOMPRESSIBILITY COM}
40     \end{equation}
41     %
42 jfenwick 3295 The subscript comma $i$ denotes the derivative of the function with respect to $x_{i}$.
43     %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}:
44 lgraham 2191 %
45     %\begin{equation}
46 jfenwick 3295 %\sigma^{'}_{ij} = 2\eta D^{'}_{ij},
47 lgraham 2191 %\label{STRESS}
48     %\end{equation}
49     %
50 jfenwick 3295 %where the deviatoric stretching $D^{'}_{ij}$ is defined as
51 lgraham 2191 %
52     %\begin{equation}
53 jfenwick 3295 %D^{'}_{ij} = D^{'}_{ij} - \frac{1}{3}D_{kk}\delta_{ij}.
54 lgraham 2191 %\label{DEVIATORIC STRETCHING}
55     %\end{equation}
56     %
57 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$).
58 caltinay 3288 The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in
59 jfenwick 3295 the $x_{3}$ direction and is given as $f=-g\rho\delta_{i3}$.
60 caltinay 3288 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 lgraham 2191 %
68     \begin{equation}
69 caltinay 3291 C = \frac{v \delta t}{h}
70 lgraham 2191 \label{COURANT}
71     \end{equation}
72     %
73 caltinay 3288 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 lgraham 2191
78 caltinay 3288 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 caltinay 3291 on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}.
95 caltinay 3288 The fact that pressure and velocity are represented in different ways is
96     expressed by
97 gross 2793 \begin{python}
98 caltinay 3288 velocity=Vector(0., Solution(mesh))
99     pressure=Scalar(0., ReducedSolution(mesh))
100 gross 2793 \end{python}
101 caltinay 3288 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 jfenwick 3295 face of the domain. At the base fluid movement in the $x_{0}$-direction
105     is free, but fixed in the $x_{1}$-direction, and similarly at the left
106     face fluid movement in the $x_{1}$-direction is free but fixed in
107     the $x_{0}$-direction.
108 caltinay 3288 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 lgraham 2191 %
119     \begin{python}
120 caltinay 3288 from esys.escript import *
121     import esys.finley
122     from esys.escript.linearPDEs import LinearPDE
123     from esys.escript.models import StokesProblemCartesian
124 lgraham 2191
125 caltinay 3288 # physical constants
126     eta=1.
127     rho=100.
128     g=10.
129 lgraham 2191
130 caltinay 3288 # 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 lgraham 2191
138 caltinay 3288 # 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 lgraham 2191
145 caltinay 3288 # gravitational force
146     Y=Vector(0., Function(mesh))
147     Y[1] = -rho*g
148 lgraham 2191
149 caltinay 3288 # element spacing
150     h = Lsup(mesh.getSize())
151 lgraham 2191
152 caltinay 3288 # boundary conditions for slip at base
153     boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
154 lgraham 2191
155 caltinay 3288 # velocity and pressure vectors
156     velocity=Vector(0., Solution(mesh))
157     pressure=Scalar(0., ReducedSolution(mesh))
158 lgraham 2191
159 caltinay 3288 # Stokes Cartesian
160     solution=StokesProblemCartesian(mesh)
161     solution.setTolerance(tolerance)
162 lgraham 2191
163 caltinay 3288 while t <= t_end:
164     print(" ----- Time step = %s -----"%t)
165     print("Time = %s seconds"%time)
166 lgraham 2191
167 caltinay 3288 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 lgraham 2191
171 caltinay 3288 print("Max velocity =", Lsup(velocity), "m/s")
172 lgraham 2191
173 caltinay 3288 # 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 lgraham 2191 \end{python}
190 lgraham 2193 %
191 caltinay 3288 The results from the simulation can be viewed with \mayavi, by executing the
192     following command:
193 lgraham 2191 %
194 caltinay 3288 \begin{verbatim}
195     mayavi2 -d vel.00.vtu -m SurfaceMap
196     \end{verbatim}
197 lgraham 2191 %
198 caltinay 3288 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 lgraham 2193 %
205 gross 2654 \begin{figure}[ht]
206 lgraham 2191 \center
207 caltinay 3288 \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 caltinay 3279 %\includegraphics[scale=0.25]{stokes-fluid-colorbar}
218 caltinay 3288 \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 lgraham 2191 \end{figure}
225 lgraham 2193 %
226 caltinay 3288 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 gross 2371 %The Level Set Method is discussed in Section \ref{LEVELSET CHAP}.
232 caltinay 3288

  ViewVC Help
Powered by ViewVC 1.1.26