Asian Journal of Scientific Research

Year: 2014 | Volume: 7 | Issue: 3 | Page No.: 366-375
DOI: 10.3923/ajsr.2014.366.375

Studying the Effect of Preconditioner Scheme on the Solution of Convection Diffusion Equations

M. Molavi-Arabshahi

Abstract: We discuss a new method for solving unsteady convection–diffusion problems arising high order compact difference approximations. The method can investigate not only the effects of numerical error but also those of uncertainty in a physical model at the same time. In this study, we developed numerical methods by replacing the time and space derivatives by compact finite-difference approximations. Furthermore, numerical experiments are conducted to verify its high accuracy and to compare it in combination preconditioned methods for stability, convergence.

Keywords: preconditioner scheme, convection diffusion equation, compact high order scheme, Finite difference scheme and krylov subspace methods

INTRODUCTION

An attractive problem is the Diffusion-Convection equation for which various numerical methods have been suggested by a number of researchers to solve them (Golbabai and Arabshahi, 2010). Application of this equation may be seen in computational fluid dynamics, computational acoustics and computational hydraulics and electromagnetic for modeling Diffusion-Convection of quantities such as mass, energy, vorticity, heat, etc. (Dehghan and Molavi-Arabshahi, 2007; Evans, 1999; Zhao et al., 2008). This study is devoted to the numerical solution of the Diffusion-Convection equation in the form:

Subject to initial and dirichlet boundary conditions:

where, , and v₁, v₂ are constants that present the convective velocities and diffusivities with respect to x and y directions, respectively. The grid points are given by:

The convection-diffusion Eq. 1 above occurs in many practical problems in which the diffusion coefficient is very small compared to the velocity field which drives the convection, precisely the case which is most difficult to solve accurately. The solution then contains many scales composed of a complex collection of boundary and interior layers. It has been realized that if the convection becomes very large as compared to the diffusion, then second order approximations of the convection term with a large spatial step size give rise to oscillations in the computed solution. Attempts have been made to overcome this difficulty either by the use of an unrealistically small grid or by modifying the method.

Various numerical finite difference schemes have been proposed to solve convection-diffusion problems approximately. Most of these schemes are either first-order or second-order accurate in space and have poor quality for convection dominated flows if the mesh is not sufficiently refined. Higher order discretizations are generally associated with large (non-compact) stencils which increase the band-width of the resulting matrix and lead to a large number of arithmetic operations, especially for higher dimensional problems. To obtain satisfactory higher order numerical results with reasonable computational cost, there have been attempts to develop Higher Order Compact (HOC) schemes, which utilize only the grid nodes directly adjacent to the central node.

Usual finite element methods typically produce approximate solutions with large, non-physical oscillations unless either the mesh width h is globally small with respect to the diffusion coefficient or enough is known about the exact solution to generate Shishkin-like meshes which are locally small with respect to in all transition regions in a very precise sense. Thus, various stabilizations have proven to be essential computational tools.

Many scientists are interested in the convection-dominated case, which may come from a linearized Navier-Stokes equation with high Reynolds number, the drift-diffusion equations of semiconductor device modeling and the Black-Scholes equation from financial modeling. Due to the small diffusion, the solution to Eq. 1 has singularities in the form of boundary or interior layers. The standard numerical approximation, e.g., standard finite element method, on quasi-uniform grids will yield nonphysical oscillatgions unless the mesh is fine enough to capture the layers. To obtain a robust numerical approximation for convection-dominated problems, one approach is to modify the discretization of the convection term while keeping the underlying uniform or quasi uniform grids unchanged (Spotz, 1995; Golbabai and Arabshahi, 2011). Moreover the accuracy depends crucially on the uniformity of the grids which are away from the singularity. On the other hand, we further show that the accuracy depends crucially on the uniformity of the grids which are away from the singularity. The accuracy of the approximation is very sensitive to the perturbation of grid points in the region where the solution is smooth but, in contrast, it is robust with respect to the perturbation of properly adapted grid points inside the boundary or interior layers. Up to now, various numerical methods have been suggested for solving Eq. 1. Many authors are mainly interested in the case |||, ||>>v₁, v₂. Since the diffusion coefficient is often much smaller than the transport velocity, it is very difficult to simulate the Diffusion-Convection problems numerically. The central finite difference and classical Galerkin methods may introduce non-physical oscillations into numerical solutions. Many schemes have been used in the simulations of these problems and have a large number of successful flow simulations. For high Péclet numbers, numerical solutions often exhibit an oscillatory behavior, one possibly chooses a grid-size h small enough to avoid oscillations. However, this is sometimes unpractical and costs much time in the process of iteration, because one has to solve a large linear system that has too many unknowns and may be unfeasible in higher dimensional cases.

Therefore, for the Diffusion-convection problems, there is considerable interest in constructing accurate and non-diffusive schemes, which overcome or reduce numerical oscillations with as little, computing cost as possible and which preserve numerical stability. Both explicit and implicit methods for solving the Diffusion-convection equation are not efficient because the explicit schemes are usually very time consuming due to the stability restriction and implicit methods are un-conditionality stable, thus allowing a large time step. Krylov subspace methods are one of the widely used and successful classes of numerical algorithms for solving large and sparse systems of algebraic equations but the speed of these methods are slow for problems which arise from typical applications. In order to be effective and obtaining faster convergence, these methods should be combined with a suitable preconditioner. The rate of convergence generally depends on the condition number of the corresponding matrix. Since the preconditioner plays a critical role in preconditioned Krylov subspace methods, many preconditioners have been proposed and studied amongst the ADI preconditioner. In this study, we accomplish a comprehensive study for different preconditioners in combination with Krylov subspace methods for solving linear systems arising from the compact high-order approximations. The resulting block tri-diagonal linear system of equations is solved by using Krylov subspace methods. The outline of the study is as follows:

In Section 2, we discuss some physical example of convection diffusion equations and in section 3, we briefly introduce Krylov subspace methods and in Section 4, we consider some available preconditioned techniques. In Section 5, we consider Diffusion-convection problem arising from the compact high-order approximations (Jain et al., 1992). We present the results of our comparative study in the final section.

BEHAVIOR OF CONVECTION-DIFFUSION EQUATION

The convection-diffusion equations are widely used in science and engineering as mathematical models for computational simulations, such as in oil reservoir simulations, transport of mass and energy and global weather production, in which an initially discontinuous profile is propagated by diffusion and convection, the latter with a speed of .

The convection-diffusion equation is a combination of the diffusion and convection (advection) equations and describes physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to two processes: Diffusion and convection. The general equation is:

Where:

The convection-diffusion equation is a relatively simple equation describing flows, or alternatively, describing a stochastically-changing system. Therefore, the same or similar equation arises in many contexts unrelated to flows through space:

where, M is the momentum of the fluid (per unit volume) at each point (equal to the density ρ multiplied by the velocity v), μ is viscosity, P is fluid pressure and f is any other body force such as gravity. In this equation, the term on the left-hand side describes the change in momentum at a given point; the first term on the right describes viscosity, which is really the diffusion of momentum; the second term on the right describes the adjective flow of momentum; and the last two terms on the right describes the external and internal forces which can act as sources or sinks of momentum.

KRYLOV SUBSPACE METHODS

Consider the linear system Ax = b, where A is a large sparse non-symmetric matrix. Let x_o present an arbitrary initial guess to x and r_o = b-Ax_o be a corresponding residual vector. An iterative scheme for solving Eq. 3) is called a Krylov subspace method if for any choice of ω; it produces approximate solutions of the form x = x_o+ω. The subspace k_m is the Krylov subspace .

In section 4, we solve our problems with well-known Krylov subspace methods such as generalized minimal residual method GMRES (m), quasi minimal residual method (QMR), BiCG, Conjugate gradient squared method (CGS) and Bi-Conjugate Gradient stabilized (BiCGSTAB) method that for more complete explanation refer to (Barrett et al., 1994; Saad, 1995; Van der Vorst, 2003).

PRECONDITIONED TECNIQUES

The convergence rate of iterative methods highly depends on the eigenvalue distribution of the coefficient matrix. A criterion for the width of the spectrum is the Euclidean condition number, that is for SPD matrices is:

With , the distance to the exact solution x* in the i^th iteration is bounded by:

The right hand side of (3.2) increases with growing condition number. Hence, lower condition number usually accelerates the speed of convergence. Hence we will attempt to transform the linear system into another equivalent system in the sense that it has the same solution, but has more favorable spectral properties. A preconditioner is a matrix that effects such as a transformation (Bruaset, 1995). If the preconditioner be as M = M₁M₂ then the preconditioned system is as:

The matrices M₁ and M₂ are called the left and right preconditioners, respectively. Now, we briefly describe preconditioners that we use for solving linear systems and let us take A matrix arising from fourth-order approximations that is block tri-diagonal.

Preconditioner based on relaxation technique: Let A = D+L+U such that D, L and U are diagonal, lower and upper triangular block matrices, respectively. A splitting of the coefficient matrix is as A = M-N where the stationary iteration for solving a linear system is as:

If the preconditioner M is defined as M = D, then this preconditioner is called Jacobi. Also, if M is defined as:

Then, we have SOR preconditioner where for ω = 1, we have Gauss-Seidel preconditioner. If M is defined as:

We get SSOR preconditioner. In the above notation, ω is called the relaxation parameter. We have chosen matrix M in Jacobi, G-S and SOR methods as a left preconditioner and in SSOR preconditioner, we have chosen:

As a left preconditioner and M2 = D^-1 (D+ωU) as a right preconditioner. Also, we take:

(Bruaset, 1995) (Fig. 1-3).

ADI preconditioner: Let A = H+V and matrix A is in the form:

where, A_i = tridiag (a_1i, b_1i, c_1i), B_i = tridiag (a_2i, b_2i, c_2i) and C_i = tridiag (a_3i, b_3i, c_3i) of order NxN where H and V are given in the form, V = (0.5 B_i, a_1i, c_1i, a_3i, c_3i).

The alternative direct H = (0.5 B_i, b_3i, b_1i) ion implicit method (In’t Hout and Welfert, 2009; Ma and Saad, 1993; Peaceman and Rachford Jr., 1955) for solving the linear system Ax = b is in following form:

The ADI preconditioner is defined as M = (H+t_jI) (V+r_jI) and M₁ = (H+t_jI), M₂ = (V+r_jI) where Parameters r₁ and r₂ are acceleration parameters. Young (1971) and Varga (1962) proved that the optimum value for r₁ and r₂ is where α≤u_i, v_i≤β and u_i, v_i are eigen-values of matrices H and V respectively.

BLAGE preconditioner: The block alternating group explicit (BLAGE) method was originally introduced as analogue of the Alternating Group Explicit (AGE) method (Bhuruth and Evans, 1997; Mohanty, 2007). The BLAGE uses fractional splitting technique that is applied in two half steps on linear systems with block tri-diagonal matrices of order N²xN² and in the form:

where, A_i, B_i and C_i are tri-diagonal matrices of order NxN. The splitting of matrix A is sum of matrices G₁ and G₂ in which A = G₁+G₂ where G₁ and G₂ are of the form:

and:

For odd values of n where:

The BLAGE preconditioner is as M = (G₁+ω₁I)(G₂+ω₂I) that M₁ = (G₁+ω₁I) and M₂ = (G₂+ω₂I) where ω₁ and ω₂ are optimal iteration parameters. We have experimentally chosen the relaxation parameter and where and, so that we will have the minimum condition number (Fig. 4).

NUMERICAL ILLUSTRATIONS

In this section, in order the validity and effectiveness of the proposed scheme; have considered the numerical solution of Diffusion-convection Eq. 1 with using preconditioned methods (Table 1 and 2). The computations have been done on a P.C. with Core 7 Pue 3.40 GHz and 8 GB RAM. We consider Eq. 1. in the region. The initial conditions is given by:

The analytical solution is of the form:

The initial condition is a Gaussian pulse centered at (0.5, 0.5) having height 1. The boundary conditions have been taken from the analytical solution. We choose v₁ = v₂ = 1, v₁ = v₂ = 0.8, θ = 0.01 and t = 10. We discretized Eq. 1 using compact high-order approximation.

CONCLUSION

High-order approximation are designed by the need to produce more stable schemes which are efficient with respect to the operation number and that do not experience difficulties near boundaries. A high-order compact scheme in combination preconditioner was applied successfully to Diffusion- Convection equation. We study comparison of different preconditioners in combination Krylov subspace methods. The numerical results which are given in the previous section demonstrate the good accuracy of this scheme and efficiency of preconditioned Krylov subspace methods. We got to this conclusion that the ADI preconditioner is effective for model problems rather than other. So, we propose using ADI preconditioner in combination with Krylov subspace methods for solving non-symmetric systems because this preconditioner needs to less computing time and have the less iteration number than other. Also, we propose the BiCGSTAB method because of the need to less iteration number, simplicity in implementation, flat convergence and to save in computational time.

REFERENCES