Modeling of Thermal Conductivity of CVI-Densified Composites at Fiber and Bundle Level

The evolution of the thermal conductivities of the unidirectional, 2D woven and 3D braided composites during the CVI (chemical vapor infiltration) process have been numerically studied by the finite element method. The results show that the dual-scale pores play an important role in the thermal conduction of the CVI-densified composites. According to our results, two thermal conductivity models applicable for CVI process have been developed. The sensitivity analysis demonstrates the parameter with the most influence on the CVI-densified composites’ thermal conductivity is matrix cracking’s density, followed by volume fraction of the bundle and thermal conductance of the matrix cracks, finally by micro-porosity inside the bundles and macro-porosity between the bundles. The obtained results are well consistent with the reported data, thus our models could be useful for designing the processing and performance of the CVI-densified composites.


Introduction
The chemical vapor infiltration (CVI) process has been considered as one of the most important fabrication methods for advanced thermostructural materials [1][2][3][4][5]. Compared with the conventional processing techniques, CVI techniques have distinct advantages, such as minimizing damage to fibers, imposing little mechanical stress on fibers, applicable for a wide variety of ceramic matrix composite materials, producing near-net-shape composites. Despite these advantages, they also have some drawbacks, such as the long processing time, premature closing of the pores at the surface. For fabricating the thick composites by less time, the forced flow (F-CVI) and thermal gradient chemical vapor infiltration (TG-CVI) process has been developed as alternatives to isothermal-isobaric chemical vapor infiltration (ICVI) [1,6]. In TG-CVI, a thermal gradient is imposed on the preform to prevent premature sealing of the gas entrance. Thus, the detailed knowledge of the evolution of the thermal transport properties of the preform during the infiltration is essential for understanding and optimizing TG-CVI.
For now, relatively few numerical studies [7][8][9][10] have been carried out on this problem. For instance, considering the fiber and matrix as the same material and pores as non-conductive, Tomadakis and Sotirchos [7,11] used Monte Carlo algorithm to calculate the thermal transport properties of 1-D partially overlapping and randomly overlapping fiber structure. In contrast, Skamser et al. [8] considered the contribution of void phase to the thermal transport properties and provided more realistic data for 1-D SiC f /SiC composite obtained by their simulated microstructure. Vignoles et al. [9] firstly computed the thermal conductivity for a square array of fibers coated with anisotropic pyrocarbon matrix, then estimated the thermal conductivity of stitched cloth layup composites according to the microstructure geometric information gained through X-ray computerized microtomography (CMT). To the best of our knowledge, very few numerical studies have been reported for thermal transport properties of 2D woven and 3D braided composites in the bundle scale for CVI process.
Most CVI-densified composites possess high temperature capability and have been used in a number of demanding applications in high temperature environments. Hence a knowledge of the thermal transport properties of such composites is also of considerable importance and has been the subject of extensive experimental [12][13][14][15][16][17] and theoretical [18][19][20] studies. However, although a lot of thermal conductivity models in the fiber and bundle scales have been developed, the one appropriate ones for CVI-densified composites are very few, due to the existence of the dual-scale pores inside the composites. In addition, CVI-densified composites may be fabricated depending on the fiber type selected and the chosen densification route. Hence it is important, not only to generate experimental data, but also to characterize the thermal transport behavior in terms of constituent properties across the wide range of porosity.
For the previous reasons, the present work is the continuation of our previous study [21,22] and aimed at developing an efficient thermal conductivity model for CVI-densified composites. To accomplish this, the unidirectional fiber structures, 2D woven and 3D braided bundle structures are taken as research subjects. The effect of the individual component (fiber, interfacial layer and pore) on the transverse thermal conductivity of unidirectional fiber structures is first numerically investigated. Secondly, we numerically simulate the evolution of the macro-pore structures for 2D and 3D composites and investigate the effect of the macro-pore on the thermal conductivity of the simulated structures. Finally, two analytical models for these structures are developed based on our numerical data and compared to the available experimental data.

Description of the Model
For unidirectional fiber structures, the thermal conductivity parallel to the fiber direction is considered to be well described by the rule of mixtures, while the situation in the transverse direction is more complicated. Thus, the purpose of this section is to develop an analytical model to relate the transverse thermal conductivity of these structures to the ones of their four components (fiber, interfacial layer, matrix and pore).
To investigate the fiber distribution effects on transverse thermal conductivity of fiber-reinforced composites, the microstructure with random fiber distribution can be generated using a program developed in MATLAB. Specifically, the procedure starts with a specified microstructure with square cross-section and proceeds by randomly and sequentially adding the circular fiber inside the microstructure. Each added fiber is accepted and its position is recorded only if it does not overlap with any of the ones previously accepted. For high fiber volume fraction, a small amount of overlapping is allowable. In addition, the boundaries of the microstructure act like rigid walls; thus, no fibers pass through the boundaries, and all fibers are restricted inside the microstructure. The microstructures generated in this manner are primarily governed by the choice of porosity p and the allowable inter-fiber distance d min . The former is calculated as p = 1 − N f πr f 2 /L 0 2 , where N f is the number of included fibers, r f is the fiber radius, and L 0 is the length of the microstructure. The latter (d min ) is defined such that the shortest allowable distance between the centers of two fibers is larger than (2 + d min )r f . The selection values of d min are 0.1 or 0 in this study. For the case of heat transfer transverse to the axis of randomly distributed fibers, the representative volume element (RVE) may be simulated by a square window placed on the cross section of the microstructure with the following relative scale dimensions: L 0 = 1.5L, where L is the length of the RVE (as depicted in Figure 1). In addition, to assure the RVE contains a statistical representation of all local heterogeneities, the heterogeneities, the relative RVE size with respect to the scale of microstructural dimensions should be sufficiently large. The study of the scale effects indicates that the estimated thermal conductivities (as averaged over three independent simulation runs) are almost constant when the REV size is larger than 50rf, thus the sizes of RVE are set as 50rf × 50rf in this study. To characterize the isotropy of the simulated fiber structures, we consider there are N fibers located less than distance R apart from the center of the microstructure, where R = 50rf in this study, and introduce the anisotropy parameter (  ): where i  is the orientation angle (with respect to the positive x coordinate) of the line connecting the center of the unit cell and the ith fiber. The values of i  is +1 or −1 when all fibers are distributed along the x-axis or along the y-axis. When fibers structure is isotropic,  approaches zero. The absolute values of  were found typically less than 0.01 for all simulated fiber structures, which we consider to satisfy the isotropy requirement.
The resultant microstructures are shown in Figure 2a,b. It can be observed that the generated fiber structures have different degree of local non-uniformity, although all of them are statistically homogenous. Thus, prior to any analysis, it is necessary to quantify the microstructure of the fiber distribution. In this effort, the normalized average nearest inter-fiber distance, l, whose statistics differ between different fiber distributions, is introduced to characterize the degree of the non-uniformity. l is defined as where Li is the nearest distance between the ith fiber and its neighbor, ri is the radius of the ith fiber and ri nearest is the radius of the fiber which is nearest to the ith fiber. Given the same fiber volume, To characterize the isotropy of the simulated fiber structures, we consider there are N fibers located less than distance R apart from the center of the microstructure, where R = 50r f in this study, and introduce the anisotropy parameter (φ): where θ i is the orientation angle (with respect to the positive x coordinate) of the line connecting the center of the unit cell and the ith fiber. The values of θ i is +1 or −1 when all fibers are distributed along the x-axis or along the y-axis. When fibers structure is isotropic, φ approaches zero. The absolute values of φ were found typically less than 0.01 for all simulated fiber structures, which we consider to satisfy the isotropy requirement. The resultant microstructures are shown in Figure 2a,b. It can be observed that the generated fiber structures have different degree of local non-uniformity, although all of them are statistically homogenous. Thus, prior to any analysis, it is necessary to quantify the microstructure of the fiber distribution. In this effort, the normalized average nearest inter-fiber distance, l, whose statistics differ between different fiber distributions, is introduced to characterize the degree of the non-uniformity. l is defined as where L i is the nearest distance between the ith fiber and its neighbor, r i is the radius of the ith fiber and r i nearest is the radius of the fiber which is nearest to the ith fiber. Given the same fiber volume, small values of l are associated with heterogeneous patterns, while large values of l indicate homogeneous patterns, and its variation can be seen as a measure of the fiber dispersion. In Figure 3, the data for two fiber structures are compared with Spitzig's experimental measurement [23]. We can observe that, Spitzig's experimental data is consistent with our data for the simulated fiber structure with d min = 0, hence only the thermal conductivities of the simulated fiber structure with d min = 0 are chosen as the research subject. By adopting the locations of the fibers and increasing the radius of the circles, the material systems with different phase have been developed and shown in Figure 2b-e. Appling the corresponding temperature and thermal insulation boundaries to the unit cells, the transverse thermal conductivity of these material systems can be calculated using finite element software "COMSOL Multiphysics" coupled with MATLAB. Then the effect of all constituent (fiber, interfacial layer, matrix and pore) thermal conductivities on the global thermal transport properties of the material systems can be studied. The full solving details can be found in Reference [24]. small values of l are associated with heterogeneous patterns, while large values of l indicate homogeneous patterns, and its variation can be seen as a measure of the fiber dispersion. In Figure 3, the data for two fiber structures are compared with Spitzig's experimental measurement [23]. We can observe that, Spitzig's experimental data is consistent with our data for the simulated fiber structure with dmin = 0, hence only the thermal conductivities of the simulated fiber structure with dmin = 0 are chosen as the research subject. By adopting the locations of the fibers and increasing the radius of the circles, the material systems with different phase have been developed and shown in Figure 2b-e. Appling the corresponding temperature and thermal insulation boundaries to the unit cells, the transverse thermal conductivity of these material systems can be calculated using finite element software "COMSOL Multiphysics" coupled with MATLAB. Then the effect of all constituent (fiber, interfacial layer, matrix and pore) thermal conductivities on the global thermal transport properties of the material systems can be studied. The full solving details can be found in Reference [24].

The Analytical Model
In 1892, Rayleigh [25] investigated the conductivity of an uniform medium composed of cylindrical obstacles (fibers) arranged in rectangular order, applied two boundary conditions (radial heat flux and temperature potential continuity) at the fiber/matrix interface and solved the steady-state heat conduction (mathematically the potential equation) to give approximate function to predict transverse ETC [26,27]: where B is a dimensionless variable defined as Equation (5)

The Analytical Model
In 1892, Rayleigh [25] investigated the conductivity of an uniform medium composed of cylindrical obstacles (fibers) arranged in rectangular order, applied two boundary conditions (radial heat flux and temperature potential continuity) at the fiber/matrix interface and solved the steady-state heat conduction (mathematically the potential equation) to give approximate function to predict transverse ETC [26,27]: where B is a dimensionless variable defined as where V f , K T e , K m , K T f and k are the volume fraction of fiber, the transverse thermal conductivity of the one-dimensional composite, the thermal conductivity of the isotropic matrix, the transverse thermal conductivity of isotropic fiber and the fiber-to-matrix conductivity ratio, respectively. When the fiber volume fraction is low, the higher order terms of Equation (3) can be ignored, then Rayleigh's model reduces to the so called self-consistent formula [28]: Another famous model is developed by Hasselman and Johnson [29], who considered the imperfect interface between fibers and matrix and introduced an effective interfacial conductance H i , their model is expressed as: where r f is the fiber radius. If we introduced an interfacial Biot number, B i (B i = H i r f /K T f ) and Equation (5) into Hasselman and Johnson's model, then H-J model can be easily rearranged into the following form: Comparing Equation (6) with (8), it is clearly found that the self-consistent model is just the special case of H-J model when 1/B i →0, the composites with imperfect fiber/matrix interface can be practically seen as the ones with perfect conductance, in which the effect of imperfect interface is lumped into fiber-to-matrix conductivity ratio k. Typical to the cylindrically orthotropic carbon fibers-reinforced composites, Later, Hasselman [30] used the same method to develop the following formula to estimate the transverse thermal conductivity: where K θ and K r are the radial and tangential thermal conductivity of the cylindrically orthotropic fiber, respectively. Obviously, Equation (9) agrees to Equation (7) if K θ = K r = K T f . However, above models are developed on the assumption that the heat flux through the single fiber is not affected by the others, thus these developed formulae are valid only for dilute fiber concentrations. At higher volume-fractions of fiber, these formulae are approximate or even erroneous, for the local temperature fields around the fibers will interact. The differences between the calculated transverse thermal conductivity of two-phase model (as shown in Figure 2b) and the prediction of self-consistent formula with the fiber-to-matrix conductivity ratio are shown in Figure 4. The relative differences (RD) between both are illustrated by the following function: where the subscripts MOD and FEM denote the K e values are calculated by the model and the simulation, respectively. It is found that the self-consistent formula can give a satisfactory prediction for all fiber structures when 0.01 < K T f K m < 100. When K T f K m < 0.01 or K T f K m > 100, the following fitting function can be used: where B has been defined by Equation (4). In Figure 4, Equation (11) is shown to give a satisfactory agreement to our data. In most cases, the fiber-matrix ratio is at the range of 10 −2 -10 2 , thus the self-consistent formula (or Hasselman and Johnson's model) is satisfactory for estimation. Comparing Equation (6) with (8), it is clearly found that the self-consistent model is just the special case of H-J model when 1/Bi→0, the composites with imperfect fiber/matrix interface can be practically seen as the ones with perfect conductance, in which the effect of imperfect interface is lumped into fiber-to-matrix conductivity ratio k. Typical to the cylindrically orthotropic carbon fibers-reinforced composites, Later, Hasselman [30] used the same method to develop the following formula to estimate the transverse thermal conductivity: where Kθ and Kr are the radial and tangential thermal conductivity of the cylindrically orthotropic fiber, respectively. Obviously, Equation (9) agrees to Equation (7) if However, above models are developed on the assumption that the heat flux through the single fiber is not affected by the others, thus these developed formulae are valid only for dilute fiber concentrations. At higher volume-fractions of fiber, these formulae are approximate or even erroneous, for the local temperature fields around the fibers will interact. The differences between the calculated transverse thermal conductivity of two-phase model (as shown in Figure 2b) and the prediction of self-consistent formula with the fiber-to-matrix conductivity ratio are shown in Figure 4. The relative differences (RD) between both are illustrated by the following function: where the subscripts MOD and FEM denote the K e values are calculated by the model and the simulation, respectively. It is found that the self-consistent formula can give a satisfactory prediction for all fiber structures when 0.01 < following fitting function can be used: where B has been defined by Equation (4). In Figure 4, Equation (11) is shown to give a satisfactory agreement to our data. In most cases, the fiber-matrix ratio is at the range of 10 −2 -10 2 , thus the self-consistent formula (or Hasselman and Johnson's model) is satisfactory for estimation.  (6)) and the fitted (Equation (11)) equations.  (6)) and the fitted (Equation (11)) equations.
To improve the interface between matrix and fiber, the fibers are coated by the pyrolytic carbon (PyC) interfacial layer before the deposition of the matrix. In self-consistent (Equation (6)) and H-J model (Equation (7)), a fiber is in the middle of a cylindrical matrix and the matrix is contained within an effectively homogeneous medium, which is very close to the situation of CVI processing if the thickness of the interfacial layer is not too large. Then the coated fiber with interface layer can be approximately treated as a homogeneous solid phase, whose volume fraction is (V f + V i ), where V i is the volume fraction of the interfacial layer. In addition, the equivalent transverse thermal conductivity of the solid phase, K T f i can be estimated according to H-J model: where K i is the thermal conductivity of the interfacial layer. Substituting Equation (12) into (6), the transverse thermal conductivity of the three-phase system can be estimated. To validate this model, by setting K T f = 3 W·m −1 ·K −1 , K m = 20 W·m −1 ·K −1 , V f = 0.5 and V i = 0.05, the transverse thermal conductivity of the three-phase system (as shown in Figure 2c) with K i values are calculated and shown in Figure 5. It is clearly shown that the presented model gives a good fit to our simulation results. To improve the interface between matrix and fiber, the fibers are coated by the pyrolytic carbon (PyC) interfacial layer before the deposition of the matrix. In self-consistent (Equation (6)) and H-J model (Equation (7)), a fiber is in the middle of a cylindrical matrix and the matrix is contained within an effectively homogeneous medium, which is very close to the situation of CVI processing if the thickness of the interfacial layer is not too large. Then the coated fiber with interface layer can be approximately treated as a homogeneous solid phase, whose volume fraction is (Vf + Vi), where Vi is the volume fraction of the interfacial layer. In addition, the equivalent transverse thermal conductivity of the solid phase, T fi K can be estimated according to H-J model: where Ki is the thermal conductivity of the interfacial layer. Substituting Equation (12) into (6), the transverse thermal conductivity of the three-phase system can be estimated. To validate this model, by setting and Vi = 0.05, the transverse thermal conductivity of the three-phase system (as shown in Figure 2c) with Ki values are calculated and shown in Figure 5. It is clearly shown that the presented model gives a good fit to our simulation results. Considering that the fiber, matrix and interfacial layer have the same thermal conductivity values, Ks, above a certain pore concentration, the transverse thermal conductivities of the systems are the same order of magnitude as the thermal conductivity of the pore, which indicates that the pore phase is the dominant pathway for heat transfer. As the solid volume fraction is increased and larger than the 'critical solid concentration', c s V , the matrix deposited on the surfaces of the different fibers contact with each other and form a continuous path for heat transfer. Noted a parameter called the critical pore concentration has been frequently used to estimate the gas transport properties of the densified microstructure for CVI process, then the value of c s V can be determined by subtracting the critical pore concentration for 1-D partially overlapping fiber structure from one, the following expression can be given from our simulation results: Considering that the fiber, matrix and interfacial layer have the same thermal conductivity values, K s , above a certain pore concentration, the transverse thermal conductivities of the systems are the same order of magnitude as the thermal conductivity of the pore, which indicates that the pore phase is the dominant pathway for heat transfer. As the solid volume fraction is increased and larger than the 'critical solid concentration', V c s , the matrix deposited on the surfaces of the different fibers contact with each other and form a continuous path for heat transfer. Noted a parameter called the critical pore concentration has been frequently used to estimate the gas transport properties of the densified microstructure for CVI process, then the value of V c s can be determined by subtracting the critical pore concentration for 1-D partially overlapping fiber structure from one, the following expression can be given from our simulation results: A phase interchange theorem is developed by Keller [31] for two-dimensional statistically homogeneous, isotropic structures, which simply states that, the K T e /K m value of this structure by interchanging the thermal conductivity values of fiber and matrix is equal to the reciprocal of the K T e /K m value of the original structure. When V s is less than V c s , than Tomadakis and Sotirchos [32] applied Keller's theorem and related the transverse thermal transport of 1-D fiber structure when the fiber are high-conductive and the gas transport of 1-D overlapping fiber structure when the fiber are non-conductive: where is the reduce transverse thermal conductivity of the systems when the solid thermal conductivity is highly larger then the one of pore, τ is the Fick diffusion tortuosity of the systems [22]. Due to the local non-uniformity of the fiber structure, the coated fibers will overlap each other and form the local conductive chains. If the local conductive chains can be approximated as the fibers with elliptic cross-section, then the effective thermal conductivity of this material system can be estimated by Tsai-Halpin equations [33,34]: for elliptic fibers with aspect ratios e, ξ = √ 3 ln(e). Then the reduced transverse thermal conductivity of the systems can be obtained by Tsai-Halpin equations when the solid is high-conductive: Combining Equations (14) and (16), we can obtain Finally, the transverse thermal conductivity of the systems can be estimated by Equations (15) and (17) when V s is less than V c s . When V s is larger than V c s , the pore can be seen as the isolative phase inside the solid media. Because the thermal conductivity of the pores is always less than the one of the solid by two orders, the pores can be practically seen as non-conductive. Fitting from the simulation data, 3 expression is made to put forth a porosity correction to thermal conductivity: where β is a pore shape factor. The obtained β values are plotted in Figure 6. It is found that the fluctuation of β values is approximately increased with porosity throughout the entire porosity range. The values of β are estimated as β = 70V p with R 2 ≈ 0.88. Such high level of dispersion demonstrates that the thermal conductivity of the two-phase system could considerably be affected by the randomness of the geometry of the pores. According to the previous discussions, the complete analytical model for four-phase system (as depicted in Figure 2f) is developed as follow. The fiber coated with two layers of materials can be seen as the equivalent homogeneous solid, whose thermal conductivity is estimated by using Hasselman and Johnson's model twice, first to replace the thermal conductivity of fiber coated with interface layer with the homogenous value, and second to replace the thermal conductivity of fiber coated with two layers of materials with the homogenous value, then the four-phase systems reduced to the two-phase systems. The complete analytical expressions are listed in Table 1.  According to the previous discussions, the complete analytical model for four-phase system (as depicted in Figure 2f) is developed as follow. The fiber coated with two layers of materials can be seen as the equivalent homogeneous solid, whose thermal conductivity is estimated by using Hasselman and Johnson's model twice, first to replace the thermal conductivity of fiber coated with interface layer with the homogenous value, and second to replace the thermal conductivity of fiber coated with two layers of materials with the homogenous value, then the four-phase systems reduced to the two-phase systems. The complete analytical expressions are listed in Table 1.

Expression
When To validate the presented model, by assuming the perfect interface conductance and setting K T f = 1 W·m −1 ·K −1 , K i = 40 W·m −1 ·K −1 , K m = 10 and 100 W·m −1 ·K −1 , K p = 0.01 W m −1 ·K −1 , V f = 0.5 and V i = 0.05, the transverse thermal conductivity of the densified microstructure with the solid volume fraction during the CVI process are calculated and compared with the presented and Youngblood's expressions [35] in Figure 7. In Youngblood's model, the thermal conductivity of matrix containing pores can be estimated as the one of homogeneous media, K m by the Maxwell-Eucken expression: Then the effective thermal conductivity of a composite with a well-bonded fiber coating can be estimated by "3-cylinder" Markworth's model [36]: where C and D are given by where r int is the thickness of the interface layer. It is clearly shown that the present model is consistent with the simulation results throughout the entire range of solid fraction, while Youngblood's model agrees well with the simulation results by setting β value as 10 when porosity is less than 20%.

The Effect of Matrix Cracking on the Global Thermal Conductivity
Due to the matrix cracking and interfacial de-bonding caused by the coefficient of thermal expansion (CTE) mismatch [37] after the processing, the effective thermal conductivity values predicted by the model will be larger than the actual values in practice. In the transverse direction, the matrix cracking plays little effect on the thermal conductivity if the cracking density is not very high. Thus, the H-J model which takes the interfacial de-bonding into account is sufficient, and the effective interfacial conductance in the model is considered as an average over all fibers and a distribution of de-bonding [16]. In the longitudinal direction, Lu and Hutchinson [38] pointed out

The Effect of Matrix Cracking on the Global Thermal Conductivity
Due to the matrix cracking and interfacial de-bonding caused by the coefficient of thermal expansion (CTE) mismatch [37] after the processing, the effective thermal conductivity values predicted by the model will be larger than the actual values in practice. In the transverse direction, the matrix cracking plays little effect on the thermal conductivity if the cracking density is not very high. Thus, the H-J model which takes the interfacial de-bonding into account is sufficient, and the effective interfacial conductance in the model is considered as an average over all fibers and a distribution of de-bonding [16]. In the longitudinal direction, Lu and Hutchinson [38] pointed out that the reduction of thermal conductivity due to matrix cracking can be as large as 50%-70% if the fiber conductivity is comparable to that of matrix and the fiber volume fraction is at the range of 0.3-0.5, and derived a model with matrix cracking density based on shear-lag analysis: where K L 0 is the initial effective conductivity, K L f is the longitudinal thermal conductivity of the fiber, l c is the spacing between matrix cracks, and H c is the thermal conductance of the matrix cracks. The estimated microstructural quantities used in the model are listed in Table 2.

Description of the Model
Recently, we used the nonconservative level set finite-element method to simulate the evolution of the boundaries of the solid phase for the complex bundle structures and calculate the gas transport properties for CVI process [21,22]. In this study, the same method was used to simulate the macro-pore structure of 2D woven and 3D braided bundle structures. The details of this method can be found in Reference [41] and are summarized briefly here. In this method the interface between solid and void phases is represented by a smooth step function of space and time, called level set function φ. φ equals zero in a domain which represents the void phase and one in the other which represents the solid phase. Across the interface, there is a smooth transition from zero to one, and it is defined by the 0.5 isocontour.
The level set method is first used to simulate the growth of the boundaries of solid phase. Then the conduction within the simulated structures is governed by the following form of Laplace's Equation: where T denotes temperature, H a material index which identifies the sold and void phases, and is defined as the function of φ: Appling the corresponding temperature and thermal insulating boundaries, the effective thermal conductivity K e can be calculated [24,42].
The simulated structure for 2D woven composites is shown in Figure 8a. In the constructed unit cell, we assume 2D woven composites are composed of the fiber tows with elliptical section, the tows follow a sinusoidal path and ten laminates are randomly stacked in 0 • /90 • sequence. By fixing the ratio of the gap width to the tow width between two adjacent tows in the longitudinal direction as 0.20 and varying the ratio of tow width to tow thickness, e from 5 to 10, the initial macro-porosities of the simulated structures are 0.28-0.32, which is consistent with the experimental data [43]. Thus, the simulated structures adopting these geometric parameters are considered to be appropriate for the estimation of the thermal transport properties. For convenience of illustration, the longitudinal and through-thickness directions are considered as the x and z directions, respectively. The simulated structure for 3D braided composites is shown in Figure 8b. In the constructed unit cell, Chen and Choy's analytical model [44] is adopted and the braiding angle is setting as 20 • . For convenience of illustration, the transverse and longitudinal directions are considered as the x and z directions, respectively.
ppling the corresponding temperature and thermal insulating boundaries, the effective l conductivity Ke can be calculated [24,42]. e simulated structure for 2D woven composites is shown in Figure 8a. In the constructed unit e assume 2D woven composites are composed of the fiber tows with elliptical section, the llow a sinusoidal path and ten laminates are randomly stacked in 0°/90° sequence. By fixing io of the gap width to the tow width between two adjacent tows in the longitudinal direction and varying the ratio of tow width to tow thickness, e from 5 to 10, the initial macro-porosities simulated structures are 0.28-0.32, which is consistent with the experimental data [43]. Thus, ulated structures adopting these geometric parameters are considered to be appropriate for imation of the thermal transport properties. For convenience of illustration, the longitudinal rough-thickness directions are considered as the x and z directions, respectively. The ted structure for 3D braided composites is shown in Figure 8b. In the constructed unit cell, nd Choy's analytical model [44] is adopted and the braiding angle is setting as 20°. For ience of illustration, the transverse and longitudinal directions are considered as the x and z ons, respectively.

The Thermal Conductivity Model for Bundle Structures
Assuming bundle and matrix are isotropic and the thermal conductivity of both are eq and considering the thermal conductivity of the solid is much larger than that of void, the conductivity of the composites can be expressed by Maxwell-Eucken expression: where β is a pore shape factor. The obtained β values by our simulation are plotted in Figu the estimated  values for different macro-pore structures are summarized in Table 3. It that β values generally increase with the average angle of the bundle orientation with resp global heat transfer direction, in other words, the reduction of the thermal conductivity with the angle of the bundle orientation with respect to the global heat transfer direction. by the different geometries of the macro-pores, β values may vary by about 60% or less. interesting finding is that β values of the random stacking are not between the ones of the sy stacking and the ones of the opposition stacking. Here the symmetry stacking means all exactly one over each other, while the opposition stacking means the layer is shifted a in-plane direction to make the raised portion of the layer insert into the hole of the adjacent foundation means the thermal conductivity of the random stacking is not between the on symmetry stacking and the one of the opposition stacking at the same macro-porosity, opposite of the acknowledgement for gas transport of 2D woven composites [45].

The Thermal Conductivity Model for Bundle Structures
Assuming bundle and matrix are isotropic and the thermal conductivity of both are equal to K s , and considering the thermal conductivity of the solid is much larger than that of void, the thermal conductivity of the composites can be expressed by Maxwell-Eucken expression: where β is a pore shape factor. The obtained β values by our simulation are plotted in Figure 9 and the estimated β values for different macro-pore structures are summarized in Table 3. It is found that β values generally increase with the average angle of the bundle orientation with respect to the global heat transfer direction, in other words, the reduction of the thermal conductivity increase with the angle of the bundle orientation with respect to the global heat transfer direction. Affected by the different geometries of the macro-pores, β values may vary by about 60% or less. Another interesting finding is that β values of the random stacking are not between the ones of the symmetry stacking and the ones of the opposition stacking. Here the symmetry stacking means all layers lie exactly one over each other, while the opposition stacking means the layer is shifted along the in-plane direction to make the raised portion of the layer insert into the hole of the adjacent one. This foundation means the thermal conductivity of the random stacking is not between the one of the symmetry stacking and the one of the opposition stacking at the same macro-porosity, which is opposite of the acknowledgement for gas transport of 2D woven composites [45].  In Tables 4 and 5, we test the sensitivity of the model to β values. The results show that the dispersion of β values has a small effect, order of 10% or less, on the thermal conductivities of 2D woven composite along the in-plane direction and the ones of 3D braided composites in the longitudinal direction, while the thermal conductivities of 2D woven composite along the through-thickness direction and the ones of 3D braided composites in the transverse direction vary by about 30% and 15% with β values, respectively. The residual macro-porosity is about 10% for the composites densified by CVI, then the uncertainty in the thermal conductivities in the transverse direction due to the macro-pore geometry is on the order of ~10%. Table 4. Sensitivity analysis of the pore shape factor for 2D woven and 3D braided composites.  In Tables 4 and 5, we test the sensitivity of the model to β values. The results show that the dispersion of β values has a small effect, order of 10% or less, on the thermal conductivities of 2D woven composite along the in-plane direction and the ones of 3D braided composites in the longitudinal direction, while the thermal conductivities of 2D woven composite along the through-thickness direction and the ones of 3D braided composites in the transverse direction vary by about 30% and 15% with β values, respectively. The residual macro-porosity is about 10% for the composites densified by CVI, then the uncertainty in the thermal conductivities in the transverse direction due to the macro-pore geometry is on the order of~10%. Table 4. Sensitivity analysis of the pore shape factor for 2D woven and 3D braided composites.  In the actual case, the thermal conductivity of bundle is not equal to that of matrix, the bundle and matrix may show different heat transfer capabilities in the transverse and axial directions. Hence our model should be revised to consider this heterogeneity. Considering the matrix are deposited on the surface of the bundles, then it generally follows the same direction of the bundle, hence the average thermal conductivity of the solid phase can also be approximated by the function related to the average angle of the bundle orientation with respect to the global heat transfer direction, θ:

2D, in-Plane
where K || b , K ⊥ b are the thermal conductivities of the bundle in the axial and transverse direction, respectively. For 2D woven composites along the in-plane direction, only a half of the total bundles is parallel to the heat transfer direction, while the other half is transverse to the heat transfer direction, thus θ = 45 • . For 2D woven composites along the through-thickness direction, the reported mean crimp angles are between 1 • and 5 • , thus we choose 2 • in this study and θ = 90 • − 2 • = 88 • . For 3D woven composites in the longitudinal direction, the bundles are laid along the transfer direction axis at a small braiding angle (~20 • ), and thus θ = 20 • . For 3D woven composites in the transverse direction, a half of the total bundles is laid along the heat transfer direction axis at an angle of (90 • − 20 • ), while the other half is totally transverse to the heat transfer direction, thus θ = (90 • + 70 • )/2 = 80 • .
Combining Equations (26) and (27), we can estimate the thermal conductivity of CVI-densified composites with anisotropic bundle and matrix at the different stage of densification. Because of the difficulties of determining the local conductivity orientation in the numerical mode, the model is only validated by our simulation by considering the bundles and the matrix as isotropic, and the comparison with the predictions of our model and numerical results are shown in Figure 10. It is found that the model agrees well with the simulation results when K b /K m is larger than 0.5, and slightly under-estimates the thermal conductivity when K b /K m is less than 0.5. As discussed in the Section 2, the matrix cracking may play an important role in the parallel direction of thermal transport. Then Lu and Hutchinson's model [38] (Equation (23)) should be adopted to estimate the reduction of thermal conductivity of the macro-pore structures, thus Equation (27) could be revised as: r is the equivalent radius of the bundle and can be seen as the geometric mean of the half length of the major and minor axis of the bundle cross-section. For the ease of estimation, the values of lc and Hc listed in Table 2 are also used for macro-pore structures.

Sensitive Analysis of Components' Parameters
In the previous sections, the models which relate the thermal conductivities of the CVI-densified composites with their constituent properties have been presented. However, the reliable estimation of thermal conductivity of these composites also depends on the precise thermal data and volume fractions of their components. However, most of the time these parameters may be not precisely known. Thus, before comparing with our models and the experimental data, the sensitivities of the predicted thermal conductivities to these parameters are assessed. The adopted thermal conductivities and volume fractions of the components are determined first and listed in Tables 6 and 7, respectively, which are obtained from the reliable sources or by the reasonable estimations. For instance, Taylor et al.'s experimental data [46] exhibits the axial thermal conductivity of T-300 carbon fiber would not changed after any heat treatment, if the temperature of the heat treatment is below 1400 K. Since the material are usually processed at around 1300 K, the thermal conductivity of T-300 carbon fiber at room temperature, 8 W·m −1 ·K −1 is safely retained [47]. According to Jumel et al.'s measurements [48], pan-based carbon fiber is transverse isotropic, the ratio of the axial thermal conductivity to the transverse one is 10.5 after any heat treatment, hence the transverse thermal conductivity of T-300 fiber is considered to be about 0.8 W·m −1 ·K −1 in this study. As discussed in the Section 2, the matrix cracking may play an important role in the parallel direction of thermal transport. Then Lu and Hutchinson's model [38] (Equation (23)) should be adopted to estimate the reduction of thermal conductivity of the macro-pore structures, thus Equation (27) could be revised as: where r b is the equivalent radius of the bundle and can be seen as the geometric mean of the half length of the major and minor axis of the bundle cross-section. For the ease of estimation, the values of l c and H c listed in Table 2 are also used for macro-pore structures.

Sensitive Analysis of Components' Parameters
In the previous sections, the models which relate the thermal conductivities of the CVI-densified composites with their constituent properties have been presented. However, the reliable estimation of thermal conductivity of these composites also depends on the precise thermal data and volume fractions of their components. However, most of the time these parameters may be not precisely known. Thus, before comparing with our models and the experimental data, the sensitivities of the predicted thermal conductivities to these parameters are assessed. The adopted thermal conductivities and volume fractions of the components are determined first and listed in Tables 6 and 7, respectively, which are obtained from the reliable sources or by the reasonable estimations. For instance, Taylor et al.'s experimental data [46] exhibits the axial thermal conductivity of T-300 carbon fiber would not changed after any heat treatment, if the temperature of the heat treatment is below 1400 K. Since the material are usually processed at around 1300 K, the thermal conductivity of T-300 carbon fiber at room temperature, 8 W·m −1 ·K −1 is safely retained [47]. According to Jumel et al.'s measurements [48], pan-based carbon fiber is transverse isotropic, the ratio of the axial thermal conductivity to the transverse one is 10.5 after any heat treatment, hence the transverse thermal conductivity of T-300 fiber is considered to be about 0.8 W·m −1 ·K −1 in this study. PyC interface layer displays a cylindrical symmetry in the transverse direction, hence the average values estimated by Katoh et al. [49] and Youngblood et al. [16] for the PyC interface layer have been retained. The estimated thermal data of CVI SiC matrix by Yamada et al. [50] is adopted. For the ease of the estimation, the spacing and thermal conductance of micro-cracks and macro-cracks are assumed to be equal and listed in Table 2. By inputting the listed data, changing only one parameter and fixing others, a detailed sensitivity study to model's parameters has been performed.  Table 7. Estimated volume fractions of the components for 1D, 2D and 3D composites.

Parameter Value
Volume fraction of fiber in 1D composites or in the bundles of 2D and 3D composites, V f 1D 55% Volume fraction of PyC interface layer in 1D composites or in the bundles of 2D and 3D composites, V i 1D 5% Volume fraction of micro-pores in 1D composites or in the bundles of 2D and 3D composites V p mic 8% Volume fraction of bundles in 2D and 3D composites V b 73% Volume fraction of macro-pores in 2D and 3D composites V p mac 5% We had performed the sensitive analysis for 1D and 2D composites reinforced by two fibers. It can be found that composites' thermal conductivities are not very sensitive to components' parameters, composites' thermal conductivities vary with these parameters only by less than 10%. The uncertainties of these parameters do not destroy the robustness of the model, their values can be utilized safely. The only exceptions are V p mic , V p mac , V b and K m , as shown by Tables 8-12.
These results demonstrate the uncertainties of these parameters play a significant influence on the composites' thermal conductivities, the fluctuation of the calculated thermal conductivities will be on the order of about 20%. Another foundation is the variation of composites' thermal conductivities will almost linearly vary with the relative difference of these parameters, since the reduced sensitivity index (RSI) do not vary very much with the relative difference of these parameters.    Except for the above parameters, some other parameters (e.g., the parameters of the matrix cracking) may vary in a much larger range, due to the uncontrollability of the processing. Hence, the sensitivity analysis in a larger range for these parameters are performed and shown in Figure 11. The results in Figure 11 demonstrate that the spacing between matrix cracks l c plays the most influence on the composites' thermal conductivities. Reducing l c by 90% will cause the reduction of composites' thermal conductivities by 30%-40%. However, increasing l c by one order will cause the increases of composites' thermal conductivities by 5%-10% and 20%-40% for 1D and 2D composites, respectively. These results demonstrate that the variance and contradiction of the thermal conductivity-residual porosity relation between different reports may partly due to the difficulty to distinguish the residual micro-porosity and macro-porosity, since the contributions of both to the thermal conductivity differ greatly. The second reason is the ignorance of the influence of the matrix cracking on the matrix thermal conductivity. thermal conductivity differ greatly. The second reason is the ignorance of the influence of the matrix cracking on the matrix thermal conductivity.

Comparison with Experimental Data
In Table 13, we summarize the available experimental data with the explicit geometrical information of various CVI-densified composites from various resources [13,15,16,35,53,54], adopt the thermal data of components listed in Table 6 and estimate the thermal conductivity of CVI-densified SiC matrix by our model. For the ease of estimation, the volume fraction of bundle and the volume fraction of the residual micro-porosity inside the bundle are utilizing the listed values in Table 7. It is noted that the measured 2D woven composites in Reference [13] stacked in 0°/30°/60°, which is slightly different from our geometrical model. We assume the difference of the stacking sequence on the thermal transport properties is not very large, and consider our model is still available. It is found that the estimated values by the present models are well consistent with the

Comparison with Experimental Data
In Table 13, we summarize the available experimental data with the explicit geometrical information of various CVI-densified composites from various resources [13,15,16,35,53,54], adopt the thermal data of components listed in Table 6 and estimate the thermal conductivity of CVI-densified SiC matrix by our model. For the ease of estimation, the volume fraction of bundle and the volume fraction of the residual micro-porosity inside the bundle are utilizing the listed values in Table 7. It is noted that the measured 2D woven composites in Reference [13] stacked in 0 • /30 • /60 • , which is slightly different from our geometrical model. We assume the difference of the stacking sequence on the thermal transport properties is not very large, and consider our model is still available. It is found that the estimated values by the present models are well consistent with the measured values. According to the discussion in Section 4, the departures may mainly attribute to the uncertainties of residual porosity in transverse direction and of the matrix cracking's density.

Conclusions
The numerical computations of the thermal conductivities of unidirectional, 2D woven and 3D braided composites at various stages of CVI densification are performed. The effects of the complex dual-scale pore structures on the thermal transport properties have been studied and two analytical thermal conductivity models with some empirical coefficients have been developed.
In general, there are seven adjustable parameters are uncertain or not available through experimental measurements and reported literatures: (i) Micro-porosity inside the bundle for the 2D and 3D composites, V p mac (ii) Macro-porosity between the bundle for the 2D and 3D composites, V p mac (iii) Volume fraction of the bundle for the 2D and 3D composites, V b (iv) Pore shape factor for the 2D and 3D composites, β (v) Effective fiber-interface and interface-matrix interfacial conductance, H 1,2 (vi) Thermal conductance of the matrix cracks, H c (vii) Spacing between matrix cracks, l c Through these parameters, the effective thermal conductivity of the CVi-densified composites have been found to be the most sensitive to the spacing between matrix cracks or matrix cracking's density, followed by volume fraction of the bundle and thermal conductance of the matrix cracks, finally by micro-porosity inside the bundles and macro-porosity between the bundles.
Although our model is developed in a simplified geometry, while the structures of actual composites are certainly more complex than these idealized structures, the presented models can provide insights into the magnitude of composites' thermal conductivities, provide an effective tool to assess the effects of all parameters on composites' thermal conductivities, and are useful for designing material processing and performance.