85 1
Ind. Eng. Chem. Res. 1991,30,851-864 Panayiotou, C.; Vera, J. H. Statistical Thermodynamics of r-Mer Fluids and Their Mixtures. Polym. J . 1982a,14,681-694. Panayiotou, C.; Vera, J. H. An Improved Lattice-Fluid Equation of State for Pure Component Polymeric Fluids. Polym. Eng. Sci. 1982b,22,345-348. Panayiotou, C.; Vera, J. H. On the Fluid Lattice and Gibbs-DiMarzio Theories. J. Polym. Sci., Polym. Lett. Ed. 1984,22, 601-606. Prigogine, I.; Defay, R. Chemical Thermodynamics; Longmans: London, 1954. Quach, A.; Simha, R. Pressure-Volume-Temperature Properties and Transitions of Amorphous Polymers; Polystyrene and Poly(orthomethylstyrene). J. Appl. Phys. 1971,42,4592-4606. Quach, A.; Simha, R. Statistical Thermodynamics of the Glass Transition and the Glassy State of Polymers. J. Phys. Chem. 1972,76,416-421. Raucher, D.; Sefcik, M. D. Sorption and Transport in Glassy Polymers. Gas-Polymer-Matrix Model. Industrial Gas Separations; Whyte, T. E., Yon, C. M., Wagener, E. H., Eds.; ACS Symposium Series 223;American Chemical Society: Washington, DC, 1983; Chapter 6. Rehage, G. Thermodynamics of the Glassy State of Polymers. J. Macromol. Sci., Phys. 1980,B18,423-443. Roe, R.-J. Thermodynamics of the glassy state with multiple order parameters. J. Appl. Phys. 1977,48,4085-4091. Sada, E.; Kumazawa, H.; Yakushiji, H.; Bamba, Y.; Sakata, K.; Wang, S. Sorption and Diffusion of Gases in Glassy Polymers. Ind. Eng. Chem. Res. 1987,26,433-438. Sanchez, I. C.; Lacombe, R. H. Statistical Thermodynamics of Polymer Solutions. Macromolecules 1978,11, 1145-1156. Sand, M.L. Method for Impregnating a Thermoplastic Polymer. US. Patent 4,598,006,1986. Sanders, E. S. Penetrant-Induced Plasticization and Gas Permeation in Glassy Polymers. J. Membr. Sci. 1988,37, 63-80.
Sefcik, M. D. Dilation and Plasticization of Polystyrene by Carbon Dioxide. J . Polym. Sci., Polym. Phys. Ed. 1986,24, 957-971. Simha, R. Configurational Statistical Thermodynamics of Polymer Liquids. Ann. N.Y. Acad. Sci. 1976,279,2-14. Spalding, B. J. Foamed polymers toughen up their act. Chem. Week 1988,142,34. Staverman, A. J. Thermodynamic Aspects of the Glass-Rubber Transition. Rheol. Acta 1966,5, 283-292. Stookey, D. J.; Patton, C. J.; Malcom, G. L. Membranes Separate Gases Selectively. Chem. Eng. h o g . 1986,82,36-40. Vieth, W. R.; Alcalay, H. H.; Frabetti, A. J. Solution of Gases in Oriented Poly(ethy1ene Terephthalate). J. Appl. Polym. Sci. 1964,8,2125-2138. Vieth, W. R.; Tam, P. M.; Michaels, A. S. Dual Sorption Mechanisms in Glassy Polystyrene. J. Colloid Interface Sci. 1966,22,360-370. Vrentas, J. S.;Vrentas, C. M. Volumetric Behavior of Glassy Polymer-Penetrant Systems. Macromolecules 1989,22, 2264-2266. Wang, W. V.; Kramer, E. J.; Sachse, W. H. Effects of High-pressure C02 on the Glass Transition Temperature and Mechanical Properties of Polystyrene. J . Polym. Sci., Polym. Phys. Ed. 1982,20, 1371-1384. Wissinger, R. G. Thermodynamic Behavior of Glassy Polymer-Carbon Dioxide Systems at Elevated Pressures. Ph.D. Dissertation, University of Delaware, 1988. Wissinger, R. G.; Paulaitis, M. E. Swelling and Sorption in Polymer-C02 Mixtures at Elevated Pressures. J.Polym. Sci., Polym. Phys. Ed. 1987,25,2497-2510. Wissinger, R. G.; Paulaitis, M. E. Glass Transitions in Polymer-CO, Mixtures at Elevated Pressures. J. Polym. Sci., Polym. Phys. Ed. 1991,in press. Received for review June 20, 1990 Accepted December 21, 1990
Analysis of Anomalous Diffusion and Relaxation in Solid Polymers Nikolaos S. Kalospiros, Raffaella Ocone,?Gianni Astarita,* and Jerry H. Meldon' Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716-0001
A model of the diffusion in polymeric solids is presented, which incorporates a description of the intrinsic kinetics of swelling and the possibility of relaxation in the constitutive equation for mass flux. With the simplest possible constitutive assumptions, the final governing equations contain only two adjustable parameters. Nonetheless, the equations can account for all anomalous experimental observations reported in the literature on mass transfer in polymers. In the Appendix, the architecture of the numerical procedure used to solve the governing equations is outlined. The procedure is discussed in general terms because it is in principle applicable to quasilinear hyperbolic problems arising in a variety of contexts including gas dynamics and adsorption.
Introduction Diffusion of low molecular weight solutes in glassy polymers is known to exhibit peculiarities that cannot be adequately described by any formulation of classical Fickian theory of diffusion. In particular, there are four major types of experimental results that deserve attention: sample size effects (SSEs),case 2 transport (C2T),twostage sorption (2SS), and sorption overshoot (SO). (1)SSEs manifest themselves in the nonoverlap of plots of W / W ( - ) , where W is the weight sorbed and W ( w )is its equilibrium value, versus t / X 2 ,where t is the time and X is the sample half-thickness [contrary to a well-known theorem of classical Fickian diffusion theory (Crank,
* T o whom correspondence should be addressed. Permanent address: Dipartimento di Ingegneria Chimica, Universitl d i Napoli, 80125 Naples, Italy. t Permanent address: Chemical Engineering Department, Princeton University, Princeton, N J 08540. t Permanent address: Chemical Engineering Department, Tufts University, Medford, MA 02155. 0888-5885/91/2630-0851$02.50/0
1975)], and in fact even in qualitative differences among different sample sizes (Enscore et al., 1977); models describing this have been presented (e.g., Astarita and Joshi, 1978). (2)C2T manifests itself in that dW/dt is initially finite, again contrary to classical diffusion theory [in fact W often ; models in the varies linearly with time up to W ( w ) ]again, literature describe C2T (Windle, 1985; Astarita and Sarti, 1978). (3) 2SS is sometimes observed in sequential experiments, wherein the gas-phase fugacity of the solute is increased to a series of constant values; in each step, W grows until an apparent equilibrium is reached, only to subsequently grow at a much smaller rate to a significantly larger final equilibrium value (Bagley and Long, 1955). Again models have been presented that adequately describe 2SS (Joshi and Astarita, 1979; Long and Richmann, 1960). (4) The most striking experimental observation is "overshoot": the amount of solute sorbed increases beyond its equilibrium value, reaches a maximum, and then decreases toward its final equilibrium value. In some in@ 1991 American Chemical Society
852 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991
stances (Overbergh et al., 1975; Kambour et al., 1966; Waywood and Durning, 1987; Ware et al., 1981; Mercier et al., 1967; Titow et al., 1974), the effect has been attributed to crystallization: the solute reaches equilibrium with swollen polymer prior to crystallization, but when crystals form they “squeeze out” the solute. However, in well-documented cases (Peppas and Urdahl, 1986,1988; Vrentas et al., 1984; Faulkner et al., 1977; Smith and Peppas, 19851, overshoot is observed in the absence of crystallization. While some of these results may be attributable to morphological changes other than crystallization (Durning, 1989), the data of Vrentas et al. cannot be so explained, since the authors observed overshoot in a repeat experiment on the very same sample (Duda, 1989). Furthermore, the data of Baird et al. (1971), though exhibiting only minor overshoot, defy explanation by any kind of morphological change. Apparently, no currently available model predicts the SO phenomenon in the absence of crystallization or some other, irreversible morphological change. Available models are all ad hoc, and while they may describe one or two anomalous behaviors, none account for all of them. A truly robust model would elucidate the relationship between all observed behaviors and between values of relevant parameters (Lustig et al., 1989). Among published ad hoc models, that proposed by Astarita and Sarti (1978) can be used to correlate both C2T and SSE. However, its success rests on the assumed existence of a discontinuity of the state of the polymer (a glassy-swollen interface) and, under nonequilibrium conditions, of a discontinuity of the solute chemical potential a t that interface. Consequently, a robust model should p r e d i c t t h a t , under appropriate conditions, discontinuities-or perhaps quasidiscontinuities-of the polymer state and chemical potential may develop. The analogous heat-transfer problem has been analyzed by Astarita and Kenny (1987) and, in more general form, by Ocone and Astarita (1987). Notably, the latter formulation-which involves a viscoelastic term in the constitutive equation for heat conduction and finite kinetics of morphological change-allows for temperature discontinuities. However, its extension to mass transfer is not trivial, since a fundamental nonlinearity is thereby introduced that has no counterpart in the heat-transfer case. In this paper, we present such an extension and show that the resulting equations are indeed capable of describing all observed behaviors, including SO. In the next section, we detail the physical rationale for the new model, since the resulting governing equations are rather unusual. For constitutive functions, we use Occam’s razor and choose the simplest possible form, so that the model is tested severely: its final formulation contains only two adjustable parameters. We focus attention on a polymer slab of finite thickness 2 X , exposed at both surfaces from time zero onward to a gas containing a soluble component at fixed fugacity, not in equilibrium with the initial polymer. Most published data on anomalous sorption behavior have been obtained under such conditions. As will be seen, the practical finite slab case requires much more mathematical structure than that of the seminfinite body considered by Astarita and Meldon (1989), since traveling discontinuities are reflected back and forth between the midplane and exposed surfaces of the finite slab. Since the model’s governing equations are hyperbolic, their numerical solution presents substantial difficulties, particularly when discontinuities are reflected back and forth between the two bounds (Sod, 1978). We have de-
veloped a numerical procedure that essentially eliminates this difficulty and present its architecture in the Appendix.
Physical Model To formulate the problem, three questions need to be addressed. First, since the focus is upon mass transfer, the variables whose distributions in space (and possibly time) determine the rate of mass transfer must be identified. Second, since the polymer may exist in different states (at the very least a glassy and a swollen state), a measure of the local state of the polymer is necessary. Finally, constitutive equations that govern respectively the rate of evolution of polymer state and the rate of mass transfer must be established. The first point is perhaps the easiest to discuss, since the chemical potentials of the two species (polymer and solute) are undoubtedly the relevant variables. As long as the system is dilute (say the mass fraction of the solute is always > 1, one recovers 2SS, as discussed earlier. If ER m, swelling could never be observed experimentally, and the governing equations would reduce to Fick’s law. A t the other extreme (with 4 still zero), when tR 0, one obtains a modified form of the corresponding heattransfer analysis (Astarita and Kenny, 1987), which in this limit predicts a discontinuity of state (a glassy-swollen interface) in the absence, however, of a chemical potential discontinuity across the same interface. In the heattransfer case, this represents the classical Stefan (1891) problem of the freezing of water. As already acknowledged, our model cannot be expected to produce sensible results in the limit R a. Considering finite R only, the case in which t 0 is uninteresting, because swelling would simply not occur. When 1 0, U 2 < P,which is sometimes (perhaps improperly) called the “entropy condition”; i.e., the secondary discontinuity lags behind the primary one. Its speed of propagation will in general not be constant, because the factor (m)/(h*P]will depend on the chemical potential to the left of the secondary discontinuity. [Note that this was assumed a priori in the successful ad hoc model of Astarita and Sarti (1978). Hence, in the limit 0 = 0, one reproduces that model, with the discontinuity’s speed of propagation varying with the local chemical potential in a manner determined by the
form of h*(P). Note, too, that a variable speed of propagation can also be achieved, without imposing a variable degree of swelling, by simply allowing the relaxation time 4 to depend on chemical potential (Sarti, 1989).] The following summarizes our conclusions about the possible existence of discontinuities: When 4 = 0, no discontinuity of chemical potential can develop (Astarita and Ocone, 1989). However, a discontinuity of state may arise when tR = O/T 0 (Astarita and Kenny, 1987). With finite 4, a primary discontinuity of chemical potential arises, but the state is continuous across this discontinuity. However, as t 0, the primary discontinuity becomes vanishingly small, and a secondary one develops across which both state and chemical potential jump. Again considering t 0, if one expresses h*(/3) by the power law h* = PK,the speed of propagation will always be constant, U 2 = R/(K + l),which is a consequence of the dimensional peculiarity of power law constitutive equations (Astarita, 1985). In the particular case considered here, K = 1 and U 2 = R/2. For purposes of further analysis, it is useful to define the quantity q = P2, so that the governing equations and jump conditions (eqs 41,42,44 and 46) reduce to an/& + a(q/2)/ay = -n (47) R anlay + aq/ar = o (48)
-
-
-
141 = 2 m 1
(49)
(50) 2U2 = R The secondary discontinuity is, of course, also reflected whenever it reaches either the midplane or the exposed surface, and a numerical procedure (see Appendix) can again be used to solve the asymptotic case e = 0, using the jump conditions to couple the two domains of integration. The differential equations for (q)and (n)are particularly simple for the secondary discontinuity: 2dlql/dr = -I41 (51) and an entirely analogous equation for (n). Both jumps decay exponentially in time, except at the reflection points when jumps are imposed on them.
Numerical Results Details of the numerical procedure implemented to solve the model’s governing equations are given in the Appendix; in this section, we present numerical results and, when appropriate, compare them with experimental evidence. All the numerical results have been obtained for r = 2. In particular, we present plots of dimensionless weight uptake w vs time. Defined as the actual dimensional weight uptake W divided by X p a , , w is calculated from either of the following integrals: w = &l(ph - 1) dy’ = R&‘n(O,
7’)
dr’
(52)
We examine first the “balanced” case where R = t = 1; plots of n, P, and h at three times are presented in Figure 1. The flux curves in Figure l a show cleady the traveling primary discontinuity, which is located at y = r before the first reflection, at y = 1 - ( 7 - 1) between the first and second reflections, and at y = r - 2 between the second and the third ones. The amplitude steadily decreases, except at reflections at the exposed surface, where the discontinuity actually changes sign. Conversely, as shown in Figure lb, the potential discontinuity changes sign at reflections on the midplane. As expected, the distributions shown in Figure ICdo not exhibit discontinuities, although
Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 857 I
1.2 I
2 T
=OB
..
0.0
0.2
0.6
0.4
0.8
1.0
0.0
0.2
0.6
0.4
0.8
1.0
Y
Y
Figure 2. R = 1, distributions of n values. Solid curves, complete solution with c = 0.001; dashed curves, asymptotic solution for c = 0. 5
0.0
0.2
126
0.4
0.6
-
0.8
1.0
Y
h
0.0
0.2
0.4
0.6
the gradient of h is discontinuous. While the latter discontinuity is clear at 7 = 0.8, it becomes very small after the first reflection, as the two upper curves indicate. Next consider the asymptotic case where t 0, which is the most interesting one. Figure 2 gives the flux distributions at three different times as calculated from the general set of equations for R = 1 and t = 0.001 (solid curves) and from the asymptotic analysis for R = 1 and t = 0 (broken curves). The solid curves show clearly that the primary discontinuity has become vanishingly small (e.g., the curve for T = 0.8 seems to approach the abscissa axis asymptotically, while there effectively remains a vanishingly small dicontinuity at y = 0.8) and that a secondary quasidiscontinuity lags behind the primary one. Except for the comparatively narrow region where the solid curves exhibit a marked sigmoid shape, the solid and broken curves in Figure 2 clearly coincide. This indicates that the sigmoid part of the solid curves is a moving boundary layer of the type discussed by Astarita and Kenny (1987), which is sandwiched between “inner” and “outer” solutions obtained from the asymptotic analysis. Note also that the inflection points in the solid curves are located where the asymptotic analysis pinpoints the discontinuity. Similar considerations apply to the distributions of potential and state presented in Figure 3; here, notably, even the solid curves for /3 and h are undistinguishable, indicating that the approximation in eq 40 is justified. Finally, Figure 4 gives the flux at the exposed surface as a function of time. Solid and broken curves coincide except in a small region near the discontinuity. Note that at time zero the solid curve takes the value of 1, which is quite different from the asymptotic value. This is not surprising, since the asymptotic analysis is not valid in the vicinity of T < E; however, the solid and broken curves essentially coincide at all T > e. Figure 5 gives the dimensionless weight uptake vs time for the original “balanced” case R = 6 = 1. Notably, the curve has the typical shape associated with C2T,with the initial slope being finite and in fact w(7) being approximately linear up to about 80% of saturation. Since R = 1, 7 is the ratio of the actual time to the diffusion time. Saturation is reached at T of order 4; i.e., sorption is slower than would be calculated on the basis of Fickian diffusion, consistent with C2T.
0.8
1 .O
Y Figure 1. Balanced case, R = 1,c = 1. (a, top) Flux n; (b, middle) fugacity 6;(c, bottom) swelling h.
858 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 2.2
4
2.c 3 1.e
P,h
w1z)
1.e
2
I
\J
7=0.8 1.4
I
1
1.2
1.c 0.0
C 0.2
0.4
0.8
0.6
5
1.0
Y Figure 3. R = 1. Distributions of and h (which are undistinguishable). Solid curves, complete solution with c = 0.001;dashed curves, asymptotic solution for c = 0.
0
1
2
3
4
10
15
0
z Figure 5. Dimensionless weight uptake vs dimensionless time, R = c = 1. Case 2 transport is predicted; the initial slope is finite, and w ( i ) is linear up to T = 2.
O
I 0 100 200 300 400 500 600 700
z Figure 4. Flux at exposed surface vs time. Continuous curve, complete solution for c = 0.001;the isolated point at T = 0 belongs to this curve. Dashed curve, calculated from asymptotic analysis for c = 0.
Figure 6. w(7) for R = 0.01, c = 1. Pseudo-Fickian sorption is predicted; the initial slope is infinite.
Figure 6 gives the dimensionless weight uptake vs time for the case R = 0.01 and t = 1. Here, behavior is quite similar to that of ordinary diffusion: the initial slope is infinity, and equilibration is reached at 7 of order 200, which correspond to =2 diffusion times. We now come to Figure 7, which presents the results for cases where swelling is slow; i.e., t >> 1. Saturation is expected to be approached when dimensional time is of order 8, corresponding to 7 = O(t). Two-stage sorption is anticipated when swelling time far exceeds diffusion time; Le., tR >> 1. Figure 7a shows one such case of 2SS (tR = 25): a (small) region of positive curvature is observed, due to interaction of the tendencies of sorption overshoot ( R = 1) and slow post-diffusion swelling. Saturation is indeed approached at 7 of order t, Figure 7b shows a case where, since eR = 10, 2SS is only approximated. Here, an initial fast uptake into the unswollen polymer takes place over a time period of order T (corresponding to 7 = 1/R = IO), followed by a slow approach to saturation, again at 7 of order e . Finally, Figure 7c shows a more typical 2SS case:
an initial rapid approach to saturation of the unswollen polymer, followed by a very slow approach to final saturation. Figure 8 gives the dimensionless weight uptake for the most interesting case of R = 1 and t = 0.001; again, the dashed curve is obtained from the asymptotic analysis for R = 1 and t = 0. Two points are of importance. First, regarding weight uptake, the asymptotic analysis can be applied with confidence, since the two curves essentially overlap. Second, and most important, sorption overshoot is predicted. It is interesting to note that, of the model’s two parameters, R depends on the sample thicknes but t does not. Hence, sample size effects must be explained by considering changes in R. This is indeed possible, since transition from case 2 transport to apparent Fickian diffusion, which has been observed experimentally as a sample size effect (Astarita and Joshi, 1978), is predicted when R changes, and indeed in the right direction (apparently Fickian diffusion in larger samples: compare Figures 5 and 6). On
Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 859
.....................
3 ..-....I....-----.............I
"
,
0
.
.
50
.
,
.
,
.
150
100
,
200
5
0
15
10
z
z
Figure 8. W ( T ) for R = 1. Full curve, complete solution for c = 0.001; dashed curve, asymptotic analysis for c = 0. Sorption overshoot is predicted.
3.0 0
1
2
3
4
5
6
7
8
9
10
R Figure 9. Maximum value of w vs R for t = 0.001. Any value in excess of 3 is an overshoot.
3
1 0
.. ...
.........................
100
200
300
...............
400
z Figure 7. Two-stage sorption behavior. (a, top) R = 0.25,t = 100. (b, middle) R = 0.1, t = 100. (c, bottom) R = 0.04,t = 1000.
the other hand, transition from C2T to 2SS behavior, which requires a change in e, has never been reported in the literature. It is also important to notice that, all else being equal, the model equations predict that the extent of sorption overshoot increases with R, see Figure 9. An increase in R corresponds to a decrease in sample thickness, and indeed the experimental data of Vrentas et al. (1984) show that the amount of overshoot increases with decreasing sample thickness. The latter data significantly represent the best documented cases of overshoot, which cannot be attributed to either crystallization or some other, irreversible structural change. Conclusions A model containing only two adjustable parameters has been presented, which is capable of describing qualitatively all reported experimental observations of anomalous diffusion phenomena in solid polymers. While agreement with experimental data is not quantitative, this is presumably due to the assumption of oversimplified constitutive equations with constant parameters.
860 Ind. Eng. Chem. Res., Vol. 30, No. 5,1991
T = diffusion time, s : t = time, s
t’ = dummy time variable, s U = speed of secondary discontinuity u = any variable V = speed of primary discontinuity V * = value of V at reflection instant W = weight sorbed, kmol w = dimensionless weight sorbed X = sample half-thickness, m x = distance from exposed surface, m y = dimensionless distance from exposed surface Y = thickness of ED z = dimensionless distance from discontinuity a = I / & - 1) /3 = dimensionless fugacity p’() = defined in eq A8 p’() = defined in eq A9 c =
1
0
Y Figure 10. Strip over which the solution is sought, and the two domains in which it is divided by ~ ( 7 ) .
Nomenclature A = constant in linearized /3 distribution B = coefficient in linearized /3 distribution b = coefficient in linearized Q distribution C = coefficient in linearized /3 distribution c = concentration, kmol/m3 c* = equilibrium concentration, kmol/m3 C() = defined by eqs A4 and A5 C*O = C(hb,r),/3b,7))= C*(Y,d d = coefficient in linearized n distribution
D = diffusivity, m2/s f = q / p , see eq A32 F ( ) = function in eq 3, kmol/(m3 atm s) g = value of h at primary discontinuity
H = c / p , kmol/(m3atm)
Ho= value of H at time zero, kmol/(m3 atm) H* = equilibrium value of H, kmol/(m3 atm) h = dimensionless local state h* = equilibrium value of h h’ ( ) = defined in eq A10 k ( ) = see eq A2 k* = value of k() at reflection instant L = mobility, kmol/(m atm s) L’ = lower limit mobility, kmol/(m atm s) 1 = dimensionless mobility m = defined in eq 43 m ( ) = see eq A1 m* = value of m( ) at reflection instant M = number of grid points along y N = flux, kmol/(m2 s) n = dimensionless flux n’() = defined in eq A6 n ” ( ) = defined in eq A7 n* = value of n at reflection instant p = fugacity, atm po = initial value of fugacity, atm p ( ) = see eq A1 P = defined in eq A21 PC = value of P at reflection instant 9 = P2
q ( ) = see eq A 2 Q = defined in eq A22 Q* = value of Q at reflection instant R = d/T r = imposed fugacity ratio r ( ) = see eq A3 s = instantaneous position of discontinuity
0/+
el, t2 = discrepancies 8 = swelling time, s p T T*
= kernel in eq 6, kmol/(m s atm) = dimensionless time = time of first reflection
= time of Nth reflection C#J = relaxation time, s 4’ = retardation time, s T~
Subscripts L = to the left of the discontinuity R = to the right of the discontinuity Operators [ ] = jump across primary discontinuity { } = jump across secondary discontinuity Acronyms C2T = case 2 transport ED = expanding domain PDP = primary discontinuity problem SD = shrinking domain SDP = secondary discontinuity problem SO = sorption overshoot SSE = sample size effects 2SS = two-stage sorption Appendix In the body of this paper, a hyperbolic set of governing equations was obtained. Here we discuss the numerical procedure that we have developed to deal with this type of problem. The procedure is believed to be of rather wide applicability, in particular to problems of the type arising in gas dynamics, adsorption, and relaxation phenomena. We present it in some generality, concentrating on the formulation that arises in relaxation phenomena. We limit the presentation to the architecture of this procedure, omitting discussion of the numerical details, since the latter are quite standard. The simplest original formulation was developed by Ocone (1987) for a heat-transfer problem, in which context the equations are essentially linear. The formulation was discussed concisely by Ocone and Astarita (1987). Subsequent extension to a mass-transfer problem is given by Astarita and Meldon (1989);in this, significantnonlinearity is introduced. Both Ocone (1987) and Astarita and Meldon (1989) considered only a semiinfinite slab, where there is only one traveling wave to account for and no reflections; here we extend the analysis to the more realistic case of a finite slab. We also generalize the equations (although we do not claim to have cast them in the most general form, still patient of numerical solution along the con-
Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 861 ceptual lines presented) to include rather strong nonlinearities. In particular, we allow for nonconstant speed of propagation of discontinuities, as opposed to the simple cases considered in the body of this paper. We consider the following quasilinear hyperbolic set of differential equations: w a y + P ( ~ , Pw/a7 ) + W,P) = o (AI) a n l a 7 + s ( h , ~ap/ay )
+ Mn,P) = o
ah/a7 = r ( h , ~ )
(A2) (A3)
Functions p and q are positive valued, and r is such that r = 0 has a unique solution h = h*(P),r > 0 when h < h*(P),and r < 0 when h > h*(@).It is helpful to think of y as position, 7 as time, n as flux, P as potential, and h as a yardstick of the local state; eq A1 is then seen to be a balance equation, eq A2 a constitutive equation for the flux which includes relaxation (without which the set would not be hyperbolic), and eq A3 a kinetic equation for the temporal evolution of the state. Now consider the function C(h,@), which, to within an arbitrary additive constant, is defined by
ac/ao = P ( ~ , P )
the constitutive equation (when this is feasible), one can establish both the discontinuity’s speed of propagation and the coupling conditions across it which couple the solutions in the two domains, thus transforming the problem to a well-formulated free boundary one. The primary discontinuity problem (PDP) is represented by eqs 23-25. This can be seen to be a special case of the set (Al)-(A3), where the functions of interest (omitting the additional assumption that l(h) = h) are m(h,P) = Pr(h,P)/R (All) p(h,P) = h / R
(-412)
q(h,P) = l(h)
(A131 k(n,P) = n (A14) The secondary discontinuity problem (SDP) is an asymptotic limit of the first one and is described by eqs 41 and 42. This is again a special case of set (Al)-(A3), with the following choice of functions (where h*’ is the ordinary derivative of h*): P(h,P) = (h*(P) + Ph*’(P))/R (AIS)
=0 (A161 dh,P) = l(h*(P)) (A171 aC/ah = m(h,P)/r(h,P) (A5) Mn,P) = n (Ala) This definition reduces eq A1 to anlay dC/& = 0, and Note that, in the SDP, eq A3 plays no role and the hence the value of C may be identified with concentration. variable h drops out. Function C relates concentration to the potential and local The Kotchine procedure is essentially as follows (Asstate. tarita and Ocone, 1989): Let y = s ( ~be ) the position of The boundary conditions imposed are of the following V = ds/dr be its speed of propagation. a discontinuity and type: The equations are first transformed with the new pair of (A61 nb,O) = n’b) independent variables z = y - s and 7. One then integrates each term in the equations over -6 < z < e and finally takes n(1,7) = n”(7) (A7) the limit of the resulting expression as e 0. P(Y,O) = p’b) (A81 If u denotes any dependent variable that may appear in the equations, the Kotchine procedure yields the folP(0,7) = V(7) (A9) lowing results: hb,O) = h’b) (A101 (i) when u itself appears, it transforms to zero provided u is bounded everywhere, as we assume to be the case for where the terms on the right are known, smooth functions the three functions m, k, and r. of their arguments. (ii) When aulay appears, it transforms to the jump of A jump is imposed on the system provided that P’(0) is u over the discontinuity, which will be denoted by [ u ] ;if different from j3”(0) or n’(1) is different from n”(0) or both. uL and uR are the limits of u as the discontinuity is apAlthough the procedure discussed below can be applied proached from the left and from the right, [u]= uL - u p whatever the forcing functions, the formulation in the body (iii) When aula7 appears, it transforms to -V[u]. of this paper is restricted to the simple case where n’ = We first apply the procedure to eq A3, which immedin” = 0,P’ = h’ = 1 (sample initially at equilibrium, and ately yields an impervious midplane), and pl’ = r, so that only the simple constant parameter r appears in the boundary [h]= 0 (AH) conditions, and a jump of 0is imposed at time zero. The This means that the state variable h suffers no disconsolution is sought over the semiinfinite strip 0 I y I 1 , ~ tinuities, and hence, one may define, without ambiguity, > 0. the quantity g as The essential basis of our numerical procedure is implementation of the Kotchine (1926) procedure to transg = h(S(7),7) (A20) form the problem, as formulated above, to a variation of The second terms in eqs A1 and A2 are not amenable a free boundary problem (Crank, 1984; Fasano and Primto the Kotchine procedure as they now stand. However, icerio, 1983). This is done because numerical integration eq A18 makes possible the transformation to a more useful of hyperbolic sets is complicated by potential discontinuform by defining the quantities P and Q as follows (as will ities in the domain of integration (Sod, 1978). Transforbe seen, the lower limit K is arbitrary and plays no role mation to a free boundary problem obviates this difficulty, in the subsequent analysis): since integration is then carried out separately in domains bounded by the traveling discontinuity. Thus, the target Pk,P) = L p @ , p l )do’ (A21) solution is smooth within such domains. In ita classical formulation, the Kotchine theorem (1926) is applied to the balance equation only and thereby would Qk,P)= l’qk,pl) K dF (A22) only produce a relationship between the jumps of the flux and of concentration across the discontinuity. By exSince the Kotchine procedure only implies consideration tending the procedure underlying Kotchine’s theorem to of a neighborhood about the discontinuity, h may be taken (A4)
+
-
862 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991
to be g throughout this procedure, so that the second terms in eqs A1 and A2 become derivatives of P and Q. One may now apply Kotchine’s procedure to eqs A1 and A2 to obtain [nl = V[PI (A23) V[nI = [QI (~24) V2 = [Ql/[Pl Equation A23 is the jump balance condition emerging from the classical Kotchine (1926) theorem, which is restricted to balance equations (indeed, the jump [PI is the jump of concentration, since eqs A5 and A19 imply that C may jump only because P itself jumps). Equation A24 is of the type referred to by Astarita and Ocone (1989) as a jump constitutive equation. It can be obtained because the constitutive equation is only mildly nonlinear. Equation A25 gives the absolute value of the speed of propagation of the discontinuity, the sign being then determined by either eq A23 or eq A24. By providing both the speed of propagation and a coupling condition between the solutions on both sides of the discontinuity, eqs A23-A25 provide the basis for transformation to a free boundary problem. For the PDP and the SDP, the quantities [PI and [Q] are given by PDP [PI = g[PI/R (-426)
181 = W [ P 1
(A271
[PI = [Ph*(P)I/R
(-428)
[QI = JBLNh*(i3’))dD’
(A29)
SDP
BR
With the additional assumption that l(h)= h, the PDP produces a constant speed of propagation, V2 = R. If, as we have done, one also assumes that h*(P)= P, the SDP also produces a constant speed of propagation, V2 = R / 2 . Indeed, a constant speed is obtained for a milder assumption for the SDP, namely, that l(h*(P))= PK, which yields V2 = R/(K + 1). The strip over which the solution is sought is sketched in Figure 10. The position of the traveling discontinuities S ( T ) is also sketched; the S ( T ) curve is, of course, not known in advance (except in those special cases when the speed of propagation is constant). Let T ~ N, = 1, 2, ...,be the times at which s is either unity or zero (see Figure 10); these are the times at which reflections of the traveling wave take place. It is useful to separately consider the respective distributions of n, 0, and h between two consecutive reflections; PN(y,7) will denote the distribution of P for T~ < T < T ~ and + ~ similarly for n and h. Over each interval of time, there is both a shrinking domain (SD) and an expanding domain (ED): when N is odd, the SD is on the right and the ED on the left and vice versa when N is even. At the end of an odd-N interval, a reflection occurs at y = 1. This is due to the fact that nN(l, T~ - 0) will in general differ from n”(TN),and hence, a jump of n will be imposed by boundary condition (BC) (A7). Conversely, at the end of an even-N interval, a reflection occurs at y = 0, because PN(O, TN - 0) will in general differ from p 1 ) ( ~ ~ ) , and hence BC (A9) will impose a jump of P. Given the jump of n imposed at an odd-N reflection, the corresponding P jump can be calculated from eq (A23). This will also yield the jump of p, since p > 0 and hence eq (A21) is invertible for p. Conversely, given the jump of
imposed at an even-N reflection, the corresponding P jump may be calculated from eq (A21) and the corresponding n jump from eq (A23). We now focus on some time T other than one of the reflection instants. Suppose the spatial distributions of all quantities of interest, as well as the position of the discontinuity, are known at this time. It is then possible to calculate the instantaneous speed of propagation from eq (A25), with V > 0 for odd-N intervals and V C 0 for even-N intervals. Let 67 be the time step chosen. One may now integrate forward in time (with any of the many known techniques for doing so) to obtain first the new position S ( T + 67) and then the new distributions. This is done separately for the SD and the ED and is discussed in some detail below for an odd-N interval (the case of an even-N interval is entirely analogous). In the SD, forward integration yields directly the values at the SD side of the discontinuity, e.g., for an odd-N interval, the values to the right of the discontinuity. In the ED, forward integration yields the new distributions only up to y = 4 7 ) . Hence, tentative values of Q, P, and n to the left of the discontinuity, identified by an asterisk, are first obtained by extrapolation. Since these values will in general not satisfy eqs (A23) and (A24) exactly, one has v(PL* - P R ) = nL* - nR + (A30) (A30 QL* - QR = V(nL* - nR) + Adjustments are thus needed in the values of QL, PL,and nL, and eqs A30 and A31 are clearly insufficient to determine them. However, 6QL and 6PLare not independent of one another, since eqs A21 and A22 imply that ~ Q L= qk,PL)8PL/Pk,PL) = f 6 P L (A32) Equation A32 furnishes the additional relationship that unequivocally determines the adjustments in terms of el and c2: 6nL = 6PL =
w e 2 (62
- f e l ) / ( V 2- f )
(A33)
V€,)/(V2 - f )
(A34)
-
The above equations identify completely the adjusted values of nL and PL, except in the case where the speed of propagation is constant (which is in fact the case in both the PDP and the SDP). If V2 were indeed constant at some value R, eq A25 would imply that [Q] = R[P], which can be satisfied for arbitrary jumps of P and Q (see eqs A21 and A22) only if q(g,P) = Rp(g,P). However, the latter result implies that the right-hand sides of eqs A33 and A34 reduce to the indeterminate form 010. Hence, the procedure sketched above cannot be applied to those cases where V2 is constant. However, for the constant V2case, a different procedure can be adopted. Differentiation of eq A23 with respect to time yields d[n]/ds = V d[P]/dr (A35) and the two sides of eq A35 may be expanded as d [ n ] / d ~= [ d n / d ~ ]+ V[an/ay] (A36) V d [ P ] / d ~= V[aP/ar]
+ V2[dP/dy]
(A37)
If one now applies the [ ] operator to all terms in eqs A1 and A2 and substitutes the result into eqs A35-A37, one obtains 2 d [ n ] / d ~= -[k] - [m]/V (A38) Since the right-hand side of eq A38 can be calculated explicitly in terms of the values of 0and n on the two sides of the discontinuity, eq A38 is an explicit differential
Ind. Eng. Chem. Res., Vol. 30, No. 5,1991 863 equation for the time evolution of [ n ]and , it can therefore be used instead of eqs A33 and A34 to obtain the value of [n](and thus also of [PI) at time 7 t 87. In the particularly simple case of the SDP, [m]= 0 and [ k ] = [n], so that eq A38 provides a very simple explicit differential equation for [ n ] . A differential equation for [n]can also be obtained in the general case. However, it contains additional jumps of the derivatives of the dependent variables. We next consider dealing with the instants of reflection. Consider the beginning of an odd-N interval (the case of an even-N interval is again analogous). The values of all quantities of interest to the left of the discontinuity at T~ + 0 may be calculated as discussed in the previous section; let these values be denoted by an asterisk. While there is no problem integrating forward in the SD (on the right) , ED initially has a zero thickness, starting at time T ~ the and hence, a difficulty arises. This can be circumvented as follows. Consider time T N + t, where E is an appropriately small time interval. One can derive a linearized solution, which satisfies the differential equations and boundary conditions to within O(t). This can always be done, as the following simple example shows. Consider the case where Y ( T=) r. If one writes for the (linearized) distribution of 4 at time TN + t P = A B y + CT (A391
+
one immediately sees that C = 0 and A = r, since otherwise the BC could not be satisfied. It is easy to convince oneself that, even when pl’ is not constant, only the dolay derivative, Le., the constant B, remains arbitrary. Since C = 0, eq A1 now implies that anlay = -m*,and hence, the linearized distribution of n introduces only one more arbitrary constant: n = n* - m*y + d7 (~40)
To within O(c), the fact that 4 is a linear function of y implies that so, too, is Q; Le.,
Q = Q*
+ by
(~41)
Equation A2 now implies that
d+b+k*=O (A42) An additional condition relating constants b and d is obtained by enforcing eq A24, noting that the values to the left of the discontinuity are nL = n* + ( d - m*V)t (A43) QL
= Q* + bVt
(A44)
The constants b and d are now determined, and it takes tedious but straightforward algebra to establish that eq A21 is also satisfied to within O(t). Equations A38 and A39 are now used as the initial distributions in the forward integration of the equations in the ED. The procedure sketched above works in principle for all free boundary problems but is particularly useful in the case of hyperbolic sets of equations. The linearized solution holds at time T~ t, when the thickness of the ED is V*c. This thickness must be divided into a finite number of grid points, e.g., 11,so that the initial space step is 6y = V*t/lO. Should the condition of numerical stability for explicit solutions to parabolic problems (as arises with ordinary diffusion) hold, one would need to choose an initial time step b~ no larger than 6y2/2, i.e., of the order of t2/100. This would make forward numerical integration very slow indeed. However, as was observed by Ocone (19871, the parabolic condition for numerical stability may be grossly violated with hyperbolic problems, and one
+
could in fact use t itself as the initial time step-though we have not done so (see below). Let Y ( T )be the thickness of the ED. For an odd-N interval, Y = s; for an even-N interval, Y = 1 - s. The thickness of the SD is of course 1 - Y. At each time step, the value of Y changes, and therefore, if it is divided into a fixed number M of steps, 6y, new grid points need be established. The values of the variables at the new grid points are calculated by interpolation from the values at the old grid points, and this imposes the restriction that 67 should be less than Y / M ) q . Strictly speaking, there is no such restriction in the SD, but we have chosen to apply it as well; i.e., 67 must also be less than (1 - Y)/M(w. Therefore, 87 begins at the value €/it4q, 1 increases until s = 0.5, and then decreases once again. When it falls below 6 , the last time step can be chosen to end exactly at the instant of reflection. Finally, an internal check of the consistency of the numerical procedure is always possible. The function C * ( ~ , T ) can always be calculated by substituting ~ ( Y , T and ) /~(Y,T) into C(h,P),and C*(y,O) can similarly be determined. One may first integrate eq A1 from y = 0 to y = 1 and then integrate the result from zero to time T . The result is JTn(O,it)d7’ -
n”(7’) dT’ = soT
x ’ C * ( y ’ , ~ dy’) slC*(y’,O) 0 dy’ (A45) The two sides of the equation (which, in the physical example of diffusion in polymers, represent the total weight sorbed) can be calculated separately, and they must agree within numerical accuracy.
Literature Cited Aifantis, E. C. On the Problem of Diffusion in Solids. Acta Mech. 1980, 37, 265. Astarita, G. Scaleup: Overview, Closing Remarks, and Cautions. In Scaleup of Chemical Processes; Bisio, A., Kabel, A. L., Eds.; J. Wiley: New York, 1985. Astarita, G. Thermodynamics. An Advanced Textbook for Chemical Engineers; Plenum Publishing Co.: New York, 1989. Astarita, G.; Joshi, S. R. Sample Dimension Effects in the Sorption of Solvents in Polymers. A Mathematical Model. J. Membr. Sci. 1978, 4, 165. Astarita, G.; Kenny, J. M. The Stefan and Deborah Numbers in Polymer Crystallization. Chem. Eng. Commun. 1987,53, 69. Astarita, G.; Marrucci, G. Principles of Non-Newtonian Fluid Mechanics; McGraw Hill: New York, 1974. Astarita, G.; Meldon, J. H. Modeling Anomalous Diffusion Behavior in Polymers. Latin Am. Appl. Res. 1989, in press. Astarita, G.; Ocone, R. Discontinuities and Interfaces in Transport Phenomena Coupled with Kinetics and Relaxation. Adu. Tramp. Processes 1989, in press. Astarita, G.; Sarti, G. C. A Class of Mathematical Models for Sorption of Swelling Solvents in Glassy Polymers. Polym. Eng. Sci. 1978, 18, 388. Astarita, G.;Paulaitis, M. E.; Wissinger, R. G. Thermodynamics of the Glass Transition. J. Polym. Sci. B 1989, 27, 2105. Bagley, E.; Long, F. A. Two-Stage Sorption and Desorption of Or1955, 77, ganic Vapors in Cellulose Acetate. J . Am. Chem. SOC. 2172. Baird, B. R.; Hopfenberg, H. B.; Stannett, V. T. The Effect of Molecular Weight and Orientation on the Sorption of n-Pentane by Glassy Polystyrene. Polym. Eng. Sci. 1971, 11, 274. Billovitz, G. F.; Durning, C. J. Viscoelastic Diffusion in the Polystyrene-Ethylbenzene System. AIChE Meeting, San Francisco, Nov 1989; Paper 210e. Bowen, R. M. Theory of Mixtures. In Continuum Physics, III; Eringen, A. C., Ed.; Academic Press: New York, 1976. Cattaneo, C. Sulla Conduzione de Calore. Atti Semin. Mat. Fis. Univ. Modena 1948, 3, 83. Cattaneo, C. Sur une Forme de l’Equation de la Chaleur Eliminant le Paradoxe d‘une Propagation Instantanee. C.R. Acad. Sci. 1958, 247, 431.
864 Ind. Eng. Chem. Res., Vol. 30,No. 5, 1991 Cohen, D. S. Mathematical Models for non-Fickian Sorption and Diffusion in Polymers. AIChE Meeting, San Francisco, Nov 1989; Paper 210c. Crank, J. A Theoretical Investigation of the Influence of Molecular Relaxation and Internal Stress on Diffusion in Polymers. J. Polym. Sci. 1953, l l , 151. Crank, J. Mathematics of Diffusion; Clarendon Press: Oxford, 1975. Crank. J. Free and Moving- Boundarv Problems; Clarendon Press: Oxford, 1984. De Socio, L. M.; Gualtieri, G. A Hyperbolic Stefan Problem. Q. Appl. Math. 1983,42, 253. Doghieri, F.; Carbonell, R. G.; Sarti, G. C. Stress, Deformation and Concentration Profiles in the Diffusion through Solid Polymers. AIChE Meeting, San Francisco, Nov 1989; Paper 210g. Duda, J. L. Personal communication, 1989. Durning, C. J. Differential Sorption in Viscoelastic Fluids. J. Polym. Sci., Polym. Phys. Ed. 1985, 23, 1831. Durning, C. J. Personal communication, 1989. Durning, C. J.; Spencer, J. L. Differential Sorption and Permeation in Viscous Media. J. Polym. Sci., Polym. Phys. Ed. 1985,23, 171. Durning, C. J.; Tabor, M. Mutual Diffusion in Concentrated Polymer Solutions under a small Driving Force. Macromolecules 1986, 19, 2220.
Enscore, D. J.; Hopfenberg, H. B.; Stannett, V. T. Effect of Particle Size on the Mechanism Controlling n-Hexane Sorption in Glassy Polystyrene Microspheres. Polymer 1977, 18, 1105. Fasano, A,, Primicerio, M., Eds. Free Boundary Problems: Theory and Applications; Pitman: Houston, 1983. Faulkner, D. L.; Hopfenberg, H. B.; Stannett, V. T. Residual Solvent Removal and n-Hexane Sorption in Blends of Atactic and Isotactic Polystyrene. Polymer 1977, 18, 1130. Fick, A. Uber Diffusion. Pogg. Ann. Phys. Chem. 1855, 94, 58. Frisch, H. L. Isothermal Diffusion in Systems with Glasslike Transitions. J. Chem. Phys. 1964, 41, 3679. Frisch, H. L. Sorption and 'hansport in Glassy Polymers-A Review. Polym. Eng. Sci. 1980, 20, 2. Goddard, J. D.; Schulz, J. S.; Bassett, R. J. On Membrane Diffusion with Near Equilibrium Reaction. Chem. Eng. Sci. 1970,25,665. Greenberg, J. M. A Hyperbolic Heat Transfer Problem with Phase Changes. IMA J. Appl. Math. 1987, 38, 1. Joshi, S. R.; Astarita, G. Diffusion-Relaxation Coupling in Polymers which show Two-Stage Sorption Phenomena. Polymer 1979,20, 455.
Kambour, R. P.; Karasz, F. E.; Daane, J. H. Kinetic and Equilibrium Phenomena in the System Acetone Vapor and Polycarbonate Film. J. Polym. Sci. A-2 1966, 4, 327. Kotchine, P. Sur la Theorie des Ondes de Choc dans un Fluide. R. Circ. Matem. Palermo 1926,50, 305. Lax, P.; Wendroff, B. Systems of Conservation Laws. Comm. Pure Appl. Math. 1960, 13, 217. Long, F. A.; Richmann, D. Concentration Gradients for Diffusion of Vapors in Glassy Polymers and Their Relation to Time Dependent Diffusion Phenomena. J. Am. Chem. SOC.1960,85,513. Lustig, S. R.; Caruthers, J. M.; Peppas, N. A. Development and Implementation of a Model for Penetrant Transport in Glassy Polymers. AIChE Meeting, San Francisco, Nov 1989; Paper 210f. Maxwell, J. C. On the Dynamical Theory of Gases. Phil. Trans. R. SOC.London 1867, 157,49. Mercier, J. P.; Groeninckx, G.; Lesne, M. Some Aspects of Vapor Induced Crvstallization of Polvcarbonate of Bisphenol A. J . Polym. Sci. C-1967, 16, 2059. Mueller, I.; Villaggio, P. Conditions of Stability and Wave Speeds for Fluid Mixtures. Meccanica 1976, 1I, 191. Neogi, P. Anomalous Diffusion of Vapors through Solid Polymers. I, Irreversible Thermodynamics of Diffusion and Solution Priocesses. AZChE. J. 1983a, 29,829; 11, Anomalous Sorption. AIChE J . 1983b, 29, 833.
Norman, K. N. Stress Assisted Diffusion of Perfect Fluids in Simple Solids with Fading Memory. AIChE Meeting, San Francisco, Nov 1989; Paper 210d. Ocone, R. Cristallizzazione di Polimeri Caratterizzati da un Tempo di Rilassemento del Flusso Termico. Chem. Eng. Thesis, University of Naples, 1987. Ocone, R.; Astarita, G. Continuous and Discontinuous Models for Transport Phenomena in Polymers. AIChE J. 1987, 33, 423. Overbergh, N.; Berghmans, H.; Smets, G. Crystallization of Isotactic Polystyrene Induced by Organic Vapors. Polymer 1975,16,703. Ozisik, M. N. Heat Conduction; J. Wiley: New York, 1980. Peppas, N. A,; Urdahl, K. G. Anomalous Penetrant Transport in Glassy Polymers. VII, Overshoots in cyclohexane uptake in crosslinked polystyrene. Polym. Bull. 1986, 16, 201. Peppas, N. A.; Urdahl, K. G. Anomalous Penetrant Transport in Glassy Polymers. VIII, Solvent Induced Cracking in Polystyrene. Eur. Polym. J. 1988, 24, 13. Petropoulos, J. H.; Roussis, P. P. Study of Non-Fickian Diffusion Anomalies Through Time Lags. J. Chem. Phys. 1967,47,1491. Petropoulos, J. H.; Roussis, P. P. Anomalous Diffusion of Good and Poor Solvents or Swelling Agents in Amorphous Polymers. J. Polym. Sci. 1969, C22, 917. Sadd, M. H.; Didlake, J. E. Non-Fourier Melting of a Semiinfinite Solid. J. Heat Trans., Trans. ASME 1977, 99, 25. Salomon, A. D.; Alexiades, V.; Wilson, D. G.; Drake, J. On the Formulation of Hyperbolic Stefan Problems. Q.Appl. Math. 1985, 43, 295.
Sarti, G. C. Personal communication, 1989. Sarti, G. C.; Carbonell, R. Paper presented a t the AIChE Meeting, New York, 1987. Showalter, R. E.; Walkington, N. J. A Hyperbolic Stefan Problem. Q. Appl. Math. 1987,44, 769. Smith, M. J.; Peppas, N. A. Effect of the Degree of Crosslinking on Penetrant transport in Polystyrene. Polymer 1985, 26, 569. Sod, G. A. A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws. J. Comput. Phys. 1978, 27, 1. Stefan, J. Ueber die Theorie der Eisbildung insbesondered ueber die Eisbildung in Polarmeere. Ann. Phys. Chem., N.F. 1891,42,261. Titow, W. V.; Braden, M.; Currell, B. R.; Loneragan, R. J. Diffusion and some Structural Effects of Two Chlorinated Hydrocarbon Solvents in Bisphenol A Polvcarbonate. J. Appl. _ _ Polym. Sci. 1974, 18, 867.
Truesdell, C. A. Mechanical Basis for Diffusion. J. Chem. Phys. 1962.37. 2336.
Truesdell, C. A. Rational Thermodynamics, 2nd ed.; Springer Verlag: New York, 1985. Vrentas, J. S.; Duda, J. L. Diffusion in Polymer-Solvent Systems. 111, Construction of Deborah Number Diagrams. J. Polym. Sci., Polym. Phys. Ed. 1977, 15,441. Vrentas, J. S.; Jarzebski, C. M.; Duda, J. L. A Deborah Number for Diffusion in Polymer-Solvent Systems. AIChE J. 1975,21,894. Vrentas, J. S.; Duda, J. L.; Hou, A. C. Anomalous Sorption in Polyethyl Methacrylate. J. Appl. Polym. Sci. 1984, 29, 399. Vrentas, J. S.; Duda, J. L.; Huang, W. J. Regions of Fickian Diffusion in Polymer-Solvent Systems. Macromolecules 1986, 19, 1718. Ware, R. A.; Tirtowidjojo, S.; Cohen, C. Diffusion and Induced Crystallization in Polycarbonate. J. Appl. Polym. Sci. 1981,26, 2975.
Waywood, W. J.; Durning, C. J. Solvent Induced Crystallization of a Compatible Polymer Blend. Polym. Eng. Sci. 1987, 27, 1265. Windle, A. H. Case-I1 Sorption, in Polymer Permeability. J.Comyn Ed.; Elsevier: London, 1985. Received for review April 4, 1990 Revised manuscript received July 13, 1990 Accepted August 1, 1990
