Chemical
Engineering
Science,
1974, Vol. 29, pp. 1457-1464.
Pergamon
Press.
Printed
in Great Britain
TRANSPORT IN COMPOSITE MATERIALS: REDUCTION TO A SELF ADJOINT FORMALISM D. RAMKRISHNA Indian Institute of Technology,
and NEAL
R. AMUNDSON
Kanpur 16, U. P. India and University Minnesota 55455, U.S.A.
of Minnesota, Minneapolis,
(Received 24 September 1973; accepted 20 November 1973) Abstract-A general formalism is presented for the solution of boundary value problems arising in the description of transport of mass or energy through multiphase media, which consists in identifying physical conditions of interfaces with criteria for self-adjointness of the operators concerned. Several examples have been presented which demonstrate the advantages of the formalism. 1.INTRODUCTION The transient
multilayered
2.THEORY
Consider the Sturm-Liouville
analysis of heat conduction through a composite solid medium has long been
of interest to engineers [l-2]. The traditional approach to the solution of such problems has generally been through the use of Laplace transformation techniques [3]. In finite geometry, however, the problem submits neatly to solution by what is popularly known as the method of finite Fourier transforms which is essentially an expansion of the solution in terms of eigenfunctions arising from associated self-adjoint operators. That this is so has not been recognized, or if recognized has not been adequately exploited by engineers. Of course this procedure is entirely equivalent to the method of Laplace transforms since the poles of the Laplace transform of the solution are directly related to the eigenvalues of the self-adjoint operators referred to above in the method of finite Fourier transforms. The method of eigenfunction expansion has the merit of being more direct, however. Besides, this approach has led the authors to provide elegant solutions for a variety of problems whcih are interesting both from the engineering and mathematical standpoints[4]. We will present the theory of onedimensional second Sturm-Liouville order operators with discontinuous coefficients, the extension to higher dimensions requiring no special imagination. The restriction to second order operators is of course due to the nature of the engineering problems under discussion.
operator
l d [p(x)&]+q(x)u, O<x{U’(X)U(X)- u(x)u’(xNl: = 0.
(2)
Ordinarily Eq. (2) would suffice for the symmetry of L, i.e. when p(x) is continuous in (0,l) and the inner product of two functions u and v is defined by I
(u, u) = / r(x)u(x)u(x)dx.
(3)
0 where the restriction to real valued functions is obvious. Suppose however that p(x) is everywhere continuous except at a number of points xl, XZ, . . . , x. in (0, l).t Then it is readily shown that (Lu, v) - (4 Lv) =
z
x {u’(x)u(x)
[P(X) -
u(x)u’(x)}l:j:oo
(4)
Thus if we impose the conditions tThe restrictions on r(x) and q(x) are not stated but are found in any standard theoryt51.
treatment
U(Xi +o)=
U(Xi -0).
(5)
of Sturm-Liouville
p(xi +o)u’(xi +o) = P(Xi _O)U’(Xi -0). 1457
(6)
1458
D.
RAMKRISHNA
and NEAL R. AWUNDSON
on u and v then clearly the operator L is symmetric.? Conditions (5) and (6) together with homogeneous boundary conditions at x = 0, and x = 1, and absolute continuity of u(x) and pu’ define a dense linear subspace of % [O, 11, denoted 9. For u, v& we have from the symmetry of L (Lu, v) = (a, Lv).
=-T(x)%
ld[
“I
p(x)-&u
+q(x)u
-
5 [ai+,p (xi + O){u(xi t- O)v’(xi + 0)
i=l
- U’(Xi +o)v(xi +o)} -aip(xi -O){u(xi -O)v’(xi -0) - u’(x, -O)v(x, -O)}l.
(7)
It is not difficult to show the existence of a symmetric Green’s function to represent the inverse of L. Since the integral of the Green’s function is a selfadjoint operator on Zz[O, l] we deduce from the spectral theorem[6] that L has a countably infinite number of eigenvalues and a complete orthonormal family of eigenfunctions forming a basis in Zz[o, 11. Also that the eigenvalues are non-negative for nonnegative-valued functions q(x) is inferred from the positivity of (Lu, u) for ME% The characteristic equation for the eigenvalues is obtained by solving the differential equation Lu
(Lu, v)-(26 Lv)= a,+lp(l)[u(l)v'(l) - u'(l)v(l)l
=Au.
(10)
On careful examination of (10) it will become evident that subject to boundary conditions at x = 0 and x = 1 for u and v such that a”+,P(l)[u(l)v’(l)-
u’(l)uU)l - a,p(O)[u(O)v’(O)- u’(O)v(O)] = 0.
(11)
and discontinuity conditions a,+,u(x, +0) = aiu(xi -0).
(12)
p(xi+o)u’(xi+o)=p(xi-O)U’(Xi--0).
(13)
i=l,2,...,n
(8)
the operator
L is symmetric,
i.e.
separately in each of the intervals
(ximl,xi) and substituting the result into the boundary conditions and conditions (5) and (6). From this procedure will ensue a set of linear homogeneous algebraic equations in the coefficients associated with the linear combination of independent solutions of (8) so that the characteristic equation is merely a statement of the requirement that the above homogeneous equations have no trivial solutions. The foregoing analysis fits naturally into the solution of conduction problems in composite solid media. Prior to demonstrating this through an example, we would like to consider a small variation of the above Sturm-Liouville problem which appears necessary to consider problems of mass diffusion in multiphase media. Thus we consider the operator (1) with the inner product
(u, v) = g ai [ r(x)u(x)v(x)dx. XI?,
(9)
(Lu, v) = (u, Lv)
whenever u, v belong to the subspace defined by conditions (1 l), (12) and (13). As before, we recognize that these conditions are sufficient for the symmetry of L but not necessary. The simplest way to realize (11) is to have unmixed boundary conditions at Oand 1. Thus, through the existence of a self-adjoint completely continuous inverse of L, we infer that L has infinitely many eigenvalues and a set of orthonormal functions complete in Y*[O, 11. The orthogonality of these eigenvectors must of course be considered in the light of the inner product (9). As before the eigenvalues are non-negative for nonnegative q(x). We consider some examples illustrating the use of the aforementioned formalism for the solution of boundary value problems. 3. HEAT
wherexo = 0, x.+, = 1 and the ai’s are strictly positive quantities. The positivity of ai and r(x) suffices to admit (9) as a valid inner product. It is readily established that tit should be noted that conditions (5) and (6) suffice for symmetry of L and are not as such necessary.
(14)
COMPOSITE
CONDUCTION IN A SLAB OF ” LAYERS
For the sake of simplicity, we restrict ourselves to one-dimensional unsteady state conduction in Cartesian coordinates. The composite medium is of unit length consisting of n layers of materials of varying thermal properties, the interfaces occurring at x,, x2,. . . , x,. The material between xi_, and xi
1459
Transport in composite materials: Reduction to a self adjoint formalism has thermal conductivity ki and heat capacity per unit volume ci. The unsteady state heat conduction equation is written as:
1
-$-& [ k(x)Z =g;,
ocx o T=O,
x=1.
(17)
The interface conditions shown to imply that
(18j
A=
(5) and (23) are easily
and the initial condition T(x, 0) = To,
0 5 x 5 1.
Since the temperature must be continuous of the interfaces xi we have T(xi +O, t) = T(x, -0,
The continuity
t).
at each
(19)
of heat flux at the i”’ interface gives
ki+,$$x,
+ 0, t) = k$xi
- 0, t).
cos J($)Xi
+ [sin J(i)
cos J($+
-+J(z)
cos J(-+-)xi
cos
sin J(&)xi]Bi+,. (27)
Bi= [cos J($)xi
sin J(&)xi sin J&xi
cos J(&)xi]A+,
(21)
+[cos J(-$ cos J(-J&+
u(l) = 0. -0).
ki+,u’(xi +O) = ki-,u’(xi -0).
(23)
From the previous section, it is clear that L is a symmetric operator with a countably infinite number of real positive eigenvalues {A.} and eigenfunctions that are orthogonal with respect to the inner product I (U,U)=lc(x)U(x)V(x)dx=~ci 0
+?J($)
%d(z)
U’(0) = 0.
+o>= u(x,
sin J(&)xi
(20)
L = - [l/c(x)] operator (d::) x z(x;(ed/dxa); operating on functions u belonging to a subspace of JZS[O,11 defined
U(Xi
[sin d($)xi
?J(z)
sin j/&xi
sin J(&)xi]Bi+b (28) i = 1,2 ,....
n.
Written in matrix form one obtains [$]
= G*J$;;].
(291
1 u(x)v(x)dx. (24)
where the coefficients of the (2 + 2) matrix I&i are evident from Eqs. (27) and (28). (For a single
D. RAMKRISHNA
1460
and NEAL R. AMUNDSON
medium rui = ai+,, ki = k,,, and Kh.i is readily seen to be 5 which is a partial verification of (27) and (28)). Repeated application of (29) gives
Before applying the specific boundary conditions (16) and (17) we pause for a generalization. Mixed boundary conditions at x = 0 and x = 1 would together produce another equation
The evaluation of the coefficients KII,Aand K12,1 of the matrix would involve some complex algebra but they are not beyond representation by compact algebraic expressions for a limited number of discontinuities. In any case for numerical evaluation of eigenvalues such expressions are not necessary Once the eigenvalues are determined the eigenfunctions are readily available from (26) through (29). We represqnt the normalized kth eigenfunction, uI, by
I
x IPk, xi-,
(31) Combination tic equation
of (30) and (31) yields the characteris-
where \I$ - i$) represents the determinant of the matrix (I$* - I$). Since the purpose of the present work is to expound a convenient formalism for the solution of the boundary value problems under discussion, we will avoid discussion of computational aspects such as evaluation of eigenvalues from Eq. (32). For the unmixed boundary conditions (16) and (17) the characteristic equation is obtained as follows: From (16), it is clear that Bi = 0 or Ku,, A.+1 + Km* B.+I = 0.
(33)
while Eq. (17) implies that sin
-&
J(
>
(34)
(37)
where
4
+ 2%(hk) cos
(38)
To solve the boundary value problem (15)-(20) we rewrite the problem in the operator formt LT=-;iiT,
d
T(l)&,
t >o
(39)
T(0) = To
(40)
where 9 is the domain of L which includes the boundary and other conditions associated with the points of discontinuity. Taking the inner product of Eqs. (39) and (40) with Uk one gets (LT,
A,+, +cos
i=l,2,...n+l.
x I xi
5
uk)
=
(T,
bk)
=
hk(T,
Uk)
(T(O), uk)
=
=
(To,
-&T>
uk)
(41) (42)
uk)
The pair is readily solved to obtain so that the eigenvalues
are to be evaluated
form (T,
I
I
(35)
or tan
(36)
tin using the operator form, we deliberately suppress the argument x in writing the dependent variable T; it is for the same reason that the derivative w.r.t. time is written as an ordinary derivative.
uk)
=
(To, Uk)‘+’
(43)
Thus the temperature is given by
profile at any time in the slab
T(x, t)= To$sg
ci/XX+ [A,(Ak)sin J(l)& I I
+ Bi(Ak) cos Ai (At) sin
J(
$
,>
X,-l 5x
i=l,2,...,n+l
x+B;(A,.)cos IX,
4
A,, - x , (Yi>I (44
Transport
in composite
materials: Reduction
It is easy to see that the solution may also be represented in terms of coefficients A,+, and B,,, by using the matrices Kh.i through Eq. (29). We conclude this discussion with the observation that the problem in which the thermal conductivity varies continuously as a function of position may be approximately solved by substituting for k(x) a suitable linear combination of step functions such as $ k(xi)H(x
-xi) where H(x) = [;I
1461
to a self adjoint formalism
Besides at x = x, we assume that c(x,-0,
t) = c(x, +o, t).
-0, t) = - DzA+(x,
- D,A,$(x,
+o, t).
(48) (49)T
The differential expression in (45) is not formally self-adjoint. In any given section one may write
,” z,”
By reducing the number of discontinuities quick estimates of eigenvalues may be obtained for the original operator.
introduce the transformation ~e’“*~‘*~) = c, xi-, < x <xi, into Eq. (50) which becomes
We
4. AXIAL DISPERSION OF A SOLUTE IN A TUBE OF VARYING CROSS-SECTION
At first it might seem somewhat redundant to consider this example separately. In fact, while it does share several common features with the previous example, it brings out one interesting variation which relates to the fact that a further modified inner product is required in the function space for a self-adjoint problem. To simplify the discussion, we assume that the tube has two sections of differing uniform cross-sectional area in each of which the axial dispersion model holds (with of course different values of dispersion coefficients in the two sections) for the distribution of a solute introduced with the feed. The mass conservation for the solute is described by the following differential equation in solute concentration c. -- 1
a
A(X) ax
D(x)A(x)$-
v(x)A(x)c
1
= $.
(45)
(51) The boundary
conditions
- D,$O,
t)++O,t)
D+,
t)+%u(l,
The conditions
to = u,cf.
t) =O.
(52)
(53)
at x = x1 become
u (x, - 0, t)e’“PWl = u (x, + 0, t)e’~‘*W~. (54) so that the function u is actually discontinuous at x = x,. The continuity of the mass flux transforms to - D,A, $x,
- 0, t)e’“~‘2D~‘“~ = - D,A,$(x,
where
transform
+ 0, t)e’@4”1.
(55)
Consider the operator
The Danckwerts boundary x = 1 are given by
-+O,
conditions
t)+u,c(~,t)=u,c~.
at x = 0 and
(W
xi-, < x < xi.
(56)
with a domain 9 that is a linear subspace Y2 [0, 11, comprising functions which satisfy - D,f’(O) + Tf(O) = 0.
E(l,
t) = 0.
of
(57)
(47) DJ’(l)+?f(l)=O.
tin writing (49) we have used (48) and incompressibility of the fluid which implies u,A2 = u2AZ.
which originate from the boundary
(58) conditions
(52)
D. RAMKRISHNA and NEALR. AMUNDSON
1462
and (53). Further conditions
f must satisfy the discontinuity
f(xl-o)e”lXJZD~ = f(x, + 0)e”2”1’*4. D,A$(x,
The boundary value problem (51)-(S) be rewritten in the operator form
(59)
du
Lu =-Z’
may now
t >o.
(65)
- O)e”@D~ = D2Azf’(x, + O)e”+1’*4. (60)
which arise from (54) and (55) respectively. If one defines the inner product between functions f and g as
two
where ZJ(t) is not an element of the domain of symmetry of L. However it is readily established that for every fe 9 (Lu, f)=(u,
1
Lf)-Ale”l”@lvlc~f(0).
(66)
=I
Cf,8) = j- A(x)f(x)g(x)dx 0
= j- A,f(x)g(x)dx
Eq. (65) is to be solved subject to an initial condition
0
I
+
I
Azf(xMx)dx.
(61)
u(0) =
XI then it is easy to see that L is not symmetric with respect to this inner product. It is however easily verified that the use of the modified inner product
h (x)e-‘“JZDt”,
I
~(x)e-‘~‘*4’“,
0< X XI <X