Introduction
Genuinely multidimensional methods, known as Residual Distribution (RD) or Fluctuation Split-
ting (FS) schemes, are a well-established technique for solving hyperbolic problems on triangular
meshes, see [3] and references therein, and [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. On the other hand,
RD methods for quadrilateral grids, which, for example, are the optimal choice for discretizing
boundary layer regions, have not yet reached the same level of maturity and reliability. The
origin of residual distribution methods for quadrilateral cells can be traced back to the work
of Ni [14], who proposed a Lax–Wendroff-type scheme with second-order accuracy in space and
time. Then, several upwind schemes have been designed, such as the ones provided in [15]
and [16], where the residual is distributed to downstream nodes with respect to the advection
velocity. Furthermore, the RD schemes successfully developed for triangular cells, such as the
N, LDA, LW, and PSI schemes, have been extended to cell-vertex quadrilateral grids in [17]. A
few years ago, a method has been proposed to apply the RD schemes to arbitrary conservation
laws [18] and finite elements [19] in the absence of a conservative linearization [20]. In spite
of these efforts, a thorough analysis of the accuracy and stability properties of RD schemes
on rectangular grids is still lacking and is worth pursuing. More recently, in [21] an extension
of the Lax-Wendroff like theorem for RD schemes on quadrilaterals has been proposed, along
with the N-modified scheme with the addition of an SUPG based artificial dissipation term to
improve convergence. In [22] a WENO reconstruction is employed to obtain the nodal values
and integrate the fluctuation, thus naturally leading to high-order accuracy. Using the WENO
reconstruction in the framework of distribution schemes has advantages in terms of computa-
tional cost over WENO finite volume schemes, and in terms of mesh-regularity requirements
over WENO finite difference schemes.
The motivations of this work are manifold. Firstly, it tries to address some of the issues aris-
ing when extending RD schemes from triangular grids to quadrilateral ones, like the marginal
stability of high-wavenumber modes for LP schemes. From another and related point of view,
it tries to address the problem of a consistent discretization of the convective and the diffusive
terms in an RD framework, an argument to which also other researchers are devoting more
2and more effort, see for instance [23, 24, 25]. The Fourier and truncation error analyses on a
structured mesh have been used, in order to get more insight into FS schemes and to provide
a strategy for the design of second-order-accurate stable schemes. In particular, the proposed
spectral analysis is multidimensional in the sense that it is not restricted to the case of Fourier
modes aligned with the advection velocity, as done for example in [26, 27, 28], and all of the
grid-resolved two-dimensional wavenumbers are analyzed for a given advection angle.
The final aim of this thesis is to make a step towards the design of stable and accurate FS
schemes for the solution of the Euler and Navier-Stokes equations on unstructured hybrid meshes,
which are commonly used in commercial and industrial codes.
This thesis is organized as follows. Chapter 1 introduces the theory of hyperbolic systems of
conservation laws, focusing in particular on the scalar conservation law and on the system of
Euler equations, which are the subject of numerical tests in subsequent sections. Chapter 2
presents the procedure to obtain an FS discretization assuming the scalar conservation law
as the model partial differential equation (PDE), both on triangular and quadrilateral meshes.
In particular, the two main steps of RD schemes are described: i) evaluating the fluctuation,
namely, the flux balance over the cell; ii) distributing the fluctuation to its vertices. Concerning
the advection term, the cell residual is evaluated using a contour integral of the flux along the
cell faces, which assumes linear variation of the unknown and thus avoids the need for a conser-
vative linearization. For the distribution step, the linearity preserving LDA and LW schemes [17]
have been used as a starting point. It is noteworthy that for the case of scalar advection, the
extension to Cartesian grids of the LW scheme for triangles recovers Ni’s scheme [14]. Then, a
conservative SUPG method for quadrilateral cells is derived. Concerning the diffusion term, it
is discretized either by a standard Galerkin finite element (FE) scheme or by a residual based
approach. In the latter case, the diffusive fluctuation is evaluated by a contour integral, whereby
the gradients at the nodes are computed by a Green-Gauss reconstruction. For the distribution
of the diffusive fluctuation, one could choose between an upwind or a centred strategy. The two
possibilities have been discussed thereafter. Chapter 3 addresses the issue of accuracy in the
discretization of advection and advection-diffusion equations using the truncation error analysis
based on 2D Taylor series expansions. Relations between FS, FEand Finite Difference (FD)
schemes are underlined, and the conditions for a hybrid Linearity Preserving (LP) (for advec-
tion) - finite element Galerkin (FE-GAL) (for diffusion) approach to be second-order-accurate
are defined, both on triangles and quadrilaterals. The Fourier analysis is applied in chapter 4,
which shows that LP schemes are marginally stable when applied to quadrilateral cells, thus
requiring an additional numerical dissipation to compute advection problems. In this thesis,
Introduction 3
LP schemes are effectively stabilized by the addition of on the bias term of the SUPG scheme,
without affecting their order of accuracy. A thorough analysis of second-order-accurate dis-
cretizations of the advection-diffusion equation is then performed. This has been accomplished
developing a mathematical framework for the multidimensional spectral analysis, both of dis-
persive and dissipative properties of FD, FE , and FS schemes. In chapter 5 the truncation
error and Fourier analyses are used to design a second-order-accurate minimum-dispersion error
RD linearity preserving scheme which is suitable for the discretization of the convective term
in advection-diffusion problems. The theoretical findings are finally validated in chapter 6 by a
series of numerical experiments carried out for the scalar conservation law. Chapter 7 is devoted
to the Euler equations. Matrix schemes are introduced and some test-cases are used to validate
the method. The computation of a turbine cascade flow is proposed using an hybrid grid, to
show the flexibility of the method and its applicability in a near future in commercial codes.
Finally, chapter 8 addresses some of the problems arising when trying to adopt higher-order
polynomial interpolation for the solution reconstruction (to achieve higher-order accuracy). In a
successive step, the Fourier analysis is used to design and optimize multidimensional (implicit)
compact schemes, which are very suitable for aeroacoustic and Large Eddy Simulation (LES)
for their excellent dispersive properties.
Chapter 1
Hyperbolic systems of conservation
laws
The governing equations of fluid dynamics are usually obtained applying a conservation state-
ment for an intensive property of the flow over a given space domain. The conservation laws for
matter (continuity equation), momentum (corresponding to the Newton’s law), and energy (first
principle of thermodynamics) represent the fundamental starting point to obtain the Navier–
Stokes and the Euler equations of gas-dynamics. In this chapter the general theory on the
numerical approximation of hyperbolic systems of conservation laws is presented, with a partic-
ular focus on the system of Euler equations, model of inviscid compressible flow with negligible
heat conductivity, and on the scalar advection equation, which is the simplest model of hyper-
bolic equation. In depth theory can be found in reference [29], to which this chapter is largely
inspired.
1.1 Conservation laws: mathematical formulation
Consider p conserved variables qm depending on space and time, and gathered into a vector-
valued function Q:
qm = qm(x, t), for m = 1, . . . , p;
Q =
q1
.
.
.
qp
, with Q : Rd × [0,∞[ → Ω ,
where d is the number of space dimensions, x = (x1, . . . , xd), t is the time variable, and Ω is
an open subset of Rp, called set of states. We call qm conserved variables since they represent
intensive1 quantities for which we can state a conservation principle. Given an arbitrary domain,
1For intensive we mean for unit-volume; e.g. for the mass conservation the intensive quantity is the mass for
unit volume, i.e. the density ρ [Kg/m3], whereas the associated extensive quantity is, obviously, the mass itself.
6 Hyperbolic systems of conservation laws
D ⊂ Rd, consider D as a control volume in which we define gDm as the (scalar) extensive quantities
associated to each qm:
gDm(t) =
∫
D
qm(x, t) dx for m = 1, . . . , p .
Let nˆ = (n1, . . . , nd)T be the outward pointing unit vector normal to the boundary ∂D of D.
Stating a conservation principle over D, we say that the variation of each qm in time is only due
to the flux through the boundary ∂D, namely:
d
dt
∫
D
Q dx +
∮
∂D
F · nˆ dS = 0 (1.1)
where the tensor F is defined as
F = (F1, . . . ,Fd),
and Fi are the flux vector function in the i-th direction, which contains the p scalar flux-functions
fm,i:
Fi =
f1,i
.
.
.
fp,i
, with Fi : Ω → Rp for i = 1, . . . , d ;
the scalar flux-function fm,i defines the specific flux of qm in the i-th direction.
Remark 1.1. The flux-functions depend explicitly on Q and implicitly on x and t through Q :
F = F(Q) = F(Q(x, t)).
In case production-destruction phenomena are taken into account, the integro-differential
system of equations expressed by (1.1) is enriched by the so called source term:
d
dt
∫
D
Q dx
︸ ︷︷ ︸
time variation of conserved
quantities over D
+
∮
∂D
F · nˆ dS
︸ ︷︷ ︸
flux through the boundary ∂D
=
∫
D
S dx
︸ ︷︷ ︸
production-destruction term
which requires the definitions of the vector-valued function S.
Remark 1.2. In the more general case the source terms depend on the state Q, position x, and
time t : S = S(Q,x, t).
From now on we will always refer to model equations like (1.1), i.e. without source terms.
1.1 Conservation laws: mathematical formulation 7
1.1.1 Differential forms
Starting from the system of p conservation laws expressed by (1.1) in integro-differential form,
we derive the conservative differential form, also referred to as the divergence form. Applying
the Green-Gauss theorem to the flux term we can write
∮
∂D
F · nˆ dS =
∫
D
∇·F dx ,
while for the first term, since D does not depend on time, we can transfer the time derivative
onto Q as follows:
d
dt
∫
D
Q dx =
∫
D
∂Q
∂t
dx ,
and finally we get
∫
D
[
∂Q
∂t
+ ∇·F
]
dx = 0 .
Since the last equation holds for any arbitrary domain D′ ⊂ D, in the limit of vanishing D′ (i.e.
|D| → 0) we obtain the differential form, valid for all x ∈ D:
∂Q
∂t
+ ∇·F = 0 conservative or divergence form (1.2)
Remark 1.3. A differential equation is said to be in divergence or conservative form when the
coefficients multiplying the derivatives do not contain variables which appear in the derivatives.
Remark 1.4. The differential form (1.2) requires, with respect to the integral one (1.1), more
coercive requirements on the functions Q and F : indeed in the latter the integrability of Q
on D and F on ∂D is required, while, in the former, Q and F should be both C 1 on D (i.e.
continuous and differentiable).
1.1.2 Hyperbolic systems
In order to define the hyperbolic character of a system of conservation law, let us start defining
the Jacobian matrices Ai collected into A:
A = (A1, . . . ,Ad) with Ai(Q) = ∂Fi∂Q for i = 1, . . . , d,
such that another form of the differential equation, called quasi-linear form, is obtained from (1.2):
∂Q
∂t
+ A ·∇Q = 0 quasi-linear form . (1.3)
Definition 1.1. The system (1.3) is (strictly) hyperbolic if, for any Q ∈ Ω and any ω =
(ω1, . . . , ωd) ∈ Rd the matrix A(Q,ω) = ω ·A has p real (and distinct) eigenvalues, λm, and p
linearly independent eigenvectors, rm, such that A rk = λk rk for all k = 1, . . . , p.
8 Hyperbolic systems of conservation laws
1.1.3 Cauchy or initial value problem
The Cauchy problem, also referred to as initial value problem (IVP), consists in finding the
function Q(x, t) satisfying
1. equation (1.2) ∀ (x, t) ∈ Rd × [0,∞[ and
2. the initial condition
Q(x, 0) = Q0(x) ∀ x ∈ Rd (1.4)
where Q0 : Rd → Ω is a given function.
The Riemann problem
In one space dimension, the Cauchy problem with initial condition given by
Q0(x) =
{
Q`, for x > 0
Qr, for x < 0
is called the (one-dimensional) Riemann problem, and is of great importance both in the exper-
imental and numerical gas-dynamics.
1.1.4 Initial boundary value problem
The initial boundary value problem (IBVP) is a generalization of the Cauchy problem to the
case of conservation laws applied to a space domain D ⊂ Rd. The function Q(x, t) is sought
which satisfies:
1. equation (1.2) ∀ (x, t) ∈ D × [0,∞[ ;
2. the initial condition Q(x, 0) = Q0(x) ∀ x ∈ D;
3. the boundary conditions
Q(x, t) = QB(x, t) ∀ x ∈ ∂D and t > 0 ,
where QB (x, t) : ∂D × [0,∞[ → Ω is a given function.
Remark 1.5. A fundamental issue related to initial boundary value problems is the way in
which boundary conditions affect the well-posedness of the problem itself, i.e. there may be no
solution or one that does not depend in a continuous way on the initial or boundary data or it
may be non-unique.
1.2 Classical and weak solutions 9
1.2 Classical and weak solutions
Definition 1.2. A classical solution to the Cauchy problem is a C 1 function Q which satis-
fies (1.2) and (1.4) point-wisely.
Unfortunately for a system of conservation laws there does not exist, in general, a classical
solution beyond some finite time interval, even for initial conditions given by a smooth function
Q0. Indeed discontinuities may arise. Here we introduce weak solutions, a wider class which
embraces both continuous and discontinuous solutions.
Let C 10 (Rd× [0,∞[) denote the space of C 1 functions with compact support in Rd× [0,∞[ , with
the vector function ϕ = (ϕ1, . . . , ϕp) made of p scalar functions ϕm ∈ C 10 (Rd× [0,∞[) for m =
1, . . . , p. Let us start from equation (1.2), multiplying by ϕ and integrating on Rd× [0,∞[ ; then
by Green’s theorem (or integrating by parts) we get
∫ +∞
0
∫
Rd
[
∂Q
∂t
+∇·F
]
·ϕ dxdt = −
∫ +∞
0
∫
Rd
[
Q · ∂ϕ
∂t
+ F ·∇ϕ
]
dxdt
+
∫ +∞
0
∫
Rd
[
∂(Q · ϕ)
∂t
+ ∇· (F · ϕ)
]
dxdt.
The second term at the right-hand-side can be further simplified noticing that
∫ +∞
0
∫
Rd
∂
∂t
(Q · ϕ) dxdt =
∫
Rd
∣
∣
∣Q · ϕ
∣
∣
∣
+∞
0
dx = −
∫
Rd
Q(x, 0) ·ϕ(x, 0) dx ,
where, since ϕ has compact support in Rp it turns out that lim
t→∞ϕ(x, t) = 0 ∀x ∈ R
d
. As a
consequence, the term for t→ ∞ vanishes. Moreover, the other term related to the flux vector
vanishes since:
∫ +∞
0
∫
Rd
∇· (F · ϕ) dxdt =
∫ +∞
0
[∮
∂D→∂Rd
ϕ · (F · nˆ) dS
]
dt = 0 ,
exploiting, once more, the compact support of ϕ in Rp. Finally the following equivalence is
obtained:
−
∫ +∞
0
∫
Rd
[
∂Q
∂t
+ ∇·F
]
·ϕ dxdt =
=
∫ +∞
0
∫
Rd
[
Q · ∂ϕ
∂t
+ F ·∇ϕ
]
dxdt +
∫
Rd
Q(x, 0) ·ϕ(x, 0) dx ,
and it is a simple matter to check that a classical solution of the Cauchy problem satisfies both
10 Hyperbolic systems of conservation laws
the following equations for all ϕ:
∫ +∞
0
∫
Rd
[
∂Q
∂t
+ ∇·F
]
· ϕ dxdt = 0 , (1.5)
∫ +∞
0
∫
Rd
[
Q · ∂ϕ
∂t
+ F ·∇ϕ
]
dxdt +
∫
Rd
Q(x, 0) ·ϕ(x, 0) dx = 0 . (1.6)
Remark 1.6. Equations (1.5), (1.6) are equivalent for a classical solution; for a discontinuous
function Q, the derivatives of Q and F do not exist for all x and t making equation (1.5) useless;
on the other hand, equation (1.6) involves only derivatives of ϕ (which is C 1).
Definition 1.3. A function Q is a weak solution of the Cauchy problem if Q(x, t) ∈ C 1 a.e.
(in the sense of distribution) and satisfies (1.6) for any function ϕ.
We observe that:
i) a classical solution is, by construction, also a weak solution;
ii) a weak solution is a classical solution in any domain where Q is C 1;
iii) weak solutions are piecewise smooth, namely there exists a finite number of smooth ori-
entable surfaces Σi in the (x, t) plane, outside of which Q is a C 1 function.
1.2.1 Rankine-Hugoniot conditions
Across a discontinuity surface Σ, Q has a jump discontinuity but, since each qm obeys a conser-
vation principle, not every discontinuity is admissible for a given system of conservation laws.
Let nΣ = (nt, nx1 , . . . , nxd)T = (nt, nˆΣ)T be the vector normal to Σ at a point P ∈ Σ, in which
nt is the time component and nˆΣ ∈ Rd is the unit vector normal to Σ at each time t > 0. With
Q+ and Q− we denote the limits of the function Q approaching to Σ from both sides:
Q± = lim
ε→0+Q ( (x, t)± εnΣ ) .
Consider a small region D ⊂ Rd × [0,∞[ containing P, and crossed by a discontinuity surface
(i.e. Σ). We denote by D± the two open components in which Σ divides D, assuming, for
instance, that nΣ points toward D+. Denoting with wfvec the vector valued function (with p
elements) chosen in the space of infinitely differentiable functions with compact support in D,
namely ϕ ∈ C∞0 (D)p, we obtain from (1.6):
∫
D
[
Q · ∂ϕ
∂t
+ F ·∇ϕ
]
dxdt =
∫
D+
[•] dxdt +
∫
D−
[•] dxdt = 0 ,
Ricevi informazioni sui nostri servizi, sulle offerte e non perdere news e consigli su università e lavoro.
