# A Cell-Centered Semi-Lagrangian Finite Volume Method for Solving Two-Dimensional Coupled Burgers’ Equations

Authors: Ilham Asmouh; Mofdi El-Amrani; Mohammed Seaid; Naji Yebari
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

A cell-centered finite volume semi-Lagrangian method is presented for the numerical solution of two-dimensional coupled Burgers’ problems on unstructured triangular meshes. The method combines a modified method of characteristics for the time integration and a cell-centered finite volume for the space discretization. The new method belongs to fractional-step algorithms for which the convection and the viscous parts in the coupled Burgers’ problems are treated separately. The crucial step of interpolation in the convection step is performed using two local procedures accounting for the element where the departure point is located. The resulting semidiscretized system is then solved using a third-order explicit Runge-Kutta scheme. In contrast to the Eulerian-based methods, we apply the new method for each time step along the characteristic curves instead of the time direction. The performance of the current method is verified using different examples for coupled Burgers’ problems with known analytical solutions. We also apply the method for simulation of an example of coupled Burgers’ flows in a complex geometry. In these test problems, the new cell-centered finite volume semi-Lagrangian method demonstrates its ability to accurately resolve the two-dimensional coupled Burgers’ problems.


---

## Body

## 1. Introduction

In this paper, given a two-dimensional bounded domainΩ in R2 with Lipschitz boundary Γ and a time interval 0,T, we are interested in solving the two-dimensional unsteady nonlinear coupled Burgers’ equations
(1)∂u∂t+u∂u∂x+v∂u∂y−ν∂2u∂x2+∂2u∂y2=0,(2)∂v∂t+u∂v∂x+v∂v∂y−ν∂2v∂x2+∂2v∂y2=0,where u=u,vΤ is the velocity field, with u the velocity in x-direction and v the velocity in y-direction, and ν=1/Re is a coefficient of viscous diffusion with Re is the Reynolds number. It should be stressed that the Reynolds number is a nondimensional number which controls the relative dominance of convective part with respect to the diffusive part in the equation (1). To provide a well-posed mathematical problem, equation (1) is equipped with suitable boundary and initial conditions. In practice, these conditions are problem-dependent, and their discussion is postponed for the section on numerical results. The coupled Burgers’ equation (1) arises in many physical phenomena such as turbulent flows, fluid mechanics, gas dynamics, nonlinear acoustic waves, and traffic flow, see [1–3] among others. The problem is widely viewed as a simplified version of the incompressible Navier-Stokes equations and is therefore an important benchmark problem towards the development of robust numerical methods for complicated incompressible flows. For many applications in incompressible flows, the convection terms are distinctly more important than the diffusion terms particularly when the Reynolds number reaches high values. In this case, solutions of coupled Burgers’ equations (1) exhibit steep fronts, sharp gradients, shock discontinuities, and boundary layers among other difficulties that most of Eulerian methods fail to resolve accurately. In addition, the nonlinear nature of the equation (1) is also a challenge for most Eulerian-based methods which require linearization techniques at each time step resulting in computationally demanding solvers.Many contributions have been published in the literature to solve the two-dimensional coupled Burgers’ equation (1). For instance, finite difference methods have been investigated in [4–7] among others. An operator compact method has been proposed in [8] for the numerical solution of the coupled Burgers’ equations. This technique leads to a system of nonlinear equations which is solved using the Newton method at each time step. A fourth-order finite difference method has also been proposed in [9] for solving the coupled Burgers’ equations. The authors in this reference eliminated the convective part by applying the Hopf-Cole transformation which results in a class of heat equations to be discretized using a fourth-order finite difference scheme. Apparently, the main drawback in these finite difference methods is the extension to irregular geometries with complex boundaries. The same Hopf-Cole transformation has been adopted in [10] using the finite volume discretization of the coupled Burgers’ equations. A serious weakness in this method is the necessity to approximate the initial and boundary conditions via a quadrature rule which may create accumulation of errors during the simulations. In [11], the nonlinear advection part is separated from the diffusion part using a time splitting operator to solve the coupled Burgers’ equations. The nonlinear part is linearized based on a special Taylor expansion, whereas the diffusion part is discretized using a finite element method. The authors in [12] have introduced a MacCormack approach based on a predictor-corrector formulation to solve the coupled Burgers’ equations. However, for highly convective problems on unstructured meshes, this method can involve numerical instability specially at the areas with strong gradients within the computational domain. In the framework of meshless methods, the authors in [13] have combined the advantages of differential quadratures and local multiquadric radial basis functions to solve the coupled Burgers’ equations. The main advantage of these meshless methods lies on their potential to handle complex geometries. However, these methods require solutions of ill-conditioned linear systems which become computationally demanding for coupled Burgers’ equations at high Reynolds numbers. The modified cubic B-spline method has also been applied in [14] to the coupled Burgers’ equations. The idea in this method consists of approximating the spatial derivatives using the weighted sum of the functional values at certain discrete points with the weighting coefficients are determined using spline functions. Recently, a strong form meshless collocation method has been investigated in [15] for solving the coupled Burgers’ equations. The main idea of this method is to split the solution into a primary approximation which is treated as the homogeneous solution satisfying the considered boundary conditions and correcting the approximation by the series of modified radial basis functions. However, time truncation errors in these Eulerian methods dominate their numerical solutions and are subject to the Courant-Friedrichs-Lewy (CFL) stability conditions, which put severe restrictions on the size of time steps taken in the numerical simulations of the coupled Burgers’ equations.Semi-Lagrangian methods employ the modified method of characteristics and have been used to solve a wide variety of physical and engineering applications. Indeed, finite element semi-Lagrangian approaches have been used for example in [16] for convection-diffusion problems, in [17] for incompressible Navier-Stokes equations, in [18] for tidal flows, and in [19] for natural and mixed convection flows. The central idea in these finite element semi-Lagrangian methods consists on rewriting the governing equations in terms of Lagrangian coordinates as defined by the characteristics associated with the problem under study. The time derivative and the advection terms are associated in a directional derivative manner along the characteristics leading to a characteristic time-stepping procedure. The time truncation errors are greatly reduced due to the Lagrangian treatment, see, for instance, [20–23]. In addition, the semi-Lagrangian method offers the possibility of using time steps that go beyond those allowed by the CFL stability conditions in Eulerian finite element methods for convection-dominated flows. The main objective of the current study is to devise a numerical approach able to accurately approximate solutions of the coupled Burgers’ equations on unstructured triangular grids. The aim is to develop a powerful numerical algorithm which is robust, easy to implement, and accurately solves the coupled Burgers’ equations without relying on highly demanding solvers. The proposed cell-centered finite volume semi-Lagrangian algorithm can be interpreted as a fractional-step scheme in which we treat the convective and diffusive parts in a separate manner. In [24], a conservative semi-Lagrangian finite volume approach for convection-dominated diffusion equations has been presented. The method is limited to the linear scalar convection-diffusion equations. To achieve an efficient numerical algorithm for nonlinear problems, the present study is carried out for the coupled Burgers’ equations. This is reflected by the differences in the calculation of departure points in [24] with the procedure proposed in the present study. Since it was necessary to decouple the system of equations, we decided to design explicit procedures to eliminate iterative procedures which are not possible with the method reported in [24]. The use of explicit procedures in the cell-centered finite volume semi-Lagrangian method improves consistently the performance of the new numerical tools for solving nonlinear problems. In addition, a third-order explicit Runge-Kutta scheme is implemented in the present work for the time integration. This feature yields a high precision in the numerical solution for the selected Reynolds numbers, and therefore, it makes the method more efficient and robust. The interpolation stage in the semi-Lagrangian is carried using two procedures, namely, the least square (LS) method [25, 26] and radial basis functions (RBF) [27, 28], using the host element where the departure point is allocated and its neighboring cells. It should be stressed out that the presented method should not be confused with splitting techniques widely used in the literature to split convection and diffusion parts in the time direction. Here, we integrate the system (1) along the characteristics using the method of characteristics to deal with the convection part, and no splitting is used in our method. The remaining diffusion part is dealt with using the Runge-Kutta scheme while the convection term is solved using the modified method of characteristics. The proposed method should be interpreted as time stepping scheme along the characteristic curves rather than in the time direction only. This is a well-established technique used in the semi-Lagrangian methods applied to transport and convection-diffusion problems. Finally, a numerical comparison to the finite element semi-Lagrangian approach is also presented in this study. This comparison is particularly interesting in the nonlinear case where the strong gradients are expected in the computed results at high Reynolds numbers. For the proposed cell-centered finite volume semi-Lagrangian method, increasing the Reynolds number in the simulations does not deteriorate the performance of the method. Evaluation of the interpolation procedures has also been discussed using both qualitative and quantitative comparisons. For the considered examples, it has been found that the RBF method is more accurate than the LS method, and it provides highly accurate algorithm for all examples. Moreover, the new cell-centered finite volume semi-Lagrangian approach combined with RBF interpolation has the potential to produce efficient numerical solutions free from spurious oscillations, and it guarantees both strong stability and high accuracy for problems containing sharp gradients. The performance of the proposed cell-centered finite volume semi-Lagrangian method is demonstrated for several test examples of coupled Burgers’ equations. Combining a cell-centered finite volume discretization with the semi-Lagrangian method, to the best of our knowledge, is reported for the first time. Numerical results presented in this study show that an interesting feature of the semi-Lagrangian method is to allow large time steps without deteriorating accuracy of the computed solutions.This paper is organized as follows. Formulation of the cell-centered finite volume semi-Lagrangian method for the convective part is presented in Section2. Two main interpolation techniques are also introduced in this section. In Section 3, we present the spatial discretization of the viscous part of the coupled Burgers’ problems. The resulting semidiscrete equations are then integrated using a third-order explicit Runge-Kutta scheme. The numerical performance of the proposed method is examined in Section 4 using different test examples of coupled Burgers’ problems. The obtained results confirm the high accuracy and the strong stability achieved by our cell-centered finite volume semi-Lagrangian method. Section 5 includes concluding remarks and future works.


## 2. Cell-Centered Finite Volume Semi-Lagrangian Method

To formulate the cell-centered finite volume semi-Lagrangian method, we consider the pure convection part of the problem (1) reformulated using the material derivative as
(3)DuDt≔∂u∂t+u·∇u=0,t,x∈0,T×Ω,where Du/Dt is the total derivative which measures the rate of change of the solution u following the trajectories of the flow particles. Notice that the central idea behind the semi-Lagrangian approaches consists on imposing a regular mesh at the new time level and backtracking the flow trajectories to the previous time level. The quantities needed at the old time level are approximated by interpolation from their known values on a regular mesh. The formulation of main steps of the cell-centered finite volume semi-Lagrangian proposed in the present study is detailed in what follows.


### 2.1. Approximation of Departure Points

For the spacial discretization, we discretize the domainΩ¯=Ω∪Γ into a set of conforming triangles Ki such as Ω¯=∪i=1MKi where Th denotes a partition of Ω¯ into conforming triangles, with M denotes the total number of triangles in the computational mesh. Here, each triangle represents a control volume where the variables are located at the geometric centers of the cells. We also divide the time interval into subintervals tn,tn+1 with a length Δt and denote by wn the value of a generic function w at time tn and win to denote the average value of w in the cell Ki at time tn as
(4)win=1Ki∫Kiwx,tndx,where Ki denotes the area of the control volume Ki. Hence, the characteristic curves Xitn;xi,tn+1 associated with the advection problem (3) are calculated for each control volume Ki, i=1,⋯,M by solving the backward differential equations
(5)dXiτ;xi,tn+1dτ=uiXiτ;xi,tn+1,τ,∀τ∈tn,tn+1,(6)Xitn+1;xi,tn+1=xi,where Xiτ;xi,tn+1=Xiτ;xi,tn+1,Yiτ;xi,tn+1Τ is the departure point defined at time τ of a particle that will reach xi=xi,yiΤ the center of the control volume Ki at time tn+1. It is worth noting that the cell-centered finite volume semi-Lagrangian method does not follow the flow particles forward in time, as a Lagrangian method does, instead it traces backwards the position at time tn of particles that will reach the points of a fixed mesh at time tn+1, see Figure 1 for an illustration. Therefore, the cell-centered finite volume semi-Lagrangian method avoids the grid distortion difficulties that the conventional Lagrangian schemes have.Figure 1
A schematic diagram showing departure points for the cell-centered finite volume semi-Lagrangian method.Note that accurate approximation of the characteristic curvesXiτ;xi,tn+1 is crucial to the overall accuracy of the cell-centered finite volume semi-Lagrangian method. Some authors approximate the solutions of (3) using a second-order explicit Runge-Kutta scheme, which is not accurate enough to maintain a particle on its curved trajectory, see, for instance, [19]. In [22, 29], a second-order extrapolation based on the midpoint rule is used to approximate the solution of (5), but this method involves an iterative procedure which may become computationally demanding. In the present study, we consider the third-order explicit Runge-Kutta method introduced in [30]. Thus, the procedure to approximate the solution of the ordinary differential equation (5) is achieved by
(7)Ki1=xi−Δtuin,(8)Ki2=34xi+14Ki1−14Δtuin,(9)Xitn;xi,tn+1=13xi+23Ki2−23Δtuin.It should be stressed that since the departure pointsXitn;xi,tn+1 and the stages Ki1 and Ki2 would not necessarily lie on a mesh point in the computational domain, the solution at the departure points must be obtained by interpolation from known values at the gridpoints of the element where the points Xitn;xi,tn+1 and the stages Ki1 and Ki2 are localized. In this study, a class of interpolation techniques are presented using the host element associated with the departure points and its neighboring control volumes as detailed below.


### 2.2. Procedures for Interpolation Stage

Let us use the notationXin≔Xitn;xi,tn+1 to denote the characteristics feet calculated in (7). Hence, the finite volume solution of (3) in the control volume Ki at instant tn+1 is given by
(10)uin+1=uXin,tn.Since in general the departure pointXin does not coincide with a mesh point, the numerical solution uXin,tn is computed by interpolation from known values at the control volume K^i where Xin resides and its neighboring elements. In the current work, we perform this step using the following interpolation techniques.


#### 2.2.1. Least Square (LS) Interpolation

This method is usually used to find the best polynomial approximation to a given set of input data, compare [25, 26]. Let us denote by Iin the set of indices of control volumes which are close to the host control volume K^i where the departure point Xin belongs at time t=tn, and let Un be the vector of solutions ujn at mesh points xj with j∈Iin. Hence, the finite volume solution (10) is obtained using the LS interpolation as
(11)uXin,tn=∑k=0mβknϕkXin,where βkn are the fitting coefficients, ϕkx,y is the polynomial basis functions, and m is the total number of fitting data. Notice that we assume that m is less than the dimension of the set Iin. To solve the problem above, we use a linear regression to calculate the coefficients βkn. To this end, we define the merit function F(12)F2=∑j∈Iinujn−∑k=0mβknϕkxj2.It should be noted that ifβn minimizes the functional (12), then the LS procedure precisely fits the vector βn=β0n,β1n,⋯,βmnΤ for a given data set. Hence, we obtain the fitting coefficients βkn by solving the m normal equations
(13)∂F2∂βkn=0,k=0,1,⋯,m.Using equation (12), the normal equations yield
(14)∑j∈Iinujn−∑l=0mβlnϕlxjϕkxj=0,k=0,1,⋯,m,which we can rewrite for each departure point Xin as a linear system of the form
(15)Aβn=r,where the matrix A is given by A=BΤB with B denotes the matrix with entries ϕlxj, 1≤l≤m, j∈Iin. The right-hand side vector r is given by r=BΤf in which f refers to the vector with entries ujn, j∈Iin. Thereby, we define a linear approximation for the function unx,y as follows:
(16)unx,y=β0n+β1nx−Xin+β2ny−Yin.To calculate the coefficientsβ0n, β1n, and β2n, we first evaluate (16) at the departure point Xin=Xin,YinΤ to obtain β0n; then, we solve the linear system (15) to obtain β1n and β2n. For instance, the inverse matrix A−1 is defined as
(17)A−1=1γi∑j∈Iinyj−Yi2−∑j∈Iinxj−Xiyj−Yi−∑j∈Iinxj−Xiyj−Yi∑j∈Iinxj−Xi2,where
(18)γi=∑j∈Iinxj−Xi2∑j∈Iinyj−Yi2−∑j∈Iinxj−Xiyj−Yi2.A similar procedure is used for the componentv in the finite volume solution (10). In the current study, the number m is related to the total number of neighbor elements of the host control volume K^i. Thus, m can be either 2, 3, or 4, see Figure 2 for an illustration. It is worth noting that the quadratic interpolant in the least square framework is more accurate than linear interpolant and would improve the numerical results. In general, the least square techniques have advantages of excellent computational and applicability properties. However, these techniques suffer from the oversmoothing properties, see [31] among others. It is also well known that the least square approach can always achieve a perfect fit to the measurements by selecting the polynomial degree equals to m−1, where m is the total number of data measurements. In the present study, we have considered the cases with m=2,3, or 4. More precisely, m is at most, equals to 2 at the boundaries. Hence, high degree least square polynomial would likely be unacceptable at the boundaries and would introduce nonphysical oscillations. Moreover, the polynomials tend to behave erratically near the boundaries. Notice also that in highly stiff problems, as the case of nonlinear Burgers’ equations at high Reynolds numbers, it may be desirable to have a lower degree least square polynomial in steep gradient regions in order to avoid any overfitting, see, for instance, [32, 33].Figure 2
The host control volumeK^i where the departure point Xin resides and its neighboring control volumes used for the LS approach (a) and the set Iin of points surrounded by the dashed circle used for RBF interpolation (b).


#### 2.2.2. Radial Basis Function (RBF) Interpolation

Interpolation using the RBF surface splines is a typical focus in the theory of spline, see, for instance, [27, 28]. Here, the solution u in (10) is calculated using the RBF interpolation as
(19)uXin,tn=∑j∈IinζjXin−xjlogXin−xj+γ0+γ1xi+γ2yi,where Iin denotes the set with dimension Ni introduced previously, compare Figure 2 for an illustration. Notice that Iin denotes the set of closing neighbors to the cell in question which varies depending on the interpolation procedure used whereas, m is at most 4 in the case of least square technique. In (19), the coefficients ζjj∈Iin, γ0, γ1, and γ2 are obtained by solving the linear system
(20)ψ11ψ12⋯ψ1Ni1x1y1ψ21ψ22⋯ψ2Ni1x2y2⋮⋮⋱⋮⋮⋮⋮ψNi1ψNi2⋯ψNiNi1xNiyNiζ1⋮ζNiγ0γ1γ2=u1u2⋮uNi,where ψij are the radial basis functions defined as
(21)ψij=Xin−xjlogXin−xj.By introducing the following matrices(22)ψ=ψ11ψ12⋯ψ1Niψ21ψ22⋯ψ2Ni⋮⋮⋱⋮ψNi1ψNi2⋯ψNiNi,L=1x1y11x2y2⋮⋮⋮1xNiyNi,and the following vectors
(23)ζ=ζ1⋮ζNi,U=u1u2⋮uNi,γ=γ0γ1γ2,the linear system (20) can be rewritten in a more compact form as
(24)ψLζγ=U.It should be stressed that the linear system (24) contains Ni equations with Ni+3 unknowns. Consequently, to construct a nonsingular system, we should add the natural constraints
(25)∑j=1Niζj=0,∑j=1Niζjxj=0,∑j=1Niζjyj=0.Hence, the coefficientsζjj∈Iin, γ0, γ1, and γ2 in (19) are obtained by solving the block linear system
(26)ψLLT0ζγ=U0.It should be noted that other radial basis functions can also be used in (19) without major conceptual modifications. Note that a similar system to (26) is solved for the component v in the finite volume solution (10).


## 2.1. Approximation of Departure Points

For the spacial discretization, we discretize the domainΩ¯=Ω∪Γ into a set of conforming triangles Ki such as Ω¯=∪i=1MKi where Th denotes a partition of Ω¯ into conforming triangles, with M denotes the total number of triangles in the computational mesh. Here, each triangle represents a control volume where the variables are located at the geometric centers of the cells. We also divide the time interval into subintervals tn,tn+1 with a length Δt and denote by wn the value of a generic function w at time tn and win to denote the average value of w in the cell Ki at time tn as
(4)win=1Ki∫Kiwx,tndx,where Ki denotes the area of the control volume Ki. Hence, the characteristic curves Xitn;xi,tn+1 associated with the advection problem (3) are calculated for each control volume Ki, i=1,⋯,M by solving the backward differential equations
(5)dXiτ;xi,tn+1dτ=uiXiτ;xi,tn+1,τ,∀τ∈tn,tn+1,(6)Xitn+1;xi,tn+1=xi,where Xiτ;xi,tn+1=Xiτ;xi,tn+1,Yiτ;xi,tn+1Τ is the departure point defined at time τ of a particle that will reach xi=xi,yiΤ the center of the control volume Ki at time tn+1. It is worth noting that the cell-centered finite volume semi-Lagrangian method does not follow the flow particles forward in time, as a Lagrangian method does, instead it traces backwards the position at time tn of particles that will reach the points of a fixed mesh at time tn+1, see Figure 1 for an illustration. Therefore, the cell-centered finite volume semi-Lagrangian method avoids the grid distortion difficulties that the conventional Lagrangian schemes have.Figure 1
A schematic diagram showing departure points for the cell-centered finite volume semi-Lagrangian method.Note that accurate approximation of the characteristic curvesXiτ;xi,tn+1 is crucial to the overall accuracy of the cell-centered finite volume semi-Lagrangian method. Some authors approximate the solutions of (3) using a second-order explicit Runge-Kutta scheme, which is not accurate enough to maintain a particle on its curved trajectory, see, for instance, [19]. In [22, 29], a second-order extrapolation based on the midpoint rule is used to approximate the solution of (5), but this method involves an iterative procedure which may become computationally demanding. In the present study, we consider the third-order explicit Runge-Kutta method introduced in [30]. Thus, the procedure to approximate the solution of the ordinary differential equation (5) is achieved by
(7)Ki1=xi−Δtuin,(8)Ki2=34xi+14Ki1−14Δtuin,(9)Xitn;xi,tn+1=13xi+23Ki2−23Δtuin.It should be stressed that since the departure pointsXitn;xi,tn+1 and the stages Ki1 and Ki2 would not necessarily lie on a mesh point in the computational domain, the solution at the departure points must be obtained by interpolation from known values at the gridpoints of the element where the points Xitn;xi,tn+1 and the stages Ki1 and Ki2 are localized. In this study, a class of interpolation techniques are presented using the host element associated with the departure points and its neighboring control volumes as detailed below.


## 2.2. Procedures for Interpolation Stage

Let us use the notationXin≔Xitn;xi,tn+1 to denote the characteristics feet calculated in (7). Hence, the finite volume solution of (3) in the control volume Ki at instant tn+1 is given by
(10)uin+1=uXin,tn.Since in general the departure pointXin does not coincide with a mesh point, the numerical solution uXin,tn is computed by interpolation from known values at the control volume K^i where Xin resides and its neighboring elements. In the current work, we perform this step using the following interpolation techniques.


### 2.2.1. Least Square (LS) Interpolation

This method is usually used to find the best polynomial approximation to a given set of input data, compare [25, 26]. Let us denote by Iin the set of indices of control volumes which are close to the host control volume K^i where the departure point Xin belongs at time t=tn, and let Un be the vector of solutions ujn at mesh points xj with j∈Iin. Hence, the finite volume solution (10) is obtained using the LS interpolation as
(11)uXin,tn=∑k=0mβknϕkXin,where βkn are the fitting coefficients, ϕkx,y is the polynomial basis functions, and m is the total number of fitting data. Notice that we assume that m is less than the dimension of the set Iin. To solve the problem above, we use a linear regression to calculate the coefficients βkn. To this end, we define the merit function F(12)F2=∑j∈Iinujn−∑k=0mβknϕkxj2.It should be noted that ifβn minimizes the functional (12), then the LS procedure precisely fits the vector βn=β0n,β1n,⋯,βmnΤ for a given data set. Hence, we obtain the fitting coefficients βkn by solving the m normal equations
(13)∂F2∂βkn=0,k=0,1,⋯,m.Using equation (12), the normal equations yield
(14)∑j∈Iinujn−∑l=0mβlnϕlxjϕkxj=0,k=0,1,⋯,m,which we can rewrite for each departure point Xin as a linear system of the form
(15)Aβn=r,where the matrix A is given by A=BΤB with B denotes the matrix with entries ϕlxj, 1≤l≤m, j∈Iin. The right-hand side vector r is given by r=BΤf in which f refers to the vector with entries ujn, j∈Iin. Thereby, we define a linear approximation for the function unx,y as follows:
(16)unx,y=β0n+β1nx−Xin+β2ny−Yin.To calculate the coefficientsβ0n, β1n, and β2n, we first evaluate (16) at the departure point Xin=Xin,YinΤ to obtain β0n; then, we solve the linear system (15) to obtain β1n and β2n. For instance, the inverse matrix A−1 is defined as
(17)A−1=1γi∑j∈Iinyj−Yi2−∑j∈Iinxj−Xiyj−Yi−∑j∈Iinxj−Xiyj−Yi∑j∈Iinxj−Xi2,where
(18)γi=∑j∈Iinxj−Xi2∑j∈Iinyj−Yi2−∑j∈Iinxj−Xiyj−Yi2.A similar procedure is used for the componentv in the finite volume solution (10). In the current study, the number m is related to the total number of neighbor elements of the host control volume K^i. Thus, m can be either 2, 3, or 4, see Figure 2 for an illustration. It is worth noting that the quadratic interpolant in the least square framework is more accurate than linear interpolant and would improve the numerical results. In general, the least square techniques have advantages of excellent computational and applicability properties. However, these techniques suffer from the oversmoothing properties, see [31] among others. It is also well known that the least square approach can always achieve a perfect fit to the measurements by selecting the polynomial degree equals to m−1, where m is the total number of data measurements. In the present study, we have considered the cases with m=2,3, or 4. More precisely, m is at most, equals to 2 at the boundaries. Hence, high degree least square polynomial would likely be unacceptable at the boundaries and would introduce nonphysical oscillations. Moreover, the polynomials tend to behave erratically near the boundaries. Notice also that in highly stiff problems, as the case of nonlinear Burgers’ equations at high Reynolds numbers, it may be desirable to have a lower degree least square polynomial in steep gradient regions in order to avoid any overfitting, see, for instance, [32, 33].Figure 2
The host control volumeK^i where the departure point Xin resides and its neighboring control volumes used for the LS approach (a) and the set Iin of points surrounded by the dashed circle used for RBF interpolation (b).


### 2.2.2. Radial Basis Function (RBF) Interpolation

Interpolation using the RBF surface splines is a typical focus in the theory of spline, see, for instance, [27, 28]. Here, the solution u in (10) is calculated using the RBF interpolation as
(19)uXin,tn=∑j∈IinζjXin−xjlogXin−xj+γ0+γ1xi+γ2yi,where Iin denotes the set with dimension Ni introduced previously, compare Figure 2 for an illustration. Notice that Iin denotes the set of closing neighbors to the cell in question which varies depending on the interpolation procedure used whereas, m is at most 4 in the case of least square technique. In (19), the coefficients ζjj∈Iin, γ0, γ1, and γ2 are obtained by solving the linear system
(20)ψ11ψ12⋯ψ1Ni1x1y1ψ21ψ22⋯ψ2Ni1x2y2⋮⋮⋱⋮⋮⋮⋮ψNi1ψNi2⋯ψNiNi1xNiyNiζ1⋮ζNiγ0γ1γ2=u1u2⋮uNi,where ψij are the radial basis functions defined as
(21)ψij=Xin−xjlogXin−xj.By introducing the following matrices(22)ψ=ψ11ψ12⋯ψ1Niψ21ψ22⋯ψ2Ni⋮⋮⋱⋮ψNi1ψNi2⋯ψNiNi,L=1x1y11x2y2⋮⋮⋮1xNiyNi,and the following vectors
(23)ζ=ζ1⋮ζNi,U=u1u2⋮uNi,γ=γ0γ1γ2,the linear system (20) can be rewritten in a more compact form as
(24)ψLζγ=U.It should be stressed that the linear system (24) contains Ni equations with Ni+3 unknowns. Consequently, to construct a nonsingular system, we should add the natural constraints
(25)∑j=1Niζj=0,∑j=1Niζjxj=0,∑j=1Niζjyj=0.Hence, the coefficientsζjj∈Iin, γ0, γ1, and γ2 in (19) are obtained by solving the block linear system
(26)ψLLT0ζγ=U0.It should be noted that other radial basis functions can also be used in (19) without major conceptual modifications. Note that a similar system to (26) is solved for the component v in the finite volume solution (10).


## 2.2.1. Least Square (LS) Interpolation

This method is usually used to find the best polynomial approximation to a given set of input data, compare [25, 26]. Let us denote by Iin the set of indices of control volumes which are close to the host control volume K^i where the departure point Xin belongs at time t=tn, and let Un be the vector of solutions ujn at mesh points xj with j∈Iin. Hence, the finite volume solution (10) is obtained using the LS interpolation as
(11)uXin,tn=∑k=0mβknϕkXin,where βkn are the fitting coefficients, ϕkx,y is the polynomial basis functions, and m is the total number of fitting data. Notice that we assume that m is less than the dimension of the set Iin. To solve the problem above, we use a linear regression to calculate the coefficients βkn. To this end, we define the merit function F(12)F2=∑j∈Iinujn−∑k=0mβknϕkxj2.It should be noted that ifβn minimizes the functional (12), then the LS procedure precisely fits the vector βn=β0n,β1n,⋯,βmnΤ for a given data set. Hence, we obtain the fitting coefficients βkn by solving the m normal equations
(13)∂F2∂βkn=0,k=0,1,⋯,m.Using equation (12), the normal equations yield
(14)∑j∈Iinujn−∑l=0mβlnϕlxjϕkxj=0,k=0,1,⋯,m,which we can rewrite for each departure point Xin as a linear system of the form
(15)Aβn=r,where the matrix A is given by A=BΤB with B denotes the matrix with entries ϕlxj, 1≤l≤m, j∈Iin. The right-hand side vector r is given by r=BΤf in which f refers to the vector with entries ujn, j∈Iin. Thereby, we define a linear approximation for the function unx,y as follows:
(16)unx,y=β0n+β1nx−Xin+β2ny−Yin.To calculate the coefficientsβ0n, β1n, and β2n, we first evaluate (16) at the departure point Xin=Xin,YinΤ to obtain β0n; then, we solve the linear system (15) to obtain β1n and β2n. For instance, the inverse matrix A−1 is defined as
(17)A−1=1γi∑j∈Iinyj−Yi2−∑j∈Iinxj−Xiyj−Yi−∑j∈Iinxj−Xiyj−Yi∑j∈Iinxj−Xi2,where
(18)γi=∑j∈Iinxj−Xi2∑j∈Iinyj−Yi2−∑j∈Iinxj−Xiyj−Yi2.A similar procedure is used for the componentv in the finite volume solution (10). In the current study, the number m is related to the total number of neighbor elements of the host control volume K^i. Thus, m can be either 2, 3, or 4, see Figure 2 for an illustration. It is worth noting that the quadratic interpolant in the least square framework is more accurate than linear interpolant and would improve the numerical results. In general, the least square techniques have advantages of excellent computational and applicability properties. However, these techniques suffer from the oversmoothing properties, see [31] among others. It is also well known that the least square approach can always achieve a perfect fit to the measurements by selecting the polynomial degree equals to m−1, where m is the total number of data measurements. In the present study, we have considered the cases with m=2,3, or 4. More precisely, m is at most, equals to 2 at the boundaries. Hence, high degree least square polynomial would likely be unacceptable at the boundaries and would introduce nonphysical oscillations. Moreover, the polynomials tend to behave erratically near the boundaries. Notice also that in highly stiff problems, as the case of nonlinear Burgers’ equations at high Reynolds numbers, it may be desirable to have a lower degree least square polynomial in steep gradient regions in order to avoid any overfitting, see, for instance, [32, 33].Figure 2
The host control volumeK^i where the departure point Xin resides and its neighboring control volumes used for the LS approach (a) and the set Iin of points surrounded by the dashed circle used for RBF interpolation (b).


## 2.2.2. Radial Basis Function (RBF) Interpolation

Interpolation using the RBF surface splines is a typical focus in the theory of spline, see, for instance, [27, 28]. Here, the solution u in (10) is calculated using the RBF interpolation as
(19)uXin,tn=∑j∈IinζjXin−xjlogXin−xj+γ0+γ1xi+γ2yi,where Iin denotes the set with dimension Ni introduced previously, compare Figure 2 for an illustration. Notice that Iin denotes the set of closing neighbors to the cell in question which varies depending on the interpolation procedure used whereas, m is at most 4 in the case of least square technique. In (19), the coefficients ζjj∈Iin, γ0, γ1, and γ2 are obtained by solving the linear system
(20)ψ11ψ12⋯ψ1Ni1x1y1ψ21ψ22⋯ψ2Ni1x2y2⋮⋮⋱⋮⋮⋮⋮ψNi1ψNi2⋯ψNiNi1xNiyNiζ1⋮ζNiγ0γ1γ2=u1u2⋮uNi,where ψij are the radial basis functions defined as
(21)ψij=Xin−xjlogXin−xj.By introducing the following matrices(22)ψ=ψ11ψ12⋯ψ1Niψ21ψ22⋯ψ2Ni⋮⋮⋱⋮ψNi1ψNi2⋯ψNiNi,L=1x1y11x2y2⋮⋮⋮1xNiyNi,and the following vectors
(23)ζ=ζ1⋮ζNi,U=u1u2⋮uNi,γ=γ0γ1γ2,the linear system (20) can be rewritten in a more compact form as
(24)ψLζγ=U.It should be stressed that the linear system (24) contains Ni equations with Ni+3 unknowns. Consequently, to construct a nonsingular system, we should add the natural constraints
(25)∑j=1Niζj=0,∑j=1Niζjxj=0,∑j=1Niζjyj=0.Hence, the coefficientsζjj∈Iin, γ0, γ1, and γ2 in (19) are obtained by solving the block linear system
(26)ψLLT0ζγ=U0.It should be noted that other radial basis functions can also be used in (19) without major conceptual modifications. Note that a similar system to (26) is solved for the component v in the finite volume solution (10).


## 3. Finite Volume Solution of Viscous Terms

Let us consider coupled Burgers’ problem (1) reformulated using the total derivative as
(27)DuDt−νΔu=0,(28)DvDt−νΔv=0.It should be stressed that, in contrast to most of numerical methods discussed in the introduction where a large system of nonlinear equations should be solved using the Newton method at each time step, in our method, equation (27) is decoupled and can be solved separately for each solution component. From a computational view point, this is very advantageous because the nonlinearity is dealt with using the semi-Lagrangian method, and thus, no need to use nonlinear solvers which are computationally demanding. Here, for brevity in the presentation, the method is formulated only for the component u, and the same procedure is carried out for the component v. Hence, integrating equation (27) over a control volume Ki∈Th and using the divergence theorem, it yields
(29)∫KiDuDtdx−∑j∈Ni∮σijν∇u·nijdσ=0,where σij is the edge between Ki and Kj, j∈Ni with Ni is the set of all neighboring elements, i.e., control volumes which have a common one-dimensional face with Ki. In (29), nij denotes the unit outer normal vector to σij with respect to Kj. Introducing the approximation
(30)ui=1Ki∫Kiudx,where ui denotes the value of u at centroid of the cell Ki, equation (29) becomes
(31)DuiDt=1Ki∑j∈Ni∮σijν∇u·nijdσ.It is clear that equation (31) can be rewritten as
(32)DuiDt=1Ki∑j∈NisijFijσij,where sij=nij·nσ, with nσ denotes the normal to the edge σij, see, for example, [34]. In (32), the flux function Fij is defined by
(33)Fij=1σij∮σijν∇u·nσdσ.In the current work, to discretize the diffusive fluxes, we implement the Green-Gauss diamond reconstruction as presented [34, 35] among others. This reconstruction is demonstrated to be efficient and second-order accurate with no serious restrictions on the angles of control volumes as it can be extended to general unstructured grids with large deformation [34, 35]. Hence, connecting the barycenters of the control volumes that share the edge σ and its endpoints as shown in Figure 3, a diamond-cell denoted by Dh is reconstructed as
(34)1Dh∫Dh∇udx=1Dh∮∂Dhu·nextdσ,where next is the outward normal vector to ∂Dh. Following [34], the discretization of (34) yields the following approximation of the gradient over σ(35)Gσ=1Dh∑γ∈∂Dh12uP1γ+uP2γγnγ,where nγ denotes the outward normal to the edge γ whereas P1γ and P2γ are the endpoints of an edge γ of ∂Dh, see Figure 3 for an illustration. Here, uE and uW are solutions at cell centers while uN and uS are the solutions at the vertices xN and xS expressed in terms of cell center values uK using the RBF interpolation technique discussed in the previous section. Thus, the reconstructed gradient can be written as follows
(36)Gσ=uE−uWhσ−θσuN−uSσnσ+uN−uSσtσ,where hσ=lWE·nσ and σ=lSN·tσ, with tσ is the unit tangential vector to nσ. Here, lWE is the vector connecting the points W and E, and lSN is the vector connecting the points S and N, and
(37)θσ=tannσ,lWE=lWE·tσlWE·nσ.Figure 3
Illustration of a covolume reconstructed from two adjacent control volumes used to discretize the viscous terms in the present study.Using the approximated gradientGσ, the numerical diffusive flux in (33) can be then written as
(38)Fij=1σij∮σijνGσij·nσijdσ=νGσij·nσij.Reformulating the semidiscrete system (31) in a compact form, the equations lead to a system of ordinary differential equations as
(39)DUDt=RU,where U is solution vectors with entries ui and RU=SU, with S denotes the stiffness matrix with entries given above. In the current work, to solve the system (39), we consider the third-order explicit Runge-Kutta method [30]. Thus, the procedure to advance the solution of an equation of the structure (39) from the time tn to the next time tn+1 can be achieved by
(40)U1=U~n+ΔtRU~n,(41)U2=34U~n+14U1+14ΔtRU1,(42)Un+1=13U~n+23U2+23ΔtRU2,where U~n is the solution vector with entries uXin,tn calculated by interpolation as described in the previous section. It should be pointed out that the main advantage of this scheme lies on the fact that (40) is a convex combination of the first-order Euler steps which exhibits strong stability properties. Therefore, the scheme (40) is TVD, third-order accurate in time, and stable under the Courant-Friedrichs-Lewy (CFL) condition given for example in [36]
(43)νΔth2≤12.Note that other time-stepping schemes can also be used for solving the system (39).


## 4. Numerical Results

To examine the performance of the new cell-centered finite volume semi-Lagrangian method proposed in this study, we consider a number of numerical examples of two-dimensional coupled Burgers’ problems. Since analytical solutions for these examples are known, the relativeL1-error and L2-error at time tn are evaluated as
(44)L1−error=∫Ωuin−uexactxi,tndx∫Ωuexactnxi,tndx,L2−error=∫Ωuin−uexactnxi,tn2dx1/2∫Ωuexactnxi,tn2dx1/2,where uin and uexactxi,tn are, respectively, the numerical and exact solutions at time tn in the control volume with centroid xi. For all results reported in this section, the radius of the circle used in the RBF procedure is set to 2h, with h is the space stepsize. The total number of points within this circle are used for the interpolation in the RBF procedure. For comparison reasons, we also compare the results obtained using our cell-centered finite volume semi-Lagrangian method to those obtained using the finite element semi-Lagrangian method. For completeness, the reformulation of this method is briefly described in the Appendix. All the computations were performed on an Intel Core(TM) i7-7500 U @ 2.70 GHz with 16 GB of RAM.


### 4.1. Example 1

In this example, we solve the coupled Burgers’ equations (1) in the squared domain Ω=0,1×0,1 with initial and boundary conditions obtained from the following analytical solution
(45)ux,y,t=34−141+exp−4x−4y+tRe/32,(46)vx,y,t=34+141+exp−4x−4y+tRe/32.This example has been widely used in the literature to examine numerical methods for solving the two-dimensional Burgers’ equations, see, for instance, [9–12, 15]. The exact solution (45) is used to quantify the errors in the proposed cell-centered finite volume semi-Lagrangian method for solving the Burgers’ system for two different Reynolds numbers Re=100 and Re=1000 at time t=1 using a time step Δt=0.01. Notice that the considered values of the Reynolds numbers yield two different flow regimes with different solution features.In Table1, we summarize the relative errors, convergence rates, and CPU times for the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques on different uniform structured meshes with spatial step h. For comparison purposes, we also include in Table 1 the computational results obtained using the finite element semi-Lagrangian method (FEM). It is clear that refining the mesh results in a decrease of the relative L1-error and L2-error and an increase in the computational cost for all considered methods. In addition, better convergence rates are obtained for all considered methods for the simulations using Re=100 than using Re=1000. However, using the RBF interpolation technique, the proposed cell-centered finite volume semi-Lagrangian method is more accurate than the LS interpolation technique and the FEM. To have a fair comparison between the LS and RBF interpolation, the same number of points used in the LS interpolation is used for the RBF procedure, and the obtained results are presented in Table 1. It is clear that using the same number of points in the interpolation stage, the RBF (4) scheme exhibits better results than the LS scheme for all considered meshes and Reynolds numbers. In terms of errors, the RBF (4) procedure achieves higher convergence rates than its LS counterpart with a slightly higher computational cost. For both values of the Reynolds numbers, the RBF scheme exhibits a higher order accuracy than the LS scheme and the FEM. In terms of computational time, Table 1 reveals that the RBF scheme requires almost the triple of CPU time required for the LS scheme and about 1.5 times larger than the FEM.Table 1
Results for Example 1 obtained by the LS and RBF procedures using different meshes and two Reynolds numbers at timet=1. Computational results obtained using the finite element method are also included. The listed CPU times are in seconds.
Re=100Re=1000hL1−errorRateL2−errorRateCPUL1−errorRateL2−errorRateCPULS185.62809E-02—6.66197E-02—1.584.13981E-01—4.10072E-01—1.851161.92205E-021.552.43842E-021.654.171.38469E-011.581.37162E-011.504.961326.12442E-031.658.50242E-031.7115.04.35143E-021.674.46236E-021.6216.31641.79573E-031.772.69050E-031.8263.21.28475E-021.761.36396E-021.7185.6RBF (4)182.51041E-02—3.82963E-02—2.512.32246E-01—2.90329E-01—2.691166.23066E-032.011.07381E-021.836.116.10053E-021.938.10952E-021.877.161321.45335E-032.102.47027E-032.1219.231.47318E-022.051.93136E-022.0721.221642.87052E-042.345.05109E-042.29813.13500E-032.234.29168E-032.1789RBF181.55020E-02—1.96313E-02—3.312.01827E-01—2.84819E-01—4.251161.74416E-033.152.55042E-032.949.592.42007E-023.064.82984E-022.5610.61321.88486E-043.213.25501E-042.9737.51.96833E-033.625.55544E-033.1227.61641.98122E-053.254.01275E-053.021761.44283E-043.774.87645E-043.51239FEM184.24047E-02—5.09578E-02—3.013.10254E-01—3.21647E-01—3.761161.08899E-021.961.37588E-021.897.997.81031E-021.998.44095E-021.938.761322.74141E-031.993.55840E-031.9528.81.92569E-022.022.02081E-022.0220.91646.62007E-042.058.83457E-042.011254.39940E-032.134.68917E-032.15169Figure4 depicts the results obtained for the solution u by the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques on a uniform structured mesh with h=1/64. For comparison reason, we also include the analytical solution in the same figure. Note that since the solution for this example satisfies u+v=3/2, only the solution u is presented. For both considered Reynolds numbers, it is clear that the LS technique involves a numerical diffusion which is more pronounced than in those obtained using the RBF interpolation technique. At Re=1000, strong gradients are formed in the solution which has been well captured by the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique. The LS interpolation technique fails to correctly capture this sharp solution, and the numerical diffusion can clearly be seen in its results presented in Figure 4. To further demonstrate these features, one-dimensional cross-sections of the computed solutions at the off-diagonal line y=1−x are shown in Figure 5. The results computed using the FEM and the analytical solutions are also included in Figure 5 for comparison reasons. For Re=100, the physical diffusion dominates the flow solution, and both results computed using the FEM and RBF methods coincide with the analytical solution. However, the results obtained using the LS interpolation technique are less accurate than those obtained using the FEM and RBF methods. Increasing the Reynolds number to Re=1000, the results obtained using the FEM are oscillatory and do not preserve monotonicity, observe the overshoots and undershoots in the FEM results in Figure 5. From the same figure, we can see that the results obtained using the RBF method are monotone and free from any nonphysical oscillations. It is also evident that the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique preserves the high-order accuracy for selected values of the Reynolds number.Figure 4
Results for the solutionu in Example 1 at time t=1 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) on a mesh with h=1/64 and the Reynolds number Re=100 (1st row) and Re=1000 (2nd row).Figure 5
One-dimensional cross-sections of the computed solutionu at y=1−x for Example 1 at time t=1 on a mesh with h=1/64 using Re=100 (a) and Re=1000 (b).


### 4.2. Example 2

We solve the Burgers’ system (1) in a squared domain Ω=0,1×0,1 subject to initial and boundary conditions obtained from the following analytical solution
(47)ux,y,t=−4πexp−5π2t/Recos2πxsinπyRe2+exp−5π2t/Resin2πxsinπy,(48)vx,y,t=−2πexp−5π2t/Resin2πxcosπyRe2+exp−5π2t/Resin2πxsinπy.It is easy to verify that the pair solutionu,v satisfies the relation u=2v, and therefore, only results obtained for the solution u are presented for this problem. Compared to the previous test example, the solution of this problem is smother with a fast decay as the time progresses. This example has also been investigated in many references, see [9, 12] among others. Similarly, we present numerical results obtained for two different Reynolds numbers Re=100 and Re=1000 at time t=2 using a time step Δt=0.01. Table 2 summarizes the relative errors, convergence rates, and CPU times for the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques on different uniform structured meshes varying the spatial step h. In Table 2, we also include the computational results obtained using the finite element semi-Lagrangian method (FEM) and the RBF scheme (RBF (4)) using the same number of interpolant points as the LS scheme. A decrease in the relative L1- and L2-error is clearly seen in all the results as the mesh is refined with a faster convergence in the results obtained using the RBF interpolation technique compared to the other schemes. It is also worth remarking that higher convergence accuracy is achieved for Re=100 in all considered methods compared to the results for Re=1000. For this test example, the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique is the most accurate compared to the other methods. For instance, at Re=100 on the mesh with h=1/64, convergence rates in the L2-error for the LS, RBF (4), RBF, and FEM schemes are 1.52, 2.31, 2.55, and 2.12, respectively. At Re=1000, these rates become 1.45, 2.23, 2.49, and 2.01 for the LS, RBF (4), RBF, and FEM, respectively. Under the considered flow conditions, the computational time required for the simulations using the RBF procedure is higher than for the LS and FEM. Using the same number of points in the RBF (4) procedure as in the LS scheme, it also demonstrates the better convergence rates achieved in the RBF method for this example.Table 2
Results for Example 2 obtained by the LS and RBF procedures using different meshes and two Reynolds numbers at timet=2. Computational results obtained using the finite element method are also included. The listed CPU times are in seconds.
Re=100Re=1000hL1−errorRateL2−errorRateCPUL1−errorRateL2−errorRateCPULS181.22884E-02—1.54575E-02—3.561.44367E-02—2.00704E-02—4.551164.62427E-031.415.98037E-031.379.325.74246E-031.338.09481E-031.319.891321.65775E-031.482.25048E-031.4128.22.16095E-031.413.17553E-031.3530.61645.62228E-041.567.84711E-041.521157.64011E-041.501.16231E-031.45122RBF(4)185.51041E-03—6.83363E-03—4.516.52321E-03—7.13223E-03—5.691161.46628E-031.911.76865E-031.8310.411.73578E-031.931.93771E-031.8812.111323.71687E-042.104.36075E-042.1233.534.36963E-042.054.48865E-042.1135.221648.55054E-052.128.79391E-052.311389.84533E-052.159.56795E-052.23145RBF184.49127E-04—5.66690E-04—8.185.18864E-04—5.95316E-04—10.61168.80945E-052.351.19132E-042.2522.31.03912E-042.321.22564E-042.2824.11321.56814E-052.492.18024E-052.4570.51.91493E-052.442.52357E-052.4178.31642.69630E-062.543.72287E-062.552993.33855E-062.524.49211E-062.49324FEM189.34059E-03—1.05866E-02—7.369.88733E-03—9.91181E-03—9.631162.55534E-031.872.81701E-031.9118.22.74267E-041.852.74946E-031.8120.11326.52258E-041.967.29088E-041.9553.17.09847E-041.957.26557E-041.9260.21641.57509E-042.051.67724E-042.122071.73809E-042.031.80385E-042.01231Results obtained for the solutionu using the proposed cell-centered finite volume semi-Lagrangian method and the analytical solution are illustrated in Figure 6 for Re=100 and Re=1000 at time t=2 on the mesh with h=1/64. For both values of the Reynolds number, there is no visible difference between the results obtained using the RBF interpolation technique and the exact solutions. On the other hand, the results obtained using the LS interpolation technique are less accurate with excessive numerical diffusion particularly for the simulations using Re=1000. To further demonstrate these features, we illustrate in Figure 7 the one-dimensional cross-sections of the computed solutions at the diagonal y=x. For comparison, we depict in the same figure the computational results obtained using the finite element semi-Lagrangian method. As expected, a smooth behavior can be seen in the solution u for this test example at both Reynolds numbers Re=100 and Re=1000. Obviously, the LS interpolation introduces a noticeable numerical diffusion, and it becomes more intense at Re=1000, compare its failure to resolve sharp gradients in Figure 7. From the same figure, we can see that the results obtained using the RBF interpolation technique are slightly more accurate than those obtained using the FEM. From the results presented, it is clear that the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique performs very well for this test example.Figure 6
Results for the solutionu in Example 2 in a squared domain at time t=2 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) on a mesh with h=1/64 and the Reynolds number Re=100 (1st row) and Re=1000 (2nd row).Figure 7
One-dimensional cross-sections of the computed solutionu at y=x for Example 2 in a squared domain at time t=2 on a mesh with h=1/64 using Re=100 (a) and Re=1000 (b).Our next objective is to test and examine the performance of the present cell-centered finite volume semi-Lagrangian for solving the Burgers’ equation (1) in complex domains using unstructured meshes. To this end, we solve this example on the complex domain shown in Figure 8 for which the initial and boundary conditions are obtained from the analytical solution (47). Here, the computational domain is defined by its parametric curved boundary
(49)Γ=x,y∈2:xθ=Rθcosθ,yθ=Rθsinθ,where
(50)Rθ=2+12sin7θ,θ∈0,2π.Figure 8
Computational mesh used for Example 2 in a complex domain.Note that similar geometry has been considered in [15, 37] among others. An unstructured mesh with 4056 elements is used in our simulations, and the obtained results are presented for two different Reynolds numbers Re=100 and Re=1000 at time t=1 using a fixed time step Δt=0.01. In Figure 9, we display the numerical results for the solution u obtained by the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques along with the analytical solution. In contrast to the previous simulation in the squared domain for which one period in the results is obtained, the results in this complex domain exhibit several periods. As it can be seen in the results obtained using the LS approach, significant differences between these results and the exact solution are observed especially for the convection-dominated regime at Re=1000. Indeed, the computational treatment of such complex patterns on complex geometries often requires efficient numerical algorithms, and clearly, the LS approach is not suitable for these situations. However, the RBF method shows a very good behavior and matches the exact solution at both Re=100 and Re=1000. An examination of the one-dimensional cross-sections of the computed results at y=x in Figure 10 confirms these claims. For a comparison reason, we depict in the same figure the computational results of the finite element semi-Lagrangian approach. Again, the LS method produces results with large numerical diffusion, and it fails to capture the exact solutions. The FEM is in reasonable agreement with exact solution although it substantially manifests small distortion particularly in regions near sharp gradients. From the same figure, the accuracy in the results obtained using the RBF method is superior to those obtained using the FEM and LS methods. It should also be stressed that the complex nature of the geometry in this example reflects the ability of the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique to resolve the complex patterns.Figure 9
Results for the solutionu in Example 2 in a complex domain at time t=1 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) using Re=100 (1st row) and Re=1000 (2nd row).Figure 10
One-dimensional cross-sections of the computed solutionu at y=x for Example 2 in a complex domain at time t=1 using Re=100 (a) and Re=1000 (b).


### 4.3. Example 3

Our final test example consists in solving the Burgers’ equation (1) in the complex domain shown in Figure 11 defined by its parametric boundary
(51)Γ=x,y∈2:xθ=Rθcosθ,yθ=Rθsinθ,where
(52)Rθ=53cos3θ+8−sin3θ31/3,θ∈0,2π.Figure 11
Computational mesh used for Example 3 in a complex domain.A similar geometry has been considered in [15] for the coupled Burgers’ equations. To obtain the corresponding initial and boundary conditions, we use the following exact solution
(53)ux,y,t=121−tanhx4ν+y4ν−t4ν,vx,y,t=121−tanhx4ν+y4ν−t4ν.This test example has also been studied in [38, 39] for a squared domain. In our simulations, we consider an unstructured mesh with 4883 elements, and numerical results are presented at time t=1 for two different Reynolds numbers Re=10 and Re=1000 using a time step fixed to Δt=0.01. Note that to reduce the computational cost, the time steps Δt are chosen as large as possible in our simulations. This makes most explicit Eulerian-based methods noncompetitive, since they are subject to stability restriction conditions for the convection terms. Therefore, the criteria of choosing time steps in our simulations performed in this study was mainly based on accuracy considerations and stability condition (43) required for the solution of viscous parts.In Figure12, we present results obtained for the solution u using the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques compared to the analytical results. The one-dimensional plots in Figure 13 correspond to cross-section solutions at the off-diagonal y=1−x of these results including those obtained using the FEM. A visual comparison of the results in these figures shows severe numerical dissipation in the solutions obtained using the LS method. This numerical dissipation is more pronounced for large values of Re, and it reduces as the mesh is refined. Under the considered flow conditions, the FEM and LS methods exhibit substantially large numerical diffusion particularly at the feet of areas where the gradient of the solution is very sharp. However, such numerical dissipation is completely absent in the results attained using the RBF method. It is also clear that the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique performs best for this test example. It should be mentioned that the proposed cell-centered finite volume semi-Lagrangian method is typically designed to solve this class of convection-dominated flow problems using time steps ten to twenty times larger than its Eulerian counterparts. It should be noted that the proposed cell-centered finite volume semi-Lagrangian method performs very well as the solution remains stable and accurate even in the case of coarse meshes with no need for nonlinear solvers or small time steps in the simulations.Figure 12
Results for the solutionu in Example 3 at time t=1 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) using Re=10 (1st row) and Re=1000 (2nd row).Figure 13
One-dimensional cross-sections of the computed solutionu at y=1−x for Example 3 at time t=1 using Re=10 (a) and Re=1000 (b).


## 4.1. Example 1

In this example, we solve the coupled Burgers’ equations (1) in the squared domain Ω=0,1×0,1 with initial and boundary conditions obtained from the following analytical solution
(45)ux,y,t=34−141+exp−4x−4y+tRe/32,(46)vx,y,t=34+141+exp−4x−4y+tRe/32.This example has been widely used in the literature to examine numerical methods for solving the two-dimensional Burgers’ equations, see, for instance, [9–12, 15]. The exact solution (45) is used to quantify the errors in the proposed cell-centered finite volume semi-Lagrangian method for solving the Burgers’ system for two different Reynolds numbers Re=100 and Re=1000 at time t=1 using a time step Δt=0.01. Notice that the considered values of the Reynolds numbers yield two different flow regimes with different solution features.In Table1, we summarize the relative errors, convergence rates, and CPU times for the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques on different uniform structured meshes with spatial step h. For comparison purposes, we also include in Table 1 the computational results obtained using the finite element semi-Lagrangian method (FEM). It is clear that refining the mesh results in a decrease of the relative L1-error and L2-error and an increase in the computational cost for all considered methods. In addition, better convergence rates are obtained for all considered methods for the simulations using Re=100 than using Re=1000. However, using the RBF interpolation technique, the proposed cell-centered finite volume semi-Lagrangian method is more accurate than the LS interpolation technique and the FEM. To have a fair comparison between the LS and RBF interpolation, the same number of points used in the LS interpolation is used for the RBF procedure, and the obtained results are presented in Table 1. It is clear that using the same number of points in the interpolation stage, the RBF (4) scheme exhibits better results than the LS scheme for all considered meshes and Reynolds numbers. In terms of errors, the RBF (4) procedure achieves higher convergence rates than its LS counterpart with a slightly higher computational cost. For both values of the Reynolds numbers, the RBF scheme exhibits a higher order accuracy than the LS scheme and the FEM. In terms of computational time, Table 1 reveals that the RBF scheme requires almost the triple of CPU time required for the LS scheme and about 1.5 times larger than the FEM.Table 1
Results for Example 1 obtained by the LS and RBF procedures using different meshes and two Reynolds numbers at timet=1. Computational results obtained using the finite element method are also included. The listed CPU times are in seconds.
Re=100Re=1000hL1−errorRateL2−errorRateCPUL1−errorRateL2−errorRateCPULS185.62809E-02—6.66197E-02—1.584.13981E-01—4.10072E-01—1.851161.92205E-021.552.43842E-021.654.171.38469E-011.581.37162E-011.504.961326.12442E-031.658.50242E-031.7115.04.35143E-021.674.46236E-021.6216.31641.79573E-031.772.69050E-031.8263.21.28475E-021.761.36396E-021.7185.6RBF (4)182.51041E-02—3.82963E-02—2.512.32246E-01—2.90329E-01—2.691166.23066E-032.011.07381E-021.836.116.10053E-021.938.10952E-021.877.161321.45335E-032.102.47027E-032.1219.231.47318E-022.051.93136E-022.0721.221642.87052E-042.345.05109E-042.29813.13500E-032.234.29168E-032.1789RBF181.55020E-02—1.96313E-02—3.312.01827E-01—2.84819E-01—4.251161.74416E-033.152.55042E-032.949.592.42007E-023.064.82984E-022.5610.61321.88486E-043.213.25501E-042.9737.51.96833E-033.625.55544E-033.1227.61641.98122E-053.254.01275E-053.021761.44283E-043.774.87645E-043.51239FEM184.24047E-02—5.09578E-02—3.013.10254E-01—3.21647E-01—3.761161.08899E-021.961.37588E-021.897.997.81031E-021.998.44095E-021.938.761322.74141E-031.993.55840E-031.9528.81.92569E-022.022.02081E-022.0220.91646.62007E-042.058.83457E-042.011254.39940E-032.134.68917E-032.15169Figure4 depicts the results obtained for the solution u by the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques on a uniform structured mesh with h=1/64. For comparison reason, we also include the analytical solution in the same figure. Note that since the solution for this example satisfies u+v=3/2, only the solution u is presented. For both considered Reynolds numbers, it is clear that the LS technique involves a numerical diffusion which is more pronounced than in those obtained using the RBF interpolation technique. At Re=1000, strong gradients are formed in the solution which has been well captured by the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique. The LS interpolation technique fails to correctly capture this sharp solution, and the numerical diffusion can clearly be seen in its results presented in Figure 4. To further demonstrate these features, one-dimensional cross-sections of the computed solutions at the off-diagonal line y=1−x are shown in Figure 5. The results computed using the FEM and the analytical solutions are also included in Figure 5 for comparison reasons. For Re=100, the physical diffusion dominates the flow solution, and both results computed using the FEM and RBF methods coincide with the analytical solution. However, the results obtained using the LS interpolation technique are less accurate than those obtained using the FEM and RBF methods. Increasing the Reynolds number to Re=1000, the results obtained using the FEM are oscillatory and do not preserve monotonicity, observe the overshoots and undershoots in the FEM results in Figure 5. From the same figure, we can see that the results obtained using the RBF method are monotone and free from any nonphysical oscillations. It is also evident that the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique preserves the high-order accuracy for selected values of the Reynolds number.Figure 4
Results for the solutionu in Example 1 at time t=1 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) on a mesh with h=1/64 and the Reynolds number Re=100 (1st row) and Re=1000 (2nd row).Figure 5
One-dimensional cross-sections of the computed solutionu at y=1−x for Example 1 at time t=1 on a mesh with h=1/64 using Re=100 (a) and Re=1000 (b).


## 4.2. Example 2

We solve the Burgers’ system (1) in a squared domain Ω=0,1×0,1 subject to initial and boundary conditions obtained from the following analytical solution
(47)ux,y,t=−4πexp−5π2t/Recos2πxsinπyRe2+exp−5π2t/Resin2πxsinπy,(48)vx,y,t=−2πexp−5π2t/Resin2πxcosπyRe2+exp−5π2t/Resin2πxsinπy.It is easy to verify that the pair solutionu,v satisfies the relation u=2v, and therefore, only results obtained for the solution u are presented for this problem. Compared to the previous test example, the solution of this problem is smother with a fast decay as the time progresses. This example has also been investigated in many references, see [9, 12] among others. Similarly, we present numerical results obtained for two different Reynolds numbers Re=100 and Re=1000 at time t=2 using a time step Δt=0.01. Table 2 summarizes the relative errors, convergence rates, and CPU times for the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques on different uniform structured meshes varying the spatial step h. In Table 2, we also include the computational results obtained using the finite element semi-Lagrangian method (FEM) and the RBF scheme (RBF (4)) using the same number of interpolant points as the LS scheme. A decrease in the relative L1- and L2-error is clearly seen in all the results as the mesh is refined with a faster convergence in the results obtained using the RBF interpolation technique compared to the other schemes. It is also worth remarking that higher convergence accuracy is achieved for Re=100 in all considered methods compared to the results for Re=1000. For this test example, the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique is the most accurate compared to the other methods. For instance, at Re=100 on the mesh with h=1/64, convergence rates in the L2-error for the LS, RBF (4), RBF, and FEM schemes are 1.52, 2.31, 2.55, and 2.12, respectively. At Re=1000, these rates become 1.45, 2.23, 2.49, and 2.01 for the LS, RBF (4), RBF, and FEM, respectively. Under the considered flow conditions, the computational time required for the simulations using the RBF procedure is higher than for the LS and FEM. Using the same number of points in the RBF (4) procedure as in the LS scheme, it also demonstrates the better convergence rates achieved in the RBF method for this example.Table 2
Results for Example 2 obtained by the LS and RBF procedures using different meshes and two Reynolds numbers at timet=2. Computational results obtained using the finite element method are also included. The listed CPU times are in seconds.
Re=100Re=1000hL1−errorRateL2−errorRateCPUL1−errorRateL2−errorRateCPULS181.22884E-02—1.54575E-02—3.561.44367E-02—2.00704E-02—4.551164.62427E-031.415.98037E-031.379.325.74246E-031.338.09481E-031.319.891321.65775E-031.482.25048E-031.4128.22.16095E-031.413.17553E-031.3530.61645.62228E-041.567.84711E-041.521157.64011E-041.501.16231E-031.45122RBF(4)185.51041E-03—6.83363E-03—4.516.52321E-03—7.13223E-03—5.691161.46628E-031.911.76865E-031.8310.411.73578E-031.931.93771E-031.8812.111323.71687E-042.104.36075E-042.1233.534.36963E-042.054.48865E-042.1135.221648.55054E-052.128.79391E-052.311389.84533E-052.159.56795E-052.23145RBF184.49127E-04—5.66690E-04—8.185.18864E-04—5.95316E-04—10.61168.80945E-052.351.19132E-042.2522.31.03912E-042.321.22564E-042.2824.11321.56814E-052.492.18024E-052.4570.51.91493E-052.442.52357E-052.4178.31642.69630E-062.543.72287E-062.552993.33855E-062.524.49211E-062.49324FEM189.34059E-03—1.05866E-02—7.369.88733E-03—9.91181E-03—9.631162.55534E-031.872.81701E-031.9118.22.74267E-041.852.74946E-031.8120.11326.52258E-041.967.29088E-041.9553.17.09847E-041.957.26557E-041.9260.21641.57509E-042.051.67724E-042.122071.73809E-042.031.80385E-042.01231Results obtained for the solutionu using the proposed cell-centered finite volume semi-Lagrangian method and the analytical solution are illustrated in Figure 6 for Re=100 and Re=1000 at time t=2 on the mesh with h=1/64. For both values of the Reynolds number, there is no visible difference between the results obtained using the RBF interpolation technique and the exact solutions. On the other hand, the results obtained using the LS interpolation technique are less accurate with excessive numerical diffusion particularly for the simulations using Re=1000. To further demonstrate these features, we illustrate in Figure 7 the one-dimensional cross-sections of the computed solutions at the diagonal y=x. For comparison, we depict in the same figure the computational results obtained using the finite element semi-Lagrangian method. As expected, a smooth behavior can be seen in the solution u for this test example at both Reynolds numbers Re=100 and Re=1000. Obviously, the LS interpolation introduces a noticeable numerical diffusion, and it becomes more intense at Re=1000, compare its failure to resolve sharp gradients in Figure 7. From the same figure, we can see that the results obtained using the RBF interpolation technique are slightly more accurate than those obtained using the FEM. From the results presented, it is clear that the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique performs very well for this test example.Figure 6
Results for the solutionu in Example 2 in a squared domain at time t=2 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) on a mesh with h=1/64 and the Reynolds number Re=100 (1st row) and Re=1000 (2nd row).Figure 7
One-dimensional cross-sections of the computed solutionu at y=x for Example 2 in a squared domain at time t=2 on a mesh with h=1/64 using Re=100 (a) and Re=1000 (b).Our next objective is to test and examine the performance of the present cell-centered finite volume semi-Lagrangian for solving the Burgers’ equation (1) in complex domains using unstructured meshes. To this end, we solve this example on the complex domain shown in Figure 8 for which the initial and boundary conditions are obtained from the analytical solution (47). Here, the computational domain is defined by its parametric curved boundary
(49)Γ=x,y∈2:xθ=Rθcosθ,yθ=Rθsinθ,where
(50)Rθ=2+12sin7θ,θ∈0,2π.Figure 8
Computational mesh used for Example 2 in a complex domain.Note that similar geometry has been considered in [15, 37] among others. An unstructured mesh with 4056 elements is used in our simulations, and the obtained results are presented for two different Reynolds numbers Re=100 and Re=1000 at time t=1 using a fixed time step Δt=0.01. In Figure 9, we display the numerical results for the solution u obtained by the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques along with the analytical solution. In contrast to the previous simulation in the squared domain for which one period in the results is obtained, the results in this complex domain exhibit several periods. As it can be seen in the results obtained using the LS approach, significant differences between these results and the exact solution are observed especially for the convection-dominated regime at Re=1000. Indeed, the computational treatment of such complex patterns on complex geometries often requires efficient numerical algorithms, and clearly, the LS approach is not suitable for these situations. However, the RBF method shows a very good behavior and matches the exact solution at both Re=100 and Re=1000. An examination of the one-dimensional cross-sections of the computed results at y=x in Figure 10 confirms these claims. For a comparison reason, we depict in the same figure the computational results of the finite element semi-Lagrangian approach. Again, the LS method produces results with large numerical diffusion, and it fails to capture the exact solutions. The FEM is in reasonable agreement with exact solution although it substantially manifests small distortion particularly in regions near sharp gradients. From the same figure, the accuracy in the results obtained using the RBF method is superior to those obtained using the FEM and LS methods. It should also be stressed that the complex nature of the geometry in this example reflects the ability of the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique to resolve the complex patterns.Figure 9
Results for the solutionu in Example 2 in a complex domain at time t=1 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) using Re=100 (1st row) and Re=1000 (2nd row).Figure 10
One-dimensional cross-sections of the computed solutionu at y=x for Example 2 in a complex domain at time t=1 using Re=100 (a) and Re=1000 (b).


## 4.3. Example 3

Our final test example consists in solving the Burgers’ equation (1) in the complex domain shown in Figure 11 defined by its parametric boundary
(51)Γ=x,y∈2:xθ=Rθcosθ,yθ=Rθsinθ,where
(52)Rθ=53cos3θ+8−sin3θ31/3,θ∈0,2π.Figure 11
Computational mesh used for Example 3 in a complex domain.A similar geometry has been considered in [15] for the coupled Burgers’ equations. To obtain the corresponding initial and boundary conditions, we use the following exact solution
(53)ux,y,t=121−tanhx4ν+y4ν−t4ν,vx,y,t=121−tanhx4ν+y4ν−t4ν.This test example has also been studied in [38, 39] for a squared domain. In our simulations, we consider an unstructured mesh with 4883 elements, and numerical results are presented at time t=1 for two different Reynolds numbers Re=10 and Re=1000 using a time step fixed to Δt=0.01. Note that to reduce the computational cost, the time steps Δt are chosen as large as possible in our simulations. This makes most explicit Eulerian-based methods noncompetitive, since they are subject to stability restriction conditions for the convection terms. Therefore, the criteria of choosing time steps in our simulations performed in this study was mainly based on accuracy considerations and stability condition (43) required for the solution of viscous parts.In Figure12, we present results obtained for the solution u using the proposed cell-centered finite volume semi-Lagrangian method using the LS and RBF interpolation techniques compared to the analytical results. The one-dimensional plots in Figure 13 correspond to cross-section solutions at the off-diagonal y=1−x of these results including those obtained using the FEM. A visual comparison of the results in these figures shows severe numerical dissipation in the solutions obtained using the LS method. This numerical dissipation is more pronounced for large values of Re, and it reduces as the mesh is refined. Under the considered flow conditions, the FEM and LS methods exhibit substantially large numerical diffusion particularly at the feet of areas where the gradient of the solution is very sharp. However, such numerical dissipation is completely absent in the results attained using the RBF method. It is also clear that the proposed cell-centered finite volume semi-Lagrangian method using the RBF interpolation technique performs best for this test example. It should be mentioned that the proposed cell-centered finite volume semi-Lagrangian method is typically designed to solve this class of convection-dominated flow problems using time steps ten to twenty times larger than its Eulerian counterparts. It should be noted that the proposed cell-centered finite volume semi-Lagrangian method performs very well as the solution remains stable and accurate even in the case of coarse meshes with no need for nonlinear solvers or small time steps in the simulations.Figure 12
Results for the solutionu in Example 3 at time t=1 obtained using the LS approach (1st column), RBF approach (2nd column), and exact solution (3rd column) using Re=10 (1st row) and Re=1000 (2nd row).Figure 13
One-dimensional cross-sections of the computed solutionu at y=1−x for Example 3 at time t=1 using Re=10 (a) and Re=1000 (b).


## 5. Conclusions

We have proposed a new cell-centered finite volume semi-Lagrangian method for the numerical solution of two-dimensional coupled Burgers’ problems on unstructured triangular meshes. The new method incorporates the combination of a semi-Lagrangian scheme with a cell-centered finite volume for space discretization and an explicit Runge-Kutta scheme for time integration. In this study, we have assessed two different interpolation techniques accounting for the element where the departure points are located. The proposed cell-centered finite volume semi-Lagrangian method is suitable for complex geometries, independent of the sizes, and arrangement of the mesh elements, and it can be performed using time steps larger than those required for its Eulerian counterparts. In addition, the proposed method is simple and stable and eliminates the numerical difficulties related to the discretization of nonlinear convective terms in the considered problems. A comparison to the conventional finite element semi-Lagrangian method shows the practicability of the current cell-centered finite volume semi-Lagrangian approach combined with the RBF interpolation, to solve the two-dimensional coupled Burgers’ problems. Verification of the proposed cell-centered finite volume semi-Lagrangian method has been carried out using several test problems with known analytical solutions. Comparison between two different interpolation methods reveals that the LS method is about three times faster than the RBF method. However, incorporating the RBF interpolation in the Lagrangian stage results in a highly accurate algorithm exhibiting good capturing of sharp gradients, high accuracy, and strong stability for all flow regimes considered. On the other hand, the LS method introduces excessive numerical diffusion and fails to capture steep gradients in the computed solutions. The presented results demonstrate the capability of the cell-centered finite volume semi-Lagrangian method that can provide insight to complex coupled Burgers’ features. The current study has been devoted to the numerical computations of two-dimensional coupled Burgers’ problems. However, our future interest consists on the use of the effective high-order cell-centered finite volume semi-Lagrangian method for coupled Burgers’ problems in three space dimensions using parallel processing. Our current effort is therefore to extend these techniques to coupled Burgers’ problems and incompressible Navier-Stokes equations in three space dimensions using unstructured meshes.



---

*Source: 8192192-2022-02-13.xml*