1012
Ind. Eng. Chem. Res. 1996, 35, 1012-1023
Convective-Diffusive Mass Transfer with Chemical Reaction in Squeezing Flows† A. S. Nagarkar, P. Arce,* and B. R. Locke Department of Chemical Engineering, FAMU-FSU College of Engineering, 2525 Pottsdamer Street, Tallahassee, Florida 32310-6046
A model for convective-diffusive mass transfer in squeezing flow is formulated and solved in the analysis reported in this contribution. A spatial averaging technique is applied to reduce the 2D problem to a 1D problem. It is assumed that the hydrodynamics is decoupled from the mass transfer and the model is solved using a sequential coupling between the convective velocity field and concentration field equations. To account for the loss of information due to averaging, a closure scheme is developed for the deviations of average values of velocity and concentration. This involves derivation of the equations for the deviations of the respective fields. The effective transport parameters are evaluated in terms of the molecular diffusion coefficient, the average velocity, and the dimensions of the system. The associated eigenvalue problem is solved by using a combination of exponentials and Kummer functions. The concentration profile is then found analytically by a Fourier series expansion of the initial condition. Two specific cases are considered in this paper. These include (1) no reaction in the bulk and (2) a linear homogeneous reaction in the bulk. The analysis presented in this paper promotes understanding of the behavior of the mass-transfer problem in, for example, polymer reaction molding. 1. Introduction The problem under investigation here is that of a squeezing flow between two parallel plates which contain a Newtonian fluid that is assumed incompressible and under isothermal conditions. The flow between the two parallel plates is produced by action of a constant force that has been applied normally to the plate (see Figure 1). In addition to the basic hydrodynamic aspects, the present study also considers the analysis of the convective-diffusive transport with and without chemical reaction. This type of problem is of relevance to several important engineering applications, which include, for example, the following: (1) The basic flow configuration is used in an industrial testing device that has been called a “parallel-plate plastometer” (Dienes and Klem, 1946; Oka, 1960), a “parallel-plate plastimeter” (Scott, 1931, 1932), a “parallel-plate viscometer” (Gent, 1960), a “traverse-flow viscometer” (van Wazer et al., 1963), and a “compression plastometer” (Mooney, 1958). This device is useful for investigating materials with high viscosity at relatively high temperatures and/or high shear rates. In spite of the fact that this process is time dependent, the equations for the hydrodynamic field have usually been used under the quasi-steady-state assumption (Denn, 1980). (2) Squeezing flows have also been studied in the past in association with lubrication devices (Metzner, 1968a,b; Tanner, 1985; Bird et al., 1987) where the lubrication approximation has been used to analyze the hydrodynamics of cases with walls moving in relative lateral motion [see, for example, Langlois, 1962]. Comparatively, less attention has been given to the analysis of situations where the pressure generated in the fluid film is due to the relative normal motion of the two surfaces or in journal bearings (Gross, 1960). In addition, most of the published work has considered only the hydro† Paper submitted by invitation to the Ind. Eng. Chem. Res. special issue honoring Professor J. King on the occasion of his 60th birthday.
0888-5885/96/2635-1012$12.00/0
Figure 1. Geometrical sketch of the squeezing device with the coordinate system.
dynamics of the system, with virtually no analysis of the basic convective-diffusive transport of mass in the system. (3) The pure hydrodynamic problem associated with squeezing flows is of intrinsic interest to rheologists because of their interest in evaluating rheological equations under transient conditions and in developing numerical techniques for computing nonviscometric flows (Tanner, 1985; Bird et al., 1987). The flow is interesting because the hydrodynamic field depends upon both the radial and axial coordinates and also because the kinematics is not known a priori and both shearing and elongational forces are present. (4) Examples of the squeezing flow arise in plastic fabrication operations such as in stamping and injection molding (Dienes, 1947). The squeezing flow problem applies to the processing of composite materials such as fiber-reinforced matrices [see, for example, Coleman, 1990; Skartsis et al., 1992] and to the polymer reaction molding processes [see, for example, Broyer and Macosko, 1976]. Furthermore, most of the work reported on this topic has used numerical techniques. During the fabrication of fiber-reinforced plastics, the material is usually in a slowly moving molten state and it can be regarded as a fiber-reinforced Newtonial fluid. In addition, several manufacturing processes involve an important squeezing action (Coleman, 1990). Most of the applications mentioned above only considered the hydrodynamic aspects. The transport of mass by convection and diffusion is particularly important in the © 1996 American Chemical Society
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1013
last class of applications mentioned above. Furthermore, in some cases the convective-diffusive transport must be coupled with chemical reaction in order to have a meaningful characterization of the overall operation in the system. The analysis presented in this paper focuses on this particular aspect of the problem. Processes with convection, diffusion, and chemical reaction in squeezing flows found in practical applications are extremely complicated. Therefore, our aim in this paper shall be rather modest and will focus on gaining a basic understanding of the behavior of systems under no reaction conditions and when it is coupled with a firstorder chemical reaction. One of the objectives is to derive effective transport coefficients such as the effective (axial) diffusivity and the effective convective velocity in terms of basic geometrical parameters and transport properties. The effective transport coefficients and the averaged(macro)transport equations are derived by following an area averaging approach performed on the microtransport equations of the original differential model. The approach follows ideas of Whitaker (1967, 1973), Slattery (1967, 1981), and Carbonell (see, for example, Zanotti and Carbonell, 1984). This methodology applied to the analysis of the squeezing flow problem in this paper allows the derivation of effective transport coefficients (i.e., effective diffusivity and effective convective velocity) in terms of fundamental parameters and geometrical dimensions of the physical system. Because of this functionality, the equations are useful for characterizing the overall transport in the system. Eventually, they will also be useful for developing guidelines for process design. The transient averaged transport equations, featuring the effective transport coefficients, are solved by using a Fourier expansion in terms of the eigenvalues of the associated Sturm-Liouville problem. The solution to this problem can be obtained analytically in terms of the cylindrical parabolic functions (Abramawitz and Stegun, 1972) that can be calculated in terms of Kummer’s functions. Details about the model formulation, the averaging approach, and the solution aspects are given in the sections below.
of continuity and motion (Denn, 1980) is then coupled with the molar species continuity equation to determine the concentration profile. The equation of continuity for a Newtonian fluid under incompressible flow conditions in two-dimensional rectangular coordinates is given by
∂vx ∂vz + )0 ∂x ∂z
and x and z components of the equation of motion are given by
[ [
] ]
0)-
∂2vx ∂2vx ∂p +µ + 2 ∂x ∂x2 ∂z
(2)
0)-
∂2vz ∂2vz ∂p +µ + 2 ∂z ∂x2 ∂z
(3)
Equations 2 and 3 are components of the NavierStokes equation where the inertial terms have been neglected due to the creeping flow approximation. It is likely that closer to each plate, due to the applied force, the fluid moves in the z direction uniformly at all x positions. Therefore, vz can be assumed to be independent of x along each plate. Hence, assuming vz ) φ(z) and substituting in eq 1
∂vz dφ(z) )∂x dz
(4)
which gives
vx ) -
dφ(z) x + c1 dz
(5)
Hence at x ) 0, due to the fact that the problem is symmetric around the z axis, vz ) 0, which yields c1 ) 0, and
vx ) -
2. General Model Formulation Consider two flat plates moving toward each other under the influence of an applied force, F (see Figure 1). The system is assumed isothermal, and, therefore, the physical properties of the system remain constant. The y-directional variations are neglected, and it is assumed that there is no flow in the y-direction. The flow in the device is considered to be under creeping flow conditions, and therefore all the inertial terms in the equation of motion are neglected. Also, as in previous studies (Denn, 1980), the hydrodynamic field equations will be considered under a quasi-steady-state condition. This leads to the conclusion that the equations of motion feature time-dependent terms only through the velocity of the plates and the height or distance between the plates, V(t) and H(t). The model formulation presented below first introduces the differential equations that describe the hydrodynamics, and then the molar species conservation equations are given. Because of the assumptions made above, the hydrodynamic field equations are decoupled from the molar concentration species continuity equation. This fact allows the solution to the hydrodynamics equation to be derived independently of the concentration field. The velocity profile obtained by solving the equations
(1)
dφ(z) x dz
(6)
It should be noted that, in general, vz depends on t, but the explicit time dependence disappears since quasistate-state flow is assumed. The boundary conditions are given by the no-slip condition at the plate
vz ) V vz ) -V
at z ) -H
(7)
at z ) H
(8)
The fluid adheres to the plate and hence moves with the same velocity as the plates. The negative sign is utilized at z ) H because the bottom plate moves in a direction opposite to the top plate. Since the system is assumed symmetric with respect to the coordinate x and because of the squeezing flow conditions described above, the kinematics of the system implies that
vx ) 0
at x ) 0, z ) H, z ) -H
(9)
The dilute solution molar species continuity equation in two-dimensional rectangular coordinates neglecting homogeneous chemical reactions is given by [Bird et al., 1960]
1014
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996
[
]
∂cA ∂cA ∂2cA ∂2cA + vx )D + 2 ∂t ∂x ∂z2 ∂x
(10)
where D is the diffusion coefficient of species A. In eq 10 the convective transport in the z direction has been neglected under the assumption that vx . vz in the system. The initial state is assumed to be uniform with respect to both the x and z coordinates. Thus,
cA ) c0
at t ) 0
(11)
The assumption of impermeable plates gives rise to the no-flux boundary conditions, for the solute species which are given by
∂cA )0 ∂z
at z ) H and z ) -H
(12)
The symmetry boundary conditions for the solute molar species concentration is given by
∂cA )0 ∂x
at x ) 0
at x ) L
(14)
An alternative to this condition would be a flux or Danckwerts type of boundary condition (Nagarkar, 1994). The analysis reported in this paper can be extended to include this or other more general boundary conditions. 3. Solution to the Hydrodynamic Problem In order to solve the molar species continuity equation, written above, the velocity profiles are required in eq 10. This is the task performed in this section below. Based on the fact that the hydrodynamic problem is sequentially coupled to the concentration field, the hydrodynamics, i.e., v[vx,vz], field can be computed from the hydrodynamic model given in section 2. The velocity component in the z direction is given by
vz )
{
}
z3 V 3z + 3 2 H H
(15)
The velocity component in the z direction is a function only of z coordinate of the system. The velocity component in the x direction, vx, is, however, a function of both x and z and is given by
vx )
{
}
3Vx z2 1- 2 2H H
(16)
The time dependency of the solution is accounted for through H(t) and V(t), as was pointed out in the previous section. The pressure profile is obtained from the equation of motion given in section 2 to produce [Denn, 1980]
[(
)
Ft 1 1 ) 2+ H (t) H0 µWL3 2
(18)
where H0 is the spacing between plates at t ) 0. Equation 18 shows that the variation of H(t) with time follows a function that starts at H0 and decays asymptotically toward zero when time increases. The function shows faster decay when larger values of the force, F, are applied. Also, fluids with increasing values in the viscosity, µ, will delay the decay of the values of the function H(t). In order for the quasi-steady-state approximation to hold, the Reynolds number should be less than 1, which gives rise to a constraint on the force (Denn, 1980)
F,
µ WL3 F H4
(19)
(13)
The system is assumed to interact with a reservoir with a constant concentration of the reactant species at x ) L. Therefore, a Dirichlet boundary condition is imposed at the external edge of the plate to give
cA ) c∞
coordinates. The Stefan equation, which gives the variation of the height, H(t), of the plates with time is obtained by integrating the pressure profile over the x direction and using a substitution V(t) ) -dH/dt to yield
]
L2 - x2 3µV z2 1 + p) 2H H2 H2
(17)
The pressure field is a function of both the x and z
This constraint can also be derived by an order of magnitude analysis performed directly on the NavierStokes equation (Denn, 1980). The veracity of eq 18 has been tested experimentally for other geometries and for different types of fluids. The agreement between theoretical predictions and experimental measurements is fairly good. These experiments are relevant for the analysis performed in this paper because they provide an indirect way of checking the quasi-steady-state assumption. As stated in the Introduction section, one of the purposes of the analysis performed here is to derive area-average transport equations and identify equations for the effective parameters. The area-averaging approach used on the microscopic molar species conservation equation (see section 4) requires an equation for the averaged velocity in the x direction (i.e., vx velocity component). In order to accomplish this task, first, for a given value of the axial coordinate, x of the parallelplate channel, the following function u(z) associated with vx (see eq 16) can be identified:
u(z) ≡
(
)
3V z2 1- 2 2H H
(20)
and using the definition of an area average gives H u(z) dz V ∫-H 〈u(z)〉 ≡ ) H H dz ∫-H
(21)
which yields the value V/H for the area average of the function u(z). Now, by recalling the form of the velocity profile vx given by eq 16 and by using eq 21, the following value of the area-average velocity profile is obtained
〈vx〉 ) (V/H)x
(22)
The averaged value of the velocity profile, i.e., 〈vx〉, can be used to determine the deviation of the point velocity field from the average, i.e., by Gray’s decomposition (see Gray, 1975). The difference of function u(z) identified above and its average value given by eq 1
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1015
〈cAu〉 ) 〈cA〉〈u〉 + 〈cˆ Auˆ 〉
yields
uˆ ) u(z) - 〈u〉
(23)
Integrating eq 30 into the averaged equation (27) gives
∂〈cˆ Auˆ 〉 ∂〈cA〉 ∂2〈cA〉 ∂〈cA〉 + x〈u〉 +x )D ∂t ∂x ∂x ∂x2
By replacing the results given by eq 21 and eq 20 into eq 23, the following expression is obtained for the deviation field, uˆ
(
uˆ ) 〈u〉
) (
3z2 3z2 V 1 1 ) 2 2H2 H 2 2H2
)
(24)
The usefulness of the relationships derived in this section will be clearly demonstrated in the area-averaging process of the molar species continuity equation. This is the subject matter of the next section. 4. Averaged Species Continuity Equation This section deals with the derivation of the averaged species continuity equation by applying the area-average technique following the approach proposed by Whitaker (1967, 1973), Slattery (1967, 1981), and Carbonell (see Zanotti and Carbonell, 1984). The central equation in the analysis is the species continuity equation introduced in section 2. However, results derived in the previous section and related to the velocity field are needed because of the sequential coupling between the species continuity equation and the hydrodynamic field equations. The analysis begins with the definition of the area average of the concentration profile and follows with the application of such procedures to the molar species continuity, the development of a closure problem, the identification of the effective transport parameters, and the formulation of the averaged equations for the concentration field. Defining the averaged concentration in the z direction as H cA dz ∫-H 〈cA〉 ≡ H dz ∫-H
(25)
∂cˆ A ∂〈cA〉 ∂cˆ A ∂〈cA〉 ∂〈cA〉 + + x〈u〉 + xuˆ + x(〈u〉 + uˆ ) ) ∂t ∂t ∂x ∂x ∂x 2 2 2 ∂ 〈cA〉 ∂ cˆ A ∂ cˆ A D + 2 + 2 (32) ∂x2 ∂x ∂z
[
]
The equation for determining the deviation of the concentration can be obtained by subtracting the above equation from the averaged molar species continuity equation (31) to give
∂cˆ A ∂t
+ xuˆ
∂cˆ A ∂〈uˆ cˆ A〉 ∂〈cA〉 + x(〈u〉 + uˆ ) -x ) ∂x ∂x ∂x 2 ∂ cˆ A ∂2cˆ A D + 2 (33) ∂x2 ∂z
[
]
In order to simplify the equation for the deviation field (eq 32), an order of magnitude analysis is carried out (Howes and Whitaker, 1985; Zanotti and Carbonell, 1984). Assuming that the value of this average concentration is much larger than the deviation value of the concentration, i.e., 〈cA〉 . cˆ A, gives
uˆ
∂cA ∂cA ∂2cA ∂2cA + vx )D + ∂t ∂x ∂x2 ∂z2
uˆ
(26)
(31)
This equation contains a mixture of averaged and deviation field variables that need to be removed to yield an equation that describes the transport process solely in terms of the averaged quantities. This task is accomplished by a closure scheme discussed below. 4.1. Closure Problem. A closure scheme is developed to resolve the terms containing the deviation field of the concentration and velocity profiles. Substituting eqs 30, 28, and 23 into eq 10 gives
and applying this average to the molar species continuity equation (10) gives
〈 〉 〈 〉 [〈 〉 〈 〉]
(30)
∂〈cA〉 ∂x
(34)
∂〈cA〉 ∂〈uˆ cˆ A〉 . ∂x ∂x
(35)
∂cˆ A ∂x
, uˆ
By using vx ) xu(z), where u(z) is the function identified in eq 20 and the no-flux conditions at z ) +H and Z ) -H, eq 26 gives
Also assuming a long channel approximation and the quasi-steady-state approximation for the deviation field of the concentration profile, the following equation can be written:
∂〈cA〉 ∂〈cAu〉 ∂2〈cA〉 +x )D ∂t ∂x ∂x2
2 2 1 ∂ cˆ A 1 ∂ cˆ A , 2 2 L2 ∂x2 H ∂z
(27)
(36)
Using Gray’s decomposition (Gray, 1975) for the velocity profile as indicated by eq 23 and for the concentration profile
Finally after neglecting the higher order terms, the equation for the deviation field, i.e., eq 33, of the concentration becomes
cA ) 〈cA〉 + cˆ A
(28)
∂2cˆ A ∂〈cA〉 )D 2 xuˆ ∂x ∂z
(29)
Following the approach of Paine et al. (1983), a solution of the type
the following equation can be derived
〈cAu〉 ) 〈(〈cA〉 + cˆ A)(〈u〉 + uˆ )〉
where cˆ A and uˆ are the deviations in concentration and velocity, respectively. Expansion of eq 29 noting that the averages of the derivations are zero results in
cˆ A ) g(z)
d〈cA〉 x dx
(37)
(38)
1016
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996
can be proposed. Substituting this solution into the deviation equation (37) leads to
Dg′′(z) - uˆ ) 0
(39)
The following constraints for the deviation field, cˆ A, can be obtained using the boundary condition (eq 12)
∂cˆ A ∂z
)0
at z ) H and z ) -H
(40)
and also for the field g(z)
g′(z) ) 0
at z ) H and z ) -H
(41)
Integration of eq 39 requires the specification of an additional condition since the no-flux condition g′(z) ) 0 at both ends allows the calculation of only one integration constant. An additional constraint on the function g(z) arises from the fact that the average of the deviation is zero and is given by H g(z) dz ) 0 ∫-H
〈g(z)〉 )
(42)
evaluate effective transport parameters the proposed solution (eq 47) is used in the averaged molar species continuity equation. By substituting this proposed solution in the averaged molar species continuity equation, it is possible to write
〈(
(44)
where
)
Equation 49 is useful for defining the effective diffusivity and the effective convective velocity of the system. Thus, the effective convective velocity is given by Veff ) xUeff where
Ueff ≡ 〈u〉 + 〈uˆ g(z)〉
( )
1 1 〈u〉, R ) 〈u〉 4D 8DH2
(45)
7H (120D )〈u〉
H uˆ g(z) dz ∫-H 〈uˆ g(z)〉 ) H dz ∫-H
(46)
γ)-
The function g(z) is necessary for expressing the deviation of the concentration field, cˆ A, as (see eq 38)
(
) (
)]
d〈cA〉 1 2 1 7H z z4 〈u〉x cˆ A(z,x) ) 2 4D 120D dx 8DH (47) 2
which is given entirely in terms of known functions and the first derivative of the average concentration profile, 〈cA〉. Equation 47 provides an interesting way to interpret the physical meaning of the deviation of the concentration field, cˆ A. This variable can be viewed as a modified convective mass flux in the axial direction, x, of the device. The coefficient proceeding d〈cA〉/dx in eq 47 is not the convective effective velocity of the system since this must be determined by using the averaged molar species continuity equation (see eq 31). This is the subject matter of section 4.2. 4.2. Effective Parameters. Effective parameters play a role in characterizing the overall macrotransport in the system and also in regaining the information lost in the averaging approach. These effective transport parameters take into account the respective deviation fields analyzed in the previous section. In order to
(51)
These two equations show the relationships between such fundamental parameters as the molecular diffusivity, D, with variables related to the deviation field, i.e., uˆ and g(z). Since these fields have been already determined in the previous sections, it is possible to obtain the expression for the term 〈uˆ g(z)〉. The deviation component 〈uˆ g(z)〉 is then found by averaging over the z direction
2
)
(50)
Deff ≡ D - x2〈uˆ g(z)〉
and γ is computed by using eq 43 to be
[(
∂x2 (48)
∂〈cA〉 ∂2〈cA〉 ∂〈cA〉 + x[〈u〉 + 〈uˆ g(z)〉] ) [D - 〈x2uˆ g(z)〉] ∂t ∂x ∂x2 (49)
(43)
g(z) ) Rz2 + β + γ
(
∂2〈cA〉
and the effective diffusivity by
〈u〉 3〈u〉 2 z 2D 2DH2
The solution to eq 43 can now be readily computed to be
β) -
)D
A reorganization of the terms of the above equation leads to
In order to solve eq 39, the function uˆ (z) given by eq 24 must be used to yield
g′′(z) )
)〉
d〈cA〉 ∂ uˆ g(z) x ∂〈cA〉 ∂〈cA〉 dx + x〈u〉 +x ∂t ∂x ∂x
(52)
By substituting the appropriate terms into eq 52, the following result is obtained:
〈uˆ g(z)〉 )
1 2H
[
2
]
H 1 3z [Rz2 + βz4 + γ] dz ∫-H 2 2H2
(53)
which, after the integration, leads to
〈uˆ g(z)〉 )
[ (
)
H5 H3 1 〈u〉 2R + 2β + 2γ 2H 2 3 5 3〈u〉 H7 H3 H5 + 2β + 2γ 2R 5 7 3 2H2
(
)] (54)
Substituting the values of R, β and γ previously determined in the section above leads to
〈uˆ g(z)〉 )
(
)
-2H2 〈u〉2 105D
(55)
By invoking eqs 50 and 51 and using the equation above, the effective parameters can be related to the fundamental and/or geometrical quantities of the system by the following two equations:
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1017
a Figure 3. Effective velocity vs time and with the Pe number as a parameter (Pe ) 5, Pe ) 10, Pe ) 20).
b Figure 2. (a) Variation of effective velocity with respect to the axial position, x (Pe ) 5). (b) Variation of effective velocity with respect to the axial position, x (Pe ) 10).
[
]
〈u〉 H x 105D
(56)
x2〈u〉2H2 105D
(57)
2
Veff ) 〈u〉 -2 Deff ) D + 2
2
In order to illustrate the behavior of the effective transport parameters, several figures are constructed. For example, parts a and b of Figure 2 show that Veff increases with the coordinate x and with time. For the calculations Pe has been defined as VL/D. Calculations of Figure 2a were performed for Pe ) 5, while those of Figure 2b were performed for Pe ) 10. The general trend for the coefficients are qualitatively very similar for both cases; however, for the case of Pe ) 10 the values of the effective velocity are somewhat smaller than those for Pe ) 5. A comparison for three different values of the Peclet number and as a function of time is shown in Figure 3. The difference between, for example, Pe ) 5 and Pe ) 20 of the values of the effective velocity (Veff/V) is not very significant. Figure 4 shows the variation in the effective diffusivity as a function of the axial coordinate, x, and for different values of the Peclet number (i.e., Pe ) 5, Pe ) 10 and Pe ) 15 and Pe ) 20). The different cases meet at the value of 1 for the origin of the coordinate (x ) 0) as should be, and they increase with the increases in values of x for a given value of the Peclet number. This trend is more definite when the values of the Peclet
Figure 4. Variation of effective diffusivity in the x direction.
number increases. For example, the ratio Deff/D is below the value 2 for Pe ) 5; however, the same ratio goes to approximately 8.5 for the value Pe ) 20. Substituting the values of effective parameters in the deviation equation yields the following averaged equation that describes the convective-diffusive transport process in the axial direction of the plates.
∂〈cA〉 ∂〈cA〉 ∂2〈cA〉 + Veff ) Deff ∂t ∂x ∂x2
(58)
Note that this is a multiple-time-scale model. It is seen that 〈u〉 ) V(t)/H(t) is a function of time. It is assumed that the variation of V(t) with t is negligible as compared to the variation of 〈cA〉 with the same variable. Comparing the above equation to the Taylor-Aris results for flow in channels and tubes, it is seen that the effective parameters depend qualitatively on the axial parameter in a somewhat analogous fashion to that of the radial parameter in the Taylor-Aris result (Aris, 1956). The result for the squeezing flow, however, shows a quadratic dependence with the axial variable, x, because of the nature of the hydrodynamic field. 5. Solution to the Transient Problem The area average equation (58) can be written in a compact nondimensional fashion by using the definition of the effective transport coefficient, Veff and Deff,
1018
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996
identified previously and by using the following definition for parameters and variables:
x s≡ ; L
〈cA〉 〈c〉 ≡ ; c0
( )
UeffL2 Pe ≡ Deff (59)
UeffL τ) t; Pe
Equation 58 is transformed into the following parabolic differential equation with variable coefficients.
where φ0 is a constant to be determined from initial conditions. The ordinary differential equation associated with the function ψn(s) in terms of the axial variable, s, is given by
[
(60)
This equation may be written in a self-adjoint form (Ramkrishna and Amundson, 1985) by using the weighting factor
2
ds
s [-Pe 2 ]
2
[
]
〈c〉 - c∞ W≡ c∞
(62)
d2ψ′n(s) -
ds2
1 ∂W ∂W ∂ ) r(s) ∂τ ∂s Pe r(s) ∂s
]
(63)
Next, the separation of variables technique is used to obtain the function W from two independent functions φn(τ) and ψn(s) by the superposition principle (Haberman, 1987). There are two ordinary differential equations that can be separated and solved for each of the functions just mentioned. Let us denote Wn by the following product:
Wn ) φn(τ) ψn(s)
]
dψn(s) dφn(τ) d ψn(s) ) r(s) φn(τ) (65) dτ ds ds
Separating the functions φn(τ) and ψn(s) gives two differential equations
]
φ′n(τ) dψn(s) 1 d r(s) ) -λn2 ) ds φn(τ) r(s) ψn(s) ds
(66)
The solution of the transient ordinary differential equation, which is in terms of the function φn(τ), leads to an exponential function
φn(τ) ) φ0 exp(-λn τ)
(67)
]
s2Pe2 1 - Pe - λn2 ψ′n(s) ) 0 4 2
(71)
d2y 1 - z2 + a y ) 0 2 4 dz
(
)
(72)
where
Pe - λn2 2
a)-
(73)
and
z ) s Pe
(74)
The solution of the differential equation for the functions ψn(s) is now given by
(64)
where the function φn(τ) will feature the transient part of the problem and the function ψn(s) will show the spatial variation of the problem and the SturmLiouville problem for the system. By using the separation of variable principle (eq 63), the molar species continuity equation becomes
[
Comparing eq 71 with the parabolic cylinder functions (Abramowitz and Stegun, 1972), the following equations can be derived
Equation 62 is further transformed into the following equation after using eq 62
[
(70)
eq 69 is written as
(61)
In order to yield homogeneous boundary conditions, the following transformation is proposed:
2
(69)
(s 4Pe)
∂〈c〉 ∂〈c〉 1 ∂ ) r(s) ∂τ ∂s Pe r(s) ∂s
[
dψn(s) + λn2ψn(s) ) 0 ds
ψn(s) ) ψ′n(s) exp
to yield
r(s) Pe
- s Pe
The above differential equation is transformed into a parabolic cylinder equation whose solution can be expressed in terms of hypergeometric functions. Using the transformation
r(s) ≡ exp
[
(68)
which gives a second-order differential equation
d2ψn(s)
∂〈c〉 ∂2〈c〉 ∂〈c〉 Pe ) - sPe 2 ∂τ ∂s ∂s
]
dψn d r(s) ) -λn2ψn(s) ds ds
ψ′n(s) ) C1y1 + C2y2
(75)
where y1 is a fundamental solution of eq 72 and can be expressed in terms of the Kummer functions as
(
)
(
1 y1 ) exp - s2Pe2 M 4
-
)
Pe - λn2 1 1 s2Pe2 2 + , , 2 4 2 2
(76)
and y2 is the other fundamental solution of eq 72 and is given by
(
)
(
1 y2 ) exp - s2Pe2 M 4
-
)
Pe - λn2 5 3 s2Pe2 2 + , , 2 4 2 2
(77)
The boundary condition dψ′n(0)/ds ) 0 gives C2 ) 0; hence, the general solution to eq 71 is given by
(
)
(
1 ψ′n(s) ) exp - s2Pe2 M 4
-
)
Pe - λn2 1 1 s2Pe2 2 + , , 4 4 2 2 (78)
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1019
which leads to the solution of eq 69 ψn(s) that can be written as
ψn(s) )
[
](
2
)
1 1 s2Pe2 1 Pe λn exp - s2(Pe2 - Pe) M + , , 4 4 2 4 2 2 (79) The solution to the differential equation (69) given by eq 79 is an analytical solution that, to the best of the authors’ knowledge, has not been reported in the literature. This new aspect of such a solution is very important from the computational point of view in solving the Sturm-Liouville problem of the model. In this regard, the solution to the eigenvalues problem (see section 5.1) reduces to simply finding the roots of a polynomial type function. Also, the computation of the eigenfunctions of the problem can be obtained by simply evaluating the functions ψn(s) given in eq 79. Details about the aspects can be found in section 5.1. 5.1. Eigenvalue Problem. The calculation of eigenvalues, λn, is an important step in solving the transient problem for the concentration field. In order to obtain the values of λn, the characteristic equation of the problem must be derived. The characteristic equation can be obtained by using the analytical solution given by functions ψn(s) (see eq 79) and the remaining boundary conditions at the external edge of the channel, i.e., at s ) 1. In terms of the variables of interest this condition is given by
ψn(s) ) 0
at s ) 1
(80)
Figure 5. First five eigenfunctions for the case without reaction, Pe ) 5. ∞
W(s,τ) )
(
)
∞
W(s,τ) )
(81)
where M(an, bn, λn) is given by the series (Slater, 1960; Abramowitz and Stegun, 1972)
M(an, bn, λn) ) 1 +
|
(a)iλni
∑ i)1 (b) i! i
∑0 An-2〈w°,ψn(s)〉e-λ
2τ
n
ψn(s)
(86)
where the constant φ° has been normalized with the normalization factor An which is given by
An2 )
∫01r(s) ψn2(s) ds
(87)
and
1 1 Pe2 Pe λn M+ , , )0 4 2 4 2 2
∞
(85)
Equation 85 leads to the general solution of eq 64 by invoking the transient solution φn(τ) to yield
which upon substitution of eq 79 yields the characteristic equation, in terms of Kummer’s function as 2
∑0 An-2〈w°,ψn(s)〉φn(τ) ψn(s)
χn)(λnF2)
〈W°,ψn(s)〉 )
∫01W°(s) r(s) ψn(s) ds
(88)
with r(s) being the weighting function defined in a section above. 6. Homogeneous Reaction Case
(82)
and
ai ) a(a + 1)(a + 2)...(a + i - 1)
(83)
bi ) b(b + 1)(b + 2)...(b + i - 1)
(84)
The eigenvalues can be found by solving for the roots of the characteristic equation by using numerical root finders such as the subroutine ZEROIN from the AT&T research lab (Netlib). Figure 5 shows the behavior of the eigenfunctions ψn(s) with respect to the axial coordinate, s, for five different values of λn. It shows that the symmetric boundary and Dirichlet boundary condition are satisfied at the respective ends (i.e., at s ) 0 and s ) 1, respectively). The dimensionless concentration profile is given by the Fourier series expansion of the initial condition using the function φn(τ) and ψn(s) and the linear superposition principle to yield
The central focus of this section is the study of the case of a homogeneous irreversible reaction with kinetics of first order (for this case, the same spatial averaging approach used in the no reaction case). The molar species continuity equation is given by
[
]
∂cA ∂cA ∂2cA ∂2cA + vx )D + 2 - kcA ∂t ∂x ∂z2 ∂x
(89)
The deviation equation for the concentration field, after the order of magnitude analysis has been performed, is given by
∂2cˆ A ∂〈cA〉 xuˆ ) D 2 - kcˆ A ∂x ∂z
(90)
Following the approach of Paine, Carbonell, and Whitaker (Paine et al., 1983), the proposed solution to eq 89 is given by
cˆ A ) g(z)
d〈cA〉 x dx
(91)
1020
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996
Substituting the proposed solution in the deviation equation gives
Dg′′ - kg(z) - uˆ ) 0
Table 1. First Eight Eigenvalues of the Stu 1 rm-Liouville Problem for the Reaction Case Pe ) 2, Da ) 2, Da ) 5
(92)
The boundary conditions for the deviation field can be obtained from the original boundary conditions which are given by no-flux conditions at the plates, and they lead to
∂cˆ A ∂z
)0
at z ) H and z ) -H
at z ) H and z ) -H
(94)
As before (case without reaction), an additional constraint, 〈g(z)〉 ) 0, is required to obtain the solution to eq 91. Effective transport parameters: The derivation of the general transport effective coefficients can be obtained in a way similar to the one used in section 4. A preliminary idea about the behavior of the area average model (i.e., averaged equations in the z variable with effective transport coefficients) can be obtained by assuming that the magnitude of the reaction rate constant, k, multiplied by the field g(z) (see eq 92) is relatively small compared to the other terms and, therefore, the same function g(z) derived in section 4 can be used for analyzing the case where the system undergoes a chemical reaction of first order. The assumption just mentioned will lead to the conclusion that the effective diffusivity, Deff, and the effective convective velocity, Veff, are the same as those for the case without a chemical reaction. The average equations, however, will feature an effective source term given by k〈cA〉 to yield
∂〈cA〉 ∂〈cA〉 ∂2cA + Veff(x) ) Deff(x) 2 - k〈cA〉 ∂t ∂x ∂x
(95)
Equation 95 is a transient parabolic differential equation whose solution will be considered in the next section of the paper. 6.1. Solution to the Concentration Profile. The averaged unsteady molar species continuity equation for the case of convective-diffusive transport coupled with a chemical reaction was derived in the previous section, and it is given by eq 95. By using the same nondimensional variables and constants defined in the previous section, the dimensionless version of the averaged equation (see eq 95) is then given by
[ ]
∂〈cA〉 ∂2〈c〉 ∂〈c〉 ) - Da〈c〉 - s Pe 2 ∂τ ∂s ∂s
Da ) 2
Da ) 5
1 2 3 4 5 6 7 8
3.40995 9.51818 15.76421 22.03136 28.30562 34.58312 40.86237 47.14267
3.82462 9.67449 15.85908 22.09934 28.35857 34.62647 40.89906 47.17448
[
(93)
This conditions lead to the following equations for the variable g(z) at the plates of the channel
g′(z) ) 0
n
]
1 ψn(s) ) exp - s2(Pe2 - Pe) × 4 2 Da 1 1 s2Pe2 Pe λn M+ + , , (98) 4 2 2 4 2 2
(
)
where the parameter Da is the Damko¨hler number given by
Da ≡ kL2/D
(99)
The coefficient an in Kummer’s function M(an, bn, z) has an additional term, Da (i.e., the Damko¨hler number), which shifts in eigenvalues. The characteristic equation for the reaction case can be obtained by evaluating ψn(1) ) 0 which yields
[
M-
]
2
Da 1 1 Pe2 Pe λn + + , , )0 4 2 2 4 2 2
(100)
Equation 100 reduces to the case of no reaction (see eq 81) by taking Da ) 0. The first eight eigenvalues or, alternatively, the roots of eq 100 are shown in Table 1 for the case of Pe ) 2 and Da ) 2 and Da ) 5. Figures 6 and 7 show that the Dirichlet boundary condition at the end and the symmetry condition at the center is satisfied by the eigenfunctions for different Pe numbers. Also, for higher eigenvalues, more oscillations are observed in the eigenfunction, as is expected. The dimensionless concentration profile is given by the Fourier series expansion of the initial condition and the eigenfunctions by the following equation: ∞
W(s,τ) )
∑0 An-2[〈W°,ψn(s)〉φn(τ) ψn(s) + 〈1,ψn(s)〉Znψn(s) (φn(τ) - 1)] (101)
where
Zn )
Da Da + λn2
(102)
An is the normalization factor
An2 )
(96)
∫01r(s) ψn2(s) ds
(103)
and which can also be written in terms of the variable W, defined in the previous section, to give
[ ] 2
∂W ∂W ∂W ) - Da(W + 1) - s Pe ∂τ ∂s ∂s2
(97)
Using the method of separation of variables, as used previously for the case without reaction, the eigenfunctions for this case are given by
〈W°,ψn(s)〉 )
∫01W°(s) r(s) ψn(s) ds
(104)
where r(s) is the weighting function (Weinberger, 1965) given by the same equation identified in the previous section. The concentration profiles for the squeezing flow problem with a homogeneous first-order reaction are shown in Figures 8 and 9. Calculations in Figure 8 were
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1021
Figure 6. First five eigenfunctions for the Dirichlet boundary condition (Pe ) 2, Da ) 2). Figure 9. Transient concentration profiles for the case with reaction (Pe ) 0.5, Da ) 5).
Figure 7. First five eigenfunctions for the Dirichlet boundary condition (Pe ) 5, Da ) 2).
for τ ) 0.05 the concentration at s ) 0 is higher than the value 1.5; however, at the same position such values decrease to a value close to 1.25 for Da ) 5 (see Figure 5). Also, for τ > 0.10 calculations performed (Nagarkar, 1994) have shown that the concentration profile does not change appreciably, which indicates that the steady state of the system has been reached. The particular concentration profile shows values of concentration inside the system that are smaller than the concentration of the environment. This implies that there is some solute entering the system which is balanced by the convective transport and reaction processes in the system. 7. Concluding Remarks
Figure 8. Transient concentration profiles for the case with reaction (Pe ) 0.5, Da ) 2).
performed for the case of Pe ) 0.5 and Da ) 2, and they were based on the transient solution given by eq 101. From the figure, it is seen that the symmetry condition at the position s ) 0 is satisfied for the values of the time used in the analysis. Also, the Dirichlet boundary condition at the external edge of the system (i.e., s ) 1) is also satisfied. The concentration decreases in the system for all the time values shown in the figure. This is consistent with the fact that the reaction is irreversible and first order. The qualitative situation of Figure 9 for Pe ) 0.5 and Da ) 5 is similar to that of Figure 8. The shifting effect of the Da number on the eigenvalues of the system is clearly seen in the Figure. For example,
In the analysis reported in this work the method of spatial averaging was applied to study convectivediffusive mass transfer in a squeezing flow process between two flat plates with and without homogeneous reaction. The resulting effective transport parameters were identified and evaluated. Furthermore, the area average transport equation for the field of concentration of the system was derived. This equation features the effective velocity and the effective diffusivity which have been shown to be functions of the fundamental parameters of the system. Also, the effective velocity is an implicit function of time. This is because the plates move under the effect of an external applied force, F. The effective convective-diffusive problem (with and without reaction) has been solved analytically by the method of separation of variables. An interesting aspect of the analysis presented in this paper is the fact that the eigenvalue problem has been solved analytically by using cylindrical hyperbolic functions which are related to Kummer’s functions. Although both types of boundary conditions at the external edge of the plates, i.e., the Dirichlet and Danckwerts (von Neumann) boundary conditions, can be used, the problem was solved for the case of Dirichlet boundary conditions. The extension to include Danckwerts boundary conditions is not a difficult task (Nagarkar, 1994). The characteristic equation resulting from the eigenvalue problem has been solved by a numerical subroutine, and this is the only part of the problem which requires a numerical procedure. For the case with a linear homogenous reaction in the bulk, the eigenvalues obtained are higher than those
1022
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996
for the case without reaction. This is due to the shift in the eigenvalues because of the effect of the Damko¨hler number. Also, as the Damko¨hler number, Da, increases, the values of the eigenvalues increase. This implies that for the case of higher Da more conversion occurs and the steady state is reached faster. These observations are consistent with what is expected from our analysis based on physical grounds. The final concentration fields are given in terms of a Fourier series expansion of the initial condition in terms of eigenvalues and the eigenfunctions of the problem. The methodology presented in this paper can be extended to include nonlinear reactions by following the integral-spectral approach (Arce and Locke, 1994) where all the aspects related to the eigenvalue problem presented can be used to produce the expansion of Green’s function. An integral equation which features this Green’s function is then used to compute the concentration field. This analysis will be reported in future communications. Some aspects of interest to be studied in future efforts include the following: (1) A generalized boundary condition for well-mixed end conditions (Locke and Arce, 1993) should replace the Dirichlet and Danckwerts boundary conditions. This will help in developing a general model for the convective-diffusive processes, i.e., the squeezing flow processes. The various other boundary conditions can thus be solved as special cases of this model. (2) Computation of the numerical solution of the problem and comparison with the averaged analytical system. This type of comparison may be useful for further generalization and application of the model to other situations of interest in, for example, polymer processing. (3) A case with a reaction at the wall of the plates. This analysis will involve inclusion of the reaction rate constant in the boundary conditions. This case will give rise to a different eigenvalue problem and hence a different characteristic equation if a linear reaction is included as part of the eigenvalue problem. However, the eigenvalue problem reported in this paper will remain as it is if an integral-spectral approach is used (Arce and Locke, 1994). This case will find interesting applications such as in membrane reactors and the various biochemical separation processes. (4) A case of a reaction with nonlinear kinetics. This analysis will involve the integral-spectral method (Arce and Locke, 1994). In this case, an integral equation which features a Green’s function is used to compute the concentration profile. (5) A study of multicomponent cases. This study will deal with the solution to the molar species continuity equation of each component. This case will find applications in industrial squeezing flow processes where more than one component is involved. Acknowledgment A.S.N. was a holder of a Presidential Graduate Fellowship from Florida State University during the course of this research. P.A. is indebted to Dr. W. Langlois for references related to the hydrodynamics of squeezing flow systems. A preliminary version of this work was presented at the Annual AIChE Meeting, St. Louis, MO, Nov 1993. Early versions of this manuscript were typed by Brian Scott and the latest revisions by Heidy Booeshaghi. Nomenclature An ) normalization factor of Fourier expansions c1 ) molar concentration of species A
cˆ ) molar concentration deviation field D ) molecular diffusivity Da ) Damko¨hler number Deff ) effective diffusivity F ) applied force on plates g(z) ) scalar field as defined in eq 38 H(t) ) distance between the two plates k ) kinetic constant L ) length of plates p ) pressure field Pe ) Peclet number r(s) ) weighting function as defined before eq 61 s ) nondimensional axial coordinates Pe/L t ) time (dimensional) veff ) convective velocity vx, vz ) components of the velocity field v(t) ) velocity of plates w ) concentration variable (nondimensional) x ) axial coordinate (dimensional) z ) lateral coordinate (dimensional) zn ) transportation variable (see eq 102) Greek Letters R ) integration constant β ) integration constant γ ) integration constant ψn(s) ) eigenfunctions of the Stu¨rm-Liouville problem φn(τ) ) time-dependent solution given by eq 67 γn ) eigenvalue of the Stu¨rm-Liouville problem τ ) time variable (nondimensional) µ ) Newtonian viscosity of fluid Superscript ° ) initial condition
Appendix. Solution to the Steady-State Problem without Reaction This appendix is designed to check that the solution to the steady-state problem satisfies the result expected in the case where no reaction occurs. The differential equation describing the steady-state problem in the system is simply
d2〈c〉
- s Pe
2
ds
d〈c〉 )0 ds
(1)
This equation can be obtained by dropping the transient term in eq 62 of section 5 of the paper. This equation can be written in the following fashion:
[
]
d〈c〉 d exp(-s2Pe2) )0 ds ds
(2)
A first and second integration of eq 2 yields
)C (-s2Pe) d〈c〉 ds 2
exp
(3)
1
and
∫
〈c〉 ) C1 exp
(-s2Pe) ds + c 2
2
(4)
Now, by using the transformation s2Pe/s ) -t2 and performing the integration indicated in eq 4, the following result is obtained:
x(2π/Pe)
〈c〉 ) -C1
2
erf(s) + C2
(5)
Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1023
The symmetric boundary condition of the original differential problem is, in terms of the nondimensional variables, given by
d〈c〉 )0 ds
at s ) 0
(6)
By using this boundary condition in eq 5, the following equation can be written:
d〈c〉 x(2π/Pe) 2 ) -C1 H (s) exp(-s2) ds 2 π 0
(7)
where H0(s) is the Hermitian polynomial [see Abramowitz and Stegun, 1972]. By expressing the Hermitian polynomial in terms of Kummer’s functions (Abramowitz and Stegun, 1972), the following result is obtained:
(
)
d〈c〉 x(2π/Pe) 2 1 s2 ) -C1 M 0, , exp(-s2) dt 2 π 2 2
(8)
Now by invoking the boundary condition at s ) 0 and by noting that M(s)0) * 0, the following result is found:
C1 ) 0
and
〈c〉 ) C2
(9)
The Dirichlet boundary condition at s ) 1 is given by
〈c〉 ) c∞
(10)
which produces the result for C2 ) c∞ and, therefore, the steady-state solution satisfies 〈c〉 ) c∞, as is expected from an analysis on physical grounds. Literature Cited Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Functions; Wiley: New York, 1972. Arce, P.; Locke, B. R. Transport and Reaction: An Integral Equation Approach, Mathematical Formulation and Computational Approaches. Chemical Engineering Research Trends; Council of Scientific Research Integration: Velayil Gardens, Trivandrum, India, 1994; Vol. 2, pp 89-158. Arce, P.; Cassano, A. E.; Irazoqui, H. A. The Tubular Reactor with Laminar Flow Regime: An Integral Equation ApproachsI: Homogeneous Reaction with Arbitrary Kinetics. Comput. Chem. Eng. 1988, 12, 1103. Arce, P.; Locke, B. R.; Trigatti, I. M. B. An Integral-Spectral approach for Convective-Diffusive Transport and Reaction in Tubular Fouseuille Flows: Mathematical Formulation, Computational Strategies and Selected Applications. AIChE J., in press. Aris, R. On the Dispersion of a Solute in a Fluid Flowing Through a Tube. Proc. R. Soc. London 1956, A235, 67. Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; Wiley: New York, 1960. Bird, R. B.; Hassager, O.; Armstrong, R. C. Dynamics of Polymeric Liquids: Vol. I, Fluid Mechanics; John Wiley and Sons: New York, 1987. Broyer, E.; Macosko, C. W. Heat Transfer and Curing in Polymer Reaction Molding. AIChE J. 1976, 22 (No. 2), 268. Coleman, C. J. On the Squeezing Flow of a Fibre-Reinforced Fluid. J. Non-Newtonian Fluid Mech. 1990, 37, 379-385. Danckwerts, P. V. Continuous Flow Systems. Distribution of Residence Times. Chem. Eng. Sci. 1953, 2, 1. Davis, E. J. Exact Solution for a Class of Heat and Mass Transfer Problems. Can. J. Chem. Eng. 1973, 51, 562.
DelGatto, J. V. Boom in Castable Polyurethane. Rubber World 1973, 169, 43. Denn, M. M. Process Fluid Mechanics; Prentice-Hall: Englewood Cliffs, NJ, 1980. Dienes, G. J. J. Colloid Sci. 1947, 2, 131. Dienes, G. J.; Klemm, H. F. J. Appl. Phys. 1946, 17, 458. Gent, A. N. Br. J. Appl. Phys. 1960, 11, 85. Gray, W. G. A Derivation of the Equations of Multiphase Transport. Chem. Eng. Sci. 1975, 30, 229. Gross, W. A. Unsteady Bearing Films and Bearing Systems. IBM Research Report RJ-117-8; 1960. Haberman, R. Elementary Applied Partial Differential Equations; Prentice-Hall: Englewood Cliffs, NJ, 1987. Howes, F. A.; Whitaker, S. Spatial Averaging Theorem Revisited. Chem. Eng. Sci. 1985, 40, 1387. Langlois, W. E. Isothermal Squeeze Films. Q. Appl. Math. 1962, 20 (2), 131-150. Lauwerier, H. A. Poiseuille Functions. Appl. Sci. Res. 1951, A3, 58. Locke, B. R.; Arce, P. Applications of Self-Adjoint Operators to Electrophoretic Transport, Enzyme Reactions, and Microwave Heating Problems in Composite MediasI: General Formulation. Chem. Eng. Sci. 1993, 48, 1675. Metzner, A. B. Trans. Soc. Rheol. 1968a, 12, 57. Metzner, A. B. J. Lubr. Technol. 1968b, 90, 531. Mooney, M. In Rheology; Eirich, F. R., Ed.; Academic Press: New York, 1958; Vol. 2, Chapter 5, pp 195-196. Nagarkar, A. Diffusive-Convective Transport with Chemical Reaction in Squeezing Flows. MS Thesis, The Florida State University, Tallahassee, FL, 1994. Oka, S. In Rheology; Eirich, F. R., Ed.; Academic Press: New York, 1960; Vol. 3, Chapter 2, pp 43-75. Paine, M. A.; Carbonell, R. G.; Whitaker, S. Dispersion in Pulsed SystemssI: Heterogeneous Reactions and Reversible Adsorption in Capillary Tubes. Chem. Eng. Sci. 1983, 38, 1781. Phan-Thein, N.; Graham, A. L. The Squeezing Flow of a Model Suspension Fluid. Rheol. Acta 1990, 29, 433. Ramkrishna, D.; Amundson, N. R. Linear Operator Methods in Chemical Engineering; Prentice-Hall: Englewood Cliffs, NJ, 1985. Scott, J. R. Trans., Inst. Rubber Ind. 1931, 7, 169. Scott, J. R. Trans., Inst. Rubber Ind. 1932, 8, 481. Skartsis, L.; Khomani, B.; Kardos, S. L. Resin Flow Through Fiber Beds During Composite Manufacturing Process. Part I. Review of Newtonian Flow Through Fiber Beds. Polym. Eng. Sci. 1992, 32 (4), 221. Slater, L. J. Confluent Hypergeometric Functions; Cambridge University Press: Cambridge, England, 1960. Slattery, J. C. Flow of Viscoelastic Fluids Through Porous Media. AIChE J. 1967, 13, 1066-1071. Slattery, J. C. Momentum, Energy, and Mass Transfer in Continua; McGraw-Hill: New York, 1981. Tanner, R. Engineering Rheology; Oxford University Press: New York, 1985. Van Wazer, J. R.; Lyons, J. N.; Kim, K. Y.; Colwell, R. E. Viscosity and Flow Measurement; Wiley: New York, 1963; pp 292-295. Weinberger, M. D. A First Course in Partial Differential Equations; Blaisdell: New York, 1965. Whitaker, S. Diffusion and Dispersion in a Porous Media. AIChE J. 1967, 13, 920-927. Whitaker, S. The Transport Equations for Multiphase Systems. Chem. Eng. Sci. 1973, 28, 139-147. Wood, A. S. RIM New High Speed Processing Option. Mod. Plast. 1947, 51 (No. 6), 54. Zanotti, F.; Carbonell, R. G. Development of Transport Equations for Multiphase SystemssI. Chem. Eng. Sci. 1984, 39 (2), 263.
Received for review April 10, 1995 Revised manuscript received October 26, 1995 Accepted December 12, 1995X IE950235X X Abstract published in Advance ACS Abstracts, February 15, 1996.
