Fourth order compact schemes for variable coefficient parabolic problems with mixed derivatives

# Fourth order compact schemes for variable coefficient parabolic problems with mixed derivatives

Shuvam Sen Department of Mathematical Sciences, Tezpur University,
Tezpur 784028, Assam, INDIA
###### Abstract

In this article, we have developed a higher order compact numerical method for variable coefficient parabolic problems with mixed derivatives. The finite difference scheme, presented here for two-dimensional domains, is based on fourth order spatial discretization. The time discretization has been carried out using using second order Crank-Nicolson. The present scheme shows good dispersion relation preserving property and has been thoroughly investigated for stability. The discrete Fourier analysis shows that the method is unconditionally stable. The fact that the method has been particularly developed for parabolic equations with mixed derivatives makes it suitable for solving incompressible Navier-Stokes (N-S) equations in irregular domains. To verify the proposed method, several problems with exact and benchmark solutions has been investigated. The proposed compact discretization has been extended to tackle flows of varying complexities governed by the two-dimensional unsteady N-S equations in domain beyond rectangular. The results show good agreement for all the problems considered.

###### keywords:
2D parabolic differential equations, mixed derivatives, Navier-Stokes equations, irregular domains
journal: Journal of Computational Physics

## 1 Introduction

### 1.1 Problem formulation

Let be a rectangular domain in and be the time interval with . The paper describes a compact fourth order numerical method for solving parabolic partial differential equation (PDE) of the form

 ⎧⎨⎩∂tϕ(X,t)+Aϕ(X,t)=s(X,t),(X,t)∈Ω×(0,T]ϕ(X,0)=ϕ0(X),X∈Ωb1(X,t)ϕ+b2(X,t)∂nϕ=g(X,t),X∈∂Ω,t∈(0,T] (1)

for the unknown transport variable defined over . Here is a variable coefficient partial differential operator defined as

 [Aϕ](X,t) = −α1(X,t)[∂xx]ϕ(X,t)−β(X,t)[∂xy]ϕ(X,t)−α2(X,t)[∂yy]ϕ(X,t) (2) +c1(X,t)[∂x]ϕ(X,t)+c2(X,t)[∂y]ϕ(X,t)+d(X,t)ϕ(X,t)

with . Further the coefficients , , , , , and forcing function together with and are assumed to be sufficiently smooth. Here the only additional restriction we have on the Eq. (1) is the positive definiteness of the diffusion matrix. This is equivalent to , and . and are arbitrary coefficients describing the boundary condition as a Dirichlet, Neumann, or Robin type in the boundary normal direction .

The generalized convection-diffusion equations with mixed derivatives, given in Eq. (1), arise in many applications. For example Heston equation which financial mathematician use for option pricing in stochastic volatility models hou_fou_10 (); dur_fou_12 (). We also note the occurrence of such differential equations in mathematical biology hou_wel_07 () and also in numerical mathematics when coordinate transformations are applied to convection-diffusion equations on non-rectangular domains mck_wal_wil_96 (); pan_kal_dal_07 (). Such transformations allow us to work on simple rectangular domains or uniform grids although the original problem is considered on complex domains or with non-uniform grids which may be necessary to account for the stiff behaviour of solution in some part of the domain.

### 1.2 Prior work

In the mathematical literature there exist various different approaches to approximate solutions of PDE. The most popular one amongst these approaches, that has been historically used by the computational fluid dynamics (CFD) community, is the finite difference (FD) method. A higher order compact (HOC) finite difference scheme is one which in two-dimensions employ a nine-point computational stencil using the eight directly adjacent nodes of the reference grid point and offer an accuracy of order four or higher. These schemes lead to a system of equations resulting in a coefficient matrix with smaller bandwidth as compared to non-compact schemes. Apart from solving convection-diffusion equation, different compact schemes have been used successfully to solve non-linear incompressible Navier-Stokes (N-S) equations in all three forms viz. the stream function-vorticity pan_kal_dal_07 (); spo_car_95 (); li_tan_for_95 (); kal_dal_das_02 (); kal_das_dal_04 (); wan_zho_zha_06 (); ge_cao_11 (); sen_13 (), the primitive variables kal_sen_07 (); tia_lia_yu_11 () and the biharmonic ben_cro_fis_05 (); gup_kal_05 (); tia_yu_11 (); sen_kal_gup_13 () formulations.

Here it is worthwhile to point out that although a plethora of compact schemes spo_car_95 (); kal_dal_das_02 (); kal_das_dal_04 (); wan_zho_zha_06 (); ge_cao_11 (); sen_13 (); gup_man_ste_84 (); spo_car_01 (); kar_zha_04 (); you_06 (); tia_ge_07 (); tia_11 () have been developed for convection-diffusion equation but few can tackle the generalized one as specified in Eq. (1). Of course a handful of contributions having linkage to the equation can be found in the works of Fournié and Karaa fou_kar_06 (), Pandit et al. pan_kal_dal_07 (), Karaa kar_07 () and Düring and Fournié dur_fou_12 (). Fournié and Karaa fou_kar_06 () in 2006 derived a fourth order compact FD scheme for a two-dimensional (2D) elliptic PDE with mixed derivative by considering the PDE itself as an auxiliary relation. But in their work they restricted and also considered to be a constant with . Application of this approach to more general problems such as Eq. (1) is not straightforward and may not be possible. This is accentuated in the work of Pandit et al. pan_kal_dal_07 () in 2007, where the authors tried to derive a compact approximation of parabolic equation. Here the authors were ultimately constrained to work with situations where the mixed derivative is absent. Karaa kar_07 (), also in 2007, proposed a fourth order compact FD scheme for solving 2D elliptic and parabolic equations with mixed derivative having variable coefficient by using polynomial approximation, but was again limited by the choice of . In 2012, Düring and Fournié dur_fou_12 () used a compact scheme, having fourth order accuracy in space and second order accuracy in time, for option pricing using Heston model. In that manuscript the authors using a variable transformation arrived at a system with . Note that such transformations are not always certain for a parabolic PDE with variable coefficients.

Although HOC approximation of Eq. (1) is yet to be established, different splitting schemes and their alternating direction implicit (ADI) implementation can be found in the literature. Among them the works of in’t Hout and Foulon hou_fou_10 (), in’t Hout and Welfert hou_wel_07 (), Mckee et al. mck_wal_wil_96 () and the references therein deserves special mention. Recently Martinsson mar_13 () has designed a composite spectral collocation scheme for the equation with smooth solutions. These works clearly establishes importance of generalized parabolic equations. Here we will like to highlight that in many applications it is desirable to use higher order numerical methods to obtain accurate solution. Of late Sen sen_13 () has developed a new family of implicit HOC schemes for unsteady convection-diffusion equation with variable convection coefficient. The schemes, where transport variable and its first derivatives are carried as the unknowns, combine virtues of compact discretization and Padé approximation. In this manuscript we generalize this philosophy to propose a HOC formulation for variable coefficient parabolic problem with mixed derivative. We are able to obtain a new compact stable scheme with truncation error of order four in space and two in time. To the best of our knowledge compact schemes for generalized parabolic 2D convection-diffusion equations is not available in the literature and in this paper we intend to address the same. The schemes thus developed have been tested for their stability by using discrete Fourier analysis and shows better phase and amplitude error properties. We also plan to augment the scheme thus developed to solve incompressible N-S equations in irregular domains.

### 1.3 Outline of this paper

The paper is organized as follows. In the first part, we present the HOC difference scheme for Eq. (1) and its extension to the 2D incompressible Navier-Stokes equations. In the second part, we apply this method to numerical computation to validate our approach. We consider some boundary layer and local singularity problems, as well as the classical impulsively started flow past circular problem as numerical examples to demonstrate the accuracy, effectiveness and efficiency of the present method. Finally, concluding remarks have been presented.

## 2 Fourth order compact schemes for parabolic problem

### 2.1 Spatial compact discretization

We begin by briefly discussing the development of HOC formulation for the steady state form of equation (1), which is obtained when , , , , , , and are independent of . Under these conditions, equation (1) becomes

 {Aϕ(X)=s(X),X∈Ωb1(X)ϕ+b2(X)∂nϕ=g(X),X∈∂Ω (3)

For simplicity we assume . In order to obtain a compact spatially fourth order accurate discretization we lay out a grid , with for and for . Consider the following approximations for second order space derivatives appearing in equation (3)

 ∂xxϕi,j=2δ2xϕi,j−δxϕxi,j+O(h4), (4)
 ∂yyϕi,j=2δ2yϕi,j−δyϕyi,j+O(k4), (5)
 ∂xyϕi,j=δxϕyi,j+δyϕxi,j−δxδyϕi,j+O(h2k2). (6)

Here , , and are usual central difference operators and denote the approximate value of at a typical grid point . Detailed derivation of the above discretizations can be found in sen_13 (). We thus obtain an approximation for equation (3) on a nine point stencil as

 Ah,kϕi,j=si,j (7)

where the discrete operator is defined as

 Ah,kϕi,j = (−2α1i,jδ2x−2α2i,jδ2y+βi,jδxδy+di,j)ϕi,j (8) +(α1i,jδx−βi,jδy+c1i,j)ϕxi,j+(α2i,jδy−βi,jδx+c2i,j)ϕyi,j.

The finite difference operator given above depends on the three grid functions , and . If we are interested in an approximation depending only on , then we need compatible fourth order approximations for space derivatives and in terms of . This is accomplished by using Padé approximations

 (I+h26δ2x)ϕxi,j=δxϕi,j, (9)

and

 (I+k26δ2y)ϕyi,j=δyϕi,j. (10)

Note that compared to standard HOC formulation kal_dal_das_02 (), we are not required to approximate the derivatives of the convection coefficients , and forcing function . The equation (7) can be viewed as a banded system with only nine non zero diagonals; of course drawback of requiring to approximate and separately using (9) and (10) respectively remain.

### 2.2 Modified wave number analysis

A detailed wave number analysis of the fourth order compact approximation for was carried out by Sen sen_13 (). Here we will like to examine the characteristic of the newly used fourth order compact approximation for the mixed derivative . We begin by considering the test equation

 ϕxy=s(x,y). (11)

Consider the trial function where and are the wave numbers corresponding to and directions respectively. Then the exact characteristic of the equation (11) is

 λExact=−κ1κ2. (12)

Using the compact approximation given in Eq. (6) in Eq. (11) we get the corresponding discretized equation as

 δxϕyi,j+δyϕxi,j−δxδyϕi,j=si,j (13)

together with the Padé approximation given in Eqs. (9) and (10). It is easy to see that the above fourth order compact discretization for mixed derivative has characteristic

 λ4OC−M=−sin(κ1h)sin(κ2k)hk[32+cos(κ1h)+32+cos(κ2k)−1]. (14)

Note that the characteristics for the second order accurate central difference approximation and fourth order accurate wide stencil approximation for the mixed derivatives are

 λ2OC=−sin(κ1h)sin(κ2k)hk (15)

and

 λ4OW=−sin(κ1h)sin(κ2k)9hk(4−cos(κ1h))(4−cos(κ2k)) (16)

respectively. To the best of our knowledge no HOC approximation for the mixed derivative is available in litarature. To get a clear idea of the dissipation errors associated with each of the above discretizations the non-dimensional characteristics as a function of corresponding to four different values of have been shown in figure 2. Note that the expressions for characteristics being symmetrical with respect to and the above set of values give comprehensive idea of associated errors. The figure clearly indicates that the fourth order compact discretization discussed here has superior wave resolution property than the other discretization procedures that can be used for mixed derivatives.

### 2.3 Implicit time discretization

The HOC approach developed for the steady case can be extended directly to the unsteady case by simply replacing by in Eq. (3). At grid point and time , the semi-discrete fourth order scheme for the parabolic equation with variable coefficients will be

 ∂tϕi,j(t)+Ah,kϕi,j(t)=si,j(t) (17)

Clearly any time integrator can be used in Eq. (17). Introducing weighted time average parameter such that for , where denote the -th time level, we obtain a family of integrators; for example, forward Euler for , backward Euler for and Crank-Nicholson for . The resulting fully discrete difference scheme for grid point at time level then becomes

 [1+ιδtAh,k]ϕ(n+1)i,j=[1−(1−ι)δtAh,k]ϕ(n)i,j+ιδts(n+1)i,j+(1−ι)δts(n)i,j. (18)

The accuracy of the scheme is , with . The second order accuracy in time is obtained for . On expansion the Eq. (18) reduces to

 [1−ιδt(2α(n+1)1i,jδ2x+2α(n+1)2i,jδ2y−β(n+1)i,jδxδy−d(n+1)i,j)]ϕ(n+1)i,j =[1+(1−ι)δt(2α(n)1i,jδ2x+2α(n)2i,jδ2y−β(n)i,jδxδy−d(n)i,j)]ϕ(n)i,j −(1−ι)δt[(α(n)1i,jδx−β(n)i,jδy+c(n)1i,j)ϕ(n)xi,j+(α(n)2i,jδy−β(n)i,jδx+c(n)2i,j)ϕ(n)yi,j−s(n)i,j] −ιδt[(α(n+1)1i,jδx−β(n+1)i,jδy+c(n+1)1i,j)ϕ(n+1)xi,j +(α(n+1)2i,jδy−β(n+1)i,jδx+c(n+1)2i,j)ϕ(n+1)yi,j−s(n+1)i,j]. (19)

It is clear that for all values of , except and , the difference stencil requires nine points in both -th and -th time levels resulting in what may be called a scheme. The compact stencil emerging in this way has been illustrated in figure 1. and schemes are obtained for and , respectively. The algorithm given by Sen sen_13 () can be used to solve the algebraic system associated with the finite difference approximation given in Eq. (2.3).

### 2.4 Stability analysis for constant coefficients

We carry out the stability analysis for the scheme given by the equation Eq. (18) by using von Neumann method. This is applicable in the special case when the coefficients , , , , and are constant. We consider a solution to the difference equation (2.3) to be

 ϕ(n)i,j=b(n)eIθxieIθyj

where is the amplitude at time level , and , are the phase angles with wavelengths , respectively. Also periodic boundary conditions are assumed. The amplification factor for all is defined as .

###### Theorem 1 (Stability)

The finite difference scheme (18) with constant coefficients is stable, in the von Neumann sense, for if and for the scheme is stable under sufficient condition .

\newproof

pfProof {@proof}[Proof.] In order to prove the above theorem we will require the following two lemmas.

###### Lemma 0

whenever .

Proof of the Lemma 1: This is clearly true if ; otherwise since we get

 |β|<2√α1α2 ⇒ |βABr|<2√α1α2|AB| ⇒ −2√α1α2|AB|<βABr ⇒ α1A2+α2B2−2√α1α2|AB|<α1A2+α2B2+βABr ⇒ (α1|A|−α2|B|)2<α1A2+α2B2+βABr ⇒ 0<α1A2+α2B2+βABr.
###### Lemma 0

The function

 R(θx,θy) = 2(8+cosθx+cosθy−cosθxcosθy)cos(θx/2)cos(θy/2)√(5+cosθx)(2+cosθx)(5+cosθy)(2+cosθy)

is a periodic function of period with respect to both the variables. This function defined over attains (i) global maxima at and and its maximum value is 1, (ii) global minima at and and its minimum value is -1.

We note that Lemma 3 can be arrived at, as a problem of maxima-minima.

Now we proceed to prove the stability theorem.

From equations (9) and (10) we see that

 ϕ(n)xi,j=I3sinθxh(2+cosθx)ϕ(n)i,j, (20)
 ϕ(n)yi,j=I3sinθyk(2+cosθy)ϕ(n)i,j, (21)

which upon substitution into the Eq. (8) yields

 Ah,kϕ(n)i,j=F(α1,α2,β,c1,c2,d,h,k,θx,θy)ϕ(n)i,j (22)

where

 F(α1,α2,β,c1,c2,d,h,k,θx,θy) ≡ (2α1(2−2cosθx)h2+2α2(2−2cosθy)k2−βsinθxhsinθyk+d) +(−α1sinθxh+βsinθyk)3sinθxh(2+cosθx) +(−α2sinθyk+βsinθxh)3sinθyk(2+cosθy) +I(c13sinθxh(2+cosθx)+c23sinθyk(2+cosθy)) ≡ FR+d+IFI.

Thus the amplification factor of the scheme given in Eq. (18) is then

 |G|=∣∣∣1−(1−ι)δtF1+ιδtF∣∣∣ (23) ⇒ |G|2=(1−(1−ι)δt(FR+d))2+(1−ι)2δt2F2I(1+ιδt(FR+d))2+ι2δt2F2I.

Now for the scheme is unconditionally stable

 if−2(FR+d)+((FR+d)2+F2I)δt(1−2ι)≤0 i.e. ifFR≥0for(1−2ι)≤0 i.e. ifFR≥0for0.5≤ι

Here,

 FR = α1h2[4(1−cosθx)−3sin2θx2+cosθx]+α2k2[4(1−cosθy)−3sin2θy2+cosθy] +βhksinθxsinθy[32+cosθx+32+cosθy−1] = α1h2[sin2θx−4cosθx+42+cosθx]+α2k2[sin2θy−4cosθy+42+cosθy] +βhksinθxsinθy[8+cosθx+cosθy−cosθxcosθy(2+cosθx)(2+cosθy)] = 2α1h2sin2(θx/2)5+cosθx2+cosθx+2α2k2sin2(θy/2)5+cosθy2+cosθy +βhksinθxsinθy[8+cosθx+cosθy−cosθxcosθy(2+cosθx)(2+cosθy)] = 2α1h2sin2(θx/2)5+cosθx2+cosθx+2α2k2sin2(θy/2)5+cosθy2+cosθy +2βhksin(θx/2)sin(θy/2)√5+cosθx2+cosθx√5+cosθy2+cosθyR.

By Lemma 3, we see that . Hence by Lemma 2, .

Next, correspond to a growth term in the governing equation (1) and the stability condition requires with independent of smith_book_86 (). From (23) we see that this is obviously true for .

This completes the proof.

### 2.5 Extension to the 2D incompressible Navier-Stokes equations

In this section we discusses applicability and effectiveness of the newly developed scheme in tackling coupled non-linear system of Navier-Stokes (N-S) equations in geometries beyond rectangular. It is worthwhile to point out that this fourth order compact scheme can be also used as a basis for discretizing unsteady incompressible N-S equations. One of the goals of this manuscript is to design schemes suited to simulate the time development of 2D viscous flows of varied physical complexities in different geometrical settings and accurately report fluid-embedded body interaction. Note that to resolve complex flow phenomena accurately, we need to develop HOC schemes for body fitted coordinate system. It must be mentioned here that the idea of constructing a rectangular grid on irregular domains has its inherent limitations; furthermore accurate imposition of boundary condition becomes an issue. This is where the motivation to develop HOC schemes on nonuniform grids comes from. Although the primitive variable form of N-S equations can be discretized we choose to work with stream function-vorticity formulation as it decouples pressure field from velocity calculation. The system of equations presented below can be found elsewhere pan_kal_dal_07 (); sen_kal_gup_13 () but we reproduce them for completeness.

In rectangular cartesian coordinates the formulation is given by

 −ψxx−ψyy=ω, (24)
 ∂tω−1Re(ωxx+ωyy)+(uωx+vωy)=0. (25)

Here is the velocity field and is the vorticity. Also is the Reynolds number based on characteristic length and characteristic velocity of the flow. We consider non-singular mapping

 x=x(ξ,η),y=y(ξ,η). (26)

to transform the physical plane into a computational plane. Under this transformation equations (24) and (25) become

 −˜a1ψξξ−˜e1ψξη−˜b1ψηη+˜c1ψξ+˜d1ψη=˜f1, (27)
 ∂tω−˜a2ωξξ−˜e2ωξη−˜b2ωηη+˜c2ωξ+˜d2ωη=0, (28)

with the velocity components expressed as

 u=1J(ψηxξ−ψξxη), (29)
 v=1J(−ψξyη+ψηyξ). (30)

Details of the coefficients appearing in the above relations are available in pan_kal_dal_07 (); sen_kal_gup_13 (). Although Pandit et al. in their work pan_kal_dal_07 () arrived at the above system of equations but ultimately restricted themselves to the development of HOC schemes for only conformably transformed equations. Note that in case of conformal transformations coefficient of the mixed derivative term in the above equations will vanish. Here in this work we are not restricted by such choices of transformations and hence hold much wider appeal. To the best of our knowledge no fourth order compact formulation to tackle the above generalized N-S system is available in the literature. Both the equations (27) and (28) satisfy positive definiteness of associated diffusion matrix and are perfectly amenable to the discretization procedures detailed in this manuscript. The steady formulation given in Eq. (7) is to be used for Eq. (27) where as the implicit scheme presented in Eq. (18) is to be used for Eq. (28). In this manuscript we restrict ourselves to to have second order time accuracy. Since the gradients of the flow variables are already available at all grid points, necessity to compute separately does not arise. For the coupled nonlinear PDEs discussed above an outer-inner iteration procedure can be adopted. Outer one is the temporal iteration and in the inner iteration transient and steady equations are repeatedly solved till convergence is achieved.

## 3 Numerical Examples

### 3.1 Problem 1: Analytical test case

We consider the two-dimensional parabolic problem

 ∂u∂t−∂2u∂x2+(1−x)(1−y)ex+y∂2u∂x∂y−∂2u∂y2+10x(1−y)∂u∂x−10y∂u∂y=f(x,y,t) (31)

defined over the domain . Following Karaa kar_07 () we use this problem to test the proposed compact scheme and to confirm the theoretical fourth order accuracy in space and second order accuracy in time. In the Eq. (31) the forcing term is selected to satisfy exact solution . The initial and Dirichlet boundary conditions are also taken from this exact solution. We consider three different mesh sizes , and and compute errors with respect to the exact solution using , and norms. The computed results at time and have been presented in table 1. A spatial order of accuracy close to four, for all cases, can be seen in this table. We estimate order of accuracy using the formula where and are the errors with the grid sizes and respectively. Temporal second order accuracy of the scheme is verified in the table 2 where we compute using three different time steps , and keeping unchanged.

### 3.2 Problem 2: Boundary layer in non-rectangular domain

This problem has been considered check the effectiveness of the spatial discretization in resolving boundary layers in non-rectangular domains. Here we tackle an elliptic problem by using the nine point scheme presented in equation (2.3) with . The domain of the problem is given by

 Ω={(x,y):0≤x≤1,0≤y≤11−0.3sin(6x)}.

The elliptic differential equation is

 −ϵ(∂2ψ∂x2+∂2ψ∂y2)−1.8cos(6x)(1−0.3sin(6x))(2−0.3sin(6x))∂ψ∂x+11+y∂ψ∂y=f, (32)

which admits a boundary layer along the top curved surface of the domain . The analytical solution of the problem is

 ψ=ey−x+(1+y)1+1ϵ[1−0.3sin(6x)2−0.3sin(6x)]1ϵ.

A nonuniform grid along -axis with clustering at the top and uniform grid along -axis has been used following the transformation

 x=ξ,y=η+λπsin(πη)1−0.3sin(6ξ)

from plane to plane. Here is a parameter controlling the stretching along the top curved boundary and was considered to be 0.9 in our computations. With this transformation the equation (32) takes the form presented in Eq. (2.3). The physical domain along with a grids have been shown in figure 3. We test our scheme using two small values of and using a fine grid of . The numerical results thus produced have been presented in figure 4(a) and (b) respectively. The results clearly show that our computation has been able to capture the boundary layers in both the cases. For maximum error where as for .

### 3.3 Problem 3: Flow past circular cylinder

Finally we intent to use the scheme developed here to solve Navier-Stokes equations and simulate flow around circular cylinder. This problem serves as an ideal test case for flow past bluff bodies on non-rectangular domains and of late, has quite often being used to validate schemes le_kho_per_06 (); ber_fal_08 (); wan_fan_cen_09 (); sen_sum_sin_10 (); pen_mit_sau_10 () for their accuracy and effectiveness. We assume the cylinder to be of unit diameter placed in an infinite domain and compute flow field for by solving the coupled system of Eqs. (27) and (28). A multi-block CFD grid has been used for computing flow field. Precisely two blocks of grids have been used. Around the cylinder a circular grid of size has been generated by using the transformation Downstream of the flow, to capture the process of vortex shedding, another block of grids of size have been numerically generated. The two sets of blocks have been point matched at the abutting boundary. The grid thus generated has been shown in the figure 5(a). The boundary conditions for the and its derivatives have been derived following the work of Kalita and Sen kal_sen_12-I (). Where as for we follow the boundary conditions given by Kalita and Ray kal_ray_09 (). Tian et al. tia_lia_yu_11 () in their work have detailed various biased finite difference schemes for estimating first derivatives at the boundary. These approximations can be used here as well to approximate derivatives of vorticity field at the boundaries. We computing using . In figure 5(b), we show the time histories of the drag and lift coefficients which establishes the eventual periodic nature of the flow. In table 3 we compare Strouhal numbers, drag and lift coefficients for the flow with the values available in the literature. A good comparison can be seen. We present the temporal evolution of streamlines and vorticity contours over one complete vortex shedding cycle of duration in figure 6. The evolution of von Krmn vortex street can be clearly seen in these figures. Two eddies are shed just behind the cylinder within each period. Figures 6(a) and 6(b) are half a vortex-shedding cycle apart and correspond to peak value of lift coefficient. Figure 6(b) is a mirror image of figures 6(a) and 6(c). The crests and troughs of the sinuous waves in the streamlines reflect the alternatively positive and negative vorticities of the eddies.

## 4 Conclusion

In this paper, we present a higher order compact numerical method for 2D parabolic problems with mixed derivatives. Apart from solving convection-diffusion equation the formulation has been used to simulate unsteady two-dimensional Navier- Stokes equations. The robustness of the scheme is highlighted when it captures not only the boundary layer of a test case but also the process of vortex shedding for incompressible viscous flows on geometries beyond rectangular. Because of its high computational efficiency, our scheme has a good potential for efficient application. To the best of our knowledge compact schemes for parabolic 2D convection-diffusion equations is not available and thus this work is an important addition to the high accuracy solution procedures available in literature.

## References

•  K. J. in’t Hout, S. Foulon, Adi finite difference schemes for option pricing in the heston model with correlation, International Journal of Numerical Analysis and Modeling 7 (2010) 303-320.
•  B. Düring, M. Fournié, High-order compact finite difference scheme for option pricing in stochastic volatility models, Journal of Computational and Applied Mathematics 236 (2012) 4462-4473.
•  K.J. in’t Hout, B.D. Welfert, Stability of ADI schemes applied to convectiondiffusion equations with mixed derivative terms, Applied Numerical Mathematics 57 (2007) 19-35.
•  S. Mckee, D. P. Wall, And S. K. Wilson, An alternating direction implicit scheme for parabolic equations with mixed derivative and convective terms, Journal of Computational Physics 126 (1996) 64-76.
•  S.K. Pandit, J.C. Kalita, D.C. Dalal, A transient higher order compact scheme for incompressible viscous flows on geometries beyond rectangular equation on non-uniform grids, Journal of Computational Physics 225 (2007) 1100-1124.
•  W.F. Spotz, G.F. Carey, High-order compact scheme for the steady stream-function vorticity equations, International Journal for Numerical Methods in Engineering 38 (1995) 3497-3512.
•  M. Li, T. Tang, B. Fornberg, A compact fourth-order finite difference scheme for the incompressible Navier-Stokes equations, International Journal for Numerical Methods in Fluids 20 (1995) 1137-1151.
•  J.C. Kalita, D.C. Dalal, A.K. Dass, A class of higher order compact schemes for the unsteady two-dimensional convection-diffusion equation with variable convection coefficients, International Journal for Numerical Methods in Fluids 38 (2002) 1111-1131.
•  J.C. Kalita, A.K. Dass, D.C. Dalal, A transformation-free HOC scheme for steady convection-diffusion on non-uniform grids, International Journal for Numerical Methods in Fluids 44 (2004) 33-53.
•  J. Wang, W. Zhong, J. Zhang, High order compact computation and nonuniform grids for stream function vorticity equations, Applied Mathematics and Computation 179 (2006) 108-120.
•  Y.B. Ge, F. Cao, Multigrid method based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems, Journal of Computational Physics 230 (2011) 4051-4070.
•  S. Sen, A new family of (5,5)CC-4OC schemes applicable for unsteady Navier-Stokes equations, Journal of Computational Physics 251 (2013) 251-271.
•  J.C. Kalita, S. Sen, The (9,5) HOC formulation for the transient Navier-Stokes equations in primitive variable, International Journal for Numerical Methods in Fluids 55 (2007) 387-406.
•  Z.F. Tian, X. Liang and P. Yu, A higher order compact finite difference algorithm for solving the incompressible Navier-Stokes equations, International Journal for Numerical Methods in Engineering 88 (2011) 511-532.
•  M. Ben-Artzi, J-P Croisille, D. Fishelov and S. Trachtenberg, A pure-compact scheme for the streamfunction formulation of Navier-Stokes equations, Journal of Computational Physics 205 (2005) 640-664.
•  M.M. Gupta, J.C. Kalita, A new paradigm for solving Navier-Stokes equations: streamfunction-velocity formulation, Journal of Computational Physicis 207 (2005) 52-68.
•  Z.F. Tian, P.X. Yu, An efficient compact difference scheme for solving the streamfunction formulation of the incompressible Navier-Stokes equations, Journal of Computational Physics 230 (2011) 6404-6419.
•  S. Sen, J.C. Kalita, M.M. Gupta, A robust implicit compact scheme for two-dimensional unsteady flows with a biharmonic stream function formulation, Computers & Fluids 84 (2013) 141-163
•  M.M. Gupta, R.P. Manohar, J.W. Stephenson, A single cell high order scheme for the convection-diffusion equation with variable coefficients, International Journal for Numerical Methods in Fluids 4 (1984) 641-651.
•  W.F. Spotz, G.F. Carey, Extension of high-order compact schemes to time-dependent problems, Numerical Methods for Partial Differential Equations 17 (2001) 657-672.
•  S. Karaa, J. Zhang, High order ADI method for solving unsteady convection-diffusion problems, Journal of Computational Physics 198 (2004) 1-9.
•  D. You, A high-order Padé ADI method for unsteady convectiondiffusion equations, Journal of Computational Physics 214 (2006) 1-11.
•  Z.F. Tian, Y.B. Ge, A fourth-order compact ADI method for solving two dimensional unsteady convection-diffusion problems, Journal of Computational and Applied Mathematics 198 (2007) 268-286.
•  Z.F. Tian, A rational high-order compact ADI method for unsteady convection-diffusion equations, Computer Physics Communications 182 (2011) 649-662.
•  M. Fournié, S. Karaa, Iterative methods and high-order difference schemes for 2d elliptic problems with mixed derivative, Journal of Applied Mathematics and Computation 22 (2006) 349-363.
•  S. Karaa, High-order difference schemes for 2D elliptic and parabolic problems with mixed derivatives, Numerical Methods for Partial Differential Equations, 23 (2007) 366-378.
•  P.G. Martinsson, A direct solver for variable coefficient elliptic PDEs discretized via a composite spectral collocation method, Journal of Computational Physics 242 (2013) 460-479.
•  G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Third Edition, Oxford University Press, New York, 1986.
•  D. V. Le, B. C. Khoo, J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, J. Comput. Phys. 220 (2006) 109-138.
•  P. A. Berthelsen, O. M. Faltinsen, A local directional ghost cell approach for incompressible viscous flow problems with irregular boundaries, J. Comput. Phys. 227 (2008) 4354-4397.
•  Z. Wang, J. Fan, K. Cen, Immersed boundary method for the simulation of 2D viscous flow based on vorticity-velocity formulations, J. Comput. Phys. 228 (2009) 1504-1520.
•  T. K. Sengupta, V. K. Suman, N. Singh, Solving NavierStokes equation for flow past cylinders using single-block structured and overset grids, J. Comput. Phys. 229 (2010) 178-199.
•  Y. F. Peng, R. Mittal, A. Sau, R. R. Hwang, Nested Cartesian grid method in incompressible viscous fluid flow, J. Comput. Phys. 229 (2010) 7072-7101.
•  J.C. Kalita, S. Sen, The biharmonic approach for unsteady flow past an impulsively started circular cylinder, Commun. Comput. Phys. 12 (2012) 1163-1182.
•  J.C. Kalita, R.K. Ray, A transformation-free HOC scheme for incompressible viscous flows past an impulsively started circular cylinder, J. Comput. Phys. 228 (2009) 5207-5236. Figure 2: Comparison of nondimensional characteristics for three different schemes with the exact one at (a) κ2k=0.5, (b) κ2k=1.0, (c) κ2k=1.5, (d) κ2k=2.0. Figure 5: Problem 3: (a) Grid used for computing flow field. (b) Variation of drag and lift coefficients for Re=300. Figure 6: Problem 3: The streamlines and vorticity contours for three successive instants of time (a) t=T, (b) t=T+T0/2, (c) t=T+T0.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters   