Asymptotic solution of a diffusion problem with mixed boundary conditions

-1

ARCHEF

_.

SdI

T;:

KONINKL. NEDERL. AKADEMIE VAN

_WETENSCHAPP'

Reprinted from Proceedings, Series B, 69, No. 2, 1966 MECHANICS

ASYMPTOTIC SOLUTION OF A DIFFUSION _PROBLEM

WITH MIXED BOUNDARY CONDITIONS

BY

L. VAN WIJNGAARDEN

(Netherlands Ship Model Ba.sin, Wageningen)

(Conmumicated by Prof. W. P. A. VAN LAMMBBEN_{at the meeting of Dec. 18, 1965)}

Summary

òc (Ö2c ò2c1

The differential equation U

_{- = D}

+ - is

_{considered with}

òx Iòx2 òy2)

mixed boundary conditions on y = O: c = CO Ofl O x 1, y = O; òc/òy= O

on x> 1, y=O and x<O, y=O. The solution can be represented by the integral equation

5

eKo{2!xel}!() d=cc0.

An asymptotic solution for large 2 is constructed by a series of solutions of similar problems for _{a semi-infinite strip. It is shown that after two}

steps the error in (òc/òy)o is of order (e/2). _{The present solution is}

compared with a solution by previous _writers. 1. _Introduction

During the last years considerable attention has been paid to the onset of cavitation in a liquid. For the first phase, the growth of small, gas and vapour containing, bubbles to macroscopic size the name "cavitation inception" has been introduced. Vaporous _{growth is marked by}

4a

Pi<Pv

-where

pj = pressure in the liquid outside the bubble Pv = vapor pressure

a = surface tension R= bubble radius.

In water this condition_{can for R =O( 1O m) not be satisfied by pressures} such as occurring in the flow through a ship's screw, and then the growth of bubbles is effected by diffusion of solved gas from the liquid into the bubble. The rate of growth by diffusion in _{a liquid at rest has been} con-sidered by various writers, _{e.g. EPSTEIN and PLESSET i . Exeerimental}

ibIoH,ek van de

£fi

t

-DDCJ:

-/

(1.9)

264

work of KERMEEN and PARKIN [2] dealt with growth in an ambient liquid with non-zero velocity relative to the bubble. In that case, apart from molecular diffusion, convective diffusion plays a rôle. To complete the experimental work an approximatetheoretical approach was developed in ref. [2] to allow for this convective diffusion. In this context the following problem was formulated in ref. [2].

The concentration of air in water c has to satisfy in the x, y* plane, U being the velocity of the water in x" direction and D the diffusion coefficient of air in water,

òc 1ò2c ò2c1

(1.2)

U- =D-.+-1,

with the boundary conditions

y=O, O<x*<L, C=Co

y=

O,x*>L;y=O,x*

(x*2+y*2)l/2 ± oo, C =

The choice of this model to represent the convective diffusion from the flow to a bubble near an impermeable wall is open to question. The present writers objection is that a difficult problem (convective diffusion to a bubble) is replaced by an almost equally difficult problem (diffusion to a strip of finite extent), for which no exact mathematical solution is known to exist at present. A different approach will be dealt with elsewhere. Here we consider approximative solutions of (1.2 )( 1.5). We introduce

x = xL, y = yL, and the Péclet number 2, based on the half length of

the strip (1.6) By substitution of (1.7) C = ue (1.2) turns into ò2u ò2u + = o.

This equation has the "source" solution or Green's function

(1.8) u = Ko{2 Vx2+y2)

where K0 is the modified Bessel function of the second kind. Introducing the "sink" distribution on the strip

f(x) =

_{2 (-)}

265

we find by virtue of (1.7) and (1.8) the following integral representation for c

(1.10) e e

eKo(2

_{V(x_E2+y2} ¡(C) d.}

By (1.10) the conditions (1.4) and (1.5) are satisfied. The boundary condition (1.3) then yields the integral equation for f(x)

11

- co = zlc = - J eKo{AIxCj}/(C) dC.

2

This integral equation can not be solved exactly. In _{the experiments of} ref.

[2] U

10 m/s; L 10 m; D=2 x 10 rn2/.s. From (1.6) we

conclude that A = 0(104). In stead of an exact solution, an approximative solution for large 2 may be attempted. In the next section _{we describe} the method of solution used in ref. [2]. In sections 3 and 4 two alternative methods are given. The results are compared with those of ref. [2]. 2. _{Solution o/the integral equation (1.11) according to Kermeenand Parkin}

The interval is divided in the parts 0x and _{x - 1. In both parts the}

Bessel function K0 is, in view of the large value of A, approximated by

il

(2.1) _Ko{AIxCI}

ez'

V_221xE!

This substition yields

/(C) 1(C)

2(2tA)'/2Lic

= ¡

(xC)1!2

dC + J e2'-

(C__x)h/2

Again in view of the large value of A, the integration in the _{second term} is replaced by

°°

/(x) _{(_x)1/2} The resulting integral equation

X

(2.2) _{2(2A)1/24c = $}

_{d +}

J/f(x)

o

(x)'/2

is subsequently solved by applying a Laplace transformation. The solution, obtained in ref. [2], then is

(2.3) _{/(x) = 4AAc e2}

_erfc(j/)

where

(3.1)

For further reference it is useful to note that for (2Âx)l/2 ' i (see e.g.

MCLACHLAN [3])

(2.5) erfc(2Ax)1/2 f i - J_. +

V22xk 42x

The Bessel function K0 behaves like the r.h.s. of (2.1) only for

aixEI

L

In the neighbourhood of x =E, however, the behaviour is like in Ix - Ei. Since this neighbourhood is expected to contribute considerably to the integral in (1.11), it is difficult to give an a priori justification of the approximation used in ref. [2] and outlined in this section.

3. Boundary layer approximation

The diffusion velocity from the strip outward is proportional to D/y * The convective velocity is U. The diffusion takes place in a layer with thickness ô, determined by the condition that the rate of increase dò/dx* is equal to the ratio of the above mentioned velocities

dâ D

dx* _ôU.

from which we infer

]¡Dx''

v-ti-.

From this expression and the definition (1.6) it follows that the Péclet number 2 can be interpreted as the square of the ratio betweenLand the average diffusion layer thickness.Since 2=0(104), we conclude that along the strip ô L. Then the variations of e in x" direction are small with respect to those in y*direction,and(1.2),may.in dimensionless coordinates, be approximated with

(3.2) _{òy2 òx}

The solution, satisfying c = Co on the x-axis, and e = C atinfinity is y i

(3.3) c

= C + (coc)

erfc

From (3.3) and (2.4) it follows that

(3.4)

-

2()

_{òy ,=o}=

2Ac.

_7X

Although this approximation fails to account for the finite extent of the strip, because e = Co on the whole positive x-axis, on the strip the solution (3.3) will be good far from x= 1. Comparison of (2.3) with (3.4) shows that with the method of section 2 one does not obtain the square root singularity, present in (3.4).

67

For 2x _{1, however, (3.4) is identical with the first term of the expan} -sion of the r.h.s. of (2.3) with the aid of (2.5). In section 4 _{we consider} the complete equation (1.2) and allow for the finite extent of the strip by a successive approximation to the exact solution.

4. _{Successive approximation, using the Wiener - Hopf technique.}

Because Ah/2 (see section 3) can be interpreted _{as the ratio between} the length of the strip and the diffusion boundary layer thickness, it suggests itself, in view of the large value of A(O(104)) to consider first a semi-infinite strip. Then (1.11) becomes

100

(4.1) Ac =

- f

eA1 _{>Ko{AJx-CI}/i(C) dC O<x<o.}

The equation (4.1), which is equivalent to the differential equation _(1.2) together with the boundary conditions for _{a semi-infinite strip, is more} difficult to solve than the simplified equation (3.2). The reason is that with the elliptic equation (1.2) the sources on x>O, y=O affect the con-centration in the whole x, y plane while in the case of the parabolic equation (3.2) the domain of influence of these_{sources is restricted to the right half} plane. In order to define (4.1) for the whole x-axis we introduce

(4.2) _{ti*(x) = /i(x) H(x)}

where H(x) is Heaviside's unit step function, and

(4.3)

(x) = H(x) 5eA(Kox_}/i(C) de.

Together with (4.2) and (4.3), (4.1) yields 1 00

(4.4)

-

_{$ e_Ko{AIx__CI}/i*(C) dC = ¿ldB(x)+çh(x).}

00

This equation can be solved with the WienerHopf technique. _We in-troduce the two sided Laplace transforms

+00

(4.5) F(p)

_{= 5 eP}

fi*(x) dx

+00

(4.6) _{(p) =}

_eP

_{(x) dx.}

Application of the two sided Laplace transform to (4.4) results, _with the aid of the convolution theorem (BREMMER and _{VAN DER POL [4],} p. 39), in

(4 7) F(p)

-

Ac_{+ (p).}

Ac(22)l/2 ¿le

Multiplying both sides with (21p)1'2, yields

(4.8) (22Lp)'I2--'(p)(2Ap)l/2 = 0.

From the definitions of F(p) and (p) it follows, that the singularities of F(p) are in the left half of the p plane, and of (p) in the right half. If we assume, that the singularities of (p) have a real part larger than or equal to 22, it follows that the first term on the l.h.s. of (4.8) is regular in Re (p)>O and the third term in Re (p) < 21, where Re is shorthand

for "real part of".

The essential step in the WienerHopf technique now is to split the second term on the 1.h.s. of (4.8) in two parts, one regular in Re (p) >0 and the other in Re (p)< 21. In the present case, inspection shows that the required partition is

{(22_p)l/2 - (22)1/2}.

P p

Using this, we sTite (4.8) as

(49) (22)1/2 Ac

(2Ap)'/2(p) +

{(21 p)112 - (22)1/2}.

The I.h.s. of (4.9) is regular in Re (p)> 0, the r.h.s. in Re (p) <21. Because these regions have a common domain, they represent the same function, which is regular in the whole p pIane and therefore is a polynomial in

p, J(p) say, so that

(4.10) F(p) = 2 (21)1/2 Ac ± J(p) l/2

p1'2

Because F(p) must be bounded for p - , we conclude J(p) = 0, and from inversion of the remainder of (4.10) we obtain

(4.11) = 24e

This is identical with the solution (3.4) for f(x) of the simplified equation (3.2). For both equations (1.2) and (3.2), the source distribution along a semi-infinite strip that gives C=C on the strip, is the same, although of course the concentrations in the x, y plane are different.

From the fact that also the r.h.s. of (4.9) is identically zero we obtain

Ac Ac(21)'/2

(4.12)

p

p(21p)'/2

For the inversion of (4.12) we make use of a theorem, given in ref. [4], p. 33:

C1=C,_C erre ( j X j o li o I zk. i X Fig. la.

270

If a(p) is the Laplace transform of a(i), regular in c <Re (p) <ß, than a(p/'v) is the Laplace transform of a(iv), regular in vß <Re (p) <rnx; y < 0. Applying this theorem for y = - I to (4.12) we obtain with the aid of a table of transforms (in this case CAMPBELL and FOSTER [5])

2H(x)Ac

°o

2

(4.13) (x) = _{_}

_f

e-z dz = Ac erfc(22IxI)"2.

1.Z {2ÀIzI)'

Because (x) (cf. 4.3) represents the value of c00c on the negative part of the x-axis we have in x<O, y=O:

(4.14) c=c erf V22Ix + e0 erfc VX[x

Where erf{}=1 erfc{}

The situation obtained with this first approximation is schematically given in fig. la.

The conditions for the full problem are represented in fig. lb.

Fig. lb. Boundary conditions of the exact solution.

To obtain a second approximation to the full problem we consider next the situation given in fig. le.

The boundary conditions for this problem are

c=0 atx<l, y=O

c=0 atVx+y2+oo

271

Fig. le. Boundary conditions of the second approximation.

We are interested in the distribution at x 1, y = O. Denoting this distribution with g(x), and shifting the origin to z= 1, we can write

(4.14) /2*(X) =

-

2H(x) 4c (2

1)

+ H(x) g(x).

If we represent the concentration on the x-axis by V(x) H(x), the

integral equation to be solved in this second approximation, is i +

(4.15) _{Ko{2IxJ} f2*() d = V(x) H(x).}

Applying a two sided Laplace transform as in the first approximation, and multiplying with p1/2, we obtain

4eP4c(22'/2) G(

(4.16)

- 7'/2p'/2(22p)'/2(pI8 eZ2dZ +

()1/2 -

x(P)

112

where G(p) and x() are the Laplace transforms of H(- z) g( - x) and V(x) H(x) respectively. The first term has to be separated in two parts, one regular in p> O, the other in p < 22. Whereas this could be done by inspection in the first problem of this section, this is not possible here.

The partition must be obtained by a contour integration, which is

given in the appendix. Here we give the result. The two functions _are

44c(22)1/2 [n Se'dz

i

/

i I Ni(p) =

-1/2 regular in p < 22. (4.17a) LP"2(22P)"2 2p (22_p)h/2 (22)1/2)] regular in p>O. (4.17b) N2(p) 24c(21)'/2 ( 1

-

(22)1/2)1 = 'I2p 2(2_p)1/2

With aid of these expressions (4.16) can be rearranged to

G(p)

(4.18) (p)p'f2_Ni(p) N2(p).

= (2,a_p)1/2

The 1.h.s. of (4.18) is regular in Re (p)>O, the r.h.s. in Re (p)<2.t. Both sides represent the same function, regular in the whole p plane. This function P(p) must be a polynomial in p. From the r.h.s. of (4.18) we obtain

G(p) = (2Ap)1/2 N2(p)+P(p) (2)p)'/2

Because G(p) must be bounded for p -* o, it follows that P(p) _ 0. Hence, by virtue of (4.17b)

2(2))h/2 Ac (1 (2)_p)1I2'\

(4.19) _°(P)

-

1f2 _{p (22)'/2p J}

Making use of the theorem mentioned before and of the table of Laplace transforms, we obtain

=

[('ia

(2')h12erfc(22xDh12]

Substituting this expression in (4.14) and returning to the original abscis where the origin is at the left edge of the strip, we get

12*(x) = - 2H(x 1) Ac +

211(1x) Ac [ e1

1/2

L(1 x)}'/2

(2A)l/2 erfc{2(1 _X)}1/2]. Together with the result (4.11), obtained as a first approximation, this gives for the "sink" strength on the x-axis

(4.20) ç (4.21) 1òc = t*(x) =/*(x) + 12*(x) '2A"112 [1

2Ac(_

(H(x)I1(x-1)}±

x

11(1x)

erfc(2A(1_x))h12}]. (22'r(1 x))"2

If we substitute the r.h.s. of (4.21) for f() in (1.10), thereby extending the integration from - oc to + oc, the resulting concentration c is an exact solution of the differential equation (1.2). Also the boundary conditions (1.3) and (1.5) are exactly satisfied. From (4.21) it follows that in our approximation òc/òy = O on z> 1, y =0. The only condition, that

273

from (4.21) and the asymptotic expansion for large argument of the complementory error function (cf. 2.5) it follows that

òy on x<O, y=O.

Because òc/òy is exponentially small on this part of the x-axis and exactly zero on z> 1, y=O, the solution (4.21) for (òc/òy)o can be considered as a very good approximation to the exact solution for large ).

5. Discussion

In this section we compare our solution (4.21) with the solution (2.3) obtained in ref. [2].

Defining ki(x) as the expression on the r.h.s. of (2.3) divided by

2(22/7r)'/2 zlc, we obtain from (2.3),

(5.1) k1(x) = (22)1/2 eV erfc(22x)'/2 with x>O.

On the strip O<x< 1, we obtain, also dividing the r.h.s. of (4.21) by the above mentioned factor

1

e21

(5.2) k2(x) ₌ +

A(1 )}1/2 erfc(2À(1x)}'/2.

We recall that in the approximation, resulting in (5.1) the Bessel

function K0 is (cf. 2.1) replaced by its asymptotic behaviour for large argument and that the strip is of semi-infinite length. Comparing (5.1) in Ox< i with (5.2), we can discern three regions

O<x.

Expanding both ki(x) and k2(x) for these values of x, we obtain

(5.3) ki(x) (2.P)1/2{i1/2 - 8(2x)1'2 + O(Ax)}

(5.4)

Because K0 {Ax ¿I} behaves for small argument as ln x j rather than as _J}1/2 the solution for 1(x) (and hence for k2(x)) has a square root singularity in x = O. The use throughout the strip of the asymptotic representation (2.1) results (cf. 5.3) in a solution which does not display this singularity and behaves also otherwise quite different from (5.4). In this region therefore the solution (5.1) is seriously in error.

Here the complementory error functions both in (5.1) and (5.2) may be replaced by their asymptotic behaviour (cf. 2.4) for large arrgument:

i

(5.5) ki(x)

'-

+ O(2'x3/2)

1

(5.6)

112

In this region both solutions coincide as far as the first term is concerned.

(1x).

In this region we obtain

(5.7) k(x) 1(1---x)+O(A-').

(5.8) kz(x)

22(1x)}1/2 + O {21/2(1_x)l/2}

Because of the finite width of the strip, k(x) possesses a square root singularity in x= I, which is not present in (5.7). Therefore, also at the right edge x= 1, the solution (5.1) is in error.

Finally we note that although /ci(x) differs considerably from k2(x) both near x = O and near x = 1, this does not lead to serious errors in the investigation of ref. [2], since there only f 1c1(x) dx is used. It is readily verified from the above results that the contributions from ki(x) and !C2(x) to the integrals S k1,2 dx are in region i) of order 2_1/2, in region iii)

of order 2' and in region ii) of order 1.

The main contribution therefore stems from region ii) where both functions coincide as far as the first term is concerned.

Appendix

4Ac(22)1/2

we have to split the function

Apart from a factor 1 /2

e,o

(Al)

K(p)

= 1/2(2) )l/2 5 eZ2 dz

in two functions, one regular in Re (p)> O, the other regular in Re (p) < 22. Consider in a complex w plane the point w=p and the lines Re (w)=d,

Re(w)=c such that O<c<Re(p)<d<22.

Let the points A and B be located in c ± ja, the points C and D in d + ia. By virtue of Cauchy's theorem

iK(w)dK()

2n wp

along ABCD. (see figure 2).

275

wp

Fig. 2. Contour integration in w plane.

Ifa - oo, the integrals along AD and BC vanish owing to the behaviour

of K(w)/w p at infinity, so that

J ( c+ioo cl-joe

K(p)=.-H $ +

S

c-joe d+ioe

We close both paths by a semi-circle and the path E1 E2 E3 along the branchcut extending from the branch point w = 22 of K(w) to infinity along the real axis.

Then by Cauchy's theorem

i c+ioe _i

5 = Residu in w p

_-

_$

-

c-joe EiE2E3

1d-oo1

The integral along E1 E2 E3 is

i e2 _{ds 5;+8)} e-z2 dz

sf(2 ± )1f2 (2 - p + s)

Because .. is large we may cf. (2.5) approximate the second integral by

e -2%

(A-5)

276

2(2) + 8)h/2

Then (A-5) reduces to

i _f"° ds "? ds i i i i

I

s'(2 + sp)

sh'2(2A + s)J - i (2_ p)'!2 (2)I2J

Evaluating (A-3) and (A-4), with the aid of this result, and multiplying 4z1 c(22)'/2

with the factor - leads to the functions Ni(p) and N2(p)

1/2 in (4.17).

REFERENCES

EPSTEIN, P. S. and M. S. PLESSET, On the stability of gas bubbles in liquid-gas

solutions. Journal of Chemical Physics, 18, 11, 1505 (1950). PARKIN, B. R. and R. W. KxiurEN, The roles of convective air diffusion and

liquid tensile stresses during cavitation inception; Proc. Jahr sympo-sium on cavitation and hydraulic machinery Sendai, Japan, (1963).

MCLACHLAN, N. W., Complex variable theory and transform calculus, Cambridge

at the University Press, (1955).

PoL, B. VAN DER and H. BREMMER, Operational calculus based on thetwo.sided

Laplace integral, Cambridge at the University Press, (1955).

CAMPBELL, G. A. and R. M. Fosvm, Fourier integrals for practical applications,