Abstract
A high-resolution, total variation diminishing (TVD) stable scheme is derived for scalar hyperbolic problems using the method of flux limiters. The scheme was constructed by combining the 1st-order upwind scheme and the 3rd-order quadratic upstream interpolation scheme (QUICK) using new flux limiter function. The new flux limiter function was established by imposing several conditions to ensure the TVD properties of the scheme. For temporal discretization, the theta method was used, and values for the parameter θ were chosen such that the scheme is unconditionally stable. Numerical results are presented for one-dimensional pure advection problems with smooth and discontinuous initial conditions and are compared to those of other known numerical schemes. The results show that the proposed numerical method is stable and of higher order than other common schemes.
1. Introduction
Development of stable numerical schemes for modeling hyperbolic equations has been a crucial topic for decades of research in the area of computational analysis and simulation. This field of study poses several challenging difficulties. One of the main difficulties is the formation of different types of discontinuities in the solution. For nonlinear equations discontinuities can appear even for smooth initial conditions. Different numerical schemes were developed, which exhibit different stability and accuracy problems near these discontinuities. First-order schemes, for example, are known to be stable near discontinuities, but they pose an accuracy problem due to the large numerical dissipations. On the other hand, higher order linear schemes are less dissipative but they are known to exhibit nonphysical oscillations near discontinuities. This will be discussed further in Section 2 of this paper. A general form of a scalar conservation law in one dimension is given byEquation (1) represents a hyperbolic problem, with a characteristic speed of . This speed represents the propagation speed of information. The solution of this problem is constant along the characteristic curves . For nonlinear problems, the slopes (being dependent on the solution) will not be the same for different characteristic curves, which implies either an intersection or divergence in the - plane. The possibility of intersection means multiple values of the solution at some points and discontinuities are most likely to appear even for smooth initial conditions. In this context it is a necessity to solve hyperbolic problems using a numerical scheme capable of handling discontinuities.
Different approaches have been adapted to capture discontinuous solutions, one of which is adding an artificial viscosity term to the scheme. The addition of artificial viscosity can be achieved either by explicitly adding additional terms or implicitly within the numerical scheme used, such as 1st-order schemes (also called monotone schemes; see Section 2). Other approaches consider adaptive schemes depending on the smoothness of the solution. One way to achieve an adaptive scheme is by considering monotone schemes only at the regions of discontinuities and higher order schemes elsewhere. Examples on this approach are the works done by Harten [1], Sweby [2] and Kadalbajoo and Kumar [3]. Other examples of adaptive schemes are the ENO and WENO schemes, where the numerical stencils are adaptive in a way that discontinuities are avoided whenever possible [4–6]. The third approach for handling discontinuous solution is Godunov’s approach, where the nonlinearity property of the problem is explicitly included in the numerical scheme. In this approach, an approximation of the solution is constructed in the discretized spatial domain, and the evolution of the solution is solved at the interfaces of the spatial cells by solving the so-called Riemann problem exactly or approximately. For more details on the Riemann problem see [7, 8]. Examples on some numerical methods based on this approach can be found in [9, 10].
In this paper we present a systematical approach of deriving an adaptive high-resolution scheme for hyperbolic equations. In Section 2, we present numerical notation and concepts useful for our derivation. We proceed with the derivation in Section 3 and present numerical results for smooth and discontinuous solutions in Section 4. Some concluding remarks are posed in Section 5.
2. Numerical Notations
To allow for discontinuities along a discretized spatial and temporal domain, we need to admit weak solutions for (1) by integrating over . The conservative form of a numerical scheme in explicit formulation comes directly from this process; the numerical scheme for (1) can be written aswhereOther forms of an explicit numerical scheme are the general formand the incremental formwhereA numerical scheme of the general form (4) is called a monotone scheme if the operator is nondecreasing for all its arguments; namely,for any in the spatial domain. Monotone schemes are desired because they have the property of not introducing numerical oscillations with new extrema near discontinuities.
The total variation (TV) is also an important notion that gives a measure for oscillations due to a numerical scheme. For a given piecewise grid function that approximates an exact solution , the total variation is given by [7]A numerical method is called Total Variation Diminishing (TVD) if the total variation does not increase due to the numerical operator in (4), mathematicallyFor a numerical scheme in the incremental form (5), this condition is satisfied if [1, 11]On the relation between monotonicity and total variation, Harten suggests that a monotone scheme has a nonincreasing TV, and a scheme with nonincreasing TV preserves monotonicity [1].
Godunov proved that when numerical schemes are used for solving the advection equation, monotonicity cannot be achieved with linear schemes of order higher than one. This is referred to as Godunov’s order barrier theorem [12]. This means that high-order linear schemes such as Lax-Wendroff [13] and the quadratic upstream interpolation scheme (QUICK) [14] cannot be TVD, and using such schemes will introduce nonphysical oscillations whenever there is a discontinuity in the solution. On the other hand, 1st-order schemes are TVD, but they exhibit excessive dissipation and tend to drastically smear the solution, which in turn masks physical discontinuities and interfaces. To overcome the drawbacks of the two families of schemes, we seek a combination of the two in the form of a nonlinear high-order monotone scheme. This is achieved by imposing additional constraints to ensure the TVD properties of the scheme. In the next section we present our approach for deriving such schemes.
3. High-Resolution Scheme
In this section we present a systematical approach for deriving a high-resolution TVD scheme based on the concept of flux limiters. This is similar to the work of Sweby [2] and Kadalbajoo and Kumar [3]. The main idea is to combine the TVD 1st-order upwind scheme and the high-order accurate QUICK scheme by using the first one near discontinuities and the second one in smooth regions. The numerical scheme in (2) is used with a numerical flux , where is the numerical flux associated with a 1st-order upwind scheme, and is the numerical flux associated with 3rd-order QUICK scheme. Both numerical fluxes depend on the direction of characteristic speeds and consequently obey the hyperbolicity property of the problem.
The limiter function depends on the smoothness of the solution in the computational domain and is chosen such that the scheme in (2) is TVD. To find the TVD region of the scheme we consider the explicit formulation (), and we apply it to the linear scalar equation for simplicity of analysis; namely,For this equation the numerical flux due to the 1st-oder upwind scheme is defined asand for the 3rd-order QUICK schemeIf we substitute the numerical fluxes (for the case ) into (11) and (2) and use , we obtain the following scheme:where is the Courant number. With additional algebraic manipulation, this can be written as:where is referred to as the smoothness parameter. This is a numerical scheme written in an incremental form (5) with the following parameters:Implementing the conditions for the scheme to be TVD (10), we obtain the following inequality:At this point we impose another restriction to our scheme by requiring . This yields the following:For our scheme to be TVD at extreme points, we must have for . This condition along with (19) is satisfied for the bounds:Equation (20) defines the TVD region shaded in Figure 1. This region is bounded by the function . The same TVD region can be obtained for the case of negative advection velocity. For the numerical scheme to be TVD, the limiter function used has to lie in the TVD shaded area of Figure 1.
It is obvious that maximizing the flux limiter increases the antidiffusivity of the scheme. This means that a flux limiter function corresponding to the upper boundary of the TVD region yields the least diffusive scheme possible. Nevertheless, choosing such a limiter does not guarantee the highest possible order of accuracy. To ensure high-order of accuracy we require one additional condition; the limiter function should be chosen such that the scheme in (16) is 3rd-order accurate whenever possible. To impose this constraint, we investigated another numerical scheme and the 3rd-order upwind scheme. For this scheme, the numerical fluxes are defined asThe discretization for the advection equation (12) due to this scheme (for ) is given byWith algebraic manipulation, it can be shown that this scheme is equivalent to the scheme in (16) when the flux limiter function is chosen to be . Hence, we define the darker shaded region in Figure 1 as the desired 3rd-order TVD region for our scheme.
To guarantee all the conditions we imposed are satisfied, we derive a family of flux limiters lying in the 3rd-order TVD region (darker shaded region) whenever possible. Any flux function in that region can be expressed as an arithmetic average of two limiter functions: that corresponds to the QUICK scheme and for the 3rd-order upwind scheme; namely,Intersection points and are found to be and (see Figure 2). We can see that for any chosen value of , the limiter function satisfies our condition for smooth regions . The arithmetic average and the intersection points are shown in Figure 2.
Based on the above, we propose the following family of flux limiters characterized by the parameter :Equation (24) shows that the proposed limiter function is continuous and bounded for any value of the smoothness parameter . It also satisfies the condition for smooth regions . The flux limiter functions in (24) lie in the TVD region for any choice of ; hence, the scheme is TVD.
In this paper, the theta method (also called the generalized Crank-Nicolson method) is used, by considering a weighted average of explicit and implicit terms in the numerical scheme of (2). The scheme can be written as Equation (25) represents a family of schemes characterized by . We obtain the explicit scheme for , the implicit scheme for , and the Crank-Nicolson scheme for . The theta method is unconditionally stable for any choice of [15, 16]. Spatially, the accuracy of the scheme is determined by the way of calculating the numerical fluxes .
4. Numerical Results and Discussion
In this section we present two sets of numerical results, one for the case of smooth solutions and the other for discontinuous solutions. Both sets of results are shown for the case of pure advection with positive velocity. Results from the proposed scheme are also compared to those from other numerical schemes.
4.1. Smooth Initial Conditions
For this set of results we consider (12) with a smooth initial condition given byBoundary conditions for this problem are periodic with a convection velocity on a 1 m domain. Total time of the simulation was set to 1 sec, so that the solution at the end of the simulation is the same as the initial condition.
Figure 3 shows a comparison of the new scheme with other classical schemes and existing high-resolution schemes for the case of implicit time discretization (). Comparison was made for a grid size of 0.05 m and a Courant number of 0.2.
Excessive dissipation of the 1st-order upwind scheme can be observed and the solution is smeared out due to the inherent numerical viscosity of the scheme. For the two higher order linear schemes (2nd-order upwind scheme and Lax-Wendroff scheme) dissipative errors were minimized, but the signal in the numerical solution was either leading or lagging due to dispersion errors. This dispersive nature of the error can be observed by studying the modified equation for the corresponding schemes. For example, the modified equation corresponding to the 2nd-order upwind scheme isIt can be seen from (27) that the 2nd-order upwind scheme poses two types of errors: dissipation errors due to temporal discretization and dispersive errors due to both temporal and spatial discretization. The latter type of error causes either leading or lagging of the numerical solution, as shown in Figure 3.
On the other hand, more satisfactory results can be observed when implementing high- resolution schemes with different flux limiters. The new scheme with resulted in a better numerical solution as compared to other schemes.
Figure 4 shows a comparison of the same set of schemes for the case of Crank-Nicolson time discretization (). By examining (27) it is clear that the dissipative error terms cease to exist for the case of . This explains the reduction in dissipation for the high-order linear schemes (2nd-order upwind and Lax-Wendroff). On the other hand, our new scheme continues to produce the best results as compared to other schemes.
The exact solution was used to determine the accuracy of different numerical schemes. Accuracy was assessed on the basis of the norm, defined aswhere is the number of cells and is the numerical solution at the cell. Convergence rate is represented by a log-log plot of the norm of the error versus mesh size. Figure 5 shows the results of the spatial convergence for the case . Convergence rates for several considered schemes are listed in Table 1. The results show competitive rate of convergence for the new scheme as compared to other schemes.
4.2. Discontinuous Initial Conditions
For this set of results we consider (12) with a discontinuous initial condition given byBoundary conditions for this problem are periodic with a convection velocity on a 1 m domain, and the time of the simulation is 1 sec.
Two sets of results are shown: implicit time discretization () is shown in Figure 6 and Crank-Nicolson time discretization is shown in Figure 7. Comparison is made for a grid size of 0.0125 m and Courant number of 0.2.
Results show significantly more dissipation for the case of implicit discretization. For the 2nd-order upwind scheme and the QUICK scheme, dispersion error appears in the form of overshoots and undershoots. For the QUICK scheme, these oscillations were damped for the case of due to dissipation. High-resolution schemes, including the new scheme, give better results capturing the discontinuity, especially for the case with where dissipation is reduced.
Figure 8 shows results of the spatial convergence study for different schemes. Results of this study are listed in Table 2. We can observe a very good rate of convergence for the new scheme compared to the other schemes.
5. Concluding Remarks
A new high-resolution stable scheme is derived by implementing a hybridization of the monotone 1st-order upwind scheme and the quadratic upstream interpolation scheme (QUICK). For temporal discretization, the generalized Crank-Nicolson method characterized by the parameter was used. The combination of the 1st-order upwind scheme and QUICK scheme was done by means of a flux limiting function, which depends on the smoothness of the solution. Constraints were imposed to make the resulting scheme TVD and high-order at smooth regions of the solution. For the one-dimensional linear case, results from the proposed scheme for different values of were compared to classical and popular high-resolution schemes, and the convergence rate for the scheme was investigated for smooth and discontinuous solutions. The proposed scheme was shown to exhibit very good results as compared to other schemes. It is shown that the high-resolution schemes based on 2nd order Lax-Wendroff discretization exhibit a convergence rate of about 1 for smooth solutions and 0.5 for discontinuous solutions. The new scheme based on 3rd-order QUICK discretization exhibits a convergence rate of about 1.3 for smooth solutions and 0.8 for discontinuous solutions. In the case of discontinuous solutions, the new method derived in this paper shows the best convergence rate of all schemes investigated here.
Appendix
Finding the TVD Region for a Negative Velocity ()
The numerical scheme is given byFrom (13) and (14), we substitute the numerical fluxes to getDefining the smoothness parameter the above equation can be written asThis is an incremental form of the numerical scheme withLet ; we implement conditions in (10) to getFor ,One condition we require for the limiter function isThis condition yields the following:From which we obtain the two inequalities:This is the same result as found for the case of a positive velocity.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.