Spectral properties of the incompressible Navier-Stokes equations

(1)

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Spectral

properties

of

the

incompressible

Navier-Stokes

equations

Fredrik Laurén

a

,

∗

,

Jan Nordström

a

,

b

a_Department_of_Mathematics,_{Computational}_Mathematics,_Linköping_University,_SE-581₈₃_Linköping,_Sweden

b_Department_of_Mathematics_and_Applied_Mathematics,_University_of_{Johannesburg,}_P.O._Box_524,_Auckland_Park_2006,_South_Africa

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received13August2020

Receivedinrevisedform17November2020 Accepted18November2020

Availableonline20November2020 Keywords:

IncompressibleNavier-Stokes Convergencetosteady-state Openboundaryconditions Dispersionrelation High-orderaccuracy Summation-by-parts

Theinfluenceofdifferentboundaryconditionsonthespectralpropertiesofthe incompres-sible Navier-Stokes equations is investigated. By using the Fourier-Laplace transform technique,wedeterminethe spectra,extract thedecayrateintime, and investigatethe dispersion relation. In contrast to an infinite domain, where only diffusion affects the convergence, weshowthatalsothepropagationspeedinfluencetherateofconvergence to steady state for a bounded domain. Once the continuous equations are analyzed, we discretize using high-order finite-difference operators on summation-by-parts form and demonstrate that the continuous analysis carries over to the discrete setting. The theoreticalresultsareverifiedbynumericalexperiments,wherewehighlightthenecessity ofhighaccuracyforacorrectdescriptionoftime-dependentphenomena.

1. Introduction

In most fluid dynamic simulations, the numerical domain is truncated with an artificial open boundary which requires boundary conditions. For the incompressible Navier-Stokes (INS) equations, several types of (more or less) non-reflecting boundary conditions have been developed to avoid non-physical effects caused by lack of appropriate data at the artificial boundary. The most common choice is the so-called natural (or “do-nothing”) condition [42], but there are many more, see for example [19,20,8,11,38].

Solvers that rapidly converge to steady state are important in many applications. Different numerical techniques such as local time-stepping [41,39,18] and multigrid [44,45,15,7] are often used to enhance the convergence rate. For truly unsteady problems, the eﬃciency of the time-integration procedure is even more important. A fundamental requirement of the un-derlying initial-boundary value problem (IBVP) is that its spectrum resides in the appropriate half of the complex plane [9,26,30]. Otherwise most time-integration and convergence acceleration techniques fail.

We will investigate the inﬂuence of different sets of open boundary conditions on the spectral properties of the INS equations and compare the results to solutions on an inﬁnite domain. Both a von Neumann and normal mode analysis [6,13] will be used to compute the spectrum, and the difference between these approaches is discussed. The normal mode analysis was used in [9,23] for hyperbolic problems and in [26,30,38] for the Navier-Stokes equations. Besides analyzing the decay rate in time, as is done in previous studies [9,23,26,30,38], we will also include the effect on the time-dependent behavior characterized by the dispersion relation. This is used to highlight the effect of a bounded domain on the convergence rate to steady state.

*

Correspondingauthor.

E-mailaddress:fredrik.lauren@liu.se(F. Laurén).

https://doi.org/10.1016/j.jcp.2020.110019

(2)

Once the continuous problem is analyzed, we proceed to the discrete approximation. High-order ﬁnite-difference opera-tors on summation-by-parts (SBP) form are used for the discretization [48,47]. The boundary conditions are imposed weakly using the Simultaneous-Approximation-Term (SAT) technique [22,4]. By using the SBP-SAT framework, the discrete approxi-mation mimics its continuous analogue and produce spectra that converge to the continuous ones. Additionally, we will use the normal mode analysis [13,21] to show that high-order accuracy is necessary for a correct and eﬃcient description of time-dependent problems.

The paper proceeds as follows. In Section 2, the IBVP is introduced. We obtain energy estimates by using the energy method and appropriate boundary conditions. The spectrum of the continuous problem is derived in Section 3, where the decay and growth in time and space as well as the dispersion relation are investigated. In Section 4, the equations are discretized in space by using the SBP-SAT technique and we derive discrete energy estimates that mimics the continuous ones. Furthermore, the discrete spectrum is computed and compared to its continuous analogue. The fully discrete scheme and numerical experiments are presented in Section 5, where the predictions from the continuous analysis are verified. We also illustrate the influence of bounded and infinite domains on steady-state computations. Conclusions are drawn in Section6.

2. Thecontinuousproblem

To produce well-deﬁned spectra, we assume periodicity of the solution in one of the spatial dimensions. Consider the IBVP for the INS equations

˜

I Ut

+

AUx

+

BUy

=

˜

I

(

Uxx

+

Uy y

),

(

x

,

y

)

∈ ,

t

≥

0

,

U

=

f

(

x

,

y

),

(

x

,

y

)

∈ ,

t

=

0

,

L0U

=

g0

(

t

,

x

,

y

),

(

x

,

y

)

∈

0

,

t

≥

0

,

L1U

=

g1

(

t

,

x

,

y

),

(

x

,

y

)

∈

1

,

t

≥

0

,

(1)

posed on the domain

= [

0

,

1

]

× [−

π

,

π

]

. In (1),

>

0 is the viscosity coeﬃcient and U

= (

u

,

v

,

p

)

is the solution vector, where u

,

v are

the velocities in the

x andy direction,

respectively, and

p is

the pressure. The matrices in (

1) are

A

=

⎛

⎝

u0 0u 10 1 0 0

⎞

⎠ ,

B

=

⎛

⎝

v0 0v 01 0 1 0

⎞

⎠ , ˜

I

=

⎛

⎝

1 0 00 1 0 0 0 0

⎞

⎠

and f , g0, g1 are smooth functions that are periodic in the y

−

direction. We denote the left and right boundaries by

0

= {(

x

,

y

)

∈

:

x

=

0

}

and

1

= {(

x

,

y

)

∈

:

x

=

1

}

, respectively. At the open boundaries

0 and

1, data is imposed through the boundary operators L0 and L1. To simplify the setup, we postulate that u

>

0, with

0as the inﬂow boundary and

1 as the outﬂow boundary, see Fig.1.

2.1. Anaprioriestimate

The energy method [13,14,21] will be used to derive estimates for the norm of the solution. For two vector functions U and V deﬁned

on

, we introduce the inner product and norm

U

,

V

=

ˆ

UV d

,

U

2

=

U

,

U

.

(3)

d dt

U

2 ˜ I

+

2

∇

U

2 ˜ I

=

B T

.

(2) In (2),

∇

U

2 ˜ I

=

´

(˜

I Ux

,

˜

I Uy

)

(

Ux

,

Uy

)

d

and the boundary terms are B T

=

ˆ

0 U

(

AU

−

2

˜

I Ux

)

dy

−

ˆ

1 U

(

AU

−

2

˜

I Ux

)

dy

,

₍₃₎

where periodicity have removed the boundary contribution at y

= ±

π

. We can now state the following result.

Proposition1.Consider

(

1) with homogeneousboundaryconditions.IfB T

≤

0 in (3), then

U

2_˜

I isbounded. Proof. If B T

≤

0, integrating (2) in time yields directly

U

2_˜

I

≤

f

2 ˜ I.

Remark1.

In the upcoming analysis, the linearized version of (

1) is considered, where A

=

A

(

u

¯

),

B

=

B

(

v

¯

),

L0

=

L0

( ¯

U

)

and

L1

=

L1

( ¯

U

)

. This formulation covers the Oseen [40] and Stokes [46] equations. We will in the following view the solution

W

= (

u

,

v

,

p

)

as a perturbation from the constant mean U

¯

= (¯

u

,

v

¯

,

p

¯

)

, such that U

= ¯

U

+

W .

3. Thecontinuousspectrum

One of the objectives with this paper is to illustrate the difference between initial value problems (IVPs), with periodic solutions (or equivalently Cauchy problems), and IBVPs, where boundary conditions are included.

3.1. VonNeumannanalysis

We consider the fully periodic IVP on the domain

[−

π

,

π

]

2_{and apply Fourier expansions in both the}_{x and} _y-directions to get W

(

t

,

x

,

y

)

=

∞

κ,ω=−∞ V

(

t

,

κ

,

ω

)

eiκxeiωy

,

(4)

where κ and ωare integers. Together with (4), the ﬁrst equation in (1) becomes

˜

I Vt

+ (

i

κ

A

+

i

ω

B

+

I

˜

(

κ

2

+

ω

)

2

)

V

=

0

.

(5) By inserting the ansatz V

= ψ

est, we ﬁnd that the Laplace transformed version of (5) with zero initial condition satisﬁes

[

s

˜

I

+

i

κ

A

+

i

ω

B

+

˜

I

(

κ

2

+

ω

2

)

]ψ =

0

.

(6)

For non-trivial solutions

ψ

, it is required that det

(

s

˜

I

+

i

κ

A

+

i

ω

B

+

˜

I

(

κ

2

+

ω

2

))

=

0

,

which yields

s∗

= −

i

(

u

¯

κ

+ ¯

v

ω

)

−

(

κ

2

+

ω

2

) .

(7)

The set of s∗in (7), that allows for a non-trivial solution, is the spectrum of the IVP. Next, we use (7) to rewrite each mode in (4) as

˜

W

= ψ

e

−(κ

2+ω2)

t decay eiκ(x−¯ut)eiω(y−¯vt)

propagation

.

(8)

Remark2.

Several observations can be made by considering (

8):

•

There is no decay or growth in space.

•

Re

(

s∗

)

= −

(

κ

2

₊

_ω

2

₎

_{in (}₇_{) determines the decay in time.}

•

All modes propagate and oscillate as eiκ(x−¯ut)_eiω(y−¯vt) _and_hence_advect_with_the_same_velocity₍_u

_¯

_,

_v

_¯

₎

_(except_for

κ

,

ω

=

0).

For the periodic case, the analysis is straightforward and the results above are known. However, as we shall see in the next section, including boundaries leads to a signiﬁcantly more complex situation.

(4)

det

(

[(

s

+ (

ω

)

−

κ

)˜

I

+

i

(

ω

)

B

+

κ

A

]) = [˜

s

+

κ

u

¯

−

κ

][(

ω

)

−

κ

] =

0

,

must hold, where s

˜

=

s

+

i

(

ω

)

v

¯

+ (

ω

)

2. Hence,

κ

1,2

= ±|

ω

|,

κ

3,4

=

1 2

¯

u

±

u

¯

2

₊

₄

_˜

_s

_.

₍₁₀₎

Unlike the fully periodic solution in (4), we now have κ

=

κ

(

s

,

ω

)

∈ C

. In the absence of coinciding roots, each mode can be written

˜

W

(

s

,

x

,

ω

)

=

4

j=1

σ

j

ψ

jest+κjx/

=

X

σ

est

,

(11)

where X

=

diag

(

eκ1x/

_,

_eκ2x/

_,

_eκ3x/

_,

_eκ4x/

₎

_{. The corresponding eigenvectors}

_ψ

j to every κj, j

=

1

,

2

,

3

,

4, are stored in the matrix

given by

=

⎡

⎣

ψ

|

1

, ψ

|

2

, ψ

|

3

, ψ

|

4

,

|

⎤

⎦ =

⎡

⎣

⎛

⎝

|

i

ω

|

−

γ

+

⎞

⎠ ,

⎛

⎝

−

|

ω

i

ω

|

γ

−

⎞

⎠ ,

⎛

⎝

ω

i

κ

3 0

⎞

⎠ ,

⎛

⎝

ω

i

κ

4 0

⎞

⎠

⎤

⎦ ,

where γ±

=

s

+

i

(

ω

)

v

¯

±|

ω

|¯

u.

Also,

σ

= (

σ

1

,

σ

2

,

σ

3

,

σ

4

)

are weights that will be determined by the boundary conditions. Remark3.

Since there are four unknown weights in (

11), four boundary conditions are necessary. Too few conditions result in a non-unique solution, while too many affect the existence.

Inserting (11) into the homogeneous form of the boundary conditions yields

E

(

s

)

σ

=

L0

X

(

0

)

L1

X

(

1

)

σ

=

0

.

(12)

In (12), E

(

s

)

is a 4

×

4 matrix with complex entries (two boundary conditions at the inﬂow boundary and two at the outﬂow boundary are imposed).

A non-trivial solution to (12) only exists if E

(

s

)

is singular. Hence, the spectrum of the IBVP [13] is given by

s∗

= {

s

∈ C :

det

(

E

(

s

))

=

0

} .

(13)

Note the difference (caused by truncating the domain of the IVP and imposing boundary conditions) between the sets of s∗ that allow for non-trivial solutions in (7) and (13). By using that κj

= (

κ

R

)

j

+

i

(

κ

I

)

j, j

=

1

,

2

,

3

,

4 and s

=

sR

+

isI, we can similarly to (8) rewrite (11) to reﬂect the growth and wave propagation characteristics as

˜

W

=

4

j=1

σ

j

ψ

j

esRt+(

κR)jx/

decay/growth ei(sIt+(κI)jx/)

propagation

.

(14)

Remark4.

Several observations can be made by considering (

14):

• κ

R determines the growth/decay in space and where the boundary conditions must be imposed [14]. In all cases, there are two growing and two decaying modes (see (10)). This will be discussed in Section3.4.

•

sR determines the growth or decay in time. If all eigenvalues reside in the left half of the complex plane, then the solution only consists of decaying modes in time, see Section3.5for details.

(5)

•

The factor ei(sIt+κIx/)_{in (}₁₄_{) is related to the propagation and oscillatory behavior of the solution, which we discuss in} Section3.6.

We realize from the bullets in Remark4and Remark2that the truncated domain complicates the situation and thereby makes time-dependent behavior of IBVPs signiﬁcantly more complex than for IVPs.

3.3. Differentsetsofopenboundaryconditions

It was proven in [33,34] that the homogeneous form of Set 1, Set 2 and Set 3 (presented below) satisfy Proposition 1

and therefore bound

U

2 ˜ I.

The ﬁrst set for the IBVP is obtained by specifying the so-called ingoing variables, see [33] for details. Explicitly, the boundary conditions are

Set 1:

λ

5u

+

p

−

ux

=

g1

,

λ

4v

−

vx

=

g2

,

at

0

,

λ

1u

+

p

−

ux

=

g3

,

λ

2v

−

vx

=

g4

,

at

1

,

where

λ

1,5

=

u 2

∓

u2 2

+

2

,

λ

2,4

=

u 2

∓

u2 2

+

1

.

The second set for the IBVP imposes Dirichlet conditions at the inﬂow boundary and natural boundary conditions at the outﬂow boundary:

Set 2:

u

=

g1

,

v

=

g2

,

at

0

,

p

−

ux

=

g3

,

−

vx

=

g4

,

at

1

.

The stabilized natural boundary conditions [1] is the third IBVP set. They have the same form both at inﬂow and outﬂow boundaries and are

Set 3: h

(

u

)

u

/

2

+

p

−

ux

=

g1

,

h

(

u

)

v

/

2

−

vx

=

g2

,

at

0and

1

,

where h

(

u

)

=

u

,

at

0

,

0

,

at

1

.

At

1, the stabilized natural boundary conditions reduce to the natural boundary conditions. For the IVP (the fully periodic case), the solution satisﬁes

Set 4: U

0

=

U

1

,

˜

I Ux

0

= ˜

I Ux

1

,

which are periodic conditions in the x-direction.

Set 4 yields

B T

=

0 and Proposition1is therefore satisﬁed. 3.4. Ill-posedboundaryconditions

From the normal mode analysis in Section3.2, the correct number and position of the boundary conditions are provided [32]. For Re

(

s

)

>

0, we have Re

(

κ

1

),

Re

(

κ

3

)

>

0 and Re

(

κ

2

),

Re

(

κ

4

)

<

0, i.e. two modes are growing and two are decaying. Hence, two energy-stable boundary conditions must be imposed at each boundary to properly bound the solution. To illus-trate this point, we present two different sets of boundary conditions that show the importance of i

)

the position where the boundary conditions are applied and ii

)

the f orm of

the boundary conditions.

The ﬁrst set is Set 5:

u

=

g1

,

at

0

,

v

=

g2

,

p

−

ux

=

g3

,

−

vx

=

g4

,

at

1

,

which is almost the same as Set 2, except that v

=

g2is imposed at

1 instead of

0. The second set is

Set 6:

p

−

u2

/

2

−

ux

=

g1

,

uv

/

2

−

vx

=

g2

,

at

0

,

p

−

ux

=

g3

,

−

vx

=

g4

,

at

1

.

The sign of the u2_{-term in Set 6 is shifted compared to the same term in Set 3. Imposing the homogeneous form of either} Set 5 or Set 6 do not result in B T

≤

0 and Proposition1is not satisﬁed.

(6)

Fig. 2. Eigenvaluesto(13) fortheboundaryconditionsinSet5andSet6. Thereareeigenvaluesintherighthalfofthecomplexplaneforbothsets, leadingtoexponentialgrowth.

Fig. 3. Theeigenvaluesfrom(7) andto(13).Set1,Set2andSet3 corre-spondtotheIBVP,whileSet4isthefullyperiodicIVP.

Table 1

The threeutmostright lying eigenval-uestotheIBVP.

Set 1 s∗₀ −3.92±2.43i s∗1 −7.33±7.80i s∗2 −10.04±12.94i Set 2 s∗0 −6.74±3.12i s∗1 −8.09±8.81i s∗2 −10.15±14.19i Set 3 s∗₀ −0.6746 s∗₁ −9.425±3.151i s∗2 −11.08±8.765i

The Secant method (with initial values covering a large domain around the origin) is used to ﬁnd the singularities s∗ in (13). Unless otherwise speciﬁed, u

¯

=

1

,

v

¯

=

0, ω

=

1 and

=

0

.

01 are used in the computations below. Fig.2 shows the utmost right lying eigenvalues to (13) for Set 5 and Set 6. Since there are eigenvalues with Re

(

s

)

>

0 for both sets, exponentially growing solutions in time will be the result.

3.5. Thedecayorgrowthintime

We now shift focus to the sets of boundary conditions that satisfy Proposition1. To start, we apply the energy method to (9). All of the boundary conditions in Set 1, Set 2, Set 3 and Set 4 leads to

d dt

V

2 ˜ I

≤ −

2

Vx

2 ˜I

−

ω

2

_V

2 ˜ I

.

(15)

Clearly the decay rate in time is affected by both and ω.

Fig.3 shows the eigenvalues for the different sets of boundary conditions in the domain s

∈ [−

15

,

0

]

×

i

[−

20

,

20

]

. All cases reside in the left half of the complex plane. The three utmost right lying eigenvalues for Set 1, Set 2 and Set 3 are given in Table1.

In Fig.4, the absolute value of the real part of the utmost right lying eigenvalue of the spectrum (denoted by

|

η

∗_C

|

) is listed. The coeﬃcient η∗_C determines the continuous decay in time. For small values of , Set 2 is the most dissipative one. As ωand increase, the diffusion mechanism in the y-direction becomes more signiﬁcant and the effect of the boundary conditions at

0,1 on

|

η

∗C

|

is relatively reduced. For ω

=

0,

|

η

∗C

|

increases signiﬁcantly for Set 2 as is reduced, while it remain approximately constant for Set 1 and Set 3. These effects are due to the presence of coinciding roots when ω

=

0 and analyzed in Appendix A.

On an infinite domain, a larger diffusion coefficient leads to faster decay, see the second bullet in Remark2. However, this is not always the case for a finite domain. In Fig.5, ω

=

1

,

3

,

5 and is varied for Set 1. The combination of small ω and small indicates a relatively fast decay. As increases, the decay gets slower. This somewhat counter-intuitive result is most easily understood by considering the solution to the advection-diffusion equation ut

+

ux

=

uxx. By studying Fig.6, we

(7)

Fig. 4. The decay rate|η∗C|for differentandω. Set 1, Set 2 and Set 3 correspond to the IBVP, while Set 4 relate to the fully periodic IVP.

Fig. 5. Thecontinuousdecayrateintime|ηC∗|fordifferentandωfortheIBVPwithboundaryconditionsSet1.Therightﬁgureisaclose-upoftheleft

ﬁgureforω=1.

realize that the slowdown of the convergence is due to the fact that the solution remains longer inside the limited domain if is large since it becomes more smeared and hence elongated. This is essentially a one-dimensional effect. For two-dimensional problems, the diffusion mechanism also acts on the solution in the transverse spatial dimension and reduces or masks this effect. If ω is increased, the diffusive effect is increased and starts to dominate, see Fig.5. This phenomenon will be further discussed in Section5.4.

3.6. Thepropagationcharacteristics

The exponent of the factor ei(sIt+κIx/)_{in (}₁₄_{) is the phase [}₅₀_{]. Thus,}

vp

=

dx dt

= −

sI

κ

I

,

κ

I

=

0

,

(16)

(8)

Fig. 6. Solutiontotheequationut+ux=uxxfordifferentvalueson.Thesolutionbecomesmoreelongatedasisincreased.Sincethedomainisﬁnite,

thepropagationcharacteristicsofthesolutioninﬂuencethedecayrate.

Fig. 7. The spectrum for Set 1 asu¯→0. The right panel is a close-up of the left panel.

Fig. 8. Thevaluesκ3(leftﬁgure)andκ4(rightﬁgure)fortheIBVPwiththeboundaryconditionsinSet1correspondingtotheutmostrightlyingrootsof thespectrum.

is the phase velocity, which is related to the advective speed of the solution.

Fig.7illustrates the spectrum for Set 1 as u

¯

→

0. The eigenvalues approach the real line (i.e.

|

s∗_I

|

→

0) as u decreases

¯

(this is in agreement with the fully periodic case, since

|

s∗_I

|

→

0 in (7) as u

¯

,

v

¯

→

0). In Fig.8, the corresponding values of

κ

3,4 are illustrated for Set 1 as a function of u.

¯

Since

|

s∗I

|

decreases, but the imaginary parts of κ3,4 are roughly unchanged,

|

vp

|

decreases and the solution propagates slower as u

¯

→

0. For u

¯

=

0, (1) becomes the Stokes equations [46] and the spectrum becomes purely real (s∗_I

=

0), which result in vp

=

0. Thus, for the Stokes equations, the only change in time is due to diffusion. There is no propagation mechanism.

(9)

Table 2

Thephasevelocities for the IBVP (computedfrom (16)) ofthe modescorrespondingtotheeigenvalueinTable1.

|vp| = |dx/dt| = |s∗I/κ3,I| = |s∗I/κ4,I| s∗ u¯=1 u¯=0.5 u¯=0.25 Set 1 s∗0 0.920 0.436 0.200 s∗1 0.860 0.398 0.179 s∗2 0.817 0.373 0.162 Set 2 s∗0 0.858 0.377 0.140 s∗₁ 0.848 0.364 0.133 s∗₂ 0.842 0.358 0.127 Set 3 s∗0 0 0 0 s∗1 0.794 0.307 0.0748 s∗2 0.780 0.296 0.0740

Table2lists the phase velocities (computed by using (16)) corresponding to the three utmost right lying eigenvalues for

¯

u

=

1

,

0

.

5

,

0

.

25. As expected, the speed decreases with u (one

¯

of the

rows in Table2 for Set 3 only contains zeros since there is one purely real root). We can conjecture that the choice of boundary conditions inﬂuence the convergence rate to steady state in two ways since it affects both the

decay in time and the propagation characteristics.

Remark5.

A

major difference compared to the fully periodic case (8) is that not all modes have the same propagation velocity, see Remark2.

3.6.1. Time-varyingmodes

For a given

ˆ

s

∈

s∗, det

(

E

(

s

ˆ

))

=

0 and E

(

s

ˆ

)

has a non-trivial nullspace. Hence,

σ

=

0 in (11) can be obtained from the nullspace of E

(

ˆ

s

)

. With all parameters (U

¯

,

ω

,

) prescribed and

ˆ

s ﬁxed,

only one free parameter remains in (

11) for speciﬁc values of t

,

x.

By inserting the eigenvalues displayed in Table

1into (11), explicit expressions for the modes generating the solution can be obtained. In the ﬁgures below, the free parameter in

σ

is determined such that max

| ˜

W

(

ˆ

s

,

x

,

ω

,

0

)

|

=

1.

Fig.9shows W

˜

(

ˆ

s

,

x

,

ω

,

t

)

corresponding to Set 1 for ω

=

1 at the initial time. In Fig.10, the time evolution of V

(

ˆ

s

,

x

,

ω

,

t

)

is illustrated for the eigenvalue

ˆ

s

= −

7

.

33

+

7

.

80i at t

=

0

,

0

.

1

,

0

.

2. The variables uand vboth decrease in time and advects in space. Since the third components (which correspond to p) of

ψ

3,4 is zero, that variable is only affected by the purely real roots κ1 and κ2. By studying (14), we realize that p does not advect since κI

=

0. As u decreases,

¯

we expect the modes to propagate slower (see Table2). In Fig.11, the time evolution of the imaginary part of vcorresponding to the ﬁrst mode for different values of u is

¯

illustrated. In agreement with the prediction, the mode propagates slower as

u decreases.

¯

4. Thediscreteproblem

To make the paper self-contained, we provide a brief introduction of the finite-difference version of the SBP-SAT technique and recommend [48,10] for extensive reviews. Several procedures such as finite-volume [30,31,28,49], spectral elements [2,3], flux reconstruction [5,17,43] and discontinuous Galerkin schemes [12,16] are included in the SBP-SAT frame-work, and the upcoming results holds with small technical modifications also for those schemes.

4.1. TheSBP-operators

Let the domain

= [

0

,

1

]

× [−

π

,

π

]

be discretized with N

+

1 and M

+

1 grid points; xi

=

i

/

N, i

=

0

,

. . . ,

N and

yj

= −

π

+

j

/

M, j

=

0

,

. . . ,

M.

The discrete approximation of the variable

u

=

u

(

x

,

y

)

, is a vector of length

(

N

+

1

)

× (

M

+

1

)

arranged as u

= (

u00

. . .

u0M

,

. . . ,

uN0

. . .

uN M

)

. Hence, the approximation of the solution vector U

= (

u

,

v

,

p

)

has the structure U

= (

u

,

v

,

p

)

.

The approximation of the spatial derivatives are given by

DxU

= (

I3

⊗

Dx

)

U

≈

Ux

,

DyU

= (

I3

⊗

Dy

)

D

≈

Uy

,

where

⊗

denotes the Kronecker product, Dx

= (

Px−1Qx

)

⊗

Iy and Dy

=

Ix

⊗ (

P−y1Qy

)

denote the difference operators on SBP form and Ux,y contain the analytical partial derivatives evaluated on the grid. The matrices Px,y are diagonal and positive deﬁnite, such that P

=

I3

⊗

Px

⊗

Py forms a quadrature rule that deﬁnes the norm

U

2P

=

UPU

≈

˜

UU d

. Furthermore, the matrices Qx,ysatisfy the SBP-property

Qx

+

Qx

=

EN

−

E0

,

Qy

+

Qy

=

EM

−

E0

,

(17) where E0

=

diag

(

1

,

0

,

0

,

. . . ,

0

)

and EN,M

=

diag

(

0

,

0

,

0

,

. . . ,

1

)

are matrices of appropriate sizes.

(10)

(11)

Fig. 10. Timeevolutionofthemodecorrespondingtos∗= −7.33+7.80i.Thearrowsindicatethedirectionofmovement.Notethatthepdoesnotadvect. Differentscalesareusedineachﬁgure.

(12)

Fig. 11. TimeevolutionoftheimaginarypartofvcorrespondingtothemostrightmodeoftheIBVPwiththeboundaryconditionsinSet1foru¯=1 and ¯

u=0.25.

4.2. Semi-discretenonlinearstability

The semi-discrete version of (1) with homogeneous boundary conditions is

˜

IUt

+

Dh

(

U

)

=

0

,

t

>

0

,

U

=

f

,

t

=

0

,

(18) where Dh

(

U

)

=

1 2

(

DxA

+

ADx

)

U

+

1 2

DyB

+

BDy

U

−

D2_x

+

D2_y

˜

IU

−

S ATx

−

S ATy

.

(19) In (19), we have introduced A

=

⎛

⎝

d

(

0u

)

d

(

0u

)

0I I 0 0

⎞

⎠ ,

B

=

⎛

⎝

d

(

0v

)

d

(

0v

)

0I 0 I 0

⎞

⎠ , ˜

I

=

⎛

⎝

0I 0 0I 0 0 0 0

⎞

⎠

where d

(

u

),

d

(

v

)

are diagonal matrices with u

,

v on

the diagonal, respectively. Moreover,

0

,

I are

the zero and the identity

matrices of appropriate sizes. In (18), the skew-symmetric form of the INS equations have been discretized, which is nec-essary for obtaining an energy estimate by using the energy method [27,33]. The terms S ATx and S ATy are penalty terms that impose boundary conditions weakly via the SAT-technique in the x- and y-direction,

respectively.

Multiplying (18) by 2UP from

the left and using the SBP-property (

17) of the operators Dx,y yield d dt

U

2 ˜ IP

+

2

∇

U

2 ˜ IP

=

BT

,

(20) where

∇

U

2 ˜ IP

= (

DxU

)

_(˜

_IP

₎₍

_D_x_U

₎

_{+ (}

_D_y_U

₎

_(˜

_IP

₎₍

_D_y_U

₎

_._{The term}_{S AT}_y _{imposes periodic boundary conditions that cancels}

all boundary terms at y

= ±

π

[29]. Hence, BT in

(

20) denotes all terms evaluated at the west and east boundaries:

BT

=

UP0

(

AU

−

2

˜

IDxU

)

−

UP1

(

AU

−

2

˜

IDxU

)

+

2UPS ATx

.

(21) In (21), P0,1

=

I3

⊗

P0,1 and P0,1

=

E0,N

⊗

Py are the quadratures corresponding to the west and east sides of the domain. We can now formulate the discrete version of Proposition1.

Proposition2.Considerthesemi-discreteapproximation

(

18). IfBT

≤

0 in (21), then

U

2_˜

IPisbounded. Proof. IfBT

≤

0, integrating (20) in time yields directly

U

2

˜ IP

≤

f

2 ˜ IP.

To illustrate the SBP-SAT procedure, we consider the boundary conditions in Set 2. Generally, for the boundary condition H U

=

g at

boundary

k (k can

be any of the sides of

), the SAT-term is of the form

(13)

S AT

=

P−1

kPk

(

HU

−

g

) ,

where

k is a penalty matrix to be determined for stability and Pk is the boundary quadrature. Speciﬁcally for Set 2,

S ATx

=

P−1

0P0

⎛

⎝

uv

−

gg21 0

⎞

⎠

HU−g

+

P−1

1P1

⎛

⎜

⎝

p

−

Dxu

−

g3

−

Dxv

−

g4 0

⎞

⎟

⎠

HU−g (22)

where P0,1

=

I3

⊗

P0,1 and P0,1

=

E0,N

⊗

Py. The penalty matrices are

0

=

⎛

⎝

−

d

(

u

)/

2

+

D x 0 0 0

−

d

(

u

)/

2

+

D_x 0

−

I 0 0

⎞

⎠ ,

1

= (

I3

⊗

I

) .

For the homogeneous form (i.e. g1

=

g2

=

g3

=

g4

=

0), multiplying (22) by 2UP results

in

2UPS ATx

= −

uP0

[

d

(

u

)

u

+

d

(

v

)

v

+

2p

−

2

Dxu

] +

2

vP0Dxv

+

2uP1

[

p

−

Dxu

] −

2

vP1Dxv

.

Hence, (21) becomes BT

= −

uP1

[

d

(

u

)

u

+

d

(

v

)

v

]

≥0

≤

0

.

(23) The last inequality holds since the right side is an outﬂow boundary (i.e. u

>

0) and Proposition 2 is satisﬁed for the homogeneous form of Set 2.

Remark6.

By imposing the homogeneous form of Set 2 into (

3), we get

B T

= −

ˆ

1

u

(

u2

+

v2

)

dy

≤

0

.

(24)

The inequality holds by the same argument as for (23). Hence, we see that (23) is the discrete analogue of (24).

It was proven in [33] that, with a correct SAT procedure, the imposition of the homogeneous form of Set 1, Set 2, Set 3 and Set 4 satisfy Proposition2and therefore bound

U

2

˜

IP. We do not obtain estimates for Set 5 and Set 6. 4.3. Thediscretespectrum

To compute the discrete spectrum, we discretize the continuous domain

[

0

,

1

]

with N

+

1 grid points: xi

=

i

/

N, i

=

0

,

. . . ,

N and

consider the semi-discrete version of the Fourier-transformed problem (

9):

(˜

I

⊗

Ix

)

Vt

+

MV

=

0

,

(25)

where

M

=

A

⊗

Dx

+ (

i

ω

B

+

ω

2

˜

I

)

⊗

Ix

−

˜

I

⊗

D2x

−

S ATx (26) and Ix is the identity matrix of size N

+

1. The penalty term S ATxis determined from the analysis in the previous section. By inserting the ansatz V

= ψ

est we ﬁnd

[

s

(˜

I

⊗

Ix

)

+

M

]ψ =

0

.

(27)

The discrete spectrum is deﬁned as s∗_D

= {

s

∈ C :

det

[

s

(˜

I

⊗

Ix

)

+

M

]

=

0

}

.

In the computations below, we compare SBP21 and SBP84 [48]. The operator SBP21 is second order accurate in the interior and ﬁrst at the boundaries, while SBP84 is eighth order accurate in the interior and fourth order accurate at the boundaries. Fig.12 illustrates the discrete spectrum for the boundary conditions in Set 1 for ω

=

1 and

=

0

.

01. As the number of grid points increases, the continuous spectrum is increasingly better approximated. Similar results are obtained for all sets of boundary conditions, including Set 5 and Set 6. Hence, if the continuous spectrum has eigenvalues in the right half of the complex plane, then a consistent numerical scheme will have that as well.

(14)

Fig. 12. Comparison between continuous and discrete spectra for Set 1. Note the poor convergence for SBP21 compared to SBP84.

Table 3

Comparisonofthethreeutmostrighteigenvaluesforthecontinuousanddiscrete spectraforallsetsofboundaryconditions.TheoperatorsSBP21andSBP84have beenusedforthespatialdiscretizationonthegridN=60.

s∗ Continuous SBP21 SBP84 Set 1 s∗₀ −3.92±2.43i −3.96±2.32i −3.92±2.43i

s∗1 −7.33±7.80i −5.17±3.51i −7.35±7.81i s∗2 −10.04±12.94i −5.22±6.93i −10.05±12.95i Set 2 s∗0 −6.74±3.12i −4.48 −6.74±3.19i

s∗1 −8.09±8.81i −4.34±2.79i −8.09±8.81i s∗₂ −10.15±14.19i −4.43±5.60i −10.14±14.19i Set 3 s∗₀ −0.675 −0.675 −0.675

s∗1 −9.43±3.151i −4.65±2.75i −9.43±3.15i s∗2 −11.08±8.765i −4.68 −11.08±8.76 Set 4 s∗0 −0.010 −0.010 −0.010

s∗1 −0.020±1.00i −0.020±0.99i −0.020±1.00i s∗₂ −0.050±2.00i −0.048±1.95i −0.050±2.00i Set 5 s∗₀ 4.63±9.12i 3.32±7.28i 4.72±9.34i

s∗1 4.20±16.85i 2.18±9,58i 3.42±15.00i s∗2 2.72 2.69 2.72 Set 6 s∗0 1.31 1.32 1.31

s∗1 −7.92±3.21i −4.733±3.51i −7.92±3.21i s∗2 −9.61±8.82i −4.65±2.755i −9.61±8.82i

In Table3, the three utmost right eigenvalues of the discrete spectrum are compared to their continuous analogues for N

=

60 for all sets of boundary conditions. Clearly, the high-order operator approximate the continuous spectrum much better than the low-order one. All cases for the imposition of Set 1, Set 2, Set 3 and Set 4 reside in the left half of the complex plane. However, as indicated by the stability analysis, using Set 5 or Set 6 lead to eigenvalues on the right side of the imaginary axis. Similarly to the continuous setting, the real part of the utmost right lying eigenvalue determines the discrete decay in time and it is denoted by η∗D. The low-order approximation is poor in all cases except for the periodic case (Set 4). This indicates that a von Neumann analysis for a low-order method can yield, or indicate, an overly positive result. For realistic problems with boundary conditions, steeper gradients will be present and high order is required.

Next, see Table 4,

|

vp

|

is computed (by using (10) and (16) for the eigenvalues in Table3) for Set 1, Set 2 and Set 3. Again, the high-order operator is superior compared to the low-order one. The purely real eigenvalues in Table3result in vp

=

0 in Table4.

Every eigenvalue in s∗ have a corresponding mode W .

˜

Similarly, every discrete eigenvalue

in s∗_D has a corresponding mode W.

˜

In Fig.

13, we compare the continuous mode corresponding to s∗

= −

3

.

92

+

2

.

43i for t

=

0 to the discrete mode corresponding to the closest eigenvalue s∗_D for SBP21 and SBP84. In Fig. 14, the time evolution of u of the same modes are illustrated. For the SBP21 operator, an oscillatory behavior is polluting the discrete mode, which is not the case for the high-order SBP84 operator.

(15)

Table 4

Thephasevelocities(computedfrom(10) and(16))ofthe modescorrespondingtotheeigenvaluesinTable3forSet1, Set2andSet3.AllpurelyrealeigenvaluesinTable3result invp=0. vp= |dx/dt| = |s∗I/κ3,I| = |s∗I/κ4,I| s∗ u¯=1 SBP21 SBP84 Set 1 s∗0 0.920 0.924 0.920 s∗1 0.860 0.913 0.860 s∗₂ 0.817 0.919 0.834 Set 2 s∗₀ 0.858 0.000 0.858 s∗1 0.848 0.911 0.849 s∗2 0.842 0.916 0.842 Set 3 s∗0 0.000 0.000 0.000 s∗1 0.794 0.905 0.793 s∗2 0.780 0.000 0.780

5. Thefullydiscreteformulation

To verify the continuous predictions of the convergence to steady state and to exemplify the different behaviors on an inﬁnite and a ﬁnite domain, we will evolve the solution in time. Since there is no time-derivative on the pressure in (1), the implicit SBP-in-time [35,24,36] procedure will be used.

5.1. Fullydiscretenonlinearstability

Following [33], consider a two-dimensional spatial grid with

(

N

+

1

)(

M

+

1

)

points and L

+

1 time levels. The fully discrete approximation of a variable u

=

u

(

t

,

x

,

y

)

is a vector of length

(

N

+

1

)(

M

+

1

)(

L

+

1

)

arranged as u

= (

u0

. . . ,

uL

)

, where each component uk, k

=

0

,

. . . ,

L holds

the solution at the

kth

time level and is of size

(

N

+

1

)(

M

+

1

)

. For the INS equations, we let U

= (

U0

,

. . . ,

UL

)

, where each Uk

= (

uk

,

vk

,

pk

)

contains the solution for all variables over the whole domain for the kth

time level. The approximation of the time derivative is

DtU

= (

Pt−1Qt

⊗

I3

⊗

Ix

⊗

Iy

)

U

≈

Ut

,

where Qt satisﬁes the SBP-property (17) and Ut is a vector containing the time-derivative of the solution evaluated on the fully discretized grid. The matrix Pt is diagonal and positive deﬁnite.

By using the deﬁnition introduced above, the fully discrete SBP-SAT approximation of (1) with homogeneous boundary conditions becomes

(

It

⊗ ˜

I

)

DtU

+

Ds

(

U

)

=

σ

t

(

Pt−1E0

⊗ ˜

I

)(

Us

−

f

) ,

(28) where Ds

(

U

)

= [

Dh

(

U0

),

. . . ,

Dh

(

UL

)

]

is a vector containing the discrete spatial operator (Dh

(

Uk

)

is deﬁned in (19)) for every time step. We assume that S ATx and S ATy in Dh

(

Uk

)

satisfy Proposition2, which is true for Set 1, Set 2, Set 3 and Set 4, but not for Set 5 and Set 6 (see Section 4.2). The right-hand side of (28) is a temporal penalty term imposing the initial data f

= [

f

,

0

. . . ,

0

]

at the ﬁrst time step. The penalty parameter σt will be determined from the following stability analysis [35].

Let Pf

=

Pt

⊗

P,

then

U

2Pf

=

U _P

fU deﬁnes

a fully discrete norm with respect to time and space. Multiplying (

28) by 2UPf from the left and using the SBP-property of Qt yield

UN

2_˜_IP

−

U0

2_˜_IP

+

D I

=

2

σ

tU0P

(

U0

−

f

) .

The term D I is

dissipative and obtained from the spatial part in (

20) at every time step

D I

=

L

i=0 pii

2

∇

Ui

2_˜_IP

−

BTi

≥

0

,

(29)

where pii

>

0 are the diagonal entires in Pt. The inequality in (29) holds for Set 1, Set 2, Set 3 and Set 4 since the spatial penalty terms satisfy Proposition2and yield BTk

≤

0, k

=

0

,

. . . ,

L.

It does not hold for Set 5 and Set 6. Adding and

subtracting

f

2

˜

IPand choosing σt

= −

1 yield

UL

2_˜_IP

≤

f

_IP_˜2

−

U0

−

f

2_˜_IP

.

(30)

Hence,

U

2_˜

IP is bounded at the ﬁnal time level for Set 1, Set 2, Set 3 and Set 4. The estimate in (30) is the discrete equivalence to the continuous estimate in the proof of Proposition1. The additional term

−

U0

−

f

2_˜_IPadds extra dissipation to the scheme.

(16)

Fig. 13. Comparisonbetweenthecontinuousanddiscretepartofthemodecorrespondingtos∗= −3.92+2.43i forN=60 andt=0 withtheboundary conditionsinSet1.AnoscillatorybehaviorpollutesthediscretemodefortheSBP21operator,whichisnotthecasefortheSBP84operator,thatproperly resolvethemodes.

(17)

Fig. 14. Comparisonofthetimeevolutionbetweenthecontinuous(solidanddashedlines)anddiscretemodecorrespondingtos∗= −3.92+2.43i (forSet 1)forN=60.Alsohere,oscillationspollutes thesolutionfromSBP21operator.Thearrowsintheﬁguresshowthedirectionofmovement.

5.2. Theorderofaccuracy

The method of manufactured solution is used to verify the implementation of (28). The time step is chosen such that the error from the time integration can be neglected and the simulations are terminated at t

=

1. A trust-region method [25], which is an optimization algorithm that can be used to solve systems of nonlinear equations, is used to evolve the solution in time via the Matlab function fsolve.

For the SBP-operators SBP21 and SBP84, the expected order of accuracy are 2 and 5 [48]. The manufactured solution is u

=

1

+

sin

(

π

x

−

0

.

01t

)

sin

(

2

π

y

−

0

.

01t

)

, v

=

p

=

sin

(

π

x

−

0

.

01t

)

sin

(

2

π

y

−

0

.

01t

)

and

=

0

.

01. The L2-error at the ﬁnal time UL is

e

2

=

ePe,

where

e is

the point-wise error and

P is

the quadrature corresponding to the SBP-operators.

Furthermore, the convergence rate is r

=

log

(

e

i

/

e

i−1

/

log

(

hi

/

hi−1

)

, where h is

the grid spacing and

i refers

to a speciﬁc

number of grid points. The error and convergence rate for the boundary conditions in Set 1 are presented in Table5and agree well with the theoretical rates.

5.3. Thedecayratetosteady-state

A small disturbance from the constant background ﬁeld

(

u

¯

,

v

¯

,

p

¯

)

should decay as eη∗Dt, which in turn should approach eη∗Ct (the predicted decay rate from the continuous analysis) for an increasingly ﬁne mesh. To illustrate that, we choose the background

(

u

¯

,

v

¯

,

p

¯

)

= (

1

,

0

,

0

)

, the initial state u

= ¯

u

+

exp

(

−(

x

−

0

.

5

)

2

_/

₀

_.

₀₅

₎

_sin

₍

_y

₎

_{and diffusion parameter}

₌

₀

_.

_{01. In} all computations below, 21 grid-points is used in the periodic y-directions.

In the left panel of Fig.

15, the decay rates for different resolutions in the x-direction

(N

=

30

,

60) are illustrated for SBP21 and SBP84 using the boundary conditions in Set 1. The theoretical decay rate is included for comparison, where the parameter η∗_C (the continuous decay in time) is the real part of s∗₀ in Table1.

(18)

Fig. 15. Left:Comparisonbetweenthecomputedandtheoreticaldecayrate(dashedline)forSet1.Right:Decayratefromsteady-statecomputationsforall setsofboundaryconditions.ThetheoreticaldecayratesforSet1,Set2andSet3(dashedlines)areincludedforcomparison.Thesolutiongrowsintime forSet5andSet6.Inthecomputations,theSBP84operatorhasbeenusedonagridofsizeN×M=100×20.

As the number of grid points increases, the numerical decay rate converges to the one predicted from the continuous analysis. Again, the high-order operator is superior to the low-order operator. In the right panel of Fig. 15, we use the SBP84 operator for all sets of boundary conditions and N

=

100. The theoretical decay rates predicted for Set 1, Set 2 and Set 3 are included and the results agree reasonably well with theory (the spectrum of Set 2 is located further to the left in Fig.3compared to the other sets and yields a more stiff problem, which is harder to resolve). As indicated by the results presented in Fig.4, the periodic case decay slower than Set 1, Set 2 and Set 3. Note that Set 5 and Set 6 lead to growth in time.

5.4. Thedifferencebetweenaﬁniteandaninﬁnitedomain

Since we are considering a ﬁnite domain, the solution can decrease in two different ways: i

)

by a diffusion mechanism in time related to esRt _and_ii

₎

_by_{a propagation mechanism related to}_ei(sIt+κIx)_._{This is in contrast to}_an_inﬁnite_domain case, where only factor i

)

is relevant. In fact, an increased diffusion coeﬃcient may lead to a slower decay on ﬁnite domains since the solution becomes more smeared and elongated, as discussed in Section3.5.

To illustrate this combined effect, we change the initial state to u

= ¯

u

+

0

.

01 exp

(

−(

x

−

0

.

5

)

2

/

0

.

045

)

sin

(

2 y

)

and compare the time-evolution of the solution for different for Set 2 on the domain

[

0

,

2

]

× [−

π

,

π

]

with N

×

M

=

100

×

20. Note that the initial state is a thin pulse that decays fast away from x

=

0

.

5.

Consider the left panel in Fig.16. For

=

0

.

025, the solution decays slowly in the beginning since most of the disturbance is still inside the domain. When the disturbance leaves the domain and interacts with the boundary, the solution decays rapidly. For

=

0

.

1, the decay rate is more evenly distributed in time. In the beginning, it decays more rapidly compared to the case with

=

0

.

025, while the situation is reversed after t

≈

1

.

8 as the pulses hit the boundary. The right panel of Fig. 16illustrates the situation at t

=

2

.

25, when most of the disturbance has left the domain. The increased diffusion coefficient results in a more elongated solution and slower decay. To compare these finite domain results to results on an infinite domain, we have performed the same experiment but used periodic boundary conditions in the x-direction.

The

only decay is now due to the diffusion mechanism since the whole pulse remains inside the computational domain. In this case, larger result in a consistently faster decay. Both types of experiments lead to identical decays during the ﬁrst part of the simulation (before the disturbance reaches the boundary).

(19)

Fig. 16. Leftfigure:decayforthesteady-statecomputationsonaboundedandaninfinitedomain.Rightfigure:Screenshotsofthesolutionfordifferent

att=2.25.Theleftpanelcorrespondsto=0.025 andtherightpanelto=0.1. 6. Summaryandconclusions

We have investigated the influence of open boundary condition for the INS equations on the spectrum of the IBVP and compared the result to the fully periodic IVP. By using a von Neumann and normal mode analysis, we calculated the spectrum of the linearized equations. It was shown that both the decay in time and the propagation characteristics influence the convergence rate of the solution. We compared the convergence rates on infinite and finite domains and found that a larger diffusion coefficient may lead to slower convergence rate on a finite domain since the solution becomes more elongated in space. Furthermore, the relation between the phase velocity and the spectrum of the IBVP was illustrated.

The equations were discretized by using the SBP-SAT technique and the discrete spectrum was derived. We verified the convergence of the discrete spectrum to its continuous analogue and computed the resulting phase velocity from the discrete spectrum. Furthermore, the theoretical predictions of the convergence rate to steady state were verified by nonlinear computations and we highlighted the importance of correctly specified boundary conditions. The sets of boundary conditions that produced spectra in the right half of the complex plane did not converge to steady state. Two different operators were used in the computations. The approximations of the high-order operator predicted the time-dependent solution without oscillatory behavior, in contrast to the low-order one. Finally, as previously suggested by the continuous analysis, the choice of boundary condition was shown to heavily impact the convergence rate to steady state.

CRediTauthorshipcontributionstatement

FredrikLaurén:

Theoretical development. Jan

Nordström:

Theoretical development.

Declarationofcompetinginterest

The authors declare that they have no known competing ﬁnancial interests or personal relationships that could have appeared to inﬂuence the work reported in this paper.

Acknowledgements

The authors acknowledge the support from The Swedish e-Science Research Centre (SeRC).