COMPUTATIONAL EVALUATION OF THERMAL BARRIER COATINGS

In the power plant industry, the turbine inlet temperature (TIT) plays a key role in the efficiency of the gas turbine and, therefore, the overall—in most cases combined—thermal power cycle efficiency. Gas turbine efficiency increases by increasing TIT. However, an increase of TIT would increase the turbine component temperature which can be critical (e.g., hot gas attack). Thermal barrier coatings (TBCs)—porous media coatings—can avoid this case and protect the surface of the turbine blade. This combination of TBC and film cooling produces a better cooling performance than conventional cooling processes. The effective thermal conductivity of this composite is highly important in design and other thermal/structural assessments. In this article, the effective thermal conductivity of a simplified model of TBC is evaluated. This work details a numerical study on the steady-state thermal response of two-phase porous media in two dimensions using personal finite element analysis (FEA) code. Specifically, the system response quantity (SRQ) under investigation is the dimensionless effective thermal conductivity of the domain. A thermally conductive matrix domain is modeled with a thermally conductive circular pore arranged in a uniform packing configuration. Both the pore size and the pore thermal conductivity are varied over a range of values to investigate the relative effects on the SRQ. In this investigation, an emphasis is placed on using code and solution verification techniques to evaluate the obtained results. The method of manufactured solutions (MMS) was used to perform code verification for the study, showing the FEA code to be second-order accurate. Solution verification was performed using the grid convergence index (GCI) approach with the global deviation uncertainty estimator on a series of five systematically refined meshes for each porosity and thermal conductivity model configuration. A comparison of the SRQs across all


INTRODUCTION
The future of increased efficiency in gas turbine-driven power plant energy generation is heavily dependent on the increased temperature of the turbine inlet temperature (TIT) [1][2][3][4]. In spite of the pursuit of higher TIT levels, virgin metal material in the turbine components are susceptible to a variety of aggressive and high-consequence degradation or failure modes. Such issues include accelerated thermal creep, material degradation due to oxidation, and cycle fatigue [5,6]. Thermalrelated corrosion issues arise at various temperature levels anywhere in the range of 650-1700 °C [2,[7][8][9][10][11].
To mitigate thermal issues, thermal barrier coatings (TBCs) are used to protect the virgin material. Such coatings, including spray technologies are typically comprised of some form of ceramic composite material [12][13][14]. These composite materials are porous in nature and provide significant protection to the metal turbine components to allow for better overall performance and efficiency.
Understanding the thermophysical properties of TBC materials is critical to the design and advancement of gas turbines for plant power generation for current and future installations. As systems continue to produce higher TIT levels, the thermal performance of advanced TBCs becomes more critical. This investigation evaluates the effective thermal transport in simplified models of TBC porous structures, specifically two-phase porous media, as an initial study into TBC material nano-structure effective thermal response.

SYSTEM DESCRIPTION
In order to analyze the thermal transport phenomena within the TBC as a part the authors' efforts, a two-dimensional system is constructed, consistent with the illustration given in Figure 1. For this study, a unit half-cell is given with height, L, and width, W, and is oriented in the Cartesian x-y plane with unit thickness out-of-plane. The bulk matrix material is defined as the first material phase with thermal conductivity k 1 , where the pore phase is defined to have thermal conductivity k 2 . The pore pattern is circular, centered in the unit cell, with radius, R. Both the x=0 and x=W boundaries are given periodic boundary conditions for system symmetry, and the y=0 and y=L boundaries have enforced hot and cold boundary temperatures, T H and T C , respectively. By enforcing Dirichlet boundary conditions at the hot and cold boundaries, heat flow-per unit length-is induced in the positive y-direction, perpendicular to each boundary, indicated by Q H and Q C , respectively.

Figure 1. Unit half-cell nano-porous structure
Of interest in this work is the effect of relative pore size and relative pore thermal conductivity on the overall effective thermal conductivity, k eff , of the TBC. Thus, the dimensionless parameter, porosity, α, is used to define the relative pore size, where The k eff value is defined by overall heat flow through the system and the enforced boundary temperatures such that where the domain is assumed to have unit thickness out-ofplane. In theory, Q H =Q C , but because this study is numerical in nature, the average of the two heat flow values is used in Eq.
(2) to define the overall induced heat flow through the domain. To account for variable material properties of the matrix and porous phases in different TBCs, the ultimate system response quantity (SRQ) of interest in this analysis is the dimensionless effective thermal conductivity, k*, which is merely the ratio of x y the effective thermal conductivity with the matrix material thermal conductivity, simply evaluated as * = / 1 .
This study will show the effects on the SRQ of varying the k 2 -to-k 1 ratio along with varying the α, using numerical verification techniques to quantify uncertainty in that value. Ratios in the two-phase thermal conductivities are evaluated from 0 to 2, and α values are varied from 5% to 65%.

DISCRETIZATION
In order to determine the SRQ, the temperature distribution, T, and resultant heat source terms, S, must be solved to satisfy the governing partial differential equation (PDE), the heat equation, as expressed in Eq. (4), where k is the thermal conductivity.
To numerically discretize the heat equation, unstructured triangular meshes were generated over the domain within each material phase. To accommodate the verification approaches used in this study, systematically refined meshes were generated for each porosity level, an example of which is shown in Figure 2a through Figure 2e, for mesh number, H, 5 through 1, respectively. Likewise, the finest meshes for a selection of porosity models of 5%, 25%, 45%, and 65% are given in Figure 3a through Figure 3d, respectively. Figure 4. Each H th mesh for a given porosity level is then described with a characteristic mesh size h H . If N H,t is the total number of triangles in mesh H, then

Each t th triangle is configured with three vertices such that the i th vertex is located at (x t,i ,y t,i ) with temperature T t,I and area A t , is shown in
Thus, systematic mesh refinement here approximately halves the characteristic mesh size from mesh H+1 to mesh H. The domain discretization is used with a second order accurate Galerkin Finite Element Method (FEM) approach which is employed to solve Eq. (4) over the entire domain [15]. In short, Eq. (4) can be solved over each element, individually, using the discretized PDE matrix equation where Ν is the linear shape function row vector used to interpolate field variables across a single element, Κ is the 2x2 thermal conductivity matrix, T is the vector of triangle vertex temperature values, q n , is the boundary normal heat flux on a triangle boundary, Ω is the element domain, and Γ is the element boundary. Furthermore, the different parts of Eq. (6) and

SOLVER
The meshes employed in this study were generated using open-source Gmsh software [16]. The custom code used for this study was written in Fortran as a generic two-dimensional FEM heat transfer solver. A basic successive over-relaxation method was used to iteratively update solutions to Eq. (9) until some residual level exit criteria was met.
The residual vector, ρ, for the system is computed as where the L ∞ norm residual, ρ L∞ , defined as A ρ L∞ value of 10 -8 was used as exit criteria for solution completion.

CODE VERIFICATION
In order to verify the formally second order accurate solution, the Method of Manufactured Solutions (MMS) was employed [17,18]. The MMS approach allows for a userdefined temperature distribution to be selected, T MMS , where the operator-in this case defined by the PDE in Eq. (4)-acts on the manufactured solution. The resulting source term, S MMS , falls out of the operator, where Thus, by applying S MMS and the boundary conditions that satisfy T MMS to the domain, the solved system of equations should converge to T MMS with mesh refinement. Here, T MMS is defined as and is applied as a boundary condition at y=0 and y=L as Dirichlet conditions. It follows that the first partial derivatives are given to be
For this study, the L ∞ norm metric is used on temperature error, ε MMS,H , over all N v vertices in mesh H, where based on the solution temperature, T i , of the i th vertex. The observed order of accuracy for mesh H, p O,H , is determined by

SOLUTION VERIFICATION
The numerical error, U num , for the SRQ, k*, is approximated in this method by performing solution verification using a modern version of a more traditional grid convergence index (GCI) [18,19] approach that incorporates a global deviation estimator for the solution [20]. This approach uses an empirical trend to scale a factor of safety, FS, for U num based on how closely the convergent observed order of accuracy is to the formal order of accuracy, thus acknowledging how close the solution is to the asymptotic solution region.
For this problem, where the SRQ is k*, the global deviation parameter, Δp, is computed as with p f being the formal order of accuracy of the solution method which, in this study, is 2. Lastly, the global order of accuracy, p*, used to compute the global deviation GCI FS is * = − such that = 3 − 1.9( * / ) 8

CODE VERIFICATION
Using the MMS solution and conditions described in Eq. (13) through Eq. (15) to perform the MMS verification, the trends of logarithmic error against logarithmic mesh size are shown in Figure 5, where the trend of a formal second order accurate solution is also shown in qualitative comparison. Likewise, Figure 6 shows the convergence of the observed order of accuracy with respect to mesh number.

SOLUTION VERIFICATION
As mentioned previously, k* is computed using the average numerically-produced values of Q H and Q C , as described by Eq. After performing the solution verification procedures described with Eq. (18) through Eq. (23), the following tables show the computed global deviation parameters and ultimate U num approximation for the finest mesh in each material configuration model. Note that the U num values given in these tables are expressed as percent values relative to the presented fine mesh k* values. All U num values are obtained to be less than 0.2% of the fine mesh k* values.

CHARACTERISTIC TRENDS
Compiling all of the k* from Table 8 through Table 14 shows a family of similar curves-roughly quadratic in formrelating k* to both α and k 2 /k 1 , as illustrated in Figure 7. In Figure 7, it is clear that an increase in the pore material thermal conductivity leads to an increase in overall effective thermal conductivity of the composite material. Also, it is apparent from Figure 7 that an increase in pore size amplifies the relative effect of k 2 on k*. All curves intersect at (k 2 /k 1 =1.0,k*=1.0), due to the fact that the condition of k 2 /k 1 =1.0 implies that the composite material behaves as a single homogeneous material from the perspective of steadystate thermal transport. Note that error bars on the k* values are indistinguishable due to their relatively low values.

FUTURE WORKS
The study performed in this work has the potential to be extended to describe the effective thermal behavior in an analytical solution with prescribed uncertainty bounds. Likewise, the effects of pore topology and geometry on effective thermal responses of composite thermal material could be investigated as well as obtaining optimal thermal conductivity for microporous insulations. Similarly, such a study could have implications on the uncertainty incurred by mesh geometry. Further investigations could be made into analysis of transient effects or the implications of more than one pore material type in a given matrix.

CONCLUSION
This paper has presented a numerical approach to capturing the effective thermal transport behavior of a two-phase porous material with respect to pore phase thermal conductivity and pore size. It was shown that the employed FEM code was roughly second order accurate using the MMS code verification approach, and numerical uncertainty was estimated using the global deviation estimator GCI method. A trend was illustrated, tying the effective thermal conductivity of the overall medium to both relative pore size and relative pore thermal conductivity.