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

A cell-centered ﬁ nite volume semi-Lagrangian method is presented for the numerical solution of two-dimensional coupled Burgers ’ problems on unstructured triangular meshes. The method combines a modi ﬁ ed method of characteristics for the time integration and a cell-centered ﬁ nite 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 veri ﬁ ed using di ﬀ erent examples for coupled Burgers ’ problems with known analytical solutions. We also apply the method for simulation of an example of coupled Burgers ’ ﬂ ows in a complex geometry. In these test problems, the new cell-centered ﬁ nite volume semi-Lagrangian method demonstrates its ability to accurately resolve the two-dimensional coupled Burgers ’ problems.


Introduction
In this paper, given a two-dimensional bounded domain Ω in R 2 with Lipschitz boundary Γ and a time interval ½0, T, we are interested in solving the two-dimensional unsteady nonlinear coupled Burgers' equations ∂v ∂t 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 num-ber. 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][2][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][5][6][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 predictorcorrector 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][21][22][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 convectiondominated 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 cellcentered 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 2 Computational and Mathematical Methods 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 wellestablished 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 cellcentered finite volume semi-Lagrangian method for the convective part is presented in Section 2. 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.

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 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.

Approximation of Departure Points.
For the spacial discretization, we discretize the domain Ω = Ω ∪ Γ into a set of conforming triangles K i such as Ω 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 ½t n , t n+1 with a length Δt and denote by w n the value of a generic function w at time t n and w n i to denote the average value of w in the cell K i at time t n as where jK i j denotes the area of the control volume K i . Hence, the characteristic curves X i ðt n ; x i , t n+1 Þ associated with the advection problem (3) are calculated for each control volume K i , i = 1, ⋯, M by solving the backward differential equations ÞÞ Τ is the departure point defined at time τ of a particle that will reach x i = ðx i , y i Þ Τ the center of the control volume K i at time t n+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 t n of particles that will reach the points of a fixed mesh at time t n+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. Note that accurate approximation of the characteristic curves X i ðτ ; x i , t n+1 Þ is crucial to the overall accuracy of the cell-centered finite volume semi-Lagrangian method.

Computational and Mathematical Methods
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 thirdorder explicit Runge-Kutta method introduced in [30]. Thus, the procedure to approximate the solution of the ordinary differential equation (5) is achieved by It should be stressed that since the departure points X i ð t n ; x i , t n+1 Þ and the stages K ð1Þ i and K ð2Þ i 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 X i ðt n ; x i , t n+1 Þ and the stages K 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.

Procedures for Interpolation Stage.
Let us use the notation X n i ≔ X i ðt n ; x i , t n+1 Þ to denote the characteristics feet calculated in (7). Hence, the finite volume solution of (3) in the control volume K i at instant t n+1 is given by Since in general the departure point X n i does not coincide with a mesh point, the numerical solution uðX n i , t n Þ is computed by interpolation from known values at the control volumeK i where X n i resides and its neighboring elements. In the current work, we perform this step using the following interpolation techniques.

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 I n i the set of indices of control volumes which are close to the host control volumeK i where the departure point X n i belongs at time t = t n , and let U n be the vector of solutions u n j at mesh points x j with j ∈ I n i . Hence, the finite volume solution (10) is obtained using the LS interpolation as where β n k are the fitting coefficients, ϕ k ðx, 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 I n i . To solve the problem above, we use a linear regression to calculate the coefficients β n k . To this end, we define the merit function F It should be noted that if β n minimizes the functional (12), then the LS procedure precisely fits the vector β n = ðβ n 0 , β n 1 , ⋯, β n m Þ Τ for a given data set. Hence, we obtain the fitting coefficients β n k by solving the m normal equations Using equation (12), the normal equations yield which we can rewrite for each departure point X n i as a linear system of the form Figure 1: A schematic diagram showing departure points for the cell-centered finite volume semi-Lagrangian method. 4 Computational and Mathematical Methods where the matrix A is given by A = B Τ B with B denotes the matrix with entries ϕ l ðx j Þ, 1 ≤ l ≤ m, j ∈ I n i . The right-hand side vector r is given by r = B Τ f in which f refers to the vector with entries u n j , j ∈ I n i . Thereby, we define a linear approximation for the function u n ðx, yÞ as follows: To calculate the coefficients β n 0 , β n 1 , and β n 2 , we first evaluate (16) at the departure point X n i = ðX n i , Y n i Þ Τ to obtain β n 0 ; then, we solve the linear system (15) to obtain β n 1 and β n 2 . For instance, the inverse matrix A −1 is defined as where A similar procedure is used for the component v 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 volumeK 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].

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 where I n i denotes the set with dimension N i introduced previously, compare Figure 2 for an illustration. Notice that I n i 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 ðζ j Þ j∈I n i , γ 0 , γ 1 , and γ 2 are obtained by solving the linear system where ψ ij are the radial basis functions defined as By introducing the following matrices the linear system (20) can be rewritten in a more compact form as It should be stressed that the linear system (24) contains N i equations with N i + 3 unknowns. Consequently, to construct a nonsingular system, we should add the natural constraints Hence, the coefficients ðζ j Þ j∈I n i , γ 0 , γ 1 , and γ 2 in (19) are obtained by solving the block linear system ψ L 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).

Finite Volume Solution of Viscous Terms
Let us consider coupled Burgers' problem (1) reformulated using the total derivative as 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 equa-tion (27) over a control volume K i ∈ T h and using the divergence theorem, it yields ð where σ ij is the edge between K i and K j , j ∈ N i with N i is the set of all neighboring elements, i.e., control volumes which have a common one-dimensional face with K i . In (29), n ij denotes the unit outer normal vector to σ ij with respect to K j . Introducing the approximation where u i denotes the value of u at centroid of the cell K i , equation (29) becomes It is clear that equation (31) can be rewritten as where s ij = n ij · n σ , with n σ denotes the normal to the edge σ ij , see, for example, [34]. In (32), the flux function F ij is defined by 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 D h is reconstructed as where n ext is the outward normal vector to ∂D h . Following [34], the discretization of (34) yields the following approximation of the gradient over σ where n γ denotes the outward normal to the edge γ whereas P 1 ðγÞ and P 2 ðγÞ are the endpoints of an edge γ of ∂D h , see 6 Computational and Mathematical Methods Figure 3 for an illustration. Here, u E and u W are solutions at cell centers while u N and u S are the solutions at the vertices x N and x S expressed in terms of cell center values u K using the RBF interpolation technique discussed in the previous section. Thus, the reconstructed gradient can be written as follows where h σ = l WE · n σ and jσj = l SN · t σ , with t σ is the unit tangential vector to n σ . Here, l WE is the vector connecting the points W and E, and l SN is the vector connecting the points S and N, and Using the approximated gradient G σ , the numerical diffusive flux in (33) can be then written as Reformulating the semidiscrete system (31) in a compact form, the equations lead to a system of ordinary differential equations as where U is solution vectors with entries u i and RðUÞ = ½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 t n to the next time t n+1 can be achieved by whereŨ n is the solution vector with entries uðX n i , t n Þ 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] ν Δt Note that other time-stepping schemes can also be used for solving the system (39).

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 twodimensional coupled Burgers' problems. Since analytical solutions for these examples are known, the relative L 1 -error and L 2 -error at time t n are evaluated as where u n i and u exact ðx i , t n Þ are, respectively, the numerical and exact solutions at time t n in the control volume with centroid x i . 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.  Figure 3: Illustration of a covolume reconstructed from two adjacent control volumes used to discretize the viscous terms in the present study.

Computational and Mathematical Methods
This example has been widely used in the literature to examine numerical methods for solving the twodimensional Burgers' equations, see, for instance, [9][10][11][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 Table 1, 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 L 1 -error and L 2 -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,  Figure 4 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 cellcentered finite volume semi-Lagrangian method using the

Computational and Mathematical Methods
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, onedimensional 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.

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  It is easy to verify that the pair solution ðu, 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 L 1 -and L 2 -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 cellcentered 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 L 2 -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.
Results obtained for the solution u 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

Computational and Mathematical Methods
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.
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 where 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, signif-icant 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   13 Computational and Mathematical Methods 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 12: Results for the solution u 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). 14 Computational and Mathematical Methods

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 where 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 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 Figure 12, 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 onedimensional 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.

Conclusions
We have proposed a new cell-centered finite volume semi-Lagrangian method for the numerical solution of twodimensional 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.