An adaptive mesh strategy for singularly perturbed convection diffusion problems

Applied Mathematical Modelling xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

An adaptive mesh strategy for singularly perturbed convection diffusion problems Vivek Kumar a,⇑, Balaji Srinivasan b a b

Delhi Technological University, Dept. of Applied Mathematics, Delhi 110042, India Indian Institute of Technology, Dept. of Applied Mechanics, Delhi 110016, India

a r t i c l e

i n f o

Article history: Received 25 October 2013 Received in revised form 9 August 2014 Accepted 1 October 2014 Available online xxxx Keywords: Central finite difference schemes Convection dominated singularly perturbed problems Layer-adaptive meshes

a b s t r a c t In this paper, a new adaptive mesh strategy has been developed for solving convection dominated, convection–diffusion singularly perturbed problems (SPP) using second order central difference schemes. Our strategy uses a novel, entropy-like variable as the adaptation parameter for convection diffusion SPP. Further, unlike the popular layer adapted meshes mainly by Bakhvalov (B-type) and Shishkin (S-type), no pre-knowledge of the location and width of the layers (boundary as well as interior) is needed. The method is completely free of arbitrary perturbation parameters  (robust) and results in oscillation free solutions to a range of convection–diffusion SPP. Numerical results for various test problems: linear (boundary layers with and without turning points), nonlinear and systems of coupled equations are presented. This method is expected to aid in robust computations of convection dominated SPP. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Singularly perturbed differential equations, which have non smooth solutions with singularities related boundary layers, and their numerical solutions have been quite an active field of research over the last few decades. Starting from late 60’s by Bakhvalov [1], a large number of numerical techniques have been given by many researchers. Well developed techniques are now available for the numerical solutions of outside the layers regions but solving the problem inside the layers region is limited and still under investigation. To deal with such problems, robust layers resolving methods: methods which generates numerical approximation at each point in the domain for arbitrarily small values of the singular perturbation parameter  are often needed in the real world applications. Layer adapted meshes have been first proposed by Bakhvalov for reaction diffusion equations and later on by Gartland [2] and others for convection diffusion equations. In early 90s, special piecewise uniform meshes have been introduced by Shishkin [3]. Because of simple structure of Shishkin meshes they have attracted much attention and are widely used for numerically approximation of SSP. The limitation with Shishkin meshes is the requirement of a priori knowledge of the location of the layer regions. A detailed survey about layer adapted meshes for convection diffusion equations can be found in [4]. It is well known that whenever central finite difference schemes (second order or higher order central schemes, nonmonotone methods) are applied to solve convection diffusion singularly perturbed problems numerically on uniform meshes, unphysical oscillations are observed in the numerical solution and their magnitude increases in the layers ⇑ Corresponding author. E-mail addresses: [email protected] (V. Kumar), [email protected] (B. Srinivasan). http://dx.doi.org/10.1016/j.apm.2014.10.019 0307-904X/Ó 2014 Elsevier Inc. All rights reserved.

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

2

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx

(boundary and interior) regions. The presence of these large unphysical oscillations in the approximated solution shows that central finite difference operators on the uniform meshes are unstable below a certain value of the singular perturbation parameter. To eliminate these oscillations while retaining the order of accuracy, one needs either a fine mesh at the layer. This maybe done either via uniformly fine meshing or via an adaptive mesh strategy. The former strategy increases significantly in computational cost as the perturbation parameter  decreases. In this paper, we exploit such oscillations indirectly to locate the layers regions and then generate the layer adaptive meshes, via an entropy production operator which is positive automatically whenever the solution is unphysical, which effectively corresponds to insufficient resolution. We restrict ourselves to the class of SPPs which correspond to second order boundary value problems with a non-zero coefficient of the convection (first derivative) term. Section 2 has been dedicated to adaptive mesh strategies and various numerical results have been discussed in Section 3 followed by the conclusion. 2. Adaptive mesh strategy 2.1. Linear convection diffusion problems We consider the following one dimensional steady-state convection diffusion problem as

LuðxÞ ¼ u00 ðxÞ  aðxÞu0 ðxÞ þ bðxÞuðxÞ ¼ f 2 ð0; 1Þ;

uð0Þ ¼ a;

uð1Þ ¼ b

ð2:1Þ

with a small positive perturbation parameter . The above equation physically represents the steady state of transport of a passive scalar by a background velocity. For instance, uðxÞ could represent the temperature and aðxÞ could represent the (possibly spatially varying) background velocity and  is the diffusion coefficient. As  becomes smaller, convective processes start dominating over the diffusive processes [5]. Without loss of generality, we consider bðxÞ P 0 and aðxÞ can change the sign (as in the case of turning point problem) in the given domain ½0; 1. For aðxÞ P c > 0, the above problem has unique solution and exhibit a boundary layer of OðÞ at x ¼ 0. It was proved by Kellogg and Tsan [6] that u and its derivatives up to an arbitrary prescribed order q (depending on the smoothness of the data) can be bounded by

juðkÞ ðxÞj 6 Cf1 þ k ecx= g;

for k ¼ 0; 1; . . . ; q and x 2 ½0; 1;

where C denotes a generic positive constant independent of  and the number of mesh points used. We now define an auxiliary equation – the entropy production equation – by multiplying with an appropriate test function. From the case of scalar hyperbolic conservation laws [7], it is well known that for scalar conservation laws, u2 is always an appropriate entropy variable and, therefore, 2uðxÞ is a suitable multiplying test function. On multiplying (component wise), we obtain

Lu  2u ¼ f  2u;

ð2:2Þ

after performing the calculation we can write the above Eq. (2.2) as 2

ðS00  2ðu0 Þ Þ  aS0 þ 2bS  2uf ¼ 0 where S ¼ u2 :

ð2:3Þ

Rewriting the above equation, we get 2

S00  aS0  2uf ¼ 2bS  2ðu0 Þ :

ð2:4Þ

We label the linear operator on the left hand side (LHS) in the above Eq. (2.4), as the entropy production operator with analogy to similar operators in hyperbolic conservation laws [7]. The continuous operator should obviously be negative for all x 2 ½0; 1 (as right hand side (RHS) is always negative for all x 2 ½0; 1).

Fig. 1. Ghost points at the left and the right boundaries.

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx

3

Main Idea: As we know that if we solve the discrete analogue of Eq. (2.1) computationally using central difference method, we get non physical oscillations inside and near the boundary layer region. Similarly if we calculate the discrete analogue of the left hand side equation in (2.4) using the same central difference operator by taking Si ¼ u2i ; where ui is the central difference computed solution for Eq. (2.1), we observe that LHS is negative wherever the solution is smooth enough and positive where we have boundary layers (or oscillation in the computed solution of Eq. (2.1)). After investigation, we have 2

found that if we write the RHS term 2bu  2ðu0 Þ2 of Eq. (2.4) in difference operator at the ith mesh point, we get

   ui  ui1 uiþ1  ui 2bi ðui ui1 Þ  2 xi  xi1 xiþ1  xi

ð2:5Þ

by applying a combination of forward and backward difference operator for u0 and ðui1 ui Þ for u2 which are of opposite sign wherever the oscillations occur and hence results in positive value of the difference operator as in Eq. (2.4) for the oscillatory numerical solution ui . We exploit this property of the entropy production operator to precisely determine regions of insufficient resolution in what follows. 2.2. Adaptive mesh generation algorithm To generate the adaptive mesh for the convection dominant SPP, we follow the following algorithm: Step 1: Calculate the entropy with a minimum number of initial uniform mesh points (minimum 3 points in our case). Step 2: Locate the mesh points where entropy S is positive and then choose the mesh point with the maximum entropy. Step 3: Add one mesh point on the left and the right side of the location found in step 2. Step 4: go to Step 1 till all the entropy becomes non positive over the complete mesh. The resulting mesh is our adapted mesh which is represented by N a . 2.3. Entropy calculation at the boundaries For a general Dirichlet boundary conditions (DBCs)

uðxL Þ ¼ uL ; uðxR Þ ¼ uR ; we introduce a ghost point say x0 at the left of the node xL and xnþ1 on the right of the node xR as shown in the Fig. 1, then we calculate

 d2 u1  a1 D0 u1 þ b1 u1 ¼ f 1 ; 0

2

 d un  an D un þ bn un ¼ f n :

ð2:6Þ ð2:7Þ

As we already know the computed values of u1 ; u2 and then we calculate u0 at the ghost point x0 . Then we can find the entropy at the node xL as usual (discussed above in Eq. (2.4)). Similarly we can calculate the entropy at the right boundary node xR . 2.4. Finite difference operators for a non-uniform mesh Let u be any mesh function defined on a given non uniform (or adaptive mesh N a ) mesh xi ; i ¼ 0; 1; . . . N; xi 2 ð0; 1Þ. The first order forward, backward and centered finite difference operators Dþ ; D ; D0 can be defined as:

uðxiþ1 Þ  uðxi Þ ; xiþ1  xi uðxi Þ  uðxi1 Þ ; D uðxi Þ ¼ xi  xi1 uðxiþ1 Þ  uðxi1 Þ D0 uðxi Þ ¼ xiþ1  xi1

Dþ uðxi Þ ¼

ð2:8Þ

and the second order centered difference operator d2 is defined by

d2 uðxi Þ ¼

2ðDþ uðxi Þ  D uðxi ÞÞ : xiþ1  xi1

ð2:9Þ

3. Numerical results In this section we solve different types of various examples to validate the proposed adaptive method. Without loss of generality, now onward we represent uðx; Þ as u. Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

4

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx

3.1. One dimensional (1D) linear convection diffusion problems Test case 1 [8]: We consider 1D linear convection–diffusion problem as

u00  2u0 ¼ 0; x 2 ð0; 1Þ; 0 <   1;

ð3:10Þ

with boundary conditions

uð0Þ ¼ 1 and uð1Þ ¼ 0;

ð3:11Þ

The exact solution is

u¼

expð2x=Þ  expð2=Þ : 1  expð2=Þ

ð3:12Þ

This problem has regular boundary layers of width OðÞ at x ¼ 0. When we apply central difference scheme on uniform mesh for this problem, we get

  h

2

þ

   1 2  1 ui1 þ 2 ui þ uiþ1 ¼ 0:  2 h h h h

ð3:13Þ

On solving the corresponding linear difference equation as given in the Eq. (3.13), we get

ui ¼

si  sN 1  ð1=NÞ ; where s ¼ : 1 þ ð1=NÞ 1  sN

ð3:14Þ

Where h ¼ 1=N with N number of mesh points. Now if we take  < h then for  sufficiently small, we have s < 0 and si is positive or negative depending on whether i is even or odd. In that case the exact solution of the difference Eq. (3.13) will show some non physical oscillations, which increases unboundedly for N even and lim  ! 0. This behavior is completely different from the exact solution as given above (3.12), which is smooth and decreases from 1 to 0. The numerical approximated solution on a uniform mesh for  ¼ 105 and N ¼ 500 is shown below in Fig. 2(a). Fig. 2(b) shows that the boundary layer is resolved very well with adaptive mesh containing only 30 points. Table 1 shows the maximum error computed for the corresponding adaptive meshes using central (CS) and upwind schemes (US) for various values of . Table 1 illustrates the uniformity in  of the method. The proposed method also satisfies  6 2c N 1 as is generally the case for discretizations of convection-dominated a problems. As discussed in [4], the central difference discrete operator say L is ðl1 ; l1 Þ stable for any mesh function u (computed solution in our case) given on any arbitrary mesh x (adaptive mesh in our case) with

Computed solution on uniform mesh

Computed solution on adaptive mesh

1

1.2 exact sol. computed sol.

0.8

adaptive mesh

1 0.6 0.8

0.4

0.6 u(x)

u(x)

0.2 0

0.4

−0.2 −0.4

0.2

−0.6 0 −0.8 −1

−0.2 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0

Fig. 2. Results for the test case 1 for

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

 ¼ 105 .

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

5

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx

0.24

1 exact sol. central

0.22

upwind

0.8 0.2

central scheme 0.6

0.18

u(x)

Max. error

upwind scheme

0.4 0.16 0.2 0.14

0

0.12 0

0.1

0.2

εN

0.3

0.4

0.5

0

0.002

0.004

0.006

0.008

0.01

x

a

Fig. 3. Results for the test case 1.

Table 1 Maximum errors for central and upwind schemes using same adaptive mesh having points N a for test case 1.



Central scheme

Upwind scheme

Na

24

0.1353352

0.2209561

07

26

0.1353352

0.2218050

09

28

0.1353352

0.2218075

11

210

0.1353352

0.2218075

13

212

0.1353352

0.2218075

15

214

0.1353352

0.2218075

17

216

0.1353352

0.2218075

19

218

0.1353352

0.2218075

21

jjujj1;x 6

a 1 81 NX hi þ hiþ1 hi j½Lui j where hi ¼ ; 4c i¼1 2

hi ¼ xi  xi1 :

ð3:15Þ

provided

    N Y  1  1   a  þ  6 1=4   i¼1 hi ai1 2 hi ai 2  and hi 6 chj for i 6 j with some constant c. We have calculated the bound as given in Eq. (3.15) for various values of  as shown in the Table 2. Fig. 3(a) shows the maximum error verses various values of  and N a using central and upwind schemes. We observe that results corresponding to central schemes are better than upwind scheme, which is of first order, using the proposed adaptive mesh. Therefore, we can claim that the proposed method is definitely better than first order but lesser than second order which is the case for uniform meshes for  ¼ 1. Reduction in order of convergence is the price one has to pay for working with nonuniform adaptive meshes. Fig. 3(b) illustrates the distribution of mesh points in the boundary layer region corresponding to different schemes. Test case 2 [9]: We consider another 1D linear convection–diffusion problem with non zero source term as

u00 þ u0 þ u ¼ cosðpxÞ; x 2 ð0; 1Þ; 0 <   1;

ð3:16Þ

with boundary conditions

uð0Þ ¼ 0 and uð1Þ ¼ 0;

ð3:17Þ

The exact solution is

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

6

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx

Table 2 Maximum errors, number of adaptive mesh points N a and bounds (3.15) for u for the test case 1.



Max. error

Na

Bound (3.15) central

2

10

0.0868

20

2.8350

9.0088

104

0.0531

29

6.4638

12.0296

105

0.0830

30

3.1460

9.2625

107

0.0505

45

6.8619

12.3680

108

0.0793

42

3.4645

9.5234

u ¼ a cosðpxÞ þ b sinðpxÞ þ A expðk0 xÞ þ B expðk1 ð1  xÞÞ:

Bounds (3.15) upwind

ð3:18Þ 2

where k0 < 0 and k1 > 0 are the real solutions of the characteristic equation k þ k þ 1 ¼ 0, and

a¼

p2 þ 1 ; p2 þ ðp2 þ 1Þ2

A ¼ a

b¼

1 þ expðk1 Þ ; 1  expðk0  k1 Þ

p ; p2 þ ðp2 þ 1Þ2 B¼a

1 þ expðk0 Þ : 1  expðk0  k1 Þ

Table 3 shows the maximum error computed for the corresponding adaptive meshes. Test case 3 [9]: We consider another 1D linear convection–diffusion problem as

u00 þ u0 þ ð1 þ Þu ¼ 0; x 2 ð0; 1Þ; 0 <   1;

ð3:19Þ

with boundary conditions

uð0Þ ¼ 1 þ expðð1 þ Þ=Þ and uð1Þ ¼ 1 þ expð1Þ;

ð3:20Þ

The exact solution is

u ¼ exp

  ð1 þ Þðx  1Þ



þ expðxÞ:

ð3:21Þ

In this problem the coefficient of reaction term and left boundary condition are the function of perturbation parameter . Table 4 shows the maximum error computed for the corresponding adaptive meshes. We have calculated the bound as given in Eq. (3.15) for various values of  as shown in the Table 4, that also varifies the uniform convergence of the proposed method. Fig. 4(a) shows the maximum error verses various values of  and N a using central and upwind schemes. We again, observe that the results corresponding to central schemes are better than upwind scheme using the proposed adaptive mesh. Fig. 4(b) illustrates the distribution of mesh points in the boundary layer region corresponding to different schemes. Test case 4: [10] We consider 1D linear convection–diffusion problem with turning point as

u00 þ 2ð2x  1Þu0 þ 4u ¼ 0; x 2 ð0; 1Þ; 0 <   1;

ð3:22Þ

with boundary conditions

uð0Þ ¼ 1 and uð1Þ ¼ 1;

ð3:23Þ

The exact solution is

u ¼ exp

  2xð1  xÞ :

ð3:24Þ



Table 3 Maximum errors and adaptive mesh points for test case 2, and starting number of uniform mesh points N = 5.



Max. error

Na

2

10

0.0308

09

103

0.0385

12

104

0.0256

16

105

0.0299

19

106

0.0376

22

107

0.0255

30

108

0.0291

31

109

0.0367

34

1010

0.0255

40

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

7

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx Table 4 Maximum errors, number of adaptive mesh points N a and bounds (3.15) for u for the test case 3.



Max. error (CS)

Max. error (US)

Bounds (3.15) (CS)

Bounds (3.15) (US)

4

2

0.0873306

0.2209864

Na 5

12.656

20.196

26

0.0968888

0.2442103

7

10.757

18.159

28

0.0996411

0.2520102

9

10.283

17.606

210

0.1003516

0.2541902

11

10.164

17.466

212

0.1005307

0.2547530

13

10.134

17.431

214

0.1005755

0.2548949

15

10.127

17.423

216

0.1005868

0.2549305

17

10.125

17.420

218

0.1005896

0.2549393

19

10.125

17.420

0.26

1.6

0.24

1.4

central upwind

0.22

1.2

0.2

central scheme

1

upwind scheme

0.18

u(x)

Max. error

exact sol.

0.8

0.16 0.6 0.14 0.4 0.12 0.2

0.1 0.08

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0 0.99

εN

0.992

0.994

0.996

0.998

1

x

a

Fig. 4. Results for the test case 3.

This test case has a turning point at x ¼ 1=2. In this case we have boundary layers at both the boundary points and the turning point does not coincide with the layer region. Such type of problems are hard to solve with the conventional methods. The proposed method gives quite impressive results as shown in the Fig. 5. Table 5 shows the maximum error computed for the corresponding adaptive meshes. 3.2. System of coupled convection diffusion problems We also consider the coupled system of singularly perturbed convection diffusion problem given as

1 u001  a1 ðxÞu01 þ b1 ðxÞu1 þ c1 ðxÞu2 ¼ f 1 ðxÞ; u1 ð0Þ ¼ a1 ; u1 ð1Þ ¼ b1 ;

ð3:25Þ

2 u002  a2 ðxÞu02 þ b2 ðxÞu2 þ c2 ðxÞu1 ¼ f 2 ðxÞ; u2 ð0Þ ¼ a2 ; u2 ð1Þ ¼ b2

ð3:26Þ

with two overlapping layers near x = 0 (assuming a1 ðxÞ; a2 ðxÞ P 0). Such problems can have perturbation parameters 1 ; 2 of different magnitude. These types of problems have been solved widely using upwind methods on Shishkin type meshes and more literature can be found in paper [11]. The corresponding entropy production operator equations can be found by taking inner product of each of the equation of the coupled systems with respect to 2u1 ðxÞ and 2u2 ðxÞ respectively as given 2

 1 ðS001  2ðu01 Þ Þ  a1 S01 þ 2b1 S1 þ 2c1 u1 u2 ¼ 2u1 f 1 ; 

00 2 ðS2



2 2ðu02 Þ Þ



a2 S02

þ 2b2 S2 þ 2c2 u1 u2 ¼ 2u2 f 2 :

ð3:27Þ ð3:28Þ

The entropies can be calculated at the boundaries points in the same way as done before in the linear case. Results by the proposed method are discussed for the following test case. Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

8

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx Computed solution on adaptive mesh

Computed solution on adaptive mesh

1.2

1.2 exact sol.

exact sol.

computed sol.

computed sol.

1

adaptive mesh

1

adaptive mesh

0.8

0.6

0.6

u(x)

u(x)

0.8

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

x

x

Fig. 5. Results for test case 4. Table 5 Maximum errors and adaptive mesh points for test case 4, and starting No. of points N = 3.



Max. error

Na

2

10

0.0825

35

103

0.1293

29

104

0.0530

53

105

0.0830

55

106

0.1248

49

107

0.0505

85

108

0.0793

79

109

0.1197

69

1010

0.0480

109

Test case 5: We consider system of coupled convection–diffusion problem as [11]

 1 u001  u01 þ 2u1  u2 ¼ expðxÞ; u1 ð0Þ ¼ 0; u1 ð1Þ ¼ 0;

ð3:29Þ

 2 u002  2u02 þ 4u2  u1 ¼ cosðxÞ; u2 ð0Þ ¼ 0; u2 ð1Þ ¼ 0;

ð3:30Þ

which shows two overlapping layers near x = 0. The exact solution of this equation is not known. The results for various values of 1 and 2 for adaptive meshes have been shown in the Figs. 6 and 7. In the Fig. 6, we have generated the adaptive mesh by taking maximum of the entropies S1 ; S2 . The results show the required adaptive mesh but some unexpected adaptive mesh points can be seen near boundary point x ¼ 1. So, we have generated the adaptive mesh by taking minimum of the entropies S1 ; S2 and the corresponding results are shown in the Fig. 7. The adaptive mesh seems better by taking choice of minimum as compared to maximum for the entropies S1 ; S2 . 3.3. Nonlinear convection diffusion problems We also consider the nonlinear singularly perturbed boundary value problem [12]

u00 ðxÞ  uðxÞðu0 ðxÞ  1Þ ¼ 0 2 ð0; 1Þ; uðaÞ ¼ a; uðbÞ ¼ b; a 6 x 6 b:

ð3:31Þ

The resulting finite difference equations give a system of non-linear equations. As suggested by [12], for efficient solution of these equations, we use Newton–Raphson iterations to determine the solution of the non-linear system. The initial guess for the solution is the approximate solution given by

uðxÞ ¼ x  x þ w0 tanhðw0 ðx  xÞ=ð2ÞÞ;

ð3:32Þ

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

9

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx Computed solution on adaptive mesh

Computed solution on adaptive mesh

1

1 sol. u

sol. u

1

1

0.9

sol. u

2

0.8

0.8

0.7

0.7

0.6

0.6

u(x)

u(x)

0.9

0.5

2

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

sol. u

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

x

0.5

0.6

0.7

0.8

0.9

1

x

Fig. 6. Results for coupled case with starting uniform mesh for N = 8, using maximum entropy condition.

Computed solution on adaptive mesh

Computed solution on adaptive mesh

1

1 sol. u

sol. u

1

1

0.9

sol. u

2

0.8

0.8

0.7

0.7

0.6

0.6 u(x)

u(x)

0.9

0.5

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

2

0.5

0.4

0

sol. u

0

0

0.1

x

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

Fig. 7. Results for coupled case with starting uniform mesh for N = 8, using minimum entropy condition.

where

x ¼ 0:5ða þ b  a  bÞ;

w0 ¼ 0:5ða  b þ b  aÞ:

This problem shows an interior layer at x ¼ 0:25 for a ¼ 0; b ¼ 1; a ¼ 1 and b ¼ 1:5. Its conservative form can be written as

u00 ðxÞ 

 2 0 u ðxÞ þ uðxÞ ¼ 0 2 ð0; 1Þ; uðaÞ ¼ a; uðbÞ ¼ b; a 6 x 6 b: 2

ð3:33Þ

There are several choices for the entropy production equation. Again, in analogy to the entropy equation in hyperbolic conservation laws, we use the conservative of the entropy equation which can be obtained as

S00 ðxÞ 

 3 0 2u ðxÞ 2 þ 2SðxÞ ¼ 2ðu0 ðxÞÞ : 3

ð3:34Þ

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

10

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx Computed sol. on adaptive mesh for ε = 0.01

Computed sol. on adaptive mesh for ε = 0.001 1.5

1

1

0.5

0.5

u(x)

u(x)

1.5

0

0

−0.5

−0.5

−1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

−1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Fig. 8. Results for nonlinear case with starting uniform mesh with N = 10.

Computed sol. on adaptive mesh for ε = 0.0001

Computed sol. on adaptive mesh for ε = 1e−05 1.5

1

1

0.5

0.5

u(x)

u(x)

1.5

0

0

−0.5

−0.5

−1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

−1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Fig. 9. Results for nonlinear case with starting uniform mesh with N = 10.

The above choice, with a conservative operator on the left, is the one that captures the correct location of the interior layer. Figs. 8 and 9 show the result of the adapted mesh obtained for different values of the perturbation parameter. As can be seen, points are clustered correctly right around the location of the interior layer. This shows the robustness of the current approach as it works on both linear and non-linear problems with equal felicity. 4. Conclusions In this paper we have presented a new type of robust layer resolving method to deal with oscillation produced by central finite difference method, which is a non-monotone  uniform, when applied to convection dominated singularly perturbed problems with the help of entropy-like variable as the adaptation parameter. We have tested our proposed method for Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019

V. Kumar, B. Srinivasan / Applied Mathematical Modelling xxx (2014) xxx–xxx

11

varieties of problems (linear problem with and without turning point, nonlinear and coupled systems of equations) and results seems worth when dealing with SPP. The proposed method can be considered as an initiative for an alternative to B-type or S-type meshes layers resolving meshes. Our hope is to ultimately extend this method to the Navier–Stokes equations. The coupled and nonlinear cases presented above were done with this in mind. We believe this extension is possible since there is an essentially unique entropy variable for the three-dimensional Navier–Stokes equations. However, several essential gaps need to be filled before this may be attempted. Most importantly, we have found the method presented above to be particularly resistant to any usual methods of convergence analysis. We believe that the encouraging results obtained herein warrant a closer look at methods that could give estimates of convergence rates and convergence behavior. This would be an important step in determining the actual potential of the method. We know from the previous experience of using central difference schemes (non monotone method), one need to choose appropriate values of the turning parameters in order to suppress non-physical oscillations in the numerical solutions using piecewise uniform meshes of S-type and choice of the parameter is also ad hoc. In the proposed method this problems has already been taken care automatically. Also the proposed method seems to be  uniform of order not less than the order of the first order upwind schemes. Only few theoretical analysis of the convergence of non monotone methods has been established so far by researchers for piecewise uniform meshes of S and B-type but is still in its infancy for a completely arbitrary (i.e. non pre-defined) meshes for SPP as discussed in this paper. Finally, it must be mentioned that, in our experience, the above method is effective only for problems with non-zero convection coefficient. There exist BVPs such as the Bratu problem [13] which do not fall under this category and require other approaches. Acknowledgments We would like to thank the reviewers for their valuable and constructive comments and suggestions. References [1] N.S. Bakhvalov, Towards optimization of methods for solving boundary value problems in the presence of boundary layers, Zh. Vychisl. Mater. Fiz. 9 (1969) 841–859 (In Russian). [2] E.C. Gartland, Graded-mesh difference schemes for singularly perturbed two-point boundary value problems, Math. Comput. 51 (184) (1988) 631657. [3] J.J.H. Miller, E. ORiordan, I.G. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [4] Torsten Linb, Layer-adapted meshes for convection diffusion problems, Comput. methods Appl. Mech. Eng. 192 (2003) 1061–1105. [5] S.V. Patankar, Numerical Heat Transfer and Fluid Flow Hemisphere, Washington, DC, 1980. [6] R.B. Kellogg, A. Tsan, Analysis of some difference approximations for a singular perturbation problem without turning points, Math. Comput. 32 (1978) 1025–1039. [7] R.J. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002. [8] P.A. Farrell, A. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin Robust, Computational Techniques for Boundary Layers, Chapman and Hall, 2000. [9] Mohan K. Kadalbajoo, Vikas Gupta, Numerical solution of singularly perturbed convection–diffusion problem using parameter uniform B-spline collocation method, J. Math. Anal. Appl. 355 (2009) 439452. [10] P. Lagerstrom, Matched Asymptotic Expansions, Springer, New York, 1988. [11] Torsten Linb, Analysis of an upwind finite-difference scheme for a system of coupled singularly perturbed convection–diffusion equations, Computing 79 (2007) 2332. [12] R.J. Leveque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007. [13] A. Mohsen, A simple solution of the Bratu problem, Comput. Math. Appl. 67 (2014) 2633.

Please cite this article in press as: V. Kumar, B. Srinivasan, An adaptive mesh strategy for singularly perturbed convection diffusion problems, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.10.019