doc/user/stokesflow.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2003-2018 by The University of Queensland
% http://www.uq.edu.au
%
% Primary Business: Queensland, Australia
% Licensed under the Apache License, version 2.0
% http://www.apache.org/licenses/LICENSE-2.0
%
% Development until 2012 by Earth Systems Science Computational Center (ESSCC)
% Development 2012-2013 by School of Earth Sciences
% Development from 2014 by Centre for Geoscience Computing (GeoComp)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Stokes Flow}
\label{STOKES FLOW CHAP}
In this section we 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.
The fluid dynamics is governed by the Stokes equation. In geophysical problems
the velocity of fluids is 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:
%
\begin{equation}
\nabla \cdot (\eta(\nabla \vec{v} + \nabla^{T} \vec{v})) - \nabla p = -f,
\label{GENERAL NAVIER STOKES}
\end{equation}
%
with the incompressibility condition
%
\begin{equation}
\nabla \cdot \vec{v} = 0.
\label{INCOMPRESSIBILITY}
\end{equation}
%
where $p$, $\eta$ and $f$ are the pressure, viscosity and body forces, respectively.
Alternatively, the Stokes equations can be represented in Einstein summation
tensor notation (compact notation):
%
\begin{equation}
-(\eta(v_{i,j} + v_{j,i})),_{j} + p,_{i} = f_{i},
\label{GENERAL NAVIER STOKES COM}
\end{equation}
%
with the incompressibility condition
%
\begin{equation}
-v_{i,i} = 0.
\label{INCOMPRESSIBILITY COM}
\end{equation}
%
The subscript comma $i$ denotes the derivative of the function with respect to $x_{i}$.
%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}:
%
%\begin{equation}
%\sigma^{'}_{ij} = 2\eta D^{'}_{ij},
%\label{STRESS}
%\end{equation}
%
%where the deviatoric stretching $D^{'}_{ij}$ is defined as
%
%\begin{equation}
%D^{'}_{ij} = D^{'}_{ij} - \frac{1}{3}D_{kk}\delta_{ij}.
%\label{DEVIATORIC STRETCHING}
%\end{equation}
%
%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$).
The body force $f$ in \eqn{GENERAL NAVIER STOKES COM} is the gravity acting in
the $x_{3}$ direction and is given as $f=-g\rho\delta_{i3}$.
The Stokes equation 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. A more detailed discussion of the class can be
found in Chapter \ref{MODELS CHAPTER}.
In order to keep numerical stability and satisfy the Courant-Friedrichs-Lewy condition (CFL condition)\index{Courant number}\index{CFL condition}, the
time-step size needs to be kept below a certain value.
The Courant number \index{Courant number} is defined as:
%
\begin{equation}
C = \frac{v \delta t}{h}
\label{COURANT}
\end{equation}
%
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$.

The following \PYTHON script is the setup for the Stokes flow simulation, and
is available in the \ExampleDirectory 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 (Preconditioned Conjugate Gradients).
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 for velocity but the elements are subdivided for the velocity.
This approach is called \textit{macro elements}\index{macro elements} and
needs to be applied to make sure that the discretized problem has a unique
solution, see~\cite{LBB} for details\footnote{Alternatively, one can use
second order elements for the velocity and first order elements for pressure
on the same element. You can set \code{order=2} in \class{esys.finley.Rectangle}.}.
The fact that pressure and velocity are represented in different ways is
expressed by
\begin{python}
  velocity=Vector(0., Solution(mesh))
  pressure=Scalar(0., ReducedSolution(mesh))
\end{python}
The gravitational force is calculated based 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_{0}$-direction
is free, but fixed in the $x_{1}$-direction, and similarly at the left
face fluid movement in the $x_{1}$-direction is free  but fixed in
the $x_{0}$-direction.
An instance of the \class{StokesProblemCartesian} class is defined for the
given computational mesh, and the solver tolerance set.
Inside the \code{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 based on the Courant-Friedrichs-Lewy condition
(CFL condition)\index{Courant number}\index{CFL 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 saved to a file
for visualization.
%
\begin{python}
  from esys.escript import *
  import esys.finley
  from esys.escript.linearPDEs import LinearPDE
  from esys.escript.models import StokesProblemCartesian
  from esys.weipa import saveVTK

  # physical constants
  eta=1.
  rho=100.
  g=10.

  # solver settings
  tolerance=1.0e-4
  max_iter=200
  t_end=50
  t=0.0
  time=0
  verbose=True

  # define mesh
  H=2.
  L=1.
  W=1.
  mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20)
  coordinates = mesh.getX()

  # gravitational force
  Y=Vector(0., Function(mesh))
  Y[1] = -rho*g

  # element spacing
  h = Lsup(mesh.getSize())

  # boundary conditions for slip at base
  boundary_cond=whereZero(coordinates[1])*[0.0,1.0]
      +whereZero(coordinates[0])*[1.0,0.0]

  # velocity and pressure vectors
  velocity=Vector(0., Solution(mesh))
  pressure=Scalar(0., ReducedSolution(mesh))

  # Stokes Cartesian
  solution=StokesProblemCartesian(mesh)
  solution.setTolerance(tolerance)

  while t <= t_end:
    print(" ----- Time step = %s -----"%t)
    print("Time = %s seconds"%time)

    solution.initialize(fixed_u_mask=boundary_cond, eta=eta, f=Y)
    velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter, \
                                     verbose=verbose)

    print("Max velocity =", Lsup(velocity), "m/s")

    # CFL condition
    dt=0.4*h/(Lsup(velocity))
    print("dt =", dt)

    # displace the mesh
    displacement = velocity * dt
    coordinates = mesh.getX()
    newx=interpolate(coordinates + displacement, ContinuousFunction(mesh))
    mesh.setX(newx)

    time += dt

    vel_mag = length(velocity)

    #save velocity and pressure output
    saveVTK("vel.%2.2i.vtu"%t, vel=vel_mag, vec=velocity, pressure=pressure)
    t = t+1.
\end{python}
%
The results from the simulation can be viewed with \mayavi, by executing the
following command:
\begin{verbatim}
mayavi2 -d vel.00.vtu -m Surface
\end{verbatim}
%
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.
\fig{FLUID OUTPUT} 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.
%
\begin{figure}[ht]
\center
\subfigure[t=1]{\label{FLOW OUTPUT 01}\includegraphics[height=5cm]{stokes-fluid-t01}}
\hspace{1.6cm}
\subfigure[t=20]{\label{FLOW OUTPUT 10}\includegraphics[height=5cm]{stokes-fluid-t10}}
\hspace{1.6cm}
\subfigure[t=30]{\label{FLOW OUTPUT 20}\includegraphics[height=5cm]{stokes-fluid-t20}}\\
\subfigure[t=40]{\label{FLOW OUTPUT 30}\includegraphics[height=5cm]{stokes-fluid-t30}}
\hspace{1cm}
\subfigure[t=50]{\label{FLOW OUTPUT 40}\includegraphics[height=5cm]{stokes-fluid-t40}}
\hspace{1cm}
\subfigure[t=60]{\label{FLOW OUTPUT 50}\includegraphics[height=5cm]{stokes-fluid-t50}}
%\includegraphics[scale=0.25]{stokes-fluid-colorbar}
\caption{Simulation output for Stokes flow. Fluid body starts off as a
rectangular shape, then progresses downwards under the influence of gravity.
Colour 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.}
\label{FLUID OUTPUT}
\end{figure}
%
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.
%The Level Set Method is discussed in Section \ref{LEVELSET CHAP}.

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