Solution of nonlinear boundary value problems—VIII. Evaluation of branching points based on shooting method and GPM technique

Solution of nonlinear boundary value problems—VIII. Evaluation of branching points based on shooting method and GPM technique

Chemical Engineering Science, 1974, Vol. 29, pp. 1695-1699. Pergamon Press. Printed in Great Britain SOLUTION OF NONLINEAR BOUNDARY VALUE PROB...

385KB Sizes 1 Downloads 38 Views

Chemical

Engineering

Science,

1974, Vol. 29, pp. 1695-1699.

Pergamon

Press.

Printed

in Great Britain

SOLUTION OF NONLINEAR BOUNDARY VALUE PROBLEMS-VIII. EVALUATION OF BRANCHING POINTS BASED ON SHOOTING METHOD AND GPM TECHNIQUE MILAN Department

of Chemical

KUBB?EK

Engineering,

and

VLADIMiR

HLAVACEK

Institute of Chemical Technology, Prague 6, Czechoslovakia

(Receiued

14 November

Suchbatorova

1903, 166 28

1973)

Abstract-General straightforward method for finding branching points of nonlinear boundary value problems is presented. The technique proposed utilizes the equations used in the GPM method. The method can be applied to problems arising in a number of physical applications, the technique is tailored to diffusion and heat conduction problems.

1. INTRODUCTION

subject to linear boundary

The nonlinear boundary value problems often may exhibit for a given set of parameters more than one solution. The kinds of problems we treat arise in wide variety of applications. In particular, some problems in chemical engineering as, for instance, axial mixing in tubular reactors, heat and mass transfer within and outside of a porous catalyst, boundary layer problems for exothermic reaction, equilibrium of neighbouring drops at different potentials etc. provide excellent motivation for our study. Of course, it is a difficult task to establish all possible solutions if more than one is expected. If the solutions are very close to each other the majority of straightforward methods fail to calculate the relevant profiles because of oscillation of a particular iteration procedure between both solutions. If both neighbouring profiles coincide a branching profile results. Because this profile corresponds in a functional space to a point the term branching point will be used in the following text. A great savings of effort can usually be made if it is possible to evaluate all branching points before calculation of solutions for particular values of governing parameters. It is the purpose of this paper to present a straightforward method for an evaluation of the branching points appearing in a two-point nonlinear boundary value problem. The method of solution will be outlined and employed to solve two practical problems. 2. DEVELOPMENT

@Y(a)

by’(b)

= 7%

(24

= ~1

(2b)

a function of the parameter initial condition

LY.On chasing the

as

y(a)

and supposing

= A

(3)

POZ 0, Eq. (2a) yields y’(a)=+-an&.

(4)

In accordance with our recent paper[l] a set of differential equations for the auxiliary variables

p(x)=w$

aY(x)

q(x)=aa

(5)

can be devloped[2]:

pfdfp +-$ p’

p(a)=

p’(a)=

1,

-2

af af q~~=-4+-rq~+-g. q(a)=q’(a)=O. ay a~ After rewritting

(7)

Eq. (2b) we have

F,(A, a) = a,y(b, A, a) + Pgy’(b, A, a) -

yl =

0 (8)

A necessary condition for a branching point yields: aF,

da zi -=---_& dA aF,

aa (1)

0

(6)

OF THE METHOD

Y, Y’, a)

+ by’(a)

a,y(b)+

2.1 One second order equation Consider the boundary value problem Y” = fk

conditions:

This equation 1695

gives, supposing

aF,laa # 0,

1696

MILAN

KUBftEK

and VLADIMIR HLAVXEK

F,(A,u)=$$=nip(b,A,n)+~,p’(b,A,a)=O Equations (8) and (9) are sufficient the unknown coordinates (A, &) of branching point. For instance, using method to solve (8) and (9) the Jacobian

(9) to establish a particular the Newton matrix r

r=

tions (3), (4) (6) (7) and (15) from x = a to x = b. (3) enumeration of F, and Ft according to (8) and (9) (4) enumeration of the elements of the Jacobian matrix (10) according to (11) and (16). (5) Evaluation of a new approximation of A and (Y using the Newton method

(10)

must be determined. The values of the first derivatives evaluated because

of F, can be

~=a,p(b,A,a)+Plp’(b,A,a) (11) $

where Ae(0, 1 > is a damping factor. (6) The calculation is stopped for ll(A, a)“‘” (A, CX)~‘~[~ < E (e is the prescribed tolerance), if not go to step 2. The method developed will be illustrated on a calculated example. Example 1. Heat and mass transfer in a porous catalyst of the plate shape is described by[5]: Y) $ =Sy exp [ 1+rP(JP(1- -Y) 1

= culq(b, A, a) + P,q’(b, A, a)

y$-0,

To calculate the first derivatives of F2 two additional auxiliary equations for the variables

(12) must be developed. These equations result after differentiation with respect to A and (Y as well:

y(l)=

1.

For y(0) = A and setting QI= 6 we can readily develop a set of differential Eqs. (6), (7), (13) and (14); for instance for r we have:

of (6) I”=

YP

6 [

(

YPY (l+P(l-y)P

(l+P(l-y))’

)

_rP(1+P(1-Y))+W2Y (1 +pu-

st&fpq +-&+P’I’+P%)+&& a’f a2

Apparently,

+-$st.

the initial conditions

are

r(u) = r’(u) = s(a) = s’(u) = 0. Now the first derivatives evaluated:

$ = alr(b,

(14)

(15)

of F2 can be easily

A, a)+ P,r’(b, A, a) (16)

g

= cx,s(b, A, a)+ P,s’(b, A, a).

A general algorithm for the calculation of Luand A is following: (1) initial guess of A and (Y (2) integration of 5 second order differential Eq. (l), (6), (7), (13) and (14) with the initial condi-

I

rP(l- Y) exp [ l+p(l-y)

I n* rP(l-Y) I exp [ l+p(l-y)

YPY

0 +-+J’+&p~+$s a(Y )

YP

+6 [ ‘-(l+p(l-y))Z

I r.

The relevant initial value problems have been integrated by virtue of Merson’s method[6]. Some results of our computations are shown in Table 1. The convergence properties of this procedure are very good. 2.2 Two second order equations We shall develop in this section an algorithm a set of two second order differential equations, general algorithm will be presented elsewhere However, the development for more equations be made in the same way. Consider a set of two differential equations: y:‘= fl(X,

Yl,

Y2,Yl, Yi, a)

y: = ji(X,

Yl,

Y2,Y’l, YL a)

for the [3]. can

(18)

Solution of nonlinear boundary value problems-VIII

1697

Table 1. Course of iteration for single second order equation (y = 20, p = 0.4, A = 1) Iteration

A

a=6

0 1 2 3

0.8000 0.7954 0.7929 0.7928

0

1 2 3 0 1 2 3 4

Ft

F,

0.1000 0.1367 0.1376 0.1376

- 0.04772 - 0~00104 0~00002

0.23304 0.01286 - owOO3

0.2000 0.2218 0.2272 0.2273

0.1000 0.0771 0.0779 0.0779

0.15636 - 0wKxO - OTlOO16

- 0.02447 ow470 0~00103

0.6000 0.8832 0.8024 0.7928 0.7928

0.1000 0.1451 0.1463 0.1372 0.1376

- 0.03936 0.01903 0.01065 - omo4o

- 0.31023 0.24825 - 0.01402 0.00165

subject to linear boundary

conditions

After inserting the values of y and y’ (which are the results of integration of (18), (20) and (21) into (19b) two nonlinear equations can be written R(A,, AZ, a) =

cukydb,AI, Azr a)

+ ~,~y;(b,A,,AZ, (Y)- ylk = 0

(24)

k = 1,2. The necessary condition for the branching point is K(A,,

ao’yi(a) + Poiy’r(a) = yoi

i = 1,2

(19a)

cul’yi(b)+ &‘y’,(b) = yli

i = 1,2.

Wb)

A~,(Y) = det WI,

here T is the Jacobian

AZ, a) ‘0,

(25)

matrix

Let us select the following two initial conditions y,(a)=A,, and supposing

yz(a)=Az

(20)

(26)

l-=

PO’Z 0, i = 1,2, we arrive at

y\(a)=L(yo’-c~o’A~) PO’

i=l,2.

(21)

In the same way as in the one-dimensional case a set of differential equations for auxiliary variables:

am)

Piik =

A!$

qi(x)

Pii T-j-&

a*Ydx) a~,a~k

qij@) =

(22)

zt.l&

This condition results from the implicit function theorem. Clearly, the branching point can be evaluated as a result of solving three nonlinear Eqs. (24) and (25) which yield the unknowns A,, A2 and (Y.If the Newton method is used to solve these nonlinear equations the terms aE/aA,, aTlao: i = 1,2,3; j = 1,2 must be evaluated. To evaluate all these terms the following relations can be developed:

I

i,j,k=l,2;

$=

kaj

can be developed [ 1,4] with the initial conditions = Pij(a)

Pii

=O,

i# j;

pii

= 1,

I

$=u,iq~(b)+Pl’q:(b) aF3 -=_

aAi

p:,(+-$ C&(U)= q:(a) = = qij(a)

pij!f(a)

= q:j(ll)

i, j, k = 1,2,

For instance,

for qij we have

=

P:jk(a)

= 0,

aFI

(23)

(27a)

i = 1,2 aF,

a2F2

a2F,

aA, m+%m-aA, ___aF2 a2F, aA, aA2aAj

cl

i, j = 1,2

cUl’pij(h)+ P,‘p:j(b)

Wb) aF,

a’F2

aA,aAj j=1,2

here, e.g. a2F2 -= aA2aAj

k 3 j.

and

a12p,,j(b) + fl1*pLj(b)r

(%a)

MILAN KUB~~EK and VLADIMIRHLAVACEK

1698

a& -=-

aa

aF, a2F2 aF, a2F2 aF, d2F, aA,ZGG+faA2aA,acvZaA,a~

---aF2 a2F, aA,

(28b)

akhaa'

here, e.g. we have

2

=a,‘q,,(b)+ P,*qL(b).

2

To illustrate the method for a set of two second order differential equations an example will be presented. Example 2. Axial heat and mass transfer in tubular reactors in described by two second-order differential equations [7,8] y - PeMy’ + PeMDaU - Y) em 0” - PeH8’ -

e

(- > 1 + le

=o

Pe&(O - 8,) +PeHDaB(l-y)exp

subject to boundary

(

&

>

=0

conditions

y’(1) = V(l) = 0 PeMy(0) - y’(0) = 0

e(l)=At,

0.1

a05

We denote y(l)=A,,

1

0.0

PeHe(0) - e’(o) = 0.

a=m

I 0.7!

Fig. 1. Dependences A,(a) and A*(a) for Pe, = PeM = 2, /3 =2.& =E =O,B = 12.

a=Da.

Table 2. Course of iteration for the values of parameters PeH=Pem=2,h=1. Iteration

A, = y(1)

A2 = e(l)

a = Da

0 1 2 3

0.2000 0.3759 0.3871 0.3871

1mOO 1.8061 1.8652 1.8652

0~1000 0.1330 0.1319 o-1319

0.00469 0~00003 0~0ooOo

0 1 2 3

0~7.500 0.7933 0.7860 0.7858

3-8000 4.0761 4.0313 4.02%

0.1150 0.1108 0.1115 0.1115

-0.01787 0.00424 OTIOO18

0 1 2 3 4 5 6

09000 0.9713 0.9533 0.9309 0.9207 0.9185 0.9184

4.0000 3.6593 36495 3.8366 3.8911 3.8990 3.8993

0.1400 0.1338 0.1206 0.1230 0.1216 0.1214 0.1214

0 1 2 3 4

0.9800 0.9677 0.970 1 0.9696 0.%96

3.8000 3.3720 3.3476 3.3527 3.3529

0.1150 0.1175 0.1154 0.1155 0.1155

F,

B = 12, p = 2, E = 8, = 0,

F*

F3

0.36868 OTIO138 0@0001

- 099197 -0aO444 -0aKKIo

- 0.02230 OdO815 OdOO56

0.45113 -0.11603 - oaO357

- 0.03548 0.23797 0.06945 0.02078 oaO393 0.00016

0.50492 244774 0.61371 0.22125 0.03630 ool145

- 1.89113 8.58409 3.94400 0.80386 0.11828 0.00421

0.47920 0.00328 - oaO113 - omOO3

3.97732 0.09019 - 0.01242 - 0.00038

20.8162 1.31692 - 040002 - 0~01564

Solution of nonlinear boundary

To find the branching points it is necessary to integrate 18 differential equations of second order (initial value problem). The course of the iteration process is reported in Table 2. The particular branching points are displayed in Fig. 1. where the dependences A,(a) and At(a) are drawn. It can be noted that the rate of convergence of the Newton iteration process is high. 3. DISCUSSION

AND CONCLUSIONS

There are numerous situations in which one is interested only in an evaluation of the branching points, i.e. the determination of boundaries where multiple solutions occur, not in the details of the whole dependences between some parameters and dependent variables. When such is the case there is no need to evaluate these dependences because the straightforward method proposed locates easily the parameters of the branching points in question. So far the location of a branching point was established only in an indirect way, i.e. it was necessary to explore the vicinity of a particular branching point by mapping or other methods (e.g. GPM technique [ 1,4,9]), however, the method suggested has proved effective because of its straight-forward character. Since obviously nonlinear elliptic partial differential equations can be converted by virtue of the ‘method of lines’ to a nonlinear boundary value problem for a set of ordinary differential equations the technique of finding the branching points can be also easily adopted. It should be

value problems-VIII

1699

noted, however, that there is a number of chemical engineering problems which cannot be integrated by shooting methods because of inherent instability of a relevant initial value problem. Clearly, the algorithm described can be also adapted to nonlinear algebraic equations; the boundary value problem can be then converted by making use of orthogonal collocation or high order difference scheme to a set of nonlinear algebraic equations which may be handled easily. The algorithm for algebraic equations will be presented elsewhere[3]. Further applications will include, probably, studies of thermal explosion equations, boundary layer structure, BCnard convection problems etc. REFERENCES

[l] Kubfcek M. and HlavaEek V., Chem. Engng Sci. 1972 27 743. [2] Coddington E. A. and Levinson N., Theory of Ordinarv Differential Eauations. McGraw-Hill, New York 1955. [3] Kubfcek M., .I. Inst. Math. Appfic., to be published. 141 Kubic’ek M. and Hlavacek V., J. Inst. Math. Applic. 1973. [5] HlavSlEek V., Marek M. and Kubfcek M., Chem. Engng Sci. I%8 23 1083. [6] Fox L., Numerical Solution of Ordinary and Partial Diflerential Equations. Pergamon Press, Oxford, (1%2). [7] HlavaEek V. and Hofmann H., Chem. Engng Sci. 1970 25 1517. [S] Hlavafek V., Hofmann H. and Kubicek M., Chem. 1973 12 287. [9] Kubicek M., Communs of the ACM, 1974.

Zusammenfassung-Eine allgkmeine direkte Methode zur Bestimmung der Verzweigungspunkte einer Die vorgeschlagene Technik verwendet die nichlinearen Randwertaufgabe wird angeftlhrt. Gleichungen, die fur die GPM Methode entwickelt wurden. Das Verfahren kann zur Liisung von vielen physikalischen Problemen verwendet werden, der Algorithmus ist zur Analyse der Diffusionsbezw. Warmeleitf5higkeitsprobleme geeignet.