Ind. Eng. Chem. Res. 1997, 36, 357-367
357
MATERIALS AND INTERFACES Multidimensional Modeling of Chemical Vapor Infiltration: Application to Isobaric CVI John Y. Ofori and Stratis V. Sotirchos* Department of Chemical Engineering, University of Rochester, Rochester, New York 14627
Multidimensional (2D and 3D) chemical vapor infiltration models are formulated and used to simulate densification of isotropic or anisotropic preforms of various geometries using SiC deposition from methyltrichlorosilane. A generalized form of the dusty-gas model for mass transport in anisotropic porous structures is used as a flux model in the pore space, and the multidimensional model equations are solved using the Galerkin finite element method. Structures consisting of freely overlapping fibers parallel to a line (one-directional), parallel to a plane (two-directional), or without preferred orientation (three-directional) are employed to model the microstructure of the preforms. The obtained results show that the effects of preform geometry can lead to partial pressure and deposition profiles in the preform that are considerably different from those suggested by one-dimensional analyses of the chemical vapor infiltration process. Introduction Chemical vapor infiltration (CVI) is a method of ceramic or carbon matrix composite fabrication in which chemical vapor deposition reactions are employed to deposit the matrix material (ceramic or carbon) on the internal surface of a porous preform. Compared to other methods of composite fabrication, such as pressing and sintering, CVI has a number of attractive advantages, the most important being its ability to fabricate nearnet-shape components and deposit a variety of matrix materials without damaging the material used for reinforcement (usually fibers) (Naslain, 1992; Warren, 1992). There are a number of different ways in which CVI can be implemented in practice, their main distinguishing features being the primary mass transport mechanism that controls transport in the pore space and the presence or absence of externally imposed temperature gradients inside the preform. The primary mass transport mechanism in isobaric CVI is diffusion, whereas in forced-flow and pulse CVI, transport occurs mainly by viscous (pressure-gradient-driven) flow. Each of these CVI variants can be operated with a thermal gradient that can be generated in a variety of ways, such as surface heating or cooling, radio frequency or microwave heating, and use of mandrels. As with any other physicochemical processes, mathematical models of chemical vapor infiltration are indispensable for the analysis and interpretation of experimental data, the design of new CVI processing schemes, and the optimization of existing processes. Several modeling studies have been presented in the literature (e.g., Gupte and Tsamopoulos, 1989; Stinton et al., 1986; Chung et al., 1992; Naslain et al., 1989; Starr and Smith, 1990; Lin and Burggraaf, 1991; Currier, 1990; Fedou et al., 1993a,b; Sheldon and Chang, 1993). Most of these models use simple mass transport models in the preform for the diffusion process (usually of the form of Fick’s law), which neglect the * To whom correspondence should be addressed. S0888-5885(96)00487-3 CCC: $14.00
interaction of the fluxes of the reactants and products in the preform. Moreover, they assume one-dimensional spatial flux dependence for the partial pressure profiles of the gases and the conversion of the pore space. With regard to the latter assumption, the only exception is a two-dimensional mathematical model by Starr and Smith (1990), in which, however, a simplified treatment of the mass transport problem in the pore space is also employed. A general mathematical model accounting for multicomponent interactions during chemical vapor infiltration was formulated by Sotirchos (1991), and its versatility was demonstrated by applying it to the analysis of isobaric CVI and pulse CVI. On the basis of this model, we developed rigorous mathematical models for various processing schemes (Ofori and Sotirchos, 1996a,b) which we employed to investigate the effects of various operating conditions and preform characteristics on the relationship between processing time and deposition uniformity. The mathematical model of Sotirchos (1991) is, in principle, a multidimensional model, but the transport model that it employs for describing the interaction of fluxes, partial pressures, and partial pressure gradients in the interior of the preform, namely, the threeparameter dusty-gas model (Jackson, 1977; Mason and Malinauskas, 1983; Sotirchos, 1989), is rigorously valid for multidimensional transport in isotropic media. For anisotropic pore structures, as is the case with several preforms in CVI, it can be used only if the principal coordinate systems of the effective transport coefficient tensors for bulk diffusion, Knudsen diffusion, and viscous flow coincide and only for unidirectional transport in the direction of one of the principal axes. For this reason, our past investigations of chemical vapor infiltration processes have been limited to one-dimensional transport in the preform. However, the assumption of a dominant direction for mass transport in the preformswhich must hold for a one-dimensional analysis to be employedscannot be justified for preforms having aspect ratios close to unity or complex shape. © 1997 American Chemical Society
358 Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997
Satisfactory description of the evolution of the density profile in the interior of such preforms requires that allowance be made for multidimensional spatial dependence of the partial pressures of the gaseous reactants and products in the pore space. Multidimensional (two-dimensional (2D) and threedimensional (3D)) chemical vapor infiltration models are formulated in the present study for anisotropic, in general, preforms of arbitrary pore structure. Tensor instead of scalar forms are employed for the parameters of the dusty-gas model in order to render it capable of describing multidimensional transport in anisotropic porous media. To the best of our knowledge, this is the first time that a rigorous multidimensional model for multicomponent mass transport (by bulk diffusion, Knudsen diffusion, and viscous flow), energy transport, reaction, and structure evolution in anisotropic media is presented. The model equations are solved using the Galerkin finite element method with quadratic basis functions. The effects of preform anisotropy, shape, and aspect ratio are investigated by applying the multidimensional (2D and 3D) models to chemical vapor infiltration of SiC through methyltrichlorosilane decomposition in preforms of various shapes that have microstructure that can be represented by one-directional (parallel to a line), two-directional (parallel to a plane), or three-directional (randomly oriented), freely overlapping fibers. Multidimensional Mathematical Model for CVI The general CVI system consists of a porous preform exposed to a multicomponent mixture of gases, with a gaseous reactant undergoing decomposition to deposit a solid product on the gas-solid interface. The multidimensional mathematical model equations for CVI are similar to those given by us in earlier studies (Sotirchos, 1991; Ofori and Sotirchos, 1996a,b), the only differences originating from the consideration of anisotropic porous media in the present study. They consist of mass balances for the gaseous species present in the reaction mixture, a mass balance for the solid material that is deposited, and an energy balance. A multicomponent mass transport flux model is required to relate the fluxes of the species with their partial pressures (or concentrations) and partial pressure gradients. A key assumption in the formulation of the mathematical model is that the porous medium can be treated as a continuum; that is, in the least statistically representative volume of the porous medium, the partial pressures and temperature can be assumed to be uniform. Transport and Reaction Equations. The mass balance for the ith gaseous species is written for the general case of a system with n gaseous species and m chemical reactions as
where RvF is the rate of the Fth homogeneous reaction (per unit volume), RsF is the rate of the Fth heterogeneous reaction (per unit area), and Se is the accessible internal surface area (per unit volume of porous medium). The mass balance for the solid produced by the deposition reaction is given in terms of the fractional filling of the pore space by
) ∂t
∂( pi/RT) ∂t
+ ∇‚Ni )
∑F νiFRvF
Cep
∂T
RvF ) SeRsF
ξ)
0
∂t
) ∇‚(λe∇T) -
0 -
(3a,b)
0
∑i NiCpi‚∇T + ∑F (-∆HF)RvF + Q˙
(4)
where λe is the effective thermal conductivity tensor of the composite, Cep and Cpi are the heat capacity per unit volume of the composite and the molar heat capacity of the ith gas species, respectively, ∆HF is the heat of reaction of the Fth reaction, and Q˙ is the energy input to the system from an external source, such as radio frequency or microwave heating. If surface heating were to be employed, its effect would be represented in the mathematical model by a term appearing in the boundary conditions for the energy equation. The effective thermal conductivity can be determined using data given by Tomadakis and Sotirchos (1993a) for bulk diffusion in randomly overlapping arrays of conductive cylinders (the solid phase) dispersed in a nonconducting medium (the gas phase), if the effect of the gas phase is neglected, or by Tomadakis and Sotirchos (1996) for conduction in structures of conductive cylinders and a conductive matrix. For isotropic media the effective medium theory (Kirkpatrick, 1973) can be employed to estimate the effective thermal conductivity from the conductivities of the gas and the solid. Dusty-Gas Model Equations. The dusty-gas model equations (Mason and Malinauskas, 1983; Jackson, 1977; Sotirchos, 1989) were used in our past studies to describe the coupling of the fluxes, partial pressures, and partial pressure gradients of the gas species in the void space of the porous medium. For an isotropic porous medium, the equations for the diffusive fluxes, ND i , are written as n
1 RT
∇pi )
∑
i)1,j*i
D xjND i - xiNj
ND i +
Deij
DeKi
(5)
(1)
where pi and Ni are the partial pressure and molar flux of the ith gaseous species, respectively, e is the accessible porosity, νiF is the stoichiometric coefficient of the ith species in the Fth reaction, R is the ideal gas constant, and T is the temperature. The reaction rate RvF (per unit volume of porous medium) of the Fth reaction is given by either of the equations
RvF ) eRF;
∑νSFRvF; F
where vS is the molar volume of the solid, 0 is the initial porosity of the preform, νSF is the stoichiometric coefficient of the solid species in the Fth reaction, and ξ is the fractional conversion of pore space to solid. The energy balance is given by
e
vS
∂ξ
(2a,b)
where xi is the mole fraction of the ith species, and ND i and ∇pi are its diffusive flux and partial pressure gradient in the direction of interest. The effective binary diffusivity Dije for the (i, j) gas pair and the e effective Knudsen diffusivity DKi for the ith gas species are given by the expressions
(p*p )(T*T )
Deij ) S1D* ij
3/2
;
(T*T )
DeKi ) S2D*Ki
1/2
(6a,b)
where D*ij is the binary diffusion coefficient for the (i, j)
Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997 359
pair computed at reference pressure and temperature p* and T*, and D*Ki is the Knudsen diffusivity of the ith species in a capillary of reference radius r* at reference temperature T*. S1 and S2 are structural parameters used to express the effects of the structure of the porous medium on the effective binary and Knudsen diffusivities, respectively. The total flux of the ith species is obtained by adding to the diffusive flux the viscous flow contribution: V Ni ) ND i + Ni
(7)
1
(S1)x
RT*
Bepi ∇p; µRT
(T*T )
µ ) µ*
1/2
Be is the effective permeability of the porous medium, µ is the viscosity of the gaseous mixture, and µ* is the viscosity of the gaseous mixture at the reference temperature T*. The dusty-gas model that is described by eq 5 is applicable for multidimensional transport only to isotropic porous structures. For anisotropic media, the effective transport coefficients for bulk diffusion, Knudsen diffusion, and viscous flow are of tensorial form. Equations 5 and 8a may be applied to an anisotropic medium only if the principal axes of the three transport coefficient tensors coincide and, even then, only for unidirectional transport along one of the principal axes. When this is the case, the effective transport coefficients e that appear in eq 5 and 8a, Dije , DKi , and Be, stand for transport parameters along the principal axes that coincide with the direction of transport. This is exactly what was done in our past studies (Ofori and Sotirchos, 1996a-c), where the multidirectional fibrous or capillary structures we employed had the same principal coordinate systems for bulk diffusion, Knudsen diffusion, and viscous flow. The extension of the three-parameter dusty-gas model to multidimensional transport in anisotropic porous structures leads to the following modified form for eqs 5 and 8a (Sotirchos, 1996):
S1 ∇pi
RT*
( ) T
T*
1/2
)
n
D pjND i - piNj
i)1,j*i
D* ijp*
∑
Bepi )-
D pj(ND i )x - pi(Nj )x
i)1,j*i
D* ijp*
)
µRT
∇p
ND i T -1 D*Ki T*
∑
(10)
with the total flux still given by eq 7. S1 is the structural parameter tensor for bulk diffusion, S2 is the corresponding tensor for Knudsen diffusion, and Be is the effective permeability tensor for viscous flow. Using S1 and S2 in place of S1 and S2 in eqs 6a and 6b gives the effective diffusivity tensors for bulk diffusion e and Knudsen diffusion, respectively, Dije and DKi . If the principal axes of the tensors S1 and S2 coincide with each other, eq 9 can be written for the diffusive fluxes in the x-direction as
(11)
where (Si)x is the xx element of the Si tensor. If Be also has the same principal axis as S1 and S2, then the viscous fluxes in the x-direction are given by:
(Be)xpi T -3/2 ∂p )µ*RT* T* ∂x
( )
(12)
Analogous equations can be written for the other axes for both the diffusive and viscous fluxes. It should be noted that most cases encountered in practice fall into the class of materials in which the principal axes for the different transport mechanisms are the same. Structure Evolution Models. Statistical structural models are usually based on the use of either hollow objects to represent the void spaces or dense objects to describe the solid phase of the preform. In this study we employ structures of randomly overlapping fibers with various directionalities one-directional, two-directional, and three-directional) to describe the structural preform. The following equations are used to describe the variation of S1, S2, and Be with conversion:
(S1)j )
; ηbj
(S2)j )
rj ; ηKj r*
(Be)j )
rj2 (13a,b,c) 8ηpj
with
2 S
rj )
(13d)
where ηbj , ηKj , and ηpj are the tortuosity factors for bulk diffusion, Knudsen diffusion, and viscous transport, respectively, in the j-direction, and rj is the average pore radius of the fiber structure. ηbj and ηKj are obtained from the correlations developed by Tomadakis and Sotirchos (1991a, 1993a) using Monte Carlo estimates of effective diffusivities in structures of uniform size fibers:
(
ηj ) ηj|)m (9)
+
(D* Ki) (S2)x T*
+
(S1)(S2)
NVi
∂x T*
(NVi )x
(8a,b)
n
1/2
(ND i )x (S1)x T
The viscous flux, NVi , is found using D’Arcy’s law:
NVi ) -
( )
∂pi T
)
m - p - p
R
(14)
The values of m, ηj|)m, and R for the various directions of diffusion in each structure are given in Table 1. p is the percolation threshold of the structure in the respective direction of diffusion. The tortuosity factor for viscous permeability, ηpj , is assumed to be equal to that for bulk diffusion (ηbj ). The surface area per unit volume of the preform is given for structures of randomly overlapping fibers by
S)-
2 ln(); r
r ) r0
( ) ln() ln(0)
1/2
(15)
where r is the radius of the fibers. The computations of Burganos and Sotirchos (1989) and Tomadakis and Sotirchos (1991b) showed that, in random fiber structures, formation of inaccessible space occurs only in a narrow range of porosities next to the percolation
360 Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997 Table 1. Transport Parameters for Fiber Structures (Tomadakis and Sotirchos, 1991, 1993)a Knudsen regime
a
fibrous structure and direction of diffusion
p
ηj|m)1
R
1D (parallel) 1D (perpendicular) 2D (parallel)
0 0.33 0.11
0.541 1.747 1.149
0 1.099 0.954
bulk regime, ηj|m)1 ) 1, R
2D (perpendicular) 3D
0.11 0.037
1.780 1.444
1.005 0.921
0 0.707 0.521 ( > 0.4) 0.872 ( < 0.4, m ) 0.4) 0.785 0.661 ( > 0.4)
The tortuosity factor for viscous flow is taken equal to that for bulk diffusion.
threshold. Therefore, for all practical purposes, e and Se can be set equal to and S, respectively, for porosities above the percolation threshold and to zero below it. We use randomly overlapping fibers to represent the structure of the porous medium because we want to make the results of the present study directly comparable to those of past studies (e.g., Ofori and Sotirchos, 1996a-c), in which randomly overlapping fiber structures were employed. A more realistic representation of the preform structure can be obtained using structures of initially nonoverlapping fibers, uniformly distributed in space according to some orientation distribution or organized into bundles (tows or threads). Simulation results for the structural properties and the effective diffusivities, as well as correlations of the form of eq 14, are available for these structures in the literature (Tomadakis and Sotirchos, 1991b, 1993b; Transvalidou and Sotirchos, 1996; Transvalidou, 1996). These results can be used in the general model we formulated in the preceding sections to study the densification behavior of nonoverlapping fiber structures. However, the results obtained in this study are applicable, qualitatively and to a certain extent quantitatively, to structures of nonoverlapping fibers or fiber bundles since the dependence of the structural and transport properties of such media on the extent of densification exhibits the same qualitative characteristics as that of overlapping fiber structures. Application to Densification by SiC Deposition Reaction and Kinetics. We apply the transport and reaction model to the case where the matrix material is SiC and is deposited through decomposition of methyltrichlorosilane (MTS). The deposition is assumed to proceed according to the overall reaction H2
CH3SiCl3 98 SiC + 3HCl; 0 6 ∆H1300 K ) 220 × 10 J/kmol (16a,b)
The rate expression for the reaction is assumed to have the form
Rs ) ks0 exp(-E/RT)pMTS
(16c)
with ks0 ) 2.290 × 10-2 kmol/(m2‚s‚atm) and E ) 120 × 106 J/kmol (Brennfleck et al., 1984). Boundary and Initial Conditions. The CVI model that was formulated in the previous section only requires an appropriate set of boundary conditions in order to describe a particular process. For the isobaric CVI application we examine here, the composition of the mixture in the gas phase surrounding the preform is assumed to be the same everywhere, with negligible mass transport limitations between the gas phase and the external surface of the preform. The composition of the gas mixture at the external surface of the preform
is therefore known and equal to that in the bulk phase, i.e.,
pi|external surface ) pib,
i ) 1, ..., n
(17a)
In cases where the densification problem inside the preform has a symmetry plane, it is not necessary to model the whole preform, but only a part of it. For this part of the preform, the standard no-flux (symmetry) boundary conditions are applied at the symmetry plane of the slab. Specifically, we have
n‚Ni ) 0
(17b)
where n is the normal vector to the symmetry plane. Of course, eq 17b can also be applied to a sealed surface. The boundary condition for the energy equation is that the temperature at the surface of the preform is fixed, i.e.
T ) Ts
(18a)
whereas at a symmetry plane or insulated boundary
n‚∇T ) 0
(18b)
The initial conversion is set equal to zero. Since we are interested in the long-time (pseudo-steady-state) behavior of the process, it is immaterial what initial conditions are used for the partial pressures and the temperature. Computational Scheme for Multidimensional CVI Application of the CVI model equations given here to describe the densification process requires that the total flux of the gaseous species, Ni, be obtained from the multidimensional multicomponent mass transport model. When the Galerkin finite element method is employed, it is not necessary to obtain the divergence of the fluxes, ∇‚Ni, appearing in eq 1 because only Ni and the spatial derivatives of the basis functions appear in the final equations. This reduces significantly the computational time in comparison to other numerical methods, such as finite differences, that requires the computation of ∇‚Ni. The anisotropic dusty-gas model equations (eq 11) can be written in matrix form for the x-direction as
-
1 C ) (B)x(ND)x RT* x
(19)
where Cx is an n-dimensional vector with the following elements
(T*T )
(Ci)x ) (S1)x
1/2
∂pi ∂x
(20)
Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997 361
and Bx is the n × n matrix:
(Bij)x )
-pi
; D*ijp*
j * i; n
(Bii)x )
pj
1 (S1)x T
+ ∑ D* i-1,j*iD*p* ij
Ki
(S2)x T*
(21)
Given the values of T, p, ∂p/∂x, and ξ, we can find (ND i )x for all n gaseous species by inverting eq 19. This procedure is repeated for the computation of (ND i )y and (ND ) . The total diffusive flux of the ith species is then i z constructed as Figure 1. Schematic diagram of a three-dimensional preform. D D D T ND i ) [(Ni )x, (Ni )y, (Ni )z]
(22)
The viscous flux in the x-direction is found from eq 12, and those in the other directions have the corresponding equations. Then the total flux of each of the n gaseous species present in the system is assembled using eq 7. The CVI model equations are discretized in their 2D or 3D forms employing computer codes that were formulated on the basis of the Galerkin finite element method. The underlying theory is given in a number of texts (Strang and Fix, 1974; Finlayson, 1980; Mori, 1986; Rao, 1989; Zienkiewicz and Taylor, 1989a,b). Quadratic Lagrangian basis functions are employed in all cases to represent the partial pressure, temperature, and conversion profiles inside the preform. In the 2D case the domain of the problem is divided into triangular elements, while in the 3D case tetrahedral elements were used. The implicit set of ordinary differential equations that is generated by the discretization is solved using an implicit algebraic-differential equation solver. Multiple integrals are evaluated using quadrature formulas given by Stroud (1971), Zienkiewicz and Taylor (1989a), and Strang and Fix (1974). The transport properties required in each direction for anisotropic structures are calculated independently and stored in a database, from which values at any conversion level are obtained by interpolation. Further details of the use of the Galerkin finite element method in solving the multidimensional CVI model equations are given by Ofori (1996). Results and Discussion We employ the 2D and 3D CVI models described earlier to simulate densification by isobaric CVI in three types of fiber structures, namely, one-directional (fibers parallel to a line), two-directional (fibers parallel to a plane), and three-directional (fibers oriented randomly). The first two structures are anisotropic, but the principal coordinate systems of the three mass transport tensors coincide. Unless otherwise specified, the direction with the largest mass transport coefficient is generally taken to be the x-direction, while the direction with the smallest coefficient is chosen as the z-direction. The tortuosity factors in the x-, y-, and z-directions for the fiber structures are presented in Table 1. In all of our calculations we use an initial fiver size of 20 µm and an initial porosity of 0.5, and the reaction conditions are fixed at 1300 K, 1 atm, and 10% MTS in H2. The effects of operating conditions on the isobaric CVI process have been examined in past studies (Ofori and Sotirchos, 1996a-c), and these papers may be consulted for further details. The dimensions of the preform are
given in each figure. The Chapman-Enskog equations (Bird et al., 1960) were used to calculate the binary diffusion coefficients and viscosities, and the Wilke equation (Bird et al., 1960) was used to estimate the gas mixture viscosity. The modified Eucken equation (Svehla, 1962) was used to calculate the gas thermal conductivities from which the gas mixture thermal conductivity was estimated using the Wilke equation. The thermal conductivity for SiC given by Fieldhouse et al. (1958) was employed. The dimensionless time is defined in such a way that it would require 10 units of it for a pore of radius 20 µm to be filled with solid product at the conditions prevailing of the external surface of the preform, that is
τ)
10vSiCks0 exp(-E/RTb)pMTS,b t r0
(23)
Effects of Preform Anisotropy. We begin by examining the effect of preform anisotropy on the progress of densification. We consider a Cartesian coordinate system placed at the center of the preform with each axis perpendicular to a pair of faces. Figure 2 shows the pseudo-steady-state MTS partial pressure profiles that are obtained along the three axes in each of the three structures at the beginning of densification ( ) 0.5) for a preform (shown in Figure 1) with equal dimensions in the x-, y-, and z-directions. The magnitude of the differences in the MTS profiles in different directions for the one-directional and two-directional structures is a reflection of the degree of anisotropy. From a comparison of the MTS profiles in the different directions in Figure 2, the one-directional structure appears to be more anisotropic than the two-directional structure. The three-directional structure is isotropic, and, thus, the MTS profile is identical along all three axes. The partial pressure in the x-direction in the onedirectional structure is higher than that in the threedirectional structure, while the partial pressures in the y- and z-directions in the one-directional structure are lower than those in the three-directional structure everywhere except close to the center of the preform. Figures 3 and 4 show the distribution of the MTS partial pressure on the z ) 0 plane for preforms composed of three-directional and one-directional fiber structures, respectively. Because of the isotropy of the threedirectional structure, the spacing of the concentration contours (a measure of the concentration gradient) in Figure 3 is the same along the two axes. This is not the case in Figure 4, where the concentration field is
362 Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997
Figure 2. Effect of preform anisotropy on MTS partial pressure profiles. pb ) 1 atm; xMTS,b ) 0.1; xHCl,b ) 0; xH2,b ) 0.9; Tb ) 1300 K; r0 ) 20 µm; 0 ) 0.5.
Figure 3. Contour plot of MTS partial pressure distribution in the preform shown in Figure 1 at the plane z ) 0 for a threedirectional fiber structure. ax ) ay ) ax ) 0.01 m; ) 0.5. Conditions the same as in Figure 2.
“stretched” in the y-direction, the direction of lower effective diffusivity. At the reaction conditions of Figures 2-4, the diffusion process is in the bulk regime, and, therefore, bulk diffusion is the principal mechanism of transport. Figure 5 shows the variation of the structural parameter for bulk diffusion in each direction for the three fiber structures. The orientationally averaged value for each structure, given by the expression
1 S1 ) ((S1)x + (S1)y + (S1)x) 3
(24)
is also shown. In agreement with the results of Figure 2, it is seen that the one-directional structure is much more strongly anisotropic than the two-directional structure. It is interesting to note that the partial pressure of MTS at the center of the preform is slightly
Figure 4. Contour plot of MTS partial pressure distribution in the preform shown in Figure 1 at the plane z ) 0 for a onedirectional fiber structure. ax ) ay ) az ) 0.01 m; ) 0.5. Conditions the same as in Figure 2.
Figure 5. Variation of directional and volume-averaged dimensionless bulk diffusivities.
higher in the one-directional structure than in the isotropic three-directional structure, even though the orientationally averaged effective bulk diffusivity of the former structure is smaller. It should be noted, though, that the average partial pressure in the preform was found to vary among the three structures in the same way as the orientationally averaged effective bulk diffusivity. Since the relative differences among the bulk diffusivity structural factors for the various structures and transport directions examined in Figure 5 vary with increasing conversion, the relative differences among the conversion profiles at conversions different from zero may be larger than those seen in Figure 2 among the MTS concentration profiles. Figure 6 presents the densification profiles along the three axes in the three fiber structures when the surface conversion has reached two different stages, ξs ) 0.32 (s ) 0.34) and ξs ) 0.58 (s ) 0.21). The conversion profile for the two-direc-
Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997 363
Figure 6. Effect of preform anisotropy on conversion profiles in the 3D model. Conditions the same as in Figure 2.
tional structure in the x-direction is not presented for ξs ) 0.32 since it almost coincides with that for the three-directional structure. The curves for ξs ) 0.32 are left unlabeled because they are in the same positions relative to each other as those for ξs ) 0.58. In agreement with Figure 2, the conversion at the center of the preform is higher in the one-directional structure. The differences among the conversion profiles of the two-directional and three-directional structures are relatively small and similar to those observed among the MTS concentration profiles at zero conversion in Figure 2. However, since the effective diffusivity in the three-directional structure becomes larger than both the z and x (or y) effective diffusivities in the two-directional structure (see Figure 5), the differences between the two structures tend to increase with increasing conversion. As Figure 5 shows, the one-directional structure reaches its percolation threshold for transport perpendicularly to the fibers at about 33% porosity, which corresponds to ξ ) 0.34. When this happens, the sides of the preform that are parallel to the fibers (i.e., the x-axis) are effectively sealed, and the concentration in the preform at the surfaces cannot remain equal to the concentration in the bulk. As a result, the conversion on the y- and z-axes at the external surface starts to deviate from that on the x-axis (see Figure 6). Effects of Preform Aspect Ratio. In this section, we examine the influence of the aspect ratio of the preform on the predictions of 2-D and 3-D CVI models under isobaric conditions. Only results for the onedirectional and three-directional fiber structures are used in the comparison since two-directional fiber structures are weakly anisotropic and thus give results close to those of the 3D fiber structure (see Figures 2 and 6). Figure 7 presents the pseudo-steady-state MTS partial pressure profiles obtained along the three axes at the beginning of the densification process ( ) 0.5) for three aspect ratios (ax/az, with ay and az fixed) using the 3D model. The corresponding results for the 2D model on the (x, y) plane are shown in Figure 8. The 2D results can also be obtained from the 3D model if the sides at z ) (az are treated as sealed. Figures 7 and 8 show that, as the aspect ratio increases, the differences between the concentration
Figure 7. Effect of preform anisotropy and aspect ratio on MTS partial pressure profiles in the 3D model. Conditions the same as in Figure 2.
Figure 8. Effect of preform anisotropy and aspect ratio on MTS partial pressure profiles in the 2D model. Conditions the same as in Figure 2.
profiles along the two axes become larger. Since there is one more direction available for transport in the 3D model, the concentrations that are computed along the x- and y-axes using this model are much higher than those found from the 2D model. For very large aspect ratios, mass transport along the direction corresponding to the largest dimension of the preform influences the concentration profile only in the vicinity of the external surfaces that are normal to this direction. This behavior is seen to occur for aspect ratios equal to 5 and 10 in Figures 7 and 8. These results point to the conclusion that for large enough aspect ratiossas is the case with many of the objects that are densified using chemical vapor infiltrationsvery accurate reactant concentration (and therefore conversion) profiles away from the edges of the preform can be obtained by using mathematical
364 Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997
Figure 9. Variation of MTS partial pressure at the center of the preform with the aspect ratio (2D and 3D models). ) 0.5. Conditions the same as in Figure 2.
models based on one-dimensional transport in the direction of the smallest dimension. A more informative picture on the effects of the aspect ratio on the concentration profile in the interior of the preform is provided by the results of Figure 9. This figure presents the variation of the center to surface MTS partial pressure ratio with the aspect ratio that is predicted by the 2D and 3D models for one-directional and three-directional structures. It is seen that the value of the aspect ratio beyond which the effects of the largest dimension on the center to surface pMTS ratio becomes negligible depends on the magnitude of the effective diffusivity in that direction relative to those of the effective diffusivities in the other directions. For example, for the one-directional structure, (pMTS)c/(pMTS)s remains practically invariant for aspect ratios greater than about 2 as the dimension in the direction exhibiting the smallest diffusivity is increased, whereas aspect ratios greater than 4 are needed for the same situation to occur when the preform is stretched in the direction of the largest diffusivity. When one dimension becomes very large, the transport and reaction problem away from the edges obviously becomes two-dimensional in the other two dimensions. This is the reason for which the limit approached by the center to surface partial pressure ratio in Figure 9 by the three-directional structure in the 3D model is the same as the center to surface partial pressure ratio in Figure 8s(pMTS)s is equal to 0.1sfor the same structure for ax/ay ) 1. Effects of Preform Shape. We use in this section preforms of various shapes to investigate the dependence of the partial pressure and conversion profiles on preform shape. The 2D model is first employed with the isotropic three-directional structure to isolate the influence of the preform shape on the behavior of the process. The 3D model is subsequently applied to both isotropic and anisotropic structures to study the combined effects of shape and anisotropy. Figure 10 presents a contour plot of the conversion field at dimensionless time τ ) 7.5 when an I-shaped 2D domain is densified using isobaric CVI. It is clearly seen in the picture that the spacing of the equal
Figure 10. Effect of preform shape on conversion pattern in an I-shaped cross section. Conditions the same as in Figure 2. Base length ) 2ax ) 0.03 m; height ) 2ay ) 0.03 m. τ ) 7.5; ξs ) 0.76.
conversion curves at a certain depth from the external surface depends strongly on the location in the preform. The spacing is the largest and the penetration of the reactant the deepest at outward pointing corners, whereas the opposite situation prevails in the neighborhood of inward pointing corners. This behavior is due to the availability of much more surface area per unit volume of the preform at outward pointing corners for transport of the reactant to the interior. Locations away from the edges, i.e., the top and bottom sides, are characterized by conversion profiles in the normal direction to the external surface that are very close to those obtained from a 1D model. If the angle of the outward pointing corner becomes smaller, the masstransfer surface area increases, and the reactant penetrates deeper into the structure. This effect may be seen in Figure 11, which presents conversion contours at τ ) 7.5 for a 2D domain having the shape of a parallelogram. The combined effects of preform geometry and anisotropy are examined in Figure 13, which presents pseudo-steady-state MTS concentration contour plots at zero conversion for a preform having the shape shown in Figure 12. It is a cube with two holes cut into it from the two sides that are perpendicular to the x-axis. The preform structure is composed of one-directional fibers with the fibers aligned in the x-direction. The concentration profiles that are shown in Figure 13 refer to the x ) 0 plane, and Figure 14 presents the corresponding concentration (partial pressure) field for the case in which the same preform, but without any holes cut into it, is reacted at the same conditions. Comparison of Figures 13 and 14 shows that the presence of holes alters dramatically the distribution of the partial pressure of MTS in the preform. Not only does it raise the minimum MTS partial pressure from 0.063 and 0.077 atm but also it shifts its location to an off-center position. It should be noted that since the fibers are oriented in the direction of the holes and the transport coefficients have their largest values in this direction (see Figure 5), the effects of the holes in Figure 13 are stronger than those that would be observed in an isotropic structure. The effects of the aspect ratio of simple preform shapes along with the effects seen in Figures 10, 11,
Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997 365
Figure 11. Effect of preform shape on conversion patterns in a parallelogram-shaped cross section. Conditions the same as in Figure 3. Base length ) 2ax ) 0.02 m; height ) 2ay ) 0.02 m. τ ) 7.5; ξs ) 0.76.
Figure 13. Contour plot of MTS partial pressure distribution in the preform shown in Figure 12 (composed of one-directional fibers) at the plane x ) 0.
acteristic size equal to the ratio of the volume of the pellet to its external surface (Aris, 1975). Obviously, this approach can also be utilized in the isobaric chemical vapor infiltration of preforms of complex shape when only a rough estimate of the average rate of deposition is needed. Figure 12. Plane and elevation views of preform with pieces cut out. ax ) ay ) az ) 0.01 m.
Summary and Conclusions
and 13 for objects of more complex geometry point to the conclusion that unless a preform has a single dominant direction of transport, that is, it has one local dimension much smaller than the others, a multidimensional (2D or 3D) model would have to be employed to obtain a satisfactory description of concentration and conversion distributions. A well-known result in the area of catalytic gas-solid reactions occurring in porous media is that a satisfactory estimate of the effectiveness factor of the catalyst, the ratio of the average reaction rate to the reaction rate at the ambient conditions, can be obtained for catalyst pellets of complex shape from a 1D model for a slab with char-
General 2D and 3D chemical vapor infiltration models were formulated to describe chemical vapor infiltration of preforms of complex shape under any reaction conditions. No restrictions were imposed on the form of the kinetics that may be employed to describe the deposition process or on the structure of the porous medium. A generalized form of the dusty-gas model was used to describe multicomponent transport in the preform in the general case of an anisotropic structure. The model was applied to SiC deposition, and anisotropic (one-directional and two-directional) and isotropic (three-directional) randomly overlapping fiber structures were used to represent the solid phase of the porous preform and
366 Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997
Acknowledgment This work was supported by a grant from the National Science Foundation. Notation
Figure 14. Contour plot of MTS partial pressure distribution in the preform shown in Figure 1 (composed of one-directional fibers) at the plane x ) 0.
to provide expressions for evolution of the structural and transport properties. The dynamic 2D and 3D CVI model equations were solved using computer codes based on the Galerkin finite element method. The domain of the problem was divided into triangular (for the 2D problem) or tetrahedral (for the 3D case) elements, and the partial pressures, temperature, and conversion profiles were represented by Lagrangian quadratic basis functions. The effects of preform anisotropy, preform aspect ratio, and irregular preform shape were examined in isobaric CVI using the 2D and 3D models. The orientationally-averaged effective diffusivity was found to be the main parameter influencing the average partial pressures and the average conversions for preforms of simple shape, but the local partial pressure and conversion values were determined by all the elements of the effective diffusivity tensor. However, even for the structure with the higher degree of anisotropy, that is, one-directional fibers, large differences among the profiles in the different directions for preforms of cubic shape appeared only after the percolation threshold is reached for transport perpendicular to the fibers. As expected, the aspect ratio had a strong influence on the model predictions, with the 3D and 2D model results away from the edges of the preform reverting to the corresponding 2D and 1D results, respectively, with increasing aspect ratio. The relative magnitude of the effective diffusivity in the direction of the increasing dimension determined how the variation of the aspect ratio affected the model predictions and the value of the aspect ratio needed to reach the behavior of the lower dimensionality model. The influence of the geometry and anisotropy of the preform on the distribution of the partial pressures and of the conversion in the preform became stronger with increasing complexity of the preform shape. These results suggest that, unless a dominant direction of transport can be identified in every section of a preform of complex geometry, a multidimensional model would have to be employed to describe satisfactorily an isobaric chemical vapor infiltration process occurring in its interior.
Symbols that do not appear here are defined in the text. Be ) tensor of effective permeabilities of the porous medium (m2) Dei ) tensor of effective diffusion coefficients of species i (m2/s) Dij, Dije ) tensors of binary diffusivity of the (i, j) pair and effective binary diffusivity, respectively (m2/s) e DKi, DKi ) tensor of Knudsen diffusivities of species i in a pore and effective Knudsen diffusivity (m2/s) N ) vector of mass transport fluxes Ni ) mass transport flux of species i (kmol/m2‚s) p ) total pressure of the mixture (atm) pi ) partial pressure of species i (atm) p ) vector of partial pressures S ) internal surface area (m2/m3) t ) time (s) T ) temperature (K) vi ) molar volume of solid i (solid S or SiC) (m3/kmol) xi ) mole fraction of species i x, y, z ) distance variables (Cartesian coordinates) (m) Greek Letters ) porosity ηj ) tortuosity factor in the j-direction λe ) effective thermal conductivity tensor (W/(m‚K)) µ ) viscosity of the gaseous mixture (kg/m‚s) νiF ) stoichiometric coefficient of species i in reaction F τ ) dimensionless time ξ ) fractional conversion of pore space to solid Subscripts b ) denotes quantities in the gas or bulk phase around the preforms i ) denotes quantities referring to species i K ) refers to Knudsen quantities S ) denotes quantities referring to the solid product 0 ) denotes quantities referring to the unreacted preform or to t ) 0 Superscripts b ) refers to bulk diffusive transport K ) refers to Knudsen diffusive transport p ) refers to viscous transport e ) denotes effective quantities * ) denotes reference quantities
Literature Cited Aris, R. The mathematical theory of diffusion and reaction in permeable catalysts; Clarendon, Oxford, U.K., 1975. Bird, R. B.; Stewart, W. E.; Lightfoot, E N. Transport phenomena; Wiley: New York, 1960. Brennfleck, K.; Fitzer, E.; Schoch, G.; Dietrich, M. CVD of SiCinterlayers and their interaction with carbon fibers and with multilayered NbN-coatings. In Proceedings of the Ninth International Conference on CVD; Robinson, M., et al., Eds.; The Electrochemical Society: Pennington, NJ, 1984; pp 649-662. Burganos, V. N.; Sotirchos, S. V. Fragmentation in random pore structures. Chem. Eng. Commun. 1989, 85, 95-112. Chung, G. Y.; McCoy, B. J.; Smith, J. M.; Cagliostro, D. E. Chemical vapor infiltration: Modeling solid matrix deposition for ceramic composites reinforced with layered woven fabrics. Chem. Eng. Sci. 1992, 47, 311-323. Currier, R. P. Overlap model for chemical vapor infiltration of fibrous yarns. J. Am. Ceram. Soc. 1990, 73, 2274-2280.
Ind. Eng. Chem. Res., Vol. 36, No. 2, 1997 367 Fedou, R.; Langlais, F.; Naslain, R. A model for the isothermal isobaric chemical vapor infiltration in a straight cylindrical pore. 1. Description of the model. J. Mater. Synth. Proc. 1993a, 1, 43-52. Fedou, R.; Langlais, F.; Naslain, R. A model for the isothermal isobaric chemical vapor infiltration in a straight cylindrical pore. 1. Application to the CVI of SiC. J. Mater. Synth. Proc. 1993b, 1, 61-73. Fieldhouse, I. B.; Hedge, J. C.; Lang, J. I.; Waterman, T. E. Thermal properties of high temperature materials. Armour Research Foundation, WADC-TR-57-487, AD-150954, Contract AF-33(616)-3701, 1958. Finlayson, B. A. Nonlinear analysis in chemical engineering; McGraw-Hill: New York, 1980. Gupte, S. M.; Tsamopoulos, J. A. Densification of porous materials by chemical vapor infiltration. J. Electrochem. Soc. 1989, 136, 555-561. Jackson, R. Transport in porous catalysts; Elsevier: New York, 1977. Kirkpatrick, S. Percolation and conduction. Rev. Mod. Phys. 1973, 45, 574. Lin, Y. S.; Burggraaf, A. J. Modelling and analysis of CVD processes in porous media for ceramic composite preparation. Chem. Eng. Sci. 1991, 46, 3067-3080. Mason, E. A.; Malinauskas, A. P. Gas transport in porous media: The dusty-gas model; Elsevier: New York, 1983. Mori, M. The finite element method and its application; Macmillan: New York, 1986. Naslain, R. In Ceramic matrix composites; Warren, R., Ed.; Chapman and Hall: Glasgow, U.K., 1992; Chapter 8. Naslain, R.; Langlais, F.; Fedou, R. The CVI processing of ceramic matrix composites. J. Phys. 1989, 50, C5-(191-207). Ofori, J. Y. Ph.D. dissertation, University of Rochester, Rochester, NY, 1996. Ofori, J. Y.; Sotirchos, S. V. Optimal pressures and temperatures for isobaric, isothermal chemical vapor infiltration. AIChE J. 1996a, 42, 2828-2840. Ofori, J. Y.; Sotirchos, S. V. Structural model effects on the predictions of chemical vapor infiltration models. J. Electrochem. Soc. 1996b, 143, 1962-1973. Ofori, J. Y.; Sotirchos, S. V. Multicomponent mass transport in chemical vapor infiltration. Ind. Eng. Chem. Res. 1996c, 35, 1275-1287. Rao, S. S. The finite element method in engineering, 2nd ed.; Pergamon Press: Oxford, U.K., 1989. Sheldon, B. W.; Chang, H. C. Minimizing infiltration times during the final stage of isothermal chemical vapor infiltration. Ceram. Trans. 1993, 42, 81-93. Sotirchos, S. V. Multicomponent diffusion and convection in capillary structures. AIChE J. 1989, 35, 1953-1961. Sotirchos, S. V. Dynamic modeling of chemical vapor infiltration. AIChE J. 1991, 37, 1365-1378. Sotirchos, S. V. Convection and diffusion flux models for anisotropic porous media. 1996, in preparation.
Starr, T. L.; Smith, A. W. 3-D modeling of forced-flow thermalgradient CVI for ceramic composite fabrication. Mater. Res. Soc. Proc. 1990, 168, 55-60. Stinton, D. P.; Caputo, A. J.; Lowden, R. A. Synthesis of fiberreinforced SiC composites by chemical vapor infiltration. Am. Ceram. Soc. Bull. 1986, 65, 347-350. Strang, G.; Fix, G. J. An analysis of the finite element method; Prentice-Hall: Englewood Cliffs, NJ, 1973. Stroud, A. H. Approximate calculation of multiple integrals; Prentice-Hall: Englewood Cliffs, NJ, 1971. Svehla, R. A. Estimated viscosities and thermal conductivities of gases at high temperatures. NASA Tech. R-132; Lewis Research Center: Cleveland, OH, 1962. Tomadakis, M. M.; Sotirchos, S. V. Effective Knudsen diffusivities in structures of randomly overlapping fibers. AIChE J. 1991a, 37, 74-86. Tomadakis, M. M.; Sotirchos, S. V. Knudsen diffusivities and properties of structures of unidirectional fibers. AIChE J. 1991b, 37, 1175-1186. Tomadakis, M. M.; Sotirchos, S. V. Transport properties of random arrays of freely overlapping cylinders with various orientation distributions. J. Chem. Phys. 1993a, 98, 616-626. Tomadakis, M. M.; Sotirchos, S. V. Effective diffusivities and conductivities of random dispersions of nonoverlapping and partially overlapping unidirectional fibers. J. Chem. Phys. 1993b, 98, 9820-9827. Tomadakis, M. M.; Sotirchos, S. V. Transport through random arrays of conductive cylinders dispersed in a conductive matrix. J. Chem. Phys. 1996, 104, 6893-6900. Transvalidou, F. Ph.D. Dissertation, University of Rochester, Rochester, NY, 1996. Transvalidou, F.; Sotirchos, S. V. Effective diffusion coefficients in square arrays of filament bundles. AIChE J. 1996, 42, 24262438. Warren, R., Ed. Ceramic matrix composites; Chapman and Hall: Glasgow, U.K., 1992. Zienkiewicz, O. C.; Taylor, R. L. The finite element method (4th Edition), Volume 1: Basic formulation and linear problems; McGraw-Hill: London, 1989a. Zienkiewicz, O. C.; Taylor, R. L. The finite element method (4th Edition), Volume 2: Solid and fluid mechanics, dynamics and nonlinearity; McGraw-Hill: London, 1989b.
Received for review August 9, 1996 Revised manuscript received October 29, 1996 Accepted November 13, 1996X IE960487D
X Abstract published in Advance ACS Abstracts, January 1, 1997.
