Numerical Analysis of Convection–Diffusion Using a Modified Upwind Approach in the Finite Volume Method

The main focus of this study was to develop a numerical scheme with new expressions for interface flux approximations based on the upwind approach in the finite volume method. Our new proposed numerical scheme is unconditionally stable with second-order accuracy in both space and time. The method is based on the second-order formulation for the temporal approximation, and an upwind approach of the finite volume method is used for spatial interface approximation. Some numerical experiments have been conducted to illustrate the performance of the new numerical scheme for a convection–diffusion problem. For the phenomena of convection dominance and diffusion dominance, we developed a comparative study of this new upwind finite volume method with an existing upwind form and central difference scheme of the finite volume method. The modified numerical scheme shows highly accurate results as compared to both numerical schemes.


Introduction
Partial differential equations (PDEs) have a vital role in describing different phenomena in real life. Such applications appear in many fields-in industry, biology, agriculture, petro-chemistry, heat transfer in draining films, the spread of solutes in liquid flowing through a tube, flow in porous media, meteorology and electric discharge [1][2][3][4][5][6]. Fluid dynamics is an important scientific area where governing equations are most often differential equations, and any analytical solutions for these equations are only possible for a restricted and limited number of cases. Due to this, several numerical techniques have been developed for the numerical approximations of convection-diffusion problems [7][8][9][10][11][12][13][14][15][16][17][18]. The most common numerical techniques are based on the finite difference method, the finite volume method and the finite element method. The comparison of finite difference and finite volume method for the numerical approximation of the convection-diffusion equation was studied in the papers [19][20][21][22]. The finite volume method is a special kind of numerical approach which is widely used for solving convection-diffusion problems in computational fluid dynamics [23][24][25][26][27] and also to simulate advection-diffusion problems with convection dominated processes [28][29][30]. The numerical simulation and the error estimate of steady-state convection-dominated diffusion problem was investigated in [31] using a cell-vertex finite volume approximation, whereas [32] used cell-centered finite difference approximations for a steady-state convection-diffusion problem to get second-order convergence that fulfilled the discrete maximum principle with discrete Sobolev spaces. By looking at the work in [33] a singular perturbed convection-diffusion problem was discretized with the help of the upwind finite volume scheme with uniform point-wise convergence over layer-adapted meshes. The authors [34] discussed the numerical solution of the unsteady-state advection-diffusion problem in two dimensions using the upwind finite volume method and an alternating-implicit, upwind finite volume method. Numerical simulations of non-Newtonian fluids were investigated using the finite volume formulation [35], where the approximation of the convection process is managed by using the QUICK differences scheme. A two-dimensional, incompressible in-viscid-free surface was simulated in [36] and they used the level set method for capturing the interface. A numerical solution of the convection-diffusion problem with variable coefficients over the nonstandard control volume grid was analyzed in [37] using the weighted upwind finite volume method. In the [38] finite volume method combined with the Lattice-Boltzmann method was studied for discretizing the convection operator over different regions, whereas [39] discussed the error estimation of the non-linear advection operator using the upwind finite volume method over a geometric corrector. In the meantime, the convergence analysis for the general initial and boundary condition with linear transport problems over the regularity property was solved in [40] using the upwind finite volume method.
Since in flow simulations it is imperative to satisfy conservation laws at all levels [41], by avoiding generation and consumption of mass, energy and momentum due to artificial terms inside the control volume, this is not guaranteed in finite difference and finite element methods. The control volume integration is a key step in the finite volume method, which ensures the conservation of relevant properties at each finite control cell level. This control volume integration gives a semi-discretize conservation law which involves fluxes at the interfaces of control volumes. The most popular forms of the finite volume method discretization are cell-centered, vertex-centered and face-centered [42]. The important contributions of the finite volume method are that it can also be formulated in the physical space on unstructured meshes and it is relatively easy to implement a variety of boundary conditions in a non-invasive manner [43].
The numerical schemes for time-dependent problems are known as explicit, implicit and fully implicit based on the chosen temporal approximation. One of the most common and widely used techniques for temporal approximation is the Crank-Nicolson method [44,45]. This technique forms a second-order implicit scheme in time. John Crank and Phyllis Nicolson constructed this well-known numerical technique in the mid 20th century [44,[46][47][48][49][50][51]. In the past, based on spatial interface approximations, several finite volume numerical schemes were developed to solve the convection-diffusion problem.
In this study, based on the upwind approach, new expressions were obtained for spatial approximation at the interface. Crank-Nicolson formulation is used for the approximation along the temporal direction, and using these new interface approximations, a numerical scheme was developed for the numerical approximation of convection-diffusion problem. This newly developed scheme is unconditionally stable and gives second-order accuracy. Some numerical experiments have been carried out to validate our newly developed theoretical algorithm for the upwind approach, and we have compared the numerical results of our modified upwind approach with the central differences scheme and the most common upwind approach [4,21], and the numerical approximation of the unsteady convection-diffusion problem.
The paper is organized as follows: we have introduced preliminary remarks and definitions in Section 2. We have also shown the discretized and formulated form of our problem in Section 3, and in Section 4, the numerical scheme is treated and analyzed with the boundary condition. In addition, our stability and convergence order analysis are discussed using the von Neumann-Fourier method in Section 5. In Section 6, numerical experiments are shown to verify our theoretical study and the conclusions are presented in Section 7.

Preliminary Remarks
Definition 1. The finite volume method is an integration-based numerical method, where after integration, we need to approximate the variables at interfaces and different choices of interface approximations give rise to different numerical schemes. The existing and most common techniques for interface approximation are central difference approximation (piece-wise linear profile assumption) and upwind approach approximation (stepwise profile assumption). Interface approximation using central difference approach are given as The upwind or donor cell approach takes into account the flow direction. This approach makes the interface approximation in terms of the upstream node, which means the stronger influencing node at the interface (see Figure 1). In the positive directional flow, the upstream node for west and east faces are respectively And for downwind or negative directional flow, the upstream nodes for interfaces are defined as In the theory of finite volume method, interface Γ is the common boundary of adjacent control volumes; it can be divided into two classes, internal faces (between two control volumes) and boundary faces (which coincide with the boundaries of the physical domain). Let the computational domain Ω be divided into n c sub-domains also called control volumes as The boundary of the control volume is express as: where n c f a denotes the number of faces or edges of a control volume. In our one dimensional study, each control volume has two edges, i.e., n c f a =2 (see Figure 2), whereas, for triangular and quadrilateral cells, there are 3 and 4 respectively. On a uniform and structured mesh, cells centers or computational nodes are indexed by subscript i or (i,j) for one and two-dimensional domains respectively. For an interface between the cell i and i + 1 the useful subscript notation is i + 1 2 (see in Figure 2).
Remark 1. In our study, by taking the transportiveness of the property or flow direction under consideration, we have obtained new expressions for the interface approximation, and for the treatment of boundary nodes we have adjusted imaginary or mirror nodes Φ 1−∆x,j and Φ n+∆x,j at a distance of ∆x 2 outside the left and right boundaries of computational domain shown (indicated in Figure 3) which are approximated by: ([40,52,53]). Let Ω = [0, L] be a bounded smooth domain and let T > 0; then the transportive linear equation with initial and boundary conditions with v > 0, k ≥ 0 being constant coefficients, has a finite and unique solution.

Discretization and Formulation of the Problem
The convection-diffusion equation can be either convection or diffusion dominated; it depends on the rate of convection and the rate of diffusion. The impact of convection on the transport of physical quantity is measured by Reynold number or the local Peclet number ( v∆x k ). For the discretization we have used the domains Ω = [0, L] and ∆x = L/M x , ∆t = T/N t , 0 ≤ i ≤ M x − 1, 0 ≤ j ≤ N t − 1. Consider the one-dimensional convection-diffusion equation with constant coefficients: and boundary condition where v > 0, k ≥ 0 represent the convective velocity and diffusion coefficient respectively, whereas Φ(x, t) is the variable convected in x-direction. v ∂Φ ∂x and k ∂ 2 Φ ∂x 2 are called convective (sometimes advective) and diffusive (sometimes viscous) terms respectively. For easy implementation of the new upwind approach of the finite volume method, we do not consider the source/sink term in our discretization of the Equation (1). The term convection means the movement of molecules within the fluid. In contrast, diffusion describes the spread of particles through the random motion from the region of high concentration to lower concentration region. In problems in which fluid flow plays important roles, we must consider the effects of convection dominance. A uniform discretization of the finite volume method is shown in Figure 2.
By integrating the problem over the control volume, Here, the first term will be approximated by using trapezoidal rule for numerical integration Here we need to make an approximation of variable and its gradients at the interface, as gradients will be approximated using central difference approximation as commonly used upwind approach considers the flow direction with the first-order accuracy. In the current study a numerical scheme based on upwind approach is proposed with second-order accuracy.
For an upwind approach, The upstream or high influential node for each interface of ith control volume is indicated in Figure 1. In the case of the upwind approach (v > 0) the flow is coming from the left side and approximates at its east face Φ i+ ∆x 2 ,j in terms of its upstream node Φ i,j will be: Here Φ i,j can be approximated by using the central difference approximation between the nodes Φ i−1,j and Φ i+1,j as follows: By putting this approximation in the above Equation (5), we have: This Equation (6), shows our expression for east interface Φ i+ ∆x 2 ,j approximation based on the upwind approach (v > 0). Now for the upwind approach (v > 0) the west face Φ i− ∆x 2 ,j will be approximated by its upstream node Φ i−1,j : Here Φ i−1,j can be approximated by central difference scheme approximation between the node Φ i−2,j and Φ i,j as follows: Substitute this approximation into (7) and it becomes: This Equation (8) shows our approximation at the west interface Φ i− ∆x 2 ,j of the ith control volume when (v > 0). Now substitute the approximations given in Equations (3), (4), (6) and (8) in above semi-discretized Equation (2); we get: After some simplification we get: After some rearrangements, the above expression becomes: Here θ method is used to evaluate the temporal integration: where For the Crank-Nicolson scheme θ= 1 2 and also divide by ∆x Let the coefficients be This Equation (13) shows our numerical scheme based on the upwind approach of the finite volume method, and it is valid for interior nodes from node 3 to n − 1 with the convergence of second-order in both space and time. For i = 1, 2 and i = n this scheme gives Φ 0,j , Φ −1,j and Φ n+1,j which we do not have in our computational domain; such cells are called ghost cells and values for these cells are called ghost values. For the special treatment of these boundary nodes (see in Section 4), our modified scheme holds the property of boundedness and also ensures the conservation of the variable using the consistent and same approximation expressions for a common interface of the two adjacent control volumes.

Remark 3.
For a downwind approach also we have define the upstream nodes (see Definition 2) and in the case of the downwind approach, i.e., v < 0, our interface approximations will be With these interface approximations for the downwind approach (v < 0) the numerical scheme will be

The Formulation of the Problem within the Analysis of Boundary Conditions
The boundary nodes 1, 2 and n need special treatment for the discretization of a governing equation over these nodes.
For node 1 the wet face Φ 1− ∆x 2 ,j is known because it coincides with the left boundary of the physical domain.
The east face of node 1 Φ 1+ ∆x 2 ,j can be approximate from Equation (6) by taking i = 1. For the approximation of east face of node 1 we have adjusted an imaginary node φ 1−∆x,j outside the boundary of the domain at a distance ∆x 2 (see in Remark 1 and Figure 3). For east face of node 1, Equation (6) becomes Here Φ 1−∆x,j is our left imaginary node; substitute its approximated value (see Remark 1) in This expression (16) shows the east face approximation of node 1. By substituting the approximations in expression (14) and (16) into Equation (2), the discretized equation for node 1 will be This is the discretized equation for node 1. As the uniqueness and major advantage of the finite volume method is to ensure the conservation of variable by using the consistent expression for a common interface between two adjacent control volumes, the west face of node 2 and east face of node 1 must have same approximation expression: The east face of node 2 can be approximated from expression (6) by taking i = 2.
Using the expressions (18) and (19) for west and east face approximations of node 2 respectively, the expression (2) after some simplification and applying the θ method gives the discretized equation for node 2 as This is the final discretized equation for node 2. The upstream node for east face Φ n+ ∆x 2 ,j is node n and can by approximated form expression (6) by taking i = n: Φ n+1 is our right side imaginary node (see Figure 3); put the approximation of imaginary node Φ n+1,j (see Remark 1) in Equation (21) and it becomes Φ R is known because it lies over the right boundary of physical domain. The west face Φ n− ∆x 2 ,j can be approximated using Equation (8) by taking i = n: By using the Equations (22) and (23) for node n the discretized equation will be

Stability Analysis
In this section, our objective is to find the stability region for the modified scheme of the finite volume method. The von Neumann technique based on Fourier mode analysis is used. Now the the stability of a scheme actually requires |ξ| ≤ 1.

Using this inequality, expression (27) becomes
From the Expression (29) it is clear that the magnitude |ξ j | is bounded, which indicates that the numerical scheme is unconditionally stable, as we expected. It means that there are no restrictions on the spatial and temporal step sizes ∆x and ∆t respectively, but we should choose ∆x and ∆t in such a way that the accuracy order of the scheme should not be degraded.
We can use the mathematical induction to prove the Theorem 2. Let j = 1 and assume |e 1 Suppose that if j ≤ n, e n i ∞ ≤ c(∆x 2 + ∆t 2 ); hold and assume j = n + 1; let |e n+1 which completes the proof.

Numerical Experiments
1. Let φ(x, t) be the property transported by means of convection and diffusion processes through one dimensional domain. Using equally spaced cells with the modified scheme of the finite volume method for the convection-diffusion equation, calculate the distribution of φ(x, t). The analytical solution of the problem given in [15] by φ(x, t) = 1 + e (−π 2 t/10) sin(π(x − t)) is considered on the finite domain [0, 1] × [0, 1]. The initial and boundary values are calculated from the exact solution by taking t = 0 and x = 0, x = 1, respectively.
For time step ∆t = 1/1500 and spatial step size ∆x = 1/320 the convection-diffusion problem is solved using a modified upwind approach of the finite volume method, obtaining numerical solution and maximum are shown in Figure 4.  The spatial and temporal convergence order of the newly constructed scheme are displayed in Tables 1 and 2, respectively. In Table 1 spatial convergence has been proved for temporal step sizes ∆t = 1/500, 1/1000 and 1/1500, whereas spatial step size ∆x taken as 1/10, 1/20, 1/40 and 1/80. In Table 2 results for the temporal convergence order have been shown. Temporal convergence order is proved for spatial step size (∆x) taken as 1/300, 1/800 and 1/1500,;the temporal step size (∆t) is 1/20, 1/40, 1/80 or 1/160. The obtained spatial and temporal convergence order mentioned in Tables 1 and 2 supports our theoretical approach about the convergence order of our scheme (see Sections 3 and 5.2). 2. Consider a homogeneous one dimensional convection-diffusion problem over a computational domain, 0 ≤ x ≤ 1, 0 ≤ t ≤ 1.; its exact solution is given in [16] with v ≥ 0 as the convection velocity to the x− direction and k > 0 as the diffusion coefficient. Now, we can find the transported property φ(x, t) for two different cases of convection dominance and diffusion dominance respectively, using the modified upwind finite volume method with equally spaced cells. The initial and boundary values are calculated from an exact solution by taking t = 0 and x = 0, x = 1, respectively. In Figure 5 numerical results for the convection dominant case (v = 1, k = 0.07) are expressed. The left figures show the numerical approximation and maximum error of the modified upwind finite volume method. In contrast, figures on the right side represent the numerical solution and maximum error using an existing upwind finite volume method for the convection-diffusion problem. In this case, the numerical value of maximum error produced by the new upwind scheme was 3.1930 × 10 −5 , but for the existing upwind scheme, the maximum error was 0.0010. This indicates new upwind scheme has the ability to deal with a dominant convection problem numerically as compared to the existing upwind scheme, and this shows consistency with our theoretical approach.  Tables 3 and 4 respectively show the maximum errors with spatial convergence order and temporal convergence order respectively, for the new scheme and existing upwind scheme of finite volume method. The spatial order of convergence in the case of convection dominance (k = 0.4 and v = 1) has been proven for temporal step sizes (∆t) 1/5000 and 1/10,000. For temporal convergence order, spatial discretization width (∆x) is taken as 1/300 and 1/400. The results shown in both tables support to our theoretical algorithm.  In Figure 6 numerical results for the diffusion dominant case (v = 0.05, k = 0.1) are expressed. Left figures show the numerical approximation and maximum error of the modified upwind finite volume method. In contrast, figures at the right side represent the numerical solution and maximum error using an existing upwind finite volume method for the convection-diffusion problem. In this case, the numerical value of maximum error produced by the new upwind scheme is 2.0803 × 10 −5 , but for the existing upwind scheme, it is 0.0012. This implies that for a dominant diffusion problem, the new upwind scheme gives a higher accuracy solution than the existing upwind FVM. Tables 5 and 6 respectively show the maximum error with spatial convergence order and temporal convergence order for the new numerical scheme and existing upwind scheme of finite volume method in the case of diffusion dominance. To prove the spatial and temporal convergence order in case of diffusion dominance, coefficients are taken as k = 0.7 and v = 0.05. For temporal convergence, the numbers of spatial nodes are taken as (M x ) = 300 and 400. For spatial convergence, time level Nt is fixed at Nt = 5000 and 10,000. The results shown in both tables support to our theoretical algorithm. From both tables it is clear that the new upwind scheme gives less maximum error; this implies that the numerical approximation is highly accurate as compared to the existing numerical scheme.  Figure 7 shows the numerical approximation (for numerical example 2) of the convection-diffusion problem using the modified upwind finite volume method with different time levels. Numerically, we have tested the modified scheme for convection dominance (see left of Figure 7) and for diffusion dominance (see right of Figure 7). Both figures indicate that the increase of the time level for a numerical solution using our proposed numerical scheme get closeness to an exact solution. 3. Consider a homogeneous one dimensional convection-diffusion problem over a computational domain, 0 ≤ x ≤ 1, 0 ≤ t ≤ 1, where its exact solution is given in [16] as: φ(x, t) = √ 20 √ t + 20 e (−(x−2−vt) 2 )/4k(t+20) with v ≥ 0 is the convection velocity to the x− direction and k > 0 is the diffusion coefficient. Now, we can find the transported property φ(x, t) for two different cases of convection dominance using new upwind finite volume method with equally spaced cells. The initial and boundary values are calculated from the exact solution by taking t = 0 and x = 0, x = 1, respectively. We took time step ∆t = 1/3000 and spatial step size ∆x = 1/320. This numerical experiment also shows a numerical comparison of the modified finite volume method and the most common form (central difference scheme) of finite volume method for the approximation of the convection-diffusion problem.
In Figures 8 and 9 numerical approximations using new form of finite volume method and the most central difference scheme of FVM are compared for different cases of convection dominance. When we compare the max-errors of both numerical schemes, in the first case (when v = 3 and k = 0.01) the max-errors produced by new FVM scheme and central (FVM) scheme are 8.0116 × 10 −7 and 1.6529 × 10 −4 respectively. For the second convection dominance case (v = 2, k = 0.001) the new scheme gives max-error 9.2800 × 10 −10 and the central scheme gives 3.9683 × 10 −8 . This implies that, for both cases, this newly constructed FVM scheme gives highly accurate results and it has the good capability to deal with convection dominance as compared to the central scheme of finite volume method. Figure 10 shows the numerical behavior of both numerical schemes for different values of the Peclet number(Pe = v∆x k ). Figure 10 left and right express the results for Pe = 6 and Pe = 14 respectively. From both figures, it is clear that the central difference scheme produces a solution that appears to oscillate about the exact solution. These oscillations are often called "wiggles" in the literature. The consistency with the analytical solution is weak as compared to the modified scheme of finite volume method.

Conclusions
In this article, an upwind approach of finite volume method is proposed to treat the convection-diffusion equation numerically with the Dirichlet-type boundary condition. With the consideration of convection dominance or flow direction, new expressions for interface approximations have been obtained for the directional flow phenomenon, and using those approximations a numerical scheme was developed with second-order accuracy in both space and time. The stability of this modified method has been discussed and it was shown that it is unconditionally stable; it has been numerically proven too. The convection diffusion equations for different values of convective velocity (v) and diffusion coefficient (k) have been numerically solved by this newly developed upwind approach, and also we have displayed the numerical results of convergence order for space and time which are consistent with our modified theoretical algorithm. The performance of this current numerical scheme for the convection-diffusion problem has been studied by measuring the maximum error and convergence order for convection and diffusion dominance phenomena. We showed a numerical comparison of our modified upwind approach of the finite volume method with an existing upwind approach of the finite volume method and an upwind approach of the finite difference method. The numerical results produced by these numerical methods indicate that our proposed numerical scheme gives numerical results which are highly accurate and near to the exact solution, so this approach is a good alternative to some existing ones for solving physical problems.