Manual
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 96
| Download | |
| Open PDF In Browser | View PDF |
The Finite Element Method in Geodynamics
C. Thieulot
January 6, 2019
Contents
1 Introduction
1.1 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Essential literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
4
4
4
2 The
2.1
2.2
2.3
2.4
2.5
physical equations of Fluid Dynamics
The heat transport equation - energy conservation equation
The momentum conservation equations . . . . . . . . . . .
The mass conservation equations . . . . . . . . . . . . . . .
The equations in ASPECT manual . . . . . . . . . . . . . .
the Boussinesq approximation: an Incompressible flow . . .
3 The building blocks of the Finite Element Method
3.1 Numerical integration . . . . . . . . . . . . . . . . .
3.1.1 in 1D - theory . . . . . . . . . . . . . . . . .
3.1.2 in 1D - examples . . . . . . . . . . . . . . . .
3.1.3 in 2D/3D - theory . . . . . . . . . . . . . . .
3.2 The mesh . . . . . . . . . . . . . . . . . . . . . . . .
3.3 A bit of FE terminology . . . . . . . . . . . . . . . .
3.4 Elements and basis functions in 1D . . . . . . . . . .
3.4.1 Linear basis functions . . . . . . . . . . . . .
3.4.2 Quadratic basis functions . . . . . . . . . . .
3.4.3 Cubic basis functions (Q3 ) . . . . . . . . . .
3.5 Elements and basis functions in 2D . . . . . . . . . .
3.5.1 The Q1 basis in 2D . . . . . . . . . . . . . . .
3.5.2 The Q2 basis in 2D . . . . . . . . . . . . . . .
3.6 The penalty approach . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
6
6
6
8
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
9
9
9
11
11
12
12
12
12
13
14
16
17
19
19
4 Additional techniques
4.1 The method of manufactured solutions . . . . . . . . .
4.2 Assigning values to quadrature points . . . . . . . . .
4.3 Matrix (Sparse) storage . . . . . . . . . . . . . . . . .
4.3.1 2D domain - One degree of freedom per node .
4.3.2 2D domain - Two degrees of freedom per node
4.3.3 in fieldstone . . . . . . . . . . . . . . . . . . . .
4.4 Mesh generation . . . . . . . . . . . . . . . . . . . . .
4.5 The value of the timestep . . . . . . . . . . . . . . . .
4.6 Tracking materials . . . . . . . . . . . . . . . . . . . .
4.7 Visco-Plasticity . . . . . . . . . . . . . . . . . . . . . .
4.8 Picard and Newton . . . . . . . . . . . . . . . . . . . .
4.9 The choice of solvers . . . . . . . . . . . . . . . . . . .
4.10 The SUPG formulation for the energy equation . . . .
4.11 Tracking materials and/or interfaces . . . . . . . . . .
4.12 Dealing with a free surface . . . . . . . . . . . . . . . .
4.13 Pressure smoothing . . . . . . . . . . . . . . . . . . . .
4.14 Pressure normalisation . . . . . . . . . . . . . . . . . .
4.15 Static condensation . . . . . . . . . . . . . . . . . . . .
4.16 Exporting data to vtk format . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
22
22
22
24
24
25
26
27
29
29
29
29
29
29
29
29
30
31
32
33
1
5 fieldstone: simple analytical solution
36
6 fieldstone: Stokes sphere
39
7 fieldstone: Convection in a 2D box
40
8 fieldstone: The lid driven cavity
8.1 the lid driven cavity problem (ldc=0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2 the lid driven cavity problem - regularisation I (ldc=1) . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3 the lid driven cavity problem - regularisation II (ldc=2) . . . . . . . . . . . . . . . . . . . . . . . . . .
42
42
42
42
9 fieldstone: solcx benchmark
44
10 fieldstone: solkz benchmark
46
11 fieldstone: solvi benchmark
47
12 fieldstone: the indentor benchmark
49
13 fieldstone: the annulus benchmark
51
14 fieldstone: stokes sphere (3D) - penalty
53
15 fieldstone: stokes sphere (3D) - mixed formulation
54
16 fieldstone: consistent pressure recovery
55
17 fieldstone: the Particle in Cell technique (1) - the effect of averaging
56
18 fieldstone: solving the full saddle point problem
59
19 fieldstone: solving the full saddle point problem in 3D
19.0.1 Constant viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19.0.2 Variable viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
62
63
20 fieldstone: solving the full saddle point problem with Q2 × Q1 elements
65
21 fieldstone: solving the full saddle point problem with Q3 × Q2 elements
68
22 fieldstone: the Busse benchmark
70
23 fieldstone: The non-conforming Q1 × P0 element
72
24 fieldstone: The stabilised Q1 × Q1 element
73
25 fieldstone: compressible flow (1)
74
26 fieldstone: compressible flow (2)
26.1 The physics . . . . . . . . . . . . . . . . . .
26.2 The numerics . . . . . . . . . . . . . . . . .
26.3 The experimental setup . . . . . . . . . . .
26.4 Scaling . . . . . . . . . . . . . . . . . . . . .
26.5 Conservation of energy 1 . . . . . . . . . . .
26.5.1 under BA and EBA approximations
26.5.2 under no approximation at all . . . .
26.6 Conservation of energy 2 . . . . . . . . . . .
26.7 The problem of the onset of convection . . .
26.8 results - BA - Ra = 104 . . . . . . . . . . .
26.9 results - BA - Ra = 105 . . . . . . . . . . .
26.10results - BA - Ra = 106 . . . . . . . . . . .
26.11results - EBA - Ra = 104 . . . . . . . . . .
26.12results - EBA - Ra = 105 . . . . . . . . . .
26.13Onset of convection . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
76
76
76
77
78
78
78
79
79
80
82
84
85
86
88
89
27 fieldstone: Rayleigh-Taylor instability (1)
90
3
1
Introduction
WARNING: this is work in progress
practical hands-on approach
as little as possible jargon
no mathematical proof
no optimised codes (readability over efficiency). avoiding as much as possible to have to look elsewhere. very
sequential, so unavoidable repetitions (jacobian, shape functions)
FE is one of several methods.
All the python scripts and this document are freely available at https://github.com/cedrict/fieldstone.
1.1
Acknowledgments
Jean Braun, Philippe Fullsack, Arie van den Berg. Lukas van de Wiel. Robert Myhill. Menno, Anne Too many BSc
and MSc students to name indivisually, although Job Mos did produce the very first version of fieldstone as part of
his MSc thesis. The ASPECT team in general and Wolfgang Bangerth and Timo Heister in particular.
1.2
Essential literature
a)
1.3
b)
c)
d)
Installation
python3.6 -m pip install --user numpy scipy matplotlib
4
e)
2
The physical equations of Fluid Dynamics
Symbol
t
x, y, z
v
ρ
η
λ
T
∇
∇·
p
ε̇(v)
α
k
Cp
H
βT
τ
σ
2.1
meaning
Time
Cartesian coordinates
velocity vector
mass density
dynamic viscosity
penalty parameter
temperature
gradient operator
divergence operator
pressure
strain rate tensor
thermal expansion coefficient
thermal conductivity
Heat capacity
intrinsic specific heat production
isothermal compressibility
deviatoric stress tensor
full stress tensor
unit
s
m
m· s−1
kg/m3
Pa· s
Pa· s
K
m−1
m−1
Pa
s−1
K−1
W/(m · K)
J/K
W/kg
Pa−1
Pa
Pa
The heat transport equation - energy conservation equation
Let us start from the heat transport equation as shown in Schubert, Turcotte and Olson [63]:
ρCp
DT
Dp
− αT
= ∇ · k∇T + Φ + ρH
Dt
Dt
with D/Dt being the total derivatives so that
DT
∂T
=
+ v · ∇T
Dt
∂t
Dp
∂p
=
+ v · ∇p
Dt
∂t
Solving for temperature, this equation is often rewritten as follows:
ρCp
Dp
DT
− ∇ · k∇T = αT
+ Φ + ρH
Dt
Dt
A note on the shear heating term Φ: In many publications, Φ is given by Φ = τij ∂j ui = τ : ∇v.
Φ
=
τij ∂j ui
=
2η ε̇dij ∂j ui
1
2η ε̇dij ∂j ui + ε̇dji ∂i uj
2
1
2η ε̇dij ∂j ui + ε̇dij ∂i uj
2
1
2η ε̇dij (∂j ui + ∂i uj )
2
2η ε̇dij ε̇ij
=
=
=
=
2η ε̇d : ε̇
1
= 2η ε̇d : ε̇d + (∇ · v)1
3
=
=
2η ε̇d : ε̇d + 2η ε̇d : 1(∇ · v)
=
2η ε̇d : ε̇d
(1)
Finally
Φ = τ : ∇v = 2η ε̇d : ε̇d = 2η (ε̇dxx )2 + (ε̇dyy )2 + 2(ε̇dxy )2
5
2.2
The momentum conservation equations
Because the Prandlt number is virtually zero in Earth science applications the Navier Stokes equations reduce to the
Stokes equation:
∇ · σ + ρg = 0
Since
σ = −p1 + τ
it also writes
−∇p + ∇ · τ + ρg = 0
Using the relationship τ = 2η ε̇d we arrive at
−∇p + ∇ · (2η ε̇d ) + ρg = 0
2.3
The mass conservation equations
The mass conservation equation is given by
Dρ
+ ρ∇ · v = 0
Dt
or,
∂ρ
+ ∇ · (ρv) = 0
∂t
In the case of an incompressible flow, then ∂ρ/∂t = 0 and ∇ρ = 0, i.e. Dρ/Dt = 0 and the remaining equation is
simply:
∇·v =0
2.4
The equations in ASPECT manual
The following is lifted off the ASPECT manual. We focus on the system of equations in a d = 2- or d = 3-dimensional
domain Ω that describes the motion of a highly viscous fluid driven by differences in the gravitational force due to
a density that depends on the temperature. In the following, we largely follow the exposition of this material in
Schubert, Turcotte and Olson [63].
Specifically, we consider the following set of equations for velocity u, pressure p and temperature T :
1
−∇ · 2η ε̇(v) − (∇ · v)1 + ∇p = ρg
in Ω,
(2)
3
∇ · (ρv) = 0
ρCp
∂T
+ v · ∇T
∂t
in Ω,
(3)
− ∇ · k∇T = ρH
1
+ 2η ε̇(v) − (∇ · v)1
3
1
: ε̇(v) − (∇ · v)1
3
(4)
+ αT (v · ∇p)
in Ω,
where ε̇(u) = 21 (∇u + ∇uT ) is the symmetric gradient of the velocity (often called the strain rate).
In this set of equations, (96) and (97) represent the compressible Stokes equations in which v = v(x, t) is the
velocity field and p = p(x, t) the pressure field. Both fields depend on space x and time t. Fluid flow is driven by the
gravity force that acts on the fluid and that is proportional to both the density of the fluid and the strength of the
gravitational pull.
Coupled to this Stokes system is equation (98) for the temperature field T = T (x, t) that contains heat conduction
terms as well as advection with the flow velocity v. The right hand side terms of this equation correspond to
• internal heat production for example due to radioactive decay;
• friction (shear) heating;
• adiabatic compression of material;
In order to arrive at the set of equations that ASPECT solves, we need to
6
• neglect the ∂p/∂t. WHY?
• neglect the ∂ρ/∂t . WHY?
from equations above.
—————————————Also, their definition of the shear heating term Φ is:
Φ = kB (∇ · v)2 + 2η ε̇d : ε̇d
For many fluids the bulk viscosity kB is very small and is often taken to be zero, an assumption known as the Stokes
assumption: kB = λ + 2η/3 = 0. Note that η is the dynamic viscosity and λ the second viscosity. Also,
τ = 2η ε̇ + λ(∇ · v)1
but since kB = λ + 2η/3 = 0, then λ = −2η/3 so
2
τ = 2η ε̇ − η(∇ · v)1 = 2η ε̇d
3
7
2.5
the Boussinesq approximation: an Incompressible flow
[from aspect manual] The Boussinesq approximation assumes that the density can be considered constant in all
occurrences in the equations with the exception of the buoyancy term on the right hand side of (96). The primary
result of this assumption is that the continuity equation (97) will now read
∇·v =0
This implies that the strain rate tensor is deviatoric. Under the Boussinesq approximation, the equations are much
simplified:
−∇ · [2η ε̇(v)] + ∇p = ρg
ρ0 C p
∂T
+ v · ∇T
∂t
in Ω,
(5)
∇ · (ρv) = 0
in Ω,
(6)
− ∇ · k∇T = ρH
in Ω
(7)
Note that all terms on the rhs of the temperature equations have disappeared, with the exception of the source term.
8
3
The building blocks of the Finite Element Method
3.1
Numerical integration
As we will see later, using the Finite Element method to solve problems involves computing integrals which are more
often than not too complex to be computed analytically/exactly. We will then need to compute them numerically.
[wiki] In essence, the basic problem in numerical integration is to compute an approximate solution to a definite
integral
Z b
f (x)dx
a
to a given degree of accuracy. This problem has been widely studied [?] and we know that if f (x) is a smooth function,
and the domain of integration is bounded, there are many methods for approximating the integral to the desired
precision.
There are several reasons for carrying out numerical integration.
• The integrand f (x) may be known only at certain points, such as obtained by sampling. Some embedded systems
and other computer applications may need numerical integration for this reason.
• A formula for the integrand may be known, but it may be difficult or impossible to find an antiderivative that is
an elementary function. An example of such an integrand is f (x) = exp(−x2 ), the antiderivative of which (the
error function, times a constant) cannot be written in elementary form.
• It may be possible to find an antiderivative symbolically, but it may be easier to compute a numerical approximation than to compute the antiderivative. That may be the case if the antiderivative is given as an infinite
series or product, or if its evaluation requires a special function that is not available.
3.1.1
in 1D - theory
The simplest method of this type is to let the interpolating function be a constant function (a polynomial of degree
zero) that passes through the point ((a + b)/2, f ((a + b)/2)).
This is called the midpoint rule or rectangle rule.
Z
b
f (x)dx ' (b − a)f (
a
a+b
)
2
insert here figure
The interpolating function may be a straight line (an affine function, i.e. a polynomial of degree 1) passing through
the points (a, f (a)) and (b, f (b)).
This is called the trapezoidal rule.
Z b
f (a) + f (b)
f (x)dx ' (b − a)
2
a
insert here figure
For either one of these rules, we can make a more accurate approximation by breaking up the interval [a, b] into
some number n of subintervals, computing an approximation for each subinterval, then adding up all the results. This
is called a composite rule, extended rule, or iterated rule. For example, the composite trapezoidal rule can be stated
as
!
Z b
n−1
b − a f (a) X
b−a
f (b)
f (x)dx '
+
f (a + k
)+
n
2
n
2
a
k=1
where the subintervals have the form [kh, (k + 1)h], with h = (b − a)/n and k = 0, 1, 2, . . . , n − 1.
a)
b)
The interval [−2, 2] is broken into 16 sub-intervals. The blue lines correspond to the approximation of the red curve
by means of a) the midpoint rule, b) the trapezoidal rule.
9
There are several algorithms for numerical integration (also commonly called ’numerical quadrature’, or simply
’quadrature’) . Interpolation with polynomials evaluated at equally spaced points in [a, b] yields the NewtonCotes
formulas, of which the rectangle rule and the trapezoidal rule are examples. If we allow the intervals between interpolation points to vary, we find another group of quadrature formulas, such as the Gauss(ian) quadrature formulas.
A Gaussian quadrature rule is typically more accurate than a NewtonCotes rule, which requires the same number of
function evaluations, if the integrand is smooth (i.e., if it is sufficiently differentiable).
An n−point Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule constructed to
yield an exact result for polynomials of degree 2n − 1 or less by a suitable choice of the points xi and weights wi for
i = 1, . . . , n.
The domain of integration for such a rule is conventionally taken as [−1, 1], so the rule is stated as
Z
+1
f (x)dx =
−1
n
X
wiq f (xiq )
iq =1
In this formula the xiq coordinate is the i-th root of the Legendre polynomial Pn (x).
It is important to note that a Gaussian quadrature will only produce good results if the function f (x) is well
approximated by a polynomial function within the range [−1, 1]. As a consequence, the method is not, for example,
suitable for functions with singularities.
Gauss-Legendre points and their weights.
As shown in the above table, it can be shown that the weight values must fulfill the following condition:
X
wiq = 2
(8)
iq
and it is worth noting that all quadrature point coordinates are symmetrical around the origin.
Since most quadrature formula are only valid on a specific interval, we now must address the problem of their use
outside of such intervals. The solution turns out to be quite simple: one must carry out a change of variables from
the interval [a, b] to [−1, 1].
We then consider the reduced coordinate r ∈ [−1, 1] such that
2
(x − a) − 1
b−a
r=
This relationship can be reversed such that when r is known, its equivalent coordinate x ∈ [a, b] can be computed:
b−a
(1 + r) + a
2
x=
From this it follows that
dx =
and then
Z
b
f (x)dx =
a
b−a
2
Z
b−a
dr
2
+1
f (r)dr '
−1
n
b−a X
wi f (riq )
2 i =1 q
q
10
3.1.2
in 1D - examples
example 1 Since we know how to carry out any required change of variables, we choose for simplicity a = −1,
b = +1. Let us take for example f (x) = π. Then we can compute the integral of this function over the interval [a, b]
exactly:
Z +1
Z +1
I=
f (x)dx = π
dx = 2π
−1
−1
We can now use a Gauss-Legendre formula to compute this same integral:
Z
+1
f (x)dx =
Igq =
−1
nq
X
wiq f (xiq ) =
iq =1
nq
X
wiq π = π
iq =1
nq
X
wiq = 2π
iq =1
| {z }
=2
where we have used the property of the weight values of Eq.(8). Since the actual number of points was never specified,
this result is valid for all quadrature rules.
example 2
Let us now take f (x) = mx + p and repeat the same exercise:
+1
Z
Z
−1
−1
Z
+1
Igq =
f (x)dx =
−1
+1
f (x)dx =
I=
nq
X
wiq f (xiq ) =
1
(mx + p)dx = [ mx2 + px]+1
−1 = 2p
2
nq
X
wiq (mxiq + p) = m
wiq xiq +p
iq =1
iq =1
iq =1
nq
X
|
nq
X
wiq = 2p
iq =1
{z
=0
}
| {z }
=2
since the quadrature points are symmetric w.r.t. to zero on the x-axis. Once again the quadrature is able to compute
the exact value of this integral: this makes sense since an n-point rule exactly integrates a 2n − 1 order polynomial
such that a 1 point quadrature exactly integrates a first order polynomial like the one above.
example 3
Let us now take f (x) = x2 . We have
+1
Z
I=
Z
+1
f (x)dx =
−1
−1
and
Z
+1
Igq =
f (x)dx =
−1
nq
X
1
2
x2 dx = [ x3 ]+1
−1 =
3
3
wiq f (xiq ) =
nq
X
iq =1
iq =1
(1)
• nq = 1: xiq = 0, wiq = 2. Igq = 0
√
√
(1)
(2)
(1)
(2)
• nq = 2: xq = −1/ 3, xq = 1/ 3, wq = wq = 1. Igq =
2
3
• It also works ∀nq > 2 !
3.1.3
in 2D/3D - theory
Let us now turn to a two-dimensional integral of the form
Z +1 Z +1
I=
f (x, y)dxdy
−1
−1
The equivalent Gaussian quadrature writes:
Igq '
nq nq
X
X
f (xiq , yjq )wiq wjq
iq =1 jq
11
wiq x2iq
3.2
The mesh
3.3
A bit of FE terminology
We introduce here some terminology for efficient element descriptions [30]:
• For triangles/tetrahedra, the designation Pm × Pn means that each component of the velocity is approximated
by continuous piecewise complete Polynomials of degree m and pressure by continuous piecewise complete Polynomials of degree n. For example P2 × P1 means
u ∼ a1 + a2 x + a3 y + a4 xy + a5 x2 + a6 y 2
with similar approximations for v, and
p ∼ b1 + b2 x + b3 y
Both velocity and pressure are continuous across element boundaries, and each triangular element contains 6
velocity nodes and three pressure nodes.
• For the same families, Pm × P−n is as above, except that pressure is approximated via piecewise discontinuous
polynomials of degree n. For instance, P2 × P−1 is the same as P2 P1 except that pressure is now an independent
linear function in each element and therefore discontinuous at element boundaries.
• For quadrilaterals/hexahedra, the designation Qm × Qn means that each component of the velocity is approximated by a continuous piecewise polynomial of degree m in each direction on the quadrilateral and likewise for
pressure, except that the polynomial is of degree n. For instance, Q2 × Q1 means
u ∼ a1 + a2 x + a3 y + a4 xy + a5 x2 + a6 y 2 + a7 x2 y + a8 xy 2 + a9 x2 y 2
and
p ∼ b1 + b2 x + b3 y + b4 xy
• For these same families, Qm × Q−n is as above, except that the pressure approximation is not continuous at
element boundaries.
• Again for the same families, Qm × P−n indicates the same velocity approximation with a pressure approximation
that is a discontinuous complete piecewise polynomial of degree n (not of degree n in each direction !)
+
or Q+
• The designation Pm
m means that some sort of bubble function was added to the polynomial approximation
for the velocity. You may also find the term ’enriched element’ in the literature.
• Finally, for n = 0, we have piecewise-constant pressure, and we omit the minus sign for simplicity.
Another point which needs to be clarified is the use of so-called ’conforming elements’ (or ’non-conforming elements’). Following again [30], conforming velocity elements are those for which the basis functions for a subset of
H 1 for the continuous problem (the first derivatives and their squares are integrable in Ω). For instance, the rotated
Q1 × P0 element of Rannacher and Turek (see section 23) is such that the velocity is discontinous across element edges,
so that the derivative does not exist there. Another typical example of non-conforming element is the Crouzeix-Raviart
element [12].
3.4
3.4.1
Elements and basis functions in 1D
Linear basis functions
Let f (r) be a C 1 function on the interval [−1 : 1] with f (−1) = f1 and f (1) = f2 .
Let us assume that the function f (r) is to be approximated on [−1, 1] by the first order polynomial
f (r) = a + br
12
(9)
Then it must fulfill
f (r = −1)
= a − b = f1
f (r = +1)
= a + b = f2
This leads to
1
1
(f1 + f2 )
b = (−f1 + f2 )
2
2
and then replacing a, b in Eq. (9) by the above values on gets
1
1
f (r) =
(1 − r) f1 + (1 + r) f2
2
2
a=
or
f (r) =
2
X
Ni (r)f1
i=1
with
3.4.2
N1 (r)
=
N2 (r)
=
1
(1 − r)
2
1
(1 + r)
2
(10)
Quadratic basis functions
Let f (r) be a C 1 function on the interval [−1 : 1] with f (−1) = f1 , f (0) = f2 and f (1) = f3 .
Let us assume that the function f (r) is to be approximated on [−1, 1] by the second order polynomial
f (r) = a + br + cr2
(11)
Then it must fulfill
f (r = −1) = a − b + c = f1
f (r = 0) = a
= f2
f (r = +1)
= a + b + c = f3
This leads to
1
1
(−f 1 + f 3) c = (f1 + f3 − 2f2 )
2
2
and then replacing a, b, c in Eq. (11) by the above values on gets
1
1
2
f (r) =
r(r − 1) f1 + (1 − r )f2 + r(r + 1) f3
2
2
a = f2
b=
or,
f (r) =
3
X
Ni (r)fi
i=1
with
N1 (r)
=
N2 (r)
=
N3 (r)
=
13
1
r(r − 1)
2
(1 − r2 )
1
r(r + 1)
2
(12)
3.4.3
Cubic basis functions (Q3 )
The 1D basis polynomial is given by
f (r) = a + br + cr2 + dr3
with the nodes at position -1,-1/3, +1/3 and +1.
= a − b + c − d = f1
b
c
d
= f2
f (−1/3) = a − + −
3 9 27
b
c
d
f (+1/3) = a − + −
= f3
3 9 27
f (+1) = a + b + c + d = f4
f (−1)
Adding the first and fourth equation and the second and third, one arrives at
f1 + f4 = 2a + 2c
f2 + f3 = 2a +
2c
9
and finally:
1
(−f1 + 9f2 + 9f3 − f4 )
16
9
c=
(f1 − f2 − f3 + f4 )
16
Combining the original 4 equations in a different way yields
a=
2b 2d
+
= f3 − f2
3
27
2b + 2d = f4 − f1
so that
1
(f1 − 27f2 + 27f3 − f4 )
16
9
(−f1 + 3f2 − 3f3 + f4 )
d=
16
b=
Finally,
f (r)
= a + b + cr2 + dr3
1
=
(−1 + r + 9r2 − 9r3 )f1
16
1
+
(9 − 27r − 9r2 + 27r3 )f2
16
1
+
(9 + 27r − 9r2 − 27r3 )f3
16
1
+
(−1 − r + 9r2 + 9r3 )f4
16
4
X
=
Ni (r)fi
i=1
where
N1
=
N2
=
N3
=
N4
=
1
(−1 + r + 9r2 − 9r3 )
16
1
(9 − 27r − 9r2 + 27r3 )
16
1
(9 + 27r − 9r2 − 27r3 )
16
1
(−1 − r + 9r2 + 9r3 )
16
14
Verification:
• Let us assume f (r) = C, then
fˆ(r) =
X
Ni (r)fi =
X
Ni C = C
i
X
Ni = C
i
so that a constant function is exactly reproduced, as expected.
• Let us assume f (r) = r, then f1 = −1, f2 = −1/3, f3 = 1/3 and f4 = +1. We then have
X
fˆ(r) =
Ni (r)fi
1
1
= −N1 (r) − N2 (r) + N3 (r) + N4 (r)
3
3
= [−(−1 + r + 9r2 − 9r3 )
1
− (9 − 27r − 9r2 − 27r3 )
3
1
+ (9 + 27r − 9r2 + 27r3 )
3
+(−1 − r + 9r2 + 9r3 )]/16
=
[−r + 9r + 9r − r]/16 + ...0...
= r
(13)
The basis functions derivative are given by
∂N1
∂r
∂N2
∂r
∂N3
∂r
∂N4
∂r
=
=
=
=
1
(1 + 18r − 27r2 )
16
1
(−27 − 18r + 81r2 )
16
1
(+27 − 18r − 81r2 )
16
1
(−1 + 18r + 27r2 )
16
Verification:
• Let us assume f (r) = C, then
∂ fˆ
∂r
=
X ∂Ni
i
= C
fi
∂r
X ∂Ni
i
=
∂r
C
[(1 + 18r − 27r2 )
16
+(−27 − 18r + 81r2 )
+(+27 − 18r − 81r2 )
+(−1 + 18r + 27r2 )]
=
0
15
• Let us assume f (r) = r, then f1 = −1, f2 = −1/3, f3 = 1/3 and f4 = +1. We then have
∂ fˆ
∂r
X ∂Ni
=
i
∂r
fi
1
[−(1 + 18r − 27r2 )
16
1
− (−27 − 18r + 81r2 )
3
1
+ (+27 − 18r − 81r2 )
3
+(−1 + 18r + 27r2 )]
1
[−2 + 18 + 54r2 − 54r2 ]
=
16
= 1
=
3.5
Elements and basis functions in 2D
Let us for a moment consider a single quadrilateral element in the xy-plane, as shown on the following figure:
Let us assume that we know the values of a given field u at the vertices. For a given point M inside the element in
the plane, what is the value of the field u at this point? It makes sense to postulate that uM = u(xM , yM ) will be
given by
uM = φ(u1 , u2 , u3 , u4 , xM , yM )
where φ is a function to be determined. Although φ is not unique, we can decide to express the value uM as a weighed
sum of the values at the vertices ui . One option could be to assign all four vertices the same weight, say 1/4 so that
uM = (u1 + u2 + u3 + u4 )/4, i.e. uM is simply given by the arithmetic mean of the vertices values. This approach
suffers from a major drawback as it does not use the location of point M inside the element. For instance, when
(xM , yM ) → (x2 , y2 ) we expect uM → u2 .
In light of this, we could now assume that the weights would depend on the position of M in a continuous fashion:
u(xM , yM ) =
4
X
Ni (xM , yM ) ui
i=1
where the Ni are continous (”well behaved”) functions which have the property:
Ni (xj , yj ) = δij
or, in other words:
N3 (x1 , y1 )
=
0
(14)
N3 (x2 , y2 )
=
0
(15)
N3 (x3 , y3 )
=
1
(16)
N3 (x4 , y4 )
=
0
(17)
The functions Ni are commonly called basis functions.
Omitting the M subscripts for any point inside the element, the velocity components u and v are given by:
û(x, y)
=
4
X
Ni (x, y) ui
(18)
Ni (x, y) vi
(19)
i=1
v̂(x, y)
=
4
X
i=1
16
Rather interestingly, one can now easily compute velocity gradients (and therefore the strain rate tensor) since we
have assumed the basis functions to be ”well behaved” (in this case differentiable):
4
˙xx (x, y)
=
∂u X ∂Ni
=
ui
∂x
∂x
i=1
˙yy (x, y)
=
∂v X ∂Ni
=
vi
∂y
∂y
i=1
˙xy (x, y)
=
1 ∂u 1 ∂v
1 X ∂Ni
1 X ∂Ni
+
=
ui +
vi
2 ∂y
2 ∂x
2 i=1 ∂y
2 i=1 ∂x
(20)
4
(21)
4
4
(22)
How we actually obtain the exact form of the basis functions is explained in the coming section.
3.5.1
The Q1 basis in 2D
In this section, we place ourselves in the most favorables case, i.e. the element is a square defined by −1 < r < 1,
−1 < s < 1 in the Cartesian coordinates system (r, s):
add corner numbering
This element is commonly called the reference element. How we go from the (x, y) coordinate system to the (r, s)
once and vice versa will be dealt later on. For now, the basis functions in the above reference element and in the
reduced coordinates system (r, s) are given by:
N1 (r, s)
=
0.25(1 − r)(1 − s)
N2 (r, s)
=
0.25(1 + r)(1 − s)
N3 (r, s)
=
0.25(1 + r)(1 + s)
N4 (r, s)
=
0.25(1 − r)(1 + s)
The partial derivatives of these functions with respect to r ans s automatically follow:
∂N1
(r, s)
∂r
∂N2
(r, s)
∂r
∂N3
(r, s)
∂r
∂N4
(r, s)
∂r
∂N1
(r, s)
∂s
∂N2
(r, s)
∂s
∂N3
(r, s)
∂s
∂N4
(r, s)
∂s
=
−0.25(1 − s)
=
+0.25(1 − s)
=
+0.25(1 + s)
=
−0.25(1 + s)
=
−0.25(1 − r)
=
−0.25(1 + r)
=
+0.25(1 + r)
=
+0.25(1 − r)
17
Let us go back to Eq.(19). And let us assume that the function v(r, s) = C so that vi = C for i = 1, 2, 3, 4. It then
follows that
v̂(r, s) =
4
X
Ni (r, s) vi = C
i=1
4
X
Ni (r, s) = C[N1 (r, s) + N2 (r, s) + N3 (r, s) + N4 (r, s)] = C
i=1
This is a very important property: if the v function used to assign values at the vertices is constant, then the value of
v̂ anywhere in the element is exactly C. If we now turn to the derivatives of v with respect to r and s:
4
4
4
4
X ∂Ni
X ∂Ni
∂v̂
(r, s) =
(r, s) vi = C
(r, s) = C [−0.25(1 − s) + 0.25(1 − s) + 0.25(1 + s) − 0.25(1 + s)] = 0
∂r
∂r
∂r
i=1
i=1
X ∂Ni
X ∂Ni
∂v̂
(r, s) =
(r, s) vi = C
(r, s) = C [−0.25(1 − r) − 0.25(1 + r) + 0.25(1 + r) + 0.25(1 − r)] = 0
∂s
∂s
∂s
i=1
i=1
We reassuringly find that the derivative of a constant field anywhere in the element is exactly zero.
If we now choose v(r, s) = ar + bs with a and b two constant scalars, we find:
v̂(r, s)
=
4
X
Ni (r, s) vi
(23)
Ni (r, s)(ari + bsi )
(24)
i=1
=
4
X
i=1
= a
4
X
Ni (r, s)ri +b
Ni (r, s)si
(25)
i=1
i=1
|
4
X
{z
|
}
r
{z
s
}
= a [0.25(1 − r)(1 − s)(−1) + 0.25(1 + r)(1 − s)(+1) + 0.25(1 + r)(1 + s)(+1) + 0.25(1 − r)(1 + s)(−1)]
+ b [0.25(1 − r)(1 − s)(−1) + 0.25(1 + r)(1 − s)(−1) + 0.25(1 + r)(1 + s)(+1) + 0.25(1 − r)(1 + s)(+1)]
= a [−0.25(1 − r)(1 − s) + 0.25(1 + r)(1 − s) + 0.25(1 + r)(1 + s) − 0.25(1 − r)(1 + s)]
+ b [−0.25(1 − r)(1 − s) − 0.25(1 + r)(1 − s) + 0.25(1 + r)(1 + s) + 0.25(1 − r)(1 + s)]
= ar + bs
(26)
verify above eq. This set of bilinear shape functions is therefore capable of exactly representing a bilinear field. The
derivatives are:
∂v̂
(r, s)
∂r
=
4
X
∂Ni
i=1
=
a
∂r
(r, s) vi
4
X
∂Ni
i=1
∂r
(r, s)ri + b
(27)
4
X
∂Ni
i=1
∂r
(r, s)si
=
a [−0.25(1 − s)(−1) + 0.25(1 − s)(+1) + 0.25(1 + s)(+1) − 0.25(1 + s)(−1)]
+
b [−0.25(1 − s)(−1) + 0.25(1 − s)(−1) + 0.25(1 + s)(+1) − 0.25(1 + s)(+1)]
a
[(1 − s) + (1 − s) + (1 + s) + (1 + s)]
4
b
[(1 − s) − (1 − s) + (1 + s) − (1 + s)]
4
a
=
+
=
(28)
(29)
Here again, we find that the derivative of the bilinear field inside the element is exact: ∂∂rv̂ = ∂v
∂r .
However, following the same methodology as above, one can easily prove that this is no more true for polynomials
of degree strivtly higher than 1. This fact has serious consequences: if the solution to the problem at hand is for
instance a parabola, the Q1 shape functions cannot represent the solution properly, but only by approximating the
parabola in each element by a line. As we will see later, Q2 basis functions can remedy this problem by containing
themselves quadratic terms.
18
3.5.2
3.6
The Q2 basis in 2D
The penalty approach
In order to impose the incompressibility constraint, two widely used procedures are available, namely the Lagrange
multiplier method and the penalty method [2, 34]. The latter is implemented in elefant, which allows for the
elimination of the pressure variable from the momentum equation (resulting in a reduction of the matrix size).
Mathematical details on the origin and validity of the penalty approach applied to the Stokes problem can for
instance be found in [13], [56] or [32].
The penalty formulation of the mass conservation equation is based on a relaxation of the incompressibility constraint and writes
p
(30)
∇·v+ =0
λ
where λ is the penalty parameter, that can be interpreted (and has the same dimension) as a bulk viscosity. It is
equivalent to say that the material is weakly compressible. It can be shown that if one chooses λ to be a sufficiently
large number, the continuity equation ∇ · v = 0 will be approximately satisfied in the finite element solution. The
value of λ is often recommended to be 6 to 7 orders of magnitude larger than the shear viscosity [21, 35].
Equation (30) can be used to eliminate the pressure in Eq. (??) so that the mass and momentum conservation
equations fuse to become :
∇ · (2η ε̇(v)) + λ∇(∇ · v) = ρg = 0
(31)
[47] have established the equivalence for incompressible problems between the reduced integration of the penalty
term and a mixed Finite Element approach if the pressure nodes coincide with the integration points of the reduced
rule.
In the end, the elimination of the pressure unknown in the Stokes equations replaces the original saddle-point
Stokes problem [3] by an elliptical problem, which leads to a symmetric positive definite (SPD) FEM matrix. This is
the major benefit of the penalized approach over the full indefinite solver with the velocity-pressure variables. Indeed,
the SPD character of the matrix lends itself to efficient solving stragegies and is less memory-demanding since it is
sufficient to store only the upper half of the matrix including the diagonal [29] . ToDo: list codes which use this
approach.
19
The stress tensor σ is symmetric (i.e. σij = σji ). For simplicity I will now focus on a Stokes flow in two dimensions.
Since the penalty formulation is only valid for incompressible flows, then ˙ = ˙ d so that the d superscript is
ommitted in what follows. The stress tensor can also be cast in vector format:
σxx
−p
˙xx
σyy = −p + 2η ˙yy
σxy
0
˙xy
˙xx + ˙yy
˙xx
= λ ˙xx + ˙yy + 2η ˙yy
0
˙xy
∂u
∂x
2 0 0
1 1 0
∂v
1 1 0 +η 0 2 0 ·
= λ
∂y
0 0 1
0 0 0
{z
}
{z
}
|
|
∂u
+ ∂v
K
Remember that
∂y
C
4
∂x
4
∂u X ∂Ni
=
ui
∂x
∂x
i=1
∂v X ∂Ni
=
vi
∂y
∂y
i=1
4
4
X
X
∂u ∂v
∂Ni
∂Ni
+
=
ui +
vi
∂y
∂x
∂y
∂x
i=1
i=1
so that
∂u
∂x
∂v
∂y
∂u
∂y
+
∂v
∂x
=
∂N1
∂x
0
∂N2
∂x
0
∂N3
∂x
0
∂N4
∂x
0
0
∂N1
∂y
0
∂N2
∂y
0
∂N3
∂y
0
∂N4
∂y
∂N1
∂y
∂N1
∂x
∂N2
∂y
∂N2
∂x
∂N3
∂y
∂N3
∂x
∂N3
∂y
∂N4
∂x
|
{z
B
·
}
|
u1
v1
u2
v2
u3
v3
u4
v4
{z
V
}
Finally,
σxx
~σ = σyy = (λK + ηC) · B · V
σxy
We will now establish the weak form of the momentum conservation equation. We start again from
∇·σ+b=0
For the Ni ’s ’regular enough’, we can write:
Z
Z
Ni ∇ · σdΩ +
Ωe
Ni bdΩ = 0
Ωe
We can integrate by parts and drop the surface term1 :
Z
Z
∇Ni · σdΩ =
Ωe
Ni bdΩ
Ωe
or,
Z
∂Ni
∂x
0
∂Ni
∂y
0
∂Ni
∂y
∂Ni
∂x
Ωe
1 We
Z
σxx
· σyy dΩ =
Ni bdΩ
Ωe
σxy
will come back to this at a later stage
20
Let i = 1, 2, 3, 4 and stack the resulting four equations on top of one another.
Z
∂N1
∂x
0
∂N1
∂y
∂N2
∂x
0
Ωe
Z
0
∂N2
∂y
∂N3
∂x
0
Ωe
Z
∂N1
∂x
∂N2
∂y
∂N2
∂x
∂N3
∂y
0
∂N3
∂y
∂N4
∂x
0
∂N3
∂x
∂N4
∂y
0
∂N4
∂y
∂N4
∂x
Ωe
Z
∂N1
∂y
0
Ωe
σxx
· σyy dΩ
σxy
σxx
· σyy dΩ
σxy
σxx
· σyy dΩ
σxy
σxx
· σyy dΩ
σxy
Z
=
N1
Ωe
Z
=
Ni
Ωe
Z
=
N3
Ωe
Z
=
N4
Ωe
bx
by
bx
by
bx
by
bx
by
dΩ
(32)
dΩ
(33)
dΩ
(34)
dΩ
(35)
We easily recognize B T inside the integrals! Let us define
NbT = (N1 bx , N1 by , ...N4 bx , N4 by )
then we can write
Z
σxx
B T · σyy dΩ =
Nb dΩ
Ωe
Ωe
σxy
Z
and finally:
Z
Z
T
B · [λK + ηC] · B · V dΩ =
Ωe
Nb dΩ
Ωe
Since V contains the velocities at the corners, it does not depend on the x or y coordinates so it can be taking outside
of the integral:
Z
Z
B T · [λK + ηC] · BdΩ · |{z}
V =
Nb dΩ
Ωe
Ωe
{z
} (8x1) | {z }
|
Ael (8×8)
Bel (8×1)
or,
Z
Z
Z
T
T
λB · K · BdΩ +
V =
Nb dΩ
ηB · C · BdΩ · |{z}
Ω
Ω
| Ωe
{z
} | e
{z
} (8x1) | e {z }
Aλ
el (8×8)
Aη
el (8×8)
INTEGRATION - MAPPING
1. partition domain Ω into elements Ωe , e = 1, ...nel .
2. loop over elements and for each element compute Ael , Bel
3. a node belongs to several elements
→ need to assemble Ael and Bel in A, B
4. apply boundary conditions
5. solve system: x = A−1 · B
6. visualise/analyse x
21
Bel (8×1)
4
Additional techniques
4.1
The method of manufactured solutions
4.2
Assigning values to quadrature points
As we have seen in Section ??, the building of the elemental matrix and rhs requires (at least) to assign a density and
viscosity value to each quadrature point inside the element. Depending on the type of modelling, this task can prove
more complex than one might expect and have large consequences on the solution accuracy.
Here are several options:
• The simplest way (which is often used for benchmarks) consists in computing the ’real’ coordinates (xq , yq , zq ) of
a given quadrature point based on its reduced coordinates (rq , sq , tq ), and passing these coordinates to a function
which returns density and/or viscosity at this location. For instance, for the Stokes sphere:
def rho(x,y):
if (x-.5)**2+(y-0.5)**2<0.123**2:
val=2.
else:
val=1.
return val
def mu(x,y):
if (x-.5)**2+(y-0.5)**2<0.123**2:
val=1.e2
else:
val=1.
return val
This is very simple, but it has been shown to potentially be problematic. In essence, it can introduce very large
contrasts inside a single element and perturb the quadrature. Please read section 3.3 of [33] and/or have a look
at the section titled ”Averaging material properties” in the ASPECT manual.
• another similar approach consists in assigning a density and viscosity value to the nodes of the FE mesh first,
and then using these nodal values to assign values to the quadrature points. Very often ,and quite logically, the
shape functions are used to this effect. Indeed we have seen before that for any point (r, s, t) inside an element
we have
m
X
fh (r, s, t) =
fi Ni (r, s, t)
i
where the fi are the nodal values and the Ni the corresponding basis functions.
In the case of linear elements (Q1 basis functions), this is straightforward. In fact, the basis functions Ni can be
seen as moving weights: the closer the point is to a node, the higher the weight (basis function value).
However, this is quite another story for quadratic elements (Q2 basis functions). In order to illustrate the
problem, let us consider a 1D problem. The basis functions are
N1 (r) =
1
r(r − 1)
2
N2 (r) = 1 − r2
N3 (r) =
1
r(r + 1)
2
Let us further assign: ρ1 = ρ2 = 0 and ρ3 = 1. Then
ρh (r) =
m
X
ρi Ni (r) = N3 (r)
i
There lies the core of the problem: the N3 (r) basis function is negative for r ∈ [−1, 0]. This means that
the quadrature point in this interval will be assigned a negative density, which is nonsensical and numerically
problematic!
use 2X Q1. write about it !
The above methods work fine as long as the domain contains a single material. As soon as there are multiple fluids
in the domain a special technique is needed to track either the fluids themselves or their interfaces. Let us start with
markers. We are then confronted to the infernal trio (a menage a trois?) which is present for each element, composed
of its nodes, its markers and its quadrature points.
22
Each marker carries the material information (density and viscosity). This information must ultimately be projected
onto the quadrature points. Two main options are possible: an algorithm is designed and projects the marker-based
fields onto the quadrature points directly or the marker fields are first projected onto the FE nodes and then onto the
quadrature points using the techniques above.
23
4.3
Matrix (Sparse) storage
The FE matrix is the result of the assembly process of all elemental matrices. Its size can become quite large when
the resolution is being increased (from thousands of lines/columns to tens of millions).
One important property of the matrix is its sparsity. Typically less than 1% of the matrix terms is not zero and
this means that the matrix storage can and should be optimised. Clever storage formats were designed early on since
the amount of RAM memory in computers was the limiting factor 3 or 4 decades ago. [58]
There are several standard formats:
• compressed sparse row (CSR) format
• compressed sparse column format (CSC)
• the Coordinate Format (COO)
• Skyline Storage Format
• ...
I focus on the CSR format in what follows.
4.3.1
2D domain - One degree of freedom per node
Let us consider again the 3 × 2 element grid which counts 12 nodes.
8=======9======10======11
|
|
|
|
| (3) | (4) | (5) |
|
|
|
|
4=======5=======6=======7
|
|
|
|
| (0) | (1) | (2) |
|
|
|
|
0=======1=======2=======3
In the case there is only a single
X
X
X
X
degree of freedom per node, the assembled FEM matrix will look like this:
X
X X
X X
X X X
X X X
X X X
X X
X X
X
X X
X X
X X
X X X
X X X
X X X
X X X
X X X
X X
X X
X X
X X
X X
X X X
X X X
X X X
X X X
X X
X X
where the X stand for non-zero terms. This matrix structure stems from the fact that
• node 0 sees nodes 0,1,4,5
• node 1 sees nodes 0,1,2,4,5,6
• node 2 sees nodes 1,2,3,5,6,7
• ...
• node 5 sees nodes 0,1,2,4,5,6,8,9,10
• ...
• node 10 sees nodes 5,6,7,9,10,11
• node 11 sees nodes 6,7,10,11
24
In light thereof, we have
• 4 corner nodes which have 4 neighbours (counting themselves)
• 2(nnx-2) nodes which have 6 neighbours
• 2(nny-2) nodes which have 6 neighbours
• (nnx-2)×(nny-2) nodes which have 9 neighbours
In total, the number of non-zero terms in the matrix is then:
N Z = 4 × 4 + 4 × 6 + 2 × 6 + 2 × 9 = 70
In general, we would then have:
N Z = 4 × 4 + [2(nnx − 2) + 2(nny − 2)] × 6 + (nnx − 2)(nny − 2) × 9
Let us temporarily assume nnx = nny = n. Then the matrix size (total number of unknowns) is N = n2 and
N Z = 16 + 24(n − 2) + 9(n − 2)2
A full matrix array would contain N 2 = n4 terms. The ratio of N Z (the actual number of reals to store) to the full
matrix size (the number of reals a full matrix contains) is then
R=
16 + 24(n − 2) + 9(n − 2)2
n4
It is then obvious that when n is large enough R ∼ 1/n2 .
CSR stores the nonzeros of the matrix row by row, in a single indexed array A of double precision numbers.
Another array COLIND contains the column index of each corresponding entry in the A array. A third integer array
RWPTR contains pointers to the beginning of each row, which an additional pointer to the first index following the
nonzeros of the matrix A. A and COLIND have length NZ and RWPTR has length N+1.
In the case of the here-above matrix, the arrays COLIND and RWPTR will look like:
COLIN D = (0, 1, 4, 5, 0, 1, 2, 4, 5, 6, 1, 2, 3, 5, 6, 7, ..., 6, 7, 10, 11)
RW P T R = (0, 4, 10, 16, ...)
4.3.2
2D domain - Two degrees of freedom per node
When there are now two degrees of freedom per node, such as in the case of the Stokes equation in two-dimensions,
the size of the K matrix is given by
NfemV=nnp∗ ndofV
In the case of the small grid above, we have NfemV=24. Elemental matrices are now 8 × 8 in size.
We still have
• 4 corner nodes which have 4 ,neighbours,
• 2(nnx-2) nodes which have 6 neighbours
• 2(nny-2) nodes which have 6 neighbours
• (nnx-2)x(nny-2) nodes which have 9 neighbours,
but now each degree of freedom from a node sees the other two degrees of freedom of another node too. In that case,
25
the number of
X X
X X
X X
X X
X X
X X
X X
X X
nonzeros has been multiplied by four and the assembled FEM matrix looks like:
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
Note that the degrees of freedom are organised as follows:
(u0 , v0 , u1 , v1 , u2 , v2 , ...u11 , v11 )
In general, we would then have:
N Z = 4 [4 × 4 + [2(nnx − 2) + 2(nny − 2)] × 6 + (nnx − 2)(nny − 2) × 9]
and in the case of the small grid, the number of non-zero terms in the matrix is then:
N Z = 4 [4 × 4 + 4 × 6 + 2 × 6 + 2 × 9] = 280
In the case of the here-above matrix, the arrays COLIND and RWPTR will look like:
COLIN D = (0, 1, 2, 3, 8, 9, 10, 11, 0, 1, 2, 3, 8, 9, 10, 11, ...)
RW P T R = (0, 8, 16, 28, ...)
4.3.3
in fieldstone
The majority of the codes have the FE matrix being a full array
a mat = np . z e r o s ( ( Nfem , Nfem ) , dtype=np . f l o a t 6 4 )
and it is converted to CSR format on the fly in the solve phase:
s o l = s p s . l i n a l g . s p s o l v e ( s p s . c s r m a t r i x ( a mat ) , r h s )
Note that linked list storages can be used (lil matrix). Substantial memory savings but much longer compute
times.
26
4.4
Mesh generation
Before basis functions can be defined and PDEs can be discretised and solved we must first tesselate the domain with
polygons, e.g. triangles and quadrilaterals in 2D, tetrahedra, prisms and hexahedra in 3D.
When the domain is itself simple (e.g. a rectangle, a sphere, ...) the mesh (or grid) can be (more or less) easily
produced and the connectivity array filled with straightforward algorithms [69]. However, real life applications can
involve etremely complex geometries (e.g. a bridge, a human spine, a car chassis and body, etc ...) and dedicated
algorithms/softwares must be used (see [71, 25, 75]).
We usually distinguish between two broad classes of grids: structured grids (with a regular connectivity) and
unstructured grids (with an irregular connectivity).
On the following figure are shown various triangle- tetrahedron-based meshes. These are obviously better suited
for simulations of complex geometries:
add more examples
Let us now focus on the case of a rectangular computational domain of size Lx × Ly with a regular mesh composed
of nelx×nely=nel quadrilaterals. There are then nnx×nny=nnp grid points. The elements are of size hx×hy with
hx=Lx/nelx.
We have no reason to come up with an irregular/ilogical node numbering so we can number nodes row by row or
column by column as shown on the example hereunder of a 3×2 grid:
8=======9======10======11
|
|
|
|
| (3) | (4) | (5) |
|
|
|
|
4=======5=======6=======7
|
|
|
|
| (0) | (1) | (2) |
|
|
|
|
0=======1=======2=======3
2=======5=======8======11
|
|
|
|
| (1) | (3) | (5) |
|
|
|
|
1=======4=======7======10
|
|
|
|
| (0) | (2) | (4) |
|
|
|
|
0=======3=======6=======9
"row by row"
"column by column"
The numbering of the elements themselves could be done in a somewhat chaotic way but we follow the numbering
of the nodes for simplicity. The row by row option is the adopted one in fieldstoneand the coordinates of the points
are computed as follows:
x = np . empty ( nnp , dtype=np . f l o a t 6 4 )
y = np . empty ( nnp , dtype=np . f l o a t 6 4 )
counter = 0
f o r j i n r a n g e ( 0 , nny ) :
f o r i i n r a n g e ( 0 , nnx ) :
x [ c o u n t e r ]= i ∗hx
y [ c o u n t e r ]= j ∗hy
c o u n t e r += 1
27
The inner loop has i ranging from 0 to nnx-1 first for j=0, 1, ... up to nny-1 which indeed corresponds to the row
by row numbering.
We now turn to the connectivity. As mentioned before, this is a structured mesh so that the so-called connectivity
array, named icon in our case, can be filled easily. For each element we need to store the node identities of its vertices.
Since there are nel elements and m=4 corners, this is a m×nel array. The algorithm goes as follows:
i c o n =np . z e r o s ( (m, n e l ) , dtype=np . i n t 1 6 )
counter = 0
f o r j in range (0 , nely ) :
f o r i in range (0 , nelx ) :
i c o n [ 0 , c o u n t e r ] = i + j ∗ nnx
i c o n [ 1 , c o u n t e r ] = i + 1 + j ∗ nnx
i c o n [ 2 , c o u n t e r ] = i + 1 + ( j + 1 ) ∗ nnx
i c o n [ 3 , c o u n t e r ] = i + ( j + 1 ) ∗ nnx
c o u n t e r += 1
In the case of the 3×2 mesh, the icon is filled as follows:
element id→
node id↓
0
1
2
3
0
1
2
3
4
5
0
1
5
4
1
2
6
5
2
3
7
6
4
5
9
8
5
6
10
9
6
7
11
10
It is to be understood as follows: element #4 is composed of nodes 5, 6, 10 and 9. Note that nodes are always stored
in a counter clockwise manner, starting at the bottom left. This is very important since the corresponding basis
functions and their derivatives will be labelled accordingly.
In three dimensions things are very similar. The mesh now counts nelx×nely×nelz=nel elements which represent
a cuboid of size Lx×Ly×Lz. The position of the nodes is obtained as follows:
x = np . empty ( nnp , dtype=np . f l o a t 6 4 )
y = np . empty ( nnp , dtype=np . f l o a t 6 4 )
z = np . empty ( nnp , dtype=np . f l o a t 6 4 )
c o u n t e r=0
f o r i i n r a n g e ( 0 , nnx ) :
f o r j i n r a n g e ( 0 , nny ) :
f o r k i n r a n g e ( 0 , nnz ) :
x [ c o u n t e r ]= i ∗hx
y [ c o u n t e r ]= j ∗hy
z [ c o u n t e r ]=k∗ hz
c o u n t e r += 1
The connectivity array is now of size m×nel with m=8:
i c o n =np . z e r o s ( (m, n e l ) , dtype=np . i n t 1 6 )
counter = 0
f o r i in range (0 , nelx ) :
f o r j in range (0 , nely ) :
f o r k in range (0 , n e l z ) :
i c o n [ 0 , c o u n t e r ]=nny∗ nnz ∗ ( i
)+nnz ∗ ( j
)+k
i c o n [ 1 , c o u n t e r ]=nny∗ nnz ∗ ( i +1)+nnz ∗ ( j
)+k
i c o n [ 2 , c o u n t e r ]=nny∗ nnz ∗ ( i +1)+nnz ∗ ( j +1)+k
i c o n [ 3 , c o u n t e r ]=nny∗ nnz ∗ ( i
)+nnz ∗ ( j +1)+k
i c o n [ 4 , c o u n t e r ]=nny∗ nnz ∗ ( i
)+nnz ∗ ( j
)+k+1
i c o n [ 5 , c o u n t e r ]=nny∗ nnz ∗ ( i +1)+nnz ∗ ( j
)+k+1
i c o n [ 6 , c o u n t e r ]=nny∗ nnz ∗ ( i +1)+nnz ∗ ( j +1)+k+1
i c o n [ 7 , c o u n t e r ]=nny∗ nnz ∗ ( i
)+nnz ∗ ( j +1)+k+1
c o u n t e r += 1
produce drawing of node numbering
28
4.5
The value of the timestep
4.6
Tracking materials
4.7
Visco-Plasticity
4.8
Picard and Newton
4.9
The choice of solvers
4.10
The SUPG formulation for the energy equation
4.11
Tracking materials and/or interfaces
4.12
Dealing with a free surface
29
4.13
Pressure smoothing
It has been widely documented that the use of the Q1 × P0 element is not without problems. Aside from the
consequences it has on the FE matrix properties, we will here focus on another unavoidable side effect: the spurious
pressure checkerboard modes.
These modes have been thoroughly analysed [31, 9, 59, 60]. They can be filtered out [9] or simply smoothed [43].
On the following figure (a,b), pressure fields for the lid driven cavity experiment are presented for both an even
and un-even number of elements. We see that the amplitude of the modes can sometimes be so large that the ’real’
pressure is not visible and that something as simple as the number of elements in the domain can trigger those or not
at all.
b)
a)
c)
a) element pressure for a 32x32 grid and for a 33x33 grid;
b) image from [19, p307] for a manufactured solution;
c) elemental pressure and smoothed pressure for the punch experiment [70]
The easiest post-processing step that can be used (especially when a regular grid is used) is explained in [70]: ”The
element-to-node interpolation is performed by averaging the elemental values from elements common to each node;
the node-to-element interpolation is performed by averaging the nodal values element-by-element. This method is not
only very efficient but produces a smoothing of the pressure that is adapted to the local density of the octree. Note
that these two steps can be repeated until a satisfying level of smoothness (and diffusion) of the pres- sure field is
attained.”
In the codes which rely on the Q1 × P0 element, the (elemental) pressure is simply defined as
p=np . z e r o s ( n e l , dtype=np . f l o a t 6 4 )
while the nodal pressure is then defined as
q=np . z e r o s ( nnp , dtype=np . f l o a t 6 4 )
The element-to-node algorithm is then simply (in 2D):
count=np . z e r o s ( nnp , dtype=np . i n t 1 6 )
f o r i e l in range (0 , nel ) :
q [ i c o n [ 0 , i e l ]]+=p [ i e l ]
q [ i c o n [ 1 , i e l ]]+=p [ i e l ]
q [ i c o n [ 2 , i e l ]]+=p [ i e l ]
q [ i c o n [ 3 , i e l ]]+=p [ i e l ]
count [ i c o n [ 0 , i e l ]]+=1
count [ i c o n [ 1 , i e l ]]+=1
count [ i c o n [ 2 , i e l ]]+=1
count [ i c o n [ 3 , i e l ]]+=1
q=q/ count
produce figure to explain this
link to proto paper
link to least square and nodal derivatives
30
4.14
Pressure normalisation
31
4.15
Static condensation
32
4.16
Exporting data to vtk format
This format seems to be the universally accepted format for 2D and 3D visualisation in Computational Geodynamics.
Such files can be opened with free softwares such as Paraview 2 , MayaVi 3 or Visit 4 .
Unfortunately it is my experience that no simple tutorial exists about how to build such files. There is an official
document which describes the vtk format5 but it delivers the information in a convoluted way. I therefore describe
hereafter how fieldstone builds the vtk files.
I hereunder show vtk file corresponding to the 3x2 grid presented earlier 4.4. In this particular example there are:
• 12 nodes and 6 elements
• 1 elemental field: the pressure p)
• 2 nodal fields: 1 scalar (the smoothed pressure q), 1 vector (the velocity field u,v,0)
Note that vtk files are inherently 3D so that even in the case of a 2D simulation the z-coordinate of the points and
for instance their z-velocity component must be provided. The file, usually called solution.vtu starts with a header:
We then proceed to write the node coordinates as follows:
NumberOfComponents= ’ 3 ’ Format= ’ a s c i i ’>
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00
These are followed by the elemental field(s):
−1.333333 e+00
−3.104414 e −10
1 . 3 3 3 3 3 3 e+00
−1.333333 e+00
8 . 2 7 8 4 1 7 e −17
1 . 3 3 3 3 3 3 e+00
Nodal quantities are written next:
0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
8 . 8 8 8 8 8 5 e −08 −8.278405 e −24 0 . 0 0 0 0 0 0 e+00
8 . 8 8 8 8 8 5 e −08 1 . 6 5 5 6 8 2 e −23 0 . 0 0 0 0 0 0 e+00
0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
1 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
1 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
1 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
1 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 e+00
2 https://www.paraview.org/
3 https://docs.enthought.com/mayavi/mayavi/
4 https://wci.llnl.gov/simulation/computer-codes/visit/
5 https://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf
33
−1.333333 e+00
−6.666664 e −01
6 . 6 6 6 6 6 4 e −01
1 . 3 3 3 3 3 3 e+00
−1.333333 e+00
−6.666664 e −01
6 . 6 6 6 6 6 4 e −01
1 . 3 3 3 3 3 3 e+00
−1.333333 e+00
−6.666664 e −01
6 . 6 6 6 6 6 4 e −01
1 . 3 3 3 3 3 3 e+00
To these informations we must append 3 more datasets. The first one is the connectivity, the second one is the
offsets and the third one is the type. The first one is trivial since said connectivity is needed for the Finite Elements.
The second must be understood as follows: when reading the connectivity information in a linear manneer the offset
values indicate the beginning of each element (omitting the zero value). The third simply is the type of element as
given in the vtk format document (9 corresponds to a generic quadrilateral with an internal numbering consistent with
ours).
0 1 5 4
1 2 6 5
2 3 7 6
4 5 9 8
5 6 10 9
6 7 11 10
4
8
12
16
20
24
9
9
9
9
9
9
The file is then closed with
The solution.vtu file can then be opened with ParaView, MayaVi or Visit and the reader is advised to find tutorials
online on how to install and use these softwares.
34
To Do:
write about impose bc on el matrix
full compressible
total energy calculations
constraints
compositions, marker chain
van keken initial value with deformed mesh
free-slip bc on annulus and sphere . See for example p540 Gresho and Sani book.
non-linear rheologies (two layer brick spmw16, tosn15)
Picard vs Newton
markers
Schur complement approach
periodic boundary conditions
open boundary conditions
consistent boundary flux (CBF)
free surface
SUPG
zaleski disk advection
Busse convection pb, compare with aspect
cvi !!!
pure elastic
including phase changes (w. R. Myhill)
discontinuous galerkin
formatting of code style
navier-stokes ? (LUKAS)
compute strainrate in middle of element or at quad point for punch?
GEO1442 code
GEO1442 indenter setup in plane ?
in/out flow on sides for lith modelling
Problems to solve:
colorscale
better yet simple matrix storage ?
35
5
fieldstone: simple analytical solution
From [19]. In order to illustrate the behavior of selected mixed finite elements in the solution of stationary Stokes
flow, we consider a two-dimensional problem in the square domain Ω = [0, 1] × [0, 1], which possesses a closed-form
analytical solution. The problem consists of determining the velocity field v = (u, v) and the pressure p such that
−ν∆v + ∇p = b
∇·v =0
v=0
in Ω
in Ω
on Γ
where the fluid viscosity is taken as ν = 1. The components of the body force b are prescribed as
bx
(12 − 24y)x4 + (−24 + 48y)x3 + (−48y + 72y 2 − 48y 3 + 12)x2
=
+(−2 + 24y − 72y 2 + 48y 3 )x + 1 − 4y + 12y 2 − 8y 3
by
(8 − 48y + 48y 2 )x3 + (−12 + 72y − 72y 2 )x2
=
+(4 − 24y + 48y 2 − 48y 3 + 24y 4 )x − 12y 2 + 24y 3 − 12y 4
With this prescribed body force, the exact solution is
Note that the pressure obeys
features
R
Ω
u(x, y)
= x2 (1 − x)2 (2y − 6y 2 + 4y 3 )
v(x, y)
= −y 2 (1 − y)2 (2x − 6x2 + 4x3 )
p(x, y)
= x(1 − x) − 1/6
p dΩ = 0
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (no-slip)
• direct solver
• isothermal
• isoviscous
• analytical solution
36
vx
1.0
vy
1.0
0.010
0.8
0.8
0.005
0.4
0.4
0.2
0.6
0.8
0.01
x
0.6
0.8
1.0
0.000010
0.8
0.000005
0.6
yy
0.8
0.000010
0.000005
0.000010
0.2
0.4
x
0.6
0.8
1.0
0.04
0.04
0.2
0.4
x
0.6
0.8
1.0
p pth
1.0
0.00010
0.00005
0.8
0.00000
0.00005
y
0.000000
0.4
0.000005
xy
0.6
y
y
0.4
0.0
0.0
0.00010
0.4
0.000005
0.2
0.000010
0.2
0.4
x
0.6
0.8
1.0
0.100000
0.15
0.00
0.0
0.0
1.0
0.8
1.0
0.02
vy tyth
1.0
0.8
0.2
0.03
0.6
0.6
0.4
0.02
x
x
0.6
0.01
0.4
0.4
0.02
0.01
0.2
0.2
0.8
0.02
0.000000
0.0
0.0
0.10
1.0
0.03
0.6
0.2
0.2
0.0
0.0
1.0
0.00
0.0
0.0
vx txth
1.0
0.8
0.2
0.03
0.4
0.6
0.4
0.02
0.2
x
0.6
0.01
0.2
0.4
y
y
0.00
0.4
0.2
0.8
0.02
0.6
0.010
1.0
0.03
0.8
0.0
0.0
0.0
0.0
1.0
xx
1.0
0.05
y
0.010
x
0.4
0.005
0.005
0.2
0.4
0.00
0.6
0.000
y
y
0.000
y
0.6
0.2
0.05
0.8
0.005
0.6
0.0
0.0
p
1.0
0.010
0.00015
0.2
0.0
0.0
0.00020
0.2
0.4
x
0.6
0.8
1.0
0.00025
velocity
pressure
x2
x1
0.010000
error
0.001000
0.000100
0.000010
0.000001
0.01
0.1
h
Quadratic convergence for velocity error, linear convergence for pressure error, as expected.
ToDo:
pressure normalisation?
different cmat, a la schmalholz
To go further:
37
1. make your own analytical solution
38
6
fieldstone: Stokes sphere
Viscosity and density directly computed at the quadrature points.
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• buoyancy-driven flow
vx
1.0
0.8
0.8
0.0005
0.6
0.4
0.0005
0.0010
0.2
0.4
x
0.6
0.8
1.0
xx
1.0
0.8
0.2
0.6
0.2
0.0
0.0
0.2
0.2
0.002
0.2
0.4
x
0.6
0.8
0.4
0.0
0.0
1.0
yy
0.2
0.010
0.0
0.0
1.0
y
0.0050
0.2
0.010
0.2
0.4
x
0.6
0.8
density
0.4
x
0.6
0.8
0.005
0.010
0.015
0.2
0.4
x
0.6
0.8
1.0
viscosity
2.0
1.0
0.8
1.8
0.8
80
0.6
1.6
0.6
60
0.4
1.4
0.4
40
0.2
1.2
0.2
20
1.0
0.0
0.0
0.0025
0.2
0.000
0.0
0.0
1.0
y
0.0125
0.0075
0.005
100
y
0.0150
0.4
0.010
0.2
1.0
0.0175
0.0100
0.015
0.005
0.2
0.6
xy
0.4
0.005
0.8
1.0
y
0.000
0.4
x
0.8
0.6
y
y
0.4
0.8
x
0.6
0.8
0.6
0.6
0.4
0.005
0.6
0.4
0.2
1.0
0.010
0.005
0.2
0.0
0.4
0.8
0.000
p
0.4
0.001
1.0
0.010
0.8
0.0
0.0
0.001
0.000
0.4
0.2
1.0
y
y
0.0000
0.002
y
0.6
0.0
0.0
vy
1.0
0.0010
0.0
0.0
0.2
0.4
x
0.6
0.8
39
1.0
0.2
0.4
x
0.6
0.8
1.0
7
fieldstone: Convection in a 2D box
This benchmark deals with the 2-D thermal convection of a fluid of infinite Prandtl number in a rectangular closed
cell. In what follows, I carry out the case 1a, 1b, and 1c experiments as shown in [5]: steady convection with constant
viscosity in a square box.
The temperature is fixed to zero on top and to ∆T at the bottom, with reflecting symmetry at the sidewalls (i.e.
∂x T = 0) and there are no internal heat sources. Free-slip conditions are implemented on all boundaries.
The Rayleigh number is given by
αgy ∆T h3 ρ2 cp
αgy ∆T h3
=
(36)
Ra =
κν
kµ
In what follows, I use the following parameter values: Lx = Ly = 1,ρ0 = cP = k = µ = 1, T0 = 0, α = 10−2 ,
g = 102 Ra and I run the model with Ra = 104 , 105 and 106 .
The initial temperature field is given by
T (x, y) = (1 − y) − 0.01 cos(πx) sin(πy)
(37)
The perturbation in the initial temperature fields leads to a perturbation of the density field and sets the fluid in
motion.
Depending on the initial Rayleigh number, the system ultimately reaches a steady state after some time.
The Nusselt number (i.e. the mean surface temperature gradient over mean bottom temperature) is computed as
follows [5]:
R ∂T
∂y (y = Ly )dx
(38)
N u = Ly R
T (y = 0)dx
Note that in our case the denominator is equal to 1 since Lx = 1 and the temperature at the bottom is prescribed to
be 1.
Finally, the steady state root mean square velocity and Nusselt number measurements are indicated in Table ??
alongside those of [5] and [65]. (Note that this benchmark was also carried out and published in other publications
[72, 1, 27, 14, 45] but since they did not provide a complete set of measurement values, they are not included in the
table.)
4
Ra = 10
Ra = 105
Ra = 106
Vrms
Nu
Vrms
Nu
Vrms
Nu
Blankenbach et al
42.864947 ± 0.000020
4.884409 ± 0.000010
193.21454 ± 0.00010
10.534095 ± 0.000010
833.98977 ± 0.00020
21.972465 ± 0.000020
Tackley [65]
42.775
4.878
193.11
10.531
833.55
21.998
Steady state Nusselt number N u and Vrms measurements as reported in the literature.
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• Boussinesq approximation
• direct solver
• non-isothermal
• buoyancy-driven flow
• isoviscous
• CFL-condition
40
0.4
20
0.2
0.2
0.4
x
0.6
0.8
40 0.0
0.0
1.0
0.4
100
0.8
x
0.6
0.8
1.0
yy
1.0
0
100
0.4
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
0.2
200 0.00.0
0.25
0.8
100
0.2
0.4
x
0.6
0.8
1.0
qy
1.0
1.00 0.4
3
0.2
1.25 0.2
2
y
0.4
x
0.6
0.8
1.0
1.50 0.0
0.0
1.75
1
0.2
0.4
x
0.6
0.8
1.0
0
0.6
0.4
0.4
y
0.6
0.0
400000 0.0
1.0
40
x
0.6
0.8
1.0
0.0
rho
0.998
0.996
0.6
0.4
0
0.994
10 0.2
0.2
0.4
x
0.6
0.8
1.0
20 0.0
0.0
30
80
7
70
6
0.2
0.4
x
0.6
0.8
1.0
0.992
5
50
40
30
4
3
20
2
10
0
0.4
0.8
20
10
0.2
0.2
1.0
30
60
5
4
0.4
0.8
xy
0.0
0.0
6
0.8
0.75 0.6
0.2
x
0.6
0.8
7
0.6
0.0
0.0
0.4
0.2
y
0.50
0.2
0.4
0
0.00
qx
1.0
0.0
0.0
0.6
0.4
100
0.8
200000
0.2
0.2
y
0.6
y
y
60
0.4
200 1.0
0.8
0.6
0
y
40
0.2
60
xx
1.0
0.6
20
0.2
T
2000000.8
Nu
0.4
0
y
y
0
0.8
20
0.6
400000
1.0
p
1.0
40
0.8
20
0.6
0.0
0.0
1.0
40
0.8
vrms
1.0
1.0
60
vy
y
60
vx
1
0.00 0.05 0.10 0.15 0.20 0.25 0.30
time
ToDo:
implement steady state criterion
reach steady state
do Ra=1e4, 1e5, 1e6
plot against blankenbach paper and aspect
look at critical Ra number
41
0.00 0.05 0.10 0.15 0.20 0.25 0.30
time
8
fieldstone: The lid driven cavity
The lid driven cavity is a famous Computational Fluid Dynamics test case and has been studied in countless publications with a wealth of numerical techniques (see [24] for a succinct review) and also in the laboratory [41].
It models a plane flow of an isothermal isoviscous fluid in a rectangular (usually square) lid-driven cavity. The
boundary conditions are indicated in the Fig. ??a. The gravity is set to zero.
8.1
the lid driven cavity problem (ldc=0)
In the standard case, the upper side of the cavity moves in its own plane at unit speed, while the other sides are fixed.
This thereby introduces a discontinuity in the boundary conditions at the two upper corners of the cavity and yields
an uncertainty as to which boundary (side or top) the corner points belong to. In this version of the code the top
corner nodes are considered to be part of the lid. If these are excluded the recovered pressure showcases and extremely
large checkboard pattern.
This benchmark is usually dicussed in the context of low to very high Reynolds number with the full Navier-Stokes
equations being solved (with the noticeable exception of [59, 60, 9, 23] which focus on the Stokes equation). In the
case of the incompressible Stokes flow, the absence of inertia renders this problem instantaneous so that only one time
step is needed.
8.2
the lid driven cavity problem - regularisation I (ldc=1)
We avoid the top corner nodes issue altogether by prescribing the horizontal velocity of the lid as follows:
u(x) = x2 (1 − x)2 .
(39)
In this case the velocity and its first derivative is continuous at the corners. This is the so-called regularised lid-driven
cavity problem [53].
8.3
the lid driven cavity problem - regularisation II (ldc=2)
Another regularisation was presented in [16]. Here, a regularized lid driven cavity is studied which is consistent in the
sense that ∇ · v = 0 holds also at the corners of the domain. There are no-slip conditions at the boundaries x = 0,
x = 1, and y = 0.
The velocity at y = 1 is given by
2
x1 − x
π)
x ∈ [0, x1 ]
1 − cos(
x1
u(x) = 1
x ∈ [x1 , 1 − x1 ]
2
1
x − (1 − x1 )
u(x) = 1 −
1 − cos(
π)
x ∈ [1 − x1 , 1]
4
x1
u(x)
=
1−
1
4
Results are obtained with x1 = 0.1.
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• isothermal
• isoviscous
42
(40)
A 100x100 element grid is used. No-slip boundary conditions are prescribed on sides and bottom. A zero vertical
velocity is prescribed at the top and the exact form of the prescribed horizontal velocity is controlled by the ldc
parameter.
200
ldc0
ldc1
ldc2
150
pressure
100
50
0
-50
-100
-150
-200
0
0.1
0.2
0.3
0.4
0.5
x
43
0.6
0.7
0.8
0.9
1
9
fieldstone: solcx benchmark
The SolCx benchmark is intended to test the accuracy of the solution to a problem that has a large jump in the
viscosity along a line through the domain. Such situations are common in geophysics: for example, the viscosity in a
cold, subducting slab is much larger than in the surrounding, relatively hot mantle material.
The SolCx benchmark computes the Stokes flow field of a fluid driven by spatial density variations, subject to a
spatially variable viscosity. Specifically, the domain is Ω = [0, 1]2 , gravity is g = (0, −1)T and the density is given by
ρ(x, y) = sin(πy) cos(πx)
(41)
Boundary conditions are free slip on all of the sides of the domain and the temperature plays no role in this benchmark.
The viscosity is prescribed as follows:
1
f or x < 0.5
µ(x, y) =
(42)
106 f or x > 0.5
Note the strongly discontinuous viscosity field yields a stagnant flow in the right half of the domain and thereby yields
a pressure discontinuity along the interface.
The SolCx benchmark was previously used in [22] (references to earlier uses of the benchmark are available there)
and its analytic solution is given in [77]. It has been carried out in [42] and [28]. Note that the source code which
evaluates the velocity and pressure fields for both SolCx and SolKz is distributed as part of the open source package
Underworld ([50], http://underworldproject.org).
In this particular example, the viscosity is computed analytically at the quadrature points (i.e. tracers are not
used to attribute a viscosity to the element). If the number of elements is even in any direction, all elements (and
their associated quadrature points) have a constant viscosity(1 or 106 ). If it is odd, then the elements situated at the
viscosity jump have half their integration points with µ = 1 and half with µ = 106 (which is a pathological case since
the used quadrature rule inside elements cannot represent accurately such a jump).
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• analytical solution
44
1.0
0.0010
0.8
0.002
0.8
0.6
0.001
0.6
0.6
y
0.0000
0.4
y
0.0005
0.0005
0.0010
0.2
0.0
0.0
0.0015
0.0
0.0
0.4
x
0.6
0.8
1.0
xx
1.0
0.0075
0.8
0.0025
y
0.0000
0.4
0.0
0.0
0.0050
0.0075
0.2
0.4
x
0.6
0.8
1.0
vx txth
1.0
0.4
x
0.6
0.8
0.0100
yy
0.0050
0.0075
0.4
x
0.6
0.8
1.0
vy tyth
0.6
0.8
1.0
0.005
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
p pth
0.010
0.020
0.4
0.00000200.4
0.005
0.0000025
0.2
0.0000030
0.0
0.0
0.010
0.010
0.000
y
0.0
0.0
0.6
0.000
0.2
0.005
0.000003 0.0
0.0
x
0.0100
0.005
0.4
0.6
0.0000015
0.000002 0.2
0.4
0.010
0.6
0.2
0.2
0.020
0.015
y
y
0.000001
0.025
0.0000005
0.8
0.0000010
0.000000
0.4
1.0
xy
0.0000000 1.0
0.000002 0.8
0.000001
0.8
0.6
0.0025
0.2
x
0.6
0.015
0.0025
0.2
0.4
0.8
0.0050
0.4
0.2
1.0
0.0100
0.0000
0.0
0.0
0.2
0.0
0.0
0.0075
0.000003 1.0
0.8
0.1
1.0
0.6
0.0025
0.2
0.2
y
0.6
0.0
0.4
0.002
0.8
0.0050
0.1
0.2
1.0
0.0100
0.2
0.001
0.2
0.2
0.000
0.4
p
1.0
0.003
y
0.8
vy
0.0015
y
vx
1.0
0.2
0.4
x
0.6
0.8
What we learn from this
45
1.0
0.015
0.2
0.4
x
0.6
0.8
1.0
0.020
10
fieldstone: solkz benchmark
The SolKz benchmark [57] is similar to the SolCx benchmark. but the viscosity is now a function of the space
coordinates:
µ(y) = exp(By) with B = 13.8155
(43)
It is however not a discontinuous function but grows exponentially with the vertical coordinate so that its overall
variation is again 106 . The forcing is again chosen by imposing a spatially variable density variation as follows:
ρ(x, y) = sin(2y) cos(3πx)
(44)
Free slip boundary conditions are imposed on all sides of the domain. This benchmark is presented in [77] as well and
is studied in [22] and [28].
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
y
0.6
0.4
0.2
0.2
0.4
x
0.6
0.8
vx txth
1.0
0.4
0.6
0.8
0.0
0.0
0.2
101
x
0.6
0.8
1.0
0.8
0.4
0.2
x
error
y
0.2
0.4
x
0.6
0.8
1.0
0.0008
0.0006
0.0004
0.0002
0.0000
0.0002
0.0004
0.0006
0.0008
0.6
0.8
10
4
10
5
0.4
10
6
0.2
10
7
0.0
0.0
1.0
0.00010000
0.00001000
0.00000100
0.00000010
0.01
p pth
0.6
0.00100000
0.00000001
1.0
0.8
velocity
pressure
x2
x1
0.01000000
0.8
0.0004
0.0003
0.0002
0.0001
0.0000
0.0001
0.0002
0.0003
0.0004
1.0
0.6
0.10000000
0.0
0.0
0.75
0.50
0.25
0.00
0.25
0.50
0.75
0.4
0.6
0.2
density
0.2
x
0.4
1.0
0.8
0.0
0.0
0.4
0.6
1.0
0.6
0.2
0.8
0.5
x
1.0
xy
1.0
0.5
0.4
0.8
0.2
1e 7
0.2
0.6
0.4
0.0
0.0
0.0
x
0.6
1.0
0.2
0.4
0.8
1.0
y
y
102
0.4
0.8
0.4
0.0
0.0
103
0.2
0.6
0.6
104
0.0
0.0
x
vy tyth
1.0
0.2
0.4
0.2
1.0
0.00075
0.00050
0.00025
0.00000
0.00025
0.00050
0.00075
0.8
105
0.4
1.0
yy
1.0
1.0
0.6
0.8
0.2
1.0
0.8
0.6
0.4
2
4
6
0.2
x
0.0
0.0
x
y
y
0.6
0.4
0.00006
0.6
1e 7
0.8
0.2
0.2
0.0
0.0
0.8
1.0
6
4
2
0
0.0
0.0
0.00004
0.4
0.075
0.050
0.025
0.000
0.025
0.050
0.075
0.4
0.00002
0.2
p
0.6
0.2
y
0.8
0.0
0.0
0.00000
0.4
1.0
0.00075
0.00050
0.00025
0.00000
0.00025
0.00050
0.00075
0.8
0.00002
0.6
xx
1.0
1.0
0.00004
y
0.4
0.8
0.00006
y
y
0.6
vy
1.0
y
0.8
0.000100
0.000075
0.000050
0.000025
0.000000
0.000025
0.000050
0.000075
0.000100
y
vx
1.0
0.1
h
46
0.2
0.4
x
0.6
0.8
1.0
11
fieldstone: solvi benchmark
Following SolCx and SolKz, the SolVi inclusion benchmark solves a problem with a discontinuous viscosity field, but
in this case the viscosity field is chosen in such a way that the discontinuity is along a circle. Given the regular nature
of the grid used by a majority of codes and the present one, this ensures that the discontinuity in the viscosity never
aligns to cell boundaries. This in turns leads to almost discontinuous pressures along the interface which are difficult to
represent accurately. [62] derived a simple analytic solution for the pressure and velocity fields for a circular inclusion
under simple shear and it was used in [17], [64], [22], [42] and [28].
Because of the symmetry of the problem, we only have to solve over the top right quarter of the domain (see Fig.
??a).
The analytical solution requires a strain rate boundary condition (e.g., pure shear) to be applied far away from the
inclusion. In order to avoid using very large domains and/or dealing with this type of boundary condition altogether,
the analytical solution is evaluated and imposed on the boundaries of the domain. By doing so, the truncation error
introduced while discretizing the strain rate boundary condition is removed.
A characteristic of the analytic solution is that the pressure is zero inside the inclusion, while outside it follows the
relation
µm (µi − µm ) ri2
cos(2θ)
(45)
pm = 4˙
µi + µm r2
where µi = 103 is the viscosity of the inclusion and µm = 1 is the viscosity of the background media, θ = tan−1 (y/x),
and ˙ = 1 is the applied strain rate.
[17] thoroughly investigated this problem with various numerical methods (FEM, FDM), with and without tracers,
and conclusively showed how various averagings lead to different results. [22] obtained a first order convergence for
both pressure and velocity, while [42] and [28] showed that the use of adaptive mesh refinement in respectively the
FEM and FDM yields convergence rates which depend on refinement strategies.
47
vx
1.0
vy
1.0
0.0
0.2
0.8
0.6
0.6
0.6
0.4
0.6
0.4
0.4
0.4
0.6
0.4
0.2
0.2
0.2
0.8
0.2
0.0
0.0
0.0
0.0
0.0
0.8
1.0
xx
1.0
0.8
1.75
1.0
1.50
0.8
1.25
0.6
0.25
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
0.00
vx txth
1.0
0.000
0.4
x
0.6
1.50
0.2
0.0
0.0
1.75
0.0
0.0
0.8
0.8
1.0
vy tyth
103
0.6
0.8
0.8
102
y
0.6
y
0.6
0.4
0.4
101
0.2
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
100
0.0
0.0
10.00000000
0.2
0.4
x
0.6
0.8
1.0
velocity
pressure
x2
x1
1.00000000
0.10000000
error
0.01000000
0.00100000
0.00010000
0.00001000
0.00000100
0.00000010
0.00000001
0.01
0.0
0.0
0.1
h
48
1.100
1.075
1.050
1.025
1.000
0.975
0.950
0.925
0.900
0.4
x
0.6
0.8
1.0
p pth
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
15
10
5
0
5
0.2
1.0
density
1.0
0.8
x
0.2
15
1.0
0.4
0.000
0.4
0.8
0.6
0.005
0.2
0.6
0.8
0.015
0.010
x
xy
1.0
0.020
0.0
0.0
1.0
1.0
0.6
y
y
0.2
0.2
x
0.4
0.4
1.25
0.4
0.2
0.6
0.75
0.2
10
0.8
0.50
0.2
0.020
0.0
0.0
0.25
1.00
5
1.0
0.00
0.4
0.015
0.2
0.0
0.0
1.0
0.6
0.010
0.4
0.8
0.8
0.005
0.6
0.6
yy
1.0
0.8
x
0.4
0.50
0.2
0.4
y
y
0.75
0.4
0.2
0.6
1.00
0
y
0.6
5
y
x
10
10
0.2
0.4
x
0.6
0.8
1.0
1.0
15
100
0.8
0.6
10
1
10
2
10
3
y
0.4
15
y
0.8
y
0.8
y
0.8
0.2
p
1.0
0.4
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
12
fieldstone: the indentor benchmark
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.4
x
0.6
0.8
1.0
1.0
10.0
7.5
5.0
2.5
0.0
2.5
5.0
1.0
20.0
17.5
15.0
12.5
10.0
7.5
5.0
2.5
xx
0.2
0.4
0.2
0.4
x
x
0.6
0.6
0.8
0.8
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.2
vy
0.0
0.2
0.4
y
y
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.6
0.2
0.4
x
0.6
0.8
1.0
0.4
x
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.6
0.8
1.0
5.0
2.5
0.0
2.5
5.0
7.5
10.0
0.5
0.4
0.3
0.2
0.1
0.0
0.0
6
p
5
4
3
2
0.2
0.4
1.0
yy
0.2
0.8
y
0.2
0.3
0.2
0.1
0.0
0.1
0.2
0.3
x
0.6
0.8
1.0
xy
0.2
0.4
x
0.6
0.8
1.0
103
0.5
0.4
0.3
0.2
0.1
0.0
0.0
102
101
100
0.2
0.4
x
0.6
0.8
1.0
y
0.5
0.4
0.3
0.2
0.1
0.0
0.0
vx
y
0.5
0.4
0.3
0.2
0.1
0.0
0.0
y
y
y
y
The punch benchmark is one of the few boundary value problems involving plastic solids for which there exists an
exact solution. Such solutions are usually either for highly simplified geometries (spherical or axial symmetry, for
instance) or simplified material models (such as rigid plastic solids) [39].
In this experiment, a rigid punch indents a rigid plastic half space; the slip line field theory gives exact solutions
as shown in Fig. ??a. The plane strain formulation of the equations and the detailed solution to the problem were
derived in the Appendix of [70] and are also presented in [26].
The two dimensional punch problem has been extensively studied numerically for the past 40 years [80, 79, 11,
10, 36, 76, 6, 55] and has been used to draw a parallel with the tectonics of eastern China in the context of the
India-Eurasia collision [67, 49]. It is also worth noting that it has been carried out in one form or another in series of
analogue modelling articles concerning the same region, with a rigid indenter colliding with a rheologically stratified
lithosphere [52, 15, 38].
Numerically, the one-time step punch experiment is performed on a two-dimensional domain of purely plastic
von Mises material. Given that the von Mises rheology yield criterion does not depend on pressure, the density of
the material and/or the gravity vector is set to zero. Sides are set to free slip boundary conditions, the bottom to
no slip, while a vertical velocity (0, −vp ) is prescribed at the top boundary for nodes whose x coordinate is within
[Lx /2 − δ/2, Lx /2 + δ/2].
The following parameters are used: Lx = 1, Ly = 0.5, µmin = 10−3 , µmax = 103 , vp = 1, δ = 0.123456789 and the
yield value of the material is set to k = 1.
The analytical solution predicts that the angle of the shear bands stemming from the sides of the punch is π/4,
that√the pressure right under the punch is 1 + π, and that the velocity of the rigid blocks on each side of the punch is
vp / 2 (this is simply explained by invoking conservation of mass).
10
1
10
2
10
3
0.5
0.4
0.3
0.2
0.1
0.0
0.0
1
20
15
10
5
0
5
10
15
20
5
p (nodal)
4
3
2
0.2
0.4
x
0.6
0.8
1.0
1
6
0.0
5
0.2
4
0.00
p
0.2
0.05
v
u
0.4
0.10
0.4
3
0.6
0.05
2
0.8
0.10
1
1.0
0.0
0.2
0.4
x
0.6
0.8
1.0
0.0
0.2
0.4
x
0.6
49
0.8
1.0
0.0
0.2
0.4
x
0.6
0.8
1.0
ToDo: smooth punch
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (no-slip)
• isothermal
• non-isoviscous
• nonlinear rheology
50
13
fieldstone: the annulus benchmark
This benchmark is based on Thieulot & Puckett [Subm.] in which an analytical solution to the isoviscous incompressible
Stokes equations is derived in an annulus geometry. The velocity and pressure fields are as follows:
vr (r, θ)
=
g(r)k sin(kθ),
(46)
vθ (r, θ)
=
f (r) cos(kθ),
(47)
p(r, θ)
=
kh(r) sin(kθ),
(48)
ρ(r, θ)
=
ℵ(r)k sin(kθ),
(49)
with
f (r)
g(r)
h(r)
ℵ(r)
A
B
= Ar + B/r,
B
C
A
r + ln r + ,
=
2
r
r
2g(r) − f (r)
=
,
r
g0
f
f0
g
= g 00 − − 2 (k 2 − 1) + 2 + ,
r
r
r
r
2(ln R1 − ln R2 )
,
= −C 2
R2 ln R1 − R12 ln R2
R22 − R12
= −C 2
.
R2 ln R1 − R12 ln R2
(50)
(51)
(52)
(53)
(54)
(55)
The parameters A and B are chosen so that vr (R1 ) = vr (R2 ) = 0, i.e. the velocity is tangential to both inner
and outer surfaces. The gravity vector is radial and of unit length. In the present case, we set R1 = 1, R2 = 2 and
C = −1.
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions
• direct solver
• isothermal
• isoviscous
• analytical solution
• annulus geometry
• elemental boundary conditions
51
10
velocity
pressure
x2
x1
1
error
0.1
0.01
0.001
0.0001
0.01
h
52
14
fieldstone: stokes sphere (3D) - penalty
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• 3D
• elemental b.c.
• buoyancy driven
resolution is 24x24x24
53
15
fieldstone: stokes sphere (3D) - mixed formulation
This is the same setup as Section 15.
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• 3D
• elemental b.c.
• buoyancy driven
54
16
fieldstone: consistent pressure recovery
p = −λ∇ · v
q1 is smoothed pressure obtained with the center-to-node approach.
q2 is recovered pressure obtained with [78]. R
All three fulfill the zero average condition: pdΩ = 0.
y
0.000
0.4
0.2
0.4
x
0.6
0.8
vx txth
1.0
0.000002
x
0.6
0.8
0.8
0.4
0.2
0.10
0.0
0.0
0.15
x
0.6
0.8
1.0
0.6
y
0.6
0.4
0.2
0.2
0.4
x
0.6
0.8
1.0
0.10
0.2
0.15
0.2
0.4
x
0.6
0.8
1.0
q2 pth
1.0
0.00002
0.00000
0.00002
0.00004
0.00006
0.2
0.4
0.004
0.8
0.003
1.0
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
1.0
0.8
0.002
0.6
0.4
0.001
0.4
0.2
0.000
0.2
0.2
0.4
x
0.6
0.8
1.0
0.001
0.0
0.0
velocity
pressure (el)
pressure (q1)
pressure (q2)
x2
x1
0.010000
error
0.001000
0.000100
0.000010
0.01
0.8
0.8
0.6
0.0
0.0
0.100000
0.000001
x
0.6
y
0.8
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000
0.05
0.4
q1 pth
1.0
0.00
0.0
0.0
0.15
1.0
0.05
y
y
0.05
0.4
q2
1.0
p pth
0.0
0.0
1.0
0.8
0.00
0.6
x
0.8
0.2
0.000002
0.8
0.6
0.4
0.000001
0.6
x
0.6
0.000000
0.4
0.4
0.8
0.000001
0.2
0.2
1.0
0.000002
1.0
0.05
0.2
vy tyth
0.0
0.0
1.0
0.0
0.0
1.0
0.2
q1
1.0
0.8
0.4
0.000002
0.4
0.6
0.6
0.000001
0.2
x
y
y
0.000000
0.4
0.4
0.8
0.000001
0.6
0.2
0.2
0.10
0.2
0.010
1.0
0.8
0.0
0.0
0.0
0.0
1.0
0.05
0.4
0.005
0.2
0.010
0.00
0.6
0.000
0.4
0.005
0.2
0.0
0.0
0.6
0.05
0.8
0.005
y
0.6
0.0
0.0
0.8
0.005
p
1.0
0.010
y
0.8
vy
1.0
0.010
y
vx
1.0
0.1
h
In terms of pressure error, q2 is better than q1 which is better than elemental.
QUESTION: why are the averages exactly zero ?!
TODO:
• add randomness to internal node positions.
• look at elefant algorithms
55
17
fieldstone: the Particle in Cell technique (1) - the effect of averaging
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (no-slip)
• isothermal
• non-isoviscous
• particle-in-cell
After the initial setup of the grid, markers can then be generated and placed in the domain. One could simply
randomly generate the marker positions in the whole domain but unless a very large number of markers is used, the
chance that an element does not contain any marker exists and this will prove problematic. In order to get a better
control over the markers spatial distribution, one usually generates the marker per element, so that the total number
of markers in the domain is the product of the number of elements times the user-chosen initial number of markers
per element.
Our next concern is how to actually place the markers inside an element. Two methods come to mind: on a regular
grid, or in a random manner, as shown on the following figure:
In both cases we make use of the basis shape functions: we generate the positions of the markers (random or
regular) in the reference element first (rim , sim ), and then map those out to the real element as follows:
xim =
m
X
Ni (rim , sim )xi
yim =
i
m
X
Ni (rim , sim )yi
i
where xi , yi are the coordinates of the vertices of the element.
When using active markers, one is faced with the problem of transferring the properties they carry to the mesh on
which the PDEs are to be solved. As we have seen, building the FE matrix involves a loop over all elements, so one
simple approach consists of assigning each element a single property computed as the average of the values carried by
the markers in that element. Often in colloquial language ”average” refers to the arithmetic mean:
n
hφiam =
1X
φi
n
k
where < φ >am is the arithmetic average of the n numbers φi . However, in mathematics other means are commonly
used, such as the geometric mean:
!
n
Y
hφigm =
φi
i
and the harmonic mean:
n
hφihm =
1X 1
n i φi
!−1
Furthermore, there is a well known inequality for any set of positive numbers,
hφiam
≥
hφigm
56
≥
hφihm
which will prove to be important later on.
Let us now turn to a simple concrete example: the 2D Stokes sphere. There are two materials in the domain, so
that markers carry the label ”mat=1” or ”mat=2”. For each element an average density and viscosity need to be
computed. The majority of elements contains markers with a single material label so that the choice of averaging does
not matter (it is trivial to verify that if φi = φ0 then hφiam = hφigm = hφihm = φ0 . Remain the elements crossed by
the interface between the two materials: they contain markers of both materials and the average density and viscosity
inside those depends on 1) the total number of markers inside the element, 2) the ratio of markers 1 to markers 2, 3)
the type of averaging.
This averaging problem has been studied and documented in the literature [61, 17, 68, 54]
Nodal projection. Left: all markers inside elements to which the green node belongs to are taken into account.
Right: only the markers closest to the green node count.
The setup is identical as the Stokes sphere experiment. The bash script runall runs the code for many resolutions,
both initial marker distribution and all three averaging types. The viscosity of the sphere has been set to 103 while
the viscosity of the surrounding fluid is 1. The average density is always computed with an arithmetic mean. Root
mean square velocity results are shown hereunder:
random markers, elemental proj
regular markers, elemental proj
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
0.009
0.008
0.007
vrms
0.006
0.005
0.004
0.003
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.01
0.009
0.008
0.007
0.006
vrms
0.01
0.005
0.004
0.003
0.002
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.002
0.001
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0
0.01
0.02
0.03
0.04
h
0.05
0.06
0.07
0.08
0.09
0.1
h
random markers, nodal proj
regular markers, nodal proj
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
0.009
0.008
0.007
vrms
0.006
0.005
0.004
0.003
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.01
0.009
0.008
0.007
0.006
vrms
0.01
0.005
0.004
0.003
0.002
0.002
0.001
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0
0.01
0.02
0.03
0.04
h
0.05
0.06
0.07
0.08
0.09
0.1
h
random markers, nodal proj
regular markers, nodal proj
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
0.009
0.008
0.007
vrms
0.006
0.005
0.004
0.003
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.01
0.009
0.008
0.007
0.006
vrms
0.01
0.005
0.004
0.003
0.002
0.002
0.001
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0
h
0.01
0.02
0.03
0.04
0.05
h
57
0.06
0.07
0.08
0.09
0.1
Conclusions:
• With increasing resolution (h → 0) vrms values seem to converge towards a single value, irrespective of the
number of markers.
• At low resolution, say 32x32 (i.e. h=0.03125), vrms values for the three averagings differ by about 10%. At
higher resolution, say 128x128, vrms values are still not converged.
• The number of markers per element plays a role at low resolution, but less and less with increasing resolution.
• Results for random and regular marker distributions are not identical but follow a similar trend and seem to
converge to the same value.
• At low resolutions, elemental values yield better results.
• harmonic mean yields overal the best results
64 markers per element
0.01
0.009
am, reg, el
gm, reg, el
hm, reg, el
am, rand, el
gm, rand, el
hm, rand, el
am, reg, nod1
gm, reg, nod1
hm, reg, nod1
am, rand, nod1
gm, rand, nod1
hm, rand, nod1
am, reg, nod2
gm, reg, nod2
hm, reg, nod2
am, rand, nod2
gm, rand, nod2
hm, rand, nod2
0.008
0.007
vrms
0.006
0.005
0.004
0.003
0.002
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
h
58
0.08
0.09
0.1
18
fieldstone: solving the full saddle point problem
The details of the numerical setup are presented in Section 5.
The main difference is that we no longer use the penalty formulation and therefore keep both velocity and pressure
as unknowns. Therefore we end up having to solve the following system:
K G
V
f
·
=
or,
A · X = rhs
P
h
GT 0
Each block K, G and vector f , h are built separately in the code and assembled into the matrix A and vector rhs
afterwards. A and rhs are then passed to the solver. We will see later that there are alternatives to solve this approach
which do not require to build the full Stokes matrix A.
Each element has m = 4 vertices so in total ndof V × m = 8 velocity dofs and a single pressure dof, commonly
situated in the center of the element. The total number of velocity dofs is therefore N f emV = nnp × ndof V while
the total number of pressure dofs is N f emP = nel. The total number of dofs is then N f em = N f emV + N f emP .
As a consequence, matrix K has size N f emV, N f emV and matrix G has size N f emV, N f emP . Vector f is of size
N f emV and vector h is of size N f emP .
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (no-slip)
• direct solver (?)
• isothermal
• isoviscous
• analytical solution
• pressure smoothing
59
0.010
vx
1.0
0.010
vy
1.0
0.005 0.8
0.005 0.8
0.6
0.6
0.6
0.4
0.0
0.0
0.4
0.0050.2
0.2
0.2
0.4
x
0.6
0.8
1.0
0.00
0.05
0.10
y
y
y
0.000
0.05
0.0
0.010 0.0
0.4
0.15
0.0050.2
0.2
0.4
x
0.6
0.8
1.0
0.00
q
1.0
0.05
0.8
0.10
0.6
y
0.8
0.000
p
1.0
0.0
0.010 0.0
0.20
0.2
0.4
x
0.6
0.8
1.0
0.15
0.4
0.20
0.2
0.25 0.00.0
0.2
0.4
x
0.6
0.8
1.0
0.06
0.8
0.01
0.6
0.02
0.8
0.01
0.6
0.00
y
y
0.00
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.4
x
0.6
0.8
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.0000010
1.0
vx txth
1.0
0.03
1.0
0.4
x
0.6
0.8
1.0
0.03
0.0
0.0
0.0000005
0.2
0.4
x
0.6
0.8
1.0
0.0
0.0
0.0000010
0.2
0.4
x
0.6
0.8
1.0
0.4
x
0.6
0.8
1.0
0.04
0.6
0.03
0.0
0.0
0.0000010
0.4
0.02
0.2
0.04 0.0
0.0
0.2
0.00
p pth
0.02
0.04
0.06
0.08
0.0000005
0.2
0.2
0.05
0.8
0.4
0.06
y
0.6
0.0000000
0.4
y
0.6
0.0000000
0.4
y
0.6
0.2
0.02
0.2
0.8
0.0000005
0.0
0.0
0.00
0.4
0.8
0.0000005
0.2
0.02
0.6
0.8
0.4
0.04 1.0
0.8
0.0000010
1.0
vy tyth
0.06
y
0.02
xy
x
0.6
0.8
1.0
q pth
1.0
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
0.02
0.8
0.04
0.6
0.06
0.4
0.08
0.2
0.0
0.0
0.01
0.00
y
1.0
0.03
yy
y
0.03
xx
0.25
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
Unlike the results obtained with the penalty formualtion (see Section 5), the pressure showcases a very strong
checkerboard pattern, similar to the one in [20].
Left: pressure solution as shown in [20]; Right: pressure solution obtained with fieldstone.
Rather interestingly, the nodal pressure (obtained with a simple center-to-node algorithm) fails to recover a correct
pressure at the four corners.
60
19
fieldstone: solving the full saddle point problem in 3D
When using Q1 ×P0 elements, this benchmark fails because of the Dirichlet b.c. on all 6 sides and all three components.
However, as we will see, it does work well with Q2 × Q1 elements. .
This benchmark begins by postulating a polynomial solution to the 3D Stokes equation [18]:
x + x2 + xy + x3 y
y + xy + y 2 + x2 y 2
v=
(56)
−2z − 3xz − 3yz − 5x2 yz
and
p = xyz + x3 y 3 z − 5/32
(57)
While it is then trivial to verify that this velocity field is divergence-free, the corresponding body force of the Stokes
equation can be computed by inserting this solution into the momentum equation with a given viscosity µ (constant
or position/velocity/strain rate dependent). The domain is a unit cube and velocity boundary conditions simply use
Eq. (56). Following [7], the viscosity is given by the smoothly varying function
µ = exp(1 − β(x(1 − x) + y(1 − y) + z(1 − z)))
(58)
One can easily show that the ratio of viscosities µ? in the system follows µ? = exp(−3β/4) so that choosing β = 10
yields µ? ' 1808 and β = 20 yields µ? ' 3.269 × 106 .
We start from the momentum conservation equation:
˙ =f
−∇p + ∇ · (2µ)
The x-component of this equation writes
fx
∂
∂
∂
∂p
+
(2µ˙xx ) +
(2µ˙xy ) +
(2µ˙xz )
∂x ∂x
∂y
∂z
∂p
∂
∂
∂
∂µ
∂µ
∂µ
= −
+ 2µ ˙xx + 2µ ˙xy + 2µ ˙xz + 2 ˙xx + 2 ˙xy + 2 ˙xz
∂x
∂x
∂y
∂z
∂x
∂y
∂z
= −
Let us compute all the block separately:
˙xx
=
1 + 2x + y + 3x2 y
˙yy
=
1 + x + 2y + 2x2 y
˙zz
= −2 − 3x − 3y − 5x2 y
2˙xy
=
(x + x3 ) + (y + 2xy 2 ) = x + y + 2xy 2 + x3
2˙xz
=
(0) + (−3z − 10xyz) = −3z − 10xyz
2˙yz
=
(0) + (−3z − 5x2 z) = −3z − 5x2 z
In passing, one can verify that ˙xx + ˙yy + ˙zz = 0. We further have
∂
2˙xx
∂x
∂
2˙xy
∂y
∂
2˙xz
∂z
∂
2˙xy
∂x
∂
2˙yy
∂y
∂
2˙yz
∂z
∂
2˙xz
∂x
∂
2˙yz
∂y
∂
2˙zz
∂z
=
2(2 + 6xy)
=
1 + 4xy
= −3 − 10xy
=
1 + 2y 2 + 3x2
=
2(2 + 2x2 )
= −3 − 5x2
= −10yz
=
0
=
2(0)
61
(59)
(60)
∂p
∂x
∂p
∂y
∂p
∂z
= yz + 3x2 y 3 z
(61)
= xz + 3x3 y 2 z
(62)
= xy + x3 y 3
(63)
Pressure normalisation Here again, because Dirichlet boundary conditions are prescribed on all sides the pressure
is known up to an arbitrary constant. This constant can be determined by (arbitrarily) choosing to normalised the
pressure field as follows:
Z
p dΩ = 0
(64)
Ω
This is a single constraint associated to a single Lagrange
K
G
GT
0
0 CT
multiplier λ and the global Stokes system takes the form
0
V
C P
λ
0
In this particular case the constraint matrix C is a vector and it only acts on the pressure degrees of freedom because
of Eq.(64). Its exact expression is as follows:
Z
X
XZ X p
X X Z
XZ
p
Ni dΩ pi =
Ce · pe
p dΩ =
Ni pi dΩ =
p dΩ =
Ω
e
Ωe
e
Ωe
e
i
i
Ωe
e
where pe is the list of pressure dofs of element e. The elemental constraint vector contains the corresponding pressure
basis functions integrated over the element. These elemental constraints are then assembled into the vector C.
19.0.1
Constant viscosity
Choosing β = 0 yields a constant velocity µ(x, y, z) = exp(1) ' 2.718 (and greatly simplifies the right-hand side) so
that
∂
µ(x, y, z) = 0
∂x
∂
µ(x, y, z) = 0
∂y
∂
µ(x, y, z) = 0
∂z
and
fx
=
=
=
fy
=
=
=
fz
=
=
∂p
∂
∂
∂
+ 2µ ˙xx + 2µ ˙xy + 2µ ˙xz
∂x
∂x
∂y
∂z
−(yz + 3x2 y 3 z) + 2(2 + 6xy) + (1 + 4xy) + (−3 − 10xy)
−
−(yz + 3x2 y 3 z) + µ(2 + 6xy)
∂
∂
∂
∂p
+ 2µ ˙xy + 2µ ˙yy + 2µ ˙yz
−
∂y
∂x
∂y
∂z
−(xz + 3x3 y 2 z) + µ(1 + 2y 2 + 3x2 ) + µ2(2 + 2x2 ) + µ(−3 − 5x2 )
−(xz + 3x3 y 2 z) + µ(2 + 2x2 + 2y 2 )
∂p
∂
∂
∂
−
+ 2µ ˙xz + 2µ ˙yz + 2µ ˙zz
∂z
∂x
∂y
∂z
−(xy + x3 y 3 ) + µ(−10yz) + 0 + 0
= −(xy + x3 y 3 ) + µ(−10yz)
Finally
yz + 3x2 y 3 z
2 + 6xy
f = − xz + 3x3 y 2 z + µ 2 + 2x2 + 2y 2
xy + x3 y 3
−10yz
Note that there seems to be a sign problem with Eq.(26) in [7].
62
(65)
(66)
(67)
0.01000
velocity
pressure
x3
x2
error
0.00100
0.00010
0.00001
0.1
h
19.0.2
Variable viscosity
The spatial derivatives of the viscosity are then given by
∂
µ(x, y, z)
∂x
∂
µ(x, y, z)
∂y
∂
µ(x, y, z)
∂z
=
−(1 − 2x)βµ(x, y, z)
=
−(1 − 2y)βµ(x, y, z)
=
−(1 − 2z)βµ(x, y, z)
and thr right-hand side by
yz + 3x2 y 3 z
2 + 6xy
f = − xz + 3x3 y 2 z + µ 2 + 2x2 + 2y 2
(68)
xy + x3 y 3
−10yz
2˙xx
2˙xy
2˙xz
−(1 − 2x)βµ(x, y, z) 2˙xy − (1 − 2y)βµ(x, y, z) 2˙yy − (1 − 2z)βµ(x, y, z) 2˙yz
2˙xz
2˙yz
2˙zz
2 3
yz + 3x y z
2 + 6xy
= − xz + 3x3 y 2 z + µ 2 + 2x2 + 2y 2
xy + x3 y 3
−10yz
2
2 + 4x + 2y + 6x y
x + y + 2xy 2 + x3
−3z − 10xyz
−3z − 5x2 z
− (1 − 2x)βµ x + y + 2xy 2 + x3 − (1 − 2y)βµ 2 + 2x + 4y + 4x2 y − (1 − 2z)βµ
2
−3z − 5x z
−4 − 6x − 6y − 10x2 y
−3z − 10xyz
Note that at (x, y, z) = (0, 0, 0), µ = exp(1), and at (x, y, z) = (0.5, 0.5, 0.5), µ = exp(1 − 3β/4) so that the
63
maximum viscosity ratio is given by
µ? =
exp(1 − 3β/4)
= exp(−3β/4)
exp(1)
By varying β between 1 and 22 we can get up to 7 orders of magnitude viscosity difference.
features
• Q1 × P0 element
• incompressible flow
• saddle point system
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• 3D
• elemental b.c.
• analytical solution
64
20
fieldstone: solving the full saddle point problem with Q2 × Q1 elements
The details of the numerical setup are presented in Section 5.
Each element has mV = 9 vertices so in total ndofV × mV = 18 velocity dofs and ndofP ∗ mP = 4 pressure dofs.
The total number of velocity dofs is therefore N f emV = nnp × ndof V while the total number of pressure dofs is
N f emP = nel. The total number of dofs is then N f em = N f emV + N f emP .
As a consequence, matrix K has size N f emV, N f emV and matrix G has size N f emV, N f emP . Vector f is of size
N f emV and vector h is of size N f emP .
renumber all nodes to start at zero!! Also internal numbering does not work here
The velocity shape functions are given by:
N0V
N1V
N2V
N3V
N4V
N5V
N6V
N7V
N8V
1
1
r(r − 1) s(s − 1)
2
2
1
1
=
r(r + 1) s(s − 1)
2
2
1
1
=
r(r + 1) s(s + 1)
2
2
1
1
=
r(r − 1) s(s + 1)
2
2
1
= (1 − r2 ) s(s − 1)
2
1
=
r(r + 1)(1 − s2 )
2
1
= (1 − r2 ) s(s + 1)
2
1
=
r(r − 1)(1 − s2 )
2
= (1 − r2 )(1 − s2 )
=
65
and their derivatives:
∂N0V
∂r
∂N1V
∂r
∂N2V
∂r
∂N3V
∂r
∂N4V
∂r
∂N5V
∂r
∂N6V
∂r
∂N7V
∂r
∂N8V
∂r
∂N0V
∂s
∂N1V
∂s
∂N2V
∂s
∂N3V
∂s
∂N4V
∂s
∂N5V
∂s
∂N6V
∂s
∂N7V
∂s
∂N8V
∂s
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
1
1
(2r − 1) s(s − 1)
2
2
1
1
(2r + 1) s(s − 1)
2
2
1
1
(2r + 1) s(s + 1)
2
2
1
1
(2r − 1) s(s + 1)
2
2
1
(−2r) s(s − 1)
2
1
(2r + 1)(1 − s2 )
2
1
(−2r) s(s + 1)
2
1
(2r − 1)(1 − s2 )
2
(−2r)(1 − s2 )
1
1
r(r − 1) (2s − 1)
2
2
1
1
r(r + 1) (2s − 1)
2
2
1
1
r(r + 1) (2s + 1)
2
2
1
1
r(r − 1) (2s + 1)
2
2
2 1
(1 − r ) (2s − 1)
2
1
r(r + 1)(−2s)
2
1
(1 − r2 ) (2s + 1)
2
1
r(r − 1)(−2s)
2
(1 − r2 )(−2s)
features
• Q2 × Q1 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (no-slip)
• isothermal
• isoviscous
• analytical solution
66
0.0100000
velocity
pressure
x3
x2
0.0010000
error
0.0001000
0.0000100
0.0000010
0.0000001
0.1
h
67
21
fieldstone: solving the full saddle point problem with Q3 × Q2 elements
The details of the numerical setup are presented in Section 5.
Each element has mV = 16 vertices so in total ndofV × mV = 32 velocity dofs and ndofP ∗ mP = 9 pressure dofs.
The total number of velocity dofs is therefore N f emV = nnp × ndof V while the total number of pressure dofs is
N f emP = nel. The total number of dofs is then N f em = N f emV + N f emP .
As a consequence, matrix K has size N f emV, N f emV and matrix G has size N f emV, N f emP . Vector f is of size
N f emV and vector h is of size N f emP .
60===61===62===63===64===65===66===67===68===70
||
||
||
||
50
51
52
53
54
55
56
57
58
59
||
||
||
||
40
41
42
43
44
45
46
47
48
49
||
||
||
||
30===31===32===33===34===35===36===37===38===39
||
||
||
||
20
21
22
23
24
25
26
27
28
29
||
||
||
||
10
11
12
13
14
15
16
17
18
19
||
||
||
||
00===01===02===03===04===05===06===07===08===09
Example of 3x2 mesh. nnx=10, nny=7, nnp=70, nelx=3, nely=2, nel=6
12===13===14===15
||
||
||
||
08===09===10===11
||
||
||
||
04===05===06===07
||
||
||
||
00===01===02===03
06=====07=====08
||
||
||
||
||
||
03=====04=====05
||
||
||
||
||
||
00=====01=====02
Velocity (Q3)
Pressure (Q2)
(r,s)_{00}=(-1,-1)
(r,s)_{01}=(-1/3,-1)
(r,s)_{02}=(+1/3,-1)
(r,s)_{03}=(+1,-1)
(r,s)_{04}=(-1,-1/3)
(r,s)_{05}=(-1/3,-1/3)
(r,s)_{06}=(+1/3,-1/3)
(r,s)_{07}=(+1,-1/3)
(r,s)_{08}=(-1,+1/3)
(r,s)_{09}=(-1/3,+1/3)
(r,s)_{10}=(+1/3,+1/3)
(r,s)_{11}=(+1,+1/3)
(r,s)_{12}=(-1,+1)
(r,s)_{13}=(-1/3,+1)
(r,s)_{14}=(+1/3,+1)
(r,s)_{15}=(+1,+1)
(r,s)_{00}=(-1,-1)
(r,s)_{01}=(0,-1)
(r,s)_{02}=(+1,-1)
(r,s)_{03}=(-1,0)
(r,s)_{04}=(0,0)
(r,s)_{05}=(+1,0)
(r,s)_{06}=(-1,+1)
(r,s)_{07}=(0,+1)
(r,s)_{08}=(+1,+1)
The velocity shape functions are given by:
N1 (r)
=
(−1 + r + 9r2 − 9r3 )/16
N2 (r)
=
(+9 − 27r − 9r2 + 27r3 )/16
N3 (r)
=
(+9 + 27r − 9r2 − 27r3 )/16
N4 (r)
=
(−1 − r + 9r2 + 9r3 )/16
68
N1 (t)
=
(−1 + t + 9t2 − 9t3 )/16
N2 (t)
=
(+9 − 27t − 9t2 + 27t3 )/16
N3 (t)
=
(+9 + 27t − 9t2 − 27t3 )/16
N4 (t)
=
(−1 − t + 9t2 + 9t3 )/16
N0 1(r, t)
= N1 (r)N1 (t)
(69)
N0 2(r, t)
=
N2 (r)N1 (t)
(70)
N0 3(r, t)
=
N3 (r)N1 (t)
(71)
N0 4(r, t)
=
N4 (r)N1 (t)
(72)
N0 5(r, t)
=
N1 (r)N2 (t)
(73)
N0 6(r, t)
=
N2 (r)N2 (t)
(74)
N0 7(r, t)
=
N3 (r)N2 (t)
(75)
N0 8(r, t)
=
N4 (r)N2 (t)
(76)
N0 9(r, t)
=
N1 (r)N3 (t)
(77)
N1 0(r, t)
=
N2 (r)N3 (t)
(78)
N1 1(r, t)
=
N3 (r)N3 (t)
(79)
N1 2(r, t)
=
N4 (r)N3 (t)
(80)
N1 3(r, t)
=
N1 (r)N4 (t)
(81)
N1 4(r, t)
=
N2 (r)N4 (t)
(82)
N1 5(r, t)
=
N3 (r)N4 (t)
(83)
N1 6(r, t)
=
N4 (r)N4 (t)
(84)
and their derivatives:
Write about 4 point quadrature.
features
• Q3 × Q2 element
• incompressible flow
• mixed formulation
• isothermal
• isoviscous
• analytical solution
0.0001
velocity
pressure
x4
x3
1x10-5
error
1x10-6
1x10-7
1x10-8
1x10-9
1x10-10
0.1
h
velocity error rate is cubic, pressure superconvergent since the pressure field is quadratic and therefore lies into the
Q2 space.
69
22
fieldstone: the Busse benchmark
This three-dimensional benchmark was first proposed by [8]. It has been subsequently presented in [65, 72, 1, 51, 14, 42].
We here focus on Case 1 of [8]: an isoviscous bimodal convection experiment at Ra = 3 × 105 .
The domain is of size a × b × h with a = 1.0079h, b = 0.6283h with h = 2700km. It is filled with a Newtonian fluid
characterised by ρ0 = 3300kg.m−3 , α = 10−5 K−1 , µ = 8.0198×1023 Pa.s, k = 3.564W.m−1 .K−1 , cp = 1080J.K−1 .kg−1 .
The gravity vector is set to g = (0, 0, −10)T . The temperature is imposed at the bottom (T = 3700◦ C) and at the
top (T = 0◦ C).
The various measurements presented in [8] are listed hereafter:
• The Nusselt number N u computed at the top surface following Eq. (38):
RR
∂T
dxdy
z=L ∂y
N u = Lz R R z
T
dxdy
z=0
• the root mean square velocity vrms and the temperature mean square velocity Trms
• The vertical velocity w and temperature T at points x1 = (0, 0, Lz /2), x2 = (Lx , 0, Lz /2), x3 = (0, Ly , Lz /2)
and x4 = (Lx , Ly , Lz /2);
• the vertical component of the heat flux Q at the top surface at all four corners.
The values plotted hereunder are adimensionalised by means of a reference temperature (3700K), a reference
lengthscale 2700km, and a reference time L2z /κ.
0.0009
2124
0.0008
2123
2122
0.0007
2121
vrms
0.0006
0.0005
0.0004
2119
2118
0.0003
2117
0.0002
0.0001
2120
2116
0
2115
2x108 4x108 6x108 8x108 1x109 1.2x1091.4x1091.6x109
0
2x108
4x108
time/year
0
0.0015
-0.5
0.001
1x109 1.2x109 1.4x109 1.6x109
-1
0.0005
heat flux
w at Lz/2
8x108
time/year
0.002
0
-0.0005
-1.5
-2
-2.5
-0.001
-3
-0.0015
-0.002
6x108
0
-3.5
1x108
2x108 4x108 6x108 8x108 1x109 1.2x1091.4x1091.6x109
time/year
1.5x108
2x108
2.5x108
time/year
70
3x108
3.5x108
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• 3D
• elemental b.c.
• buoyancy driven
ToDo: look at energy conservation. run to steady state and make sure the expected values are retrieved.
71
fieldstone: The non-conforming Q1 × P0 element
23
features
• Non-conforming Q1 × P0 element
• incompressible flow
• mixed formulation
• isothermal
• non-isoviscous
• analytical solution
• pressure smoothing
try Q1 mapping instead of isoparametric.
0.005 0.8
0.6
0.4
0.6
0.000
0.0
0.0
0.4
0.0050.2
0.2
0.2
0.4
x
0.6
0.8
1.0
0.05
0.00
0.005 0.8
y
y
0.000
p
1.0
0.05
0.6
0.10
y
0.8
0.010
vy
1.0
0.0
0.010 0.0
0.4
0.15
0.0050.2
0.2
0.4
x
0.6
0.8
1.0
0.00
q
1.0
0.05
0.8
0.10
0.6
y
0.010
vx
1.0
0.0
0.010 0.0
0.20
0.2
0.4
x
0.6
0.8
1.0
0.15
0.4
0.20
0.2
0.25 0.00.0
0.2
0.4
x
0.6
0.8
1.0
0.06
0.01
0.6
0.02
0.8
0.01
0.6
0.00
y
y
0.00
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.4
x
0.6
0.8
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.0000010
1.0
vx txth
1.0
0.03
1.0
0.4
x
0.6
0.8
1.0
0.03
0.0
0.0
0.0000005
0.2
0.4
x
0.6
0.8
1.0
0.0
0.0
0.0000010
0.2
0.4
x
0.6
0.8
1.0
0.4
x
0.6
0.8
1.0
0.04
0.6
0.03
0.0
0.0
0.0000010
72
0.4
0.02
0.2
0.04 0.0
0.0
0.2
0.00
p pth
0.02
0.04
0.06
0.08
0.0000005
0.2
0.2
0.05
0.8
0.4
0.06
y
0.6
0.0000000
0.4
y
0.6
0.0000000
0.4
y
0.6
0.2
0.02
0.2
0.8
0.0000005
0.0
0.0
0.00
0.4
0.8
0.0000005
0.2
0.02
0.6
0.8
0.4
0.04 1.0
0.8
0.0000010
1.0
vy tyth
0.06
y
0.02
0.8
xy
x
0.6
0.8
1.0
q pth
1.0
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
0.02
0.8
0.04
0.6
0.06
0.4
0.08
0.2
0.0
0.0
0.01
0.00
y
1.0
0.03
yy
y
0.03
xx
0.25
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
24
fieldstone: The stabilised Q1 × Q1 element
The details of the numerical setup are presented in Section 5.
73
25
fieldstone: compressible flow (1)
We first start with an isothermal Stokes flow, so that we disregard the heat transport equation and the equations we
wish to solve are simply:
1
−∇ · 2η ε̇(v) − (∇ · v)1 + ∇p = ρg
3
∇ · (ρv) = 0
in Ω,
(85)
in Ω
(86)
The second equation can be rewritten ∇ · (ρv) = ρ∇ · v + v · ∇ρ = 0 or,
1
∇ · v + v · ∇ρ = 0
ρ
Note that this presupposes that the density is not zero anywhere in the domain.
We use a mixed formulation and therefore keep both velocity and pressure as unknowns. We end up having to
solve the following system:
K
G
V
f
·
=
or,
A · X = rhs
P
h
GT + Z 0
Where K is the stiffness matrix, G is the discrete gradient operator, GT is the discrete divergence operator, V the
velocity vector, P the pressure vector. Note that the term ZV derives from term v · ∇ρ in the continuity equation.
Each block K, G , Z and vectors f and h are built separately in the code and assembled into the matrix A and
vector rhs afterwards. A and rhs are then passed to the solver. We will see later that there are alternatives to solve
this approach which do not require to build the full Stokes matrix A.
Remark: the term ZV is often put in the rhs (i.e. added to h) so that the matrix A retains the same structure as
in the incompressible case. This is indeed how it is implemented in ASPECT. This however requires more work since
the rhs depends on the solution and some form of iterations is needed.
In the case of a compressible flow the strain rate tensor and the deviatoric strain rate tensor are no more equal
(since ∇ · v 6= 0). The deviatoric strainrate tensor is given by6
1
1
˙ = (v)
˙
˙
− (∇ · v)1
˙ d (v) = (v)
− T r()1
3
3
In that case:
˙dxx
=
˙dyy
=
2˙dxy
=
∂u 1 ∂u ∂v
2 ∂u 1 ∂v
−
+
−
=
∂x 3 ∂x ∂y
3 ∂x 3 ∂y
1 ∂u ∂v
1 ∂u 2 ∂v
∂v
−
+
=−
+
∂y 3 ∂x ∂y
3 ∂x 3 ∂y
∂u ∂v
+
∂y
∂x
and then
˙ d (v) =
2 ∂u
3 ∂x
−
1 ∂v
3 ∂y
1 ∂u
2 ∂y
+
1 ∂v
2 ∂x
1 ∂u
2 ∂y
+
1 ∂v
2 ∂x
− 13 ∂u
∂x +
2 ∂v
3 ∂y
(87)
(88)
(89)
From ~τ = 2η~d we arrive at:
d
˙xx
τxx
2/3
τyy = 2η ˙dyy = 2η −1/3
τxy
0
˙dxy
−1/3
2/3
0
0
0 ·
1/2
∂u
∂x
∂v
∂y
∂u
∂v
∂y + ∂x
4/3
= η −2/3
0
or,
~τ = Cη BV
6 See
the ASPECT manual for a justification of the 3 value in the denominator in 2D and 3D.
74
−2/3
4/3
0
0
0 ·
1
∂u
∂x
∂v
∂y
∂u
∂v
∂y + ∂x
In order to test our implementation we have created a few manufactured solutions:
• benchmark #1 (ibench=1)): Starting from a density profile of:
ρ(x, y) = xy
(90)
We derive a velocity given by:
vx (x, y) =
With gx (x, y) =
1
x
Cx
Cy
, vy (x, y) =
x
y
and gy (x, y) = y1 , this leads us to a pressure profile:
4Cx
4Cy
p = −η
+ 2 + xy + C0
3x2
3y
(91)
(92)
This gives us a strain rate:
−Cx
−Cy
˙yy = 2
˙xy = 0
x2
y
In what follows, we choose η = 1 and Cx = Cy = 1 and for a unit square domain [1 : 2] × [1 : 2] we compute C0
so that the pressure is normalised to zero over the whole domain and obtain C0 = −1.
˙xx =
• benchmark #2 (ibench=2): Starting from a density profile of:
ρ = cos(x) cos(y)
(93)
We derive a velocity given by:
vx =
With gx =
1
cos(y)
and gy =
Cx
Cy
, vy =
cos(x)
cos(y)
(94)
1
cos(x) ,
this leads us to a pressure profile:
!
4Cx sin(x) 4Cy sin(y)
p=η
+
+ (sin(x) + sin(y)) + C0
3 cos2 (x)
3 cos2 (y)
(95)
sin(y)
sin(x)
˙yy = Cy
˙xy = 0
cos2 (x)
cos2 (y)
We choose η = 1 and Cx = Cy = 1. The domain is the unit square [0 : 1] × [0 : 1] and we obtain C0 as before
and obtain
1
C0 = 2 − 2 cos(1) + 8/3(
− 1) ' 3.18823730
cos(1)
(thank you WolframAlpha)
˙xx = Cx
• benchmark #3 (ibench=3)
• benchmark #4 (ibench=4)
• benchmark #5 (ibench=5)
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (no-slip)
• isothermal
• isoviscous
• analytical solution
• pressure smoothing
ToDo:
• pbs with odd vs even number of elements
• q is ’fine’ everywhere except in the corners - revisit pressure smoothing paper?
• redo A v d Berg benchmark (see Tom Weir thesis)
75
26
26.1
fieldstone: compressible flow (2)
The physics
Let us start with some thermodynamics. Every material has an equation of state. The equilibrium thermodynamic
state of any material can be constrained if any two state variables are specified. Examples of state variables include
the pressure p and specific volume ν = 1/ρ, as well as the temperature T .
After linearisation, the density depends on temperature and pressure as follows:
ρ(T, p) = ρ0 ((1 − α(T − T0 ) + βT p)
where α is the coefficient of thermal expansion, also called thermal expansivity:
1 ∂ρ
α=−
ρ ∂T p
α is the percentage increase in volume of a material per degree of temperature increase; the subscript p means that
the pressure is held fixed.
βT is the isothermal compressibility of the fluid, which is given by
1
1 ∂ρ
βT =
=
K
ρ ∂P T
with K the bulk modulus. Values of βT = 10−12 −10−11 Pa−1 are reasonable for Earth’s mantle, with values decreasing
by about a factor of 5 between the shallow lithosphere and core-mantle boundary. This is the percentage increase
in density per unit change in pressure at constant temperature. Both the coefficient of thermal expansion and the
isothermal compressibility can be obtained from the equation of state.
The full set of equations we wish to solve is given by
−∇ · 2η ˙ d (v) + ∇p = ρ0 ((1 − α(T − T0 ) + βT p) g
in Ω
1
∇ · v + v · ∇ρ = 0
in Ω
ρ
∂T
∂p
d
d
ρCp
+ v · ∇T − ∇ · k∇T = ρH + 2η ˙ : ˙ + αT
+ v · ∇p
∂t
∂t
(96)
(97)
in Ω,
(98)
Note that this presupposes that the density is not zero anywhere in the domain.
26.2
The numerics
We use a mixed formulation and therefore keep both velocity and pressure as unknowns. We end up having to solve
the following system:
K
G+W
V
f
·
=
or,
A · X = rhs
P
h
GT + Z
0
Where K is the stiffness matrix, G is the discrete gradient operator, GT is the discrete divergence operator, V the
velocity vector, P the pressure vector. Note that the term ZV derives from term v · ∇ρ in the continuity equation.
As perfectly explained in the step 32 of deal.ii7 , we need to scale the G term since it is many orders of magnitude
smaller than K, which introduces large inaccuracies in the solving process to the point that the solution is nonsensical.
This scaling coefficient is η/L. After building the G block, it is then scaled as follows: G0 = Lη G so that we now solve
K
G0 + W
V
f
·
=
P0
h
G0T + Z
0
After the solve phase, we recover the real pressure with P = Lη P 0 .
adapt notes since I should scale W and Z too. h should be caled too !!!!!!!!!!!!!!!
Each block K, G , Z and vectors f and h are built separately in the code and assembled into the matrix A and
vector rhs afterwards. A and rhs are then passed to the solver. We will see later that there are alternatives to solve
this approach which do not require to build the full Stokes matrix A.
Remark 1: the terms ZV and WP are often put in the rhs (i.e. added to h) so that the matrix A retains the
same structure as in the incompressible case. This is indeed how it is implemented in ASPECT, see also appendix A
of [44]. This however requires more work since the rhs depends on the solution and some form of iterations is needed.
7 https://www.dealii.org/9.0.0/doxygen/deal.II/step
32.html
76
Remark 2: Very often the adiabatic heating term αT (v · ∇p) is simplified as follows: If you assume the vertical
component of the gradient of the dynamic pressure to be small compared to the gradient of the total pressure (in
other words, the gradient is dominated by the gradient of the hydrostatic pressure), then −ρg ' ∇p and then
αT (v · ∇p) ' −αρT v · g. We will however not be using this approximation in what follows.
We have already established that
~τ = Cη BV
The following measurements are carried out:
• The root mean square velocity (vrms):
s
vrms =
1
V
Z
1
V
Z
v 2 dV
V
• The average temperature (Tavrg):
< T >=
T dV
V
• The total mass (mass):
Z
M=
ρdV
V
• The Nusselt number (Nu):
Nu = −
1 1
Lx ∆T
Lx
Z
0
∂T (x, y = Ly )
dx
∂y
• The kinetic energy (EK):
Z
EK =
V
1 2
ρv dV
2
• The work done against gravity
Z
< W >= −
ρgy vy dV
V
• The total viscous dissipation (visc diss)
Z
< Φ >=
1
ΦdV =
V
Z
2η ε̇ : ε̇dV
• The gravitational potential energy (EG)
Z
ρgy (Ly − y)dV
EG =
V
• The internal thermal energy (ET)
Z
ET =
ρ(0) Cp T dV
V
Remark 3: Measuring the total mass can be misleading: indeed because ρ = ρ0 (1 − αT ), then measuring the
total mass amounts to measuring a constant minus the volume-integrated temperature, and there is no reason
why the latter should be zero, so that there is no reason why the total mass should be zero...!
26.3
The experimental setup
The setup is as follows: the domain is Lx = Ly = 3000km. Free slip boundary conditions are imposed on all four
sides. The initial temperature is given by:
Ly − y
πx
πy
T (x, y) =
− 0.01 cos( ) sin( ) ∆T + Tsurf
Ly
Lx
Ly
with ∆T = 4000K, Tsurf = T0 = 273.15K. The temperature is set to ∆T + Tsurf at the bottom and Tsurf at the top.
We also set k = 3, Cp = 1250, |g| = 10, ρ0 = 3000 and we keep the Rayleigh number Ra and dissipation number Di
as input parameters:
αg∆T L3 ρ20 Cp
αgL
Ra =
Di =
ηk
Cp
77
DiCp
gL ,
From the second equation we get α =
Ra =
which we can insert in the first one:
DiCp2 ∆T L2 ρ20
ηk
or,
DiCp2 ∆T L2 ρ20
Ra k
η=
For instance, for Ra = 104 and Di = 0.75, we obtain α ' 3 · 10−5 and η ' 1025 which are quite reasonable values.
26.4
Scaling
Following [40], we non-dimensionalize the equations using the reference values for density ρr , thermal expansivity
αr , temperature contrast ∆Tr (refTemp), thermal conductivity kr , heat capacity Cp , depth of the fluid layer L and
viscosity ηr . The non-dimensionalization for velocity, ur , pressure pr and time, tr become
ur =
kr
ρr Cp L
(refvel)
pr =
ηr kr
ρr Cp L2
(refpress)
tr =
ρr Cp L2
kr
(reftime)
In the case of the setup described hereabove, and when choosing Ra = 104 and Di = 0.5, we get:
alphaT 2.083333e-05
eta 8.437500e+24
reftime 1.125000e+19
refvel 2.666667e-13
refPress 7.500000e+05
26.5
Conservation of energy 1
26.5.1
under BA and EBA approximations
Following [44], we take the dot product of the momentum equation with the velocity v and integrate over the whole
volume8 :
Z
Z
[−∇ · τ + ∇p] · vdV =
ρg · vdV
V
V
or,
Z
−
Z
Z
(∇ · τ ) · vdV +
V
∇p · vdV =
V
ρg · vdV
V
Let us look at each block separately:
Z
Z
Z
Z
Z
− (∇ · τ ) · vdV = − τ |v{z
· n} dS +
τ : ∇vdV =
τ : ε̇dV =
ΦdV
V
S
V
V
V
=0 (b.c.)
which is the volume integral of the shear heating. Then,
Z
Z
Z
∇p · vdV =
p v
· n} dS −
| {z
V
S
V
=0 (b.c.)
∇
| {z· v}
pdV = 0
=0 (incomp.)
which is then zero in the case of an incompressible flow. And finally
Z
ρg · vdV = W
V
which is the work against gravity.
Conclusion for an incompressible fluid: we should have
Z
Z
ΦdV =
ρg · vdV
V
8 Check:
V
this is akin to looking at the power, force*velocity, says Arie
78
(99)
This formula is hugely problematic: indeed, the term ρ in the rhs is the full density. We know that to the value of ρ0
corresponds a lithostatic pressure gradient pL = ρ0 gy. In this case one can write ρ = ρ0 + ρ0 and p = pL + p0 so that
we also have
Z
Z
[−∇ · τ + ∇p0 ] · vdV =
ρ0 g · vdV
V
V
which will ultimately yield
Z
Z
ΦdV =
V
Z
ρ0 g · vdV =
(ρ − ρ0 )g · vdV
V
(100)
V
Obviously Eqs.(99) and (100) cannot be true at the same time. The problem comes from the nature of the (E)BA
approximation: ρ = ρ0 in the mass conservation equation but it is not constant in the momentum conservation
equation, which
R is of course inconsistent. Since the mass conservation equation is ∇ · v = 0 under this approximation
then the term V ∇p · vdV is always zero for any pressure (full pressure p, or overpressure p − pL ), hence the paradox.
This paradox will be lifted when a consistent set of equations will be used (compressible formulation). On a practical
note, Eqs.(99) is not verified by the code, while (100) is.
In the end:
Z
Z
(ρ − ρ0 )g · vdV
ΦdV =
V
{z
}
| {z } |
(101)
V
visc diss
26.5.2
work grav
under no approximation at all
Z
Z
∇p · vdV
=
V
Z
p v
· n} dS −
| {z
S
=0 (b.c.)
Z
=
V
∇ · v pdV = 0
(102)
V
1
v · ∇ρ pdV = 0
ρ
(103)
(104)
ToDo:see section 3 of [44] where this is carried out with the Adams-Williamson eos.
26.6
Conservation of energy 2
Also, following the Reynold’s transport theorem [48],p210, we have for a property A (per unit mass)
Z
Z
Z
d
∂
AρdV =
(Aρ)dV +
Aρv · ndS
dt V
V ∂t
S
Let us apply to this to A = Cp T and compute the time derivative of the internal energy:
Z
Z
Z
Z
Z
d
∂
∂ρ
∂T
ρCp T dV =
(ρCp T )dV +
Aρ v
·
n
dS
=
C
T
dV
+
ρCp
dV
p
{z
}
|
dt V
∂t
∂t
∂t
V
S
V
V
{z
} |
{z
}
|
=0 (b.c.)
I
(105)
II
In order to expand I, the mass conservation equation will be used, while the heat transport equation will be used
for II:
Z
I=
V
Z
II =
V
∂ρ
Cp T dV
∂t
∂T
ρCp
dV
∂t
Z
= −
Z
Cp T ∇ · (ρv)dV = −
V
Z
=
−ρCp v · ∇T
V
Z
=
−ρCp v · ∇T
V
Z
=
−ρCp v · ∇T
V
Z
=
−ρCp v · ∇T
V
V
Z
Cp T ρ v
· n} dS +
| {z
=0 (b.c.)
ρCp ∇T · vdV
(106)
V
∂p
+ ∇ · k∇T + ρH + Φ + αT
+ v · ∇p dV
∂t
Z
∂p
+ ρH + Φ + αT
+ v · ∇p dV +
∇ · k∇T dV
∂t
V
Z
∂p
+ ρH + Φ + αT
+ v · ∇p dV +
k∇T · ndS
∂t
S
Z
∂p
+ ρH + Φ + αT
+ v · ∇p dV −
q · ndS
∂t
S
Finally:
79
(107)
(108)
(109)
(110)
d
I + II =
dt
Z
|V
ρCp T dV
{z
}
Z
Z
∂p
=
+ v · ∇p dV −
ρH + Φ + αT
q · ndS
∂t
V
S
(111)
ET
Z
=
Z
Z
∂p
αT v · ∇pdV −
αT dV +
ΦdV +
∂t
V
V
V
{z
} |
{z
}
| {z } |
Z
ρHdV +
V
visc diss
extra
adiab heating
Z
q · ndS
| {z }
(112)
S
heatflux boundary
This was of course needlessly complicated as the term ∂ρ/∂t is always taken to be zero, so that I = 0 automatically.
The mass conservation equation is then simply ∇ · (ρv) = 0. Then it follows that
Z
Z
Z
ρCp ∇T · vdV
(113)
0 =
Cp T ∇ · (ρv)dV = −
Cp T ρ v
· n} dS +
| {z
V
V
=0 (b.c.)
V
Z
ρCp ∇T · vdV
=
(114)
V
so that the same term in Eq.(110) vanishes too, and then Eq.(112) is always valid, although one should be careful
when computing ET in the BA and EBA cases as it should use ρ0 and not ρ.
26.7
The problem of the onset of convection
[wiki] In geophysics, the Rayleigh number is of fundamental importance: it indicates the presence and strength of
convection within a fluid body such as the Earth’s mantle. The mantle is a solid that behaves as a fluid over geological
time scales.
The Rayleigh number essentially is an indicator of the type of heat transport mechanism. At low Rayleigh numbers
conduction processes dominate over convection ones. At high Rayleigh numbers it is the other way around. There is
a so-called critical value of the number with delineates the transition from one regime to the other.
This problem has been studied and approached both theoretically and numerically [73, e.g.] and it was found that
the critical Rayleigh number Rac is
Rac = (27/4)π 4 ' 657.5
in setups similar to ours.
VERY BIG PROBLEM
The temperature setup is built as follows: Tsurf is prescribed at the top, Tsurf + ∆T is prescribed at the bottom.
The initial temperature profile is linear between these two values. In the case of BA, the actual value of Tsurf is of no
consequence. However, for the EBA the full temperature is present in the adiabatic heating term on the rhs of the hte,
and the value of Tsurf will therefore influence the solution greatly. This is very problematic as there is no real way to
arrive at the surface temperature from the King paper. On top of this, the density uses a reference temperature T0
which too will influence the solution without being present in the controlling Ra and Di numbers!!
In light thereof, it will be very difficult to recover the values of King et al for EBA!
80
features
• Q1 × P0 element
• compressible flow
• mixed formulation
• Dirichlet boundary conditions (no-slip)
• isoviscous
• analytical solution
• pressure smoothing
Relevant literature: [4, 37, 66, 44, 40, 45, 46, 33]
ToDo:
• heat flux is at the moment elemental, so Nusselt and heat flux on boundaries measurements not as accurate as
could be.
• implement steady state detection
• do Ra = 105 and Ra = 106
• velocity average at surface
• non dimensional heat flux at corners [5]
• depth-dependent viscosity (case 2 of [5])
81
26.8
results - BA - Ra = 104
These results were obtained with a 64x64 resolution, and CFL number of 1. Steady state was reached after about
1250 timesteps.
7x10-6
-7.73687x1023
7.67188x1022
6x10-6
-7.73687x1023
7.67188x1022
5x10-6
-7.73687x1023
7.67188x1022
4x10-6
-7.73687x1023
3x10-6
EG
ET
EK
7.67188x1022
7.67188x1022
-7.73687x1023
7.67188x1022
1x10-6
0
(a)
-7.73687x1023
7.67188x1022
0
10000 20000 30000 40000 50000 60000 70000 80000 90000100000
time (Myr)
7.67188x1022
(b)
-7.73687x1023
0 100002000030000400005000060000700008000090000
100000
time (Myr)
2273.15
-7.73687x1023
-7.73687x1023
7.67188x1022
2x10-6
-7.73687x1023
4500
-7.73687x1023
(c)
time (Myr)
0.0015
min(T)
max(T)
4000
0 100002000030000400005000060000700008000090000
100000
2273.15
min(u)
max(u)
min(v)
max(v)
0.001
3500
2273.15
2273.15
3000
0.0005
velocity
Temperature
2273.15
2500
2000
1500
0
-0.0005
1000
2273.15
-0.001
500
2273.15
(d)
0
100002000030000400005000060000700008000090000100000
time (Myr)
0
(e)
time (Myr)
-10000
80000
80000
70000
70000
-40000
-50000
-60000
60000
50000
40000
30000
20000
-70000
10000
-80000
0
(g)
work against gravity
90000
-30000
0
10000 20000 30000 40000 50000 60000 70000 80000 90000100000
time (Myr)
(h)
60000
50000
40000
30000
20000
10000
0
10000 20000 30000 40000 50000 60000 70000 80000 90000100000
time (Myr)
0
(i)
0.07
8
2x10-8
0.06
7
-2x10-8
(j)
0
10000 20000 30000 40000 50000 60000 70000 80000 90000100000
time (Myr)
8x10-7
0
(k)
2
1
0
10000 20000 30000 40000 50000 60000 70000 80000 90000100000
time (Myr)
90000
0
HF
dETdt
6x10-7
4
3
0.01
-4x10-8
-5x10-8
5
0.03
0.02
-3x10-8
4.88
0.04
Nu
vrms (cm/yr)
-1x10-8
10000 20000 30000 40000 50000 60000 70000 80000 90000100000
6
0.05
0
0
time (Myr)
3x10-8
1x10-8
0 100002000030000400005000060000700008000090000
100000
time (Myr)
90000
viscous dissipation
adiabatic heating
(f)
0
-20000
heat flux
0 100002000030000400005000060000700008000090000
100000
-0.0015
0
(l)
0
10000 20000 30000 40000 50000 60000 70000 80000 90000100000
time (Myr)
VD
WG
80000
70000
4x10-7
60000
2x10-7
50000
0
40000
30000
-2x10-7
20000
-4x10-7
-6x10-7
(l)
10000
0
0 1x1010
2x1010
3x1010
4x1010
5x1010
6x1010
7x1010
8x1010
9x1010
1x1011
0
1x10102x10103x10104x10105x10106x10107x10108x10109x10101x1011
time (Myr)
(m)
AH: adiabatic heating, VD: viscous dissipation, HF: heat flux, WG: work against gravity
time (Myr)
Eq.(112) is verified by (l) and Eq.(101) is verified by (m).
82
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
83
26.9
results - BA - Ra = 105
These results were obtained with a 64x64 resolution, and CFL number of 1. Steady state was reached after about
1250 timesteps.
0.0007
7.67188x1022
0.0006
7.67188x1022
0.0005
-7.73687x1023
-7.73687x1023
-7.73687x1023
7.67188x1022
-7.73687x1023
7.67188x1022
EG
ET
EK
0.0004
0.0003
7.67188x1022
0.0002
0
(a)
5000
10000
15000
20000
25000
30000
35000
time (Myr)
7.67188x1022
(b)
2273.150000001
4500
2273.150000001
4000
Temperature
-7.73687x1023
(c)
2273.150000000
2273.150000000
2273.150000000
2000
1500
1000
500
2273.149999999
0
min(u)
max(u)
min(v)
max(v)
0
-0.005
-0.01
0
5000
10000 15000 20000 25000 30000 35000
time (Myr)
(e)
0
5000
-0.015
10000 15000 20000 25000 30000 35000
(f)
time (Myr)
0
-100000
800000
800000
-200000
700000
700000
-400000
-500000
-600000
600000
500000
400000
300000
-700000
200000
-800000
100000
-900000
0
5000
10000
15000
20000
25000
30000
35000
time (Myr)
(h)
0.8
2x10-6
0.7
vrms (cm/yr)
1x10-6
5x10-7
0
10000
15000
20000
25000
30000
35000
0
5000
10000
15000
20000
25000
30000
35000
1x10-5
(k)
5000
10000
15000
20000
25000
30000
35000
18
10.53
14
12
0.4
0.3
10
8
6
4
2
0
5000
0
10000 15000 20000 25000 30000 35000
(l)
time (Myr)
900000
700000
600000
2x10-6
500000
0
400000
-2x10-6
300000
-4x10-6
200000
-6x10-6
100000
0
5x109 1x1010 1.5x1010 2x1010 2.5x1010 3x1010 3.5x1010
0
5000
10000
15000
20000
25000
30000
time (Myr)
VD
WG
800000
4x10-6
0
0
16
6x10-6
-8x10-6
300000
time (Myr)
193.2150
0
HF
dETdt
8x10-6
400000
0
0.5
0
500000
(i)
0.1
time (Myr)
(l)
5000
0.2
-5x10-7
600000
100000
0
0.6
1.5x10-6
10000 15000 20000 25000 30000 35000
200000
time (Myr)
2.5x10-6
(j)
work against gravity
900000
viscous dissipation
900000
-300000
5000
time (Myr)
0
-1x10-6
10000 15000 20000 25000 30000 35000
0.005
2500
2273.149999999
(g)
5000
time (Myr)
3000
2273.149999999
0
0
0.015
min(T)
max(T)
Nu
10000 15000 20000 25000 30000 35000
3500
2273.150000000
adiabatic heating
5000
0.01
2273.150000000
(d)
-7.73687x1023
0
time (Myr)
2273.150000001
heat flux
-7.73687x1023
velocity
0
-7.73687x1023
-7.73687x1023
7.67188x1022
0.0001
-7.73687x1023
0
5x109
1x1010 1.5x1010 2x1010 2.5x1010 3x1010 3.5x1010
time (Myr)
(m)
AH: adiabatic heating, VD: viscous dissipation, HF: heat flux, WG: work against gravity
time (Myr)
Eq.(112) is verified by (l) and Eq.(101) is verified by (m).
84
35000
26.10
results - BA - Ra = 106
85
results - EBA - Ra = 104
26.11
These results were obtained with a 64x64 resolution, and CFL number of 1. Steady state was reached after about
2500 timesteps
8x10-7
7.7x1022
7.6x1022
6x10-7
-7.74x1023
7.55x1022
5x10-7
7.5x1022
-7.745x1023
7.45x1022
EG
4x10-7
ET
EK
-7.735x1023
7.65x1022
7x10-7
7.4x1022
3x10-7
-7.75x1023
7.35x1022
2x10-7
7.3x1022
1x10-7
-7.755x1023
7.25x1022
0
0
(a)
50000
100000
150000
200000
250000
time (Myr)
7.2x1022
(b)
0
50000
100000
150000
200000
-7.76x1023
250000
(c)
time (Myr)
2280
4500
2200
2180
2140
(d)
50000
100000
150000
200000
250000
2000
(e)
-0.0003
0
50000
100000
150000
200000
-0.0004
250000
(f)
10000
10000
9000
9000
8000
8000
-3000
7000
-5000
-6000
-7000
6000
5000
4000
3000
-8000
2000
-9000
1000
-10000
0
0
50000
100000
150000
200000
250000
time (Myr)
work against gravity
0
-2000
-4000
(h)
50000
100000
150000
200000
250000
200000
250000
6000
5000
4000
3000
0
0
50000
100000
150000
200000
250000
time (Myr)
2.4
2.2
0.02
2
7000
5000
4000
3000
1.8
0.015
Nu
vrms (cm/yr)
6000
0.01
1.2
0.005
1000
1.6
1.4
2000
1
0
0
50000
100000
150000
200000
250000
time (Myr)
0
(k)
1000
0
50000
100000
-1000
200000
0.8
250000
(l)
0
50000
100000
150000
200000
time (Myr)
10000
VD
WG
9000
8000
-2000
7000
-3000
6000
-4000
5000
-5000
4000
-6000
3000
-7000
2000
-8000
1000
-9000
-10000
150000
time (Myr)
0
AH+VD+HF
dETdt
0
(l)
150000
7000
(i)
0.025
8000
(j)
100000
1000
0
9000
-1000
50000
2000
time (Myr)
10000
0
time (Myr)
-1000
(g)
heat flux
0
-0.0002
time (Myr)
viscous dissipation
adiabatic heating
time (Myr)
min(u)
max(u)
min(v)
max(v)
-0.0001
1500
0
250000
0.0001
2500
500
0
200000
0.0002
3000
1000
2160
150000
0.0003
velocity
Temperature
2220
100000
0.0004
3500
2240
50000
time (Myr)
min(T)
max(T)
4000
2260
0
0
50000
100000
150000
200000
0
250000
0
50000
100000
150000
200000
250000
time (Myr)
(m)
AH: adiabatic heating, VD: viscous dissipation, HF: heat flux, WG: work against gravity
time (Myr)
Eq.(112) is verified by (l) and Eq.(101) is verified by (m).
86
250000
a)
b)
c)
d)
e)
f)
g)
h)
i)
j)
k)
l)
m)
n)
87
results - EBA - Ra = 105
26.12
These results were obtained with a 64x64 resolution, and CFL number of 1. Simulation was stopped after about 4300
timesteps.
9x10-5
7.7x1022
8x10-5
7x10-5
6x10-5
0
20000
40000
60000
80000
100000
120000
time (Myr)
7.4x1022
(b)
EG
0
20000
40000
60000
80000
100000
5000
Temperature
2240
2230
2220
2210
-7.75x1023
20000
40000
60000
80000
100000
120000
time (Myr)
3500
0.002
3000
2500
2000
-0.002
1000
-0.003
-0.004
0
20000
40000
60000
80000
100000
-0.005
120000
(f)
time (Myr)
0
viscous dissipation
100000
80000
60000
40000
-100000
20000
-120000
0
0
20000
40000
60000
80000
100000
120000
(h)
0
20000
40000
60000
80000
100000
120000
80000
100000 120000
60000
40000
0
0
20000
40000
60000
80000
100000
0.3
6
5.5
5
50000
4.5
0.2
30000
20000
10000
4
Nu
vrms (cm/yr)
40000
0.15
3.5
3
2.5
0.1
2
0
1.5
0.05
-10000
1
0
20000
40000
60000
80000
100000
120000
time (Myr)
20000
0
(k)
0
20000
40000
60000
80000
100000
120000
time (Myr)
0.5
(l)
0
20000
40000
60000
80000
100000
time (Myr)
120000
0
AH+VD+HF
dETdt
10000
VD
WG
100000
0
80000
-10000
-20000
60000
-30000
40000
-40000
20000
-50000
-60000
(l)
120000
time (Myr)
0.25
(j)
60000
80000
(i)
60000
-20000
40000
20000
time (Myr)
70000
20000
time (Myr)
120000
time (Myr)
heat flux
0
1500
100000
(g)
120000
0.001
120000
-80000
100000
-0.001
0
-60000
80000
min(u)
max(u)
min(v)
max(v)
-20000
-40000
60000
0.004
0.003
(e)
40000
0.005
500
0
20000
time (Myr)
4000
0
0
(c)
velocity
2260
2250
120000
min(T)
max(T)
4500
2270
-7.748x1023
time (Myr)
2280
adiabatic heating
-7.746x1023
7.45x1022
1x10-5
work against gravity
EK
ET
7.5x1022
2x10-5
(d)
7.55x1022
-7.744x1023
3x10-5
2200
-7.74x1023
-7.742x1023
4x10-5
(a)
-7.738x1023
7.6x1022
5x10-5
0
-7.736x1023
7.65x1022
0
20000
40000
60000
80000
0
100000 120000
0
20000
40000
60000
80000
100000
120000
time (Myr)
(m)
AH: adiabatic heating, VD: viscous dissipation, HF: heat flux, WG: work against gravity
time (Myr)
88
120000
26.13
Onset of convection
The code can be run for values of Ra between 500 and 1000, at various resolutions for the BA formulation. The value
vrms (t) − vrms (0) is plotted as a function of Ra and for the 10 first timesteps. If the vrms is found to decrease, then
the Rayleigh number is not high enough to allow for convection and the initial temperature perturbation relaxes by
diffusion (and then vrms (t) − vrms (0) < 0. If the vrms is found to increase, then vrms (t) − vrms (0) > 0 and the system
is going to showcase convection. The zero value of vrms (t) − vrms (0) gives us the critical Rayleigh number, which is
found between 775 and 790.
8x10-15
1.5x10-16
16x16
32x32
48x48
64x64
6x10-15
5x10-17
0
vrms trend
4x10-15
vrms trend
16x16
32x32
48x48
64x64
1x10-16
2x10-15
0
-5x10-17
-1x10-16
-1.5x10-16
-2x10-16
-2.5x10-16
-2x10-15
-3x10-16
-4x10-15
-3.5x10-16
770
1000
Ra
775
780
Ra
89
785
790
27
fieldstone: Rayleigh-Taylor instability (1)
This numerical experiment was first presented in [74]. It consists of an isothermal Rayleigh-Taylor instability in a
two-dimensional box of size Lx = 0.9142 and Ly = 1.
Two Newtonian fluids are present in the system: the buoyant
layer
is placed at the bottom of the box and the
πx
interface between both fluids is given by y(x) = 0.2 + 0.02 cos Lx The bottom fluid is parametrised by its mass
density ρ1 and its viscosity µ1 , while the layer above is parametrised by ρ2 and µ2 .
No-slip boundary conditions are applied at the bottom and at the top of the box while free-slip boundary conditions
are applied on the sides.
In the original benchmark the system is run over 2000 units of dimensionless time and the timing and position
of various upwellings/downwellings is monitored. In this present experiment only the root mean square velocity is
measured at t = 0: the code is indeed not yet foreseen of any algorithm capable of tracking deformation.
Another approach than the ones presented in the extensive literature which showcases results of this benchmark
is taken. The mesh is initially fitted to the fluids interface and the resolution is progressively increased. This results
in the following figure:
0.01000
vrms
0.00100
0.00010
0.00001
0.01
0.1
hx
The green line indicates results obtained with my code ELEFANT with grids up to 2000x2000 with the exact same
methodology.
90
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• isothermal
• numerical benchmark
91
References
[1] M. Albers. A local mesh refinement multigrid method for 3D convection problems with strongly variable viscosity.
J. Comp. Phys., 160:126–150, 2000.
[2] K.-J. Bathe. Finite Element Procedures in Engineering Analysis. Prentice-Hall, 1982.
[3] M. Benzi, G.H. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta Numerica, 14:1–137,
2005.
[4] D. Bercovici, G. Schubert, and G.A. Glatzmaier. Three-dimensional convection of an infinite Prandtl-number
compressible fluid in a basally heated spherical shell. J. Fluid Mech., 239:683–719, 1992.
[5] B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch,
G. Marquart, D. Moore, P. Olson, H. Schmeling, and T. Schnaubelt. A benchmark comparison for mantle
convection codes. Geophys. J. Int., 98:23–38, 1989.
[6] H.H. Bui, R. Fukugawa, K. Sako, and S. Ohno. Lagrangian meshfree particles method (SPH) for large deformation
and failure flows of geomaterial using elasticplastic soil constitutive model. Int. J. Numer. Anal. Geomech.,
32(12):1537–1570, 2008.
[7] C. Burstedde, G. Stadler, L. Alisic, L.C. Wilcox, E. Tan, M. Gurnis, and O. Ghattas. Large-scale adaptive mantle
convection simulation. Geophy. J. Int., 192:889–906, 2013.
[8] F.H. Busse, U. Christensen, R. Clever, L. Cserepes, C. Gable, E. Giannandrea, L. Guillou, G. Houseman, H.-C.
Nataf, M. Ogawa, M. Parmentier, C. Sotin, and B. Travis. 3D convection at infinite Prandtl number in Cartesian
geometry - a benchmark comparison. Geophys. Astrophys. Fluid Dynamics, 75:39–59, 1993.
[9] J.S. Chen, C. Pan, and T.Y.P. Chang. On the control of pressure oscillation in bilinear-displacement constantpressure element. Comput. Methods Appl. Mech. Engrg., 128:137–152, 1995.
[10] Edmund Christiansen and Knud D. Andersen. Computation of collapse states with von mises type yield condition.
International Journal for Numerical Methods in Engineering, 46:1185–1202, 1999.
[11] Edmund Christiansen and Ole S. Pedersen. Automatic mesh refinement in limit analysis. International Journal
for Numerical Methods in Engineering, 50:1331–1346, 2001.
[12] M. crouzeix and P.A. Raviart. Conforming and non-conforming finite element methods for solving the stationary
Stokes equations. RAIRO, 7:33–76, 1973.
[13] C. Cuvelier, A. Segal, and A.A. van Steenhoven. Finite Element Methods and Navier-Stokes Equations. D. Reidel
Publishing Company, 1986.
[14] D.R. Davies, C.R. Wilson, and S.C. Kramer. Fluidity: A fully unstructured anisotropic adaptive mesh computational modeling framework for geodynamics. Geochem. Geophys. Geosyst., 12(6), 2011.
[15] P. Davy and P. Cobbold. Indentation tectonics in nature and experiment. 1. experiments scaled for gravity.
Bulletin of the Geological Institutions of Uppsala, 14:129–141, 1988.
[16] J. de Frutos, V. John, and J. Novo. Projection methods for incompressible ow problems with WENO nite
difference schemes. J. Comp. Phys., 309:368–386, 2016.
[17] Y. Deubelbeiss and B.J.P. Kaus. Comparison of Eulerian and Lagrangian numerical techniques for the Stokes
equations in the presence of strongly varying viscosity. Phys. Earth Planet. Interiors, 171:92–111, 2008.
[18] C.R. Dohrmann and P.B. Bochev. A stabilized finite element method for the Stokes problem based on polynomial
pressure projections. Int. J. Num. Meth. Fluids, 46:183–201, 2004.
[19] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. 2003.
[20] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. John Wiley & Sons, 2003.
[21] Jean Donea and Antonio Huerta. Finite Element Methods for Flow Problems. John Wiley & Sons, 2003.
[22] T. Duretz, D.A. May, T.V. Gerya, and P.J. Tackley. Discretization errors and free surface stabilisation in the finite
difference and marker-in-cell method for applied geodynamics: A numerical study. Geochem. Geophys. Geosyst.,
12(Q07004), 2011.
92
[23] R. Eid. Higher order isoparametric finite element solution of Stokes flow . Applied Mathematics and Computation,
162:1083–1101, 2005.
[24] E. Erturk. Discussions on Driven Cavity Flow. Int. J. Num. Meth. Fluids, 60:275–294, 2009.
[25] P.J. Frey and P.-L. George. Mesh generation. Hermes Science, 2000.
[26] M. Gerbault, A.N.B. Poliakov, and M. Daignieres. Prediction of faulting from the theories of elasticity and
plasticity: what are the limits? Journal of Structural Geology, 20:301–320, 1998.
[27] Taras Gerya. Numerical Geodynamic Modelling. Cambridge University Press, 2010.
[28] T.V. Gerya, D.A. May, and T. Duretz. An adaptive staggered grid finite difference method for modeling geodynamic Stokes flows with strongly variable viscosity. Geochem. Geophys. Geosyst., 14(4), 2013.
[29] G.H. Golub and C.F. van Loan. Matrix Computations, 4th edition. John Hopkins University Press, 2013.
[30] P.M. Gresho and R.L. Sani. Incompressible flow and the Finite Element Method, vol II. John Wiley and Sons,
Ltd, 2000.
[31] D. Griffiths and D. Silvester. Unstable modes of the q1-p0 element. Technical Report 257, University of MAnchester/UMIST, 1994.
[32] M. Gunzburger. Finite Element Methods for Viscous Incompressible Flows: A Guide to Theory, Practice and
Algorithms. Academic, Boston, 1989.
[33] T. Heister, J. Dannberg, R. Gassmöller, and W. Bangerth. High Accuracy Mantle Convection Simulation through
Modern Numerical Methods. II: Realistic Models and Problems. Geophy. J. Int., 210(2):833–851, 2017.
[34] T.J.R. Hughes. The Finite Element Method. Linear Static and Dynamic Finite Element Analysis. Dover Publications, Inc., 2000.
[35] T.J.R. Hughes, W.K. Liu, and A. Brooks. Finite element analysis of Incompressible viscous flows by the penalty
function formulation. J. Comp. Phys., 30:1–60, 1979.
[36] Hoon Huh, Choong Ho Lee, and Wei H. Yang. A general algorithm for plastic flow simulation by finite element
limit analysis. International Journal of Solids and Structures, 36:1193–1207, 1999.
[37] J. Ita and S.D. King. Sensitivity of convection with an endothermic phase change to the form of governing
equations, initial conditions, boundary conditions, and equation of state. J. Geophys. Res., 99(B8):15,919–15,938,
1994.
[38] L. Jolivet, P. Davy, and P. Cobbold. Right-lateral shear along the Northwest Pacific margin and the India-Eurasia
collision. Tectonics, 9(6):1409–1419, 1990.
[39] L.M. Kachanov. Fundamentals of the Theory of Plasticity. Dover Publications, Inc., 2004.
[40] S. King, C. Lee, P. van Keken, W. Leng, S. Zhong, E. Tan, N. Tosi, and M. Kameyama. A community benchmark
for 2D Cartesian com- pressible convection in the Earths mantle. Geophy. J. Int., 180:7387, 2010.
[41] J.R. Koseff and R.L. Street. The Lid-Driven Cavity Flow: A Synthesis of Qualitative and Quantitative Observations. J. Fluids Eng, 106:390–398, 1984.
[42] M. Kronbichler, T. Heister, and W. Bangerth. High accuracy mantle convection simulation through modern
numerical methods . Geophy. J. Int., 191:12–29, 2012.
[43] R. Lee, P. Gresho, and R. Sani. Smooting techniques for certain primitive variable solutions of the Navier-Stokes
equations. . Int. Journal for Numerical Methods in Engineering, 14:1785–1804, 1979.
[44] W. Leng and S. Zhong. Viscous heating, adiabatic heating and energetic consistency in compressible mantle
convection. Geophy. J. Int., 173:693–702, 2008.
[45] W. Leng and S. Zhong. Implementation and application of adaptive mesh refinement for thermochemical mantle
convection studies. Geochem. Geophys. Geosyst., 12(4), 2011.
[46] X. Liu and S. Zhong. Analyses of marginal stability, heat transfer and boundary layer properties for thermal
convection in a compressible fluid with infinite Prandtl number. Geophy. J. Int., 194:125–144, 2013.
93
[47] D.S. Malkus and T.J.R. Hughes. Mixed finite element methods - reduced and selective integration techniques: a
unification of concepts. Comput. Meth. Appl. Mech. Eng., 15:63–81, 1978.
[48] L.E. Malvern. Introduction to the mechanics of a continuous medium. Prentice-Hall, Inc., 1969.
[49] P. Molnar and P. Tapponnier. Relation of the tectonics of eastern China to the India-Eurasia collision: Application
of the slip-line field theory to large-scale continental tectonics. Geology, 5:212–216, 1977.
[50] L. Moresi, S. Quenette, V. Lemiale, C. Mériaux, B. Appelbe, and H.-B. Mühlhaus. Computational approaches to
studying non-linear dynamics of the crust and mantle. Phys. Earth. Planet. Inter., 163:69–82, 2007.
[51] C. O’Neill, L. Moresi, D. Müller, R. Albert, and F. Dufour. Ellipsis 3D: a particle-in-cell finite element hybrid
code for modelling mantle convection and lithospheric deformation. Computers and Geosciences, 32:1769–1779,
2006.
[52] G. Peltzer and P. Tapponnier. Formation and evolution of strike-slip faults, rifts, and basins during the india-asia
collision: an experimental approach. J. Geophys. Res., 93(B12):15085–15177, 1988.
[53] A. Pinelli and A. Vacca. Chebyshev collocation method and multidomain decomposition for the incompressible
Navier-Stokes equations. International Journal for numerical methods in fluids, 18:781–799, 1994.
[54] A.E. Pusok, B.J.P. Kaus, and A.A. Popov. On the Quality of Velocity Interpolation Schemes for Marker-in-Cell
Method and Staggered Grids. Pure and Applied Geophysics, pages doi:10.1007/s00024–016–1431–8, 2016.
[55] T. Rabczuk, P.M.A. Areias, and T. Belytschko. A simplified mesh-free method for shear bands with cohesive
surfaces . Int. J. Num. Meth. Eng., 69:993–1021, 2007.
[56] J.N. Reddy. On penalty function methods in the finite element analysis of flow problems. Int. J. Num. Meth. Fluids, 2:151–171, 1982.
[57] J. Revenaugh and B. Parsons. Dynamic topography and gravity anomalies for fluid layers whose viscosity varies
exponentially with depth. Geophysical Journal of the Royal Astronomical Society, 90(2):349–368, 1987.
[58] Y. Saad. Iterative methods for sparse linear systems. SIAM, 2003.
[59] R.L. Sani, P.M. Gresho, R.L. Lee, and D.F. Griffiths. The cause and cure (?) of the spurious pressures generated
by certain FEM solutions of the incompressible Navier-Stokes equations: part 1. Int. J. Num. Meth. Fluids,
1:17–43, 1981.
[60] R.L. Sani, P.M. Gresho, R.L. Lee, D.F. Griffiths, and M. Engelman. The cause and cure (?) of the
spurious pressures generated by certain fem solutions of the incompressible navier-stokes equations: part 2.
Int. J. Num. Meth. Fluids, 1:171–204, 1981.
[61] H. Schmeling, A.Y. Babeyko, A. Enns, C. Faccenna, F. Funiciello, T. Gerya, G.J. Golabek, S. Grigull, B.J.P.
Kaus, G. Morra, S.M. Schmalholz, and J. van Hunen. A benchmark comparison of spontaneous subduction
models - Towards a free surface. Phys. Earth. Planet. Inter., 171:198–223, 2008.
[62] D.W. Schmid and Y.Y. Podlachikov. Analytical solutions for deformable elliptical inclusions in general shear.
Geophy. J. Int., 155:269–288, 2003.
[63] G. Schubert, D.L. Turcotte, and P. Olson. Mantle Convection in the Earth and Planets. Cambridge University
Press, 2001.
[64] J. Suckale, J.-C. Nave, and B.H. Hager. It takes three to tango: 1. Simulating buoyancy-driven flow in the
presence of large viscosity contrasts. J. Geophys. Res., 115(B07409), 2010.
[65] P. Tackley. Three-dimensional models of mantle convection: Influence of phase transitions and temperaturedependent viscosity. PhD thesis, California Institute of Technology, 1994.
[66] E. Tan and M. Gurnis. Compressible thermochemical convection and application to lower mantle structures.
J. Geophys. Res., 112(B06304), 2007.
[67] Paul Tapponnier and Peter Molnar. Slip-line field theory and large-scale continental tectonics. Nature, 264:319–
324, November 1976.
[68] M. Thielmann, D.A. May, and B.J.P. Kaus. Discretization errors in the Hybrid Finite Element Particle-In-Cell
Method. Pure and Applied Geophysics, 2014.
94
[69] C. Thieulot. GHOST: Geoscientific Hollow Sphere Tesselation. Solid Earth, 9(1–9), 2018.
[70] C. Thieulot, P. Fullsack, and J. Braun. Adaptive octree-based finite element analysis of two- and three-dimensional
indentation problems. J. Geophys. Res., 113:B12207, 2008.
[71] J.F. Thompson, B.K. Soni, and N.P. Weatherill. Handbook of grid generation. CRC press, 1998.
[72] R.A. Trompert and U. Hansen. On the Rayleigh number dependence of convection with a strongly temperaturedependent viscosity. Physics of Fluids, 10(2):351–360, 1998.
[73] D.L. Turcotte and G. Schubert. Geodynamics, 2nd edition. Cambridge, 2012.
[74] P.E. van Keken, S.D. King, H. Schmeling, U.R. Christensen, D. Neumeister, and M.-P. Doin. A comparison of
methods for the modeling of thermochemical convection. J. Geophys. Res., 102(B10):22,477–22,495, 1997.
[75] H. Xing, W. Yu, and J. Zhang. In Advances in Geocomputing, Lecture Notes in Earth Sciences. Springer-Verlag,
Berlin Heidelberg, 2009.
[76] X. Yu and F. Tin-Loi. A simple mixed finite element for static limit analysis. Computers and Structures, 84:1906–
1917, 2006.
[77] S. Zhong. Analytic solutions for Stokes flow with lateral variations in viscosity. Geophys. J. Int., 124:18–28, 1996.
[78] O. Zienkiewicz and S. Nakazawa. The penalty function method and its application to the numerical solution of
boundary value problems. The American Society of Mechanical Engineers, 51, 1982.
[79] O.C. Zienkiewicz, M. Huang, and M. Pastor. Localization problems in plasticity using finite elements with
adaptive remeshing. International Journal for Numerical and Analytical Methods in Geomechanics, 19:127–148,
1995.
[80] O.C. Zienkiewicz, C. Humpheson, and R.W. Lewis. Associated and non-associated visco-plasticity and plasticity
in soil mechanics . Géotechnique, 25(4):671–689, 1975.
95
Index
Pm × Pn , 12
Pm × P−n , 12
Q1 × P0 , 36, 42, 44, 50, 56, 59, 75, 81, 91
Q2 × Q1 , 12, 61, 66
Q3 × Q2 , 69
Qm × P−n , 12
Qm × Qn , 12
Qm × Q−n , 12
second viscosity, 7
static condensation, 32
structured grid, 27
analytical solution, 36, 44, 59, 66, 69, 72, 75, 81
arithmetic mean, 56
weak form, 20
work against gravity, 78
thermal expansion, 76
trapezoidal rule, 9
unstructured grid, 27
basis functions, 16
Boussinesq, 8
bubble function, 12
bulk modulus, 76
bulk viscosity, 7
checkerboard, 30
Compressed Sparse Column, 24
Compressed Sparse Row, 24
compressibility, 76
compressible flow, 81
conforming element, 12
connectivity array, 28
convex polygon, 27
CSC, 24
CSR, 24
dynamic viscosity, 7
Gauss quadrature, 10
geometric mean, 56
harmonic mean, 56
incompressible flow, 44, 50, 56, 59, 66, 69, 72, 75, 91
isothermal, 36, 42, 44, 50, 56, 59, 66, 69, 72, 75, 91
isoviscous, 36, 42, 59, 66, 69, 75, 81
Legendre polynomial, 10
midpoint rule, 9
mixed formulation, 59, 66, 69, 72, 75, 81, 91
Newton-Cotes, 10
non-conforming element, 12
non-isoviscous, 44, 50, 56, 72
nonconforming Q1 × P0 , 72
nonlinear rheology, 50
numerical benchmark, 91
particle-in-cell, 56
penalty formulation, 19, 36, 42, 44, 50, 56
piecewise, 12
pressure normalisation, 62
pressure smoothing, 30, 59, 72, 75, 81
quadrature, 10
rectangle rule, 9
96
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.5 Linearized : No Page Count : 96 Page Mode : UseOutlines Author : Title : Subject : Creator : LaTeX with hyperref package Producer : pdfTeX-1.40.18 Create Date : 2019:01:06 22:10:33+01:00 Modify Date : 2019:01:06 22:10:33+01:00 Trapped : False PTEX Fullbanner : This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) kpathsea version 6.2.3EXIF Metadata provided by EXIF.tools