1 Introduction

In the development of computational methods for compressible fluid flows, there are a lot of successful achievements in the finite difference, finite volume, and discontinuous Galerkin methods and their variants, e.g., see Ref. [1] and the references therein. One of the central ingredients in the design of these methods is the construction of numerical fluxes. The genuine multi-dimensionality becomes a key criterion evaluating the resulting scheme for practical applications. In the literature, there are many direct genuinely multidimensional methods^[2-7], e.g., the kinetic flux vector splitting (KFVS) method. The difficulty in the design of genuinely multidimensional methods lies at least in the following aspects:

(ⅰ) The balance laws for compressible fluid flows tell that the change rates of the total mass, momentum, and energy over any control volume are equal to the corresponding flux across the boundary of the control volume. The laws on their own do not explain the role of the transversal effect in the flow patterns even though the derivative velocity field can have longitudinal and transversal decomposition. In parallel, the direct numerical expression of the balance laws is the finite volume formula. As far as the resulting partial differential equations are adopted, the Gauss-Green theorem provides the longitudinal (normal) flux, but the transversal variation relative to the computational cell boundaries is not included.

(ⅱ) For a single advection equation, the multidimensional upwinding, as pointed out in Ref. [8], is useful to provide all information for the flux. In contrast, for a fluid dynamical system, all information propagates along the characteristic cones (not only bicharacteristics). Therefore, the information on the boundary of the control volumes should be tracked backward along the characteristic cones, as done in Refs. [3] and [4], which is not easy to achieve for nonlinear systems, particularly when discontinuities are present in the solutions.

(ⅲ) In practical computations, the dimensional operator splitting technique^[9] is often used to achieve the multi-dimensionality. However, this technique has limitation when complex geometry domains are involved, and high order extension is hard to be carried out.

This paper aims to show, using the classical wave system, how the transversal variation plays a role in the computation of compressible fluid flows. We adopt the multi-step Runge-Kutta (RK) line method^[10] and the acoustic generalized Riemann problem (GRP) method to compare the differences. It turns out from the numerical results that, even superficially with the same convergence rate, numerical errors are substantially improved when the transversal effect is reflected precisely in the computation.

We organize the remanent paper as follows. In Section 2, we introduce mutidimensional numerical schemes for a linear advection equation and a wave system. We focus on the differences between the two-dimensional (2D) GRP method and the multi-step method, and highlight the transversal effect on the numerical flux. We display the numerical performance of the GRP method with the 2D GRP solver in Section 3 by comparing with the RK method and the GRP method with the original planar one-dimensional (1D) GRP solver. In Section 4, we summarize the conclusions that we have obtained.

2 Numerical schemes for multidimensional conservation laws 2.1 A brief review

There are a lot of multidimensional numerical schemes for hyperbolic conservation laws^[2-4]. Examples are as follows:

(ⅰ) Flux splitting methods

With the Strang splitting approach^[9], a multidimensional problem is decomposed into several 1D problems, and then 1D numerical methods are used. Particularly, 1D Godunov-type flux solvers are adopted with great success^[11]. However, such a splitting often brings errors for practical problems. For example, we solve the problem with

Then, this problem can be solved exactly if and only if the operators can commute, i.e., ^[8]. As far as the computational domains with complex geometry are involved, such methods need great modification to suit for practical computations. There is no intrinsic extension into higher order methods yet.

(ⅱ) Unsplit methods

It is natural to design unsplit methods for computing compressible fluid flows^[2], thanks to the physical background or integral form (balance law). Most of those methods assume that information is exchanged between the neighbor cells by 1D waves traveling normal to the interface, subject to given data at each initial time step t=t_n. Due to the ignorance of the transversal effect on the cell interfaces, we cannot see their advantages over dimensional splitting methods when computing shear waves, even though unstructured meshes are used^[8].

In this paper, we will demonstrate how important the transversal effect is included in the flux.

2.2 A linear advection equation

Consider a linear advection equation

(1)

over any computational control volume, where a is a constant vector. Here, we restrict to the 2D case with rectangular meshes. Our arguments are certainly applied to other cases, even with unstructured meshes. Let

Then, we let

be a computational cell, Δx and Δy be the spatial increments, and Δt be a time increment. The average of ϕ at t_n over Ω_ij is denoted as follows:

(2)

This value can be updated by

(3)

where the numerical fluxes are numerical fluxes across the interfaces and during the temporal interval [t_n, t_n+1], respectively. As in Ref. [2], we recast this into the finite volume framework with the following two steps:

(ⅰ) Solution evolution

Given the data , we solve Eq. (1) approximately to obtain the solution ϕ(x, t) for t_n < t < t_n+1 and incorporate it properly into Eq. (3), where is a persistent functional space and is often chosen to be a space of piecewise polynomials.

(ⅱ) Data reconstruction

Use a certain reconstruction technology, e.g., the weighted essentially non-oscillatory (WENO) scheme^[12], to construct the data at the next time step t=t_n+1 as follows:

(4)

In this presentation, we take ϕⁿ⁺¹(x) as a piecewise linear distribution of the form

(5)

for the second-order accurate schemes.

We will not discuss the data reconstruction step, but concentrate on the solution evolution step. For this purpose, it is natural to make approximations for numerical fluxes, i.e.,

(6)

One way is, as shown in Ref. [2], that, within the second-order accuracy, the numerical flux is taken as follows:

(7)

Then, we proceed to obtain

(8)

The instantaneous values and should be taken upwind. Here, we assume a>0 and b>0 for simplicity. Then, we have

(9)

(10)

This is consistent with the direct method of characteristics

(11)

by recalling the linear distribution of data in Eq. (5). The transversal effect is reflected by including the term in the flux.

An alternative way is to adopt a multi-step method, e.g., the RK method, to include the transversal effect. To be more precise, the flux is approximated as follows:

(12)

where is a Gaussian point, and ω_ℓ is the corresponding weight. The value is obtained by some exact or approximate Riemann solvers

(13)

We find that the transversal terms are canceled out as follows:

(14)

It turns out that the transversal effect has to be reflected through the temporal iteration. This is very unlike the single step method presented above. This fact is highlighted for the fluid dynamical system. In the following content, we will discuss the linear wave system to show how to include the transversal effect.

2.3 Wave system

Consider the linear wave system

(15)

This system is equivalent to the wave equation

(16)

Regarding (u, v) as the velocity field, the vorticity ω=v_x - u_y is preserved, i.e.,

(17)

Therefore, a good scheme should preserve such a property.

For the clarity of presentation below, we denote U =(p, u, v)^T and u=(u, v). Then, the fluxes are

The system (15) is written in the conservation form as follows:

(18)

The numerical flux across any interface at any time t is

(19)

where is the normal of the interface . Due to the rotational invariance of Eq. (15), we can assume

for the simplicity of presentation but without loss of generality. Then, we have

(20)

Relative to the interface

the transversal effect is exhibited in terms of the tangential component of the velocity v. It is clear that the transversal effect cannot be directly included, no matter how any first-order exact or approximate Riemann solver is adopted. Therefore, let us think of high order flux solvers, and use the GRP solver, e.g., by assuming the piecewise initial data at t=t_n,

(21)

Then, we use the governing equations in Eq. (15) to see

(22)

so that the v-variation is included in the flux.

Practically, we write the system (18) in the form as follows:

(23)

with

(24)

Then, a 2D GRP solver is stated as follows:

The values satisfy

(25)

(26)

where λ_l (l=1, 2, 3) are the eigenvalues of A, L is the left eigenmatrix of A, and

The justification of this proposition is put in Appendix A for completeness. Here, it is observed that the temporal derivative of p is

(27)

Even though there is no obvious term for v in (see Eq. (20)), the numerical flux by the GRP solver is expressed as follows:

(28)

in which the p-term is

(29)

However, due to the discontinuity of v_y across the interface, we cannot use Eq. (27) directly. Instead, we approximate the midpoint value by solving the following Riemann problem on the left side in the second order. We solve the Riemann problem at as follows:

(30)

It turns out that the mid-point value for p is

(31)

We see that the transversal effects are included.

The same as the advection case, if multistage time evolution approaches are used to devise high-order numerical methods, it is common to use the Gauss numerical integrals for approximating numerical fluxes in order to guarantee the second-order accuracy in space. In detail, we evaluate the flux by

(32)

By the Riemann solver at the interface , we have

(33)

The transversal effect is not reflected for the wave system, i.e.,

(34)

The numerical examples in the next section show that the transversal effect at the cell interfaces will affect the numerical solution.

3 Numerical results

In this section, we present some examples for the wave system and the Euler equations. For all examples, the computational grid is the structured grid with the cell width Δx=Δy=h. The time step Δt is restricted by Δt = νh/c with the Courant-Friedrichs-Lewy (CFL) number ν = 0.5. Then, we present the numerical results obtained by using the RK method and the GRP methods. We abbreviate RK2 for the two-stage second-order RK method, GRP2D for the second-order GRP method with the 2D GRP solver, and GRP1D for the second-order GRP method with the original planar 1D GRP solver^[13]. All these methods consist of two steps, i.e., solution evolution and data reconstruction. In the following subsection, we will determine how the data are reconstructed.

3.1 Data reconstruction

For the second-order case, we use the piecewise reconstruction with the minmod limiter as introduced in Ref. [14]. For the GRP2D, we use the reconstruction with the piecewise constant slope calculated by

(35)

with

(36)

(37)

where α ∈ [1, 2], and it is set as 1.9 here. The generalized minmod function is

(38)

Since the term in Eq. (27) is discontinuous at , the limiting values on the left-hand and right-hand sides of the interface are different, which leads to the fact that

(39)

Therefore, we need to update the slope .

3.2 Numerical examples for the wave system

Here, we devise two examples to test the accuracy of schemes and the transversal effect. The constant sonic speed c₀ ≡ 1. The following numerical examples show that the GRP2D performs the best.

Example 1 Periodic waves

The first example is an initial value problem in Ref. [15]. We take the initial data as follows:

(40)

In this example, the exact solutions are

(41)

(42)

The computational domain [-2, 2] × [-2, 2] is divided into N× N cells for N=40, 80, 160, 320, 640, and is extended periodically in both directions.

The L₁ errors of the numerical results are displayed in Table 1. This table shows clearly that the solutions of the GRP2D and RK2 reach second-order accuracy, while the GRP1D cannot reach second-order accuracy. According to the results in Table 1, the GRP2D performs better than the RK2 and GRP1D in terms of the L₁ errors of U, and has the best rate of convergence. The errors by the GRP2D are almost 1/9 of those by the RK2.

For 80× 80 cells, the local numerical solutions of p and u at the center-line y = 0 are shown in Fig. 1. We can see that the solution of the GRP2D is the closest to the exact solution. This example illustrates that the 2D GRP solver with the transversal effect can greatly improve the convergence error on the solution evolution step.

Example 2 2D Riemann problem

This example is a 2D Riemann problem in Ref. [4]. We consider the initial data including four pieces of discontinuous data as follows:

(43)

The computational domain [-1, 1] × [-1, 1] is divided into 400× 400 cells, whose boundaries are non-reflective.

The solutions obtained by the GRP2D, RK2, and GRP1D are at the final time t=0.4. The two discontinuities at the line x=-y propagate in the positive and negative normal directions of this line with the speed of c₀. For the exact solution of the 2D Riemann problem outside the subsonic domain at the origin, we can refer to Ref. [16]. In Fig. 2, we show the comparison of the final solutions by the GRP2D, RK2, and GRP1D. We can see clearly that the GRP2D resolves the discontinuities, propagating at the characteristic velocity ±c₀, better than the RK2, and resolves the stationary discontinuity, better than the GRP1D.

More minutely, we present the local solution of p at y=0 in Fig. 3. We can observe that the GRP2D produces no oscillation and has better resolution of discontinuities than the RK2 and GRP1D. This example illustrates that the 2D GRP solver with the transversal effect can improve the resolution of discontinuities intersecting with the cell interfaces.

4 Conclusions

In this paper, we use the 2D GRP solver to highlight the importance of the transversal effect for the design of numerical schemes for fluid dynamical problems. From the comparison of the RK approach and the original planar 1D GRP solver, we can see that the 2D GRP solver performs the best in the accuracy test and resolution of discontinuity. This work is done for the wave equation system. We think that the conclusion is instructive for other systems of hyperbolic conservation laws including the velocity field, e.g., the Euler equations and the Navier-Stokes equations. In particular, as far as turbulent flows are simulated, the transversal effect becomes extremely prominent^[17-24]. In this direction, we will further carry out our research in the near future.

Appendix A

The Riemann solution is the local solution to the following Riemann problem of Eq. (18):

(A1)

Diagonalizing Eq. (23), we obtain

(A2)

where L is the left eigenmatrix of A defined by

(A3)

and Λ is the diagonal matrix with entries of eigenvalues, i.e.,

(A4)

According to Ref. [23], the l-th element of the Riemann solution vector is

(A5)

In the limit of t→t_n+0, we have

(A6)

Thus,

(A7)

where A⁺=L^-1Λ⁺L, and A^-=L^-1Λ^-L.

Taking the spatial derivative of Eq. (A2), we have

(A8)

The solution of these linearized equations, corresponding to the 2D GRP at the cell interfaces , is obtained by solving the Riemann problem as follows:

(A9)

Similarly, we have

(A10)

Since λ₂=0, the second element of , is discontinuous at the interface . We have , and thus