*Article* **Characterization of Thermochemical and Thermomechanical Properties of Eyjafjallajökull Volcanic Ash Glass**

#### **Rebekah I. Webster 1, Narottam P. Bansal 2, Jonathan A. Salem 2, Elizabeth J. Opila <sup>1</sup> and Valerie L. Wiesner 3,\***


Received: 23 December 2019; Accepted: 19 January 2020; Published: 23 January 2020

**Abstract:** The properties of a volcanic ash glass obtained from the Eyjafjallajökull eruption of 2010 were studied. Crystallization experiments were carried out on bulk and powdered glass samples at temperatures between 900 and 1300 ◦C. Iron oxides, Fe3O4 and Fe2O3, and a silicate plagioclase, (Na,Ca)(Si,Al)4O8, were observed. Bulk samples remained mostly amorphous after up to 40 h at temperature. Powdered glass samples showed increased crystallinity after heat treatment compared to bulk samples. The average coefficient of thermal expansion of the glass was 7.00 <sup>×</sup> 10−<sup>6</sup> K−<sup>1</sup> over 25–720 ◦C. The Vickers hardness of the glass was 6–7 GPa and the indentation fracture toughness, 1–2 MPa <sup>√</sup>m Values for density, elastic modulus, and Poisson's ratio were 2.52 g/cm3, 75 GPa, and 0.24, respectively. The viscosity of the glass was determined experimentally and compared to three common models from the literature. The implications for the deposition of volcanic ash on hot section components of aircraft turbine engines are discussed.

**Keywords:** CMAS; EBCs; TBCs; volcanic ash; crystallization; thermal expansion; hardness; indentation fracture toughness; viscosity

#### **1. Introduction**

Particulates comprised of mainly calcium, magnesium, aluminum, and silicon oxides (CaO–MgO–Al2O3–SiO2 or CMAS) stand as a challenge in the development of coating systems for next-generation aircraft turbine engines. CMAS originates as siliceous debris, such as sand, volcanic ash, or runway dust, which can be ingested into the turbine on takeoff, landing, or during flight and deposit on hot-section components [1,2]. CMAS melts at ~1200 ◦C and can infiltrate the protective barrier coatings (thermal and environmental barrier coatings, T/EBCs) needed for engine components.

Currently, superalloy components can reach a maximum temperature of ~1200–1300 ◦C due to the presence of a thermal barrier coating (TBC) which acts as a thermally protective topcoat [3]. Silicon carbide (SiC)-based ceramic matrix composites (CMCs) are proposed to replace some conventional Ni-base superalloys as turbine engine materials due to their higher operating temperatures and lighter weight, resulting in increased efficiency [4]. EBCs are required as a topcoat for CMCs to prevent thermally grown silicon oxide (SiO2) from volatilizing in water vapor, a species resulting from the combustion of jet fuel in the engine. In addition to stability in water vapor and high-temperature environments, among other requirements [5], EBCs must be resistant to CMAS attack. The target operating temperature for CMCs approaches 1500 ◦C [6]. The temperatures experienced by both T/EBC systems are typically above CMAS melting temperatures [7].

It has been proposed in the literature that inducing crystallization of CMAS at the T/EBC-melt interface is a viable strategy in mitigating CMAS attack by limiting the extent of glass ingress [8,9]. Many studies have been conducted on different coating materials wherein the ability of the T/EBC to promote crystallization of CMAS was probed and the evolution of coating/glass phases was monitored as a function of temperature and/or time [8,10–13]. Conversely, few reports have isolated the crystallization behavior of CMAS alone. CMAS composition can vary widely, depending on its source and location [14]. As such, the intrinsic properties of these glasses may also be expected to differ. In addition to crystallization behavior, CMAS composition influences its melting temperature and viscosity, both important parameters when considering coating infiltration. Mechanical properties (i.e., coefficient of thermal expansion, hardness, toughness) of CMAS can also vary and are required to fully describe coating degradation.

Zaleski et al. used differential scanning calorimetry (DSC) to probe the melting and crystallization behavior of ternary CAS, CMAS, and Fe-containing CFAS and CMFAS synthetic deposits [15]. In general, it appeared that most of the melts under study were resistant to crystallization upon cooling at 10 ◦C/min. This behavior was most prominent in CAS compositions. The addition of Fe, however, promoted the crystallization of phases including hematite (Fe2O3), esseneite (CaFeAlSiO6), and Ca-ferrite (CaFe2O4), as indicated by x-ray diffraction (XRD). Isothermal DSC measurements on a Ca13Fe10Al18Si59 (mol %) glass showed significantly increased crystallization kinetics compared to the ternary Ca15Al20Si65 having a similar Ca:Si ratio and Al content. The presence of Fe is likely important when considering volcanic ash glasses, as bulk compositions have been previously reported with up to 16 at % Fe [16].

The objective of the current work was to characterize the intrinsic properties of a volcanic ash glass containing CMAS and Fe2O3, TiO2, Na2O, and K2O. The crystallization behavior of the volcanic ash glass was determined in bulk and powdered form. Thermal and mechanical properties including melting (*T*m) and glass transition (*T*g) temperature, coefficient of thermal expansion (CTE), hardness, indentation fracture toughness, elastic modulus (*E*), and Poisson's ratio (ν) were also measured. Glass viscosity was determined experimentally and compared to calculated values based on composition using several viscosity models described in the literature. The intrinsic properties of the volcanic ash glass were compared to those of previously investigated synthetic sand and desert sand CMAS glasses [17–19]. The significance of results with respect to T/EBC design are discussed.

#### **2. Experimental Procedure**

Ash obtained from the 2010 eruption of Iceland's Eyjafjallajökull volcano was melted at 1500 ◦C for 1 h in a platinum crucible followed by quenching in water. The composition of the resulting glass was determined by inductively coupled plasma atomic emission spectroscopy (ICP-AES) performed by a commercial test laboratory and can be found in Table 1. Table 1 also includes the composition of a synthetic sand glass [17,18] and a desert sand glass [19] for comparison.

Differential thermal analysis (DTA) was carried out using a Netzsch STA 409 (Burlington, MA, USA) on bulk and powder glass samples to determine melting and crystallization behavior. Scans were performed in flowing air from room temperature to 1500 ◦C at a ramp rate of 10 ◦C/min. A platinum pan was used to contain each sample. Prior to each run, a background scan was performed under identical conditions and was subsequently subtracted from the sample scan.

Isothermal heat treatments were carried out in a stagnant air box furnace on bulk glass pieces weighing approximately 250–350 mg. Samples were placed on a platinum foil and were either exposed for 1, 10, or 20 h at 900 and 1000 ◦C (low-temperature furnace; Neytech Vulcan 3-550; Bloomfield, NJ, USA) or for 1, 10, 20, or 40 h at 1100, 1150, 1200, and 1300 ◦C (high-temperature furnace; Carbolite HTF 18/27; Hope Valley, UK). Samples were cooled at either 10 ◦C/min or quenched in air. Some experiments were carried out using glass powder instead of bulk specimens; in these cases, bulk pieces weighing 250–350 mg were ground with mortar and pestle to a fine grit prior to exposure. After heat treatment, samples were either crushed for x-ray diffraction (XRD) phase analysis (Bruker D8 Advance; Billerica, MA, USA) or mounted and polished in cross-section for imaging by scanning electron microscopy (SEM—Phenom-World Phenom Pro; Eindhoven, The Netherlands). Mounted and polished specimens were coated with a thin layer of platinum prior to imaging. Quantitative XRD was performed by mixing 20–25 wt % α-Al2O3 (corundum) with the powdered samples after heat treatment. Powders were mixed thoroughly by grinding with a mortar and pestle. Quantitative XRD scans were run from a 2θ of 10◦ to 120◦ with a step size of 0.02◦ and a step time of 2.5 s/step. Sample peaks were referenced to the Al2O3 standard in the Whole Pattern Fitting (WPF) module of Jade 2010 software from Materials Data Inc. (MDI; Livermore, CA, USA).

**Table 1.** Composition (in mol % and wt %, determined by inductively coupled plasma atomic emission spectroscopy—ICP-AES) of the volcanic ash glass prepared in this study compared to previously investigated synthetic sand and desert sand glasses [17–19]. Estimated uncertainty for values in the range of 1.0−3.0% = ±10% of actual value, for values in the range of 3.0−10.0% = ±5% of actual value, for values in the range of 10.0−25.0% = ±2% of actual value, for values in the range of 25.0−75.0% = ±1% of actual value.


Sample discs were prepared for mechanical testing and dilatometry using a Centorr mini hot press (Nashua, NH, USA). Glass powder was loaded into a graphite die and treated at 600 ◦C and 2 ksi (13.8 MPa) for 10 min under vacuum. Specimens were cooled to room temperature at a rate of 10 ◦C/min after the applied pressure was removed. Graphite foil was used to cover the die pieces in contact with the glass during pressing, and any residual foil was removed from sample surfaces by polishing.

Glass density was calculated from the measured mass and volume of the hot-pressed glass disc. The impulse excitation method described in ASTM C1259 [20] was used to determine the elastic modulus and Poisson's ratio of the glass at room temperature. An Audio Technica ATM350 condenser microphone with an M-Audio DMP preamplifier was used to acoustically measure and amplify the natural frequency of the glass disc upon mechanical excitation in the desired mode. The acoustic signal acquisition hardware and Sound & Vibration Toolset software from National Instruments (Austin, TX, USA) were used to determine the disc's frequency, which was used to calculate the elastic modulus and Poisson's ratio.

A Struer's DuraScan (Cleveland, OH, USA) was used to apply a Vickers diamond indent onto the glass disc surface, which was polished according to ASTM C1327-08 [21]. Three indentations were made on the polished disc at each individual load of 1.96, 2.94, 4.9, and 9.8 N applied for 15 s. The diagonal (2*a*) and crack (2*c*) lengths of each indentation were measured.

A glass bar with length of 2.5 cm was evaluated using a Netzsch differential dilatometer model 402 C (Burlington, MA, USA) interfaced with a computerized data acquisition and analysis system to measure the glass transition temperature (*T*g), softening point (*T*d), and linear coefficient of thermal expansion (CTE; α). The glass bar was heated at a rate of 5 ◦C/min in the air from room temperature to 1000 ◦C.

Viscosity measurements were performed using an Orton RSV-1600 viscometer furnace setup (Westerville, OH, USA) equipped with a Brookfield (Middleboro, MA, USA) HA-DV2T viscosity measuring unit and platinum spindle. A 50 mL platinum crucible was filled with approximately 80 g of volcanic ash glass for a molten glass volume of approximately 30 mL. Viscosity measurements were performed between about 1300 and 1500 ◦C. The viscometer furnace was held at 1500 ◦C for at least

20 min to equilibrate the melt. Two methods were utilized for collecting data. In the first, the spindle was lowered into the melt and allowed to rotate continuously as the furnace cooled at 2 ◦C/min. Viscosity data was collected during cooling. In the second, the spindle was lowered into the melt and allowed to rotate. The furnace was cooled at 50 ◦C intervals and held at each temperature for 20 min. The viscosity of the glass was determined at each temperature after a stable reading was obtained.

The viscosity of the glass melt was determined by measuring the percent spindle torque required to maintain the spindle rotation at the desired speed. The viscosity of the melt (η; in centipoise) is related to the spindle torque (τ) and rotational speed (RPM) by the following equation, found in the Orton RSV-1600 instrument manual:

$$
\eta = \frac{100}{\text{RPM}} \times \text{TK} \times \text{SMC} \times \tau \tag{1}
$$

where TK and SMC are the viscometer torque constant and the spindle multiplier constant, respectively. These values were provided in the viscometer software. The viscosity of a borosilicate glass standard (NIST SRM 717a, Sigma Aldrich; St. Louis, MO, USA) was determined experimentally between 1000–1325 ◦C, confirming that the supplied constants gave accurate data in the temperature region of interest (Figure 1). The borosilicate glass standard was not run above 1325 ◦C due to volatility concerns. The calculated viscosity curve (Figure 1) was determined using viscosity constants provided by NIST.

**Figure 1.** Calculated and experimental viscosity curves for a borosilicate glass standard.

#### **3. Results**

#### *3.1. Crystallization Behavior—DTA*

DTA heating and cooling curves at a ramp rate of 10 ◦C/min are given in Figure 2 for the bulk and powdered glass. A list of observed DTA thermal events, described for bulk and powdered glasses, is given in Table 2. On heating, the bulk glass exhibited exothermic peaks at approximately 900 (Figure 2—1B) and 1100 ◦C (2B), which corresponded to crystallization. Each peak is likely attributed to the crystallization of a distinct phase. There appeared to be a melt endotherm with a peak temperature of approximately 1350 ◦C; however, the peak itself was very broad. A dip in the baseline was evident at about 1050 ◦C (3B), suggesting that this glass started to melt around 1050 ◦C and might not have a distinct melting temperature. Indeed, in separate box furnace experiments glass specimens were observed to bead by 1100 ◦C and wet the platinum holder substrate by 1150 ◦C. Bulk glass specimens after heat treatment at 1100 and 1150 ◦C are shown in Figure 3. Bulk glass heat treated at 1300 ◦C for 1 h and air quenched was almost entirely amorphous by XRD analysis, suggesting that by this temperature, the glass was nearly completely molten. On cooling, an exothermic peak was also observed around

900 ◦C (4B) suggesting that a similar phase that crystallized around 900 ◦C on heating (1B) might have also evolved upon cooling at 10 ◦C/min.

**Figure 2.** Bulk (bold line) and powder (dashed line) volcanic ash glass differential thermal analysis (DTA) curves with a ramp rate of 10 ◦C/min. Labels 1B, 2B, 1P, 2P, and 3P correspond to crystallization in the bulk (B) and powder (P) samples on heating. Labels 3B and 4P show the onset of melting. Labels 4B and 5P indicate solidification of molten glass upon cooling.

**Table 2.** List of DTA thermal events presented in Figure 2 for the volcanic ash glass.


**Figure 3.** Bulk volcanic ash glass samples before heat treatment (**a**) and after heat treatment at 1100 (**b**) and 1150 ◦C (**c**). The solidified glass samples are the dark droplets on Pt foil.

The heating curve for powdered glass was quite different from that seen for the bulk, as shown in Figure 2. There was a broad exothermic hump that centered around 900 ◦C. This peak was not as distinct as for the bulk glass and had a comparatively lower onset temperature (1P). In addition, instead of having a second exothermic peak around 1100 ◦C, there were two very large peaks at approximately 1225 (2P) and 1275 ◦C (3P). It is somewhat difficult to distinguish background information from actual thermal events for this curve, however, the melting behavior appeared similar to that for the bulk—the onset (4P) and peak melting temperatures being shifted to slightly lower values. On cooling, there was an exothermic peak at approximately 900 ◦C (5P), echoing that observed for the bulk sample.

#### *3.2. Crystallization Behavior—Box Furnace Experiments*

XRD plots for bulk glass pieces treated for 1 h at 900–1200 ◦C can be seen in Figure 4. All scans were collected for samples after being heated and cooled at 10 ◦C/min. XRD patterns for specimens that were placed in and removed from the furnace at temperature (i.e., quenched in air) did not show differences in phase or amount of phase present when compared to samples that were ramped to and from temperature. Consequently, XRD results for quenched specimens were excluded from the rest of this report. In Figure 4, the XRD scan for glass not exposed to temperature is also shown (labeled "Control") for comparison. Only two phases were discerned—Fe3O4 (magnetite) and Fe2O3 (hematite). After heat treatment at 900 and 1000 ◦C, broad, shallow peaks for magnetite were present, but the glass appeared to have remained mostly amorphous. At 1100, 1150, and 1200 ◦C, the peak for magnetite appeared larger and sharper. At 1150 and 1200 ◦C, the 100% intensity peak for hematite began to emerge.

**Figure 4.** XRD plots for bulk volcanic ash glass heat treated for 1 h at 900, 1000, 1100, 1150, and 1200 ◦C (ramp rate 10 ◦C/min). The "Control" spectrum is for volcanic ash glass not exposed to temperature.

As the amount of time held at temperature was increased for 1100–1200 ◦C treatments, the magnetite was converted to hematite (Table 3). The transition from magnetite to hematite occurred more rapidly as the temperature increased, such that there was a higher relative amount of hematite after 10 h at 1200 ◦C compared to after 10 h at 1100 ◦C (Table 3). This is likely due to the enhanced diffusion of oxygen through the bulk glass at higher temperatures [22]. The reaction of Fe3O4 to Fe2O3 is given by [23]

$$4\text{Fe}\_3\text{O}\_4 + \text{O}\_2(\text{g}) \leftrightarrow 6\text{Fe}\_2\text{O}\_3\tag{2}$$

The evolution of the microstructure of samples held at 1150 ◦C can be seen in Figure 5. Energy dispersive spectroscopy (EDS) was unable to differentiate Fe3O4 from Fe2O3. The bright (by backscattered electron imaging, BSE) Fe-based precipitates appear as both faceted particulates and rod-like structures. Faceted particles are mostly submicron in size after 1 h at temperature but grow to a maximum diameter of approximately 10 μm with increasing time. In some areas, the surface

of the glass had a thin layer or "shell" of Fe-based particulates (somewhat visible in Figure 5c,d, more apparent in Figure 6).


**Table 3.** The phase composition of bulk volcanic ash glass samples after heat treatment using quantitative XRD analysis.

**Figure 5.** Cross-section backscattered electron imaging (BSE) images of bulk volcanic ash glass after exposure at 1150 ◦C for 1 h (**a**), 10 h (**b**), 20 h (**c**), and 40 h (**d**).

**Figure 6.** Cross-section BSE images for bulk volcanic ash glass samples treated at 1100 (**a**,**c**) and 1150 ◦C (**b**,**d**) for 10 h.

At 1100 ◦C, a third phase, likely a solid solution of anorthite (CaAl2Si2O8) and albite (NaAlSi3O8), was detected after holding for 10 h (see Figure 6). This phase (hereto referred to as plagioclase) may also be present in small amounts at 1000 ◦C.

SEM BSE cross-section images are presented in Figure 6 for samples held at 1100 and 1150 ◦C for 10 h. The results given by SEM supported those of XRD. After 10 h at 1100 ◦C, there appeared to be two phases growing at the surface of the glass (Figure 6a). EDS confirmed a Fe-rich phase (lighter contrast by BSE) and a Si-rich phase (darker by BSE). Iron oxides (likely mostly Fe3O4 based on XRD) were also present throughout the bulk of the sample (Figure 6c) but were generally smaller than those in the 1150 ◦C sample (Figures 5b and 6d). A thin layer of concentrated Fe precipitates can be seen at the surface of both the 1100 and 1150 ◦C (Figure 6c,d) samples. This "shell" is likely the Fe2O3 phase, as it has been previously observed in magnetite iron ore pellets oxidized at 800–1200 ◦C [24]. The "shell" was not always continuous across the surface of samples, which may be due to the removal of material when separating the bulk glass from its platinum substrate. As expected from XRD, no plagioclase phase was observed in the 1150 ◦C sample (Figure 6b).

When powdered samples were heated in the box furnace at 1100 and 1150 ◦C, the prominent crystallized phases were Fe2O3 and plagioclase (Figure 7). The amounts of each phase are given in Table 4. It is expected that the increased amount of crystallized glass compared to bulk samples is due to the increased surface area of the powder.

FactSage free energy minimization calculations were performed for both the volcanic ash glass and the synthetic sand glass compositions (Table 1), using the Equilibrium Module and FToxide database, including SLAGA [25], to determine the expected phase distributions at a given temperature based on thermodynamics alone. Respective amounts of the different CMAS oxides (in mol %, Table 1) were used as FactSage inputs. Calculations were performed at 1 atm. As was observed experimentally, FactSage indicated that Fe2O3 and albite (NaAlSi3O8) and anorthite (CaAl2Si2O8) were among the expected phase formations. Between 900–1150 ◦C, SiO2, MgSiO3, KAlSi3O8, CaSiTiO5, and CaMgSi2O6 (diopside) were among other phases also predicted to appear. Albite was not seen above 1050 ◦C and Fe2O3 was the only expected solid-phase above ~1200 ◦C. The amount of the observed phases as a function of temperature is given in Figure 8. A Scheil–Gulliver cooling profile was also performed, with a start and stop temperature of 1300 and 900 ◦C, respectively, and a step size of 25 ◦C. This method subtracts precipitated phase constituents from the overall glass composition as they form. At and above 1200 ◦C, only Fe2O3 is expected. Below 1200 ◦C, anorthite forms; at ≤1100 ◦C, MgSiO3, diopside, and CaSiTiO5 appear in addition to Fe2O3 and anorthite. Albite is only discerned at ≤1000 ◦C.

**Figure 7.** XRD spectra for bulk (solid line) and powder (dashed line) volcanic ash glass samples exposed to 1100 ◦C for 1 h.

**Figure 8.** Expected phase formations in the volcanic ash glass as a function of temperature, calculated using the FactSage Equilibrium module.

The same processes were performed for the synthetic sand glass. Wollastonite (CaSiO3) and diopside were among the expected phase formations between 900–1300 ◦C. Other predicted phases (between 900–1000 ◦C) included albite, SiO2, and Na2Ca3Si6O16. Figure 9 plots the expected phase formations in the synthetic sand glass as a function of temperature.


**Table 4.** The phase composition of powder volcanic ash glass samples after heat treatment using quantitative XRD analysis.

**Figure 9.** Expected phase formations in the synthetic sand glass as a function of temperature, calculated using the FactSage Equilibrium module.

#### *3.3. Thermal Behavior—Dilatometry*

Figure 10 shows the dilatometric thermal expansion curve for the hot-pressed volcanic ash glass bar heated at 5 ◦C/min from RT to 1000 ◦C. Dilatometric values for the volcanic ash glass, synthetic sand glass, and desert sand glass are found in Table 5. The *T*<sup>g</sup> and *T*<sup>d</sup> of the volcanic ash glass were observed to be 741 and 844 ◦C, respectively. Comparatively, the synthetic and desert sand glasses referred to in Table 1 had reported *T*<sup>g</sup> values of 694 and 706 ◦C, respectively, and *T*<sup>d</sup> values of 751 and 764 ◦C, respectively. The average linear CTE for the bulk volcanic ash glass was 7.00 <sup>×</sup> 10−<sup>6</sup> K−<sup>1</sup> between 25–720 ◦C. This value was lower than the reported ~9–10 <sup>×</sup> 10−<sup>6</sup> K−<sup>1</sup> for the synthetic and desert sand glasses over the same temperature range.

**Figure 10.** Dilatometric thermal expansion curve on heating for the volcanic ash glass at a heating rate of 5 ◦C/min in air. Glass transition temperature (*T*g) and glass softening temperature (*T*d) are indicated.



\* Values obtained from [17,18]; ± Values obtained from [19].

#### *3.4. Density and Mechanical Properties*

The bulk density, Elastic Modulus (*E*), Poisson's ratio (ν), Vickers hardness (*H*V), and indentation fracture toughness (evaluated by three different equations) for the volcanic ash glass are reported in Tables 5 and 6. Density was found to be 2.52 g cm−<sup>3</sup> and E and ν to be 75 GPa and 0.24, respectively. These values were comparable to those reported for the synthetic and desert sand glasses, though E for the synthetic and desert sand glasses was approximately 84 and 92 GPa, respectively, indicating that the volcanic ash glass was less stiff.

Vickers microhardness was determined at each load using the following equation [21]:

$$H\_{\rm V} = 1.8544 \left[ \frac{P}{\left(2a\right)^2} \right] \tag{3}$$

where *P* was the applied load and 2*a* was the indentation diagonal length. Vickers hardness values determined at loads of 1.96, 2.94, 4.9, and 9.8 N are reported in Table 6.

Indentation fracture toughness (*K*C) of the volcanic ash glass was calculated from indentation load and crack and diagonal lengths using relationships found in the literature. Miyoshi et al. [26] used the following:

$$K\_{\mathbb{C}} = 0.026 a E^{0.5} P^{0.5} \text{ c}^{-1.5} \tag{4}$$

where *E* was the Young's Modulus, *P* was the indentation load, and *a* and *c* were half the indent length and crack length, respectively. Marshall and Evans [27] evaluated indentation fracture toughness using the following relation:

$$K\_{\mathbb{C}} = 0.036 E^{0.4} \, P^{0.6} \, a^{0.8} \, c^{-1.5} \tag{5}$$

Anstis et al. [28] presented:

$$K\_{\rm C} = \left(0.016 \pm 0.004\right) P \left(\frac{E}{H\_{\rm V}}\right)^{0.5} c^{-1.5} \tag{6}$$

to determine the indentation fracture toughness of glass specimens. *K*<sup>C</sup> values calculated using Equations (4)–(6) can be found in Table 6 for a range of indentation loads. The loads were chosen to prevent spalling while giving a reasonably long crack (*c* > 2*a*). The hardness values for the volcanic ash glass are comparable to those of the synthetic and desert sand glasses, but the volcanic ash glass is significantly tougher (~1–2 MPa <sup>√</sup>m). The toughness of the synthetic and desert sand glasses is akin to that of many glasses, however, that of the volcanic ash glass more closely matches a glass-ceramic [29]. Hot pressing the glass at 600 ◦C resulted in the crystallization of some Fe3O4 precipitates in the sample (data not shown), which may explain the increased toughness.

**Table 6.** Vickers hardness and indentation fracture toughness at various indent loads for the volcanic ash glass and synthetic sand and desert sand glasses.


#### *3.5. Viscosity*

Experimental viscosity values for the volcanic ash glass, obtained by both continuous cooling and isothermal holds, are plotted in Figure 11 as a function of temperature. EDS analysis of the glass following viscosity testing suggested that the composition of the glass did not change during testing.

The experimental viscosities range from ~15–105 Pa·s (log values ~1.25–2) with a decreasing temperature. These values are considerably lower than expected from model predictions. Glass viscosity was predicted using three different viscosity models (Giordano et al., Fluegel, and the FactSage Viscosity Module with Melt Database). The Giordano et al. model [30] calculates the non-Arrhenian temperature dependence of viscosity for naturally occurring silicate melts by connecting experimentally obtained viscosity profiles with VFT (Vogel–Fulcher–Tamman) constants. The VFT equation is given by [31–33]

$$
\log \eta = \mathbf{A} + \frac{\mathbf{B}}{T - \mathbf{C}} \tag{7}
$$

where A, B, and C are constants. The Fluegel model [34] also correlates VFT constants to experimental viscosity data, but for commercial silicate-based glasses. On the other hand, the FactSage method [25] of calculation is non-empirical, instead relating viscosity to melt structure. The Modified Quasichemical Model is utilized, along with thermodynamic values, to calculate viscosity for a given composition. The viscosity, based off of the Giordano et al., Fluegel, and FactSage models, is given as a function of temperature for the volcanic ash glass in Figure 11, alongside experimental values. It is assumed that the glass is fully molten at temperatures 1325–1500 ◦C based on the results discussed previously.

**Figure 11.** Viscosity curves for the volcanic ash glass used in this study based on experimental and calculated values. Calculated values were determined using the FactSage Viscosity Module and Melt Database (black squares), the Giordano et al. model (blue circles), and the Fluegel model (red triangles). Experimental data are given in green. The continuous curve is for the experimental method in which the glass was cooled at a constant 2 ◦C/min while the separate diamond-shaped data points are for the experimental method in which the glass was equilibrated at 50 ◦C intervals.

It can be seen from Figure 11 that the FactSage and Fluegel models give values that are in good agreement, while the Giordano model viscosity is greater by up to about one order of magnitude. Wiesner et al. performed similar calculations on the synthetic sand glass composition given in Table 1 and also saw that the FactSage and Fluegel models showed reasonable agreement. Experimental measurements confirmed that the FactSage and Fluegel models were more accurate than the Giordano et al. model for that particular composition [35,36]. This result did not translate to the volcanic ash glass, which had an experimental viscosity lower than the FactSage and Fluegel models by about half an order of magnitude.

It is important to note that the input composition for the Giordano model is slightly different from that reported in Table 1 due to model constraints—notably, that the input was FeO instead of Fe2O3. The composition of the glass for this model was 5.1CaO-2.5MgO-15.2Al2O3-59.9SiO2-9.3FeO-1.7TiO2-4.6Na2O-1.7K2O (wt %). The Fe2O3 reported in Table 1 was converted to FeO (assuming 2 mol FeO for 1 mol Fe2O3 based on Fe alone) and trace oxides were not taken into account due to their total weight being unknown.

Wiesner and Bansal [18] estimated the amount of time (t) it would take for their synthetic sand CMAS to infiltrate a 200 μm TBC, barring CMAS crystallization and/or other thermochemical interactions with the coating, using the following

$$t \sim \left[\frac{k\_t}{8D\_c} \left(\frac{1-\omega}{\omega}\right)^2 L^2\right] \frac{\eta}{\sigma\_{\rm LV}}\tag{8}$$

where a tortuosity (*k*t) value of ~3 was considered for the coating [5], the capillary diameter (*D*c) was 1 μm [37], ω ≈ 0.1 was the pore fraction open to flow [37], and surface tension (σLV) was assumed to be ~0.4 J m<sup>−</sup>2, calculated using an approach for silicate glass melts at 1400 ◦C [38]. Assuming these values for the current system and using the viscosities calculated at 1400 ◦C by each model, the infiltration time was determined to be ~29 min for the Giordano model and ~7–8 min for both the FactSage and Fluegel models. These values are compared to the calculated infiltration times for the synthetic sand glass in Table 7. Given that the experimental viscosity for the volcanic ash glass is about half an order of magnitude less than that calculated by the FactSage and Fluegel approaches, the expected time to complete infiltration at 1400 ◦C is actually around 2 min. Experimental viscosities and infiltration times at ~1400 ◦C are also given in Table 7 for the volcanic ash glass and the synthetic sand glass reported previously [35].

**Table 7.** Viscosity (η) values (in Pa·s) as a function of temperature, model, and glass composition. The predicted time to complete the infiltration of a 200 μm thick thermal barrier coating (TBC) (dependent on viscosity) at 1400 ◦C is given for synthetic sand [35,36] and volcanic ash glasses. The expected infiltration time based on experimental viscosity data (using isothermal holds) is also given.


#### **4. Discussion**

The crystallized phases observed in this study are similar to those reported elsewhere. Mechnich et al. previously investigated a synthetic volcanic ash glass, prepared by sol-gel, with a composition approximating that of an actual Eyjafjallajökull deposit [39]. Calcination of the base glass powder at 800 ◦C resulted in the formation of Fe2O3, as well as some pseudobrookite (Fe2TiO5) and aegirine (NaFeSi2O6), within the glass. After heat treatment of a powdered mixture of the volcanic ash glass with TBC yttria-stabilized zirconia (YSZ) at temperatures of 900–1200 ◦C, the hematite phase was still discerned and plagioclase was evident ≥1000 ◦C. An electron beam-physical vapor deposited (EB-PVD) YSZ coating exposed to the glass for 1 h at 1200 ◦C also showed the formation of plagioclase and an Fe-rich oxide at its surface. Jang et al. characterized an as-received volcanic ash deposit from the Eyjafjallajökull 2010 eruption and determined the presence of anorthite/albite (plagioclase), augite (Ca(Mg,Fe,Al)(Si,Al)2O6), and analcime (Na(Al,Si2O6)H2O) [40]. Interestingly, they did not report any Fe-based phases in the as-received sample despite it being composed of nearly 10 wt % Fe2O3. The heat treatment history of their sample is unknown.

The volcanic ash glass showed considerably different crystallization behavior compared to a previously investigated synthetic sand glass (Table 1). Wiesner and Bansal [17] saw the formation of wollastonite (CaSiO3) and aluminum diopside (Ca(Mg,Al)(Si,Al)2O6) in bulk CMAS samples at temperatures as low as 925 ◦C. For the glass in this study, while DTA indicated crystallization at 900 ◦C, and FactSage calculations suggested the same, furnace tests revealed that only a very small amount of Fe3O4 formed at this temperature. The synthetic sand CMAS glass was fully crystalline after 20 h at 925 ◦C or 5 h at 960 ◦C. Based on the results of this study, it is clear that bulk volcanic ash glass pieces would not fully crystallize at these temperature/time scales. It is interesting to note that the powderized glass showed significantly increased crystallization kinetics (Table 4) compared to the bulk (Table 3) at 1100 and 1150 ◦C. It is likely that the glass powder would also show increased crystallization at 900, 1000, and 1200 ◦C compared to the bulk. It is possible, and apparently more desirable in terms of crystallization, that a protective coating (T/EBC) in service will come into contact with fine particles instead of bulk pieces.

Volcanic ash glass of this composition in contact with a T/EBC at temperatures nearing 1300 ◦C will likely soften/melt and penetrate the coating via defects such as open channels (in a TBC) or grain boundaries and pores/cracks (in an air plasma spray (APS)-deposited EBC). Barring any chemical interaction between the glass/coating, it is unlikely that infiltration will be halted by intrinsic glass crystallization, even with the assumption that there is a thermal gradient in the coating. A 250 μm thick EB-PVD YSZ TBC was nearly completely infiltrated after only 1 h at 1200 ◦C in contact with an artificial volcanic ash glass [39]. Another study by Mechnich et al. investigated an alternative TBC material, Gd2Zr2O7 (GZO), in contact with the artificial volcanic ash and saw slowed infiltration (50 μm after 1 h at 1200 ◦C) compared to YSZ [41]. This slowing was likely due to the dissolution of the coating in the glass to form newly crystallized reaction products, including an oxyapatite phase, Ca2Gd8(SiO4)6O2. In the study by Jang et al., the authors investigated the reaction between volcanic ash and a dense EBC material, Yb2SiO5, at 1400 ◦C [40]. Infiltration was halted by the formation of a thin layer of Yb2Si2O7. Yb2SiO5, in contact with a CMAS composition more closely approximating that of the synthetic sand glass (Table 1; increased CaO content, decreased SiO2 content), formed a much thicker reaction layer containing Yb apatite, Ca2Yb8(SiO4)6O2.

TBC materials YSZ and GZO have higher CTEs than the bulk volcanic ash glass, with values of ~11–12 [42,43]. If a TBC is infiltrated and then cooled, strain from a thermal mismatch between the coating and glass can lead to cracking and delamination. EBC materials, rare earth monosilicates (Y2SiO5, Yb2SiO5), have CTEs ~6–7.5 [6,44,45], which are similar to that of the bulk volcanic ash glass. These materials, however, will have a relatively high thermal mismatch with an underlying SiC-based CMC. SiC has a CTE of about 4.5–5.5 <sup>×</sup> 10−<sup>6</sup> ◦C−<sup>6</sup> [6]. Other proposed EBC materials, rare earth (RE) disilicates (Y2Si2O7, Yb2Si2O7), have CTEs that more closely match that of SiC, ranging ~3.5–5 <sup>×</sup> <sup>10</sup>−<sup>6</sup> ◦C−<sup>6</sup> [46,47]. The CTE of the bulk glass is greater than that of the CMC/disilicate system and the CTEs of the intrinsically crystallized phases in the glass are considerably so. CTE values (to 1000 ◦C) for the crystalline products of Fe3O4 and Fe2O3 are around 9–15 and 9–12 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>K</sup>−<sup>1</sup> [48], respectively, and 3 and 1.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>K</sup>−<sup>1</sup> for albite and anorthite, respectively [49]. Similar to that expected in a TBC system, CTE mismatch between the glass/crystallized glass and an EBC will likely result in spallation at the coating/glass interface.

In addition to CTE, mechanical properties such as elastic modulus, hardness, and toughness of the coated CMC can be compromised by CMAS attack. The bulk elastic moduli of TBC materials YSZ and GZO are on the order of 200–250 GPa [42,50–52]. However, the typical in-plane modulus for an EB-PVD TBC is about 30 GPa [13]. Upon infiltration, the coating is stiffened, with an E reaching at least that of the glass (75 GPa), and the in-plane compliance of the TBC is degraded [13]. For current EBCs of γ-Y2Si2O7 and β-Yb2Si2O7, E is around 155 and 168 GPa [53,54], respectively, while Y2SiO5 and Yb2SiO5 have values of 123 and 149 GPa, respectively [55,56]. The effect of changes in the elastic modulus is likely not as important in EBCs because they are relatively dense compared to TBCs.

Hardness and toughness values for YSZ and GZO are between ~12–14 and 1–2, respectively [50,52]. The hardness of the volcanic ash glass is much lower than these TBC materials, while the toughness is similar. For the Yb- and Y-based silicates, *H*<sup>V</sup> and *K*<sup>C</sup> values are 6.4 and 2.3, respectively, for Yb2SiO5 [56], 7.3 and 2.8 for Yb2Si2O7 [54], 5.3 and 2.2 for Y2SiO5 [55], and 6.2 and 2.1 for Y2Si2O7 [53]. Hardness values reported here for the solidified volcanic ash are on the order of those for the RE silicates; more concerning is the apparent decrease in toughness, which could lead to diminished resistance of the coating to fracture due to CMAS infiltration. It is noted that the volcanic ash glass has a greater toughness than that reported for the synthetic and desert sand glasses.

The viscosity of the volcanic ash glass determined experimentally in this study was much lower than predicted by model calculations. Wiesner et al. reported that the FactSage and Fluegel models better represented experimental viscosity values for the synthetic sand glass composition (Table 1) than the Giordano et al. model [35,36]. While the FactSage and Fluegel model predictions more closely matched experimental data for the volcanic ash glass, compared to the Giordano et al. model, they were still higher by about half an order of magnitude. Preliminary Fourier-transform infrared spectroscopy (FTIR) measurements on the volcanic ash glass suggest that H2O was incorporated into its structure. The presence of water in glass can significantly lower glass viscosity [57], which is a

possible explanation for the lower measured viscosity values. It is important to consider the effect of glass bonding with water when modeling CMAS viscosity, especially when considering an actual engine environment wherein high water vapor partial pressures are expected. It is also possible that none of the models investigated are adequate to describe viscosity for this particular composition. Additionally, the valence state of Fe in the glass is uncertain and can vary as a function of composition, temperature and pO2, and its role in the glass structure. The valence state's effect on glass structure is related to the overall viscosity of the melt. The available iron oxide inputs for the Fluegel and Giordano et al. models are limited to Fe2O3 and FeO, respectively. The FactSage model can incorporate both Fe2O3 and FeO as input compositions. When substituting FeO for Fe2O3 in the volcanic ash glass composition given by Table 1, the predicted viscosity was lower (Table 8), but not to the degree given by experimental data. It is proposed that all Fe2<sup>+</sup> acts as a glass network modifier while Fe3<sup>+</sup> acts as a network former when Fe3<sup>+</sup>:Fe2<sup>+</sup> > 1 [58], explaining the observed changes in calculated viscosity.

**Table 8.** Model FactSage viscosity (η) values using Fe2O3 or FeO as inputs. Using FeO in place of Fe2O3, the composition was the same as that for the Giordano et al. model (5.1CaO–2.5MgO–15.2Al2O3–59.9SiO2–9.3FeO–1.7TiO2–4.6Na2O–1.7K2O wt %).


The infiltration times calculated for the synthetic sand glass and the volcanic ash glass determined from experimental viscosity values are not very different. The expected time to complete infiltration (determined using Equation (7) and experimental viscosity values; Table 7) is only slightly longer for coatings exposed to the volcanic ash glass (1.6 min at 1410 ◦C) compared to the synthetic sand glass (0.41 min at 1411 ◦C). This small difference in infiltration time is likely due to a slightly higher viscosity of the volcanic ash glass resulting in part from a greater SiO2 and lesser CaO content.

Equation (7) has limitations in its ability to predict the infiltration rate in T/EBCs. The competing effect of crystallization in the determination of infiltration time should be considered for a TBC such as GZO. This coating material dissolves to some degree in contact with CMAS and contributes metal ions to the melt, inducing crystallization of oxyapatite. In addition, this equation is not well suited to predict infiltration in APS-deposited EBCs. APS EBCs are nominally dense and do not contain open channels for CMAS to penetrate (although grain boundaries and pores are susceptible). However, regardless of type, a coating material that reacts quickly with the volcanic ash glass to form newly crystallized phases will likely slow penetration. The low viscosities (Table 7) and sluggish crystallization reported for the volcanic ash glass studied here can lead to fast infiltration via coating channels and defects.

#### **5. Summary and Conclusions**

The intrinsic thermal and mechanical properties of a volcanic ash glass were studied using a variety of characterization techniques. It was determined that the glass had a low propensity to crystallize in bulk form, with ≤20 wt % transforming to Fe3O4, Fe2O3, and/or plagioclase after up to 40 h at temperatures of 900–1200 ◦C. This was in contrast to a previously investigated synthetic sand glass, which was more easily able to crystallize and formed wollastonite (CaSiO3) and diopside (Ca(Mg,Al)(Si,Al)2O6). The ability of the volcanic ash glass to crystallize was improved when exposed in powder form.

The melting temperature (*T*m), glass transition temperature (*T*g), and glass softening temperature (*T*d) of the volcanic ash glass were greater than those of sand glass compositions. There was a broad melting range for the volcanic ash glass, determined by DTA; after heating to 1300 ◦C at 10 ◦C/min and quenching in air, the glass was mostly amorphous by XRD analysis. The coefficient of thermal expansion (CTE) and Young's modulus (*E*) of the volcanic ash glass were lower than those reported for sand glasses. Hardness values were similar, but the indentation fracture toughness of the volcanic ash glass was about twice that of the sand glasses, likely due to the presence of some Fe3O4 crystallites in the tested sample.

The viscosity of the volcanic ash glass was lower than expected from model predictions. This is unlike that reported for the synthetic sand glass composition, which showed relatively good agreement with Fluegel and FactSage viscosity models. The discrepancies in this study could be due to the difference in the glass composition, or possibly the incorporation of water in the volcanic ash glass structure. Experimental viscosities for the volcanic ash glass were higher than for the synthetic sand glass, in part likely due to a greater SiO2 and lesser CaO content.

In conclusion, it has been shown that the chemical and mechanical properties of the Eyjafjallajökull volcanic ash CMAS glass can vary significantly from sand-based CMAS glass compositions. In assessing the potential of new coating (T/EBCs) materials for use over a wide range of operating conditions, the variation of CMAS properties with composition must be understood.

**Author Contributions:** Conceptualization, V.L.W.; methodology, V.L.W. and R.I.W.; validation, V.L.W., J.A.S. and N.P.B.; formal analysis, R.I.W.; investigation, R.I.W.; resources, V.L.W., J.A.S., N.P.B. and E.J.O.; data curation, V.L.W., J.A.S., N.P.B. and R.I.W.; writing—original draft preparation, R.I.W.; writing—review and editing, V.L.W., J.A.S., N.P.B. and E.J.O.; visualization, V.L.W.; supervision, V.L.W., J.A.S. and N.P.B.; project administration, V.L.W.; funding acquisition, V.L.W. and E.J.O. All authors have read and agree to the published version of the manuscript.

**Funding:** This research was supported by NASA's Transformative Tools and Technologies (TTT) Project within the Transformative Aeronautics Concept Program (TCAP) and the Lewis' Educational and Research Collaborative Internship Project (LERCIP) at NASA Glenn Research Center.

**Acknowledgments:** The authors are grateful to Gustavo Costa, Bryan Harder, Nathan Jacobson, Dereck Johnson, and Richard Rogers for training and helpful discussion, and Marie Kløve Keiding, Magnus Gudmundsson, and Tinna Jónsdóttir for assistance in obtaining volcanic ash from the Eyjafjallajökull eruption used in this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Flow Kinetics of Molten Silicates through Thermal Barrier Coating: A Numerical Study**

#### **Mohammad Rizviul Kabir 1,\*, Anil Kumar Sirigiri 2, Ravisankar Naraparaju <sup>2</sup> and Uwe Schulz <sup>2</sup>**


Received: 16 April 2019; Accepted: 21 May 2019; Published: 23 May 2019

**Abstract:** Infiltration of molten calcium–magnesium–alumina–silicates (CMAS) through thermal barrier coatings (TBCs) causes structural degradation of TBC layers. The infiltration kinetics can be altered by careful tailoring of the electron beam physical vapor deposition (EB-PVD) microstructure such as feather arm lengths and inter-columnar gaps, etc. Morphology of the feathery columns and their inherent porosities directly influences the infiltration kinetics of molten CMAS. To understand the influence of columnar morphology on the kinetics of the CAMS flow, a finite element based parametric model was developed for describing a variety of EB-PVD top coat microstructures. A detailed numerical study was performed considering fluid-solid interactions (FSI) between the CMAS and TBC top coat (TC). The CMAS flow characteristics through these microstructures were assessed quantitatively and qualitatively. Finally, correlations between the morphological parameters and CMAS flow kinetics were established. It was shown that the rate of longitudinal and lateral infiltration could be minimized by reducing the gap between columns and increasing the length of the feather arms. The results also show that the microstructures with long feather arms having a lower lateral inclination decrease the CMAS infiltration rate, therefore, reduce the CMAS infiltration depth. The analyses allow the identification of key morphological features that are important for mitigating the CMAS infiltration.

**Keywords:** TBCs; CMAS; infiltration; microstructure; modeling; finite element

#### **1. Introduction**

Hot section engine components of gas turbines are insulated with thermal barrier coatings (TBCs) to protect them against excessive heat. The TBC layers fabricated by electron beam physical vapor deposition (EB-PVD) show superior structural stability and long-life compared to the TBCs produced by the atmospheric plasma spray (APS) method due to their feathery columnar microstructure [1–6]. However, state-of-the art standard 7YSZ (7 wt.% yttria stabilized zirconia) coatings are vulnerable against environmental attacks, i.e., sand, volcanic ash, etc., that comes from dusty air intake during engine operation. These silicate-based particles melt under high engine temperature forming calcium–magnesium–alumina–silicates (CMAS), deposit on coated parts, and infiltrates through the porous channels of the TBC layers [7,8]. The contamination of CMAS eventually degrades the structural properties of TBCs, causing spallation failure under thermal cycles [9–13].

In past years, considerable attempts were made on the improvement of TBC systems aiming at reducing the infiltration and contamination by molten silicates [14–19]. Several methodologies were suggested to enforce the mitigation mechanisms against CMAS, which fall under two broad categories. In the first approach, new TBC compositions were proposed that react with CMAS, and recrystallize

into new stable phases, which restrict further CMAS infiltration [20–23]. A long-term resistance against the CMAS flow cannot be guaranteed due to the continuous reaction and recession of the TBC material. In the second approach, microstructures of the top coat were custom-tailored for distributed porosities, through which the deposition of CMAS occurs locally over the top coat. After solidification of CMAS inside these pores, further infiltration was restricted, thus, the total infiltration depth was reduced.

The latter approach was systematically analyzed in some previous experimental works [24–26], where the CMAS infiltration behavior was compared for different top coat microstructures. It was found out that by using the optimized columnar microstructure the CMAS infiltration kinetics could be changed, where infiltration mainly proceed due to the capillary forces. The results suggest that custom-tailored well-defined porosities of the TBC microstructure may potentially impede the infiltration of CMAS. However, the knowledge provided by the previous works is not sufficient for identifying such an optimized morphology of TBCs. A systematic analysis to understand the microstructure dependent infiltration behavior is highly required. This goal can be achieved by conducting additional experiments with a number of columnar structures obtained via EB-PVD, which involves many trial-and-error modifications of the top coat morphology. Additionally, advanced numerical analyses can be suggested to predict the microstructural influence on infiltration kinetics, which will obviously reduce the experimental efforts, and additionally, will provide an in-depth understanding of the infiltration processes.

In past decades, the numerical investigations of TBCs were performed focusing mostly on the problems of structural stability, life-cycle reduction, and failure [27–30]. The effects of deposited CMAS were investigated in terms of structural property degradation. The CMAS infiltration depth was assumed based on the experimental observation. As the infiltration depth is dependent on the morphological factors of the TBC, it is desirable to develop models that firstly describe the infiltration kinetics, and secondly, predict the infiltration depth considering the morphological obstacles.

In one of the recent works by the authors [24], the infiltration depth of CMAS through different porous EB-PVD top coats was modeled analytically, where the longitudinal porous channels were assumed as open or concentric pipes, through which CMAS infiltrated under the action of capillary forces. Microstructure obstacles imposed on the CMAS flow was considered by a parameter "tortuosity factor". The model predicted the microstructure dependency of the infiltration depth by fitting tortuosity factors for different morphologies. A good approximation of the infiltration time compared to depth was obtained which closely predicts the experimental observation [31]. However, interest on direct microstructure modeling incorporating the flow behavior is growing tremendously, as this approach allows describing realistic models based on actual microstructural statistics, and also allows fluid structure interactions during liquid flow. Thus, no "tortuosity" assumptions are required to be imposed on the moving liquid.

Microstructure based models also give an opportunity to analyze a number of equivalent real microstructure numerically, which give an easy access to the details of molten CMAS flow kinetics that are not easily accessible by experimental means. The know-how may provide qualified information for custom-tailoring the morphology of the EB-PVD top coat that potentially acts as a barrier to the CMAS infiltration.

So far, no detailed modeling approach was found that captures microstructure-sensitivity of the CMAS flow kinetics, however, the demand of an advanced modeling framework is increasing. In the present work, an effort has been made to establish a preliminary computational and numerical framework for such analysis. The following main objectives are focused:


insight into the flow behavior with respect to the individual morphological parameter will be obtained.

#### **2. Materials and Methods**

#### *2.1. Material and Microstructure*

Thermal barrier coating (TBC) systems are generally comprised of three layers: (a) a bond coat layer (BC) of a Ni-rich nickel aluminide or MCrAlY on Ni-based superalloy substrate; (b) thermally grown alumina oxide layer (TGO); and (c) a top coat (TC) typically of 7YSZ. Figure 1 shows an EB-PVD fabricated 7YSZ-TBC on a Ni-based superalloy component. The top coat consists of a very fine feathery columnar microstructure. The thickness of the columns decreases from top to bottom, and the size of the porous channel reduces accordingly. At the bottom, very dense and fine thin columns are generated. Very fine nanopores can be identified at this zone.

**Figure 1.** Cross-section of electron beam physical vapor deposition (EB-PVD) fabricated thermal barrier coating (TBC) showing layered structure with deposited calcium–magnesium–alumina–silicates (CMAS) on top.

Variation of top coat microstructures can be obtained by custom tailoring the columnar morphology with controlled EB-PVD process parameters [25,32]. Several tailored microstructures with different feathery columnar lengths are shown in Figure 2. The morphological details can be identified from the top and cross-sectional views.

**Figure 2.** Tailored TBC microstructures fabricated by EB-PVD. Top view (upper A1, B1 and C1) and cross-sectional views (bottom A2, B2 and C2) are shown for: (**a**) Fine columnar microstructure; (**b**) coarse columnar microstructure; and (**c**) columnar microstructure with feathery arms (high porosity). Reprinted with permission from [25]; Copyright 2006 Elsevier.

Here, only three exemplary morphologies, named Fine, Course, and Feathery are shown, which vary with respect to the columnar microstructure. The cross-sectional views A2, B2, and C2 show numerous single crystalline columns and nano-sized feather arms along the complete column's periphery. The top views A1, B1, and C1 show different columnar tips, which also vary with respect to size and shape. The overall porosity of the top coat can be attributed to the gaps between columns, feathery arms, and gaps in the feather tips. Due to the morphological variations in Fine, Course, and Feathery columns, the overall porosity has been altered.

#### *2.2. Microstructure Modeling of EB-PVD Top Coat*

#### 2.2.1. Microstructure Idealization

Geometries of the top coat columns and feathers are highly complex. Moreover, the microstructure of each column varies slightly with each other. For modeling purposes, the microstructure variations need to be idealized in a statistical manner providing that the main structural features remain consistent. Therefore, a statistically representative microstructure unit cell (MUC) needs to be defined with much flexibility to incorporate necessary geometrical details as found in different types of EB-PVD top coats. Details of the MUC are given in Section 2.2.2.

To generate such a MUC, the microstructure parameters are identified from a high magnification SEM micrograph of an EB-PVD TBC column as shown in Figure 3a. The total column diameter is *D*TC, and the columnar diameter excluding the feather arm length is the reduced column diameter, *D*C, as shown in the figure. The "reduced column diameter" will be termed simply as "column diameter" in the following sections unless otherwise stated. Similarly, for the sake of brevity, the reduced column radius, which is the half of the reduced column diameter, will be termed as "column radius".

**Figure 3.** SEM image of a cross-sectional view of an EB-PVD 7YSZ column. The microstructures are shown for: (**a**) columnar geometry of the top coat (TC) with microstructural parameters identified for modeling purposes; and (**b**) closer look of the top view of pyramidal column tip, side view of feathers with irregular feather arrangements (top), and tortuous columnar gaps due to irregular feather arms (bottom). Reprinted with permission from [25]; Copyright 2006 Elsevier.

Other morphological parameters are total length of the column, the column length (*L*C), and gap between two columns, the columnar gap (*G*C). The thickness of each feather arm was named as the feather arm width or feather thickness (*T*F), the horizontal length of the feather arm along the column width is the feather arm length (*L*F), and the gap between individual feather arms is considered as void between feather arms (*V*F), and will be termed shortly the feather gap. Further, the inclination of feather arms with horizontal axis was considered as the feather inclination angle (θF), the vertical angle of the pyramidal tip is the feather tip angle (θT).

Figure 3b shows additional microstructural details with morphological variations of columnar tips, feather arms, and tortuosity due to feather arms. The pyramidal tip angle is around 30◦–45◦. The facets of the tips depict staggered arrangements causing gaps between pyramidal tips. The feather arms are inclined at a 40◦–60◦ angle. The length of the feather arm depends on this angle of inclination.

Bottom pictures of the Figure 3b show two cross-sectional views of the tortuous columnar gaps due to uneven staggered arrangements of feathery arms. The feather arms and feather tips are irregular in shape and size along this columnar gap, i.e., they vary at left and right columnar edge as well as from top to bottom of a column. The degree of irregularities causes high or low tortuosity along columnar gaps. The columnar gaps may show almost regular width along the column length with a low tortuosity over the column height, or they might be highly irregular with high tortuosity contour along the column length. It must be noted that all microstructural features are three-dimensional, and they may vary along the circumference of each column considerably. The microstructures analyzed here are two-dimensional cuts; however, several orientations of cuts can be taken into account.

From the geometrical details stated above, the most relevant morphological parameters influencing the CMAS flow behavior were identified as illustrated in Figure 4a. A python code was written to handle these variables in a parametric way for generating a number of synthetic (virtual) columnar microstructures resembling real EB-PVD columnar structures. Several exemplary models generated by this code are shown in Figure 4b. A number of possible morphological variations can be obtained by cross-mixing these microstructure parameters.

**Figure 4.** Model representation of the EB-PVD columnar microstructure idealized from real microstructural SEM images: (**a**) definition of morphological parameters for computer generated parametric model of EB-PVD columns; and (**b**) automated model generation using morphological parameter variations. Column structures with different feather geometries were obtained.

In these models, the feather arm length *(L*F*)* can be made short or long with sharp or rounded feather tips to obtain irregularities in feather arrangements as observed in the microstructure picture of Figure 3a,b. A sharp feather tip was defined by an angular cut at the feather arm. The top of the column with pyramidal tip was idealized as a triangular shape keeping the width of the pyramid equal to the columnar diameter. According to the microstructure picture, the pyramidal width increases gradually from top to bottom. This geometry was obtained by defining two gradual angles along the pyramidal edges. The width of the pyramid bottom can be shortened to obtain a large CMAS entry area over the columnar gap and enlarged for obtaining a smaller CMAS entry area.

In this parametric model, the feather arm angle is taken either 40◦ or 60◦. It should be noted that due to trigonometric functions in the parametric model, a higher feather inclination angle results in a thinner feather gap, and a lower inclination angle yields a wider father gap. In reality, the feather gaps are arbitrary, and are not directly correlated with feather inclination angles.

Total porosity of the model can be calculated adding the columnar gaps and voids between feather arms (feather gaps).

#### 2.2.2. Microstructure Unit Cell (MUC) Model

The columnar microstructural setting of the top coat was modeled by a periodic microstructure unit cell (MUC) consisting of symmetric parts of the columnar sub-structures arranged side-by-side, facing each other along the feathery arms (Figure 5). It was assumed that each microstructure is symmetric along the mid-axis of a column. This assumption conforms with our observation that the real EB-PVD TBC microstructure is also nearly symmetric at least along the two main directions, i.e., parallel to the rotational axis and perpendicular to it. For a 2D-simplification, only two EB-PVD columns having a particular columnar gap, *C*G, were considered as a statistical feathery unit of columns, as shown in Figure 5a. The 2D-MUC model can be arranged periodically in *X*-direction to obtain a complete EB-PVD top coat layer. In *Z*-direction, the columns are assumed to be infinite. For a 3D-simplification, four 1/8th symmetric parts of the feathery columns can be re-arranged in both *X*- and *Z*-directions, keeping the symmetry plane at the boundaries, and with a face-to-face setting of the feather arms separated by a columnar gap, *C*G. A more realistic EB-PVD top coat columnar microstructure can be obtained by this 3D-MUC model taking idealized symmetric microstructure in both *X*- and *Z*-directions, as shown in Figure 5c. In 3D-MUC, the tortuosity in the feathery structures can be modeled taking further geometrical complexity into account. For the present work, the investigations were done on 2D-MUCs only to obtain fundamental understanding of the flow behavior, however, advance analyses with 3D-MUC models are open for the future.

**Figure 5.** Microstructure unit cell (MUC) for infiltration analysis: (**a**) a two-dimensional (2D) representation of feathery columns with sideway periodic arrangements; (**b**) fine details of the feather tips idealized form real tortuous feathery structures; and (**c**) a three-dimensional (3D) representation of MUC consisting of four one-eighth feathery columns, arranged with alternate zigzag patterns of feathery structures.

For modeling simplicity, the columnar gap was assumed to be vertically straight, and has a constant radius over the length (Figure 3b). The tortuosity along the columnar gap was obtained by defining a particular arrangement of feather arm and feather tip as described in the following.

At first, the feathery arm patterns with non-uniformities were idealized as periodic irregularities, as can be seen in Figure 5a,c along the columnar gaps. A single set of opposite feathers is shown in Figure 5b with enlargement. Here, the feather-tips of the left column were made slightly longer and the feather-tips of the right column were made shorter such that the tips of the longer feathers extend slightly into the columnar gap, and the tips of the shorter feathers stay at the columnar gap. In the left column, a periodic arrangement of longer and shorter feather-arms was done over the length of columnar gaps (from top to bottom) to obtain a realistic irregular (tortuous) feathery microstructure. For the right column, the arrangement is with an opposite order, one feather was made shorter, and the next was longer, and they are successively arranged along top to bottom. Note that the feathery

structures of the opposite columns are not the mirror image along the columnar gap, which is further illustrated in the figures of the following sections.

For the numerical analysis of infiltration kinetics, only a 2D-MUC was used. The geometries of the top coat columns were varied in the MUC in order to obtain different types of EB-PVD morphologies.

#### *2.3. Numerical Approach*

#### 2.3.1. Modeling Approach for Molten CMAS Flow

The molten CMAS flow through EB-PVD top coats involves multi-physics Fluid (Liquid)-Solid Interactions (FSI), which were modeled using a coupled Eulerian–Lagrangian (CEL) approach available in the commercial software ABAQUS [33]. The top coat was defined as a solid structure in the Lagrangian framework allowing structural deformation and rotation under the elastoplastic material behavior. On the other hand, the molten CMAS was assumed to be a highly dense, compressible, viscous fluid. To model the severe deformation and distortion of the molten CMAS, an Eulerian approach was adopted. Further, the flow of CMAS was assumed to be Newtonian, as the infiltration of the CMAS through TBC occurs only by capillary actions in the presence of the low strain-rate, i.e., strain rate effects on the viscous flow was neglected. Only the temperature dependency on viscosity was considered.

To capture multi-physics interactions along the solid-liquid domain boundaries (here, the domain between top coat and CMAS), the free surface evolution of the liquid CMAS flow was modeled using a volume-of-fluid (VOF) method. In this approach, the prescribed liquid material is contained within a defined domain of Eulerian elements bounded by a large domain. The amount of liquid material that is filled or emptied in an Eulerian element is defined by the Eulerian Volume Fraction (EVF) for that element. An element is fully filled with a material if EVF = 1.0, and this element is empty (void element) if EVF = 0. EVF ranges within 0.25, 0.5, and 0.75 for partially filled material. The free surface of a liquid flow is approximated from the boundaries of fully filled, partially filled or voided elements. During the material flow, the free surface of the moving material is tracked, and approximated as a continuous flow-front.

Interactions between Lagrangian solids, and Eulerian liquids were modeled by defining particular contact interactions among these two domain boundaries. During the liquid flow, the evolving Eulerian surface exerts pressure on the solid surfaces according to the flow density and pre-defined friction behavior, and both the liquid-solid interface share a common deformation behavior. However, this approach requires intensive computational power. For the present work, some computational efforts were reduced by assuming the EB-PVD columnar structures as rigid bodies with no translational and rotational degree of freedom. This assumption was acceptable for the proposed analysis, as the goal was primarily focused on the infiltration kinetics of the molten CMAS, but not the deformation and stress evolution in the TBCs.

In Figure 6, the microstructure unit cell (MUC) was re-arranged for FSI modeling in a CEL framework. To define the Eulerian CMAS domain, an isosceles trapezoidal shape of the CMAS domain was identified according to the observation of the real microstructure (Figure 6a). Accordingly, an isosceles trapezoid sub-domain was defined over the void area of the MUC model as shown in Figure 6b. This CMAS sub-domain was assumed to be in the molten state, which has already been moved (or infiltrated) up to the column tips. To capture the primary effects of morphological influences on the infiltration behavior, it is sufficient to model the CMAS flow through the columnar gaps taking different feather morphologies.

**Figure 6.** Identification of the CMAS initial position for modeling: (**a**) trapezoidal shape of the CMAS domain is isolated at the top of the columns; and (**b**) microstructure unit cell (MUC) with CMAS sub-domain and Eulerian domains for material flow. Reprinted with permission from [25]; Copyright 2006 Elsevier.

For this simulation, an additional Eulerian domain with void elements need to be defined such that the domain occupies the columnar as well as feathery gaps, and it extends slightly beyond the Lagrangian domains by a small clearance. This clearance was required to ensure proper contacts in the fluid-solid interaction model. During simulations, the molten CMAS sub-domain flows through the Eulerian domain creating a moving free surface of the CMAS flow.

#### 2.3.2. Boundary Conditions

The applied boundary conditions are summarized in Figure 7. The columns are kept undeformed by applying rigid body constraints. The symmetric boundary condition was applied at the sides, and the fixed boundary condition at the bottom (Figure 7a).

**Figure 7.** Boundary conditions for the fluid-solid interaction model (FSI): (**a**) boundary conditions for the rigid top coat model; (**b**) applied capillary pressure at the bottom surface of the molten CMAS; and (**c**) velocity boundary condition for the Eulerian domain.

To simulate the infiltration due to the capillary pressure along the feather gaps, a suction pressure (acting downwards) was imposed on the Eulerian subdomain (Figure 7b). An approximate value of the capillary pressure was calculated using Equation (1), as employed by [8] for an analytical infiltration analysis of the CMAS flow.

$$
\Delta P = \frac{\sigma \cos \theta}{r\_{\rm eff}} \tag{1}
$$

Here, θ is the wetting angle of the molten CMAS, σ is the surface tension of the liquid CMAS, and *r*eff is the effective hollow radius between the columnar gaps. Assuming the fully wetted condition of CMAS, we obtained a zero wetting angle and cosθ = 1. The value of σ was taken as 0.438 N/m

from the approximated value calculated for CMAS in [18], and *r*eff was 18 μm from the microstructure analysis (Section 2.2.1).

Symmetry conditions of the Eulerian domain were implemented with velocity constraints in the required directions as follows: On *X*- and *Z*-planes the velocity is zero, i.e., *VX* = 0 and *Vz* = 0. In the *Y*-direction the velocity was kept unconstraint allowing an equilibrium liquid flow through the top (inlet) and bottom (outlet) surfaces (Figure 7c). The contact surfaces between the columns and Eulerian domains are assumed to be slightly rough where the interactions between these surfaces were defined by a Coulomb friction model taking a small friction coefficient 0.1. This assumption was taken due to the lack of sufficient knowledge that characterizes the CMAS-TBC contact surface interactions during the CMAS flow. As the friction coefficient is small, the flow behavior will be mainly controlled by the viscosity parameter of the CMAS liquid. For all the models the friction coefficient was kept the same, therefore, the results from these models were comparable.

2.3.3. Constitutive Behavior and Material Parameters of CMAS

• Equation of State (EOS)

The properties of liquid silicate were described by an equation of state (EOS). With the EOS a quantitative functional relationship between the volumetric properties of a phase, and its specified composition, temperature (*T*) and pressure (*P*) can be established. In the present work, the EOS properties of the molten CMAS were approximated from the phase relationship of molten silicate: 'Diopside (Ca2Mg2SiO6)–Anorthite (CaAl2Si2O8) eutectoid liquid' due to the fundamental similarities of chemical compositions and physical behavior of CMAS [34,35].

For this particular molten CMAS, the complex relationships of their state variables, e.g., pressure, volume, and temperature were described by a Mie–Gruneisen type EOS material model as given in the ABAQUS software [33]. The constitutive model assumes the pressure as a function of the current density and internal energy per unit mass [*p* = *ƒ(*ρ*,E*m*)*], and takes the following form:

$$p = \frac{\rho c\_0^2 \eta}{\left(1 - s\eta\right)^2} \left(1 - \frac{\Gamma\_0 \eta}{2}\right) + \Gamma\_0 \rho E\_{\rm m} \tag{2}$$

Here, *p* is the pressure, ρ is the reference density, η is the nominal volumetric compressive strain, *E*<sup>m</sup> is the specific energy, and Γ<sup>0</sup> is the Grüneisen ratio. The parameter *c*<sup>0</sup> is the bulk speed of sound in the material and *s* is a parameter determined from the shock compression test by relating the shock velocity (*Us*) and the particle velocity (*Up*). Plotting the *Up* in *X*-axis and *Us* in *Y*-axis, a linear relationship of *c*<sup>0</sup> and *s* can be defined, which is known as the Hugoniot equation, and written as follows:

$$
\mathcal{U}\_{\mathbb{S}} = \mathcal{c}\_0 + s\mathcal{U}\_{\mathbb{P}} \tag{3}
$$

For the present work, the Hugoniot parameters for CMAS are approximated from the molten silicate data [36–38], which for the parameters *c0*, *s,* and Γ<sup>0</sup> are: 3060 m/s, 1.36, and 1.52 respectively.

#### • Density and Viscosity of Molten CMAS

The density and viscosity of CMAS vary with their chemical compositions and are temperature dependent [39]. The density of laboratory produced CMAS in this work has been taken from a similar CMAS composition [40], where the density is reported to be in the range of 2.6–2.8 g/cc at 1300 ◦C. The approximated value for the present work is 2.69 g/cc at 1260 ◦C.

On the other hand, viscosity for the investigated CMAS was estimated by two methods: Experimental and theoretical. For a particular CMAS chemical composition (here, CMAS1: 41.68(SiO2)–24.56(CaO)–12.43(MgO)–11.05(Al2O3)–8.71(FeO)–1.57(TiO2)–0(SO3) at.%) in the temperature range of 1200–1400 ◦C, the obtained viscosity shows a rapid decrease with increased temperature [26]. Due to the complicated data measurement techniques at high temperature, only one

data was taken at each temperature. In a second approach, the viscosity was estimated theoretically using a well-established Giordano (GRD) empirical viscosity model that takes into account the compositional variation of CMAS and temperature dependency [41]. This model was verified on a large database of chemical compositions of molten silicates that were experimentally obtained.

In literature, the Giordano viscosity model has been widely used to estimate the temperature dependent viscosity of CMAS for a wide range of chemical compositions [42–45]. For the present work, the estimated GRD viscosity for the CMAS1 shows a lower viscosity compared to the experimental obtained viscosity data [31]. In Figure 8 both the experimentally estimated, and GRD model based viscosity is plotted for different temperature ranges. One can find that the viscosity prediction of the GRD model for CMAS1 at 1260 ◦C is ~0.5 Pa s, and the experimental value is close to 5.8 Pa s [26], which is almost 11 times higher.

**Figure 8.** Viscosities for CMAS1 measured experimentally and using the Giordano (GRD) model at different temperatures.

For the numerical simulations, we have considered two arbitrary viscosity parameters, representative to the GRD (lower) and experimental values (higher), to obtain a qualitative understating of the influence of viscosities on the infiltration behavior. Both the parameters are shown in Figure 8 at the vertical axis. Several preliminary studies of the infiltration behavior were performed. Calculations with the lower viscosity parameter were numerically faster. Therefore, most of the analyses were performed with the low viscosity parameter, from which a fundamental comparative understanding of infiltration kinetics with respect to the EB-PVD morphology was obtained.

#### 2.3.4. EB-PVD Microstructure Model for Infiltration Analysis

Geometric models having 40◦ and 60◦ feather inclinations with varying feather arm lengths (short and long), and columnar gaps (narrow and wide) were considered for the investigation of infiltration kinetics. For this purpose, six top coat model geometries were generated using a Python-based parametric model, where the morphological parameters were systematically varied. The models are shown in Figure 9. For each feather inclination angle (θF) two feather-arm gaps, *G*F, and two columnar gaps, *G*<sup>C</sup> were considered. Thus, changing the columnar gaps from narrow to wide and keeping the length of feather arms from short to long, the overall void vol.% (equivalent to porosity of the TBC) has been varied.

**Figure 9.** EB-PVD microstructure unit cell models including the TBC inter columnar gaps, molten CMAS domain and an Eulerian domain, generated using the Python script: (**a**) 40◦ inclined feather arm model with short arm and narrow columnar gap; (**b**) with short arm and wide columnar gap; (**c**) with long arm and wide columnar gap; (**d**) 60◦ inclined feather arm model with short arm and narrow columnar gap; (**e**) with short arm and wide columnar gap; and (**f**) with long arm and wide columnar gap.

The six models are named according to the variables used for this analysis: θ(Feather angle)–(Feather gap/void)*G*F–(Columnar gap)*G*C. Thus, for the 40◦ feather inclination model with double length feather arm and half columnar gap is denoted as: θ40–2*G*F–0.5*G*C. Figure 9a–c shows the variants of 40◦ feather inclination models, and Figure 9c–e shows the variants of 60◦ feather inclination models.

It should be mentioned that the higher feather angle (60◦) reduces the overall gaps between the feather arms, as can be seen in Figure 9b,e. During the model generation, the vertical gap parameter (*V*F) was kept the same for 40◦ and 60◦ inclinations. However, due to the angular decomposition of *V*F, the diameter of the feather gap for 60◦ inclined feathers (*D*60) was reduced to 0.94*D*40, where *D*<sup>40</sup> is the feather gap diameter at 40◦ inclination.

To study the CMAS flow behavior, data were evaluated along the columnar gap, represented with I1, I2, I3, and I4 sections, and inside the void between feather arms, represented with 1, 2, 3, and 4 numbers, in Figure 10.

**Figure 10.** Defined sections at the columnar gap and void between feather arms for numerical data evaluation during infiltration kinetics study

The time occupied by the molten CMAS to fill one Eulerian element with EVF = 1 located at the flow front was taken to be the infiltration time for that element. The distance from this element to the top surface was taken as the depth (or height) of infiltration. The total infiltration time is the time required for filling up the empty Eulerian elements along any direction, which includes longitudinal (vertical) infiltration in the columnar gaps and lateral (sideway) infiltration in the feather arm voids. It should be noted that the columnar height is the same for all the models but due to different feather morphology (e.g., feather inclination and length) the distance travelled by the molten CMAS within the feathery voids are not the same.

To obtain equivalent reliable numerical solutions from all the FE models, some pre-analyses were performed to optimize the mesh size, domain shape, and initial CMAS domain size. Only the optimized models were used for data evaluation and comparative analysis of the CMAS flow behavior. All simulations were performed using the low viscosity parameter until the CMAS infiltrated into the I4 section. For mutual comparison of infiltration depths, the data were evaluated after 80 μs of infiltration time.

In Figure 11, the infiltration along the longitudinal and lateral direction was shown for six models measured at 80 μs of infiltration time. Differences in the CMAS flow pattern can be well identified. The systematic analyses with the results are presented in the next section, where the infiltration kinetics with respect to the morphological variations is compared.

**Figure 11.** Infiltrated CMAS height at time of 80 μs simulated for six top coat morphologies. Infiltration depths are shown for: (**a**) Model\_40a; (**b**) model\_40b; (**c**) model\_40c; (**d**) model\_60a; (**e**) model\_60b; and (**f**) model\_60c.

#### **3. Results**

#### *3.1. Morphological Influences on Infiltration Kinetics of CMAS*

#### 3.1.1. Influence of Feather Inclination and Feather Arm Length on CMAS Infiltration

Infiltration depths along the columnar gaps in the longitudinal (vertical) direction and along the void between feather arms (feather gaps) for lateral (sideway) infiltration were analyzed. The estimated total time for longitudinal infiltration involves the time required to fill up the vertical gaps in columns, preceded by filling up the sideway gaps between feathers. In Figure 12, the longitudinal infiltrations are summarized for four models, with feather arm angles at 40◦ and 60◦, and with short and long feather gap. The columnar gap was the same for all models to gain a comparative overview of the morphological influences on the infiltration kinetics.

The columnar infiltration depth increased almost logarithmically with time, i.e., the infiltration was initially faster, and became slower as the CMAS approached towards the bottom sections. This effect came from the gradual increase of friction effects from the increased contact surfaces during material flow. Friction between the Eulerian material and column boundaries offers continuous resistance during contact, as a result the non-linear decrease of the flow was found. It can be seen that for the models with longer feather arms (or gaps), the overall vertical infiltration decreased and a lower infiltration depth was obtained for a particular infiltration time (for example, at 75 μs time). Comparing Figure 12a,b it can also be found that the infiltration depth is higher for 60◦ feather inclination compared to 40◦ inclined feathers. This result is compared with horizontal lines at an infiltration time of 75 μs in Figure 12a,b. The prediction indicates that a higher inclination angle of feather arms facilitates an overall easy flow of CMAS in the vertical direction.

**Figure 12.** Vertical infiltration in columnar gaps with respect to geometrical factors variations, e.g., feather angle, feather arm length: (**a**) comparison for 40◦ inclined feathers with short and long arm; and (**b**) comparison for 60◦ inclined feathers with short and long arm.

The lateral (sideway) infiltration through long and short feather gaps was also analyzed in order to obtain a comparative view of local flow characteristics inside different types of feather gaps. The results are summarized in Figure 13 in terms of the feather position from the top (feather height) and time of infiltration. Data points were taken at the end of feather gaps in feather sets (FS), marked by 1, 2, 3, 4 in Figure 10. In these simulations, primary infiltration occurred vertically in the columnar gap, and in the course of the CMAS flow the secondary infiltration occurred sideways along the feather gaps. In Figure 13a,b, the results were compared for infiltration through feather gaps for the third feather-set FS3 for 40◦ and 60◦ feather angles. The apparent relationship between the infiltration time and infiltration depth inside the feather gaps are marked for both the cases. As shown in the plots, the CMAS flow inside the longer feather gaps requires higher time to fill-up. For the longer feather arm model, the dotted trend lines for data points (FS1–FS3) lies slightly above the trend line of the short feather arm model. The trend lines indicate that the infiltration depths at a particular time are slightly higher for longer feather gaps. This result is due to the slight differences in the contact surfaces for the long and short feather gaps. The short feather models exhibit slightly more tortuous flow path, which consists of several feather tips with complex geometry. The frictional effect on the flow front is affected by this geometry; as a result, slightly lower infiltration height was achieved in this model.

**Figure 13.** Lateral (sideway) infiltration in feather gaps with respect to geometrical factors variations, e.g., feather angle, feather arm length. Data are measured at four feather sets FS1–FS4: (**a**) comparison for 40◦ inclined feathers with short and long arm; (**b**) comparison for 60◦ inclined feathers with short and long arm.

#### 3.1.2. Influence of Columnar Gap (Narrow and Wide) on CMAS Infiltration

The influence of the columnar gap on the infiltration behavior was studied for narrow (*G*<sup>C</sup> = 0.9 μm) and wide (*G*<sup>C</sup> = 1.8 μm) columnar gaps. Narrow columnar gap models were obtained simply by halving the inter-columnar gaps from the previous models. Models with short feather arms taking 40◦ and 60◦ feather inclinations were used. In Figure 14a,b, the infiltration depths with respect to time were summarized. For a particular infiltration depth large differences in the infiltration time were obtained between narrow and wide columnar gaps.

**Figure 14.** CMAS infiltration behavior for narrow and wide columnar gaps for short feather arm model: (**a**) Comparison for 40◦ inclined feather arm models; (**b**) comparison for 60◦ inclined feather arm models.

In Figure 14a, the infiltration time was almost 3.5 times higher for narrow columnar gaps to travel an infiltration depth of ca. 14 μm. By normalizing the data for total infiltration time, it can be shown that the infiltration time was increased by 250% in narrow columnar gaps (*G*<sup>C</sup> = 0.9). The feather inclination was 40◦ for these models. Similarly, it can be shown from Figure 14b that the infiltration time was increased by ca. 200% in the narrow columnar gap for an infiltration depth of ca. 14 μm considering 60◦ feather inclination. The results showed that the reduction of the columnar gap width by a factor of 2.0 could reduce the infiltration time significantly. Comparing the results for feather inclinations one may also find that the models with 40◦ feather inclination showed reduced overall flow, which means, this model microstructure was more resistive to the molten CMAS infiltration compared to the 60◦ feather inclination model.

#### 3.1.3. Microstructure Influence on the Infiltration Rate and Velocity

The morphological influence on the local infiltration kinetics can be understood clearly by analyzing the instantaneous rate of infiltration through columnar gaps for different feather morphologies. The instantaneous rate of infiltration was defined as the incremental change (Δ*h*) in the columnar infiltration depth of the molten CMAS with respect to the incremental infiltration time (Δ*t*). The infiltration rate was measured considering the morphological variations of feathery columns obtained by changing the feather arm size (for short or long *G*<sup>F</sup> = 1 or 2) and feather inclination angle (θ<sup>F</sup> = 40◦ or 60◦).

In Figure 15, the infiltration rate with respect to the total time of infiltration was summarized for these models. From the plots (a–d) exponential decreasing trends of the curves were observed. The curves showed very high infiltration rate at the top of the columnar gaps just below the throat area. At this location only a little contact surface between TBC and CMAS was present, therefore, less frictional drag was activated allowing faster CMAS flow rate.

**Figure 15.** Change of infiltration rate through columnar gaps (vertical infiltration) influenced by different feather morphologies. The results are shown for: (**a**) short feathers with 40◦ feather inclination; (**b**) long feathers with 40◦ feather inclination; (**c**) short feathers with 60◦ feather inclination; and (**d**) long feathers with 60◦ feather inclination.

The rate of flow drops drastically when the contact surface was successively increased with the moving flow-front of CMAS. After filling up several feather gaps, the slope of the curves decreased gradually due to the decreasing flow rate and finally the rate became the slowest at the bottom of the columnar gap.

The rates of infiltration showed locally up and down zigzag patterns in the curves. These undulations were due to the interactions of CMAS at each feather-arm tips it encountered. At this location the primary flow path at the vertical direction put forth a new flow-path along lateral (sideway) direction through the feather gaps. Eventually, new contact surfaces came gradually in contact with the sideway flow which reduced the overall flow kinetics. As a resulting effect, the vertical infiltration became temporary stagnant or progressed with a minimal rate. Once the feather gaps were filled with CMAS, the vertical flow became faster again. The zigzag pattern was the result of slow-fast infiltration rates along the columnar gaps due to the branching of the CMAS flow path. The slow rate was induced by the lateral movement of CMAS through feather gaps having different morphological variations.

For a comparative study of the overall rate of the vertical infiltration with respect to the variations of feather arm length and inclination, the previous flow curves were reorganized with respect to the trend-line behavior of infiltration kinetics. The trend-lines basically approximate the curve behavior by fitting polynomial functions based on the available data points, thus the zigzag data patterns were idealized. In Figure 16, the morphological influences on the flow behavior can be recognized more clearly, from which two conclusive results can be drawn: (i) changing the feather inclination angle from 40◦ to 60◦, the vertical infiltration rate (overall infiltration) at column gaps became faster; and (ii) increasing the feather length from short to long the vertical infiltration rate becomes slower.

For the first case, the 60◦ feather inclination narrowed down the feather gap diameter slightly (can be seen in Figure 9). As a consequence, the frictional drag per unit area was increased; therefore, the CMAS flow rate was reduced inside the feather gaps. However, the overall porosity is reduced, therefore, during less lateral flow in the feather gaps, more liquid flows into the vertical gaps, resulting in a faster flow rate (on average) compared to the 40◦ feather angle model.

**Figure 16.** Simulated infiltration rates in columnar gaps compared for short and long feather arm models considering 40◦ and 60◦ feather arm inclination.

For the second case, the longer feather arm increased the feather gaps. The frictional contact surface that the flow-path encounters becomes longer as well. In this model, during the CMAS flow the vertical flow-path subdivided into wider and longer branches along the feather gaps, reducing the effective flow rate along the vertical direction.

To understand the overall rate of lateral (sideway) infiltration, the infiltration velocity was computed along the feather gaps for each feather set (FS). An average infiltration velocity was estimated by measuring the required time for distance travelled by CMAS at the feather gaps. In Figure 17, the results were summarized for the four models. As the distances of the feather arm sets (FS) from the top position were different, different infiltration depth for FS1–FS4 was obtained. To compare the results, the trend-line behavior of the velocity data was analyzed.

**Figure 17.** Lateral (sideway) infiltration in feather gaps showing velocity reduction in feather arm sets (FS) and the change of velocity gradient: (**a**) Comparison for 40◦ inclined feathers; (**b**) comparison for 60◦ inclined feathers.

The trend-lines showed that the flow was faster at the first feather gaps, and with the course of time, it was reduced successively for the second and third feather gaps. The reduction rate of the velocity depends on the frictional surface areas against which the CMAS interacts during flowing. For longer feather gaps the velocity reduction is higher, as the CMAS has to travel longer distances by overcoming longer frictional flow paths.

It should be noted that the high velocity (faster flow) at the FS1 is influenced by the wide CMAS entry area, which is slightly different for each model (recall Figure 11). The feather set FS3 is away from the CMAS entry area, therefore, at this position the flow becomes steady. This steadiness with respect to the feather gap size and inclination is compared with the velocity gradient (*V*<sup>G</sup> = Δ*v*/Δ*h*) data, as shown in Figure 17.

#### *3.2. Influence of Viscosity on CMAS Infiltration*

The influence of viscosity on the flow behavior for long feather arms and 40◦ feather inclination (model: θ40–2*G*F–*G*C) was investigated by using two viscosity values: Low and high as stated in Section 2.3.3. The characteristics of the flow kinetics were distinguished by comparing the flow rate compared to the infiltration depth relationship, as shown in Figure 18. Here, the instantaneous rate of infiltration (Δ*h*/Δ*t*) was plotted, which captured the local change of the infiltration rate during traveling through the columnar gaps. The instantaneous rate gives a better understanding of the flow kinetics than the average rate of change, as for the earlier case the local variation can be better explored and understood.

**Figure 18.** Decrease in the infiltration rate with respect to the infiltration depth measured in columnar gaps based on low and high viscosity values for 40◦ feather inclination, 1.8 μm columnar gap and long feather arm condition.

As seen in Figure 18, the low viscous CMAS entered the column with a high flow rate and then the flow was gradually subsided while filling up the columnar gap. On the other hand, the high viscous fluid entered the columnar gap with an almost 8–9 times slower flow rate, as shown at a distance of 2 μm infiltration depth in the figure. After filling the columnar and feathery gaps at a distance of 20 μm, the higher viscos CMAS almost became sluggish, but the low viscous CMAS continued further with a slower flow rate.

For the high viscosity parameter a higher flow resistance was imposed on the molten CMAS, which eventually reduces the flow velocity inside the columnar and feathery gaps.

A comparison of the flow velocity for low and high viscosity parameters along the columnar gaps shows a similar decreasing pattern as given in Figure 18, and therefore, not presented here.

Further, the lateral (sideway) infiltration behavior for low and high viscosity parameters was studied for two different morphologies: One is for short and the other is for long feather arms (model: θ40–*G*F–*G*<sup>C</sup> and θ40–2*G*F–*G*C), with the feather inclination angle at 40◦. The infiltration time required for filling up the void between feather arms considering short and long feather arms were compared. In Figure 19, the total infiltration time was recorded for filling up each feather set (FS) in the θ40–*G*F–*G*<sup>C</sup> and θ40–2*G*F–*G*<sup>C</sup> model. This infiltration time also includes the time of vertical infiltration, as the vertical and horizontal infiltration occurs successively.

In Figure 19a, the short feather arm model shows a slower infiltration for a high viscous CMAS compared to a lower viscous CMAS, as expected. The dotted lines indicate that the infiltration time increased exponentially for the high viscous CMAS while filling the longitudinal and lateral voids (in columnar and feathery gaps), while a moderate increase of time was observed for the low viscous CMAS. For the low viscous CMAS (Figure 19b), the infiltration time remained qualitatively the same for both short and long feather gaps. For the high viscous CMAS, the infiltration time increased remarkably for long feather gaps compared to the shorter one.

**Figure 19.** Comparison of CMAS infiltration into the feather gaps for different feather sets using a low and a high viscosity parameter simulated in a model with 40◦ feather inclination and 1.8 μm columnar gap: (**a**) comparison of infiltration time to fill up three short feather gaps; and (**b**) comparison of infiltration time to fill up three long feather gaps.

#### *3.3. Qualitative Comparison with Experimental Infiltration Behavior*

The CMAS infiltration kinetics was studied experimentally in previous researches [19,24] for the case of coarse and feathery morphologies (see Figure 2b,c). The overall porosity (gaps between columns and feathers) in the feathery structure was almost 2.5 times higher than that of the coarse structure. The main goal of the experimental investigation was to identify which microstructure resists CMAS infiltration depending on the geometrical factors. During the experiment infiltration was carried out for different times at a constant temperature. Then for both microstructures, the infiltration depth was measured by means of SEM imaging with the help of the energy-dispersive X-ray spectroscopy (EDS) analysis. The maximum infiltration depth of CMAS is normalized for both microstructures for a comparative study. If the maximum infiltration depth of CMAS is considered to be 100% (for a coating thickness of 400 μm), and the maximum infiltration time is 100%, then the infiltration rates for coarse and feathery structures can be qualitatively compared, which are shown in Figure 20a.

**Figure 20.** Comparison of infiltration depths and infiltration times for 'feathery' and 'coarse' morphology: (**a**) experimental results (normalized) at two columnar distances; and (**b**) simulation based prediction (normalized) at two columnar distances.

It was found that after 50% and 100% infiltration time, the difference in infiltration depth between feathery to coarse microstructure was found to be 50%. At any point of the infiltration time the feathery morphology exhibits almost 50% lesser infiltration depth.

For the numerical study, an equivalent microstructure for coarse and feathery columns was obtained from the θ40–*G*F–*G*<sup>C</sup> and θ40–2*G*F–0.5*G*<sup>C</sup> model. For both the models, simulations were performed only up to 15 μm infiltration depths, which cover full infiltration beyond four feather-sets. The respective variations of infiltration depths for two intervals of the infiltration time similar to the experimental case were compared. The results are depicted in Figure 20b. The difference of infiltration depth was found to be 46.8% at 100% infiltration time, which is close to the experimental result. At 50% of the maximum infiltration time the difference was found to be 48.21%. In both experimental and simulation results the infiltration depth was found to be ~50% lower in the feathery (feathers with double length) microstructure.

The numerical prediction of infiltration depth at 100% infiltration time predicts the experimental observation reasonably well, but for 50% infiltration time the model showed higher infiltration depth compared to experiments. This discrepancy can be explained on the basis of the melting nature of the CMAS melt. It is to be noted that at 1225 ◦C, the CMAS melt was in semi molten state (crystalline + glass) and longer heating times were required to make it fully glassy in nature which was the case of 20 h infiltration results (100% infiltration at 1225 ◦C). For the model simplicity, it was assumed that the whole CMAS was in a glassy and flowing condition, which in turn has resulted in higher infiltration depth than in the experiments. However, it was shown that the comparative trend of infiltration kinetics matches very well with the experimental observation.

#### **4. Discussion**

The CMAS infiltration kinetics was assessed for an EB-PVD TBC microstructure with varying morphology, which includes geometrical features for feather arm length, feather inclination angle, feather gaps, columnar gaps, etc. Based on the model parameters the infiltration behavior was quantitatively estimated. Some qualitative understanding of the infiltration kinetics with respect to the morphological variations was also gained. Some important results are discussed as following.

#### *4.1. Influence of Columnar Gap width on Infiltration Behavior*

Columnar gaps were arranged vertically in the models as described previously in Figures 3 and 4. The tortuosity along the column height was modeled by arranging the longer and shorter feather arms alternatively. The negative pressure force (equivalent to a constant capillary force) was active in this gap. The infiltration flow was obstructed by the long and short feather arms with sharp feathery tips, and branched by the feather gaps with different inclination. The contact surfaces imposed frictional drag during flow progression. For a wide columnar gap the unit frictional surface area with respect to the total surface area is reduced because the ratio of the friction surface length and the columnar diameter (*L*/*D*C) is less compared to a narrow columnar gap. That means, for wide columnar gaps a unit volume of liquid is exposed to less frictional drag during infiltrating. For a narrow columnar gap the overall porosity is reduced, and the frictional effect is significantly higher for a unit volume of liquid to enter, which obviously increases the infiltration time (Figure 14).

#### *4.2. Influence of Feather Length (Gap) on Total Infiltration Behavior*

The vertical infiltration time estimated in the models is the time required to fill-up all the voids, e.g., columnar gaps arranged in vertical directions and feather gaps arranged laterally (sideway). Therefore, the vertical infiltration time is considered as the time of total infiltration. It was shown that for a constant columnar gap the variations in feather morphology, e.g., feather arm lengths and feather inclination influence the infiltration kinetics remarkably.

The feather arm length basically defines the length of feather gaps between two feathers. For a set of long feather arms the void between feather arms increases. Infiltration time is expected to be higher for this type of long feather model. Similar behavior was also observed in the experiments comparing coarse and feathery structures of the top coat columns [19,24]

As the contact surface areas between the CMAS flow path and the top coat boundaries retard the flow rate due to frictional effects, for an increased contact surface area slow infiltration was obtained. Therefore, a large feather arm model containing large feather gaps and large contact surface areas showed reduced overall flow rate.

Another morphological factor, the feather inclination, also affects the infiltration kinetics, basically due to the narrow feather gaps obtained during model generation. Narrow feather gaps increase the frictional drag effects per unit area as discussed in Section 3.1.3. As a result, the CMAS flow rate inside a more inclined feather gap was reduced. For 40◦ inclined feathers the gap was wider, and thus larger, therefore, the acted frictional drag on CMAS was less, and inside the feather gap flow was faster.

However, for the 60◦ inclined feather model the overall porosity is smaller compared to the 40◦ inclined feather as the both have the same inter-columnar gaps. As a result, the vertical flow rate, which is the resultant of the flow along inter-columnar and feather gaps, was higher for the 60◦ inclined model.

It should be noted that the narrow feather gap for a higher feather angle is only due to the trigonometric conversation of data in the parametric model output. In reality, an inclined feather can possess any arbitrary feather gaps, not necessarily a narrower feather gap.

#### *4.3. Influence of Morphology on the Flow Velocity and Rate*

For all of the models the rate of infiltration was higher at the very beginning, and it was reduced gradually with respect to the vertical infiltration along the columnar gaps. This behavior was obvious as the CMAS entry area offered less frictional contacts, and for the progressive CMAS flow increased frictional contact areas were activated on the fly.

It was shown that the flow velocity in feather gaps and columnar gaps are dependent on the gap diameter, which determines the vol.% of porosity and effective frictional surface encountered by the flow-front. For both the columnar and feather gaps the flow velocity was reduced for narrower gaps due to large frictional effects per unit area. For a more vertically inclined feather arms (e.g., 60◦ inclined), the flow at the narrow feather gaps was reduced but the vertical infiltration occurred continuously showing an overall increased flow rate (Section 3.1.3). Therefore, we observed higher infiltration depth in the 60◦ feather angle model at a particular infiltration time (Figure 12).

#### *4.4. Influence of Viscosity*

The known effects of viscosity on the CMAS flow kinetics in TBCs were once again proven with the developed model showing that the CMAS with less viscosity has the faster infiltration rates. For a high viscous flow, the effects of frictional drag increase exponentially with the distance travelled by the CMAS flow. For the low viscosity parameter, the frictional drag increases only slightly during continuous flow along the frictional flow path.

No linear relationship between the flow rates can be found comparing the results of low and high viscosity values. That means, a linear increase of viscosity parameter will not increase the infiltration time linearly; rather an exponential increase of infiltration time can be expected.

#### *4.5. Suggestions and Remarks toward Future Direction*

The presented simulations captured a quantitative picture of microstructure sensitivity of the CMAS infiltration behavior. Valuable information can be extracted from the predictive analysis, from which the critical morphological features that need to be customized for mitigating CMAS infiltration can be identified.

For example, the results obtained from the six idealized 2D models suggest that EB-PVD top coats may comprise some morphological features, like narrow inter-columnar gaps, long feather arms, lower feather inclination angles, and above all higher lateral porosity, in order to reduce the CMAS intake through infiltration. Simulation-based knowledge can thus be used for advanced TBC design against environmental attack.

It should be noted that the predictability of the present model is restricted due to several idealizations in flow behavior and morphology. These limitations can be improved by employing detailed physics-based modeling incorporating complex 3D representations of the TBC microstructures. In particular, a thermodynamic equilibrium based melt/coating interaction model taking the thermal gradient as well as the temperature dependent viscous flow of CMAS, considering crystallization dynamics during infiltration needs to be developed. Realization of such thermo-chemical interaction modeling requires an enhancement of the presented multi-physical computational approach, which is currently in progress and remains beyond the scope of this work.

The computational method presented here is not limited to the CMAS-TBC interaction analysis only. This modeling approach can be extended for the analysis of CMAS-EBC (environmental barrier coating) interactions with few modifications in the microstructure model.

#### **5. Conclusions**

In this work, a numerical approach was developed for studying the infiltration kinetics of molten CMAS considering the morphological variation of EB-PVD TBCs. For a numerical representation of the top coat TBC microstructure, a microstructure model was established by rational idealization of morphological features from the SEM images. Six models were generated for the TBC-CMAS interaction analysis. Detailed analyses of the flow kinetics were performed, and then the results were verified by experimental infiltration depths of the feathery and coarse microstructures. Based on the predicted results the flow characteristics regarding infiltration depth, rate, and velocity were evaluated with respect to the columnar microstructural features of the top coat, enabling a qualitative comparison of the microstructure sensitive flow kinetics of CMAS. An in-depth understanding of the CMAS flow behavior was gained, and several morphological features of the EB-PVD TBC controlling infiltration kinetics were identified. The results are summarized as follows:


In future, the microstructure sensitive flow kinetics of molten CMAS will be assessed using 3D microstructure models considering extended geometrical features of the TBC columns. Furthermore, the model can be used for the thermomechanical stress analysis to predict the CMAS assisted stress evolution and fracture in advanced TBC layers.

**Author Contributions:** Conceptualization, M.R.K. and R.N.; Methodology, M.R.K.; Validation, M.R.K. and N.R.; Formal Analysis, A.K.S.; Investigation, M.R.K., A.K.S. and R.N.; Writing–Original Draft Preparation, M.R.K.; Writing–Review and Editing, A.K.S., R.N. and U.S.; Supervision, R.N. and U.S.; Project Administration, R.N. and U.S.; Funding Acquisition, R.N. and U.S.

**Funding:** This research was funded by Deutsche Forschungsgemeinschaft (DFG) (No. Schu 1372/5-1).

**Acknowledgments:** The authors acknowledge the supports from ICAMS, Ruhr-University Bochum for providing short-term computer facilities.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Crack Healing in Mullite-Based EBC during Thermal Shock Cycle**

#### **Hyoung-IL Seo 1, Daejong Kim <sup>2</sup> and Kee Sung Lee 1,\***


Received: 9 August 2019; Accepted: 10 September 2019; Published: 17 September 2019

**Abstract:** Crack healing phenomena were observed in mullite and mullite + Yb2SiO5 environmental barrier coating (EBC) materials during thermal shock cycles. Air plasma spray coating was used to deposit the EBC materials onto a Si bondcoat on a SiCf/SiC composite substrate. This study reveals that unidirectional vertical cracks (mud cracks) formed after several thermal shock cycles; however, the cracks were stable for 5000 thermal shock cycles at a maximum temperature of 1350 ◦C. Moreover, the crack densities decreased with an increasing number of thermal shock cycles. After 3000 thermal shock cycles, cracks were healed via melting of a phase containing SiO2 phase, which partially filled the gaps of the cracks and resulted in the precipitation of crystalline Al2O3 in the mullite. Post-indentation tests after thermal shock cycling indicated that the mullite-based EBC maintained its initial mechanical behavior compared to Y2SiO5. The indentation load–displacement tests revealed that, among the materials investigated in the present study, the mullite + Yb2SiO5 EBC demonstrated the best durability during repetitive thermal shocks.

**Keywords:** crack healing; environmental barrier coating; thermal shock; Hertzian indentations; mechanical behavior; CMC composite

#### **1. Introduction**

Gas turbines that operate at high temperatures are being developed for and installed in stationary and aviation applications. Higher operating temperatures can increase the energy efficiency and contribute to a cleaner environment. However, because materials that can withstand temperatures greater than 1000 ◦C are required, studies to develop a ceramic-based material beyond a single-crystalline superalloy are needed. In particular, all-ceramic materials are likely to be used for hot gas components such as combustor liners, blades, and nozzles, which require durable materials that can withstand high temperatures.

Silicon carbide exhibits both excellent heat resistance and high strength at high temperatures in addition to exceptional thermal shock resistance [1,2]. Silicon carbide also exhibits high hardness at room temperature and excellent wear resistance. The problem of brittleness can be overcome through the use of silicon carbide fibers (i.e., SiC-fiber-reinforced SiC (SiCf/SiC) composites). Even if the silicon carbide matrix cracks, the fibers with relatively high mechanical strength support the mechanical load and suppress crack propagation, thereby preventing complete failure [3,4]. SiCf/SiC composites are also expected to resist high-temperature creep and fatigue.

However, high-temperature oxidation and hot corrosion remain unsolved problems impeding the use of SiCf/SiC composite materials in gas turbines. Because gas turbines operate in high-temperature environments for extended periods, the oxide film initially formed on their surface by the oxidation of silicon carbide is expected to function as a protective film that prevents further oxidation. However, in an environment where heating and cooling cycles are repeated, the entire system can be damaged

because of inadequate protection by the oxide film. In particular, when exposed to a steam atmosphere at high temperature, SiCf/SiC composite materials can undergo mass reduction due to high-temperature corrosion. This behavior is attributed to Si in a silicon carbide reacting with water vapor (H2O) in the atmosphere to generate hydroxides such as Si(OH)*x* and Si–O*x*–H*y*. The hydroxide escapes into a gaseous state, possibly causing damage where the gas exits the matrix. Therefore, mass reduction due to high-temperature reactions with water vapor should be prevented to improve the operating lifetime of turbines [5,6].

To prevent the corrosion of Si-based ceramics at high temperatures, studies on the environmental barrier coatings (EBCs) have been extensively investigated [7–13]. The application of an oxide ceramic layer with a thickness of several tens to several hundreds of micrometers onto the silicon carbide surface can prevent high-temperature oxidation and corrosion by blocking the oxygen and steam penetration. Mullite and rare-earth silicates, which have excellent chemical resistance at high temperatures, are well known to be suitable materials for EBCs. Kang Lee first introduced a method to inhibit the reduction of mass of Si-based materials by applying an EBC [6]. This approach is feasible because the EBC exhibits sufficient thermal and chemical resistance and good mechanical properties. Specifically, the EBC is required to have a low thermal mismatch with the sublayer, high chemical resistance to steam, low oxygen permeability, good phase stability, and high impact and wear resistance. However, delamination of the coating layer occurs during multiple heating cycles of a gas turbine system, where heating and cooling are repeated. This delamination is due to differences in thermal expansion coefficients and elastic moduli between the coating and sublayer materials. As the heating and cooling are repeated, the expansion and contraction of the coating layer constrained by the sublayer material is restricted, inevitably resulting in the accumulation of stress in the coating layer. The literature reports numerous accounts of cracks forming in EBCs [14–16], where mud-cracks form in EBCs [17]; however, they do not grow to cause interface delamination because the cracks are unidirectional and vertical. Nonetheless, if these cracks can grow along the interface during continuous thermal cyclic loading, they may lead to interface delamination.

Crack healing of materials is an interesting subject. In particular, crack healing of brittle materials such as concrete or engineering ceramics is an interesting approach to prolonging their life [18–20]. In the case of concrete, crack healing occurs extrinsically. A crack healing agent in the form of filled capsules or particles is added to the concrete matrix, where the capsules or particles are ruptured upon crack penetration if cracking occurs. The healing agent is then released to fill and heal the crack. However, in the case of engineering ceramics, crack healing occurs intrinsically via a high-temperature reaction such as oxidation. For example, an additive material can react with atmospheric oxygen at a high temperature and fill the crack when it is exposed to high temperatures in air for a prolonged period. The initial strength has been reported to be recovered by crack healing [18,21,22].

In this study, mullite-based EBCs were applied by air plasma spray (APS) onto a Si bondcoat on a SiCf/SiC composite. For comparison, a coating layer of Y2SiO5 was also prepared by APS. The prepared coating layer was subjected to repeated thermal shock cycles from 1350 ◦C to room temperature to investigate the effects of thermal shock cycles on crack density. The mechanical behavior of indentation load–displacement, as evaluated by Hertzian indentation, was subsequently analyzed to investigate the change in mechanical properties of the EBCs after the thermal shock cycles.

#### **2. Experimental Procedures**

#### *2.1. Material Preparation*

The starting commercial raw powders used in this study were Yb2O3 (3N, Kojundo Chemical Laboratory Co., Ltd., Saitama, Japan, submicron) or Y2O3 (5N, Hebei Pejin International Co., Ltd., Hebei, China, particle size of 3 μm) and SiO2 (Kojundo Chemical Laboratory Co., Ltd., Saitama, Japan). The Yb2SiO5 or Y2SiO5 powder was mixed in molar composition of Yb2O3:SiO2 = 1:1 or Y2O3:SiO2 = 1:1, respectively, with zirconia balls in a ball mill. Each powder was dispersed in alcohol in the ball mill

and mixed for 24 h. After mixing, they were dried in a fume hood for 24 h in air, for 4 h or more in a drying oven, and then sieved using a 60-μm sieve. The ball-milled powders were first heated at 600 ◦C for 1 h and then heated at 1400 ◦C for 20 h to synthesize Yb2SiO5 or Y2SiO5. The heat-treated powders were crushed and then sieved again.

The synthesized Yb2SiO5 powder was dispersed in alcohol with mullite powder (DURAMUL 325F, Washington Mills, NY, USA) and then mixed homogeneously in a ball mill with zirconia balls. The weight ratio between mullite and the synthesized Yb2SiO5 powder was 88:12. The mixed powders were subjected to the same conditions of drying and sieving previously described. The 100% mullite powder was also prepared from different mullite powders with an average diameter of 43 μm (3Al2O3·SiO2, Yichang Kebo Refractories Co., Ltd., Yichang, China, KB-M-003).

The organic additives, which included a binder (PVA, Sigma-Aldrich, St. Louis, MO, USA) and a dispersant (San NOPCO, Tokyo, Japan, CERASRERSE 5468CF), were added to each synthesized or mixed powder with distilled water and zirconia balls in a ball mill for 24 h. At this time, 5 wt % of the binder and 4 wt % of the dispersant were added to the powder mixture.

The slurry was fed at a constant rate of 10 L/h and then sprayed as a soft agglomerated powder using the rotating atomizer of a spray dryer (Sewon Hardfacing Co., Mokpo, Korea). At this time, the inlet temperature was controlled at 170–210 ◦C, the outlet temperature at 70–110 ◦C, and the rotation speed of the atomizer at 6000 to 10,000 rpm.

Heat treatment was carried out at 600–800 ◦C in air to remove the organic additives in the spherical agglomerates. The agglomerates were continuously heat treated at 1250 ◦C for 2 h to impart strength to the powders in air. The heat-treated powders were sieved to various sizes and then granules within the range of 10 to 63 μm were classified and used as a powder for coating.

Liquid silicon infiltration (LSI) SiC-fiber-reinforced SiC composite material was prepared at the Korea Institute of Energy Research by reactive sintering after Si liquid infiltration with a diameter of 1 inch (25.4 mm) and a thickness of approximately 3 mm as a substrate for the EBC. The preforms were prepared by laminating silicon carbide fibers (Ube Industry, SA grade, Japan) with [0◦, 90◦] stacking and then filling them with phenolic resin to prepare fiber-reinforced plastics. The fiber fraction was maintained at 40 vol %. The prepared preform was cured at a constant temperature of 180 ◦C while being compressed through vacuum bagging. Afterwards, the cured product was subjected to carbonization of the phenolic resin at 1000 ◦C under a nitrogen gas atmosphere. The SiC-fiber-reinforced composite material was then prepared by an LSI process in which elemental silicon was heated to 1450 ◦C in a vacuum atmosphere to melt and infiltrate.

Silicon bondcoat having a thermal expansion coefficient intermediate between the coating material and the substrate material was prepared from Si powder (Si, particle size 40 μm, Saint Gobain, Anderlecht, Belgium).

Sand blasting was performed with alumina sand to improve bonding strength between the substrate and the bondcoat and also between the bondcoat and the topcoat. The bondcoat material was coated with a high-speed flame coating method (HVOF) to a thickness of approximately 400 μm, and a topcoat was applied by APS (9 MB, Sulzer Metco Holding AG, Winterthur, Switzerland) to a thickness of approximately 400 μm.

Mullite, mullite + Yb2SiO5, and Y2SiO5, prepared as previously described, were coated on the bondcoat. The coating powders were supplied to the plasma stream at a constant feed rate of 10 g/min with a mixed gas (mixing ratio of Ar/He = 55/5), and the current and voltage were kept constant at 500–530 A and 99 V, respectively. The distance between the spray gun and the substrate to be coated was controlled to 100–120 mm, and the sprayed velocity was controlled at a rate of 1000 mm/s to obtain a coating layer of constant thickness.

#### *2.2. Thermal*/*Mechanical Characterization*

The surface of the coating layer was ground and then polished with 15, 6, 3 and 1 μm diamond suspensions to observe the cracks in the EBCs. The thickness was measured at each polishing step so that the degree of polishing was controlled to be the same for all samples.

Thermal shock tests were performed to evaluate the durability of the EBC against repetitive thermal fatigue. The thermal shock was induced using a laboratory-made horizontal thermal shock machine in which the chamber could be moved horizontally by a motor. The schematic diagram of the thermal shock test is shown in Figure 1a. The test samples were placed on a fixed hotplate such that the substrate was positioned to receive heat from the lower part of the hot plate, and the coating layer was exposed to the top to receive heat from the chamber. First, the test samples placed on the hotplate were maintained at a temperature of 1100 and 1350 ◦C in the upper coating layer after the heating chamber was closed; the heating rate was 5 ◦C/min. After the target temperature had been maintained for 10 min, the upper chamber was moved to expose the sample to room temperature for 4 min as shown in Figure 1b, and then the chamber was moved back to expose the sample at 1350 ◦C again. This cycle was defined as one cycle, and repetitive thermal shocks were applied for a maximum of 5000 cycles.

**Figure 1.** (**a**) Schematic of the thermal shock test from 1350 ◦C to room temperature; (**b**) temperature profiles during the thermal shock tests in this study.

We observed the surface using an optical microscope (Olympus, GX51, Tokyo, Japan) after thermal shock tests of 1000, 2000, 3000, 4000, and 5000 cycles. Cracks and damages were observed with an optical microscope, and the density of cracks was measured with an image analyzer. For this purpose, the surface was observed at 200× magnification using an optical microscope. The crack density was defined by measuring the number of cracks intercepting the baseline after the baseline had been marked at positions of 1/4, 1/2, and 3/4 of each line of the horizontal, vertical, and diagonal length, respectively. We defined the crack density as the number of cracks divided by the length of the baseline and expressed it in terms of the number per inch. The numbers of cracks were determined from 10 pictures, and the mean values were plotted. The crack morphology of the surface was also observed with a high-resolution optical microscope (VHX-5000, Keyence Co., Itasca, IL, USA) after thermal shock tests of 100 and 5000 cycles.

Scanning electron microscopy (SEM) was used to observe the cracks in detail after crack healing. The electron microscope (JSM-6010F, JEOL Co., Tokyo, Japan) was operated at an acceleration voltage of 10 kV. Energy-dispersive X-ray spectroscopy (EDS) was used to determine the element species with an acceleration voltage of 10 kV and a spot size of 3.0.

The phase of the coating layer was analyzed by X-ray diffraction (XRD, Rigaku Co., Ltd., Tokyo, Japan) to analyze the change in the phase before and after the crack healing. The XRD analyses were

conducted with Cu Kα radiation generated at 45 kV and 200 mA with a scan rate of 2◦/min over the 2θ range from 20◦ to 80◦.

Spherical indentations were induced to compare and evaluate the mechanical behavior before and after the thermal shock test. The schematic diagram of spherical indentation is shown in Figure 2. Indentation load–displacement curves were obtained during indentation tests. A tungsten carbide ball with a radius of *r* = 3.18 mm was attached to the jig of a universal testing system (Instron 5567, Boston, MA, USA) and the pre-load was applied to the coating layer to *P* = 10 N at a rate of 3 N/min. The load was then applied to *P* = 500 N at a rate of 0.1 mm/min. The load was removed at the same rate as soon as the load had reached *P* = 500 N. The displacement was measured with an extensometer during load changes. Tests were conducted on more than five areas of the coating layer of the test specimens, and the indentation load and displacement results were averaged and plotted. The relative hardness was calculated from relative displacement after unloading in each indentation load–displacement and compared [23–26].

**Figure 2.** Schematic of a spherical indentation using a WC ball with a radius of 3.18 mm in this study.

#### **3. Results and Discussion**

The cracks formed on the surface and cross section of the EBC coating layer by thermal shock were observed by optical microscopy, as shown in Figure 3. Thermal shock tests were carried out at 1350 ◦C on the coating surface of the sample and at 1100 ◦C on the bottom of the sample because a gas turbine operating at 1600 ◦C is cooled simultaneously by the cooling gas through the flow path. Mud-type cracks formed on the surface. After 100 thermal shock cycles, cracks developed in the EBC layer fabricated by APS, as shown in Figure 3a. We considered that these cracks developed because of the difference in thermal expansion coefficients between the EBC coating and the metal sublayer. Several researchers have observed and reported mud cracks in EBCs [8]. However, cross-sectional observations of these cracks reveal that they formed in a direction perpendicular to the interface, not parallel to the interface. The vertical cracks shown in Figure 3b propagated through 400 μm of thickness; their propagation stopped at the top of the bondcoat layer. Such cracks can be stable if they do not induce horizontal cracking along the interface between the EBC and the sublayer, but rather tolerate strain caused by stress, thereby contributing to preventing delamination of the coating layer [27,28].

**Figure 3.** Optical micrographs of the (**a**) surface and (**b**) cross section of EBCs after 100 thermal shock cycles at 1350 ◦C.

The results of the measurements of the crack density developed on the surface of the EBC by thermal shock are shown in Figure 4. The periodic repetition of heating and cooling in turbine systems, and especially an abrupt stop of the system, will exacerbate crack formation. Therefore, the stability of cracks was evaluated by thermal shock tests. The mullite-based EBCs (i.e., the mullite and mullite + Yb2SiO5 EBCs) were not delaminated when the thermal shock cycle of exposure to room temperature was repeated for 5000 cycles while the coating material was maintained at 1350 ◦C and the SiCf/SiC substrate material was maintained at 1100 ◦C. Even though a crack is generated initially, interface delamination does not occur because vertical cracks are formed. That is, an initially developed crack in the EBC layer is very stable. Although the Y2SiO5 EBC was stable for 4000 cycles, interface delamination occurred after 5000 thermal shock cycles. The crack density on the surface of a sample after each 1000 cycles in the thermal shock test was measured and plotted in Figure 4. As shown in the graph, the changes in crack density depend on the EBC material. The crack density did not substantially change compared with the initial crack density in the case of Y2SiO5, whereas the crack density of the mullite-based EBC increased slightly and then drastically decreased after 3000 cycles. Furthermore, it exhibited a substantial decrease in crack density after 5000 cycles.

**Figure 4.** Plot of crack density of EBCs after each 1000 cycles of the thermal shock test at 1350 ◦C.

The optical microscope images in Figure 5 show the surfaces of the mullite, mullite + Yb2SiO5, and Y2SiO5 after 2000 and 4000 cycles of thermal shock, indicating crack healing. Although obvious cracks are observed after 2000 cycles of thermal shock irrespective of the EBC material (Figure 5a–c), surprisingly, the surface cracks were not observed after 4000 cycles of thermal shock (Figure 5d–f). In the case of Y2SiO5, relatively large macrocracks developed and were maintained without the formation of microcracks; thus, the crack density did not substantially change. By contrast, the population of microscale cracks increased at 2000 cycles and new phases covered the surface completely at 4000 cycles in the mullite-based EBC.

**Figure 5.** Optical micrographs of the surface of the EBCs after the thermal shock test at 1350 ◦C: (**a**) mullite, 2000 cycles; (**b**) mullite + Yb2SiO5, 2000 cycles; (**c**) Y2SiO5, 2000 cycles; (**d**) mullite, 4000 cycles; (**e**) mullite + Yb2SiO5, 4000 cycles; and (**f**) Y2SiO5, 4000 cycles.

The cracks on the surface of the mullite-based EBC after 5000 cycles of thermal shock were observed with a high-resolution optical microscope, as shown in Figure 6. The optical microscope images of the surfaces of the mullite-based EBC after 100 cycles are also presented in the figure. Although mud cracks were noticeable after 100 cycles of thermal shock, the crack gaps were clearly filled after 5000 cycles of thermal shock.

**Figure 6.** High-resolution optical micrographs of the surface of the EBCs after the thermal shock test at 1350 ◦C: (**a**) mullite, 100 cycles; (**b**) mullite + Yb2SiO5, 100 cycles; (**c**) mullite, 5000 cycles; and (**d**) mullite + Yb2SiO5, 5000 cycles.

Figure 7 shows the cracks of the mullite-based EBC after 100 and 5000 thermal shock cycles, as observed at high magnification by SEM. Only cracks and grains are observed after 100 cycles of thermal shock; cracks are observed around the grain. By contrast, after 5000 thermal shock cycles, new phases are formed to fill the gaps and cover the cracks. The results of the EDS composition analysis of these phases indicate that they contain only Si and Al, which are the constituents of mullite. We therefore speculate that the mullite itself underwent a change when repeatedly exposed to high temperatures.

**Figure 7.** SEM micrographs showing the microstructure on the surface of the EBC after thermal shock test at 1350 ◦C: (**a**) mullite, 100 cycles; (**b**) mullite + Yb2SiO5, 100 cycles; (**c**) mullite, 5000 cycles; (**d**) mullite + Yb2SiO5, 5000 cycles.

Figure 8 shows the cracks on the cross-section of the mullite, mullite + Yb2SiO5, and Y2SiO5 after 5000 thermal shock cycles, as observed by SEM. Each arrow in the mullite, mullite + Yb2SiO5 indicate crack healed region. On the other hand, each arrow in the Y2SiO5 indicate cracked region. Only cracks are observed after 5000 cycles of thermal shock in the case of Y2SiO5 (Figure 8c). By contrast, after 5000 thermal shock cycles, new phases are formed to fill the gaps and cover the cracks, and this result is more pronounced in the mullite + Yb2SiO5 (Figure 8b).

**Figure 8.** SEM micrographs showing the microstructure on the cross-section of the EB after thermal shock test at 1350 ◦C: (**a**) mullite, 5000 cycles; (**b**) mullite + Yb2SiO5, 5000 cycles; and (**c**) Y2SiO5, 5000 cycles.

XRD analysis was performed on the EBC materials after 5000 thermal shock cycles for the analysis of the phases formed around the cracks, as shown in Figure 9. In the case of the mullite material, an Al2O3 phase was detected in addition to the mullite phase. Therefore, the newly formed phases in Figures 7 and 8 are attributed to mullite that melted under repeated thermal shock cycles at high temperature. We speculate that the low-melting-point phases containing SiO2 migrated into the cracks, filling the gaps of the cracks, and that the Al2O3 phases precipitated as crystalline phase during cooling. In particular, Al2O3 phases are considered desirable because they exhibit excellent oxidation and corrosion resistance at high temperatures and can therefore resist oxidation and corrosion of EBC materials after crack healing. Similarly, an Al2O3 phase was detected along with the mullite, Yb2SiO5, and Yb2Si2O7 phases in the case of the mullite + Yb2SiO5 material, as shown in Figure 9b. However, in the case of Y2SiO5, Y2O3, and SiO2 phases were detected along with Y2SiO5 and Y2Si2O7 phases, as shown in Figure 9c. The Y2SiO5 phase melts at high temperatures and is decomposed into crystalline Y2O3 and SiO2, similar to the process that occurs in the mullite-based EBC.

**Figure 9.** XRD patterns of (**a**) mullite, (**b**) mullite + Yb2SiO5, and (**c**) Y2SiO5 after 5000 thermal shock cycles at 1350 ◦C.

Rare-earth silicate materials such as Y2SiO5 and Yb2SiO5 have been reported to exhibit less mass loss than mullite in a steam atmosphere at relatively higher temperatures [11,17,29,30]. However, if cracks occur in the coating layer because of thermal shock, the crack healing ability of Y2SiO5 and Yb2SiO5 is considered to be relatively poor. Therefore, studies are needed to develop a eutectic phase or new crack healing agent. The mullite + Yb2SiO5 material suggested in this study exhibits high water vapor resistance because of the Yb2SiO5 phase and exhibits crack healing behavior because of the mullite phase. Therefore, the mullite + Yb2SiO5 can be used as a high-temperature EBC material with crack healing ability.

The surfaces of the EBC materials were contacted with a WC ball with a radius of 3.18 mm to a load *P* of 500 N and then unloaded to determine the mechanical behavior after crack healing, as shown in Figure 10. The graphs show the test results for the sample before the thermal shock and the sample after 5000 thermal shock cycles. The Y2SiO5 EBC shows a large change from elastic to plastic after 5000 thermal shock cycles, as shown in Figure 10c because the test was carried out on the delaminated coating samples. By contrast, the load–displacement behavior of the mullite-based EBCs, which did not exhibit interfacial delamination even after 5000 thermal shock cycles, did not substantially differ from their initial indentation load–displacement behavior, as shown in Figure 10a,b. In the case of mullite, the residual displacement related to the hardness of the material was slightly increased after thermal shock, and the tangential slope of the stress–strain at unloading, which is related to the modulus of elasticity, decreased after thermal shock [10,25,26,31]. Thus, the hardness and elastic modulus of the sample decreased by 5000 thermal shock cycles. The decrease in the mechanical properties is attributed to cracks even though the sintering of the material during repetitive exposure to high temperatures occurs [32,33]. However, the hardness and elastic modulus of the mullite + Yb2SiO5 did not substantially change even after 5000 thermal shock cycles. Therefore, the initial mechanical properties of the mullite + Yb2SiO5 were recovered by crack healing.

$$\mathbf{a})$$

(**c**)

**Figure 10.** Plot of indentation load–displacement curves of (**a**) mullite, (**b**) mullite + Yb2SiO5 and (**c**) Y2SiO5 using a WC ball with a radius of 3.18 mm at load *P* = 500 N.

Figure 11 shows the change in the relative hardness of the material, as revealed by analysis of the residual displacement of the indentation load–displacement curves in Figure 10. The hardness value before the thermal shock test is assumed to be 100%, and the hardness value after 5000 cycles of thermal shock is plotted in percent. As shown in the graph, the mullite + Yb2SiO5 maintains its initial hardness value. These results indicate that the mullite + Yb2SiO5 is a potential EBC material because gas turbine coating materials exposed to repeated thermal cycles must not only not exhibit interfacial delamination due to the thermal shock but also resist impact and wear due to external FOD.

**Figure 11.** Hardness changes of EBC materials during thermal shock tests.

#### **4. Conclusions**

The purpose of this study was to introduce an EBC to protect SiCf/SiC composites from high-temperature oxidation and corrosion and to investigate the durability of the coating when exposed to high temperatures and repeated thermal shocks. Crack formation, crack stability, and crack healing were investigated after specimens were subjected to repeated heating to 1350 ◦C followed by cooling. In addition, spherical post-indentation tests were performed to determine the change in mechanical behavior after thermal shock. As a result, the following conclusions were obtained:


around the gap of the cracks. A crystalline phase of Al2O3 was detected by XRD analysis. We considered that the mullite phase melted at high temperatures and filled the gaps of the cracks and that a crystalline phase of Al2O3 was precipitated during cooling, resulting in crack healing. Crystalline SiO2 and Y2O3 phases were observed in the case of Y2SiO5.

(4) The mechanical indentation load–displacement behaviors of the EBC samples were characterized using a spherical indentation method after the specimens were subjected to 5000 thermal shock cycles. The mullite-based EBC showed no substantial change in mechanical behavior even after 5000 thermal shock cycles, whereas the Y2SiO5 exhibited a change from elastic to plastic behavior because of interfacial delamination. In particular, the mullite + Yb2SiO5 EBC maintained its initial mechanical behavior even after 5000 thermal shock cycles because of its crack healing behavior.

**Author Contributions:** Conceptualization, K.S.L., H.S.; Methodology, H.S., D.K.; Software, D.K., H.S.; Validation, D.K.; Formal Analysis, H.S.; Investigation, H.S.; Resources, K.S.L.; Data Curation H.S., D.K.; Writing—Original Draft Preparation, K.S.L., H.S.; Writing—Review and Editing, K.S.L., H.S. and D.K.; Visualization, K.S.L., D.K.; Supervision, K.S.L.; Project Administration, K.S.L.; Funding Acquisition, K.S.L., D.K.

**Funding:** This research was funded by the IT R&D program of MOTIE/KEIT [NO. 20000192, Development on the crack self-healing technology of CMC composite and coating material for gas turbine hot gas component].

**Acknowledgments:** We wish to thanks to Sewon Hardfacing Company and KIER for coating and the substrate preparation.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Crack Initiation Criteria in EBC under Thermal Stress**

#### **Emi Kawai 1, Hideki Kakisawa 2, Atsushi Kubo 3, Norio Yamaguchi 1, Taishi Yokoi 4, Takashi Akatsu 5, Satoshi Kitaoka <sup>1</sup> and Yoshitaka Umeno 3,\***


Received: 12 September 2019; Accepted: 21 October 2019; Published: 24 October 2019

**Abstract:** For design of multi-layered environmental barrier coatings (EBCs), it is essential to assure mechanical reliability against interface crack initiation and propagation induced by thermal stress owing to a misfit of the coefficients of thermal expansion between the coating layers and SiC/SiC substrate. We conducted finite element method (FEM) analyses to evaluate energy release rate (ERR) for interface cracks and performed experiment to obtain interface fracture toughness to assess mechanical reliability of an EBC with a function of thermal barrier (T/EBC; SiC/SiAlON/mullite/Yb-silicate gradient composition layer/Yb2SiO5 with porous segment structure) on an SiC/SiC substrate under thermal stress due to cooling in fabrication process. Our FEM analysis revealed that a thinner SiAlON layer and a thicker mullite layer are most suitable to reduce ERRs for crack initiation at the SiC/SiAlON, SiAlON/mullite and mullite/Yb2Si2O7 interfaces. Interface fracture tests of the T/EBC with layer thicknesses within the proposed range exhibited fracture at the SiC/SiAlON and SiAlON/mullite interfaces. We also estimated the approximate fracture toughness for the SiC/SiAlON and SiAlON/mullite interfaces and lower limit of fracture toughness for the mullite/Yb2Si2O7 interface. Comparison between ERR and fracture toughness indicates that the fabricated T/EBC possesses sufficient mechanical reliability against interface crack initiation and propagation.

**Keywords:** fabrication process; layer thickness design; interface crack initiation/propagation; fracture toughness; energy release rate; finite element method analysis

#### **1. Introduction**

Silicon carbide (SiC) fiber reinforced SiC matrix (SiC/SiC) composites are one-third lighter and have approximate 100–200 K higher heat resistance than current heat-resistant super alloys (Nickel-based super alloys) [1–3]. Application of SiC/SiC to advanced hot-section components in airplane engines is expected for improving fuel consumption and curbing emission of carbon dioxide [2,3]. SiC/SiC composites react with oxygen to form silica, which hinders degradation of the composites at high temperature over 1373–1473 K [4,5]. In the condition of high temperature and humidity (i.e., combustion gas environment), however, the silica reacts with water vapor to form silicic acid (Si(OH4)) gas [4,5], leading to wall thinning of the composites. Thus, application of the SiC/SiC composites to the hot-section components inevitably requires environmental barrier coatings (EBCs).

To achieve superior environmental shielding performance and thermomechanical durability, EBCs consist of several layers with different shielding functions [6,7]. For using EBCs at about 1673 K, which is the durable temperature of the heat-resistant SiC fiber, it is of immense importance to select a material with the following properties as the water vapor shielding layer; significant wall-thinning resistance in the combustion gas environment, and coefficients of thermal expansion (CTEs) being close to that of the SiC/SiC. Klemm conducted burner-rig tests for various heat-resistant materials under the condition simulating the actual combustion gas environment to quantitatively evaluate relationship between CTEs and wall-thinning rate [8]. The result indicated that ytterbium silicate (Yb2SiO5; YbMS and Yb2Si2O7; YbDS) has prominent wall-thinning resistance. He reported that the resistance of YbMS is particularly outstanding and that the CTE of YbMS is slightly larger than that of the SiC/SiC substrate, which is almost equal to that of YbDS. Yb-silicate is unlikely to fracture by phase transition during thermal cycle because phase transition causing the large volume change does not appear in the material throughout a wide temperature range from room temperature to high temperature. On the basis of these features, Yb-silicate is considered to be an attractive material as the water vapor shielding layer. By the way, YbMS is expected to be applied as the thermal barrier coating (TBC) for SiC/SiC substrate because the thermal conductivity of YbMS [9] is smaller than that of Y2O3-ZrO2 (YSZ), which is used as TBC for Nickel-based super alloys.

Under the combustion gas environment, EBCs are exposed to a steep gradient of oxygen chemical potential (*d*μO). It is known that, in an oxide layer under the environment of high temperature and *d*μO, the oxide ion migrates from high oxygen partial pressure (*P*O2 ) to low *P*O2 while the cation migrates in the opposite direction following the Gibbs-Duhem equation. This suggests that oxygen permeation occurs because of migration of cations as well as oxide ions. In addition, the cation migration may promote phase decomposition and fracture of the multi-layer component. For improving the phase stability of EBCs, therefore, it is of primary importance to prevent the migration of cations. Kitaoka et al. conducted oxygen permeability tests for representative heat-resistant oxides (Al6Si2O13 (mullite) and Yb-silicate) under high temperature [10–12]. They clarified that the oxygen shielding of mullite was superior to that of Yb-silicate and that the mullite was decomposed due to the motion of Al ions from the low-*P*O2 to high-*P*O2 side producing deficiency of Al ions in the layer on the low-*P*O2 side. Kitaoka et al. demonstrated that the phase of mullite is stabilized by depositing a SiAlON layer on the low-*P*O2 side of the mullite layer because it feeds Al ions [13].

On the basis of the abovementioned studies, we propose an EBC with a function of thermal barrier as shown in Figure 1, which we call T/EBC hereafter. The structure consists of layers as follows (in order from the side of SiC/SiC substrate); the SiC layer for smoothing the surface of the substrate, the SiAlON layer for bonding with the substrate and ensuring structural stability of the mullite, the mullite dense layer for oxygen shielding, the Yb-silicate dense layer for water vapor shielding, and the layer of YbMS with porous segment structure for thermal shock resistance. Here, for decreasing thermal stress that occurs by the difference in CTEs of the coating layers and the substrate, the Yb-silicate dense layer is formed with gradient composition to have YbDS and YbMS on the substrate and top sides, respectively. Hereafter it is referred to as the Yb-silicate gradient composition layer. The topmost YbMS layer with a porous segment structure (topcoat) is employed for thermal shock resistance so as to relax strain due to rapid temperature change in the thickness direction associated with start and stop of the engine.

T/EBC is fabricated by depositing these layers on a SiC/SiC substrate at high temperature and cooling it to room temperature. During the cooling process, the thermal stress may cause interface crack initiation/propagation and delamination along interfaces between the layers with the function of environmental barrier. This leads to deterioration of the mechanical reliability of T/EBC. Thus, it is essential to design T/EBC in such a way that crack initiation and propagation are prevented during the cooling process.

**Figure 1.** Schematic of thermal/environmental barrier coating (T/EBC) consisting of SiAlON (SA; bonding layer), mullite (Mu; oxygen shielding), layer of gradient composition of Yb2SiO5 (YbMS) and Yb2Si2O7 (YbDS) (Yb-silicate gradient composition layer; water vapor shielding) and top YbMS with porous segment structure (thermal shock resistance).

In general, fracture of brittle materials follows the Griffith theory; i.e., a crack initiates or propagates when the energy release rate (ERR) exceeds fracture toughness for the corresponding fracture mode. While fracture toughness is a material property that should be evaluated in experiments, ERR can be calculated by numerical simulations through, for example, finite element method (FEM) calculations. Unlike conventional TBCs having simple structures, the mechanical behavior of T/EBC is expected to be complicated due to its complex structure consisting of multiple layers on a substrate. Numerical simulation can be an efficient tool to evaluate ERR even in such a complex structure.

In designing T/EBC structure with sufficient mechanical reliability, we must optimize a wide range of parameters, including conditions of multi-layer deposition processes and structural dimensions such as layer thicknesses. Layer thicknesses should affect much the mechanical state and play a crucial role in the mechanical stability of T/EBC under thermomechanical loading, which the component should undergo in operation. It is of importance to understand how such structural parameters influence the mechanical reliability. The aim of this study is to propose the condition of layer thicknesses for preventing interface crack initiation and propagation during the cooling process in the fabrication of the T/EBC. We perform FEM analysis for a T/EBC model with varying thicknesses of the coating layers to investigate the dependence of the layer thicknesses on ERR of crack initiation at interfaces, where thermal stress concentrates significantly. The analysis result facilitates the determination of layer thicknesses to prevent interface crack initiation. Then, the T/EBC is fabricated within the proposed thicknesses of the coating layers. We calculate ERR for interface crack initiation and propagation by conducting FEM analysis to a model of the fabricated T/EBC. Comparing the ERRs and fracture toughness obtained by interface fracture tests, we assess the validity of the proposed layer thicknesses.

#### **2. Experiment and Simulation Procedures**

#### *2.1. Fabrication Procedure of T*/*EBC*

In this study, the substrate was SiC/SiC (50 mm square and 3 mm in thickness) whose matrix was formed by chemical vapor infiltration and subsequent melt infiltration. Before T/EBC deposition, a SiC layer with an adequate thickness was deposited by chemical vapor deposition in order to smoothen the surface to be coated.

Each constituent layer of T/EBC shown in Figure 1 is a complex oxide or an oxynitride containing multiple cations. When these compounds are heated to a high temperature, the composition of vapor containing cation differs substantially from that of compounds in solid phase (incongruent vaporization). Thus, deposition methods accompanying the melting process of raw materials, such as plasma spraying, bring about the incongruent vaporization of the materials during deposition, resulting in significantly deviated compositions of coating layers.

In this study, therefore, dual electron beam-physical vapor deposition (EB-PVD) was employed to prepare layers of compounds which were evaporated incongruently. Two target raw materials located in the deposition chamber were evaporated using two electron guns regulated individually, which enabled independent control of vapor pressures of gases generated from the targets. Compound layers with arbitrary compositions were deposited on a substrate placed at a proper distance from the targets. During deposition, small amount of source gas was introduced in the coating chamber, and the coated surface was heated up to about 1473 K by irradiation of a direct diode laser (wavelength of 915 nm) to promote full crystallization and densification of each layer. After the deposition finished, the specimen was cooled from the deposition temperature to room temperature at the average cooling rate of 20 K/min.

Table 1 lists the target materials (Ingots A and B) and source gases used in the deposition of the layers. During deposition of the Yb-silicate gradient composition layer, the ratio of two electron beam powers irradiated to two targets was gradually changed with coating time to control the composition toward the thickness direction. The YbMS topcoat was deposited on a rotating sample for out-of-surface to form the porous segmented structure owing to the shadowing effect.


**Table 1.** Raw material targets and source gases for preparation of T/EBC.

The cross-sectional microstructure of T/EBCs was observed using a scanning electron microscope (SEM; SU8000 SEM, HITACHI, Tokyo, Japan) with energy dispersive spectroscopy (EDS).

#### *2.2. ERR Calculation by FEM Analysis*

For ensuring the mechanical reliability of T/EBC in the fabrication procedure, we need to determine its layer thicknesses with which ERR for interface crack initiation due to the thermal stress, <sup>G</sup>th init, can be reduced so that interface crack is not likely to initiate. In general, it is difficult to fabricate T/EBC without any initial interface crack such as microcrack. Thus, under the proposed layer thicknesses for preventing interface crack initiation, evaluating the behavior of interface crack propagation is also required. This means that we need to calculate ERR for interface crack propagation due to the thermal stress, <sup>G</sup>th prop, and compare with it and interface fracture toughness.

In this study, in order to investigate the effect of layer thicknesses on <sup>G</sup>th init, we conducted two-step FEM analyses using a T/EBC model with varying layer thickness; (1) Thermal stress analysis: FEM analysis for evaluating the thermal stress state after the cooling process and (2) Interface crack introduction analysis: FEM analysis for calculating ERR by introducing interface crack in the thermal stress state that was obtained in Step (1). Then, for obtaining <sup>G</sup>th prop, we conducted the two-step FEM analyses to the T/EBC model with the optimized layer thickness to prevent interface crack initiation. ABAQUS (ver.6.14.6, Dassault Systemes, Vélizy-Villacoublay, France) was used in the calculations.

#### 2.2.1. FEM Model of T/EBC

The state of thermal stress in T/EBC is affected by the thicknesses of the coating layers. In particular, the difference in CTEs of the SiAlON and mullite layers is large compared to those at other interfaces (CTEs are explained in Section 2.2.2), and thus the thicknesses of the SiAlON and mullite layers are expected to significantly affect the state of thermal stress after cooling. We examined the simulation models mimicking T/EBC with various thicknesses of the SiAlON and mullite layers as shown in

Figure 2. Note that the FEM model used in this analysis is a model of the vicinity of interface edge of T/EBC deposited sample, which is indicated by a blue rectangle in Figure 2.

**Figure 2.** Schematic illustration of finite element method (FEM) analyses model and boundary condition.

In this study, the Yb-silicate gradient composition layer was modeled with five layers with different ratios of the YbDS and YbMS as follows; YbDS:YbMS = 100:0, 80:20, 50:50, 20:80, and 0:100. For computational efficiency and simplicity, the segmented YbMS layer was regarded as a homogenized material with an anisotropic (in-plane isotropic) mechanical property equivalent to the segmented layer [14]; hereafter it is referred to as the YbMS anisotropic homogeneous layer (Figure 2). The SiC/SiC substrate was treated as an in-plane isotropic material owing to the orientation of SiC fibers. T/EBC structure was represented by a 2-D simulation model, where displacements of the left and bottom edges of the model along the *x* and *y* directions, respectively, (shown in Figure 2) were constrained. A plane strain state was assumed to the *z* direction. In this analysis, the thicknesses of the SiAlON and mullite layers were changed (5, 15 and 25 μm), while the thicknesses of other layers were fixed as below; the SiC, Yb-silicate gradient composition and topcoat YbMS layers were 25, 100 (20 μm times five layers) and 200 μm, respectively. The thickness and width of the SiC/SiC substrate were 3000 and 4000 μm, respectively. The regions near all interfaces, where strong stress concentration is expected, were divided into a finer mesh (0.5 μm × 0.5 μm).

#### 2.2.2. Thermal Stress Analysis

Since the experimental temperature is decreased gradually during the cooling process in fabrication of the T/EBC as explained in Section 2.1, the temperature dependence of the elastic properties (Young' modulus and Poisson's ratio) and the thermal property (CTE) must be taken into account for the thermal stress analysis. In addition, it is necessary to take account of the effect of creep on the mechanical state for the conditions of the experimental temperature and the cooling rate in the present deposition process. Thus, in this study, a series of the thermal stress analyses was carried out to evaluate the thermal stress distribution in T/EBC, where the creep effect was considered. Thermal stress is accumulated in T/EBC throughout the cooling process; i.e., stress is assumed to be null at the high temperature in the deposition process, and we evaluate the thermal stress state in T/EBC after cooling to room temperature from the high temperature. The thermal stress is dominated by stress in the *x*-axis direction in Figure 2 mainly, and thus the loading parallel to the interface, in other words, mode II-rich thermal loading is applied to T/EBC.

The temperature condition for this analysis was determined from the experimental condition of the cooling process after deposition of T/EBC (cooling from a high temperature in the deposition process to room temperature at cooling rate of 20 K/min). A uniform temperature distribution in T/EBC was assumed during the cooling process. Although the surface temperature of the coated layers is up to 1473 K in the deposition process as described in Section 2.1, in this analysis we set the deposition temperature (the initial temperature) to be 1673 K, which is an assumed operation temperature of T/EBC and provides an estimation on the safer side. In addition, the cooling rate of 20 K/min can cause a considerable effect of creep during the cooling process. In general, it is known that creep emerges significantly at temperatures higher than half the melting point [15]. Thus, in this analysis, the stress was set to be null at the initial state of 1673 K, and the development of stress distribution was calculated at a cooling rate of 20 K/min while creep was considered within the temperature range of 1673–1173 K. Creep was assumed to be negligible for SiC and SiC/SiC substrate. Since creep was not taken into account for all materials (i.e., the analysis is independent of time) below 1173 K, the system was cooled immediately from 1173 to 303 K.

As we consider a wide temperature range, we adopted temperature-dependent material properties; Young's modulus, *E*, Poisson's ratio, ν, CTE (reference temperature: 1673 K), α', and creep properties (creep coefficient, *C*, and creep index, *n*) listed in Tables 2–8, respectively. The elastic properties (i.e., *E* and ν) of the YbMS anisotropic homogeneous layer and the SiC/SiC substrate were regarded as anisotropic from the microstructure and the fiber orientation, respectively. Also, *E* and ν of the mullite layer were regarded as anisotropic on the basis of the experimental result of the in-plane and out-of-plane loading tests on the deposited mullite layer (along the *x* and *y* axes in Figure 2, respectively). The other layers were dealt with as isotropic materials. CTE and the creep properties of each layer and substrate were assumed to be isotropic. The temperature-dependent material properties were calculated from experimental measurements and literature (see Appendix A).


**Table 2.** Young's moduli (*E*) of isotropic elastic layers.


**Table 3.** Poisson's ratios (ν) of isotropic elastic layers.

<sup>1</sup> No temperature dependence is assumed.


**Table 4.** Elastic properties (*E*, ν, and *G*) of SiC/SiC substrate.

<sup>1</sup> No temperature dependence is assumed.

**Table 5.** Elastic properties (*E*, ν*,* and *G*) of mullite layer.



**Table 5.** *Cont.*

<sup>1</sup> No temperature dependence is assumed.

**Table 6.** Elastic properties (*E*, ν, and *G*) of YbMS anisotropic homogenous layer.


<sup>1</sup> No temperature dependence is assumed.

**Table 7.** Coefficients of thermal expansion (CTEs) at reference temperature of 1673 K (α') of substrate and coating layers.



**Table 8.** Creep properties (*C* and *n*) of coating layers.

<sup>1</sup> The unit of *C* depends on the corresponding creep index *n*. <sup>2</sup> No temperature dependence is assumed.

#### 2.2.3. Interface Crack Introduction Analysis

ERR owing to the thermal stress, <sup>G</sup>th, is defined as the reduction in strain energy, <sup>Π</sup>, stored in a structure per increased crack area. When the crack area, *A*, increases to *A* + d*A*, ERR is evaluated as follows:

$$\mathcal{G}^{\text{th}} = -\frac{\text{d}\Pi}{\text{d}A}.\tag{1}$$

Thus, it is necessary to obtain the strain energies stored in T/EBC with different crack lengths to evaluate <sup>G</sup>th init- <sup>=</sup> <sup>G</sup>th *A*→0 and <sup>G</sup>th prop. We carried out FEM analyses for T/EBC models with various crack length under the stress condition obtained from the thermal stress analysis, and the strain energy in T/EBC was evaluated as a function of crack length.

The geometry, the boundary conditions and the mesh sizes of the model in this analysis were set in the same way as in the thermal stress analysis. The initial stress condition was derived from the stress distribution of the uncracked T/EBC model obtained by the thermal stress analysis. The temperature was kept constant at 303 K (temperature after the cooling process) during the analyses, and thus we used the material properties *E* and ν at 303 K shown in Tables 2–6. Interface cracks were introduced from the right edge of the simulation model shown in Figure 2. According to a preliminary FEM analysis, the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interfaces were chosen as the objective interfaces to evaluate <sup>G</sup>th init and <sup>G</sup>th prop (the detail of the preliminary analysis is found in Appendix B).

#### 2.2.4. Evaluation of ERR for Interface Crack Initiation and Propagation

We evaluated <sup>G</sup>th init and <sup>G</sup>th prop from the strain energy in T/EBC with different crack lengths obtained by the FEM analyses in Section 2.2.3 in the following ways.

#### Evaluation of <sup>G</sup>th init

<sup>G</sup>th init is defined as ERR in the limit when *A* approaches 0, as follows:

$$\mathcal{G}\_{\text{init}}^{\text{th}} = -\lim\_{A \to 0} \frac{\text{d}\Pi}{\text{d}A}.\tag{2}$$

Note that <sup>G</sup>th init is equal to the first-order coefficient of the Taylor expansion of Π(*A*) around *<sup>A</sup>* <sup>=</sup> <sup>0</sup> with the inverted sign. Therefore, we approximated Π(*A*) as a polynomial from the FEM analyses for the T/EBC models with several different crack lengths, and <sup>G</sup>th init was evaluated as the first-order coefficient of the polynomial with the inverted sign. Two cases are possible according to the sign of <sup>G</sup>th init: for <sup>G</sup>th init <sup>&</sup>gt; 0, crack initiation occurs if the interface fracture toughness, <sup>Γ</sup>, is smaller than <sup>G</sup>th init; for <sup>G</sup>th init < 0, no crack initiation is expected. From the result of the preliminary analysis (Appendix B), we examined crack length *<sup>a</sup>* <sup>=</sup> 2.5, 5.0 and 7.5 <sup>μ</sup>m to evaluate <sup>G</sup>th init at the objective interfaces.

#### Evaluation of <sup>G</sup>th prop

ERR when a crack propagates from *A* to *A* + Δ*A* is given by numerical differentiation as follows:

$$\mathcal{G}\_{\text{prop}}^{\text{th}}\left(A + \frac{\Delta A}{2}\right) = -\frac{\text{d}\Pi}{\text{d}A}\Big|\_{A + \Delta A/2} \approx -\frac{\Pi(A + \Delta A) - \Pi(A)}{\Delta A}.\tag{3}$$

From the result of the preliminary analysis (Appendix B), the *a* for evaluation of <sup>G</sup>th prop was determined as 2.5, 5.0, 7.5, 10, 50, 100, 250, ..., 3500 μm in increments of 250 μm between 250 and 3500 μm.

#### *2.3. Interface Fracture Test*

Interfacial toughness was measured by an interface fracture test. The test method was a modification of the one designed for EBCs on SiC/SiC composites having a weak inter-laminar strength [16]. Figure 3 shows schematic of the test set-up.

**Figure 3.** Schematics of experimental set-up of interface fracture test: (**a**) Cross-sectional side view; (**b**) Front view.

A block of specimen with a height, *L*, of 4 mm and a width, *w*, of 3 mm was prepared using a diamond blade saw in such a way that one of the fiber directions in the woven fabric of the SiC/SiC substrate was oriented to the longitudinal side of the specimen. The surface was as cut without polishing to avoid the coating delamination during the processing. The cut surface was smooth enough for us to distinguish the interface in the cross-sectional observation by conventional optical microscopy.

The specimen was set in a steel jig like a horizontal die, the upper part of which was adjustable. The specimen was put in the die and slid horizontally from the back side by a punch as indicated by the dashed arrow 1 in Figure 3 so that only the coating protruded from the die while the substrate was buried. Then the upper part was adjusted in the vertical direction (the dashed arrow 2 in Figure 3) to fix the specimen. The jig with the fixed specimen was placed under a universal testing machine (load cell capacity: 500 N, EZ-LX, Shimadzu Corporation, Kyoto, Japan) having a steel pushing plate with a thickness of 0.5 mm and a width of 5 mm.

The edge of the coating was compressively loaded by the plate with a constant crosshead speed of 0.024 mm/min till the coating was delaminated from the substrate; thus, the mode II-rich loading at the interface was achieved with little loading to the substrate in the inter-laminar shear direction. The shear stress ought not to concentrate at only an interface between SiC/SiC substrate and T/EBC

(SiC/SiAlON interface); i.e., it occurs along some interfaces in T/EBC. The pushing load, *P*, and the crosshead displacement, *u*, were measured during the test. The coating near the loading point was observed in the cross-sectional direction of the specimen by optical microscopy to determine when a crack was initiated at the loaded edge and the corresponding load, *P*init. The number of the tests was six.

The energy release rate, Ginit, associated with the initiation of an interface crack at the loaded edge under an applied load of *P*init is given by

$$\mathcal{G}\_{\rm init} = \frac{\left(P\_{\rm init}\right)^2}{2w^2h\_{\rm coat}} \left\{ \frac{1}{E\_{\rm coat}'} - \frac{1}{E\_{\rm sub}'} \left( \frac{1}{D} + \frac{\left(h\_{\rm cost}/2 + h\_{\rm sub} - \delta\right)^2}{I\left(h\_{\rm cost}\right)^2} \right) \right\}. \tag{4}$$

Definition of parameters and details of deriving Equation (4) are shown in Appendix C. The Ginit is defined as fracture toughness for the interface crack initiation, Γinit. Note that, Γinit is a nominal fracture toughness value of interface crack initiation including the effect of thermal residual stress.

After the tests, four of the specimens were mounted in resin and their cross-sections of the delaminated coating and substrate were observed; the cross-section was polished to mirror finish. SEM observation (VE–7800, Keyence corporation, Osaka, Japan) and energy dispersive X-ray (EDX) analysis (EX74120 attached to field emission SEM (JSM6500F, JEOL Ltd., Tokyo, Japan)) of the cross section were done to identify the fractured interface. The remaining two specimens were served for the analyses of fracture surfaces; optical and SEM observations and EDX analysis were done.

#### **3. Results of T**/**EBC Fabrication**

#### *3.1. Layer Thickness Condition for Preventing the Objective Interface Crack Initiation*

Figure 4a–c show the relationships between the SiAlON layer thickness, *h*SA, and <sup>G</sup>th init at the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interfaces, respectively, with various mullite layer thicknesses, *h*Mu. From Figure 4a, the ERR for interface crack initiation at the SiC/SiAlON interface, <sup>G</sup>th init, SC/SA, is nearly zero for any combination of layer thicknesses. Thus, a crack is unlikely to initiate at the SiC/SiAlON interface during the cooling process after deposition if thicknesses of the SiAlON and mullite layers are both within the range of 5–25 μm.

The ERR for interface crack initiation at the SiAlON/mullite interface, <sup>G</sup>th init, SA/Mu, depends on the thicknesses of *h*SA and *h*Mu as shown in Figure 4b, which is found to be expressed as the following equation within the examined range of *h*SA and *h*Mu (5–25 μm):

$$\mathcal{G}\_{\text{init,SA/Mu}}^{\text{th}}(h\_{\text{SA}}, h\_{\text{Mu}}) = \sum\_{k=0}^{2} \sum\_{l=0}^{2} c\_{kl} (h\_{\text{SA}})^k (h\_{\text{Mu}})^l \tag{5}$$

with the coefficients, *ckl*, listed in Table 9. Figure <sup>5</sup> shows the contour chart of <sup>G</sup>th init, SA/Mu as a function of *<sup>h</sup>*SA and *<sup>h</sup>*Mu shown in Equation (5). The level of <sup>G</sup>th init, SA/Mu is indicated by the depth of color (violet); a deeper-colored region corresponds to a lower <sup>G</sup>th init, SA/Mu and thus a better combination of layer thicknesses. This result suggests that one of the SiAlON or mullite layers must be thick and the other be thin in order to decrease <sup>G</sup>th init, SA/Mu and hinder crack initiation at the SiAlON/mullite interface.

**Figure 4.** Relationship between SiAlON layer thickness and ERR for interface crack initiation due to thermal stress at (**a**) SiC/SiAlON, (**b**) SiAlON/mullite and (**c**) mullite/YbDS interfaces with various mullite layer thicknesses.



<sup>1</sup> Unit of *ckl*: μm(*k*+*l*).

**Figure 5.** Contour chart of <sup>G</sup>th init,SA/Mu(*h*SA, *<sup>h</sup>*Mu) interpolated with Equation (5).

As shown in Figure 4c, it is found that a thin mullite layer (*h*Mu = 5 μm) results in a rise in the ERR for interface crack initiation at the mullite/YbDS interface, <sup>G</sup>th init, Mu/YD, with increasing *h*SA, while a thicker mullite layer (*h*Mu = 15 and 25 <sup>μ</sup>m) reduces <sup>G</sup>th init, Mu/YD to almost zero regardless of *h*SA. Therefore, a thick mullite layer is necessary to decrease <sup>G</sup>th init, Mu/YD, i.e., to suppress crack initiation at the mullite/YbDS interface.

From the above results, it is found that a thin SiAlON layer and a thick mullite layer are required for reducing <sup>G</sup>th init and preventing crack initiation at the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interfaces, where strong stress concentration is expected during the cooling process in fabrication of T/EBC as shown in Figure 1.

#### *3.2. Phase Structure of T*/*EBC*

In Section 3.1, we found that thinning the SiAlON layer together with thickening the mullite layer was effective in suppressing interface crack formation during cooling where a large thermal stress is expected, i.e., the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interfaces. In this study, therefore, the thickness of the SiAlON layer and the mullite layer were determined to be 5 and 20 μm, respectively, and T/EBC was prepared by the dual EB-PVD process described in Section 2.1. We confirmed no delamination at all the interfaces of the fabricated T/EBC. Cross-sectional SEM images and elemental mapping images of T/EBC were shown in Figure 6a,b, respectively, where we found the thicknesses of constituent layers to be; SiAlON: about 5 μm, mullite: about 20 μm, Yb-silicate gradient composition layer: about 100 μm, and porous segmented YbMS layer: about 200 μm. That is, T/EBC was successfully fabricated nearly as designed.

**Figure 6.** (**a**) Cross-sectional SEM images of T/EBC and (**b**) elemental mapping images focused on SiC/SiAlON/mullite layers.

#### **4. Behavior of Interface Crack Initiation and Propagation in T**/**EBC**

To investigate the behavior of crack initiation and propagation during the cooling process in fabrication of T/EBC shown in Section 3.2, calculated ERRs (Gth init, <sup>G</sup>th prop) must be compared with interface fracture toughness obtained by experiment. Thus, we carried out the interface fracture toughness tests (Section 2.3) to obtain <sup>Γ</sup>, which is to be compared with <sup>G</sup>th init and <sup>G</sup>th prop with *h*SA = 5 μm and *h*Mu = 20 μm.

#### *4.1.* G*th init and* <sup>G</sup>*th prop for Objective Interface Cracks in Fabricated T*/*EBC*

Table <sup>10</sup> lists <sup>G</sup>th init of the objective interface cracks for the T/EBC model with *h*SA = 5 μm and *h*Mu = 20 <sup>μ</sup>m. From this result, the SiAlON/mullite interface is found to possess the highest <sup>G</sup>th init of the objective interfaces.

**Table 10.** <sup>G</sup>th init for T/EBC model with *h*SA = 5 μm and *h*Mu = 20 μm.


Figure <sup>7</sup> shows the relationships between <sup>G</sup>th prop and *a* at the objective interfaces in the T/EBC with *h*SA = 5 <sup>μ</sup>m and *h*Mu = 20 <sup>μ</sup>m. In the <sup>G</sup>th prop–*a* relationships at all the objective interfaces, we observe three stages with increasing *a*: (A) <sup>G</sup>th prop shows a relatively rapid increase; (B) <sup>G</sup>th prop decreases; (C) <sup>G</sup>th prop increases again.

**Figure 7.** *Cont.*

**Figure 7.** Relationships between ERR for interface crack propagation due to the thermal stress and crack length at objective interfaces in T/EBC with *h*SA = 5 μm and *h*Mu = 20 μm; (**a**) SiC/SiAlON interface: (**a–1**) 0 ≤ *a* ≤ 3500 μm and (**a–2**) 0 ≤ a ≤ 100 μm, (**b**) SiAlON/mullite interface: (**b–1**) 0 ≤ a ≤ 3500 μm and (**b–2**) 0 ≤ a ≤ 200 μm and (**c**) mullite/YbDS interface: (**c–1**) 0 ≤ a ≤ 3500 μm and (**c–2**) 0 ≤ a ≤ 100 μm.

The reason for the decrease in <sup>G</sup>th prop at Stage B is explained by the distribution of out-of-plane thermal stress in the vicinity of the interface edge. Figure 8a shows the σ*<sup>y</sup>* distributions in the coating layers along the *y* axis, which are obtained at the positions *rx* distant from the right interface edge in the *x* direction (*rx* = 0, 5, 10, 20, 30 and 40 μm). The highest level of σ*<sup>y</sup>* was observed at the interface edge (*rx* = 0 μm), and the σ*<sup>y</sup>* distribution reduces immediately to a negligible level within *rx* = 30 μm. Thus, the effect of σ*<sup>y</sup>* on the coating layers is significant but extremely localized in the vicinity of interface edge, while that is negligible inside the simulation model. Note that the effect of <sup>σ</sup>*<sup>x</sup>* on <sup>G</sup>th prop is marginal for sufficiently small *a* because we find σ*<sup>x</sup>* ~ 0 at the interface edge (*rx* = 0 μm) directly from the balance of stress component in the *x* direction. Figure 8b shows the σ*<sup>x</sup>* distributions in the coating layers along the *y* axis at various *rx*. These results indicate that <sup>σ</sup>*<sup>y</sup>* has a dominating effect on <sup>G</sup>th prop for very short cracks. The mechanism of the second rise in <sup>G</sup>th prop at Stage C is attributed to the increase in the level of σ*<sup>x</sup>* inside the simulation model as shown in Figure 8b. Therefore, the three-stage behavior in the <sup>G</sup>th prop–*<sup>a</sup>* relationship is explained as follows: in the beginning, <sup>G</sup>th prop rises as a crack propagates owing to the strong <sup>σ</sup>*<sup>y</sup>* near the interface edge (Stage A); then <sup>G</sup>th prop decreases because of a steep drop in σ*<sup>y</sup>* inside the model (Stage B); for a sufficiently long crack, <sup>G</sup>th prop is mainly affected by σ*<sup>x</sup>* and rises again with increasing σ*<sup>x</sup>* (Stage C).

**Figure 8.** *Cont.*

*x*

SiAlON SiC

Mullite

SiC/SiC sub.

*rx*

0

*ry*

**Figure 8.** (**a**) σ*y* distributions, and (**b**) σ*x* distributions in coating layers along the *y* axis at the positions *rx* distant from right SiC/substrate interface edge in the *x* direction.

#### *4.2. Interface Fracture Toughness of T*/*EBC*

Figure 9 shows a typical *P*–*u* curve during the tests. By the cross-sectional optical microscopy, the onset of interface fracture near the loaded edge was identified at Point A in *P*–*u* curve before the maximum load. The load at Point A is defined as *P*init. The crack then propagated down through the interface to the midway till the load reached the maximum; finally it grew unstably to the complete delamination at the maximum load. Buckling or compressive fracture of the coating was not observed during the test. The measured *P*init values for respective tests are listed in Table 11.

**Figure 9.** Typical load-crosshead displacement curve during test (Specimen No. 2).

**Table 11.** Crack initiation load (*P*init) obtained by interface fracture test.


Examples of the cross-sectional SEM images of the specimens after the tests were shown in Figure 10. At the loaded edge where the interface crack started to propagate, fracture occurred below or above the SiAlON layer, i.e., either at the interface between the SiC and SiAlON layers (Figure 10a) or at the interface between the SiAlON and mullite layers (Figure 10b). The crack path was along either of the two interfaces, and it was switched to each other by crack kinking across the SiAlON layer. However, crack propagation within a layer, such as a crack passing through the inside of layer parallel to the coating plane, was not observed.

**Figure 10.** Examples of cross-sectional SEM images of specimens after interface fracture test; Fracture at (**a**) SiC/SiAlON interface and (**b**) SiAlON/mullite interface.

Figure 11a shows a fracture surface near the loaded edge after the tests observed by optical microscopy. The fracture surface was classified into three according to their colors: the largest black area denoted by Region A, the gray area by Region B, and the small white area by Region C as illustrated in Figure 11a. EDX maps of a selected area in the fracture surface (indicated by a dashed rectangle in Figure 11a) are shown in Figure 11b. Mapping of element distributions shows the concentration of Si in Region A. In Region B the existence of Si, Al and N was detected. Region C was characterized by Yb concentration. Apparently Al also looks rich in Region C, but it is due to misdetection because the Mα line of Yb is very close to Kα line of Al which was used for the mapping of Al distribution. These results suggest that Region A with the largest area corresponds to the fracture surface created by the interface fracture between the SiC and SiAlON layers observed in Figure 10a; Region B showing the second largest area corresponds to the surface created by the SiAlON/mullite interface fracture (Figure 10b). We could not observe the fractured interface corresponding to Region C in the cross-sectional SEM observation of the T/EBC and substrate (Figure 10) because of the limited area of Region C. To determine the exact location of the fracture surface in Region C, we cut after the tests one of the substrate in the cross-section including Region C. The result is shown in Figure 12 suggesting that the fracture occurred in the Yb-silicate gradient composition layer.

To calculate the line fractions of Regions A, B and C along the loaded edge, the fracture surfaces near the edge were analyzed using image processing software (Image J 1.48v, National Institute of Health, Bethesda, MD, USA) (Table 12). As shown in Table 12, fracture surface was dominantly composed of Regions A and B; the ratio of Region C was relatively small. Thus the crack initiation at the edge was supposed to occur primarily either at the interface between the SiC and SiAlON layers or at the interface between the SiAlON and mullite layers.

**Figure 11.** (**a**) Top view of fracture surface of specimen after interface fracture test and (**b**) EDX maps of selected area in fracture surface indicated by a dashed rectangle in (**a**).

**Figure 12.** Cross-sectional SEM observation including Region C.

**Table 12.** Line fractions of Regions A, B and C at loaded edge.


The energy release rate for the edge crack initiation along the SiC/SiAlON interface under *P*init (in Table 11) was calculated from Equations (A14)–(A18) and (4), regarding the SiC/SiC plate and the SiC layer as "substrate" and the rest of the layers as "coating" in Equation (A14); the mean value of 6.4 J/m<sup>2</sup> was obtained. When the SiAlON layer is counted as part of the substrate in addition to the SiC/SiC and SiC, the mean energy release rate for the crack along the SiAlON/mullite interface becomes 6.5 J/m2. These values of energy release rates are very close to each other because the thickness of the SiAlON layer lying between the two interfaces is quite small compared to the total thickness of the system. Under the assumption that the edge crack was initiated simultaneously through the edge (i.e., in the direction of depth) at *P*init regardless of the fractured interfaces, we can estimate that the SiC/SiAlON interface fracture and the SiAlON/mullite interface fracture both occur with significant proportions at almost the same energy release rates. This suggests that the interface toughness of these two interfaces are very close; both the interfaces should have an interface toughness of ~6.4 J/m2. The mean energy release rate at *P*init when the edge crack was initiated along the mullite/YbDS interface, which was supposed to have large thermal stress according to the FEM analysis, was calculated to be 6.8 J/m2. Since no fracture occurred at this interface at *P*init, we can expect that the toughness of the interface should be larger than 6.8 J/m2. The abovementioned results indicate that the approximate Γinit for the SiC/SiAlON and SiAlON/mullite interfaces are 6.4 J/m<sup>2</sup> and the lower limit of Γinit for the mullite/YbDS interface is 6.8 J/m2.

Table 13 shows comparison between <sup>G</sup>th init of the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interface cracks and <sup>Γ</sup>init of the corresponding interfaces. <sup>G</sup>th inits of every objective interface crack are smaller than Γinit. This result suggests that the effect of thermal residual stress on the intrinsic fracture toughness for interface crack initiation is sufficiently small. Therefore, in this study, we regard the nominal value as the fracture toughness of interface crack initiation. Table 13 also suggests that the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interface cracks are not likely to initiate by cooling in the T/EBC fabrication process.

**Table 13.** Comparison between <sup>G</sup>th init of target interface cracks and fracture toughness for interface crack initiation (Γinit) of corresponding interfaces.


<sup>1</sup> Approximate value. <sup>2</sup> Lower limit value.

Figure 13a–c show comparisons between <sup>G</sup>th prop and fracture toughness for interface crack propagation, Γprop, at the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interfaces. Here, we assumed that fracture surfaces formed by interface crack initiation and propagation were the same and regarded <sup>Γ</sup>prop as equal to <sup>Γ</sup>init. As shown in Figure 13a,c, <sup>G</sup>th prop of the SiC/SiAlON and mullite/YbDS interface cracks are smaller than Γprop of the corresponding interfaces. The results suggest that these interface cracks are not likely to propagate regardless of initial interface crack length. Figure 13b suggests that, if the T/EBC has an initial length of the SiAlON/mullite interface crack of about 1200 μm, the crack should propagate instantly because <sup>G</sup>th prop exceeds Γprop at the crack length of 1200 μm. However, delamination along the SiAlON/mullite interface was not observed in the experiment (see Section 3.2). These results suggest that there was no initial crack exceeding the length threshold at the SiAlON/mullite interface.

The results shown in this section reveal that the T/EBC with the proposed layer thicknesses can be fabricated without delamination along interfaces by cooling in the fabrication process. This is confirmed by the result that no delamination along interface is observed as described in Section 3.2.

#### *4.3. Future Prospects*

In this study, we mainly focused on the mechanical reliability of T/EBC during the cooling process in fabrication. From the practical viewpoint, however, it is also essential to evaluate the mechanical reliability under an operating condition. In general, an EBC is exposed to a thermal cycle condition with a high humidity in operation, which can induce a microstructural change and chemical transformation of the coating layers through a reaction with heated oxygen and water vapor. Consequently, changes in the material properties of each layer due to microstructural change and chemical transformation are crucial to the mechanical state in T/EBC. Thus, it is necessary to evaluate a time dependent ERR by the FEM analysis where changes in material properties over time are incorporated. Furthermore, the interface fracture toughness changes presumably over time owing to microstructural change and chemical transformation, which should also be taken into account for a reliable design of T/EBC. Nevertheless, it is beyond the scope of this paper and will be considered in our future work.

**Figure 13.** Comparison between ERR for crack propagation due to thermal stress and fracture toughness for crack propagation at (**a**) SiC/SiAlON interface, (**b**) SiAlON/mullite interface and (**c**) mullite/YbDS interface.

#### **5. Conclusions**

To assess the mechanical reliability of T/EBC (SiC/SiAlON/mullite/Yb-silicate gradient composition layer/YbMS with porous segment structure) against interface crack initiation and propagation due to thermal stress induced in fabricated process, we carried out FEM analysis to evaluate ERR for interface cracks and performed experiment to obtain interface fracture toughness. Our FEM calculations revealed that <sup>G</sup>th init of the objective interfaces decreases by making the SiAlON layer thinner and the mullite layer thicker. We fabricated the T/EBC with layer thicknesses within the proposed range (*h*SA = 5 μm and *h*Mu = 20 μm) and confirmed no delamination along the interfaces. In the interface fracture test for the fabricated T/EBC, fracture surfaces were found at the SiC/SiAlON and SiAlON/mullite interfaces. We estimated the approximate fracture toughness for the SiC/SiAlON and SiAlON/mullite interfaces and minimum limit of fracture toughness for the mullite/YbDS interface. Comparison between the fracture toughness by the experiment and calculated <sup>G</sup>th init and <sup>G</sup>th prop indicated that the fabricated T/EBC possesses sufficient mechanical reliability against interface cracks.

Of course, the present study considers only the stress state just after the fabrication process and does not assure the mechanical reliability in operation under thermal cycle with a high humidity, which is regarded as our future work.

**Author Contributions:** Conceptualization, E.K. and Y.U.; Fabrication of T/EBC, N.Y. and T.Y.; Measurement of material properties, N.Y., T.Y., H.K. and T.A.; FEM analysis, E.K. and A.K.; Interface fracture test, H.K.; Writing–original draft preparation for fabrication process and result, N.Y. and T.Y.; Writing–original draft preparation for FEM analysis and result, E.K. and A.K.; Writing–original draft preparation for interface fracture test procedure and result, H.K.; Writing–review and editing, S.K. and Y.U.; Supervision, H.K., T.A., S.K. and Y.U.; Project administration, S.K.; Funding acquisition, H.K., T.A., S.K. and Y.U.

**Funding:** This work was supported by Council for Science, Technology and Innovation (CSTI), Cabinet Office, Japan, through Cross-ministerial Strategic Innovation Promotion Program (SIP), "Structural Materials for Innovation" (Funding agency: JST, Japan Science and Technology Agency).

**Acknowledgments:** The authors would like to thank Kaori Sakata, Shoko Taniguchi and Takahumi Kawano for technical assistance.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Procedure for Calculating the Temperature-Dependent Material Properties**

*Appendix A.1. Young's Modulus, Poisson's Ratio and Shear Modulus*

In this study, Young's moduli *E* at 303–1673 K of the SiC and SiAlON layers were considered to be identical to that of their bulk, *E*bulk. *E* at the same temperature range of the mullite, YbDS and YbMS layers was estimated from the measured values of coating specimens deposited by the dual EB-PVD process with the same conditions described in Section 2.1, *E*film. The outline of the estimation method of *E* of each layer is shown below;

Step 1: Derivation of an approximation formula of relation between temperature, *T*, and *E*bulk of measured values or literature values.

Step 2: Estimation of *E*bulk at the temperature range from 303 to 1673 K based on the approximation formula derived in Step 1.

Step 3: Measurement of *E*film at 303 K of a film/substrate layered specimen.

Step 4: Estimation of temperature dependence of the *E*, *E*(*T*), of the coating layers using the temperature dependence of *E*bulk and the ratio of *E*film(303) to *E*bulk(303) shown in Equation (A1);

$$E(T) = E\_{\text{bulk}}(T) \frac{E\_{\text{film}}(\text{303})}{E\_{\text{bulk}}(\text{303})}.\tag{A1}$$

*E*s of the SiC and SiAlON layers were estimated through Steps 1 and 2. The values of *E* of the mullite, YbDS and YbMS layers were estimated through Steps 1–4. Tables A1 and A2 summarize the estimation procedure for each layer at Step 1 and Steps 2–3, respectively. For the Yb-silicate gradient composition layer, *E* was obtained from Young's moduli of the YbMS and YbDS layers by the following Equation [17]:

$$E\_{\beta:1-\beta} = (E\_{\rm YM})^{\beta}(E\_{\rm ND})^{1-\beta}, \beta = 0, 0.2, 0.5, 0.8, \text{ and } 1. \tag{A2}$$

Here, β and 1 − β denote volume content of YbMS and YbDS, respectively.


**Table A1.** Estimation procedure of *E* of SiC, SiAlON, mullite, YbDS and YbMS layers (Step 1).

<sup>1</sup> CVD-SiC wafer, ADMAP Inc.. <sup>2</sup> Raw powder: BSI3–001B, AG Materials Inc.. Sintering condition: 2023 K for two hours, applied pressure 40MPa, *P*N2 = 0.6MPa. <sup>3</sup> Raw powder: KM-101, KCM Corporation. Sintering condition: 1573 K for fifty hours in air → 2023 K for five hours in air. <sup>4</sup> Raw powder: self-build by ultrasonic spray pyrolysis. Sintering condition: 1673 K for two hours in air → Jet mill grinding → 1773 K for five hours in air. <sup>5</sup> Resonant method with cantilever bending and torsion mechanism. Dimensions of a specimen are 2 mm × 10 mm × 60 mm.


**Table A2.** Detail of estimation procedure of *E* of SiC, SiAlON, mullite, YbDS and YbMS layers (Steps 2 and 3).

<sup>1</sup> Reference [19], Optimization of *E*0, *F*, and *T*0. <sup>2</sup> Oliver and Pharr method.

The values of ν of the SiC, SiAlON and YbDS layers were considered to be identical to those of their bulk. The outline of the estimation procedure of ν of each layer at Step 1 is shown in Table A3. Step 1: Measurement of ν of the bulk, νbulk in some temperature range.

Step 2: Obtainment of ν in the temperature range from 303 to 1673 K by the assumption that νbulk outside the temperature range in Step 1 was identical to that measured at adjacent temperature.



ν*ij* of the anisotropic mullite layer in the temperature range from 303 to 1673 K was estimated using ν*ij* of the bulk specimen described in Table A1 and *E*film of the mullite layer derived from the procedure shown in Table A2. ν of the YbMS used in this estimation was literature value [9] and temperature-independent. For the Yb-silicate compositional gradient layers, ν was obtained from Young's moduli and bulk moduli, *B*, of YbMS and YbDS by the following equations:

$$\nu\_{\beta:1-\beta} = \frac{1}{2} \left( 1 - \frac{E\_{\beta:1-\beta}}{3B\_{\beta:1-\beta}} \right) \cdot B\_{\beta:1-\beta} = (B\_{\text{YM}})^{\beta} (B\_{\text{YD}})^{1-\beta} \,. \tag{A3}$$

The YbMS anisotropic homogeneous layer is characterized by anisotropic Young's modulus, *E*ani *<sup>i</sup>* , Poisson's ratio, νani *ij* , and shear modulus, *<sup>G</sup>*ani *ij* (*i*, *j* = *x*, *y*, *z*; *i j*). A method of calculations for these parameters will be discussed in Reference [14]. Shear modulus of mullite layer was estimated by Equation (A4).

$$G\_{ij} = \frac{E\_i E\_j}{E\_i + E\_j + 2E\_j \nu\_{ij}}.\tag{A4}$$

#### *Appendix A.2. CTE (Reference Temperature: 1673 K)*

α was calculated from CTE at the reference temperature of 303 K (α; evaluated as a function of temperature from experimental measurements) as follows:

$$\alpha'(T) = \frac{a(T) \times (T - 303) - \alpha(1673) \times (1673 - 303)}{(T - 1673) \times (a(1673) \times (1673 - 303) + 1)}.\tag{A5}$$

α of the SiC and SiAlON layers was considered to be identical to that of their bulk. The values of α of the sintered bodies and the wafer in the temperature range from 573 to 1673 K were measured based on JIS R1618. Then, α of the sintered bodies and the wafer at 323–1673 K was estimated by extrapolation of the relationship between α and *T*. α of the mullite, YbDS and YbMS layers at given temperature was estimated using measured values of film/substrate layered specimens deposited with the same conditions as in Section 2.1. The outline of the calculation procedure of α of each layer is described below;

Step 1: Derivation of relationship between strain, ε, and *T* for each film/substrate layered specimen (described in Table A2), based on digital image correlation method.

Step 2: Approximation of the ε–*T* relationship derived in Step 1 by Equation (A6).

$$
\varepsilon = N\_1 (T - 303) + N\_2 (T - 303)^2. \tag{A6}
$$

Step 3: Calculation of α in the temperature range from 323 to 1673 K using the approximation formula derived in Step 2.

The value of α of the Yb-silicate gradient composition layers was obtained by the law of mixture as follows [20]:

$$\alpha\_{\beta:1-\beta} = \frac{\alpha\_{\text{YM}}\beta B\_{\text{YM}} + \alpha\_{\text{YD}}(1-\beta)B\_{\text{YD}}}{\beta B\_{\text{YM}} + (1-\beta)B\_{\text{YD}}}. \tag{A7}$$

In this study, α of the YbMS anisotropic homogeneous layer was regarded as identical to that of the YbMS dense layer.

#### *Appendix A.3. Creep Properties*

In this study, we assumed that the creep property of each layer follows the strain-hardening law expressed as .

$$
\dot{\varepsilon}^{\mathcal{C}} = \mathbb{C}\sigma^{u},
\tag{A8}
$$

where . ε <sup>c</sup> denotes creep strain rate.

As described in Section 2.2.2, we assumed that creep was negligible for the SiC layer. Creep properties of the SiAlON, mullite, YbDS and YbMS layers were estimated by the method described below;

Step 1: Estimation of *n* at a given temperature and the average of *n* within the temperature range from the graph of the true stress and the strain rate derived from literature values and/or measured values. Noted that the average of *n* was used in the FEM analysis.

Step 2: Estimation of *C* at a given temperature using the average of *n*, true stress and strain rate derived in Step 1.

Step 3: Estimation of *C* in the temperature range from 1173 to 1673 K using *C*∗ and *Q*/*R* in Equation (A9) based on *C* versus inverse of *T* plot derived in Step 2.

$$\mathcal{C} = \mathcal{C}^\* \exp\left(-\frac{\mathcal{Q}}{RT}\right). \tag{A9}$$

The estimation method of the creep properties of each layer at Step 1 is shown in Table A4. For the Yb-silicate gradient composition layers, *C* and *n* were obtained as follows:

$$\mathbf{C}\_{\beta:1-\beta} = (\mathbf{C}\_{\text{YM}})^{\beta} (\mathbf{C}\_{\text{YD}})^{1-\beta},\tag{A10}$$

$$n\_{\beta:1-\beta} = (n\_{\rm YM})^{\beta} (n\_{\rm ND})^{1-\beta}. \tag{A11}$$

The creep property of the YbMS anisotropic homogeneous layer was regarded as identical to that of the YbMS dense layer.


**Table A4.** Method of creep property estimation of the SiAlON, mullite, YbDS and YbMS layers (Step 1).

<sup>1</sup> Raw powder: self-build by ultrasonic spray pyrolysis. Sintering condition: 1873 K for two hours in air → Jet mill grinding → 1823 K for 0.5 hours applied pressure 50 MPa in Ar, 1873 K for five hours in air. <sup>2</sup> Dimensions of a specimen are 3 mm × 2 mm × 2 mm.

#### **Appendix B. Preliminary FEM Analyses: Determination of Objective Interfaces and Crack Lengths**

As was explained in Introduction, an interface crack in the T/EBC is likely to occur at interfaces with strong thermal stress concentration due to the difference in the CTEs from the SiC/SiC substrate. Here we determined the objective interfaces and crack lengths from the stress distribution obtained by a preliminary FEM analysis. The thermal stress analysis was carried out in the way shown in Section 2.2.2 for the uncracked T/EBC model with the SiAlON and mullite layers both being 25 μm thick.

Figure A1 shows the distribution of σ*<sup>x</sup>* near the interface edge in the constituent layers of the T/EBC after the cooling process. Thermal stress is found to be concentrated particularly at the SiC/SiAlON, SiAlON/mullite and mullite/YbDS interfaces. Thus, we focused on crack initiation and propagation at these objective interfaces in this study.

From the balance of the stress component in the *x* direction, we deduce σ*<sup>x</sup>* ~ 0 at the interface edge and thus a dominating effect of <sup>σ</sup>*<sup>y</sup>* on <sup>G</sup>th init. Therefore, we determined the crack lengths examined to evaluate <sup>G</sup>th init, based on the σ*<sup>y</sup>* distribution in the vicinity of the interface edge. Figure A2 shows the σ*<sup>y</sup>* distribution in the coating layers along the *y* axis at the positions *rx* distant from the interface edge (*rx* = 0, 5, 10, 60, 100, 300, 500 and 1000 μm). The highest level of σ*<sup>y</sup>* was observed at the interface edge (*rx* = 0 μm) while the σ*<sup>y</sup>* distribution reduces to a negligible level inside the model. On the basis of the results, we examined the crack lengths of *a* = 2.5, 5.0 and 7.5 <sup>μ</sup>m to evaluate <sup>G</sup>th init. For evaluation of <sup>G</sup>th prop, the crack length *a* was varied as *a* = 2.5, 5.0, 7.5, 10, 50, 100, 250, ..., 3500 μm in increments of 250 μm between 250 and 3500 μm.

**Figure A1.** Thermal stress distribution near right interface edge of T/EBC deposited on SiC/SiC substrate without interface cracks.

**Figure A2.** σ*y* distributions in coating layers along the *y* axis at the positions *rx* distant from right SiC/substrate interface edge in *x* direction.

#### **Appendix C. Derivation Process of Energy Release Rate under Applying Mechanical Loading**

Figure A3 shows the loads and moments applied to the specimen. *Pb* and *Mb* (*b* = 1, 2) are the load and moment from surrounding jigs, respectively, and *P* is the applied load by the pushing plate. Equilibria of forces and moments are given by the following equations, respectively [16,23,24]:

$$P - P\_1 - P\_2 = 0,\tag{A12}$$

and

$$P\left(\frac{h\_{\text{coat}}}{2} + h\_{\text{sub}} - \delta\right) - M\_1 - M\_2 = 0,\tag{A13}$$

where *h*sub and *h*coat are the thicknesses of the substrate and coating, respectively, and δ is the distance from the bottom of the substrate to the neutral axis of the specimen. When *L* is small compared to *h*sub, large *P*<sup>2</sup> and *M*<sup>2</sup> can be generated, but in this case these are negligible. The specimen is a multi-layer composite of ten layers including the substrate with different Young's moduli, *Ep*, and thicknesses, *hp* (*p* = 1, 2, ··· , 10), so δ is given in a quite complicated form. Fortunately, in the loading condition for this test, the load applied to the coating is almost purely compressive parallel to the coating plane and the contribution of moment to the strain energy seems relatively small, thus the Young's modulus of the coating, *E*coat, and the substrate, *E*sub, may be simply approximated with the Voigt model in the rule of mixture, respectively,

$$E\_{\rm sub} = \frac{1}{h\_{\rm sub}} \sum E\_{\rm s} h\_{\rm s} E\_{\rm coat} = \frac{1}{h\_{\rm coat}} \sum E\_{\rm t} h\_{\rm l}, \rm s + t = 10,\tag{A14}$$

where *s*, *t*, *h*sub and *h*coat are varied according to the fractured interface. Thus δ can be obtained as the neutral axis position of a simple two-layer composite (single equivalent coating layer with *E*coat and substrate with *E*sub):

$$
\delta = \Lambda l\_{\rm coat} \, \Lambda = \frac{1 + 2\Sigma \eta + 2\Sigma \eta^2}{2\eta (1 + \Sigma \eta)},
\tag{A15}
$$

where

$$\Sigma = \frac{E\_{\text{cout}}'}{E\_{\text{sub}}'}, \eta = \frac{h\_{\text{cout}}}{h\_{\text{sub}}}, E\_{\text{cout}}' = \frac{E\_{\text{cout}}}{1 - \left(\nu\_{\text{cout}}\right)^2}, E\_{\text{sub}}' = \frac{E\_{\text{sub}}}{1 - \left(\nu\_{\text{sub}}\right)^2}. \tag{A16}$$

νcoat and νsub are Poisson's ratios of coating and substrate, respectively.

Energy release rate for the edge crack under an applied load of *P* is denoted by G and expressed as [24]

$$\mathcal{G} = \frac{1}{2w^2 E\_{\text{cout}}'} \left( \frac{P^2}{h\_{\text{cout}}} \right) + \frac{1}{2w^2 E\_{\text{sub}}'} \left\{ -\frac{\left(P\_1\right)^2}{Dh\_{\text{cout}}} - \frac{\left(M\_1\right)^2}{I \left(h\_{\text{cout}}\right)^3} + \frac{\left(P\_2\right)^2}{h\_{\text{sub}}} + 12 \frac{\left(M\_2\right)^2}{\left(h\_{\text{sub}}\right)^3} \right\}. \tag{A17}$$

Here, *D* is the dimensionless cross-sectional area and *I* is the moment of inertia of area written as

$$D = \frac{1}{\eta} + \Sigma,\\ I = \Sigma \left\{ \left(\Lambda - \frac{1}{\eta}\right)^2 - \left(\Lambda - \frac{1}{\eta}\right) + \frac{1}{3} \right\} + \frac{\Lambda}{\eta} \left(\Lambda - \frac{1}{\eta}\right) + \frac{1}{\eta^3}.\tag{A18}$$

*I* is again calculated assuming that the material system is a two-layer composite with a coating and substrate. Substituting Equations (A12), (A13) and *P* = *P*init into (A17) and ignoring *P*<sup>2</sup> and *M*2, Equation (4) is derived.

**Figure A3.** Loads and moments applied to specimen.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Coatings* Editorial Office E-mail: coatings@mdpi.com www.mdpi.com/journal/coatings

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18