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

Revision 3325 - (hide annotations)
Fri Oct 29 00:30:51 2010 UTC (11 years, 9 months ago) by caltinay
File MIME type: application/x-tex
File size: 9361 byte(s)
Section 6.1.

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