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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6928 - (hide annotations)
Mon Dec 9 05:29:21 2019 UTC (16 months, 1 week ago) by acodd
File MIME type: application/x-tex
File size: 9829 byte(s)
more updates to userguide 

1 lgraham 2191
2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 6651 % Copyright (c) 2003-2018 by The University of Queensland
4 jfenwick 3989 % http://www.uq.edu.au
5 caltinay 3325 %
6     % Primary Business: Queensland, Australia
7 jfenwick 6112 % Licensed under the Apache License, version 2.0
8     % http://www.apache.org/licenses/LICENSE-2.0
9 caltinay 3325 %
10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 % Development 2012-2013 by School of Earth Sciences
12     % Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3989 %
14     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 caltinay 3325
16 lgraham 2191 \section{Stokes Flow}
17     \label{STOKES FLOW CHAP}
18 acodd 6928 In this section we look at Computational Fluid Dynamics (CFD) to simulate
19 caltinay 3288 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 lgraham 2191 %
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 caltinay 3288 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 lgraham 2191 %
44     \begin{equation}
45 acodd 6928 -(\eta(v_{i,j} + v_{j,i})),_{j} + p,_{i} = f_{i},
46 lgraham 2191 \label{GENERAL NAVIER STOKES COM}
47     \end{equation}
48     %
49     with the incompressibility condition
50     %
51     \begin{equation}
52 jfenwick 3295 -v_{i,i} = 0.
53 lgraham 2191 \label{INCOMPRESSIBILITY COM}
54     \end{equation}
55     %
56 jfenwick 3295 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 lgraham 2191 %
59     %\begin{equation}
60 jfenwick 3295 %\sigma^{'}_{ij} = 2\eta D^{'}_{ij},
61 lgraham 2191 %\label{STRESS}
62     %\end{equation}
63     %
64 jfenwick 3295 %where the deviatoric stretching $D^{'}_{ij}$ is defined as
65 lgraham 2191 %
66     %\begin{equation}
67 jfenwick 3295 %D^{'}_{ij} = D^{'}_{ij} - \frac{1}{3}D_{kk}\delta_{ij}.
68 lgraham 2191 %\label{DEVIATORIC STRETCHING}
69     %\end{equation}
70     %
71 sshaw 4552 %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 caltinay 3288 The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in
73 jfenwick 3295 the $x_{3}$ direction and is given as $f=-g\rho\delta_{i3}$.
74 caltinay 3288 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 caltinay 5296 In order to keep numerical stability and satisfy the Courant-Friedrichs-Lewy condition (CFL condition)\index{Courant number}\index{CFL condition}, the
79 caltinay 3288 time-step size needs to be kept below a certain value.
80 gross 3379 The Courant number \index{Courant number} is defined as:
81 lgraham 2191 %
82     \begin{equation}
83 caltinay 3291 C = \frac{v \delta t}{h}
84 lgraham 2191 \label{COURANT}
85     \end{equation}
86     %
87 caltinay 3288 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 lgraham 2191
92 caltinay 3288 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 caltinay 5296 solver used is PCG (Preconditioned Conjugate Gradients).
101     The mesh is defined as a rectangle to represent the body of fluid.
102 caltinay 3288 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 caltinay 3291 on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}.
109 caltinay 3288 The fact that pressure and velocity are represented in different ways is
110     expressed by
111 gross 2793 \begin{python}
112 caltinay 3288 velocity=Vector(0., Solution(mesh))
113     pressure=Scalar(0., ReducedSolution(mesh))
114 gross 2793 \end{python}
115 caltinay 3288 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 jfenwick 3295 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 caltinay 3288 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 sshaw 4552 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 caltinay 3288 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 lgraham 2191 %
134     \begin{python}
135 caltinay 3288 from esys.escript import *
136     import esys.finley
137     from esys.escript.linearPDEs import LinearPDE
138     from esys.escript.models import StokesProblemCartesian
139 caltinay 3348 from esys.weipa import saveVTK
140 lgraham 2191
141 caltinay 3288 # physical constants
142     eta=1.
143     rho=100.
144     g=10.
145 lgraham 2191
146 caltinay 3288 # 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 lgraham 2191
154 caltinay 3288 # 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 lgraham 2191
161 caltinay 3288 # gravitational force
162     Y=Vector(0., Function(mesh))
163     Y[1] = -rho*g
164 lgraham 2191
165 caltinay 3288 # element spacing
166     h = Lsup(mesh.getSize())
167 lgraham 2191
168 caltinay 3288 # boundary conditions for slip at base
169 jfenwick 6678 boundary_cond=whereZero(coordinates[1])*[0.0,1.0]
170     +whereZero(coordinates[0])*[1.0,0.0]
171 lgraham 2191
172 caltinay 3288 # velocity and pressure vectors
173     velocity=Vector(0., Solution(mesh))
174     pressure=Scalar(0., ReducedSolution(mesh))
175 lgraham 2191
176 caltinay 3288 # Stokes Cartesian
177     solution=StokesProblemCartesian(mesh)
178     solution.setTolerance(tolerance)
179 lgraham 2191
180 caltinay 3288 while t <= t_end:
181     print(" ----- Time step = %s -----"%t)
182     print("Time = %s seconds"%time)
183 lgraham 2191
184 caltinay 3288 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 lgraham 2191
188 caltinay 3288 print("Max velocity =", Lsup(velocity), "m/s")
189 lgraham 2191
190 gross 3379 # CFL condition
191 caltinay 3288 dt=0.4*h/(Lsup(velocity))
192     print("dt =", dt)
193    
194     # displace the mesh
195     displacement = velocity * dt
196     coordinates = mesh.getX()
197 jfenwick 3914 newx=interpolate(coordinates + displacement, ContinuousFunction(mesh))
198     mesh.setX(newx)
199 caltinay 3288
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 lgraham 2191 \end{python}
208 lgraham 2193 %
209 caltinay 3288 The results from the simulation can be viewed with \mayavi, by executing the
210     following command:
211     \begin{verbatim}
212 caltinay 5296 mayavi2 -d vel.00.vtu -m Surface
213 caltinay 3288 \end{verbatim}
214 lgraham 2191 %
215 caltinay 3288 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 lgraham 2193 %
222 gross 2654 \begin{figure}[ht]
223 lgraham 2191 \center
224 caltinay 3288 \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 caltinay 3279 %\includegraphics[scale=0.25]{stokes-fluid-colorbar}
235 caltinay 3288 \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 lgraham 2191 \end{figure}
242 lgraham 2193 %
243 caltinay 3288 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 gross 2371 %The Level Set Method is discussed in Section \ref{LEVELSET CHAP}.
249 caltinay 3288

  ViewVC Help
Powered by ViewVC 1.1.26