carbon composites by a film-boiling chemical vapor infiltration process

Chemical Engineering Science 61 (2006) 7509 – 7527 www.elsevier.com/locate/ces

Analytical and numerical study of the densification of carbon/carbon composites by a film-boiling chemical vapor infiltration process Nathalie Nadeau a,b , Gérard L. Vignoles a,∗ , Claude-Michel Brauner b a Laboratoire des Composites ThermoStructuraux (UMR 5801 CNRS-SNECMA-CEA-University Bordeaux 1), 3, Allée La Boëtie, Université Bordeaux 1,

F 33600 Pessac, France b Laboratoire de Mathématiques Appliquées (UMR 5466 CNRS-University Bordeaux 1), 351, Cours de la Libération, Université Bordeaux 1,

F 33405 Talence Cedex, France Received 21 November 2005; received in revised form 11 July 2006; accepted 8 August 2006 Available online 25 August 2006

Abstract The film-boiling chemical vapor infiltration (CVI) process is a fast process developed for composite material fabrication, and especially carbon/carbon composites. In order to help define optimal conditions, a local 1D model has been developed to study the densification front which establishes itself during the processing of a carbon/carbon fibrous preform. The model features heat conduction, precursor gas diffusion, densification reactions and structural evolution of the porous medium. The effects of total mass flux, Thiele modulus, porous medium geometry on front behavior parameters such as width, velocity and residual porosity are presented as semi-analytical correlations. An existence criterion is produced, which involves a minimal heat flux. Comparison between process-scale experiments and simulation is then possible by inserting the semi-analytical results achieved in the local study of the front into a light numerical model involving the entire preform. The model has been validated with respect to previous experimental and numerical data. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: C/C composites; Chemical vapor infiltration (CVI); Densification front; Numerical modelling

1. Introduction Ceramic-matrix composites (CMC), and particularly carbon/carbon composites (C/C), are more and more frequently used in high-temperature applications, such as rocket nozzle heat shields, airplane brakes, and furnace components. The current development of these materials is limited by economical and technical constraints. One of the most attractive methods for high-quality CMC production is chemical vapor infiltration (CVI) (Naslain and Langlais, 1990; Besmann et al., 1991) which makes it possible to obtain composite materials with excellent mechanical and thermal properties, but at the expense of high fabrication costs. This method is based on

∗ Corresponding author. Tel.: +33 5 56 84 47 34; fax: +33 5 56 84 12 25.

E-mail addresses: [email protected] (N. Nadeau), [email protected] (G.L. Vignoles), [email protected] (C.-M. Brauner). 0009-2509/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.08.027

the infiltration of a gaseous precursor inside a preform, i.e., a porous arrangement of fibers which will reinforce the composite. When the preform is heated, the gases decompose and chemical deposition of the matrix occurs. For C/C composite materials, the preform made of carbon fibers is densified by a pyrocarbon deposit originated in the cracking of gaseous pyrocarbons, usually at high temperature (ca. 1000–1300 K) and at low pressures. During the last 10 years, many studies have been devoted to a better understanding of CVI of carbon-based materials. The main goal is to obtain better control and a high efficiency for preparing C/C composites. The classical implementations of CVI are isothermal-isobaric (I-CVI), forced (F-CVI) and thermal gradient CVI (TG-CVI). I-CVI allows to process batches of numerous pieces with complicated shapes, but requires long processing times. F-CVI only works in practice when a thermal gradient is also present. Thus, the key of processing time reduction lies in the control of a thermal gradient. This is particularly true when using a promising modification of TG-CVI technique : the “Film-Boiling” CVI

7510

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

Safety valve Heating device

Goretex® Preform

Condenser

Film-boiling zone

Diphasic zone

Boiling precursor

Goretex®textile sheets Densified material

90 mm

Resistor (hollow)

Densification zone

Preform 8

Inner flask

10

25

25.5

r (mm)

Current generator 20 mm

50 mm

Fig. 1. Descriptive scheme of the film-boiling CVI process.

(also called “Kalamazoo”) (Houdayer et al., 1984). In this technique, a central susceptor or resistor brings heat to the inner (“hot”) face of the fibrous preform and the thermal gradient is obtained by exposing the outer (“cold”) face to the liquid precursor (e.g. cyclohexane for pyrocarbon matrices) in which the preform is directly immersed. The precursor is partially vaporized in situ; a densification front is created by the strong temperature gradient inside the porous preform and moves from the hot face towards the cold face (Delhaès, 2003). Fig. 1 is a schematical illustration of this setup. One of the most important issues in “Kalamazoo” control is to find the optimal values of technological parameters, providing both minimal infiltration times and required quality. A previous numerical study (Lines et al., 2005; Vignoles et al., 2006) has successfully accounted for the process behavior; a discussion of the results has given a particular importance to the notion of densification front. Indeed, this notion had already been evidenced before in various studies of CVI with thermal gradients (see e.g. Melkote and Jensen, 1990); however, there is a need for a more thorough understanding of its nature and behavior. This is the aim of the present paper. A simplified one-dimensional model for Kalamazoo and other TG-CVI processes is proposed: it is subdivided in a local front model and a global process model that makes use of the results from the local model. They are based on balances of heat and masses of precursor gas and solid. In a first part, the local model is built in a form suited to mathematical analysis. Analytical and numerical results are then presented and discussed. The more global process model is presented in a second part. It merges the laws obtained for the local front behavior into the whole preform. The local front behavior is described using only three characteristic parame-

ters: its width, its velocity, and the residual porosity left behind it. Making use of these quantities, mass and heat balances are then solved for the preform zones outside the preform, allowing to provide a very light numerical routine in a 1D geometry. Numerical results are given, and a validation with respect to the previous, more detailed numerical approach and to experimental data is made. 2. Local front modelling This part is devoted to the presentation and discussion of a one-dimensional model restricted to the front region, where various simplifications are made with respect to the previous model (Vignoles et al., 2006), in order to make an analytical study possible. 2.1. Model set-up The densification process involves the coupling of heat and mass transfer, as well as structural changes of the porous medium, which occur simultaneously during the process. A satisfactory starting point for the model set-up is to consider the deposition reaction rate. After homogenization on a small porous medium volume, it reads R =
v ().k(T ).C  ,

(1)

where  is the reaction order. The heterogeneous rate k(T ) follows an Arrhenius law
 Ea , (2) k(T ) = k0 exp − RT

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

7511

Table 1 Balance equations and constitutive laws used in the local model Subset

Balance equation

Constitutive laws  conv J = Cv g = −CK −1 ∇p J diff = −D(, T )∇C  cond q = −(, T )∇T q conv = c(, T )vg T

jC + ∇ · (J conv + J diff ) = −R jt jT + ∇ · (q cond + q conv ) = 0 c(, T ) jt j − = Vm R jt 

Gas mass balance Heat balance Solid mass balance

R =
v () · k(T ) · C 

Table 2 List of provided parameter values Symbol

Meaning



Thermal conductivity

Expression (1 − )0 (1 − )c0

Specific heat capacity

4 × 106 J m−3 K−1

c0

Solid phase capacity



8.2 W m−1 K−1

0

Solid phase thermal conductivity

c

2 × 103 kg m−3

Deposit density

Vm

Deposit molar volume

D

Effective gas diffusivity

Ms


v

6.0 × 10−6 m3 mol−1

s

D0 m+1

Pure gas diffusion coefficient

9 × 10−5 m2 s−1

D0 A1/n (1 − B )

Internal surface area Parameters

4.0 × 105 m−1

A B




Kinetic constant

Ea k0 exp − RT k0

Activation energy

Ea



Deposition rate order

1

q

Heat flux

k

Chemical deposition rate constant

Ea . RTref

(3)

Thus, (2) can be written



Tref −1 , k(T ) = k(Tref ) f (T ) = k(Tref ) exp − T



1/2

1.1 × 108 m s−1 3.26 × 105 J mol−1 5 × 105 W m−2

where Ea is an activation energy. By choosing a reference temperature, a dimensionless activation energy may be defined =

Value/Unit

(4)

where k(Tref ) = k0 exp(−) is a suited reference value for the reaction rate constant. Eq. (1) implies that densification is controlled essentially by three factors: (1) The temperature T, through its influence on the heterogeneous reaction rate k. (2) The precursor gas concentration C. (3) The pore volume fraction, or porosity , through its influence on the internal surface area
v . Any mathematical model of densification has to feature at least balance equations for three quantities, which are chosen to be T, C, and . They are, respectively, an energy balance, a precursor

gas mass balance, and a solid mass balance. The equations are described in Table 1. Associated properties are listed in Table 2. The equation used for mass transfer is fairly simple, since the transport and reaction of only one species (the precursor) has been considered. Indeed, the counter-flow of gaseous byproducts is only weakly coupled to the precursor consumption, i.e., it has no direct influence on it. It is also assumed that the gas velocity follows Darcy’s law, and experimental evidence shows that the pressure drop is low (under 10 Pa). This has also been confirmed by the preceding numerical study (Vignoles et al., 2006). Indeed, the mass Péclet number is
 K
p 1 J conv  = P e´m  . . (5) diff J   D The permeability is of order 10−12 –10−13 m2 , and the viscosity is around 10−5 Pa s, so the estimates of P e´m lie around 10−3 –10−4 . Hence, convective flow may be discarded. Similarly, in the heat equation, the convective flux is found negligible since the thermal Péclet number is small: q conv  K()
p c(, T ) = P e´th   1. cond q   

(6)

7512

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

Moreover, the heat source term arising from chemical deposition is negligible (Vignoles et al., 2006). The most important structural parameters are then the effective diffusivity (D), the internal surface area (
v ), the volume heat c and the thermal conductivity . Since the variations of D with T are small with respect to the variations of k with T, it is chosen to neglect them. Also, it has been proved (Vignoles et al., 2006) that neglecting Knudsen diffusion (rarefied gas transport) does not introduce strong errors in the model. Exponent m = j ln D/j ln  − 1 is a traduction of the notion of tortuosity, and m + 1 is known as Archie’s law exponent. If m were 0, this would mean that a simple law of mixtures is valid for diffusivity, as e.g. when the pores are perfectly parallel at all densification stages. Usually, in whatever actual porous medium, m is strictly positive, denoting the existence of some tortuosity effect (Tomadakis and Robertson, 2005). Values of up to 10 have been reported, but usual values range between 1 and 2 (Tomadakis and Sotirchos, 1993). Exponent n = (j ln
v /j ln )−1 relates the evolution of the specific surface with porosity. For example, it is shown by straightforward area and volume computations that the value of n is 2 in a medium formed by parallel cylindrical pores, and 23 in a medium formed by spherical pores. In a random pixelized binary image, exponent n is 1 (Burganos, 1998). In some other random porous media, logarithmic dependencies may also occur, e.g.
v ∝ − ln  when porosity is close to zero for random parallel fibers (Tomadakis and Sotirchos, 1991; Vignoles et al., 2006). Rikvold and Stell have developed a model for partially overlapping families of spheres and cylinders (Rikvold and Stell, 1985a,b; Tomadakis and Sotirchos, 1991; Coindreau and √ Vignoles, 2005); in that case it can be shown that
v ∝  − ln  when porosity is close to zero. In this work, it will be assumed that values of n ranging between 1 and 2 are representative of actual porous media. The volume heat c and thermal conductivity  are assumed to follow a linear law with respect to porosity (i.e., a law of mixtures), where the solid phase only has been accounted for, since the gas-phase related properties are negligible with respect to the solid-phase ones. Media with bimodal pore size distributions (e.g. 3D weave or cloth lay-up) or continuous pore size distributions (i.e., fractal media) may not be correctly described with the structural parameter laws described here, principally because they are only functions of porosity; in that case, more elaborate modelling has to be performed, since the introduction of pore size distribution affects the whole approach. The system has been further simplified by making use of the following assumptions: (i) The well-defined moving front observed in experiments leads us to consider a 1D moving coordinate frame attached to the front velocity, v. As in most combustion problems, stationary solutions are sought in that frame. (ii) Introducing the reference quantities listed in Table 3, a study of characteristic times, summarized in Table 4, allows to neglect the transient times that appear in heat and gas mass balances, which are very short with respect to

Table 3 List of reference quantities Symbol

Meaning

Value

Tref Lref Cref

Reference temperature Reference length Reference concentration

1373 K 7.8 × 10−3 m 8.9 mol m−3

Table 4 List of characteristic time values Symbol

Meaning

 2

Densification Chemical reaction

3

Gas diffusion

4

Heat diffusion

Definition

Value −1

(Vm Cref Ak(Tref )) (Ak(Tref ))−1 L2ref D0 c0 L2ref

0

9.1 × 102 s 0.06 s 0.41 s 15.0 s

the porous medium evolution time. Accordingly, the introduced velocity (v) appears only in the densification equation (8). (iii) An ignition temperature Ti below which the reaction rate is considered negligible, is introduced. This allows to modify the Arrhenius law into an expression of k(T ) where k equals strictly 0 for all temperatures lower than Ti . There are different possible choices for f (T ) in Eq. (4). In this study, the following choice is retained: 

⎧ Tref ⎪ ⎪ ⎪ exp − T − 1 ⎨

 Tref (7) f (T ) = ⎪ − exp − −1 if T > Ti , ⎪ ⎪ Ti ⎩ 0 otherwise. Since there is no heat source or sink, the heat equation is immediately integrated once, and the heat flux q that traverses the front appears as a prime integral. The system is now ⎧ dT ⎪ = q, −0 (1 − ) ⎪ ⎪ dx ⎪  ⎨ d
dC −m+1 D0 = −A(1 − B)k(Tref )f (T )C, ⎪ dx dx ⎪ ⎪ ⎪ ⎩ d v = Vm A(1 − B)k(Tref )f (T )C. dx

(8)

This problem statement has to be completed by appropriate conditions on the boundaries of the local study domain. 2.1.1. Boundaries The front’s boundaries are the locus of points for which the deposition rate R tends to zero. From Eqs. (1) and (7), and Table 2, this is possible in three cases • when T < Ti , a “cold side” is defined, • when C → 0, a “depletion side” is defined, • when  → 0, a “complete densification side” is defined.

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

ψ

1 εc

ε

γ

7513

ψ ε

γ

θ

1 εc

θ

ψh

εh

θi

0 − ξf

0

0 ξ

θi −ξf

0

ξ

Fig. 2. Depletion front: profiles of , , and  in [−f , 0].

Fig. 3. Total densification front: profiles of , , and  in [−f , 0].

Quite obviously, if there is always a “cold side”, the opposite side (say the “hot side”) may be either a “depletion side” or a “complete densification side”. In a 1D geometry, the local analysis of the front will thus be restricted to a bounded interval I = [−xf , 0]. By convention the origin (x = 0) is the cold side, where the deposition rate is zero and the porous medium is not yet densified. It is described by the following Cauchy data:  C(x = 0) = Cc , (9) T (x = 0) = Ti , (x = 0) = c .

width depends on temperature through the reaction rate. In this way, the dimensionless front width is expected to be of order unity. Using this reference length, the traditional expression for a Thiele modulus, representing the competition between reaction and diffusion, is (Thiele, 1939)
v k = Lref . (13) D

The cold side receives a positive mass flux coming from regions where the porous medium is not densified, defined as D0 Jc = J (x = 0) = m+1 c

dC (x = 0). dx

(10)

The hot side lies at x =−xf . At this point, one has to distinguish two cases (Nadeau et al., 2004), as sketched at Figs. 2 and 3 (a) If the densification is partial  T (−xf ) = Th , (−xf ) = h = 0, C(−xf ) = 0, where the residual porosity is noted h . (b) The densification is total  T (−xf ) = Th , (−xf ) = 0, C(−xf ) = Ch = 0.

(11)

(12)

In this case, there is a residual reactant concentration noted Ch . 2.1.2. Dimensionless problem The last step of the model construction consists in rewriting the problem in dimensionless form. Reference quantities, dimensionless numbers and variables are listed in Table 5. It has been chosen to use a thermal-transfer related reference length including the scaled activation energy since the front

Making use of the expressions listed in Table 2, this quantity is in our case 2 = 20 g()f (T ),

(14)

where 20 = L2ref

Ak(Th ) D0 (Th )

and g() is an “infiltrability function”. The dimensionless system is ⎧ d

1 ⎪ =− , ⎪ ⎪ ⎪ d 1 −   ⎪ ⎨ d
d −m+1 = − 20 1/n (1 − B)f ( ) , ⎪ d d ⎪ ⎪ ⎪ ⎪ ⎩ d = 1 1/n (1 − B)f ( ) d

(15)

(16)

with the following expression for the Arrhenius-like term


/ − 1 f ( ) = exp Ze − exp(−), (17) (1 + ( / − 1)) where  = (Th − Ti )/Th is a thermal expansion term and Ze =  is a Zeldovich number. This quantity is usually large, due to the high activation energy value. However, it will not be intended to perform a traditional asymptotic analysis in the limit of Ze → ∞. The sum of the last two equations of system (16) leads to
 d d −m+1 + 20  = 0. (18) d d

7514

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

Table 5 List of dimensionless quantities, numbers and variables Symbol

Meaning

Reference quantities Th

Temperature

qref

Heat flux

Definition

Value (unit)

1373 K 5 × 104 W m2

Lref

Length

R0 Th2

Cref

Concentration

Cc

7.8 × 10−3 m

Ea q

8.9 mol m−3 −1



Time

vref

Velocity

Jref

Mole flux

 Cc vref

c

Raw preform porosity

–

Ti

Ignition temperature

(Vm Cc Ak(Th )) Lref

9 × 102 s 0.87 × 10−5 m s−1 7.7 × 10−5 mol m−2 s−1 0.94 1000 K

Dimensionless numbers (estimated values)



Condensation ratio



Activation energy



Thermal expansion ratio

Ze

Zeldovich number

Vm Cc Ea RTh Th − Ti Th

6.3 × 10−5



7.8

x Lref C Cc T − Ti  Th − Ti

[−f , 0]

28.6 0.27

Dimensionless variables and parameters



Space variable



Concentration



Temperature

Mole flux

0

Thiele modulus

d − 20 . d

[0; +∞)

Jc j = Jref  20 j Ak(Th ) Lref D0 −

It is of practical interest to introduce a constant K which is the prime integral of (18). K is the total mass flux in the moving reference frame (which makes solid deposition appear as a convective-like term): K = m+1

[0, ]

v vref

Velocity

J˜c

[0, 1]

(19)

Note that K and J˜c are related to each other through the relation 

1 K J˜c = (20)

c + 2 .  0 The system is now reduced to a set of three ordinary differential equations (ODEs) ⎧ d

1 ⎪ ⎪ =− , ⎪ ⎪ d 1− ⎪ ⎨ d (21) = K−(m+1) + 20 −m , ⎪ d ⎪ ⎪ ⎪ ⎪ ⎩ d = 1 1/n (1 − B)f ( ) d

1+m

with Cauchy data at  = 0 (cold side) 

(0) = 0, (0) = 1, (0) = c

[0; +∞) [0; +∞)

(22)

and two alternative sets of boundary conditions on the hot side, for the two kinds of possible fronts at  = −f : ⎧

(−f ) = , ⎪ ⎨ (−f ) = 0, (−f ) = h (reactant depletion front), ⎪ or ⎩ (−f ) = h , (−f ) = 0 (total densification front). (23) 2.2. Analytical results 2.2.1. Upper bound for K Since the front progression exclusively relies on irreversible deposition, the total mass flux traversing the front in the moving frame can not be negative. However, it is not necessarily zero,

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

and it is of interest to find an estimate of its highest possible value. The slope of concentration is

2 d K = m+1 + m 0 . d  

(24)

At  = −f , it is convenient to define d . Υh = Υ (−f ) = (−f ). d

(25)

In the case of a depletion front, since (−f ) = 0 and that cannot be negative, then Υh has to be positive. In the case of a full densification front, this quantity tends to infinite as −m as  tends to −f . It will be shown later that K has to be exactly zero in this case, so the rest of this discussion will focus on the depletion case. At  = 0, one has  20

d J˜ = (0), m+1 c d c

(26)

where J˜c is the scaled impinging gas flux, counted positive when the gas is brought to the front from the cold side, which is always the case. From Eq. (20), is seen that the value of K at  = 0 is K = 20 (J˜c − c ).

(27)

As is also a positive value, an upper bound for K appears K
20 J˜c .

(28)

In other words, the total mass flux in the moving frame cannot exceed the reactant mass flux. Note that it has not to be equal to it, because of the possibility of a residual porosity or concentration: this renders the study of this front much more elaborate than for e.g. a smoldering front (Dosanjh et al., 1987). 2.2.2. Criterion for full densification The second point is the determination of the conditions for which a total densification front is possible. This depends precisely on the competition between diffusion and reaction at the hot side of the front: if reaction is too strong, then only a depletion front will be possible. The “infiltrability function” has the following expression: g() = 1/n−1−m (1 − B).

(29)

Since usual values of n are superior to 1 and usual values of m are positive, then it is expected that the Thiele modulus will increase when porosity goes to zero, rendering more effective the diffusional limitations in this limit. So one expects that the best infiltrable material will be the one with the lowest possible values of m and n. A local study of the hot side of the front allows to produce, among other results, a more precise criterion. Choosing  as a new independent variable, the autonomous 3-variable differential system (21) is reduced to a non-autonomous 2-variable differential system for ( (), (), ). The porosity is a monotonous, increasing function of , since there is only

Values for “usual” porous media NO FRONT m (diffusion exponent)

Υ =

7515

5.

(0) Degenerate front (1) Depletion front

A

2. (0) Degenerate front (1) Depletion front (2) Total densification front

1. 0. 1.

2.

B

3. n (internal surface exponent)

Fig. 4. The different possible fronts in (m, n) parameter subspace.

growth and no erosion of the solid phase, so few information is lost when performing this step: ⎧ j

⎪ =− , ⎪ 1/n (1 − B) f ( ) ⎪ j (1 − ) ⎪ ⎪ ⎨ j K = m+1/n+1 (30) j  (1 − B) f ( ) ⎪ ⎪ ⎪ 2 2 ⎪

0 ⎪ ⎩ + m+1/n .  (1 − B) f ( ) When the porosity is close to h , the temperature is close to . So, the non-autonomous equation for becomes in that region d 2 /2

∼ ( 20  + K)−(m+1/n+1) . d f ()

(31)

From this relation, it follows clearly that: (a) If m + 1/n > 1, then  may not go to zero, i.e., h has to be strictly positive. The only way to have a well-defined hot side is to ensure h = 0. This is indeed possible and is a depletion front, as sketched in Fig. 2. In that case, K may only be positive. (b) If m + 1/n < 1, then (i) if K = 0, then h
0 and h
0, excluding the existence of the front, (ii) or K = 0 and three sub-cases are possible: • h
0, h = 0 (depletion front), • h = 0, h = 0 (limiting case), • h = 0, h
0 (full densification front, sketched at Fig. 3). (c) If m + 1/n + 1 = 1, the study is more algebraically involved because of the logarithm that shows up in the integration; however, it is of poor practical significance. This leads to an existence domain for two kinds of front in the (m, n) parameter space (see Fig. 4), as well as conditions on K: m + 1/n − 1 > 0: only a depletion front is possible (see Fig. 2). Here, K belongs to [0,  20 J˜c ]. m + 1/n − 1 < 0: a total densification is also possible (see Fig. 3). Here, K is exactly zero.

7516

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

However, when n = 2, the calculus is simpler √ √ 2

lim = 2(F () − F (0))(arctanh( c ) − arctanh( h )).

Unlike expected from the study of the g function, it is seen that larger values of n will favor full densification. However, this criterion confirms the key role of the evolution of diffusivity and internal surface with porosity. The usual values for m and n in porous media (1
m
10, 1
n
2) exclude a total densification front. The obtained criterion is extremely severe; however, an ideal medium like a bundle of parallel cylindrical fibers would satisfy it. On the other hand, a depletion front with extremely low values of h is possible in many cases; in practice, it is very difficult to distinguish from the rigorously total densification case. Since in the case of a total densification front, the prime integral K has to be zero, the expression of the velocity is straightforward to obtain

These expressions show that the velocity has a finite limit when h tends to zero, and that this limit is of order unity in dimensionless units, due to the proper choice for the definition of

and f. The second region is characterized by ≈  and  ≈ h , noting that h may be arbitrarily close to zero. Let us introduce Y defined by

˜

= −1 c Jc .

j j j = · jY jY j

(32)

This makes this kind of front relatively similar to a smoldering front. On the other hand, if K is not equal to zero, then the front is a depletion front.

 = h Y,

(36)

Y ∈ [1, ∞].

(37)

Using the non-autonomous equation (31) with (37), one has

= − h

K

m+1 Y m+1 h

20

+

1/n h Y m



m+1/n 1/n f ()h Y

. (38)

2.2.3. Asymptotic behavior of the depletion front for low residual porosities The depletion front may be analyzed through an asymptotic approach. The principle of the asymptotic development lies in the decomposition of the concentration profile into two regions:

Thus, integrating (38) on [1, Yc ] and taking the limit Yc → ∞ leads to 20 1 1 1

K

2 + . = m+1/n m+1/n−1 2 m + 1/n f ()  m + 1/n − 1 f ()  h

h

(39)

• A region, close to the cold side, where decreases (when  decreases) very slowly from the right-hand Cauchy value of 1, • A region close to the hot side, where decreases (when  decreases) sharply to zero, while is close to .

This is a relation between K, 0 and h . According to it, when either K or 0 are small, then the residual porosity may be n also small. In this relation, the expression of lim can be used if h tends to zero. Two cases may be distinguished, depending on the relative importance of the two terms in the right-hand side of Eq. (39):

In the first region, since ≈ 1, the non-autonomous equation for temperature becomes

• if K  20 h , then

d

≈− . d (1 − )f ( )1/n (1 − B)

n

(33)

F () − F (0) = c . 1/n h d/(1 − )(1 − B)

(34)

Let us remark that this velocity value depends on n and the calculus of the integral is not explicit in the general case. For sake of simplicity, parameter B will be set to zero in the rest of this section; this does not change the qualitative nature of the results. An expression available for all n > 0 is ⎛

 ∞ k   −1/n c n n+

= (F () − F (0)) ⎝c k − 1/n k=0 ⎛ ⎞⎞ j ∞  h −1/n ⎝ ⎠⎠ . − h n+ (35) j − 1/n j =0

,

• in the converse case, then

The primitive function of f ( ) will be noted F ( ) and Eq. (33) is integrated from ( = h ) =  to ( = c ) = 0, giving the following relation on n

2/(m+1/n−1)

h  (2( lim )2 /f ()(m + 1/n − 1))1/(m+1/n−1) 0 n

h  (2 lim /f ()(m + 1/n))1/(m+1/n) K1/(m+1/n) . The criterion may be rewritten using the obtained expressions for h under the form 2/(1−1/(m+1/n))

• if K  0

, then

n

2/(m+1/n−1)

h  (2( lim )2 /f ()(m + 1/n − 1))1/(m+1/n−1) 0

,

• in the converse case, then n

h  (2 lim /f ()(m + 1/n))1/(m+1/n) K1/(m+1/n) . Note that in any case, h → 0 implies that K → 0. So, the mass flux J˜c tends to a lower bound lim J˜c = lim c .

h →0

We will see later a physical interpretation of this result.

(40)

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

7517

, Φ0 fixed

εh ~ εc

1 εc

ε

Choose initial guess for ω

γ ψ θ

Input

θi

0 − ξf

I. C.

0 ξ

Fig. 5. Degenerate front: profiles of , , and  in [−f , 0].

Going back to the autonomous system (21), and noticing that  is close to c in a large subset of [−f ; 0], one has the following useful lower bound for f f (1 − c )( − i ).

(ξ = ξ + ∆ξ) ODE Solver for ε, ψ, θ

Increase ω

The last four items may be deduced from the first one by a straightforward asymptotic analysis. Such a case is sketched in Fig. 5. By setting   c in system (16), integration of the

and equations yields

no

ε<0?

ε (-ξf) ψ (-ξf)

Return ξf and

θ (-ξf)

<

θ (-ξf) > γ ?

>

= Return ω, ξf, εh or ψh

κJc = ω εc +

(42)

This provides an existence condition on m since h has to be positive ln(f J˜c −2 0 ) m< . ln c

Decrease ω

(41)

• A residual porosity h which is close to c , that is, densification is negligible. • Linear temperature and concentration profiles. • A dimensionless front width close to , i.e., a large but finite quantity. • A divergence of front velocity as (c − h )−2 . • A divergence of the associated mass flux as 1 .

f = (1 − c )( − i ), −(m+1) . h = 1 − f (K + 20 c )c

no

ψ<0?

2.2.4. Degenerate front Along with the two kinds of front already mentioned, a “degenerate” solution is also possible. It is characterized by



ε (0) = εc ψ (0) = 1 θ (0) = θi

2

Φ0

Fig. 6. The shooting algorithm.

2.3. Numerical method (43)

This means that for m above this limit—depending on J˜c —no front at all is possible. Such a front has not been observed in practice, since it would require some violent phenomenon (detonation or else) to show up. Such conditions would be also in contradiction with the hypotheses of incompressible flow and of slow gas convection that underlie the whole model set-up. Accordingly, no special attention will be payed to this solution.

System (21) constitutes a nonlinear eigenvalue problem, since there are more independent Cauchy data (treated as boundary conditions) than independent differential relationships, but some parameters are yet unknown, like the front velocity, residual porosity, and width. The system is studied numerically with a combination of shooting and continuation methods. The shooting method is used to determine solutions for a given fixed set of input parameters. Then, continuation is applied in order to follow the solutions with input parameter variations.

7518

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

Table 6 Summary of the different studied cases

m

n

Variable Variable Fixed Fixed Fixed

Variable Variable Fixed Fixed Variable

Fixed Fixed Variable Variable Fixed

Fixed Fixed Fixed Fixed Fixed

1e+06

No front

10000

Front 100

1

0.01 0.001

0.01

0.1

1

Dimensionless inverse heat flux or reference Thiele modulus

Fig. 7. Domain of front existence in ( 0 , log(J˜c )) parameter space.

io

n

flu

x

100

al

di

ffu s

10 1

Front

0.1 w Lo y

ea

sit

rr fo

o or

lp

ua sid re

No front

0.01

x flu

0.1

1

10

on

n

0.001

gi re

io ct

2.4.1. Case 1 This case is the most directly significant for the practical implementation of the process, since it yields results on the front behavior as a function of control parameters. Indeed, parameter 0 relies on the reference length Lref = 0 RTh2 /Ea · 1/q and may thus be thought of as a reciprocal heat flux for given porous medium and gas properties and given reference temperature. The second control parameter is indeed a mass flux parameter, and we have chosen to use J˜c = Jc /(Cc k(Th ) · ALref ). Results have been produced for m = 2, n = 2 and m = 41 , n = 2; they are very similar, so only results for m = 2, n = 2 are shown. The domain of front existence in the ( 0 , J˜c ) parameter space is shown in Fig. 7. In logarithmic coordinates, is has a roughly

1e+08

ry sa

2.4. Results and interpretation

or or h = 0/

es ec

These cases are summarized in Table 6.

h = 0/ h = 0/

h , f h , f h , f h , f h , f

N

• The first one is a search for the dimensionless variables , h and f for the two kinds of front, (i.e., for two appropriate sets of parameters m and n), K and 0 being the two continuation parameters. • The second part is a description of the front behavior when the diffusion exponent m varies, while all other parameters are fixed. • The last part is a study of the behavior of h , , and f when the densification is total. In this case, one has to state that h = 0, K = 0. Suitable fixed values have been chosen for m and n; so, only 0 is the continuation parameter.

,

,

,

,

,

im

The shooting algorithm is described in Fig. 6. For sake of commodity, parameter K has been used instead of the more physically tractable parameter J˜c , because it allows to avoid the presence of the velocity in the concentration equation, thus making the numerical integration easier. However, J˜c is computed as a post-processing result and used for the graphical representations shown below. The “AUTO 2000” software (Doedel, 1981) has been used to compute parameter-dependent families of solutions to nonlinear equations using numerical continuation (Doedel et al., 2002). This method allows to follow the solutions even through folds in parameter-solution space. In order to study the problem under different angles and to scan as many cases as possible, the numerical approach has been divided in three parts.

Conditions and regions

ax

3

0

Front behavior variables

M

2

K

Dimensionless mass flux

1

Continuation parameters

Scaled mass flux

Case

100

Scaled heat flux

Fig. 8. Domain of front existence in scaled fluxes parameter space.

triangular shape. There is a maximal admissible value for 0 close to unity, which means that there exists a minimum heat flux below which no front exists:  0 RTh2 Ak(Th ) . (44) q D0 Ea

Residual porosity εh

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

7519

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

1e+09 1e+07 100000 1000

0.001

1 0.1 0.01

0.01 0.1

Φ0

1

10

Fig. 9. Residual porosity h vs. 0 and J˜c .

1 1.0 0.9 0.8

0.7 0.6

0.5

0.4 0.3 0.2

Residual porosity εh

0.1

0.1

0.01 Φ0 = 0.01

0.001

0.0001 0.0001

0.01

1

100

10000

1e+06

1e+08

1e+10

~

Dimensionless mass flux Jc Fig. 10. Residual porosity h vs. J˜c at different fixed values of 0 .

This means that the critical flux decreases • when the activation energy increases, • when the gas diffusion coefficient increases (with a pressure decrease unless a Knudsen regime is attained), • when the reference temperature decreases, • when the solid phase conductivity decreases, • when the initial pore diameter dp,0 increases (since A ≈ 1/dp,0 ) (Vignoles et al., 2006). Two bounds for the mass flux appear: Cc k(Th ) · ALref Jc 

Cc D0 . Lref

(45)

The upper bound is of poor practical utility since it corresponds to the degenerate front. It is related to the gas diffusion. On the other hand, the lower bound is only related to the chemical reaction rate and to the porous medium structure (remember that A is an inverse pore diameter size). When both bounds are equal, one comes back to the existence condition on q given in Eq. (44). This suggests an interesting formulation for the existence criterion: the mass flux has to be at least equal to what the chemical reaction needs and is at most equal to what diffusion is able to bring within the front, as illustrated by Fig. 8 where the coordinates are now scaled heat and mass fluxes. The evolutions of residual porosity, dimensionless front velocity and width, as computed with this model, are shown in

7520

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

Dimensionless velocity

100

10

1

1e+09

0.1

1e+07 100000 1000 0.001 0.01 0.1 1

Φ0

1 0.1 0.01

10

Fig. 11. Dimensionless front velocity vs. 0 and J˜c .

100

Dimensionless velocity ω

0.1

0.2

Φ0 = 0.01

10 0.3 0.4 0.5 0.6 0.7

1

0.8 0.9 1.0

0.1 0.0001

0.01

1

100 10000 Dimensionless mass flux

1e+06

1e+08

1e+10

Fig. 12. Dimensionless front velocity vs. J˜c at different fixed values of 0 .

Figs. 9–14. They confirm the two distinct front behaviors which were to be expected:

tends to zero when the mass flux approaches its lower bound, as predicted in the preceding section (see Fig. 10).

• the degenerate front with large residual porosity (h  c ) and high velocity ( → +∞). • the depletion front with a dimensionless velocity of order unity and a residual porosity that may be small. Indeed, it

The latter case is indeed the one that is realized in practice. The asymptotic study of this case, performed in the preceding section, is confirmed by these results. Indeed, Eq. (36) yields 2

lim ≈ 0.3, in close agreement with the numerical value. From

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

7521

2

Dimensionless front width

1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 0.001 0.01 0.01

0.1

1

10

100

0.1 1000 10000 100000 1e+06 1e+07

Φ0

1 1e+08

1e+09

Fig. 13. Dimensionless front length f vs. 0 and J˜c .

2

Dimensionless front width ξf

1.8

1.6

1.4

1.0

1.2

1 0.0001

0.01

0.9 0.8 0.7 0.6 0.5 0.4

0.3 0.2

0.1

1

100

10000

Φ = 0.01

1e+06

1e+08

1e+10

Dimensionless mass flux Jc

Fig. 14. Dimensionless front width f vs. J˜c at different fixed values of 0 .

the definition of vref in Table 5 it can be seen that the front moves slower when the reference temperature decreases and −1 when the initial pore diameter (recalling that A ∝ dp0 ), the activation energy, or the heat flux increase. Fig. 14 shows that the dimensionless width of the front has a small variation, lying between 1.2 and 2 for the whole parameter domain of existence. This confirms that the proposed choice for a reference length is appropriate with respect to the front width. From the definition of the reference length (see Table 5), it appears now that the front broadens when the hotside temperature and the solid phase conductivity increase, and when the activation energy and the heat flux decrease. Again,

the numerical values of f are satisfactorily close to the analytical lower bound f 1.15 from Eq. (41). 2.4.2. Case 2 Variations of the front properties have been investigated as a function of m and n, which are characteristic exponents describing the geometrical complexity of the porous medium. Fig. 15 is a plot of the front behavior as a function of m for n = 2, 0 = 0.3. Here, the three kinds of front are obtained: the total densification front exists for K = 0 and m lying below a critical value (here m
21 ). Above this value, a depletion front arises, and the residual porosity begins to differ from zero.

7522

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527 0.94 0.90 Degenerate front

0.80 0.70

K = 0. 0.60

Depletion front

A

B K = 0.

0.50

Residual porosity

0.40 0.30 0.20 0.10 0.00 0.

1.

2.

3.

4. 5. 6. m (diffusion exponent)

7.

8.

9.

10.

Total densification front

Fig. 15. Residual porosity h vs. exponent m at fixed n.

Residual concentration ψh

2.4.3. Case 3 The last numerical run is an investigation of the front behavior in conditions of total densification. Exponents m and n have been chosen in order to be in that case (see criterion arising from Eq. (31)). The most important difference with case 1 is that K has to be null, thus reducing by one the number of independent parameters acting on the system. For example, since J˜c = c there is no need to vary J˜c as a parameter, since its value is fixed as soon as the velocity is known. Another difference is that now h is tracked as a system response instead of h . Again, there appears a critical value of 0 , close to unity, above which no front is possible. As 0 decreases below this critical value, there is a continuous increase in h , but on the other hand dimensionless velocity and width do not vary strongly. Both are again close to unity, leading to the same conclusions as in case 1 (Figs. 16–18).

1

0.6

0.4

0.2

0

(a)

0.1

0.2

0.3 0.4 0.5 0.6 0.7 Reference Thiele modulus Φ0

0.8

0.3 Dimensionless mass flux Jc

0.35

1

2.5. Conclusions

0.8

0.6

0.4

0.2

The simplified formulation that has been proposed for the local model has shown great power in: • providing an existence criterion for the front, • providing suitable estimates of the magnitudes of the front characteristics, that is, width, velocity, and residual porosity,

0.8

0

Residual concentration ψh

A degenerate branch is also possible, as before, and it merges into the depletion front branch for a critical value of exponent m (here approximately 9.5), above which no front is possible at all. This matches the prediction from the brief analysis of the degenerate front in the preceding section.

0 0.25 (b)

Fig. 16. Total densification front: residual concentration plots with respect to (a) 0 , and (b) J˜c .

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

7523

1.94

0.44

Dimensionless front width ξf

Dimensionless front velocity ω

1.93

0.42 0.4 0.38 0.36 0.34 0.32

1.92 1.91 1.9 1.89 1.88 1.87 1.86 1.85 1.84

0.3

1.83

0 (a)

0.1

0.2 0.3 0.4 0.5 0.6 Reference Thiele modulus Φ0

0.7

0.8

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Reference Thiele modulus Φ0

(a) 1.94

0.44 Dimensionless front width ξf

Dimensionless front velocity ω

1.93

0.42 0.4 0.38 0.36 0.34 0.32

1.91 1.9 1.89 1.88 1.87 1.86 1.85 1.84 1.83

0.3 0.25 (b)

1.92

0.3 Dimensionless mass flux

0.25

0.35 ~ Jc

Fig. 17. Total densification front: velocity plots with respect to (a) 0 , and (b) J˜c .

• performing rapidly (i.e., at low CPU time costs) various parametric studies. Two kinds of physically possible fronts have been found: a partial densification front, and a total densification front. The total front is related to diffusion and surface area laws which are not likely to be realized in classical porous media. However, it will be difficult to distinguish between a low-residual porosity partial densification front (which occurs in practice when the mass flux is just above the minimal required value for existence) and a total densification front. The discussion on front existence allows to provide two guidelines for an optimal process control (Vignoles et al., 2005, 2006): (i) it is important to bring enough heat flux at the very beginning of the process, and (ii) it is unnecessary to increase the heat flux after the front has started to settle. Also, it is seen that exponents m and n play an important role, because they control the ease with which the fibrous medium is infiltrated, as translated by an “infiltrability function”. This suggests future work in the direction of optimizing fibrous architectures for infiltration, using an infiltrability functional, i.e., a generalization of g() in Eq. (14).

(b)

0.3

0.35

Dimensionless mass flux

Fig. 18. Total densification front: front width plots with respect to (a) 0 , and (b) J˜c .

3. Process modelling using the front behavior results The local front results may be used in a more global system which represents the actual process, so a “light” numerical scheme has been designed to model the whole preform containing the front. This makes it possible to perform a direct comparison between the numerical predictions and experimental data, as well as results from the more detailed numerical simulation built before (Vignoles et al., 2006), that has already been validated with respect to the same experimental data. 3.1. Construction of the model Fig. 19 shows the 1D resolution domain spanning the whole preform width LT , which is suited to the description of the experiment. If criterion (44) is satisfied, a front exists and propagates towards the cold side. It is assumed that its build-up time is relatively short with respect to the total processing time. The front velocity is considered slow with respect to the transient times associated to heat and gas mass transport. Then, at any moment, the preform is considered to be divided into three zones: a hot zone, with parameters associated to the densified

7524

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527 Table 7 List of provided parameter values for the global model

Ts

Graphite heater

Th

Ce

Ti Cc

Densified preform

Densification front

Raw preform

Te C=0

h

Tb

Symbol

Meaning

Value/Unit

LT h Tb Ts Ptot

Total preform width Heat transfer coefficient Precursor boiling point Hot-side temperature Total pressure

15 mm 200 W m−2 K−1 353 K 1373 K 1.013 × 105 Pa

m n

Exponent for diffusion Exponent for surface area Dimensionless front velocity Dimensionless front width Dimensionless mass flux

1

f J˜c

3 2

0.3 1.2 0.6

xf

Lh LT

Lc

Fig. 19. Scheme of the global 1D model featuring the front.

material, a front zone, the characteristic properties of which are given by the results of the local front study, and a cold zone, with parameters associated to the raw preform. The preform outer surface is held at a given temperature Te and concentration Ce , by maintaining it in contact with the boiling precursor. The heat transfer coefficient with the exterior biphasic medium has been identified in the previous model. On the other side, the graphite heating device has been regulated so as to maintain a constant temperature Ts . At any time t = 0, the front will be lying between a raw preform zone (porosity c ) and a densified preform zone (porosity h , considered as close to zero). The balance equations for heat and precursor gas are very simple in the zones which do not contain the front, since there is no source or sink term for any of them. They are then integrated straightforwardly. In the raw (cold) preform zone, the following relations are given: ⎧ Lc = Lt − Lh − xf , ⎪ ⎪ ⎪ c (Ti − Te ) ⎪ ⎪ q= , ⎪ ⎪ ⎨ Lc (46) D0 m+1 (Ce − Cc ) c ⎪ J = , c ⎪ ⎪ L ⎪ c ⎪ ⎪ ⎪ ⎩ Ce = Ptot /(RTe ), Te = Tb + q/ h, where c is the conductivity on the cold side, (1 − c )0 . Inside the front, the following equations hold: ⎧ R ⎪ , ⎪ Jc = (J˜c )Cc Ak(Th )0 Th2 ⎪ ⎪ E aq ⎪ ⎨ R v = Vm Cc Ak(Th )0 Th2 , (47) ⎪ E aq ⎪ ⎪ ⎪ ⎪ ⎩ xf = f 0 T 2 R . h Ea q Values for the dimensionless parameters are provided by the preceding local front study with m = 1 and n = 23 ; they are listed in Table 7.

On the dense (hot) side, when lh = 0, i.e., when t = 0, the computation has been initialized with the following equations: ⎧T 0 = T , s h ⎪ ⎪ ⎪ 0 = L − x0 , ⎪ L ⎪ t c ⎪ f ⎪ ⎪ ⎨  (T − Te ) c i q0 = , (48) ⎪ Lt − xf0 ⎪ ⎪ ⎪ ⎪ ⎪ D m+1 (Ce − Cc ) ⎪ ⎪ ⎩ Jc0 = 0 c , Lt − xf0 where the superscript 0 indicates initial values. At this initial step, Cc0 , v0 , and xf0 verify Eqs. (46) and (47), with fluxes (q, Jc ) set to (q 0 , Jc0 ). Then, a time integration is performed over N time intervals of size
t, and in an explicit fashion. The ith integration step is performed on Lh i−1 Lih = Li−1
t. h +v

(49)

All other values are updated by the resolution of two subsystems: ⎧ i Te = Tb + q i / h, ⎪ ⎪ ⎪ ⎪ ⎪ Ti − Tei ⎪ i = ⎪ q , ⎪ c ⎪ ⎪ Lt − xfi − Lih ⎨ (50) i 2 ⎪ i = 0 R(Th ) ·  , ⎪ x ⎪ f ⎪ f ⎪ Ea q i ⎪ ⎪ ⎪ i i ⎪ l q ⎪ ⎩ T i = Ts − h , h h ⎧ i Ce = Ptot /(RTei ), ⎪ ⎪ ⎪ ⎪ i 2 ⎪ ⎨ i 0 R(Th ) Vm (Cci )2 Ak(Thi ) · J˜c , Jc = (51) Ea q i ⎪ ⎪ i ⎪ ⎪ Jc ⎪ ⎩ Cci = Cei − (lt − xfi − lhi ). D0 m+1 c The last step in the ith iteration is to evaluate the new value of vi from Eq. (47). This scheme has the advantage of describing the process either in 1D cartesian, 1D cylindrical, or 1D radial coordinates, if the curvature radius is large compared to the thickness of the front. If this last hypothesis is not verified, one has to modify

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

7525

14

20

12

Full model

18

Exptl.

10 Quantities

Radial Position (mm)

Simplified model

16

8 6 4

14

2 0

12

0

1

2

3

4

5

Time (h)

10

Fig. 21. Computed front width, front velocity and heat flux, as compared to Vignoles et al. (2006).

8 0

1

2

3

4

5

Time (h) Fig. 20. Front position vs. time: simulation results and experimental values.

the local front study in order to incorporate explicitly the effect of curvature in the evaluations. Since the scheme is explicit in essence, the stability of the computations with respect to the time step size is always checked. 3.2. Results and validation The numerical integration has been performed and compared to experimental data, using the front position at different times as a validation parameter. Fig. 20 is a comparison plot showing reasonable agreement between experiments, this model and the previous model (Lines et al., 2005; Vignoles et al., 2006), which featured much more details (multicomponent diffusion, Knudsen effects, dependence of gas diffusion and solid heat conduction on temperature, etc., …) and no assumption about the existence of a front. The fact that the simplified model even looks better than the full model is purely incidental. Indeed, either the simplified model or full model could be tuned up for a better agreement, but such a systematic optimization has not been looked for. Fig. 21 is a comparison of front width, front velocity, and heat flux estimations performed with both models. There is a quantitative agreement between both models concerning velocity; the width is 15% overestimated in the simplified model, but with a totally similar evolution in time. Also, the heat flux is in reasonable agreement but its time evolution is too fast in the simplified model—possibly because it does not incorporate the effect of the front curvature. All models diverge from the actual experimental data in the end of the run, because the heat transfer coefficient drops, as film-boiling begins to occur around the whole preform, instead of being localized only in

the interior of the preform. However, this shows that the approach presented in this paper is a successful reduction of the more detailed model even if it is not able to describe the onset of the front.

4. Conclusion The fabrication of carbon/carbon composites by the “filmboiling” or “Kalamazoo” technique, derived from thermalgradient CVI, has been modelled using the concept of densification front. Simplified one-dimensional equations have allowed to perform a local study of the front, yielding useful results: • Existence criteria for three types of front as a function of process parameters and of the porous preform properties; • Description of front behavior in terms of velocity, width, and residual porosity or concentration left behind it, as a function of all parameters; One should notice that the described densification front is somewhat different in nature and behavior from some well-known flame fronts, since it features at least three variables, and the role of temperature is totally distinct: here, this quantity acts more like a reactant in combustion. Of the three types of front, one is achieved in practice, with a complete depletion of reactant and a moderate residual porosity. The presence of this latter quantity which is a priori unknown implies the use of a shooting algorithm in a light numerical (or semi-analytical) scheme. The asymptotic limits in parameter space have been also studied analytically. The semi-analytical results for the front behavior are reused in a larger-scale model describing the infiltration of the whole preform, and the computed results compare favorably with experimental data and with a previous, more detailed, numerical simulation.

7526

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527

Notation A c C D Ea f F g h J J˜ k k0 K K L m n p q r R R t T v vg Vm x Y

characteristic value for surface area, m−1 specific heat, J kg−1 K −1 precursor gas-phase concentration, mol m−3 precursor gas mass diffusion coefficient, m2 s−1 activation energy, J mol−1 thermal dependence of reaction rate, dimensionless primitive of function f “infiltrability function”, dimensionless heat transfer coefficient, W m−2 K −1 precursor gas mole flux, mol m−2 s−1 scaled precursor gas mole flux, dimensionless deposition rate constant, m (mol m−3 )(1−) K − s−1 pre-exponential term of Arrhenius equation, m (mol m−3 )(1−) s−1 scaled total mole flux through the front, dimensionless permeability, m2 length, m Archie’s law exponent, dimensionless surface area variation exponent, dimensionless total pressure, Pa heat flux, W m−2 K −1 radial coordinate, m perfect gas constant, J mol−1 K −1 precursor consumption rate, mol m−3 s−1 time, s temperature, K front velocity, m s−1 local gas velocity, m s−1 deposit molar volume, m3 mol−1 space coordinate, m densification parameter, dimensionless

Greek letters   

     
v  Υ

reaction rate order, dimensionless scaled activation energy, dimensionless porosity, dimensionless scaled temperature, dimensionless condensation ratio, dimensionless thermal conductivity, W m−1 K −1 viscosity, Pa s thermal expansion parameter, dimensionless scaled space coordinate, dimensionless density, kg m−3 internal surface area, m−1 characteristic time, s scaled concentration gradient, dimensionless Thiele modulus, dimensionless scaled concentration, dimensionless scaled front velocity, dimensionless

Subscripts 0 0 b

refers to initial state (raw material, pure gas, …) refers to initial time in large-scale simulation boiling

c f g h i lim ref s t

refers to the cold side of the front refers to the front gas refers to the hot side of the front ignition limiting value reference value value at susceptor location total

Acknowledgments The authors wish to thank Région Aquitaine for financial support to N.N. through a Ph.D. grant, and Jan-Bouwe van den Berg (Vrije Universiteit Amsterdam) for support in utilizing AUTO2000 software. Numerous interesting discussions with G. Damamme, P. David, and C. Gachet (CEA, Le Ripault) are also gratefully acknowledged. References Besmann, T.M., Sheldon, B.W., Lowden, R.A., Stinton, D.P., 1991. Vapor-phase fabrication and properties of continuous-filament ceramic composites. Science 253, 1104–1109. Burganos, V.N., 1998. Gas diffusion in random binary media. Journal of Chemical Physics 109 (16), 6772–6779. Coindreau, O., Vignoles, G.L., 2005. Assessment of geometrical and transport properties of a fibrous C/C composite preform as digitized by X-ray computed micro-tomography. Part I: image acquisition and geometrical properties. Journal of Materials Research 20, 2328–2339. Delhaès, P., 2003. Chemical vapor deposition and infiltration processes of carbon materials. Carbon 40 (5), 641–657. Doedel, E.J., 1981. Auto, a program for automatic bifurcation analysis of autonomous systems. Congressus Numerantium 30, 265–384. Doedel, E.J., Paffenroth, R.C., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Yu.A., Oldeman, B.E., Sandstede, B., Wang, X., 2002. Auto 2000: continuation and bifurcation software for ordinary differential equations (with HomCont). Technical report, CalTech. Dosanjh, S.S., Pagni, P.J., Fernandez-Pello, A.C., 1987. Forced cocurrent smoldering combustion. Combustion and Flame 68, 131–142. Houdayer, M., Spitz, J., Tran Van, D., 1984. US Patent no. 4 472, 454. Lines, J.-F., Vignoles, G.L., Goyhénèche, J.-M., Puiggali, J.-R., 2005. Thermal modelling of a carbon/carbon composite material fabrication. In: Denis, S. (Ed.), Proceedings of the International Conference on Thermal Process Modelling and Computer Simulation 2003, Les Ulis, France, EDP Sciences, Journal de Physique IV France 120, 291–297. Melkote, R.R., Jensen, K.F., 1990. A model for chemical vapor infiltration of fibrous substrates. In: Besmann, T.M., Gallois, B.M. (Eds.), Chemical Vapor Deposition of Refractory Metals and Ceramics, Pittsburgh, Pennsylvania, Materials Research Society, Materials Research Society Symposium Proceedings 168, 67–72. Nadeau, N., Brauner, C.-M., Vignoles, G.L., 2004. Modélisation de front de densification de milieux poreux. In: Thovert, J.-F., Quintard, M. (Eds.), Combustion, Pyrolyse des milieux poreux (Actes de la journée du 28 Janvier 2004 à Paris), VandWuvre, France, SFT. Naslain, R., Langlais, F., 1990. Fundamental and practical aspects of the chemical vapor infiltration of porous substrates. High Temperature Science 27, 221–235. Rikvold, P.A., Stell, G., 1985a. d-dimensional interpenetrable-sphere models of random two-phase media: microstructure and an application to chromatogaphy. Journal of Colloid Inter Science 108, 158–162. Rikvold, P.A., Stell, G., 1985b. Porosity and specific surface for interpenetrable-sphere models of two-phase random media. Journal of Chemical Physics 31, 1014–1020.

N. Nadeau et al. / Chemical Engineering Science 61 (2006) 7509 – 7527 Thiele, E.W., 1939. Relation between catalytic activity and size of particules. Industrial Engineering Chemistry 31, 916. Tomadakis, M.M., Robertson, T.J., 2005. Viscous permeability of random fiber structures: comparison of electrical and diffusional estimates with experimental and analytical results. Journal of Composite Materials 39 (2), 163–188. Tomadakis, M.M., Sotirchos, S.V., 1991. Knudsen diffusivities and properties of structures of unidirectional fibers. A.I.Ch.E. Journal 37, 1175–1186. Tomadakis, M.M., Sotirchos, S.V., 1993. Ordinary, transition and Knudsen regime diffusion in random capillary structures. Chemical Engineering Science 48 (19), 3323–3333.

7527

Vignoles, G.L., Nadeau, N., Brauner, C.-M., Lines, J.-F., Puiggali, J.-R., 2005. The notion of densification front in CVI processing with temperature gradients. In: Zhu, D., Lara-Curzio, E., Kriven, W.M. (Eds.), Mechanical Properties and Performance of Engineering Ceramics and Composites, Westerville, OH, The American Ceramic Society, Ceramic Engineering and Science Proceedings 26, 187–195. Vignoles, G.L., Goyhénèche, J.-M., Sébastian, P., Puiggali, J.-R., Lines, J.-F., Lachaud, J., Delhaès, P., Trinquecoste, M., 2006. The filmboiling densification process for C C composite fabrication: from local scale to overall optimization. Chemical Engineering Science 61, 5336–5353.