# Diff of /trunk/doc/user/stokesflow.tex

revision 2792 by gross, Tue Sep 8 07:11:12 2009 UTC revision 2793 by gross, Tue Dec 1 06:10:10 2009 UTC
# Line 2  Line 2
2  \section{Stokes Flow}  \section{Stokes Flow}
3  \label{STOKES FLOW CHAP}  \label{STOKES FLOW CHAP}
4
5  In this section we will look at Computational Fluid Dynamics (CFD) to simulate the flow of fluid under the influence of gravity. The StokesProblemCartesian class will be used to calculate the velocity and pressure of the fluid.  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  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:  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
# Line 49  The subscript comma $i$ denotes the deri Line 49  The subscript comma $i$ denotes the deri
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$).  %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}$.  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  The Stokes equations is a saddle point problem, and can be solved using a Uzawa scheme. A class called 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}.  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  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:  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
# Line 59  C = \frac{v \delta t}{h}. Line 59  C = \frac{v \delta t}{h}.
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.  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  The following PYTHON script is the setup for the Stokes flow simulation, and is available in the example directory as 'fluid.py'. It starts off by importing the classes, such as the 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. 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 of the mesh; fluid movement in the x-direction is free, but fixed in the y-direction. An instance of the 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.  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  %  %
72  \begin{python}  \begin{python}
73  from esys.escript import *  from esys.escript import *
# Line 84  verbose=True Line 92  verbose=True
92  H=2.0  H=2.0
93  L=1.0  L=1.0
94  W=1.0  W=1.0
95  mesh = esys.finley.Rectangle(l0=L, l1=H, order=2, n0=20, n1=20)  mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20)
96  coordinates = mesh.getX()  coordinates = mesh.getX()
97
98  #gravitational force  #gravitational force
# Line 95  Y[1]=-rho*g Line 103  Y[1]=-rho*g
103  h=Lsup(mesh.getSize())  h=Lsup(mesh.getSize())
104
105  #boundary conditions for slip at base  #boundary conditions for slip at base
106  boundary_cond=whereZero(coordinates[1])*[0.0,1.0]  boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
107
108  #velocity and pressure vectors  #velocity and pressure vectors
109  velocity=Vector(0.0, ContinuousFunction(mesh))  velocity=Vector(0.0, Solution(mesh))
110  pressure=Scalar(0.0, ContinuousFunction(mesh))  pressure=Scalar(0.0, ReducedSolution(mesh))
111
112  #Stokes Cartesian  #Stokes Cartesian
113  solution=StokesProblemCartesian(mesh)  solution=StokesProblemCartesian(mesh)
# Line 146  Figures \ref{FLUID OUTPUT1} and \ref{FLU Line 154  Figures \ref{FLUID OUTPUT1} and \ref{FLU
154  %  %
155  \begin{figure}[ht]  \begin{figure}[ht]
156  \center  \center
157  \subfigure[t=1]{\label{FLOW OUTPUT 01}\includegraphics[scale=0.25]{figures/stokes-fluid-t01}}  \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.25]{figures/stokes-fluid-t10}}  \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.25]{figures/stokes-fluid-t20}}  \subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[scale=0.2]{figures/stokes-fluid-t20}}
160  \includegraphics[scale=0.25]{figures/stokes-fluid-colorbar}  %\includegraphics[scale=0.25]{figures/stokes-fluid-colorbar}
161  \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.}  \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}  \label{FLUID OUTPUT1}
163  \end{figure}  \end{figure}
164  %  %
165  \begin{figure}[ht]  \begin{figure}[ht]
166  \center  \center
167  \subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[scale=0.25]{figures/stokes-fluid-t30}}  \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.25]{figures/stokes-fluid-t40}}  \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.25]{figures/stokes-fluid-t50}}  \subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[scale=0.2]{figures/stokes-fluid-t50}}
170  \includegraphics[scale=0.25]{figures/stokes-fluid-colorbar}  %\includegraphics[scale=0.25]{figures/stokes-fluid-colorbar}
171  \caption{Simulation output for Stokes flow.}  \caption{Simulation output for Stokes flow.}
172  \label{FLUID OUTPUT2}  \label{FLUID OUTPUT2}
173  \end{figure}  \end{figure}

Legend:
 Removed from v.2792 changed lines Added in v.2793