The numerical oscillation suppression method for the large time step scheme

The large time step scheme, which is commonly used in hydrodynamic numerical simulation, adopts the multielement reconstruction method in order to improve both the computational efficiency and accuracy; however, this results in numerical oscillation and calculation divergence. Therefore, it is necessary to develop methods to suppress the occurrence of numerical oscillation near the discontinuities. The numerical oscillation that is induced can be suppressed by using the random choice method (RCM) with the van der Corput random number sequence. In order to optimize the suppression performance, error analysis has been carried out in this paper for the calculation results from the RCM with the popular van der Corput random sequence combinations for selecting the optimal sequence combination. In order to further eliminate any short-wave numerical oscillations near the discontinuities, a method that combines the RCM with a nonlinear filter, which often used in aeroacoustic calculations, has been proposed in this paper to suppress short-wave oscillations. Specifically, the RCM and the nonlinear filter method were sequentially implemented in order to process the calculation results at each time step. The results showed that this method had a good suppression effect on the numerical oscillations near the discontinuities, which was demonstrated by the numerical experiments.The large time step scheme, which is commonly used in hydrodynamic numerical simulation, adopts the multielement reconstruction method in order to improve both the computational efficiency and accuracy; however, this results in numerical oscillation and calculation divergence. Therefore, it is necessary to develop methods to suppress the occurrence of numerical oscillation near the discontinuities. The numerical oscillation that is induced can be suppressed by using the random choice method (RCM) with the van der Corput random number sequence. In order to optimize the suppression performance, error analysis has been carried out in this paper for the calculation results from the RCM with the popular van der Corput random sequence combinations for selecting the optimal sequence combination. In order to further eliminate any short-wave numerical oscillations near the discontinuities, a method that combines the RCM with a nonlinear filter, which often used in aeroacoustic calculations, has been proposed in th...


I. INTRODUCTION
The shallow water equation is one of the macrogoverning equations that are commonly used in hydraulic simulation. The equation is mainly applied to the simulation of water flow in channels, dam-breaks, dike breaches, estuaries, shallow seas, and other water bodies. [1][2][3] The most significant difficulty in the use of the shallow water equation simulation is the presence of the discontinuity of physical quantities in shallow water flow. 4,5 The upwind method has been proven to be a method suitable for the discretization of shallow water equations. It is able to predict the water profile as well as the flow rate in hydraulic models, and it is widely used in the field of hydraulic engineering. 6,7 However, the traditional upwind scheme is limited by the Courant-Friedrichs-Lewy (CFL) condition in order to maintain numerical stability during the application process. 8 A mathematical model with an excessively small CFL number is impractical for practical applications. The instability of flow is studied by many scholars. 9,10 Some scholars have tried to increase the CFL number in order to overcome this difficulty, under the condition of stability for the numerical simulation.
Leveque first proposed the large time step (LTS) scheme; in his method, the characteristic wave is able to travel through multiple cells in one-time step using the traveling wave method, thus enabling the scheme to take a larger time step than the general explicit scheme and improving the calculation efficiency. 11 In terms of the functional properties of reducing the calculation time and maintaining reasonable accuracy, the LTS scheme has become increasingly popular in industrial applications. [12][13][14] Recently, Morales-Hernandez et al. 15 applied the approximate Riemann solution to the shallow water equation via the rarefaction wave approximation. As with other high-precision schemes, the LTS scheme with multielement reconstruction is able to improve the precision but results in numerical oscillation. 16 In order to overcome this difficulty, Xu et al. 17 explored the random choice method (RCM) and found that this method acts to damp the oscillation. At the ARTICLE scitation.org/journal/adv same time, it also generates dispersion and a shift in the resultant discontinuous position. Satisfactory calculation results of the RCM can only be achieved when the CFL number is small; as the CFL number increases, there will be some numerical oscillation at the discontinuity. The nonlinear filter technique is a shock wave capturing method that was proposed by Engquist et al. 18 It adopts a nonlinear filter factor to process the numerical solution at each time step, and it only acts on part of the grid points; therefore, the computational efficiency is relatively high. This method only works on short waves and only has a small effect on long waves. It can achieve a decent shock suppression effect on numerical dispersion in aeroacoustic calculations. 19 In this paper, the RCM and nonlinear filter from aeroacoustics have been combined to achieve the effect of suppressing the numerical oscillation at the discontinuity generated by the LTS scheme. Before this combination of methods was conducted, the random sequence combination in the RCM was statistically analyzed, and error analysis was performed to optimize the combination of sequences that is appropriate for the RCM.

II. USING THE LTS SCHEME TO SOLVE THE GOVERNING EQUATION OF THE ONE-DIMENSIONAL SHALLOW WATER EQUATION
A. One-dimensional shallow water model equation The parameters in Eq. (1) are as follows: where h, u, and g are the water depth, flow velocity, and acceleration due to gravity, respectively. t and x denote time and spatial distance, respectively. It should be noted that since the effects of bottom slope and the riverbed roughness on the water model equation are negligible, they were ignored in this study.

B. The LTS scheme for one-dimensional shallow water equations
The traditional method of solving nonlinear equations is to assume that the space variables are piecewise constants distributed in discrete grids and then update each cell i according to the contributions of the left and right interfaces, e.g., where λ is the wave speed of the characteristic wave, m/s; Δt is the time step, s; and Δx is the space step size, m. When the value of CFL = (λ j,R ) i− stability condition (Fig. 1). When |λj| Δt Δx > 1, the waves from the i − 1/2 interface and the i + 1/2 interface not only affect the ith cell but more interfaces are involved to affect the ith cell.
In order to realize the large time step (LTS) scheme in the onedimensional shallow water equation, as shown in Fig. 2, Xu combined the traveling wave method, based on Leveque's research, and Riemann's exact solution. An exact Riemann solver was used at the interface discontinuity between cells i and i + 1, and a new region S was derived between the left and right states Ui and U i+1 , where the physical variables (e.g., depth h and velocity u) were constants. In this way, the boundary between the S region and the right state generated a right-propagating wave with a wave velocity of λ i+ 1 2 ,R (λ i+ 1 2 ,R > 0), and the boundary between the S region and the left state generated a left-propagating wave with a wave velocity of λ i+ 1 With the framework of Leveque's method, the scheme can be extended to a larger time step where each wave propagates independently of the other waves. The strengths of the left-propagating and right-propagating waves that are generated at the interface discontinuity between cells i and i + 1 can be written as The number of cells passed by a left-propagating wave is equal The number of cells through which the right-propagating wave passes is equal to μ i+ 1 2 ]. There are μ i+ 1 2 ,L + 1 cells on the left side of the interface i + 1/2 that need to be updated in order to ensure that the μ i+ 1 2 ,L cell wave has passed completely. Therefore, this can be updated as follows: When the last cell wave has passed through, this can be updated to There are μ i+ 1 2 ,R + 1 cells on the right side of the interface i + 1/2 that need to be updated to ensure that the μ i+ 1 2 ,R cell wave has passed completely. Therefore, Eq. (6) can be updated as follows: When the last cell wave has passed through, it can be updated as follows:

III. NUMERICAL EXPERIMENTS
The large time step (LTS) scheme is able to accurately capture shock waves, while rarefaction waves appear in the Riemann problem, which requires the multiwave approximation. In order to fully understand how rarefaction problems are handled by the LTS scheme, a one-dimensional shallow water equation with a length of 100 m, a width of 1 m, an initial velocity of 0, and a discontinuous water surface was used as a test case. The initial conditions were as follows: The calculation range x ∈ [0, 100 m] was divided into 200 grids with a space step dx = 0.5 m. The initial time of dam-break is shown in Fig. 3(a). At the initial time, the water surface discontinuity (x = 50 m) generated a shock wave propagating to the right and a rarefaction wave propagating to the left [ Fig. 3(b)]. In the LTS scheme, rarefaction waves must be approximated by a finite number of shock waves [ Fig. 3(c)]. Higher accuracy can be achieved for the numerical solution, if approximated wave numbers are involved. The ten-wave approximation is able to obtain good computational results, which has been demonstrated by Xu, 17 therefore taking into consideration the accuracy and efficiency of calculation, the tenwave approximation is a good fit for the rarefaction waves in this paper.

IV. SUPPRESSION OF NONPHYSICAL OSCILLATION IN THE LTS SCHEME
A. The effect of grid number on numerical oscillation Xu found that the magnitude of the numerical oscillation at the plateau became larger as the CFL number increased.
However, the influence of time and grid number on the numerical oscillation was not considered in his paper. Based on his work, the calculation results of the RCM in the LTS scheme for different times and grid numbers were obtained, as shown in Fig. 4; from which it can be seen that both time and grid number affected the accuracy of the calculation. As time and grid number both increased, it was found that the RCM could not completely eliminate the numerical oscillation. In particular, when the number of grids was reduced, the resolution began to drop, especially at the discontinuities. A larger number of grids can improve the resolution but will also result in an increase in numerical oscillation at the plateau. Therefore, a compromise must be reached between the resolution and numerical oscillation; in this example, it was more appropriate to select 200 grids.

B. The random choice method
The LTS scheme is able to improve the calculation efficiency by increasing the time step; however, this results in undesired numerical oscillation. In order to suppress the numerical oscillation that is caused by the LTS scheme, Leveque adopted the Random Choice Method (RCM); this method selects a random number at each time step. If the random number is smaller than the fractional part of the CFL number that is transmitted to the cell, then the whole wave strength needs to be updated, otherwise it is discarded. For the last cell, the unified solution can be written as follows: Δx − μ) ≥ β n and β is the random number. In the RCM, the quality of the solution relies heavily on the choice of the random number β. A better choice of random number should be a sequence of numbers that is able to quickly reach an even distribution. By comparing the effects of different sampling schemes from aerodynamics on the performance of the RCM, Elperin 20 found that the best random sequence for application in the RCM is the van der Corput sequence. The van AIP Advances der Corput random sequence 21 is now generally applied to Leveque's RCM in order to analyze the influence of the parameters of the van der Corput random sequence on the results of the calculation. The van der Corput random number sequence {β n } is composed of k 1 and k 2 , which are two coprime and positive numbers, where k 1 > k 2 . The van der Corput random number sequence, composed of (k 1 , k 2 ), can be defined as where Ai = k 2 ai(mod k 1 ), n = m ∑ i=0 aik 1 i is the number of sequences, k 1 is the decimal number of n, and a i is the coefficient of n when it is expressed as k 1 ary. Toro et al. 4 investigated the statistics of the average valuex, the standard variance e, and the chi-square statistic i 2 for the commonly   Table I. Through comparison of the various (k 1 , k 2 ) combinations shown in Table I, it can be seen that the difference for the average value and the standard variance for each combination (k 1 , k 2 ) was relatively small compared to the variance in the chi-squared value. It was also found that there was no combination of values, which was able to optimize all three statistical characteristics. Therefore, their performance can only be assessed by trial and error. In order to further illustrate the effect of the simulation, the error function E was then introduced as follows: where N is the number of cells, CExact is the analytical solution of cell i, and C Numerical is the numerical solution of cell i. Table II shows the errors of the results of the calculation for various combinations of (k 1 , k 2 ) when CFL = 10 and time = 7 s. Figure 5(a) shows the results of the RCM-free calculation. Moreover, Figs. 5(b) and 5(c) show the results of the calculation for the RCM with the random number sequence combination (2, 1) and (7, 3), respectively. As shown in Table II and Fig. 5, among all the combinations of (k 1 , k 2 ), the combination (2, 1) resulted in the smallest calculation error, whereas the combination (7, 3) produced the largest calculation error. If the RCM was not utilized, there would be a large number of numerical oscillations at the plateau which would be located near the discontinuity. On the other hand, if the RCM with the combinations (2, 1) or (7, 3) was used, then the oscillation in the plateau was basically eliminated, although there was some residual numerical oscillation near the discontinuity. The oscillation amplitude for the combination (2, 1) was smaller than that for the combination (7,3). Although the chi-square statistic i 2 varied greatly for each combination, the results showed that it had little influence on the results of the calculation. The analysis of the error suggested that the highest accuracy could be achieved if the combination of the van der Corput's random number sequence in the RCM was (2, 1), which is consistent with the results of Olivier's research on the application of the RCM for a two-dimensional dam-break. 22 Table II shows the errors of the calculation results for each combination (k 1 , k 2 ) when CFL = 10 and time = 7 s.

C. Nonlinear filter method
The oscillation at the discontinuity is caused by numerical dispersion, which consists of a large number of short waves, which lowers the quality of the numerical solution. 23 Numerical oscillations that are caused by a short wave can be effectively filtered by the nonlinear filter method in aeroacoustics; however, it only works on short waves and has no effect on long waves. In this paper, the nonlinear filter method has been extended to the LST solution for shallow water equations. The nonlinear filter method and the RCM have been used to optimize the calculation results in order to achieve the effect of suppressing any numerical oscillation at the discontinuities.
Figures 6(a) and 6(b) show the original signals that were composed of a long wave as well as a short wave and the signals after the nonlinear filter, respectively. As shown in this figure, the nonlinear filter had a good suppression effect on the short wave but only had a small effect on the long wave. Consequently, the random selection method was combined with the nonlinear filter which was employed in sequence at each time step, as shown in Fig. 7. The results of the calculation that were obtained, when the CFL number was 10 and time was 7 s. The results after the RCM and nonlinear filter are shown in Fig. 8. The results after nonlinear filter are shown in Fig. 9.
From Figs. 5, 8, and 9, it can be seen that the simultaneous application of the random selection method and the nonlinear filter method in the LTS scheme was able to further inhibit the residual short wave numerical oscillation at the discontinuity and improve the accuracy of the numerical solution when compared with their individual application. In order to further verify the shock suppression effect of this method, the initial condition of dam-break was then changed to a one-dimensional cylindrical dam-break [Eq. (13)] and the results are shown in Fig. 10, As shown in Fig. 10, when the LTS scheme was adopted in order to calculate the one-dimensional cylindrical dam-break, it can be seen that there were a large number of numerical oscillations near the discontinuity. Although the RCM was able to suppress the oscillation in the calculation process, short wave oscillation still remained at the discontinuity. The use of the RCM and the nonlinear filter method in sequence at each time step was able to effectively eliminate the numerical oscillation at the discontinuity and boost the accuracy of the numerical solution.
The performance in three different scenarios, including the case where no numerical oscillation suppression methods were used as well as the cases for the RCM framework with and without nonlinear filter, was investigated for the relationship mapping between time and the error function E, as shown in Fig. 11. The root mean square error is shown in Table III. The results showed that the resultant error of the water depth increased slightly with time, and the water velocity error gradually decreased as time increased. This was because the initial water velocity was very close to zero, and the denominator of the error function E was small, which resulted in a larger value of E. Without the use of any numerical oscillation suppression methods, the resulting error function and root mean square error were the largest and it was the smallest when the RCM  and nonlinear filter method were combined. Therefore, it can be concluded that after the implementation of the RCM and the nonlinear filter method, the numerical error was able to be efficiently lowered.

V. CONCLUSION
Based on the research of Leveque, Xu, and others, the results in this paper have shown that the number of grids in the simulation should be within a reasonable range when the LTS scheme is implemented. If the grid is too sparse, the resolution of the numerical solution will be lower, especially at the discontinuity. On the other hand: while an over-dense grid will improve the resolution, this will result in the residual short-wave numerical oscillation being strengthened.
Although Leveque utilized the RCM in his research in order to suppress any oscillation in the LTS scheme, residual oscillation still occurred at the discontinuity. The RCM can be optimized by introducing the van der Corput random series, the performance of which was evaluated in this paper using error analysis of the calculation results; the results showed that the combination (2, 1) offered the highest accuracy.

ARTICLE scitation.org/journal/adv
In order to avoid any residual short-wave numerical oscillation at the discontinuity in the simulation, the nonlinear filtering method from aeroacoustics was integrated into the LTS scheme. At each time step, the RCM and the nonlinear filter method were executed sequentially in order to effectively eliminate the numerical oscillation at the discontinuity. Through the integration of the RCM and the nonlinear filter method, it was found that the overall calculation error can be effectively reduced and the accuracy of the numerical solution was improved through the implementation of relationship mapping between time and the error function E.