' $
Department of Mathematics
A Well-posed and Stable Stochastic Galerkin
Formulation of the Incompressible
Navier-Stokes Equations with Random Data
Per Pettersson, Jan Nordstr¨om and Alireza Doostan
LiTH-MAT-R--2015/06--SE
Department of Mathematics
Link¨oping University
A Well-posed and Stable Stochastic Galerkin
Formulation of the Incompressible Navier-Stokes
Equations with Random Data
Per Petterssona,∗, Jan Nordstr¨omb, Alireza Doostanc
aUni Research, N-5007 Bergen, Norway.
bDepartment of Mathematics, Computational Mathematics, Link¨oping University,
SE-58183 Link¨oping, Sweden.
cAerospace Engineering Sciences, University of Colorado Boulder, CO 80309, USA.
Abstract
We present a well-posed stochastic Galerkin formulation of the incompressible Navier-Stokes equations with uncertainty in model parameters or the initial and boundary conditions. The stochastic Galerkin method involves represen-tation of the solution through generalized polynomial chaos expansion and projection of the governing equations onto stochastic basis functions, result-ing in an extended system of equations. A relatively low-order generalized polynomial chaos expansion is sufficient to capture the stochastic solution.
We derive boundary conditions for an energy estimate that leads to zero divergence of the velocity field. In other words, the incompressibility con-dition is not imposed directly in the problem formulation but is instead a consequence of the combination of the partial differential equations and the boundary conditions.
Based on the analysis of the continuous equations, we present a semi-discretized system where the spatial derivatives are approximated using finite difference operators with a summation-by-parts property. With a suitable choice of dissipative boundary conditions imposed weakly through penalty terms, the semi-discrete scheme is shown to be stable. Numerical experiments corroborate the theoretical results and we obtain high-order accurate results
∗Corresponding author
Email addresses: per.pettersson@uib.no (Per Pettersson),
jan.nordstrom@liu.se (Jan Nordstr¨om), alireza.doostan@colorado.edu (Alireza Doostan)
for the solution variables and the velocity divergence converges to zero as the mesh is refined.
Keywords: Uncertainty quantification, Incompressible Navier-Stokes
equations, Summation-by-parts operators, Stochastic Galerkin method, Boundary conditions
1. Introduction
Complex fluid flow problems are frequently dependent on model param-eters that rely on approximations. Accurate and reliable numerical simula-tions therefore require that we properly account for uncertainties in param-eters, and in initial and boundary conditions. Forward uncertainty quan-tification deals with systematic and quantitative means of accounting for input uncertainty and methods to propagate it through the numerical model of interest. This process in turn requires numerical methods that produce accurate solutions under uncertainty [1].
A framework for efficient representation and propagation of uncertainty is offered by polynomial chaos methods, e.g., the stochastic Galerkin method where the governing equations are projected onto stochastic basis functions [2, 3]. The aim of this paper is to present and investigate a new stochastic Galerkin formulation of the incompressible Navier-Stokes equations and a corresponding finite difference formulation. We are mainly interested in the case of uncertainty in the viscosity, but provide a framework that may include uncertainty in other model parameters as well as uncertainty in the initial and boundary conditions.
In numerical computations, the governing equations are approximated within the order of accuracy of the numerical method. Neither of them is satisfied exactly on a finite computational grid, but we will derive a method that leads to high-order convergence in all solution variables. By enforcing an advection-diffusion equation for the velocity divergence in combination with suitable boundary conditions, the divergence of the velocity is converging to zero with the order of the scheme.
The stochastic Galerkin formulation introduces an error due to the trun-cated gPC representation of the problem, and with a low-order stochastic approximation, this error may dominate the discretization error in both the solution variables and the divergence. However, even in the absence
gence for a provably stable solution of high order of accuracy where the error in the solution variables and the divergence can be made arbitrarily small through mesh refinement.
The incompressible Navier-Stokes equations with uncertain viscosity was previously investigated in [4], where a stochastic Galerkin formulation was presented and implemented with a fractional step method resulting in an advection-diffusion problem decoupled from the solution of the pressure and divergence system. The method was generalized to the weakly compress-ible Navier-Stokes equations and random process input uncertainty in [5]. A similar stochastic Galerkin method for the incompressible Navier-Stokes equations was devised in [6], and evaluated for different uncertain bound-ary conditions. Efficient representation of uncertainty was ensured by the choice of stochastic basis, i.e., orthogonal polynomials with the inner prod-uct weight function matching the probability density function of the bound-ary conditions. The stochastic Galerkin formulation of the incompressible Navier-Stokes equations was also employed for fluid-structure interaction to study vortex-induced vibrations in [7], and for flow under stochastic non-Boussinesq conditions in [8].
The purpose of this work is not a direct comparison with the previous incompressible Navier-Stokes stochastic Galerkin solvers listed above; that is the goal of a future investigation. Instead, we focus on the development of a new stochastic Galerkin solver and perform detailed analysis to prove relevant properties of the continous formulation and the proposed numerical method. The stochastic Galerkin projection results in an extended system of equations that is dependent on both the sources of uncertainty and the stochastic rep-resentation. In order to design general, robust and stable numerical schemes that can handle different forms of uncertainty and different stochastic rep-resentations, we use finite-difference operators with a summation-by-parts property [9, 10]. This allows for the design of methods that are provably time stable [11] in combination with boundary conditions that are imposed weakly, i.e., through discrete penalty terms that vanish with spatial mesh refinement [12].
The analysis of well-posedness and stability for the stochastic Galerkin formulation presented in this paper extends the analysis presented for the deterministic incompressible Navier-Stokes problem in [13]. In that analysis, the condition of exact enforcement of the incompressibility condition was re-laxed in order to gain overall stability and accuracy of the numerical method. We provide a stochastic Galerkin formulation that leads to an energy estimate
for the velocity divergence. In combination with certain data requirements, this estimate shows that the zero divergence condition indeed is satisfied. Fi-nally, by mimicking that procedure for the semi-discretized problem, we can prove stability and high-order convergence. The result is a stable high-order method that handle stochastic, non-periodic and time-dependent boundary conditions for the incompressible Navier-Stokes equations.
To appreciate both the similarities and the differences between the deter-ministic and the stochastic Galerkin formulation, we reproduce some of the results from [13] along with new results. The deterministic analysis serves as a roadmap for the stochastic Galerkin analysis, and the analysis of the contin-uous problem serves as a guide for the analysis of the semi-discrete problem. By proceeding step-wise, from continuous to semi-discrete formulation and from deterministic to stochastic Galerkin formulation, we derive a numerical method for the solution of the stochastic problem that may otherwise be difficult to derive.
In Section 2, we present the generalized polynomial chaos framework, and the stochastic Galerkin incompressible Navier-Stokes equations are presented in Section 3. The continuous formulation is analyzed in Section 4. The deterministic analysis is reproduced from [13], followed by analysis of well-posedness for the stochastic Galerkin formulation. The continuous analysis provides boundary conditions and a recipe for the analysis of stability of the spatially discretized problem in Section 5. Numerical results are presented in Section 6 for a Taylor vortex and a Poiseuille flow problem with stochastic viscosity. Finally, conclusions are discussed in Section 7.
2. Representation of Uncertainty
Following the seminal work in [2, 14], we represent stochastic functions by using the polynomial chaos framework. To this end, consider a probability space (Ω, F , P) with the set of elementary events Ω and probability measure
P defined on the σ-algebra F . Let ξ = (ξ1, . . . , ξd)T be a d-dimensional
independent random vector defined on this space. The inner product between two stochastic functions f (ξ) and g(ξ) is defined by
hf gi = Z
Ω
f (ξ)g(ξ)dP(ξ).
span-probability space. Possible choices of basis functions include the orthogonal polynomials from the Askey scheme [15] and, for functions with non-smooth behavior, stochastic multiwavelets [16]. Multi-dimensional basis functions are generated by tensor products of one-dimensional basis functions. The basis functions are assumed to be orthonormal with respect to the inner product h·i with the measure P, i.e., they satisfy
hψiψji = δij, i, j ∈ N,
where δij is the Kronecker delta. Any second order random field u(x, t, ξ)
(typically the solution of a partial differential equation) can be expressed as a gPC expansion in the form
u(x, t, ξ) =
∞
X
i=1
ui(x, t)ψi(ξ), (1)
where the coefficients ui(x, t) are defined by the projections
ui(x, t) = hu(x, t, ξ)ψi(ξ)i, i ∈ N. (2)
The solution statistics are functions of the gPC coefficients and can easily be obtained by post-processing of the solution. For instance, the expectation and variance are given by, respectively,
E(u) ≡ hui = u1, V ar(u) =
∞
X
i=2
u2i. (3)
In practice, the infinite expansion in (1) is truncated to a finite number of
terms such that, for instance, the total polynomial order of each ψi(ξ) is at
most p. The total number of gPC terms is then prescribed a priori, unless an adaptive method is used, in that case the number of basis functions change over time [17].
3. Problem Formulation
For ease of presentation, we consider the two-dimensional incompressible Navier-Stokes equations, while the proposed developments are applicable to
the three-dimensional setting. Let x = (x, y) and consider the spatial domain D = [0, 2] × [0, 2]. We have
ut+ uux+ vuy+ px− [(µux)x+ (µuy)y] = 0, (4)
vt+ uvx+ vvy + py − [(µvx)x+ (µvy)y] = 0, (5)
ux+ vy = 0. (6)
In (4)-(6), u and v are the velocity components in the x- and y-directions, respectively, p is the pressure, and µ is the viscosity. The subscripts x, y
and t refer to the spatial and temporal derivatives. Setting φ = ux+ vy and
taking the divergence of (4)-(5) without imposing (6), gives
φt+ uφx+ vφy − (µφx)x− (µφy)y = F, (7)
where
F = (µxux)x+ (µxuy)y+ (µyvx)x+ (µyvy)y− u2x− 2uyvx− vy2− pxx− pyy. (8)
In the case of a constant µ, which was considered in [13], we have F =
−u2
x− 2uyvx− vy2− pxx− pyy.
In the numerical calculations we approximate the relations (4)-(7) within the order of accuracy of the numerical method. Neither of them is satisfied exactly on a finite computational grid, but we will derive a method that leads to high-order convergence in all solution variables u, v, p and the divergence of the velocity converging to zero. By enforcing (7) in combination with suitable boundary conditions and a Poisson equation for the pressure, we also satisfy (6) without explicitly imposing it.
3.1. Stochastic Galerkin Formulation
We next consider the case where the viscosity µ in (4)-(5) as well as the boundary and initial conditions are uncertain and characterized by an independent random vector ξ. This implies that the velocity components u, v and the pressure p are also random functions, represented by the truncated
gPC series u(x, t, ξ) = M X i=1 ui(x, t)ψi(ξ), (9) v(x, t, ξ) = M X i=1 vi(x, t)ψi(ξ), (10) p(x, t, ξ) = M X i=1 pi(x, t)ψi(ξ). (11)
The number of basis functions M should be chosen sufficiently large to ac-curately represent the stochastic variation in the output, which could be unknown a priori. Spectral stochastic projection is performed by
multiply-ing (4) and (5) with ψk for k = 1, . . . , M and integrating w.r.t. the measure
P. Following the stochastic Galerkin formulation, we require the residual of
the equation to be orthogonal to each and every ψk, which leads to the
fol-lowing coupled systems of partial differential equations (PDEs) for the gPC expansion coefficients. (uk)t+ M X i=1 M X j=1 [ui(uj)x+ vi(uj)y] hψiψjψki + (pk)x − M0 X i=1 M X j=1 [(µi(uj)x)x+ (µi(uj)y)y] hψiψjψki = 0, (12) (vk)t+ M X i=1 M X j=1 [ui(vj)x+ vi(vj)y] hψiψjψki + (pk)y − M0 X i=1 M X j=1 [(µi(vj)x)x+ (µi(vj)y)y] hψiψjψki = 0. (13)
The viscosity µ is represented by a gPC expansion of order M0, where M0 ≥
M is chosen sufficiently large for maximum accuracy of the projection of (4)
and (5) on the finite space spanned by {ψk}Mk=1.
Let uM = (u
1, . . . , uM)T and vM = (v1, . . . , vM)T denote the vectors of
gPC coefficients of u and v, respectively. Let the matrix A(uM
by [A(uM )]kj = M X i=1 uihψiψjψki , j, k = 1, . . . , M.
By the definition of A(.), for any vectors aM, bM
∈ RM, it holds that [A(aM )bM ]k= M X j=1 M X i=1 hψiψjψki aibj = M X j=1 M X i=1 hψiψjψki biaj = [A(bM)aM]k,
so we can interchange the vector valued argument of A(.) with the vector the matrix is multiplied with, i.e. the order of the two factors of the product does not matter. Thus, we arrive at the expression
A(aM
)bM
= A(bM
)aM
, (14)
that is frequently used in our derivations. By using (14), the stochastic Galerkin formulation of the momentum equations (12)-(13) can be written
uM t + A(u M )uM x + A(v M )uM y + p M x − [A(µ)u M x]x+ [A(µ)uMy ]y = 0 (15) vM t + A(u M )vM x + A(v M )vM y + p M y − [A(µ)v M x]x+ [A(µ)vMy ]y = 0, (16) where µ = (µ1, . . . , µM0)T and pM = (p1, . . . , pM)T.
Remark 1. The order of the gPC expansion of µ is chosen sufficiently large to guarantee that A(µ) is positive definite. This is explained and analyzed in [18] for a single stochastic dimension.
Performing stochastic Galerkin projection of the incompressibility condi-tion (6) onto the stochastic basis yields
uM
x + v
M
y = 0. (17)
Similarly, stochastic Galerkin projection of (7) gives
φM t + A(u M )φM x + A(v M )φM y − (A(µ)φ M x)x− (A(µ)φMy )y = FM, (18) where FM = (A(µx)uM x)x+ (A(µx)u M y )y+ (A(µy)v M x)x+ (A(µy)v M y )y− A(uMx)u M x
Remark 2. The stochastic Galerkin version of the incompressibility condi-tion is (17), i.e., a zero-divergence condicondi-tion for each gPC coefficient of the velocity components.
Remark 3. Note that by setting M = 1 in (15)-(19), we recover the deter-ministic case (4)-(8).
4. Analysis of Well-posedness
As was indicated above, the incompressibility condition (17) is not explic-itly enforced. As we will see, (17) is enforced indirectly by solving (18) with
FM = 0 and a specific set of boundary conditions. To ensure
incompress-ibility, we will derive an energy estimate that shows that incompressibility is indeed satisfied. The energy estimates in the following sections rely on time integration of expressions for the temporal derivative of the norm of the quantity of interest. Without subscripts specifying a weight function, k.k denotes the standard Euclidean norm.
4.1. An Energy Estimate for Divergence Preservation
We start by first left-multiplying the advection-diffusion equation (18) for
the divergence gPC coefficients by 2(φM
)T
. Next, we use that (φM
)T A(uM ) = (uM)TA(φM ) and (φM )TA(vM) = (vM)TA(φM ) by (14). Finally, we inte-grate over the spatial domain to obtain the following expression for the time derivative of the norm of the velocity divergence.
d dtkφ M k2+ Z D (uM )T [A(φM )φM ]x+ (vM)T[A(φM)φM]ydxdy − 2 Z D (φM )T [A(µ)φM x]x+ (φM)T[A(µ)φMy ]ydxdy = 2 Z D (φM )T FM dxdy. Integration by parts and rearrangement of terms yield
d dtkφ Mk2 + 2 Z D (φM x) T A(µ)φM x + (φ M y ) T A(µ)φM y dxdy + I ∂D (uM n) T A(φM )φM − 2(φM )T A(µ)∂φ M ∂n ds = Z D (φM )T A(φM )φM + 2(φM )T FM dxdy, (20)
where n = (n1, n2)T is the outward pointing normal, and we have adopted the notation uM n ≡ n1uM + n2vM, (21) ∂φM ∂n ≡ n1φ M x + n2φ M y . (22)
In order to arrive at an energy estimate for the norm of the divergence coefficients, the boundary integral terms in (20) will be rewritten in a form that allows us to identify growth terms, which can then be controlled by
subscription of boundary conditions. Assuming that A(uM
n) is invertible
everywhere along the boundary ∂D implies that the boundary terms in (20) can be rewritten as (φM )T A(uM n)φ M − 2(φM )T A(µ)∂φ M ∂n = A(uM n)φ M − A(µ)∂φ M ∂n T [A(uM n)] −1 A(uM n)φ M − A(µ)∂φ M ∂n − A(µ)∂φ M ∂n T [A(uM n)] −1 A(µ)∂φ M ∂n . (23)
Let the eigenvalue decomposition [A(uM
n)] −1
= XΛ−1XT be given, where X
and Λ−1 may both vary in space. Then the boundary terms are given by
XT A(uM n)φ M − A(µ)∂φ M ∂n T Λ−1 XT A(uM n)φ M − A(µ)∂φ M ∂n − XT A(µ)∂φ M ∂n T Λ−1 XT A(µ)∂φ M ∂n . (24) To get a bound on (23), we assign boundary conditions to prescribe the growth terms, XT j A(uM n)φ M − A(µ)∂φM ∂n if [Λ]jj < 0 XT jA(µ) ∂φM ∂n if [Λ]jj ≥ 0 ) = gj, (25)
data. By using (25) in (24), we get the estimate − I ∂D (uM n) T A(φM )φM − 2(φM )T A(µ)∂φ M ∂n ds ≤ I ∂D gT|Λ−1 |g − [XT A(uM n)φ M − g]T |Λ−1 | [XT A(uM n)φ M − g] ds, (26)
for the boundary integral in (20).
Remark 4. If A(uM
n) is not invertible but has at least one nonzero
eigen-value, we may write its eigenvalue decomposition as
A(uM n) = X X0 Λ 0 0 0 XT XT 0 ,
where X contains the basis vectors corresponding to the nonzero eigenvalues
(Λ) and X0 contains the vectors of the nullspace. Then, since this implies
A(uM
n) = XΛX
T
, we may use the latter reduced eigenvalue decomposition and proceed as indicated above to prescribe boundary conditions and obtain an estimate for the boundary integral terms.
For the volume terms on the right hand side of (20), we use the inequality
2(φM )T FM ≤ η(φM )T φM + (FM )T FM
/η, for any real η > 0. (27)
Consider sufficiently short times t ∈ [0, T ] and let δ = max
t,x eig(A(u M x + v M y )). It holds that Z D (φM )T A(uM x +v M y )φ M +2(φM )T FM dxdy ≤ (δ +η) kφMk2 +1 ηkF Mk2 . (28) Inserting the estimates (26) and (28) into (20), multiplying by the factor e−(δ+η)t and integrating over time from 0 to T , we obtain the energy estimate
kφM k2+ Z T 0 h 2 kφM xk 2 A(µ)+ 2 φM y 2 A(µ) i e−(δ+η)(t−T )dt + Z T 0 kXT A(uM n)φ M − gk2|Λ−1|e−(δ+η)(t−T )dt ≤ kφM (t = 0)k2e(δ+η)T + Z T 0 kgk2|Λ−1|+ 1 ηkF Mk2 e−(δ+η)(t−T )dt, (29)
where k.k|Λ−1| is the L2 vector norm weighted by |Λ −1
|. We impose zero
initial and boundary data, φM
(t = 0) = g = 0, and require that FM
= 0.
Then the argument is repeated for a sufficiently short time interval [T, T0],
and we obtain kφM
(t = T0)k2 ≤ kφM
(t = 0)k2. Iterating this argument for
sufficiently short time intervals, we conclude that φM
= 0 for all times. Thus, the zero divergence condition is satisfied for all times.
The boundary conditions (25) are not the only possible dissipative choice for the divergence of the velocity. A similar but slightly less dissipative energy estimate is obtained for the boundary condition
φM
= 0, x, y ∈ ∂D. (30)
4.2. An Energy Estimate for the Momentum Equations
The derivation of well-posedness for the momentum equations and the resulting energy estimate is a generalization of the deterministic case. For completeness and appreciation of the similarities and differences imposed by the stochastic Galerkin formulation, in Section 4.2.1 we start with the deter-ministic derivation published in [13]. With a roadmap given by that anal-ysis, we then proceed with the analysis of well-posedness for the stochastic Galerkin formulation of momentum equations in Section 4.2.2.
4.2.1. The Deterministic Case
Multiplying the momentum equations (4)-(5) from the left with the ve-locities and integrating over the spatial domain D, we get
d dt kuk 2 + kvk2 + 2 Z D u2ux+ uvuy + uvvx+ v2vy | {z } 1 2u(u2+v2)x+ 1 2v(u2+v2)y dxdy +2 Z D (upx+vpy)dxdy−2 Z D
(u(µux)x+u(µuy)y+v(µvx)x+v(µvy)y)dxdy = 0.
Integration by parts of (31) yields d dt kuk 2 + kvk2 + I ∂D (u2+ v2)(un1+ vn2)ds − Z D (ux+ vy)(u2+ v2)dxdy | {z } Advective terms + 2 I ∂D (un1+ vn2)pds − 2 Z D (ux+ vy)pdxdy | {z } Pressure terms −2 I ∂D µ(uux+ uuy+ vvx+ vvy)ds + 2 k∇uk2µ+ k∇vk2µ | {z } Viscous terms = 0, (32)
where n = (n1, n2)T is the outward pointing normal to ∂D. It was shown
in Section 4.1 that φM
= uM + vM = 0 for all M including M = 1 which
is equivalent to the deterministic case. Thus, the volume integrals related to the advective and pressure terms in (32) vanish. The boundary integral terms can be rewritten as
I
∂D
(u2+ v2)(un1+ vn2) + 2(un1+ vn2)p − 2µ(uux+ uuy+ vvx+ vvy)ds
= I ∂D wT Bwds, where w = un us p µ∂un ∂n µ∂us ∂n , B = un 0 1 −1 0 0 un 0 0 −1 1 0 0 0 0 −1 0 0 0 0 0 −1 0 0 0 . (33)
The eigenvalues of B are
λ1,2 = un 2 ± r un 2 2 + 1, λ3 = 0, λ4,5 = un 2 ± r un 2 2 + 1. (34)
Boundary conditions are imposed on the characteristics corresponding to the
negative eigenvalues λ2 and λ5, i.e. (XTw)j = gj for j = 2, 5 with X being
the matrix of eigenvectors of B. This ensures that the correct number of boundary conditions is imposed.
4.2.2. The Stochastic Galerkin Case
Multiply (15)-(16) by [(uM)T, (vM)T] and integrate over the spatial
do-main. The resulting expression can be written d
dt
kuM
k2+ kvM
k2+advective terms+pressure terms+diffusive terms = 0.
(35) For clarity, we define and treat the terms groupwise below.
For the advective terms, we use that A(.) is linear in its argument to rewrite the spatial derivatives. Then we perform integration by parts,
2 Z D (uM )T A(uM )uM x + (v M )T A(uM )vM x + (u M )T A(vM )uM y dxdy + 2 Z D (vM )T A(vM )vM y dxdy = Z D (uM )T [A(uM )uM + A(vM )vM ]x+ (vM )T [A(uM )uM + A(vM )vM ]ydxdy = I ∂D (uM n1+ vMn2)T[A(uM)uM + A(vM)vM] ds − Z D (φM )T [A(uM )uM + A(vM )vM ] dxdy.
In the first equality, we used that 2A(uM)uM
x = A(u M)uM x + A(u M x)u M = [A(uM )uM
]x, and similar expressions for vM and the y-derivatives.
For the pressure terms, we perform integration by parts, 2 Z D (uM )T pM x +(v M )T pM y dxdy = 2 I ∂D (uM n1+vMn2)TpMds−2 Z D (φM )T pM dxdy. For the diffusive terms, we perform integration by parts and rewrite the
resulting volume integral, assuming that A(µ) defines a norm, i.e. that
positive almost surely (in the probabilistic sense). We obtain 2 Z D (uM )T [A(µ)uM
x]x+ (uM)T[A(µ)uMy ]y + (vM)T[A(µ)vMx]x dxdy
+ 2 Z D (vM )T [A(µ)vM y ]ydxdy = 2 I ∂D n1(uM)TA(µ)uMx + n2(uM)TA(µ)uMy ds + 2 I ∂D n1(vM)TA(µ)vxM + n2(vM)TA(µ)vMy ds −2 Z D (uM x) T A(µ)uM x + (u M y ) T A(µ)uM y + (v M x) T A(µ)vM x + (v M y ) T A(µ)vM y dxdy | {z } kuM x k 2 A(µ)+kvMx k 2 A(µ)+kuMy k 2 A(µ)+kv M y k 2 A(µ) .
Assembling the advective, pressure and diffusive terms, the expression (35) becomes d dt kuM k2+ kvM k2+ 2 h kuM xk 2 A(µ)+ kv M x k 2 A(µ)+ uM y 2 A(µ)+ vM y 2 A(µ) i = −2 I ∂D (uM n1+ vMn2)T[A(uM)uM/2 + A(vM)vM/2 + pM] ds + 2 I ∂D
n1(uM)TA(µ)uMx + n2(uM)TA(µ)uMy + n1(vM)TA(µ)vMxds
+2 I ∂D n2(vM)TA(µ)vMy ds+2 Z D (φM )T [A(uM )uM /2 + A(vM )vM /2 + pM ] dxdy. (36)
In order to rewrite (36) to obtain an energy estimate, we use s = (n2, −n1)T
to denote the vector tangential to the spatial boundary ∂D. Let
uM s ≡ n2uM − n1vM, (37) ∂uM n ∂n ≡ n1(u M n)x+ n2(uMn)y, (38) ∂uM s ∂n ≡ n1(u M s )x+ n2(uMs )y. (39)
By combining the property (14) of matrix A(.) with (21) and (37), we get the identity A(uM n)u M n + A(u M s )u M s = [n1A(uM) + n2A(vM)] [n1uM + n2vM] + [n2A(uM) − n1A(vM)] [n2uM − n1vM] = (n21+ n 2 2)A(u M )uM + (n21+ n22)A(vM )vM = A(uM )uM + A(vM )vM . (40)
Also, using (21) and (37)-(39), we have (uM n) T A(µ)∂u M n ∂n + (u M s ) T A(µ)∂u M s ∂n = (n1uM + n2vM)TA(µ)[n21u M x + n1n2vMx + n1n2uMy + n22v M y ] + (n2uM − n1vM)TA(µ)[n1n2uMx − n 2 1v M x + n 2 2u M y − n1n2vMy ] = n1(n21+ n 2 2)(u M )T A(µ)uM x + n2(n21+ n 2 2)(u M )T A(µ)uM y + n1(n21 + n 2 2)(v M )T A(µ)vM x + n2(n21+ n 2 2)(v M )T A(µ)vM y
= n1(uM)TA(µ)uMx+n2(uM)TA(µ)uMy +n1(vM)TA(µ)vMx+n2(vM)TA(µ)vyM.
(41) Finally, by combining (40), (41) and (36), we obtain,
d dt kuMk2 + kvMk2 + 2hkuM xk 2 A(µ)+ kv M x k 2 A(µ)+ uM y 2 A(µ)+ vM y 2 A(µ) i = − I ∂D (wM )T BM wM ds+2 Z D (φM )T [A(uM )uM /2 + A(vM )vM /2 + pM ] dxdy, (42) where wM = uM n uM s pM A(µ)∂uMn ∂n A(µ)∂uMs ∂n , BM = A(uM n) 0 I M −IM 0 0 A(uM n) 0 0 −I M IM 0 0 0 0 −IM 0 0 0 0 0 −IM 0 0 0 . (43) Note that (43) mimics the expression (13) in [13] and (33) in this paper. Just
as in the deterministic case, we need the sign of the eigenvalues of BM to be
able to impose boundary conditions.
4.2.3. Eigenvalues of BM and Imposition of Boundary Conditions
To find the eigenvalues of BM
, consider the equivalent problem of finding
the eigenvalues of a similar matrix. Let Xunbe the matrix of eigenvectors and
A(uM
n) = XunΛunX
T
un. Then, the matrix B
M Λ, defined by BM Λ = Λun 0 I M −IM 0 0 Λun 0 0 −I M IM 0 0 0 0 −IM 0 0 0 0 0 −IM 0 0 0 = (I5 ⊗ XT un)B M (I5⊗ X un), is similar to BM. Thus, 0 = det(λI5M− BM ) = det λIM − Λ un 0 −I M IM 0 0 λIM − Λ un 0 0 I M −IM 0 λIM 0 0 IM 0 0 λIM 0 0 IM 0 0 λIM ,
where I5M is the identity matrix of size 5M × 5M . We will use the fact that
for any matrices ˆA ∈ Rm×m, ˆB ∈ Rm×n, ˆC ∈ Rn×m and ˆD ∈ Rn×n it holds
that det ˆ A Bˆ ˆ C Dˆ = det I Bˆ 0 Dˆ " ˆ A − ˆB ˆD−1C 0ˆ ˆ D−1Cˆ I #! = det( ˆD) det( ˆA − ˆB ˆD−1C), (44)ˆ
assuming that ˆD is invertible [19]. Let
ˆ A = (λIM − Λ un) ⊗ I 2, B = ˆˆ CT = −I M IM 0 0 0 IM , D = λIˆ 3M . Then, by (44) 0 = det(λI5M − BM ) = = det λIM 0 0 0 λIM 0 0 0 λIM det λIM − Λ un− 2λ −1IM 0 0 λIM − Λ un− λ −1IM = λ3M det(diag(λIM − Λ un − 2λ −1 IM , λIM − Λ un− λ −1 IM )) = λM det(diag(λ2IM − λΛ un− 2I M , λ2IM − λΛ un − I M )). (45)
The superscripts on λ denote powers, but superscripts on matrices (boldface) always denote the order of gPC. The solutions λ of (45) are
λ1,j = λun,j 2 + s λ2 un,j 4 + 1, j = 1, . . . , M, (46) λ2,j = λun,j 2 − s λ2 un,j 4 + 1, j = 1, . . . , M, (47) λ3,j = 0, multiplicity M, j = 1, . . . , M, (48) λ4,j = λun,j 2 + s λ2 un,j 4 + 2, j = 1, . . . , M, (49) λ5,j = λun,j 2 − s λ2 un,j 4 + 2, j = 1, . . . , M. (50)
The subscript un refers to the velocity in the normal direction n. As was
mentioned above, the stochastic Galerkin eigenvalues (46)-(50) are direct generalizations of the deterministic eigenvalues (34).
Similarly to the deterministic case, the negative eigenvalues are always
λ2,j and λ5,j, j = 1 . . . , M , independent of both the mean flow direction and
the sign of higher order coefficients. The number of boundary conditions therefore remains exactly 2M .
Remark 5. The fixed number of boundary conditions will be a great advan-tage for the numerical solution of the problem since in this case no special considerations are needed for characteristic curves changing direction. This situation is far from self-evident, and this is in fact not always the case even for scalar linear problems. For example, c.f. [20] for a stochastic Galerkin formulation of an advection equation where the number of boundary condi-tions depends on the input uncertainty.
The boundary conditions are prescribed on the characteristics correspond-ing to the negative eigenvalues. We have
XunΛ2X T unu M s − A(µ) ∂uM s ∂n = g M 2 , ∂uM (51)
where Λ2 = diag(λ2,1, . . . , λ2,M) and Λ5 = diag(λ5,1, . . . , λ5,M).
We showed in Section 4.1 that φM
= 0. Imposing zero divergence and integrating (42) in time using the boundary conditions (51), leads to the energy estimate kuMk2 + kvMk2 ≤ kuM (t = 0)k2+ kvM (t = 0)k2+ Z T 0 kgMk2 ∂Ddt − 2 Z T 0 kuM xk 2 A(µ)+ kv M x k 2 A(µ)+ uM y 2 A(µ)+ vM y 2 A(µ)dt.
The energy estimate is the first building block in an analysis of well-posedness. If the problem was linear, the energy estimate would directly lead to unique-ness and by choosing the correct number of boundary conditions, existence would be guaranteed. An energy estimate, uniqueness and existence imply well-posedness for linear problems. For nonlinear problems, the situation is similar but more technically involved. For more details on uniqueness and existence, we refer the reader to [21] where the relations between linear and nonlinear problems are discussed.
5. Numerical Method
So far, we have derived continuous results for arbitrary spatial domains. Henceforth, for clarity and simplicity we will present the numerical method for the square domain depicted in Figure 1. We introduce a uniform mesh
(xi, yj), for i = 1, . . . , mx and j = 1, . . . , my on the physical domain D =
n s north (n) n s south (s) n s east (e) n s west (w) 0 2 2 x y
Figure 1: Physical domain D with outward normal (n) and tangent (s) vec-tors indicated for each boundary (w, e, s, n).
We use finite-difference operators with a summation-by-parts (SBP) prop-erty [9] for the semi-discrete problem, i.e., the spatial discretization of the momentum equation (15)-(16) with the time derivative left continuos. The same SBP operators will be used for the discrete update of the pressure field. Summation-by-parts is the discrete equivalent of integration-by-parts, which we frequently used for the continuous energy estimates in Section 4. By mimicking the continuous analysis, we will derive energy estimates for the semi-discrete momentum equations.
The boundary conditions are enforced weakly through a simultaneous approximation term (SAT) technique [12] that penalizes the deviation from the exact boundary data g. The matrices
Em = diag(0, . . . , 0, 1), E1 = diag(1, 0, . . . , 0), m = mx, my,
The first derivative SBP operator was introduced in [9, 10] and is
approx-imated by ux ≈ P−1Qu, where
Q + QT
= diag(−1, 0, . . . , 0, 1) = Em− E1 ≡ ˜B. (52)
P is a positive definite diagonal matrix and defines a matrix norm. Operators of order 2n, n ∈ N, in the interior of the domain are combined with bound-ary closures of order of accuracy n. It is possible to design operators with higher order accuracy at the boundary, but this requires P to have nonzero off-diagonal entries which can lead to stability problems [22]. We restrict ourselves to diagonal matrices P since the proofs of stability to be presented later rely on this condition. Subscripts x and y on the operators are used to distinguish between the discretizations in the two coordinate directions, whereas subscripts on vectors denote partial derivatives.
For the approximation of the second derivative SBP operator, introduced in [23, 24], we can either use the first derivative operator twice, or use the
compact approximation uxx ≈ P−1(−M + ˜BD)u, where M + MT ≥ 0, ˜B is
given by (52), and D is a first-derivative approximation at the boundaries, i.e., D = 1 ∆x d1 d2 d3 . . . 1 . .. 1 . . . −d3 −d2 −d1 .
Here, di, i = 1, 2, 3, . . . , are scalar values leading to a consistent
first-derivative approximation at the boundaries.
As a convenient computational tool, we will employ the Kronecker prod-uct, denoted by ⊗ and defined for two matrices B and C, as
B ⊗ C = [B]11C . . . [B]1nC .. . . .. ... [B]m1C . . . [B]mnC .
The following notation for the discretization of the derivative operators and matrices of gPC parameters of the momentum and pressure equations is
used: Dx = P−1x Qx⊗ Iy⊗ IM, Dy = Ix⊗ P−1y Qy⊗ IM, Dxx = P−1x (−Mx+ ˜BxDx) ⊗ Iy⊗ IM, Dyy = Ix⊗ P−1y (−My + ˜ByDy) ⊗ IM, ˜ A(uM ) = diag(A(uM i=1,j=1), . . . , A(u M i=1,j=my), . . . , A(u M i=mx,j=my)), ˜ A(vM ) = diag(A(vM i=1,j=1), . . . , A(v M i=1,j=my), . . . , A(v M i=mx,j=my)), ˜ A(uM x) = diag(A([DxuM]1:M), . . . , A([DxuM](mx+1)(my−1)M+1:mxmy M)), ˜ A(vM y ) = diag(A([DyvM]1:M), . . . , A([DyvM](mx+1)(my−1)M+1:mxmy M)), D(skew)x = 1 2 ˜ A(uM )Dx+ 1 2Dx ˜ A(uM ) − 1 2 ˜ A(DxuM), D(skew)y = 1 2 ˜ A(vM )Dy + 1 2Dy ˜ A(vM ) −1 2 ˜ A(DyvM), ˜
A(µ) = diag(A(µi=1,j=1), . . . , A(µi=mx,j=my)),
where subscripts x and y denote matrices of size mx and my, respectively.
Note that the matrices involving ˜A(.) are block diagonal (but not diagonal).
5.1. Discrete Momentum Equations
The stochastic Galerkin momentum equations (15)-(16) are semi-discretized in space, uM t + D (skew) x u M +D(skew)y uM + DxpM
= DxA(µ)D˜ xuM + DyA(µ)D˜ yuM + SAT(u), (53)
vM t + D (skew) x v M +D(skew)y vM + DypM
= DxA(µ)D˜ xvM + DyA(µ)D˜ yvM + SAT(v), (54)
where the penalty terms SAT(u) and SAT(v) for the boundary conditions will
be given below. The application of the first derivative twice for the viscous second derivative terms yields a wide stencil. Alternatively, one may use the second derivative operator for varying coefficients, see [25], to obtain a more narrow stencil. In order to construct discrete boundary conditions for (53) and (54), the block-diagonal eigenvalue decompositions
are introduced. Let the diagonal matrices containing the eigenvalues (47) and (50) evaluated on the four boundaries (w, e, s, n), be given by
Λw5 = −1 2Λu− 1 4Λ 2 u+ 2I M 1/2 , Λe5 = 1 2Λu− 1 4Λ 2 u+ 2I M 1/2 , Λs5 = −1 2Λv− 1 4Λ 2 v+ 2I M 1/2 , Λn5 = 1 2Λv− 1 4Λ 2 v+ 2I M 1/2 , Λw2 = −1 2Λu− 1 4Λ 2 u+ I M 1/2 , Λe2 = 1 2Λu− 1 4Λ 2 u+ I M 1/2 , Λs2 = −1 2Λv− 1 4Λ 2 v+ I M 1/2 , Λn2 = 1 2Λv− 1 4Λ 2 v+ I M 1/2 . The difference in sign in the above expressions is due to whether the outward normal vector is aligned with the positive velocity direction or not.
The discrete boundary conditions corresponding to (51) are given by
−(XuΛ5XTuu M )w + pMw − A(µ)(DxuM)w = gw(u), (55) (XuΛ5XTuu M )e+ pMe − A(µ)(DxuM)e = ge(u), (56) (XvΛ2XTvu M )s+ A(µ)(DyuM)s = gs(u), (57) −(XvΛ2XTvu M )n+ A(µ)(DyuM)n = gn(u), (58) −(XuΛ2XTuv M )w− A(µ)(DxvM)w = gw(v), (59) (XuΛ2XTuv M )e− A(µ)(DxvM)e = ge(v), (60) −(XvΛ5XTvv M )s+ pMs − A(µ)(DyvM)s = gs(v), (61) (XvΛ5XTvv M )n+ pMn − A(µ)(DyvM)n = gn(v), (62)
for the left, right, lower and upper boundaries respectively. Note the dif-ference in sign in (55)-(62) depending on the orientation of the normal and tangential vectors as displayed in Figure 1.
The boundary conditions are imposed weakly by choosing appropriate
penalty coefficients τα(w) (α = w, e, s, n and w = u, v) on the deviation of
discretized momentum equations (53)-(54) are given by
SAT(u) = τw(u) P−1x E1⊗ Iy⊗ IM
−XuΛw5X T uu M + pM − ˜ A(µ)DxuM − gw(u) + τe(u) P−1x Emx ⊗ Iy ⊗ I M XuΛe5X T uu M + pM − ˜ A(µ)DxuM − ge(u) + τs(u) Ix⊗ P−1y E1⊗ IM XvΛs2X T vu M + ˜A(µ)DyuM − g(u)s + τn(u) Ix⊗ P−1y Emy⊗ I M −XvΛn2X T vu M + ˜A(µ)DyuM − g(u)n , SAT(v) = τw(v) P−1x E1⊗ Iy ⊗ IM −XuΛ w 2XTuv M − ˜ A(µ)DxvM − gw(v) + τe(v) P−1x Emx ⊗ Iy⊗ I M XuΛe2X T uv M − ˜ A(µ)DxvM − ge(v) + τs(v) Ix⊗ P−1y E1⊗ IM −XvΛs5X T vv M + pM − ˜ A(µ)DyvM − gs(v) + τn(v) Ix⊗ P−1y Emy ⊗ I M XvΛn5X T vv M + pM − ˜A(µ)DyvM − gn(v) . 5.2. The Discretized Poisson Equation for the Pressure
The momentum equations require the pressure to be known. The SBP discretization of the stochastic Galerkin pressure equation (19) (temporarily ignoring the boundary conditions) is given by
(Dxx+ Dyy) pM = Dx(A(µx)DxuM)+Dy(A(µx)DyuM)+Dx(A(µy)DxvM)
+Dy(A(µy)DyvM)−A(DxuM)DxuM−2A(DyuM)DxvM−A(DyvM)DyvM.
We use a weak imposition of boundary conditions and solve a linear system of equations for the pressure update in each time-step. The penalty parameters are chosen to ensure that the resulting matrix operator is non-singular. The boundary conditions are given in the form
pM
= g(p)k , (x, y) ∈ ∂Dk, k = w, e, s, n.
The boundary data g(p)k are manufactured by solving (55), (56), (61) and (62)
for the pressure and substituting the DxuM derivative using the dissipative
divergence boundary conditions (30) in discrete form. The result is
g(u)w − A(µ)(DyvM)w+ (XuΛ5XTuu M )w = g(p)w , ge(u)− A(µ)(DyvM)e− (XuΛ5XTuu M )e= g(p)e , g(v)s − A(µ)(DxuM)s+ (XvΛ5XTvv M )s= g(p)s ,
5.3. Time Stability of the Numerical Method
To get a semi-discrete energy estimate for the momentum equations, we
multiply (53) from the left by (uM)TP where P = (P
x⊗Py⊗IM), and add the
transpose of the result to itself. By using P ˜A = ˜AP and the SBP property
(52), we get d dtku Mk2 P+ +1 2(u M )T˜ A(uM )[ ˜Bx⊗ Py⊗ IM]uM + 1 2(u M )T [ ˜Bx⊗ Py⊗ IM] ˜A(uM)uM | {z }
Boundary terms from D(skew)x
+1 2(u M )T˜ A(vM )[Px⊗ ˜By ⊗ IM]uM + 1 2(u M )T [Px⊗ ˜By ⊗ IM] ˜A(vM)uM | {z }
Boundary terms from D(skew)y
−(uM )T P ˜A(uM x ) + ˜A(v M y ) uM | {z } Residual term −2(uM )T (QT x ⊗ Py⊗ IM)pM | {z }
Pressure PDE term
+2(uM )T (Emx ⊗ Py⊗ I M )pM − 2(uM )T (E1⊗ Py ⊗ IM)pM | {z }
Pressure boundary terms
= −2(uM
)T
(QT
x ⊗ Iy⊗ IM)(P−1x ⊗ Py⊗ IM) ˜A(µ)(Qx⊗ Iy ⊗ IM)uM
| {z }
Viscous PDE x-term
−2(uM )T (Ix⊗ QTy ⊗ I M )(Px⊗ P−1y ⊗ I M ) ˜A(µ)(Ix⊗ Qy ⊗ IM)uM | {z }
Viscous PDE y-term
+2(uM )T [(Emx ⊗ Py⊗ I M ) − (E1⊗ Py ⊗ IM)] ˜A(µ)DxuM | {z }
Viscous boundary terms east and west
+2(uM )T(P x⊗ Emy ⊗ I M ) − (Px⊗ E1⊗ IM) ˜ A(µ)DyuM | {z }
Viscous boundary terms north and south
+ SAT. (63)
We group the terms in (63) and the corresponding expression for vM into
SAT/boundary terms, dissipative terms, and residual terms, and analyze them separately.
5.3.1. Boundary Terms and SAT
The matrix (P−1x ⊗ Py⊗ IM) ˜A(µ) = ˜A(µ)(P−1x ⊗ Py⊗ IM) is symmetric
west boundary and the corresponding SAT with zero data, we have BTw(u) = (uM )T (E1⊗ Py ⊗ IM) ˜A(uM)uM + 2(uM)T(E1⊗ Py⊗ IM)pM − 2(uM )T (E1 ⊗ Py ⊗ IM) ˜A(µ)DxuM + 2τw(u)(uM )T (E1 ⊗ Py ⊗ IM) −XuΛw5X T uu M + pM − ˜ A(µ)DxuM .
With the choice τw(u)= −1, the west boundary terms simplify to
BTw(u) = (uM )T (E1⊗ Py ⊗ IM) ˜A(uM) + 2XuΛw5X T u uM = −(uM )T (E1⊗ Py⊗ IM) Xu Λ2u+ 8I 1/2 XT uu M , with − (E1⊗ Py ⊗ IM) Xu Λ2u+ 8I 1/2 XT
u being negative semi-definite.
The analysis is similar for the remaining three boundaries. For the east
boundary and SAT with ge= 0, we obtain
BTe(u) = −(uM )T (Emx⊗ Py⊗ I M ) ˜A(uM )uM − 2(uM )T (Emx⊗ Py⊗ I M )pM + 2(uM )T (Emx⊗ Py ⊗ I M )A(µ)DxuM + 2τe(u)(uM )T (Emx ⊗ Py⊗ I M ) (XuΛ e 5X T uu M + pM − A(µ)D xuM) .
With τe(u) = 1, the east boundary terms become
BTe(u) = (uM )T (Emx⊗ Py ⊗ I M ) − ˜A(uM ) + 2XuΛe5X T u uM = −(uM )T (Emx ⊗ Py⊗ I M ) Xu Λ2u+ 8I 1/2 XT uu M , which again is negative semi-definite.
The south and north boundaries involve no pressure terms. With τs(u) = 1,
we obtain BTs(u) = −(uM )T (Px⊗ E1⊗ IM) Xv Λ2v+ 4I 1/2 XT vu M ,
and τn(u) = −1 yields
BTn(u) = −(uM )T Px⊗ Emy⊗ I M X v Λ2v + 4I 1/2 XT vu M .
gener-τi(w) = ±1 (w = u, v, i = w, e, s, n) and thus do not change with the order
M of gPC expansion. A similar analysis for the vM momentum equation (see
appendix A) yields τw(v) = τs(v) = −1 and τe(v)= τn(v) = 1 for stability.
The boundary terms are collected together and we get
BT ≡ BTw(u)+BTe(u)+BTn(u)+BTs(u)+BTw(v)+BTe(v)+BTn(v)+BTs(v). (64)
5.3.2. Residual Terms The residual terms
RT ≡ (uM )T P ˜A(uM x) + ˜A(v M y ) uM + (vM )T P ˜A(uM x) + ˜A(v M y ) vM + 2 (DxuM + DyvM) T PpM (65) are proportional to the discrete divergence but indefinite. The impact of these terms is assumed to be small, and that will be numerically investigated in Section 6.
5.3.3. Dissipative Terms The viscous PDE terms are
DT ≡ −2 (DxuM) T ˜ A(µ)PDxuM − 2 (DyuM) T ˜ A(µ)PDyuM − 2 (DxvM) T ˜ A(µ)PDxvM − 2 (DyvM) T ˜ A(µ)PDyvM. (66)
The relation (66) is negative semi-definite if ˜A(µ)P is positive definite.
˜
A(µ)P is block diagonal and thus positive definite if and only if all of its diagonal blocks are positive definite. The matrix block corresponding to the
spatial point (xi, yj) is given by [Px]ii[Py]jjA(µ(xi, yj)), for i = 1, . . . , mx
and j = 1, . . . , my. The diagonal entries of Px and Py are all positive, and
A(µ(xi, yj)) is positive definite for all (xi, yj). Hence, ˜A(µ)P is positive
definite.
Assembling the boundary, residual and dissipative terms, the semi-discrete energy estimate becomes
d dt kuMk2 P+ kv Mk2 P = BT + RT + DT,
where BT is defined by (64), RT by (65), and DT by (66). The boundary terms BT and the dissipative terms DT are negative and lead to decay of the
norm of the velocity. The residual terms RT are proportional to the discrete
divergence and expected to be small. Thus, we expect the norm of the
discrete velocity to decay over time and the solution to be time-stable. The discrete energy estimate is analogous to the continuous energy estimate (42). It is important to note that the continuous analysis of well-posedness served as a roadmap to obtain a discrete energy estimate by replacing integration-by-parts by summation-integration-by-parts and using discrete versions of the continuous boundary conditions.
6. Numerical Results
We will consider two exact solutions to the incompressible Navier-Stokes equations and verify the accuracy of the numerical approximation. We are interested in the numerical convergence of statistical properties f of the
so-lution variables and measure the discrete `2 errors. We define the error in f
as (f )` 2 = kf − frefk`2 = ∆x∆y mx X i=1 my X j=1 (f (xi, yj, t) − fref(xi, yj, t))2 !1/2 ,
where a quantity of interest f can be calculated as a function of the gPC coefficients (for instance, if f = V ar(u), then we use the approximation
f ≈ PM
k=2u 2
i). If f is taken to be vector valued, i.e., the vector of gPC
coefficients, we first compute the L2 norm of f . In that case,
(f )
`2 denotes a
double norm with some abuse of notation. The error in the statistics is due to both the discretization error and the truncation of the gPC series.
6.1. Taylor Vortex
Consider the stochastic Taylor vortex model with components exactly satisfying the incompressible Navier-Stokes equations,
u(x, y, t, ξ) = − cos(α) sin(β) exp(−2π2µ(ξ)t) + u∞cos(θ),
v(x, y, t, ξ) = sin(α) cos(β) exp(−2π2µ(ξ)t) + u∞sin(θ),
p(x, y, t, ξ) = −1
4[cos(2α) + cos(2β)] exp(−4π
2µ(ξ)t),
α = π(x − x0− u∞cos(θ)t), β = π(y − y0− u∞sin(θ)t).
The gPC coefficients of the exact solution (67) are given by (9)-(11) and the reference solution is calculated through numerical quadrature approx-imation of the stochastic integrals for the gPC coefficients defined by (2). The analytical solution is deterministic at t = 0 but stochastic for all t > 0. Reference values for statistics, e.g. mean and variance, are computed via numerical integration of functions of the analytical solution. Even though the test case has a periodic solution, we use the more general non-periodic boundary conditions derived previously.
The number of gPC coefficients needed to obtain a given accuracy of the solution is in general not known in advance. We will use relatively low-order gPC representations (M = 3, 4) for the numerical experiments. For the Taylor vortex problem, the gPC coefficients can be calculated up to the accuaracy of the numerical quadrature rule. Figure 2 shows the decay of the spatial norms of the gPC coefficients. The Gauss-Legendre quadrature rule integrates up to tenth-order polynomials exactly. Most of the variance is captured by 3-4 gPC terms due to the fast decay of the coefficients.
1 2 3 4 5 6 10−10 10−8 10−6 10−4 10−2 100 Order of gPC u, µ p, µ u, µ(x) p, µ(x)
Figure 2: Decay of the gPC coefficients of u and p for the Taylor vortex prob-lem. Constant µ ∼ U [0.03, 0.07] and spatially varying µ(x) ∼ U [0.005, 0.05].
6.1.1. Spatially uniform µ
We assume µ = 0.05 + 0.02ξ, ξ ∼ U [−1, 1] and compare the results to the deterministic case µ = 0.05. The gPC coefficients decay quickly for this diffusive problem. The low-order expansion M = 3 leads to sufficiently good stochastic resolution in the sense that the error from the spatial discretization is dominant and we may check the spatial order of convergence by numerical
experiments. The means and standard deviations of the stochastic solution are shown in Figure 3.
0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y 0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y 0.005 0.01 0.015 0.02 0.025 0.03 0.035 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
(a) Mean (left) and standard deviation (right) of velocity component u.
0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y 0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 0.005 0.01 0.015 0.02 0.025 0.03 0.035
(b) Mean (left) and standard deviation (right) of velocity component v.
0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y 0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.005 0.01 0.015 0.02 0.025 0.03
(c) Mean (left) and standard deviation (right) of pressure p.
Figure 3: Means and standard deviations at t = 0.2, mx = my= 121, order
of gPC M = 3, and µ ∼ U [0.03, 0.07].
Table 1 displays the spatial convergence with mesh refinement for the deterministic case. Table 2 shows the convergence for the stochastic Galerkin formulation of the stochastic problem with M = 3 gPC terms. A
fourth-temporal error is negligible compared to the spatial error. We use spatial SBP operators with 8th order of accuracy in the interior and 4th order of accuracy on the boundary. According to the analysis in [26], we may expect 5th order of accuracy globally for the incompressible Navier-Stokes equations and the given choice of spatial operators. The numerical results indicate a slightly higher order of accuracy.
∆x, ∆y (u) `2 q (u) (p) `2 q (p) (φ) `2 q (φ) 2/60 2.2 · 10−5 5.9 5.9 5.9 5.9 1.1 · 10−4 5.7 5.8 5.9 5.9 3.7 · 10−4 6.1 6.2 6.2 6.2 2/80 4.1 · 10−6 2.1 · 10−5 6.4 · 10−5 2/100 1.1 · 10−6 5.8 · 10−5 1.6 · 10−5 2/120 3.7 · 10−7 2.0 · 10−6 5.2 · 10−6 2/140 1.5 · 10−7 7.9 · 10−7 2.0 · 10−6
Table 1: Numerical convergence for the deterministic vortex problem, t = 0.2. `2 error () and order of convergence (q).
We investigate the errors in expectations, standard deviations and the vectors of gPC cofficients. The results for the different statistics are not independent but provide more accessible information. The error in the stan-dard deviation is affected by the truncation error resulting from truncation of the infinite sum in (3) to M terms. For statistics that are independent of the truncated terms, the truncation error will still have some impact on the re-maining terms since the truncation implies that we solve a modified problem. The lower-order terms are expected to be less affected by the stochastic trun-cation. With sufficient stochastic resolution, we can thus expect the error in the expectation (i.e., the first gPC coefficient) to be dominated by the spatial discretization error. Thus, we expect the observed rate of convergence for the expectations in Table 2 to be the most accurate indicator of the actual spa-tial order of convergence since the other statistics are somewhat polluted by the stochastic truncation error. Note also the remarkably high order of con-vergence for the velocity dicon-vergence (which consists of derivatives) in Table 1 and 2.
∆x, ∆y (E(u)) `2 q (E(u)) (Std(u)) `2 q (Std(u)) (u) `2 q (u) 2/40 2.2 · 10−4 5.6 5.9 5.9 5.8 3.2 · 10−5 5.5 6.1 6.2 6.1 2.2 · 10−4 5.6 5.9 5.9 5.7 2/60 2.2 · 10−5 3.5 · 10−6 2.2 · 10−5 2/80 4.1 · 10−6 6.1 · 10−7 4.2 · 10−6 2/100 1.1 · 10−6 1.5 · 10−7 1.2 · 10−6 2/120 3.8 · 10−7 5.1 · 10−8 4.0 · 10−7 ∆x, ∆y (E(p)) `2 q (E(p)) (Std(p)) `2 q (Std(p)) (p) `2 q (p) 2/40 7.5 · 10−4 5.6 5.8 5.9 5.9 7.2 · 10−5 5.8 6.0 5.9 5.0 7.5 · 10−4 5.6 5.8 5.8 5.6 2/60 7.9 · 10−5 7.0 · 10−6 7.9 · 10−5 2/80 1.5 · 10−5 1.3 · 10−6 1.5 · 10−5 2/100 4.0 · 10−6 3.4 · 10−7 4.1 · 10−6 2/120 1.4 · 10−6 1.4 · 10−7 1.5 · 10−6
∆x, ∆y (E(φ))`2 q(E(φ)) (Std(φ))
`2 q (Std(φ)) (φ)`2 q(φ) 2/40 4.0 · 10−3 5.7 6.1 6.2 6.2 1.1 · 10−3 5.5 6.1 6.3 6.3 4.1 · 10−3 5.7 6.1 6.2 6.2 2/60 3.9 · 10−4 1.2 · 10−4 4.1 · 10−4 2/80 6.8 · 10−5 2.1 · 10−5 7.2 · 10−5 2/100 1.7 · 10−5 5.2 · 10−6 1.8 · 10−5 2/120 5.6 · 10−6 1.7 · 10−6 5.8 · 10−6
Table 2: Numerical convergence for the stochastic vortex problem, M = 3, t = 0.2. Error and order of convergence q for mean value (E(.)), standard deviation (Std(.)), and the vector of gPC coefficients (u/p/φ).
6.1.2. Spatially varying µ
To investigate the performance of the present formulation for problems with spatially varying viscosity, we next consider the solution to (4)-(6),
where the viscosity is a random field µ(x, ξ) = ¯µ+σµk−1/2cos(kπx) cos(kπy)ξ.
We set ¯µ = 0.0275, and σµ= 0.0225 (i.e., µ varies in the interval [0.005, 0.05]).
Using the method of manufactured solutions [27, 28], (4)-(5) and (8) are mod-ified by adding appropriate source terms to satisfy the solution (67). Results on spatial convergence are shown in Table 3.
tic space and we set k = 1, but the extension to multiple random variables is straightforward by using a series of terms with different integer k. This
implies a Karhunen-Lo`eve-like expansion with, respectively, eigenvalues and
eigenfunctions of the form k−1 and cos(kπx) cos(kπy) for the covariance
function of µ(x, ξ). ∆x, ∆y (E(u)) `2 q (E(u)) (Std(u)) `2 q (Std(u)) (u) `2 q (u) 2/40 3.1 · 10−4 5.7 5.9 5.9 5.9 7.8 · 10−5 4.5 5.0 5.2 5.2 3.2 · 10−4 5.6 5.7 5.7 5.6 2/60 3.0 · 10−5 1.3 · 10−5 3.3 · 10−5 2/80 5.5 · 10−6 3.1 · 10−6 6.4 · 10−6 2/100 1.5 · 10−6 9.7 · 10−7 1.8 · 10−6 2/120 5.1 · 10−7 3.8 · 10−7 6.5 · 10−7
∆x, ∆y (E(p))`2 q(E(p)) (Std(p))
`2 q (Std(p)) (p)`2 q(p) 2/40 9.0 · 10−4 5.5 5.8 5.9 5.9 1.5 · 10−4 4.7 5.3 5.5 5.5 9.1 · 10−4 5.4 5.7 5.8 5.9 2/60 9.8 · 10−5 2.2 · 10−5 1.0 · 10−4 2/80 1.9 · 10−5 4.8 · 10−6 1.9 · 10−5 2/100 5.0 · 10−6 1.4 · 10−6 5.2 · 10−6 2/120 1.7 · 10−6 5.1 · 10−7 1.8 · 10−6 ∆x, ∆y (E(φ)) `2 q (E(φ)) (Std(φ)) `2 q (Std(φ)) (φ) `2 q (φ) 2/40 7.3 · 10−3 5.5 5.9 6.0 6.1 1.8 · 10−3 4.5 5.1 5.4 5.6 7.5 · 10−3 5.4 5.8 5.9 6.0 2/60 7.9 · 10−4 3.0 · 10−4 8.4 · 10−4 2/80 1.4 · 10−4 6.9 · 10−5 1.6 · 10−4 2/100 3.8 · 10−5 2.1 · 10−5 4.3 · 10−5 2/120 1.3 · 10−5 7.5 · 10−6 1.5 · 10−5
Table 3: Numerical convergence for the stochastic vortex problem, M = 4, t = 0.2. Error and order of convergence for mean value (E(.)), standard deviation (Std(.)), and the vector of gPC coefficients (u/p/φ). q denotes the observed order of convergence.
The decay of the velocity divergence with mesh refinement verifies that the residual terms (65) proportional to the discrete divergence are indeed
small. As in the case of spatially uniform µ, remarkably high-order of con-vergence is observed for the velocity dicon-vergence also for the case of spatially varying µ.
6.2. Poiseuille Flow
We next consider laminar flow between two plates where the fluid is ini-tially at rest, u = v = 0. A constant pressure gradient is applied in the horizontal direction with a parabolic inflow profile at the west boundary. The steady state solution with unit density and deterministic boundary con-ditions is the parabolic Poiseuille profile,
u(x, y, t, ξ) = u∞(1 − (y − 1)2),
v(x, y, t, ξ) = 0,
p(x, y, t, ξ) = 1 − 2u∞µ(ξ)x.
The zero Neumann condition at the outlet is replaced with the east boundary conditions (56) and (60) evaluated with steady state data.
Figure 4 shows the determinitic solution at t = 10 on a coarse mesh,
mx = my = 24. The numerical solution is close to the steady state solution
and the difference is on the order 10−6 for both the solution variables and
the divergence. For this problem, the boundary data and the initial function are not consistent and the divergence is initially non-zero but converges to zero with time. The decay of the velocity divergence with time, even with inconsistent initial data, verifies that the residual terms (65) proportional to the discrete divergence are indeed small. As before, we do not impose zero divergence strongly in the interior but provide dissipative boundary data. As a result, the error in the divergence has been reduced by seven orders of magnitude at t = 10 for the coarse mesh problem of Figure 4.
For the stochastic Poiseuille flow with µ = 0.05 + 0.02ξ, ξ ∼ U [−1, 1], we observe a similar convergence to a divergence-free solution, but for the case M = 3, the error is an order of magnitude larger. The errors in the standard deviation of the solution variables and the divergence are shown in Figure 5. The decay of the velocity divergence over time is shown in Figure 6 for two different mesh sizes, ∆x = ∆y = 0.1 and ∆x = ∆y = 0.05. The decay in the
norm of all the included gPC coefficients (φ)`
2 is plotted as a function of time
(a) Velocity component u. (b) Velocity component v. (c) Pressure p. 0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 x 10−5 (d) Divergence φ. Figure 4: Deterministic Poiseuille solution at t = 10, mx= my= 24.
(a) Velocity component u. 0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y 2 4 6 8 x 10−6 (b) Velocity component v. (c) Pressure p. 0 0.5 1 1.5 2 0 0.5 1 1.5 2 x y 2 4 6 8 10 12 14 16 x 10−5 (d) Divergence φ.
Figure 5: Error in the standard deviations for Poiseuille flow at t = 10, mx= my = 24. Stochastic µ = 0.05 + 0.02ξ, ξ ∼ U [−1, 1].
0 5 10 15 20 10−10 10−8 10−6 10−4 10−2 100 Time ∆x = 0.05 ∆x = 0.1
Figure 6: Decay of the discrete `2norm of the gPC coefficients of the velocity
divergence for the Poiseuille flow problem over time. For each spatial point we compute the L2 norm of the local vector of gPC coefficients. Constant
µ ∼ U [0.03, 0.07].
The Poiseuille flow problem and the Taylor vortex problem are comple-mentary. The initial velocity field of the Poiseuille flow problem is incon-sistent with the boundary data, which leads to relatively large values of the divergence field. Asymptotically, the numerical solution converges to the true steady-state solution with zero velocity divergence. For the Taylor vor-tex problem, we provided consistent initial data and boundary data satisfying the zero divergence condition. In that case, we could verify high-order spatial convergence for any given time.
7. Conclusions
We have presented a numerical method for the stochastic Galerkin for-mulation of the incompressible Navier-Stokes equations that enforces incom-pressibility weakly. The analysis relies on the use of summation-by-parts op-erators for different orders of accuracy. The result is a stable and high-order accurate numerical method where the divergence is zero to the approximation order of the scheme.
Stability of numerical methods is a prerequisite for convergence to the true solution, but non-trivial to ensure for non-linear problems. For stochas-tic Galerkin formulations, the situation is further complicated and choosing boundary conditions based on physical arguments is not sufficient to ensure
a numerically stable method. The boundary conditions are stochastic, non-periodic in space, and time-dependent. The number of boundary conditions depends on the order of the gPC expansion, but does not change over time. In order to derive a stable scheme for the semi-discrete stochastic Galerkin incompressible Navier-Stokes system, the continuous formulation and the de-terministic scheme serve as a roadmap. A stable stochastic Galerkin method is highly non-trivial to find without first performing the analysis for the deterministic and continuous problems. Once the appropriate general for-mulation of the numerical method has been found, design parameters are straightforward to choose.
The method presented here leads to high-order accuracy for all variables including the velocity divergence. The high-order convergence for the di-vergence is notable since it would normally be lower for quantities that are derivative approximations of the solution variables. Numerical results are presented with discretization operators of eighth-order accuracy in the inte-rior spatial domain, and fourth-order accuracy on the boundaries. With a sufficiently accurate stochastic representation, we observe fifth-order global convergence.
8. Acknowledgements
This material is based upon work of Alireza Doostan supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, under Award Number DE-SC0006402. The authors would like to thank Anna Nissen for analysis of stability leading to the choice of penalty parameters for the pressure updates.
Appendix A Stability Analysis
A.1 Energy estimate for the vM momentum equation
The energy estimate for vM is analogous to the energy estimate for uM.
Multiply (54) from the left by (vM)TP to obtain
d dtkv M k2P+ (A.1) +1 2(v M )T˜ A(uM )[ ˜Bx⊗ Py ⊗ IM]vM + 1 2(v M )T [ ˜Bx⊗ Py⊗ IM] ˜A(uM)vM | {z }
Boundary terms fromD(skew)x
+1 2(v M )TA(v˜ M )[Px⊗ ˜By ⊗ IM]vM + 1 2(v M )T [Px⊗ ˜By⊗ IM] ˜A(vM)vM | {z }
Boundary terms from D(skew)y
−(vM )T P ˜A(uM x) + ˜A(v M y ) vM | {z } Residual term −2(vM )T (Px⊗ QTy ⊗ I M )pM | {z }
Pressure PDE term
+2(vM )T (Px⊗ Emy⊗ I M )pM − 2(vM )T (Px⊗ E1⊗ IM)pM | {z }
Pressure boundary terms
= −2(vM
)T
(QTx ⊗ Iy ⊗ IM)(P−1x ⊗ Py⊗ IM) ˜A(µ)(Qx⊗ Iy ⊗ IM)vM
| {z }
Viscous PDE x-term
−2(vM )T(Ix⊗ QTy ⊗ I M )(Px⊗ P−1y ⊗ I M ) ˜A(µ)(Ix⊗ Qy ⊗ IM)vM | {z }
Viscous PDE y-term
+2(vM
)T [(Emx ⊗ Py⊗ I
M
) − (E1 ⊗ Py ⊗ IM)] ˜A(µ)DxvM
| {z }
Viscous boundary terms east and west
+2(vM )T (Px⊗ Emy ⊗ I M ) − (Px⊗ E1⊗ IM) ˜ A(µ)DyvM | {z }
Viscous boundary terms north and south
Set g(v)w = g(v)e = g(v)s = gn(v)= 0, and collect the terms boundary-wise. For
the west boundary, we have
BTw(v) = (vM )T (E1⊗ Py⊗ IM) ˜A(uM)vM − 2(vM )T (E1⊗ Py⊗ IM) ˜A(µ)DxvM + 2τw(v)(v M )T (E1⊗ Py⊗ IM) × × −Xu " −1 2Λu− 1 4Λ 2 u+ I 1/2# XT uv M − ˜ A(µ)DxvM ! . (A.2)
With the choice τw(v)= −1, (A.2) simplifies to
BTw(v) = −(vM )T (E1⊗ Py⊗ IM) Xu Λ2u+ 4I 1/2 XT uv M . For the east boundary,
BTe(v) = −(vM )T (Emx ⊗ Py⊗ I M ) ˜A(uM )vM + 2(vM )T (Emx ⊗ Py⊗ I M ) ˜A(µ)DxvM + 2τe(v)(v M )T (Emx ⊗ Py ⊗ I M ) × × Xu " 1 2Λu− 1 4Λ 2 u+ I 1/2# XT uv M − ˜ A(µ)DxvM ! , which simplifies to BTe(v)= −(vM )T (Emx ⊗ Py ⊗ I M ) Xu Λ2u+ 4I 1/2 XT uv M with τe(v)= 1.
The south boundary terms are given by
BTs(v) = (vM )T (Px⊗ E1⊗ IM) ˜A(vM)vM + 2(vM )T (Px⊗ E1⊗ IM) h pM − ˜ A(µ)DyvM i + 2τs(v)(vM )T (Px⊗ E1⊗ IM) × × −Xv " −1 2Λv− 1 4Λ 2 v+ 2I 1/2# XT vv M + pM − ˜ A(µ)DyvM ! ,
which for τs(v) = −1 becomes
Finally, the north boundary terms are BTn(v) = −(vM )T Px⊗ Emy ⊗ I M ˜ A(vM )vM +2(vM )T (Px⊗Emy⊗I M )h−pM + ˜A(µ)DyvM i +2τn(v)(vM )T Px⊗ Emy⊗ I M × × Xv " −1 2Λv− 1 4Λ 2 v+ 2I 1/2# XT vv M + pM − ˜ A(µ)DyvM ! . Setting τn(v) = 1 yields BTn(v) = −(vM )T Px⊗ Emy ⊗ I M X v Λ2v+ 8I 1/2 XT vv M . References
[1] O. P. Le Maˆıtre, O. M. Knio, Spectral Methods for Uncertainty Quan-tification, 1st Edition, Springer, 2010.
[2] R. G. Ghanem, P. D. Spanos, Stochastic finite elements: a spectral approach, Springer-Verlag, New York, 1991.
[3] D. Xiu, Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton University Press, 2010.
[4] O. P. Le Maˆıtre, O. M. Knio, H. N. Najm, R. G. Ghanem, A stochastic projection method for fluid flow: I. basic formulation, J. Comput. Phys. 173 (2) (2001) 481–511. doi:10.1006/jcph.2001.6889.
[5] O. P. Le Maˆıtre, M. T. Reagan, H. N. Najm, R. G. Ghanem, O. M. Knio, A stochastic projection method for fluid flow: II. random process, J. Comput. Phys. 181 (1) (2002) 9–44. doi:10.1006/jcph.2002.7104. [6] D. Xiu, G. E. Karniadakis, Modeling uncertainty in flow simulations via
generalized polynomial chaos, J. Comput. Phys. 187 (1) (2003) 137–167. doi:10.1016/S0021-9991(03)00092-5.
[7] D. Xiu, D. Lucor, C. H. Su, G. E. Karniadakis, Stochastic modeling of Flow-Structure interactions using generalized polynomial chaos, J. Fluids Eng. 124 (1) (2002) 51–59. doi:10.1115/1.1436089.
[8] O. P. Le Maˆıtre, M. T. Reagan, B. J. Debusschere, H. N. Najm, R. G. Ghanem, O. M. Knio, Natural convection in a closed cavity under stochastic non-Boussinesq conditions., SIAM J. Scientific Computing 26 (2) (2004) 375–394. doi:10.1137/S1064827503422853.
[9] H.-O. Kreiss, G. Scherer, Finite element and finite difference methods for hyperbolic partial differential equations, in: Mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, New York, 1974, pp. 179–183.
[10] B. Strand, Summation by parts for finite difference
approxi-mations for d/dx, J. Comput. Phys. 110 (1) (1994) 47–67.
doi:10.1006/jcph.1994.1005.
[11] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time dependent problems and difference methods, Pure and applied mathematics, J. Wiley, New York, Chichester, Brisbane, 1995. doi:10.1002/9781118548448.
[12] M. H. Carpenter, D. Gottlieb, S. Abarbanel, Time-stable boundary con-ditions for finite-difference schemes solving hyperbolic systems: method-ology and application to high-order compact schemes, J. Comput. Phys. 111 (2) (1994) 220–236. doi:http://dx.doi.org/10.1006/jcph.1994.1057.
[13] J. Nordstr¨om, K. Mattsson, C. Swanson, Boundary conditions
for a divergence free velocity-pressure formulation of the
Navier-Stokes equations, J. Comput. Phys. 225 (1) (2007) 874–890.
doi:10.1016/j.jcp.2007.01.010.
[14] D. Xiu, G. E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2) (2002) 619–644. doi:http://dx.doi.org/10.1137/S1064827501387826.
[15] R. Askey, J. A. Wilson, Some Basic Hypergeometric Orthogonal Poly-nomials That Generalize Jacobi PolyPoly-nomials, no. 319 in Memoirs of the American Mathematical Society, American Mathematical Society, 1985.
[16] B. K. Alpert, A class of bases in L2 for the sparse
representa-tions of integral operators, SIAM J. Math. Anal. 24 (1993) 246–262. doi:10.1137/0524016.
[17] X. Wan, G. E. Karniadakis, An adaptive multi-element generalized poly-nomial chaos method for stochastic differential equations, J. Comput. Phys. 209 (2005) 617–642. doi:10.1016/j.jcp.2005.03.023.
[18] P. Pettersson, A. Doostan, J. Nordstr¨om, On stability and monotonicity
requirements of finite difference approximations of stochastic conserva-tion laws with random viscosity, Comput. Method. Appl. M. 258 (2013) 134–151. doi:10.1016/j.cma.2013.02.009.
[19] J. Silvester, Determinants of block matrices, Math. Gaz. 84 (501) (2000) 460–467. doi:10.2307/3620776.
[20] D. Gottlieb, D. Xiu, Galerkin method for wave equations with uncertain coefficients, Comm. Comput. Phys. 3 (2) (2008) 505–518.
[21] H.-O. Kreiss, J. Lorenz, Initial-Boundary Value Problems and
the Navier-Stokes Equation, Classics in Applied
Mathemat-ics, Society for Industrial and Applied Mathematics, 2004.
doi:http://dx.doi.org/10.1137/1.9780898719130.
[22] M. Sv¨ard, On coordinate transformations for summation-by-parts
opera-tors, J. Sci. Comput. 20 (1) (2004) 29–42. doi:10.1023/A:1025881528802.
[23] M. H. Carpenter, J. Nordstr¨om, D. Gottlieb, A stable and conservative
interface treatment of arbitrary spatial accuracy, J. Comput. Phys. 148 (1999) 341–365. doi:10.1006/jcph.1998.6114.
[24] K. Mattsson, J. Nordstr¨om, Summation by parts operators for finite
dif-ference approximations of second derivatives, J. Comput. Phys. 199 (2) (2004) 503–540. doi:10.1016/j.jcp.2004.03.001.
[25] K. Mattsson, Summation by parts operators for finite difference approx-imations of second-derivatives with variable coefficients, J. Sci. Comput. 51 (2012) 650–682. doi:10.1007/s10915-011-9525-z.
[26] M. Sv¨ard, J. Nordstr¨om, On the order of accuracy for difference
approx-imations of initial-boundary value problems, J. Comput. Phys. 218 (1) (2006) 333–352. doi:10.1016/j.jcp.2006.02.014.
[27] P. J. Roache, Verification of codes and calculations, AIAA Journal 36 (5) (1988) 696–702. doi:10.2514/2.457.
[28] K. Salari, P. Knupp, Code verification by the method of manufactured solutions, Tech. Rep. SAND 2000-1444, Sandia National Laboratories, Albuquerque, NM (2000).
