Reach Us +44-1474-556909
High Order Numerical Solutions to Convection Diffusion Equations with Different Approaches | OMICS International
ISSN: 2168-9679
Journal of Applied & Computational Mathematics
Make the best use of Scientific Research and information from our 700+ peer reviewed, Open Access Journals that operates with the help of 50,000+ Editorial Board Members and esteemed reviewers and 1000+ Scientific associations in Medical, Clinical, Pharmaceutical, Engineering, Technology and Management Fields.
Meet Inspiring Speakers and Experts at our 3000+ Global Conferenceseries Events with over 600+ Conferences, 1200+ Symposiums and 1200+ Workshops on Medical, Pharma, Engineering, Science, Technology and Business
All submissions of the EM system will be redirected to Online Manuscript Submission System. Authors are requested to submit articles directly to Online Manuscript Submission System of respective journal.

High Order Numerical Solutions to Convection Diffusion Equations with Different Approaches

Don Liu* and Yifan Wang

Haibo Zhang Mathematics and Statistics, Louisiana Tech University Ruston, LA 71272, USA

*Corresponding Author:
Don Liu
Haibo Zhang Mathematics and Statistics
Louisiana Tech University Ruston, LA 71272 USA
Tel: (+1) 318 257 4670
E-mail: [email protected]

Received Date: January 22, 2015; Accepted Date: February 22, 2015; Published Date: March 10, 2015

Citation: Liu D, Wang Y (2015) High Order Numerical Solutions to Convection Diffusion Equations with Different Approaches. J Appl Computat Math 4:208. doi: 10.4172/2168-9679.1000208

Copyright: © 2015 Liu D, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Visit for more related articles at Journal of Applied & Computational Mathematics


Spectral element method; Finite element method; Compact difference method; Convection diffusion equation; Lagrange interpolation


Numerically solving nonlinear PDE has been on the center stage for computational sciences, especially in computational fluid dynamics (CFD), a mature interdisciplinary area with many applications. The key component of CFD is obtaining reliable and accurate numerical solutions to governing equations such as Burgers’ equations, convection diffusion equations, and Navier-Stokes equations under a variety of boundary and initial conditions in different dimensions. Convection diffusion equation and Burgers’ equation have been widely used to test numerical algorithms as they are simple nonlinear PDE describing convection and dif-fusion of a fluid, gas dynamics, shock waves, acoustic waves, and traffic flows etc. Handling nonlinear terms is central to solving convection diffusion equations and Burgers’ equations.

Previous studies on convection diffusion and Burgers’ equations are mainly analytical and numerical. Courant et al. [1] proposed the method of characteristics to solve Burgers’ equation. But it has a limitation and would not work for the inviscid Burgers’ equation because the characteristics may intersect at certain time and may result in discontinuities in solution. Hopf et al. [2] introduced a nonlinear transformation to convert an original PDE into a linear heat equation which can be solved with less effort. Zhu et al. [3] reported a discrete Adomain decomposition method, which derived a fully implicit finite difference scheme for Burgers’ equation, and nonlinear terms was defined by the infinite series of Adomian’s polynomials. Srivastava et al. [4,5] used finite difference method to discretize spatial derivatives and adopted Newton’s root-finding method [6,7] to solve a system of nonlinear equations without linearizing Burgers’ equation. Although this method has the advantage of solving Burgers’ equation implicitly, generating Jacobian matrix of the system is computational expensive.

Nonlinear terms could be numerically solved either with a linearization or without a linearization such as using conservative forms. Linearization could relieve the constraint on the size of time steps due to stability condition because a linearized PDE could be solve implicitly with much larger time step than the time step of an explicit method. For very large systems of nonlinear PDE, linearization is especially useful in decoupling systems of PDE and making them possible to be solved individually and sequentially. For nonlinear PDE in multiple dimensions, linearization also facilitates decoupling and makes it possible to solve one component ahead of time. On the other hand, a relatively recent and successful approach is to cast a nonlinear term into a conservative form and use numerical fluxes to solve conservation laws [8-13]. This approach is very effective for hyperbolic equation where convection is dominant and information propagates along certain direction.

There are many methods to solve convection diffusion equations. Compact difference method (CDM) [14-18], due to its simple nature and convenience in implementa-tion, offers high order algebraic convergence and uses smaller stencils than finite difference method. But the linear system is usually twice or even three times bigger than the one from finite difference. Finite element method (FEM) is known for handling geometrical complex- ity. Fletcher et al. obtained finite element solution to convection diffusion equations [19,20]. However, FEM is limited in acquiring higher order accuracy as basis functions become less orthogonal mutually once the order of basis functions increases beyond fifth and mass ma- trices become ill-conditioned. For high accuracy resolution, although one could use h-type discretization refinement, i.e., smaller and more elements, rounding errors would accumulate and even defeat accuracy at certain point. Spectral element method (SEM), first appeared in [21], could be a good alternative. It achieves very high accuracy (e.g. 15th) by using or- thogonal polynomial bases and zeros of orthogonal polynomials as quadrature points. SEM is capable of hp-refinement, and especially the p-type refinement which enhances resolution without extra numbers of elements [22-24]. Due to the capability of handling complexity geometry and both hp-type refinement for high accuracy, SEM has successfully appeared in computational fluid dynamics [21,25-34], especially in simulating microfluidic devices [35-41]. SEM could also be used to model advanced microfluidic biosensors [42] with applications in magnetic labeling, sorting, medical diagnostics, and weak magnetic field detection [43].

Another relatively new approach, discontinuous Galerkin finite element method (DG- FEM), has been successful in handling nonlinear PDE in which nonlinear terms are written in conservative forms [22,44-54]. Using discontinuous basis functions and appropriate numerical fluxes at boundaries, DG-FEM is free from the C0 constraint on basis functions and is capable of capturing discontinuity in the solution. How to choose numerical fluxes is crucial to some strong hyperbolic equations.

In this paper, we discuss four different approaches, namely Method I: Taylor expansion, Method II: Cole-Hopf transformation, Method III: Compact Difference Preprocessing, and Method IV: Conservative Form, to handle nonlinearity. With these approaches, we produce numerical solutions to one and two dimensional convection diffusion equations using a blend of CDM, FEM, and SEM. Convergence rates are demonstrated in h-type and/or p-type refinement tests. Pros and cons of different approaches are discussed to provide a reference for solving nonlinear PDE in multiple dimensions and systems of nonlinear PDE.

Methodology and Numerical Examples

We consider convection diffusion equations in one or two dimensional domains:

Equation      (1)

where μ is the dynamic viscosity. In some occasions, for the convenience of discussion, we treat the viscosity as a small value and drop the diffusion term; therefore, the convection diffusion becomes the inviscid Burgers’ equation. In order to verify the accuracy of numerical solutions, we need an exact solution to compute the norm of errors at all quadrature points. Therefore, the convection diffusion equation was modified with the diffusion term replaced by a known function of space and time. In some examples, we solve the following modified convection diffusion equations instead:

Equation      (2)

In the following subsections, four different approaches are discussed to handle nonlinearity in one dimensional (1D) or two dimensional (2D) situations. Numerical solutions are provided with SEM, FEM, CDM or a blend of the above.

Linearization with Taylor expansion (Method I)

We start with the following simplest nonlinear PDE:

Equation      (3)

Equation      (3a)

Equation      (3b)

and use two different discretization methods in time and space in subsequent discussions.

Backward Euler and finite element method

Discretizing Equatation. (3) implicitly in time, we obtain:

Equation      (4)

Due to the nonlinear term, we could not directly treat the entire convective term implicitly, because after a Galerkin projection, there will be three different indices for the nonlinear term. Therefore, we attempted to linearize Equatation. (4). The first approach to handle nonlinearity, denoted as Method I, applies the Taylor expansion to the coefficient of un+1 in time and repetitively replace time derivatives with space derivatives using the exact information in the governing PDE. This linearization is similar to the general idea of Lax-Wendroff scheme [55-57], except we use first or higher order Taylor expansion [58] to achieve higher accu- racy and implement a finite element or high order difference method to compute the spatial derivatives. The CFL number here is no more than 1. The first order approximation for the nonlinear term has the second order truncation error in time:

Equation      (5)

Substitute Equatation. (5) into Equatation. (4), we have:

Equation      (6)

Using the exact information from the PDE: un = f n − unun, we substitute it into Equatation. (6):

Equation      (7)

Let (un)*=un+Δtfn−Δtunun, which is known at the time step n, then Equatation. (7) becomes a linearized ordinary differential equation:

Equation      (8)

Next we apply the nodal Galerkin projection to form a linear system A Uˆ=B . We divide the domain into N elements. Within each element Ωe, we expand u and ux in kth order Lagrangian basis functions over k+1 uniform points in order to simplify the evaluation of (un)*, although an alternative way could use Lagrangian basis function over Gauss-Lobatto- Legendre points. In a typical element Ωe:

Equation      (9)

Equation      (10)

Substituting Equs. (9) and (10) into Equatation. (8) and applying Galerkin projection, we obtain the weak form of Equatation. (8) within one element Ωe:

Equation      (11)

where (Un)* denotes the expansion term for (un)* at the point of i. After global assembling, the matrix-vector form of Equatation. (8) is:

Equation      (12)

Equation      (13)

where, M is the global mass matrix; Fn+1 is obtained by global assembly; the matrix K is obtained by the scalar multiplies global advection matrix. If we choose a linear expansion basis, the matrix K has a tridiagonal form as below, while if quadratic and cubic basis function and pentadiagonal and seven-diagonal matrix for using respectively:

In which, Ul=Uˆn+ΔtF n−Δt(Uˆ nUˆ n)l, l stands for the numbering of global points, and each Klm corresponds to the entry of global advection matrix at l, m. Here, l and m is from 2 to N, since boundary points at 1 and N+1 are known. The key idea here is to treat (un)* explicitly and take it as a scalar multiplier for un+1, which is treated implicitly, then include it into the global advection matrix to form matrix K. Therefore, we obtain a set of discretized equations with only two indices which are represent Table 1 in a linear system.

Order α β A B C
4th 1/4 0 3/2 0 0
6th 1/3 0 14/2 1/9 0
8th 4/9 1/36 40/27 25/54 0

Table 1: Where the coefficients are.

Finally, Equatation. (12) becomes:

Equation      (14)

To examine the accuracy of this approach, we set up the problem as below:

Equation      (15)

subject to the boundary and initial conditions:

Equation      (16)

The analytic solution to this well-posed problem is:

Equation      (17)

Since this problem is time-dependent, time integrations are crucial in numerical solutions. We attempted both implicit and explicit time integration. The CFL condition (for the latter) and the diffusion condition were satisfied by restricting the time step to be small enough:

Equation      (18)

Equation      (19)

where Umax is the maximum of absolute phase velocity, Δt and Δx are time step and distance between quadrature points, respectively, and μ is the diffusion coefficient Table 2.

Order α β γ A B C D E G H
4th 1 3 0 -17/6 3/2 3/2 -1/6 0 0 0
5th 1 4 0 -37/12 2/3 3 -2/3 1/12 0 0
6th 1 5 0 -197/60 -5/12 5 -5/3 5/12 -1/20 0
7th 1/10 1 1 -227/600 -13/12 7/6 1/3 -1/24 1/300 0
8th 1/12 1 5/4 -79/240 -77/60 55/48 5/9 -5/48 1/60 -1/720

Table 2: Where the coefficients are.

When we choose linear, quadratic, and cubic bases in space and backward Euler in time, the error estimations are order of O(Δt1+Δx2), O(Δt1+Δx3), O(Δt1+Δx4), respectively. We compute the Euclidean norm of errors (L2 Error) at all nodes. Figure 1 presents the L2 Error versus the number of elements in the log-log scale at a fixed time step Δt=10−4. The convergence rates are basically one order higher than the order of expansion functions. To investigate the order of accuracy in time, we fix element number N=10, and vary Δt=0.01, 0.005, 0.0025, 0.00125. Using a first order Taylor expansion in the log-log scale, the first order convergence in time is observed in Figure 2. When we switch it to the second order Taylor expansion, a second order accuracy in time is obtained as shown in Figure 3.


Figure 1: Convective Equ. (15) linearized with Method I: Algebraic convergence in space at a fixed time step Δt=10−4.


Figure 2: Convective Equ. (15) linearized with Method I: Temporal convergence at a fixedspatial resolution of using 10 elements for a total time of 0.2 dimensionless unit.


Figure 3: The convective equation (15) linearized with Method I: Convergence in space at afixed time step.

Crank-Nicholson and compact difference method

We test the same approach for nonlinearity but using the Crank- Nicholson scheme in time and compact method in space for Equatation. (3):

Equation      (20)

To linearize the term un+1 uxn+1, we use the Taylor expansion in time again to express un+1 with some terms in the time level n, such that un+1 is replaced with u*. Therefore, Equatation. (20) becomes:

Equation      (21)

In order to keep the 2nd order temporal accuracy, we use a few more terms in the expansion such that

Equation      (22)

After truncating the 2nd order term, we have:

Equation      (23)

in which, utn can be obtained from the PDE exactly:

Equation      (24)

Hence, Equatation. (20) becomes linearized as below:

Equation      (25)

Rearrange the above equation, we have:

Equation   (26)

We use the following 4th order Pa´de scheme to compute un+1 and uxn+1 simultaneously:

Equation      (27)

To examine the spatial and temporal accuracy of using Crank- Nicholson and Compact Difference, we use the same problem as before, Equatations. (15) and (16), with the same exact solution Equ. (17). For the spatial error only, we fix Δt=0.0001, compute u(x, 0.0001), and compare it with the exact solution in L2 error. We notice that the number of elements N=10, 20, 30, 40 increases, the error goes down consistently. We observed that the order of convergence rate is close to 4th in space, as shown in Figure 3.

Figure 4 illustrates a temporal convergence rate close to the second order at the total time of 1. In this Figure 5, we set the time step to be Δt=0.001, 0.0005, 0.00025, 0.000125, use only 21 elements in space. Since the Crank-Nicholson was used in time and 4th order Pa´de schemes in space, the second order convergence in time is anticipated.


Figure 4: The convective equation (15) linearized with Method I: Convergence in time at afixed spatial resolution.


Figure 5: Evolution of numerical solution of 1D convection diffusion equation (28) linearizedwith Method II: Cole-Hopf transformation.

In summary, the idea of using Taylor expansion to express an unknown of time level n+1 explicitly to approximate nonlinearity is convenient and stable, provided that CFL condition is obeyed. This idea is rather simple and could be used as long as an explicit expression for time derivative is possible or other approaches are difficult to use.

Linearization with Cole-Hopf transformation (Method II)

One-dimensional convection diffusion equation: Now we include the diffusion and consider the one dimensional convection diffusion equation:

Equation      (28)

with the initial condition and two Dirichlet boundary conditions below:

u(x, 0) = f (x),0 ≤ x ≤ 1,      (28a)

u(0, t) = α(t),u(1, t) = β(t),t ≥ 0.      (28b)

With the Cole-Hopf transformation [2]:

Equation     (29)

the original equation becomes a linear heat equation of in w(x, t),

Equation     (30)

with new initial and boundary conditions as below:

Equation     (30a)

Equation     (30b)

Equation     (30c)

We solve the transformed linear equation, then use the inverse transformation, i.e., Equatation. (29), to acquire the solution to the original problem.

Two-dimensional convection diffusion equation: For twodimensional domains, the Cole-Hopf transformation [59] is given below:

Equation     (31)

Equation     (32)

We could transform the vector form of convection diffusion equation, Equatation. (1), into a scalar equation,

Equation     (33)

To obtain numerical solutions, we compute w in Equatation. (33) with the nodal spectral element method using high order Lagrangian interpolants. Then we use high order finite difference schemes to compute wx and wy from w. Finally, u and v could be obtained from w using Equatations. (31) and (32), respectively. When we use the kth order Lagrangian interpolants as the basis functions, in order to guarantee the exponential convergence, we used k + 1st order difference schemes to compute wx and wy .

Numerical results: One-dimensional Example:

To validate Method II in one dimension (1D), we solve Equatation. (28) with μ=0.02, u(x,0)=sin(πx), 0 ≤ x ≤ 1, and α(t)=β(t)=0. After the Cole- Hopf transformation, we have the heat equation in the new variable w(x) with new initial and boundary conditions:

Equation     (33)

Equation     (34)

Since the original nonlinear PDE becomes linear, we use implicit treatment in time and finite element method in space. After we obtain the solution in w(x), we use Cole-Hopf inverse transformation of to obtain the numerical solution to the original problem. Figure 5 shows the numerical solution to the 1D problem at different time. The number of elements is fixed to be 10. As in the figure, the wave propagates towards right and forms a wave front, and the magnitude of the wave decreases with the time due to the viscous dissipation. We used cubic basis functions, as the number of elements varying from 10, 20, 30, 40, to 50, the 4th order convergence rate in space was obtained as shown in Figure 6, where the time step was Δt=0.0001. For the following two dimensional (2D) situations, we test two examples using different initial and boundary conditions.


Figure 6: Spatial convergence for the convection diffusion equation (28) linearized withMethod II: Cole-Hopf transformation.

Two-dimensional Example I:

We set the viscosity μ=0.02 and the initial and boundary conditions for Equatation. (1) as below:

u(x, y, 0)=sin(πx) sin(πy), v(x,y,0)=y,Ω : 0 ≤ x, y ≤ 1, t ≥ 0.      (36)

u(x,y,t)=0,v(x, y, t) = y,x, y ∈ ∂Ω,      (37)

We perform the Cole-Hopf transformation to reduce a vector PDE in u(x,y,t) into a scalar one in w(x,y,t). The domain was divided into 4 elements. Nodal SEM with 5th order basis function and Lagrangian interpolants over Gauss-Lobatto-Legendre points were used to solve the 2D scalar heat equation. We still use Lagrangian interpolation and difference method to compute wx and wy. By Equatations. (31) and (32), we found the solutions to u(x,y) for the original problem. Figure 7a and b show the contour lines for the velocity component u and v at t=0.1, respectively. Figure 7c illustrates the elements and quadrature points. Figure 7d shows the vector field of u=(u,v). Because velocity components u and v were asymmetric and coupled together, the vector field convects in both x and y directions.


Figure 7: The two-dimensional convection diffusion equation (28) (Example I) linearized withMethod II: Velocity components u, v, the computational elements and quadrature points.

Two-dimensional Example II:

In the second test, we change the domain into Ω={(x, y): 0 ≤ x, y ≤ 2π}, set μ=0.005, and use new initial conditions

u(x,y,0)=e(1−(x−π)2−(y−π)2),      (38)

v(x,y,0)=sin(x)sin(y),      (39)

and new boundary conditions on the boundary ∂Ω

u(x,y,t)=0,v(x,y,0)=0,x,y ∈ ∂Ω.      (40)

The domain was divided into 16 elements. We used the nodal SEM with 5th order basis function in each direction in each element to solve the 2D heat equation. Then, we used Lagrange interpolation and compact schemes to compute the first derivatives. Finally, solu- tions for u and v were obtained from Equatations. (31) and (32). Numerical results were shown at time t=0.5. Figure 8a and 8b show contour lines for u and v, respectively. Figure 8c illustrates 16 elements and quadrature points. Figure 8d is the vector field for u=(u, v).


Figure 8: The two-dimensional convection diffusion equation (28) (Example II) linearizedwith Method II: Velocity components u, v, the mesh with 16 elements and quadrature points.

Generally speaking, Cole-Hopf transformation offers two advantages: reducing a vector equation into a linear scalar one and allowing an implicit temporal treatment in which the time step is unrestricted by the stability condition and time integration is unconditionally stable. However, Cole-Hopf transformation could only be applied to certain nonlinear PDE [58], such as hyperbolic PDE which contains second or higher spatial derivatives [60] and the Korteweg-de Vries equation (KdV equation) [61]. Chu et al. [62] listed a system of PDE with second order spatial derivatives that Cole-Hopf transformation is applicable.

Compact Difference Preprocessing (Method III)

Consider the 1D nonlinear problem in Equatations. (3), (3a), (3b), and (3c), we linearize it by treating ux explicitly as the coefficient of u. This is called Method III in this paper, which is an explicit and nonconservative method [58]. The linearized equation below could be solved with an explicit scheme:

Equation      (41)

or implicit scheme:

Equation      (42)

The success of Method III relies on the accuracy of evaluating ux. We compute the nodal values of ux with high order compact difference [15-17,63] schemes before multiplying u. To demonstrate the idea of Method III, we compute ux at interior points with sixth order compact difference scheme as shown below:

Equation      (43)

where h=1 . For boundary points, we use Dirichlet boundary conditions (for simplicity) and the following sixth order scheme:

Equation      (44)

Once we computed the values of (ux)in at all points, Equatations. (41) and (42) are linearized.We divide the whole domain into a total of N elements (Ωe) and express un and un+1 with expansions of basis functions (for example, linear) on each element:

Equation      (45)

Equation      (46)

Then we use the Galerkin projection to acquire the weak form of Equatation. (42):

Equation      (47)

After global assembling procedure, the linear system of Equatation. (42) in vector form is:

Equation      (48)

in which, Equationis the global mass matrix; Equation is the assembled RHS of the Equatation. (47) over all element. The matrix Equation consists of the products of the global mass matrix and the explicit values of Equation given by Equatations. (43) and (44). The matrix M* has the following form:

For the convenience of time integration, Equatation. (48) can be written as:

Equation      (49)

Method III can readily be a high order scheme by adopting higher order compact schemes to evaluate ux on equispaced grids and use higher order polynomials as basis functions in a Galerkin finite element method. Using cubic or quartic basis functions could provide solu- tions with decent algebraic convergence. Beyond fifth order, the convergence rate starts to deteriorate due to the Runge’s phenomenon and a spectral element method using orthogo- nal bases is a better alternative. We demonstrate this by using Lagrangian interpolation for ux from a uniform finite difference grid to Gauss-Lobatto-Legendre quadrature points. For two dimensional situation, the same idea works since Lagrangian interpolation in 2D can be utilized to calculate first derivatives of u at Gauss-Lobatto-Legendre points.

To test the accuracy of Method III in 1D, we consider the previous problem with the initial and boundary conditions in Equatations. (15) and (16). We choose linear basis functions in space and forward Euler method in time, and use the 2nd order central difference scheme to compute value of ux. Therefore, the estimated error is of order O(Δt1 + Δx2). The CFL condition is obeyed with Δt=10−4. By varying the spatial resolution, the number of elements from N=5 to 35, we present the spatial convergence in the log-log scale in Figure 9. The order of convergence is almost 2.


Figure 9: The 1D convective equation linearized with Method III: Convergence in space ata fixed time step.

To examine the order of accuracy in time, we fix N=10 and decrease Δt=0.01, 0.005, 0.0025, 0.00125. The time convergence rate is almost one, as shown in the log-log scale in Figure 10. This convergence rate is expected because the forward Euler method was used.


Figure 10: The 1D convective equation linearized with Method III: Convergence in time ata fixed spatial resolution.

For two-dimensional convection diffusion equations, we denote velocity components with u, v along x and y and let f, g be corresponding forcing terms. We perform linearization and set up the following implicit schemes:

Equation      (50)

Equation      (51)

The derivative terms, ux, uy , vx and vy , at interior points, were computed from u and v at time level n using compact difference method given below [64] and were treated as known coefficients in Equatations. (50) and (51):

Equation      (52)

Equation      (53)

After a Galerkin projection, Equatations. (50) and (51) could form a 2N by 2N system:

Equation      (54)

where M is global mass matrix. M1, M2, M3 and M4 are constructed by following the same pattern as the products of global mass matrix and corresponding spatial differentiation (ux)n, (uy )n, (vx)n and (vy)n respectively.

To facilitate validation, we choose exact solutions of u, v as following:

u = sin(x-t)sin( y-t),      (55)

Equation      (56)

Equation      (57)

Equation      (58)

and the initial and boundary conditions are acquired from Equatations (55) and (56). After we computed the first derivatives of velocity components with respect to x and y, we use two dimensional Lagrangian interpolants to compute their values at Gauss-Lobatto-Legendre quadrature points. Figure 11a and 11b show the contour lines of u and v at time t=0.5. Figure 11c shows four elements and the interior quadrature points. Figure 11d shows the velocity vector field. For all elements, the highest order of polynomial is 8. We pre-compute ux, uy, vx and vy with compact difference schemes of 8th order accuracy. Then we interpolate values of u and v on uniform grids to Gauss-Lobatto- Legendre points and use the nodal SEM with 8th order basis expansions, which are Lagrangian interpolants on Gauss-Lobatto-Legendre points. The values of u and v could be computed with high accuracy.


Figure 11: Two dimensional convection diffusion equation linearized with Method III: Velocity components u, v and computational mesh.

Figure 12 presents exponential convergence of L2 norm of errors of u as the polynomial order varying from 2 to 8. We only derived and tested 8th compact difference schemes for pre-processing first derivatives. It has been observed that Method III works well for convection diffusion equations. For irregular domains, we have to use a mapping and precompute the first derivatives on a mapped grid or using full inclusion of metrics [17] which could add difficulty to Method III.


Figure 12: Two dimensional convection diffusion equations linearized with Method III: p-Refinement.

In summary, the approach of Compact Difference Preprocessing uses previously obtained values in prior time steps to compute unknowns explicitly. This way of handling nonlinearity provides convenience in decoupling PDE and especially for systems of PDE or multiple dimensions. It opens up the possibility of using an implicit method in time in which a large time step is feasible and is especially good for long time integration.

Using conservative form (Method IV): The one dimensional convection diffusion equation, Equatation. (28), could be written in conservative form [65], since it is always conservative in one dimensional situation (See the Appendix in Section 4). The approaching of casting it into conservative form is called Method IV.

Equation      (59)

We use the forward Euler scheme in time for the above equation:

Equation      (60)

and discretize the whole domain into N elements. Within a typical element Ωe, we expand terms such as un, (u2)n and un with the kth order Lagrangian interpolants φ on the k+1 2 xxx Gauss-Lobatto- Legendre points:

Equation      (61)

Equation      (62)

Equation      (63)

We substitute the above expansions into Equatation. (60) and apply a Galerkin projection to obtain the weak form of Equatation. (59) within Ωe:

Equation      (64)

Note that within each time step, (u2)n is obtained explicitly from un. Assembling Equatation. (64) over all elements, we obtain the global form of Equatation. (60):

Equation      (65)

In which, M , K and L are global mass, convection and stiffness matrices, respectively.

Two dimensional convection diffusion equation: For two dimensional convection diffusion equations, Equatation. (1), subject to some initial and boundary conditions, we assume u and v are velocity components along x and y directions, respectively, and write Equatation. (1) in the component form:

Equation      (66)

Equation      (67)

For Equatations. (66) and (67), Method IV is only conditionally conservative [58]. The proof of this property is given in the Appendix in Section 4. Now we change Equatations. (66) and (67) into the partial conservation form:

Equation      (68)

Equation      (69)

Unless we know that u=(u, v) is a conservative vector field, i.e., uy=vx, or incompressible flow field, i.e., ux+vy=0, we cannot take it for granted to change coupled terms vuy and uvx into conservative forms. In Equatation. (68), the coupled term vuy remains nonlinear. We could linearize it by treating v explicitly as known coefficient for uy in Method IV. Similarly, the coupled term uvx in Equatation. (69) is linearized by treating u explicitly as the known coefficient for vx.

In terms of spatial discretization, we first use rectangular elements (irregular elements are discussed later) and choose the kth order Lagrangian interpolants, Φ, on Gauss-Lobatto-Legendre points as the basis functions in both x and y directions, although different orders of bases could be used in x and y. In a typical element Ωe, we expand u(x, y) as below:

Equation      (70)

Similarly for v. The derivatives of u with respect to x and y, respectively, are:

Equation      (71)

Equation      (72)

We perform a Galerkin projection with the test function w:

Equation      (73)

to acquire a linear algebraic system about global variables.

Numerical examples one dimensional example: Consider the same 1D problem, Equatations. (15) and (16) as in Method I, we use Method IV to solve it. We discretize the whole domain into 2 elements and use the forward Euler method in time. To achieve exponential convergence in space, we use SEM and perform a p-refinement by varying the order of polynomial basis functions. Figure 13 shows the exponential convergence with the order of polynomial basis functions versus the L2 norm of point-wise errors. Notice that, when polynomial order is less than 4, the exponential convergence of L2 norm of errors is not reached, because the approximation space have not been well set up yet. When the polynomial order is greater than 16, the L2 norm of error will not decrease any more and oscillate around 10−13 due to double decision rounding error.


Figure 13: One dimensional convection diffusion equation (28) with Method IV: p-Refinement.

Two Dimensional Rectangular Domain

The two dimensional computational domain is set to be Ω={(x, y) : 0 ≤ x, y ≤ 1, t ≥ 0}. Equatation. (1) is the governing equation and the initial and boundary conditions are generated from the exact solution of u and v which are given below:

Equation      (74)

Equation      (75)

We use Method IV to solve the 2D problem on a rectangular domain, which was divided into 4 and 16 elements, respectively. Numerical results on 4 elements are shown in Figure 14 and on 16 elements are shown in 15. In both figures, the top plots ((a) and (b)) show contour lines of the velocity component u and v, separately at t = 0.5. The plot (c) illustrates the elements and the quadrature points. The plot (d) shows the velocity vector field. The Reynolds number is about 200; therefore, the magnitude of solution waves decrease with time due to viscous dissipation. Figure 15 gives more details about the convective flow since it has higher resolution than Figure 14 using 4 elements.


Figure 14: Two dimensional convection diffusion equation, Equ. (1) with Method IV: Velocity components u, v with 4 elements N=4 and the polynomial order Po=10.


Figure 15: Two dimensional convection diffusion equation, Equ. (1) with Method IV: Velocity components u, v with 16 elements N=16 and the polynomial order Po=10.

The exponential convergence in space is shown in Figure 16. Polynomial order were up to 40. The time step Δt was chosen to satisfy the CFL and diffusion conditions. Specifically, we chose Δt=0.001 and compared numerical result of u at time t=0.001 with the exact solution in the infinity norm. Using 16 elements gives faster convergence rate than using 4 elements and the error reaches machine zero at the polynomial order of 24.


Figure 16: Two dimensional convection diffusion equation, Equ. (1) with Method IV: hp-refinement in rectangular domain

Two-dimensional Irregular Domain

For the same problem as given in Equs. (66) and (67) but defined on an irregular domain with similar initial and boundary conditions that we designed here, we map the physical coordinates of (x, y) of a quadrilateral element Ωe to the coordinates (ξ1, ξ2) of a standard element (square) Ωst with the Jacobian:

Equation      (76)

For any convex quadrilateral element, we denote its vertices as A, B, C and D, and their coordinates as Ax, Ay , Bx, By , Cx, Cy , Dx and Dy. The components of the Jacobian are:

Equation      (77)

Equation      (78)

Equation      (79)

Equation      (80)

With a known mapping and Jacobian, we could construct the mass, convection, and stiff- ness matrices for any convex quadrilateral element. Similar procedures were performed and Method IV was used to obtain numerical solutions.

Figure 17 shows the velocity component u, v at time t=0.5, quadrature points, and the velocity vector field. In this figure, the domain was divided into 4 elements and the spectral nodal element method with polynomials of order Po=10 in each direction and each element was used as in the simulation. We completed a very high order solution with nodal basis polynomial functions of order Po=35 using 36 nodes, i.e., Gauss-Legendre-Lobatto points in each direction. This high accuracy solution was used as our “exact” solution to compute nodal errors for lower order runs. Figure 18 shows the L∞ norm of errors at all nodes versus the order of expansion polynomial order. A spectral convergence rate was shown in this figure. The solutions are consistent as the order increases.


Figure 17: Two dimensional convection diffusion equation, Equ. (1) with Method IV: Velocity components u, v and quadrature points in irregular domain


Figure 18: Two dimensional convection diffusion equation, Equ. (1) with Method IV: p-refinement in irregular domain

In summary, Method IV is capable of capturing discontinuity developed in the solution and preserves good convergence rates which is unaffected by irregularity in geometry. This approach is very effective for strong hyperbolic equations provided that the conservative form exists.


In this paper, four different approaches for solving convection diffusion equations in 1D and 2D are discussed. Method I, II and III use linearization procedures which give the advantage of using an implicit time discretization and a large time step without the stability constraint. For the linearized problem, an implicit method in time such as Crank-Nicholson or backward Euler is preferred especially for long time integration. Spatial discretization methods used in this paper include spectral modal and nodal element, finite element, and compact difference. Of course other methods such as discontinuous Galerkin finite element or finite volume etc. are possible.

Since Method I uses a Taylor expansion in time to linearize a convective term, the trun- cation rror in time is consistent with the temporal discretization method. This method is simple and effective for strongly nonlinear terms such as (un+1)2un+1 as long as the first temporal derivative could be expressed in terms of the rest terms in the original PDE. The drawback of Method I is that it requires an explicit expression for the time derivative which could be difficult to derive for PDE with second or higher order temporal derivative and may introduce complexity in linearized equation.

Method II is an excellent approach for convection diffusion equations and possibly for other nonlinear equations if a Cole-Hopf transformation is applicable. The success of this method relies on performing a nonlinear transformation, which could reduce the original nonlinear vector form of PDE into a linear one. This is the major advantage of Method.

II. Apparently, both explicit and implicit time schemes can be used afterwards, although implicit is better in time integration. Nevertheless, the Cole-Hopf transformation is difficult for higher order spatial derivatives. Depending on the specific equation, it could be difficult to find analytical forms of transformed initial and boundary conditions. In general, it is not a trivial work to find a suitable Cole- Hopf transformation for a nonlinear equation.

Method III is a good idea in this paper which uses compact difference method to pre-compute first order spatial derivatives so that a convective term could be approximated by the pre-computed value multiplying the unknown which could be treated implicitly. Basically, Method III uses previously obtained values to compute a part of the unknown explicitly. Handling nonlinearity in this approach provides convenience in decoupling PDE and especially for systems of PDE or multiple dimensions. It gives the convenience to use an implicit method in time in which a large time step is feasible and is especially good for long time integration. However, for irregular domains, Method III becomes complicated due to the spatial discretization and a mapping may be necessary. If SEM is used for numerical solutions, Method III requires deriving difference schemes on uniform grids to compute the first derivatives and then interpolate them to Gauss-Lobatto- Legendre nodes. Method III could be improved by using schemes of full inclusion of metrics [16,17] on non-uniform grids or an efficient mapping between grids to compute the first derivatives.

Compared with the above approaches, in terms of convenience and efficiency, Method IV could be superior in dealing with strong nonlinearity and complex geometry. It is straight forward for general weighted residual methods and is simple for multiple dimensions, since generating conservative form does not rely on mesh grids. However, the disadvantage is that using conservative forms to treat nonlinear terms could make the state vector and the matrix system much larger than without using it. This could rely on more memory in computation. In general, for nonlinear PDE, if possible, one could use Method II and design a Cole-Hopf Transformation, to reduce a multi-dimensional nonlinear PDE into a scalar linear PDE so as to bypass nonlinearity. Depending on the dimension and problem, alternatively, one could select an appropriate linearization with controlled errors to open an option for implicit time treatment, which relax the constraint on time step size and could be beneficial to long time integration. For solutions develops discontinuity or there is discontinuous initial or boundary conditions, Method IV, using conservative form is very effective and capable of capturing discontinuity. The accuracy is usually unaffected by complexity of geometry or discontinuity in the domain.


For Method IV, the one dimensional convection diffusion equation, Equatation. (28), is always conservative. To prove this claim, we use the chain rule backward for the convective term:

Equation      (81)

then we introduce new variables F and G, respectively:

Equation      (82)

Equation      (83)

then we change Equatation. (81) into the following:

Equation      (84)

which is conservative. Therefore, the one dimensional convection diffusion equation is always conservative.

However, the two dimensional convection diffusion equations are not always conservative [58] because of the existence of coupled terms. Again, we use the chain rule backward in the two dimensional component-wise convection diffusion equation for the convective terms to obtain the partial conservative forms:

Equation      (85)

Equation      (86)

Focusing on two coupled terms above, we use the chain rule to form below equations:

(vu)y=vy u + vuy,      (87)

(uv)x=uxv + uvx.      (88)

If the irrotaional condition is valid for the vector field which is defined as u=(u, v), i.e., conservation of the total pressure in an irrotational flow [66], then we have the following:

Equation      (89)

which is equivalent to the vorticity-free condition in fluid mechanics:

∇×u=0,      (90)

By Equatation. (89), we replace vuy with vvx in Equatation. (85) and put vvx in conservative form:

Equation      (91)

By Equatation. (89), we replace uvx with uuy in Equatation. (86) to and put uuy in conservative form:

Equation      (92)

Now Equatations. (91) and (92) are both conservative.

Alternatively, if the incompressible flow condition is satisfied by the vector field which is defined as u = (u, v), i.e., conservation of mass, then we have the following equations:

Equation      (93)

which is equivalent to the divergence-free condition in fluid mechanics:

∇·u=0,      (94)

then we replace −vy with ux in the variation of the convective term vuy in Equatation. (85) so that vuy becomes:

Equation      (95)

then Equatation. (85) becomes:

Equation      (96)

Now the new form of Equatation. (85) is already in conservative form. Similarly, we replace ux with −vy in the variation of the convective term uvx in Equatation. (86) so that uvx becomes:

Equation      (97)

then Equatation. (86) becomes:

Equation      (98)

Now the new form of Equatation. (86) is also in conservative form. To summarize these proofs:

1. The one dimensional convection diffusion equation is conservative.

2. If ∇×u=0 is valid, then the dimensional convection diffusion equation is in conservative form:

Equation      (99)

Equation      (100)

If ∇ · u = 0 is valid, then the two dimensional Burger’s equation is in conservative form:

Equation      (101)

Equation      (102)


This study was supported by National Science Foundation under grants DMS-1115546, DMS-1115527, and Louisiana Board of Regents grant LEQSF (2007-10)-RDA-22. We thank the support from the Louisiana Optical Network Initiative (LONI) and Center for Computation and Technology, Louisiana State University. We are especially grateful to the encouragement from Professor Chi-Wang Shu (Brown University) and his continued constructive comments and insightful advice to improve the quality of this research.


Select your language of interest to view the total content in your interested language
Post your comment

Share This Article

Article Usage

  • Total views: 13768
  • [From(publication date):
    April-2015 - Sep 22, 2019]
  • Breakdown by view type
  • HTML page views : 9820
  • PDF downloads : 3948