# Numerical Modeling Chemical Vapor Infiltration of SiC Composites

Authors: Yaochan Zhu; Eckart Schnack
Journal: Journal of Chemistry
 (2013)

---

## 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 H2 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.


---

## Body

## 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.Scanning electron microscope (SEM) image of the ceramic infiltrated with SiC by pyrolysis of MTS [5].
(a)
(b)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 i is represented by phase field order parameters (ϕi,i=1…N,N is the number of phases or species) to be 1, otherwise 0 in the rest phases, and 0<ϕi<1 in the interfaces between ith and the other phases. The sum of these phase field order parameters should be conversed as the following equation:
(1)∑i=1N‍ϕi=1,
here ϕ1 denotes the multi-component gaseous phase and ϕ2, ϕ3, ϕ4 denote silicon, silicon carbide, and carbon, respectively.Schematic illustration of the definition of the phase field parameterϕ.
(a)
Scheme in the fiber axis
(b)
Scheme in the fiber vertical directionBenefited from the equivalent ratio of silicon to carbon, methyltrichlorosilane (CH3SiCl3, 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 (H2) 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 H2/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−1 to 400 kJmol−1 [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 CH3 and C2 species are assumed to be the carbon sources, while SiCl3 and HSiCl3 are assumed to be the silicon sources for the deposition of C, Si, and SiC.Dissociation reaction of MTS is as follows:(2)CH3SiCl3+H2⟷SiCl3+CH3+H2
C-bearing species is as follows:
(3)CH3+H2⟷CH4+H(4)CH3+CH3⟷C2H4+H2
Si-bearing species is as follows:
(5)SiCl3⟷SiCl2+Cl(6)SiCl3+H2⟷HSiCl3+H
Formation of the depositions (* indicates the adsorption of the intermediate species on free surface site) is as follows:
(7)CH3*+3H⟷C*+3H2(8)C2H4*+4H⟷3C*+4H2(9)SiCl3*+3H⟷Si*+3HCl(10)HSiCl3*+2H⟷Si*+3HCl(11)Si*+C*⟷SiC
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
(12)ΔGr=∑ziΔGf·i0+RTlnQ,
where zi is the stoichiometric coefficient of species i in reaction (negative for reactants, positive for products). R (J·mol-1·K-1) is the ideal gas constant. Q is the reaction quotient which is related to the activity and partial pressure of species i; Gf·i0 is the standard free energy of formation of species i at temperature T(K) and 1 atm as
(13)ΔGf·i0=ΔHi+TΔSi,
where
(14)Hi=RT(a1+a22T+a33T3+a44T4+a55T5+a6T),(15)Si=R(a1lnT+a2T+a32T2+a43T4+a54T4+a7),
with coefficient ai,i=1…7 being obtained for each species from the thermodynamic database in the NASA equilibrium code. The standard thermodynamic data, that is, the enthalpies H and the entropies S 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(16)F=∫Ω∑j>i∑i=1N[{(∑i=1Nϕi-1)}-εij22∇ϕi∇ϕj+Wijϕiϕj+fe+λL(∑i=1Nϕi-1)]dΩ,
here εij is the gradient energy coefficient and Wij is the height of the energy barrier. εij 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 εij on the orientation [24] (ϕi, ∇ϕi). The Lagrange multiplier λL 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. fe is the chemical free energy density of local bulk phase, which can be interpolated by free energy density of the coexisting phases Gi(c,T) using the weight function p(ϕi):
(17)fe=∑iNp(ϕi)ΔGri
here the weight function takes the form p(ϕi)=ϕi to ensure the high numeric efficiency. Thereby any point (x,t) in the system is assumed to be a mixture of N phases with fraction of ϕi for i phase. In previous work [16, 18], the Gi is phenomenologically expressed by Gi=AcMTS+BcHCl for gas phase and Gi=CcSiC for solid phase, where A, B, and C 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:(18)∂ϕi∂t=∑j≠in-2sijMijn[(∂F∂ϕi-∇∂F∂∇ϕi)-(∂F∂ϕj-∇∂F∂∇ϕj)]=∑j≠in-2sijMijn×[∑j≠in[εij22(∇2ϕj-∇2ϕi)+Wij(ϕj-ϕi)]+ΔGi-ΔGj{∑j≠in}],
where sij is the identification function to locate interface between i phase and j phase. By interface field theory [25], sij=sisj where the step function si(x,t)=1 if ϕi>0 and si(x,t)=0 otherwise. The interface field theory assumes that the junction is composed of the multipairwise interfaces and thereby the Lagrange multiplier limit is canceled. Mij is the interface mobility coefficient which includes kinetics information of the reaction should come from experiments. n is the number of phases really coexisting in a given point and can be expressed as
(19)n(x,t)=∑iNsi(x,t).


### 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:(20)∂tci=-ϕ1(u·∇)ci+∇·[D∇ci]+rieff,i∈{gaseousspecies},
where ci (mol·m-3) is the molar concentration of gaseous species i. i∈{gaseousspecies} indicates that the diffusion of each species in solid phases is disregarded due to its low diffusivity. u (m·s-1) is the velocity vector. D (m2·s-1) is the diffusivity which is a function of effective binary diffusion coefficient Di in gas. ri (mol·m-3s-1) is the effective chemical reaction rate. In this paper the Milke equation [26] is used to calculate the effective binary diffusion coefficient Di:
(21)Di=1∑j≠iγj/Di,j,
where γj is the molar fraction of component j. Di,j is the binary diffusion coefficient of species i and j, which can be calculated by Chapman-Enskog equation:
(22)Di,j=3162(RT)3π(1Mi+1Mj)1Npσij2ΩDij
here Mi (g·mol-1) is the molecular weight of the species i, N (mol-1) is the Avogadro number, p(Pa) is the pressure of the gas mixture, σij is the collision diameter of species i and j, and ΩDij is the dimensionless collision integral for species i, j. Using mixture rule, then the diffusivity D can be expressed as:
(23)D=∑iNϕiDi,
where an artificial small diffusivity is assumed for solid phase in order to prevent D 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 raterieff can be regarded as a combination of the homogeneous gaseous reaction and heterogeneous surface reaction with the reaction intensities (ϕ1) in gas phase and (1-ϕ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 δ(ϕi) is introduced in present work:
(24)δ(ϕi)={0ϕi≥0.99,i≠11othercases.
Based on the above function, we can describe the effective chemical reaction rate in the diffuse domain of species i as
(25)rieff=ϕ1rihom+(1-ϕ1)δ(ϕi)rihet,
with rihom as the homogeneous reaction rate (gaseous reaction rate) of species i and rihet as the heterogeneous reaction rate (surface reaction rate) of species i.


### 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:
(26)∇·(ϕ1u)=0,(27)ρ∂(ϕ1u)∂t=-ρ(ϕ1u·∇)u+ϕ1∇·(-p)I+η∇·[∇(ϕ1u)+(∇(ϕ1u))T]+Bη(ϕ12(1-ϕ1)u),
here u (m·s-1) is the velocity. ρ (kg·m-3) is the density (the same density assumed for different gaseous species for simplicity), and η (Pa·s) is the viscosity. I is the unit metric and T is the inversion sign. B 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
(28)∇·(ϕ1u)=0,ρ∂(ϕ1u)∂t=F1(ϕ,u,p)=-ρ(ϕ1u·∇)u+ϕ1∇·(-p)I+η∇·[∇(ϕ1u)+(∇(ϕ1u))T]+Bη(ϕ12(1-ϕ1)u),∂ϕi∂t=F2(ϕ,c)=∑j≠in-2sijMijn×[∑j≠in[εij22(∇2ϕj-∇2ϕi)+Wij(ϕj-ϕi)]+ΔGi-ΔGj{∑j≠in}],∂tci=F3(ϕ,c)=-ϕ1(u·∇)ci+∇·[D∇ci]+rieff,i∈{gaseousspecies},
where vector ϕ=(ϕ1,…,ϕN) and c=(c1,…,cM).


### 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:(29)∇·(ϕ1nun+1)=0,ρϕ1nun+1-unΔt=F1(ϕn,un+1,pn+1),ϕin+1-ϕinΔt=F2(ϕn+1,cn),cin+1-cinΔt=F3(ϕn,cn+1).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 ϕ~i∈Vm, c~l∈Vm, p~∈Vm1, and u~∈Vmd:
(30)∫Ωp~∇·ϕ1udv=0∫Ωρu~·∂t(ϕ1u)dv+∫Ωu~:ρϕ1u·∇udv+∫Ω(-ϕ1pI+η(∇(ϕ1u)+∇(ϕ1u)T)):∇u~dv-∫Γ(-ϕ1pI+η(∇(ϕ1u)+∇(ϕ1u)T)):u~·nda-∫ΩBηϕ12(1-ϕ1)u·u~dv=0,(31)∫Ωϕ~i∂tϕidv-∫Ω∑j≠in2sijMijn[∑j≠in(εij22∇ϕi)∇ϕ~i]dv+∫Γ∑j≠in2sijMijn[∑j≠in(εij22∇ϕi)ϕ~i·n]da+∫Ω∑j≠in2sijMijn×[∑j≠in(εij22∇2ϕj+Wij(ϕj-ϕi))+ΔGi-ΔGj{∑j≠in(εij22∇2ϕj+Wij(ϕj-ϕi))}]ϕ~idv,(32)∫Ω∂tcic~idv+∫Ω(ϕ1(u·∇)ci-rieff)c~idv-∫ΩD∇ci·∇c~idv+∫ΓD∇cic~i·nda=0,
with n as the outward unit normal.Here the spatial discretization of (32) is proposed as an example. Defining the approximation ci=NTci, where NT=(N1,…,Nn) is vector composed by shape functions, ci=(ci1,…,cin), n is the node number in element. Choosing shape functions as test functions and assuming the boundary conditions ∂ci/∂n=q, then we have
(33)K1∂tci+K2ci=fK1=∫ΩeNNTdvK2=∫Ωe[D(∂N∂x∂NT∂x+∂N∂y∂NT∂y)+(ϕ1(u∂N∂x+v∂N∂y)-rieff)NT]dxdyf=-∫ΓDNqdΓ.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)Ax=f.Notice that the coefficient matrixA is large, sparse, and nonsymmetric. Here the parallel generalized minimal residual method (GMRES) is employed as an iterative linear solver.


## 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 i is represented by phase field order parameters (ϕi,i=1…N,N is the number of phases or species) to be 1, otherwise 0 in the rest phases, and 0<ϕi<1 in the interfaces between ith and the other phases. The sum of these phase field order parameters should be conversed as the following equation:
(1)∑i=1N‍ϕi=1,
here ϕ1 denotes the multi-component gaseous phase and ϕ2, ϕ3, ϕ4 denote silicon, silicon carbide, and carbon, respectively.Schematic illustration of the definition of the phase field parameterϕ.
(a)
Scheme in the fiber axis
(b)
Scheme in the fiber vertical directionBenefited from the equivalent ratio of silicon to carbon, methyltrichlorosilane (CH3SiCl3, 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 (H2) 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 H2/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−1 to 400 kJmol−1 [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 CH3 and C2 species are assumed to be the carbon sources, while SiCl3 and HSiCl3 are assumed to be the silicon sources for the deposition of C, Si, and SiC.Dissociation reaction of MTS is as follows:(2)CH3SiCl3+H2⟷SiCl3+CH3+H2
C-bearing species is as follows:
(3)CH3+H2⟷CH4+H(4)CH3+CH3⟷C2H4+H2
Si-bearing species is as follows:
(5)SiCl3⟷SiCl2+Cl(6)SiCl3+H2⟷HSiCl3+H
Formation of the depositions (* indicates the adsorption of the intermediate species on free surface site) is as follows:
(7)CH3*+3H⟷C*+3H2(8)C2H4*+4H⟷3C*+4H2(9)SiCl3*+3H⟷Si*+3HCl(10)HSiCl3*+2H⟷Si*+3HCl(11)Si*+C*⟷SiC
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
(12)ΔGr=∑ziΔGf·i0+RTlnQ,
where zi is the stoichiometric coefficient of species i in reaction (negative for reactants, positive for products). R (J·mol-1·K-1) is the ideal gas constant. Q is the reaction quotient which is related to the activity and partial pressure of species i; Gf·i0 is the standard free energy of formation of species i at temperature T(K) and 1 atm as
(13)ΔGf·i0=ΔHi+TΔSi,
where
(14)Hi=RT(a1+a22T+a33T3+a44T4+a55T5+a6T),(15)Si=R(a1lnT+a2T+a32T2+a43T4+a54T4+a7),
with coefficient ai,i=1…7 being obtained for each species from the thermodynamic database in the NASA equilibrium code. The standard thermodynamic data, that is, the enthalpies H and the entropies S 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(16)F=∫Ω∑j>i∑i=1N[{(∑i=1Nϕi-1)}-εij22∇ϕi∇ϕj+Wijϕiϕj+fe+λL(∑i=1Nϕi-1)]dΩ,
here εij is the gradient energy coefficient and Wij is the height of the energy barrier. εij 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 εij on the orientation [24] (ϕi, ∇ϕi). The Lagrange multiplier λL 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. fe is the chemical free energy density of local bulk phase, which can be interpolated by free energy density of the coexisting phases Gi(c,T) using the weight function p(ϕi):
(17)fe=∑iNp(ϕi)ΔGri
here the weight function takes the form p(ϕi)=ϕi to ensure the high numeric efficiency. Thereby any point (x,t) in the system is assumed to be a mixture of N phases with fraction of ϕi for i phase. In previous work [16, 18], the Gi is phenomenologically expressed by Gi=AcMTS+BcHCl for gas phase and Gi=CcSiC for solid phase, where A, B, and C 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:(18)∂ϕi∂t=∑j≠in-2sijMijn[(∂F∂ϕi-∇∂F∂∇ϕi)-(∂F∂ϕj-∇∂F∂∇ϕj)]=∑j≠in-2sijMijn×[∑j≠in[εij22(∇2ϕj-∇2ϕi)+Wij(ϕj-ϕi)]+ΔGi-ΔGj{∑j≠in}],
where sij is the identification function to locate interface between i phase and j phase. By interface field theory [25], sij=sisj where the step function si(x,t)=1 if ϕi>0 and si(x,t)=0 otherwise. The interface field theory assumes that the junction is composed of the multipairwise interfaces and thereby the Lagrange multiplier limit is canceled. Mij is the interface mobility coefficient which includes kinetics information of the reaction should come from experiments. n is the number of phases really coexisting in a given point and can be expressed as
(19)n(x,t)=∑iNsi(x,t).


## 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:(20)∂tci=-ϕ1(u·∇)ci+∇·[D∇ci]+rieff,i∈{gaseousspecies},
where ci (mol·m-3) is the molar concentration of gaseous species i. i∈{gaseousspecies} indicates that the diffusion of each species in solid phases is disregarded due to its low diffusivity. u (m·s-1) is the velocity vector. D (m2·s-1) is the diffusivity which is a function of effective binary diffusion coefficient Di in gas. ri (mol·m-3s-1) is the effective chemical reaction rate. In this paper the Milke equation [26] is used to calculate the effective binary diffusion coefficient Di:
(21)Di=1∑j≠iγj/Di,j,
where γj is the molar fraction of component j. Di,j is the binary diffusion coefficient of species i and j, which can be calculated by Chapman-Enskog equation:
(22)Di,j=3162(RT)3π(1Mi+1Mj)1Npσij2ΩDij
here Mi (g·mol-1) is the molecular weight of the species i, N (mol-1) is the Avogadro number, p(Pa) is the pressure of the gas mixture, σij is the collision diameter of species i and j, and ΩDij is the dimensionless collision integral for species i, j. Using mixture rule, then the diffusivity D can be expressed as:
(23)D=∑iNϕiDi,
where an artificial small diffusivity is assumed for solid phase in order to prevent D 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 raterieff can be regarded as a combination of the homogeneous gaseous reaction and heterogeneous surface reaction with the reaction intensities (ϕ1) in gas phase and (1-ϕ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 δ(ϕi) is introduced in present work:
(24)δ(ϕi)={0ϕi≥0.99,i≠11othercases.
Based on the above function, we can describe the effective chemical reaction rate in the diffuse domain of species i as
(25)rieff=ϕ1rihom+(1-ϕ1)δ(ϕi)rihet,
with rihom as the homogeneous reaction rate (gaseous reaction rate) of species i and rihet as the heterogeneous reaction rate (surface reaction rate) of species i.


## 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:
(26)∇·(ϕ1u)=0,(27)ρ∂(ϕ1u)∂t=-ρ(ϕ1u·∇)u+ϕ1∇·(-p)I+η∇·[∇(ϕ1u)+(∇(ϕ1u))T]+Bη(ϕ12(1-ϕ1)u),
here u (m·s-1) is the velocity. ρ (kg·m-3) is the density (the same density assumed for different gaseous species for simplicity), and η (Pa·s) is the viscosity. I is the unit metric and T is the inversion sign. B 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
(28)∇·(ϕ1u)=0,ρ∂(ϕ1u)∂t=F1(ϕ,u,p)=-ρ(ϕ1u·∇)u+ϕ1∇·(-p)I+η∇·[∇(ϕ1u)+(∇(ϕ1u))T]+Bη(ϕ12(1-ϕ1)u),∂ϕi∂t=F2(ϕ,c)=∑j≠in-2sijMijn×[∑j≠in[εij22(∇2ϕj-∇2ϕi)+Wij(ϕj-ϕi)]+ΔGi-ΔGj{∑j≠in}],∂tci=F3(ϕ,c)=-ϕ1(u·∇)ci+∇·[D∇ci]+rieff,i∈{gaseousspecies},
where vector ϕ=(ϕ1,…,ϕN) and c=(c1,…,cM).


## 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:(29)∇·(ϕ1nun+1)=0,ρϕ1nun+1-unΔt=F1(ϕn,un+1,pn+1),ϕin+1-ϕinΔt=F2(ϕn+1,cn),cin+1-cinΔt=F3(ϕn,cn+1).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 ϕ~i∈Vm, c~l∈Vm, p~∈Vm1, and u~∈Vmd:
(30)∫Ωp~∇·ϕ1udv=0∫Ωρu~·∂t(ϕ1u)dv+∫Ωu~:ρϕ1u·∇udv+∫Ω(-ϕ1pI+η(∇(ϕ1u)+∇(ϕ1u)T)):∇u~dv-∫Γ(-ϕ1pI+η(∇(ϕ1u)+∇(ϕ1u)T)):u~·nda-∫ΩBηϕ12(1-ϕ1)u·u~dv=0,(31)∫Ωϕ~i∂tϕidv-∫Ω∑j≠in2sijMijn[∑j≠in(εij22∇ϕi)∇ϕ~i]dv+∫Γ∑j≠in2sijMijn[∑j≠in(εij22∇ϕi)ϕ~i·n]da+∫Ω∑j≠in2sijMijn×[∑j≠in(εij22∇2ϕj+Wij(ϕj-ϕi))+ΔGi-ΔGj{∑j≠in(εij22∇2ϕj+Wij(ϕj-ϕi))}]ϕ~idv,(32)∫Ω∂tcic~idv+∫Ω(ϕ1(u·∇)ci-rieff)c~idv-∫ΩD∇ci·∇c~idv+∫ΓD∇cic~i·nda=0,
with n as the outward unit normal.Here the spatial discretization of (32) is proposed as an example. Defining the approximation ci=NTci, where NT=(N1,…,Nn) is vector composed by shape functions, ci=(ci1,…,cin), n is the node number in element. Choosing shape functions as test functions and assuming the boundary conditions ∂ci/∂n=q, then we have
(33)K1∂tci+K2ci=fK1=∫ΩeNNTdvK2=∫Ωe[D(∂N∂x∂NT∂x+∂N∂y∂NT∂y)+(ϕ1(u∂N∂x+v∂N∂y)-rieff)NT]dxdyf=-∫ΓDNqdΓ.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)Ax=f.Notice that the coefficient matrixA 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 Figure3. 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 ϕ1=0.5.Figure 3
Schematic illustrations of unidirectional fiber-reinforced substrates.For the multiphase field order parametersϕi, there is no flux at the inlet, outlet, and wall, which can be described by the Neumann boundaries:
(35)n·∇ϕi=0on∂Ωinlet,∂Ωoutlet,∂Ωwall.For the species variableci, 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:
(36)cMTS=cMTS0,(uci-D∇ci)·n=0on∂Ωinlet(uci-D∇ci)·n=0on∂Ωwall(D∇ci)·n=0on∂Ωoutlet.HerecMTS0 is the molar concentration of MTS in the imported gas. It can be calculated by the ideal gas law (pMTS is the partial pressure of MTS):
(37)cMTS0=pMTS/(RT).For the variables (u,p) in Navier-Stokes equations, the velocity u 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 p is fixed as a constant. It is assumed there are no viscous effects hence no boundary layer develops on the wall. Then we have
(38)u=u0on∂Ωinletu=0on∂Ωg/sp=p0on∂Ωoutletn·u=0,t·[pI-η(∇u+(∇u)T)]n=0on∂Ωwall,
where u0=(0,0,u0) and t is a tangential vector to the boundary.The initial conditions for variables are straightforward:ϕi|t=0=1 if ith phase exists; ci|t=0=0 for not MTS species; u|t=0=u0 and p|t=0=p0.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 Figure3. Therefore the corresponding boundary conditions are n·∇ϕi=0 and n·∇ci=0. Notice that it is not necessary to set up boundary condition for the variables in Navier-Stokes equations (u,p), because the actual boundary for these variables come to be the interface with ϕ1=0.5, namely, gas/solid boundary ∂Ωg/s. The parameters related to model are listed in Table 1.Table 1
Parameters related to model.
Parameter
Value
ϵ
i
j
(
Jm
-
1
)
1
/
2
5.8e-4
p (Pa)
1.01325e5
T (K)
1123
ρ
g
(Kg·m-3)
0.067
D
CH
3
SiCl
3
,
H
2
(
m
2
·
s
-
1
)
7.26e-4
γ
H
2
0.9
l
2
D
(
m
)
2e-5
η
s
(
Pa
·
s
)
2.6
W
i
j
(
Jm
-
1
)
1.32e7
M
i
j
(
ms
-
1
)
2.09e-7
u
x
(
ms
-
1
)
3e-2
ρ
s
(
Kg
·
m
-
3
)
3200.0
D
HCl
,
H
2
(
m
2
·
s
-
1
)
1.47e-3
γ
MTS
0.1
l
3
D
(
m
)
1e-5
η
g
(
Pa
·
s
)
2.62e-5


### 3.2. 2D and 3D Simulation

The numerical simulation is performed based on the software “COMSOL Multiphysics” with the relative tolerance1e-2 and absolute tolerance 1e-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).Figure4 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 H2 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 cMTS along both these two sections decreases; the gradient decreasing of cMTS along section A is obvious which is mainly caused by chemical reaction as being explained later; there barely exists change of cMTS 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)t=0.1; (b) t=0.5; (c) t=1; (d) t=1.5.Figure 5
Time dependentcMTS at section A and B.By changing the reaction diffusion equation, the effect of fluid flow and chemical reaction oncMTS 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 oncMTS.Under appreciated infiltration conditions, the single phase SiC is deposited. Figure7(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 z-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.SiC deposited as grains.
(a)
SiC deposited as grains
(b)
SiC grain in 3DFigure 8
Blockage formation during SiC fiber deposited with SiC.


## 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 Figure3. 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 ϕ1=0.5.Figure 3
Schematic illustrations of unidirectional fiber-reinforced substrates.For the multiphase field order parametersϕi, there is no flux at the inlet, outlet, and wall, which can be described by the Neumann boundaries:
(35)n·∇ϕi=0on∂Ωinlet,∂Ωoutlet,∂Ωwall.For the species variableci, 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:
(36)cMTS=cMTS0,(uci-D∇ci)·n=0on∂Ωinlet(uci-D∇ci)·n=0on∂Ωwall(D∇ci)·n=0on∂Ωoutlet.HerecMTS0 is the molar concentration of MTS in the imported gas. It can be calculated by the ideal gas law (pMTS is the partial pressure of MTS):
(37)cMTS0=pMTS/(RT).For the variables (u,p) in Navier-Stokes equations, the velocity u 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 p is fixed as a constant. It is assumed there are no viscous effects hence no boundary layer develops on the wall. Then we have
(38)u=u0on∂Ωinletu=0on∂Ωg/sp=p0on∂Ωoutletn·u=0,t·[pI-η(∇u+(∇u)T)]n=0on∂Ωwall,
where u0=(0,0,u0) and t is a tangential vector to the boundary.The initial conditions for variables are straightforward:ϕi|t=0=1 if ith phase exists; ci|t=0=0 for not MTS species; u|t=0=u0 and p|t=0=p0.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 Figure3. Therefore the corresponding boundary conditions are n·∇ϕi=0 and n·∇ci=0. Notice that it is not necessary to set up boundary condition for the variables in Navier-Stokes equations (u,p), because the actual boundary for these variables come to be the interface with ϕ1=0.5, namely, gas/solid boundary ∂Ωg/s. The parameters related to model are listed in Table 1.Table 1
Parameters related to model.
Parameter
Value
ϵ
i
j
(
Jm
-
1
)
1
/
2
5.8e-4
p (Pa)
1.01325e5
T (K)
1123
ρ
g
(Kg·m-3)
0.067
D
CH
3
SiCl
3
,
H
2
(
m
2
·
s
-
1
)
7.26e-4
γ
H
2
0.9
l
2
D
(
m
)
2e-5
η
s
(
Pa
·
s
)
2.6
W
i
j
(
Jm
-
1
)
1.32e7
M
i
j
(
ms
-
1
)
2.09e-7
u
x
(
ms
-
1
)
3e-2
ρ
s
(
Kg
·
m
-
3
)
3200.0
D
HCl
,
H
2
(
m
2
·
s
-
1
)
1.47e-3
γ
MTS
0.1
l
3
D
(
m
)
1e-5
η
g
(
Pa
·
s
)
2.62e-5


## 3.2. 2D and 3D Simulation

The numerical simulation is performed based on the software “COMSOL Multiphysics” with the relative tolerance1e-2 and absolute tolerance 1e-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).Figure4 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 H2 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 cMTS along both these two sections decreases; the gradient decreasing of cMTS along section A is obvious which is mainly caused by chemical reaction as being explained later; there barely exists change of cMTS 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)t=0.1; (b) t=0.5; (c) t=1; (d) t=1.5.Figure 5
Time dependentcMTS at section A and B.By changing the reaction diffusion equation, the effect of fluid flow and chemical reaction oncMTS 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 oncMTS.Under appreciated infiltration conditions, the single phase SiC is deposited. Figure7(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 z-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.SiC deposited as grains.
(a)
SiC deposited as grains
(b)
SiC grain in 3DFigure 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 systemAx=f 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.



---

*Source: 836187-2012-12-23.xml*