Numerical Modeling Chemical Vapor Infiltration of SiC Composites

Abstract

The multiphase field model for chemical vapor infiltration (CVI) of SiC/SiC composites is developed in this study, thereby to reproduce the microstructure evolution during CVI process and to achieve better understanding of the effect of process parameters (e.g., temperature, pressure, etc.) on the final product. In order to incorporate the thermodynamics of methyltrichlorosilane (MTS) pyrolysis into phase field model framework, the reduced chemical reaction mechanism is adopted. The model consists of a set of nonlinear partial differential equations by coupling Ginzburg-Landau type phase field equations with mass balance equations (e.g., convection-diffusion equation) and the modified Navier-Stokes equations which accounts for the fluid motion. The microstructure of preferential codeposition of Si, SiC under high ratio of H₂ to MTS is simulated and the potential risk of blockage of the premature pores during isothermal CVI process is predicted. The competitive growth mechanism of SiC grains is discussed and the formation process of potential premature pore blockage is reproduced.

1. Introduction

The lightweight SiC fiber-reinforced SiC matrix composite (SiC/SiC) is of great importance in industry because of its excellent thermophysical and mechanical properties, for example, high strength, light weight, and excellent chemical-thermal resistance [1, 2]. It is expected to be served as structural and functional components, especially in extreme environments, such as for gas turbine, spacecrafts, and nuclear fusion reactors. This composite is now industrially widely produced by chemical vapor infiltration (CVI) process, a unique method for producing continuous fiber ceramic composite which belongs to chemical vapor deposition (CVD) on the internal surface of a porous preform. Among the different precursors for the deposition of SiC, methyltrichlorosilane (MTS) is the most frequently used precursor due to the equivalent ratio of Si to C, wide deposition temperature range, and low toxicity [3, 4]. If MTS is the only precursor for manufacturing of the SiC matrix by infiltration, SiC with excess Si, SiC with excess C, and stoichiometric SiC can be deposited by pyrolysis of MTS as shown in Figure 1.

(a)

(b)

Figure 1 

Scanning electron microscope (SEM) image of the ceramic infiltrated with SiC by pyrolysis of MTS [5].

Extensive experimental works have been devoted to investigate the influence of the infiltration conditions on SiC deposition, for example, the temperature, pressure, and the kind of precursor. Many studies also concentrated on the morphology of the deposited matrix. Most of them are focused on the composition of the deposited matrix in terms of the C/Si atomic ratio caused by macroscopic infiltration conditions [5–7]. Considering that the experimental investigation is time consuming and costly, more and more studies are now focused on numerically modeling chemical vapor infiltration process with the development of the computational technology and advanced numerical algorithms. Based on the conception of porosity evolution, both the chemical reaction mechanism and mass transport are taken into account in the model [8–10].

Besides the macroscale porosity, the morphology of the deposited matrix and microscale components play a critical role in determining the final properties of the composites. There are mainly two kinds of methodologies to describe the microstructure evolution on mesoscale: sharp interface method and diffusive interface method. The sharp interface method describes the interface between different phases or species explicitly by the mathematical analytical model, which always results in the complicated refinement of the mesh near free moving boundary during numerical discretization and is nearly impractical 3D cases. As an example of diffusive interface method, the level-set method uses the convection equation to implicitly capture the evolution of microstructures, which may be numerically instable. Recently, the phase field method has become an important and extremely versatile technique for modeling microstructure evolution on mesoscale [11]. The phase field method inherits the advantages of the level-set method and the diffusive region is introduced to replace the sharp interface in the sharp interface method. The phase field order parameter is introduced to denote the thermodynamics state there. On the other hand, the convection diffusion equation is applied to describe the evolution of the interface to avoid the numerical problems in the level-set method. In the past two decades, the phase-field method has exhibited its potential in different applications.

Until now, few studies investigating the morphological phenomena during the CVI process have been published. Palmer [12, 13] has studied the morphological instabilities during CVD process. Jin et al. [14, 15] tried to apply the level-set method to predict the evolution of the microstructures during CVI process. There is still a challenge in the morphological study of the CVI process. In the our previous research [16, 17], the phase field model has been successfully developed and applied in predicting the microstructure evolution during CVI process of C/C and SiC/SiC composites in two dimensions. Due to the codeposition of condensed solid phases 𝛽-SiC(s) with Si(s) and/or C(graphite), a multiphase field model is formulated in this paper to model the CVI process of SiC fiber-reinforced SiC matrix composite (SiC/SiC).

2. Multiphase Field Model of the Chemical Vapor Infiltration

2.1. Phase Field Equations

In phase field model, the spatial diffusive distribution of phase field order parameter 𝜙 is used to denote the gas-solid interface implicitly instead of the explicit sharp interface (shown in Figure 2). Therefore the spatiotemporal behavior of phase field order parameter reproduces the deposited SiC microstructure evolution. Due to the fact that the codeposition of condensed solid phases occurs under specific processing conditions during the SiC CVI process, the multiphase field order parameters are adopted to describe not only the gas-solid but also the solid-solid interfaces. The existence of different bulk phase 𝑖 is represented by phase field order parameters (𝜙_𝑖,𝑖=1…𝑁,𝑁 is the number of phases or species) to be 1, otherwise 0 in the rest phases, and 0<𝜙_𝑖<1 in the interfaces between 𝑖th and the other phases. The sum of these phase field order parameters should be conversed as the following equation: ^𝑁_𝑖=1𝜙_𝑖=1, (1) here 𝜙₁ denotes the multi-component gaseous phase and 𝜙₂, 𝜙₃, 𝜙₄ denote silicon, silicon carbide, and carbon, respectively.

(a) Scheme in the fiber axis

(b) Scheme in the fiber vertical direction

Figure 2 

Schematic illustration of the definition of the phase field parameter 𝜙.

Benefited from the equivalent ratio of silicon to carbon, methyltrichlorosilane (CH₃SiCl₃, MTS) is adopted as precursor in the isothermal chemical vapour infiltration (ICVI) process of SiC/SiC composites. During the process, MTS as precursor and hydrogen (H₂) as catalytic gas and as diluting gas are pumped into the reactor. The deposition production (SiC with Si and/or C) from MTS depends on the temperature, pressure and the ratio of H₂/MTS. In the previous work [16, 18], the lumped reaction is adopted to predict the dissociation of MTS using the Arrhenius’ equation for the dependence of the reaction rate on temperature. For this lumped reaction, it still seems not being well understood which is indicated by the following facts: the reaction order with respect to pyrolysis of MTS has been reported varying from −3 to 2.5 and activation energies from 40 kJmol⁻¹ to 400 kJmol⁻¹ [19] depending on the infiltration condition. Now experimental research increasingly focuses on the reduced mechanism of thermal decomposition of MTS. Instead of simplified morphological system free energy description under lumped reaction in previous research, the reduced mechanism permits a detailed description which is necessary for constructing phase field model. On the other side, only several species are considered in the reduced mechanism instead of hundreds of species in equilibrium gas-phase reaction thermodynamics calculation of CVD/CVI system.

In the following, based on the work by Zhang and Huttinger [19], Allendof and Melius [20], and Fitzer and Kehr [21] a model of both homogeneous gas phase chemistry and heterogeneous surface chemistry is presented. Here CH₃ and C₂ species are assumed to be the carbon sources, while SiCl₃ and HSiCl₃ are assumed to be the silicon sources for the deposition of C, Si, and SiC.

Dissociation reaction of MTS is as follows: CH₃SiCl₃+H₂⟷SiCl₃+CH₃+H₂ (2) C-bearing species is as follows: CH₃+H₂⟷CH₄+HCH₃+CH₃⟷C₂H₄+H₂ (3)(4) Si-bearing species is as follows: SiCl₃⟷SiCl₂+ClSiCl₃+H₂⟷HSiCl₃+H (5)(6) Formation of the depositions (∗ indicates the adsorption of the intermediate species on free surface site) is as follows: CH₃^∗+3H⟷C^∗+3H₂C₂H₄^∗+4H⟷3C^∗+4H₂SiCl₃^∗+3H⟷Si^∗+3HClHSiCl^∗₃+2H⟷Si^∗+3HClSi^∗+C^∗⟷SiC (7)(8)(9)(10)(11) The above chemical reactions are governed by thermodynamics and kinetics. The thermodynamics defines the driving force which indicates the potential direction that the reaction will proceed. The kinetics defines the transport process and subsequently determines the rate-control step during chemical vapor deposition. Both of these two aspects of information should be introduced into the multiphase field model. The free energy changing for a reaction is defined as Δ𝐺_𝑟=𝑧^𝑖Δ𝐺⁰_𝑓⋅𝑖+𝑅𝑇ln𝑄, (12) where 𝑧^𝑖 is the stoichiometric coefficient of species 𝑖 in reaction (negative for reactants, positive for products). 𝑅 (J⋅mol⁻¹⋅K⁻¹) is the ideal gas constant. 𝑄 is the reaction quotient which is related to the activity and partial pressure of species 𝑖; 𝐺⁰_𝑓⋅𝑖 is the standard free energy of formation of species 𝑖 at temperature 𝑇(K) and 1 atm as Δ𝐺⁰_𝑓⋅𝑖=Δ𝐻_𝑖+𝑇Δ𝑆_𝑖, (13) where 𝐻_𝑖=𝑅𝑇𝑎₁+𝑎₂2𝑎𝑇+³3𝑇³+𝑎₄4𝑇⁴+𝑎₅5𝑇⁵+𝑎₆𝑇,𝑆_𝑖=𝑅𝑎₁ln𝑇+𝑎₂𝑎𝑇+³2𝑇²+𝑎₄3𝑇⁴+𝑎₅4𝑇⁴+𝑎₇, (14)(15) with coefficient 𝑎_𝑖,𝑖=1…7 being obtained for each species from the thermodynamic database in the NASA equilibrium code. The standard thermodynamic data, that is, the enthalpies 𝐻 and the entropies 𝑆 at 298 K can be found in JANAF thermochemical tables [22]. The thermodynamics properties of individual substances can be refered to [23].

The Ginzburg-Landau type of system free energy can be expressed as 𝐹=_Ω^𝑁_{𝑗>𝑖𝑖=1}−𝜀²_𝑖𝑗2∇𝜙_𝑖∇𝜙_𝑗+𝑊_𝑖𝑗𝜙_𝑖𝜙_𝑗+𝑓_𝑒+𝜆_𝐿^𝑁_𝑖=1𝜙_𝑖−1𝑑Ω, (16) here 𝜀_𝑖𝑗 is the gradient energy coefficient and 𝑊_𝑖𝑗 is the height of the energy barrier. 𝜀_𝑖𝑗 is dependent on the interface tensor and the interface thickness 𝜉 [11]. Anisotropy of the interface energy is incorporated in the model by the dependence of 𝜀_𝑖𝑗 on the orientation [24] (𝜙_𝑖, ∇𝜙_𝑖). The Lagrange multiplier 𝜆_𝐿 is introduced to account for that the sum of phase field in the system is conserved. The first item with the gradient of phase fields accounts for energy penalty for the formation of interfaces/boundaries. The second item is the energy barrier to prevent phase transition occurring randomly. 𝑓_𝑒 is the chemical free energy density of local bulk phase, which can be interpolated by free energy density of the coexisting phases 𝐺^𝑖(𝐜,𝑇) using the weight function 𝑝(𝜙_𝑖): 𝑓_𝑒=^𝑁_𝑖𝑝𝜙_𝑖Δ𝐺^𝑖_𝑟 (17) here the weight function takes the form 𝑝(𝜙_𝑖)=𝜙_𝑖 to ensure the high numeric efficiency. Thereby any point (𝑥,𝑡) in the system is assumed to be a mixture of 𝑁 phases with fraction of 𝜙_𝑖 for 𝑖 phase. In previous work [16, 18], the 𝐺_𝑖 is phenomenologically expressed by 𝐺_𝑖=𝐴𝑐_MTS+𝐵𝑐_HCl for gas phase and 𝐺_𝑖=𝐶𝑐_SiC for solid phase, where 𝐴, 𝐵, and 𝐶 are phenomenological constants. These two expressions also indicate that the driving force of chemical reaction is change of the Gibb’s free energy which is accordant with thermodynamic principle; however, they fail to indicate that that Gibb’s free energy is highly dependent on temperature and pressure, not only on the concentration of species.

The evolution of phase order parameter fields follows the time-dependent Ginzburg-Landau equations (TGDL) assuming minimization of total system free energy as follow: 𝜕𝜙_𝑖=𝜕𝑡_{𝑛𝑗≠𝑖}−2𝑠_𝑖𝑗𝑀_𝑖𝑗𝑛𝜕𝐹𝜕𝜙_𝑖−∇𝜕𝐹𝜕∇𝜙_𝑖−𝜕𝐹𝜕𝜙_𝑗−∇𝜕𝐹𝜕∇𝜙_𝑗=^𝑛_𝑗≠𝑖−2𝑠_𝑖𝑗𝑀_𝑖𝑗𝑛×^𝑛_𝑗≠𝑖𝜀²_𝑖𝑗2∇²𝜙_𝑗−∇²𝜙_𝑖+𝑊_𝑖𝑗𝜙_𝑗−𝜙_𝑖+Δ𝐺_𝑖−Δ𝐺_𝑗, (18) where 𝑠_𝑖𝑗 is the identification function to locate interface between 𝑖 phase and 𝑗 phase. By interface field theory [25], 𝑠_𝑖𝑗=𝑠_𝑖𝑠_𝑗 where the step function 𝑠_𝑖(𝑥,𝑡)=1 if 𝜙_𝑖>0 and 𝑠_𝑖(𝑥,𝑡)=0 otherwise. The interface field theory assumes that the junction is composed of the multipairwise interfaces and thereby the Lagrange multiplier limit is canceled. 𝑀_𝑖𝑗 is the interface mobility coefficient which includes kinetics information of the reaction should come from experiments. 𝑛 is the number of phases really coexisting in a given point and can be expressed as 𝑛(𝑥,𝑡)=^𝑁_𝑖𝑠_𝑖(𝑥,𝑡). (19)

2.2. Mass Transport

In the following, the mass transport equation in the multiphase field model will be deduced. The modified mass transport of the species in the process of ICVI can be stated in the compact form of a time-dependent convection-diffusion problem: 𝜕_𝑡𝑐_𝑖=−𝜙₁(𝐮⋅∇)𝑐_𝑖+∇⋅𝐷∇𝑐_𝑖+𝑟^eff_𝑖,𝑖∈gaseousspecies, (20) where 𝑐_𝑖 (mol⋅m⁻³) is the molar concentration of gaseous species 𝑖. 𝑖∈{gaseousspecies} indicates that the diffusion of each species in solid phases is disregarded due to its low diffusivity. 𝐮 (m⋅s⁻¹) is the velocity vector. 𝐷 (m²⋅s⁻¹) is the diffusivity which is a function of effective binary diffusion coefficient 𝐷_𝑖 in gas. 𝑟_𝑖 (mol⋅m⁻³s⁻¹) is the effective chemical reaction rate. In this paper the Milke equation [26] is used to calculate the effective binary diffusion coefficient 𝐷_𝑖: 𝐷_𝑖=1∑_𝑗≠𝑖𝛾_𝑗/𝐷_𝑖,𝑗, (21) where 𝛾_𝑗 is the molar fraction of component 𝑗. 𝐷_𝑖,𝑗 is the binary diffusion coefficient of species 𝑖 and 𝑗, which can be calculated by Chapman-Enskog equation: 𝐷_𝑖,𝑗=3162(𝑅𝑇)³𝜋1𝑀_𝑖+1𝑀_𝑗1𝑁𝑝𝜎²_𝑖𝑗Ω_𝐷^𝑖𝑗 (22) here 𝑀_𝑖 (g⋅mol⁻¹) is the molecular weight of the species 𝑖, 𝑁 (mol⁻¹) is the Avogadro number, 𝑝(Pa) is the pressure of the gas mixture, 𝜎_𝑖𝑗 is the collision diameter of species 𝑖 and 𝑗, and Ω_𝐷^𝑖𝑗 is the dimensionless collision integral for species 𝑖, 𝑗. Using mixture rule, then the diffusivity 𝐷 can be expressed as: 𝐷=^𝑁_𝑖𝜙_𝑖𝐷_𝑖, (23) where an artificial small diffusivity is assumed for solid phase in order to prevent 𝐷 from being zero in the diffusive interface near solid phase side.

Two kinds of reactions occur during the ICVI process of MTS: homogeneous gaseous reaction and heterogeneous surface reaction. The effective chemical reaction rate 𝑟^eff_𝑖 can be regarded as a combination of the homogeneous gaseous reaction and heterogeneous surface reaction with the reaction intensities (𝜙₁) in gas phase and (1−𝜙₁) in the diffuse interface, respectively. It is assumed that the heterogeneous process ends in the full-deposited solid region Ω_solid (shown in Figure 2). To account for this property of the process, a step function 𝛿(𝜙_𝑖) is introduced in present work: 𝛿𝜙_𝑖=0𝜙_𝑖≥0.99,𝑖≠11othercases. (24) Based on the above function, we can describe the effective chemical reaction rate in the diffuse domain of species 𝑖 as 𝑟^eff_𝑖=𝜙₁𝑟^hom_𝑖+1−𝜙₁𝛿𝜙_𝑖𝑟^het_𝑖, (25) with 𝑟^hom_𝑖 as the homogeneous reaction rate (gaseous reaction rate) of species 𝑖 and 𝑟^het_𝑖 as the heterogeneous reaction rate (surface reaction rate) of species 𝑖.

2.3. Fluid Motion

In the previous work [16, 18], the solid phase is assumed to be Newtonian fluid with a large viscosity. Thereby the system consists of two Newtonian fluids: gas fluid and solid fluid, while, in this work, the solid-fluid system is assumed to be rigid and stationary. Assigning a zero-velocity to the solid, the calculation for solid phase is avoided instead of still solving the equations in previous work. Here a modified Navier-Stokes equations is employed to describe the evolution of the gaseous fluid motion: ∇⋅𝜙₁𝐮𝜌𝜕𝜙=0,₁𝐮𝜙𝜕𝑡=−𝜌₁𝐮⋅∇𝐮+𝜙₁𝐈∇𝜙∇⋅−𝑝+𝜂∇⋅₁𝐮+∇𝜙₁𝐮^𝑇+𝐵𝜂𝜙²₁1−𝜙₁𝐮, (26)(27) here 𝐮 (m⋅s⁻¹) is the velocity. 𝜌 (kg⋅m⁻³) is the density (the same density assumed for different gaseous species for simplicity), and 𝜂 (Pa⋅s) is the viscosity. 𝐈 is the unit metric and 𝑇 is the inversion sign. 𝐵 is the parameter depending on the gradient energy coefficient 𝜀. The last item in (27) is the force item, which is chosen such that the no slip condition can be accurately reproduced regardless of the diffusive interface thickness [27].

2.4. Multiphase Field Model of CVI

Combining (18), (20), (26), and (27), the formulation of the multiphase field model of CVI can be expressed as ∇⋅𝜙₁𝐮𝜌𝜕𝜙=0,₁𝐮𝜕𝑡=𝐹₁𝜙𝝓,𝐮,𝑝=−𝜌₁𝐮⋅∇𝐮+𝜙₁𝐈∇𝜙∇⋅−𝑝+𝜂∇⋅₁𝐮+∇𝜙₁𝐮^𝑇+𝐵𝜂𝜙²₁1−𝜙₁𝐮,𝜕𝜙_𝑖𝜕𝑡=𝐹₂()=𝝓,𝐜_{𝑛𝑗≠𝑖}−2𝑠_𝑖𝑗𝑀_𝑖𝑗𝑛×^𝑛_𝑗≠𝑖𝜀²_𝑖𝑗2∇²𝜙_𝑗−∇²𝜙_𝑖+𝑊_𝑖𝑗𝜙_𝑗−𝜙_𝑖+Δ𝐺_𝑖−Δ𝐺_𝑗,𝜕_𝑡𝑐_𝑖=𝐹₃(𝝓,𝐜)=−𝜙₁(𝐮⋅∇)𝑐_𝑖+∇⋅𝐷∇𝑐_𝑖+𝑟^eff_𝑖,𝑖∈gaseousspecies, (28) where vector 𝝓=(𝜙₁,…,𝜙_𝑁) and 𝐜=(𝑐₁,…,𝑐_𝑀).

2.5. Finite Element Implementation for Multiphase Field Model

The numerical discretization for the phase field formulation consists of two parts: temporal discretization and spatial discretization. In this work, we use the implicit Euler finite difference approximation for temporal discretization and Galerkin weighted residual method for spatial discretization. The implicit Euler format discretization is of the advantage of stability and loose limit on time step: ∇⋅𝜙^𝑛₁𝐮^𝑛+1=0,𝜌𝜙^𝑛₁𝐮^𝑛+1−𝐮^𝑛Δ𝑡=𝐹₁𝝓^𝑛,𝐮^𝑛+1,𝑝^𝑛+1,𝜙_𝑖^𝑛+1−𝜙^𝑛_𝑖Δ𝑡=𝐹₂𝝓^𝑛+1,𝐜^𝑛,𝑐_𝑖^𝑛+1−𝑐^𝑛_𝑖Δ𝑡=𝐹₃𝝓^𝑛,𝐜^𝑛+1. (29)

In the following, we will propose the weak formulations for (28). Based on the Green formula, the solutions of phase field model should satisfy the following equations for the arbitrary test functions 𝜙_𝑖∈𝑉_𝑚, 𝑐_𝑙̃∈𝑉_𝑚, 𝑝∈𝑉_𝑚1, and 𝐮^∼∈𝑉^𝑑_𝑚: _Ω𝑝∇⋅𝜙₁𝐮𝑑𝑣=0_Ω𝜌𝐮^∼⋅𝜕_𝑡𝜙₁𝐮𝑑𝑣+𝐮_Ω^∼∶𝜌𝜙₁𝐮⋅∇𝐮𝑑𝑣+_Ω−𝜙₁𝑝𝐈+𝜂∇𝜙₁𝐮+∇𝜙₁𝐮^𝑇∶∇𝐮^∼−𝑑𝑣_Γ−𝜙₁𝑝𝐈+𝜂∇𝜙₁𝐮+∇𝜙₁𝐮^𝑇∶𝐮^∼−⋅𝐧𝑑𝑎_Ω𝐵𝜂𝜙²₁1−𝜙₁𝐮⋅𝐮^∼𝑑𝑣=0,_Ω𝜙_𝑖𝜕_𝑡𝜙_𝑖𝑑𝑣−_Ω^𝑛_𝑗≠𝑖2𝑠_𝑖𝑗𝑀_𝑖𝑗𝑛^𝑛_𝑗≠𝑖𝜀²_𝑖𝑗2∇𝜙_𝑖∇𝜙_𝑖+𝑑𝑣_{Γ𝑛𝑗≠𝑖}2𝑠_𝑖𝑗𝑀_𝑖𝑗𝑛^𝑛_𝑗≠𝑖𝜀²_𝑖𝑗2∇𝜙_𝑖𝜙_𝑖+⋅𝐧𝑑𝑎_{Ω𝑛𝑗≠𝑖}2𝑠_𝑖𝑗𝑀_𝑖𝑗𝑛×^𝑛_𝑗≠𝑖𝜀²_𝑖𝑗2∇²𝜙_𝑗+𝑊_𝑖𝑗𝜙_𝑗−𝜙_𝑖+Δ𝐺_𝑖−Δ𝐺_𝑗𝜙_𝑖𝑑𝑣,_Ω𝜕_𝑡𝑐_𝑖𝑐_𝑖̃𝑑𝑣+_Ω𝜙₁(𝐮⋅∇)𝑐_𝑖−𝑟^eff_𝑖𝑐_𝑖̃−𝑑𝑣_Ω𝐷∇𝑐_𝑖⋅∇𝑐_𝑖̃𝑑𝑣+_Γ𝐷∇𝑐_𝑖𝑐_𝑖̃⋅𝐧𝑑𝑎=0, (30)(31)(32) with 𝐧 as the outward unit normal.

Here the spatial discretization of (32) is proposed as an example. Defining the approximation 𝑐_𝑖=𝐍^𝑇𝐜_𝐢, where 𝐍^𝑇=(𝑁₁,…,𝑁_𝑛) is vector composed by shape functions, 𝐜_𝐢=(𝑐¹_𝑖,…,𝑐^𝑛_𝑖), 𝑛 is the node number in element. Choosing shape functions as test functions and assuming the boundary conditions 𝜕𝑐_𝑖/𝜕𝑛=𝑞, then we have 𝐊₁𝜕_𝑡𝐜_𝑖+𝐊₂𝐜_𝑖𝐊=𝐟₁=_Ω𝑒𝐍𝐍^𝑇𝐊𝑑𝑣₂=_Ω𝑒𝐷𝜕𝐍𝜕𝑥𝜕𝐍^𝑇+𝜕𝑥𝜕𝐍𝜕𝑦𝜕𝐍^𝑇+𝜙𝜕𝑦₁𝑢𝜕𝐍𝜕𝑥+𝑣𝜕𝐍𝜕𝑦−𝑟^eff_𝑖𝐍^𝑇𝑑𝑥𝑑𝑦𝐟=−_Γ𝐷𝐍𝑞𝑑Γ. (33)

The same procedure can be applied to the rest of equations. Thereby with damped Newton iteration [28, 29], an overall linear system is obtained: 𝐴𝐱=𝐟. (34)

Notice that the coefficient matrix 𝐴 is large, sparse, and nonsymmetric. Here the parallel generalized minimal residual method (GMRES) is employed as an iterative linear solver.

3. Result and Discussion

3.1. Initial and Boundary Conditions

In this paper, the evolution of the microstructure in unidirectional fibre-reinforced composites is studied. The geometric model is shown in Figure 3. For three-dimensional case, the one-fiber or the multifiber model can be used. Due to the symmetry along the radial section of the fiber, this three-dimensional model often can be reduced to a 2D model as illustrated in Figure 3. To accomplish the multiphase field model equations system, appropriate initial and boundary conditions should be supplemented. For the one-fiber or multifiber model in 3D case, the boundaries are composed of (a) the inlet face 𝜕Ω_inlet and outlet face 𝜕Ω_outlet which are normal to the fiber direction; (b) the artificial wall faces 𝜕Ω_wall with same physical properties which are parallel to the fiber direction; (c) the gas/solid boundary 𝜕Ω_g/s which is especially required by fluid motion for no slip boundary condition is defined by 𝜙₁=0.5.

Figure 3 

Schematic illustrations of unidirectional fiber-reinforced substrates.

For the multiphase field order parameters 𝜙_𝑖, there is no flux at the inlet, outlet, and wall, which can be described by the Neumann boundaries: 𝐧⋅∇𝜙_𝑖=0on𝜕Ω_inlet,𝜕Ω_outlet,𝜕Ω_wall. (35)

For the species variable 𝑐_𝑖, they are limited by mass conservation which leads to either a specification of the concentration boundary conditions with Dirichlet type or a prescription of the mass fluxes with generalized Neumann type. The precursor MTS imported at the inlet is assumed to be a constant (Dirichlet type) and there are no mass fluxes for the other species (Neumann type). Compared with convection, the influence of the diffusion at the outlet can be neglected (Neumann type). On artificial wall, there is no mass flux across the boundaries due to the close system, the symmetric boundary conditions are implemented (Neumann type). Then the boundary conditions for species take the following forms: 𝑐_MTS=𝑐_MTS0,𝐮𝑐_𝑖−𝐷∇𝑐_𝑖⋅𝐧=0on𝜕Ω_inlet𝐮𝑐_𝑖−𝐷∇𝑐_𝑖⋅𝐧=0on𝜕Ω_wall𝐷∇𝑐_𝑖⋅𝐧=0on𝜕Ω_outlet. (36)

Here 𝑐_MTS0 is the molar concentration of MTS in the imported gas. It can be calculated by the ideal gas law (𝑝_MTS is the partial pressure of MTS): 𝑐_MTS0=𝑝_MTS/(𝑅𝑇). (37)

For the variables (𝐮,𝑝) in Navier-Stokes equations, the velocity 𝐮 along fiber direction on the inlet is constant while those normal to the fiber direction are zero. The no slip boundary is applied on the gas/solid boundary. On the outlet, the total pressure 𝑝 is fixed as a constant. It is assumed there are no viscous effects hence no boundary layer develops on the wall. Then we have 𝐮=𝐮₀on𝜕Ω_inlet𝐮=0on𝜕Ω_g/s𝑝=𝑝₀on𝜕Ω_outlet𝐧⋅𝐮=0,𝐭⋅𝑝𝐈−𝜂∇𝐮+(∇𝐮)^𝑇𝐧=0on𝜕Ω_wall, (38) where 𝐮₀=(0,0,𝑢₀) and 𝐭 is a tangential vector to the boundary.

The initial conditions for variables are straightforward: 𝜙_𝑖|_𝑡=0=1 if 𝑖th phase exists; 𝑐_𝑖|_𝑡=0=0 for not MTS species; 𝐮|_𝑡=0=𝐮₀ and 𝑝|_𝑡=0=𝑝₀.

The similar procedure can be applied to the 2D case, except that the bottom wall in 3D case now is occupied by fiber as shown in the section of Figure 3. Therefore the corresponding boundary conditions are 𝐧⋅∇𝜙_𝑖=0 and 𝐧⋅∇𝑐_𝑖=0. Notice that it is not necessary to set up boundary condition for the variables in Navier-Stokes equations (𝐮,𝑝), because the actual boundary for these variables come to be the interface with 𝜙₁=0.5, namely, gas/solid boundary 𝜕Ω_g/s. The parameters related to model are listed in Table 1.

Table 1 

Parameters related to model.

3.2. 2D and 3D Simulation

The numerical simulation is performed based on the software “COMSOL Multiphysics” with the relative tolerance 1𝑒−2 and absolute tolerance 1𝑒−4. The distributed parallelization calculation is implemented on high performance computing (HPC) clusters “InstitutsCluster” in Karlsruhe Institute of Technology. Each node contains two Quad-core Intel Xeon processors which run at a clock speed of 2.667 GHz and each node has 16 GB of main memory, 4 local disks (250 GB each).

Figure 4 shows the microstructure evolution and corresponding MTS concentration field during SiC deposition. It can be found that the growth of solid phase nearby precursor inlet side is much faster than that nearby outlet side, which implies that the risk of blockage of the premature pores increases. Also we can find that the C (graphite) phase tends to be eliminated by SiC phase because of its overgrowth. This is accordant with experiments [30, 31] that only codeposition of Si(s) and SiC(s) exists while C (graphite) phase is absent in under the high ratios of H₂ to MTS (>104) as in this study. As time passes by, concentration of MTS nearby the inlet side is much higher than that of the rest, and this distribution favours the growth of the solid phase nearby inlet side. The concentration of MTS at section A and B at different moment is shown in Figure 5. The step increasing in the curves along section B indicates the location of gas/solid interface, with left part in solid phase and the right part in gas phase. It shows that as time goes on, the 𝑐_MTS along both these two sections decreases; the gradient decreasing of 𝑐_MTS along section A is obvious which is mainly caused by chemical reaction as being explained later; there barely exists change of 𝑐_MTS along section B because the effect of diffusion is weak comparing to that caused by chemical reaction.

Figure 4 

Phase field evolution and gas composition distribution at scaled time: (a) 𝑡=0.1; (b) 𝑡=0.5; (c) 𝑡=1; (d) 𝑡=1.5.

Figure 5 

Time dependent 𝑐MTS at section A and B.

By changing the reaction diffusion equation, the effect of fluid flow and chemical reaction on 𝑐_MTS distribution is investigated. Three cases are set up: (a) only diffusion is considered; (b) constant fluid flow is considered in mass transport equation based on case a; (c) chemical reaction in mass transport equation is considered based on case a. Simulation shows that the fluid flow in case b changes MTS concentration field as in case a and tends to reduce the concentration gradient as shown in Figure 6. This homogenization of MTS is in favour of the even deposition to reduce the potential pore blockages which cause the porosity in matrix. After chemical reaction is introduced into the diffusion equation as in case c the concentration nearby inlet side is quickly decreased to small value nearby outlet side of the reactor, which indicates that MTS is highly reactive and is prone to decomposing. This inhomogeneous MTS distribution impedes the even deposition to achieve highly densified production. Because the decomposition of MTS is highly dependent on temperature; the thermogradient CVI (TGCVI) is frequently employed in practice to improve the SiC deposition environment. The phase field modeling TGCVI process will be a part of our further research as scheduled in renewal proposal. From this simulation, it can be found that the flow is in favour of the homogeneous distribution of MTS while this concentration distribution is mainly dominated by chemical reaction.

Figure 6 

The effect of fluid and chemical reaction on 𝑐MTS.

Under appreciated infiltration conditions, the single phase SiC is deposited. Figure 7(a) shows the growth of multi-SiC grains with different initial seed size. It is found that as deposition proceeds, the SiC grains with small size tend to be eliminated by those surrounding grains by the competitive growth. The coarse SiC grains in final microstructure are accordant with experimental results which show that the tiny grain on the deposition surface is gradually coarsen. The 3D simulation of SiC grains growth with uniform seeds is shown in Figure 7(b). The SiC deposition on SiC fiber surface under the constant fluid is simulated as shown in Figure 8. Herein the fibers are unidirectionally preformed as bundle. The gas flow is along 𝑧-axis direction, from left to right. The slices of fiber at different moment show clearly how the premature pore blockage comes into being with narrow entrance and wide end.

(a) SiC deposited as grains

(b) SiC grain in 3D

Figure 7 

SiC deposited as grains.

Figure 8 

Blockage formation during SiC fiber deposited with SiC.

4. Summary and Conclusion

The multiphase field model is formulated in this paper to predict the microstructures evolution during the CVI process of SiC. By this model the codeposition behavior of the multisolid phases is reproduced. The model consists of phase field equations, mass conservation equations, and the modified Navier-stokes equations. Thereby the complex coupling is accounted among the interprocesses of chemical reaction, deposition, diffusion, and fluid flow. A new form of the system free energy density is constructed based on chemical reaction free energy changing, instead of phenomenological description of the dependence of free energy on the concentration of the gaseous reactants in previous research. Both homogeneous process and heterogeneous process are included in the model through process intensities. The modified Navier-Stokes equations for fluid motion is proposed by using phase field order parameter as coupling variable, which successfully avoids the assumption for solid phase as Newtonian fluids; thereby potential computation load is reduced. Based on the analysis of linear system 𝐴𝐱=𝐟 which arises from the finite element discretization of the model, the parallelized implement is applied. The microstructure simulation shows the potential risk of blockage of the premature pores during CVI process. It can be found that the diffusion effect tends to cause the homogeneous distribution of MTS and meanwhile the flow is in favour of it; however, this concentration distribution of MTS is mainly dominated by chemical reaction. Because the pyrolysis of MTS is highly sensitive to temperature, the thermogradient CVI (TGCVI) is expected to improve the SiC deposition environment aiming at the full densification of SiC preform. The quick decreasing of MTS from inlet side to the outlet side is accordant with the experiment investigation. The competitive growth mechanism is reproduced during SiC grain deposition and result shows that the SiC grains with small size tend to be eliminated by those surrounding grains. Its 3D simulation illustrates the formation process of the premature pore blockage characterized by narrow entrance and wide end.

Acknowledgment

The financial support by the German Research Foundation (DFG Schn 245/31-1) is gratefully acknowledged.