Numerical study of non-linear waves for one-dimensional planar, cylindrical and spherical flow using B-spline finite element method

: In a recent study, an evolution equation is found for waves' behavior at far-field with relaxation mode of molecules. An analytical technique was used to solve this evolution problem, which is a generalized Burger equation. The analytical approach has limitations and requires a very accurate initial guess by a trial method. In this paper, the evolution equation for one-dimensional planar, cylindrical, and spherical flow in the presence of relaxation mode is solved using a collocation approach with a cubic B-spline function. The numerical results are graphed and compared with the exact solution for planar flow. The obtained numerical results match the exact solution quite well and show that the technique is quite reliable and can deal with the nonlinearity involved in the present problem. Results have also been obtained for cylindrical and spherical flow at the far-field. The obtained numerical results show that the present approach with the cubic B-spline function works well and accurately. Fourier stability analysis is used to investigate the stability of the cubic B-spline collocation method.


Introduction
The Euler system of equations governs the conservation laws for compressible flows. It's still challenging to find an analytical solution for this highly non-linear system of equations. The case of planar flow is less complicated than that of cylindrical and spherical flow. It is due to the source term in geometric form when Euler's system is expressed in polar coordinates. It shows that there is a singularity for = 0 in the source term. The singularity at = 0 makes it challenging to find the approximate solutions of the cylindrically and spherical flow. Therefore, studying the behavior of converging flows is still a challenge. These types of flows are quite important for researchers due to their usefulness in detonations. Mainly, during World War II, many studies were performed to know the dynamics of these non-linear waves [1]. However, these studies have faced a serious challenge as the numerical solutions blow up at the axis of convergence for various flow variables. Finding the accurate numerical solution at the axis of convergence for imploding flows is still an open problem. Many groups, therefore, have worked to give approximate numerical solutions in the vicinity of the axis of convergence.
To make the problem simpler, researchers have first considered the stationary flows. This problem has been achieved by taking two concentric cylinders. The gas is injected in the bounded region from one end and extracted at the opposite end. The rate of injection and extraction of gas were made equal. It creates stationary flows. Chen and Glim [2] have studied spherical stationary shocks at the point of extraction. However, the problem of non-stationary flows is much more complex and still an open problem. Guderley [3] has found the solution of imploding shocks passing through a channel with a non-uniform cross-section. Whitham [4] has proposed geometrical shock dynamics to find an approximate analytical solution for symmetric flows. His theory was based on the characteristic equation entering the shock being used to approximate the shock front. Geometrical Shock Dynamics (GSD) is a bit simpler as it makes the problem non-singular. Many research groups have studied similarity solutions to imploding and exploding shocks. The method is useful as it converts the onedimensional flow represented by partial differential equations to ordinary differential equations. In the context of similarity solution, we would like to highlight the work of Siddiqui et al. [5], Zakeri [6], Sharma and Arora [7], Sedov [8], Taylor [9], and Sakurai [10,11].
Due to the above difficulties in finding approximate solutions, asymptotic solutions of the problem are quite useful. These solutions may also be used to verify the accuracy of the numerical code. To find the behavior of these wave propagation at the far-field is a problem that researchers are very interested in solving. It gives valuable information about flow variables at far places from a source. To deal with these types of problems, a number of numerical and analytical approaches have been developed [12][13][14][15][16][17][18][19][20][21][22]. Many researchers have studied this problem. Among them, the authors would like to highlight the work of Siddiqui et al. [23], Sharma and Jena [24], and Manickam & Radha [25]. These studies have presented an asymptotic equation governing the behavior of waves. A non-linear modified Burger's equation served as the asymptotic equation. The perturbation approach has been suggested by Sharma and Jena [24] to obtain the solution to the modified Burger's equation. The authors concluded that perturbation techniques do not produce a general solution and are only helpful in certain situations. In a recently published article by Siddiqui et al. [23], Homotopy Analysis Method [26,27] is used. This method is quite suitable for weakly non-linear problems but faces serious difficulty for highly non-linear problems. Since the evolution equation obtained by asymptotic analysis is highly non-linear and governs all characteristics of the system, some numerical technique should be used to provide the general solution.
In this study, modified Burger's equation for planar, cylindrical, and spherical flow in the presence of relaxation mode is solved using a cubic B-spline collocation approach. The proposed problem is solved as a Cauchy problem. Planar flows have an analytical solution to the evolution equation. We compared our numerical results to the analytical solutions for planar flow that is currently available. In addition, the flow profiles for the cylindrical and spherical flows have also been plotted. This approximate solution is challenging to find for symmetric cylindrical and spherical flow. By providing a more general numerical scheme that can be used with any initial condition, the authors have attempted to fill the gap in this article.
The following is the paper's structure. Section 1 introduces the background and importance of non-linear wave propagation in an ideal medium. Asymptotic analysis is used in Section 2 to describe the mathematical formulation of the governing equation. Section 3 describes the methodology used in the paper. The spatial and time discretization of the problem is discussed in this section. Section 4 presents the stability analysis of the method. The results and discussion are presented in Section 5. The figures are plotted for several values of the specific gas ratio. Finally, the study's conclusion is written in Section 6.

Governing evolution equation
In non-ideal and relaxing gas, the set of equations that govern planar, cylindrical, and spherical flow is given by [23][24][25].
The density, velocity, and pressure of the gas are denoted by , , and , respectively. and represent the vibrational energy and rate of change of , respectively. The time variable is , and the spatial variable is . The speed of sound is represented by = √ (1− 0 ) . denotes the ratio of specific heats under constant pressure and volume and 0 is van der Waals excluded volume. The following formula yields as: where ̅ = 0 + ( − 0 ). Relaxation time, gas constant, and temperature are denoted by , R and T, respectively. The ratio of vibrational specific heat to the specific gas constant is . The initial values of , and are 0 , 0 and 0 respectively. The state equation for a real gas is The fast variables = , = 2 , are significantly less than 1, and are introduced to explore non-linear wave behavior at a distance from the source. where = − . Following some algebra, the evolution equation is derived by applying asymptotic expansion of flow variables in a set of equations expressed in stretched variables [23] which govern the flow behavior of waves in the farfield.
where (1) denotes to first approximation to the pressure and where the gas pressure is denoted by . = , = 2 are the stretched variable of spatial variable and time-variable respectively. represents a small parameter. is the time taken by gas molecules from an exciting state to a relaxing state. γ denotes the heat capacity ratio of the gas. represents the ratio of specific heat.
is used for the particular geometry of the flow. Planar, cylindrical and spherical waves are represented by = 0, 1, 2, respectively.
Consider (1) = . The Eq (7) can be rewritten as: The initial and boundary conditions are: The Eq (10) will be solved as a Cauchy problem based on cubic B-spline functions to observe the far-field behavior of non-linear waves. The flow profiles' time evolution is shown.
For > 0, Using the de Boor recursion formula (15), the third degree ( = 3) B-splines termed as cubic Bsplines ,3 ( ), = 0, … , are defined as where Explicitly, cubic B-splines can be written as Because each cubic B-spline ,3 ( ) covers 4 intervals or elements, each finite interval [ , +1 ] is covered by 4 splines. Approximate solution ( , ) to the exact solution of Eq (10) with initial condition (11) and boundary conditions (12) and (13) can be defined as cubic B-spline function in the following form where ( ) represent the time-dependent coefficients. The unknown parameters ( ) will be determined using by requiring that satisfies Eq (4) at + 1 collocation points and boundary conditions.
The approximation ( , ) over typical element [ , +1 ] is given by According to the local support properties of the basis function, there are only three nonzero basis functions that will be included for the evaluation at each . Using m-th B-spline basis function (17) and B-spline function (18), we get the values of and its various order derivatives ′ , ′′ at the knots in terms of as follows: where dashes " ′, ′′ " represent the first and the second derivatives with respect to , respectively. To use the collocation approach, first, choose collocation points that correspond with knots, then substitute the values of and its derivatives from Eqs (20)- (22) in evolution equation (10). This gives the first-order ordinary differential equations system as: The dot (⋅) is used to show the derivative with respect to and = −1 + 4 + +1 is a non-linear term of Eq (23).
Derivative terms ̇ in Eq (23) is discretized using an upwind scheme as follows: Crank-Nicolson scheme is used for to get where upper index represents the time layer.
Values of ̇ and are substituted in Eq (23). The following non-linear system of equations is obtained.
In a simplified form, the above equation can be written as where 1 = 4ℎ 2 + ℎ 2 Δ − 6ℎ Δ − 12 Δ , 2 = 16ℎ 2 + 4 ℎ 2 Δ + 24 Δ , System (27) has ( + 1) equations in ( + 3) unknown. Two unknown parameters −1 +1 and +1 are eliminated by using boundary conditions (12 and 13) to solve this system of equations. The system (23) is solved using a Thomas algorithm [29]. The initial vectors 0 are required to start the iteration of the system (23), which can be determined using the initial and boundary conditions: where all parameters 0 are determined. must satisfy the following relations at points : A tridiagonal matrix system is produced using the condition mentioned above and the initial vectors 0 are obtained from the following system: . (31) The Thomas algorithm [29] is used to solve the tridiagonal system (31). The obtained initial vectors will serve as the starting point for our scheme in (27).
The following approach is used at each time steps for the element parameters to deal with the nonlinearity of the system (27).
Once we find the approximation n+1 utilizing the systems (27), the new approximation to n+1 is obtained using the iteration approach (32). Before moving to the next time step solution n+1 , the iteration should be performed two or three times to refine n+1 .

Fourier stability analysis
In this section, Fourier stability analysis is used to investigate the stability of the cubic B-spline collocation method as applied to a linear differential equation. For this, the non-linear term is considered as a linear term by taking as locally constant g in Eq (10). We obtain the following linear system The following linear system of equations is obtained by utilizing the collocation method, as mentioned in section 3.
The simplified form of Eq (34) is Substituting the Fourier mode = ℎ , into linearized Eq (35), where represents the wavenumber and ℎ is the element width.
We get: where 1 = (8ℎ 2 − 2 ℎ 2 Δ + 24 Δ ) cos( ℎ) + 2 , 2 = (8ℎ 2 + 2 ℎ 2 Δ − 24 Δ ) cos( ℎ) + 2 , The Von-Neumann stability condition for the above system is that the maximum modulus of the eigenvalues has to be less than or equal to one. Since 2 > 1 , we obtain the modulus | | < 1 . This means the proposed scheme (27) is unconditionally stable according to von Neumann's stability analysis. Previous studies indicated that there is only a marginal loss in accuracy when such linearization is adopted.

Results and discussion
This section presents the approximate results of a modified generalized Burgers' equation (10). The proposed method's accuracy is assessed using discrete ∞ error norms by: where, and represent the approximate solution and exact solution respectively. For the first test example, we consider the Shock-like solution of Eq (10) with the exact solution when = 0 is: where 0 = exp (1/8 ). The initial condition is: To test performance of the proposed scheme, the test example is solved with the given values ℎ = 0.005 and = 0.01, = 0.005. The simulation is run at various times. The numerical results of the current scheme and the exact solution at = 2.5 are compared in Table 1. The obtained numerical solutions reflect exact solutions within very acceptable limits. The proposed scheme's results are noticeably better than earlier results [30]. Figure 1 shows the graphical results of the present scheme at different times. Figure 1 also shows a comparison between numerical and exact solutions.
The numerical solutions for viscosity coefficients = 0.0005 at various times are graphed in Figure 2. It is noticed that the propagation front is steeper at lower viscosity values.  Second, test problem, the exact solution of Eq (10) for planar flow when =0 is given by where This exact solution is obtained by using Cole Hopf transformation. The proposed scheme's results have been calculated for various values of spatial variables and time as tabulated in Table 2. The maximum error norms are obtained relatively very small, and the flow profiles of exact and numerical solutions are almost identical.   It is noticed that even with the smooth initial profile, the problem develops a discontinuity as time progresses for planar flow. This discontinuity converts into a shock at some critical time, as seen in Figure 4. Moreover, further pressure increases very rapidly after the initial phase with the progress in time. It, therefore, suggests that the shock that develops with progress in time for planar flow. To investigate the cylindrical flow at the far-field, we take = 1 in equation (10). Flow profiles for different times are obtained by the scheme, as displayed in Figure 5. Again, the initial smooth profile culminates into a shock after a particular time. Pressure (W) increases rapidly after the initial phase with the increment in time, and a shock wave develops after a critical time value. However, the formation of shock waves is faster in planar flow than in cylindrical flow. The behavior of spherical waves in the far-field is studied, taking = 2 in the evolution equation (10) that has been obtained from the conservation laws. Numerical solutions for different time steps are calculated using the cubic B-spline collocation method. Figure 6 displays that even a smooth initial profile develops shocks quickly after some initial phase. Further, the pressure increases with time, suggesting that the developed shock will strengthen with time. Numerical solutions of the evolution equation (10) are computed at different periods in the flow profiles to analyze the various possible geometry of the flow. The values of the other physical parameters are chosen. It is evident from Figure 7 that there is not any noticeable difference in the flow profiles in the initial phase. However, the flow profiles are influenced by the geometry after the initial phase, as shown in Figure 7. The pressure decreases faster with the increment in . It was also noticed that the pressure is more for planar flow than cylindrical flow.
Similarly, it is noted from Figure 7 that the cylindrical flow pressure is more than the pressure for the spherical flow. However, any flow profile corresponding to any geometry develops a shock after some critical value. Planar flow produces shock more quickly than cylindrical and spherical flow. Similarly, shock arises faster in cylindrical flow than in spherical flow.

Conclusions
This paper presents an efficient technique, namely the cubic B-spline collocation method as an alternative to existing methods in partial differential equations. It has been used to solve the modified Burger's equation for planar, cylindrical, and spherical flow in the presence of relaxation with great success. Examples demonstrate that the proposed technique is a useful and powerful mathematical tool for dealing with non-linear problems of this nature. A cubic B-spline collocation method is presented in this article to study the non-linear wave propagation at the far-field. This has been incorporated in the problem by the introduction of fast variable = , = / 2 , where = 1. The flow profiles for planar ( = 0), cylindrical ( = 1), and spherical ( = 2) geometries of non-linear waves were determined using the proposed scheme. The proposed scheme's results have been calculated for various values of spatial variables and times. The maximum error norms are obtained relatively very small. It is observed that the flow culminated into a shock as time progresses and the formation of the shock is faster in planar flow than in corresponding cylindrical or spherical flow. A good choice of collocation points for the spline difference operator on the piecewise uniform mesh ensures the discrete maximum principle. As a result, the scheme is stable in the maximum-norm case. The simplicity of the cubic B-spline collocation method is an advantage over the other numerical methods to obtain the numerical solution of non-linear differential equations. Due to the presence of a B-spline function in the proposed scheme, the desired accuracy of the results is obtained by taking a cubic B-spline function along with a collocation approach. Due to the properties of a uniform Bspline basis, this method can only be used on a structured grid. As a result, this method is incapable of dealing with arbitrary geometry or irregular node distribution. This necessitates a new line of research that incorporates the non-uniform B-spline rather than the uniform B-spline. As a result, the proposed cubic B-spline collocation method can be used to approximate two-dimensional nonlinear differential equations.