- Email: [email protected]
Finite element modeling of reactive liquid silicon infiltration
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
Contents lists available at ScienceDirect
Journal of the European Ceramic Society journal homepage: www.elsevier.com/locate/jeurceramsoc
Original Article
Finite element modeling of reactive liquid silicon infiltration Peter Josef Hofbauera, , Friedrich Raethera, Edda Rädleinb ⁎
a b
Fraunhofer Center for High Temperature Materials and Design, 95448 Bayreuth, Germany Department of Inorganic-Nonmetallic Materials, Technische Universität Ilmenau, 98684 Ilmenau, Germany
ARTICLE INFO
ABSTRACT
Keywords: Reactive melt infiltration Porous carbon preform Liquid silicon
The LSI process, i.e. the infiltration of molten silicon into porous structures, is one of the most economical techniques for the production of C/C-SiC and C/SiC ceramics. However, despite decades of development, the infiltration behavior affected by phenomena at the infiltration front has not been understood sufficiently. In the present work, a numerical model, based on the finite element method, was developed to simulate the infiltration process. The 3D model includes the penetration of silicon into the porous preform as well as the exothermal reactions at the infiltration front caused by the growth of SiC layers. For model validation, a special measuring furnace was used, enabling in situ optical inspection and weight measurement during liquid silicon infiltration into C/C-preforms in a controlled atmosphere. For the first time, a numerical model could be established which provides a tool to simulate the infiltration kinetics as well as the thermal processes during the LSI process in three dimensions. The model enables the optimization of melt infiltration processes with complex components within reasonable computer times.
1. Introduction In the recent work of the authors, a mesoscopic model for the infiltration of silicon into carbon capillaries was presented [1]. In contrast to many other publications, constant infiltration rates were found in model experiments using gap capillaries as well as porous C/C preforms. Basically, this was attributed to an evaporation – adsorption mechanism as the rate controlling step. Silicon vapor is formed at the infiltration front and adsorbed at the carbon pore walls close to it. This silicon reacts with carbon forming a thin layer of SiC, which is why the process can also be called chemisorption. Only after this SiC layer is covering the carbon, wetting of the pore walls by silicon becomes sufficient to allow proceeding of the infiltration front. This leads to a constant infiltration rate, as long as boundary conditions are steady [1]. These boundary conditions are: (i) constant temperature and pressure, (ii) an open pore – or micro crack system with an average pore diameter between 20 μm and 300 μm, (iii) infiltration heights smaller than the maximum capillary rise due to gravity, (iv) potential infiltration rate due to the capillary effect greater than the actual infiltration rate. Details of the experiments and the model as well as a survey on the literature on the LSI process are given in [1]. Regarding computer simulations of the LSI process, only few publications exist [2–7]. Most of them are based on the well-known Washburn equation [8], which implies a dependence of infiltration
⁎
height on the square root of time. Since the Washburn infiltration kinetics is different from real kinetics, a good match with experiments cannot be expected and indeed was not found indeed [4,7]. Moreover, existing models suffer from long computer times, which is attributed to very fine finite element (FE) meshes required. Therefore, simulations presented so far have been performed either in one or two dimensions [4,5,7], which does not allow simulations of infiltration processes with complex 3D components. In the following sections, we will present experimental results on the infiltration kinetics and a FE model which allows fast simulations considering steady infiltration rates. 2. Infiltration process The infiltration set-up to be simulated consisted of a special designed device for the infiltration of cylindrical porous C/C samples with a diameter of 3.2 cm and height of 3.0 cm, see Fig. 1. The samples were produced by machining larger C/C plates made by warm pressing of carbon short fiber bundles. For this, carbon fiber rovings were pressed and chopped to get short fiber bundles with a length of 6 mm, a width of 1 mm and a thickness of 0.2 mm [9]. Afterwards the short fiber bundles were mixed with phenolic resin as binder. Then the mixture was filled into a mould where it was uniaxially pressed at 150 °C to achieve the desired green density. The compact was then heated up to 200 °C to cure the phenolic resin [10]. A horizontal blind hole with a
Corresponding author. E-mail address: [email protected] (P.J. Hofbauer).
https://doi.org/10.1016/j.jeurceramsoc.2019.09.041 Received 29 May 2019; Received in revised form 18 September 2019; Accepted 24 September 2019 0955-2219/ © 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: Peter Josef Hofbauer, Friedrich Raether and Edda Rädlein, Journal of the European Ceramic Society, https://doi.org/10.1016/j.jeurceramsoc.2019.09.041
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
Fig. 2. Simplified CAD model of TOM_ac with the built-in experiment set-up. (1) Suspension rod. (2) Muffle. (3) Sample cage. (4) Optical beam path. (5) Graphite insulation. (6) Steel vessel.
atmosphere, oxygen respectively CO impurities were below 1 pm – typically 0.1 pm. In the examples shown below, the furnace atmosphere was argon with a purity of 99.999% and a pressure of 100 mbar. After reaching the target temperature, the weight sensor together with the graphite cage and the sample were lowered, until the bottom side of the wick was immersed into the liquid silicon. The point of immersion was registered by the weight sensor due to the buoyancy of the wick in the silicon melt. Thereafter, the weight change during silicon infiltration was monitored in dependence on time. In addition, the thermo-optical measuring arrangement also allowed an optical inspection with a camera during the reactive infiltration. Brightness measurements on the cylinder jacket and in the core of the C/C sample enabled local temperature recording. The effect of the number of drillings in the wick on infiltration rate was investigated in several infiltration experiments. It turned out, that infiltration kinetics did not change until the number of drillings was reduced below two holes. From that, it was concluded that in the present infiltration experiment with 43 drillings infiltration kinetics was controlled by the silicon transport within the C/C preform and not by the supply of silicon melt through the wick. During the silicon infiltration, the displacement of the infiltration front could be monitored at the lateral sample surface due to additional heat radiation resulting from fast exothermal silicon carbide formation. The starting point of infiltration was set at the moment when the silicon melt had passed the wick.
Fig. 1. Infiltration test set-up of a C/C specimen with silicon. – The drilling in the middle of the cylindrical C/C sample serves for the measurement of the temperature inside the specimen. (1) Suspension rod. (2) Sample cage made of graphite. (3) C/C sample with central drilling Ø 6 mm up to the cylinder axis. (4) Graphite wick with 43 cylindrical capillary drillings with Ø 0.6 mm in diameter. (5) Silicon granules. (6) Crucible made of graphite.
diameter of 6 mm ending in the center of the samples was prepared to monitor heat radiation during the infiltration (see Fig. 1). Temperature calibration was done by measuring the intensity of the heat radiation of C/C preforms via an IR camera in the same arrangement at different furnace temperatures without a reactive infiltration process. In the same way, temperature at the shell surface of the samples was measured and calibrated. The sample was supported by a cylindrical graphite wick with a length of 27 mm. In axial direction, 43 holes were drilled with a diameter of 0.6 mm each to supply the silicon melt via capillary forces. At the top side, the wick had a 0.2 mm deep, circular pocket to avoid throttling effects at the interface between the test sample and the wick. The wick was mounted in a graphite cage, see Fig. 1, which had a solid bottom plate to prevent undesired deposition of silicon vapor on the sample, before infiltration experiments were started. The graphite cage was fixed to the bottom side of a long rod suspended to the weight sensor at its upper side. A graphite crucible filled with silicon granules was placed some centimeters below the wick, as can be seen in Fig. 1. After installing the set-up in a special thermo-optical measuring furnace, named TOM_ac, the furnace was heated up under vacuum or in inert atmosphere to the desired infiltration temperature above the melting point of silicon (1413 °C) [11]. TOM_ac was equipped with graphite heaters and graphite insulation in a vacuum tight water-cooled vessel of stainless steel, see Fig. 2 in Section 2. Horizontal windows enabled an optical inspection of samples during the heat treatment. In addition, a weight sensor was mounted in a vacuum tight chamber above the furnace. This chamber was connected with the furnace using vacuum bellows. The graphite cage was hang-up via a suspension rod at the weight sensor. This made gravimetric measurements possible on the one hand, and the sample could be moved up and down during the experiments with a linear guiding on the other hand. The vacuum chamber was evacuated through the bellows during evacuation of the vessel. Maximum temperature of the furnace was 2200 °C. In inert
3. Finite element method 3.1. Infiltration kinetics Since the infiltration rate was constant in the present experiments, it was mathematically described as a chemical reaction with spatial propagation, which is expressed by the Fisher equation, a special case of the reaction diffusion equation [12]: 2u u = D 2 + k u (C t x
u),
(1)
with D denoting the diffusion coefficient, k the rate constant for the reaction, C the maximum concentration and u the time dependent local concentration [12]. Eq. (1) was originally used by Luther to describe autocatalytic, isotherm reactions in homogeneous media. He found out that the spatial spreading of the reaction is given by a superposition of diffusion and chemical reaction, where the first term on the right hand 2
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
side accounts for the diffusion and the second term for the autocatalytic reaction [13,14]:
du = k u (C dt
u)
(2)
The solution of Eq. (2), also called the logistic equation, shows a sigmoidal course in the direction perpendicular to the reaction front. Due to catalyst formation, initially the concentration curve increases exponentially and then flattens on account of the material consumption from the reaction [14,15]. To describe the propagation of the reaction front, the well-known formulation of displacement x − v t can be used [16]. However, this leads to the fact that stable solutions exist only for infiltration rates greater or equal to the minimum rate v with
v
(3)
4k C D.
The problem of the minimum spreading rate was already addressed by Luther and has so far been unsolved [13,17]. If a few assumptions are made, the Fisher equation can also be used to describe the infiltration front of liquid silicon in porous carbon preforms. For this purpose, it is assumed, that the concentration c behaves analogously to the degree of saturation α of the pores. Thus, the maximum concentration, c = 1, corresponds to the infiltrated region and the concentration c = 0 to the not infiltrated region. The reaction constant k describes the sharpness of the transition from the unsaturated to the saturated domain. To solve the aforementioned problem of the smallest possible propagation respectively infiltration velocity, Eq. (3) was solved for the diffusion coefficient and used in Eq. (1):
t
=
2 vinf 4k
+k
(1
).
Fig. 4. Convergence behavior and computing time as a function of mesh quality and k-value. – An infiltration rate of 10−3 m/s was used for the computation and the α-value was determined after 100 s at an infiltration height of 100 mm.
distribution results from twisted pore channels in the preform leading to locally different infiltration path lengths and – from a macroscopic view – to a smearing out of the infiltration front. Numerically the kvalue is simply adapted to the mesh quality since the real distribution is unknown. Very fine meshes as customarily used for the simulation of infiltration processes are not necessary since the overall kinetics does not change significantly. Thus, it is possible to calculate the infiltration behavior of complex components in a short time. Fig. 4 illustrates the convergence of α as a function of the element size, which was used in a uniform mesh. Numerical calculations were performed using Intel i73820 processor, with 4 cores at a frequency of 2.70 GHz. And finally, the solution function α(x, y, z, t), which returns the value 1 for the infiltrated area and the value 0 for the uninfiltrated area, can be used to assign the material properties E for the different regions of the preform:
(4)
Using a stability analysis, it can be shown that stable solutions exist for the new formulated infiltration Eq. (4) at any rate vinf . The procedure of stability analysis can be found in mathematical text books, e.g. in [16]. Fig. 3 shows the numerically calculated solutions for the one-dimensional case of Eq. (4) using different infiltration rates and times. The width of the sigmoidal function decreases with increasing parameter k. With t = x / vinf Eq. (2) gives the slope k / vinf of the sigmoidal curve at the inflection point. From a physical point of view, parameter k describes the fuzziness of the silicon distribution at the infiltration front. Mainly, this
Esample =
Einfiltrated + (1
) Enon-infiltrated.
(5)
3.2. Time and temperature dependent material parameter Thermal diffusivity of the C/C preforms used in the infiltration experiments was measured up to temperatures of 2000 °C using a thermooptical measuring device called TOM_1 which combines a special laserflash technique with optical dilatometry [18]. In addition, thermal diffusivity was measured at lower temperatures with a commercial laser-flash device (Netzsch LFA 457). Due to the manufacturing process, these preforms have transversely isotropic material properties: equal properties in the horizontal plane and different properties in the z-direction – the direction of warm pressing. Therefore, two independent laser-flash measurements with the heat flow in z-direction and in xydirection were performed. Results shown in Fig. 5 reflect this anisotropy. Thermal expansion was measured by optical dilatometry [11]. From the temperature dependent change of length, the coefficient of longitudinal thermal expansion for an orthotropic material can be calculated by a symmetrical tensor with its components αij(T): ij (T )
=
ij (T )
T
,
(6)
with the strain ϵij(T) and for the temperature change ΔT [19,20]. Using the transversal isotropy α11 = α22, the volumetric thermal expansion coefficient γ(T) and the density ρ(T) were calculated:
Fig. 3. Degree of saturation α versus coordinate x for one-dimensional numerical solutions of Eq. (4) with different rates v and parameter k at different times t.
(T ) = 3
V =2 T
11
+
33
(7)
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
thermodynamic databases and c AR was determined to be 56.3 × 103 m2/m3 by pore analysis of an μ X-ray computer tomograph (Diondo 225 kV microfocus tube) of the C/C sample. With increasing reaction time, the SiC crystals within the SiC layer coarsen due to Ostwald Ripening. Thus, within the SiC layer the ratio of grain boundary diffusion to lattice diffusion decreases. The interdiffusion of Si and C species was considered by a time-dependent apparent diffusion coefficient Dapp(t) which was defined by: SiC
t
=
Dapp (t ) (16)
SiC
Dapp (t ) = D0,app (t ) exp
EA (t ) RT
(17)
where both the preexponential factor D0,app(t) and activation energy EA(t) depend on time. The two functions were obtained by evaluation of literature and measurement data:
(T ) =
exp
(
EA (t ) = 260
)
T
(T ) dT .
T0
(8)
(9)
= aij(T ) (T ) Cp (T ),
Ph 1 EPh 1
+
Ph 2
EPh 2 +
+
Ph n
EPh n .
C/C (T )
=
C/C
Cp,C/SiC (T ) = C/SiC (T )
1
C/C
C/C
Si
Cp,C/C (T ) + C/C (T )
=
+
1
+
Si
Si
+
Si (T )
SiC
Cp,Si (T ) +
Si (T )
+
SiC
SiC (T )
SiC
Cp,SiC (T )
SiC (T )
(19)
First of all, a simplified geometry derived from a computer aided design (CAD) model of the furnace and infiltration set-up was generated using the commercial CAD software AutoCAD Inventor. Then, the threedimensional geometry of the sample was divided into smaller sections to generate subdomains using the FE software COMSOL Multiphysics. In this way, the subdomains could be meshed with quadrangular elements at the subdomains surface. The meshed geometry is shown in Fig. 7. Within the subdomains a finer meshed submodel with tetrahedral elements was generated. Care was taken, to ensure that the tetrahedral elements were grouped in such a way, that their edge nodes coincided with the corresponding master nodes of the elements of the subdomain. Thereby the tetrahedral element groups were coupled to the master
(10)
Where EComp is the material property of the composite material, EPh n is the property of the particular phase and φPh n is the corresponding volume fraction of this phase. Thus, the material properties are given by: C/SiC (T )
kJ 110s exp . mol t + 185s
3.3. Discretization
where Cp(T) is the specific heat capacity of the sample [21,19,20]. The change of material properties during infiltration of silicon and formation of silicon carbide were estimated using the rule of mixture [22–24]:
EComp
(18)
Results of the numerical integration of Eq. (16) are shown in Fig. 6. A good agreement was found with the data from Zhou and Singh [25] as well as Favre et al. [26] at T = 1510 °C respectively T = 1600 °C.
With the measured thermal diffusivity aij(T), the thermal conductivity λij(T) was calculated: ij (T )
cm2 1700s exp s t + 86s
D0,app (t ) = 4.09
Fig. 5. Thermal diffusivity a T measurement of an orthotropic C/C material with transversal isotropy.
(11) (12) (13)
with ρC/C and λC/C as bulk density and thermal diffusivity, respectively, as well as the porosity Φ. Literature data were used for specific heat of all phases as well as density and thermal conductivity of silicon and silicon carbide. The heat Q generated due to the exothermal reaction of carbon and silicon to silicon carbide is given by:
Q = H c AR
SiC
t
s(
SiC ),
(14)
where ΔH and c AR are the enthalpy of reaction and the reactive surface per volume respectively. SiC is the increase of the thickness of the SiC t layer between carbon and silicon melt per time and s(δSiC) is a step function taking into account that the reaction stops, when either no carbon or no silicon are available in the corresponding pore channel:
c AR SiC 2 s ( SiC) = c AR SiC 1 for 2 0 else 1 for
0, Si
, 0, C/C
(15) Fig. 6. Numerical integration of Eq. (16) compared with literature data from Zhou and Singh [25] as well as Favre et al. [26].
where φ0,Si and φ0,C/C are the volume fractions of silicon and carbon before the formation of the SiC layer. ΔH was obtained from 4
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
4. Computational results and validation 4.1. Infiltration kinetics Fig. 9 shows the simulation of the propagating infiltration front with a constant velocity of 1 × 10−3 m/s based on Eq. (4). The shape and movement of the infiltration front reflect the infiltration behavior which was experimentally observed. Using Eqs. (11) and (16), the increase of mass during infiltration was calculated. The comparison with the gravimetrically measured data shows a good correlation between both curves, see Fig. 10. However, there is a small deviation at point 2 in Fig. 10, when the measured infiltration is faster than the simulation. This is attributed to a slight asymmetry in the set-up, which leads to variations of the rate of mass increase in the initial state of infiltration, when the hemispherical infiltration front arrives at the lateral surface of the cylinder. 4.2. Temperature fields Heat generation depends on the region which is already infiltrated and on the position of the infiltration front. It is calculated according to Eqs. (14)–(19). Due to the moving infiltration front, a dynamic interplay between heat generation and heat conduction within the sample and heat radiation at its surfaces occurs, which was considered in the FE model. Fig. 11 shows an example of local heat generation and temperature distribution after 18 s of silicon infiltration. As already observed experimentally, heat generation is especially strong at the infiltration front and less pronounced in the infiltrated region, see Fig. 11a, which reflects the decreasing formation rate of SiC with increasing layer thickness, see Eq. (16). Due to heat transfer to cooler regions, the temperature distribution within the sample is smoother than the heat generation rate, compare Fig. 11b. Fig. 12 shows the measured and simulated temperatures at two fixed points during infiltration: the cylinder jacket and at the core of the sample. The latter was measured in the blind hole drilled into the sample, compare Fig. 1. To check the reproducibility, two infiltration experiments were performed. It can be seen that the core temperature rises from the start temperature at 1550 °C to a maximum at 2150 °C within 20 s and decreases slowly thereafter. Even at the surface of the sample, the temperature increases to 1990 °C. There is some scatter in the experimental curves which is attributed to the coarse structure of the C/C samples. Considering this scatter, the correlation between calculated and measured temperature curves is excellent.
Fig. 7. Meshed FE model used in the simulations. – By using a symmetry plane, the furnace as well as the experimental set-up can be meshed in the half-model.
nodes of the subdomain elements, also called superelements [27,28]. By that, heat radiation could be calculated in the global model and the infiltration equation as well as the moving heat source in the submodel. Since both, the infiltration of the sample by the silicon melt and the transport of the silicon from the crucible into the sample were simulated, moving and deformed meshes were present in the submodels. Using the superelements reduced the recalculation of view factors during the simulations and saved computing time. Both, the global model, meshed with superelements and the submodel, meshed with tetrahedral elements, are shown in Fig. 8.
Fig. 8. Meshed global model and submodel of the sample. – Fig. 8a shows an enlarged section of Fig. 7. 5
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
Fig. 9. Infiltration front after 13 s, 20 s, 27 s and 34 s based on a constant infiltration velocity of 1 × 10−3 m/s. – The grid lines show the edges of the subdomains for discretization of the geometry provided by the CAD system. – The figure shows the meshed sample in Fig. 8b.
4.3. Convergence and sensitivity analysis For an estimation of the convergence properties of the FE model, computations were performed with different mesh sizes, time steps and tolerances for the iteration errros. Infiltration velocity was set to 1 × 10−3 m/s and k-value to 0.5. The maximum temperature increase due to the exothermic reaction in the sample was used for comparison. Using a mesh size of 0.8 mm, a time step of 0.1 s and a tolerance of 1 × 10−3 °C for the iteration error results in a variation of the maximum temperature of ± 2 °C around its mean value of 2248 °C, see also Fig. 4. Only the mesh size of the preform in the submodel was refined, the rest of the model remained unchanged, see Fig. 8. To verify how the model reacts to small changes in the measured infiltration rate and in the reactive surface, a sensitivity analysis was performed. The maximum temperature increased by 5 °C at a 10% increase in the infiltration rate from 1.0 mm/s to 1.1 mm/s and increased by 26 °C at a 10% increase in the reactive surface from 56.3 × 103 m2/m3 to 61.8 × 103 m2/m3. However, as the model validation shows, the uncertainty of the numerical computations are in the same range as the uncertainty resulting from the reproducibility of the measurements, see Fig. 12.
Fig. 10. Increase of mass during the infiltration of a C/C preform with liquid silicon. – The insert shows a cross section through the sample indicating the infiltration front at different times which are also marked on the curve.
6
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
Fig. 11. Heat generation and temperature field during specimen infiltration after 18 s.
in a short computer time. This is especially important, when industrial infiltration processes are optimized and heating schedules and stacks have to be varied many times. Using the present model, the LSI process for the production of ceramic break disks was successfully optimized. For that, other effects like the melting and solidifying of silicon during heating up and cooling of the charge also had to be integrated into the model. Principles of the entire model will be presented in a further work. Acknowledgement The authors gratefully acknowledge the help of P. Döppmann and J. Baber with the measurements and the financial support of the Bayerische Forschungsstiftung within the project ISE-LSI. References [1] P.J. Hofbauer, E. Rädlein, F. Raether, Fundamental mechanisms with reactive infiltration of silicon melt into carbon capillaries, Adv. Eng. Mater. 1900184 (2019) 1–11. [2] W.B. Hillig, Melt infiltration approach to ceramic matrix composites, J. Am. Ceram. Soc. 71 (2) (1988) C-96–C-99. [3] E.O. Einset, Capillary infiltration rates into porous media with applications to silcomp processing, J. Am. Ceram. Soc. 79 (2) (1996) 333–338. [4] E.O. Einset, Analysis of reactive melt infiltration in the processing of ceramics and ceramic composites, Chem. Eng. Sci. 53 (5) (1998) 1027–1039. [5] E.S. Nelson, P. Colella, Parametric Study of Reactive Melt Infiltration, Tech. Rep. E12107, National Aeronautics and Space Administration, 2000. [6] P. Sangsuwan, S.N. Tewari, J.E. Catica, M. Singh, R. Dickerson, Reactive infiltration of silicon melt through microporous amorphous carbon preforms, Metall. Mater. Trans. B 30 (5) (1999) 933–944. [7] S.P. Yushanov, J.S. Crompton, K.C. Koppenhoefer, Simulation of Manufacturing Process of Ceramic Matrix Composites, (2008). [8] E.W. Washburn, The dynamics of capillary flow, Phys. Rev. 17 (3) (1921) 273–283. [9] A. Kienzle, I. Kratschmer, Friction-Tolerant Disks Made of Fiber-Reinforced Ceramic, (2015). [10] A. Kienzle, H. Jäger, Carbon-Fiber-Reinforced Silicon Carbide: A New Brake Disk Material, (2015). [11] F.G. Raether, Current state of in situ measuring methods for the control of firing processes, J. Am. Ceram. Soc. 92 (S1) (2009) 146–152. [12] R.A. Fisher, The wave of advance of advantageous genes, Ann. Hum. Genet. 7 (4) (1937) 355–369. [13] R. Luther, Räumliche Fortpflanzung chemischer Reaktionen, Z. für Elektrochem. Angew. Phys. Chem. 12 (32) (1906) 569–600. [14] F.T.K. Meinecke, Räumliche Fortpflanzung chemischer Reaktionen, Inauguraldissertation, 1908. [15] R.K. Hobbie, B.J. Roth, Intermediate Physics for Medicine and Biology, SpringerVerlag, New York, 2007. [16] J.D. Murray, Mathematical Biology: I. An Introduction, 3rd Edition, Vol. 17 of Interdisciplinary Applied Mathematics, Springer-Verlag Berlin Heidelberg, 2002. [17] E. Brunet, B. Derrida, An exactly solvable travelling wave equation in the FisherKPP class, J. Stat. Phys. 161 (4) (2015) 801–820. [18] F. Raether, R. Hofmann, G. Müller, H.J. Sölter, A novel thermo-optical measuring system for the in situ study of sintering processes, J. Therm. Anal. Calorim. 53 (3) (1998) 717–735. [19] C. Dudescu, J. Naumann, M. Stockmann, S. Nebel, Characterisation of thermal
Fig. 12. Measured (two infiltration experiments) and calculated temperature increase at two different points during silicon infiltration into a C/C preform. – The insert in the upper right corner shows the points of temperature detection.
5. Conclusions A FE model was developed for the simulation of the infiltration kinetics and the corresponding thermal effects during infiltration of silicon melt into porous C/C preforms. A good match to the characteristics of the infiltration process obtained from in situ measurements with the measuring furnace TOM_ac was obtained. This good correlation was achieved using material data from own measurements and from literature as a basis for the simulation. Important parameters, which are specific for any type of C/C preform, are the reactive surface area, thermal conductivity and infiltration rate. Note that the latter also depends on furnace parameters like gas pressure and temperature [1]. It is recommended to measure these parameters individually for each used type of C/C material. If the reactive infiltration into other carbonaceous preforms or pyrolized wood is to be simulated, Eqs. (17) and (19) have to adapted as well. The good correlation between experiment and simulation was also caused by the use of time dependent diffusion coefficients based on these two equations. In contrast to previous models using capillary or diffusion equations for description of infiltration kinetics, a differential equation for autocatalytic processes was developed for this purpose. This new infiltration equation can be adapted to the mesh quality of FE models and thus enables 3D simulation of infiltration processes of complex components 7
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
[20] [21] [22] [23] [24]
expansion coefficient of anisotropic materials by electronic speckle pattern interferometry, Strain 42 (3) (2006) 197–205. O. Zienkiewicz, R. Taylor, J. Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013. Z. Hashin, S. Shtrikman, On some variational principles in anisotropic and nonhomogeneous elasticity, J. Mech. Phys. Solids 10 (4) (1962) 335–342. M.F. Ashby, Materials Selection in Mechanical Design, 3rd ed., ButterworthHeinemann, Oxford, 2005. J.K. Farooqi, M.A. Skeikh, Finite Element Modeling of Thermal Transport in Composite Unit Cell Using Abaqus/CAE, (2005). Z. Fu, L. Schlier, N. Travitzky, P. Greil, Three-dimensional printing of sisic lattice
truss structures, Mater. Sci. Eng. A 560 (2013) 851–856. [25] H. Zhou, R.N. Singh, Kinetics model for the growth of silicon carbide by the reaction of liquid silicon with carbon, J. Am. Ceram. Soc. 78 (9) (1995) 7. [26] A. Favre, H. Fuzellier, J. Suptil, An original way to investigate the siliconizing of carbon materials, Ceram. Int. 29 (2003) 235–243. [27] J. Plaza, M. Abasolo, I. Coria, J. Aguirrebeitia, I. Fernández de Bustos, A new finite element approach for the analysis of slewing bearings in wind turbine generators using superelement techniques, Meccanica 50 (6) (2015) 1623–1633. [28] T. Panczak, S. Ring, M. Welch, A cad-based tool for fdm and fem radiation and conduction modeling, SAE Trans. 107 (1) (1998) 363–372.
8
Contents lists available at ScienceDirect
Journal of the European Ceramic Society journal homepage: www.elsevier.com/locate/jeurceramsoc
Original Article
Finite element modeling of reactive liquid silicon infiltration Peter Josef Hofbauera, , Friedrich Raethera, Edda Rädleinb ⁎
a b
Fraunhofer Center for High Temperature Materials and Design, 95448 Bayreuth, Germany Department of Inorganic-Nonmetallic Materials, Technische Universität Ilmenau, 98684 Ilmenau, Germany
ARTICLE INFO
ABSTRACT
Keywords: Reactive melt infiltration Porous carbon preform Liquid silicon
The LSI process, i.e. the infiltration of molten silicon into porous structures, is one of the most economical techniques for the production of C/C-SiC and C/SiC ceramics. However, despite decades of development, the infiltration behavior affected by phenomena at the infiltration front has not been understood sufficiently. In the present work, a numerical model, based on the finite element method, was developed to simulate the infiltration process. The 3D model includes the penetration of silicon into the porous preform as well as the exothermal reactions at the infiltration front caused by the growth of SiC layers. For model validation, a special measuring furnace was used, enabling in situ optical inspection and weight measurement during liquid silicon infiltration into C/C-preforms in a controlled atmosphere. For the first time, a numerical model could be established which provides a tool to simulate the infiltration kinetics as well as the thermal processes during the LSI process in three dimensions. The model enables the optimization of melt infiltration processes with complex components within reasonable computer times.
1. Introduction In the recent work of the authors, a mesoscopic model for the infiltration of silicon into carbon capillaries was presented [1]. In contrast to many other publications, constant infiltration rates were found in model experiments using gap capillaries as well as porous C/C preforms. Basically, this was attributed to an evaporation – adsorption mechanism as the rate controlling step. Silicon vapor is formed at the infiltration front and adsorbed at the carbon pore walls close to it. This silicon reacts with carbon forming a thin layer of SiC, which is why the process can also be called chemisorption. Only after this SiC layer is covering the carbon, wetting of the pore walls by silicon becomes sufficient to allow proceeding of the infiltration front. This leads to a constant infiltration rate, as long as boundary conditions are steady [1]. These boundary conditions are: (i) constant temperature and pressure, (ii) an open pore – or micro crack system with an average pore diameter between 20 μm and 300 μm, (iii) infiltration heights smaller than the maximum capillary rise due to gravity, (iv) potential infiltration rate due to the capillary effect greater than the actual infiltration rate. Details of the experiments and the model as well as a survey on the literature on the LSI process are given in [1]. Regarding computer simulations of the LSI process, only few publications exist [2–7]. Most of them are based on the well-known Washburn equation [8], which implies a dependence of infiltration
⁎
height on the square root of time. Since the Washburn infiltration kinetics is different from real kinetics, a good match with experiments cannot be expected and indeed was not found indeed [4,7]. Moreover, existing models suffer from long computer times, which is attributed to very fine finite element (FE) meshes required. Therefore, simulations presented so far have been performed either in one or two dimensions [4,5,7], which does not allow simulations of infiltration processes with complex 3D components. In the following sections, we will present experimental results on the infiltration kinetics and a FE model which allows fast simulations considering steady infiltration rates. 2. Infiltration process The infiltration set-up to be simulated consisted of a special designed device for the infiltration of cylindrical porous C/C samples with a diameter of 3.2 cm and height of 3.0 cm, see Fig. 1. The samples were produced by machining larger C/C plates made by warm pressing of carbon short fiber bundles. For this, carbon fiber rovings were pressed and chopped to get short fiber bundles with a length of 6 mm, a width of 1 mm and a thickness of 0.2 mm [9]. Afterwards the short fiber bundles were mixed with phenolic resin as binder. Then the mixture was filled into a mould where it was uniaxially pressed at 150 °C to achieve the desired green density. The compact was then heated up to 200 °C to cure the phenolic resin [10]. A horizontal blind hole with a
Corresponding author. E-mail address: [email protected] (P.J. Hofbauer).
https://doi.org/10.1016/j.jeurceramsoc.2019.09.041 Received 29 May 2019; Received in revised form 18 September 2019; Accepted 24 September 2019 0955-2219/ © 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: Peter Josef Hofbauer, Friedrich Raether and Edda Rädlein, Journal of the European Ceramic Society, https://doi.org/10.1016/j.jeurceramsoc.2019.09.041
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
Fig. 2. Simplified CAD model of TOM_ac with the built-in experiment set-up. (1) Suspension rod. (2) Muffle. (3) Sample cage. (4) Optical beam path. (5) Graphite insulation. (6) Steel vessel.
atmosphere, oxygen respectively CO impurities were below 1 pm – typically 0.1 pm. In the examples shown below, the furnace atmosphere was argon with a purity of 99.999% and a pressure of 100 mbar. After reaching the target temperature, the weight sensor together with the graphite cage and the sample were lowered, until the bottom side of the wick was immersed into the liquid silicon. The point of immersion was registered by the weight sensor due to the buoyancy of the wick in the silicon melt. Thereafter, the weight change during silicon infiltration was monitored in dependence on time. In addition, the thermo-optical measuring arrangement also allowed an optical inspection with a camera during the reactive infiltration. Brightness measurements on the cylinder jacket and in the core of the C/C sample enabled local temperature recording. The effect of the number of drillings in the wick on infiltration rate was investigated in several infiltration experiments. It turned out, that infiltration kinetics did not change until the number of drillings was reduced below two holes. From that, it was concluded that in the present infiltration experiment with 43 drillings infiltration kinetics was controlled by the silicon transport within the C/C preform and not by the supply of silicon melt through the wick. During the silicon infiltration, the displacement of the infiltration front could be monitored at the lateral sample surface due to additional heat radiation resulting from fast exothermal silicon carbide formation. The starting point of infiltration was set at the moment when the silicon melt had passed the wick.
Fig. 1. Infiltration test set-up of a C/C specimen with silicon. – The drilling in the middle of the cylindrical C/C sample serves for the measurement of the temperature inside the specimen. (1) Suspension rod. (2) Sample cage made of graphite. (3) C/C sample with central drilling Ø 6 mm up to the cylinder axis. (4) Graphite wick with 43 cylindrical capillary drillings with Ø 0.6 mm in diameter. (5) Silicon granules. (6) Crucible made of graphite.
diameter of 6 mm ending in the center of the samples was prepared to monitor heat radiation during the infiltration (see Fig. 1). Temperature calibration was done by measuring the intensity of the heat radiation of C/C preforms via an IR camera in the same arrangement at different furnace temperatures without a reactive infiltration process. In the same way, temperature at the shell surface of the samples was measured and calibrated. The sample was supported by a cylindrical graphite wick with a length of 27 mm. In axial direction, 43 holes were drilled with a diameter of 0.6 mm each to supply the silicon melt via capillary forces. At the top side, the wick had a 0.2 mm deep, circular pocket to avoid throttling effects at the interface between the test sample and the wick. The wick was mounted in a graphite cage, see Fig. 1, which had a solid bottom plate to prevent undesired deposition of silicon vapor on the sample, before infiltration experiments were started. The graphite cage was fixed to the bottom side of a long rod suspended to the weight sensor at its upper side. A graphite crucible filled with silicon granules was placed some centimeters below the wick, as can be seen in Fig. 1. After installing the set-up in a special thermo-optical measuring furnace, named TOM_ac, the furnace was heated up under vacuum or in inert atmosphere to the desired infiltration temperature above the melting point of silicon (1413 °C) [11]. TOM_ac was equipped with graphite heaters and graphite insulation in a vacuum tight water-cooled vessel of stainless steel, see Fig. 2 in Section 2. Horizontal windows enabled an optical inspection of samples during the heat treatment. In addition, a weight sensor was mounted in a vacuum tight chamber above the furnace. This chamber was connected with the furnace using vacuum bellows. The graphite cage was hang-up via a suspension rod at the weight sensor. This made gravimetric measurements possible on the one hand, and the sample could be moved up and down during the experiments with a linear guiding on the other hand. The vacuum chamber was evacuated through the bellows during evacuation of the vessel. Maximum temperature of the furnace was 2200 °C. In inert
3. Finite element method 3.1. Infiltration kinetics Since the infiltration rate was constant in the present experiments, it was mathematically described as a chemical reaction with spatial propagation, which is expressed by the Fisher equation, a special case of the reaction diffusion equation [12]: 2u u = D 2 + k u (C t x
u),
(1)
with D denoting the diffusion coefficient, k the rate constant for the reaction, C the maximum concentration and u the time dependent local concentration [12]. Eq. (1) was originally used by Luther to describe autocatalytic, isotherm reactions in homogeneous media. He found out that the spatial spreading of the reaction is given by a superposition of diffusion and chemical reaction, where the first term on the right hand 2
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
side accounts for the diffusion and the second term for the autocatalytic reaction [13,14]:
du = k u (C dt
u)
(2)
The solution of Eq. (2), also called the logistic equation, shows a sigmoidal course in the direction perpendicular to the reaction front. Due to catalyst formation, initially the concentration curve increases exponentially and then flattens on account of the material consumption from the reaction [14,15]. To describe the propagation of the reaction front, the well-known formulation of displacement x − v t can be used [16]. However, this leads to the fact that stable solutions exist only for infiltration rates greater or equal to the minimum rate v with
v
(3)
4k C D.
The problem of the minimum spreading rate was already addressed by Luther and has so far been unsolved [13,17]. If a few assumptions are made, the Fisher equation can also be used to describe the infiltration front of liquid silicon in porous carbon preforms. For this purpose, it is assumed, that the concentration c behaves analogously to the degree of saturation α of the pores. Thus, the maximum concentration, c = 1, corresponds to the infiltrated region and the concentration c = 0 to the not infiltrated region. The reaction constant k describes the sharpness of the transition from the unsaturated to the saturated domain. To solve the aforementioned problem of the smallest possible propagation respectively infiltration velocity, Eq. (3) was solved for the diffusion coefficient and used in Eq. (1):
t
=
2 vinf 4k
+k
(1
).
Fig. 4. Convergence behavior and computing time as a function of mesh quality and k-value. – An infiltration rate of 10−3 m/s was used for the computation and the α-value was determined after 100 s at an infiltration height of 100 mm.
distribution results from twisted pore channels in the preform leading to locally different infiltration path lengths and – from a macroscopic view – to a smearing out of the infiltration front. Numerically the kvalue is simply adapted to the mesh quality since the real distribution is unknown. Very fine meshes as customarily used for the simulation of infiltration processes are not necessary since the overall kinetics does not change significantly. Thus, it is possible to calculate the infiltration behavior of complex components in a short time. Fig. 4 illustrates the convergence of α as a function of the element size, which was used in a uniform mesh. Numerical calculations were performed using Intel i73820 processor, with 4 cores at a frequency of 2.70 GHz. And finally, the solution function α(x, y, z, t), which returns the value 1 for the infiltrated area and the value 0 for the uninfiltrated area, can be used to assign the material properties E for the different regions of the preform:
(4)
Using a stability analysis, it can be shown that stable solutions exist for the new formulated infiltration Eq. (4) at any rate vinf . The procedure of stability analysis can be found in mathematical text books, e.g. in [16]. Fig. 3 shows the numerically calculated solutions for the one-dimensional case of Eq. (4) using different infiltration rates and times. The width of the sigmoidal function decreases with increasing parameter k. With t = x / vinf Eq. (2) gives the slope k / vinf of the sigmoidal curve at the inflection point. From a physical point of view, parameter k describes the fuzziness of the silicon distribution at the infiltration front. Mainly, this
Esample =
Einfiltrated + (1
) Enon-infiltrated.
(5)
3.2. Time and temperature dependent material parameter Thermal diffusivity of the C/C preforms used in the infiltration experiments was measured up to temperatures of 2000 °C using a thermooptical measuring device called TOM_1 which combines a special laserflash technique with optical dilatometry [18]. In addition, thermal diffusivity was measured at lower temperatures with a commercial laser-flash device (Netzsch LFA 457). Due to the manufacturing process, these preforms have transversely isotropic material properties: equal properties in the horizontal plane and different properties in the z-direction – the direction of warm pressing. Therefore, two independent laser-flash measurements with the heat flow in z-direction and in xydirection were performed. Results shown in Fig. 5 reflect this anisotropy. Thermal expansion was measured by optical dilatometry [11]. From the temperature dependent change of length, the coefficient of longitudinal thermal expansion for an orthotropic material can be calculated by a symmetrical tensor with its components αij(T): ij (T )
=
ij (T )
T
,
(6)
with the strain ϵij(T) and for the temperature change ΔT [19,20]. Using the transversal isotropy α11 = α22, the volumetric thermal expansion coefficient γ(T) and the density ρ(T) were calculated:
Fig. 3. Degree of saturation α versus coordinate x for one-dimensional numerical solutions of Eq. (4) with different rates v and parameter k at different times t.
(T ) = 3
V =2 T
11
+
33
(7)
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
thermodynamic databases and c AR was determined to be 56.3 × 103 m2/m3 by pore analysis of an μ X-ray computer tomograph (Diondo 225 kV microfocus tube) of the C/C sample. With increasing reaction time, the SiC crystals within the SiC layer coarsen due to Ostwald Ripening. Thus, within the SiC layer the ratio of grain boundary diffusion to lattice diffusion decreases. The interdiffusion of Si and C species was considered by a time-dependent apparent diffusion coefficient Dapp(t) which was defined by: SiC
t
=
Dapp (t ) (16)
SiC
Dapp (t ) = D0,app (t ) exp
EA (t ) RT
(17)
where both the preexponential factor D0,app(t) and activation energy EA(t) depend on time. The two functions were obtained by evaluation of literature and measurement data:
(T ) =
exp
(
EA (t ) = 260
)
T
(T ) dT .
T0
(8)
(9)
= aij(T ) (T ) Cp (T ),
Ph 1 EPh 1
+
Ph 2
EPh 2 +
+
Ph n
EPh n .
C/C (T )
=
C/C
Cp,C/SiC (T ) = C/SiC (T )
1
C/C
C/C
Si
Cp,C/C (T ) + C/C (T )
=
+
1
+
Si
Si
+
Si (T )
SiC
Cp,Si (T ) +
Si (T )
+
SiC
SiC (T )
SiC
Cp,SiC (T )
SiC (T )
(19)
First of all, a simplified geometry derived from a computer aided design (CAD) model of the furnace and infiltration set-up was generated using the commercial CAD software AutoCAD Inventor. Then, the threedimensional geometry of the sample was divided into smaller sections to generate subdomains using the FE software COMSOL Multiphysics. In this way, the subdomains could be meshed with quadrangular elements at the subdomains surface. The meshed geometry is shown in Fig. 7. Within the subdomains a finer meshed submodel with tetrahedral elements was generated. Care was taken, to ensure that the tetrahedral elements were grouped in such a way, that their edge nodes coincided with the corresponding master nodes of the elements of the subdomain. Thereby the tetrahedral element groups were coupled to the master
(10)
Where EComp is the material property of the composite material, EPh n is the property of the particular phase and φPh n is the corresponding volume fraction of this phase. Thus, the material properties are given by: C/SiC (T )
kJ 110s exp . mol t + 185s
3.3. Discretization
where Cp(T) is the specific heat capacity of the sample [21,19,20]. The change of material properties during infiltration of silicon and formation of silicon carbide were estimated using the rule of mixture [22–24]:
EComp
(18)
Results of the numerical integration of Eq. (16) are shown in Fig. 6. A good agreement was found with the data from Zhou and Singh [25] as well as Favre et al. [26] at T = 1510 °C respectively T = 1600 °C.
With the measured thermal diffusivity aij(T), the thermal conductivity λij(T) was calculated: ij (T )
cm2 1700s exp s t + 86s
D0,app (t ) = 4.09
Fig. 5. Thermal diffusivity a T measurement of an orthotropic C/C material with transversal isotropy.
(11) (12) (13)
with ρC/C and λC/C as bulk density and thermal diffusivity, respectively, as well as the porosity Φ. Literature data were used for specific heat of all phases as well as density and thermal conductivity of silicon and silicon carbide. The heat Q generated due to the exothermal reaction of carbon and silicon to silicon carbide is given by:
Q = H c AR
SiC
t
s(
SiC ),
(14)
where ΔH and c AR are the enthalpy of reaction and the reactive surface per volume respectively. SiC is the increase of the thickness of the SiC t layer between carbon and silicon melt per time and s(δSiC) is a step function taking into account that the reaction stops, when either no carbon or no silicon are available in the corresponding pore channel:
c AR SiC 2 s ( SiC) = c AR SiC 1 for 2 0 else 1 for
0, Si
, 0, C/C
(15) Fig. 6. Numerical integration of Eq. (16) compared with literature data from Zhou and Singh [25] as well as Favre et al. [26].
where φ0,Si and φ0,C/C are the volume fractions of silicon and carbon before the formation of the SiC layer. ΔH was obtained from 4
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
4. Computational results and validation 4.1. Infiltration kinetics Fig. 9 shows the simulation of the propagating infiltration front with a constant velocity of 1 × 10−3 m/s based on Eq. (4). The shape and movement of the infiltration front reflect the infiltration behavior which was experimentally observed. Using Eqs. (11) and (16), the increase of mass during infiltration was calculated. The comparison with the gravimetrically measured data shows a good correlation between both curves, see Fig. 10. However, there is a small deviation at point 2 in Fig. 10, when the measured infiltration is faster than the simulation. This is attributed to a slight asymmetry in the set-up, which leads to variations of the rate of mass increase in the initial state of infiltration, when the hemispherical infiltration front arrives at the lateral surface of the cylinder. 4.2. Temperature fields Heat generation depends on the region which is already infiltrated and on the position of the infiltration front. It is calculated according to Eqs. (14)–(19). Due to the moving infiltration front, a dynamic interplay between heat generation and heat conduction within the sample and heat radiation at its surfaces occurs, which was considered in the FE model. Fig. 11 shows an example of local heat generation and temperature distribution after 18 s of silicon infiltration. As already observed experimentally, heat generation is especially strong at the infiltration front and less pronounced in the infiltrated region, see Fig. 11a, which reflects the decreasing formation rate of SiC with increasing layer thickness, see Eq. (16). Due to heat transfer to cooler regions, the temperature distribution within the sample is smoother than the heat generation rate, compare Fig. 11b. Fig. 12 shows the measured and simulated temperatures at two fixed points during infiltration: the cylinder jacket and at the core of the sample. The latter was measured in the blind hole drilled into the sample, compare Fig. 1. To check the reproducibility, two infiltration experiments were performed. It can be seen that the core temperature rises from the start temperature at 1550 °C to a maximum at 2150 °C within 20 s and decreases slowly thereafter. Even at the surface of the sample, the temperature increases to 1990 °C. There is some scatter in the experimental curves which is attributed to the coarse structure of the C/C samples. Considering this scatter, the correlation between calculated and measured temperature curves is excellent.
Fig. 7. Meshed FE model used in the simulations. – By using a symmetry plane, the furnace as well as the experimental set-up can be meshed in the half-model.
nodes of the subdomain elements, also called superelements [27,28]. By that, heat radiation could be calculated in the global model and the infiltration equation as well as the moving heat source in the submodel. Since both, the infiltration of the sample by the silicon melt and the transport of the silicon from the crucible into the sample were simulated, moving and deformed meshes were present in the submodels. Using the superelements reduced the recalculation of view factors during the simulations and saved computing time. Both, the global model, meshed with superelements and the submodel, meshed with tetrahedral elements, are shown in Fig. 8.
Fig. 8. Meshed global model and submodel of the sample. – Fig. 8a shows an enlarged section of Fig. 7. 5
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
Fig. 9. Infiltration front after 13 s, 20 s, 27 s and 34 s based on a constant infiltration velocity of 1 × 10−3 m/s. – The grid lines show the edges of the subdomains for discretization of the geometry provided by the CAD system. – The figure shows the meshed sample in Fig. 8b.
4.3. Convergence and sensitivity analysis For an estimation of the convergence properties of the FE model, computations were performed with different mesh sizes, time steps and tolerances for the iteration errros. Infiltration velocity was set to 1 × 10−3 m/s and k-value to 0.5. The maximum temperature increase due to the exothermic reaction in the sample was used for comparison. Using a mesh size of 0.8 mm, a time step of 0.1 s and a tolerance of 1 × 10−3 °C for the iteration error results in a variation of the maximum temperature of ± 2 °C around its mean value of 2248 °C, see also Fig. 4. Only the mesh size of the preform in the submodel was refined, the rest of the model remained unchanged, see Fig. 8. To verify how the model reacts to small changes in the measured infiltration rate and in the reactive surface, a sensitivity analysis was performed. The maximum temperature increased by 5 °C at a 10% increase in the infiltration rate from 1.0 mm/s to 1.1 mm/s and increased by 26 °C at a 10% increase in the reactive surface from 56.3 × 103 m2/m3 to 61.8 × 103 m2/m3. However, as the model validation shows, the uncertainty of the numerical computations are in the same range as the uncertainty resulting from the reproducibility of the measurements, see Fig. 12.
Fig. 10. Increase of mass during the infiltration of a C/C preform with liquid silicon. – The insert shows a cross section through the sample indicating the infiltration front at different times which are also marked on the curve.
6
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
Fig. 11. Heat generation and temperature field during specimen infiltration after 18 s.
in a short computer time. This is especially important, when industrial infiltration processes are optimized and heating schedules and stacks have to be varied many times. Using the present model, the LSI process for the production of ceramic break disks was successfully optimized. For that, other effects like the melting and solidifying of silicon during heating up and cooling of the charge also had to be integrated into the model. Principles of the entire model will be presented in a further work. Acknowledgement The authors gratefully acknowledge the help of P. Döppmann and J. Baber with the measurements and the financial support of the Bayerische Forschungsstiftung within the project ISE-LSI. References [1] P.J. Hofbauer, E. Rädlein, F. Raether, Fundamental mechanisms with reactive infiltration of silicon melt into carbon capillaries, Adv. Eng. Mater. 1900184 (2019) 1–11. [2] W.B. Hillig, Melt infiltration approach to ceramic matrix composites, J. Am. Ceram. Soc. 71 (2) (1988) C-96–C-99. [3] E.O. Einset, Capillary infiltration rates into porous media with applications to silcomp processing, J. Am. Ceram. Soc. 79 (2) (1996) 333–338. [4] E.O. Einset, Analysis of reactive melt infiltration in the processing of ceramics and ceramic composites, Chem. Eng. Sci. 53 (5) (1998) 1027–1039. [5] E.S. Nelson, P. Colella, Parametric Study of Reactive Melt Infiltration, Tech. Rep. E12107, National Aeronautics and Space Administration, 2000. [6] P. Sangsuwan, S.N. Tewari, J.E. Catica, M. Singh, R. Dickerson, Reactive infiltration of silicon melt through microporous amorphous carbon preforms, Metall. Mater. Trans. B 30 (5) (1999) 933–944. [7] S.P. Yushanov, J.S. Crompton, K.C. Koppenhoefer, Simulation of Manufacturing Process of Ceramic Matrix Composites, (2008). [8] E.W. Washburn, The dynamics of capillary flow, Phys. Rev. 17 (3) (1921) 273–283. [9] A. Kienzle, I. Kratschmer, Friction-Tolerant Disks Made of Fiber-Reinforced Ceramic, (2015). [10] A. Kienzle, H. Jäger, Carbon-Fiber-Reinforced Silicon Carbide: A New Brake Disk Material, (2015). [11] F.G. Raether, Current state of in situ measuring methods for the control of firing processes, J. Am. Ceram. Soc. 92 (S1) (2009) 146–152. [12] R.A. Fisher, The wave of advance of advantageous genes, Ann. Hum. Genet. 7 (4) (1937) 355–369. [13] R. Luther, Räumliche Fortpflanzung chemischer Reaktionen, Z. für Elektrochem. Angew. Phys. Chem. 12 (32) (1906) 569–600. [14] F.T.K. Meinecke, Räumliche Fortpflanzung chemischer Reaktionen, Inauguraldissertation, 1908. [15] R.K. Hobbie, B.J. Roth, Intermediate Physics for Medicine and Biology, SpringerVerlag, New York, 2007. [16] J.D. Murray, Mathematical Biology: I. An Introduction, 3rd Edition, Vol. 17 of Interdisciplinary Applied Mathematics, Springer-Verlag Berlin Heidelberg, 2002. [17] E. Brunet, B. Derrida, An exactly solvable travelling wave equation in the FisherKPP class, J. Stat. Phys. 161 (4) (2015) 801–820. [18] F. Raether, R. Hofmann, G. Müller, H.J. Sölter, A novel thermo-optical measuring system for the in situ study of sintering processes, J. Therm. Anal. Calorim. 53 (3) (1998) 717–735. [19] C. Dudescu, J. Naumann, M. Stockmann, S. Nebel, Characterisation of thermal
Fig. 12. Measured (two infiltration experiments) and calculated temperature increase at two different points during silicon infiltration into a C/C preform. – The insert in the upper right corner shows the points of temperature detection.
5. Conclusions A FE model was developed for the simulation of the infiltration kinetics and the corresponding thermal effects during infiltration of silicon melt into porous C/C preforms. A good match to the characteristics of the infiltration process obtained from in situ measurements with the measuring furnace TOM_ac was obtained. This good correlation was achieved using material data from own measurements and from literature as a basis for the simulation. Important parameters, which are specific for any type of C/C preform, are the reactive surface area, thermal conductivity and infiltration rate. Note that the latter also depends on furnace parameters like gas pressure and temperature [1]. It is recommended to measure these parameters individually for each used type of C/C material. If the reactive infiltration into other carbonaceous preforms or pyrolized wood is to be simulated, Eqs. (17) and (19) have to adapted as well. The good correlation between experiment and simulation was also caused by the use of time dependent diffusion coefficients based on these two equations. In contrast to previous models using capillary or diffusion equations for description of infiltration kinetics, a differential equation for autocatalytic processes was developed for this purpose. This new infiltration equation can be adapted to the mesh quality of FE models and thus enables 3D simulation of infiltration processes of complex components 7
Journal of the European Ceramic Society xxx (xxxx) xxx–xxx
P.J. Hofbauer, et al.
[20] [21] [22] [23] [24]
expansion coefficient of anisotropic materials by electronic speckle pattern interferometry, Strain 42 (3) (2006) 197–205. O. Zienkiewicz, R. Taylor, J. Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013. Z. Hashin, S. Shtrikman, On some variational principles in anisotropic and nonhomogeneous elasticity, J. Mech. Phys. Solids 10 (4) (1962) 335–342. M.F. Ashby, Materials Selection in Mechanical Design, 3rd ed., ButterworthHeinemann, Oxford, 2005. J.K. Farooqi, M.A. Skeikh, Finite Element Modeling of Thermal Transport in Composite Unit Cell Using Abaqus/CAE, (2005). Z. Fu, L. Schlier, N. Travitzky, P. Greil, Three-dimensional printing of sisic lattice
truss structures, Mater. Sci. Eng. A 560 (2013) 851–856. [25] H. Zhou, R.N. Singh, Kinetics model for the growth of silicon carbide by the reaction of liquid silicon with carbon, J. Am. Ceram. Soc. 78 (9) (1995) 7. [26] A. Favre, H. Fuzellier, J. Suptil, An original way to investigate the siliconizing of carbon materials, Ceram. Int. 29 (2003) 235–243. [27] J. Plaza, M. Abasolo, I. Coria, J. Aguirrebeitia, I. Fernández de Bustos, A new finite element approach for the analysis of slewing bearings in wind turbine generators using superelement techniques, Meccanica 50 (6) (2015) 1623–1633. [28] T. Panczak, S. Ring, M. Welch, A cad-based tool for fdm and fem radiation and conduction modeling, SAE Trans. 107 (1) (1998) 363–372.
8