A new high-resolution and genuinely multidimensional numerical method for solving conservation laws is being developed. It was designed to avoid the limitations of the traditional methods, and was built from ground zero with extensive physics considerations. Nevertheless, its foundation is mathematically simple enough that one can build from it a coherent, robust, efficient and accurate numerical framework.

Two basic beliefs that set the new method apart from the established methods are at the core of its development. The first belief is that, in order to capture physics more efficiently and realistically, the modeling focus should be placed on the original integral form of the physical conservation laws, rather than the differential form. The latter form follows from the integral form under the additional assumption that the physical solution is smooth, an assumption that is difficult to realize numerically in a region of rapid change, such as a boundary layer or a shock. The second belief is that, with proper modeling of the integral and differential forms themselves, the resulting numerical solution should automatically be consistent with the properties derived from the integral and differential forms, e.g., the jump conditions across a shock and the properties of characteristics. Therefore a much simpler and more robust method can be developed by not using the above derived properties explicitly.

Specifically, to capture physics as fully as possible, the method requires that: (i) space and time be unified and treated as a single entity; (ii) both local and global flux conservation in space and time be enforced; and (iii) a multidimensional scheme be constructed without using the dimensional-splitting approach, such that multidimensional effects and source terms (which are scalars) can be modeled more realistically.

To simplify mathematics and broaden its applicability as much as possible, the method attempts to use the simplest logical structures and approximation techniques. Specifically, (i) it uses a staggered space-time mesh such that flux at any interface separating two conservation elements can be evaluated internally in a simpler and more consistent manner, without using a separate flux model; (ii) it does not use many well-established techniques such as Riemann solvers, flux splittings and monotonicity constraints such that the limitations and complications associated with them can be avoided; and (iii) it does not use special techniques that are not applicable to more general problems.

Furthermore, triangles in 2D space and tetrahedrons in 3D space are used as the basic building blocks of the spatial meshes, such that the method (i) can be used to construct 2D and 3D non-dissipative schemes in a natural manner; and (ii) is compatible with the simplest unstructured meshes.

Note that while numerical dissipation is required for shock capturing, it may also result in annihilation of small disturbances such as sound waves and, in the case of flow with a large Reynolds number, may overwhelm physical dissipation. To overcome this difficulty, two different and mutually complementary types of adjustable numerical dissipation are introduced in the present development.

Since its inception in 1991 [1], the space-time conservation element and solution element method [1-32] has been used to obtain highly accurate numerical solutions for flow problems involving shocks, rarefaction waves, acoustic waves, vortices, ZND detonation waves, shock/acoustic waves/vortices interactions, dam-break and hydraulic jump. This article is the first of a series of papers that will provide a systematic and up-to-date description of this new method (hereafter it may be referred to abbreviatedly as the space-time CE/SE method or simply as the CE/SE method). To answer frequently-asked questions and clarify possible misconceptions, we shall begin this paper with (i) an overall view of the CE/SE method and its capabilities, and (ii) an extensive comparison of the basic concepts used by the CE/SE method with those used by other methods.

Currently, the field of computational fluid dynamics (CFD) represents a diverse collection of numerical methods, with each of them having its own limitations. Generally speaking, these methods were originally introduced to solve special classes of flow problems. Development of the CE/SE method is motivated by a desire to build a brand new, more general and coherent numerical framework that avoids the limitations of the traditional methods.

The new method is built on a set of design principles given in [2]. They include: (i) To enforce both local and global flux conservation in space and time, with flux evaluation at an interface being an integral part of the solution procedure and requiring no interpolation or extrapolation; (ii) To unify space and time and treat them as a single entity; (iii) To consider mesh values of dependent variables and their derivatives as independent variables, to be solved for simultaneously; (iv) To use only local discrete variables rather than global variables like the expansion coefficients used in spectral methods; (v) To define conservation elements and solution elements such that the simplest stencil will result; (vi) To require that, as much as possible, a numerical analogue be constructed so as to share with the corresponding physical equations the same space-time invariant properties, such that numerical dissipation can be minimized [5,10,24]; (vii) To exclude the use of characteristics-based techniques (such as Riemann solvers); and (viii) To avoid the use of ad hoc techniques as much as possible.

Moreover, the development of the CE/SE method is also guided by two basic beliefs that set it apart from the established methods. The first belief is that, in order to capture physics more efficiently and realistically, the modeling focus should be placed on the original integral form of the physical conservation laws, rather than the differential form. The latter form follows from the integral form under the additional assumption that the physical solution is smooth, an assumption that is difficult to realize numerically in a region of rapid change, such as a boundary layer or a shock. The second belief is that, with proper modeling of the integral and differential forms themselves, the resulting numerical solution should automatically be consistent with the properties derived from the integral and differential forms, e.g., the jump conditions across a shock and the properties of characteristics. In other words, a much simpler and more robust method can be developed by not using the above derived properties explicitly.

With the exception of the Navier-Stokes solver, all the 1D schemes described in [2] have been extended to become their 2D counterparts [9-11,14]. A more complete account of these new 2D schemes and their applications will be given in this and the following papers [3,4]. It will be shown in Sec. 3 that the spatial meshes used in these schemes are built from triangles-in such a manner that the resulting meshes are completely different from those used in the finite element method. As a result, these schemes are (i) compatible with the simplest unstructured meshes [31], and (ii) constructed without using the dimensional-splitting approach, i.e., without applying a 1D scheme in each coordinate direction. The dimensional-splitting approach is widely used in the construction of multidimensional upwind schemes. Unfortunately, this approach is flawed in several respects [33]. In particular, because a source term is not aligned with a special direction, it is difficult to imagine how this dimensional-splitting approach, in a logically consistent manner , can be used to solve a multidimensional problem involving source terms, such as those modeling chemical energy release.

Moreover, as will be shown shortly, because the CE/SE 2D schemes share with their 1D versions the same design principles, not only is the extension to 2D a straightforward matter, each of the new 2D schemes also shares with its 1D version virtually identical fundamental characteristics.

At this juncture, note that monotonicity conditions are not observed by general flow fields, e.g., those involving ZND detonation waves [21]. As a result, techniques involving monotonicity constraints are not used in the present development.

To give the reader, in advance, a concrete example that demonstrate the validity of the two basic beliefs referred to earlier, a self-contained Fortran program is listed in Appendix A. It is a CE/SE solver [23] for an extended Sod's shock tube problem that is the original Sod's problem [38] with the additional complication of imposing a non-reflecting boundary condition at each end of the computational domain. Note that the flow under consideration contains discontinuities and, relative to the computational frame, is subsonic throughout. It is well known that implementing a non-reflecting boundary condition for a subsonic flow is much more difficult than doing the same for a supersonic flow. This difficulty is further exacerbated by the fact that the traditional non-reflecting boundary conditions, e.g., the characteristic, the radiation (asymptotic), the buffer-zone, and the absorbing boundary conditions [39-44] are all based on an assumption that is not valid for the present case, i.e., that the flow is continuous. In spite of the fact that solving the present extended Sod's problem is substantially more difficult than the original Sod's problem, the main loop in the program listed herein contains only 39 Fortran statements. Not only is it very small in size, this program has a very simple logical structure. With the exception of a single ``if'' statement used to identify the time levels at which the non-reflecting boundary conditions must be imposed, it contains no conditional Fortran statements or functions such as ``if'', ``amax'', or ``amin'' that are often used in programs implementing high-resolution upwind methods. The small size of the listed program reflects the simplicity of the techniques employed by the CE/SE method to capture shock waves. It also results from the fact that the non-reflecting boundary conditions used in the present solver are the simple extrapolation conditions Eqs. (2.66) and (2.67) given in Sec. 2. They are much simpler than the traditional non-reflecting boundary conditions. On the other hand, the absence of Fortran conditional statements is a result of avoiding the use of ad hoc techniques. In spite of its small size and simple logical structure, according to the numerical results generated by the listed program (presented here as Figs. 1(a)-(c), with the numerical results and the exact solutions denoted by triangles and solid lines, respectively; see also [23]), the present solver is capable of generating nearly perfect non-reflecting solutions using the same time-step size from t = 0. Note that, at t = 10, all the waves have exited the computational domain, i.e., the exact solution is constant within it. The theoretical values of density, velocity, and pressure are approximately 0.4262000, 0.9277462 and 0.3030000, respectively. The maximum magnitudes of the errors in the numerically computed values of density, velocity, and pressure at t = 10 are approximately 0.0004, 0.0007, and 0.0004, respectively.

Note that Eqs. (2.66) and (2.67) represent only one of many sets of simple and robust non-reflecting boundary conditions developed especially for the CE/SE method [23]. Behind this development is a radical new concept based entirely on an assumption about the space-time flux distribution in the neighborhood of a spatial boundary. As a result, implementation of these CE/SE non-reflecting boundary conditions does not require the use of characteristics-based techniques.

To further demonstrate the nontraditional nature of the CE/SE method, the numerical results generated using the steady-state non-reflecting boundary conditions that were introduced and rigorously justified in [23] will also be presented here. Consider an alternate CE/SE solver that differs from the above CE/SE solver only in the fact that the steady-state boundary conditions Eq. (2.68) given in Sec. 2 are now taking the place of Eqs. (2.66) and (2.67). At t = 0.2, the waves generated in the interior of the computational domain have not yet reached the boundaries. In this case, with the given initial conditions (i.e., two different uniform states separated by a discontinuity located at the dead center of the domain), each of the above two solvers yield the same uniform solution in the vincinity of the right or left boundary. As a result, at t = 0.2, the numerical results generated by the alternate solver are identical to those shown in Fig. 1(a). The numerical results of the alternate solver at t = 0.4 are shown in Fig. 2(a). It is seen that, by this time, the shock wave has passed cleanly through the right boundary. There is good agreement between the numerical solution and the exact solution everywhere in the interior except for a slight disagreement in the vicinity of the right boundary. Note that the right boundary values, which do not vary with time , are discontinuous with respect to the neighboring interior values. The numerical results at t = 0.6 are shown in Fig. 2(b). As seen from the density profile, by this time, the contact discontinuity has also passed through the right boundary. Agreement between the numerical solution and the exact solution continue to be good in the interior. However, both left and right boundary values are now discontinuous with respect to the neighboring interior values.

Note that several recent applications [13,16,17,26,28] of the CE/SE method to 2D aeroacoustics problems reveal that: (i) the trivial nature of implementing CE/SE non-reflecting boundary conditions is manifested even for 2D problems; (ii) accuracy of the numerical results for nonlinear Euler problems is comparable to that of a 4-6th order compact difference scheme, even though nominally the CE/SE solver used is only of 2nd-order accuracy; and (iii) most importantly, the CE/SE method is capable of accurately modeling both small disturbances and strong shocks, and thus provides a unique tool for solving flow problems where the interactions between sound waves and shocks are important, such as the noise field around a supersonic over- and under-expanded jet. The fact listed in item (i) is more fundamental in nature, and will be further discussed in a separate paper. The following comments pertain to items (ii) and (iii):

(a)

Assuming the same order of accuracy, generally speaking, the accuracy of a scheme that enforces the space-time flux-conservation property is higher than that of a scheme that does not. A compact scheme generally does not enforce the flux-conservation property of the nonlinear Euler equations. On the contrary, not only is the present scheme flux-conserving, its accuracy in nonlinear calculations is enhanced by its surprisingly small dispersive errors [2,8,13,16,17]. Moreover, the nominal order of accuracy of an Euler solver is determined assuming a linearized form of the Euler equations. Thus its significance with respect to a highly nonlinear solution of the Euler equations may be questionable.

(b)

while numerical dissipation is required for shock resolution, it may also result in annihilation of small disturbances such as sound waves. Thus, a solver that can handle both small disturbances and strong shocks must be able to overcome this difficulty. It will be explained shortly that the CE/SE method is intrinsically endowed with this capability. On the other hand, a high-resolution upwind scheme that focuses only on shock resolution may introduce too much numerical dissipation [45].

Next we shall review briefly the inviscid version of the a-m scheme described in [2]. In addition to providing a historical perspective, the review will remove, once and for all, any lingering doubt from the reader's mind that the CE/SE method indeed differs substantially in both concept and methodology from the well-established methods. In particular, it will give in advance answers to questions such as: (i) is there any difference between the space-time elements used here and those used in the finite element method? and (ii) what are the key differences between the CE/SE method and other finite volume methods?

To proceed, consider an initial-value problem involving the PDE

At this juncture, the reader may wonder what the merit is of constructing a neutrally stable scheme. After all, it is well known that its nonlinear extensions generally are unstable. To address this question, the significance of constructing such a scheme and the critical role it plays in the development of the CE/SE method will be discussed immediately.

To proceed, note that there are several explicit and implicit extensions [2,12,25] of the a scheme which are solvers for

The above property is important because of the following observation: with a few exceptions, the numerical solution of a time-marching problem generally is contaminated by numerical dissipation. For a nearly inviscid problem, e.g., flow at a large Reynolds number, numerical dissipation may overwhelm physical dissipation and cause a complete distortion of the solution. To avoid such a difficulty, ideally a CE/SE solver for Eq. (1.2) or for the Navier-Stokes equations should possess the above special property. Obviously the development of such a solver must be preceded by that of a neutrally stable solver of Eq. (1.1).

The problem of physical dissipation being overwhelmed by numerical dissipation does not exist for a pure convection problem. However, as explained in the earlier discussion about the delicate nature of simulating small disturbances in the presence of shocks, numerical dissipation must still be handled carefully in this case.

Note that numerical dissipation traditionally is adjusted by varying the magnitude of added artificial dissipation terms. However, after being stripped of these added artificial dissipation terms, almost every traditional scheme (such as the Lax-Wendroff scheme) is still not free from inherent numerical dissipation. Hence, numerical dissipation generally cannot be avoided completely using the traditional approach.

This completes the discussion about the roles of non-dissipative schemes in the current development. To proceed further, the construction of the 1D a scheme will be described briefly.

Let x₁ = x, and x₂ = t be considered as the coordinates of a two-dimensional Euclidean space E₂. By using Gauss' divergence theorem in the space-time E₂, it can be shown that Eq. (1.1) is the differential form of the integral conservation law

Let W denote the set of all mesh points (j,n) in E₂ (dots in Fig. 4(a)) with n being a half or whole integer, and (j-n) being a half integer. For each (j,n) Î W, let the solution element SE(j,n) be the interior of the space-time region bounded by a dashed curve depicted in Fig. 4(b). It includes a horizontal line segment, a vertical line segment, and their immediate neighborhood. For the discussions given in this paper, the exact size of this neighborhood does not matter. However, in case the conservation law Eq. (1.3) takes a more complicated form in which the right side is a volume integral involving a source term, the SEs must fill the entire computational domain such that the volume integral can be modeled properly [21,22]. A SE that fulfills this requirement is depicted in Fig. 4(c).

For any (x,t) Î SE(j,n), let u(x,t) and [h\vec](x,t), respectively, be approximated by u^* (x,t ;j,n) and [h\vec]^* (x,t ;j,n) which we shall define shortly. Let

We shall require that u = u^* (x,t ;j,n) satisfy Eq. (1.1) within SE(j,n). As a result,

Let E₂ be divided into non-overlapping rectangular regions (see Fig. 4(a)) referred to as conservation elements (CEs). As depicted in Figs. 4(d) and 4(e), the CE with its top-right (top-left) vertex being the mesh point (j,n) Î W is denoted by CE_-(j,n) (CE₊(j,n)). The discrete approximation of Eq. (1.3) is then

According to Fig. 4(d), S(CE_-(j,n)), i.e., the boundary of CE_-(j,n), is formed by four line segments. Among them, [`AB] and [`AD] lie within SE(j,n). As a result, the flux leaving CE_-(j,n) through these two line segments will be evaluated using Eqs. (1.6) and (1.7) with the assumption that any point (x,t) on them belongs to SE(j,n). On the other hand, because [`CB] and [`CD] lie within SE(j-1/2,n-1/2), the flux leaving CE_-(j,n) through them will be evaluated assuming any point (x,t) on them belongs to SE(j-1/2,n-1/2).

According to Eq. (1.8), the total flux of [h\vec]^* leaving the boundary of any conservation element is zero. Because the surface integration over any interface separating two neighboring CEs is evaluated using the information from a single SE, obviously the local conservation relation Eq. (1.8) leads to a global flux conservation relation, i.e., the total flux of [h\vec]^* leaving the boundary of any space-time region that is the union of any combination of CEs will also vanish.

From the above discussions, it becomes obvious that the space-time element used in the finite element method differs from the current space-time SE and CE in both concept and the roles they serve. In particular, the former is not introduced to enforce flux conservation. In contrast to this, in the CE/SE method, flux conservation transmits information between neighboring SEs, and no global smoothness requirements are made on the solution to link neighboring SEs. This strategy enables the accurate capturing of traveling multidimensional solution discontinuities, e.g., moving multidimensional shock waves.

Furthermore, the CE/SE method is also fundamentally different from the traditional finite-volume methods such as the high-resolution upwind methods [46,47] and the discontinuous Galerkin method [48] in one important respect, i.e., because of the space-time staggering nature of its solution elements, the present method has a much simpler and consistent procedure to evaluate the flux at an interface. The key features of CE/SE flux-evaluation that distinguish it from those of the traditional methods are discussed in the following remarks:

(a)

Because an interface separating two neighboring CEs lies within a SE, the flux at this interface is evaluated without interpolation or extrapolation. Furthermore, the SE to which a particular interface belongs is determined by a rule that is independent of the local numerical solution. In other words, the concept of special upwind treatments and the complications that arise from these treatments are entirely foreign to the CE/SE method. To be more specific, consider the flux at the interface [`AD] depicted in Fig. 4(d). It is completely determined by u_jⁿ and (u_x)_jⁿ, two numerical variables associated with the predetermined mesh point (j,n), i.e., point A.

(b)

Flux evaluation is straightforward and it requires only simple integration involving the first-order Taylor's expansion. No complicated techniques such as the characteristics-based techniques are ever needed.

Finally, we also want to emphasize that the concepts used in the construction of the a scheme are fundamentally different from several schemes introduced by Nessyahu and Tadmor[49], and Sanders and Weiser [50] except that the meshes used by the a scheme and the latter schemes are all staggered in time. The key features of the a scheme that distinguish it from the latter schemes include: (i) the mesh values of both the dependent variable and its spatial derivative are considered as the independent variables, to be solved for simultaneously; and (ii) no interpolation or extrapolation techniques are used in the construction of the a scheme. Note that the differences between the latter schemes and an extension of the a scheme were also clearly spelled out by Huynh [51].

This section is concluded with the following remarks:

(a)

The a scheme can be constructed from a different perspective in which both CEs and SEs have the shape of a rhombus [2]. In this alternative construction, the differential condition Eq. (1.5) is not assumed. Instead it becomes a result of a local flux conservation condition and Eq. (1.4). In other words, the a scheme can be constructed entirely from flux conservation conditions and the assumption that u^* (x,t ;j,n) is linear in x and t.

(b)

The a scheme has many non-traditional features. They were discussed in great detail in [2].

(c)

Because there are two independent marching variables at each mesh point Î W, two amplification factors appear in the von Neumann stability analysis of the a scheme [2]. It happens that these two factors are identical to those of the Leapfrog scheme [52] if the latter factors arise from a ``correct'' von Neumann analysis [2]. Note that the main Leapfrog scheme (excluding its starting scheme which relates the mesh variables at the first two time levels), the Lax scheme [52], and the main DuFort-Frankel scheme [52] share one special property, i.e., a solution to any one of these schemes is formed by two decoupled solutions. Traditionally the von Neumann analysis for these schemes is performed without taking into account this decoupled nature. It is explained in [2] why such an erroneous analysis will result in a dispersive property prediction that makes the dispersion appear worse than it really is. Moreover, because (i) the a scheme and the Leapfrog scheme share the same amplification factors, and (ii) the a scheme is a two-level scheme while the Leapfrog scheme is a three-level scheme, the a scheme can be considered as a more advanced and compact Leapfrog scheme.

   

     The fact that the amplification factors of the a scheme are related to those of a celebrated classical scheme is only one among a string of similar unexpected coincidences encountered during the development of the CE/SE method. As it turns out [2,12,25], the amplification factors of the Lax, the Crank-Nicolson, and the DuFort-Frankel schemes also are related to those of some of the extensions of the a scheme.

In this section, we shall (i) review and reformulate the 1D schemes described in [2], and (ii) fill a gap in the derivation of Eq. (4.28) in [2]. Not only does the reformulation enable the reader to see more clearly the structural similarity between the 1D solvers of Eq. (1.1) and their Euler counterparts, it also makes it easier for him to appreciate the consistency between the construction of the 1D CE/SE solvers and that of the 2D solvers to be described in the later sections.

2.1. The a Scheme

As the first step, the marching procedure of the a scheme will be cast into a form slightly different from that given in [2]. To proceed, let the Courant number n = aD t/D x. Also let

Because both (s₊)_j+1/2^n-1/2 and (s_-)_j-1/2^n-1/2 are explicit functions of the marching variables at the (n-1/2)th time level, Eqs. (2.7) and (2.8) form the explicit marching procedure for the a scheme. Note that these equations can be obtained from the inviscid form of the a-m scheme, i.e., Eq. (2.14) in [2]. Also note that the superscript symbol ``a'' in the parameter (u^a+_x)_jⁿ is introduced to remind the reader that Eq. (2.8) is valid for the a scheme.

2.2. The a-e Scheme

In the a-e scheme [2], CE₊(j,n) and CE_-(j,n), which are depicted in Figs. 4(d) and 4(e), respectively, are not considered as conservation elements, i.e., Eq. (1.8) is no longer applicable. Instead, one assumes that

It was explained in [2] that Eq. (2.7) follows directly from Eq. (2.9). As a result, the former is also valid in the a-e scheme. The a-e scheme is formed by Eq. (2.7) and another equation that differs from Eq. (2.8) only in the expression on the right side. To show more clearly the similarity of the 1D schemes and their 2D versions to be described in the later sections, in the following, the counterpart of Eq. (2.8) in the a-e scheme will be rederived from a perspective different from that presented in [2].

Let (j,n) Î W. Then (j±1/2,n-1/2) Î W. Let

According to Eq. (2.10), u_j±1/2^¢ n can be interpreted as a first-order Taylor's approximation of u at (j±1/2,n). Thus

is a central-difference approximation of ¶u / ¶x at (j,n), normalized by the same factor D x/4 that appears in Eq. (2.1). Note that the superscript ``c'' is used to remind the reader of the central-difference nature of the term (u^c+_x)_jⁿ. In the a-e scheme, Eq. (2.8) is replaced by

At this juncture, note that, at each mesh point (j,n) Î W, Eqs. (2.7) and (2.8) are the results of two conservation conditions given in Eq. (1.8). Because Eq. (2.13) does not reduce to Eq. (2.8) except in the special case e = 0, at each mesh point (j,n) Î W, generally the a-e scheme satisfies only the single conservation condition Eq. (2.9) rather than the two consevation conditions Eq. (1.8). However, because (u^a+_x)_jⁿ generally is present on the right side of Eq. (2.13), the a-e scheme generally will still be burdened with the cost of solving two conservation conditions at each mesh point. The exception occurs only for the special case e = 1/2, under which Eq. (2.13) reduces to (u⁺_x)_jⁿ = (u^c+_x)_jⁿ.

Note that the first part of the expression on the right side of Eq. (2.13), i.e., (u^a+_x)_jⁿ, emerges from the development of the non-dissipative a scheme. As a result, it is the non-dissipative part. On the other hand, the second part, whose magnitude can be adjusted by the parameter e, represents numerical dissipation introduced by the difference between the central difference term (u^c+_x)_jⁿ and the non-dissipative term (u^a+_x)_jⁿ. Thus one may heuristically conclude that the numerical dissipation associated with the a-e scheme can be increased by increasing the value of e. It was shown in [2] that this conclusion is indeed valid in the stability domain of the a-e scheme, i.e.,

According to Eqs. (2.4)-(2.6), (2.11) and (2.12), both (u^c+_x)_jⁿ and (u^a+_x)_jⁿ are explicitly dependent on n (and therefore explicitly dependent on D t). However, (u^c+_x-u^a+_x)_jⁿ is not dependent on n. As a matter of fact, it can be shown that

Note that, in the original development [2], (du_x)_jⁿ was introduced to break the symmetry of the stencil of the a scheme with respect to space-time inversion. This symmetry breaking results in the a-e scheme that was originally defined by the matrix equation Eq. (3.6) of [2]. Its two component equations are Eq. (2.7) and

2.3. The Euler a Scheme

For a reason that will soon become obvious to the reader, reformulation of the inviscid (m = 0) version of the Navier-Stokes solver described in Section 5 of [2] will precede that of the Euler solvers described in Section 4 of [2]. Because the inviscid version is also an Euler solver and, like the a scheme, is free of numerical dissipation if it is stable, it will be referred to as the Euler a scheme.

To proceed, consider the Euler equations [2]

To proceed, let (i)

For any (x,t) Î SE(j,n), u_m(x,t), f_m(x,t) and [h\vec]_m(x,t) are approximated by

For all (j,n) Î W, we assume that

For each (j,n) Î W, let (i)

Note that the flux conservation conditions Eqs. (2.2) and (2.3), and its Euler counterparts, i.e., Eqs. (2.28) and (2.29) share the same algebraic structure. As a matter of fact, the former pair will become the latter pair if the symbols 1, n, u and u_x⁺ are replaced by I, F⁺, [u\vec] and [u\vec]_x⁺, respectively. As a result, Eqs. (2.28) and (2.29) will be solved by a procedure similar to that used earlier to extract Eqs. (2.7) and (2.8) from Eqs. (2.2) and (2.3). However, because (i) matrix multiplication is not commutative and (ii) the matrix (F⁺)_jⁿ is a function of (u_m)_jⁿ, m = 1,2,3, while n is a simple constant, as will be shown shortly, the algebraic structure of the solution to Eqs. (2.28) and (2.29) is more complex than that of Eqs. (2.7) and (2.8).

To proceed, let (j,n) Î W and

Equation (2.32) represents the first part of the solution to Eqs. (2.28) and (2.29). To obtain the second part, one must assume the existence of the inverse of the matrix [I-(F⁺)²]_jⁿ for all (j,n) Î W. In the following, we shall briefly discuss the significance of this assumption.

Let v and c be the fluid speed and sonic speed, respectively. They are known functions of u_m, m = 1,2,3 [2]. For each (j,n) Î W, let v_jⁿ and c_jⁿ, respectively, denote the values of v and c when u_m, m = 1,2,3, respectively, assume the values of (u_m)_jⁿ, m = 1,2,3. Let

Let (j,n) Î W. Let the marching variables at the (n-1/2)th time level be given. Then [u\vec]_jⁿ can be evaluated using Eq. (2.32). Because [I±F⁺]_jⁿ is a function of [u\vec]_jⁿ, it follows that

To obtain the second part of the solution to Eqs. (2.28) and (2.29), they are multiplied from the left by

It has been shown by numerical experiments that the Euler a scheme is neutrally stable in the interior of the computational domain up to at least a thousand time steps when n_jⁿ < 1 for all (j,n) Î W. In these numerical experiments involving a shock-tube problem, the computational domain was allowed to grow with time, so that the undisturbed fluid state could always be prescribed at the computational boundaries as the exact solution. As a matter of fact, by using an analysis similar to that given at the end of Sec. 6 in [7], one can show that the linearized form of the Euler a scheme is neutrally stable when n_jⁿ < 1 for all (j,n) Î W.

The parameters ([S\vec]₊)_jⁿ and ([S\vec]_-)_jⁿ can be evaluated by using Eqs. (2.37) and (2.38) directly. This direct evaluation involves inverting two 3×3 matrices which is computationally costly. In the following, we shall describe a more efficient approach.

According to Eqs. (2.37) and (2.38), ([S\vec]₊)_jⁿ and ([S\vec]_-)_jⁿ are the solutions to

2.4. The Simplified Euler a scheme

In implementing the Euler a scheme, two systems of linear equations must be solved for each (j,n) Î W. As a result, the Euler a scheme is locally implicit [1, p.22]. In this subsection we shall develop a simplified version that is completely explicit.

To proceed, the expressions

Generally CE_±(j,n), (j,n) Î W, are not conservation elements in the simplified scheme. However, because Eq. (2.32) is equivalent to the conservation condition [2]

Note that by replacing the symbols s₊, s_-, u^a+_x, u, u_x⁺, 1 and n in Eqs. (2.4)-(2.8) by [s\vec]₊, [s\vec]_-, [u\vec]^a¢+_x, [u\vec], [u\vec]_x⁺, I and F⁺, respectively, these equations will become Eqs. (2.30), (2.31), (2.44), (2.32) and (2.45), respectively. In other words, the a scheme and the simplified Euler a scheme share the same algebraic structure.

The simplified Euler a scheme generally is unstable. However, as will be shown shortly, this scheme can be extended to become the simplified Euler a-e scheme which does have a large stability domain.

2.5. The Euler a-e Scheme

The process by which the a-e scheme was constructed from the a scheme will be used to construct the Euler a-e scheme from the Euler a scheme.

In the Euler a-e scheme, the conservation conditions given in Eq. (2.46) are assumed. Because Eq. (2.32) is equivalent to Eq. (2.46), the former is also a part of of the Euler a-e scheme. The Euler a-e scheme is formed by Eq. (2.32) and another equation that differs from Eq. (2.40) only in the expression on the right side.

To proceed, let (j,n) Î W and

2.6. The Simplified Euler a-e Scheme

According to Eq. (2.50), excluding the special case e = 1/2 , implementation of the Euler a-e scheme also requires the evaluation of ([u\vec]^a+_x)_jⁿ and therefore (see Eqs. (2.37)-(2.39)) the solution of Eqs. (2.41) and (2.42). Thus the Euler a-e scheme is locally implicit if e ¹ 1/2. A totally explicit variant, referred to as the simplified Euler a-e scheme, is defined by Eq. (2.32) (or, equivalently, Eq. (2.46)) and

Note that by replacing the symbols s₊, s_-, u^a+_x, u, u_x⁺, u¢, u^c+_x, 1 and n in Eqs. (2.4)-(2.7) and (2.11)-(2.13) by [s\vec]₊, [s\vec]_-, [u\vec]^a¢+_x, [u\vec], [u\vec]_x⁺, [u\vec]¢, [u\vec]^c+_x, I and F⁺, respectively, these equations will become Eqs. (2.30), (2.31), (2.44), (2.32), (2.48), (2.49) and (2.52) respectively. In other words, the a-e scheme and the simplified Euler a-e scheme share the same algebraic structure.

It has been shown numerically that the simplified Euler a-e scheme is stable if

According to Eqs. (2.30), (2.31), (2.44), (2.48) and (2.49), both ([u\vec]^c+_x)_jⁿ and ([u\vec]^a¢+_x)_jⁿ are explicitly dependent on the the matrices (F⁺)_j+1/2^n-1/2 and (F⁺)_j-1/2^n-1/2 (and therefore explicitly dependent on D t). However, ([u\vec]^c+_x-[u\vec]^a¢+_x)_jⁿ is free from this dependency. Let (i) (du_mx)_jⁿ be the parameter defined by Eq. (4.26) in [2], and (ii) (d[u\vec]_x)_jⁿ be the column matrix formed by (du_mx)_jⁿ, m = 1,2,3. Then it can be shown that

With the above preliminaries, we are now ready to provide a proof for Eq. (4.28) in [2]. Note that the last equation was introduced in [2] simply as a ``natural generalization'' of Eq. (3.10) in [2].

To proceed, note that Eq. (2.47) is the matrix form of Eq. (4.27) in [2], i.e., [u\vec]_j±1/2 ^¢ n is the column matrix formed by (u_m¢)_j±1/2ⁿ, m = 1,2,3, which were introduced in the latter equation. As a result, with the aid of Eqs. (2.27), (2.49) and (2.54), Eq. (2.52) can be rewritten as

Because Eqs. (4.24) in [2] are equivalent to Eq. (2.32), the Euler scheme defined by Eqs. (4.24) and (4.28) in [2] is identical to the simplified Euler a-e scheme.

2.7. The a-e-a-b Scheme and Its Euler Versions

Consider the a-e scheme defined by Eqs. (2.7) and (2.13). If discontinuities are present in a numerical solution, the above scheme is not equipped to suppress numerical wiggles that generally appear near these discontinuities. In the following, we shall describe a remedy for this deficiency.

Let

The expression on the right side of Eq. (2.60) contains three parts. The first part is a non-dissipative term (u^a+_x)_jⁿ. The second part is the product of 2e and the difference between the central difference term (u^c+_x)_jⁿ and the non-dissipative term (u^a+_x)_jⁿ. The third part is the product of b and the difference between a weighted average of (u^c+_x+)_jⁿ and (u^c+_x-)_jⁿ and their simple average. Numerical dissipation introduced by the second part generally is effective in damping out numerical instabilities that arise from the smooth region of a solution. But it is less effective in suppressing numerical wiggles that often occur near a discontinuity. On the other hand, numerical dissipation introduced by the third part is very effective in suppressing numerical wiggles. Moreover, because the condition |(u^c+_x+)_jⁿ| = |(u^c+_x-)_jⁿ| more or less prevails and thus the weighted average is nearly equal to the simple average (see the comment given immediately following Eq. (2.58)) in the smooth region of the the solution, numerical dissipation introduced by the third part has very slight effect in the smooth region.

From the above disscusion, one concludes that there are two different types of numerical dissipation associated with the a-e-a-b scheme and they complement each other. As a result, the a-e-a-b scheme can handle both small disturbances and sharp discontinuies simultaneously if the values of e, a and b are chosen properly (usually e = 1/2, a = 1,2 and b = 1). Also note that, to give the CE/SE method more flexibility in controlling local numerical dissipation, the parameters e and b can even be considered as functions of local dynamical variables and mesh parameters (see Sec. 8).

Similarly, the Euler a-e scheme and the simplified Euler a-e scheme can be modified to become the Euler a-e-a-b scheme and the simplified Euler a-e-a-b scheme, respectively, by simply replacing Eqs. (2.50) and (2.52) with

2.8. The 1D CE/SE Shock-Capturing Scheme

Let e = 1/2 and b = 1. Then the Euler a-e-a-b scheme and the simplified Euler a-e-a-b scheme reduce to the same scheme. The reduced scheme is defined by Eq. (2.32) and

The above scheme is one of the simplest among the Euler solvers known to the authors. The value of a is the only adustable parameter allowed in this scheme. Because it is totally explicit and has the simplest stencil, the scheme is also highly compatible with parallel computing. Furthermore, it has been shown that the scheme can accurately capture shocks and contact discontinuities with high resolution and no numerical oscillations. For these distinctive features and for convenience of future reference, the reduced scheme will be given a special name, i.e., the 1D CE/SE shock-capturing scheme. Note that this scheme with a = 1 is implemented in the two shock-tube solvers referred to in Sec. 1. Consider only the case that all spatial boundary points (j,n) Î W are at the time levels n = 0,1,2,¼ (see Fig. 4(a)). The non-reflecting boundary conditions used in the first solver, i.e., the one listed in Appendix A, are: (i)

In Sec. 2, it was established that, for each 1D CE/SE solver, there were 2M independent marching variables per mesh point with M being the number of conservation laws to be solved. Because M conservation conditions are imposed over each CE, two CEs were introduced at each mesh point such that both the 1D a scheme and the 1D Euler a scheme can be constructed by solving, at each mesh point (j,n) Î W, for the 2M variables using the 2M conservation conditions imposed over CE_-(j,n) and CE₊(j,n).

As will be shown in the following sections, for each 2D CE/SE solver, there are 3M independent marching variables per mesh point. As a result, construction of the 2D a scheme and the 2D Euler a scheme demands that three CEs be defined at each mesh point. In this section, only the basic geometric structures of these CEs will be described.

Consider a spatial domain formed by congruent triangles (see Fig. 5). The center of each triangle is marked by either a hollow circle or a solid circle. The distribution of these hollow and solid circles is such that if the center of a triangle is marked by a solid (hollow) circle, then the centers of the three neighboring triangles with which the first triangle shares a side are marked by hollow (solid) circles. As an example, point G , the center of the triangle DBDF, is marked by a solid circle while points A, C and E, the centers of the triangles DFMB, DBJD and DDLF, respectively, are marked by hollow circles. These centers are the spatial projections of the space-time mesh points used in the 2D CE/SE solvers.

To specify the exact locations of the mesh points in space-time, one must also specify their temporal coordinates. In the 2D CE/SE development, again we assume that the mesh points are located at the time levels n = 0,±1/2,±1,±3/2,¼ with t = n D t at the n th time level. Furthermore, we assume that the spatial projections of the mesh points at a whole-integer (half-integer) time level are the points marked by hollow (solid) circles in Fig. 5.

Let the triangles depicted in Fig. 5 lie on the time level n = 0. Then those points marked by hollow circles are the mesh points at this time level. On the other hand, those points marked by solid circles are not the mesh points at the time level n = 0. They are the spatial projections of the mesh points at half-integer time levels.

Points A, C and E, which are depicted in Figs. 5 and 6(a), are three mesh points at the time level n = 0. Point G¢, which is depicted in Fig. 6(a), is a mesh point at the time level n = 1/2. Its spatial projection at the time level n = 0 is point G. Because point G is not a mesh point, it is not marked by a circle in the space-time plots given in Figs. 6(a) and 6(c). Hereafter, only a mesh point, e.g., point G¢, will be marked by a solid or hollow circle in a space-time plot.

The conservation elements associated with point G¢ are defined to be the space-time quadrilateral cylinders GFABG¢F¢A¢B¢, GBCDG¢B¢C¢D¢, and GDEFG¢D¢E¢F¢ that are depicted in Fig. 6(a). Here (i) points B, D and F are the vertices of the triangle with point G as its center (centroid) (see also Fig. 5), and (ii) points A¢, B¢, C¢, D¢, E¢ and F¢ are on the time level n = 1/2 with their spatial projections on the time level n = 0 being points A, B, C, D, E and F, respectively.

Point G¢ is a mesh point at a half-integer time level. For a mesh point at a whole-integer time-level, the conservation elements associated with it can be constructed in a similar fashion. As an example, consider Fig. 6(b). Here points B¢( B¢¢), I¢( I¢¢), J¢( J¢¢), K¢( K¢¢), D¢( D¢¢), G¢( G¢¢) and C¢( C¢¢) are on the time level n = 1/2 (n = 1) with their spatial projections on the time level n = 0, respectively, being the points B, I, J, K, D, G and C that are depicted in Fig. 5. Point C¢¢ is a mesh point at the time level n = 1. By definition, the conservation elements associated with point C¢¢ are the quadrilateral cylinders C¢J¢K¢D¢C¢¢J¢¢K¢¢D¢¢, C¢D¢G¢B¢C¢¢D¢¢G¢¢B¢¢ and C¢B¢I¢J¢C¢¢B¢¢I¢¢J¢¢. The relative space-time positions of the six CEs associated with mesh points G¢ and C¢¢ are depicted in Fig. 6(c).

Recall that, in the development of the 1D a scheme, a pair of diagonally opposite vertices of each CE_±(j,n) (see Figs. 4(d) and 4(e)) are assigned as mesh points. Furthermore, the boundary of each CE_±(j,n) is a subset of the union of the SEs associated with the two diagonally opposite mesh points of this CE. In the 2D development, as seen from Figs. 6(a)-(c), two diagonally opposite vertices of each CE are also assigned as mesh points. In Sec. 4, we shall define the SEs such that even in the 2D case, the boundary of a CE is again a subset of the union of the SEs associated with the two diagonally opposite mesh points of this CE.

As a preliminary to the derivation of several equations to be given in Sec. 4, this section is concluded with a discussion of several geometric relations involving point G and the vertices of the hexagon ABCDEF that are depicted in Fig. 5. By using the facts that (i) points A, C, E and G are the geometric centers of four neighboring congruent triangles \triangle FMB, \triangle BJD, \triangle DLF and \triangle BDF, respectively; and (ii) any two of the above four triangles form a parallelogram (note: two congruent triangles sharing one side may not form a parallelogram), it can be shown that:

(a)

[`CD], [`GE], [`BG] and [`AF] are parallel line segments of equal length.

(b)

[`AB], [`GC], [`FG] and [`ED] are parallel line segments of equal length.

(c)

[`BC], [`GD], [`AG] and [`FE] are parallel line segments of equal length.

(d)

Point G is the geometric center of the hexagon ABCDEF and the triangle ACE.

Note that the line segments [`GA], [`GC], [`GE] [`AC], [`CE] and [`EA] are not shown in Fig. 5. Also note that, because the hexagon BIJKDG (depicted in Fig. 5) is congruent to the hexagon ABCDEF, a set of geometric relations similar to those listed above also exists for the vertices and the center of the hexagon BIJKDG.

In this section, we consider a dimensionless form of the 2D convection equation, i.e.,

In the following analysis, the nontraditional space-time mesh that was sketched in Sec. 3 will be rigorously defined. To proceed, the spatial projections of the mesh points depicted in Fig. 5 are reproduced in Fig. 7. Note that the dashed lines that appear in Fig. 7 are the spatial projections of the vertical interfaces (see Fig. 6(a)-(c)) that separate different CEs. Also note that, as a result of the geometric relations listed at the end of Sec. 3, any dashed line can point only in one of three different fixed directions. We assume that the congruent triangles depicted in Fig. 5 are aligned such that one of the above fixed directions is the x-direction.

Each mesh point marked by a solid or hollow circle is assigned a pair of spatial indices (j,k) according to the location of its spatial projection. Obviously, a mesh point can be uniquely identified by its spatial indices (j,k) and the time level n where it resides. According to Figs. 8 and 9, the spatial projections of the mesh points that share the same value of j (k) lie on a straight line on the x-y plane with this straight line pointing in the direction of the k- (j-) mesh axis.

Let

Each mesh point (j,k,n) Î W is associated with (i) three conservation elements (CEs), denoted by CE_r(j,k,n), r = 1,2,3 (see Figs. 10(a) and 11(a)); and (ii) a solution element (SE), denoted by SE(j,k,n) (see Figs. 10(b) and 11(b)). Each CE is a quadrilateral cylinder in space-time while each SE is the union of three vertical planes, a horizontal plane, and their immediate neighborhoods. Note that the CEs and the SE associated with a mesh point (j,k,n) Î W₁ differ from those associated with a mesh point (j,k,n) Î W₂ in their space-time orientations.

By using the geometric relations listed at the end of Sec. 3, one can conclude that the spatial coordinates of the vertices of the hexagon ABCDEF, which appears in both Figs. 10(a) and 11(a), are uniquely determined by three positive parameters w, b and h (see Fig. 12(a)) if (i) one assumes that [`DA] is aligned with the x-direction, and (ii) the spatial coordinates of point G (the centroid of the hexagon) are given. Note that w, b and h, respectively, are the lengths of the line segments [`DM], [`MH] and [`BH] with (i) [`DM] being a median of the triangle \triangle BDF, and (ii) points G, M and H being on the line segment [`DA]. Also note that a dashed line in Fig. 7 may appear in other figures as a solid line.

According to Fig. 7, E₃ can be filled with the CEs defined above. Moreover, it is seen from Figs. 10(a), 10(b), 11(a), and 11(b) that the boundary of a CE is formed by the subsets of two neighboring SEs .

Let the space-time mesh be uniform, i.e., the parameters D t, w, b, and h are constants. Let x_j,k and y_j,k be the x- and y- coordinates of any mesh points (j,k,n) Î W. Let x_0,0 = 0 and y_0,0 = 0. Then information provided by Figs. 12(a) and 12(b) implies that

For any (x,y,t) Î SE(j,k,n), u(x,y,t) and [h\vec](x,y,t), respectively, are approximated by

where u_j,kⁿ, (u_x)_j,kⁿ, (u_y)_j,kⁿ, and (u_t)_j,kⁿ are constants within SE(j,k,n). The last four coefficients, respectively, can be considered as the numerical analogues of the values of u, ¶u / ¶x, ¶u / ¶y, and ¶u / ¶t at (x_j,k,y_j,k,tⁿ). As a result, the expression on the right side of Eq. (4.7) can be considered as the first-order Taylor's expansion of u(x,y,t) at (x_j,k,y_j,k,tⁿ). Also note that Eq. (4.8) is the numerical analogue of Eq. (4.3).

We shall require that u = u^*(x,y,t;j,k,n) satisfy Eq. (4.1) within SE(j,k,n). As a result,

As a result of Eq. (4.11), the total flux leaving the boundary of any CE is zero. Because the flux at any interface separating two neighboring CEs is calculated using the information from a single SE, the flux entering one of these CEs is equal to that leaving another. It follows that the local conservation conditions Eq. (4.11) will lead to a global conservation condition, i.e., the total flux leaving the boundary of any space-time region that is the union of any combination of CEs will also vanish .

In the following, several preliminaries will be given prior to the evaluation of Eq. (4.11). To proceed, note that a mesh line with j and n being constant or a mesh line with k and n being constant is not aligned with the x-axis or the y-axis. We shall introduce a new spatial coordinate system (z,h) with its axes aligned with the above mesh lines (see Fig. 12(c)).

Let [e\vec]_x and [e\vec]_y be the unit vectors in the x- and the y- directions, respectively. Let [e\vec]_z and [e\vec]_h be the unit vectors in the directions of DF and DB (i.e., the j- and the k- directions-see Figs. 12(a)-(c)), respectively. It can be shown that

With the aid of Eqs. (4.5), (4.18), and (4.20), it can be shown that the coordinates (z,h) of any mesh point (j,k,n) Î W are given by

Next we shall introduce several coefficients that are tied to the coordinate system (z,h). Let

(a)

(a_z,a_h) are the contravariant components with respect to the coordinates (z,h) for the spatial vector whose x- and y- components are a_x and a_y, respectively.

(b)

((u_z)ⁿ_j,k,(u_h)ⁿ_j,k) are the covariant components with respect to the coordinates (z,h) for the spatial vector whose x- and y- components are (u_x)_j,kⁿ and (u_y)_j,kⁿ, respectively.

(c)

Because the contraction of the contravariant components of a vector and the covariant components of another is a scalar, Eq. (4.9) can be rewritten as

(ut)j,kn = -[az(uz)nj,k+ah(uh)nj,k]     (4.24)

(d)

Under the linear coordinate transformation defined by Eqs. (4.17) and (4.18), (z-jDz,h-kDh) are the contravariant components with respect to the coordinates (z,h) for the spatial vector whose x- and y- components are x-x_j,k and y-y_j,k, respectively. Using the same reason given in (c), Eq. (4.10) implies that

u*(x,y,t;j,k,n) = u*(z,h,t;j,k,n)     (4.25)

where

u*(z,h,t;j,k,n) def =   uj,kn + (uz)nj,k[(z-jDz)-az(t-tn)] + (uh)nj,k[(h-kDh)-ah(t-tn)] (4.26)

Note that Eqs. (4.24) and (4.25) can also be verified directly using Eqs. (4.18), (4.20), (4.22), and (4.23).

Next, let (i)

Furthermore, let

(a)

Each of Eqs. (4.29)-(4.46) represents two equations. One corresponds to the upper signs while the other, to the lower signs.

(b)

The definitions given in Eqs. (4.29)-(4.37) will be used in the first marching step of the 2D a scheme; while those given in Eqs. (4.38)-(4.46) will be used in the second marching step. It is seen that the expressions on the right sides of the former can be converted to those of the latter, respectively, by reversing the ``+'' and ``-'' signs. Moreover, for every pair of r and s (r,s = 1,2,3), s^(1)-_rs and s^(2)-_rs are converted to s⁽²⁾⁺_rs and s⁽¹⁾⁺_rs, respectively, if n_z, and n_h are replaced by -n_z, and -n_h, respectively.

(c)

We have

s(q)±11+s(q)±21+s(q)±31 = 3,       q = 1,2     (4.47)

and

s(q)±12+s(q)±22+s(q)±32 = s(q)±13+s(q)±23+s(q)±33 = 0,       q = 1,2 (4.48)

To simplify the following development, let

Equation (4.11) is evaluated in Appendix B. Let (j,k,n) Î W_q with q = 1,2. Then, for any r = 1,2,3, the result of evaluation can be expressed as:

According to Eqs. (4.29)-(4.31), s^(1)±₁₁, s^(1)±₁₂, and s^(1)±₁₃ contain a common factor (1-n_z-n_h). Similarly, each of three consecutive pairs of coefficients defined in Eqs. (4.32)-(4.46) also contain a common factor. As a result, if one assumes that (i) 1-n_z-n_h ¹ 0, (ii) 1+n_z ¹ 0, (iii) 1+n_h ¹ 0, (iv) 1+n_z+n_h ¹ 0, (v) 1-n_z ¹ 0 and (vi) 1-n_h ¹ 0, i.e.,

respectively. Here

Note that the expressions within the brackets in Eqs. (4.53)-(4.55) and (4.59)-(4.61), respectively, can be converted to those in Eqs. (4.56)-(4.58) and (4.62)-(4.64) by reversing the ``+'' and ``-'' signs.

It can be shown that Eqs. (4.53)-(4.55) are equivalent to

At this juncture, it should be emphasized that Eqs. (4.65) and (4.68) can be derived directly from Eq. (4.51). As a matter of fact, with the aid of Eqs. (4.47) and (4.48), we can obtain Eq. (4.65) (Eq. (4.68)) by summing over the three equations with q = 1 (q = 2) and r = 1,2,3 in Eq. (4.51).

The 2D a scheme is formed by repeatedly applying the two marching steps defined by Eqs. (4.65)-(4.67) and Eqs. (4.68)-(4.70), respectively. It has been shown numerically that it is of second order in accuracy for u_j,kⁿ, (u_z)ⁿ_j,k and (u_h)ⁿ_j,k assuming that n_z and n_h are held constant in the process of mesh refinement (note: as a result of Eq. (4.28), the 2D a scheme is third order accurate for (u⁺_z)ⁿ_j,k and (u⁺_h)ⁿ_j,k). Note that the superscript symbol ``a'' in (u^a+_z)_j,kⁿ and (u^a+_h)_j,kⁿ is introduced to remind the reader that Eqs. (4.66), (4.67), (4.69) and (4.70) are valid for the 2D a scheme. Although the 2D a scheme is constructed using a procedure very much parallel to that used to construct the 1D a scheme, the former is more complex than the latter in many aspects. One key difference between these two schemes is that the 2D a scheme is formed by two distinctly different marching steps while the 1D a scheme is formed by repeatedly applying the same marching step defined by the inviscid version of Eq. (2.14) in [2]. It is this difference that, in the 2D case, makes it necessary to divide the mesh points into two sets W₁ and W₂.

As a preliminary for the stability analysis of the 2D a scheme given in Sec. 6, for any (j,k,n) Î W, let

Furthermore, let the six 3×3 matrices Q^(q)_r, q = 1,2, and r = 1,2,3, respectively, be the special cases of those defined in Eqs. (5.18)-(5.23) (see Sec. 5) with e = 0. Then Eqs. (4.65)-(4.70) can be expressed as

Combining Eqs. (4.72) and (4.49a)-(4.50c), one has (i)

The 2D a scheme has several nontraditional features. They are summarized in the following comments:

(a)

As in the case of the 1D a scheme, the 2D a scheme also has the simplest stencil possible, i.e., in each of their two marching steps, the stencil is a tetrahedron in 3D space-time with one vertex at the upper time level and the other three vertices at the lower time level.

(b)

As in the case of the 1D a scheme, each of the six flux conservation conditions associated with the 2D a scheme., i.e., those given in Eq. (4.51) (q = 1,2 and r = 1,2,3), represents a relation among the marching variables associated with only two neighboring SEs .

(c)

As in the case of the 1D a scheme, the 2D a scheme also is non-dissipative if it is stable. It is shown in Sec. 7 that the 2D a scheme is neutrally stable if

|nz| < 1.5,    |nh| < 1.5,   and    |nz+nh| < 1.5     (4.75)

As depicted in Fig. 14, the domain of stability defined by Eq. (4.75) is a hexagonal region in the n_z-n_h plane. Moreover, it will also be shown later that Eq. (4.75) can be interpreted as the requirement that the physical domain of dependence of Eq. (4.1) be within the numerical domain of dependence. Note that the points on the n_z-n_h plane that violate Eq. (4.52) form the boundary of a hexagonal region which is entirely within the stability domain defined in Eq. (4.75). As was emphasized earlier, the 2D a scheme applies even at these points.

(d)

It is shown in [9] that the 2D a scheme has the following property, i.e., for any (j,k,n) Î W,

® q    (j,k,n+1)® ® q    (j,k,n)    as     D t ® 0     (4.76)

if a_x, a_y, w, b, and h are held constant. The 1D a scheme also possesses a similar property, i.e., Eq. (2.19) in [2]. The above property usually is not shared by other schemes that use a mesh that is staggered in time, e.g., the Lax scheme [52].

(e)

As in the case of the 1D a scheme, the 2D a scheme is also a two-way marching scheme. In other words, Eqs. (4.53)-(4.58) can also be used to construct the backward time-marching version of the 2D a scheme. More discussions on this subject are given in [9].

This section is concluded with the following remarks:

(a)

the 2D a scheme is only a special case of the 2D a-m scheme described in [9]. It is a solver for the 2D convection-diffusion equation

¶u ¶t +ax  ¶u ¶x +ay  ¶u ¶y -m æ ç è ¶2 u ¶x2 + ¶2 u ¶y2 ö ÷ ø = 0     (4.77)

where a_x, a_y, and m ( ³ 0) are constants. Note that this solver, as in the case of its 1D counterpart, is unconditionally stable if a_x = a_y = 0.

(b)

It should be emphasized that, with the aid of Eqs. (4.17)-(4.20), (4.22), and (4.23), the 2D a scheme can also be expressed in terms of the marching variables and the coefficients tied to the coordinates (x,y). In other words, the coordinates (z,h) are introduced solely for the purpose of simplifying the current development. The essence of the 2D a scheme, and the schemes to be introduced in the following sections, is not dependent on the choice of the coordinates in terms of which these schemes are expressed.

The 2D a scheme is non-dissipative and reversible in time. It is well known that a non-dissipative numerical analogue of Eq. (4.1) generally becomes unstable or highly dispersive when it is extended to model the 2D unsteady Euler equations. It is also obvious that a scheme that is reversible in time cannot model a physical problem that is irreversible in time, e.g., an inviscid flow problem involving shocks. As a result, the 2D a scheme will be extended to become the dissipative 2D a-e and a-e-a-b scheme before it is extended to model the Euler equations. As will be shown, the 2D extensions are carried out in a fashion completely parallel to their 1D counterparts.

5.1. The 2D a-e Scheme

To proceed, note that the CEs for the 2D a-e scheme generally are not those associated with the 2D a scheme. Here only a single CE is associated with a mesh point (j,k,n) Î W. This CE, denoted by CE(j,k,n), is the union of CE_r(j,k,n), r = 1,2,3. In other words,

Moreover, because of Eq. (5.1), Eq. (5.2) must be true if Eq. (4.11) is assumed. As a matter of fact, a direct evaluation of Eq. (5.2) reveals that it is equivalent to Eq. (4.65) (Eq. (4.68)) if (j,k,n) Î W₁ ((j,k,n) Î W₂). As a result, Eqs. (4.65) and (4.68) are shared by the 2D a scheme and 2D a-e scheme. Recall that Eq. (2.7) is also shared by the 1D a and a-e schemes. In this section, using a procedure similar to that which was used to extend the 1D a scheme to become the 1D a-e scheme, the two marching steps that form the 2D a-e scheme will be constructed by modifying the other equations in the 2D a scheme, i.e., Eqs. (4.66), (4.67), (4.69), and (4.70). As a prerequisite, first we shall provide a geometric interpretation of the procedure by which the second equation of the 1D a scheme, i.e., Eq. (2.8), was extended to become the second equation of the 1D a-e scheme, i.e., Eq. (2.13).

The key step in extending the 1D a scheme to the 1D a-e scheme is the construction of a central difference approximation of ¶u / ¶x at the mesh point (j,n). The approximation is given as the fraction within the parentheses on the extreme right side of Eq. (2.12). Consider a line segment in the x-u space joining the two points (x_j-1/2, u_j-1/2^¢ n) and (x_j+1/2, u_j+1/2^¢ n). It is obvious that the above central-difference approximation is the value of the slope du/dx of this line segment. In the following modification, instead of considering a line segment in the x-u space joining two points, we begin with the construction of a plane in the z-h- u space that intersects three given points.

To proceed, for any (j,k,n) Î W_q, q = 1,2, let