An optimization technique in analyzing the Burgers equation

Abstract

Keywords: Finite element method; α-family of approximation; advection-diffusion process; Burgers equation; optimization; B-spline.

1. Introduction

Many physical processes encountered in physical environment are represented by differential equations. Most of the processes are used to model physical flows in various fields of sciences such as wave propagation, convection–diffusion processes, biological waves etc. One of those physical models is Burgers equation attracting much attention in dealing with evolution equations constituting different models [1]. Computation of the Burgers equation is an important first step towards developing methods. Under certain conditions, uniqueness and existence of solutions to the Burgers equation were shown and discussed [2]. The Burgers equation is the nonlinear model equation for diffusive waves in fluid dynamics. The corresponding equation has also many applied areas including theory of shock waves, sound waves in a viscous medium, mathematical modeling of turbulent fluid and so on. Much effort has been spent in solving the Burgers equation for last couple of decades. Since some exact solutions fail for small kinematic viscosity values [3], 𝜀 < 0.01, many authors [1,4– 20] have suggested various methods such as finite element, finite difference, boundary element in

Corresponding Author/Sorumlu Yazar: e-mail/e-ileti: sarim@yildiz.edu.tr, tel: +90 (212) 383 43 60

M. Sari, H. Tunc / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

numerically analysing the processes represented by the model problem. The Burgers equation were also solved exactly by using Hopf-Cole transformation [21,22]. Model of advection mechanisms, and diffusion transports leading to the Burgers equation is the nonlinear model equation for diffusive waves in fluid dynamics. Yet, the corresponding equation has many applied areas including some wave mechanisms in a viscous medium, modelling of turbulent fluid etc. This study proposes to establish a Galerkin type finite element method (FEM) in which a strong form of the considered equation is preferred rather than the weak form. Because, the use of the strong form of the FEM in analysing the advection-diffusion processes represented by the Burgers equation has some superiorities comparison to the latter. The weak form of the equation needs more complicated computers codes, more computational time than the first one. Note that the weak form and strong form are mathematically equivalent to each other but computationally it is not the case. As the weak form of the equation requires additional matrices for residual term of the integration, this gives rise to excessive computational time and may therefore lead to loss of accuracy. It is noticeable that there are few studies [23,24] mainly used quadratic B-splines Galerkin method in solving the Burgers equation. In reference [24], they considered the fist-order splitted direct Crank-Nicolson method. The present article prefers to use the considered equation itself and general time approximation. This is the main difference between the present approach and their approach. In order to find the first approximate solution, an additional equation is required for both the present approach and for theirs. As will be seen in Section 4, the derivative of the initial condition is used here to find one additional condition whilst they used homogeneous condition. Some results have been given to explain the developed procedure in Section 6. Note also that the procedure followed here is fully different from Aksan’s work [23]. The proposed technique produces accurate and unconditionally stable results. The computational cost of the technique is acceptable. Even when the advection dominated cases of the Burgers equation are considered, the present approach produces oscillation-free results. The B-spline basis functions cover the spatial domain and for each of the elements of the spatial domain the well approximated solutions can be obtained. In addition to the spatial superiorities, the present time approximation is very suitable for the solution of the dynamical process of the considered equation. Thus, the 𝛼-family of time approximation is flexible and unconditionally stable for suitable choices of the parameter 𝛼. Here a hybrid numerical approach is proposed to solve the Burgers equation. To the best of the authors’ knowledge, this approach has not previously been suggested. The present method has been shown to be unconditionally stable for α ≥ 0.5. Behaviour of many processes arising in various fields of science leads to the Burgers equation problem is considered into the following form 𝑢𝑡 + 𝑢𝑢𝑥 = 𝜀𝑢𝑥𝑥 , 𝑎 ≤ 𝑥 ≤ 𝑏

with the boundary conditions 𝑢(𝑎, 𝑡) = 𝑓1 (𝑡), 𝑡 > 0, 𝑢(𝑏, 𝑡) = 𝑓2 (𝑡), 𝑡 > 0

where 𝜀 is viscosity constant for 𝜀 > 0 and 𝑓1 , 𝑓2 and 𝑔 are known functions. The subscripts 𝑥 and 𝑡 represent differentiations with respect to 𝑥 and 𝑡 respectively. The outline of this paper is as follows. The proposed technique is explained in Section 2. The implementation to the model equation is given in Section 3. The stability of the method is analysed in Section 4. Some numerical illustrations are presented in Section 5. Section 6 consists of some concluding remarks. This article proposes a numerical approach to solve the Burgers equation (1) with a set of boundary and initial conditions given by equations (2) and (3).

An Optimization Technique in Analyzing the Burgers … / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

3. Implementation

Time dependent quantities {𝛿}𝑙 must be found out using (14) and (15) to evaluate globally approximated solution (5). The following procedures are implemented to find {𝛿}𝑙 : (i) The first approximation {𝛿}0l is obtained from the initial condition (3). (𝑁 + 1) equations are found from the initial condition as well as one additional condition found from Table 1 using first derivative of the approximate function as follows 𝑢̃(𝑥𝑙 , 0) = 𝑔(𝑥𝑙 ) = {𝛿}0l−1 +{𝛿}0l 2 0 0 ′ (𝑥 ′ 𝑢̃𝑁 𝑙 , 0) = 𝑔 (𝑥𝑙 ) = ({𝛿}l − {𝛿}l−1 )

Thence the algebraic system will be (𝑁 + 2) × (𝑁 + 2) and is solved using the Thomas algorithm. (ii) To find second approximation {𝛿}1l , the above approximation {𝛿}0l is used at the righthand side of equation (14). The time dependent matrix 𝑅1∗ is required correction of {𝛿}0l . To ∗ approximate 𝑅1∗ first, {𝛿}1l = {𝛿}0l is chosen and then {𝛿}1l is found. Thus, for the rest of iterations, the coming combination is used about 10 times, to approximate {𝛿}1l . {𝛿}1l = 𝛼{𝛿}0l + (1 − 𝛼){𝛿}1l

(iii) For the other iterations, {𝛿}Il for 𝑖 = 2, … , 𝐽 (𝐽 is the chosen time element number), case (ii) is applied except that here the refinement is done only about 5 times.

4. Stability Analysis Of The Hybrid Method

Theorem: The proposed method for solving Equation (1) is unconditionally stable for α ≥ 0.5. Proof: To realize limitations of the computed solution under the consideration of the proposed hybrid approach (14), the stability analysis has been carried out with von Neumann theory taking Fourier growth factor defined by 𝛿sn = 𝛿̃ n 𝑒 𝑖𝑠𝑘ℎ (18) where 𝑘 is mode number and ℎ is the spatial element size selected for recursive approximation (14). To obtain a typical row of (14), values of δs+1 and δs in the time dependent ∗ matrices 𝑅𝑠+1 and 𝑅𝑠∗ are taken to be locally constant and equal to 𝑝. A typical row of pentadiagonal system (14) can thus be stated as n+1 n+1 + 𝑐 δn+1 + 𝑐 δn+1 𝑐1 δn+1 4 s+1 5 s+2 s−2 + 𝑐2 δs−1 + 𝑐3 δs n n = 𝑐6 δs−2 + 𝑐7 δs−1 + 𝑐8 δns + 𝑐9 δns+1 + 𝑐10 δns+2

M. Sari, H. Tunc / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

where 𝑐1 = 𝑟1 − 𝑟2 , 𝑐2 = 26𝑟1 − 2𝑟2, 𝑐3 = 66𝑟1 + 6𝑟2 𝑐4 = 26𝑟1 − 2𝑟2 , 𝑐5 = 𝑟1 − 𝑟2 , 𝑐6 = 𝑟1 + 𝑟3 𝑐7 = 26𝑟1 + 2𝑟3 , 𝑐8 = 66𝑟1 − 6𝑟3 , 𝑐9 = 26𝑟1 + 2𝑟3, 𝑐10 = 𝑟1 + 𝑟3 ℎ 2𝜀𝛼△t 2𝜀(1−𝛼)△t 𝑟1 = , 𝑟2 = , 𝑟3 = . 30

Substitution of (18) into (19) and use of Euler expansion for exponential terms give rise to 𝑔δ̃n+1 = 𝑔* δ̃n (21) where 𝑔 = (𝑟1 − 𝑟2 )cos(2𝑘ℎ) + (26𝑟1 − 2𝑟2 ) cos(𝑘ℎ) + 33𝑟1 + 3𝑟2 𝑔∗ = (𝑟1 + 𝑟3 )cos(2𝑘ℎ) + (26𝑟1 + 2𝑟3 ) cos(𝑘ℎ) + 33𝑟1 − 3𝑟3

Scheme (14) is stable if and only if |𝑧| < 1. Thus the following inequality must be fulfilled 𝑔∗

When 𝛼 ≥ 0.5, the inequality 𝑟3 ≤ 𝑟2 is automatically satisfied. In expression (22), since 𝑟1 is depends only parameter ℎ and 𝑟1 ≪ 1, the absolute value of the term | 𝑟1 − 𝑟2 | behaves like | 𝑟2 | and the term | 𝑟1 + 𝑟3 | behaves like | 𝑟3 |. Notice that when 𝛼 ≥ 0.5, |33𝑟1 + 3𝑟2 | > |33𝑟1 − 3𝑟3 |. Hence under the condition of being 𝛼 ≥ 0.5, |𝑔| > |𝑔∗ | and (24) is satisfied. Then the proposed hybrid approximation is unconditionally stable under the consideration of the aforementioned cases. The other selections of 𝛼 values lead to a conditionally stable approximation.

5. Numerical Experiments

To figure out the effect of Galerkin FEM in strong form and α-family of approximation over numerical solutions of the Burgers equation, let us consider the following two test problems. To produce accurate results by dealing with meaningfully different values of the kinematic viscosity constant 𝜀 is considered. Example 1 [1] Let us consider homogeneous problem with initial condition 𝑢(𝑥, 0) = 𝑔(𝑥) = sin 𝜋𝑥, 0 < 𝑥 < 1

and homogenous Dirichlet boundary conditions 𝑢(0, 𝑡) = 0 , 𝑡 > 0

A smooth exact solution of (1) under the consideration of the cases (25)-(27) given by Cole [22] is ∑∞

𝑛=1 𝑛 𝑢(𝑥, 𝑡) = 2𝜋𝜀 𝑎 +∑ ∞ 𝑎 exp(−𝑛2 𝜋 2 𝜀𝑡)𝑐𝑜𝑠(𝑛𝜋𝑥) 0

An Optimization Technique in Analyzing the Burgers … / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

𝑎𝑛 = 2 ∫0 𝑒𝑥𝑝{−(2𝜋𝜀)−1 [1 − cos(𝜋𝑥)]}cos(𝑛𝜋𝑥)𝑑𝑥 . In the present example, comparison of the produced results with the results of the literature [1,9,12,26-28] and exact solutions has been carried out in Tables 2-4, for various spatial points at both small and large times. As realized from the tables, the produced results here are more accurate and more economical, even with less number of time elements, comparison to the taken results from the literature. The computed results revealed that use of far less number of time elements for the proposed method is capable of catching better accuracy than the compared results [1,9,12,26-28]. Figure 1 shows the effect of parameter 𝛼 under the consideration of 𝐿2 error. Now it is time to deal with the smaller kinematic viscosity constants. Comparison of the currently produced solutions has been done with the literature and the exact solution for various values of the physical factors, with challenging kinematic viscosity values as seen in Tables 3-4. The computed results revealed that, even with the use of far less number of time elements, one can find similar or sometimes more accurate results than the literature [1,9,12]. The calculated solutions are depicted at various values of parameters 𝜀, 𝑑𝑡, ℎ at different times in Figures 2-7. The present method is observed to be very effective on capturing the steep behavior of the solution function as seen in Figures 2-3. Figures 8-9 show effects of the optimization factor on the computed results. As seen in Figure 8 under the fixed parameters 𝜀 = 1, ℎ = 0.025 and 𝑑𝑡 = 0.02, the optimal value of 𝛼 is equal to 0.525. Under the consideration of parameters 𝜀 = 1, ℎ = 0.01 and 𝑑𝑡 = 0.05, the optimal value of 𝛼 is seen to be 0.550 (see Figure 9). The present technique has also been compared with the work of Dag et al. [24] numerically and the current results are seen to be more accurate than theirs. The absolute errors are 1.6𝐸 − 06 and 0.6𝐸 − 05 (their Table 3), respectively, for the same physical parameters. Their result is taken from their Table 3. Note that behaviour of the solution of the Burgers equation for various values of viscosity constant 𝜀 were discussed in the literature, e.g. Dag et al. [24], Sari and Gurarslan [1]. Physical behaviour of the nonlinear advection-diffusion process calculated by the proposed technique in terms of the viscosity constant exhibits the expected physical characteristics of the problem. Comparison of the present method with the differential quadrature based methods [29-31] has also been done in Table 5. Even if their time approximation technique is the RK4 and conditionally stable, the present results are the same as with their results or more accurate than theirs as seen in Table 5.

M. Sari, H. Tunc / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Figure 1. Comparison of 𝐿2 error norms presented in Table 1

Figure 2. Numerical solution of the problem at different times produced for the parameters 𝜀 = 0.001, ℎ = 0.0016, 𝛼 = 0.50 and 𝑑𝑡 = 0.1

An Optimization Technique in Analyzing the Burgers … / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Figure 3. Numerical solution of the problem at different times produced for the parameters 𝜀 = 0.0005, ℎ = 0.0014, 𝛼 = 0.50 and 𝑑𝑡 = 0.1

Figure 4. Numerical solution of the problem produced for the parameters 𝜀 = 1, ℎ = 0.0125, 𝛼 = 0.50 and 𝑑𝑡 = 0.01 377

M. Sari, H. Tunc / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Figure 5. Numerical solution of the problem produced for the parameters 𝜀 = 0.1, ℎ = 0.0125, 𝛼 = 0.50 and 𝑑𝑡 = 0.01

Figure 6. Numerical solution of the problem produced by 𝜀 = 0.01, ℎ = 0.0125, 𝛼 = 0.50 and 𝑑𝑡 = 0.01

An Optimization Technique in Analyzing the Burgers … / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Figure 7. Numerical solution of the problem produced by 𝜀 = 0.001, ℎ = 0.0025, 𝛼 = 0.50 and 𝑑𝑡 = 0.1

Figure 8. Numerical and analytical solutions of the problem at 𝑡 = 0.5 produced by 𝜀 = 1, ℎ = 0.025 and 𝑑𝑡 = 0.02 where a) 𝛼 = 0.500 b) 𝛼 = 0.525 c) 𝛼 = 0.550 d) 𝛼 = 0.575 e) 𝛼 = 0.600 f) Exact.

M. Sari, H. Tunc / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Figure 9. Numerical and analytical solutions of the problem at 𝑡 = 0.5 produced by 𝜀 = 1, ℎ = 0.01 and 𝑑𝑡 = 0.05 where a) 𝛼 = 0.500 b) 𝛼 = 0.525 c) 𝛼 = 0.550 d) 𝛼 = 0.575 e) 𝛼 = 0.600 f) Exact Table 2. Comparison of the results produced with 𝜀 = 0.05, ℎ = 0.01 for various values of the parameter 𝛼 and different selection of 𝑑𝑡.

An Optimization Technique in Analyzing the Burgers … / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Table 3. Comparison of the results produced with 𝜀 = 0.003, ℎ = 0.01 for various values of parameters 𝛼 and 𝑑𝑡. 𝑥

Table 4. Comparison of maximum error norms of various schemes for 𝜀 = 0.01, ℎ = 0.0125. 𝑥

Eefdm

M. Sari, H. Tunc / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Table 5. Comparison of the produced results for 𝛼 = 0.50 and 𝜀 = 0.1. Present N=20

Example 2 [1] Let us take now then Burgers equation (1) with initial condition 𝑢(𝑥, 0) = 𝑔(𝑥) = 4𝑥(1 − 𝑥),

The exact solution of (1) under the consideration of the cases (22)-(24) given by Cole [21] as in (28) but with the Fourier coefficients 1

𝑎𝑛 = 2 ∫0 𝑒𝑥𝑝{−𝑥 2 (3𝜀)−1 (3 − 2𝑥)}cos(𝑛𝜋𝑥)𝑑𝑥. Table 6 includes comparison of currently produced solutions with exact solution and the literature [16,28] for different values of physical parameters under the consideration of various α values. The computed results for less number of time elements, comparison to some works carried out [16,28], are seen to be more accurate or of the same accuracy with the corresponding literatures. By taking into account various values of the processes of interest, all the presented numerical solutions in Table 7 are compared with exact solution and the literature [8,14]. The numerical behaviour at 𝑡 = 0.5 for challenging values of the kinematic viscosity is exhibited both quantitatively and qualitatively (see Table 6 and Figure 10). The present work has required less number of time elements for more accurate results as compared to the studies done in references [8,14]. As previously stated, the optimum values of the parameter α are affected from the number of the time elements although the corresponding optimum values are seen to be not seriously affected from the number of spatial elements (see Figure 11).

An Optimization Technique in Analyzing the Burgers … / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Figure 10. Numerical solution of Example 2 at 𝑡 = 0.5 with different kinematic viscosity constants: a) 𝜀 = 1 b) 𝜀 = 0.5 c) 𝜀 = 0.1 d) 𝜀 = 0.01 e) 𝜀 = 0.005 f) 𝜀 = 0.001

Figure 11. Numerical and analytical solutions of the Example 2 at 𝑡 = 0.5 produced by 𝜀 = 1, ℎ = 0.025 and 𝑑𝑡 = 0.02 where a) 𝛼 = 0.500 b) 𝛼 = 0.525 c) 𝛼 = 0.550 d) 𝛼 = 0.575 e) 𝛼 = 0.600 f) Exact 383

M. Sari, H. Tunc / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

Table 6. Comparison of the results produced with 𝜀 = 1 for various values of parameters 𝛼, 𝑑𝑡 and ℎ in Example 2

Present α=0.50 ℎ = 0.0125, 𝑑𝑡 = 0.0002 0.4262850 0.2614793 0.1614772 0.0610873

Present α=0.51 ℎ = 0.0125, 𝑑𝑡 = 0.0002 0.4262904 0.2614845 0.1614819 0.0610903

Table 7. Comparison of the results produced with 𝜀 = 0.1, ℎ = 0.0125 for various values of parameters 𝛼 and 𝑑𝑡 in Example 2 𝑥

5. Conclusions And Recommendation

This paper has proposed a hybrid numerical technique to deal with the Burgers equation. The designed technique has then been shown to be unconditionally stable for parameter α ≥ 0.5. The technique has been illustrated to have great potentiality in analysing the Burgers equation with less number of elements used in both time and space even for challenging cases of the problem. The computed results have also been realized to have the same level of accuracy or to be better

An Optimization Technique in Analyzing the Burgers … / Sigma J Eng & Nat Sci 35 (3), 369-386, 2017

than existing results in the literature, even with the use of a modest time approximation. Note that the suggested approach is seen to be a very good alternative to achieve a high degree of accuracy while analysing the advection-diffusion processes. In addition to both time and space discretization, the choice of parameter α is observed to be effective on capturing the behavior of the problem. The presented results have validated all cases of α-optimization. Advantages of the current algorithm are clearly seen especially from the steep behavior of the produced results. As is the case for all versions of the finite element based methods, storage drawbacks of the proposed methods are expected to come out in two/three dimensional cases of the current models for very large domains. Future studies can focus on designing of the current technique to physical processes represented by more involved time-dependent models.

Share and Cite