Abstract
In the numerical solution for nonlinear hyperbolic equations, numerical oscillation often shows and hides the real solution with the progress of computation. Using wavelet analysis, a dual wavelet shrinkage procedure is proposed, which allows one to extract the real solution hidden in the numerical solution with oscillation. The dual wavelet shrinkage procedure is introduced after applying the local differential quadrature method, which is a straightforward technique to calculate the spatial derivatives. Results free from numerical oscillation can be obtained, which can not only capture the position of shock and rarefaction waves, but also keep the sharp gradient structure within the shock wave. Three model problems—a one-dimensional dam-break flow governed by shallow water equations, and the propagation of a one-dimensional and a two-dimensional shock wave controlled by the Euler equations—are used to confirm the validity of the proposed procedure.
1. Introduction
Due to the nonlinearity, most problems governed by hyperbolic PDEs in fluid dynamics engineering have to be solved numerically. The main difficulty is that the solution of hyperbolic PDE are bound to develop discontinuities in finite time. A conventional numerical scheme, such as finite difference or finite volume, is directly applied to solve shock wave problem. Two serious problems may appear simultaneously or separately: (1) smearing out the shock wave gradually; (2) polluting the shock wave by numerical oscillation. In the recent several decades, people have developed many numerical schemes to maintain the sharp gradient structure of shock wave, while avoiding the numerical oscillation. Godunov [1] was credited with introducing the first Riemann solver for the Euler equations, by extending the previous Courant-Isaacson-Reeves (CIR) method to nonlinear systems of hyperbolic conservation law. To increase computation’s effectiveness, Roe [2] proposed a second order Roe solver through taking average of square root of density on double sides of a cell. Subsequent works were HLL and HLLC schemes [3]. Efforts in this field led to the proposal of total variance diminishing (TVD) [3, 4] and weighted essentially nonoscillating (WENO) schemes [5, 6]. Besides these sophisticated schemes, Shyy proposed a nonlinear filtering algorithm to eliminate numerical oscillation from second order central or upwind differencing in calculation of shock wave [7]. The filter has proved to be very effective in suppressing oscillation of short wavelength. However, the effect in removing the oscillation of long wavelength is not so promising. Kang introduced a multiresolution analysis (MRA) for increasing computation’s efficiency with preserving the high order numerical accuracy of a conventional solver [8]. Inspired by their work, we discovered wavelet’s application in suppressing numerical oscillation around the shock wave. Wavelet analysis is characterized by decomposing the signal to be analyzed into multiscale coefficients; high frequency component is described by coefficients on small scale and low frequency component is described by coefficients on large scale [9, 10]. In this way, shock wave may be maintained and the numerical oscillation around the shock wave may be removed after some special treatment for wavelet coefficients.
In this paper, we propose a dual wavelet shrinkage procedure to suppress numerical oscillation from a straightforward numerical scheme, named localized differential quadrature (LDQ) method, to calculate shock wave problem. LDQ, proposed by Zong, is a high order accurate numerical method. The LDQ method for spatial discretization combined with Runge-Kutta method for the time integration was used to solve Riemann problem [11, 12]; the wave profiles are smooth and nonbreaking. However, if LDQ method is directly used, high oscillation emerged. A dual wavelet transformation is then applied to process the highly oscillatory results. Three shock wave propagation problems governed by shallow water equations and Euler equations are presented for demonstration. The results are compared to their analytical solutions and very well agreement is achieved.
2. A Brief Introduction of LDQ
Hereinafter, only the main formulas of LDQ are introduced. Details of LDQ and their applications can be referred to in Zong [13]. The first step to LDQ method is to locate the neighborhood of a grid point of interest . We use , , to denote the distance between any two points in the solution domain. We find the permutation to satisfy .
It is clear that the points falling in the neighborhood of i-th point xi are the first m points. Denoting , , Si defines the neighborhood of the grid point of interest. We may get the first and second spatial derivatives at point xi of function by its neighborhood Si as follows [13].whereMultidimensional formulas share the same form as one-dimensional formulas because they are independent in each direction.
3. A Dual Wavelet Shrinkage Procedure
Wavelet analysis, introduced by Grossmann & Morlet in the early 1980s, has significantly impacted the signal and image processing [10] and numerical solving in nonlinear partial differential equations [14–16] as well as its application on turbulence [17–20]. Wavelet provides compactly supported, orthogonal or biorthogonal basis functions with adjustable smoothness. Due to the compactly supported basis, wavelet coefficients contain local structure information; hence wavelet analysis is a proper mathematical tool to process signal with local structure. Furthermore, wavelet analysis enables us to obtain multiscale information of the analyzed signal by introducing scale dilating.
Vanishing moment is an important property of wavelet analysis, which is directly related to wavelet basis functions’ smoothness. Wavelet function is called M order vanishing moment if the following relation is satisfied [16].Daubechies wavelet, one of the most important orthogonal wavelet categories, is widely used in many applications [21]. The function’s smoothness of Daubechies wavelet is adjustable by the vanishing moment order M. Let sup be scaling function and sup be wavelet function’s supported range. For example, the first or second order of Daubechies wavelet, denoted as DB1 or DB2, is shown in Figure 1.
(a)
(b)
Wavelet bases provide us with a multiscale resolution to express function. We assume function f(x) can be totally approached on the finest scaling index J; i.e., , and then it can be decomposed in terms of the sum of a series of wavelet bases functions as follows [19]. The coefficients are obtained by wavelet decomposition formula.For function with most part is well smooth, most wavelet coefficients will be small. Consequently, we can retain a good approximation even discarding a large number of wavelets with small coefficients. More precisely, if we rewrite the approximation as a sum of two termswhere is a prescribed particular threshold andfrom Vasilyev and Paolucci [15], the approximation error caused by the significant wavelets, coefficient amplitude above a given threshold , which will be introduced by adaptive wavelet shrinkage procedure subsequently, is bounded by the following restriction: and the number of significant wavelet coefficients is bounded by and wavelet’s vanishing moment as , in (14) and (15) depend on wavelet vanishing moment and function . Threshold has two effects: making approximation adaptively and controlling the approximation error globally. The similar situation can be simply extended to multidimensional space by tensor product.
For the signal to be analyzed, including unknown noise and unknown smoothness structures, how to remove the noise and keep the unknown structure is a complicated problem. Donoho and Johnstone proposed a wavelet shrinkage procedure to extract the structure from noisy sampled data [22]. We view numerical oscillation as the noise and discontinuity such as shock and rarefaction wave as the signal’s structure. Based on this understanding, the wavelet shrinkage procedure is applied to suppress the highly numerical oscillation obtained from the LDQ method. Indeed, wavelet shrinkage procedure can extract the unknown smoothness structures contaminated by heavy noise. Donoho and Johnstone also reported that the reconstruction is essentially as smooth as the mother wavelet [22]. This indicates that the reconstructed signal’s smoothness is closely related to the adopted wavelet basis functions’ smoothness. The reconstructed signal is locally similar to the adopted wavelet basis function, so the optimal wavelet basis should be in the same order smoothness as the real physical signal. Based on the consideration, we propose a dual wavelet shrinkage procedure using DB1 and DB2 for suppressing numerical oscillation obtained from LDQ method. DB1 and DB2 are, respectively, suitable for shock and rarefaction wave’s reconstruction, because DB1 is a sharp jump function and DB2 is a function with one-order smoothness. The dual wavelet shrinkage procedure is proposed as follows.(1)Decompose the numerical result with highly oscillation obtained from LDQ method via the discrete wavelet transform using DB1 by Mallat algorithm [19]; then the wavelet coefficients are obtained, where , is scale index, and , represent the largest and smallest scale, respectively.(2)Set as the threshold value at scale j, where is the standard deviation of coefficients and is the number of wavelet coefficients , i.e., . A threshold is assigned to each dyadic resolution level for threshold estimates, which is adaptive.(3)Threshold wavelet coefficients to get revised coefficients by the following.(4)Reconstruct denoised data via the wavelet reconstruction transform using the revised wavelet coefficients .(5)Repeat the above procedure using DB2. At next time step, the dual wavelet shrinkage procedure is complemented again.
The addition of this computational effort of the overall procedure is only of order as a function of sample size by Mallat algorithm, which is a fast transform similar to fast Fourier transform, bringing little extra computation.
4. Numerical Tests
4.1. One-Dimensional Dam-Break Flow
One-dimensional shallow water equations (SWEs) are a typical Riemann problem [23].SWEs describe the flow at time at point , where is the water depth, is the average horizontal velocity, and g is the gravitational acceleration. A wide variety of physical phenomena are governed by the SWEs, such as tidal flows in coastal water regions, bore wave propagation, flood waves in rivers, surges, and dam-break modeling [24]. Here we use SWEs to model dam-break problem in a channel of length L = 2000 m and the initial condition isThe dam collapses at t = 0 and the resulting bow consists of a shock wave traveling downstream and a rarefaction wave traveling upstream. LDQ is applied to the calculation of shock wave tube. In this research, the LDQ method is employed for spatial derivative calculation, where m = 5. Uniformly spaced nodes N = 256 are set in the channel length. Four-order Runge-Kutta method with time step Δt = 0.05 s is used for time integration. In order to check the dual wavelet shrinkage procedure’s effect, results without wavelet shrinkage procedure at time t = 25 s and t = 50 s are shown in Figure 2. Highly numerical oscillatory results are obtained only by LDQ method, especially after the shock wave front. Furthermore, the oscillation develops along with the time. Real physical solutions are hidden in the highly oscillatory numerical results. The task of the wavelet shrinkage procedure is to reconstruct them by suppressing the numerical oscillation adaptively.
(a)
(b)
The dual wavelet shrinkage procedure is subsequently implemented after LDQ method to remove the numerical oscillation. Figure 3 shows the numerical result of water depth along the channel at time t = 25 s and t = 50 s, which is also compared to the analytical result referred to by Stoker [23]. The highly numerical oscillation is removed out by this dual procedure. Furthermore, it remains the sharp gradient structure near the shock wave front simultaneously. However, small disturbance is introduced in the domain of rarefaction wave, which could be ignored.
(a)
(b)
Flow velocity to examine the efficiency of the proposed approach to remove unphysical velocity oscillations near the shock wave front is compared with that of existing exact Riemann solvers (Figure 4). The numerical oscillation is eliminated and the accuracy is good, which proves the validity of the dual wavelet method.
In order to examine the dual wavelet shrinkage’s superiority over single wavelet shrinkage, we use only DB1 and DB2 in the procedure, respectively. Figure 4 is the result of water depth only using DB2 at time t = 50 s. As shown in Figure 5(a), DB2 can extract the rarefaction wave’s structure well. However, DB2 introduces obvious disturbance with similar shape as DB2 function in the solution. Besides, DB2 makes the shock front excessively smooth, which does not comply with the physical truth. The result obtained from the procedure only using DB1 is shown in Figure 5(b). DB1 can reconstruct shock wave front with sharp gradient structure, but it also introduces many stepped shape errors in rarefaction wave and before shock wave front.
(a)
(b)
The dual procedure combines the respective advantages of DB1 and DB2. DB1 is a discontinuity function, which is the best candidate to capture the shock wave front. DB2, the second order vanishing moment function, is the optimal function to reconstruct the rarefaction wave. The dual wavelet procedure using DB1 and DB2 is a proper combination for processing highly numerical oscillatory results obtained from LDQ method in Riemann problem with shock wave.
4.2. One-Dimensional Shock Tube Problem
The Euler equations for one-dimensional unsteady ideal gas flow without heat conduction are given in conservation form. where , , , is gas’ density, velocity, energy per unit volume, and pressure, respectively. We need one more equation, i.e., the state equation, to close the system.Let ratio of specific heats γ = 1.4 and spatial range . Initial conditions are specified as follows.This is a shock tube problem, setting uniformly spaced nodes with number N = 512 in the tube length and m = 5 in LDQ method for spatial derivatives’ calculation. Fourth-order Runge-Kutta method is used for time integration with time step Δt = 0.005 s.
In order to reveal the effect of wavelet shrinkage procedure, the numerical results obtained from pure LDQ method are shown in Figure 6. Intense oscillation appears as in Figure 2.
After applying the dual wavelet shrinkage procedure to remove the numerical oscillation, the gas’s density, velocity, and pressure at time t = 5 s are shown in Figure 6, also including the analytical solutions. As shown in Figure 7, numerical results are very close to their analytical counterparts as referred by Toro [3].
4.3. Two-Dimensional Shock Wave Propagation
In order to confirm the procedure’s effect in two-dimensional shock wave’s computation, we consider two-dimensional Euler equations as follows.The initial conditions are specified as below. Similar calculation parameters are set as in one-dimensional case, i.e., uniformly 512 × 512 nodes in x and y direction, respectively; m = 5 in LDQ method; and fourth-order Runge-Kutta method for time with time step Δt = 0.005 s.
A pressure contour at time s is shown in Figure 8. Shock wave spreads outward with uniform speed and numerical results keep sharp interface. Pressure keeps the same at the two sides of the contact discontinuity and varies slowly in rarefaction wave region whose position can be seen obviously in Figure 9. It can be clearly seen that the wave spreads uniformly and symmetrically. However, slightly distorted deformation is introduced in radial direction by the effects of the grid which cannot represent a circle by square grid.
In order to get a clear version of density and pressure, we display their profiles at position of at time s in Figure 8. Symmetry is maintained in the calculation and results with nonoscillation are obtained.
This numerical model indicates that the dual wavelet shrinkage procedure can be extended to 2-D Riemann problem with shock wave and the effect on numerical suppressing is as promising as in 1-D case.
5. Conclusion
A dual wavelet shrinkage procedure is proposed to suppress the numerical oscillation for nonlinear hyperbolic equations, known as Riemann problem, with shock wave in fluid dynamics engineering in this paper. Three numerical tests, i.e., one-dimensional water dam breaking problem and one-/two-dimensional air shock wave propagating problems, are used to verify the procedure’s performance. High quality results are obtained both in capturing discontinuity and in suppressing the numerical oscillation. It is demonstrated that the dual wavelet procedure is a proper combination for processing highly numerical oscillation in Riemann problem.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
We thank the National Natural Science Foundation of China (Nos. 51679021, 51609031) and China Postdoctoral Science Foundation (No. 2016M601294).