Journal of Computational Engineering

Volume 2015 (2015), Article ID 575380, 10 pages

http://dx.doi.org/10.1155/2015/575380

## Development of High-Resolution Total Variation Diminishing Scheme for Linear Hyperbolic Problems

^{1}Nuclear Engineering Department, Jordan University of Science and Technology, P.O. Box 3030, Irbid 22110, Jordan^{2}Department of Nuclear, Plasma and Radiological Engineering, University of Illinois at Urbana-Champaign. 216 Talbot Laboratory, 104 S. Wright Street, Urbana, IL 61801, USA

Received 15 March 2015; Revised 7 June 2015; Accepted 21 June 2015

Academic Editor: Jim B. W. Kok

Copyright © 2015 Rabie A. Abu Saleem and Tomasz Kozlowski. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.