A numerical approximation for the Navier-Stokes€¦ · testing, automobile design, biodigestor...
Transcript of A numerical approximation for the Navier-Stokes€¦ · testing, automobile design, biodigestor...
A numerical approximation for the Navier-Stokes
equations using the Finite Element Method
Joao Francisco Fragao Rocha Marques
Dissertacao para a obtencao do grau de mestre em
Matematica e Aplicacoes
Orientador: Prof. Doutor Jorge Filipe Duarte Tiago
Juri
Presidente: Prof. Doutora Ana Leonor Mestre Vicente Silvestre
Vogais: Prof. Doutora Marılia da Conceicao Valente de Oliveira Pires
Prof. Doutor Jorge Filipe Duarte Tiago
July 2016
Abstract
In this thesis we describe the numerical approximation of the Navier-Stokes equations, for incom-
pressible Newtonian fluids, using the Finite Element Method. In order to be able to deal with different
typical situations, a Matlab numerical implementation was done and tested for the linear case (Stokes
equations), the stationary convection dominated Navier-Stokes equations and also for the evolutive case.
This included the possibility of dealing with both Dirichlet and Neumann boundary conditions. The
inf-sup stable P2-P1 FEM were used for the spatial discretization. In the case of convection dominated
problems, a stabilization technique was added at the linear iterative step of Newton’s method. The fully
implicit Euler finite difference scheme was used for the time discretization. The starting point for the
Matlab implementation was the Q2-Q1 implementation for the Stokes equations, with Dirichlet boundary
conditions, made available in [9], [8] and [7].
Keywords: Navier-Stokes equations, Numerical approximation, Newton-Raphson Method, Stream-
line Upwind Petrov-Galerkin (SUPG) Method, Finite Element Method (FEM).
i
Resumo
Nesta tese descrevemos abordagens numericas para aproximar as equacoes de Navier-Stokes (de fluidos
incompressıveis e Newtonianos) tendo como base o metodo dos elementos finitos. Para se conseguir lidar
com diferentes situacoes tıpicas de modelacao de escoamentos de fluidos, foram implementados esquemas
numericos em Matlab considerando o caso linear (equacoes de Stokes), as equacoes estacionarias de
Navier-Stokes e finalmente o caso nao estacionario. Estes esquemas numericos incluem a capacidade de
lidar com condicoes de fronteira de Dirichlet e Neumann. Para a discretizacao espacial foram usados
elementos finitos P2-P1 que respeitam a condicao de estabilidade inf-sup discreta. No caso de problemas
onde o termo convectivo e dominante, foi introduzida uma tecnica de estabilizacao no metodo de Newton.
Para a discretizacao temporal foi usado o metodo das diferencas finitas de Euler implıcito. O ponto de
partida para implementar estes algoritmos foi o codigo que usa elementos finitos Q2-Q1 para resolver o
problema de Stokes com condicoes de Dirichlet, disponıvel em [9], [8] e [7].
Palavras Chave: Equacoes de Navier-Stokes, Aproximacao numerica, Metodo de Newton-Raphson,
Metodo Streamline Upwind Petrov-Galerkin (SUPG), Metodo dos Elementos Finitos (FEM).
ii
List of Abbreviations
1. FEM Finite Element Method.
2. SUPG Streamline upwind Petrov-Galerkin.
3. P2-P1 biquadratic, bilinear triangular elements.
4. Q2-Q1 biquadratic, bilinear quadrangular elements.
5. IFISS Incompressible Flow Iterative Solution Software.
iii
Contents
1 Introduction 1
2 Mathematical Models 3
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Non-Stationary Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.3 Stationary Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3 Stokes Equations 6
3.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.2 Weak Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.3 Approximation by the FEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.4 Stability of the Discrete System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.5 P2-P1 Algorithm Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5.1 Domain Discretization and Labeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5.2 Elementwise Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5.3 Global Matrix and Assembling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.5.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.5.5 Added Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.6 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.6.1 Poiseuille Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.6.2 Cavity Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.6.3 Custom Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4 Stationary Navier-Stokes Equations 33
4.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.2 Newton Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3 Numerical Results for Newton Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.3.1 Poiseuille Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.3.2 Cavity Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.4 Newton SUPG Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.5 Numerical Results for SUPG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.5.1 Poiseuille Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.5.2 Cavity Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
iv
5 Non-Stationary Navier-Stokes Equations 52
5.1 Numerical Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.2 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.2.1 Test Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.2.2 Cavity Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.2.3 Womersley Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.4 Obstacle Flow Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6 Conclusion and Future Work 61
7 Appendix 62
7.1 Stokes P2-P1 code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.1.1 P2-P1 Poiseuille Stokes example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.1.2 P2-P1 Cavity Stokes example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.1.3 P2-P1 Custom Stokes example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.1.4 P2-P1 Custom Stokes example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.2 Newton method P2-P1 code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.2.1 P2-P1 Newton Poiseuille Navier-Stokes example . . . . . . . . . . . . . . . . . . . . . 63
7.2.2 P2-P1 Newton Cavity Navier-Stokes example . . . . . . . . . . . . . . . . . . . . . . . 63
7.2.3 P2-P1 Newton Obstacle Navier-Stokes example . . . . . . . . . . . . . . . . . . . . . . 63
7.3 Newton SUPG method P2-P1 code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
7.3.1 P2-P1 Newton SUPG Poiseuille Navier-Stokes example . . . . . . . . . . . . . . . . . 63
7.3.2 P2-P1 Newton SUPG Cavity Navier-Stokes example . . . . . . . . . . . . . . . . . . . 63
7.3.3 P2-P1 Newton SUPG Obstacle Navier-Stokes example . . . . . . . . . . . . . . . . . . 63
7.4 Non-Stationary Newton method P2-P1 code . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
7.4.1 P2-P1 Non-Stationary Navier-Stokes Test Example . . . . . . . . . . . . . . . . . . . . 64
7.4.2 P2-P1 Non-Stationary Navier-Stokes Cavity Example . . . . . . . . . . . . . . . . . . 64
7.4.3 P2-P1 Non-Stationary Navier-Stokes Womersly Example . . . . . . . . . . . . . . . . . 64
7.4.4 P2-P1 Non-Stationary Navier-Stokes Obstacle Example . . . . . . . . . . . . . . . . . 64
v
List of Figures
1 P2-P1 macro-element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2 P2-P1 macro-element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3 Node numbering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4 Element numbering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
5 Local node numbering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6 Element assembling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
7 Domain perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
8 Poiseuille solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
9 Log-Log graph with linear regression for example 1 . . . . . . . . . . . . . . . . . . . . . . . . 28
10 Cavity solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
11 Cavity with streamlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
12 Test solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
13 Stokes Poiseuille problem with obstacle solution . . . . . . . . . . . . . . . . . . . . . . . . . . 32
14 Stokes Poiseuille problem with obstacle ux solution graph . . . . . . . . . . . . . . . . . . . . 32
15 Navier-Stokes solution of Poiseuille flow (Newton) . . . . . . . . . . . . . . . . . . . . . . . . 39
16 Log-Log graph with linear regression for example 1 . . . . . . . . . . . . . . . . . . . . . . . . 40
17 Navier-Stokes cavity solution streamlines (Newton, ν = 1/100 m2/s) . . . . . . . . . . . . . . 41
18 Navier-Stokes cavity solution (Newton) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
19 Navier-Stokes Poiseuille with obstacle solution (Newton) . . . . . . . . . . . . . . . . . . . . . 42
20 Navier-Stokes Poiseuille with obstacle solution streamlines and enlarged ux graph (Newton) . 42
21 Log-Log graph between ||eu1|| and h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
22 Navier-Stokes Poiseuille solution (Newton SUPG) . . . . . . . . . . . . . . . . . . . . . . . . . 48
23 Navier-Stokes cavity solution (Newton SUPG) with ν = 1/1000 m2/s . . . . . . . . . . . . . . 49
24 Navier-Stokes cavity solution streamlines using SUPG for ν = 1/1000 m2/s and ν = 1×10−13
m2/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
25 Navier-Stokes Poiseuille with obstacle solution (Newton SUPG) . . . . . . . . . . . . . . . . . 51
26 Navier-Stokes Poiseuille with obstacle solution streamlines and enlarged ux graph (Newton
SUPG) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
27 Log-Log graph with linear regression for example 1 . . . . . . . . . . . . . . . . . . . . . . . . 56
28 Non-stationary Navier-Stokes cavity simulation (fully implict) . . . . . . . . . . . . . . . . . . 57
29 Non-stationary Navier-Stokes cavity streamline simulation (fully implicit) . . . . . . . . . . . 57
30 Womersley flow solution with ν = 1/10 m2/s . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
31 Different ’water pipe’ slope functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
32 Different time solution plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
33 ||u||2 on different time frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
vi
List of Tables
1 Poiseuille Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2 Test Flow Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3 Navier-Stokes Poiseuille Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4 Cavity Convergence Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5 Navier-Stokes Poiseuille Benchmark (Fully implicit Newton-SUPG) . . . . . . . . . . . . . . . 47
6 Cavity (SUPG fully implicit Newton) Convergence Benchmark . . . . . . . . . . . . . . . . . 49
7 NSNS Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
vii
1 Introduction
The Navier-Stokes equations have been among the most famous mathematical modeling tolls for a few decades
now. Its ability to describe many different types of fluids allows its application to day life type of problems
that span a variety of research and industrial fields such as weather forecast, ocean stream modeling, airplane
testing, automobile design, biodigestor prototyping or, more recently, blood flow simulations. Although in
the past 50 years many advances have been made in terms of the mathematical analysis of the Navier-Stokes
equations in its different versions ( see, for instance, [12] and [32]), for the type of applications mentioned
above, such thing as analytical solutions must obviously be excluded. In this scenario, efficient computational
tools must be used instead, in order to obtain reliable numerical approximations of the real solution. Due to
the nonlinear nature of these equations, typical methods for the numerical approximation of partial differential
equations ([15]), which perform well for linear problems, when kept in their original form, fail to capture the
complexity of the Navier-Stokes equations ([24]) needing to be modified in not so intuitive ways, in order
to be considered stable. This modifications, called stabilization techniques (see [21]), must be chosen and
adjusted according to the specific type of fluid that we are modeling.
At the present stage many commercial and also some free packages can be found to compute numerical
solutions for the Navier-Stokes equations (Open Foam, Comsol Multiphysics, Fenics Project, LifeV, just to
mention a few). However, a robust understanding of the inherent methods, as it is required for research pur-
poses in fluid dynamics, requires a mastering of the implemented algorithm that can hardly be accomplished
by just using such type of ”ready-to-use” packages. In this sense, as part of a major strategy for future
local research activities on this subject, a numerical algorithm to solve the time dependent Navier-Stokes
equations for incompressible Newtonian fluids was implemented. The time discretization was done using an
implicit Euler formula for finite differences, which allows to deal with time changing boundary conditions
including pulsatile flow in a stable way ([24]). The space discretization approach is based on the FEM ([8],
[11] and [21]). To allow the simulation of flow associated to high Reynolds number, stabilization terms were
added to the weak formulation, corresponding to the so called Streamline Upwind Petrov-Galerkin stabiliza-
tion (SUPG). Although different approaches can be used in the spirit of SUPG, here the approach of [11]
was followed. Therefore, the Newton-Raphson (NR) method was used to deal with the nonlinear feature of
the momentum equation, and the SUPG terms were added in each linear iteration of NR method. A sepa-
rated code devoted to the stationary case was also prepared. It was drawn for both high Reynolds number
(convection dominated problem) or low Reynolds number (Stokes flow).
In Chapter 2 the mathematical models are presented: the derivation of the Navier-Stokes equations, its sta-
tionary counterpart and the Stokes linearization. In Chapter 3, the case of the Stokes equations is presented.
This includes the definition of its weak form and the application of the Finite Element Method (FEM) for the
numerical approximation of its solution. The discrete version of the inf-sup stability condition is shown to be
verified for the choice of P2 type elements for the velocity and P1 for the pressure. A numerical implemen-
tation based on an equivalent version for Q2-Q1 elements, given in [8], is described. Written in Matlab, this
1
will be the starting point for sequential generalizations described in the following sections. Some reference
problems ([14]) are used for numerical experiments. In Chapter 4 the application of Newton-Rhapson method
to the stationary case is described. The SUPG is introduced and its stability properties are proved. The
application of FEM is made and numerical results are shown using a generalization of the implementation
described in the previous section.
In Chapter 5 the time discretization is described and a stability estimate is proved. The implementation
is applied to several reference problems including the Womersley pulsatile profile and the flow around an
obstacle inside a pipe.
The thesis ends with some conclusions and additional guiding material for code usage.
2
2 Mathematical Models
2.1 Introduction
First we will give a brief explanation about the origin of the time-dependent Navier-Stokes equations. Next
we will focus on the case of the stationary case of the Navier-Stokes equations (as a particular case of
the non-stationary case). Finally we give some insight about the relation between the non-stationary case
of the Navier-Stokes equations and the Stokes problem. The Navier-Stokes equations are the basis for
computational modeling of the flow of an incompressible Newtonian fluid, such as water. A wide range
of engineering applications is related with formulating and solving such a problem (for example: weather
prediction or airplane designing).
2.2 Non-Stationary Navier-Stokes Equations
The Navier-Stokes equations are derived assuming not only that the velocity and pressure are at least weakly
differentiable, but also that the fluid itself respects the following principles: conservation of mass and con-
servation of momentum. These principles seem fairly acceptable to be enforced in formulating a fluid-like
behavior. Next we will break down each one of these principles and the role each one plays in the final
equation. However we will not be getting in too much detail about the derivation from the principle itself to
the final form as this would be fairly out of context. In essence the mass conservation principle states that
if we consider a generic volume domain (that may change over time) V (t), then we expect that the amount
of mass inside this volume remains the same over time. We would then write it as :
d
dt
∫V (t)
ρdx = 0,∀t ≥ 0
But if we use Reynolds transport theorem, then we are able to state that the previous is equivalent to
∂ρ∂t +∇ · (ρ~u) = 0 and 1
ρDρDt +∇ · (~u) = 0, where D
Dt is the material derivative on a given function. Assuming
the fluid being analyzed is incompressible (DρDt = 0), these two principles entail then, that ∇ · ~u = 0, which
is the second equation of the Navier-Stokes equations. The next principle we rely on is the conservation of
momentum, that is, we assume Newton’s Second law of motion which states that F = m.a. The way we
apply this principle to fluids is very simple. First, we define the forces that we consider to be on the current
fluid volume ball, V (t). They consist of the sum of the internal and external forces applied on each point:
∫V (t)
ρ~g +
∫A(t)
T · ~n (2.2.1)
where g is some external force (such as gravity) being applied to the fluid particles themselves, V (t) is an
3
arbitrary volume ball, A(t) its correspondent surface area, ρ is the fluid density and T is a tensor that is
associated with the internal forces on the fluid itself. The ma part of the Newton law can be translated to
the following expression:
d
dt
∫V (t)
ρ~u (2.2.2)
Hence, if we apply the Reynolds transport theorem once again, we get:
d
dt
∫V (t)
ρ~u =
∫V (t)
D
Dt(ρ~u) + ρ~u∇ · ~u =
∫V (t)
ρD~u
Dt
which by the second law of Newton, must be equal to the expression (2.2.1). Using the latter we can extract
the following formula:
∫V (t)
ρD~u
Dt=
∫V (t)
ρ~g +
∫A(t)
T · n (2.2.3)
By now we have enough information to assemble together and derive the Navier-Stokes equations. Using
what we now from the momentum conservation:
∫V (t)
ρD~u
Dt=
∫V (t)
ρ~g +
∫A(t)
T · n =
∫V (t)
ρ~g +
∫V (t)
∇ · T
We can also extract the corresponding differential form:
ρD~u
Dt= ρ~g +∇ · T (2.2.4)
And if we apply the constitutive law of an incompressible Newton fluid, we get the incompressible Navier-
Stokes equations:
ρ(∂~u
∂t+ (~u · ∇)~u)− µ∆~u+∇p = ρ~g (2.2.5)
∇.u = 0 (2.2.6)
The formulation we have adopted in this thesis is the following (which is obtained merely by dividing equation
(2.2.5) by ρ).
4
∂~u
∂t+ (~u · ∇)~u− ν∆~u+∇p = ~f in Ω× (0,∞) (2.2.7)
∇.u = 0 in Ω× (0,∞) (2.2.8)
Where Ω is the physical bound of our problem, and (0,∞) is the time domain. The interested reader can
find additional information about the Navier-Stokes equations, in sections 7 and 8 of [8] and section 1.3.1
from [17]. Work addressing existence, uniqueness and data dependence of the Navier-Stokes equations can
be seen in [20], [33] and [32].
2.3 Stationary Navier-Stokes Equations
In this subsection we describe the formulation adopted for the stationary case of the Navier-Stokes Equations.
These equations follow simply from assuming that ∂~u∂t = 0 from equation (2.2.7). The following is , therefore,
a possible formulation for the Navier-Stokes equations:
−ν∆~u+ (~u · ∇)~u+∇p = ~f, in Ω (2.3.1)
∇ · ~u = 0, in Ω (2.3.2)
where Ω is the considered domain, ~f is the force applied on the system under Ω, ~u : R2 → R2 for our case
and ν ∈ R+. In addition, both Dirichlet and Neumann boundary conditions may be enforced:
~u = ~g, in ΓΩD(2.3.3)
ν∂~u
∂~n− ~np = ~s, in ΓΩN
(2.3.4)
with ΓΩD∩ ΓΩN
= ∅ (2.3.5)
2.4 Stokes Equations
The generalized form of the Stokes equations can be obtained as a limit case of equation (2.3.1). If we assume
that ν∆~u >> (~u · ∇)~u, then one can approximate (~u · ∇)~u ≈ α~u (see page 435 from [23]). This simplification
origins the generalized form of the Stokes equations:
α~u− ν∆~u+∇p = ~f , in Ω (2.4.1)
∇ · ~u = 0, in Ω (2.4.2)
5
using the same boundary conditions as (2.3.3), (2.3.4) and (2.3.5). In our methods we always consider that
α = 0 and ν = 1 (for simplification reasons). Our formulation for the Stokes problem is therefore:
−∆~u+∇p = ~f (2.4.3)
∇ · ~u = 0 (2.4.4)
where ~u : Ω→ R2 is the fluid velocity, p is the pressure and f is an optional external force to be applied on
the fluid particles. Using as boundary conditions:
∂~u
∂n− n~p = ~s, in ΓΩN
(2.4.5)
~u = ~g, in ΓΩD(2.4.6)
where ΓΩNrefers to the domain (Ω) part that is associated with Neumann boundary conditions and ΓΩD
is
associated with Dirichlet boundary conditions. In addition we enforce that ΓΩN∩ΓΩD
= ∅ and ΓΩN∪ΓΩD
=
Ω. Since we are assuming that ν∆~u >> (~u · ∇)~u, this model is usually good when we are modeling an
incompressible, very viscous flow.
3 Stokes Equations
3.1 Description
In this section we solve numerically the problem of Stokes, presented in sub-section (2.4). First we walk
through the weak formulation for formulas (2.4.3) and (2.4.4). Next we solve the said weak formulation by
making use of FEM (using P2-P1 elements). A stability result is given for our problem afterwards, and finally
a few examples of solutions produced by our algorithms are presented and analyzed.
3.2 Weak Formulation
One way of solving the Stokes PDE, is to discretize it by using the FEM. In this section we derive the weak
form. But before that, we first need to define a few functional spaces where we will be working on. We start
with a very important, Sobolev Space, that we will be using from now on (with the norm defined as usual):
H1(Ω) = ~u : Ω→ R|~u, ∂~u∂x,∂~u
∂y∈ L2(Ω)
Next we define the space where the solution will lie in:
6
H1E(Ω) = ~u ∈ [H1(Ω)]2 |~u = ~g,in ΓΩD
With L2(Ω) being the space with the functions regular enough so that they verify∫
Ω|v|2 <∞. Additionally,
we also need to make use of test functions, which are an essential ingredient of the weak formulation. These
functions live in a different space:
H1E0
(Ω) = ~u ∈ [H1(Ω)]2 |~u = 0, in ΓΩD
As far as the pressure is concerned, we consider both the solution p and the test functions q to be in L2(Ω).
Now we have enough ingredients to proceed with the weak formulation process. As usual, we multiply the
first equation of system (2.4.3), both sides by an arbitrary test function ~v ∈ H1E0
and integrate over Ω:
−∫
Ω
∆~u · ~v +
∫Ω
∇p · ~v =
∫Ω
f · ~v (3.2.1)
However, by appling Green’s formula, we know that:
∫Ω
∆~u · ~v +
∫Ω
∇~u : ∇~v =
∫ΓΩ
(~n · ∇~u) · ~v (3.2.2)
Using the later, along with the fact that ~v must be zero on the Dirichlet boundary, we can get another way
of writing∫
Ω∆~u · ~v:
∫Ω
∆~u · ~v = −∫
Ω
∇~u : ∇~v +
∫ΓΩD
(~n · ∇~u) · ~v +
∫ΓΩN
(~n · ∇~u) · ~v = −∫
Ω
∇~u : ∇~v +
∫ΓΩN
(~n · ∇~u) · ~v (3.2.3)
Replacing this on equation (3.2.1), we get:
∫Ω
∇~u : ∇~v −∫
ΓΩN
(~n · ∇~u) · ~v +
∫Ω
∇p · ~v =
∫Ω
~f · ~v (3.2.4)
We can simplify the equation further, by applying Green’s Formula again on the term∫
Ω∇p · ~v:
7
∫Ω
∇p · ~v +
∫Ω
p · (∇ · ~v) =
∫ΓΩN
(p · ~n) · ~v (3.2.5)
Combining this equation with (2.2.4) (and taking into account the expression for the Neumann boundary
conditions):
∫Ω
∇~u : ∇~v −∫
Ω
p · (∇ · ~v) =
∫Ω
~f · ~v +
∫ΓΩN
(−p · ~n+ (~n · ∇~u)) · ~v =
∫Ω
~f · ~v +
∫ΓΩN
(−p · ~n+ (∂~u
∂~n) · ~v
(3.2.6)
=
∫Ω
~f · ~v +
∫ΓΩN
~s · ~v
(3.2.7)
To get the weak formulation of the second equation, one simply multiplies ∇ · ~u = 0 by a pressure test
function q ∈ L2(Ω) and integrate over Ω. Thus, getting: q · ∇ · u = 0. Finally, linking all the pieces together,
we obtain the weak formulation of the Stokes problem. Find ~u ∈ H1E(Ω), p ∈ L2(Ω) such that:
∫Ω
∇~u : ∇~v −∫
Ω
p · (∇ · ~v) =
∫Ω
~f · ~v +
∫ΓΩN
~s · ~v (3.2.8)∫Ω
q · ∇ · ~u = 0 (3.2.9)
For all ~v ∈ H1E0
(Ω) , q ∈ L2(Ω). We do not address here the existence and unicity of the Stokes flow solution,
such results can be found, for instance, in [15] and [12].
3.3 Approximation by the FEM
To solve system (3.2.8) and (3.2.9) numerically, we need to discretize it. This is first accomplished by the
discretization of our domain Ω in several P2 triangular elements (our mesh) and then by considering a finite
dimensional subspace Xh of the original function space H1E . Generally one considers that such space, Xh,
is generated by a base φi0≤i≤n for some n (so that our basis functions work as an interpolation over the
said mesh). Generally one denotes the set of all elements in a mesh by τ . To give information about how
fine is our grid we use a subscript h on τh. This h means that the biggest diameter among the triangles of
our mesh has value h. Similarly one also considers a finite-basis subspace Mh for the pressure whose basis
is ψi0≤i≤m (these basis functions only interpolate on half of the nodes used for velocity interpolation and
are linear). We also define here Xh0 to be the set of functions in Xh that have zero value on the Dirichlet
boundary. The choices on the number of functions that form either the basis of Xh or Mh is directly related
8
to the discretization of the domain Ω. Naturally, these choices will have direct impact on the results of the
computed solution. We will now proceed to formulate the discrete version of the weak formulation of the
Stokes equations that were obtained in the previous subsection. Find ~uh ∈ Xh, ph ∈Mh such that:
∫Ω
∇~uh : ∇~vh −∫
Ω
ph · (∇ · ~vh) =
∫Ω
~f · ~vh +
∫ΓΩN
~s · ~vh (3.3.1)∫Ω
qh · ∇ · ~uh = 0 (3.3.2)
for all ~vh ∈ Xh0 , qh ∈Mh. Now the key observation is that the previous holds if the following also holds:
∫Ω
∇~uh : ∇~φi −∫
Ω
ph · (∇ · ~φi) =
∫Ω
~f · ~φi +
∫ΓΩN
~s · ~φi , i = 1, ..., n (3.3.3)∫Ω
qh · ∇ · ~uh = 0 (3.3.4)
where φi are basis functions of the space Xh. Now in order to have a full linear model that we can later
solve, we need first to consider the expansion of the aproximated solution ~uh under the space Xh:
~uh =
nin∑j=1
uj~φj +
nin+nd∑j=nin+1
uj~φj (3.3.5)
where nin is the number of interior nodes on the domain (which are associated with the velocity) and nd is
the number of Dirichlet boundary nodes, which we assume to be known.By replacing the expansion of uh in
(3.3.3) we get the following system:
nin∑j=1
uj
∫Ω
∇~φj : ∇~φi −∫
Ω
ph.(∇ · ~φi) =
∫Ω
~f · ~φi +
∫ΓΩN
~s · ~φi −nin+nd∑j=nin+1
uj
∫Ω
∇~φj : ∇~φi (3.3.6)
Due to the discretized pressure ph and the discrete pressure test function qh, we don’t have a yet a full linear
system. In order to achieve that we essentially use the same idea that we used with the velocity. First we
consider a expansion of the pressure on the function space Mh:
ph =
m∑k=1
ψk.pk (3.3.7)
And we assume that if (3.3.2) holds for every qh then it should also hold for every ψk. We have then the
almost complete linear system:
9
nin∑j=1
uj
∫Ω
∇~φj : ∇~φi −m∑k=1
pk
∫Ω
ψk.(∇ · ~φi) =
∫Ω
~f · ~φi +
∫ΓΩN
~s · ~φi −nin+nd∑j=nin+1
uj
∫Ω
∇~φj : ∇~φi (3.3.8)∫Ω
ψk · ∇ · ~uh = 0 (3.3.9)
We still need to take one additional step. Note that in equation (3.3.9) we have all ~u nodes on the left hand
side of the equation, which does not make sense since we are only solving for nodes whose value is unknown
(we already know the Dirichlet boundary nodes with respect to ~u). So we do:
−nin∑j=1
uj
∫Ω
ψk · ∇ · φj =
nin+nd∑j=nin+1
uj
∫Ω
ψk · ∇ · φj (3.3.10)
which can be re-formulated to the following:
A BT
B O
.
up
=
fg
(3.3.11)
where:
Aij =
∫Ω
∇~φi : ∇~φj (3.3.12)
Bkj = −∫
Ω
ψk.(∇ · ~φi) (3.3.13)
fi =
∫Ω
~f · ~φi +
∫ΓΩN
~s · ~φi −nin+nd∑j=nin+1
uj
∫Ω
∇~φj : ∇~φi (3.3.14)
gk =
nin+nd∑j=nin+1
uj
∫Ω
ψk · ∇ · φj (3.3.15)
We have, in the meanwhile, an additional concern that is related to the velocity basis functions which is the
fact that we have a bi-dimensional space as a co-domain for these functions. This means that in some way,
we will have to encode both the x-axis and y-axis velocity information on the previous system. We do this
in the following way:
10
A O BTx
O A BTy
Bx By O
.
ux
uy
p
=
fx
fy
g
(3.3.16)
(3.3.17)
where :
Aij =
∫Ω
∇φi · ∇φj (3.3.18)
Bx,kj = −∫
Ω
ψk ·∂φi∂x
(3.3.19)
By,kj = −∫
Ω
ψk ·∂φi∂y
(3.3.20)
One should note that φi are no longer vector basis functions. Instead, we have twice the number of basis
functions but from now on, scalar functions.
3.4 Stability of the Discrete System
In this section we provide a proof that the P2-P1 finite element method approximation for the Stokes problem
satisfies a certain discrete inf-sup condition. Before introducing the result we define the following notion of
equivalence between macro-elements.
Definition 3.4.1 A macro-element M is said to be equivalent to a reference macro-element M if there is a
mapping FM : M →M that satisfies the following conditions:
1. FM is continuous and bijective.
2. FM (M) = M .
3. If M =⋃mj=1 Kj, where Kjj = 1, 2, ...,m, are the triangles or quadrilaterals in M , then Kj = FM (Kj),
j = 1, 2, ...,m, are the triangles or quadrilaterals in M .
4. FM|Kj= FKj
F−1
Kj, j = 1, 2, ...m, where FKj
and FKjare the affine or bilinear mappings from the
reference triangle (with vertexes (0,0), (0,1) and (1,0)) or unit square onto Kj and Kj, respectively.
The following is the discrete inf-sup condition that we will be proving in our further on (for additional
theorems regarding discrete inf-sup and their consequences, see [2], in addition, a proof that is not based in
macro-elements is present in [8]).
11
Theorem 3.4.2 Let εM be a class of equivalent triangular macro-elements. Then there is a positive constant
βM = β(M) such that the condition:
infph∈Mh
0 ,||ph||0,M=1sup
vh∈Xh0 ,|vh|1,M=1
(∇ · vh, ph)M ≥ βM > 0. (3.4.1)
holds for every M ∈ εM , in the case of the discrete spaces Xh0 , Mh
0 used in the P2-P1 FEM formulation.
Proof: First consider the following space, NM :
NM = qh ∈Mh|∫
Ω
qh∇ · ~vh = 0, ∀~vh ∈ Xh0 (3.4.2)
Our aim is to prove that this space is one-dimensional (this will later come in handy to prove our result by
applying another theorem ). NM space is just isomorphic to the null space of BT , this observation can be
made by considering the discrete version of NM :
NhM = q ∈ Rm| < BT q, v >= 0, ∀v ∈ Rnu (3.4.3)
This is what we will now prove, that is, we will prove that the null space of BT in our case will have dimension
1. The proof will follow by induction: we will first be considering a trivial domain with three elements and
then we will prove the induction step (assuming a valid domain, we will have a valid domain as well if we
add another element). First lets then assume that M is given by the collection of three triangles as in the
figure 1
Figure 1: P2-P1 macro-element
Note that we can expand the pressure by naturally doing p(x, y) = p1l1(x, y)+p2l2(x, y)+p3l3(x, y)+p4l5(x, y),
where each li is the bilinear function associated with the pressure on the node i. Another aspect to be taken
into account is the configuration of the space XM . In our case XM =< q6ey, q6ex, q7ex, q7ey >, where q6
12
and q7 are the bi-quadratic interpolation functions associated with nodes 6 and 7, respectively. If we expand
(∇ · v, p)M :
(∇ · v, p)M = −(v,∇p)M = −3∑j=1
(v,∇p)Tj(3.4.4)
which is obtained from integration by parts, furthermore:
(v,∇p)Tj=
∫Tj
v · ∇p (3.4.5)
However,∫Tjv · ∇p does not pose a real challenge to calculate because, like we did before, we can change the
integration domain to a reference element Tj and since the integrand function will not have a degree higher
than 2 we can integrate by Gauss quadrature. We will omit some of the integration details for the sake of size.
First lets calculate∫T1vj · ∇p (we must take into account that we defined v1 = q6ex, v2 = q6ey, v3 = q7ex
and v4 = q7ey):
∫T1
v1 · ∇p =1
6((y2 − y3)p1 + (y3 − y1)p2 + (y1 − y2)p3) (3.4.6)∫
T1
v2 · ∇p = −1
6((x2 − x3)p1 + (x3 − x1)p2 + (x1 − x2)p3) (3.4.7)
v3 and v4 are omitted because they have zero value on T1. Next we show similar results for T2:
∫T2
v1 · ∇p =1
6((y3 − y5)p1 + (y5 − y1)p3 + (y1 − y3)p5) (3.4.8)∫
T2
v2 · ∇p = −1
6((x3 − x5)p1 + (x5 − x1)p3 + (x1 − x3)p5) (3.4.9)∫
T2
v3 · ∇p =1
6((y3 − y5)p1 + (y5 − y1)p3 + (y1 − y3)p5) (3.4.10)∫
T2
v4 · ∇p = −1
6((x3 − x5)p1 + (x5 − x1)p3 + (x1 − x3)p5) (3.4.11)
Similarly, we present the equivalent results for the T3 triangle:
∫T3
v3 · ∇p =1
6((y4 − y5)p3 + (y5 − y3)p4 + (y3 − y4)p5) (3.4.12)∫
T3
v4 · ∇p = −1
6((x4 − x5)p3 + (x5 − x3)p4 + (x3 − x4)p5) (3.4.13)
13
v1 and v2 are avoided because they evaluate to zero on T3. The final matrix BT in our case looks like:
(y2 − y5) (−y1 + y3) (−y2 + y5) 0 (y1 − y3)
(x2− x5) (−x1 + x3) (−x2 + x5) 0 (x1− x3)
(y3− y5) 0 (−y1 + y4) (−y3 + y5) (y1− y4)
(x3− x5) 0 −x1 + x4 −x3 + x5 x1− x4
(3.4.14)
and if we compute the null space of this matrix using NullSpace[] from Mathematica we get the return
value of 1, 1, 1, 1, 1. Which not only means that the kernel dimension of BT is one but that ph will
only vary by a constant. Now we need to proceed to the induction step, that is, we assume that we have a
stable macro-element M such that dim(ker(BT )) = 1 and we add a new element (the new macro-element
configuration will be called M). We must consider two cases when adding a new element to the existing macro-
element: (i) the new element shares two edges with M or (ii) the new element only shares one edge with
M . Lets denote by BT the matrix associated with NhM and by BT the matrix associated with Nh
M. Naturally,
because the rank must be allways the number of columns minus the dimension of the kernel, we know that
rank(BT ) = nMh− 1, where nMh
is the dimension of the discrete pressure space. The row dimension of BT
is naturally m = 2 × (number of interior edges in M) ≥ nMh− 1. Now lets analyse BT : the row dimension
of this matrix is m ≥ nMh+ 1 since we are assuming that we have gained two new interior edges. Notice
that ∀vi ∈ XM we have∫Mvi · ∇p = 0 for any p a constant function, hence the kernel dimension must
be at least 1, therefore: rank(BT ) ≤ nMh− 1. Since rank(BT ) = nMh
− 1 that means we have nMh− 1
linearly independent rows on BT , but since these rows also exist on BT then BT has at least nMh−1 linearly
independent rows, hence rank(BT ) = nMh− 1 which means that the dimension of Nh
Mis 1. The part (ii) of
this proof is when we consider that only one additional interior edge was gained by adding the new element.
In this new case, however, we have only gained one additional interior node, but we have gained two new
exterior edges. This means that we have an additional dimension on the pressure space. Again, as we have
that ∀vi ∈ XM
∫Mvi ·∇p = 0, for an p a constant function, the kernel dimension of BT must be at least one,
therefore, rank(B) ≤ nMh. Furthermore, since for v ∈ XM , v|M/M = 0, we can tell how BT looks like:
BT =
BT 0
0 b
(3.4.15)
where the number of rows of the matrix B is two, that correspond to the basis functions v1, v2, associated
with the new shared triangle edge vertex. Figure 2 illustrates the case we are working on.
14
Figure 2: P2-P1 macro-element
Some additional calculations show that b = [−y2 + y3,−x2 + x3]. As R2 6= R3, the number of independent
rows in BT must be greater than the independent number of rows in BT . Hence rank(BT ) = nMhwhich
means that the dimension of NhM
is 1. Now we consider the following result, present in [31].
Theorem 3.4.3 Let εM be a class of equivalent macro-elements. Suppose that NM is one-dimensional,
consisting of functions that are constant on M. Then there is a positive constant βM = β(M) such that the
condition:
supvh∈Xh0 ,vh 6=0
(∇ · vh, ph)M|vh|1,M
≥ βM ||ph||0,M , ∀ph ∈Mh0 (3.4.16)
holds for every M ∈ εM .
For the proof, see page 12 of [31]. As we have seen before we are under the conditions of the previous theorem,
which means that formula (3.4.16) entails, and this formula is equivalent to (3.4.1) (see the demonstration
of the previous theorem on page 12 of [31]).
3.5 P2-P1 Algorithm Description
In this section we will explain on how the actual finite element method using P2-P1 elements is implemented
according to our code. As we walk through the implementation, we will also present the code in which such
implementation details are present. From now on, we are assuming that the approximation spaces explained
in the previous section, namely Xh and Mh, are piecewise quadratic and linear, respectively.
3.5.1 Domain Discretization and Labeling
The first part will address the discretization of the domain Ω and the different aspects of the discretization
labeling convention. Just to be clear, after what has been discussed in the previous section, we are interested
in calculating Ai,j , and Bi,j . Thats the whole point of the discretization process, so that we have a very
good estimate for these integrals, and (using Gauss quadrature) expectedly, achieve great results as far as
15
the solution accuracy is concerned. But first lets take a look at a few basic things, such as how the velocity
nodes, pressure nodes and elements are labeled. The labeling mechanism follows similar standards, which is
basically left-right and bottom-up (except for local node labeling). So, if we wanted to label each velocity
node with a number, it would be something like in figure 1.
Similar labeling convection is used for the elements in the domain, and the pressure nodes. Naturally, since
we are assuming P1 for the pressure nodes, we only consider 3 nodes per element as opposed to the 6 velocity
nodes per element. Take a look at figures 1.b and 2, for an example on the labeling conventions used on the
pressure and elements, respectively.
(a) Velocity node numbering (b) Pressure node numbering
Figure 3: Node numbering
Figure 4: Element numbering
The local node numbering convention (inside each element) follows a different numbering convention. Figure
5 illustrates this convention.
16
Figure 5: Local node numbering
where nodes 1,2 and 3 are for pressure and 1, 2, 3, 4, 5 and 6 are used for velocity. In order to access any
information on the nodes or elements, usually one uses an array (or matrix) and chooses as an index the
node or element label. Take the nodes coordinates as an example: a matrix xyi,j (this matrix is actually
constructed inside the code) is defined and one would choose as i the node label as the previous process
showed, and j as 1 or 2 to get the i-th node x or y coordinate respectively. In addition to node and element
labeling, there is also a boundary labeling on the domain Ω. Ω is defined as the [−1, 1] × [−1, 1] square.
Boundary labeling is easy and very straight forward:
Boundary number one: [−1, 1]× −1
Boundary number two: 1 × [−1, 1]
Boundary number three: [−1, 1]× 1
Boundary number four: −1 × [−1, 1]
This is the general logic regarding accessing and setting attributes of the nodes and/or elements.
3.5.2 Elementwise Integration
The calculation of Ai,j poses a small challenge -∫
Ω∇φi · ∇φj is hard to deal with on a first sight because we
need a both efficient and accurate way to calculate what lies inside the integral in an automated way. The
solution is essentially breaking down the calculation of the integral by the several different elements and then
assembling back everything together on a more general matrix. However, we still need to be able to calculate
the quantity∫4∇φj · ∇φi, and this is done by changing the coordinates to a reference element, so that the
calculation is almost independent of the shape of the current element. In order to do that we need several
tools, in particular we need a local labeling mechanism on each element of the domain. Naturally we also
need to have a correspondence from a local node, to its general correspondent node, on the domain scale. In
our code, a P2 macro-element mapping matrix is created for this purpose. This matrix is often denoted by
mv and its size is ne× 6 (where ne is the number of elements in the entire domain). The entry mvi,j gives
the global index of the local node j on the element i. The same needs to be done for the P1 nodes, therefore
a matrix often denoted by mp serves this purpose. Notice that this matrix only has 3 rows instead of six
17
because each element will only contain 3 nodes. Certainly, we need a data structure to store the quantities∫4∇φi · ∇φj for each element and each pair of local nodes i and j (so that later we are able to assemble
everything). Our code uses a tensor to store that information that is named ae and whose dimensions are
ne× 6× 6. Similar variables are created for the matrix B (with the dimensions set accordingly). The next
step is to change the calculation of∫4k∇φi ·∇φj to
∫4∗∇φi ·∇φj , where 4∗ is the reference element and 4k
is an arbitrary domain element. The reference element chosen is the right triangle with both sides in [0, 1].
So a simple change of coordinates is what it takes to do this, but first we need to define a linear mapping
from the 4∗ to 4k. For triangles, the following interpolation polynomial is adopted:
let s, t ∈ 4∗, and x, y ∈ 4k (3.5.1)
x(s, t) = x1X1(s, t) + x2X2(s, t) + x3X3(s, t) (3.5.2)
y(s, t) = y1X1(s, t) + y2X2(s, t) + y3X3(s, t) (3.5.3)
(3.5.4)
where
X1(s, t) = 1− s− t (3.5.5)
X2(s, t) = s (3.5.6)
X3(s, t) = t (3.5.7)
where xi and yi follow the aforementioned labeling convection for P1 nodes on the local element k.This
defines a linear mapping between our reference element and the current element k, lets call this mapping
Tk(s, t) = (x(s, t), y(s, t)). Now we can just change the integration to the reference element by doing:
∫4k
∇φi(x, y) · ∇φj(x, y)dV =
∫4∗∇φi(s, t) · ∇φj(s, t)|JTk
|dV (3.5.8)
where JTkis the Jacobian matrix of the coordinate change mentioned above that will be calculated next.
Since:
Tk(s, t) = (
3∑i=1
xiXi(s, t),
3∑i=1
xiYi(s, t)) (3.5.9)
we can immediately calculate JTk:
18
JTk=
∑3i=1 xi
∂Xi
∂s
∑3i=1 yi
∂Xi
∂s∑3i=1 xi
∂Xi
∂t
∑3i=1 yi
∂Xi
∂t
(3.5.10)
thus making the calculation of |JTk| trivial. However we still need to calculate dφ
dx ,dφdy . There is one trick we
can use to ease our calculations a little bit. First we know that:
dφdx
dφdy
=
∂s∂x
∂t∂x
∂s∂y
∂t∂y
.
∂φ∂s
∂φ∂t
(3.5.11)
secondly we also know that:
∂φ∂s
∂φ∂t
=
∂x∂s
∂y∂s
∂x∂t
∂y∂t
.
dφdx
dφdy
(3.5.12)
where ∂x∂s , ∂x
∂t , ∂y∂s and ∂y
∂t are values that we can specifically calculate. The values ∂s∂x ,
∂t∂x ,
∂s∂y and ∂t
∂y are
unknown but we can know these (and consequentially dφdx and dφ
dy ) simply by inverting the matrix in (3.5.12).
After calculating the inverse in (3.5.12), we find out that values of the entries in (3.5.11) are:
∂s
∂x=
1
|JTk|∂y
∂t(3.5.13)
∂t
∂x= − 1
|JTk|∂y
∂s(3.5.14)
∂s
∂y= − 1
|JTk|∂x
∂t(3.5.15)
∂t
∂y=
1
|JTk|∂x
∂s(3.5.16)
Therefore, now we have enough ingredients to proceed with the local integrations we had discussed previously.
19
For the case of aek,i,j , we have:
aek,i,j =
∫4k
dφidx· dφjdx
+dφidy· dφjdy
=
∫4∗
dφidx· dφjdx
+dφidy· dφjdy|JTk| =
(3.5.17)∫4∗
((∂s
∂x· ∂φi∂s
+∂t
∂x· ∂φi∂t
) · ( ∂s∂x· ∂φj∂s
+∂t
∂x· ∂φj∂t
) + (∂s
∂y· ∂φi∂s
+∂t
∂y· ∂φi∂t
) · (∂s∂y· ∂φj∂s
+∂t
∂y· ∂φj∂t
))|JTk| =
(3.5.18)∫4∗
(1
JTk
∂y
∂t· ∂φi∂s− 1
JTk
∂y
∂s· ∂φi∂t
) · ( 1
JTk
∂y
∂t· ∂φj∂s− 1
JTk
∂y
∂s· ∂φj∂t
)|JTk|+
(3.5.19)∫4∗
(− 1
JTk
∂x
∂t· ∂φi∂s
+1
JTk
∂x
∂s· ∂φi∂t
) · (− 1
JTk
∂x
∂t· ∂φj∂s
+1
JTk
∂x
∂s· ∂φj∂t
)|JTk| =
(3.5.20)∫4∗
(∂y
∂t· ∂φi∂s− ∂y
∂s· ∂φi∂t
) · (∂y∂t· ∂φj∂s− ∂y
∂s· ∂φj∂t
) +
∫4∗
(−∂x∂t· ∂φi∂s
+∂x
∂s· ∂φi∂t
) · (−∂x∂t· ∂φj∂s
+∂x
∂s· ∂φj∂t
)1
|JTk|
(3.5.21)
The last step in calculating the integral in (3.5.21) is to approximate the integration via a Gauss quadrature.
The Gauss quadrature method is element shape dependent, and it basically means that we approximate the
calculation of a single integral of a given f function, by making a weighted sum of the values f has on fixed
points on the element. We could apply this principle to a generic integral over the reference triangle element
4∗:
∫4∗
f(s, t)dsdt =
k max∑k=1
wk.f(sk, tk) (3.5.22)
In our case, quantities wk, sk and tk are calculated from [25]. The final calculation for expression (3.5.21),
using the Gauss quadrature is:
k max∑k=1
wk(∂y
∂t· ∂φi∂s− ∂y
∂s· ∂φi∂t
) · (∂y∂t· ∂φj∂s− ∂y
∂s· ∂φj∂t
)|(sk,tk)1
|JTk|+ (3.5.23)
k max∑k=1
(−∂x∂t· ∂φi∂s
+∂x
∂s· ∂φi∂t
) · (−∂x∂t· ∂φj∂s
+∂x
∂s· ∂φj∂t
)|(sk,tk)1
|JTk|
(3.5.24)
And this is the method used by IFISS software to calculate the entries in aek,i,j . For the sake of completeness,
we will also leave a brief explanation of the formula used to calculate the entries of bxek,i,j , which is the ae
analogue for the matrix Bx,kj presented in (3.3.19). By,kj follows the exact same logic and therefore it will
20
not be covered. Using what we already know, we see that:
bxek,i,j =
∫4k
−ψidφjdx
dxdy =
∫4∗−ψi
dφjdx|JTk|dsdt =
∫4∗−ψi(
∂s
∂x· ∂φj∂s
+∂t
∂x· ∂φj∂t
)|JTk|dsdt =
(3.5.25)∫4∗−ψi(
1
JTk
∂y
∂t· ∂φj∂s
+− 1
JTk
∂y
∂s· ∂φj∂t
)|JTk|dsdt =
∫4∗−ψi(
∂y
∂t· ∂φj∂s
+−∂y∂s· ∂φj∂t
)dsdt
(3.5.26)
And by applying the Gauss quadrature method once more we obtain:
bxek,i,j = −k max∑k=1
wkψi(∂y
∂t· ∂φj∂s
+−∂y∂s· ∂φj∂t
)|(sk,tk) (3.5.27)
We end this subsection hoping to have given a clear insight about the implementation of the algorithm used
by IFISS to calculate the elementwise integration matrices.
3.5.3 Global Matrix and Assembling
Once we have managed to build the local (elementwise) integration matrices, we need to assemble them back
together in once piece to form what is known as the stiffness matrix. The process of assembling can be
well summed up by the following pictures. Suppose you have 4 distinct local matrices (integrating the same
function) as it is shown in figure 6.
21
(a) Element 4
(b) Element 3
(c) Element 1(d) Element 2
Figure 6: Element assembling
Then we want to assemble them in a way such that the original stiffness matrix is recovered (see Figure 7).
Figure 7: Domain perspective
After the tables introduced in section 3.5.2 this is actually a very easy task. All we need to do is take
advantage of the translation table already introduced, mv, and build a new stiffness matrix Anv×nv (where
nv is a variable that denotes the number of P2 nodes on a domain perspective). The way IFISS software
works, is by first allocating such matrix and then adding integration information from the various local
22
matrices. One would express it as (in the case of matrix A):
Amvk,i,mvk,j= Amvk,i,mvk,j
+ aek,i,j , ∀ 0 ≤ i, j ≤ 6, k ∈ τh (3.5.28)
The same principle can be applied to Bx and By. However the dimensions of these matrices are different,
namely they are both np×nv (where np is the number of P1 nodes on the domain level), and consequently,
we need to apply different node label translation on the rows. To be more precise the implementation is done
using the following formula:
lBx,mpk,i,mvk,j= Bx,mpk,i,mvk,j
+ bxek,i,j (3.5.29)
By,mpk,i,mvk,j= By,mpk,i,mvk,j
+ byek,i,j (3.5.30)
∀ 0 ≤ j ≤ 6, 0 ≤ i ≤ 3, k ∈ τh (3.5.31)
This concludes the stiffness matrix assembling process. To recover the left hand side matrix of the system in
(3.3.16), standard built-in Matlab concatenation functions are used.
3.5.4 Boundary Conditions
So far we have computed with success the left hand side of (3.3.16). But what about the right hand side?
And what changes need to be done to the left hand side in order to do that as well? First we will address
the Dirichlet boundary conditions relative to the right hand side of (3.3.14) and (3.3.15). These quantities
can be calculated by multiplying a vector that is the result of applying ~g to all points that belong to ΓΩD
and multiply it by A and B and thus subtracting these values to the right hand side of (3.3.11) (so that we
recover∑nin+nd
j=nin+1 uj∫
Ω∇~φj : ∇~φi and
∑nin+nd
j=nin+1 uj∫
Ω~ψk : ∇ · ~φj) . However we need to enforce that we
already know that
ui = gi in ΓΩD(3.5.32)
The way we do this is simply by zeroing out the columns of the A and B matrices corresponding to the nodes
i ∈ ΓΩDso that we can be sure that no nodes are both in the Dirichlet boundary and unknowns at the same
time. In order to keep the same number of unknowns and equations, we also substitute the row associated
with the zeroed column and replace it by the i-th line of the identity matrix, I. If we also replace fi for
i ∈ ΓΩD, by the respective value of the Dirichlet function ~g, we naturally ensure (3.5.32). Naturally as we
saw before we also need a way to calculate∫
ΓΩN
~s · ~φi, and we do that by using a 1D Gaussian Quadrature
method:
23
∫ 1
−1
g(ε)dε =
N∑i=1
wig(εi) (3.5.33)
where g above will be ~s · ~φi in our case. The precision (N) of our method will be an available option inside
the function that actually computes this (N can go from 1 to 5).
3.5.5 Added Code
Here we just discuss what code was added to implement what has been described in the previous sections
(too see the full description of the original IFISS software, see [9] and [7]). Keep also in mind that the original
IFISS software solves our linear system derived from discretization by using matlab linearsolve function.
For other solving techniques, see [6], [19], [29], [28], [13] and [22]. One should note that all the code added
was added on top of the previous IFISS Q2-Q1 implementation. As usual, the full code regarding the P2-P1
solver for the Stokes problem can be found in the appendix. Before anything else, we store the boundary
types on a vector called outbc query. This is a 4 × 1 vector that stores a number which will identify the
i-th boundary with boundary conditions type outbc queryi. In main channel.m we have a few preset
types for the boundary:
FLOW CHANNEL NATURAL = [ 0 , 2 , 0 , 3 ] ;
FLOW CHANNEL = [ 0 , 2 , 0 , 2 ] ;
FLOW FULL DIR = [ 3 , 0 , 3 , 0 ] ;
FLOW CUSTOM = [ 2 , 0 , 2 , 0 ] ;
In addition to that, we also have a new variable called elements which will decide between choosing Q2-Q1(if
its set to 0) or P2-P1(if its set to 1). The code inside generate channel stk() is slightly modified to use
the new possible discretization type:
fpr intf ( ’ bu i l d i ng e lements . . . \ n ’ ) ;
[ x , y , xy , xyp ,mp, map ] = p2p1gridxx (x , y , xy ,mv, bound ) ;
fpr intf ( ’DOF: %d\n ’ , length ( xy ( : , 1 ) ) ) ;
fpr intf ( ’ completed .\n ’ ) ;
fpr intf ( ’ bu i ld system matrix , f s i d e . . . \ n ’ ) ;
[A,B, Bx , By , f , g ] = stk p2p1 ( xy , xyp ,mv,mp) ;
fpr intf ( ’ completed .\n ’ ) ;
These have identical meanings as in the Q2-Q1 case, therefore we won’t explain them any further. Except
for function stk p2p1 which is a bit different from stk q2q1 in the sense that now we are also implementing
the Neumann boundary conditions and this is achieved by function neumann right side 2 called inside
24
stk p2p1. As it was already told in the implementation sub-section, we can also define the how many nodes
we want for the boundary integral on∫
ΓΩN~s · ~φi - this is controlled by variable nngpt.
nngpt = 3 ; %depends on Neumann boundary f u n c t i o n .
i f nngpt==1
s (1 ) = . 0 ; wt (1 )=2 .0 ;
e l s e i f nngpt==2
s (1 ) = .577350269 ; wt (1 ) = 1 . 0 ;
s (2 ) = − .577350269; wt (2 ) = 1 . 0 ;
e l s e i f nngpt==3
s (1 ) = . 0 ; wt (1 ) = .888888889 ;
s (2 ) = .774596669 ; wt (2 ) = .555555556 ;
s (3 ) = − .774596669; wt (3 ) = .555555556 ;
e l s e i f nngpt==4
s (1 ) = .339981044 ; wt (1 ) = .652145155 ;
s (2 ) = − .339981044; wt (2 ) = .652145155 ;
s (3 ) = .861136312 ; wt (3 ) = .347854845 ;
s (4 ) = − .861136312; wt (4 ) = .347854845 ;
e l s e i f nngpt==5
s (1 ) = . 0 ; wt (1 ) = .568888889 ;
s (2 ) = .538469310 ; wt (2 ) = .478628670 ;
s (3 ) = − .538469310; wt (3 ) = .478628670 ;
s (4 ) = .906179846 ; wt (4 ) = .236926885 ;
s (5 ) = − .906179846; wt (5 ) = .236926885 ;
end
The other part of the code is very straight forward and it just implements the Gauss Quadrature integra-
tion method. The Neumann functions for each boundary, wich should be implemented by the user, are
NEUMANN F1, NEUMANN F2, NEUMANN F3 and NEUMANN F4. They refer to boundary
1,2,3 and 4 respectively. New functions such as :
[ jac , inv jac , phi , dphidx , dphidy ] = q d e r i v a t 2 ( s i gpt , t i gpt , x l v , y l v ) ;
and
[ ps i , dpsidx , dpsidy ] = qqde r iva t 2 ( s igpt , t i gpt , x l v , y l v ) ;
which are similar to what we have seen in the Q2-Q1 case in the sense that the first one should calculate
∂ψ∂x , ∂ψ∂y , ψ, |Jk| and |Jk|−1 for each element when we are processing Gauss Quadrature point in the reference
25
element (sk, tk). The same idea applies for the second function, with the difference that now we are dealing
with the bi-quadratic velocity functions φ ∈ L2. Other functions that are worth noting that have been added
are shap triangles 3nodes and q2shape triangles 6nodes which are the Q1 an Q2 shape functions,
respectively.
3.6 Numerical Results
In this section we shall give some examples of solutions obtained with the P2-P1 code along with some
Neumann-Dirichlet conditions enforced on the boundary.
3.6.1 Poiseuille Flow
We will follow with some tests on different kinds of reference problems. First we start with the standard
Poiseuille flow with full Dirichlet boundary conditions using as domain the square [−1, 1] × [−1, 1]. The
solution for the Poiseuille equation is the following:
ux = 1− y2, uy = 0, p = −2x+ cte
Picture 8 shows the result of our computed solution (using 66049 DOFs).
26
(a) u1 (b) u2
(c) Vector map
(d) Pressure
Figure 8: Poiseuille solution
But in order to evaluate a proper benchmark, we need to take into account several aspects: such as the time
each discrete system takes to be solved and the actual error that was obtained when compared with the real
solution. The errors will be computed using the following formula:
||eu|| =||u− u||||u||
(3.6.1)
where ||.|| is the euclidean norm.
In Figure 9 we show the relation of log10(h) used in the domain discretization and the log10 of the generated
relative error.
27
Table 1: Poiseuille Benchmark
DOF ||eu1 || time (seconds)25 4.7× 10−10 1.32× 10−2
81 1.6× 10−10 2.01× 10−3
289 4.59× 10−11 8.9× 10−3
1089 1.20× 10−11 4.04× 10−2
4225 3.09× 10−12 2.28× 10−1
16641 7.84× 10−13 1.39× 100
66049 2.01× 10−13 7.67× 100
Figure 9: Log-Log graph with linear regression for example 1
Figure 9 shows a plot of the results in table (1) using a ’Log-Log’ scale with a linear regression for this data
and a curve with slope 2. We can observe that both curves (green and red) have visually the same slope.
This is not in deed far from reality since we have observed that the calculated regression curve slope is 1.98.
This means that our error is near O(h2) as expected.
3.6.2 Cavity Flow
The next example we test is the problem of the enclosed flow inside a ’cavity’. In this example problem we
consider zero velocity on all the walls, except the upper one where we have lateral velocity (ex direction)
only. Once more we formulate this problem assuming Dirichlet boundary conditions only.
gy = 0 on ΓΩ (3.6.2)
gx = 0 on ΓΩ1,ΓΩ2
,ΓΩ4(3.6.3)
gx = 1− x4 on ΓΩ3 (3.6.4)
28
using as domain Ω = [−1, 1] × [−1, 1]. The Dirichlet boundary conditions imposed are present in (3.6.2),
(3.6.3) and (3.6.4). However, now we don’t have a way to test how accurate our solution is because we don’t
have an exact solution. However we can still try to see if the streamlines of our solution are as expected or
not. It is known that if we get enough discretization we should be able to see some flow recirculation on the
bottom corners.
(a) u1(b) u2
(c) Vector map(d) Pressure
Figure 10: Cavity solution
Figure 10 shows our solution at the cavity problem. But in order to see the recirculation we must choose
appropriated streamlines, such streamlines are visible in figure 11.
29
Figure 11: Cavity with streamlines
Since we don’t have the analytic expression for our solution this time, we can’t have an accurate way to test
if our solution is correct. However, we can compare it to previous works such as [8], pag. 220, and realize
the solutions look similar.
3.6.3 Custom Flow
The last flow will test our implementation of the Neumann boundary conditions. We will use the following
test flow:
ux = 20xy3, uy = 5x4 − 5y4, p = 60x2y − 20y3 + cte (3.6.5)
We will use Dirichlet boundary conditions on boundaries number 1 and 3 and the remaining will be Neumann
boundary conditions. Using the definition of the Neumann boundary function ~s, we know that:
~s2 = (−60x2y + 40y3, 20x3) (3.6.6)
~s4 = (60x2y − 40y3,−20x3) (3.6.7)
furthermore we use as domain Ω = [−1, 1]× [−1, 1]. Figure 12 shows the computed solution for this problem.
30
(a) u1
(b) u2
(c) Vector map (d) Pressure
Figure 12: Test solution
Once again, we compute the relative errors for the solution obtained.
Table 2: Test Flow Benchmark
DOF ||eu1|| ||eu2
|| time (seconds)25 5.0× 10−2 1.2× 10−1 4.14× 10−2
81 6.7× 10−3 1.19× 10−2 2.29× 10−3
289 5.84× 10−4 1.0× 10−4 1.094× 10−2
1089 4.45× 10−5 7.72× 10−5 4.13× 10−2
4225 3.11× 10−6 5.40× 10−6 2.12× 10−1
16641 2.04× 10−7 3.70× 10−7 1.249× 100
66049 1.36× 10−8 2.45× 10−8 7.35× 100
31
In addition to the previously presented tests, we will also add an obstacle based test. In this particular test,
we will be using full Dirichlet boundary conditions (for the Poiseuille flow), except for boundary layer number
2 where we will be using homogeneous Neumann conditions.
Figure 13: Stokes Poiseuille problem with obstacle solution
(a) ux
Figure 14: Stokes Poiseuille problem with obstacle ux solution graph
Figure 13 reveals the complete solution details, while figure 14 reveals only an enlarged version of the solution
the ux solution part.
32
4 Stationary Navier-Stokes Equations
4.1 Description
In this section we address the problem described in sub-section 2.3. Our main equations will be (2.3.1) and
(2.3.2) using as constraints (2.3.3), (2.3.4) and (2.3.5). Two different approaches were taken regarding solving
this problem. First we explain and apply the Newton method to our problem. Next we use our algorithms
on a few standard test cases and analyze its results. Next we apply a stabilization technique (SUPG) to
our problem and a stability result is given and proven. Finally we test this stabilization technique on more
standard test cases and analyze the results.
4.2 Newton Method
To solve the Navier-Stokes Equations we try to proceed as we did before, that is, first we write the cor-
respondent weak formulation of the Navier-Stokes problem (however, for different Navier-Stokes equation
formulations, one can get different weak formulations - for additional information see [16], [33]). This process
is very similar to what we did before for the Stokes problem, therefore we avoid showing all the intermediate
steps and we just show the final form (for details see [1], [10], [15], [24] and [30]), which is: find ~u ∈ H1E and
p ∈ L2(Ω):
ν
∫Ω
∇~u : ∇~v − ν∫
ΓΩN
(~n · ∇~u) · ~v +
∫Ω
(~u · ∇)~u · ~v +
∫Ω
p(∇ · ~v) =
∫Ω
~f · ~v, ∀~v ∈ H1E0
(4.2.1)∫Ω
q(∇ · ~u) = 0, ∀q ∈ L2(Ω) (4.2.2)
Unfortunately we no longer obtain a linear system that can easily be solved, and this happens because of the
term (~u · ∇)~u. In order to surpass this challenge, we have chosen to apply Newton’s method to the problem.
In this section we will describe Newton’s method, and how it can be applied to this problem. First, consider
the following operator F : H1E × L2(Ω)→ R:
F (~u, p) = −ν∆~u+ (~u · ∇)~u+∇p− ~f (4.2.3)
Now we can calculate its Gateaux derivatives:
33
DF (~u, p) · (~v, q) = limε→0
1
ε(F (~u+ ε~v, p+ εq)− F (~u, p)) =
limε→0
1
ε(−ν∆(~u+ ε~v) + ν∆(~u) + ((~u+ ε~v) · ∇)(~u+ ε~v)− (~u · ∇)~u+∇(p+ εq)−∇p− ~f + ~f) =
−ν∆~v + (~v · ∇)~u+ (~u · ∇)~v +∇q
Therefore we know that the derivative of operator F on points (~u, p) under the ’direction’ (~v, q) is:
DF (~u, p) · (~v, q) = −ν∆~v + (~v · ∇)~u+ (~u · ∇)~v +∇q (4.2.4)
The point is that we are seeking solutions for F (~u, p) = 0. Therefore by, applying Newton’s Method, we
have:
F (~uk+1, pk+1) = 0 u DF (~uk, pk) · (δ~uk, δpk) + F (~uk, pk)
where δ~uk = uk+1 − uk and δ~pk = pk+1 − pk. This entails the following numerical method:
DF (~uk, pk) · (δ~uk, δpk) = −F (~uk, pk) (4.2.5)
~uk+1 = ~uk + δ~uk (4.2.6)
pk+1 = pk + δpk (4.2.7)
Notice now that system (4.2.5) is potentially linear. From now on, we refer to the right hand side of (4.2.5)
as ’residual’, since as Newton is converging this member is supposed to have less and less contribution for
the final solution, and for this reason we will use it as stopping criteria further on. Now we do again a weak
formulation for (4.2.5) in order to extract a linear system that we can later on solve iteratively. We start
with the residual part:
Rei =
∫Ω
F ( ~uk, pk) · ~φi
where ~φi are base functions for the space H1E0
. And by expanding it:
34
Rei =
∫Ω
−ν∆~uk · ~φi + ~uk · ∇~uk · ~φi +∇pk · ~φi − ~f · ~φi (4.2.8)
Integration by parts on∫
Ω∆~uk · ~φi and
∫Ω∇pk · ~φi yields:
∫Ω
∆~uk · ~φi =
∫ΓΩN
~n · ∇~u · ~φi −∫
Ω
∇~uk : ∇~φi (4.2.9)∫Ω
∇pk · ~φi =
∫ΓΩN
~n.pk · ~φi −∫
Ω
pk · ∇ · ~φi (4.2.10)
By substituting into (4.2.8), we obtain:
Rei = ν
∫Ω
∇~uk : ∇~φi +
∫Ω
~uk · ∇~uk · ~φi −∫
Ω
pk · ∇~φi −∫
Ω
~f · ~φi − ν∫
ΓΩN
~n · ∇~uk · ~φi +
∫ΓΩN
~n.pk · ~φi =
ν
∫Ω
∇~uk : ∇~φi +
∫Ω
~uk · ∇~uk · ~φi −∫
Ω
pk · ∇~φi −∫
Ω
~f · ~φi −∫
ΓΩN
~s · ~φi
And therefore we obtain the weak form for (4.2.8):
Rei = ν
∫Ω
∇~uk : ∇~φi +
∫Ω
~uk · ∇~uk · ~φi −∫
Ω
pk · ∇~φi −∫
Ω
~f · ~φi −∫
ΓΩN
~s · ~φi (4.2.11)
In order to have a complete system we need to repeat the process but for the left hand side of (4.2.5). We
start then by expanding∫
ΩDF (~uk, pk) · (δ~uk, δpk) · ~φi:
∫Ω
DF (~uk, pk) · (δ~uk, δpk) · ~φi = −ν∫
Ω
∆δ~uk · ~φi +
∫Ω
δ~uk · ∇~uk · ~φi +
∫Ω
~uk · ∇δ~uk · ~φi +
∫Ω
∇δpk · ~φi =
−ν(−∫
Ω
∇δ~uk : ∇~v) +
∫Ω
δ~uk · ∇~uk · ~φi +
∫Ω
~uk · ∇δ~uk · ~φi −∫
Ω
δpk · ∇ · ~φi
This allows us to obtain the following discrete system:
νA+N +W BT
B O
.
δuδp
=
fg
(4.2.12)
where A and B are defined as in the Stokes discrete formulation (3.3.12) and (3.3.13). The rest are defined
as:
35
fi = Rei (4.2.13)
gk =
∫Ω
ψk(∇ · ~uk) (4.2.14)
Nij =
∫Ω
(~uk · ∇φj) · ~φi (4.2.15)
Wij =
∫Ω
(~φj · ∇~uk) · ~φi (4.2.16)
As was previously done in (3.3.16), we must take also into account the two dimension version of (4.2.12),
which is:
νA+N +Wxx Wxy BTx
Wyx νA+N +Wyy BTy
Bx By O
.
δux
δuy
δp
=
fx
fy
g
(4.2.17)
where:
Nij =
∫Ω
(~uk · ∇φj)φi (4.2.18)
Wxy,ij =
∫Ω
∂ux∂y
φiφj (4.2.19)
However, in order to apply stabilization we need a slightly different formulation of the previously described
Newton method. The reason for this is that stabilization algorithms are usually applied to a bilinear con-
tinuous coercive operator a(., .) linked to a problem of the kind a(uk+1, pk+1) = fk (in our case we simply
have a(δuk, δpk) = fk). So instead of adapting the stabilization operator for this Newton formulation,
we have decided to adapt this Newton formulation for the SUPG operator on Navier-Stokes. This means
we are simply reformulating the previous method to have something like A · (u, p) = f . First notice that
DF (~uk, pk) · (δ~uk, δpk) is a linear operator as it would be expected. For this reason:
DF (~uk, pk) · (δ~uk, δpk) = DF (~uk, pk) · (~uk+1, pk+1)−DF (~uk, pk) · (~uk, pk) (4.2.20)
Thus, we can reformulation Newton iteration method as:
DF (~uk, pk) · (~uk+1, pk+1) = −F (~uk, pk) +DF (~uk, pk) · (~uk, pk) (4.2.21)
36
And by expanding the right-hand side we obtain:
−F (~uk, pk) +DF (~uk, pk) · (~uk, pk) = (4.2.22)
−(−ν∆~uk + ~uk · ∇~uk +∇pk − f) + (−ν∆~uk + ~uk · ∇~uk + ~uk · ∇~uk +∇pk) = f + ~uk · ∇~uk (4.2.23)
Now we have enough ingredients to expand (4.2.21) and attempt a weak formulation for this formula.
−ν∫
Ω
∆~uk+1 · φi +
∫Ω
~uk+1 · ∇~uk · φi +
∫Ω
~uk · ∇~uk+1 · φi +
∫Ω
∇pk+1 · φi =
∫Ω
f · φi +
∫Ω
~uk · ∇~uk · φi
(4.2.24)
Using the following formulas:
∫Ω
∆~uk+1 · φi +
∫Ω
∇~uk+1 : ∇~φi =
∫ΓΩN
~n · ∇~uk+1 · φi (4.2.25)∫Ω
∇pk+1 · φi +
∫Ω
pk+1 · ∇ · φi =
∫ΓΩN
~n.pk+1 · φi (4.2.26)
We obtain the final weak formulation:
ν
∫Ω
∇~uk+1 : ∇φi +
∫Ω
~uk+1 · ∇~uk · φi +
∫Ω
~uk · ∇~uk+1 · φi +
∫Ω
∇pk+1 · φi = (4.2.27)∫Ω
f · φi +
∫Ω
~uk · ∇~uk · φi +
∫ΓΩN
~s · φi (4.2.28)
∀φi ∈ H1E0
(4.2.29)
After the discretization process is complete, we obtain the following linear system:
νA+N +Wxx Wxy BTx
Wyx νA+N +Wyy BTy
Bx By O
.
ux,k+1
uy,k+1
pk+1
=
Fx
Fy
g
+
Wxx Wxy O
Wyx Wyy O
O O O
.
ux,k
uy,k
pk
(4.2.30)
where :
37
Fx =
∫Ω
fx · φi +
∫ΓΩN
~sx · φi (4.2.31)
Fy =
∫Ω
fy · φi +
∫ΓΩN
~sy · φi (4.2.32)
where φi is as usual a basis for the discrete space Xh0 adopted. And the other matrix components are the
same as before. Naturally we need to agree on a stopping criteria for these methods as well. In our case we
will simply control the quantity:
||δ ~uk||H1E(Ω) = (||∇δ ~uk||2L2(Ω) + ||δ ~uk||2L2(Ω))
1/2 (4.2.33)
We also need to adopt a discrete version for these norms as well. We will simply calculate (diffT ·A ·diff −
diffT · r · diff)1/2, where diff = uk+1 − uk and r is a matrix such that ri,j =∫
Ωφi · φj . This will be our
stopping criteria to control the iteration error.
4.3 Numerical Results for Newton Method
4.3.1 Poiseuille Flow
For the fully implicit Newton method we choose for the first example the typical Poiseuille flow (with low
Reynolds number) where we shall construct the same table as we did before for the case of Stokes. More
specifically we are using Dirichlet boundary conditions on boundary 1,3 and 4 and neutral Neumann condi-
tions on boundary 2. The expression used on the Dirichlet boundary data is:
ux = 1− y2 (4.3.1)
uy = 0 (4.3.2)
Additional data used on our Newton method was:
1. Starting point: null solution
2. Stopping criteria: iteration error < 1× 10−5
3. ν = 1/10 (m2/s)
A solution with high-resolution is visible in figure 15.
38
Figure 15: Navier-Stokes solution of Poiseuille flow (Newton)
Notice how the pressure slope is around 0.2 as it was supposed to be in the Navier-Stokes case of the Poiseuille
flow (p = −2xν). Next we present a table with the benchmarks for the several domain subdivision levels for
the previous problem.
Table 3: Navier-Stokes Poiseuille Benchmark
DOF ||eu1|| time (seconds)
25 4.48× 10−10 1.1681 1.72× 10−10 1.30289 5.80× 10−11 1.201089 1.94× 10−11 1.374225 6.60× 10−12 2.2216641 2.282× 10−12 7.1766049 8.29× 10−13 38.74
In a similar way to what has been done before for the Poiseuille flow of the Stokes problem, figure 16 shows
a plot of the data in table (3), a linear regression for this data and a linear function with slope 2 (for
comparison). The linear regression function we calculated has a slope value of 1.7 which is relatively close to
the expected value of 2 for a convergence of O(h2) (however since Newton method requires several iterations
a worse error convergence is expected when compared with FEM for the Stokes problem due to possible error
propagation between iterations).
39
Figure 16: Log-Log graph with linear regression for example 1
4.3.2 Cavity Flow
The next example is a bit more interesting: we will be testing the full Dirichlet cavity fluid (the same
conditions we enforced in the Stokes examples) for different values of ν and testing whether we are able to
get convergence for each case and under which conditions (no convergence is when we cannot get a solution
within the error bonds of < 1× 10−5 under 100 iterations with no more than 16641 elements). This will be
a valuable resource when comparing with the stabilized (SUPG) version of the next section. Keep also in
mind that we begin our attempts with 1089 degrees of freedom.
Table 4: Cavity Convergence Benchmark
ν (m2/s) convergence (DOF, iterations)10 (1089, 4)1 (1089, 4)1/10 (1089, 5)1/100 (1089, 7)1/1000 -
Figure 17 shows a high definition solution streamlines for the Navier-Stokes problem of the previous cavity
formulation (with ν = 1/100 (m2/s)).
40
Figure 17: Navier-Stokes cavity solution streamlines (Newton, ν = 1/100 m2/s)
The respective full solution can be seen in Figure 18. For other experiments on similar problems and its
results see [14].
Figure 18: Navier-Stokes cavity solution (Newton)
The last example includes another Poiseuille instance with an obstacle similar to what we did for the Stokes
case. Since we are not using any stabilization, we will just use ν = 1/100 m2/s.
41
Figure 19: Navier-Stokes Poiseuille with obstacle solution (Newton)
(a) Streamlines (b) ux
Figure 20: Navier-Stokes Poiseuille with obstacle solution streamlines and enlarged ux graph (Newton)
Figure 19 refers to the full solution of this problem while 20 refers to the streamlines plus a enlarged version
of a ux plot.
4.4 Newton SUPG Method
In order to reduce the interpolation errors, one often uses a stabilization term on the weak problem, more on
this topic can be found in [17], section 2 and [24] sections 8.3.1 and 8.3.2. For stabilization on other problems,
42
such as streamline-diffusion or advection-diffusion, see [18], [34], [26] or [3]. In our case what we will do is to
consider first the weak formulation of the Newton method for the stationary Navier-Stokes equations: Find
~uk+1 ∈ H1E , p
k+1 ∈ L2:
a(~uk+1, pk+1, ~v) = F (~v), ∀~v ∈ H1E0
(4.4.1)
where:
a(~uk+1, pk+1, ~v) = −ν(∆~uk+1, ~v) + (~uk+1 · ∇~uk, ~v) + (~uk · ∇~uk+1, ~v) + (∇pk+1, ~v) (4.4.2)
F (~v) = (~f,~v) + (~uk · ∇~uk, ~v) (4.4.3)
The SUPG method consists of adding a stabilization operator to formulation (4.4.1), L(~w, p), which will be
zero if we evaluate it on the actual solution. The following is the adopted formulation for our stabilization
operator (SUPG case), for other possible stabilization formulas, see [17]:
L(~w, p) = δ∑k∈τh
(L(~w, p)−~b, h2K .LSS(~v, q))K (4.4.4)
where:
L(~v, q) =
−ν∆~v + (~uk · ∇)~v + (~v · ∇)~uk +∇q
∇ · ~v
(4.4.5)
~b =
~f + (~uk · ∇)~uk
0
(4.4.6)
And LSS is the skew-symmetric associate of the operator L:
LSS(~v, q) =
−ν(~uk · ∇)~v − 12 (∇ · ~uk)~v +RSS~v +∇q
∇ · ~v
(4.4.7)
where RSS~v = 12 [(~v · ∇) · ~uk − (∇~uk) · ~v]. Notice that the original weak formulation (4.4.1) can be re-written
as L(~uk+1, pk+1)−~b = 0, therefore, we can finally put together the stabilized version of (4.4.1):
43
a(~uk+1, pk+1, ~v) + L(~uk+1, pk+1) = F (~v), ∀~v ∈ H1E0
(4.4.8)
Notice that if the pair (~uk+1, pk+1) happens to be the real solution, then operator L is zero. Therefore, for the
real solution, we recover (4.4.1). The challenge is to find a good value for δ in formula (4.4.4). The optimal
value for the δ seems to vary significantly according to the problem type and the amount of discretization
used. However if we choose a value for δ that is too far from the optimal value, SUPG seems to take too
long to converge (and does not converge in a realistic amount of time). We have failed to come up with an
heuristic that behaves optimally for all problems and all discretization levels, however we did find an heuristic
that seems to behave properly for our examples (but it has limitations) as we will see in the next section.
Our heuristic (found only via trial-and-error) is:
δ =6.4× 10−4
ν.h2K
(4.4.9)
however, other heuristics to calculate δ can be found in [23]. Next we will prove a stability estimate for our
previous method. In this proof, we adapt the idea used to prove the stability of the stabilized method for
the convection-diffusion problem described in [23], page 321 (for more information on convection-diffusion
problems see [27]). For additional theory on stability addressing Navier-Stokes equations, see [5]. Before
introducing main result, we introduce a custom norm called SUPG norm. It is defined as follows:
Definition 4.4.1 Norm ||.||SUPG:
||~u||SUPG = ν||∇~u||2L2 + (~b · ∇~u, ~u) + (~u∇~b, ~u) + δ∑k∈τh
(L(~u, p), h2K .LSS(~u, p))K + (∇p, ~u)
for ~u a certain function in H1E.
The following is our stability result for the Newton SUPG method.
Theorem 4.4.2 Let ~u be the ~uk+1 ∈ H1E0
(simplification) solution in (4.4.8), ~b the previous ~uk solution, ~v a
test function ∈ H1E0
, q ∈ L2 and finally LSS and L are the operators defined in (4.4.5) and (4.4.7). Assume
that ~uT · ∇~b.~u ≥ 0 and that (LSS(~v, q), LS(~v, q))K is also 0. Then:
||~u||SUPG ≤4
3(ε1 + 1)||~g||2L2 (4.4.10)
where ε is a positive constant and ~g = ~f +~b · ∇~b.
Proof: By the antisymmetry property of the non-linear operator, we know that: (~b · ∇~u, ~u) = 0. And since,
by assumption (LSS(~v, q), LS(~v, q))K = 0 , we can simplify ||~u||SUPG expression into:
44
||~u||SUPG = ν||∇~u||2L2 + (~u∇~b, ~u) + δ∑k∈τh
(LSS(~u, p), h2K .LSS(~u, p))K + (∇p, ~u) (4.4.11)
Since ~u is a solution to our problem and assuming ~u ∈ H1E0
, we can know that: (∇p, ~u) = 0 - by using
integration by parts. Our formula therefore simplifies to:
||~u||SUPG = ν||∇~u||2L2 + (~u∇~b, ~u) + δ∑k∈τh
(LSS(~u, p), h2K .LSS(~u, p))K (4.4.12)
By definition of our SUPG method, we know that, furthermore:
||~u||SUPG =
∫Ω
~g.~u+∑k∈τh
(~g, h2K .LSS(~u, p))K (4.4.13)
where ~g = ~f +~b · ∇~b. We can also observe the following:
∫Ω
~g.~u ≤ ||~g||2L2 .||~u||2L2 ≤ ε1||~g||2L2 +1
4 · ε1||~u||2L2 (4.4.14)
By using Cauchy-Schwartz and Young’s inequalities. Furthermore we can additionally use Poincare’s in-
equality to conclude that:
||~u||2L2 ≤ C.||∇~u||2L2 (4.4.15)
Therefore, expression (4.4.14) can be simplified to :
∫Ω
~g.~u ≤ ||~g||2L2 .||~u||2L2 ≤ ε1||~g||2L2 +C
4 · ε1||∇~u||2L2 (4.4.16)
The other term of the right-hand side of (4.4.13), we can also majored:
∑k∈τh
(~g, h2K .LSS(~u, p))K ≤ δ||~g||2L2 .||h2
K .LSS(~u, p)||2L2 ≤ ε2||~g||2L2 +1
4 · ε2||h2
K .LSS(~u, p)||2L2 (4.4.17)
By choosing ε1 such that C4·ε1 = ν
4 and ε2 = 1, we have:
45
||~u||SUPG ≤ ε1||~g||2L2 +ν
4||∇~u||2L2 + ||~g||2L2 +
1
4||h2
K .LSS(~u, p)||2L2 (4.4.18)
Since, by assumption, (~u∇~b, ~u) ≥ 0, we can simply add 14 (~u∇~b, ~u) to our previous expression, getting:
||~u||SUPG ≤ ε1||~g||2L2 +ν
4||∇~u||2L2 + ||~g||2L2 +
1
4||h2
K .LSS(~u, p)||2L2 +1
4(~u∇~b, ~u) (4.4.19)
By substitution, we have finally:
||~u||SUPG ≤ (ε1 + 1)||~g||2L2 +1
4||~u||SUPG (4.4.20)
3
4||~u||SUPG ≤ (ε1 + 1)||~g||2L2 (4.4.21)
||~u||SUPG ≤4
3(ε1 + 1)||~g||2L2 (4.4.22)
which proves our point.
4.5 Numerical Results for SUPG
4.5.1 Poiseuille Flow
The examples for this section will essentially follow the same structure of the previously explored sections,
therefore we will start with a Poiseuille flow whose configuration is Dirichlet boundary conditions on boundary
layers 1,3 and 4 and neutral Neumann conditions on boundary layer 2. Similarly to what we have also done
before, we will also present a benchmark table that takes into account execution times and relative errors
with respect to the real solution. The additional settings regarding to these tests will be very similar to what
we had set up in the Numerical Methods section regarding the Newton method: we will be using as starting
point the null solution, we will also be using the same iteration error 1 × 10−5 as stopping criteria but we
will be using ν = 1/1000 instead. The table (5) shows our results for this case.
Aditional data relating h used in the discretization process and the generated relative error is plotted in
Figure 21.
46
Table 5: Navier-Stokes Poiseuille Benchmark (Fully implicit Newton-SUPG)
DOF ||eu1 || time (seconds)81 3.3× 10−4 1.07289 5.31× 10−4 1.661089 6.24× 10−4 1.634225 6.57× 10−4 2.916641 6.65× 10−4 10.0766049 6.64× 10−4 51.7
Figure 21: Log-Log graph between ||eu1 || and h
The δ values chosen are generated from the previously explained heuristic. Figure 21 shows that our error
results are nothing like we would expect them to be. One possible justification for this fact is that our δ
choice is probably not the best for this particular problem. We also provide a high resolution solution for the
previous problem (using SUPG fully implicit Newton) on figure 22.
47
Figure 22: Navier-Stokes Poiseuille solution (Newton SUPG)
4.5.2 Cavity Flow
The next example will follow the same pattern as the Newton numerical results section, that is, we will
first present a table similar to table (3) but in this case using SUPG stabilization (and guessing δ with the
previously presented heuristic). The iteration stopping criteria and starting point will remain the same as
in the previous example. The problem will be formulated with full Dirichlet boundary conditions respecting
the following constraints:
ux = 1− x2, uy = 0, on boundary layer 3 (4.5.1)
ux = uy = 0, on any other boundary layer (4.5.2)
Table (6) shows our results.
We must note that although the last table entry didn’t converge according to our convergence definition, we
got a error margin of 0.99 on the last iteration and it had just been starting to fastly converge. The results
seem very positive however for the table values of ν ≤ 1/1000 (m2/s) we only get convergence for DOF=1089.
This may seem odd, but the reader should be remembered that our heuristic depends on the macro-elements
as well, which means that the chosen heuristic cannot handle all situations properly. In addition we also
provide a high definition solution for the previous cavity problem with ν = 1/1000 m2/s and a low definition
48
Table 6: Cavity (SUPG fully implicit Newton) Convergence Benchmark
ν (m2/s) convergence (DOF, iterations)1× 101 (1089, 4)1× 100 (1089, 5)1× 10−1 (1089, 6)1× 10−2 (1089, 9)1× 10−3 (1089, 13)1× 10−4 (1089, 52)1× 10−5 (1089, 46)1× 10−7 (1089, 45)1× 10−10 (1089, 48)1× 10−13 (1089, 47)1× 10−17 -
streamline solution for ν = 1× 10−13 m2/s which was the last solution achieved in table (6).
Figure 23 refers to the full solution for cavity (using SUPG, naturally) when ν = 1/1000 m2/s. Figure 24
shows the two streamline solutions we just discussed.
Figure 23: Navier-Stokes cavity solution (Newton SUPG) with ν = 1/1000 m2/s
49
(a) ν = 1× 10−3 m2/s (b) ν = 1× 10−13 m2/s
Figure 24: Navier-Stokes cavity solution streamlines using SUPG for ν = 1/1000 m2/s and ν = 1× 10−13 m2/s
One of the reasons why we are not getting any recirculation on the streamlines for the cavity solution with the
lower ν is that probably we do not capture enough information from the solution due to the bad discretization
used. The last test we will perform is , as previously, a Poiseuille flow with an obstacle inside the domain. The
formulation adopted is again Dirichlet conditions on boundary layers 1,3 and 4 followed by neutral Neumann
conditions on layer number 2. For this test we attempted to compute a solution with ν = 1/1000 m2/s, but
unfortunately we had a small setback which was that we couldn’t find a suitable value for δ on SUPG for a
high discretization detail (DOF = 66049).However, we were able to calculate a solution for DOF=16641.
50
Figure 25: Navier-Stokes Poiseuille with obstacle solution (Newton SUPG)
(a) Streamlines (b) ux
Figure 26: Navier-Stokes Poiseuille with obstacle solution streamlines and enlarged ux graph (Newton SUPG)
Figure 25 shows the full solution details, while Figure 26 shows the respective streamlines along with an
enlarged version ux plot.
51
5 Non-Stationary Navier-Stokes Equations
In this section, we will explain our method to solve equation (2.2.7) numerically: its weak formulation and
our numerical methods. Farther on, we will also present some examples. As what has been exposed in
sub-section 2.2, we know that the non-stationary form of the Navier-Stokes Equations is:
∂~u
∂t+ (~u · ∇)~u− ν∆~u+∇p = ~f in Ω× (0,∞) (5.0.1)
∇.u = 0 in Ω× (0,∞) (5.0.2)
5.1 Numerical Approximation
Our approach will be based on what has been explained in [4], [23] and [24], more precisely, we will use a
fully implicit finite differences approximation scheme in order to time. In our simulation we will consider a
simulation interval I = [0, T ] divided in n subintervals with measure ∆t = tn+1 − tn. We use the following
first order approximation for the ~u derivative:
∂~u(tn+1)
∂t=~u(tn+1)− ~u(tn)
∆t+O(∆t) (5.1.1)
Therefore assuming the previous approximation for the derivative ∂~u∂t and substituting into the original non-
stationary Navier-Stokes equations we get the following system:
un+1 − un
∆t+ (un+1 · ∇)un+1 +∇pn+1 − ν∆un+1 = fn+1 in Ω (5.1.2)
∇.un+1 = 0 in Ω (5.1.3)
where uk = ~u(tk), pk = p(tk) and fk = ~f(tk) .In each iterations we need to solve the previous (non-linear)
system and we also need to define what solution will be u0. As far as u0 is concerned, we will use the
stationary solution for the Navier-Stokes problem. In what regards solving the non linear system (5.1.2) on
each iterations we will be using the previously explained Newton method. Here we also present a stability
estimate for this discretization. Consider the original (continuous) problem:
δ~u
δt− ν∆~u+ (~u · ∇~u) +∇p = ~f (5.1.4)
∇.~u = 0 (5.1.5)
~u = ~g on ΓΩD, ~u = ~u0 in Ω, t = 0 (5.1.6)
52
Now consider a discrete version of the previous problem:
1
∆t(~un+1 − ~un)− ν∆~un+1 +∇pn+1 + (~un+1 · ∇~un+1) = ~fn+1 (5.1.7)
∇.~un+1 = 0 (5.1.8)
assuming, of course, the following approximation for the time derivative:
∂~u
∂t(tn+1) u
1
∆t(~un+1 − ~un) (5.1.9)
Our weak formulation for the previous discrete version reads as follows: ∀n ≥ 0, find ~un+1 ∈ H1E(Ω),
pn+1 ∈ L20(Ω), for any ~v ∈ H1
E0(Ω) and q ∈ L2(Ω):
1
∆t(~un+1, ~v) + a(~un+1, ~v) + b(~v, pn+1) + c(~un+1, ~un+1, ~v) = (~fn+1, ~v) +
1
∆t(~un, ~v) (5.1.10)
b(~un+1, q) = 0 (5.1.11)
We obtain the stability estimate by setting v = ~un+1 and q = pn+1, and subtracting the continuity equation
from the momentum one:
1
∆t(~un+1 − ~un, ~un+1) + a(~un+1, ~un+1) + c(~un+1, ~un+1, ~un+1) = (~fn+1, ~un+1) (5.1.12)
As stated in pag. 346 of [11], the non-linear term c(., ., .) is skew-symmetric. Therefore since c(~u, ~z, ~w) =
−c(~z, ~u, ~w), we know that c(~un+1, ~un+1, ~un+1) = 0, which allows us to remove the non-linear term from the
previous formula. Since we know the bilinear form a is coercive, we can (also by applying Young inequality)
get the following expression from the one above:
1
∆t(~un+1 − ~un, ~un+1) + α||~un+1||2V ≤
1
4ε||F ||2V ′ + ε||~un+1||2V (5.1.13)
Denoting by Fn+1(~v) = (fn+1, ~v). The following result is known and can be seen in section 5 of [11]:
2(~un+1 − ~un, ~un+1) = ||~un+1||2L2 + ||~un+1 − ~un||2L2 − ||~un||2L2 (5.1.14)
By setting ε = α/2 we get:
53
1
∆t(||~un+1||2L2 + ||~un+1 − ~un||2L2 − ||~un||2L2) + α||~un+1||2V ≤
1
4ε||F ||2V ′ (5.1.15)
Multiplying by ∆t, and summing the time index n = 0, ...N − 1 (where N is the number of discrete time
points):
N−1∑n=0
(||~un+1||2L2 − ||~un||2L2) = ||~uN ||2L2 − ||~u0||2L2 (5.1.16)
and taking advantage of the telescopic sum we get at last:
||~uN ||2L2 + α∆t
N−1∑n=0
||~un+1||2V ≤∆t
α
N−1∑n=0
||Fn+1||2V ′ + ||~u0||2L2 (5.1.17)
which is a stability estimate for the original time dependent Navier Stokes problem.
5.2 Numerical Results
The first experiment we will consider is a very specific test flow we designed only for testing. The idea is to
first choose a ~u solution that respects ∇·~u = 0 and secondly calculate ~f from the previously set ~u expression
(and p as well). The expression chosen for ~u is :
5.2.1 Test Example
ux = sin(π.y).sin(π.t/2) (5.2.1)
uy = sin(π.x).sin(π.t/2) (5.2.2)
Therefore, for t = 0 we have no flow inside, and for t = 1 we will have the maximum amplitude for both ux
and uy. As for the pressure we set:
p = 2x (5.2.3)
Naturally, by doing the math (in this case we use Mathematica for the tedious symbolic calculations), we
obtain the following formula for ~f :
54
fx = 2 + pi ∗ cos(pi ∗ y) ∗ (sin(pi ∗ t/2)2) ∗ sin(pi ∗ x)+
0.5 ∗ pi ∗ cos(pi ∗ t/2) ∗ sin(pi ∗ y) + (pi2) ∗ v ∗ sin(pi ∗ t/2) ∗ sin(pi ∗ y);
fy = 0.5 ∗ pi ∗ (cos(pi ∗ t/2) ∗ sin(pi ∗ x) + 2 ∗ sin(pi ∗ t/2) ∗ (pi ∗ v ∗ sin(pi ∗ x)
+cos(pi ∗ x) ∗ sin(pi ∗ t/2) ∗ sin(pi ∗ y)));
We also need to define a way of calculating the error. We will be adopting the following norm:
||δ~u||L(0;T,H1E(Ω)) = (
∫ t=T
t=0
||δ~u(t)||2H1E(Ω))
12 u (
t=T∑t=0,t=t+∆t
∆t.||δ~u(t)||2H1E(Ω))
12 (5.2.4)
which means that, in our situation, we are having a nested combination of approximations for calculating
the error relative do the real solution.
The following table expresses the errors we have obtained along with the corresponding time discretization.
∆t denotes the time discretization interval between t = 0 and t = 1 and GRE denotes the general relative
error given by ||~u− ~ur||L(0;T,H1E(Ω))/||~ur||L(0;T,H1
E(Ω)), where ~ur and ~u denote the real and the approximated
solution, respectively.
Table 7: NSNS Benchmark
∆t GRE1/10 0.05011/20 0.01221/40 0.00301/80 7.36× 10−4
1/160 1.81× 10−4
Figure 27 shows the graph for relating the log10(∆t) versus the log10 of the GRE value. Additionally a best
fit linear regression and a linear function with slop 2 have been plotted.
55
Figure 27: Log-Log graph with linear regression for example 1
Visually, by observing Figure 27, we can see that the error is very close to O(h2). This is a bit better then
what was expected since we are only using a first order approximation in formula (5.1.1) (in fact the slope
for the linear regression curve generated for our error points is around 2.02).
5.2.2 Cavity Test
Another experiment we did was to include a non-stationary cavity flow, whose intensity varies from 0 to 1
(at t = 0.5) and later on it goes back to 0 again (at t = 1) using a sin() function. Interestingly enough, we
noticed that if we use enough discretization we can achieve this simulation even for values of ν = 1/1000
m2/s (which is this case). The following are the details of the simulation:
uy = 0, on ΓΩ (5.2.5)
ux = 0, on ΓΩ1,ΓΩ2
,ΓΩ4(5.2.6)
ux = (1− x2)|sin(πt
2)|, on ΓΩ3 (5.2.7)
We considered, as previously, the starting solution u0 (u evaluated at t = 0) to be the solution generated
by the associated stationary problem (which we assume to be the zero solution since no external forces are
applied). In this simulation t ∈ [0, 1], ∆t = 140 and each stationary NS problem was solved using 66049
degrees of freedom (which corresponds to a subdivision level of 8). Figure 28 is the final solution obtained
at t = 1.
56
Figure 28: Non-stationary Navier-Stokes cavity simulation (fully implict)
(a) t = 0.25(b) t = 0.5
(c) t = 0.75 (d) t = 1
Figure 29: Non-stationary Navier-Stokes cavity streamline simulation (fully implicit)
5.2.3 Womersley Test
Additionally, other known test we did was to solve a flow known as ’Womersley’. This flow is very known in
the simulation scene. The original formulation of the flow is described as mixed boundary problem (Dirichlet
57
/ Neumann) as follows:
uy = 0 as a Dirichlet boundary condition on boundaries 1 and 3 (5.2.8)
s = 2νsin(w.t)n as a Neumann boundary condition on boundary 4 (5.2.9)
s = 0 as a Neumann boundary condition on boundary 2 (5.2.10)
which is the exact formulation we will use in our simulation (previously we attempted to use the derivation
presented in page 380 of [11] however, we encountered performance difficulties regarding the approximation
of the analytic ux solution). In the following simulation results, we use ν = 1/10 m2/s and w = 3.14.
(a) t = 1.01 (b) t = 1.40
(c) t = 1.55
(d) t = 2.01
Figure 30: Womersley flow solution with ν = 1/10 m2/s
However, one can’t simply expect to see the exact same results as in [11] because the domain is different
(here we use [−1, 1]× [−1, 1] instead of [0, 1]× [0, 1]).
58
5.2.4 Obstacle Flow Test
In this section we perform a more complicated (and much more interesting) test. This time we are attempting
to see fluid recirculation on a obstacle based non-stationary Navier Stokes setup. Our domain is the same as
before (the [−1, 1]× [−1, 1] square), however this time, we set up a small square obstacle sufficiently close to
the left hand side boundary (ΓΩ4). To see some fluid recirculation we must have a sufficiently high Reynolds
number, and to accomplish this, we have combined both high speed entry on the boundary and low ν. To
be more exact, our boundary types are: Dirichlet in ΓΩ1 ,ΓΩ3 and ΓΩ4 and null Neumann on ΓΩ2 . Our ν is
1/1000 (m2/s) and our t goes from 0 to 4. The following is our Dirichlet boundary condition expression:
ux(t) = (1− y4) · ( tn
tn + 1) · 7 (5.2.11)
uy(t) = 0 (5.2.12)
(5.2.13)
With n = 6 in our particular case. With this expression, we attempt to simulate a water pipe that seems
closed at first and then as t grows the amount of water tends to be the modified parabolic profile 1− y4. In
formula (5.2.11), by modifying the n, we can get different step intensities for the moment that simulates the
’opening’ of the said water tube.
(a) t2
t2+1(b) t6
t6+1
Figure 31: Different ’water pipe’ slope functions
To solve this simulation we used 66049 degrees of freedom for the velocity and a ∆t = 0.0031. The simulation
was slow, but in the end we were able to extract and plot the expected solutions.
59
(a) t = 1.5
(b) t = 4.0
Figure 32: Different time solution plots
There are, however, better/more efficient methods to deal with turbulent flows on non-stationary Navier
Stokes (see [21]). At this point, we can try different plots for the solutions we have obtained - the last plot
will contain graphics for our solution in the form of ||u||2, where ||.|| is the Euclidean norm. We have chosen
to square it, so that differences in the values of our solution would be more explicit.
60
(a) t = 1.00 (b) t = 1.50
(c) t = 1.75 (d) t = 4.00
Figure 33: ||u||2 on different time frames
6 Conclusion and Future Work
We have shown that the proposed numerical implementation can handle with different benchmark idealized
problems in the two dimensional frame. It is, of course, far from being completed, if one aims to use it
for realistic 3D simulations. Several improvements can be introduced such as adaptive meshes for both the
domain level and the time discretization. Furthermore, parallelization techniques should be considered in
order to manage the tremendous computational effort required for realist simulations. One interesting work
to follow, would be to port those parallelized algorithms to a platform such as CUDA (from NVidia), and
take advantage of the huge parallelization power offered by the GPU.
61
7 Appendix
7.1 Stokes P2-P1 code
The code that we have implemented was based on the IFISS SolvingStokesq2q1.zip link made available in [8],
for the Stokes problem with Dirichlet boundary conditions. Here we give some guidelines about how to use
the equivalent implementation for the case of P2-P1 finite elements including Neumann boundary conditions.
To use this code, you can customize it for several different situations. Next we will describe what parameters
you are supposed to modify as a user and what files you are supposed to run. However we also provide
the preconfigured package for each and every example so that you only need to run main channel.m to
get the same results we did. In the following files you should run main channel.m to execute the main
Stokes equation solving algorithm. To configure the boundary conditions you should set outbc query with
a vector that describes the kind of boundary condition on each boundary layer. Several types are already
made available for you such as FLOW CHANNEL NATURAL or FLOW FULL DIR. If you want to
provide custom implementation for your own Dirichlet boundary conditions or your own Neumann boundary
conditions you should check dir bd flow.m (this function is called upon calculating Dirichlet boundary
conditions) or NEUMANN F1.m, NEUMANN F2.m, NEUMANN F3.m or NEUMANN F4.m to
implement Neumann boundary functions on boundary layer number 1,2,3 or 4 respectively. Furthermore
elements flag will alternate between Q2-Q1 and P2-P1 elements. nodesp variable measures the subdivision
level of our grid.
7.1.1 P2-P1 Poiseuille Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvva0pFTkFiTllRaEU
7.1.2 P2-P1 Cavity Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvOFhaZkF1NkZaMmc
7.1.3 P2-P1 Custom Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvRjI3WHBGcFZ2OHM
7.1.4 P2-P1 Custom Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvSm5DS0ZMRl84MVU
7.2 Newton method P2-P1 code
The following is the code we developed to implement our fully implicit Newton method simulations. The
examples explored in section 3.3 will be provided here as well. To run any of the following examples one
62
simply needs to execute main newton stabilized.m. To change viscosity of the fluid simply change v
variable. To change the boundary conditions applied on the fluid, simply call channel dom p2p1() with
the second argument set according to our previous definitions. To set the stopping criteria of our Newton
method, configure the stopping error with e tol. Variable it max controls the maximum number of cycles
we allow our Newton iteration to have. Naturally, as before, subdiv defines your grid subdivision level.
delta is the actual δ variable used in SUPG method. If delta is set to zero, the original Newnton method is
used, otherwise SUPG is used with δ set according to delta. The following are packages that reproduce any
of the examples explored in section 3.3
7.2.1 P2-P1 Newton Poiseuille Navier-Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvMnJfT21ReEI5QjA
7.2.2 P2-P1 Newton Cavity Navier-Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvOU1HQ3BsZE1kZWc
7.2.3 P2-P1 Newton Obstacle Navier-Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvdUNRQThDMXBNSWM
7.3 Newton SUPG method P2-P1 code
For SUPG we use the exact same code as for fully implicit Newton. The only difference is that, now, we
are not setting anywhere delta to be zero, instead, we set this variable to compute the heuristic (4.4.9).
Naturally, main newton stabilized.m needs to be run in order to achieve the same results as we did. The
following packages reproduce the same results we obtained in section 3.5.
7.3.1 P2-P1 Newton SUPG Poiseuille Navier-Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvaE53Skt5RWpwam8
7.3.2 P2-P1 Newton SUPG Cavity Navier-Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvbWVES042blY4eE0
7.3.3 P2-P1 Newton SUPG Obstacle Navier-Stokes example
https://drive.google.com/open?id=0B7jM6LDiQSvvcFBPeENFdEdsSTA
63
7.4 Non-Stationary Newton method P2-P1 code
The following code is the code used to implement our method for solving non-stationary Navier-Stokes equa-
tions. Naturally, all our examples from section 4.2 use this code. Therefore, in order to get the same results
as we did in the aforementioned section, one should execute main time.m. In order to alter the boundary
specific conditions, one needs to alter variable chn type (also altering the appropriate boundary Neuman-
n/Dirichlet functions). time init and time max mark the beginning and ending of our time dependent
simulation, respectively. Variable k defines the number of time intervals in wich we do our simulation in.
This value is calculated with 20k and stored in variable frames. delta t variable defines our ∆t value for
our simulation. v, as usual, is our variable for defining viscosity of our fluid. The following packages produce
our examples from section 4.2.
7.4.1 P2-P1 Non-Stationary Navier-Stokes Test Example
https://drive.google.com/open?id=0B7jM6LDiQSvvaFI3WUdIQ2dXd2s
7.4.2 P2-P1 Non-Stationary Navier-Stokes Cavity Example
https://drive.google.com/open?id=0B7jM6LDiQSvvZGxhcTdMQVBwZEk
7.4.3 P2-P1 Non-Stationary Navier-Stokes Womersly Example
https://drive.google.com/open?id=0B7jM6LDiQSvvWWZ2QVZSeEI4ZFU
7.4.4 P2-P1 Non-Stationary Navier-Stokes Obstacle Example
https://drive.google.com/open?id=0B7jM6LDiQSvvWVlFNm1ITkVEdzQ
64
References
[1] D. Braess. Finite elements: Theory, fast solvers, and applications in solid mechanics. Cambridge
University Press, Cambridge, United Kingdom, 2007.
[2] F. Brezzi and M. Fortin. Mixed and hybrid finite element methods, volume 15. Springer Science &
Business Media, Berlin, Germany, 2012.
[3] A. Brooks and T. Hughes. Streamline upwind/Petrov-Galerkin formulations for convection dominated
flows with particular emphasis on the incompressible Navier-Stokes equations. Computer methods in
applied mechanics and engineering, 32(1):199–259, 1982.
[4] A. de Moura. Analise numerica e simulacao das equacoes de Navier-Stokes 2d para fluidos incom-
pressıveis. Trabalho Final de Curso, Instituto Superior Tecnico, pages 51–100, 2001.
[5] C. Doering and J. Gibbon. Applied analysis of the Navier-Stokes equations, volume 12. Cambridge
University Press, Cambridge, United Kingdom, 1995.
[6] H. Elman. Preconditioning for the steady-state Navier-Stokes equations with low viscosity. SIAM Journal
on Scientific Computing, 20(4):1299–1316, 1999.
[7] H. Elman, A. Ramage, and D. Silvester. Algorithm 866: Ifiss, a matlab toolbox for modelling incom-
pressible flow. ACM Transactions on Mathematical Software (TOMS), 33(2):14, 2007.
[8] H. Elman, D. Silvester, and A. Wathen. Finite elements and fast iterative solvers. JOURNAL OF
FLUID MECHANICS, 557(1):474–475, 2006.
[9] H. C. Elman, A. Ramage, and D. J. Silvester. Incompressible flow iterative solution software (ifiss)
installation & software guide version 3.2, 2012, http://www.manchester.ac.uk/ifiss/ifiss guide 3.2.pdf.
[10] A. Ern and J. Guermond. Theory and practice of finite elements, volume 159. Springer Science &
Business Media, Berlin, Germany, 2013.
[11] L. Formaggia, F. Saleri, and A. Veneziani. Solving numerical PDEs: Problems, applications, exercises.
Springer Science & Business Media, Berlin, Germany, 2012.
[12] G. Galdi. An introduction to the mathematical theory of the Navier-Stokes equations: Steady-state
problems. Springer Science & Business Media, Berlin, Germany, 2011.
[13] A. Gauthier, F. Saleri, and A. Veneziani. A fast preconditioner for the incompressible Navier-Stokes
equations. Computing and Visualization in science, 6(2-3):105–112, 2004.
[14] U. Ghia, K. Ghia, and C. Shin. High-re solutions for incompressible flow using the Navier-Stokes
equations and a multigrid method. Journal of computational physics, 48(3):387–411, 1982.
65
[15] V. Girault and P. Raviart. Finite element methods for Navier-Stokes equations: theory and algorithms,
volume 5. Springer Science & Business Media, Berlin, Germany, 2012.
[16] P. Gresho and R. Sani. Incompressible Flow and the Finite Element Method, Volume 2, Isothermal
Laminar Flow. John Wiley & Sons, Chichester, 2000.
[17] L. John. Stabilized finite element methods for Dirichlet boundary control problems in fluid mechanics.
na, 2011.
[18] C. Johnson and J. Saranen. Streamline diffusion methods for the incompressible Euler and Navier-Stokes
equations. Mathematics of Computation, 47(175):1–18, 1986.
[19] D. Kay, D. Loghin, and A. Wathen. A preconditioner for the steady-state Navier-Stokes equations.
SIAM Journal on Scientific Computing, 24(1):237–256, 2002.
[20] O. Ladyzhenskaya and R. Silverman. The mathematical theory of viscous incompressible flow, volume 76.
Gordon and Breach, Oxfordshire, United Kingdom, 1969.
[21] R. Taylor O. Zienkiewicz and P. Nithiarasu. The Finite Element Method for Fluid Dynamics.
Butterworth-Heinemann, Oxford, United Kingdom, 2005.
[22] M. Olshanskii. An iterative solver for the Oseen problem and numerical solution of incompressible
Navier-Stokes equations. Numerical linear algebra with applications, 6(5):353–378, 1999.
[23] A. Quarteroni. Numerical models for differential problems, volume 2. Springer Science & Business Media,
Berlin, Germany, 2010.
[24] A. Quarteroni and A. Valli. Numerical approximation of partial differential equations, volume 23.
Springer Science & Business Media, Berlin, Germany, 2008.
[25] H. Rathod, K. Nagaraja, B. Venkatesudu, and N. Ramesh. Gauss Legendre quadrature over a triangle.
Journal of the Indian Institute of Science, 84(5):183, 2013.
[26] H. Roos, M. Stynes, L. Tobiska, and R. Kellogg. Numerical methods for singularly perturbed differential
equations. SIAM Review, 39(3):535, 1997.
[27] H. Roos, M. Stynes, L. Tobiska, and R. Kellogg. Numerical methods for singularly perturbed differential
equations. SIAM Review, 39(3):535, 1997.
[28] Y. Saad and M. Schultz. Gmres: A generalized minimal residual algorithm for solving nonsymmetric
linear systems. SIAM Journal on scientific and statistical computing, 7(3):856–869, 1986.
[29] D. Silvester, H. Elman, D. Kay, and A. Wathen. Efficient preconditioning of the linearized Navier-Stokes
equations for incompressible flow. Journal of Computational and Applied Mathematics, 128(1):261–279,
2001.
66
[30] O. Steinbach. Numerical approximation methods for elliptic boundary value problems: finite and boundary
elements. SSpringer Science & Business Media, Berlin, Germany, 2007.
[31] R. Stenberg. Analysis of mixed finite elements methods for the Stokes problem: a unified approach.
Mathematics of computation, 42(165):9–23, 1984.
[32] R. Temam. Navier-Stokes equations and nonlinear functional analysis, volume 66. Siam, Pennsylvania,
United States, 1995.
[33] R. Temam. Navier-Stokes equations: theory and numerical analysis, volume 343. American Mathematical
Soc., Rhode Island, United States, 2001.
[34] L. Tobiska and G. Lube. A modified streamline diffusion method for solving the stationary Navier-Stokes
equation. Numerische Mathematik, 59(1):13–29, 1991.
67