mills mixed bc paper

Solution of Mixed Boundary Value Problems by Integral Equations and Methods of Weighted Residuals with Application to Heat Conduction and Diffusion-Reaction Systems Author(s): P. L. Mills and M. P. Dudukovic Source: SIAM Journal on Applied Mathematics, Vol. 44, No. 6 (Dec., 1984), pp. 1076-1091 Published by: Society for Industrial and Applied Mathematics Stable URL: http://www.jstor.org/stable/2101396 Accessed: 25/04/2010 22:44 Your use of the JSTOR archive indicates your acceptance of JSTOR's Terms and Conditions of Use, available at http://www.jstor.org/page/info/about/policies/terms.jsp. JSTOR's Terms and Conditions of Use provides, in part, that unless you have obtained prior permission, you may not download an entire issue of a journal or multiple copies of articles, and you may use content in the JSTOR archive only for your personal, non-commercial use. Please contact the publisher regarding any further use of this work. Publisher contact information may be obtained at http://www.jstor.org/action/showPublisher?publisherCode=siam. Each copy of any part of a JSTOR transmission must contain the same copyright notice that appears on the screen or printed page of such transmission. JSTOR is a not-for-profit service that helps scholars, researchers, and students discover, use, and build upon a wide range of content in a trusted digital archive. We use information technology and tools to increase productivity and facilitate new forms of scholarship. For more information about JSTOR, please contact [email protected]. Society for Industrial and Applied Mathematics is collaborating with JSTOR to digitize, preserve and extend access to SIAM Journal on Applied Mathematics. http://www.jstor.org SIAMJ.Apri.Mari: (©194SotforndaandAppidMatern NatNo,December964 om SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS BY INTEGRAL EQUATIONS AND METHODS OF WEIGHTED RESIDUALS WITH APPLICATION TO HEAT CONDUCTION AND DIFFUSION-REACTION SYSTEMS* P.L.MILLSt AND M.P.DUDUKOVICE Abstract. Integral equation methods areconsidered asanalternate approach forthesolution ofmixed boundary value problems inheat conduction ordiffusion with chemical reaction that can bedescribed by dual ortriple series equations. The approach isbased upon reduction ofthedual ortripe series toa single ora setofFredholm integral equations ofthe frst kind whose kernel K(x,y) and forcing function f(x) are presented interms ofanappropriate infinite series which canbederived from theseparation-of-variables solution. The solution oftheintegral equation fortheunknown function g(y) isobtained byamodified ‘quadrature method which accounts forthepresence ofa logarithmic singularity inthekemel for x=ywhich isacommon occurrence inproblems ofthis type. The approach iillustrated bysolving selected example problems involvingeitherdiffusionandeactionorheatconduction.Comparisons aremadeoexactsolutions where available and also toother approximate solutions based upon themethod ofweighted residuals (MWR). The results ofvarious numerical experiments suggest that theintegral equation approach yields results ofthesame orsuperior accuracy and with anorder-of-magnitude less effort than those based upon MWR. ey words. integral equations, method ofweighted residuals, heat conduction, diffusion, reaction, dual 1,Introduction, Certain types ofheat ormass transfer processes which areoften associated with chemical ormechanical engineering analysis can result inmixed boundary value problems. These problems arise whenever thedifferential boundary operator contains adiscontinuity due tothespecification oftwo ormore different types ofboundary conditions along acommon boundary. Some particular applications which suggest thevariety ofsystems that have recently been examined include: (i) diffusion and reaction inporous catalyst pellets which have both partial internal and external wetted surfaces [16}, [18], [32], (35), (371, [421, [51], (ii)diffusion ofreactants toasolid surface which contains anonuniform distribution ofcatalyst [30], (iii) determination ofeffective diffusion coefficients forspherical catalystpelletsbyinterpre-tation ofeither dynamic [10] orsteady-state [34] tracer gasresponse measurements, (iv) analysis ofconduction heat transfer between solids which have partial contact along acommon boundary [11], [43}{45], and (v)fluid mechanical modeling ofa particular type ofStokes flow phenomena [50]. Problems ofasimilar nature arealso common inother branches ofscience andengineering such asfracture mechanics [31], [41], [55]and electrostatics [7],toname afew. Additional references regarding the origin and application ofvarious other mixed boundary value problems aregiven by Kelman [27] Despite thewidespread application andgrowing interest inmixed boundary value problems (hereafter represented asmixed BVP) that has developed over thepast eighteen years since theappearance ofSneddon’s monograph [48], nounified method ofsolution tothevarious subclasses ofproblems, such asdual ortriple orthogonal series and dual ortriple integral equations, seems toexist. Inthecase ofmixed BVP *Received bytheeditors December 29,1982, andinfinal revised form February 6,1984 +Corporate Research Laboratory, Monsanto Company, St.Louis, Missouri 63167 +Chemical Reaction Engineering Laboratory, Department ofChemical Engineering, Washington Uni- versity, St.Louis, Missouri 63130 1076 SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS 1077 which aredescribed bydual ortriple orthogonal series, forexample, closed-form expressions fortheseries coefficients exist only forafewspecial forms oftheserieskernelandmodifier(cf.KelmanandKoper[24]forthedefinition ofthis terminology). These special cases apply generally toproblems defined onsemi-infinite domains and arenot applicable tomany ofthemore recent physical application problems cited above from thechemical ormechanical engineering literature which have been posed formainly finite domains. Inthefew cases where thephysical problem isdescribed bydual orthogonal series which have aclosed-form solution, special algorithms such asthose developed byKelman andco-workers [25], (26], [28]must beapplied toobtain numerical results. The need to obtain numerical results for various mixed BVP which occur in chemical andmechanical engineering applications hasresulted inanumberofalternate approaches. These include theapplication ofmore orless standard finite difference methods (16), [30], [32], [42], [45], theartificial interface method (9),[20}{22], and theleast squares method [24], [35]. The least squares method hasalso been treated as asubset ofthevarious methods ofweighted residuals formixed BVP described by dual ortriple orthogonal series [38], [39]. The extensive numerical testing oftheleast squares method byKelman and Koper [24] onvarious dual orthogonal series aswell astheconvergence theorems established byFeinerman and Kelman [12] suggest that thismethodisperhapsthemostreliableoneusedtodate.Acomparison oftheresultsproduced bythevarious weighted residual methods (Galerkin, collocation and least squares) for anumber ofdual and triple orthogonal series arising incertain chemical engineering applications supports thisconclusion (38), (39). Akeyresult oftheextensive testing cited above fortheleast squares method on diffusion-reaction problems isthat rather large linear systems ofequations (N= 100 insome cases) must besolved before anaccurate approximation tothedual ortriple orthogonal series coefficients isobtained. Although very good approximations for certain problems has been demonstrated forsmall values ofN,forexample, N10 (cf.Kelman (24), (27}), more precise results foraquantity, such asafluxevaluated at theboundary forsmall values oftheseries discontinuity contheinterval 0=x=1, often require larger values ofN. Our recent growing interest indevelopment ofanalternate approach forsolving some particular dual orthogonal series fordiffusion-reaction systems asameans of verifying theresults obtained bytheleast squares method hasbeen onemotivation for thepresent study. This approach isbased upon reduction ofthedual orthogonal seriestoasingleintegralequationofthe first kind whose solution may beobtained numerically byapplication ofaquadrature technique. The success ofthe method relies upon summation oftheparticular kernel series which arises ineach application toasuitable closed-form. Itisespecially worth noting that these kernel series are not the type encountered inclassic dual orthogonal series. This isbecause they contain certain key dimensionless groups which appear inthekernel series asaresult ofreaction source terms orthe specification ofmixed Dirichlet, Neumann, orRobin type boundary conditions. Evidence isgiven, based upon theresults ofseveral example problems each ofwhich has aphysical significance, that this integral equation approach iscapableofyielding more accurate results with lesscomputational effort when compared tothemethod ofweighted residuals. The merits ofusing thisapproach forthesolution oftriple orthogonal series isalso given. The possibilities forimplementation ofthe methods onsmaller computing devices such ashand programmable calculators or microprocessors which has already been performed previously forthe least squares method [24] aresuggested. 1078 P.L.MILLS AND M. P.DUDUKOVIC 2.Formulation oftheproblem and integral equations. The example problems to betreated later inthepaper arebased generally upon atransport equation which has thefollowing form: Qn Vu-cu=0 inD where Disaregion enclosed bythesmooth boundary and (2.2) aD=¥aD, isused torepresent theboundary along which themixed ordiscontinuous boundary conditions occur. The parameter mmay assume values ofeither m=2 orm=3, corresponding todual ortriple orthogonal series respectively. The parameter cmay beassigned zero, positive, orimaginary values. Laplace’s equation isobtained forc=0 which iscommonly used todescribe steady-state diffusion ofheat ormass inisotropic media. Thecase where cisimaginary yields theHelmholtz equation where |c|represents thewave number. Thecase! where c>0 istheonewhich iscommonly encountered inthetheory ofdiffusion and reaction and isused todescribe theconcentration ofa species which disappears byafirst-order reaction atany location inDofacatalyst. ‘The quantity wtypically isused todenote anappropriately defined dimensionless quantity such asadimensionless temperature w=T/T, inthecase ofheat conduction, oradimensionless concentration u=C/C, inthe case ofmass transfer orreaction. The subscript 6denotes bulk conditions such asthose conditions encountered inthe absence ofany finite heat ormass transport resistance which otherwise might exist on theboundary aD. The mixed ordiscontinuous boundary conditions which arecommonly encoun- tered inapplications involving equation (2.1) areassumed inthis work toconsist of acombination ofDirichlet, Neumann, orRobin boundary conditions which canbe represented ingeneral as: 23) Aw+ y=, onaD, ‘an The parameters A,,vjand w,arerelated tocertain physical quantities and aredefined ‘on(0,00)with each having adifferent value, some possibly zero, oneach @D, Details ofthemeaning ofthese various parameters and their physical significance isomitted forbrevity butcanbeobtained bytheinterested reader inmonographs such asthose ofOzisik [40] forheat conduction, orAris [2]fordiffusion-reaction, toname afew. Specific examples of(2.1) and theassociated mixed boundary conditions expressed by(2.3) forthecase ofdual orthogonal series aregiven byKelman [20], Collins (6), and Mills and Dudukovié (37] corresponding toc=0, c=ia,and ¢>0 asdescribed above. The approach used byKelman and Koper [24] inapplying the least squares method todual orthogonal series, and thepresent authors forboth dual [38] ortriple [36], [39] orthogonal series when using one ofthemethods ofweighted residuals, is toutilize theform ofthesolution which satisfies thepartial differential equation (2.1) and allbutthemixed boundary conditions (2.3). Such aform can often befound by analytical methods such asseparation-of-variables andcanberepresented byaninfinite "The appropriate designation of(2.1) forthiscase hasbeen thesubject ofsome debate [29]. Ourpreferenceistodenoteitasaformofthe generalized Emden-Fowler equation asgiven inMehta andAris (33, p.612}, where thechoice offlu)=uismade. SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS 1079 series written generally as (2.4) CE0)=EaaXaln)Wal6) where X,() and w,(€) aretheeigenfunctions associated with thesolution uin dimensionless (usually) coordinates (n,£).The a,arethedual ortriple orthogonal series coefficients and aretheunknowns tobedetermined. Application ofthemixed boundary conditions given byequation (2.3) toequation (2.4) leads tothefollowing setofeither dual (m=2) oftriple (m=3) orthogonal series which must besolved for the a,: (2) EdatnthalE)=H(E)s §-1<E<E j=1m Intermsofthenomenclature suggested in[24],w,(€)istheserieskernels, f(€)are theprescribed functions on@D,while the,;aretheseries modifiers. Forthecase outlined above, these are obtained byevaluation ofX,(n) and itsderivatives as indicated by(2.3) fortheappropriate values of7along theboundary aD.Appropriate dimensionalization of£ensures that £=0 and &,,=1sothat theseries kernel w,(£) isorthogonal over (0,1)with respect tosome weighting function w(é), butnotover asubinterval (g,, &),ingeneral. Formation oftheresidual ¢from (2.5) followed bydefining theappropriate inner product allows thevarious methods ofweighted residuals (Galerkin, collocation, and least squares) tobeapplied. The development ofthefinal setofworking equations forthese methods isgiven elsewhere [24], [38], towhich theinterested reader isreferred fordetails. Allofthemethods lead tothefollowing infinite system ofalgebraic equations whose unknowns aretheseries coefficient aj: (2.6) LaAy= Fy7=0,1,-++, 00. i Expressions fortheinner products andmatrix elements Ayand F,foreach method aresummarized inTable |forreference. Closed-form expressions fortheAyandF, can beusually obtained byperforming apiecewise integration over (£o,gm). These have been tabulated foravariety oftrigonometric series kernels, such ascosnax and sinnzx and Legendre polynomials P,(cos @)[27], [38], [39]. Computer software pack- ages which apply theleast squares method [13] or,more generally, oneofthemethods ofweighted residuals (Galerkin, collocation, orleast squares), toauser specified series kernel and series modifier [38], [39] arealso available. This latter software package was used bythepresent authors toobtain some oftheexample problem results which aregiven later inthepaper. Tame 1 Summary ofmatrix expression forMWR. a a Leastsquares 9a,0"° jWCE)raBenp(EWE)AEjfWEE)MhiENE)AE Galerkin ke,(OO) j*weptEG(6)dEJ*wl)(E(€)aE CollocationKe,w(6)8(E- &)) “Hui) KE) 1080 P.L.MILLS AND M.P,DUDUKOVIC ‘Analternate approach, which istheoneofprime interest inthispaper, isto convert thedualortriple series expressed by(2.5)toasingle orsetof‘integral equations, respectively. Thedevelopment given below corresponds specifically todual series equations withtreatment oftriple series being similar. Returning toequation (2.5), it isassumed that oneoftheseries corresponding toaparticular value ofjcanbe ‘extended intotheremaining interval denoted byindex kbyintroduction ofanunknown function gi(€) (tobedetermined) asillustrated below: s Gr G1<8<b 21 eannpite {50 & en Zeta), bei<é<byK=12kAD Formation oftheinner product (¥,(£), Wa(€)) over(£0,2)yields thefollowing Fred- holm integral equation fortheunknown function ACH a(28) JKy(éa6)d=fl)hy) Theappropriate integral equation kernel Ky(é, €')andforcing function Ay(€) are expressed asinfinite series composed oftheseries kernel y,(¢) andseries modifiers tyandjyasshown below: nits )n(€) 29) Kyl66)=Ye ae WEED FTE WRe)de 6 (2.10) hy)=jGENK&€)aE. en ‘Thedualorthogonal series coefficients a,canbeobtained fromthefollowing, relation provided thattheintegral equation expressed by(2.8)canbesolved toyield gx(é): eu all’ »agf° ae (2.11 Gn=—FE Gbde. (En€aesBCE)UnE)aE Po207s 3.Limitations ofintegral equation methods. Thebasic notions ofsolving dualortripleorthogonal seriesbycontinuation ofoneoftheseriesintotheremaining portion oftheinterval byintroduction ofanunknown function isoutlined bySneddon [48].Thisapproach oftenreliesuponaningenious selection fortheformofthecontinuation function g,(é) sothattheresulting kernel series canbesummed toaclosed formand theintegral equation hasanexact inverse. Ourexperience inapplying themethods outlined inSneddon’s monograph [48]tocertain mixed boundary value problems for finite, two-dimensional domains having discontinuous boundary conditions hasshown thattheresulting integral equations arenotoftheclassical type. Forthese cases, closed-form inverses arenotreadily obtained. However, forproblems where asepar- ation-of-variables solution having theform of(2.4) exists, theprocedure outlined in§2canbeapplied toyieldtheFredholm integral equation givenby(2.8).Ifthekernel series given by(2.9) canbesummed inclosed form, ofsummed numerically toa specified accuracy byanappropriate seriessummation acceleration method [46],[47], [57],thentheintegral equation canusually beinverted numerically toobtain discrete solution values forg,(€). Cases where thekernel series cannot benumerically evaluated toaspecified accuracy represent anapparent limitation ofthemethod initspresent form. Forthese cases, alternate methods such asfinite differences, finite elements, or themethod ofweighted residuals, asmentioned earlier, must beused. SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS, 1081 Finite element methods, such asthose described byFinlayson [15], become necessary forproblems described onirregular geometry orforproblems having non- linear source terms [16], (32]. These situations, among others, represents another limitation ofintegral equation methods. Acomparison between the finite element method and integral equation method foralinearoperatorproblemofpracticalinterest inheterogeneous catalysis [30] hasshown that theintegral equation ismore efficient. Most oftheapplications cited previously, including theexample problems intro- duced inalater section, represent two-dimensional engineering type problems that have practical applications and aphysical significance. These problems, including others notcited here forbrevity, often involve regular geometry and linear operators having separation-of-variables solutions sothat theintegral equation method asout- lined here canbeapplied. Extension tothree dimensional problems ispossible also but has not been studied indetail. 4,Solutionoftheintegralequation.Theintegralequationgivenbyequation(2.8) hasbeen solved inthecurrent work byaquadrature method and found togive reliable results. Limitations ofthis method and potential problems associated with solving Fredholm integral equations ofthefirst kind arewell known (3),[19] and have been thesubject ofarecent symposium [49] sothat they will notbementioned here. Various quadrature rules, such astrapezoidal, Simpson's, and Gaussian quadrature, have been incorporated into ourintegral equation computer software forcomparison ofperform- ance, Numerical studies reported inthefirst example given in§5show that theerror between the exact and numerical values for the dual series coefficients decreases when more accurate quadrature rules areapplied. This generally agrees with results given byBaker (3).Afairlycommon featureisthatthekernelexpression oftencontainsalogarithmicsingularity atx=y.Either quadrature rules designed forintegrals with alogarithmic singularity canbeused asdescribed inDavis and Rabinowitz 8],ormodifiedquadrature [3],[19] canbeused. Forthelatter method, thefollowing system oflinear equations areobtained whose unknowns aretheintegral equation solution values gi(m) atthe quadrature abscissas &,: (4.1) YS.Brake=ln2=1,2,00°,N where (420) afe3tmKiy(EmEm)» =IM, WK EmEms n#m, (4.2b) n=SiEn)—hy(En) ra(4.2c) AlEn)=iKulm&')dé’. The quadrature indicated by(4.2c) represents anintegrated form ofthekernel. This hasbeen performed inclosed form fortheexample problems treated inthenext section and must betreated onaproblem-to-problem basis. Additional aspects related tothe properties thatthekernel Ky(é é')andsolution gx(£) should possess formodified quadrature toapply arediscussed intheabove cited references and areomitted here forbrevity. 1082 PLMILLSANDM.P,DUDUKOVIE 5.Numerical examples. Inthis section and §6,two example problems arenow considered thataretypical ofmixed boundary value problems inchemical engineering that lead todual orthogonal series. The basic objective isto illustrate theapplication ofboth themethod ofweighted residuals andtheintegral equation method tothese problems anddetermine some measure oftheir performance bycomparison ofnumeri-calresultsforafewkeyparameters ofinterest. Example 1.Heat conduction inapartially insulated plate. The first problem con- siders heat conduction inasemi-infinite plate. Thefaces atx*=/ andaty*->co are maintained ataprescribed temperature T=T.Thelower edge aty*=0ismaintained ataprescribed surface temperature T=T,for0x*<c*ands insulated forc*<x"5 1.The lineofsymmetry isatx*=0. Ifheat generation isabsent and theplate is isotropic, then setting c=0 in(2.1) gives theLaplace orpotential equation from which thesteady-state temperature profile intheplate canbeobtained. ‘Thegoverning equation andassociated boundary conditions aregiven below in terms ofthedimensionless temperature difference 6=(T—T)/(T,~ Ts), dimension- lesscoordinates x=y*/L, y=-7y*/, anddimensionless discontinuity ¢=mc*/: #0, #05.1 oohe0, 6) axay ae 2) =0, a0, yo, (5.2)x=0, 50, yE0, (53) x=, 020, yz0, (54a) O=1, 0sx<« =0 (5.40) 770) 00, ccxem ay (55) lim{0}=0, 0sxsm Asolution which satisfies (5.1) andtheboundary conditions given by(5.2), (5.3) and (5.5) is: 2 1 (5.6) (x,y)=Dbaei"?cos[(n+4):} Application ofthemixed boundary conditions (5.4a) and(5.4b) anddefining anew setofcoefficients a, =b,(n+1/2) leads tothefollowing dual series equations for determination ofthe series coefficients a4: (57a) F—*<c08|(n+4)x]=1, osx<*eont1/2 2)*J>e 6 2 P (5.76) Ea,cos[(n+4)x]=0 céxsm ‘These dual series contain one ofthesimplest forms fortheseries modifier encountered invarious applications and possess aclosed-form solution fortheseries coefficients. This example provides, therefore, animportant benchmark since adirect comparison canbemade between results produced bytheclosed-form solution and both theleast squares and integral equation methods SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS 1083 Thedual-series given by(5.7a) and(5.7b) canbereduced toasingle integral equation byfollowing theprocedure outlined in§2.The final result is (38) jfK(xyaa)ax=F, where thekernel K(x, »)isgiven bythefollowing infinite series which can bereadily summed using aresult given byGradshteyn and Ryzhik [17, eq.1.442-2]: _©cos{(n+1/2)x]cos[(n+1/2)9] Kisy)=2) n+l/2 (5.9) alinl=(x/2)+08{yf 2°" Joos(x/2)=cos (»/2) The kernel becomes singular asx+y, sothat determination ofg(x) bymodified quadrature requires evaluation of(4.2c). The final result is ‘ ce c+ (5.10)a=K(x,y)a(S 2)(Z)), o where the function fisgivenby[52,p.149,eq.16]: _&sin[(2n+1)z] Se)=3 2n+l ip (5.11) -4/In(tan3)dy weinz-24 §eZ ETpstmtgle 2 kokepE ™ The B,,, intheusual notation, denotes theBernoulli numbers oforder 2k. Stretched coordinates w=x/c and &=y/c areintroduced intheabove formulas forK(x,y)andA(y)sothattheintervalofintegration conveniently spans(0,1).Once theg(cw,) aredetermined, thea,areevaluated from thefollowing expression forthe series coefficients a,using Filon quadrature [8],[14]: sm a=2{'goa(ne!)s] ac Aclosed-form solution forthedual-series coefficients a,hasbeenobtained byTranterusingtwoalternate approaches [53],[54].Settingf(x)=Iandusingthesimpler approach [54] yields thefollowing closed-form expression forthedual series coefficients a4: P,(cos e) zr IKcatipy? Marte G13)*K(cos' (6/2) " where P,(cos c)andK(cos* c/2) refer totheLegendre polynomials oforder nand complete elliptic integral ofthefirst kind, respectively. Figure 1shows acomparison between the dimensionless temperature profile calculated byboth theexact solution, integral equation method, and theleast squares method forc=7/180 and y=mfor0x5 =.Little improvement intheleast squares approximation isobtained asthe number ofterms Nisincreased from N=5 to 1084 P.L.MILLSANDM.P.DUDUKOVIC exacrSOLUTION # i~_ F} * .i Ps a° \§ \ z *,RADIANS Fo. 1.Camprson between mensions emperaurerfcalculatedbytheeastsquaresmethodand exact solution forafixed value ofyandsmall¢(c=1/180). 'N=150. Even athigher values ofN,e.g., N=200, thedeviation isstill substantial. The same comparison ismade inFig. 2between theexact solution and theintegral equation method outlined above foradditional values ofy.The agreement between results was within 3-4 significant digits sothat thedifferences cannot bedetected on 12401 souron wo ee ° £: i i “oO we8 ~Ey a8 7 X,RADIANSFc.2.ComparitnBenendimensionlesstomerprofescalculatedbytheineralequationmethod and exc soir fortron eau @f)andsmallc(e= #180, SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS 1085 the figure. This represents asignificant improvement when compared tothe least squares results given ontheprevious figure. Theapproximate integral equation results were obtained using N=10Gaussian quadrature points sothat theamount ofcomputa- tional effort, when compared totheleast squares method atN=150, was significantly less and can bereadily performed onamicrocomputer. Table 2gives acomparison between numerical values forthedual series coefficients calculated bytheintegral equation method, theleast square method, and theexact solution forc=7/180. Forcomparison purposes, trapezoidal, Simpson's },and Gaussian quadrature rules with N=9 quadrature points were used fortheintegral equation method. The least square results, obtained using amatrix size ofN=150, still deviate significantly from theexact solution asindicated earlier inFig. 1.The integral equation results arealways ingood agreement with theexact solution, irrespec- tive ofthequadrature rule used. Further reduction intheerror between theexact solution and integral equation results isobtained forthemore accurate quadrature rules. Results obtained using theGaussian quadrature rule have theleast error among those tested which agrees with those reported byBaker [3].This indicates that itis nottheaccuracy ofthequadrature rule, butthenature oftheintegral equation approach that ismost significant. Tapte 2 Dual-series coeficients a,calculated byvarious methods for ¢=7/180. Integral equation method (N=9) Least squares Exact ‘trapezoidal Simpson’s}_—Gaussian method (N=150) solution °0.152369 0.163612 0.163058 0.330007x10"! 0.1631910.152346 0.163586 0.163034 0.201636%10-! 0.1631720.152301 0.163535 0.162985 0.17593910"! 0.1631230.152232 0.163458 0.162911 0.164900x10"! 0.1630440.152140 0.163385 0.162813 0.158740x10" 0.16294s0.152026 0.163227 0.162650 0.13479210"! 0.162826 0.151888 0.163073, 0.162543 0.152030 %10-! 0.16267 7 0.151728 0.162898 0.162371 0.149975 x10! 0.16250 8 o.1sisas 0.162690 o.162i75 0.148374 10-" 0.1623090.151339 0.162459 0.161954 0.187079x10" 0.16207 6.Example 2.Catalyst effectiveness factor inspherical ink-bottle pores. The second example isadopted from theinvestigation ofChuandChon [5]which dealt with thegeometry ofcatalyst pores and itseffect oncatalyst performance. The present case considers theproblem ofdiffusion and surface reaction inaspherical ink-bottle pore whose surface ismaintained atconcentration C,for05<a butundergoes first-order surface reaction r,=k,C, for a<0 7.The validty ofassuming spherical geometry versus one having anirregular shape represents anidealization, butexperimental data obtained byBroekhoff and DeBoer [4]were satisfactorily interpreted using this geometry asabasis. More recent applications involving zeolite orrelated cage-type Porous materials also support this assumption. With this supporting evidence, the governing equation and associated boundary conditions are: a( sau) 1 af. au) 1086 P.L.MILLS AND M. P.DUDUKOVIC (62a) p=0, pritao, ap (6.2b) u=l, 050<a, pat dusd, < (6.2¢)Daap t=O@<0S%, (6.24) o=0,7, 420, ospsi. 20 ‘The parameter Darefers totheDamkdhler number [2]and isgiven byDa=k,R/D, where k,isthereaction rate constant, Risthepellet radius, and D,istheeffective diffusivity. Solving (6.1) subject totheboundary condition given by(6.2a) and (6.24) gives thefollowing solution forthedimensionless concentration u=C/C,: (63) u(p,0)= 5a,p"P,(cos 8) where P,(cos 6)aretheLegendre polynomials. Application of(6.2b) and (6.2c) lead tothefollowing dual Legendre series: (64a) a,P,(cos 8)=1, 050<a, S n 6.40) a,(1+)P,(cos @)=0, a<O5m, (6.40) a-(+2)(cos)=0,a<0. These dual Legendre series can beconverted todual cosine series using Mehler's formula fortheLegendre polynomials (17, eq.8.913] followed byinversion ofthe resulting Abel integral equations. The final results are: ©,ata (1)J (x (65a) Ett cosnt5)x]=pge0s(5), OSx<a, 2 1 (6.5b)rb,cos[(n-+4)x], a<xSa, where thea,and 6,arerelated bythefollowing expression: _)Ltn/Da (66)by a Integration of(6.5a) over (0,x)andfollowingtheprocedureoutlinedin§2leadsto thefollowing Fredholm integral equation ofthefirst kind: (6.7) fK(xy)a(»)dy=f(x) where thekernel K(x, y)and forcing function f(x) aregiven by (68a) K(x,y)=S(Da, x-y)+5(Da,x+y) Qn, (x (6.8) Six)=Fsin(2). SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS 1087 The function Sintroduced inequation (6.8a) isgiven bythefollowing Fourier series: &sin[(n+1/2)z] 9) =5Salty] . Tolstov's method [52] can befollowed toobtain aformula forS(a,z) where the remainder isO(n~*) which issuitable fornumerical evaluation toaspecified accuracy. ‘The quantity ofinterest isthecatalyst effectiveness factor [2],[5]which isdefined astheratio oftheactual rate ofreaction occurring onthesurface oftheactive pore area totherate ofreaction that would exist ifthepore surface was maintained atthe surface concentration C,,The actual rate ofreaction can beevaluated byintegration ofthe flux orbyintegration oftheconcentration profile asshown below where the subscripts fand careused todistinguish between thetwo expressions: 1 *aw =-—1__|"sinoao 1 feo 6.10)=Daan 2,te!Pa(x)dx (6.10) Darter Bf ) Le, anti=Seay Lb(Pa i(c0sa)—cosaP,(cosa}, 21Fe0sa)EyPnDalPani(o0sa)—008aPa(c0sa] tot1“ats! u(1,6)sin60 tos poe 6.11 -_1_ Sa (6.1)Treo ofPa(x)dx bo, Da, _2ntt aS Pa cos Py |.3tagona)2Pataray Pasi(608a)£08aa(20sa) The results reported below arebased upon summation (6.11) using thefirst 100terms since numerical experiments showed that the calculated effectiveness factors were unaffected within 4significant digits when additional terms were used. Equation (6.11) ispreferred over equation (6.10) since theseries coefficients decrease bytheadditional factor ofn~'.The dual series coefficients b,were calculated from thecalculated values for g(y) byapplication ofFilon quadrature [14] tothe following equation which follows from (6.5b): 2[° r (6.12) bo=Z|,g(r)cos|(m+5)y| dy. Acomparison between effectiveness factor values calculated byboth theintegral equation method and themethod ofweighted residuals fora=7/180 and Da= 0.005 isgiven inFig.3.Fortheintegral equation method, thenumber ofquadrature points was varied between 2to18(ordinate values corresponding tothesecondary x-axis), while the matrix size Nwas varied between 5to150 forthe method ofweighted residuals (ordinate values corresponding totheprimary x-axis). Numerical values obtained atthesmallest and largest matrix sizes foreach method arealso given on thefigure forcomparison. The Galerkin method yields results substantially below those fortheother methods. The results fortheleast squares method decrease byonly 0.001inabsolutemagnitude asthematrixsizeNisincreased from5to150andoverestimatethe results produced bytheintegral equation and collocation methods. The results 1088 P.L.MILLS AND M.P,DUDUKOVIE 3a os § ‘COLLOCATION weri00, fi| |se 8OsINTEGRAL EQUATIONMETHOD 0878 B| fe ian byere _|= + MATRIX SIZE, N“ NUMBER OFQUADRATURE POWTS, N Fio.3.Comparison between effectiveness factor values calculated byvarious methods asafunction ofthe‘matrixsizeforanapertureanglea=7/180. produced bytheselattertwomethodsdifferbyonly0.056whichissmallwhencomparedtotheothers. Forthecollocation method, amatrix sizeN=137wasnecessary before acollocation point occurred intherange between 0to77/180 corresponding tothe range ofthefirstdual series given byequation (6.5a). Considering thatonly afew collocation points occur overtherange 0to77/180 atN=150,theordinary collocation method appears toconverge toward theintegral equation results. Values forthislatter method arerelatively constant forN=4 which represents asignificant reduction in thelinear system ofequations when compared toN=150.Chu andChon [5]used theGalerkin method andamatrix sizeofN=48 atthisaperture angle which, inview oftheabove results, explains their unreasonably lowvalues fortheeffectiveness factor inFig. |oftheir original publication. The results obtained byrepeating theabove exercise atalarger aperture angle (a= 1/18) aregiven inFig.4.Theagreement between effectiveness factors calculated bythevarious methods represents asignificant improvement when compared tothe results given inFig.3.Aninteresting observation which applies toboth Figs. 3and4 isthat effectiveness factors calculated bytheintegral equations, least squares, and collocation methods decrease with increasing values ofN,while those fortheGalerkin ‘methodincreasewithincreasing valuesofN.Thisseemstoprovideameansofassigning ‘anupper andlower bounds ontheexact value using theresults fortheintegral equation method and theGalerkin method atthelargest matrix size asthebasis. Forthecase where a=77/18, these bounds would be0.9078-0.9150. Theassessment oftheeffective- ness factor with thisaccuracy ismore than sufficient forchemical engineering applica-tions.Thevalidityofthe above bounds hasnotyetbeen rigorously established, however, and should beviewed asasemi-empirical result. Application ofthecollocation method asdescribed above uses collocation points6,intheinterval050zthatarerootsofthe series kernel, namely, Py..1(¢os @,)=0 fori=0,1,-++, N—1 where Nisthenumber ofcollocation points. This method for SOLUTION OFMIXED BOUNDARY VALUE PROBLEMS 1089 axelwTEGRALEQUATIONWETHODLS8.) See FH Icouoeanonwerioo aeTb.eseeewMATRIX SIZE |N maa NUMBER OFQUADRATURE POINTS, N Fic.4.Comparison berween efectveness factor values calculated byvarious methods asafunction ofthe aur sizeforanaperture angle =7/18 choosing thecollocation points becomes unsuitable when thediscontinuity angle approaches oneoftheinterval endpoints because thesmaller portion oftheentireintervalwillnotcontainasufficientnumberofabscissas toyieldasatisfactory approxi-mation,Theuseoffinite-element collocation [15]would divide theinterval (0,)into NE subintervals, with each subinterval containing N~2 interior collocation points. Initssimplest form fordual series, theinterval (0,)could bedivided into two subintervals, with thefirst one being defined over (0,a)and thesecond one over (a,77).This approach should givesuperior results when compared toordinary colloca- tionover theentire interval, butthenecessary development asapplied todual series equations liesoutside thescope ofthecurrent work. 7.Conclusions. An alternate method for solution ofdual and triple series equations hasbeen proposed which isbased upon reduction oftheseries toaFredholm integral equation ofthefirstkind. This equation issolved byanumerical quadrature method since closed-form inverses arenotreadily obtained oravailable formany problems involving practical applications. Theproposed method hasbeen used tosolve theexample problems selected from applications inheat conduction anddiffusion with reaction. Anexact solution was available forthefirst example and was used forcomparison totheleast squares and integral equation method. Itwasshown that theproposed integral equation method produced results that were ingood agreement with theexact solution, while those obtained bytheleast squares method hadsignificant error. Similar behavior wasfound foramoredifficultexamplewhereclosed-form solutionscannotbeobtainedbycurrentmethods, Forthisexample, results obtained byboth theintegral equation method and MWR were inbetter agreement forlarger values oftheboundary discontinuity which istypical forthese and related problems examined sofar. The results ofvarious numerical experiments suggest that theintegral equation‘methodrequiresthesolutionofamuchsmallersystemoflinearequations incomparisontoMWR toachieve thesame result. Itwas shown bythevarious examples that the 1090 P.L.MILLSANDM.P,DUDUKOVIC integral equation method yields infinite series forthekernel andforcing function whose translation toaclosed-form must betreated onaproblem-to-problem basis. Application ofMWRtoagivenproblem requires lesseffortbycomparison, especially once generalized computer software isprepared. REFERENCES [1]M.Abramowitz ANDL.StEGUN,eds,HandbookofMathematical Functions,Dover,NewYork,1970.[2]R.ARis,TheMathematical TheoryofDifusionandReactioninPermeableCatalysts,ClarendonPres,Oxford, 1975 [3]CT. H.BAKER, TheNumerical Treatment ofIntegral Equations, Clarendon Press, Oxford; 1977[4]TC.P.BROEKHOFF ANDT.H.DeBOER,Studiesonporesystemsincatalysts,J.Catal,10(1968),pp.153-165 [5]C.CHU AND K.CHON, Catalyst effectiveness factor ininkbotle pores,3.Catal,17(1970),pp.71-18 [6]W.D.Coctins,Somescalardifractonproblemsforasphericalcap,Arch.RatMech,Anal,10(1962), pp. 249-272 (71—~, Onsome triple series equations andtheirapplications, Arch. Rat.Mech. Anal, 11(1962), pp 12-137. [8]P.J.Davis AND P.RABINOWITZ, Methods ofNumerical Integration, Academic Press, NewYork, 1975. [9]M.P. Dupuxovie AND P.L,Mitts, Catalyst efectivenss factors intrckle-bed reactors, thInterna- tional Symposium onChemical Reaction Engineering, Adv. Chem. Ser. 65(1978), pp.387-399. [10] M.P.Dupuxovic, PL, Mitts ANDS.P.WALDRAM, Difusoninsinglecatalystpllesofspherical form,Paper133,TihInternational CongressofChemicalEngineering,ChemicalEquipmentDesign and Automation, CHISA 'I, Prague, Czechoslovakia, Aug. 31-Sept. 4,1981.[11]4.DUNDURSANDC.PANEK,HeatconductionBenweenbodieswithwavysurfaces,It.J.HeatMassTransfer, 19(1976), pp. 731-736. [12] R.P.FEINERMAN AND R.B.KELMAN, Theconvergence ofleast squares approximations fordual- drthogonal series, Glasgow Math J,18(1974), pp.82-84 3] RP. FEINERMAN, R.B.KELMAN ANDC.A.KOPER,JR,Dualorthogonalseries:Acasestudyof ‘heinfluence ofcomputing upon mathematical theory, Proc. Symposium Applied Mathematics, 20, ‘American Mathematica Society, Providence, RI,1974, pp.129-134 [14]L.N.G.FLON, Onaquadrature formula fortrigonometric integrals, Proc. Royal Soc. Edinburgh, 49 (1928/29), pp.38-47.[15]B.A.FINLAYSON, Non-linearAnalysisinChemicalEnginering,McGraw-Hill, NewYork,1980[16]S.Goro,A.LAKOTAANDJ.LevEC,Effectiveness factorsofnthorderkineticsintrckle-bedrectors,36(1981),pp.157-162. [17] LS. GRaDsHTEYN AND J.M. RYZHK, TablesofIntegrals,SeriesandProducts,AcademiePress,New York, 1965 [18]M.Herskowrrz, R.G.CARBONELL AND J.M.SMITH, Effectiveness factors andmass transfer in Irickle-bed reactors, AICHE 1,28(1979), pp.272-283.[19]LV.KaNToROVICH ANDV.KRYLOV,Approximate MethodsofHigherAnalysisthiededition,JohnWiley, New York, 1964. [20]R.B. KeLMAN, Thesteady temperatures incylinders withmixed Dirichlet andradiation conditions inthe lateral wallJ.Math.andPhys.46(1967),pp.333-342 [21]——,, Thepotential inacylindrical electrode partially submerged inaperfectly conducting medium, 4.Comput. Phys,2(1967),pp.120-128. [22]—— Harmonie mixed boundary-cae problems inthe plane, Q.Mech. Appl. Math, 22(1970), pp. ‘48-566,[23]——Separated-variables solutionforsteadytemperatures inrectangleswithbrokenboundarycondition,‘Trans. ASME, Ser. C,J.HeatTrans,98(1973),pp.130-132. [24] R.BLKELMAN AND C.A.KOPER, Jn,Leastsguares approximations fordual trigonometric series, Glasgow Math J,14(1973), pp.111-118. [25]RB. KeLwan AND J.T. SIMPSON, Analgrihm forsolvingadual-cosineseries,Comp.Math,with Appls,1(1978),pp.193-199. [26]RB,RELWAN,J.P.MAHEFFYANDJ.T.Sisson,Algorithmsforlscdualtrigonometric equations,Comp. Math. with Apps. 3(1977), pp.203-218.[27]R.B.Kewan,Least-squares Fourierseriessolutionstoboundaryvaleproblems,SIAMRev.21(1979), pp. 329-338 SOLUTION OF MIXED BOUNDARY VALUE PROBLEMS 1091 [28] RB. KeLMAn, Algorithm fortheclassic dual cosine equation, Appl. Math. Comp, 19(1980), pp.7-19. [29] ——, private communication, November 5,1981 [30] D.G.LorFLeR AND L.D.SCHMIDT, Catalytic activity and selectivity onheterogeneous surfaces with ‘mass transfer, AICHE J.,21(1975),pp.786-791. [31] M.LoweNnGrus, Apairofcoplanarcracksattheinterfaceoftwobondeddissimilarelastichalfplanes, Int.J.Engng. Sei, 13(1978), pp.731-741 [52] 0.M.MarTiNEZ, G.F,BARRETOANDN.O.LEMCOFF, Effectiveness factorofacatalystpelletin@ trickle-bed reactor, Chem. Engng, Sc., 36(1981), pp. 901-907. (33]B.N.MentaANDR.ARIS,AnoteonaformoftheEmden-Fowler equation,1.Math.Anal.Appl,36(1971), pp. 611-621 [34] W.W.Mever, LL. HEDGEDUS AND R.ARIS, Geometric correction factors fortheWeis: diffusivity cell, J.Catal, 42(1976), pp.135-138.[35]P.L.MILLSANDM.P.DUDUKowIG, Adual-sriessolution fortheeffectiveness factorofpartiallywettedcatalysts intickle-bed reactors, Ind. Engng. Chem., Fundam., 18(1979), pp.139-149. [36] —, Tripleseries equations: Solution bythemethod ofweighted residuals, Proc. 2nd InternationalConference onMathematical Modeling, St.Louis,July11-13,1979(X.J.R.Avula,ed.),Univ,Missouri-Rolla Press, 1980, pp.231-243.[37]—,Analysisofcatalysteffectiveness intrickle-bedreactorsprocessingvolatileornoncolatile reactants,‘Chem. Engng. Sci., 35(1980), pp.2267-2279. [38] ——, Application ofthemethod ofweighted residuals 1omixed boundary value problems, Dual-sries ‘relations, Chem. Engng. Sci, 35(1980), pp. 1557-1570. [39] —,, Application ofthemethod ofweighted residuals 10mixed boundary value problems. Triple-sries relations, Comput. Chem. Engng. 6(1982), pp.141-154[40]M.N.Oz1sik,Boundary-Value ProblemsofHeatConduction, International TextbookCo.,Scranton,PA, 1968.[41]K.S.PARIHAR, Sometripletrigonometrcal seriesandtheirapplication, Proc.Roy.Soc.Edin.Ser.A.,69(1971), pp. 255-265. [42] P.A.RAMACHANDRAN AND J.M.SMITH, Effectiveness factors intrickle-bed reactors, AICHE J.,2 (1979), pp. $38-$42.[43]S.S.SADHAL,Transientthermalresponseoftwosolidsincontactoveracirculardisk,Int.J.HeatMass ‘Transfer, 23(1980), pp.731-733.[44]——,Unsteadyheatflowbetweensolidswithpartiallycontactinginterface,1.HeatTransfer,Trans.ASME, 103(1981), pp.32-35. [43] G.E,SCHNEIDER,A.B,STRONGANDM.M.YOVANOVICH, Transientthermalresponseoftwobodies ‘communicatingthroughsmallcircularcontactarea,Int.J.HeatMassTransfer,20(1977),pp.301-308, [46]D.SHANKS,Non-linear transformations ofdivergentandslowlyconvergent sequences, 1.Math. Phys., 34.1955), pp.1-42.[47]D.A.SMITHANDW.F,FORD,Numericalcomparisons ofnon-linearconvergence accelerators, Math.Comp, 38(158), pp.481-499,[48]LN.SNEDDON, MixedBoundaryValueProblemsinPotentialTheory,North-Holland, Amsterdam,1966. [49] Society forIndustrial and Applied Mathematics 1981 Fall Meeting, Numerical Solution ofIntegral Equations Symposium, Cincinnati, Ohi. [50] H.TAKAG!, Slow motion ofasphericalcapinaviscousfluid.1.Uniformtranslationandrotation,J. Phys. Soc. Japan, 37(1974), pp. 229-236. [31] C.8.TaN AND J.M.SMrTH, Catalyst particle effetiveness with unsymmetrical boundary conditions, Chem. Engng. Sci. 35(1980), pp.1601-1609. [52] G.P.Torstov, Fourier-Series, R.A. Silverman, trans, Prentice-Hall, Englewood Clifs, NJ, 1962.[53]C.J.TRANTER, Dualtrigonometrical series,Proc.GlasgowMath.Assoc.4(1959),pp.49-57. [54] ——, Animproved method fordual trigonometrcal series, Proc. Glasgow Math. Assoc., 6(1964), pp. 136-140.[55]T.P.TsarANDW.S.FU,Aclassofasymmetric mixedboundaryvalueproblemsinthermoelasticity, IntJ.Engng. Sci, 13(1975), pp.343-352.[56]A.D.WHEELON, TablesofSummable SeriesandIntegralsIncolvingBesselFunctions,Holden-Day,San Francisco, 1968. [ST] P.WYNN, Onadevice forcomputing theén(S,) transformation, Math. Tables Aids Comp. 10(1956), pp.91-96.