Damage Assessment and Fracture Resistance of Functionally Graded Advanced Thermal Barrier Coating Systems: Experimental and Analytical Modeling Approach

Enhancement of stability, durability, and performance of thermal barrier coating (TBC) systems providing thermal insulation to aero-propulsion hot-section components is a pressing industrial need. An experimental program was undertaken with thermally cycled eight wt.% yttria stabilized zirconia (YSZ) TBC to examine the progressive and sequential physical damage and coating failure. A linear relation for parameterized thermally grown oxide (TGO) growth rate and crack length was evident when plotted against parameterized thermal cycling up to 430 cycles. An exponential function thereafter with the thermal cycling observed irrespective of coating processing. A phenomenological model for the TBC delamination is proposed based on TGO initiation, growth, and profile changes. An isostrain-based simplistic fracture mechanical model is presented and simulations carried out for functionally graded (FG) TBC systems to analyze the cracking instability and fracture resistance. A few realistic FG TBCs architectures were considered, exploiting the compositional, dimensional, and other parameters for simulations using the model. Normalized stress intensity factor, K1/K0 as an effective design parameter in evaluating the fracture resistance of the interfaces is proposed. The elastic modulus difference between adjacent FG layers showed stronger influence on K1/K0 than the layer thickness. Two advanced and promising TBC materials were also taken into consideration, namely gadolinium zirconate and lanthanum zirconate. Fracture resistance of both double layer and trilayer hybrid architectures were also simulated and analyzed.


Introduction
Thermal barrier coating (TBC) systems provide thermal insulation to gas turbine engine components exposed to high temperature and oxidizing environments for aero-propulsion. Conventionally, a TBC system is a two-layered coating consisting of eight percent YSZ (TBC), and a bond coat (BC) enriched in aluminum over a Ni base superalloy substrate. During the thermal cycling operation, a layer of oxidation product known as thermally grown oxide (TGO-α Al 2 O 3 ) forms and grows with time in between YSZ (top coat) and bond coat (BC) layer under the influence of mechanical and thermal stress cycles [1][2][3]. The predominant failure mechanisms for TBC include microcrack nucleation at the TBC/BC interface layer and coalescence and propagation of cracks to cause buckling/delamination over a large part and finally to spalling. TBC failures often occur under biaxial TBC systems continue to be of major concern. In recent times, considerable enhancement in reliability of TBC has been achieved through improvements in the coating processes, particularly in large part because of maintaining coating consistency and through understanding the failure mechanism(s) in relation to the compositional change [2,4,7]. A transition from continuous, and presumably protective, alumina scale to a chromia and spinel occurs under continued thermal exposure leading to TBC cracking as a result of aluminum depletion and large volumetric change. The delamination failure suggests that the interfacial fracture resistances of the TBC/Cr 2 O 3 and the BC/spinel interfaces are lower than that of the TBC/Al 2 O 3 interface.
In order to improve the performance and fracture resistance of the BC/TBC interface by suitable modifications, the functionally graded material (FGM) concept has been researched over the past decade [7][8][9][10][11][12][13][14][15][16][17]. In functionally gradient materials (FGMs), weak and sharp interfaces between two dissimilar materials are substituted by a number of graded layers. These layers vary in chemistry, properties, and microstructure, and progressively reduce mechanical and thermal stresses at the interface and thus enhance the performance. The graded transition across an interface of two materials can essentially reduce the thermal stresses and stress concentration at intersections with free surfaces. The stress intensity factor at the crack tip can be altered by varying the gradient properties across the interface. The metal-ceramic FGMs tend to exhibit higher fracture resistance due to crack bridging in a graded volume fraction. Varying thermal expansion in graded layers induces residual stress and affects the crack growth mode. With an increase in the number of graded layers, Sudarshan and Kokini et al. have shown that maximum surface temperature also increases at the event of cracking, while the crack length tends to decrease [8,9]. Comparing the cracking resistance between FGM TBC and conventional TBC, they confirmed that for similar thermal resistance FGM TBC offers higher resistance to interface cracking. An economical method for producing FGM TBC with three layers using the wet powder (slurry) and sintering method recently showed good quality coatings with little or no spallation [13]. The optimum time, heat flux, and applied pressure level for sintering were deduced for 40-45% ceramic, 4% binder, and 0.4% dispersant composition, and sintered for 30 min with an applied pressure of 30 MPa. Superior thermal stress relaxation in graded TBC was also demonstrated to explain better oxidation resistance and longer lifetime. High-velocity oxygen flame (HVOF) TBC with multilayers successfully withstood thermal cycling and prevented visible delamination or transverse cracks [11].
Thermo-mechanical response and fracture behavior of FG BC-TBC systems under high heat flux thermal loads carried out earlier evidenced surface cracks as well as TBC-BC interface cracks [8,9]. While time-dependent (viscoplastic) deformations led to cracking of FGMs, both crack initiation and propagation resistance were shown to increase with number of gradations. Computational fracture mechanics and thermal shock results showed that graded TBC is much more durable than the conventional TBC. Researchers at University of California, Davis, recently invented functionally graded nano-structured ceramic-metallic composite between the substrate and the top coat for enhanced corrosion, wear, and oxidation characteristics as well as extended life. The layer is formed by introducing nano-particles of aluminum oxides and nitrides into the bond coat to accelerate the formation of a uniform and protective TGO. The coating lifetime is claimed to be around double by the research group than the commonly observed coating lifetime [18].
Our current research focuses on physics-based modeling of the damage and failure process in thermally cycled TBC of gas turbine engine blades. Experimental observations of physical damages Coatings 2020, 10, 474 3 of 26 facilitate in establishing a phenomenological model. A simplistic analytical modeling based on fracture mechanics principle is proposed and used for simulations to obtain deeper insights of the fracture resistance behavior of functionally graded TBC. The work presented here attempts to address the following aspects of TBC systems as applied to gas turbines: (a) Physics-based damage analysis and modeling of failure mechanism(s) involved through experimental research. Both kinetics of oxide growth and crack formation are considered as the primary and secondary damages, respectively. (b) Fracture mechanics-based analytical approach to propose a design parameter for apply FGTBC for enhancement of fracture resistance. Development of a simplistic model to characterize fracture resistance of the interfacial region and compare various conventional and advanced TBC systems. The attempt made here should be treated as the initial step towards more refined modeling for realistic situations.

Experimental Study
The conventional TBC consisting of a bond coat (BC) and yttria stabilized zirconia (YSZ) as top coat (TC) has a sharp interface between the two layers. Prior to conducting analytical instability modeling research for TBC interface, an experimental program was undertaken to gain hands-on experience of the physical damages and failure mechanisms. The experimental study aimed to get a deeper insight on the following issues: (a) Assessment of TGO oxidation kinetics under thermal cycling effects and quantitative relations with damage induced crack length; (b) Modeling of the physical damage accumulation and cracking mechanisms in presence of changing TGO geometric profiles; (c) The following sections discuss the experimental procedures followed, physical damages assessed, results obtained, empirical relations established, and phenomenological model developed.

Experimentation
YSZ coated samples were made by air plasma spray (APS) method on the inconel superalloy flat disk substrate of size 16 mm × 8 mm. Before applying YSZ (ZrO 2 -8% Y 2 O 3 ) top coat, the substrate disk was bond coated with Co-32 Ni-21 Cr-8 Al-0.5 Y (in weight percent) using plasma spray torch system (Axial III torch, Mettech corporation, Surrey, Canada). The APS technique was used under two environmental conditions for two sets of samples, namely (a) in air (samples designated as M08) and (b) in low pressure oxygen (samples designated as M07). Low pressure oxygen environment was created in quartz tubes, evacuated to a pressure of about 0.27 Pa and back filled with argon [7,19]. The as-sprayed samples were vacuum heat treated at 1080 • C for 4 h and at a pressure less than 0.07 Pa prior to thermal cycling.
The as-sprayed and heat-treated samples were then subjected to thermal cycling in air to different numbers of cycles. Each cycle consisted of 8-12 min ramping to a temperature of 1080 • C for 45 min and then cooling slowly to ambient temperature. Four thermal exposure cycles were considered for the study in order to allow growth of TGO to different thickness and these were 10, 100, 430 h, and time until failure of TBC. Time to failure and withdrawal of samples from tests was based on appearance of significantly large spalled coating area (at least around 25% to 30% of the disk surface). The samples were sectioned, cold mounted in epoxy, and prepared metallographically and examined under a scanning electron microscope (SEM, SU7000, Hitachi, Tokyo, Japan). Representative images from different regions at various magnifications were photographed digitally and transferred to computer for further analysis. Figure 1a-d displays a few representative SEM micrographs. The SEM images of TGO regions were examined to infer both qualitative and quantitative data, namely the size/thickness, nature and geometric profile of TGO layer, crack length, and configurations.

TGO Damage Assessment
TGO thickness measurements yielding basic statistical data are important as the first step for baseline information. On the other hand, crack nucleation phenomenon is linked to local defects, microstructure and damages involving TGO size and boundary profile. To get a deeper insight of the microcrack initiation mechanism and probabilistic failure time prediction, a probabilistic analysis using statistical distribution is also required. Two measurement methods using the SEM micrographs provide both the baseline information and probabilistic distribution of the physical damages. The quantitative damage assessment was carried out for TGO as primary source and induced crack growth as the secondary damages. The thickness of TGO oxide layers as well as the crack lengths associated in different samples were measured using software named Image Tool (version 3.5) [19]. Large numbers of SEM images were used for damage assessments for each sample and a large number of measurements were made from each of the images, as reported in Table 1.

TGO Damage Assessment
TGO thickness measurements yielding basic statistical data are important as the first step for baseline information. On the other hand, crack nucleation phenomenon is linked to local defects, microstructure and damages involving TGO size and boundary profile. To get a deeper insight of the microcrack initiation mechanism and probabilistic failure time prediction, a probabilistic analysis using statistical distribution is also required. Two measurement methods using the SEM micrographs provide both the baseline information and probabilistic distribution of the physical damages. The quantitative damage assessment was carried out for TGO as primary source and induced crack growth as the secondary damages. The thickness of TGO oxide layers as well as the crack lengths associated in Coatings 2020, 10, 474 5 of 26 different samples were measured using software named Image Tool (version 3.5) [19]. Large numbers of SEM images were used for damage assessments for each sample and a large number of measurements were made from each of the images, as reported in Table 1. The first method of TGO oxide layer thickness estimation is an indirect approach and is based on the entire TGO region, where the area under high magnifications was estimated using the Image Tool software in pixel units. Area-based mean TGO oxide layer thickness (T A ) is obtained as T A = A/L, where A is the overall cross-sectional area measured from an image and L is the shortest distance to cover the area along the length of the area. The area-based method is simple and fast to provide baseline information (mean) without accounting for thickness variation (i.e., range, peaks, etc.). However, at least three times for each image the TGO area was scanned and average thickness was estimated.

Point to Point Measurement
In the second method, point to point measurement technique (PPT) was followed. Image Tool allows measurement of distances between two points in pixel units, as illustrated in Figure 1. Both magnification and spatial measurements were considered for determining calibration factors. The distance between the two points for TGO thickness lie on the two boundaries, one on BC-TGO and the other on TGO-TBC layer. However, because of irregular shape of boundaries, in particular for TGO-TBC, choosing the line often becomes an arguable issue and a source of error. To maintain consistency, criterion of shortest distance was followed as far as practicable. Around 80 to 100 measurements were made from each image at intervals of 1-2 mm, and at least five images were used for a sample. Thus, enough data is generated to provide meaningful statistical analysis.

Crack Length
Estimation of crack lengths from images was done only by the direct method, i.e., by PPT technique. Cracks are often observed to be nonlinear as these form mostly at the interfacial curvilinear regions between the layers or other places. Depending on the size and curvature, the crack length may appear to be different when measured end to end (assuming linear shape) and measured along the curved profile. This is illustrated for a crack in TBC in Figure 1d. When the difference in length appeared to be unacceptably different, the measurement of crack length was made following the shape of the curve. For all small sized cracks less than 5 µm or so, measurements were made assuming crack linearity [20].

Microstructural Analysis
Sample micrographs as displayed in Figure 1 illustrate the layer arrangements, interface profile, and growth of TGO at the BC-YSZ interface after thermal cycling. All SEM micrographs displayed here are at 300× magnifications for visual comparisons of TGO size and profile changes. Clearly, with more thermal cycling the TGO thickens and the geometry changes to more complex shapes and curvatures. In other words, highly heterogeneous oxide growth is evident irrespective of the time of exposure in all samples exposed to air and limited oxygen environments. The BC-TBC system is designated as zero (0) layer system with thickness of two layers as 150 and 360 µm, respectively for BC and YSZ TBC. The conventional zero layer system has one interface with a high degree of mismatch between the two adjacent layers. High residual stress develops at the interface to facilitate cracking even with very small preexisting crack like defects. The zigzag profile of the BC-YSZ interface is held responsible for stress concentration to cause cracking at the interface and in YSZ. Complex TGO profile, extensive damage, and cracking at failure are evident from Figure 1d as compared to other thermally exposed conditions.

TGO Thickness and Crack Length
The physical damage measurement results are summarized in Table 1, which compiles mean of TGO thickness and crack length, uncertainties, and number of measurement points. The TGO mean data as obtained by both measurement methods are shown in Figure 2. The data points fall closely to an ideal line that signifies one-to-one correlation between the two sets. An excellent matching of data points with the line is an indication that both the methods are effective to characterize the TGO growth. Besides the individual sample mean data, the population mean is also marked by x in Figure 2 for the sample set. The population mean is also in excellent agreement with the idealistic line except for cycles at failures (1500 cycles), as may be seen in Figure 2. This is due to large scale cracking of TGO and TBC layers. The widespread transverse cracking mode of TGOs ( Figure 1d) seems to have introduced appreciable errors in the measurements. The points to points measurement method is relatively less influenced by these transverse cracks as compared to direct method. Figure 2 also gives an estimate of TGO growth as a function of thermal cycles and clearly signifies the linear relation between the two. An eight-fold increase in the size of TGO is evident in Figure 2 and Table 1. exposure in all samples exposed to air and limited oxygen environments. The BC-TBC system is designated as zero (0) layer system with thickness of two layers as 150 and 360 μm, respectively for BC and YSZ TBC. The conventional zero layer system has one interface with a high degree of mismatch between the two adjacent layers. High residual stress develops at the interface to facilitate cracking even with very small preexisting crack like defects. The zigzag profile of the BC-YSZ interface is held responsible for stress concentration to cause cracking at the interface and in YSZ. Complex TGO profile, extensive damage, and cracking at failure are evident from Figure 1d as compared to other thermally exposed conditions.

TGO Thickness and Crack Length
The physical damage measurement results are summarized in Table 1, which compiles mean of TGO thickness and crack length, uncertainties, and number of measurement points. The TGO mean data as obtained by both measurement methods are shown in Figure 2. The data points fall closely to an ideal line that signifies one-to-one correlation between the two sets. An excellent matching of data points with the line is an indication that both the methods are effective to characterize the TGO growth. Besides the individual sample mean data, the population mean is also marked by x in Figure  2 for the sample set. The population mean is also in excellent agreement with the idealistic line except for cycles at failures (1500 cycles), as may be seen in Figure 2. This is due to large scale cracking of TGO and TBC layers. The widespread transverse cracking mode of TGOs ( Figure 1d) seems to have introduced appreciable errors in the measurements. The points to points measurement method is relatively less influenced by these transverse cracks as compared to direct method. Figure 2 also gives an estimate of TGO growth as a function of thermal cycles and clearly signifies the linear relation between the two. An eight-fold increase in the size of TGO is evident in Figure 2 and Table 1. The scatter and uncertainty are determined by calculating standard deviation (SD) and both population mean and standard deviation for TGO thickness. These two statistical quantities are graphically displayed in Figure 3 for all samples and for both the measured parameters. Interestingly, all the measured data points fall closely around a theoretical exponential function. However, uncertainty and scatter tend to increase exponentially as the measured parameters continue to increase. Though a matter of great concern, this implies longer the duration of thermal exposure, of  The scatter and uncertainty are determined by calculating standard deviation (SD) and both population mean and standard deviation for TGO thickness. These two statistical quantities are Coatings 2020, 10, 474 7 of 26 graphically displayed in Figure 3 for all samples and for both the measured parameters. Interestingly, all the measured data points fall closely around a theoretical exponential function. However, uncertainty and scatter tend to increase exponentially as the measured parameters continue to increase. Though a matter of great concern, this implies longer the duration of thermal exposure, of more concern will be the measured damage data and higher uncertainty in prediction. In other words, probabilistic consideration becomes more meaningful than a deterministic approach for prediction of life. The best fit function that computes the uncertainty in both the parameters is expressed by where x represents the TGO thickness or crack length in µm.
Coatings 2020, 10, x FOR PEER REVIEW 7 of 25 SD = 0.70 × exp(x/5) where x represents the TGO thickness or crack length in μm.
It is worthwhile to mention here that the oxygen availability and pressure level during APS processing do not seem to have any noticeable influence on the oxidation characteristics and damage behavior (Table 1). Percentage changes in mean thickness between two classes indicate small and inconsistent variation for the same hours of thermal exposures. The differences in thickness between the two classes turn out to be around ±7% for 430 h and up to coating failure time. However, for all further analysis we will treat the two sample groups separately, even though the data may tend to overlap. Presence of mixed oxides consisting of chromia and spinels in TGO was reported earlier when thermally exposed in air. In absence of mixed oxides in low oxygen pressure, TBC is reported to have a longer lifetime [21][22][23]. The mixed oxides facilitate easy crack nucleation due to volumetric expansion of TGO. In the present analysis, the enhanced coating lifetime in air (sample M08) is found to be around 12.8% as compared to life in a low oxygen environment. The result is found to be in contrast to earlier observations [22]. Symbol ∆ represents crack length-SD data for M08 samples; symbol • represents TGO-SD data for M07 samples; and symbol ◊ represents TGO-SD data for M08 samples.

Thermal Cycling vs. TGO Growth
The TGO growth is found to have a linear relation with thermal cycling when both the parameters are rationalized. Thermal exposure is rationalized with respect to time to failure (t/tf), while TGO growth (Ttgo) is made dimensionless by dividing the TBC layer thickness (Ttgo/Tlayer). All data points for two classes of APS treatments fall close to each other to be represented by the empirical relation as: Ttgo/Tlayer = m(t/tf) + c. However, a nonlinear TGO growth is evident for both the coating processes used in the experiments.
The present work contradicts the findings of others who have reported a three-stage TGO growth and parabolic function between the minimum TGO thickness and thermal cycling [7,14,15]. The present work accounts for oxidation growth leading to the failure of coating unlike the early stage of oxidation as considered in other works. Parabolic TGO thickening kinetic was also reported where thickness per cycle (Δh) is given by (h -h0) × Δt/2t, where h0 is the initial TGO thickness, h is the TGO thickness at time t. An increment in time (Δt) is rationalized by cumulative hot time, t [16,22]. Earlier work reported enhanced tendency for heterogeneous oxide growth with increase in exposure time [7]. Dependency of TBC lifetimes on oxidation temperature and TGO growth is discussed elsewhere considering compiled literature data [1]. These data points are taken for both MCrAlY BC with APS processing and Pt-aluminide BC with electron beam physical vapor deposition (EBPVD) processing subjected to both isothermal and thermal cycling effects. However, the data points exhibit Symbol ∆ represents crack length-SD data for M08 samples; symbol • represents TGO-SD data for M07 samples; and symbol ♦ represents TGO-SD data for M08 samples.
It is worthwhile to mention here that the oxygen availability and pressure level during APS processing do not seem to have any noticeable influence on the oxidation characteristics and damage behavior (Table 1). Percentage changes in mean thickness between two classes indicate small and inconsistent variation for the same hours of thermal exposures. The differences in thickness between the two classes turn out to be around ±7% for 430 h and up to coating failure time. However, for all further analysis we will treat the two sample groups separately, even though the data may tend to overlap. Presence of mixed oxides consisting of chromia and spinels in TGO was reported earlier when thermally exposed in air. In absence of mixed oxides in low oxygen pressure, TBC is reported to have a longer lifetime [21][22][23]. The mixed oxides facilitate easy crack nucleation due to volumetric expansion of TGO. In the present analysis, the enhanced coating lifetime in air (sample M08) is found to be around 12.8% as compared to life in a low oxygen environment. The result is found to be in contrast to earlier observations [22].

Thermal Cycling vs. TGO Growth
The TGO growth is found to have a linear relation with thermal cycling when both the parameters are rationalized. Thermal exposure is rationalized with respect to time to failure (t/t f ), while TGO growth (T tgo ) is made dimensionless by dividing the TBC layer thickness (T tgo /T layer ). All data points for two classes of APS treatments fall close to each other to be represented by the empirical relation as: T tgo /T layer = m(t/t f ) + c. However, a nonlinear TGO growth is evident for both the coating processes used in the experiments. The present work contradicts the findings of others who have reported a three-stage TGO growth and parabolic function between the minimum TGO thickness and thermal cycling [7,14,15]. The present work accounts for oxidation growth leading to the failure of coating unlike the early stage of oxidation as considered in other works. Parabolic TGO thickening kinetic was also reported where thickness per cycle (∆h) is given by (h − h 0 ) × ∆t/2t, where h 0 is the initial TGO thickness, h is the TGO thickness at time t. An increment in time (∆t) is rationalized by cumulative hot time, t [16,22]. Earlier work reported enhanced tendency for heterogeneous oxide growth with increase in exposure time [7]. Dependency of TBC lifetimes on oxidation temperature and TGO growth is discussed elsewhere considering compiled literature data [1]. These data points are taken for both MCrAlY BC with APS processing and Pt-aluminide BC with electron beam physical vapor deposition (EBPVD) processing subjected to both isothermal and thermal cycling effects. However, the data points exhibit considerable scattering without any consistent effects of BC type and thermal exposure. For the temperature of 1080 • C used in the present work, the time to reach comparable TGO thickness as found in this work is around 550-600, 160, and 30 h [7]. These data are quite low as compared to corresponding values found in our work, i.e., 100, 430, and 1428-1612 h respectively for thicknesses around 2.8, 5, and 11.5 µm. The present study also reveals a wide distribution in TGO thickness data with exponential increase in standard deviation. Scatter of TBC life time and TGO growth by a factor of 10 is attributed to BC surface roughness and possible mechanisms leading to nucleation and growth of delamination cracks [1,18]. The result demonstrates the inadequacy of lifetime analysis based on a deterministic approach, i.e., either mean or minimum thickness TGO data. This constitutes the basis for further analysis of TGO growth kinetics considering statistical distribution of TGO and lifetime as discussed elsewhere [20].

Linearity vs. Nonlinearity Trend
The TGO growth trend with thermal exposure is the most important parameter in understanding and modeling of the TBC delamination process. A simplistic representation of the TGO growth is an unambiguous input for any quantitative modeling and lifetime prediction. With the limited data generated from the experimental investigation in this work, TGO thickness data are plotted in two ways as displayed in Figure 4a,b. Untreated TGO thickness tend to follow nonlinearity as a function of time of thermal exposure, as shown in Figure 4a, for both the coatings considered in the work (M07 and M08). The data points fit best with the 2-degree polynomials for both coatings with R-square values close to one. It may be noted here the coefficients of the square terms are very small, as may be seen in the fitted equations shown in Figure 4a. A different explanation of the trend may also emerge with a closer look of the how the trend is made. It reveals that the first three points (up to 430 cycles) fits linearly while the last part, i.e., 430 to failure point of the plot follows more nonlinear trend.
Coatings 2020, 10, x FOR PEER REVIEW 8 of 25 considerable scattering without any consistent effects of BC type and thermal exposure. For the temperature of 1080 °C used in the present work, the time to reach comparable TGO thickness as found in this work is around 550-600, 160, and 30 h [7]. These data are quite low as compared to corresponding values found in our work, i.e., 100, 430, and 1428-1612 h respectively for thicknesses around 2.8, 5, and 11.5 μm. The present study also reveals a wide distribution in TGO thickness data with exponential increase in standard deviation. Scatter of TBC life time and TGO growth by a factor of 10 is attributed to BC surface roughness and possible mechanisms leading to nucleation and growth of delamination cracks [1,18]. The result demonstrates the inadequacy of lifetime analysis based on a deterministic approach, i.e., either mean or minimum thickness TGO data. This constitutes the basis for further analysis of TGO growth kinetics considering statistical distribution of TGO and lifetime as discussed elsewhere [20].

Linearity vs. Nonlinearity Trend
The TGO growth trend with thermal exposure is the most important parameter in understanding and modeling of the TBC delamination process. A simplistic representation of the TGO growth is an unambiguous input for any quantitative modeling and lifetime prediction. With the limited data generated from the experimental investigation in this work, TGO thickness data are plotted in two ways as displayed in Figure 4a,b. Untreated TGO thickness tend to follow nonlinearity as a function of time of thermal exposure, as shown in Figure 4a, for both the coatings considered in the work (M07 and M08). The data points fit best with the 2-degree polynomials for both coatings with R-square values close to one. It may be noted here the coefficients of the square terms are very small, as may be seen in the fitted equations shown in Figure 4a. A different explanation of the trend may also emerge with a closer look of the how the trend is made. It reveals that the first three points (up to 430 cycles) fits linearly while the last part, i.e., 430 to failure point of the plot follows more nonlinear trend.
However, the above discussion indicates the chosen exposure time in the thermal cycling experiment in the present work was not judiciously chosen. The time of exposure should have been chosen more uniformly over the TBC failure time. In different words, two-way trends also suggest the possibilities of two-stage delamination mechanisms involved and should be a subject of future research.
This analysis is further extended when the parameterized TGO thickness and time of exposure is plotted, as displayed in Figure 4b. This confirms a linear trend for both samples with R-square values almost close to unity. The general form of the trend can be expressed as Ttgo = m × Tlayer × (t/tf) + c, where Ttgo is the TGO thickness, Tlayer is the top-coat thickness, t is the time of exposures, and tf is the time to failure.  However, the above discussion indicates the chosen exposure time in the thermal cycling experiment in the present work was not judiciously chosen. The time of exposure should have been chosen more uniformly over the TBC failure time. In different words, two-way trends also suggest the possibilities of two-stage delamination mechanisms involved and should be a subject of future research. This analysis is further extended when the parameterized TGO thickness and time of exposure is plotted, as displayed in Figure 4b. This confirms a linear trend for both samples with R-square values almost close to unity. The general form of the trend can be expressed as where T tgo is the TGO thickness, T layer is the top-coat thickness, t is the time of exposures, and t f is the time to failure.

Crack Length vs. TGO Thickness
The mean crack length appears to follow linearity initially with thermal cycling up to 430 cycles and thereafter has an exponential function for higher number of thermal cycling, irrespective of coating processing. The maximum crack length was reported earlier to lie around 50 to 300 µm for the similar TBC system [7] up to 600 h of exposure. In contrast, the present work finds the mean crack length to be around 5 to 25 µm. The crack growth behavior up to 430 cycles is of linear nature, while the dependency changes to exponential form for further exposure up to lifetime. A qualitative explanation of the exponential function can be given by considering a number of factors closely associated with secondary damaging effects. These are primarily presence of fabrication induced preexisting flaws and discontinuities at BC/TBC interface, enlargement as well as formation of new cracks during thermal exposure, growth (sidewise and axial) and coalescence of cracks, stress build-up, strain mismatch effects with TGO formation, etc. [7]. Physical changes in mechanisms of cracking are also observed with the continuation in thermal cycling beyond 430 cycles. The population of cracks is seen to be lying at the interfacial regions between TBC and TGO during early stages (up to 430 cycles). However, cracking away from the interfaces as well as the transverse cracking (cracking through the TGOs) become dominant as TGO grows with time, as displayed in Figure 5. Mechanisms of cracking in TGO at different thermal cycles are discussed in a later section.

Crack Length vs. TGO Thickness
The mean crack length appears to follow linearity initially with thermal cycling up to 430 cycles and thereafter has an exponential function for higher number of thermal cycling, irrespective of coating processing. The maximum crack length was reported earlier to lie around 50 to 300 μm for the similar TBC system [7] up to 600 h of exposure. In contrast, the present work finds the mean crack length to be around 5 to 25 μm. The crack growth behavior up to 430 cycles is of linear nature, while the dependency changes to exponential form for further exposure up to lifetime. A qualitative explanation of the exponential function can be given by considering a number of factors closely associated with secondary damaging effects. These are primarily presence of fabrication induced preexisting flaws and discontinuities at BC/TBC interface, enlargement as well as formation of new cracks during thermal exposure, growth (sidewise and axial) and coalescence of cracks, stress buildup, strain mismatch effects with TGO formation, etc. [7]. Physical changes in mechanisms of cracking are also observed with the continuation in thermal cycling beyond 430 cycles. The population of cracks is seen to be lying at the interfacial regions between TBC and TGO during early stages (up to 430 cycles). However, cracking away from the interfaces as well as the transverse cracking (cracking through the TGOs) become dominant as TGO grows with time, as displayed in Figure 5. Mechanisms of cracking in TGO at different thermal cycles are discussed in a later section. A functional relationship between the two forms of physical damages, namely crack damage and TGO oxide thickness, is shown in Figure 5. Clearly and consistently, an exponential dependency between statistical mean of crack length in different samples and the mean TGO thickness data is evident beyond 430 cycles. Up to 430 cycles, a linear dependency between the two is evident from Figure 5. However, more experimental data is needed for confirmation of this trend as only one data point is available in this work for 1612 cycle. The two plots for samples M07 and M08 lie fairly close to each other demonstrating no significant effects of oxygen pressure on the physical damage size in the TBC system as mentioned in an earlier section. A linear relation between equivalent TGO thickness and maximum crack length was reported earlier for the same TBC system [7]. This trend suggests two-stage TGO growths kinetic from early oxidation until the failure time. During early stage of TBC life (up to 430 cycles), the cracking mechanisms in TBC is predominantly crack  A functional relationship between the two forms of physical damages, namely crack damage and TGO oxide thickness, is shown in Figure 5. Clearly and consistently, an exponential dependency between statistical mean of crack length in different samples and the mean TGO thickness data is evident beyond 430 cycles. Up to 430 cycles, a linear dependency between the two is evident from Figure 5. However, more experimental data is needed for confirmation of this trend as only one data point is available in this work for 1612 cycle. The two plots for samples M07 and M08 lie fairly close to each other demonstrating no significant effects of oxygen pressure on the physical damage size in the TBC system as mentioned in an earlier section. A linear relation between equivalent TGO thickness and maximum crack length was reported earlier for the same TBC system [7]. This trend suggests two-stage TGO growths kinetic from early oxidation until the failure time. During early stage of TBC life (up to 430 cycles), the cracking mechanisms in TBC is predominantly crack nucleation and opening controlled, while the mechanisms changes to predominant propagation-controlled mode towards the later stage (exceeding 430 cycles) and until TBC failure. Two stage kinetics of oxide growth-cracking relation can be related to the changes in oxide phase changes. Alumina formation changes to mixed oxide formation consisting of Cr, Ni, Co, Al) due to depletion in Al. Growth rate for mixed oxide is higher than alumina [21]. It is confirmed that thermal cycling up to 500 cycles TGO is made of pure alumina, while the internal oxide in BC consists of both alumina and spinels. However, some uncertainties are reported to be always associated with oxide scale measurements because of higher instability of spinels as compared to alumina.

Phenomenological Modeling
The damage and cracking mechanism(s) in TBC as observed in samples subjected to different thermal cycling during experimentation are discussed with an objective to establish a phenomenological modeling for the process of degradation. General pattern of TGO growth and specific geometrical changes with thermal cycling are in particular emphasized. Figure 6a-e represents schematically the physical damaging process as well as the changing geometry of the TGO phase in different thermally exposed conditions of TBC. Oxide layer is observed to occur by two mechanisms, namely nucleation and growth. Nucleation of oxide discretely begins at the BC-TBC interface and tends to grow along both directions, i.e., in a direction transverse to the interface and along the interface following the curvature. However, both nucleation and growth rate are highly heterogeneous in nature resulting in very irregular and uneven TGO-TBC boundaries as reflected in different SEM images shown in Figure 1. The lower boundaries of TGO (i.e., TGO-BC interface) more or less maintain similar profile as that of as-sprayed BC-TBC interface (Figure 6a). In other words, as compared to the upper boundary of TGO layer, the lower edge is more even and has regular shape. This demonstrates that due to heterogeneous TGO growth, the upper edge of TGO tends to become rough and irregular with sharp profiles.
Due to non-uniform oxidation rate at different sites, often the islands of TGOs remain connected with a thin region as illustrated schematically in Figure 6b. During the early stage of oxidation, even the TGO islands remain isolated. Alternately, islands of ceramic coatings are also observed where the oxidation occurred all around. Such discrete TGO growth is only conducive because of local supporting microstructures. Diffusion of aluminum atoms from the BC and oxygen from atmosphere through TBC occurs through various microstructural paths and sites like grain boundary, dislocation network, vacancies, pores, etc. Around 10% to 15% porosities are incorporated deliberately between splats of successively deposited TBC materials by APS technique in order to enhance strain compliance and reduce thermal conductivity [10,18]. In EB-PVD coating, columnar structures also contain microscopic porosity that reduces conductivity. These defects have the segregation tendencies in the microstructure influencing the local diffusion characteristics and oxidation kinetics.
The cracking and interface separation mechanisms are influenced by shape and curvatures at the TGO boundaries. Cracks are mostly observed at, and around, sharp locations and ridges as shown in Figure 6a-e. At these locations, stresses must be higher and deformability may be less compatible between the adjacent phases with growing TGO. In APS technique, rough bond coat with roughness around 5-10 µm are employed to ensure effective mechanical locking. Tensile stress generates above peaks while compressive stress components in the valley [22][23][24][25][26][27]. BC contracts more than topcoat during the cooling segments of thermal cycles. Hence tensile stress will generate around peaks. In general, stress concentration becomes higher as the raddi of curvature reduces. Around the peaks and ridges, the raddi of curvature is less as compared to other regions. These regions are found to be the weakest among other regions. For thermal cycling up to 100, the interface separations are observed only at TGO-TBC boundary. However, longer cycling time (430 and more) induces separation and cracking also at TGO-BC interface and also around less vulnerable locations along the upper interface.
in Figure 6a-e. At these locations, stresses must be higher and deformability may be less compatible between the adjacent phases with growing TGO. In APS technique, rough bond coat with roughness around 5-10 μm are employed to ensure effective mechanical locking. Tensile stress generates above peaks while compressive stress components in the valley [22][23][24][25][26][27]. BC contracts more than topcoat during the cooling segments of thermal cycles. Hence tensile stress will generate around peaks. In general, stress concentration becomes higher as the raddi of curvature reduces. Around the peaks and ridges, the raddi of curvature is less as compared to other regions. These regions are found to be the weakest among other regions. For thermal cycling up to 100, the interface separations are observed only at TGO-TBC boundary. However, longer cycling time (430 and more) induces separation and cracking also at TGO-BC interface and also around less vulnerable locations along the upper interface.  The higher cracking and separation tendencies at the TGO-TBC interface are primarily because of lower fracture resistance as compared to the BC-TGO interface. The fracture energy of the TBC-TGO interface is around only 2-4 Mpa·m 1/2 after 100 h of oxidation. The fracture energy for BC-TGO layers has not been reported so far, but it is likely to be much higher because BC layer toughness is more than 60 Mpa·m 1/2 [12]. Furthermore, as the TGO thickness further increases with thermal cycling, the transverse cracking across the layer is also observed (Figure 6e). During the early stage of growth, the separation and cracking mostly remain at, and around, the TGO-TBC interfaces. TBC instability Figure 6. Phenomenological modeling of major damaging mechanisms in presence of changing TGO geometry and size under the influence of increasing numbers of thermal cycling (room temperature to 1080 • C). (a) As-sprayed BC-TC interface profile with roughness around 5-10 µm; (b) Thermally exposed samples for 10 cycles. 1. Oxidation around ceramic island, 2. Scattered nucleation and limited growth, 3. Small cracks around peaks, 4. Adjacent nucleation and growth connected by thin ligament; (c) Thermally exposed samples for 100 cycles, 1. Oxide islands linked by thin ligament, 2. Uneven oxidation, 3. Discontinuous oxide growth, 4. Short cracks around peaks, 5. Long cracks along interface; (d) Thermally exposed samples for 430 cycles; 1. Divorced islands of oxide, 2. Thin TGO links, 3. Massive oxide growth with fine cracks, 4. TGO-BC interface separation, 5. Series/discrete cracking at TGO-TC interface, 6. Long TC cracking, 7. Cracks around TGO peaks; (e) Thermally exposed samples at delamination, 1. TGO growth around ceramic island, 2. Complex profile of TGO at both edges, 3. Separated TGO islands with cracks at interface, 4. Series of transverse parallel cracks, 5. Long cracks in TC parallel to TGO, 6. Extensive cracking network in TC.
The higher cracking and separation tendencies at the TGO-TBC interface are primarily because of lower fracture resistance as compared to the BC-TGO interface. The fracture energy of the TBC-TGO interface is around only 2-4 Mpa·m 1/2 after 100 h of oxidation. The fracture energy for BC-TGO layers has not been reported so far, but it is likely to be much higher because BC layer toughness is more than 60 Mpa·m 1/2 [12]. Furthermore, as the TGO thickness further increases with thermal cycling, the transverse cracking across the layer is also observed (Figure 6e). During the early stage of growth, the separation and cracking mostly remain at, and around, the TGO-TBC interfaces. TBC instability analysis using fracture mechanical models earlier was based on isostrain condition during early life of TGO growth and on TGO buckling condition for crack propagation [28]. Under isostrain condition, an applied stress of magnitude around 100 MPa results in K 1 values of 2 to 4 Mpa·m 1/2 , which is just enough for initiating a crack from a defect size of 2 µm [12]. The stress level required for crack driving force (K 1 /G 1 ) to be around and larger than the fracture resistance of the TBC component materials (TBC, TGO) are found to be in the range of 50 to 500 MPa for crack length exceeding the critical length (>1 µm) in TGO. The approximate stress levels as required for both crack initiation and propagation stages and computed from two models are consistent. The present experimental work confirms the presence of cracks in the TBC system in the size range of 2-3 µm even at the onset of thermal cycling (around 10 cycles) as shown in Figure 6b and Table 1. Thus, the results obtained earlier from fracture mechanical model analysis and the experimental findings are in excellent agreement to support the TBC instability conditions. An increasing tendency for stress intensity factor from 2 to 3.2 Mpa·m 1/2 was reported earlier with TGO growth from 4 to 8 µm [29]. The elastic energy release rate, G in J/m 2 was estimated to rise from 28 to 73 for the identical TGO growth. FE model also predicts that cracking initiates at specific geometric features of the rough interface of a PS coated system, which was confirmed by metallographic examination. Phenomenological model is a description of the damaging process observed in TBC delamination under thermal cycling. This presents in qualitative terms how TGO and crack length are intimately associated with the delamination and thermal cycling stages. This physical model could be the very first step before further research towards the development of microstructure-based deterministic modeling in future research. At present this modeling approach is complex and full of uncertainties, and is insufficiently understood for a quantitative microstructure-based modeling.

Model of Layered TBC
In functionally graded TBC, the conventional YSZ top layer is substituted with a number of graded layers with varying contents in YSZ. A model for multilayer functionally graded coating system for TBC is shown in Figure 6 with bond coat and YSZ as the bottom and top layers, respectively. The instability model was developed by the authors earlier for both crack initiation and crack propagation modes considering two-layer TBC systems [28]. In this section, the essential part of the model is discussed here to help understand the necessary modifications made for applying to a multilayered system. For fracture mechanics analysis, a three-dimensional elemental block inside the FGM beam is considered as shown in Figure 7. The block contains three adjacent layers and two interfaces. A vertical plane section the block and Figure 8 shows the 2D view of the plane after sectioning. The middle layer (layer 2) is assumed to contain a preexisting through thickness crack that is subjected to uniaxial transverse stress (mode 1 loading). The crack in layer 2 tends to initiate when subjected to applied stress. The model presented here assumes isostrain condition prevails in the composite TBC beam. The advantage is its simplicity, as the stiffness of a composite layer tends to vary linearly with the volume fraction of the constituent phases of the composites as opposed to isostress condition. In order to maintain integrity in the TBC systems, strain in all layers need to be constant while stress may develop differently. However, isostrain-based modeling requires a number of assumptions to be made that are outlined in the following sections.

Assumptions
A few basic assumptions are made to simplify the formulation of the problem and analysis thereafter considering the fracture mechanical approach [28].

Assumptions
A few basic assumptions are made to simplify the formulation of the problem and analysis thereafter considering the fracture mechanical approach [28].

Assumptions
A few basic assumptions are made to simplify the formulation of the problem and analysis thereafter considering the fracture mechanical approach [28].

•
Firstly, the different layers in TBC FGM are subjected to isostrain condition under external applied force prior to the onset of crack initiation. • Secondly, the energy balance criterion is maintained during the stage of loading and straining. The work done by the external force is used up in the process of crack initiation and this may be estimated from the local stress distribution in the vicinity of crack tip.

•
Thirdly, linear elastic fracture mechanic (LEFM) situation is maintained during mode I loading around the crack tip so that stress intensity factor criterion is the appropriate solution to deal with cracking of FGM.
• Furthermore, only externally applied mechanical forces are present and no stress effect due to temperature fluctuations is considered. TBC FGM is considered as an asymmetrical composite system where only layers of identical dimensions but different characteristics are considered for fracture mechanical analysis.

•
In all the analysis presented here, the preexisting crack length as assumed in the model is considered to be a deep crack with crack length to width ratio as 0.5.

•
No thermal and TGO stress effect is considered in the analysis, and the multilayer TBC system is assumed to operate at room temperature only. Room temperature material properties are considered in the simulations. Among other material properties affecting the thermal protection performance and lifetime, thermal conductivity (TC) play a significant role. This is truer for graded TBC systems with different compositions and materials. However, for a simple approach considered in the modeling, TC is not considered, as thermal effects due to fluctuations is beyond the scope of the present research. Also, the influence of TC in altering the RT modulus is not taken into consideration, as sufficient data is not available.

Fracture Mechanics Model
The physical model of TBC FGM as described above can now be used for developing a fracture mechanics (FM) model for further analysis. Considering Westergaard's stress function and applied stress, the generic mode I (opening mode) stress intensity factor (SIF) expression is given as [ where all terms have the usual meanings. σ app stands for the applied stress, Y is the Young's moduls, and a is the crack length.
The remote applied stress, σ will influence the local stress distribution and it is necessary to estimate the local stress for energy balance approach as stated earlier in the assumptions. The local stress distribution in layer 1 will be affected by two factors, namely the presence of the preexisting crack of assumed length of 2a ( Figure 8) in layer 2 and the presence of adjacent layer and interface (layer 3). Under isostrain condition, however, stresses are unequal in the layers, i.e., σ 1 σ 2 .
The generic form of Westergaard's stress functions can now be considered to represent the local stress distribution. The local stress, σ yy as a function of distance, x is given as [30] The local stress distribution as discussed is considered for force equilibrium for the determination of stress intensity factor through a relation between the local stress and applied stress. The force estimation from the local stress distribution may be obtained from the area under the curve. A small element of width 'dx' is assumed over the stressed region under the curves and by carrying out integration within the limits, the force may be obtained. The total force, thus, from local stress over the length (a to W 1 + W 2 as shown in is given by [28]), W stands for layer thickness, z is a quantity used to replace [1 − (a/x) 2 ] −1/2 , so the dimensionless normalized SIF (K 1 /K 0 ) can be expressed as Now the terms in the integral can be rearranged for a standard integration, The integral is now in a standard form and its integration is given as (x 2 − a 2 ) 1/2 (Equation (6)) can be rewritten using the limits of integrals as Substituting m for W 1 /W and n for a/W, (Equation (7)) is simplified to where m is the crack depth to width ratio and n is the ratio of two adjacent layer widths as shown in Figure 8. For the interface between layer 1 and layer 2, n should be considered as the ratio of W 1 to W and for the other interface, i.e., between layer 2 and 3, n should be ratio between W 2 to W. E 1 , and E 2 are the elastic moduli. Equation (8) has been used for simulations in Matlab to compute the normalized SIF for various situations that can be exploited for redesigning the BC-TBC interface by optimization [31]. K 1 /K 0 provides the necessary driving force for the preexisting crack to initiate cracking. The simulation results are discussed in the following section.

Functionally Graded TBC
A Functionally Graded (FG) TBC has both the composition and the microstructure gradually change over the TC, resulting in corresponding changes in the properties of the material. This is illustrated in Figure 9b Figure 9 shows four possible architectures for TBC systems as considered for modeling, simulations, and analysis in the present work. The conventional single layer system with one interface between BC and TC may be described as the base case as displayed in Figure 9a and is subjected to various simulations. Bilayer and trilayer systems (Figure 9c,d) are discussed and analyzed in later sections. The conventional TC is redesigned following the concept of functional gradient (FG) as discussed earlier with respect to Figures 7 and 9b. The interface between BC and single-layer TC is weak for a number of reasons and is susceptible to early cracking and failure limiting the durability. FG layers can substitute the single layer TC when the layer composition can be varied in uniform and random manners from 100% TC at the top and 0% (i.e., BC) at the bottom ( Table 2 and Figure 9b). The prime design motive is to protect the substrate with maximum efficiency, high durability, and to be cost-effective. In the present analysis, the BC layer thickness is assumed to be constant as 150 μm and the topcoat thickness is taken as constant at 360 μ, as shown in Figure 9ad. These values are chosen to be consistent with our experimental data as presented and discussed in an earlier section.

Interfaces
In FG systems with 2 and more layers, it is necessary to distinguish the strength and weakness of individual interfaces. Figure 10 displays how different interfaces are designated with numbers (1,2,3…) for identification in only four FG TBC systems. K1/K0 values are computed for all the interfaces in different FG TBC systems for the comparative evaluation of interfaces. In all systems, the first interface represents the interface between the TC and BC.  Figure 9 shows four possible architectures for TBC systems as considered for modeling, simulations, and analysis in the present work. The conventional single layer system with one interface between BC and TC may be described as the base case as displayed in Figure 9a and is subjected to various simulations. Bilayer and trilayer systems (Figure 9c,d) are discussed and analyzed in later sections. The conventional TC is redesigned following the concept of functional gradient (FG) as discussed earlier with respect to Figures 7 and 9b. The interface between BC and single-layer TC is weak for a number of reasons and is susceptible to early cracking and failure limiting the durability. FG layers can substitute the single layer TC when the layer composition can be varied in uniform and random manners from 100% TC at the top and 0% (i.e., BC) at the bottom (Table 2 and Figure 9b). The prime design motive is to protect the substrate with maximum efficiency, high durability, and to be cost-effective. In the present analysis, the BC layer thickness is assumed to be constant as 150 µm and the topcoat thickness is taken as constant at 360 µ, as shown in Figure 9a-d. These values are chosen to be consistent with our experimental data as presented and discussed in an earlier section.

Interfaces
In FG systems with 2 and more layers, it is necessary to distinguish the strength and weakness of individual interfaces. Figure 10

Top Coat (TC)
Three topcoat materials are considered in this research as briefly discussed here for three architectures for the FG TBC system.

YSZ TC
Yttria-stabilized zirconia (YSZ) is the most widely used materials for TBCs owing to its excellent shock resistance, low-thermal conductivity, and relatively high coefficient of thermal expansion. An amount of 6 wt.%-8 wt.% YSZ is considered as TC and Inconel 718 superalloy as BC. Elastic modulus (E) data in GPa is taken as 40 and 190 respectively for YSZ TC and Inconel BC alloy [1,23,24]. This value is reported to be a strong function of processing technique as well as temperature. The modulus value for YSZ is reported to be as low as 20 GPa at room temperature. E for APS and EBPVD methods is reported to be 90 and 200 GPa, respectively. General rule of mixture is used to estimate the elastic modulus for all BC-YSZ intermediate FG composite layers. As for example, for the layer with equal volume fractions of the two phases, the elastic modulus is 115 GPa, while for 75BC/25YSZ graded layer the same is 152.5 GPa. The density effect is not taken into account because of small density difference between the BC and YSZ phases. The various FG TBC systems studied, layer thickness, and compositional variations are detailed in Table 2. Equation (8) has been used to compute dimensionless K1/K0 values.

Gadolinium Zirconate TC
Traditional YSZ TBC topcoat cannot always meet more demanding requirements at elevated temperatures. It is of great industrial interest to develop a new generation of TBCs to achieve greater insulation to allow higher operation temperatures. Novel ceramics like Gadolinium zirconate (GdZ, Gd2Zr2O7), is being considered as a promising candidate because of lower thermal conductivity, better sintering ability, higher melting point, and phase stability than YSZ. Despite these attractive properties, GdZ with pyrochlore structure has lower fracture toughness than YSZ [32][33][34][35], which could lead to poor thermal cyclic durability. The fracture toughness (K1C) value of GdZ is about 1 MPa·m 1/2 , while for YSZ this ia about 4 MPa·m 1/2 . The elastic modulus of GdZ is reported to lie in the range of 205 to 215 GPa and an average value is used for simulations. YSZ BC Figure 10. Designation of interfaces in three FG TBC systems (2, 4, and 5 layers) along with the base case (YSZ TC over BC). Decreasing layer thickness and the compositional changes with increase in layer numbers is also evident as the TC thickness is kept constant at 360 µm. The number ratios represent the volume fraction of two phases in the layer, namely YSZ/BC.

Top Coat (TC)
Three topcoat materials are considered in this research as briefly discussed here for three architectures for the FG TBC system.

YSZ TC
Yttria-stabilized zirconia (YSZ) is the most widely used materials for TBCs owing to its excellent shock resistance, low-thermal conductivity, and relatively high coefficient of thermal expansion. An amount of 6 wt.-8 wt.% YSZ is considered as TC and Inconel 718 superalloy as BC. Elastic modulus (E) data in GPa is taken as 40 and 190 respectively for YSZ TC and Inconel BC alloy [1,23,24]. This value is reported to be a strong function of processing technique as well as temperature. The modulus value for YSZ is reported to be as low as 20 GPa at room temperature. E for APS and EBPVD methods is reported to be 90 and 200 GPa, respectively. General rule of mixture is used to estimate the elastic modulus for all BC-YSZ intermediate FG composite layers. As for example, for the layer with equal volume fractions of the two phases, the elastic modulus is 115 GPa, while for 75BC/25YSZ graded layer the same is 152.5 GPa. The density effect is not taken into account because of small density difference between the BC and YSZ phases. The various FG TBC systems studied, layer thickness, and compositional variations are detailed in Table 2. Equation (8) has been used to compute dimensionless K 1 /K 0 values.

Gadolinium Zirconate TC
Traditional YSZ TBC topcoat cannot always meet more demanding requirements at elevated temperatures. It is of great industrial interest to develop a new generation of TBCs to achieve greater insulation to allow higher operation temperatures. Novel ceramics like Gadolinium zirconate (GdZ, Gd 2 Zr 2 O 7 ), is being considered as a promising candidate because of lower thermal conductivity, better sintering ability, higher melting point, and phase stability than YSZ. Despite these attractive properties, GdZ with pyrochlore structure has lower fracture toughness than YSZ [32][33][34][35], which could lead to poor thermal cyclic durability. The fracture toughness (K 1C ) value of GdZ is about 1 MPa·m 1/2 , while for YSZ this ia about 4 MPa·m 1/2 . The elastic modulus of GdZ is reported to lie in the range of 205 to 215 GPa and an average value is used for simulations.

Lanthanum Zirconate TC
Lanthanum zirconate (LaZ, La 2 Zr 2 O 7 ) is a pyrochlore structure ceramic material consisting of the ZrO 6 octahedra and La 3+ ions at vacancies at the La 3+ , Zr 4+ and O 2− [36,37]. Lanthanum zirconate is considered as a promising thermal barrier coating (TBC) material as it has high sintering resistance up to its melting point. The Young modulus of Lanthanum zirconate is around 175 GPa. The low thermal expansion coefficient and poor toughness are the main factors that contribute to Lanthanum zirconate having smaller thermal cycling life. The thermal cycling performance can be improved by incorporating a multi-layer system, i.e., Lanthanum zirconate on top of YSZ. LaZ coating exhibits a better bond coat oxidation resistance than YSZ. They also provide an effective protection from calcium magnesium alumino silicate (CMAS) attack and volcanic ash.

YSZ System
Firstly, we take up the simulations with YSZ systems as detailed in Table 2 and in Figure 9a,b. Three types of simulations are identified to be meaningful and carried out to investigate the impacts of FG layers. These three simulation types are with respect to: •

Number of FG layers, •
Compositions of layers, • Thickness of layers.
For monolithic TC materials, the three are interdependent under certain conditions. The simulation details and results are presented in this section.

Number of Layers
The simplest case in the FGM systems is to study the effect of increasing number of layers substituting for YSZ top layer as considered in the conventional TBC system design. For the base case with only one interface (Figures 9a and 10 and Table 2), K 1 /K 0 is computed to be 1.74 using the model. The computed values of K 1 /K 0 for higher layers are included in the last column of the table for the BC/TC interface. These values provide the crack driving force for the respective interfaces in terms of normalized SIF between the adjacent composite layers as the case may be. The higher the K 1 /K 0 value, more is the driving force for cracking, and so more probability for failure. For instance, a value of 1.41 for K 1 /K 0 is obtained for the four-layer system between BC and 75YSZ/25BC composite layer, and so lower driving force. Results of the simulations for normalized SIF (K 1 /K 0 ) with increasing layer numbers is displayed in Figure 11 for all the interfaces in different FG TBC systems. The profile of K 1 /K 0 in interfaces decreases gradually as we move from TC towards BC-TC interface shown in Figure 11. Also, the K 1 /K 0 line lowers as the number of layers increases, implying the beneficial effects of graded layers, i.e., increasing layer numbers offers less driving force. Two interesting observations can be noted here.
Firstly, in all the cases, least value of K 1 /K 0 is observed at the BC/TC interface for a FG TBC system. But the driving force to cracking gradually decreases for the successive interfaces moving away from the topmost YSZ layer. In different words, the interface cracking probability is maximum under a given applied stress for the interface adjacent to the topmost layer (interface 1; highest K 1 /K 0 ) and the probabilities to failure diminish for successive interfaces. Secondly, as the number of intermediate layers increases, the fracture resistance to cracking tend to increases ( Figure 11 and Table 2) irrespective of the interface location. However, the extent of beneficial impacts due to increased layers is much higher around smaller interface numbers (2 or 3 layers) than the same for higher numbers.
The sharp drop in K 1 /K 0 values is evident for interfaces 1 to 2, rather than between higher indexed interfaces. The drop in K 1 /K 0 factor is around 29% and 14% over the interfaces between 1 to 2 and 2 to 3, respectively for a 3-layer system. All K 1 /K 0 profiles taper off for higher indexed interfaces. Figure 11 clearly advocates for higher number of intermediate composite layers with enhanced fracture resistance, but the increasing cost of production for FGM TBC would limit such numbers.

Compositional Effects
However, increase in the number of layers in various FGTBC with fixed BC and TBC layer thicknesses introduces two apparent changes, namely changes in the composition and changes in the layer thickness. For a multilayer system with same layer thickness, the variation in driving force reflects the effect of composition only. As the elastic moduli of successive layers improves with increasing fraction of BC (lesser fraction of YSZ), the driving force given by K1/K0 also tends to enhance. As an example, K1/K0 changes from 1.35 to 1.14 for a 5-layer system as the YSZ contents in layers decreases from 80% to 20%. On the other hand, a comparison among different system is not straightforward as both layer thickness and compositional variation are involved. Further analysis is carried out to separate these two influences.
The impact of compositional gradation on the stability of graded TBC system is studied by keeping other factors constants, namely number and thickness of layers. Figure 12 illustrates the variation of K1/K0 at different five interfaces for a 5-layer FGTBC system. The base case here may be the one with uniform compositional gradation of 0.2, i.e., five layers will have YSZ content as 1-0.8-0.6-0.4 and 0.2. In other cases, the gradation is changed by the ratio of YSZ content in the adjacent layers across an interface. As for instance, for 0.5 ratios, the YSZ contents in the layers are respectively 1-0.5-0.25-0.125 and 0.06. Each graded layer thickness is the same in all cases and it is 72 μ. Interestingly, as the compositional ratio (CR) increases, the FGTBC stability and fracture resistance tend to improve in top layers (interfaces 1 to 3 in Figure 11). In fact, for CR equals to 0.8 and 0.9, the K1/K0 factor is appreciably lower than those for the base case. However, the K1/K0 values shoot up for interfaces on the other side (BC) and so instability enhances. The compositional effect in the instability model is accounted for through elastic modulus of the composite layers. The larger difference in modulus between adjacent layers tends to enhance instability in the FGTBC system. The base case situation with uniform gradation seems to be at the centerline balancing the instability variation over entire compositional gradation. Unlike the other, the uniformly graded FGTBC has the CR varying as 0.8, 0.75, 0.67, and 0.5 for respective layers.

Compositional Effects
However, increase in the number of layers in various FGTBC with fixed BC and TBC layer thicknesses introduces two apparent changes, namely changes in the composition and changes in the layer thickness. For a multilayer system with same layer thickness, the variation in driving force reflects the effect of composition only. As the elastic moduli of successive layers improves with increasing fraction of BC (lesser fraction of YSZ), the driving force given by K 1 /K 0 also tends to enhance. As an example, K 1 /K 0 changes from 1.35 to 1.14 for a 5-layer system as the YSZ contents in layers decreases from 80% to 20%. On the other hand, a comparison among different system is not straightforward as both layer thickness and compositional variation are involved. Further analysis is carried out to separate these two influences.
The impact of compositional gradation on the stability of graded TBC system is studied by keeping other factors constants, namely number and thickness of layers. Figure 12 illustrates the variation of K 1 /K 0 at different five interfaces for a 5-layer FGTBC system. The base case here may be the one with uniform compositional gradation of 0.2, i.e., five layers will have YSZ content as 1-0.8-0.6-0.4 and 0.2. In other cases, the gradation is changed by the ratio of YSZ content in the adjacent layers across an interface. As for instance, for 0.5 ratios, the YSZ contents in the layers are respectively 1-0.5-0.25-0.125 and 0.06. Each graded layer thickness is the same in all cases and it is 72 µ. Interestingly, as the compositional ratio (CR) increases, the FGTBC stability and fracture resistance tend to improve in top layers (interfaces 1 to 3 in Figure 11). In fact, for CR equals to 0.8 and 0.9, the K 1 /K 0 factor is appreciably lower than those for the base case. However, the K 1 /K 0 values shoot up for interfaces on the other side (BC) and so instability enhances. The compositional effect in the instability model is accounted for through elastic modulus of the composite layers. The larger difference in modulus between adjacent layers tends to enhance instability in the FGTBC system. The base case situation with uniform gradation seems to be at the centerline balancing the instability variation over entire compositional gradation. Unlike the other, the uniformly graded FGTBC has the CR varying as 0.8, 0.75, 0.67, and 0.5 for respective layers.

Layer Thickness
An important issue in the effective designing of FGTBC is the dimensional parameter, i.e., the individual layer thickness or fraction of total TBC thickness. The analytical model is used to study the significance of the layer dimension with regard to the instability and integrity. A 5-layers system is considered with uniform compositional gradation of 20% YSZ content across the layers to demonstrate the influence of layer thickness. Figure 13 shows the changes in the K1/K0 values across the interfaces for different layer thickness. As before, the total thickness is assumed to be 360 μm and sum of the fractional thickness (individual to total thickness) is maintained to be unity in all cases. A base case (A) here is the one with the same fractional thickness of 0.2 for all the 5 layers. A steady drop in the K1/K0 may be observed over the first 1 to 3 interfaces and later the plot tends to be flattened. However, for the base case, the first interface exhibits a relatively high value of K1/K0 and so higher crack driving force. A few arbitrary and non-uniform thickness variation (B to E, Figure 13) is assumed in the analysis and results of four such cases are shown in Figure 12. The K1/K0 is seen to drop significantly in the first and second layer with larger layer thickness (case D and E). However, if the thickness becomes too small (0.1 fraction), there is a tendency for K1/K0 to increase (interfaces 4 and 5 in Figure 13. Case E among all appears to yield most optimized results with a smooth variation of K1/K0. A different approach to analyze the effect of layer thickness on K1/K0 is shown in Figure 14 where the adjacent layer thickness ratio is plotted against K1/K0. The effect is shown separately for five

Layer Thickness
An important issue in the effective designing of FGTBC is the dimensional parameter, i.e., the individual layer thickness or fraction of total TBC thickness. The analytical model is used to study the significance of the layer dimension with regard to the instability and integrity. A 5-layers system is considered with uniform compositional gradation of 20% YSZ content across the layers to demonstrate the influence of layer thickness. Figure 13 shows the changes in the K 1 /K 0 values across the interfaces for different layer thickness. As before, the total thickness is assumed to be 360 µm and sum of the fractional thickness (individual to total thickness) is maintained to be unity in all cases. A base case (A) here is the one with the same fractional thickness of 0.2 for all the 5 layers. A steady drop in the K 1 /K 0 may be observed over the first 1 to 3 interfaces and later the plot tends to be flattened. However, for the base case, the first interface exhibits a relatively high value of K 1 /K 0 and so higher crack driving force. A few arbitrary and non-uniform thickness variation (B to E, Figure 13) is assumed in the analysis and results of four such cases are shown in Figure 12. The K 1 /K 0 is seen to drop significantly in the first and second layer with larger layer thickness (case D and E). However, if the thickness becomes too small (0.1 fraction), there is a tendency for K 1 /K 0 to increase (interfaces 4 and 5 in Figure 13. Case E among all appears to yield most optimized results with a smooth variation of K 1 /K 0 . A different approach to analyze the effect of layer thickness on K 1 /K 0 is shown in Figure 14 where the adjacent layer thickness ratio is plotted against K 1 /K 0 . The effect is shown separately for five interfaces in a 5-layer TBC. Again, the compositional gradation is uniform. The TBC-0.8TBC is the first interface while BC-0.2TBC is the fifth interface consistent with the earlier figure (Figure 13). The K 1 /K 0 values decrease with the increase in the layer thickness ratio to a minimum and thereafter increase further. It is clear that irrespective of the interface locations, the minimum K 1 /K 0 values correspond to the layer thickness ratio lying in the range of 1 to 1.2. However, the effect is more sensitive for first interface between TBC-0.8TBC layers than others and is the least sensitive for the other extreme interface. crack driving force. A few arbitrary and non-uniform thickness variation (B to E, Figure 13) is assumed in the analysis and results of four such cases are shown in Figure 12. The K1/K0 is seen to drop significantly in the first and second layer with larger layer thickness (case D and E). However, if the thickness becomes too small (0.1 fraction), there is a tendency for K1/K0 to increase (interfaces 4 and 5 in Figure 13. Case E among all appears to yield most optimized results with a smooth variation of K1/K0. A different approach to analyze the effect of layer thickness on K1/K0 is shown in Figure 14 where the adjacent layer thickness ratio is plotted against K1/K0. The effect is shown separately for five interfaces in a 5-layer TBC. Again, the compositional gradation is uniform. The TBC-0.8TBC is the  Figure 13. Effect of different layer width (fraction of total TBC width) on K 1 /K 0 for 5-layer FGM.
Coatings 2020, 10, x FOR PEER REVIEW 21 of 25 first interface while BC-0.2TBC is the fifth interface consistent with the earlier figure (Figure 13). The K1/K0 values decrease with the increase in the layer thickness ratio to a minimum and thereafter increase further. It is clear that irrespective of the interface locations, the minimum K1/K0 values correspond to the layer thickness ratio lying in the range of 1 to 1.2. However, the effect is more sensitive for first interface between TBC-0.8TBC layers than others and is the least sensitive for the other extreme interface.

Advanced TBC
A number of advanced TBC materials are now being designed and developed to replace YSZ. Alternative structural architectures have also been proposed in recent times. Double-and triple (hybrid)-layered TBC are proposed as top coats to take advantage of the high thermal stability and low thermal conductivity of these materials and at the same time to avoid the drawbacks associated with the relatively lower CTE and toughness. Thus, wide-ranging TBC architectures, designs, and advanced materials are now available for optimizations of combinations. Our current research work with YSZ TC has been further extended here to include both DL and TL TBC systems

FG GdZ System
Significance of promising gadolinium zirconate (GdZ) TBC is discussed in an earlier section. Simulations are also carried out with FG model for the GdZ system as the FG YSZ discussed in an earlier section. As presented in Table 2, YSZ data are replaced by GdZ data and the K1/K0 values are determined for all interfaces as shown in Figure 10. Results of simulations are included in Figure 10. K1/K0 profiles may be seen to be flat as the values for all GdZ FG systems lie in the close range of 1.1 to 1.2, unlike the YSZ systems where the range of K1/K0 extends around 1.3 to 3.4. This implies that the interfaces and layers in GdZ are stronger as the stress intensity levels remain within 10% to 20% of the base level. This is primarily due to higher value of elastic modulus of Gd and quite comparable to the modulus of BC. However, the general trend and impacts of number of layers are observed to be consistent between the two systems. No layer compositional and thickness effects are studied for GdZ systems.

Advanced TBC
A number of advanced TBC materials are now being designed and developed to replace YSZ. Alternative structural architectures have also been proposed in recent times. Double-and triple (hybrid)-layered TBC are proposed as top coats to take advantage of the high thermal stability and low thermal conductivity of these materials and at the same time to avoid the drawbacks associated with the relatively lower CTE and toughness. Thus, wide-ranging TBC architectures, designs, and advanced materials are now available for optimizations of combinations. Our current research work with YSZ TC has been further extended here to include both DL and TL TBC systems

FG GdZ System
Significance of promising gadolinium zirconate (GdZ) TBC is discussed in an earlier section. Simulations are also carried out with FG model for the GdZ system as the FG YSZ discussed in an earlier section. As presented in Table 2, YSZ data are replaced by GdZ data and the K 1 /K 0 values are determined for all interfaces as shown in Figure 10. Results of simulations are included in Figure 10. K 1 /K 0 profiles may be seen to be flat as the values for all GdZ FG systems lie in the close range of 1.1 to 1.2, unlike the YSZ systems where the range of K 1 /K 0 extends around 1.3 to 3.4. This implies that the interfaces and layers in GdZ are stronger as the stress intensity levels remain within 10% to 20% of the base level. This is primarily due to higher value of elastic modulus of Gd and quite comparable to the modulus of BC. However, the general trend and impacts of number of layers are observed to be consistent between the two systems. No layer compositional and thickness effects are studied for GdZ systems.

GdZ-YSZ DL System
A double layer (DL) advanced TBC system is displayed in Figure 9c with GdZ layer sat over the first layer YSZ. Two interfaces are formed between GdZ/YSZ and YSZ/BC layers. These interfaces are unique because of the combinations of two entirely different materials with different modulus value. Here again, the interfaces can be engineered with more layers with different proportions of YSZ and GdZ. Instead, two structural variations are made here with respect to width of bilayers for simulation. As shown in Figure 9c, the width for GdZ and YSZ layers are considered to 300 and 60 µm, respectively, while BC width is considered the same as 150 µ. The simulated K 1 /K 0 values for two interfaces are 0.58 and 1.32, respectively, for Gdz/YSZ and YSZ/BC layers. In spite of large difference in modulus values between YSZ and GdZ (40 and 210 GPa), the K 1 /K 0 turned out to be rather small, and this is due to the layer width difference. GdZ layer is five times the width of the YSZ layer ( Figure 9c). Another simulation is completed with the widths of GdZ and YSZ layers as 200 and 160 µm, keeping total TC width as 360 µm. The K 1 /K 0 values are 0.284 and 3.53 for YSZ/GdZ and YSZ/BC interfaces, respectively.
In order to find a more scientific and logical explanation for dependency of K 1 /K 0 values, an attempt is made here to separate out the effects of the two ratios as may be seen in Equation (8), i.e., elastic modulus (E) and adjacent layer width (W). Data collected from the simulations point to the various FG YSZ TBC cases, as discussed earlier with respect to Figure 7 and Table 2. The K 1 /K 0 values reduce sharply with increase in E ratios between adjacent layers as displayed in Figure 15. In other words, the smaller the E-ratios, the higher is K 1 /K 0 . As the ratios of two adjacent layers tend to increase or match in properties, the K 1 /K 0 values drop. Signifying the lower stress and stress intensity at the interfaces. On the other hand, increase in K 1 /K 0 ratios is evident with increase in the values of width ratio, linearly with a small slope. Figure 15 displays the impact of width ratios between layers for a small range, i.e., from 0.2 to 0.55 only. Ignoring the large number of simulation points to a width ratio of 0.50 and wide range of K 1 /K 0 . Due to the effects of E-ratios, the width effect seems to be only moderate, unlike the strong effects of E-ratios. Table 3 includes simulation results for GdZ-YSZ DL systems. It is very evident that the E ratios are primarily responsible for high values of K 1 /K 0 when the widths among the layers are nearly the same.  GdZ-LaZ-YSZ-BC 120-120-120-150 0.96 0.264 2.71 YSZ-GdZ-YSZ-BC 120-120-120-150 6.06 0. 22 2.70 ratio, linearly with a small slope. Figure 15 displays the impact of width ratios between layers for a small range, i.e., from 0.2 to 0.55 only. Ignoring the large number of simulation points to a width ratio of 0.50 and wide range of K1/K0. Due to the effects of E-ratios, the width effect seems to be only moderate, unlike the strong effects of E-ratios. Table 3 includes simulation results for GdZ-YSZ DL systems. It is very evident that the E ratios are primarily responsible for high values of K1/K0 when the widths among the layers are nearly the same. Figure 15. Impact of the ratios of E and W on the normalized stress intensity factor, K1/K0.

Hybrid System
Hybrid TBC architecture has been considered recently to replace conventional monolithic YSZ TBC. Hybrid system will offer far more flexibility in the optimization of performance. One trilayer hybrid system is displayed in Figure 9d with three materials, namely GdZ, LaZ, and YSZ. GdZ is chosen for low thermal conductivity, YSZ for high fracture toughness, and LaZ is the choice for high elastic modulus and high sintering resistance. Simulations with the FG TBC models is carried out for the hybrid system. Table 3 displays the results of simulations for a few arbitrary TL systems along with layer width. The K 1 /K 0 values are all small because of matching E ratios between GdZ and LaZ and the interface between YSZ and BC reflects a value of 2.71, as the E values are quite different. This argument is also reflected in the values of K 1 /K 0 for YSZ-GdZ-YSZ-BC system as reported in Table 3.

Conclusions
Experimental and analytical research is presented here with an objective of modeling the damaging process and characterizing fracture resistance of functionally graded YSZ TBC and other advanced TBC systems. The following conclusions emerge from the work.
An excellent agreement in baseline TGO oxide growth measurement was observed by area-based and point to point techniques. TGO growth is found to have a linear function of time under thermal exposure till failure occurs. On the other hand, the damaging and induced crack length measured follows exponential relation during the oxide growth. The primary (i.e., TGO) and secondary physical damage (i.e., cracks) in TBC system under thermal cycling are also related exponentially. The plasma sprayed coating methods with low oxygen and air show no marked difference in the oxide growth and damage behavior.
Progressive growth of both the physical damages during thermal cycling is modeled identifying the crack nucleation sites and changing TGO profiles. The cracking and interface separation mechanisms are influenced by shape and curvatures at the interface boundaries. Cracks are mostly observed at and around sharp locations and ridges where higher tensile stress is present. For thermal cycling up to 100, the interface separations are observed only at the TGO-TBC boundary. But longer cycling time (430 and more) induces separation and cracking also at TGO-BC interface, and also around less vulnerable locations. However, as the TGO thickness increases with continued thermal cycling, the extensive transverse cracking causes the damage. The experimental work confirms the presence of cracks in the TBC system in the size range of two to three micrometers, even at the onset of thermal cycling (around ten cycles).
An isostrain-based simplistic fracture mechanical model developed for conventional TBC system and extended to functionally graded TBC. The normalized stress intensity factor, K 1 /K 0 , is found to be an effective design parameter in evaluating the fracture resistance of the interfaces in FGM TBCs. A few realistic design scenarios are considered along with a hypothetical case where YSZ-YSZ interface represents the worst case and BC-YSZ conventional system as the base case. The K 1 /K 0 values are determined to be 1.47 and 1.83 for the base and worst case, respectively. An increase in fracture resistance by around 18 percent can be achieved by introducing nine graded layers of uniform thickness in between the bond coat and top YSZ coat.
Bilayer and hybrid TBC systems consisting of gadolinium zirconate, lanthanum zirconate, and YSZ are characterized using the proposed fracture mechanics model. Fracture characterizing parameter, K 1 /K 0 , lie in the close range of 1.1 to 1.2 for FG GdZ systems unlike the YSZ, where the range is 1.3 to 3.4. The elastic modulus of the layers across an interface in a TBC system are primarily responsible for values of K 1 /K 0 when the widths effects seem to have a much smaller role. Both bilayer and hybrid systems signify improved interface fracture resistance.
The fracture mechanics modeling framework presented is the first initiative and opens up new research directions towards more useful and refined modeling for more realistic situations. Some important effects need to be incorporated, namely thermal, material anisotropy, nonlinear stress, crack interaction, and stress shielding. Future modeling and simulation should involve effects of thermal cycling, thermal expansion mismatch, temperature dependent modulus, key microstructural features, interface roughness, elasoplastic nonlinear behavior for both plasma, and EBPVD processed TBCs.

Conflicts of Interest:
The authors declare no conflicts of interest.