On a difference approximation scheme for the numerical solution of parabolic type equations

ON A DIFFERENCE APPROXIMATION SCHEME FOR TIIENUMERICAL SOLUTION OF PARABOLIC TYPE EQUATIONS* I. V. PETUKHOV Moscow (Received

27

1965)

October

THEboundary conditions for a number of two-dimensional problems, described by parabolic type equations, are formulated as follows: values of the required function u are given on the characteristic x = 0 and certain boundary conditions are given on the boundaries of the region O<,Y &.Yes x>,O. The solution of such problems can be reduced by some difference approximation of the derivative &/& to a step by step solution of a boundary value problem for an ordinary differential equation (with argument y). In the present paper a new autonomous difference scheme of the second order of accuracy is described for approximating the derivative &/&. This scheme is applied, without loss either of accuracy or degree of stability of the calculation, in the case of the degeneracy (complete and partial) of the parabolic type equation in general (with respect to ,y). The partial degeneracy of the equations is typical of the problem of a boundary layer and in certain cases leads to a decrease in the effectiveness of the analogous difference scheme used earlier for approximating to the second order of accuracy. 1. We shall consider below a numerical method applicable to the solution of equations of the form a?2 = L[u],

o&y)

20,

L[u] = (mu’)’ + 1iPU+ kiU + ko, Here and in future -~

~ Zh.

u~chisl.

(1)

m(x, Y) > 0.

ti = au / dx, u’ = du / ~3y.

Mat. mat. Fiz. 6, 6, 1019 - 1028, 1966. 103

(2)

I.V.

104

The solution

of equation

04

Petukhov

(1) is considered 17:d

XL,

in the region

O=G?/<&

(3)

with certain boundary conditions on the boundaries of the region ,y = 0 = .ye and with initial conditions on the characteristic x = 0: and N a(O, Y) = ho@). In the case (1) which we are considering, a degeneracy (complete or partial) of an equation of parabolic time is generally permissible, i.e. the vanishing (in the whole or nart of the region (3)) of the coefficient a(x, y), Such a case is characteristic of boundary layer problems where the analogue of the coefficient 0: vanishes on the boundary y = 0 because of adhesion conditions, In addition, with some tr~sformations of the original equation, connected with the standardization of the region of the solution of (3), the coefficient a vanishes on the initial characteristic x = 0 [Il.

FIG. 1. Obviously a more rational method of solving equations of such a form is to reduce (by some difference approximation of the partial derivative a} the partial differential equation in general (with argument .y) to the successive step by step solution of the corresponding boundary value problem by the recursion method. At each stage the required values U=U+ on the characteristic nt = x- + hx are determined on the basis of the calculated on the preceding characteristic values u = U-, previously With such an r = r-, and on the basis of the boundary conditions. approach the method of calculation is divided into two independent elements:

On a difference

approximation

(a) the method of difference

105

scheme

approximation of the derivative

zi,

(b) the method of solving a boundary value problem for an ordinary differential equation. Here we shall confine ourselves to investigating methods of difference approximation of the derivative zi, without considering the method of solving the boundary value problem for an ordinary differential equation. 2. We shall consider three typical schemes of difference approximation of 8. Schemes of this form have been used. for instance, in [l - 31. Scheme

I

-

an approximation scheme of the first

order of accuracy:

x+=x-+h.

C-9

Here and in future rejected terms are enclosed in braces. In the case of degeneracy a = 0 the solution is determined exactly. With a = O(l), however, the solution is determined with an error of order h with values of the argument x = nh of order unity. Scheme

II

- 811 approximation scheme of the second order of accuracy:

(5) Here x0 = X- t h, X+ = ~0 + u^ =

(see Fig. 1);

h

u++ u2

’

u+=

a;-

u-.

(6)

Equation (5) is used for determining i?. The values of ut are then determined from (fi). If a = O(1) the solution is determined with an accuracy of greater order than in Scheme I. The degenerate case, a = 0, however, deserves consideration. We rewrite

(6) in another notation: UR =

xn =

2iLn xn--1 +

Un_-f, h,

and denote by &I the error in the quantity u (defined here as the

(7)

106

I.V.

Petakhou

difference between the approximate and exact solutions ential equation). From (5) it follows that if a = 0

of the differ-

where a,,is some quantity, of single type in n, of order unity. Using (6) twice we obtain 6u, = &z&-2 + h2(a, - an-i) = 6l.b2 + 0 (h3).

(9)

Since a,, is a quantity of single-type in n the accuracy of order h2 will not be lost if the number of steps is increased. Here, however, the strong influence of the initial error on the succeeding behaviour of the solution is characteristic. Repeatedly using (7) and (8) we obtain

bl+h =(--l)k

6u, + h2 i,

(-a,+p*

w

83i

Suppose that for some reason or other the error 6u, is significant. From (10) it is ‘obvious that this error (retaining the absolute value) will be imposed alternately with the sign (+) and (-), giving the solution a characteristic “saw-tooth” form. Thus if a = 0 the influence of the initial error is not damped; in other words the stability of the solution is found at the limit. Scheme

III

- an

apgroximation scheme of the second order of accuracy

a+

3u++ u--4&’

h2 . .

---

2h

-?f”{

l-J

=

L+[u+].

Here, exactly as in Scheme I, the solution for a = 0 is determined exactly. If a = O(1) the solution is determined with an error of the same order as in Scheme II. As distinct from the two preceding schemes, however, Scheme III is not autonomous: to solve equation (1) it is necessary to know u- and uo directly on two preceding characteristics. Therefore Scheme III is not applicable in the two cases:

On a difference

appr~x~~ution

(a) with departure from the original

107

scheme

characteristic;

(b) with departure from the characteristic undergoes a discontinuity.

on which the derivative

ci

3. The idea of the proposed method consists of the following. We are given some value ii = ii0 and apply Scheme I, taking into account the term in the braces, successively to the intervals cc-, noI and LO, x+1 (Fig. I) :

ao

f\ uiO; u-+-;ii0 ) = LO [u*@],

CZ+

i

z&i+- l&o h

h +

2

‘. 1 ‘lo J=

(12’) (12”)

L+ [Uif].

If a = 0 the values ~0 = aI0 and at = ul+ are calculated exactly. Let a = 0(1) so that we can put a s 1. Subtracting (12’) from (12’“) we obtain & z=z

u- + ui+ - 2uio h2

L+ [z&q=-

Lquq

(13)

h

We shall denote by Silk = Uk - u the difference between the approximate and the exact solutions of the differential - difference system (12). Let 6iio = O(1). Then obviously that

6uio = G(P),

L[6u] = O(6U).

6&

= O(hZ). We assume (14)

Then it follows from (13) that Sii1 = O(h). Repeatedly solving equations (12). in which ii1 is substituted for iio, we successively determine the 6u2* = ~(~3~. values of u20 and u2+ with error of order h3: 6~~0 = O&3), Thus with condition (14) a twofold solution of system (12) enables us to determine the solution with error of order h3. Condition (14) is satisfied, for instance, for the case of the ordinary first-order differential equation (with respect to x) au+Bu+y=O,

if the coefficients

a and @ are of order unity.

4. There is a definite analogy between the solution of a boundary value problem with initial data on the characteristic x = i? for

108

I.V.

Petukhou

equation (1). (2) and the solution of the Cauchy problem with a given initial value for u at the point x = 0 for equation (15). Let us consider, for instance, the simplest form of equation (l), (2) uti = u”,

o = const > 0,

(16)

with boundary conditions u = U,(x) ~pproximating equations

tl

if

y=o,

ZJ= ue(s)

if

y = &.

($7)

‘: we replace equation f 16) by the system of ordinary 2

o&J-f--urn= 12

2 %-z-H+ urn-i 2 3

22

Ym+i =ym+L

System (18) can be solved by iterations, taking the values of the right-hand sides from the preceding approximation. If a > 0 such a process obviously converges automatically. From (18) it is obvious that the above analogy must be made with the condition @iO

(a > (0,

UQI

where the coefficient p can be taken as large as we please. The last condition is most essential since in this case condition (14) is not satisfied. Let us consider the iterative to equation (15):

process for solving system (12) applied

Uk+- uk” h _* h +guk-i

It is easy to obtain the following k-th approximation:

*iik=

--$&+- $.

=

expression for the error 6iik in the

(6+/a+- ‘/$OhO) h

(i’ +

pOh/aO)

(I

+

p+h,a+)*Uk-i-

From (21) it follows directly that, with condition tive process always converges, where

(l§),

(21) the itera-

On a difference

isii k

I

0(%.&-i)

approximation

if

{=

O(1)

109

scheme

and if

Blower convergence takes place with $/a = @l/h). We shall evaluate the speed of convergence for this case. Assuming that $+/a+ =@*/a0 t O(l), we put $+/a+ = @*/a* = P/a : (22) The maximumof the ratio Giik/Gi~_r is reached for @/a equal to l/8:

5. For the system (16). (17) it is elementary basis for the method under denote by 6Uk = ilk - u the difference solutions of the system (12)‘ Writing we obtain

since any approximationUk satisfies

=

1 /h and is

not difficult to give a strict and consideration. As before, we shall between the approximate and exact out system (12) for equation (16)

the boundary conditions

(17) we

have 6Uk-V) = &Zk-+(&=)= 0.

dui?(O) = 6Uk0(&?)= 0,

(25)

As is well known the solution ZI of the system v”

-

_au=-f,

h can be expanded in the Fourier series

where the normslfsed characteristic

cpl.“+(L+Yl, =o, The coefficients

v(O)=

v@,)=O

functions q,.(y) satisfy

cpm= fP&?) =

the system

0.

vl. can be expressed in terms of the corresponding

I.V.

110

coefficients

f, of the function

In the case considered

Using these

results

for system

7

=

a

(&6k+)r =---h where (&a), are the Fourier fact that

we find

and the eigenvalues

f(y)

= const.)

(a/h

(6UkO)

8iik =

Petukhov

we have

(24) we obtain

-

-

ah (6&-i)r 2 a,

(6uk’)r

A+

coefficients

6nk+ - 26uk” h2

A,:

,

ah (I%k-1) ,. 2

h-

of the function

6~. Using the

h,=;+(E)2,

that (6iik-i)

Formulae (26) anki (22) are completely also analogous ]

analogous.

P-3)

r.

The result

(23) is

(6iik)rl < ‘/aI (6Gk--i)rl.

Thus the successive approximations 6Uk for system (lFj), (17) converge no more slowly than l/Sk. The first coefficients of the Fourier series converge as O(hk). 6. The method considered was approved in calculations of currents of gas in a boundary layer. These currents were described by a system of related non-linear equations of the form (1), (2). where the coefficients of the equations depended on the required functions. We note that for systems of such a form both equations (12) are solved initially for the first, then for the second equation of the system etc. Thus each

On a difference

approxi~&~ion

schene

111

iterationis carriedout. Any other order is not permissiblein the general case: we cannot first solve equation (12') alone for all the equationsof the system and then equation (12"). Some resultsare given below of calculationsfor a bounded layer on

a sphericalblunt cone dependingon the argumentx, expressedIn degrees:

I

80’

I

100

I

150

I

200 X”

FIG. 2. =s/R, x*=x- 180 n, where R is the radiusof the blunt part, and s the length of the generator,measuredfrom the forwardcriticalpoint.

x

112

1 .V.

Petukhov

This case is one of the most difficult since the derivatives with respect to x for the solution on the conical part of the body tends to infinity with the approximation from the right to the junction point of the sphere and cone [41. This circumstance leads to some difficulties in the use of Scheme II. Thus to calculate the characteristic x = x1 t Dx, following directly behind the point of junction, the solution contains a relatively large error, due to the presence of very large derivatives with respect to X. In connection with (10) this error, which is slowly damped, is imposed on the consequent solution. To eliminate this influence the first steps POXafter the junction point x1 should be very small, as was the case in [41. In Fig. 2 we quote, relative to constant values, the friction stress -r2 and the specific heat flux q, on the surface of the blunt cone:

Here the notation of [41 is used. The results of the calculations spond to a cone of semivertical angle o = loo, the temperature of surface T, = 2000“ and the number q,,, = 6. The values of -r, and qs given on a part of the surface consisting only of a small part of spherical surface to the left of the function point x1 = 80’.

correthe are the

The calculations were carried out by two methods - according to Scheme II and the proposed scheme. In this case and the other the method of solution of the corresponding boundary value problem on the characteristic x = const., put forward in [I], was used. In both cases the calculation was carried out with the same constant step DxP = 1.25O and with the same distribution of the nodes on each characteristic: q = 0.00. 0.25. 0.50, 0.75, 1.00, 1.50, 2.60, 2.50, 3,00, 3.50, 4.00, 4.50, 5.00, 5.50, 6.00, 7.00, 8.00, 9.00, 10.00. Up to the junction point, with n xl, the results of the calculation using the proposed scheme are indicated by the smooth curve, and using Scheme II by the sawtooth line. In addition the calculation was carried out by Scheme II in correspondence with the law of change of step Dx, used with [41. The first step Dxl after the junction point was taken equal to Dx, = 1.25/32’. The results of the last calculation in Fig. 2 coincide with the results of the calculation by the proposed scheme.

On a difference

In addition calculationa ent (constant) steps

approximation

were carried out by this scheme for differDx

and for a different x = const.:

113

schene

=

2”.5,

5”,

10”

(given) number of steps II on the characteristics

To the number rz = 13 there corresponds the distribution above. For the two other cases we used if

n = 9:

if

n =

q =o.o.

0.5. 1.0, 2.0, 3.0,

of nodes given

4.0. 5.0, 6.0, 5.0, 10.0;

4: q = 0, 1, 3, 5, 8.

Each case was computed for a different

number of iterations

12= 00, 2, 1. The variant k =co corresponds to the case where the iterative process is carried out completely, i.e. right up to the coincidence with the given small error E (E = 2”’ z 0.47T10-g) of the last two iterations. The actual means of the value k, corresponding to this case, vary with the limits 3 - 6. In the remaining variants exactly k iterations were carried out for each characteristic, excluding the first two and the first after the junction point: for the latter the iterative process was TABLE1 6% I Dx

I k

CC”=40

0.2816 2”.5

2” 1

5”

2” 1

loo

00

2 1

0 0 -0.0001 0 0 3

80

I

( m= 18) 90

/

ii0

0.0739

230

t40

0.0609

0.0514

0.0438

0.0001 -0.0002 -0.0001 0.0001 -0.0004 -0.0001 0 0 0.0001

0 0 0.0001

0 0 0

0 0 0

-0.0006 -0.0062 0.0001 -0.0011 -0.0002

-0.0001 -0.0001 0.0001

0.1901

0.0952

0.0001

0.0001

0.0038

0 0.0003 -0.0025 0 0.0002 -0.0026 0.0004 -0.0002 -0.0029 I

0.0001

-0 ..0005 -0.0002 -0.0003 -0.0005 -0.0001 -0.0014

-0.0001 --0.0001 0

-0.0001 -0.0001 0

-0.0002 -0.0002 -0

-0.0001 -0.0001 -0.0001

114

I.!‘,

-

-

-

-I -I

Dx

k

Petukhov

230

X0=40

0.281f i

-0.0001 --0.0001 -4.0001

0 0 0

1

--O.QOOl --0 . 000 1 0

0 0 0

2” 1

0 0 0 IO004

2”

2O.5

1 5”

2”

10”

0.0238

--Q‘OOOI

-0.0001 0

TABLE 3

DX

I

k

I 1

T 2O 5

X0=40

80

90

0.28t8

O.fWI

0.0952

1

0.0739 I

140

0.0609

j

180

j

0.05I4

230

0.0438

2” 1

--0.0011 --O.OOfI

--0.0060 --0 .0060 --0.0060

0.0028 0.0027 0.0026

““-0.0003 --0.0004 --0.0002 --0.0001 --0.0003 --0 .0004 --0 ,0002 --0 , 0001 --O.OOO2 --0.0003 --o.O+.m --o,Om

2” f

--0.0010 j--0.0011

--0.0059 --O.QO63

0.0027 9.0019 O.OOIi

--0 .0003 --0 .0004 --0.0032 -0.0001 --0.0004 --0.0004 --a, 0002 --0 ”0001 --0 . OQO2 -0.0003 --0.0002 --0,OOOl

2”

--0.0010

--0.0057 --0.0067

--0.0006

--0.0057

o.oOOB 0.0008 0.0007

--0 .oOO5 --a. 0005 --0 <0003 --0 . OOOl --0 .0007 --cl. 0005 --0 .0003 0,0001 --0 . 00 15I --0 .0003 --0 .0002 --0 . 0001

5”

f0”

ii0

1

I

always carried out completely. On each characteristic, as the zero approximation for analogs of the second derivative ii, of (B), the values computed on the preceding characteristic were taken = In this connection we note that the variant k = 1 is not equivalent to Scheme I: it leads to a more exact result (in order).

fn

Table d, - 3 results are given of calculations

in the form of

On a difference

approximation

scheme

115

differences 6-r, between the computed and exact values of 7,. The exact values of T, are given in the upper row of each table. The tables show that the variants k = a, and k = 2 give the same’ accuracy in practice. Consequently it is reasonable to carry out exactly two iterations for each characteristic (with an obvious exception). Trans

Zated

by

H.F.

Cleaves

REFERENCES 1.

PETUKHOV, I.V.

The numerical calculation In Sb. Numerical

in a boundary layer.

of two-dimensional methods

of

solving

flows differ-

ent iat and. integral equations in quadrature formulae (Chislennye metody resheniya differ, I. integr, yr-nii i kvadraturnye formuly). Izd-vo Akad Nauk SSSR, 304 - 325, Moscow, 1964.

2.

SRAILOVSKAYA, I.Yu. and CHUDOV, L.A. The solution of boundary layer In Sb. Computational Methods and equations by a difference method. Programming (Vychisl. metody 1 programmirovanie). Izd-vo Moscow State University, Moscow, 167 - 182, 1962.

3.

GROMOV,V.G. The application the solution of boundary (Series Mekh. Mashinostr.

4.

ANKUDINOV, A.L. The results blunt cones in supersonic 960 - 966, 1965.

of a three-layer difference scheme for layer equations. IZV. A&ad. Nauk SSSR ,9 - 10, 5. 124 - 133, 1965. of calculation for flow. Zh. vj;chisl.

a boundary Mat.

mat.

layer on Fiz. 5, 5,