# A Self-Adaptive Numerical Method to Solve Convection-Dominated Diffusion Problems

Authors: Zhi-Wei Cao; Zhi-Fan Liu; Zhi-Feng Liu; Xiao-Hong Wang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Convection-dominated diffusion problems usually develop multiscaled solutions and adaptive mesh is popular to approach high resolution numerical solutions. Most adaptive mesh methods involve complex adaptive operations that not only increase algorithmic complexity but also may introduce numerical dissipation. Hence, it is motivated in this paper to develop an adaptive mesh method which is free from complex adaptive operations. The method is developed based on a range-discrete mesh, which is uniformly distributed in the value domain and has a desirable property of self-adaptivity in the spatial domain. To solve the time-dependent problem, movement of mesh points is tracked according to the governing equation, while their values are fixed. Adaptivity of the mesh points is automatically achieved during the course of solving the discretized equation. Moreover, a singular point resulting from a nonlinear diffusive term can be maintained by treating it as a special boundary condition. Serval numerical tests are performed. Residual errors are found to be independent of the magnitude of diffusive term. The proposed method can serve as a fast and accuracy tool for assessment of propagation of steep fronts in various flow problems.


---

## Body

## 1. Introduction

It is well known that convection-diffusion equations arise in a variety of important science and engineering fields, for example, thermodynamics, fluid mechanics, and chemical reactions [1, 2]. In many applications recognized as convection-dominated problems, diffusion may be quite small in contrast to convection, so that the solution will develop steep moving fronts that are nearly shocks. Efficient and accurate assessment of propagations of moving fronts is important. However, standard Finite Difference Methods (FDM) tend to introduce spurious oscillations and stabilized methods often resort to upwind schemes, which introduce excessive numerical dissipations that will oversmooth the fronts [3, 4]. Consequently, an extremely fine mesh is needed to increase the resolution, making the convection-diffusion equation a most challenging one to be solved numerically.There are mainly two lines of research in numerical methods for achieving higher resolution with limited computation cost. One is the development of higher-order methods trying to reduce the numerical diffusion, such as high-order FDM with flux limiters [5–10] and the method of characteristics [3, 11–14]. The other important approach is using adaptive meshes. That is, mesh points are concentrated in the region where the solution varies steeply and is less concentrated in slowly varying regions.Different adaptive strategies, for example, Adaptive Mesh Refinement method and moving mesh methods, have shown success in solving convection-diffusion problems [15–22]. As an adaptive scheme is designed to yield highly nonuniform mesh, discretization of the governing equation on general meshes should be settled first [23]. As for how to locate mesh points adaptively, such methods often appeal to equidistributing a monitor function or solving mesh equations [19]. As a whole, adaptive methods generally consist of two parts: one is to locate mesh points adaptively and the other is to solve the governing equation on the nonuniform mesh. As a result, computation complexity is certainly increased. Moreover, a local interpolation is usually needed to project the information to new meshes. Such a remeshing operation often brings numerical diffusion. So, it is desirable to develop a numerical method that solves the problem adaptively but without complex adaptive operations.In this paper, we propose a self-adaptive numerical method to capture the steep fronts accurately. We first consider the following convection-diffusion equation to illustrate the main ideas of our method [1, 2]:(1)∂S∂t+∂fS∂x=∂∂xϵS∂S∂x,where the diffusive coefficient ϵS is positive. Instead of a traditional spatial-discretized mesh, the proposed method is based on a “range-discrete mesh” which is uniformly distributed in the value domain [24]. It is found that the range-discrete mesh has an adaptive nature because the density of mesh points is proportionate to the spatial gradient of the solution. During the computation process, positions of the mesh points are determined by the governing equation, while their values are fixed. In this way, movement of the mesh points not only solves the equation but also locates mesh points adaptively. In other words, the proposed method combines the two parts of an adaptive method, the adaptive strategy and solving the equation, in one procedure. Therefore, neither complex adaptive operations nor any interpolations are involved in the proposed method, leading to a simple and efficient method for convection-diffusion problems.It must be emphasized that, with a linear diffusive coefficient, the solution of (1) is smooth, though there exist steep gradients, which is the most concerned situation in the literature. However, with a nonlinear diffusive coefficient, though small, spatial gradient of the solution can be discontinuous or even infinite [25–27]. This singularity phenomenon is common in many applications but often smeared out by traditional methods. By simply treating it as a boundary condition, such a singular point can be maintained in our method.It should be noted that, for monotone solutions, the proposed method is similar to the method of Fayers and Sheldon [25]. Unfortunately, as the unknown is no longer single valued in their Lagrangian form, their method does not work for nonmonotone solutions. The proposed method overcomes this disadvantage by discretizing the governing equation over dynamic control volumes. Furthermore, how to deal with different kinds of boundary conditions with the proposed method is discussed.This paper is organized as follows. Section2 introduces the range-discrete strategy to obtain an adaptive mesh. Next, the numerical scheme based on the range-discrete mesh, including boundary conditions, will be illustrated in detail in Sections 3, 4, and 5. Serval numerical tests are conducted to validate the numerical method in Section 6 and conclusions are drawn in Section 7.


## 2. Range-Discrete Mesh

Considering that the unknown in (1) has bounded value, that is, Sx,t∈m,M, a series of discrete values can be obtained as m≤S0<S1<⋯<Sj<⋯<Sk≤M, where Sj+1=Sj+ΔSj. Then, one can get a series of intersection points of the line sets S=Sjj=0,1,…,k and the curve S=Sx,t, namely, Aii=0,1,… as graphed in Figure 1. Then, with the location of each intersection known, approximation of Sx,t is obtained. Thus, a mesh can be formed by treating each intersection as a mesh point.Figure 1
Description of the range-discrete mesh.For a constant piece in the solution curve, the governing equation is satisfied trivially. The solutions on both sides of this constant piece do not interact with each other until the constant piece disappears. Consequently, each endpoint of the constant piece can be treated as a boundary condition for each side, respectively. With two adjacent mesh points assigned to the same value and marked as boundary points, for example,A1 and A2 shown in Figure 1, a constant piece is sufficiently described.IfS=Sx is nonmonotone, serval intersections of the line S=Sj and the curve S=Sx exist, for example, points Aj1,Aj2,Aj3 shown in Figure 1. It is intuitive to employ all these intersections, assigned with the same discrete value Sj, as different mesh points to describe the nonmonotone feature. Moreover, a mesh point should also be set at the position of an extremum of S=Sx, for example, point Ap in Figure 1, because, without it, the solution between Aj2 and Aj3 will be identified as a constant piece.In contrast to spatial-discrete mesh used by FDM and many other methods, the discrete mesh obtained using the rules above is named as a “range-discrete” mesh. We should state that, thoughΔSj can be different from each other, we focus our attention on the equidistance situation for the sake of simplicity.If there exitn mesh points in the spatial area xa,xb, the density of mesh points can be defined as ρ=n/xb-xa. For a range-discrete mesh, variation of the solution over the area xa,xb can be defined as TVxa,xbS=∑j=1n-1Sxj+1-Sxj‍=n-1ΔS. Thus, we can obtain(2)ρ∝TVSxa,xbxa-xb,where TVSxa,xb/xa-xb can be interpreted as the spatial variation speed of the solution. Thus, we can conclude from (2) that ρ is proportionate to the spatial variation speed of the solution, so that the density of mesh points is large where the solution varies quickly and vice versa, indicating an adaptive nature, which can also be seen from Figure 1. As a special case, the density of mesh points for a constant piece is actually zero, being quite computational effort saving.In fact, aslimxa-xb→0⁡TVSxa,xb/xa-xb=Sxx, we can observe an equidistribution of the solution gradient Sxx in the range-discrete grid. This is an obvious adaptive criterion widely used in many adaptive mesh methods [19, 28, 29]. While complex adaptive operations have to be involved in these methods, it is done by the range-discrete mesh simply via a equidistance discretization in value domain.With the technique illustrated above, a range-discrete mesh can be obtained for general solutions. Next, the numerical scheme based on the range-discrete mesh for solving (1) will be present.


## 3. Numerical Scheme Based on the Range-Discrete Mesh

For the initial-boundary value problem of (1), one can get an initial range-discrete mesh using the technique illustrated in Section 2. Once a range-discrete mesh is obtained, a dynamic control volume can be formed around each mesh point. Take x1t,x2t in Figure 2 as an example, and let Sx1t,t≡S1=Si+Si-1/2 and Sx2t,t≡S2=Si+1+Si/2, making the control volume moving with the flow. Integrating (1) on this control volume and a small time interval t0,t1 where t1-t0=Δt, one can get(3)∫t0t1∫x1tx2t∂S∂tdxdt+∫t0t1∫x1tx2t∂fS∂xdxdt=∫t0t1∫x1tx2t∂∂xϵS∂S∂xdxdt,where the integration of the convective and diffusive term can be done directly as (4)∫t0t1∫x1tx2t∂fS∂xdxdt=∫t0t1fSx2t,t-fSx1t,tdt=∫t0t1fS2-fS1dt,∫t0t1∫x1tx2t∂∂xϵS∂S∂xdxdt=∫t0t1ϵSSxx2t-ϵSSxx1tdt=∫t0t1ϵS2Sx2-ϵS1Sx1dt.As for the transient term, according to Appendix A, if the solution curve Sx,t is monotone in x1t,x2t, we have (5)∫t0t1∫x1tx2t∂S∂tdxdt=-∫S1S2xS,t1dS+∫S1S2xS,t0dS.Consequently, the integration form of (1) is(6)∫S1S2xS,t1dS-∫S1S2xS,t0dS=∫t0t1fS2-ϵS2Sx2-fS1-ϵS1Sx1dt.The item on the left-hand side of (6) is the variation of the total quantity of the control volume over the time interval t0,t1; fS1-ϵS1Sx1 and fS2-ϵS2Sx2 can be interpreted as the left and right boundary flux, respectively. Thus, (6) is actually a natural result of mass conservation principle. It should be emphasized that (6) is exact as it is derived from (1) without any approximation.Figure 2
A dynamic control volume around a mesh point.Using linear approximation in the FDM manner, we can obtain first-order discrete form of (6) as(7)xit0+Δt-xit0S2-S1=fS2-ϵS2Sx2Δt-fS1-ϵS1Sx1Δt,where (8)S2=Si+1+Si2,S1=Si+Si-12,Sx2=Si+1-Sixi+1-xi,Sx1=Si-Si-1xi-xi-1.If the items on the right-hand side of (7) are the values at t, the discrete equation is explicit; and if they are values at t+Δt, it is implicit. One may notice that (7) is quite similar to the discrete equations of FDM and can be solved in the same way of FDM, for example, Newton’s method. Generally speaking, other kinds of discrete equations and higher-order forms can also be obtained in the same manner as FDM.However, we should emphasize that the unknown in (7) is xi, not Si, which is the essential difference between the proposed method and FDM. It means that the movements of the mesh points are tracked while their values are fixed throughout the computation process. When new positions of the mesh points are located, the approximation of the solution is obtained. Moreover, adaptivity is obtained simultaneously as a benefit of the range-discrete strategy.We should note that though derived in different manners, (7) is similar to that in the method of Fayers and Sheldon [25]. And also, it only fits for monotone solutions in x1t,x2t. To make the range-discrete grid work for nonmonotone solutions, adjustments are needed.


## 4. Adjustments for the Extremum Point

As has been mentioned in Section2, a mesh point, for example, Ap in Figure 1 or Figure 3, is set at the extremum if the solution is nonmonotone. For the control volume of such an extremum point, (6) does not fit and adjustments are needed.Figure 3
The dynamic control volume around an extremum point.The dynamic control volumex1t,x2t for an extremum mesh point Ap is shown in Figure 3, where S1=S2=Sp-1+ΔS/2. According to Appendix A, The items inside the brackets in (3) for the control volume of an extremum point are (9)∫t0t1∫x1tx2t∂S∂tdxdt=∫x1t1x2t1Sx,t+Δt-S2dx-∫x1t0x2t0Sx,t-S2dx. So, we have(10)∫x1t1x2t1Sx,t1-S2dx-∫x1t0x2t0Sx,t0-S2dx=∫t0t1-ϵS2Sx2--ϵS1Sx1dt.With an parabolic approximation of the solution curve betweenx1t,x2t, we obtain (11)Sx,t=Sp-4Sp-Sp-1xp+1-xp-12x-xp+1+xp-122.Consequently, we have (12)∫x1x2Sx,t-S2dx=23Sp-S1Sp-1-S1Sp-S1xp+1-xp-1,ϵS2Sx2-ϵS1Sx1=-8ϵS1Sp-S1Sp-Sp-1Sp-Sp-1xp+1-xp-1.Thus, the discrete equation for an extremum point is obtained.One can conclude from (10) that, with a positive ϵS, left-hand side of (10) always decreases as time goes on. It means that the value of the extremum Sp is approaching S1, making the extremum point a special case: its value and position are always changing and the control volume is shrinking.To avoid too small a control volume, the mesh pointsAp-1 and Ap+1 will be removed from the range-discrete mesh if Sp-S1 is less than a threshold value, for example, Sp-S1<ΔS/4. Then, let S1=Sp-2+ΔS/2 and a new control volume is formed immediately to continue the computation. In this way, the extremum point in the initial mesh is always an extremum point and represents the movement of the extremum in the solution.


## 5. Boundary Conditions

### 5.1. Dirichlet Boundary Condition

A Dirichlet boundary condition, that is,Sxd,t=Sd, can be used directly in the discrete equation in the same way as FDM. That is, one can consider a discrete point Ak, set at the boundary to provide the boundary information, with its position and value always fixed as xk=xd and Sk=Sd during the computing process.


### 5.2. Moving Boundary Condition

As has been mentioned in Section2, either endpoint of a constant piece in the solution is treated as a boundary condition. The value of this boundary is fixed, which is quite similar to the Dirichlet boundary condition, but the position of the boundary varies. Thus, we name it a moving boundary condition.At a moving boundary, we haveSλt,t=Sf and ϵSSxλt=0 (see Appendix B). Thus, it can be involved in the range-discrete mesh by a boundary point Ak with Sk=Sf and ϵSSxAk=0. The boundary condition is not involved in a discrete manner as what was done in [25], but in a more precise manner, making it more promising for high precision.One can see that the position of the moving boundary point is not really needed in the range-discrete equation. However, a constant piece in the solution might disappear and the solution on two sides will then interact with each other. To illustrate this, we need to merge the two boundary points into a normal mesh point if their distance is less than a thresholddmin. The position of the boundary point is only needed here and can be interpolated using Sx=a0x-λn by two nearby mesh points, where n can be settled using the method illustrated in Appendix B.As a benefit of such a special boundary condition, the singular point, where the gradient of the solution might be discontinuous or infinite, as indicated in AppendixB, can be maintained by the proposed method. Numerical tests in Section 6 also confirm this.


### 5.3. Neumann Boundary Condition

Given a Neumann boundary conditionSxxb=qb, a boundary mesh point Ak will also be set. The discrete value of Ak is determined dynamically both by the value of Ai, which is the point nearest the boundary, and by the sign of qb, that is, Sk=Si-sgn⁡qbΔS for the left boundary and Sk=Si+sgn⁡qbΔS for the right boundary. The position of Ak can be calculated by an linear interpolation:(13)qb=Sk-Sixk-xior other higher-order interpolations. As we can see, the position of a Neumann boundary point determined by (13) can be either in or out of the calculation area. So, we design the following rules: (1)
Ifxk is out the calculation area, computing continuous on.(2)
Ifxk is in the calculation area, this boundary point becomes a normal mesh point and a new boundary point will be inserted into the range-discrete mesh.(3)
Ifxi is out of the calculation area, Ai will be removed immediately and the information of Ak is determined by the new nearest boundary mesh point.


## 5.1. Dirichlet Boundary Condition

A Dirichlet boundary condition, that is,Sxd,t=Sd, can be used directly in the discrete equation in the same way as FDM. That is, one can consider a discrete point Ak, set at the boundary to provide the boundary information, with its position and value always fixed as xk=xd and Sk=Sd during the computing process.


## 5.2. Moving Boundary Condition

As has been mentioned in Section2, either endpoint of a constant piece in the solution is treated as a boundary condition. The value of this boundary is fixed, which is quite similar to the Dirichlet boundary condition, but the position of the boundary varies. Thus, we name it a moving boundary condition.At a moving boundary, we haveSλt,t=Sf and ϵSSxλt=0 (see Appendix B). Thus, it can be involved in the range-discrete mesh by a boundary point Ak with Sk=Sf and ϵSSxAk=0. The boundary condition is not involved in a discrete manner as what was done in [25], but in a more precise manner, making it more promising for high precision.One can see that the position of the moving boundary point is not really needed in the range-discrete equation. However, a constant piece in the solution might disappear and the solution on two sides will then interact with each other. To illustrate this, we need to merge the two boundary points into a normal mesh point if their distance is less than a thresholddmin. The position of the boundary point is only needed here and can be interpolated using Sx=a0x-λn by two nearby mesh points, where n can be settled using the method illustrated in Appendix B.As a benefit of such a special boundary condition, the singular point, where the gradient of the solution might be discontinuous or infinite, as indicated in AppendixB, can be maintained by the proposed method. Numerical tests in Section 6 also confirm this.


## 5.3. Neumann Boundary Condition

Given a Neumann boundary conditionSxxb=qb, a boundary mesh point Ak will also be set. The discrete value of Ak is determined dynamically both by the value of Ai, which is the point nearest the boundary, and by the sign of qb, that is, Sk=Si-sgn⁡qbΔS for the left boundary and Sk=Si+sgn⁡qbΔS for the right boundary. The position of Ak can be calculated by an linear interpolation:(13)qb=Sk-Sixk-xior other higher-order interpolations. As we can see, the position of a Neumann boundary point determined by (13) can be either in or out of the calculation area. So, we design the following rules: (1)
Ifxk is out the calculation area, computing continuous on.(2)
Ifxk is in the calculation area, this boundary point becomes a normal mesh point and a new boundary point will be inserted into the range-discrete mesh.(3)
Ifxi is out of the calculation area, Ai will be removed immediately and the information of Ak is determined by the new nearest boundary mesh point.


## 6. Numerical Examples

### 6.1. A Travelling Wave Solution

WithfS=S2/2 and ϵS being a constant, we have Burgers’ equation from (1):(14)St+12S2x=ϵSxx.It is easy to verify that this equation has the travelling wave solution:(15)Sx,t=2cecct-x/ϵ1+ecct-x/ϵ,where c is the speed of the wave. With c=1, (15) is the solution of (14) with an initial condition Sx,0=2e-x/ϵ/1+e-x/ϵ. Figure 4 and Table 1 show the solutions and error tables at different time using the proposed method when ϵ=10-1,10-2,10-3. It can be seen from Figure 4 that mesh points are located adaptively according to the spatial gradient. Moreover, the adaptivity follows the movement of the travelling front. As a benefit, with only a few mesh points, the front is solved with good accuracy.Table 1
Error table of the numerical solution of Burgers’ equation with different diffusive coefficients at different time.N is the number of mesh points.
ϵ
t
errors
N
=
10
20
40
80
160
1
0
-
1
0.5
L
∞ error
1.00
E
-
02
5.23
E
-
03
2.70
E
-
03
1.37
E
-
03
6.95
E
-
04
L
1 error
4.43
E
-
03
2.39
E
-
03
1.35
E
-
03
7.61
E
-
04
4.21
E
-
04
1.0
L
∞ error
1.13
E
-
02
5.83
E
-
03
2.96
E
-
03
1.50
E
-
03
7.52
E
-
04
L
1 error
5.43
E
-
03
2.77
E
-
03
1.39
E
-
03
7.17
E
-
04
3.74
E
-
04
1
0
-
2
0.5
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
04
2.94
E
-
04
1.50
E
-
04
7.59
E
-
05
3.82
E
-
05
1.0
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
04
2.94
E
-
04
1.50
E
-
04
7.59
E
-
05
3.82
E
-
05
1
0
-
3
0.5
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
05
2.94
E
-
05
1.50
E
-
05
7.59
E
-
06
3.82
E
-
06
1.0
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
05
2.94
E
-
05
1.50
E
-
05
7.59
E
-
06
3.82
E
-
06Figure 4
Numerical solutions att=0.5 and t=1 of the travelling wave with different diffusive coefficients. Solid line represents exact solutions and squares, circles, and triangles in S-x figures are numerical solutions. The t-x figures below each S-x figure show the spatial distribution of the mesh points at different time.
(a)
ϵ
=
10
-
1
(b)
ϵ
=
10
-
2
(c)
ϵ
=
10
-
3One should notice that, with the diffusive coefficientϵ decreasing, the front becomes steeper, tending to form a moving shock. While it is hard to maintain the accuracy near the front with the same number of mesh points using FDM, the accuracy of the proposed method for different diffusive coefficients is uniform.Furthermore, the calculating area of this problem is actually infinite, that is,x∈-∞,+∞. While most numerical methods approach this within a finite area, the proposed method solves this problem in an infinite area simply with two moving boundary conditions: S-∞=1, Sx-∞=0 and S+∞=0, Sx+∞=0.


### 6.2. A Nonmonotone Case

In this numerical example, the solution of Burgers’ equation is designed to be nonmonotone to test the numerical method. The initial condition and boundary condition are set to beSx,0=sin⁡πx, S0,t=0 and S1,t=0, respectively. According to Caldwell and Smith [30], the solution is (16)Sx,t=2πϵ∑n=1∞ane-n2π2ϵtnsin⁡nπxa0+∑n=1∞ane-n2π2ϵtcos⁡nπx,where (17)a0=∫01exp⁡-2πϵ-11-cos⁡πxdx,an=2∫01exp⁡-2πϵ-11-cos⁡πxcos⁡nπxdxn=1,2,3,….The numerical solutions at different times with diffusive coefficientsϵ=10-1,10-2,10-3 are given by Figure 5 and the errors are Table 2. The results also show good performance of the proposed method for nonmonotone solutions. Again, the adaptivity functions well and tracks the front in an accurate manner.Table 2
Error table of the numerical solution of the burgers’ equation with different diffusive coefficients at different time.N is the number of mesh points.
ϵ
t
errors
N
=
10
20
40
80
160
1
0
-
1
0.5
L
∞ error
7.95
E
-
03
2.94
E
-
03
1.18
E
-
03
3.64
E
-
04
1.72
E
-
04
L
1 error
2.83
E
-
03
8.55
E
-
04
2.60
E
-
04
1.04
E
-
04
7.38
E
-
05
1.0
L
∞ error
1.46
E
-
02
5.16
E
-
03
1.50
E
-
03
5.35
E
-
04
1.26
E
-
04
L
1 error
5.41
E
-
03
1.60
E
-
03
4.46
E
-
04
1.11
E
-
04
6.53
E
-
05
1
0
-
2
0.5
L
∞ error
3.65
E
-
02
7.29
E
-
03
2.39
E
-
03
6.99
E
-
04
4.64
E
-
04
L
1 error
5.83
E
-
03
1.18
E
-
03
3.28
E
-
04
9.05
E
-
05
4.18
E
-
05
1.0
L
∞ error
6.16
E
-
02
1.97
E
-
02
9.19
E
-
03
2.93
E
-
03
9.49
E
-
04
L
1 error
7.38
E
-
03
1.66
E
-
03
6.63
E
-
04
1.95
E
-
04
5.95
E
-
05
1
0
-
3
0.5
L
∞ error
2.42
E
-
02
1.10
E
-
02
4.52
E
-
03
1.52
E
-
03
1.49
E
-
03
L
1 error
5.00
E
-
03
1.46
E
-
03
4.02
E
-
04
9.82
E
-
05
9.71
E
-
05
1.0
L
∞ error
9.05
E
-
02
4.77
E
-
02
1.68
E
-
02
7.01
E
-
03
2.94
E
-
03
L
1 error
3.31
E
-
03
1.70
E
-
03
2.07
E
-
04
1.13
E
-
04
6.32
E
-
05Figure 5
Numerical solutions obtained by the proposed method att=0.5 and t=1 of Burgers’ equation with different diffusive coefficients. Solid line represents exact solutions and squares, circles, and triangles in S-x figures are numerical solutions. The t-x figures below each S-x figure show the spatial distribution of the mesh points at different time.
(a)
ϵ
=
10
-
1
(b)
ϵ
=
10
-
2
(c)
ϵ
=
10
-
3


### 6.3. Nonlinear Diffusive Term

How a nonlinear diffusive term effects the distribution of the front is tested here. Consider Burgers’ equation with a nonlinear diffusive term:(18)St+12S2x=ϵSSxxwith two Dirichlet boundary conditions S0,t=1, S1,t=0 and a lamp initial condition, which is S0≤x≤0.2,t=0=1-5x and S0.2<x≤1,t=0=0. Diffusive coefficients set as ϵS=0.1Sα with α=0.5,1,2 are tested. Solutions at t=0.5 are shown in Figure 6.Figure 6
Numerical solutions of different nonlinear diffusive coefficients: (a)Sx and (b) dSdx.
(a)
(b)It is shown by Figure6 that the gradient of the solution with α=0.5 tends to be zero at the front and thus a smooth solution. However, the gradient with α=2 tends to be infinite at the front, resulting in a continuous but not smooth solution. As a critical condition, the gradient tends to be a limited value at the front with α=1, which is also a conclusion of Appendix B.The solutions with diffusive coefficients set asϵS=0.01Sα and ϵS=0.001Sα with α=0.5,1,2 are shown in Figure 7. With a degenerating diffusive coefficient, the nonlinear term has less impact on the solution, approaching the shock solution of the corresponding hyperbolic equation. With such a small diffusive coefficient, not only its nonlinearity but also the whole diffusive term can be neglected. If the diffusive term is not so small in contract to the convective term, its nonlinearity does have a significant influence on the solution.Figure 7
Numerical solutions of different nonlinear diffusive coefficients: (a)ϵS=0.01Sα and (b) ϵS=0.001Sα.
(a)
(b)


### 6.4. Buckley-Leverett Equation

In this example, a two-phase flow problem in porous medium is solved. Considering an infinite horizontal layer, the regionx<0 is full of the wetting phase, for example, water, and the region x>0 is occupied by a nonwetting phase, for example, oil, initially. The governing equation can be obtained by Buckley-Leverett theory [31]: (19)ϕ∂S∂t+∂∂xVtfS+KKrnμnfSdPcdS∂S∂x=0,where S is the saturation of the wetting phase, Vt is the total infiltration rate, PcS represents the capillary pressure, and fS=μoKrw/μnKrw+μwKrn is the flux coefficient of the wetting phase. In this example, Vt=1 and porosity and absolute permeability of the layer are nondimensionalized to be ϕ=1 and K=1, respectively. Relative permeability of the wetting phase (w) and the nonwetting phase (n) are Krw=S4 and Krn=1+S1-S3, respectively, and their viscosity is μw=μn=1. The capillary pressure is the J-function Pc=0.1S-0.5.The convective term is typically nonconvex. Another problem is that the diffusive term is nonlinear, which will cause infinite saturation gradients at the front of each phase. Most methods will smear this singularity out because of their numerical diffusion, no matter how small. With the range-discrete mesh method, this singular feature can be maintained by dealing the front as a moving boundary condition.The numerical solutions at different time withΔS=0.05 are shown in Figure 8. It is obvious that there exists not only a steep front where saturation varies quickly but also two singular points at S→0 and S→1, where the saturation gradient tends to be infinite. As a benefit of the range-discrete mesh and the moving boundary condition, the steep front is described preciously and two singular points are maintained intactly with only 20 mesh points.Figure 8
Numerical solutions of the Buckley-Leveret equation att=0.25 and t=0.5. The exact solution shown in solid is obtained by FDM with a fine mesh.


## 6.1. A Travelling Wave Solution

WithfS=S2/2 and ϵS being a constant, we have Burgers’ equation from (1):(14)St+12S2x=ϵSxx.It is easy to verify that this equation has the travelling wave solution:(15)Sx,t=2cecct-x/ϵ1+ecct-x/ϵ,where c is the speed of the wave. With c=1, (15) is the solution of (14) with an initial condition Sx,0=2e-x/ϵ/1+e-x/ϵ. Figure 4 and Table 1 show the solutions and error tables at different time using the proposed method when ϵ=10-1,10-2,10-3. It can be seen from Figure 4 that mesh points are located adaptively according to the spatial gradient. Moreover, the adaptivity follows the movement of the travelling front. As a benefit, with only a few mesh points, the front is solved with good accuracy.Table 1
Error table of the numerical solution of Burgers’ equation with different diffusive coefficients at different time.N is the number of mesh points.
ϵ
t
errors
N
=
10
20
40
80
160
1
0
-
1
0.5
L
∞ error
1.00
E
-
02
5.23
E
-
03
2.70
E
-
03
1.37
E
-
03
6.95
E
-
04
L
1 error
4.43
E
-
03
2.39
E
-
03
1.35
E
-
03
7.61
E
-
04
4.21
E
-
04
1.0
L
∞ error
1.13
E
-
02
5.83
E
-
03
2.96
E
-
03
1.50
E
-
03
7.52
E
-
04
L
1 error
5.43
E
-
03
2.77
E
-
03
1.39
E
-
03
7.17
E
-
04
3.74
E
-
04
1
0
-
2
0.5
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
04
2.94
E
-
04
1.50
E
-
04
7.59
E
-
05
3.82
E
-
05
1.0
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
04
2.94
E
-
04
1.50
E
-
04
7.59
E
-
05
3.82
E
-
05
1
0
-
3
0.5
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
05
2.94
E
-
05
1.50
E
-
05
7.59
E
-
06
3.82
E
-
06
1.0
L
∞ error
1.16
E
-
02
5.97
E
-
03
3.03
E
-
03
1.53
E
-
03
7.66
E
-
04
L
1 error
5.64
E
-
05
2.94
E
-
05
1.50
E
-
05
7.59
E
-
06
3.82
E
-
06Figure 4
Numerical solutions att=0.5 and t=1 of the travelling wave with different diffusive coefficients. Solid line represents exact solutions and squares, circles, and triangles in S-x figures are numerical solutions. The t-x figures below each S-x figure show the spatial distribution of the mesh points at different time.
(a)
ϵ
=
10
-
1
(b)
ϵ
=
10
-
2
(c)
ϵ
=
10
-
3One should notice that, with the diffusive coefficientϵ decreasing, the front becomes steeper, tending to form a moving shock. While it is hard to maintain the accuracy near the front with the same number of mesh points using FDM, the accuracy of the proposed method for different diffusive coefficients is uniform.Furthermore, the calculating area of this problem is actually infinite, that is,x∈-∞,+∞. While most numerical methods approach this within a finite area, the proposed method solves this problem in an infinite area simply with two moving boundary conditions: S-∞=1, Sx-∞=0 and S+∞=0, Sx+∞=0.


## 6.2. A Nonmonotone Case

In this numerical example, the solution of Burgers’ equation is designed to be nonmonotone to test the numerical method. The initial condition and boundary condition are set to beSx,0=sin⁡πx, S0,t=0 and S1,t=0, respectively. According to Caldwell and Smith [30], the solution is (16)Sx,t=2πϵ∑n=1∞ane-n2π2ϵtnsin⁡nπxa0+∑n=1∞ane-n2π2ϵtcos⁡nπx,where (17)a0=∫01exp⁡-2πϵ-11-cos⁡πxdx,an=2∫01exp⁡-2πϵ-11-cos⁡πxcos⁡nπxdxn=1,2,3,….The numerical solutions at different times with diffusive coefficientsϵ=10-1,10-2,10-3 are given by Figure 5 and the errors are Table 2. The results also show good performance of the proposed method for nonmonotone solutions. Again, the adaptivity functions well and tracks the front in an accurate manner.Table 2
Error table of the numerical solution of the burgers’ equation with different diffusive coefficients at different time.N is the number of mesh points.
ϵ
t
errors
N
=
10
20
40
80
160
1
0
-
1
0.5
L
∞ error
7.95
E
-
03
2.94
E
-
03
1.18
E
-
03
3.64
E
-
04
1.72
E
-
04
L
1 error
2.83
E
-
03
8.55
E
-
04
2.60
E
-
04
1.04
E
-
04
7.38
E
-
05
1.0
L
∞ error
1.46
E
-
02
5.16
E
-
03
1.50
E
-
03
5.35
E
-
04
1.26
E
-
04
L
1 error
5.41
E
-
03
1.60
E
-
03
4.46
E
-
04
1.11
E
-
04
6.53
E
-
05
1
0
-
2
0.5
L
∞ error
3.65
E
-
02
7.29
E
-
03
2.39
E
-
03
6.99
E
-
04
4.64
E
-
04
L
1 error
5.83
E
-
03
1.18
E
-
03
3.28
E
-
04
9.05
E
-
05
4.18
E
-
05
1.0
L
∞ error
6.16
E
-
02
1.97
E
-
02
9.19
E
-
03
2.93
E
-
03
9.49
E
-
04
L
1 error
7.38
E
-
03
1.66
E
-
03
6.63
E
-
04
1.95
E
-
04
5.95
E
-
05
1
0
-
3
0.5
L
∞ error
2.42
E
-
02
1.10
E
-
02
4.52
E
-
03
1.52
E
-
03
1.49
E
-
03
L
1 error
5.00
E
-
03
1.46
E
-
03
4.02
E
-
04
9.82
E
-
05
9.71
E
-
05
1.0
L
∞ error
9.05
E
-
02
4.77
E
-
02
1.68
E
-
02
7.01
E
-
03
2.94
E
-
03
L
1 error
3.31
E
-
03
1.70
E
-
03
2.07
E
-
04
1.13
E
-
04
6.32
E
-
05Figure 5
Numerical solutions obtained by the proposed method att=0.5 and t=1 of Burgers’ equation with different diffusive coefficients. Solid line represents exact solutions and squares, circles, and triangles in S-x figures are numerical solutions. The t-x figures below each S-x figure show the spatial distribution of the mesh points at different time.
(a)
ϵ
=
10
-
1
(b)
ϵ
=
10
-
2
(c)
ϵ
=
10
-
3


## 6.3. Nonlinear Diffusive Term

How a nonlinear diffusive term effects the distribution of the front is tested here. Consider Burgers’ equation with a nonlinear diffusive term:(18)St+12S2x=ϵSSxxwith two Dirichlet boundary conditions S0,t=1, S1,t=0 and a lamp initial condition, which is S0≤x≤0.2,t=0=1-5x and S0.2<x≤1,t=0=0. Diffusive coefficients set as ϵS=0.1Sα with α=0.5,1,2 are tested. Solutions at t=0.5 are shown in Figure 6.Figure 6
Numerical solutions of different nonlinear diffusive coefficients: (a)Sx and (b) dSdx.
(a)
(b)It is shown by Figure6 that the gradient of the solution with α=0.5 tends to be zero at the front and thus a smooth solution. However, the gradient with α=2 tends to be infinite at the front, resulting in a continuous but not smooth solution. As a critical condition, the gradient tends to be a limited value at the front with α=1, which is also a conclusion of Appendix B.The solutions with diffusive coefficients set asϵS=0.01Sα and ϵS=0.001Sα with α=0.5,1,2 are shown in Figure 7. With a degenerating diffusive coefficient, the nonlinear term has less impact on the solution, approaching the shock solution of the corresponding hyperbolic equation. With such a small diffusive coefficient, not only its nonlinearity but also the whole diffusive term can be neglected. If the diffusive term is not so small in contract to the convective term, its nonlinearity does have a significant influence on the solution.Figure 7
Numerical solutions of different nonlinear diffusive coefficients: (a)ϵS=0.01Sα and (b) ϵS=0.001Sα.
(a)
(b)


## 6.4. Buckley-Leverett Equation

In this example, a two-phase flow problem in porous medium is solved. Considering an infinite horizontal layer, the regionx<0 is full of the wetting phase, for example, water, and the region x>0 is occupied by a nonwetting phase, for example, oil, initially. The governing equation can be obtained by Buckley-Leverett theory [31]: (19)ϕ∂S∂t+∂∂xVtfS+KKrnμnfSdPcdS∂S∂x=0,where S is the saturation of the wetting phase, Vt is the total infiltration rate, PcS represents the capillary pressure, and fS=μoKrw/μnKrw+μwKrn is the flux coefficient of the wetting phase. In this example, Vt=1 and porosity and absolute permeability of the layer are nondimensionalized to be ϕ=1 and K=1, respectively. Relative permeability of the wetting phase (w) and the nonwetting phase (n) are Krw=S4 and Krn=1+S1-S3, respectively, and their viscosity is μw=μn=1. The capillary pressure is the J-function Pc=0.1S-0.5.The convective term is typically nonconvex. Another problem is that the diffusive term is nonlinear, which will cause infinite saturation gradients at the front of each phase. Most methods will smear this singularity out because of their numerical diffusion, no matter how small. With the range-discrete mesh method, this singular feature can be maintained by dealing the front as a moving boundary condition.The numerical solutions at different time withΔS=0.05 are shown in Figure 8. It is obvious that there exists not only a steep front where saturation varies quickly but also two singular points at S→0 and S→1, where the saturation gradient tends to be infinite. As a benefit of the range-discrete mesh and the moving boundary condition, the steep front is described preciously and two singular points are maintained intactly with only 20 mesh points.Figure 8
Numerical solutions of the Buckley-Leveret equation att=0.25 and t=0.5. The exact solution shown in solid is obtained by FDM with a fine mesh.


## 7. Conclusion

Achievements of this paper are summarized as follows:(i)
Traditional adaptive mesh methods involve complex adaptive operations and introduce numerical dissipation by interpolation.(ii)
We have developed a self-adaptive numerical method based on uniformly distributed range-discrete mesh for convection-diffusion equation. The governing equation is discretized by integration on the dynamic control volume of each mesh point. The goal of adaptivity and solving the equation are achieved simultaneously by solving the discretized equation, which moves the mesh points but does not change their values.(iii)
The method was verified against the analytical solution of serval numerical examples. Steep fronts were tracked accurately with only a few mesh points. The residual error is found to be independent of the magnitude of diffusion.(iv)
With a moving boundary condition, singular points at the fronts were maintained by the proposed method. It is found that, despite its magnitude, a nonlinear diffusion coefficientϵS=Sαα≥1 can produce singularities where the gradient of the solution can be discontinuous or infinite.In all, the proposed method has shown advantages in solving convection-diffusion equations with high efficiency. The method can serve as a fast numerical tool for a wide range of front propagation problems, such as underground water flow, chemical floods of oil fields, and air pollution problems. On the other hand, it is obvious that, for higher dimensions, the elements in the range-discrete mesh become contour lines or contour planes and the solution can be achieved by tracking their movements. This work is in our next research scope.



---

*Source: 8379609-2017-07-17.xml*