- Email: [email protected]
The introduction of the surfing scheme for shock capturing with high-stability and high-speed convergence
Accepted Manuscript
The Introduction of the Surfing Scheme for Shock Capturing with High-Stability and High-Speed Convergence M. Mollaei , S.M. Malek Jafarian PII: DOI: Article Number: Reference:
S1007-5704(19)30175-3 https://doi.org/10.1016/j.cnsns.2019.104861 104861 CNSNS 104861
To appear in:
Communications in Nonlinear Science and Numerical Simulation
Received date: Revised date: Accepted date:
10 September 2018 20 May 2019 26 May 2019
Please cite this article as: M. Mollaei , S.M. Malek Jafarian , The Introduction of the Surfing Scheme for Shock Capturing with High-Stability and High-Speed Convergence, Communications in Nonlinear Science and Numerical Simulation (2019), doi: https://doi.org/10.1016/j.cnsns.2019.104861
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights:
AC
CE
PT
ED
M
AN US
CR IP T
The Surfing scheme has ultra-stability. This scheme has rapid convergence. This scheme captures shock waves without any extra terms like artificial viscosity. This scheme can capture shock waves with high resolution by increasing the number of nodes used for solving the Surfer part, similar to adaptive meshes. Although the number of equations is increased, the calculation time is decreased due to increasing the time step.
1
ACCEPTED MANUSCRIPT
The Introduction of the Surfing Scheme for Shock Capturing with HighStability and High-Speed Convergence M. Mollaei a, S.M. Malek Jafarian a* a
Department of Mechanical Engineering, University of Birjand, Iran.
CR IP T
Abstract
M
AN US
In this paper, a novel shock capturing method with high stability and speed of solution is presented. The method divides the dependent variables into two parts, i.e. Surf and Surfer, like surfing sport. The Surf part is calculated by applying the concept of Sobolev gradient on the dependent variables. Furthermore, it has the least difference with these variables, along with the lowest discretization error and the same equation. Hence, it is expected that this part can be solved at high stability conditions and large time steps. In addition, the Surfer part is obtained from the difference between the Surf part and the dependent variables. Therefore, only the error-maker operators (for example, a discontinuity like a shock) are in this part and accordingly it has a local nature. Due to this feature, it can be solved with fewer points compared to the Surf part. Therefore, its equation is solved as quickly as Surf equation, despite its limited stability conditions. Finally, the summation of the Surf and Surfer solutions makes the main solution. The method has been applied to the inviscid Burgers’ equation with an initial value of the step function, and one-dimensional inviscid Euler flow through a nozzle with a shock wave. The stability condition for the Burgers’ equation has been increased to a Courant number of 1010 and the solution time for the Euler flow case has been decreased by one-fifth. Comparison of the results with traditional methods indicates the ultra-high stability and ultra-speed convergence of this method.
ED
Keywords: Shock Capturing; Sobolev Gradient; Surf; Surfer
PT
1. Introduction
AC
CE
In order to obtain a mathematical model for a fluid, its properties are considered as differentiable variables in fluid mechanics. However, at supersonic flows, that upstream information does not have enough time to move from upstream to downstream; discontinuities are formed in the fluid properties and are called shock waves. Therefore, the differentiable assumption is not correct in shock regions. This causes errors that are produced and propagated in the numerical solution of the relevant equations. Many numerical methods have been developed by investigators to overcome this and capture shock waves. These methods and a historical overview of them is presented briefly in the following. Shock capturing methods are classified into two categories of classical and modern methods. Classical methods are based on the symmetric or central difference and modern methods are based on upwind schemes. Symmetric or center-based methods, unlike modern methods, do not use any information from the sign of characteristic velocity and wave propagation direction for discretization. For capturing shock waves, both classical and modern methods require a certain amount of artificial viscosity for stability to prevent non-physical oscillations. For classical methods one can point to the
*
Corresponding author. E-mail addresses: ([email protected]) ([email protected]) 2
ACCEPTED MANUSCRIPT MacCormack, Beam and Warming, Lax and Lax–Wendroff methods [1]; and for modern methods one can point to high order TVD (Total Variation Diminishing), FCT (Flux-Corrected Transport), MUSCL (Monotonic Upstream-Centered Schemes for Conservation Laws), ENO (Essentially NonOscillatory) and PPM (Piecewise Parabolic Method) [2].
AC
CE
PT
ED
M
AN US
CR IP T
In classical methods, Ohwada et al. [3] increased the efficiency of the kinetic Lax–Wendroff method using an equilibrium distribution function of gas molecules in the region of shock waves. Since classical methods cannot function well in modeling strong shock waves, modern methods have been developed. Modern methods are second-order in the continuous region and non-oscillatory or monotonic in the shock region [4]. To eliminate oscillations, these methods use artificial viscosity in nonlinear form and locally. The non-oscillatory nature of the TVD based upwind method causes that these methods to become the basis of modern methods. Among these methods, one can point to “Essentially Non-Oscillatory” methods. The ENO method, presented by Harten et al. [5] for solving the Euler equations, is the first method in this field and discretizes flow fluxes without using weight coefficients. Yang [6,7] presented and used second- and third-order ENO method for solving Euler equations. Liu et al. [8] presented WENO method by adding weight coefficients in high-order ENO to enhance interpolation accuracy of this method. Both of these methods worked very well in capturing shock waves, and they were widely used because of the precision of these methods could be determined by the user; but in these methods, small-scale shock waves may be removed and nonphysical oscillations may be generated. Jameson [9,10] presented CUSP (Convective Upwind and Split Pressure) method with high-speed convergence and completely upwind for flows with high Mach numbers that use characteristic splitting and scalar diffusive fluxes by adding a pressure difference to momentum equation. Adams and Shariff [11] combined ENO method with the compact method to pre-calculate node-base flux derivatives and used it for the simulation of the interaction between turbulent boundary layers and shocks. Wang [12] developed ENO-Padé scheme by combining the ENO algorithm with the cell-centered Padé scheme to eliminate the nonphysical oscillatory of the Padé scheme across discontinuities and to improve the performance of the ENO scheme in smooth regions. Furthermore, similar hybrid methods were developed by Deng and Maekawa [13], Deng and Zhang [14], Jiang et al. [15]. Pirozzoli [16] combined WENO and Compact methods and presented a new method which uses explicit methods to approximate fluxes. Chen [4], after studying various high-order methods, developed a method with combining WENO method and the quasi-implicit Crank-Nicolson method to increase the stability of Pirozzoli method and applied it to free-surface flows. Sigalotti et al. [17] presented a method for capturing shock waves with the combination of high-order Godunov method and SPH-equations to increase the accuracy of the results. Kawai and Lele [18] developed a simple and effective artificial diffusion method with the aim of capturing discontinuities in curvilinear and anisotropic meshes by high-order methods. Kumar and Kadalbajoo [19] studied a class of simple and efficient high-order TVD methods for conservation laws and by introducing anti-diffusive terms to first-order three-point monotone methods and improved the efficiency of these methods in capturing shock waves. Djavareshkian and Jahdi [20] combined TVD and ACM (Artificial Compression Method) methods to enhance stability and compared their results on steady inviscid flows for various Mach numbers from subsonic to supersonic. Klöckner et al. [21] presented a high-order method using the concept of artificial viscosity for detecting shock waves based on Galerkin methods and examined various issues. Jiwari [22] proposed an efficient numerical scheme based on uniform Haar wavelets and the quasilinearization process for the numerical simulation of time dependent nonlinear Burgers’ equation. Yee et al. [23] investigated the relationship between numerical dissipation and discontinuity propagation speed in high-order methods of capturing shock waves. Yu and Yan [24] developed a method based on the Galerkin discontinuity and artificial diffusion of the Kawai model to calculate high order fluxes in unstructured grids. Kumar and Pandit [25] proposed a composite numerical scheme based on finite difference and Haar wavelets for solving time dependent coupled Burgers’ equation with appropriate initial and boundary conditions. Fu et al. [26] presented the MWCS (Modified Weighted Compact 3
ACCEPTED MANUSCRIPT
AN US
CR IP T
Scheme) method for the better approximation of fluxes on Euler equations which combined WENO and WCS (Weighted Compact Scheme) methods in order to eliminate oscillations and local calculation of weight coefficients for pressure and density at each time step and showed the reduction of calculation costs compared to the WENO method. Premasuthan et al. [27,28] applied an artificial viscosity approach to spectral difference and simulated discontinuities in compressible flows with better performance on unstructured and curvilinear grids. Jiwari [29] developed a hybrid numerical scheme based on Euler implicit method, quasilinearization and uniform Haar wavelets for the numerical solutions of Burgers’ equation. Pandit et al. [30] developed a method by extending the Jiwari’s method [22] which is reliable for planar and non-planar nonlinear Burgers’ equation and yields results better than other methods and compatible with the exact solutions. Kotov et al. [31] presented a method for reducing the amount of numerical viscosity in high-order shock capturing methods for turbulent flows at low velocities. In this method, the amount of numerical viscosity, which is modeled by LES methods, were determined with a sensor automatically. Yang et al. [32] presented a high-order method to calculate fluxes to capture shock waves on unstructured dynamic grids with artificial diffusivity, and studied its results for flow on an airfoil, moving vortex flows and supersonic flow over a bump. Mazaheri and Nishikawa [33] developed second- and third-order nonoscillatory methods for capturing shock waves for triangular nonlinear grids in problems that are discontinuous. Their results were for the two-dimensional viscous Burgers’ equation. Mittal and Pandit [34] formulated a novel technique for the numerical solutions of Shock wave Burgers’ equations for planar and non-planar geometry by using robustness of wavelets generated by dilation and translation of Haar wavelets on third scale to capture the sensitivity information.
CE
PT
ED
M
As mentioned, researchers are investigating methods that reduce calculation costs and increase stability and accuracy of shock capturing methods. Surfing scheme, which is presented in this work, is a novel method with these goals. The method, similar to the surfing sport on water, splits the dependent variables into two parts Surf and Surfer. The Surf part has the least difference with the variables, along with the lowest discretization error and the same equation. Hence, its solution is an approximation of final solution and therefore it exists across the solution field. Due to its low discretization error, it is expected that this part can be solved with a very high stability condition and very large time step. The Surfer part, which is obtained from the difference between the dependent variables and the Surf part, only contains all of the shock-maker operators. Therefore, it exists only in the shock region and its equation can be solved locally by the prevalent methods. Hence, its calculation cost, similar to the Surf part, is very low. Finally, the main solution is obtained by the summation of the Surf and Surfer solutions. The results of this scheme for the inviscid Burgers’ equation with an initial value of the step function, and one-dimensional inviscid flow through a nozzle with a shock wave is presented in the Results section. 2. Numerical Scheme
AC
A popular equation in the field of shock capturing study, which is used as a model equation to investigate various solution procedures, is the one-dimensional inviscid Burgers’ equation as a classical nonlinear first-order hyperbolic equation [1]. This equation is: (1)
Where is the velocity component that is convected along the one-dimensional space ( -direction) and is time. The initial and boundary conditions are: (
)
( )
(2) ( ) (3) The initial profile ( ), which spans the entire flow domain , -, with the constant boundary value , provide the unity solution of the nonlinear Burgers’ equation. To solve the equation (1) 4
ACCEPTED MANUSCRIPT numerically, the second term is discretized with various schemes. Most of these schemes are based on central and upwind discretization methods. Depending on which type of discretization methods is used, errors are produced in the form of high-order derivatives of the dependent variables. By using the simple central differencing scheme, equation (1) is discretized as follows: (
)
(4)
The -direction error in equation (4) is obtained by using Taylor's series expansion as below: )
(
)
(
)
)
(5)
CR IP T
( (
The dominant term of error in equation (5) can be identified by using the scale-up method. For comparison the magnitude of all terms in this equation, the magnitude scales of and are selected for and respectively, which can be specified after defining a problem. Therefore, the scales of all components can be obtained by putting these magnitude scales into all terms as: )
(
)
AN US
}
(
(
(6)
)
{ According to scales in (6), they can be sorted as follow: (
)
)
(
(
)
)
(
(
)
)
(
)
(7)
M
(
CE
PT
ED
The comparison in (7) shows that the first term in error, which contains the third-order derivative of , has a large scale compared to the scales of other terms. Hence, it is the dominant error and its order is approximately the same as the order of . This error term has a significant amount where is discontinuous or has a sudden variation. Hence, if it is not controlled, it spreads to all the domain and diverges the solution. Many numerical methods add some artificial viscosity that only reduces the effects of this error locally and do not affect any other parts of the domain. It seems that if has a variation in the -direction that the value of this term can be neglected compared to the main terms of equation (4). The error does not reach its maximum amount and does not spread on the domain. The Surfing scheme uses this fact to control the propagation of errors. 3. Surfing Scheme
AC
According to the previous section, a variable like can include sections where higher-order derivatives can be partially significant. Since in a numerical simulation, the generation of discretization errors is unavoidable, the errors have significant amounts in these sections and should be controlled to achieve the convergence. In another view, one can look at any variables like , which is a combination of two parts, the first part that contains all sections which cause errors generation and propagation such as discontinuous, plus the second part that contains the remaining without any errors. Consequently, if the error-maker operators are separated from the dependent variables, the outcome can be solved with an errorless solution. In addition, the separated parts can be solved along with the outcome solution. For an example of error-maker operators, it can be mentioned as discontinuities like shock waves. Therefore, the existing discontinuity in a variable (Fig. 1 (a)) can be considered as a soft part (Surf in Fig. 1 (b)), which carries the discontinuity-maker (Surfer in Fig. 1 (b)). This is similar to the surfing sport that the surfer moves with the surf of sea (Fig. 1 (c)). 5
ACCEPTED MANUSCRIPT (a)
(b)
(c) Surf
Surfer
Soft Part Hard Part
Discontinuity
Fig. 1. Schematic diagram of the Surfing scheme.
AN US
CR IP T
Hence, this scheme splits the base variable into two Surf and Surfer parts. The Surf part is calculated by applying the concept of the Sobolev gradient in this paper. The method of calculation is such that this part has the lowest discretization error. In addition, to approximate the base variable, the form of the base equation and Surf equation are alike; and Surf value has the least difference with the value of the base variable. This part covers all the solution domain and solves with a large time step due to its stability. Also, it can be solved by using mesh-based methods (see Fig. 2 (a)) and its time is very small. The Surfer is also obtained from the difference between the base variable and the Surf part. This part exists only in the discontinuity region. The stability condition for this part is smaller than the Surf part because all the discontinuity-makers are in this part. Therefore, it cannot use a large time step for the solution. Due to the local nature of this part, it can use very few points compared to the Surf part. Consequently, the solution time is also very small for this part. It can be solved locally by point-cloud methods. If the discontinuity-makers move, this part can be solved by using methods like SPH. For an example, points in Fig. 2 (a), which are located in the shock regions that can be seen in Fig. 2 (b), would be used for the solution of the Surfer part. Eventually, the final solution is calculated from the summation of the Surf and Surfer solutions; and its solution time is decreased marvelously in comparison with the traditional methods. (b)
M
(a)
ED
Surf Mesh
Shock Regions
PT
Surfer Points
CE
Fig. 2. Influence region of Surf and Surfer parts in a two-dimensional case.
For the case of
in equation (1), it can be split into Surf and Surfer parts as follows:
AC
(8) Where and are denoted the Surf (soft part) and Surfer (hard part) respectively. By putting (8) into (1), the following equation is obtained: (9)
In order to obtain the Surf equation, it must be considered that this equation has a form similar to the base equation. With rearranging the terms in equation (9), this equation can be rewritten as: (
)
(
)
Finally, equation (10) can be split into two equations (Surf and Surfer) as follow:
6
(10)
ACCEPTED MANUSCRIPT
(
)
(11) (
)
(12)
3.1.
CR IP T
Equations (11) and (12) are the Surf and Surfer equations respectively. Where function is an arbitrary function that can control errors of the solution. By using this function, the error can be distributed between these two equations. Since the Surf part approximates and its equation is of the same form as the base equation. In this paper, the function should be considered zero. In addition, according to the Surfer ( ) exists only in the discontinuity region; equation (12) can be solved locally by point-cloud methods. If the discontinuity moves, this equation can be solved by using methods like SPH. Surf part
AN US
As mentioned, this part should have very high stability condition and never diverges. In fact, a numerical solution is called diverged when its result differs from the analytical solution. This difference occurs when the mathematical equation is discretized to be solved numerically. The reason is that errors are generated in the discrete operations. For this part, equation (11) can be discretized similar to equation (4) as follows: (
)
(13)
Similar to equation (5), the -component of error for this equation can be obtained as: ( (
)
(
)
(
)
)
(14)
Sobolev gradient
AC
3.1.1.
CE
PT
ED
M
By using the Taylor's series expansion, the only difference between the discretized equation (13) and equation (11) is the error terms which the -component of error is obtained in equation (14). If this error has negligible value, the solution of equation (13) is exactly equal to the solution of equation (11) with a high stability condition and without any divergence. Therefore, in order to ensure that the Surf part has very high stability condition, the dominant term of equation (14) should have negligible value which achieved by controlling the variation subject to . For this purpose, the maximum amounts of error terms in (14) can be determined similarly to (6) and (7). It is important that the dominant error never reaches its maximum. Thus, variation subject to must be such that it does not happen and this means that the solution never diverges. Since the error at time approximates the error at time , by controlling variation at this time, it can be achieved. For this purpose, the Sobolev gradient concept is used for determining in this paper, which is described in the next section. Sobolev gradient ̅ can be calculated from a gradient , by applying a diffusion term as follows: ̅
(
̅
)
(15)
The Sobolev gradient ̅ is smoother than the gradient and its value tends to the straight line between two boundaries with increasing as illustrated in Fig. 3. Therefore, the maximum and minimum values of ̅ move to zero if the values of the boundaries are zero (see Fig. 3 (b)).
7
ACCEPTED MANUSCRIPT
(a)
(b) x2 2 x 2 x x2
2 (3)
Sobolev Gradient (x s /3!)
Sobolev Gradient (s)
x2 2 x 2 x x2
x
x
CR IP T
Fig. 3. Sobolev gradient applied on: (a) a discontinuity of with various and (b) a zero-boundary function ( dominant error) with various .
Since the errors cannot be eliminated, it can be negligible compared to the main terms of the original equation by reducing its order. To do this, and ̅ are replaced by and dominant errors in equation (15) respectively as follows: *
)+
AN US
(
(16)
The order of the terms of the equation (16) is obtained by using scale-up as follows:
}
(
)+
(17)
M
*
PT
→
ED
{ Where is the magnitude scale of which is the same order of and shows the maximum value that can take. The maximum value of dominant error should be times smaller than the order of dominant error and consequently, the error can be negligible by choosing an optimized . It can be achieved by applying (17) to equation (16): (18)
AC
CE
Therefore, if is equal to , the maximum value of dominant error can be controlled by . After specifying the , can be calculated from by solving equation (16). This equation has a fifthorder derivative of and its solution is usually time-consuming in two and three dimensions. However, subject to the local behavior of Surfer part, this equation also can be solved locally to reduce the calculation time. In this case, the equation (16) can be integrated to reduce the order of derivatives as below: (19)
Using central difference method, equation (19) can be discretized at time step
as follows: (20)
Equation (20) is a pentadiagonal system of equations. system.
8
is obtained at time step
by solving this
ACCEPTED MANUSCRIPT 3.2.
Surfer part
Surfer part contains all the error-maker operators. In this part, the operators can be treated as a discontinuity, where the traditional discretization methods that based on differentiable variables use a certain amount of artificial viscosity to control error propagation and guarantee the convergence. In addition, to decrease the solution time similar to equation (16), the Surfer equation can be solved locally due to the local nature of this part. Equation (12) can be written as: (
)
(21)
CR IP T
The right-hand side of equation (21) has an acceleration role in the solution process. The central difference scheme is used for discretization: (22)
3.3.
AN US
The Suring scheme was explained using the Burgers’ equation as a simple case, at the present and previous sections. The scheme is applied to the one-dimensional Euler equations (system of equations containing mass, momentum, and energy equations) as a more complex case, in the following. It also should be noted that the surfing scheme can be applied to two- and three-dimensional system of equations. Euler Equations
The one-dimensional form of Euler equations are as follows:
M
(
)
(23-c)
are density, velocity, and total energy, respectively. Also,
PT
Where , , and obtained from:
(23-b)
ED
(
)
(
(23-a)
)(
)
is pressure and it is
(24)
AC
CE
Where is the ratio of specific heats of fluid. Equations (23), which span the entire flow domain , -, require two boundary conditions at the inlet and one at the outlet for a subsonic inflow and outflow [35] as below:
Where , obtained from
( ) (25-a) ( ) (25-b) ( ) (25-c) , and are inlet density, inlet pressure, and outlet pressure, respectively. is at the boundaries and by using equation (24) as follows: (26)
The initial conditions for equations (23) are as: ( (
) )
(27-a) (27-b) 9
ACCEPTED MANUSCRIPT ( ) (27-c) Where , , and are the initial values of density, velocity and pressure of the domain, respectively. In order to use the Surfing scheme, the flow variables are split into two parts, i.e. Surf and Surfer: (28-a) (28-b)
CR IP T
(28-c) (28-d) Where and are the Surf and Surfer parts of the relevant flow properties, respectively. To obtain the Surf and Surfer equations, equations (28) can be inserted into equations (23) similar to what was done in equation (9), but the mathematical calculation is complex. Since the Surf part approximates the flow variables and its equations should be the same form as Euler equations, it can be chosen similar to equations (23) as below: (29-a)
)
AN US
(
[(
)
(29-b)
]
(29-c)
]
(
(30-a) -
)
*
(
(30-b)
)
+
-
(30-c)
CE
,[
PT
,
ED
M
The initial values of the Surf variables (Eq. (29)) are obtained by applying the Sobolev gradient concept to the initial flow variables (similar to Eq. (16)). The Surfer equations are obtained by the difference between the Euler and Surf equations. Consequently, these equations are as follow:
Finally, the flow variables are calculated from (28) by the summation of the solution of equations (29) and (30).
AC
4. Results 4.1.
Burgers’ Equation
The first example is a discontinuity which is described by the functions (31) and (32). These functions are used as initial condition (2) and boundary conditions (3) for the Burgers’ equation to investigate the discontinuity propagation throughout the domain. (
)
{
(31)
( ) (32) As mentioned in the previous section, the Surf part approximates the base variable. This part at time step is obtained from equation (20). Table 1 shows the effects of various (at ) on the order of 10
ACCEPTED MANUSCRIPT error terms. As expected, these orders are reduced times compared to the error terms of ( achieved by ). Therefore, the errors can be controlled by . In this paper, the value of has been chosen equal to 10. Table 1: Effects of various (
)
(
(
4.2 0.42 0.054 0.0059
(
)
0.78 0.067 0.0076 0.00077
)
(
-0.15 -0.011 -0.0013 -0.00013
)
(
0.022 0.0016 0.00017 0.000017
)
0.0031 0.00020 0.000021 0.0000022
CR IP T
12 1.5 0.21 0.023
)
at on order of error terms.
AN US
To investigate the stability of a numerical method, the Courant–Friedrichs–Lewy (CFL) condition is ( ⁄ ) for one-dimensional case and is known as the Courant used. It is defined as number. Where is equal to one in this case, according to the initial condition (31). Fig. 4 (a) shows the Surf solutions at the first time step for various Courant numbers. As shown, the solution has never diverged for Courant numbers up to 10 10 and has a very high stability. The discretization error has not grown up. This is because of the dominant error is negligible compared to the main terms and has been controlled. Fig. 4 (b) shows Surfer solutions at first time step for various Courant numbers. These solutions modify the Surf part solutions for presenting the results as the final exact solution. As shown in this figure, the values of this part are decreased by increasing the Courant number because of negative acceleration on the right-hand side of equation (21). In addition, oscillations are introduced in the solution due to the presence of central difference dispersion error. (a)
1.2
0.1
s
0.6 0.4 0.2 0 -0.2
0
1
2
3
0.05 0
-0.05 -0.1 4
1
1.5
2
2.5
3
x
PT
x
ED
-1
C=10 0 C=10 C=102 C=104 6 C=10 8 C=10 C=1010
h
0.8
M
1
(b)
-1
C=10 0 C=10 C=102 4 C=10 6 C=10 C=108 C=1010
Fig. 4. (a) Surf and (b) Surfer solutions at first time step for various Courant numbers.
CE
To validate the Surfing scheme, the results are compared to the solutions of Lax, Lax-Wendroff, MacCormack, FCT, and first and second order TVD methods for Courant number of one in Fig. 5. As shown in these figures, the Surfing results are in a good agreement with the results of these methods. (a)
AC
1.2
(b)
1.2
0.8
0.8
0.6
0.6
u
1
u
1
0.4
0.4
0.2 0 -0.2
0.2
Lax Lax-Wendroff MacCormack Surfing
0
1
FCT (Beam and Warming) First-Order TVD Second-Order TVD Surfing
0 2
3
-0.2
4
x
0
1
2
3
4
x
Fig. 5. Comparison of Surfing solution at Courant number of one with: (a) Lax, Lax-Wendroff and MacCormack, (b) FCT, first and second order TVD solutions. 11
ACCEPTED MANUSCRIPT Fig. 6 shows the ultra-stability of the Surfing scheme compared to the above-expressed methods. The Courant number is 1000. As shown, all other methods have been diverged, while the present method has maintained its stability. This indicates that the Surfing scheme has superiority to existing methods. (a)
1.2
(b)
1.2
0.6
0.6
u
1 0.8
u
1 0.8
0.4
0.2 0 -0.2
0.2
Lax Lax-Wendroff MacCormack Surfing
0
1
First-Order TVD Second-Order TVD Surfing FCT (Beam and Warming)
0 2
3
-0.2
4
CR IP T
0.4
0
x
1
2
3
4
x
Fig. 6. Comparison of Surfing solution at Courant number of with: (a) Lax, Lax-Wendroff and MacCormack, (b) FCT, first and second order TVD solutions.
AN US
Fig. 7 (a) and (b) show solution of Burgers’ equation at five time steps for and Courant numbers with the presented method, respectively. As shown, the discontinuity moves forward with advancing in time. In addition, it moves with larger space step by increasing the Courant number (that equivalent to time step for a constant space step). The main solution has been tended asymptotically to the Surf part one, due to reducing the amount of the Surfer part. In addition, the ultra-stability of this scheme can be seen in these figures. (a)
1.2
1.2 1
0.8
0.4
Initial Step 1 Step 2 Step 3 Step 4 Step 5
0 -0.2
0
1
2
3
PT
0.2
0.6 0.4
ED
u
0.6
0.8
u
M
1
(b)
Initial Step 1 Step 2 Step 3 Step 4 Step 5
0.2 0 -0.2
4
x
0
1
2
3
4
x
CE
Fig. 7. Surfing solutions of Burgers’ equation at five time steps for: (a) Courant number of number of .
, (b) Courant
AC
In the previous results, was chosen as zero in order to the form of equation (11) and (1) can be the same and hence can approximate . In order to investigate the effects of the arbitrary function in ) *( + has chosen for a test case. Fig. 8 shows the equations (11) and (12), the function effects of various on the solutions for Courant number of at the first step. It illustrates that can control errors.
12
ACCEPTED MANUSCRIPT
(a)
(b)
1.2 1
1 0.8
0.6
u
u
0.8
0.6
0.4
0.4
0.2
0.2
0
0
-0.2
0
1
2
3
1.2
-0.2 1.8
4
1.9
2
x
2.2
2.3
2.4
x
on results of the Surfing scheme for the Courant number of magnification of the effects.
in the first step, (b)
CR IP T
Fig. 8. (a) Effects of
2.1
The second example is described by the functions (33) and (34). These functions are used as initial condition (2) and boundary conditions (3) for the Burgers’ equation. , ⁄ ( ) { (33) (
)
(a)
1.2
1.2
M
Initial Step 1 Step 2 Step 3 Step 4 Step 5
1 0.8
0
2
PT
0.2
4
6
Initial Step 1 Step 2 Step 3 Step 4 Step 5
1 0.8
u
u 0.4
(b)
0.6
ED
0.6
-0.2 0
AN US
(34) Fig. 9 (a) and (b) show solution of Burgers’ equation at five time steps for and of Courant numbers with the presented method for initial and boundary conditions of (33) and (34), respectively. As shown, a discontinuity is formed and moves forward with advancing in time. In addition, it is formed faster when the Courant number is increased. The ultra-stability of this scheme can be seen in these figures.
0.4 0.2 0 -0.2 0
8
2
x
4
6
8
x
4.2.
CE
Fig. 9. Surfing solutions of Burgers’ equation with initial condition of (33) at five time steps for: (a) Courant number of , (b) Courant number of .
Euler Equations
AC
The boundary conditions (25) and initial conditions (27), which are used for solving the equations (23), are as follow: (
)
(
)
( ⁄(
(
)
(
( (
)
(
)) )
)
(35)
)
(36-a)
(36-b) ( ) (36-c) Where , is equal to 1.4 for air. Fig. 10 (a) shows the geometry of a nozzle. The nozzle throat is located at ⁄ , where L is the length of the nozzle. As shown in Fig. 10 (b) through (d), there is choking at the nozzle throat, after that, a normal shock wave has been formed at the divergent part of 13
ACCEPTED MANUSCRIPT the nozzle. These figures show that the Surfing scheme has a good agreement with SCDS and SLAU2 [36] results. (a)
0.5
(b) SCDS SLAU2 Surfing
0.4
SCDS SLAU2 Surfing
2
Area
1.5
u
0.3
1
0.2
0
0.2
0.4
0.6
0.8
0 0
1
0.2
x
0.6
0.8
1
(d)
1.2
SCDS SLAU2 Surfing
1
0.4
x
(c)
1.2
CR IP T
0.5 0.1
SCDS SLAU2 Surfing
1 0.8
0.8
0.6
AN US
p
0.6
0.4
0.4
0.2
0.2 0 0
0
0.2
0.4
0.6
0.8
1
x
-0.2 0
0.2
0.4
0.6
0.8
1
x
M
Fig. 10. Solution of one-dimensional Euler equations: (a) geometry of the nozzle, (b) velocity distribution, (c) density distribution, (d) pressure distribution.
Table 2: Comparison of the computational times. Methods Computational Time (millisecond) Surfing 528 SCDS 2903 SLAU2 2716
AC
CE
PT
ED
Table 2 shows the comparison of computational times of the Surfing, SCDS, and SLAU2 methods. All results have been obtained by single-thread processing by using a 64-bit operating system with 2.5 GHz processor speed and 8 GB of RAM. It seems that the computational time must be increased due to the additional calculations of the Surfer solution and smoothing procedure. Since the Surf equation is solved with a very large time step compared to SCDS and SLAU2, convergence time is decreased for this part. As shown, the computation time of the presented method has been decreased to one-fifth.
5. Conclusion
The results of Surfing scheme were presented in the previous section. As mentioned, the stability condition of Burgers’ equation has been investigated up to 1010 and the results marvelously show the ultra-high stability of the presented method. According to this, the solution process could use large time steps in comparison with the traditional methods. This increased the convergence speed in the solution of equations. This was evident in solving the Euler equations that the solution time has been decreased to one-fifth. Certainly, the advantages and convergence speed of presented scheme in twoand three-dimensional cases become more obvious compared to one-dimensional case. The advantages of this method are listed as follow: 14
ACCEPTED MANUSCRIPT -
The Surfing scheme has ultra-stability. This scheme has rapid convergence. This scheme captures shock waves without any extra terms like artificial viscosity. This scheme can capture shock waves with high resolution by increasing the number of nodes used for solving the Surfer part, similar to adaptive meshes. Although the number of equations is increased, the calculation time is decreased due to increasing the time step.
6. Acknowledgement
CR IP T
The authors would like to thank Professor Earl Dowell (Department of Mechanical Engineering and Materials Science, Duke University) for his invaluable suggestions and fruitful comments which greatly improved the manuscript. It is gratefully acknowledged.
References
AC
CE
PT
ED
M
AN US
[1] K.A. Hoffmann and S.T. Chiang, “Computational Fluid Dynamics Volume I”, Fourth Edition, 2000. [2] H. Liu, “Shock capturing methods”, Advanced Topics in Numerical Methods, 2008. [3] T. Ohwada, R. Adachi, K. Xu, and J. Luo, “On the Remedy against Shock Anomalies in Kinetic Schemes”, Journal of Computational Physics, 106–129, 2013. [4] C. Chen, “High Order Shock Capturing Schemes for Hyperbolic Conservation Laws and the Application in Open Chanel Flows”, Ph.D. Dissertation, University of Kentucky, 2006. [5] A. Harten, S. Osher, B. Engquist, and S.R. Chakravarthy, “Some Results on Uniformly High-Order Accurate Essentially Nonoscillatory Schemes.”, Applied Numerical Mathematics, 347–376, 1986. [6] J.Y. Yang, “Second-Order Accurate Essentially Non-Oscillatory Schemes for the Euler Equations”, AIAA Journal, Vol.28, No.12, 2069–2076, 1990. [7] J.Y. Yang, “Third-Order Nonoscillatory Schemes for the Euler Equations”, AIAA Journal, 1611–1618, 1991. [8] X.D. Liu, S. Osher, and T. Chan, “Weighted Essentially Non-Oscillatory Schemes”, Journal of Computational Physics, 200–212, 1994. [9] A. Jameson, “Analysis and Design of Numerical Schemes for Gas Dynamics 1 Artificial Diffusion, Upwind Biasing, Limiters and Their Effect on Accuracy and Multigrid Convergence”, International Journal of Computational Fluid Dynamics, 171–218, 1995. [10] A. Jameson, “Analysis and Design of Numerical Schemes for Gas Dynamics 2 Artificial Diffusion and Discrete Shock Structure”, International Journal of Computational Fluid Dynamics, 1–38, 1996. [11] N.A. Adams and K. Shariff, “A High-Resolution Hybrid Compact-ENO Scheme for Shock-Turbulence Interaction Problems”, Journal of Computational Physics, 27–51, 1996. [12] Z.P. Wang, “An Essentially Non-Oscillatory High Order Padé-Type (ENO-Padé) Scheme”, Ph.D. Dissertation, University of Kentucky, 2002. [13] X. Deng and H. Maekawa, “Compact High-Order Accurate Nonlinear Schemes”, Journal of Computational Physics, 77–91, 1997. [14] X. Deng and H. Zhang, “Developing High-Order Weighted Compact Nonlinear Schemes”, Journal of Computational Physics, 22–44, 2000. [15] L. Jiang, H. Shan, and C. Liu, “Weighted Compact Schemes for Shock Capturing”, International Journal of Computational Fluid Dynamics, 147–155, 2001. [16] S. Pirozzoli, “Conservative Hybrid Compact-WENO Schemes for Shock Turbulence Interaction”, Journal of Computational Physics, 81–117, 2002. [17] L.D.G. Sigalotti, H. Lopez, A. Donoso, E. Sira, and J. Klapp, “A Shock-Capturing SPH Scheme Based on Adaptive Kernel Estimation”, Journal of Computational Physics, 124–149, 2006. [18] S. Kawai and S.K. Lele, “Localized Artificial Diffusivity Scheme for Discontinuity Capturing on Curvilinear Meshes”, Journal of Computational Physics, 9498–9526, 2008. [19] R. Kumar and M.K. Kadalbajoo, “A Class of High Resolution Shock Capturing Schemes for Hyperbolic Conservation Laws”, Applied Mathematics and Computation, 110–126, 2008.
15
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
[20] M.H. Djavareshkian and M.H.A. Jahdi, “Shock-Capturing Method using Characteristic-Based Dissipation Filters in Pressure-Based Algorithm”, Acta Mechanica, 2010. [21] A. Klöckner, T. Warburton, and J.S. Hesthaven, “Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method”, Mathematical Modelling of Natural Phenomena, 57–83, 2011. [22] R. Jiwari, “A Haar wavelet quasilinearization approach for numerical simulation of Burgers’ equation”, Computer Physics Communications, 2413–2423, 2012. [23] H.C. Yee, D.V. Kotov, W. Wang, and C.W. Shu, “Spurious Behavior of Shock-Capturing Methods: Problems Containing Stiff Source Terms and Discontinuities”, International Conference on Computational Fluid Dynamics, 2012. [24] J. Yu and C. Yan, “An Artificial Diffusivity Discontinuous Galerkin Scheme for Discontinuous Flows”, Computers & Fluids, 56–71, 2013. [25] M. Kumar and S. Pandit, “A composite numerical scheme for the numerical simulation of coupled Burgers’ equation”, Computer Physics Communications, 809–817, 2014. [26] H. Fu, Z. Wang, Y. Yan, and C. Liu, “Modified Weighted Compact Scheme with Global Weights for Shock Capturing”, Computers & Fluids, 165–176, 2014. [27] S. Premasuthan, C. Liang, and A. Jameson, “Computation of Flows with Shocks Using the Spectral Difference Method with Artificial Viscosity, I: Basic Formulation and Application”, Computers & Fluids, 111–121, 2014. [28] S. Premasuthan, C. Liang, and A. Jameson, “Computation of Flows with Shocks Using the Spectral Difference Method with Artificial Viscosity, II: Modified Formulation with Local Mesh Refinement”, Computers & Fluids, 122–133, 2014. [29] R. Jiwari, “A hybrid numerical scheme for the numerical solution of the Burgers’ equation”, Computer Physics Communications, 59–67, 2015. [30] S. Pandit, M. Kumar, R.N. Mohapatra and A.S. Alshomrani, “Shock waves analysis of planar and non planar nonlinear Burgers’ equation using Scale-2 Haar wavelets”, International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 27 Issue: 8, 1814–1850, 2016. [31] D.V. Kotov, H.C. Yee, A.A. Wray, B. Sjögreen, and A.G. Kritsuk, “Numerical Dissipation Control in High Order Shock-Capturing Schemes for LES of Low Speed Flows”, Journal of Computational Physics, 189– 202, 2016. [32] J. Yang, B. Zhang, C. Liang, and Y. Rong, “A High-Order Flux Reconstruction Method with Adaptive Mesh Refinement and Artificial Diffusivity on Unstructured Moving/Deforming Mesh for Shock Capturing”, Computers & Fluids, 17–35, 2016. [33] A. Mazaheri and H. Nishikawa, “High-Order Shock-Capturing Hyperbolic Residual-Distribution Schemes on Irregular Triangular Grids”, Computers & Fluids, 29–44, 2016. [34] R.C. Mittal and S. Pandit, “Sensitivity Analysis of Shock Wave Burgers’ Equation via a Novel Algorithm Based on scale-3 Haar Wavelets”, International Journal of Computer Mathematics, 601–625, 2017. [35] C.A.J. Fletcher, “Computational Techniques for Fluid Dynamics 2”, Springer Series in Computational Physics, 1988, Page 175. [36] K. Kitamura and E. Shima, “Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for AUSM-family schemes”, Journal of Computational Physics, 62–83, 2013.
16
The Introduction of the Surfing Scheme for Shock Capturing with High-Stability and High-Speed Convergence M. Mollaei , S.M. Malek Jafarian PII: DOI: Article Number: Reference:
S1007-5704(19)30175-3 https://doi.org/10.1016/j.cnsns.2019.104861 104861 CNSNS 104861
To appear in:
Communications in Nonlinear Science and Numerical Simulation
Received date: Revised date: Accepted date:
10 September 2018 20 May 2019 26 May 2019
Please cite this article as: M. Mollaei , S.M. Malek Jafarian , The Introduction of the Surfing Scheme for Shock Capturing with High-Stability and High-Speed Convergence, Communications in Nonlinear Science and Numerical Simulation (2019), doi: https://doi.org/10.1016/j.cnsns.2019.104861
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights:
AC
CE
PT
ED
M
AN US
CR IP T
The Surfing scheme has ultra-stability. This scheme has rapid convergence. This scheme captures shock waves without any extra terms like artificial viscosity. This scheme can capture shock waves with high resolution by increasing the number of nodes used for solving the Surfer part, similar to adaptive meshes. Although the number of equations is increased, the calculation time is decreased due to increasing the time step.
1
ACCEPTED MANUSCRIPT
The Introduction of the Surfing Scheme for Shock Capturing with HighStability and High-Speed Convergence M. Mollaei a, S.M. Malek Jafarian a* a
Department of Mechanical Engineering, University of Birjand, Iran.
CR IP T
Abstract
M
AN US
In this paper, a novel shock capturing method with high stability and speed of solution is presented. The method divides the dependent variables into two parts, i.e. Surf and Surfer, like surfing sport. The Surf part is calculated by applying the concept of Sobolev gradient on the dependent variables. Furthermore, it has the least difference with these variables, along with the lowest discretization error and the same equation. Hence, it is expected that this part can be solved at high stability conditions and large time steps. In addition, the Surfer part is obtained from the difference between the Surf part and the dependent variables. Therefore, only the error-maker operators (for example, a discontinuity like a shock) are in this part and accordingly it has a local nature. Due to this feature, it can be solved with fewer points compared to the Surf part. Therefore, its equation is solved as quickly as Surf equation, despite its limited stability conditions. Finally, the summation of the Surf and Surfer solutions makes the main solution. The method has been applied to the inviscid Burgers’ equation with an initial value of the step function, and one-dimensional inviscid Euler flow through a nozzle with a shock wave. The stability condition for the Burgers’ equation has been increased to a Courant number of 1010 and the solution time for the Euler flow case has been decreased by one-fifth. Comparison of the results with traditional methods indicates the ultra-high stability and ultra-speed convergence of this method.
ED
Keywords: Shock Capturing; Sobolev Gradient; Surf; Surfer
PT
1. Introduction
AC
CE
In order to obtain a mathematical model for a fluid, its properties are considered as differentiable variables in fluid mechanics. However, at supersonic flows, that upstream information does not have enough time to move from upstream to downstream; discontinuities are formed in the fluid properties and are called shock waves. Therefore, the differentiable assumption is not correct in shock regions. This causes errors that are produced and propagated in the numerical solution of the relevant equations. Many numerical methods have been developed by investigators to overcome this and capture shock waves. These methods and a historical overview of them is presented briefly in the following. Shock capturing methods are classified into two categories of classical and modern methods. Classical methods are based on the symmetric or central difference and modern methods are based on upwind schemes. Symmetric or center-based methods, unlike modern methods, do not use any information from the sign of characteristic velocity and wave propagation direction for discretization. For capturing shock waves, both classical and modern methods require a certain amount of artificial viscosity for stability to prevent non-physical oscillations. For classical methods one can point to the
*
Corresponding author. E-mail addresses: ([email protected]) ([email protected]) 2
ACCEPTED MANUSCRIPT MacCormack, Beam and Warming, Lax and Lax–Wendroff methods [1]; and for modern methods one can point to high order TVD (Total Variation Diminishing), FCT (Flux-Corrected Transport), MUSCL (Monotonic Upstream-Centered Schemes for Conservation Laws), ENO (Essentially NonOscillatory) and PPM (Piecewise Parabolic Method) [2].
AC
CE
PT
ED
M
AN US
CR IP T
In classical methods, Ohwada et al. [3] increased the efficiency of the kinetic Lax–Wendroff method using an equilibrium distribution function of gas molecules in the region of shock waves. Since classical methods cannot function well in modeling strong shock waves, modern methods have been developed. Modern methods are second-order in the continuous region and non-oscillatory or monotonic in the shock region [4]. To eliminate oscillations, these methods use artificial viscosity in nonlinear form and locally. The non-oscillatory nature of the TVD based upwind method causes that these methods to become the basis of modern methods. Among these methods, one can point to “Essentially Non-Oscillatory” methods. The ENO method, presented by Harten et al. [5] for solving the Euler equations, is the first method in this field and discretizes flow fluxes without using weight coefficients. Yang [6,7] presented and used second- and third-order ENO method for solving Euler equations. Liu et al. [8] presented WENO method by adding weight coefficients in high-order ENO to enhance interpolation accuracy of this method. Both of these methods worked very well in capturing shock waves, and they were widely used because of the precision of these methods could be determined by the user; but in these methods, small-scale shock waves may be removed and nonphysical oscillations may be generated. Jameson [9,10] presented CUSP (Convective Upwind and Split Pressure) method with high-speed convergence and completely upwind for flows with high Mach numbers that use characteristic splitting and scalar diffusive fluxes by adding a pressure difference to momentum equation. Adams and Shariff [11] combined ENO method with the compact method to pre-calculate node-base flux derivatives and used it for the simulation of the interaction between turbulent boundary layers and shocks. Wang [12] developed ENO-Padé scheme by combining the ENO algorithm with the cell-centered Padé scheme to eliminate the nonphysical oscillatory of the Padé scheme across discontinuities and to improve the performance of the ENO scheme in smooth regions. Furthermore, similar hybrid methods were developed by Deng and Maekawa [13], Deng and Zhang [14], Jiang et al. [15]. Pirozzoli [16] combined WENO and Compact methods and presented a new method which uses explicit methods to approximate fluxes. Chen [4], after studying various high-order methods, developed a method with combining WENO method and the quasi-implicit Crank-Nicolson method to increase the stability of Pirozzoli method and applied it to free-surface flows. Sigalotti et al. [17] presented a method for capturing shock waves with the combination of high-order Godunov method and SPH-equations to increase the accuracy of the results. Kawai and Lele [18] developed a simple and effective artificial diffusion method with the aim of capturing discontinuities in curvilinear and anisotropic meshes by high-order methods. Kumar and Kadalbajoo [19] studied a class of simple and efficient high-order TVD methods for conservation laws and by introducing anti-diffusive terms to first-order three-point monotone methods and improved the efficiency of these methods in capturing shock waves. Djavareshkian and Jahdi [20] combined TVD and ACM (Artificial Compression Method) methods to enhance stability and compared their results on steady inviscid flows for various Mach numbers from subsonic to supersonic. Klöckner et al. [21] presented a high-order method using the concept of artificial viscosity for detecting shock waves based on Galerkin methods and examined various issues. Jiwari [22] proposed an efficient numerical scheme based on uniform Haar wavelets and the quasilinearization process for the numerical simulation of time dependent nonlinear Burgers’ equation. Yee et al. [23] investigated the relationship between numerical dissipation and discontinuity propagation speed in high-order methods of capturing shock waves. Yu and Yan [24] developed a method based on the Galerkin discontinuity and artificial diffusion of the Kawai model to calculate high order fluxes in unstructured grids. Kumar and Pandit [25] proposed a composite numerical scheme based on finite difference and Haar wavelets for solving time dependent coupled Burgers’ equation with appropriate initial and boundary conditions. Fu et al. [26] presented the MWCS (Modified Weighted Compact 3
ACCEPTED MANUSCRIPT
AN US
CR IP T
Scheme) method for the better approximation of fluxes on Euler equations which combined WENO and WCS (Weighted Compact Scheme) methods in order to eliminate oscillations and local calculation of weight coefficients for pressure and density at each time step and showed the reduction of calculation costs compared to the WENO method. Premasuthan et al. [27,28] applied an artificial viscosity approach to spectral difference and simulated discontinuities in compressible flows with better performance on unstructured and curvilinear grids. Jiwari [29] developed a hybrid numerical scheme based on Euler implicit method, quasilinearization and uniform Haar wavelets for the numerical solutions of Burgers’ equation. Pandit et al. [30] developed a method by extending the Jiwari’s method [22] which is reliable for planar and non-planar nonlinear Burgers’ equation and yields results better than other methods and compatible with the exact solutions. Kotov et al. [31] presented a method for reducing the amount of numerical viscosity in high-order shock capturing methods for turbulent flows at low velocities. In this method, the amount of numerical viscosity, which is modeled by LES methods, were determined with a sensor automatically. Yang et al. [32] presented a high-order method to calculate fluxes to capture shock waves on unstructured dynamic grids with artificial diffusivity, and studied its results for flow on an airfoil, moving vortex flows and supersonic flow over a bump. Mazaheri and Nishikawa [33] developed second- and third-order nonoscillatory methods for capturing shock waves for triangular nonlinear grids in problems that are discontinuous. Their results were for the two-dimensional viscous Burgers’ equation. Mittal and Pandit [34] formulated a novel technique for the numerical solutions of Shock wave Burgers’ equations for planar and non-planar geometry by using robustness of wavelets generated by dilation and translation of Haar wavelets on third scale to capture the sensitivity information.
CE
PT
ED
M
As mentioned, researchers are investigating methods that reduce calculation costs and increase stability and accuracy of shock capturing methods. Surfing scheme, which is presented in this work, is a novel method with these goals. The method, similar to the surfing sport on water, splits the dependent variables into two parts Surf and Surfer. The Surf part has the least difference with the variables, along with the lowest discretization error and the same equation. Hence, its solution is an approximation of final solution and therefore it exists across the solution field. Due to its low discretization error, it is expected that this part can be solved with a very high stability condition and very large time step. The Surfer part, which is obtained from the difference between the dependent variables and the Surf part, only contains all of the shock-maker operators. Therefore, it exists only in the shock region and its equation can be solved locally by the prevalent methods. Hence, its calculation cost, similar to the Surf part, is very low. Finally, the main solution is obtained by the summation of the Surf and Surfer solutions. The results of this scheme for the inviscid Burgers’ equation with an initial value of the step function, and one-dimensional inviscid flow through a nozzle with a shock wave is presented in the Results section. 2. Numerical Scheme
AC
A popular equation in the field of shock capturing study, which is used as a model equation to investigate various solution procedures, is the one-dimensional inviscid Burgers’ equation as a classical nonlinear first-order hyperbolic equation [1]. This equation is: (1)
Where is the velocity component that is convected along the one-dimensional space ( -direction) and is time. The initial and boundary conditions are: (
)
( )
(2) ( ) (3) The initial profile ( ), which spans the entire flow domain , -, with the constant boundary value , provide the unity solution of the nonlinear Burgers’ equation. To solve the equation (1) 4
ACCEPTED MANUSCRIPT numerically, the second term is discretized with various schemes. Most of these schemes are based on central and upwind discretization methods. Depending on which type of discretization methods is used, errors are produced in the form of high-order derivatives of the dependent variables. By using the simple central differencing scheme, equation (1) is discretized as follows: (
)
(4)
The -direction error in equation (4) is obtained by using Taylor's series expansion as below: )
(
)
(
)
)
(5)
CR IP T
( (
The dominant term of error in equation (5) can be identified by using the scale-up method. For comparison the magnitude of all terms in this equation, the magnitude scales of and are selected for and respectively, which can be specified after defining a problem. Therefore, the scales of all components can be obtained by putting these magnitude scales into all terms as: )
(
)
AN US
}
(
(
(6)
)
{ According to scales in (6), they can be sorted as follow: (
)
)
(
(
)
)
(
(
)
)
(
)
(7)
M
(
CE
PT
ED
The comparison in (7) shows that the first term in error, which contains the third-order derivative of , has a large scale compared to the scales of other terms. Hence, it is the dominant error and its order is approximately the same as the order of . This error term has a significant amount where is discontinuous or has a sudden variation. Hence, if it is not controlled, it spreads to all the domain and diverges the solution. Many numerical methods add some artificial viscosity that only reduces the effects of this error locally and do not affect any other parts of the domain. It seems that if has a variation in the -direction that the value of this term can be neglected compared to the main terms of equation (4). The error does not reach its maximum amount and does not spread on the domain. The Surfing scheme uses this fact to control the propagation of errors. 3. Surfing Scheme
AC
According to the previous section, a variable like can include sections where higher-order derivatives can be partially significant. Since in a numerical simulation, the generation of discretization errors is unavoidable, the errors have significant amounts in these sections and should be controlled to achieve the convergence. In another view, one can look at any variables like , which is a combination of two parts, the first part that contains all sections which cause errors generation and propagation such as discontinuous, plus the second part that contains the remaining without any errors. Consequently, if the error-maker operators are separated from the dependent variables, the outcome can be solved with an errorless solution. In addition, the separated parts can be solved along with the outcome solution. For an example of error-maker operators, it can be mentioned as discontinuities like shock waves. Therefore, the existing discontinuity in a variable (Fig. 1 (a)) can be considered as a soft part (Surf in Fig. 1 (b)), which carries the discontinuity-maker (Surfer in Fig. 1 (b)). This is similar to the surfing sport that the surfer moves with the surf of sea (Fig. 1 (c)). 5
ACCEPTED MANUSCRIPT (a)
(b)
(c) Surf
Surfer
Soft Part Hard Part
Discontinuity
Fig. 1. Schematic diagram of the Surfing scheme.
AN US
CR IP T
Hence, this scheme splits the base variable into two Surf and Surfer parts. The Surf part is calculated by applying the concept of the Sobolev gradient in this paper. The method of calculation is such that this part has the lowest discretization error. In addition, to approximate the base variable, the form of the base equation and Surf equation are alike; and Surf value has the least difference with the value of the base variable. This part covers all the solution domain and solves with a large time step due to its stability. Also, it can be solved by using mesh-based methods (see Fig. 2 (a)) and its time is very small. The Surfer is also obtained from the difference between the base variable and the Surf part. This part exists only in the discontinuity region. The stability condition for this part is smaller than the Surf part because all the discontinuity-makers are in this part. Therefore, it cannot use a large time step for the solution. Due to the local nature of this part, it can use very few points compared to the Surf part. Consequently, the solution time is also very small for this part. It can be solved locally by point-cloud methods. If the discontinuity-makers move, this part can be solved by using methods like SPH. For an example, points in Fig. 2 (a), which are located in the shock regions that can be seen in Fig. 2 (b), would be used for the solution of the Surfer part. Eventually, the final solution is calculated from the summation of the Surf and Surfer solutions; and its solution time is decreased marvelously in comparison with the traditional methods. (b)
M
(a)
ED
Surf Mesh
Shock Regions
PT
Surfer Points
CE
Fig. 2. Influence region of Surf and Surfer parts in a two-dimensional case.
For the case of
in equation (1), it can be split into Surf and Surfer parts as follows:
AC
(8) Where and are denoted the Surf (soft part) and Surfer (hard part) respectively. By putting (8) into (1), the following equation is obtained: (9)
In order to obtain the Surf equation, it must be considered that this equation has a form similar to the base equation. With rearranging the terms in equation (9), this equation can be rewritten as: (
)
(
)
Finally, equation (10) can be split into two equations (Surf and Surfer) as follow:
6
(10)
ACCEPTED MANUSCRIPT
(
)
(11) (
)
(12)
3.1.
CR IP T
Equations (11) and (12) are the Surf and Surfer equations respectively. Where function is an arbitrary function that can control errors of the solution. By using this function, the error can be distributed between these two equations. Since the Surf part approximates and its equation is of the same form as the base equation. In this paper, the function should be considered zero. In addition, according to the Surfer ( ) exists only in the discontinuity region; equation (12) can be solved locally by point-cloud methods. If the discontinuity moves, this equation can be solved by using methods like SPH. Surf part
AN US
As mentioned, this part should have very high stability condition and never diverges. In fact, a numerical solution is called diverged when its result differs from the analytical solution. This difference occurs when the mathematical equation is discretized to be solved numerically. The reason is that errors are generated in the discrete operations. For this part, equation (11) can be discretized similar to equation (4) as follows: (
)
(13)
Similar to equation (5), the -component of error for this equation can be obtained as: ( (
)
(
)
(
)
)
(14)
Sobolev gradient
AC
3.1.1.
CE
PT
ED
M
By using the Taylor's series expansion, the only difference between the discretized equation (13) and equation (11) is the error terms which the -component of error is obtained in equation (14). If this error has negligible value, the solution of equation (13) is exactly equal to the solution of equation (11) with a high stability condition and without any divergence. Therefore, in order to ensure that the Surf part has very high stability condition, the dominant term of equation (14) should have negligible value which achieved by controlling the variation subject to . For this purpose, the maximum amounts of error terms in (14) can be determined similarly to (6) and (7). It is important that the dominant error never reaches its maximum. Thus, variation subject to must be such that it does not happen and this means that the solution never diverges. Since the error at time approximates the error at time , by controlling variation at this time, it can be achieved. For this purpose, the Sobolev gradient concept is used for determining in this paper, which is described in the next section. Sobolev gradient ̅ can be calculated from a gradient , by applying a diffusion term as follows: ̅
(
̅
)
(15)
The Sobolev gradient ̅ is smoother than the gradient and its value tends to the straight line between two boundaries with increasing as illustrated in Fig. 3. Therefore, the maximum and minimum values of ̅ move to zero if the values of the boundaries are zero (see Fig. 3 (b)).
7
ACCEPTED MANUSCRIPT
(a)
(b) x2 2 x 2 x x2
2 (3)
Sobolev Gradient (x s /3!)
Sobolev Gradient (s)
x2 2 x 2 x x2
x
x
CR IP T
Fig. 3. Sobolev gradient applied on: (a) a discontinuity of with various and (b) a zero-boundary function ( dominant error) with various .
Since the errors cannot be eliminated, it can be negligible compared to the main terms of the original equation by reducing its order. To do this, and ̅ are replaced by and dominant errors in equation (15) respectively as follows: *
)+
AN US
(
(16)
The order of the terms of the equation (16) is obtained by using scale-up as follows:
}
(
)+
(17)
M
*
PT
→
ED
{ Where is the magnitude scale of which is the same order of and shows the maximum value that can take. The maximum value of dominant error should be times smaller than the order of dominant error and consequently, the error can be negligible by choosing an optimized . It can be achieved by applying (17) to equation (16): (18)
AC
CE
Therefore, if is equal to , the maximum value of dominant error can be controlled by . After specifying the , can be calculated from by solving equation (16). This equation has a fifthorder derivative of and its solution is usually time-consuming in two and three dimensions. However, subject to the local behavior of Surfer part, this equation also can be solved locally to reduce the calculation time. In this case, the equation (16) can be integrated to reduce the order of derivatives as below: (19)
Using central difference method, equation (19) can be discretized at time step
as follows: (20)
Equation (20) is a pentadiagonal system of equations. system.
8
is obtained at time step
by solving this
ACCEPTED MANUSCRIPT 3.2.
Surfer part
Surfer part contains all the error-maker operators. In this part, the operators can be treated as a discontinuity, where the traditional discretization methods that based on differentiable variables use a certain amount of artificial viscosity to control error propagation and guarantee the convergence. In addition, to decrease the solution time similar to equation (16), the Surfer equation can be solved locally due to the local nature of this part. Equation (12) can be written as: (
)
(21)
CR IP T
The right-hand side of equation (21) has an acceleration role in the solution process. The central difference scheme is used for discretization: (22)
3.3.
AN US
The Suring scheme was explained using the Burgers’ equation as a simple case, at the present and previous sections. The scheme is applied to the one-dimensional Euler equations (system of equations containing mass, momentum, and energy equations) as a more complex case, in the following. It also should be noted that the surfing scheme can be applied to two- and three-dimensional system of equations. Euler Equations
The one-dimensional form of Euler equations are as follows:
M
(
)
(23-c)
are density, velocity, and total energy, respectively. Also,
PT
Where , , and obtained from:
(23-b)
ED
(
)
(
(23-a)
)(
)
is pressure and it is
(24)
AC
CE
Where is the ratio of specific heats of fluid. Equations (23), which span the entire flow domain , -, require two boundary conditions at the inlet and one at the outlet for a subsonic inflow and outflow [35] as below:
Where , obtained from
( ) (25-a) ( ) (25-b) ( ) (25-c) , and are inlet density, inlet pressure, and outlet pressure, respectively. is at the boundaries and by using equation (24) as follows: (26)
The initial conditions for equations (23) are as: ( (
) )
(27-a) (27-b) 9
ACCEPTED MANUSCRIPT ( ) (27-c) Where , , and are the initial values of density, velocity and pressure of the domain, respectively. In order to use the Surfing scheme, the flow variables are split into two parts, i.e. Surf and Surfer: (28-a) (28-b)
CR IP T
(28-c) (28-d) Where and are the Surf and Surfer parts of the relevant flow properties, respectively. To obtain the Surf and Surfer equations, equations (28) can be inserted into equations (23) similar to what was done in equation (9), but the mathematical calculation is complex. Since the Surf part approximates the flow variables and its equations should be the same form as Euler equations, it can be chosen similar to equations (23) as below: (29-a)
)
AN US
(
[(
)
(29-b)
]
(29-c)
]
(
(30-a) -
)
*
(
(30-b)
)
+
-
(30-c)
CE
,[
PT
,
ED
M
The initial values of the Surf variables (Eq. (29)) are obtained by applying the Sobolev gradient concept to the initial flow variables (similar to Eq. (16)). The Surfer equations are obtained by the difference between the Euler and Surf equations. Consequently, these equations are as follow:
Finally, the flow variables are calculated from (28) by the summation of the solution of equations (29) and (30).
AC
4. Results 4.1.
Burgers’ Equation
The first example is a discontinuity which is described by the functions (31) and (32). These functions are used as initial condition (2) and boundary conditions (3) for the Burgers’ equation to investigate the discontinuity propagation throughout the domain. (
)
{
(31)
( ) (32) As mentioned in the previous section, the Surf part approximates the base variable. This part at time step is obtained from equation (20). Table 1 shows the effects of various (at ) on the order of 10
ACCEPTED MANUSCRIPT error terms. As expected, these orders are reduced times compared to the error terms of ( achieved by ). Therefore, the errors can be controlled by . In this paper, the value of has been chosen equal to 10. Table 1: Effects of various (
)
(
(
4.2 0.42 0.054 0.0059
(
)
0.78 0.067 0.0076 0.00077
)
(
-0.15 -0.011 -0.0013 -0.00013
)
(
0.022 0.0016 0.00017 0.000017
)
0.0031 0.00020 0.000021 0.0000022
CR IP T
12 1.5 0.21 0.023
)
at on order of error terms.
AN US
To investigate the stability of a numerical method, the Courant–Friedrichs–Lewy (CFL) condition is ( ⁄ ) for one-dimensional case and is known as the Courant used. It is defined as number. Where is equal to one in this case, according to the initial condition (31). Fig. 4 (a) shows the Surf solutions at the first time step for various Courant numbers. As shown, the solution has never diverged for Courant numbers up to 10 10 and has a very high stability. The discretization error has not grown up. This is because of the dominant error is negligible compared to the main terms and has been controlled. Fig. 4 (b) shows Surfer solutions at first time step for various Courant numbers. These solutions modify the Surf part solutions for presenting the results as the final exact solution. As shown in this figure, the values of this part are decreased by increasing the Courant number because of negative acceleration on the right-hand side of equation (21). In addition, oscillations are introduced in the solution due to the presence of central difference dispersion error. (a)
1.2
0.1
s
0.6 0.4 0.2 0 -0.2
0
1
2
3
0.05 0
-0.05 -0.1 4
1
1.5
2
2.5
3
x
PT
x
ED
-1
C=10 0 C=10 C=102 C=104 6 C=10 8 C=10 C=1010
h
0.8
M
1
(b)
-1
C=10 0 C=10 C=102 4 C=10 6 C=10 C=108 C=1010
Fig. 4. (a) Surf and (b) Surfer solutions at first time step for various Courant numbers.
CE
To validate the Surfing scheme, the results are compared to the solutions of Lax, Lax-Wendroff, MacCormack, FCT, and first and second order TVD methods for Courant number of one in Fig. 5. As shown in these figures, the Surfing results are in a good agreement with the results of these methods. (a)
AC
1.2
(b)
1.2
0.8
0.8
0.6
0.6
u
1
u
1
0.4
0.4
0.2 0 -0.2
0.2
Lax Lax-Wendroff MacCormack Surfing
0
1
FCT (Beam and Warming) First-Order TVD Second-Order TVD Surfing
0 2
3
-0.2
4
x
0
1
2
3
4
x
Fig. 5. Comparison of Surfing solution at Courant number of one with: (a) Lax, Lax-Wendroff and MacCormack, (b) FCT, first and second order TVD solutions. 11
ACCEPTED MANUSCRIPT Fig. 6 shows the ultra-stability of the Surfing scheme compared to the above-expressed methods. The Courant number is 1000. As shown, all other methods have been diverged, while the present method has maintained its stability. This indicates that the Surfing scheme has superiority to existing methods. (a)
1.2
(b)
1.2
0.6
0.6
u
1 0.8
u
1 0.8
0.4
0.2 0 -0.2
0.2
Lax Lax-Wendroff MacCormack Surfing
0
1
First-Order TVD Second-Order TVD Surfing FCT (Beam and Warming)
0 2
3
-0.2
4
CR IP T
0.4
0
x
1
2
3
4
x
Fig. 6. Comparison of Surfing solution at Courant number of with: (a) Lax, Lax-Wendroff and MacCormack, (b) FCT, first and second order TVD solutions.
AN US
Fig. 7 (a) and (b) show solution of Burgers’ equation at five time steps for and Courant numbers with the presented method, respectively. As shown, the discontinuity moves forward with advancing in time. In addition, it moves with larger space step by increasing the Courant number (that equivalent to time step for a constant space step). The main solution has been tended asymptotically to the Surf part one, due to reducing the amount of the Surfer part. In addition, the ultra-stability of this scheme can be seen in these figures. (a)
1.2
1.2 1
0.8
0.4
Initial Step 1 Step 2 Step 3 Step 4 Step 5
0 -0.2
0
1
2
3
PT
0.2
0.6 0.4
ED
u
0.6
0.8
u
M
1
(b)
Initial Step 1 Step 2 Step 3 Step 4 Step 5
0.2 0 -0.2
4
x
0
1
2
3
4
x
CE
Fig. 7. Surfing solutions of Burgers’ equation at five time steps for: (a) Courant number of number of .
, (b) Courant
AC
In the previous results, was chosen as zero in order to the form of equation (11) and (1) can be the same and hence can approximate . In order to investigate the effects of the arbitrary function in ) *( + has chosen for a test case. Fig. 8 shows the equations (11) and (12), the function effects of various on the solutions for Courant number of at the first step. It illustrates that can control errors.
12
ACCEPTED MANUSCRIPT
(a)
(b)
1.2 1
1 0.8
0.6
u
u
0.8
0.6
0.4
0.4
0.2
0.2
0
0
-0.2
0
1
2
3
1.2
-0.2 1.8
4
1.9
2
x
2.2
2.3
2.4
x
on results of the Surfing scheme for the Courant number of magnification of the effects.
in the first step, (b)
CR IP T
Fig. 8. (a) Effects of
2.1
The second example is described by the functions (33) and (34). These functions are used as initial condition (2) and boundary conditions (3) for the Burgers’ equation. , ⁄ ( ) { (33) (
)
(a)
1.2
1.2
M
Initial Step 1 Step 2 Step 3 Step 4 Step 5
1 0.8
0
2
PT
0.2
4
6
Initial Step 1 Step 2 Step 3 Step 4 Step 5
1 0.8
u
u 0.4
(b)
0.6
ED
0.6
-0.2 0
AN US
(34) Fig. 9 (a) and (b) show solution of Burgers’ equation at five time steps for and of Courant numbers with the presented method for initial and boundary conditions of (33) and (34), respectively. As shown, a discontinuity is formed and moves forward with advancing in time. In addition, it is formed faster when the Courant number is increased. The ultra-stability of this scheme can be seen in these figures.
0.4 0.2 0 -0.2 0
8
2
x
4
6
8
x
4.2.
CE
Fig. 9. Surfing solutions of Burgers’ equation with initial condition of (33) at five time steps for: (a) Courant number of , (b) Courant number of .
Euler Equations
AC
The boundary conditions (25) and initial conditions (27), which are used for solving the equations (23), are as follow: (
)
(
)
( ⁄(
(
)
(
( (
)
(
)) )
)
(35)
)
(36-a)
(36-b) ( ) (36-c) Where , is equal to 1.4 for air. Fig. 10 (a) shows the geometry of a nozzle. The nozzle throat is located at ⁄ , where L is the length of the nozzle. As shown in Fig. 10 (b) through (d), there is choking at the nozzle throat, after that, a normal shock wave has been formed at the divergent part of 13
ACCEPTED MANUSCRIPT the nozzle. These figures show that the Surfing scheme has a good agreement with SCDS and SLAU2 [36] results. (a)
0.5
(b) SCDS SLAU2 Surfing
0.4
SCDS SLAU2 Surfing
2
Area
1.5
u
0.3
1
0.2
0
0.2
0.4
0.6
0.8
0 0
1
0.2
x
0.6
0.8
1
(d)
1.2
SCDS SLAU2 Surfing
1
0.4
x
(c)
1.2
CR IP T
0.5 0.1
SCDS SLAU2 Surfing
1 0.8
0.8
0.6
AN US
p
0.6
0.4
0.4
0.2
0.2 0 0
0
0.2
0.4
0.6
0.8
1
x
-0.2 0
0.2
0.4
0.6
0.8
1
x
M
Fig. 10. Solution of one-dimensional Euler equations: (a) geometry of the nozzle, (b) velocity distribution, (c) density distribution, (d) pressure distribution.
Table 2: Comparison of the computational times. Methods Computational Time (millisecond) Surfing 528 SCDS 2903 SLAU2 2716
AC
CE
PT
ED
Table 2 shows the comparison of computational times of the Surfing, SCDS, and SLAU2 methods. All results have been obtained by single-thread processing by using a 64-bit operating system with 2.5 GHz processor speed and 8 GB of RAM. It seems that the computational time must be increased due to the additional calculations of the Surfer solution and smoothing procedure. Since the Surf equation is solved with a very large time step compared to SCDS and SLAU2, convergence time is decreased for this part. As shown, the computation time of the presented method has been decreased to one-fifth.
5. Conclusion
The results of Surfing scheme were presented in the previous section. As mentioned, the stability condition of Burgers’ equation has been investigated up to 1010 and the results marvelously show the ultra-high stability of the presented method. According to this, the solution process could use large time steps in comparison with the traditional methods. This increased the convergence speed in the solution of equations. This was evident in solving the Euler equations that the solution time has been decreased to one-fifth. Certainly, the advantages and convergence speed of presented scheme in twoand three-dimensional cases become more obvious compared to one-dimensional case. The advantages of this method are listed as follow: 14
ACCEPTED MANUSCRIPT -
The Surfing scheme has ultra-stability. This scheme has rapid convergence. This scheme captures shock waves without any extra terms like artificial viscosity. This scheme can capture shock waves with high resolution by increasing the number of nodes used for solving the Surfer part, similar to adaptive meshes. Although the number of equations is increased, the calculation time is decreased due to increasing the time step.
6. Acknowledgement
CR IP T
The authors would like to thank Professor Earl Dowell (Department of Mechanical Engineering and Materials Science, Duke University) for his invaluable suggestions and fruitful comments which greatly improved the manuscript. It is gratefully acknowledged.
References
AC
CE
PT
ED
M
AN US
[1] K.A. Hoffmann and S.T. Chiang, “Computational Fluid Dynamics Volume I”, Fourth Edition, 2000. [2] H. Liu, “Shock capturing methods”, Advanced Topics in Numerical Methods, 2008. [3] T. Ohwada, R. Adachi, K. Xu, and J. Luo, “On the Remedy against Shock Anomalies in Kinetic Schemes”, Journal of Computational Physics, 106–129, 2013. [4] C. Chen, “High Order Shock Capturing Schemes for Hyperbolic Conservation Laws and the Application in Open Chanel Flows”, Ph.D. Dissertation, University of Kentucky, 2006. [5] A. Harten, S. Osher, B. Engquist, and S.R. Chakravarthy, “Some Results on Uniformly High-Order Accurate Essentially Nonoscillatory Schemes.”, Applied Numerical Mathematics, 347–376, 1986. [6] J.Y. Yang, “Second-Order Accurate Essentially Non-Oscillatory Schemes for the Euler Equations”, AIAA Journal, Vol.28, No.12, 2069–2076, 1990. [7] J.Y. Yang, “Third-Order Nonoscillatory Schemes for the Euler Equations”, AIAA Journal, 1611–1618, 1991. [8] X.D. Liu, S. Osher, and T. Chan, “Weighted Essentially Non-Oscillatory Schemes”, Journal of Computational Physics, 200–212, 1994. [9] A. Jameson, “Analysis and Design of Numerical Schemes for Gas Dynamics 1 Artificial Diffusion, Upwind Biasing, Limiters and Their Effect on Accuracy and Multigrid Convergence”, International Journal of Computational Fluid Dynamics, 171–218, 1995. [10] A. Jameson, “Analysis and Design of Numerical Schemes for Gas Dynamics 2 Artificial Diffusion and Discrete Shock Structure”, International Journal of Computational Fluid Dynamics, 1–38, 1996. [11] N.A. Adams and K. Shariff, “A High-Resolution Hybrid Compact-ENO Scheme for Shock-Turbulence Interaction Problems”, Journal of Computational Physics, 27–51, 1996. [12] Z.P. Wang, “An Essentially Non-Oscillatory High Order Padé-Type (ENO-Padé) Scheme”, Ph.D. Dissertation, University of Kentucky, 2002. [13] X. Deng and H. Maekawa, “Compact High-Order Accurate Nonlinear Schemes”, Journal of Computational Physics, 77–91, 1997. [14] X. Deng and H. Zhang, “Developing High-Order Weighted Compact Nonlinear Schemes”, Journal of Computational Physics, 22–44, 2000. [15] L. Jiang, H. Shan, and C. Liu, “Weighted Compact Schemes for Shock Capturing”, International Journal of Computational Fluid Dynamics, 147–155, 2001. [16] S. Pirozzoli, “Conservative Hybrid Compact-WENO Schemes for Shock Turbulence Interaction”, Journal of Computational Physics, 81–117, 2002. [17] L.D.G. Sigalotti, H. Lopez, A. Donoso, E. Sira, and J. Klapp, “A Shock-Capturing SPH Scheme Based on Adaptive Kernel Estimation”, Journal of Computational Physics, 124–149, 2006. [18] S. Kawai and S.K. Lele, “Localized Artificial Diffusivity Scheme for Discontinuity Capturing on Curvilinear Meshes”, Journal of Computational Physics, 9498–9526, 2008. [19] R. Kumar and M.K. Kadalbajoo, “A Class of High Resolution Shock Capturing Schemes for Hyperbolic Conservation Laws”, Applied Mathematics and Computation, 110–126, 2008.
15
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
[20] M.H. Djavareshkian and M.H.A. Jahdi, “Shock-Capturing Method using Characteristic-Based Dissipation Filters in Pressure-Based Algorithm”, Acta Mechanica, 2010. [21] A. Klöckner, T. Warburton, and J.S. Hesthaven, “Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method”, Mathematical Modelling of Natural Phenomena, 57–83, 2011. [22] R. Jiwari, “A Haar wavelet quasilinearization approach for numerical simulation of Burgers’ equation”, Computer Physics Communications, 2413–2423, 2012. [23] H.C. Yee, D.V. Kotov, W. Wang, and C.W. Shu, “Spurious Behavior of Shock-Capturing Methods: Problems Containing Stiff Source Terms and Discontinuities”, International Conference on Computational Fluid Dynamics, 2012. [24] J. Yu and C. Yan, “An Artificial Diffusivity Discontinuous Galerkin Scheme for Discontinuous Flows”, Computers & Fluids, 56–71, 2013. [25] M. Kumar and S. Pandit, “A composite numerical scheme for the numerical simulation of coupled Burgers’ equation”, Computer Physics Communications, 809–817, 2014. [26] H. Fu, Z. Wang, Y. Yan, and C. Liu, “Modified Weighted Compact Scheme with Global Weights for Shock Capturing”, Computers & Fluids, 165–176, 2014. [27] S. Premasuthan, C. Liang, and A. Jameson, “Computation of Flows with Shocks Using the Spectral Difference Method with Artificial Viscosity, I: Basic Formulation and Application”, Computers & Fluids, 111–121, 2014. [28] S. Premasuthan, C. Liang, and A. Jameson, “Computation of Flows with Shocks Using the Spectral Difference Method with Artificial Viscosity, II: Modified Formulation with Local Mesh Refinement”, Computers & Fluids, 122–133, 2014. [29] R. Jiwari, “A hybrid numerical scheme for the numerical solution of the Burgers’ equation”, Computer Physics Communications, 59–67, 2015. [30] S. Pandit, M. Kumar, R.N. Mohapatra and A.S. Alshomrani, “Shock waves analysis of planar and non planar nonlinear Burgers’ equation using Scale-2 Haar wavelets”, International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 27 Issue: 8, 1814–1850, 2016. [31] D.V. Kotov, H.C. Yee, A.A. Wray, B. Sjögreen, and A.G. Kritsuk, “Numerical Dissipation Control in High Order Shock-Capturing Schemes for LES of Low Speed Flows”, Journal of Computational Physics, 189– 202, 2016. [32] J. Yang, B. Zhang, C. Liang, and Y. Rong, “A High-Order Flux Reconstruction Method with Adaptive Mesh Refinement and Artificial Diffusivity on Unstructured Moving/Deforming Mesh for Shock Capturing”, Computers & Fluids, 17–35, 2016. [33] A. Mazaheri and H. Nishikawa, “High-Order Shock-Capturing Hyperbolic Residual-Distribution Schemes on Irregular Triangular Grids”, Computers & Fluids, 29–44, 2016. [34] R.C. Mittal and S. Pandit, “Sensitivity Analysis of Shock Wave Burgers’ Equation via a Novel Algorithm Based on scale-3 Haar Wavelets”, International Journal of Computer Mathematics, 601–625, 2017. [35] C.A.J. Fletcher, “Computational Techniques for Fluid Dynamics 2”, Springer Series in Computational Physics, 1988, Page 175. [36] K. Kitamura and E. Shima, “Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for AUSM-family schemes”, Journal of Computational Physics, 62–83, 2013.
16