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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26