Koren B., Vuik K. (Editors) Advanced Computational Methods in Science and Engineering

Подождите немного. Документ загружается.

Least-Squares Spectral Element Methods in Computational Fluid Dynamics

that the full potential of solving differential equations by minimizing well-posed

variational statements was recognized in mechanical engineering.

2.2 Galerkin formulations

It is not always possible to convert a differential equation to an associated mini-

mization problem. A well-known example with applications in ﬂuid dynamics is

the steady linear advection equation given by

= 0 , x ∈(0,1) , a > 0 , (21)

with boundary condition u(0) = u

. The method described by (16) is still applicable.

This gives

a(u,v) =

Ω

= 0 . (22)

This formulation where the partial differential equation is weighted by a suitably

chosen set of test functions is called the Galerkin formulation. This weak formu-

lation does not correspond to a minimization problem and generally the space of

trial solutions, u, and test functions, v, are not the same. For the integral (22) to

make sense we need to have that u ∈ H

(0,1) and v ∈ L

(0,1), so a setting in terms

of Hilbert spaces is not possible. Consider the slightly more general case given by

Seek u ∈W such that: a(u,v) = f (v) , ∀v ∈V , (23)

where W and V are function spaces equipped with the norms k·k

and k·k

respectively. W is called the solution space and V is called the test space. We assume

that the bi-linear form is continuous as deﬁned in Deﬁnition 1 and f ∈V

′

. Then the

following theorem establishes the conditions for well-posedness, [31]

Theorem 1. Let W be a Banach space and let V be a reﬂexive Banach space. Let

a ∈ L (W ×V;R) and f ∈V

′

, then problem (23) is well-posed if and only if

∃

, inf

w∈W

sup

v∈V

a(w, v)

kwk

kvk

≥

, (24)

∀v ∈V , (∀w ∈W , a(w,v) = 0) =⇒ (v = 0) . (25)

⊓⊔

This condition is referred to as the inf-sup condition, the Ladyzhenskaya-Babuˇska-

Brezzi (LBB) condition or the Banach-Neˇcas-Babuˇska condition. It is generally

quite hard for a given weak formulation to prove the existence of the positive con-

stant

, the coercivity constant.

In contrast to the minimization problem, well-posedness on the continuous level

does not imply well-posedness for discrete conforming subspaces W

⊂W and V

⊂

183

Marc Gerritsma and Bart De Maerschalck

V. This means that for every problem discrete function spaces need to be identiﬁed

such that well-posedness also holds for (W

). Different weak formulations then

require different approximating subspaces. Such may be the case with multi-physics

simulations such as ﬂuid-structure interaction problems.

In order to establish well-posedness and to develop a consistent and convergent

scheme a variety of stabilization operators are generally required.

2.3 Least-squares formulation

An alternative approach of establishing equivalence between differential equations

and weak formulations is offered by the so-called least-squares formulation. The

least-squares formulation transforms partial differential equations into a minimiza-

tion problem and in this sense extends the Rayleigh-Ritz formulation. The Rayleigh-

Ritz formulation is limited to elliptic, self-adjoint problems. The previous subsec-

tion, 2.1 and 2.2, demonstrated that it is favorable to have a weak formulation which

corresponds to a minimization problem, but that it is not always possible to convert

a partial differential equation in terms of a minimization problem. The least-squares

formulation provides a framework in which partial differential equations can be cast

into a minimization setting, thus retrieving the appealing properties of the Rayleigh-

Ritz principle.

This is established by introducing a norm equivalence between the residual and

the error. Two norms, k·k

and k·k

, are called equivalent when two positive con-

stants, C

and C

, exist such that

kuk

≤ kuk

≤C

kuk

, ∀u ∈X . (26)

If we deﬁne the linear spaces to be

X = {u|kuk

≤ ∞} and Y = {u|kuk

≤ ∞} , (27)

norm-equivalence states that both sets contain the same elements, X = Y. Equivalent

norms on a function space X deﬁne the same topology and more speciﬁcally, Cauchy

sequences in both norms are the same. For ﬁnite-dimensional spaces all norms are

equivalent.

Let the abstract ﬁrst order partial differential equation be given by

L u(x) = f(x) , ∀x ∈

Ω

, (28)

with

Ru(x) = g(x) , ∀x ∈

∂Ω

, (29)

where L is a ﬁrst order linear partial differential operator and R is a linear trace

operator which prescribes boundary values. In the least-squares formulation, we

now look for function spaces, X and Y , with associated norms such that we have the

184

Least-Squares Spectral Element Methods in Computational Fluid Dynamics

norm-equivalence

kuk

Ω

)

≤ kL uk

Ω

)

+ kRuk

∂Ω

)

≤Ckuk

Ω

)

, ∀u ∈X . (30)

The space X is called the solution space and Y is called the residual space. Due to

the linearity of the differential and trace operator we have that, if u

∈X

ku −u

Ω

)

≤ kL u − f k

Ω

)

+ kRu −gk

∂Ω

)

≤Cku −u

Ω

)

, ∀u ∈X ,

(31)

in which u

denotes the exact solution of (28-29). This equivalence tells us that

small residual norms correspond to small error norms and conversely, small error

norms correspond to small residual norms.

Based on this equivalence we can now deﬁne the so-called least-squares func-

tional

J (u) =

kL u − f k

Ω

)

+ kRu −gk

∂Ω

)

. (32)

Solving the abstract PDE (28-29) is now equivalent to minimizing the least-squares

functional J (u). So by applying the least-squares formulation, the problem has

been recast in terms of a minimization problem thus avoiding stability prerequisites

such as inf-sup conditions as discussed in the previous section.

Minimization of the least-squares functional requires setting DJ (u)v = 0 for all

v ∈ X analogous to the minimization procedure described in Section 2.1. Note also

that the boundary conditions, Ru = g, are also incorporated in the minimization

process.

Without loss of generality, we can set g = 0 in the equation Ru = g. This al-

lows us to incorporate the boundary conditions in the space X. When we do so, the

boundary conditions are strongly enforced and can be dropped from the variational

statement.

If Y is a Hilbert space this gives the weak formulation

(L u,L v)

Ω

)

= ( f ,L v)

Ω

)

, ∀v ∈ X . (33)

Note that this formulation is symmetric; interchanging u and v on the left hand side

of this equation leaves the expression unchanged. Furthermore, the basic existence

and uniqueness criterion, the norm-equivalence, (30), is inherited by conforming

subspaces.

Although in the previous discussion general function spaces X and Y were

treated, from a practical point of view it is convenient to choose function spaces

which allow for easy calculation of the residual norm in the minimization process.

So in many applications, the boundary conditions are already incorporated in the

function space X and one takes Y (

Ω

) = L

(

Ω

). This leads to the least-squares func-

tional

J (u) =

kL u − f k

(

Ω

)

. (34)

185

Marc Gerritsma and Bart De Maerschalck

Minimization is in the L

-sense, so point-wise convergence, in the L

∞

-sense, is

only attained if the exact solution is sufﬁciently regular. The next section discusses

the construction of conforming ﬁnite-dimensional subspaces of the solution space

3 Spectral element methods

Instead of seeking the minimizer over the inﬁnite-dimensional space X we restrict

our search to a conforming subspace X

⊂X by performing a domain decomposition

where the solution within each sub-domain is expanded with respect to a polynomial

basis. The domain

Ω

is sub-divided into K non-overlapping quadrilateral open sub-

domains

Ω

[

k=1

Ω

∩

Ω

= /0 , k 6= l . (35)

Each sub-domain is mapped onto the unit cube (−1, 1)

, where d = dim(

Ω

). Within

this unit cube the unknown function is approximated by polynomials. Generally a

spectral element method based on Legendre polynomials, L

(x) over the interval

[−1,1], is employed, [15, 22, 52]. We deﬁne the Gauss-Lobatto-Legendre (GLL)

nodes by the zeroes of the polynomial



1 −x



′

(x) , (36)

and the Lagrange polynomials, h

(x), through these GLL-points, x

, by

(x) =

N(N + 1)

−1)L

′

(x)

)(x −x

)

for i = 0, .. ., N , (37)

where L

′

(x) denotes the derivative of the Nth Legendre polynomial. For multi-

dimensional problems tensor products of the one-dimensional basis functions are

employed in the expansion of the approximate solution. We can therefore expand

the approximate solution in each sub-domain in terms of a truncated series of these

Lagrangian basis functions, which for d = 2 yields

(x,y) =

∑

i=0

∑

j=0

ˆu

i j

(x)h

(y), (38)

where the ˆu

i j

’s are to be determined by the least-squares method. If we have con-

verted a general higher order PDE to an equivalent ﬁrst order system, C

-continuity

sufﬁces to patch the solutions on the individual sub-domains together.

The integrals appearing in the least squares formulation, (34), are approximated

by Gauss-Lobatto quadrature

186

Least-Squares Spectral Element Methods in Computational Fluid Dynamics

−1

f (x)dx ≈

∑

i=0

f (x

, (39)

where w

are the GL weights given by

P(P + 1)

)

, i = 0,. .. ,P ≥N . (40)

It has been shown in [19] that it is beneﬁcial for non-linear equations possessing

large gradients to choose the integration order P higher than the approximation of

the solution, N.

The method is not restricted to higher order methods based on Legendre poly-

nomials. In [18] and [71] Lagrangian basis functions based on the Chebyshev poly-

nomials were used for non-linear, time-dependent hyperbolic equations and the in-

compressible Navier-Stokes equations, respectively. Other systems of orthogonal

polynomials can be easily introduced into the least-squares spectral element frame-

work.

From a practical point of view, only least-squares formulations which allow for

the use of C

-ﬁnite or spectral elements are usable. Since C

-ﬁnite or spectral el-

ement methods are based on piecewise continuously differentiable polynomials,

standard ﬁnite and spectral elements can be used, which results in a very practi-

cal method from an implementational point of view. This can be accomplished by

ﬁrst transforming the system into a ﬁrst order system and subsequently requiring

that only (scaled) L

-norms are used in the quadratic least-squares functional, see

(34). The transformation into a ﬁrst order system has two important consequences.

First of all, the continuity requirements between neighboring spectral elements are

mitigated such that C

-ﬁnite or spectral elements can be used (in case the residuals

are measured by L

-norms). Secondly, the transformation will keep the condition

number of the resulting discrete system under control [3, 14, 16].

4 Convergence and a priori error estimates

The H

− and L

-spaces are particularly suitable as the function spaces X and Y in

least-squares ﬁnite or spectral element methods resulting from the minimization of

ﬁrst order partial differential equations with strongly imposed boundary conditions.

To appreciate this, assume that we have the following norm-equivalence

kuk

(

Ω

≤kL uk

(

Ω

)

≤Ckuk

(

Ω

)

, ∀u ∈X =



u ∈ H

(

Ω

) | Ru = 0 on

∂Ω



(41)

where the space X represents the space of functions which already satisfy the homo-

geneous boundary condition and for which the function itself and its ﬁrst derivatives

are square integrable over the domain

Ω

. Based on the norm-equivalence (41), the

quadratic least-squares functional can be obtained which upon minimization yields

the weak formulation of the least-squares problem. If one uses a conforming ﬁnite-

187

Marc Gerritsma and Bart De Maerschalck

dimensional subspace X

⊂ X = H

(

Ω

), then one can approximate u

∈ X

piecewise continuously differentiable polynomials. In [68], it has been shown that

least-squares methods based on the norm-equivalence (41) yield the following error

estimate:



u −u



≤

C inf

∈X



u −v



, (42)

where the constant

C is given by

C = 1 + 2C

, with the constants

and C from

(41). Here the subscript 1 in the norm k·k

refers to the Sobolev norm k·k

(

Ω

)

while k·k

will denote k·k

(

Ω

)

= k·k

(

Ω

)

Since the interpolation of the solution u in the space X

obviously belongs to the

ﬁnite-dimensional subspace X

, the right-hand side of equation (42) can be bounded

in the following way:

inf

∈X



u −v



≤



u −

(u)



, (43)

where

(u) ∈ X

represents the interpolation of the solution u ∈ X in the space X

consisting of polynomials of degree N. The right-hand side of expression (43) can

be further bounded from above if an h− or p−reﬁnement strategy is used to obtain

a better approximation of the solution u ∈X. Since both reﬁnement strategies result

in different estimates for the interpolation error (43) and hence in a different error

estimate (42), they are discussed separately hereafter.

u ∈ H

(

Ω

) for

some s ≥ 2



u −u



≤Ch

|u|

l+1

, (44)

where l = min(N,s−1) and where N represents the approximating order of the C

−

spectral elements and where h represents a grid parameter (for example, the maxi-

mum of the square root of the area of the spectral elements), which decreases with

increasing number of spectral elements. Here |·|

denotes the semi-norm deﬁned by

|v|

∑

|=s

(

Ω

)

, (45)

where

is a multi-index given by the vector

= (

,.. .,

) , (46)

where the

are non-negative integers and

| =

+ ···+

(47)

We used the following notation for partial derivatives

188

Combining expression (43) with the interpolation error corresponding to h−

refinement, results in the following estimate if the exact solution

Least-Squares Spectral Element Methods in Computational Fluid Dynamics

∂

...

∂

. (48)

For example, if d = 2 and

= (1,2), then

u =

∂

. (49)

The rate of convergence (44) is optimal in the H

−norm since it provides the high-

est possible rate of convergence allowed by polynomials of degree N, [73]. Note

that the optimal rate of convergence depends on the polynomial degree (N) and the

regularity of the exact solution, denoted by the Sobolev exponent s. A similar rate of

convergence can be obtained for the L

−norm. In [73], it is shown that the following

−error estimate holds if the exact solution u ∈ H

(

Ω

) for some s ≥ 2:



u −u



≤Ch

l+1

|u|

l+1

. (50)

In case of p−reﬁnement, the order N of the spectral elements is increased while

keeping the number of the spectral elements constant. Consequently, the grid pa-

rameter h remains constant. As a consequence, the convergence rates (44) and (50)

are not suitable in this case. In order to obtain useful convergence rates in case of

p−reﬁnement, one can use the Legendre interpolation operator to bound the inter-

polation error u −

(u) from above. Assuming that u ∈ H

(

Ω

) for some s ≥ 2, it

can be proven [73, page 126] that



u −

(u)



≤CN

k−s

|u|

, with k = 0, 1, (51)

where

(u) represents the Legendre interpolation operator applied to the exact

solution u. Combining the latter expression with (42) and (43) results into the fol-

lowing error estimate



u −u



≤CN

k−s

|u|

, with k = 0, 1. (52)

Note that the rate of convergence is only bounded by the smoothness degree of the

solution s and not by any other grid parameter h as it occurs in the ﬁnite element

case. As a consequence, exponential convergence rates can only be obtained for

smooth problems if a p-reﬁnement strategy is used. Since this rate of convergence

results from p-reﬁnement, it is called p-convergence hereafter.

Least-squares formulations based on the norm-equivalence (41) are called fully

-coercive formulations and they have optimal h−convergence rates in the H

and L

-norms. Moreover, if the fully H

-coercive least-squares method is supple-

mented with strongly imposed boundary conditions, one can use standard ﬁnite and

spectral elements. Consequently, it is not so surprising that fully H

-coercive for-

mulations with strongly imposed boundary conditions are the preferred setting for

189

Marc Gerritsma and Bart De Maerschalck

least-squares ﬁnite and spectral element methods since these methods are both prac-

tical and optimally accurate.

Note that the above a priori estimates only depend on continuity and coercivity of

the differential operator L and interpolation estimates. In [69, 70] these error esti-

mates have been conﬁrmed numerically with the use of Legendre polynomials. But

since the estimates are independent of the type of orthogonal basis functions used,

similar results will hold for polynomials such as Chebyshev polynomials. These er-

ror estimates will play a role in hp-adaptive schemes to be discussed in Section 7.1.

Unfortunately, not all differential equations can be converted into H

-coercive

least-squares formulations. Especially compressible, inviscid ﬂows which exhibit

shocks or contact discontinuities and models like Burgers’ equation do not ﬁt in the

-coercive framework and therefore require a different treatment to be discussed

in Section 6.

5 Incompressible ﬂows

For incompressible ﬂows we can distinguish three regimes, the case where the

Reynolds number is so low that the convective terms can be neglected (Stokes ﬂow),

the case where the Reynolds number is sufﬁciently small to yield steady solutions

and for slightly higher Reynolds numbers, we have time-dependent ﬂow.

Stokes ﬂow has been thoroughly discussed by Proot, [69, 70]. All function spaces

for which norm-equivalence can be established have been identiﬁed in these papers.

In this chapter we will show some applications to steady and unsteady incompress-

ible, viscous ﬂow problems.

5.1 Governing equations

In terms of the primitive variables (u, p) the governing equations read

∂

+ (u ·∇)u = −∇p +

∆

u+ f in

Ω

, (53)

∇·u = 0 in

Ω

, (54)

where u represents the velocity vector, p the kinematic pressure, f the forcing term

per unit mass (if applicable) and Re the Reynolds number.

190

Least-Squares Spectral Element Methods in Computational Fluid Dynamics

5.2 The ﬁrst order formulation of the Navier-Stokes equations

In order to obtain an equivalent ﬁrst order formulation of the unsteady Navier-Stokes

equations, the vorticity

has been introduced as an auxiliary variable. By using the

identity ∇×∇×u = −

∆

u+ ∇(∇·u) and by using the incompressibility constraint

∇·u = 0, the governing equations subsequently read

∂

+ u ·∇u + ∇p +

∇×

= f in

Ω

, (55)

−∇ ×u = 0 in

Ω

, (56)

∇·u = 0 in

Ω

, (57)

where, in the particular case of the two-dimensional problem, u

= [u

] repre-

sents the velocity vector, p the kinematic pressure, f

= [ f

, f

] the forcing term per

unit mass (if applicable) and Re is the Reynolds number.

5.3 Linearization of the non-linear terms

Before the least-squares principles can be applied and the corresponding weak form

discretized with spectral elements, the convective term u·∇u must be linearized. To

this end, one can use a Picard (e.g. successive substitution)

u·∇u ≈ u

·∇u, (58)

or a Newton linearization

u·∇u ≈ u

·∇u + u ·∇u

−u

·∇u

. (59)

In the latter two equations, the subscript “0” indicates that the value of the corre-

sponding variable is known from the previous iteration step.

If Picard linearization is used the following linearized momentum equation is

obtained:

∂

+ (u

·∇u + ∇p +

∇×

−f) = 0 . (60)

For steady problems the time derivative can be dropped. For time-dependent calcu-

lations either a time-stepping scheme needs to be selected or space-time elements

can be used, [19, 20, 21, 26, 63].

The least-squares formulation now becomes: Find u, p,

∈ H

(

Ω

) which mini-

mizes the functional in the absence of body forces

I (u, p,

) = k∇·uk

+ k

−∇ ×uk

191

Marc Gerritsma and Bart De Maerschalck



∂

+ (u

·∇u + ∇p +

∇×

)



. (61)

Variational analysis with respect to the four unknowns u, p and

leads to the sym-

metric positive deﬁnite system.

5.4 Steady ﬂow around a cylinder at low Reynolds

One of the test cases considered for LSQSEM is the steady ﬂow around a circular

cylinder. This research was performed by De Groot, [41], with the conventional

least-squares method discussed above and Direct Minimization, see Section 7.2.

Here we are simulating the ﬂow around a circular cylinder placed perpendicularly

in a uniform parallel ﬂow, see Fig. 1.

Fig. 1 Description of the ﬂow problem

It is known that for ﬂows with Reynolds numbers Re < 45 a steady solution exists

and for Reynolds number higher than Re = 45 time-dependent behaviour sets in and

the well-known Von Karman vortex street develops. The investigation of the ﬂow

will be in the Reynolds range: 1 < Re < 45. In this range two different ﬂow types

are observed. For Re < 6 the ﬂow is completely attached to the cylinder as can be

seen in Fig. 2 (a). For Re > 6 the ﬂow starts to separate somewhere on the cylinder

forming two attached eddies behind the cylinder, Fig. 2 (b). For a benchmark be-

tween LSQSEM on the one hand and experimental data on the other, the estimation

of the Re number for which the twin vortices appear, is an important property of the

ﬂow. This Re number is from here on referred to as Re

onset

. The twin vortices grow

in length and change shape as Re increases. The accurate prediction of this relation

is a common benchmark. Fig. 3 gives a schematic overview of the ﬂow where two

steady vortices exist behind the cylinder.

192