• No results found

Numerical technique for solving convective-reaction-diffusion equation

N/A
N/A
Protected

Academic year: 2023

Share "Numerical technique for solving convective-reaction-diffusion equation"

Copied!
13
0
0

Loading.... (view fulltext now)

Full text

(1)

Pergamon Mat& Compul. Mo~e~~~n~ Vol. 22, No. 9, pp. 113-125, 1995 Copyright@1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved

0895-7177/95 $9.50 + 0.00 089%7177(95)00172-7

Numerical Technique for Solving Convective-Reaction-Diffusion Equation

P.

C.

JAIN*, R. SHANKAR AND T.

V.

SINGH Department of Mathematics, Indian Institute of Technology, Delhi

Hauz Khas, New Delhi - 110 016, India

Abstract-By using spotting-up technique and cubic splines, a numerical algorithm of second order accuracy is developed to solve the nonlinear equation

ZL~ = Re-r %!z i-

If(U + h(4,

with prescribed initial and boundary conditions. The validity of the difference scheme is tested by applying it to the nonlinear Burgers’ equation. The numerical results obtained for Burgers’ equation at low as well as for high Reynolds numbers are found to be in good agreement with the exact solutions. Then, the proposed scheme is used for solving the convective-reaction-diffusion equation

with initial and boundary conditions. Graphs have been drawn for the numerical solutions at low as well as for high Reynolds numbers by taking the values of parameters in 0 < CI 5 10, 0 < p < 2, and 0 5 E 5 2. The numerical solutions at low values of Re are in steady state for some sets of values of parameters. For intermediate values of Re, the solutions get affected considerably by convective and/or reactive forces, and for high values of Re the solutions blow up. The computed results are compared with the available results for some particular cases.

Keywords-Convective-reaction-diffusion equation, Splitting-up technique, Cubic spiines.

1. INTRODUCTION

The one-dimensional convective-diffusion equation

Ut + vu, = TUXX,

o<x<x,

t20,

(1)

where subscripts denote the partial derivatives, ~(2, t) is a scalar variable which is convected in the z-direction with a constant velocity v > 0 and is spread with constant diffusivity T > 0, with proper initial and boundary conditions, governs the phenomena which occur in many physical situations, for example, thermal pollution in a river system [I], leaching of salts in soil 121, flow in porous media [3], dispersion of dissolved salts in flowing ground water [4] and adsorption of chemicals into beds [5], t e c. A little progress has been made by Davis [6] in order to solve equation (1) analytically with arbitrary initial and boundary conditions except for a few special cases.

In equation (l), if v is replaced by the variable U, it becomes the well-known Burgers’ equation which is a standard model for discussing the phenomena of shock waves and turbulence. It serves

*Present Address: Department of Mathematics, Arts Faculty Extension Building, University of Delhi, Delhi 110 007.

Typeset by A,++S-T&X 113

(2)

as a good test example for developing numerical techniques for the solutions of nonlinear equa- tions ss its exact solution can be obtained by means of HopECole tr~sformations [7]. Further, Levine et ab. f8] have analyzed the case 1 < p 5 m and proved some general results for station~y solutions for any p, m > 1 for the convection-reaction-di~sion equation

ut = Re-iu,, + EIPlz + up. (2)

They [9] studied equation (2) to analyze the steady state problem and have proved the results about finite time blow up. In spite of this progress, most of the nonlinear partial differential equations involving diffusion, convection and/or reaction are intractable analytically. Therefore, it is necessary to use numerical methods to solve these nonlinear equations. In this direction, we [lo] have presented a numerical scheme for the solution of Burgers’ equation with a semilinear boundary condition based on splitting-up and cubic splines. In the present paper, we propose in Section 2 an algorithm for solving the following general nonlinear convective-reaction-diffusion equation:

ut = Re-‘u,, +

mdlz + s(u),

(3)

where Re is the Reynolds number, with appropriate initial and boundary conditions. In order to examine the accuracy and validity of the scheme, in Section 3 we apply it to the Burgers’

equation and compare the numerics results with the exact solutions. Then, in Section 4, we apply the scheme to solve equation (3) and discuss the behaviour of the numerical solutions with f(u) = su2/2 and g(u) = oup.

2. DIFFERENCE SCHEME

Consider the following nonlinear convective-reaction-diffusion equation:

Ut =

Re-luzz +

Wlz

+ g(u). (4

We split equation (4) as

1 -% = fm

3 (5)

?ut = Re-‘umrr 1 (8)

;u, =

g(u).

Now, appro~mating the time derivative by forward difference and the space derivative by the first order derivative of the cubic spline function S=(z) interpolating fp (i = 1,2,. . . , I), equation (5) can be writter) as

elm;+‘/3 + (1 - I+$ = f @;+r/3 - $) ) where ul x U(Xi, tn), xi = 50 + ih, tn = nk, (8)

rnp denotes the value of the first derivative of S%(z) at xi, x0 is the initial point, h is the space mesh size, k is the increment in time and @I E [O, 11, a cubic spline parameter. F’rom the condition of ~ontinui~ of the second derivative of S&(z) at zi, we have [ll]

my+I + 4rny + mr_l = ; (frj.1 - fL).

On elimination of m from equations (8) and (9), we get

(9)

(3)

115

For computational work, we remove the nonlinearity in equation (10) by using Taylor series expansion as

fty+1/3

- f; = ;(fi,: + O(k2) = A; (I$+~‘~ - u;) , where A = fu. (11)

Substituting (11) in equation (lo), we get

In the same manner, we approximate equation (6), by using forward difference for the time derivative and the space derivative by the second order derivative of the cubic spline function S:(X) interpolating up (i = 1,2,. . . , I) to get

(13) where 02 E [0, l] and MS? denotes the second derivative of S;(z) at xi. The condition of continuity for the first derivative of S;(x) at xi is given by [ll]

M;, + 4MF + MF_, = $T2ur. (14)

Eliminating M from equations (13) and (14), we get

Further, the approximation of equation (7) gives nn+l

2 - u;+~‘~ = e3kg;+l + (1 - 03)kg”+2’3, where e3 E [0, 11. We remove the nonlinearity in equation (16) by taking

hg;+’ - hgn+2/3 = B;+2/3 (uF+’ _ u:+2/3) + O(k2), B;+2/3 = (9U);+2/3.

where

Thus, equation (16) takes the form

1 - k83B”+2’3 ) IL;+’ = (1 - ke3BnC213) + kgn+2’3.

(15)

06)

(17)

(18)

Equations (12), (15) and (18) constitute a multistep implicit finite difference scheme for solving equation (4).

Local truncation errors in equations (12), (15) and (18) are given by

CUtt)i n+1/3 + k3 (A - 2) (Uttt)n+1/3 + g(utt):+1’3, and (20)

k2(~-~)(utt):‘2’3+k3(&~)((Ltlt)~+2’3,

(21)

(4)

respectively, by neglecting the terms of higher order. Therefore, equations (19), (20) and (21) give truncation errors of order 0(k2 + k2h2 + kh4), 0(k2 + kh2) and 0(k2), respectively, for 13i,82,Bs E [O,l]; whereas for 81 = 02 = 0s = l/2, equations (19), (20), and (21) provide the truncation errors of the order of O(kh4 + k3), 0(k3 + kh2) and 0(k3), respectively. Also the truncation error in calculating ~a+’ from ua for equation (4) is of order of 0(k3 + kh2). Therefore, the algorithm given above is second order accurate in space as well as in time direction.

For the stability of the algorithm, we apply the method to the linear form of equation (4) and then use the von Neumann criterion of stability and get the necessary condition of stability as

3.

NUMERICAL SOLUTION OF BURGERS’ EQUATION

As pointed out in the introduction, we are here concerned with the solution of Burgers’ equation

by using the scheme developed in Section 2. It may be noted that the equation (23) is a particular case of equation (3) for f(u) = -u2/2 and g(u) = 0. Further, we take the following initial and boundary conditions as

0.5 Ix I 1.5, u(x,O) =

& [x +

tan

(z)] ,

u(O.5,t) = (r +lRe) ~ [0.5+tan(4(tyRe))], t20, u(i.5,r) = (r +lRe) ~ [1.5+tan(4(p+Re))], t20.

The exact solution of this problem is given by

u(x,t) = -

(t

+fRe)

[x+tan(2(tYZe))] ’

(24)

(25) For computational work, we have taken h = 0.05 and r = 0.4. It is found that the computed results agree with the exact solution at low Re within an error of order 10-2, while at intermediate and high Reynolds numbers the orders of error are 10e5 and 10-12, respectively. For the sake of testing, we have also taken h = 0.1 and h = 0.025 and found that the algorithm developed in the preceding section works satisfactorily for the solution of Burgers’ equation for different values of mesh parameters in the range 1 5 Re 2 105.

We have also considered the problem of a uniformly propagating shock [7], which is given as Ut + uu, = -uzz, 1

Re x E (-co,co), with (26)

u(x, 0) = {

ui = 0.5, x < 0, Ug = 1.5, x > 0.

The exact solution of this problem is given by

u+(u2 -w)

u = [l + exp (-Rey(x - it))] ’ (27)

where U = (ui + u2)/2.

In the computational work that follows, we have restricted the domain; that is, x E (-0.5,0.5), t > 0 by taking the boundary conditions as u(-0.5,O) = 0.5 and u(O.5,0) = 1.5 to solve this problem, and the numerical results are plotted. Figure 1 shows the computed results for Re = 10000 by taking r = 4000 and h = l.O/lOOOO.O at t = 0.1. It is clear that the proposed scheme gives the numerical results which are in good agreement with the exact solution.

After getting satisfactory results for Burgers’ equation, we now apply this scheme to solve the convective-reaction-diffusion equation in the following section.

(5)

1.60 1.5

1.20

t . u

OAO-

I-

0.5-- 0.40-

,Numerical

~ Exact

0.2OL

-0.50 - 0.50 -0.io O.lb OS0 0.40

_-x-M -

Figure 1. Numerical solution of Burger’s equation for Re = lo4 at t = 0.1.

4. NUMERICAL SOLUTION OF

CONVECTIVE-REACTION-DIFFUSION EQUATION

By using the scheme developed in Section 2, we have solved equation (3) for f(u) = EU~/~ and g(u) = cruP with the boundary conditions

21(0,t) = 21(1,t) = 0, t > 0, and following two types of initial conditions

(i) u(5,O) = sin(zr), O<zll, (28)

(ii) u(z,O) = 1

102, 5 E [O.O,O.l),

-loz + 2, 2 E [0.1,0.2), (29)

0. z E [0.2,1.0].

Since the basic equation contains Re, E, o and p as parameters, we have taken some representative values of these parameters in the range 1 < Re < 10000, 0 5 E 5 2, 0 < a 5 10 and 0 < p 5 2.

Now we shall discuss the different cases of equation (3) one by one by taking t 5 5.

CASE (i). For E = 0, (Y = 0, equation (3) reduces to

it = Re-‘u,,, (30)

which is known as diffusion equation. In this case urnax; that is, maximum value of IL decreases with time for low and intermediate values of Re; that is, Re < 1000, whereas for high value of Re;

that is, 10000, the solutions are steady state within the considered time. The numerical results do not change with initial conditions. It may also be noted here that equation (30) has the exact solution:

2 u = sin(nz) exp -$-J ,

( > (31)

with the first initial condition and our numerical solutions in agreement with the exact solution.

(6)

CASE (ii). For E # 0 and o = 0, equation (3) becomes convection-diffusion equation

ut = ReV1uz, + EUU,. (32)

In this case for low Re I 100, the value of urnax decreases for all values of E with a tilt towards left (Figure 2). The effect of E is seen on the tilt of the solution profile. It is observed that the tilt increase with E. For higher values of Re(= lOOO), the numerical results depend on the initial condition. It is found that with the first initial condition, the solution profile gets decreasing, whereas for second initial condition the solution tends to infinity for all E # 0. For a high value of Re(= lOOOO>, the sofution tends to infinity irrespective of initial conditions.

0 0.1 0.2 0.3 0.4 05 0.6 0.7 0.8 a9 1.0

X

(1) Time-0.62s (2) Time- 1.25 ES) Time- I.675 (4) Time-P.5 kjl Time- 3.125 (6) Time- 3.75 (71 Time-4.375 (81 Time- 5.0

Figure 2. Behaviour of the solution for a = 0, E = 0.5 and Re = 100 with first initial condition.

CASE (iii). For E = 0, a: # 0, equation (3) is the reaction-diffusion equation

ZQ = Re‘-‘uzr + aup. (33)

In this case, we consider the following sets of values of p:

(4 0 <

P <

1, (b)

P =

1,

(c) 1 < p 5 2.

(a) For Re = 1, it is clear from Figure 3 that the solution becomes steady state. In this case, we observed an effect of o on the solution; that is, for cw < 6, the steady state is the result of decreasing trend in the solution, and for Q: > 6, the steady state is the result of increasing trend in the profiles, whereas for a! = 6, the solution gets steady state at the initial value itself. In this case, if we take higher values of Re(= lo), the solution becomes dependent on the initial condition. It is found that the numerical solution with the first initial condition decreases with time (Figure 4) for all CY, whereas with the second initial condition, solution decreases with time for (Y 5 5 (Figure 5) and increases with time for cy > 5 (Figure 6). For more high values of Re(> loo), the solution tends to infinity for both the initial conditions.

(7)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1.0 x

Figure 3. Steady state profile for cy = 6, E = 0, p = 0.5 and Re = 1 with first initial condition.

X

[I) Time-0.623 t2) Time- 1.25 (3) Time-I.873 (4) Time-Z.5 (5) Time- 3.123 (6) Time-3.75 (71 Time- 4.373(S) Time-f.0

Figure 4. Behaviour of solution for cx = 4, E = 0, p = 0.5 and lb = 10 with first initial condition.

(b) In this case, the value of urnax depends on the value of c11 for Re = 1. It is observed that the value of urnax decreases with time for a 5 5 (Figure 7) symmetrically with z, whereas the value of %nax increases with time for cy > 5 (Figure 8). For a higher value Re( 2 lo), the solution tends to infinity irrespective of the initial conditions.

(c) In this case, the solution depends on the value of cy and initial conditions for Re = 1. It is observed that with the first initial condition and for ct: 5 3, the solution decreases with time and for a > 3, the solution increases with time towards infinity. By using the second initial condition,

(8)

(I) Time - 0.62s (21 Time- t.25 13) Time- I.875 14) Time-2.5 I31 Time- 3.itS (61 Time-3.?S 171 Time- 4.375 (8) Time- 5.0 Figure 5. Behaviour of solution for a = 4, E = 0, p = 0.5 and Re = 10 with second initial condition.

I(

t

E 3

4

21

C 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 09 tl

tl) Timr-O-625 (2) Time- 1.25 (3) T.bne- I. 1373 14) time-2.3

Figure 6. Behaviour of the solution for a: = 8, E = 0, p = 0.5 and Re = 10 with second initial condition.

we find that the solution is decreasing to zero for all values of cy. For this case, Levine [l2] gave a theorem, which is as follows.

THEOREM. Let P,(N) = 1 + 2/N (N = 1 in present case) and 1 < p 5 P,(N), then the only nonnegative global (in time) solution of equation (33) is u = 0.

For details about this work, one may see [12]. Our numerical solution for p > 1 is in conformity with 1121 as we either get the blowing up or zero solution. For Re = 10, with first initial condition

(9)

a==

tf) Time- 0.625 (2) Time- 1.25 13) Tim?- 1.875 141 Time-2.5 (5) Time3.125

Figure 7. Behaviour of solution for (2 = 4, E = 0, p = 1 and Re = 1 with first initial condition.

121

0 O*I 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.0 1.0 X

(1) Tim-O.625 12) T’imc-1.25 (3) Tim- 1.67s t4)Tim-2.5

Figure 8. Behaviour of solution for QI = 6, E = 0, p = 1 and Re = I with first initial condition.

the solution blows up for all values of CY., whereas with the second initial condition we find that the solution decreases with time for a! 2 3 (Figure 9), and the solution blows up for CY > 3. For higher values of Re(> lCKI), the solution tends to infinity irrespective of initial conditions. This phenomena can be described as follows: for higher values of Re, the diffusion coefficient becomes small and the source term auP becomes dominant and causes the solution to blow up. Also, it is observed that, in the absence of convective forces, the solution is always symmetrical with respect to 2.

(10)

(1) TiaW- 0.625 42) Time- 1.23 (31 Time- 1.675 (4) Time-2.3 tS1 lime- 3.123 (6) Time- 3.73 171 Time- 4.37s (8) Time- 3.0 Figure 9. Behaviour of soiution for a = 2, E = 0, p = 2 and Re = 10 with second initial condition

(

0 0.1 0.2 0.3 Q-4 0.1 0.6 0.7 0.8 0.9 I.0

X

Figure 10. Steady state profile for 1y = 6, E = 1, p = 0.5 and Re = 10 with second initial condition.

CASE (iv). Finally, we consider the general case, that is, E # 0 and a # 0 in equation (3). In this case also, we consider three cases depending on the values of p.

(a) 0 <p < 1. In this case, u max gets steady state for low values of Re(l 100). For Re = 1, 10, the solution profile has a steep tilt towards left (Figure lo), whereas for Re = 100, the solution gets steady state oscillating profile (Figure 11). This seems to be due to the effect of the nonlinear term aup. As we go to higher values of Re, the results become unstable; that is, they tend to infinity with large oscillations. The value of u max in steady state increases with a! and decreases with E. Also, the tilt increases with cx as well as with E.

(11)

SC

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1.0

Figure 11. Solution profile (steady state)?or (Y = 4, E = 1, p = 0.5 and Re = 100 with first initial condition.

!a

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1.0

X

Figure 12. Steady state profile for cy = 6, E = 1.0, p = 1 and Re = 10 with first initial condition.

(b) p = 1. In this case for Re = 1, u tends to zero solution for (I: 5 5 symmetrically with 2, whereas for CY > 5, u attains the steady state having a tilt towards left which increases with Q as well as with E. For Re = 10, the solution reaches the steady state profile having a tilt towards left (Figure 12) for all considered values of cx and E. For higher values of Re(= loo), the solution depends on the value of QI. It is observed that the solution attains steady state for (Y I 5 (Figure 13), whereas for (Y > 5 the solution becomes infinite. The solutions for Re > 100 are found to be approaching towards infinity.

(c) p > 1. In this case for Re = 1, the solution is approaching to zero value for o 5 3, whereas for (I > 3, the solution is infinite for all E. We observed the effect of E for Re = 10 and Re = 100.

Figure 14 shows the effects of E on the value of u,,, for different values of Q: and Re = 10. It is

(12)

6-

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1.0

Figure 13. Steady state profile for a = 4, E = 1, p = 1 and FLe = 100 with first initial condition.

16 -

14-

IZ-

IO-

a-

6-

4-

2-

O-

0 0.2 0.4 0.6 0.6 1.0 1.2 1.4 1.6 1.6 2.0

e-

(1) U=I.O (2) a = 2.0 (3)a - 3.0

Figure 14. Effect of E on urnaX (Fk = 10, p = 2, and first initial condition).

(13)

clear that for LY = 1, the solution is infinite for E < 1, whereas for E > 1 the solution is tending to zero. Similar behaviour of the values of u,,,, is found for Re = 100. Thus, it is concluded here that the presence of strong convective force can balance the system which involves diffusion and reaction.

5. CONCLUSION

From the above study, we conclude that for low values of Re, the solution of the convective- reaction-diffusion equation is found to be in steady state for some sets of values of Q and p, whereas for intermediate values of Re, the solution gets affected dominantly by convective and/or reactive forces. For high values of Re, the solution becomes unstable. Also, while comparing our numerical results for various particular cases with the results available in the literature, we find that our numerical results are in good agreement.

REFERENCES

1. M.H. Chaudhry, D.E. Cuss and J.E. Edinger, Modelling of unsteady flow water temperatures, J. Hydraulic Engg. ASCE 109, 657-668 (1983).

2. E. Bresler, Simultaneous transport of solutes and water under transient unsaturated flow conditions, Water Resource Res. 9, 975-986 (1973).

3. N. Kumar, Unsteady flow against dispersion in finite porous media, J. Hydrology 63, 345-358 (1983).

4. V. Guvansson and R.E. Volker, Numerical solutions for solute transport in unconfined quaifers. Int. J.

Numer. Methods Fluids 3, 103-123 (1983).

5. L. Lapidus and N.R. Amundsom, Mathematics of adsorption in beds. VI-The effect of longitudinal diffusion in ion exchange and chromatographic columns, J. Phys. Chem. 56, 984-988 (1952).

6 G.B. Davis, A Laplace transform technique for the analytic solution of a diffusion convection equation over a finite domain, Appl. Math. ModelZing 9, 69-71 (1985).

7. P.L. Sachdev, Nonlinear D#usive Waves, Cambridge University Press, Cambridge, (1987).

8. H.A. Levine, T.F. Chen and P.E. Sacks, Analysis of a convective-reaction-diffusion equation, Nonlinear Analysis TMA 12, 1349-1370 (1988).

9. H.A. Levine, L.E. Payne, P.E. Sacks and B. Straughan, Analysis of a convective-reaction-diffusion equa- tion II, SIAM J. Math. Anal. 20 (l), 133-147 (1989).

10. P.C. Jain, R. Shankar and T.V. Singh, Cubic spline technique for solution of Burgers’ equation with a semilinear boundary condition, Comm. Appl. Num. Methods 8, 235-242 (1992).

11. J.H. Ahlberg, E.N. Nilson and J.L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York, (1967).

12 H.A. Levine, The role of critical exponents in blow up theorems, SIAM Review 32 (2), 262-288 (1990).

References

Related documents

Attempts have been made to explore the exact periodic and solitary wave solutions of nonlinear reaction diffusion (RD) equation involving cubic–quintic nonlinearity along with

In this work, we have developed an accurate finite difference solver ‘SIBIDIF’ for numerical simulation of two-dimensional drift diffusion equation (DDE) describing two-

Some of them are: the (G /G)- expansion method [1–4], the simplest equation method [5–7], the solitary wave ansatz method [8–10], the first integral method [11–14], the

In this paper, exact solutions of Benjamin–Bona–Mahony–Peregrine equation are obtained with power-law and dual power-law nonlinearities.. The Lie group analysis as well as the

In addition, the numerical analysis of the Schrödinger–KdV equation with power-law nonlinearity is studied using the variational iteration method (VIM) and homo- topy

As a result, more types of exact solutions to the generalized KdV equation with generalized evolution are obtained, which include more general single-hump solitons, multihump

In this paper some exact solutions including soliton solutions for the KdV equation with dual power law nonlinearity and the K ( m , n ) equation with generalized evolution are

In the paper, we proposed a new irrational trial equation method and used it to obtain some exact travelling wave solutions to the Burgers–KdV equation and the dissipative