Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 8379609, 13 pages
https://doi.org/10.1155/2017/8379609

Research Article

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

Zhi-Wei Cao, Zhi-Fan Liu, Zhi-Feng Liu, and Xiao-Hong Wang

Department of Thermal Science and Energy Engineering, University of Science and Technology of China, Hefei, Anhui 230026, China

Correspondence should be addressed to Zhi-Fan Liu; zfliu13@ustc.edu.cn

Received 5 April 2017; Accepted 11 June 2017; Published 17 July 2017

Academic Editor: Mohsen Sheikholeslami

Copyright © 2017 Zhi-Wei Cao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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]:

where the diffusive coefficient 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. Section 2 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, , a series of discrete values can be obtained as , where . Then, one can get a series of intersection points of the line sets and the curve , namely, as graphed in Figure 1. Then, with the location of each intersection known, approximation of 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, and shown in Figure 1, a constant piece is sufficiently described.

If is nonmonotone, serval intersections of the line and the curve exist, for example, points shown in Figure 1. It is intuitive to employ all these intersections, assigned with the same discrete value , as different mesh points to describe the nonmonotone feature. Moreover, a mesh point should also be set at the position of an extremum of , for example, point in Figure 1, because, without it, the solution between and 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 can be different from each other, we focus our attention on the equidistance situation for the sake of simplicity.

If there exit mesh points in the spatial area , the density of mesh points can be defined as . For a range-discrete mesh, variation of the solution over the area can be defined as . Thus, we can obtain

where 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, as , we can observe an equidistribution of the solution gradient 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 in Figure 2 as an example, and let and , making the control volume moving with the flow. Integrating (1) on this control volume and a small time interval where , one can get

where the integration of the convective and diffusive term can be done directly as

As for the transient term, according to Appendix A, if the solution curve is monotone in , we have

Consequently, the integration form of (1) is

The item on the left-hand side of (6) is the variation of the total quantity of the control volume over the time interval ; and 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

where

If the items on the right-hand side of (7) are the values at , the discrete equation is explicit; and if they are values at , 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 , not , 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 . To make the range-discrete grid work for nonmonotone solutions, adjustments are needed.

4. Adjustments for the Extremum Point

As has been mentioned in Section 2, a mesh point, for example, 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 volume for an extremum mesh point is shown in Figure 3, where . According to Appendix A, The items inside the brackets in (3) for the control volume of an extremum point are

So, we have

With an parabolic approximation of the solution curve between , we obtain

Consequently, we have

Thus, the discrete equation for an extremum point is obtained.

One can conclude from (10) that, with a positive , left-hand side of (10) always decreases as time goes on. It means that the value of the extremum is approaching , 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 points and will be removed from the range-discrete mesh if is less than a threshold value, for example, . Then, let 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, , can be used directly in the discrete equation in the same way as FDM. That is, one can consider a discrete point , set at the boundary to provide the boundary information, with its position and value always fixed as and during the computing process.

5.2. Moving Boundary Condition

As has been mentioned in Section 2, 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 have and (see Appendix B). Thus, it can be involved in the range-discrete mesh by a boundary point with and . 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 threshold . The position of the boundary point is only needed here and can be interpolated using by two nearby mesh points, where 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 Appendix B, can be maintained by the proposed method. Numerical tests in Section 6 also confirm this.

5.3. Neumann Boundary Condition

Given a Neumann boundary condition , a boundary mesh point will also be set. The discrete value of is determined dynamically both by the value of , which is the point nearest the boundary, and by the sign of , that is, for the left boundary and for the right boundary. The position of can be calculated by an linear interpolation:

or 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)If is out the calculation area, computing continuous on.(2)If 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)If is out of the calculation area, will be removed immediately and the information of is determined by the new nearest boundary mesh point.

6. Numerical Examples

6.1. A Travelling Wave Solution

With and being a constant, we have Burgers’ equation from (1):

It is easy to verify that this equation has the travelling wave solution:

where is the speed of the wave. With , (15) is the solution of (14) with an initial condition . Figure 4 and Table 1 show the solutions and error tables at different time using the proposed method when . 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. is the number of mesh points.

errors

20
40
80
160

error

error

error

error

error

error

error

error

error

error

error

error

(a)

(b)

(c)

Figure 4: Numerical solutions at and of the travelling wave with different diffusive coefficients. Solid line represents exact solutions and squares, circles, and triangles in figures are numerical solutions. The figures below each figure show the spatial distribution of the mesh points at different time.

One 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, . 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: , and , .

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 be , and , respectively. According to Caldwell and Smith [30], the solution is

where

The numerical solutions at different times with diffusive coefficients 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. is the number of mesh points.

errors

20
40
80
160

error

error

error

error

error

error

error

error

error

error

error

error

(a)

(b)

(c)

Figure 5: Numerical solutions obtained by the proposed method at and of Burgers’ equation with different diffusive coefficients. Solid line represents exact solutions and squares, circles, and triangles in figures are numerical solutions. The figures below each figure show the spatial distribution of the mesh points at different time.

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:

with two Dirichlet boundary conditions , and a lamp initial condition, which is and . Diffusive coefficients set as with are tested. Solutions at are shown in Figure 6.

(a)

(b)

Figure 6: Numerical solutions of different nonlinear diffusive coefficients: (a) and (b) .

It is shown by Figure 6 that the gradient of the solution with tends to be zero at the front and thus a smooth solution. However, the gradient with 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 , which is also a conclusion of Appendix B.

The solutions with diffusive coefficients set as and with 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.

(a)

(b)

Figure 7: Numerical solutions of different nonlinear diffusive coefficients: (a) and (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 region is full of the wetting phase, for example, water, and the region is occupied by a nonwetting phase, for example, oil, initially. The governing equation can be obtained by Buckley-Leverett theory [31]:

where is the saturation of the wetting phase, is the total infiltration rate, represents the capillary pressure, and is the flux coefficient of the wetting phase. In this example, and porosity and absolute permeability of the layer are nondimensionalized to be and , respectively. Relative permeability of the wetting phase (w) and the nonwetting phase (n) are and , respectively, and their viscosity is . The capillary pressure is the -function .

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 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 and , 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 at and . 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 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.

Appendix

A. Integration of the Transient Term in (3)

In (3), we have the integration of the transient term over the dynamic control volume . To understand this term, first, we know

So, integration of the transient term in (3) is

The items inside the brackets in (A.2) are

which can be understood according to Figure 9. Consider the following three situations:

(1) is monotone decreasing in , as shown in Figure 9(a). Then (A.3) is the shaded area shown in Figure 9(a). So, we have

(2) is monotone increasing in , as shown in Figure 9(b). Then (A.3) is the shaded negative area shown in Figure 9(b). So, we have

(3) is nonmonotone in , as shown in Figure 9(c). As we have , so

 is the shaded area or negative area according to the sign of .

(a)

(b)

(c)

Figure 9: Geometric description of the integration.

In short, the integration of the transient item is

B. The Moving Boundary Condition

Assume that there is a piece of constant solution at time , a moving boundary condition is formed at either end of the constant piece. In this section, without loss of generality, the constant piece is , which means that , and is the so-called moving boundary.

A dynamic control volume , where and is an arbitrary point in the constant piece, is chosen as shown in Figure 10. For , we have

and also . For this control volume, the integration equation becomes

For , we have ; thus (B.2) becomes

Therefore, (B.3) and both hold for the front point and will serve as the boundary condition for range-discrete method.

Figure 10: The description of a moving boundary condition.

However, is not a given condition but rather a result of the solution. It can be calculated by an interpolation as one of the two following cases:

(1)If when or : from (B.3), we can conclude that . Thus, using the Taylor series, the solution when is

 Dropping the items of order more than , we have an approximation as

(2)If when or : an approximation behavior of the diffusive coefficient can be obtained as when . And we assume that the solution when can be approximated by . can be calculated using the limitation that the quality flux cannot be infinite, that is, when , where

 Thus, should approach zero at the same speed as when . So, we have , and finally

 Thus, , and so we have . This means that, when ,

It can be seen from (B.5) that when , the gradient of the solution tends to be infinite when . The situation when graphed in Figure 7 shows us an example.

With the information of the two neighboring mesh points known, can be obtained approximately using (B.5) or (B.8).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (nos. 11572315 and 11202205) and CNPC-CAS Science and Technology Cooperation Project (no. 2015A-4812).

References

B. A. Finlayson, Numerical Methods for Problems with Moving Fronts, Ravenna Park Publishing, Incorporated, Seattle, WA, USA, 1992.

L. Debnath, Nonlinear Partial Differential Equations for Scientists and Engineers, Birkhäuser Boston, Inc, Boston, MA, USA, 3rd edition, 2011.

R. E. Ewing and H. Wang, “A summary of numerical methods for time-dependent advection-dominated partial differential equations,” Journal of Computational and Applied Mathematics, vol. 128, no. 1-2, pp. 423–445, 2001.

R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, UK, 2002.

A. Kurganov and E. Tadmor, “New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations,” Journal of Computational Physics, vol. 160, no. 1, pp. 241–282, 2000.

T. Nazir, M. Abbas, A. I. Ismail, A. A. Majid, and A. Rashid, “The numerical solution of advection-diffusion problems using new cubic trigonometric B-splines approach,” Applied Mathematical Modelling. Simulation and Computation for Engineering and Environmental Systems, vol. 40, no. 7-8, pp. 4586–4611, 2016.

P. K. Sweby, “High resolution schemes using flux limiters for hyperbolic conservation laws,” SIAM Journal on Numerical Analysis, vol. 21, no. 5, pp. 995–1011, 1984.

K. Alhumaizi, “Comparison of finite difference methods for the numerical simulation of reacting flow,” Computers and Chemical Engineering, vol. 28, no. 9, pp. 1759–1769, 2004.

C.-W. Shu, “High order weighted essentially nonoscillatory schemes for convection dominated problems,” SIAM Review, vol. 51, no. 1, pp. 82–126, 2009.

H. Ding and Y. Zhang, “A new difference scheme with high accuracy and absolute stability for solving convection-diffusion equations,” Journal of Computational and Applied Mathematics, vol. 230, no. 2, pp. 600–606, 2009.

H. K. Dahle, M. S. Espedal, and R. E. Ewing, “Characteristic Petrov-Galerkin subdomain methods for convection-diffusion problems,” Computer Methods in Applied Mechanics and Engineering, vol. 64, pp. 77–86, 1988.

T. F. Russell and M. A. Celia, “An overview of research on Eulerian-Lagrangian localized adjoint methods (ELLAM),” Advances in Water Resources, vol. 25, no. 8-12, pp. 1215–1231, 2002.

Z. Chen, “Characteristic mixed discontinuous finite element methods for advection-dominated diffusion problems,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 23-24, pp. 2509–2538, 2002.

M. Seaïd, “An Eulerian-Lagrangian method for coupled parabolic-hyperbolic equations,” Applied Numerical Mathematics, vol. 59, no. 3-4, pp. 754–768, 2009.

J. Caldwell, P. Wanless, and A. E. Cook, “Solution of Burgers’ equation for large Reynolds number using finite elements with moving nodes,” Applied Mathematical Modelling, vol. 11, no. 3, pp. 211–214, 1987.

J. Behrens, A. Iske, and M. Käser, “Adaptive Meshfree Method of Backward Characteristics for Nonlinear Transport Equations,” in Meshfree Methods for Partial Differential Equations SE - 2, M. Griebel and M. Schweitzer, Eds., vol. 26 of Lecture Notes in Computational Science and Engineering, pp. 21–36, Springer Berlin Heidelberg, Berlin, Germany, 2003.

T. Plewa, T. Linde, and V. G. Weirs, “Adaptive mesh refinement-theory and applications,” in Proceedings of the Chicago Workshop on Adaptive Mesh Refinement Methods, vol. 41, Chicago, Ill, USA, 2003.

X.-H. Wang, M. Quintard, and G. Darche, “Adaptive mesh refinement for one-dimensional three-phase flow with phase change in porous media,” Numerical Heat Transfer, Part B: Fundamentals, vol. 50, no. 3, pp. 231–268, 2006.

W. Huang and R. D. Russell, Adaptive Moving Mesh Methods, Springer Science & Business Media, New York, NY, USA, 2011.

H. Dong, Z. Qiao, S. Sun, and T. Tang, “Adaptive moving grid methods for two-phase flow in porous media,” Journal of Computational and Applied Mathematics, vol. 265, pp. 139–150, 2014.

S. Reboux, B. Schrader, and I. F. Sbalzarini, “A self-organizing Lagrangian particle method for adaptive-resolution advection-diffusion simulations,” Journal of Computational Physics, vol. 231, no. 9, pp. 3623–3646, 2012.

C.-C. Chueh, N. Djilali, and W. Bangerth, “An h-adaptive operator splitting method for two-phase flow in 3D heterogeneous porous media,” SIAM Journal on Scientific Computing, vol. 35, no. 1, pp. B149–B175, 2013.

O. Angelini, K. Brenner, and D. Hilhorst, “A finite volume method on general meshes for a degenerate parabolic convection-reaction-diffusion equation,” Numerische Mathematik, vol. 123, no. 2, pp. 219–257, 2013.

Z. Cao, Z. Liu, X. Wang, A. Shi, H. Luo, and B. Noetinger, “A dissipation-free numerical method to solve one-dimensional hyperbolic flow equations,” International Journal for Numerical Methods in Fluids.

F. J. Fayers and J. W. Sheldon, “The effect of capillary pressure and gravity on two-phase fluid flow in a porous medium, petroleum transactions,” AIME, vol. 216, pp. 147–155, 1959.

F. Hussain, Y. Cinar, and P. Bedrikovetsky, “A semi-analytical model for two phase immiscible flow in porous media honouring capillary pressure,” Transport in Porous Media, vol. 92, no. 1, pp. 187–212, 2012.

I. R. Goumiri, J. H. Prévost, and M. Preisig, “The effect of capillary pressure on the saturation equation of two-phase flow in porous media,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 36, no. 3, pp. 352–361, 2012.

J. Mackenzie, “Uniform convergence analysis of an upwind finite-difference approximation of a convection-diffusion boundary value problem on an adaptive grid,” IMA Journal of Numerical Analysis, vol. 19, no. 2, pp. 233–249, 1999.

N. Kopteva and M. Stynes, “A robust adaptive method for a quasi-linear one-dimensional convection-diffusion problem,” SIAM Journal on Numerical Analysis, vol. 39, no. 4, pp. 1446–1467, 2001.

J. Caldwell and P. Smith, “Solution of Burgers’ equation with a large Reynolds number,” Applied Mathematical Modelling, vol. 6, no. 5, pp. 381–385, 1982.

S. E. Buckley and M. C. Leverett, “Mechanism of fluid displacement in sands,” Transactions of the AIME, vol. 146, no. 1, pp. 107–116, 2013.