An Element-Driven Boundary Integral Treatment for Nonlinearity: A Coupled Nonlinear Two-Dimensional Burger’s Equation Test Case
Okey Oseloka Onyejekwe
Computational Science Program Addis Ababa University Arat Kilo Campus, Addis Ababa, Ethiopia
Email address:
To cite this article:
Okey Oseloka Onyejekwe. An Element-Driven Boundary Integral Treatment for Nonlinearity: A Coupled Nonlinear Two-Dimensional Burger’s Equation Test Case. International Journal of Fluid Mechanics & Thermal Sciences. Vol. 3, No. 1, 2017, pp. 1-15. doi: 10.11648/j.ijfmts.20170301.11
Received: December 11, 2016; Accepted: December 26, 2016; Published: January 24, 2017
Abstract: To further improve on the competitiveness of the boundary element method (BEM), a hybrid version of it is used for a numerical solution of two dimensional nonlinear coupled viscous Burger’s equation. Adopting this approach to a discretized 2D spatial domain, the resulting integral equations arising from the singular integral theory are applied locally to each of the elements. The resulting nonlinear discrete equations are finally solved by the Picard iteration algorithm. The simulation results obtained, not only concur with analytical solutions, but also display high accuracy and are in agreement with those available in literature.
Keywords: Two Dimensional, Coupled Nonlinear Burger’S Equation, Hybrid Boundary Element Method, Integral Equation, Singular Integral Theory, Discretization
1. Introduction
The two dimensional coupled nonlinear Burger’s equation forms an important class of a system of nonlinear parabolic and hyperbolic partial differential equations (PDEs) that models unsteady flow equations. It bears a close resemblance to the Navier-Stokes equations by incorporating the same convective and diffusive terms but differs by the absence of the pressure term. It is often used to study various aspects of physical phenomena involving among others turbulence, flow through shockwaves [1], sedimentation of particles in fluid suspensions or colloids under the effect of gravity [2], wave processes in thermo-elastic medium [3], dispersion in porous media, and vorticity transport.
The analytic solution of the coupled nonlinear Burger’s equation involving various transformation techniques are given by Fletcher [4], Abazari and Borhanifar [5] and Ablowitz and Clarkson [6]. Other powerful methods of arriving at exact solutions include the homotopy perturbation method, Hirota’s bilinear method (Hirota [7]), and the delta method (Bender et al. [8]). Kaya [9], Zhu et. Al. [10] solved the system of coupled Burger’s equation by the Adomian decomposition method. Comparison between exact nonlinear equation solutions and those obtained numerically were carried out by Fletcher [11], Abdou and Soliman [12]. A combination of the homotopy perturbation and Pade’ approximation were applied by Kelleci and Yildrim [13] to study both the homogeneous and inhomogeneous coupled Burger’s equation.
The finite element method (FEM) in combination with the spline-interpolation techniques have been used to solve the coupled nonlinear Burger’s equation (Kutluay, and Ucar [14]). Higher order accurate schemes have also been used (Esipov [3]). This is accompanied by those involving variational iteration method (Soliman 2009 [15], Abdou and Soliman [12]), Odibat and Momani [16]). Bahadir [17] applied a finite –difference fully implicit technique to the discretization of coupled nonlinear Burger’s equation and resolved the nonlinearity with the Newton’s method. Other finite difference applications can be found in Srivastava et al. [18], Srivastava et al. [19]. These followed earlier work by Wubs and Goede [20], and Goyon [21].
BEM literature reveals a paucity of information regarding the application of boundary element methods to the solution of coupled nonlinear partial differential equations especially for 2D coupled nonlinear Burger’s equation Onyejekwe [22]. It is well known that considerable numerical difficulties arise in extending BEM application to inhomogeneous, non-linear, transient, coupled nonlinear problems Grigoriev [23], Siraj-ul-islam et al. [24], Toutip [25]. The major drawback in these cases is the need to discretize the problem domain into an aggregation of cells or elements needed to deal with domain integrations. It will not be possible to go into details here. A comprehensive literature on this topic can be found in Grigoriev’s paper [23]. However it is noted that the current research trajectory of BEM still reveals an unresolved issue namely the drive to restrict or situate all domain integrals on the boundaries of the problem domain even if the physics of the problem dictates otherwise. The compelling need to make BEM a purely boundary-driven numerical technique motivates this trend. Unfortunately this is often accompanied by some numerical challenges some of which have been mentioned above and in some previous papers (Percher et al. [26], Lorinczi [27], Lorinczi et al. [28], Archer et al. [29], Archer [30], Onyejekwe [31, 32], Onyejekwe and Onyejekwe [33] Archer and Horne [34]).
Numerical experience shows that direct application of BEM theory performs poorly in the absence of a strong link between the problem domain and the boundary as is found in the Laplace equation. Hence there is always a need to modify BEM application for those problems that unambivalently incorporate the problem domain into problem formulation such as is found with body-force terms, nonlinearity, transience, heterogeneity, source and sink terms etc. Current efforts to prevail over this challenge have resulted in a variety of hybrid BEM techniques Grigoriev [23], Taigbenu [35], Onyejekwe [36], Onyejekwe [37], Taigbenu and Onyejekwe [38], Hibersek and Skerget [39], Onyejekwe [40], Grigoriev and Dargush [41], Perata and Popov [42], Portapilla and Power [43], Sladeck et al. [44], Toutip [25].
Hybrid BEM formulations, as promising simulation techniques, can either be regarded as extensions of the direct BEM approach or in most cases as specialized discrete forms of the traditional formulation that weigh more on the effective handling of the integral equations in a domain-localized sense. It ought to be mentioned as well that one of the reasons that prompted this option is the unavailability of auxiliary equations and/or fundamental solutions for many challenging and realistic problems. Even when they exist, they do so in forms that are not easily computable. More often than not, this leads to intricate formulations and systems of equations with fully populated matrix. Hence aiming at hybrid BEM formulations that mimic other domain-based numerical methods in terms of localized discretizations, and at the same time incorporate some of their unique advantages should motivate hybrid BEM research. This approach brings a lot of attractive features such as formulation simplicity, slender sparse coefficient matrices and an efficient handling of domain discretization. Hybrid BEM techniques have not only displayed competitive potential in modelling complex engineering problems (Abashar [45], Pecher et al. [26], Onyejekwe [46-48], Lorinczi et al. [27, 28] Lorinczi et al. (2010), Nyirenda [49], Bagherinezhad and Pishvaie [50], Archer and Horne [51, 52] but also have displayed a remarkable ability to accurately represent coupled nonlinear mathematical equations (Ramsuroop [53]).
2. Numerical Formulation
Consider the two-dimensional non-linear Burger’s equation:
(1a)
(1b)
where is the 2-D gradient operator with spatial variables in x and y, t is the time dimension, . Re is the Reynolds number and physically represents the ratio of inertia to viscous forces, are velocity components in x and y directions to be determined. It should be mentioned that in convection dominated flows, there is a steepening effect of the nonlinear advection term, and this effect creates noticeable numerical challenges. Equations (1a) and (1b) admit initial and boundary conditions. For a domain with a boundary , the system is subject to the following initial conditions
(2a)
(2b)
And boundary conditions:
(3a)
(3b)
where and are known functions.
3. Hybrid BEM Procedure
Details of this procedure have previously been dealt with elsewhere (Onyejekwe [54-56], Lorinczi [27, 28], Archer [30], Ramsuroop [53], Onyejekwe [54-56]); so only the most relevant will be mentioned. By Applying the Green’s theorem, the integral analogs of equations (1a) and (1b) are given as:
(4a)
(4b)
or
(5a)
where
and
(5b)
where
where is the fundamental solution of in an infinite domain, is the Euclidean distance between the field point and the source node at i, is the shape function of the nodal angle at . Although equations (4) and (5) were arrived at by exploiting the attractive features of BEM discretization, their implementation vis-à-vis the problem domain marks the point of departure between the boundary-only implementation of BEM and its hybridized form that does not consider domain integration as a disadvantage. To this end, suitable polygonal elements on which the discretized equations are numerically solved are employed to discretize the problem domain. A domain-driven-elemental approach exemplifies the limiting case of the domain decomposition technique whereby the subdomain is reduced to its barest minimum or finite elements (Taigbenu [35], Onyejekwe [36]). The key factors which are paramount at this stage of the solution process are identical to those of the finite element method (FEM) namely (a) linearization of the governing equations (b) interpolation of the functional quantities and continuity across element boundaries (c) choice of elements. For this study, we have elected to represent the distribution of scalar quantities by linear interpolation functions, that is where are the linear interpolation functions. Introducing this interpolation relationship into equations (5a) and (5b) yields the following discrete element equations which are applicable to each element.
(6a)
where
and
(6b)
where
The coefficients are given by:
All integrations can be done analytically for each of the rectangular elements. Without any loss in generality equations (6a) and (6b) can be given generically as:
(7)
The temporal derivative in equation (7) is approximated by a finite difference scheme with a weighting factor to give
(8)
where n is the iteration number, represent the current and previous times. With the introduction of known initial and boundary conditions, equation (8) can further be simplified into its coupled components to read:
(9)
(10)
The Picard algorithm is adopted for iterating the two nonlinear coupled equations and is started by:
a Storing the known scalar values obtained from prescribed initial and boundary conditions in the right hand side vectors of equations (9) and (10)
b Solving matrix equation (9) with the known values to obtain updated values of the scalar field
c These latest iterates from (b) are then used to obtain the scalar field from equation (10) i.e.
d The mean deviations from the latest iterates from the consecutive values of the dependent variables in (b) and (c) are compared and the iteration is stopped if the mean deviations are less than or equal to an a priori specified error tolerance ?. If this is not the case the procedure from steps (b) to (c) are repeated with refined iterates until convergence if finally achieved.
e When convergence is finally achieved, a time increment is implemented to obtain a new set of converged scalar field for each time level.
4. Numerical Examples and Discussion
Some test problems have been chosen to validate the formulation described herein. The presence of analytical results will help to quantify the relative performance of this scheme when compared with other solution methods available in literature. For each of the tested problems, the initial and boundary conditions are determined by their analytical solutions.
The accuracy and consistence of the scheme is measure in terms of errors defined by:
(11a)
(11b)
where represent exact and computed solutions. Equations (13a) and (13b) will also apply to the v dependent variable.
4.1. Problem 1
The exact solutions of equations (1) and (2) can be generated as shown in Fletcher (1983a):
(12)
The computational domain is a square , a uniform grid is applied with a mesh dimension of with the following problem parameters . Comparison of numerical and exact solutions shows that the proposed formulation achieves very close results for different values of Reynolds number. Figs. 1a and 1b are very close to those obtained analytically (Tamsir and Srivasta [57]). Tables 1a,1b, 1c and 1d display the absolute errors resulting from a comparison of numerical and analytical results for velocity profiles at different values of Reynolds number. A revealing information is the magnitude of absolute errors generated in representing areas of steep slope in the u and v velocity profiles. Figs. 1c and 1d display their relative magnitudes relative to those obtained in other areas of the profiles.
(x, y) | Numerical | Exact | Abs. Error |
(0.100000, 0.100000) | 0.534028 | 0.543322 | 0.009294 |
(0.500000, 0.100000) | 0.500270 | 0.500353 | 0.000083 |
(0.900000, 0.100000) | 0.500002 | 0.500002 | 0.000000 |
(0.300000, 0.300000) | 0.516242 | 0.543322 | 0.027080 |
(0.700000, 0.300000) | 0.500028 | 0.500353 | 0.000324 |
(0.100000, 0.500000) | 0.735927 | 0.742214 | 0.006287 |
(0.500000, 0.500000) | 0.534824 | 0.543322 | 0.008498 |
(0.900000, 0.500000) | 0.500269 | 0.500353 | 0.000084 |
(0.300000, 0.700000) | 0.734717 | 0.742214 | 0.007498 |
(0.700000, 0.700000) | 0.542311 | 0.543322 | 0.001011 |
(0.100000, 0.900000) | 0.749901 | 0.749946 | 0.000045 |
(0.500000, 0.900000) | 0.741808 | 0.742214 | 0.000406 |
(0.900000, 0.900000) | 0.539545 | 0.543322 | 0.003778 |
(x, y) | Numerical | Exact | Abs. Error |
(0.100000, 0.100000) | 0.965972 | 0.956678 | 0.009294 |
(0.500000, 0.100000) | 0.999730 | 0.999647 | 0.000083 |
(0.900000, 0.100000) | 0.999998 | 0.999998 | 0.000000 |
(0.300000, 0.300000) | 0.983758 | 0.956678 | 0.027080 |
(0.700000, 0.300000) | 0.999972 | 0.999647 | 0.000324 |
(0.100000, 0.500000) | 0.764073 | 0.757786 | 0.006287 |
(0.500000, 0.500000) | 0.965176 | 0.956678 | 0.008498 |
(0.900000, 0.500000) | 0.999731 | 0.999647 | 0.000084 |
(0.300000, 0.700000) | 0.765283 | 0.757786 | 0.007498 |
(0.700000, 0.700000) | 0.957689 | 0.956678 | 0.001011 |
(0.100000, 0.900000) | 0.750099 | 0.750054 | 0.000045 |
(0.500000, 0.900000) | 0.758192 | 0.757786 | 0.000406 |
(0.900000, 0.900000) | 0.960455 | 0.956678 | 0.003778 |
(x, y) | Numerical | Exact | Abs. (?) Error |
(0.100000,0.100000) | 0.571735 | 0.578513 | 0.006778 |
(0.500000,0.100000) | 0.507501 | 0.509055 | 0.001554 |
(0.900000,0.100000) | 0.500633 | 0.500769 | 0.000136 |
(0.300000,0.300000) | 0.557219 | 0.578513 | 0.021294 |
(0.700000,0.300000) | 0.505997 | 0.509055 | 0.003058 |
(0.100000,0.500000) | 0.702971 | 0.711992 | 0.009021 |
(0.500000,0.500000) | 0.566951 | 0.578513 | 0.011562 |
(0.900000,0.500000) | 0.507953 | 0.509055 | 0.001102 |
(0.300000,0.700000) | 0.700424 | 0.711992 | 0.011568 |
(0.700000,0.700000) | 0.577499 | 0.578513 | 0.001014 |
(0.100000,0.900000) | 0.745264 | 0.746374 | 0.001110 |
(0.500000,0.900000) | 0.710125 | 0.711992 | 0.001867 |
(0.900000,0.900000) | 0.578351 | 0.578513 | 0.000162 |
(x, y) | Numerical | Exact | Abs. Error |
(0.100000,0.100000) | 0.928265 | 0.921487 | 0.006778 |
(0.500000,0.100000) | 0.992499 | 0.990945 | 0.001554 |
(0.900000,0.100000) | 0.999367 | 0.999231 | 0.000136 |
(0.300000,0.300000) | 0.942781 | 0.921487 | 0.021294 |
(0.700000,0.300000) | 0.994003 | 0.990945 | 0.003058 |
(0.100000,0.500000) | 0.797029 | 0.788008 | 0.009021 |
(0.500000,0.500000) | 0.933049 | 0.921487 | 0.011562 |
(0.900000,0.500000) | 0.992047 | 0.990945 | 0.001102 |
(0.300000,0.700000) | 0.799576 | 0.788008 | 0.011568 |
(0.700000,0.700000) | 0.922501 | 0.921487 | 0.001014 |
(0.100000,0.900000) | 0.754736 | 0.753626 | 0.001110 |
(0.500000,0.900000) | 0.789875 | 0.788008 | 0.001867 |
(0.900000,0.900000) | 0.921649 | 0.921487 | 0.000162 |
4.2. Problem 2
The computational domain is given as and the initial and boundary conditions are:
(13a)
and
(13b)
(13c)
We deploy 20×20 grids for the numerical computations for and and . There are no analytical solutions to compare our results with in this case. However figs. 2a, 2b, 2c and 2d of u and v velocity profiles are very close to those of Shukla et al. [60]. In addition Tables 2a and 2b show a comparison of numerical values with those of Glka [58], Bahadir [17], Jain and Holla [59], Srivastava et al. [18]. They all appear to give very close profiles for the u and v velocities.
Coord (x, y) | Present Work | Gulkac [58] | Bahadir [17] | Jain and Holla [59] | Srivastava et al. [18] |
(0.1,0.1) | 0.971460 | 0.966595 | 0.96668 | 0.97258 | 0.976146 |
(0.3,0.1) | 1.152764 | 1.14835 | 1.14827 | 1.16214 | 1.15280 |
(0.2,0.2) | 0.863023 | 0.85918 | 0.85911 | 0.86281 | 0.86307 |
(0.4,0.2) | 0.979359 | 0.97644 | 0.97637 | 0.96483 | 0.97981 |
(0.1,0.3) | 0.663136 | 0.66026 | 0.66619 | 0.66318 | 0.66316 |
(0.3,0.3) | 0.771935 | 0.76939 | 0.76932 | 0.77030 | 0.77230 |
(0.2,0.4) | 0.581622 | 0.57974 | 0.57966 | 0.58070 | 0.58180 |
(0.4,0.4) | 0.75770 | 0.75686 | 0.75678 | 0.74436 | 0.75856 |
Coord (x, y) | Present Work | Gulkac [58] | Bahadir [17] | Jain and Holla [59] | Srivastava et al. [18] |
(0.1,0.1) | 0.098688 | 0.09832 | 0.09824 | 0.09773 | 0.09869 |
(0.3,0.1) | 0.141566 | 0.14119 | 0.4112 | 0.14039 | 0.14158 |
(0.2,0.2) | 0.167525 | 0.16689 | 0.16681 | 0.16660 | 0.16754 |
(0.4,0.2) | 0.170948 | 0.17073 | 0.17065 | 0.17397 | 0.17110 |
(0.1,0.3) | 0.263771 | 0.26269 | 0.26261 | 0.26940 | 0.26378 |
(0.3,0.3) | 0.226400 | 0.22582 | 0.22576 | 0.22463 | 0.22654 |
(0.2,0.4) | 0.328432 | 0.32754 | 0.32745 | 0.32402 | 0.32851 |
(0.4,0.4) | 0.324710 | 0.32447 | 0.32441 | 0.31822 | 0.32500 |
4.3. Problem 3
The computational domain for the coupled nonlinear Burger’s equation is specified as . It admits the following initial and boundary conditions:
Initial conditions:
(14a)
(14b)
Boundary conditions:
(14c)
(14d)
(14e)
(14f)
(14g)
(14h)
The exact solutions are specified as
(14i)
(14j)
This example tests the ability of the scheme to handle relatively high Reynolds number. Comparison of numerical and analytical profiles of the velocity profiles (Figs. 3a, 3b, 3c and 3d) and with those of Tamsir and Srivastava [57] confirm accurate results. Tables 3a and 3b give in addition to values relative and average absolute errors confirm the reliability of this formulation.
(x, y) | Numerical | Exact | Absolute Error |
(0.100000,0.100000) | -0.006594 | -0.006498 | 0.000097 |
(0.500000,0.100000) | 0.007496 | 0.007392 | 0.000103 |
(0.900000,0.100000) | -0.004835 | -0.005464 | 0.000629 |
(0.300000,0.300000) | 0.007188 | 0.008171 | 0.000982 |
(0.700000,0.300000) | 0.004253 | 0.003791 | 0.000462 |
(0.100000,0.500000) | -0.022465 | -0.024768 | 0.002303 |
(0.500000,0.500000) | 0.025169 | 0.023923 | 0.001247 |
(0.900000,0.500000) | -0.013511 | -0.013940 | 0.000429 |
(0.300000,0.700000) | 0.007188 | 0.008171 | 0.000982 |
(0.700000,0.700000) | 0.004253 | 0.003791 | 0.000462 |
(0.100000,0.900000) | -0.006594 | -0.006498 | 0.000097 |
(0.500000,0.900000) | 0.007496 | 0.007392 | 0.000103 |
(0.900000,0.900000) | -0.004835 | -0.005464 | 0.000629 |
relative error(u) =0.000992
avg absolute error(u) =0.000028
(x, y) | Numerical | Exact | Absolute Error |
(0.100000,0.100000) | -0.007530 | -0.007265 | 0.000265 |
(0.500000,0.100000) | -0.000648 | -0.000000 | 0.000648 |
(0.900000,0.100000) | 0.005427 | 0.006109 | 0.000681 |
(0.300000,0.300000) | -0.009014 | -0.009135 | 0.000121 |
(0.700000,0.300000) | 0.004326 | 0.004238 | 0.000088 |
(0.100000,0.500000) | -0.000000 | -0.000000 | 0.000000 |
(0.500000,0.500000) | -0.000000 | -0.000000 | 0.000000 |
(0.900000,0.500000) | 0.000000 | 0.000000 | 0.000000 |
(0.300000,0.700000) | 0.009014 | 0.009135 | 0.000121 |
(0.700000,0.700000) | -0.004326 | -0.004238 | 0.000088 |
(0.100000,0.900000) | 0.007530 | 0.007265 | 0.000265 |
(0.500000,0.900000) | 0.000648 | 0.000000 | 0.000648 |
(0.900000,0.900000) | -0.005427 | -0.006109 | 0.000681 |
relative error(v) =0.000995
avg absolute error(v) =0.000007
5. Conclusion
By using numerical examples with analytical solutions, we have shown the utility of the recent formulation in handling nonlinear transient coupled nonlinear 2D partial differential equations of the Burger’s or Navier-Stokes types [60, 61] in a straightforward manner. The emphasis here is to arrive at a hybrid boundary integral approach devoid of undue complexity in handling some of the major numerical complexities that degrade the competitiveness of the boundary element method. Another equally important consideration is to explore how a relatively straightforward modification of BEM theory could render it accessible to many BEM users when handling the type of problems that relate to real life applications. Accurate representation of the scalar field has been produced in those areas where direct BEM application had previously been deficient.
References