Computer-aided design and prediction of the functional properties of double-layered thermal barrier coatings NiAl-YSZ

An approach to the stochastic simulation of the formation process of the lamellar structure of double-layered thermal barrier coatings (TBCs), including high-velocity spraying of NiCrAlY bond layer and atmosphere plasma spraying of YSZ top layer, is proposed. The approach is based on the developed and experimentally verified theoretical foundations which could be used for rapid and sufficiently accurate prediction of thickness and diameter of metal and metal oxide splats depending on the key physical parameters (KPPs): temperature, velocity and size of particle, and substrate temperature. The results are presented, which concern the development of the computational algorithm and the program code for simulation the process of stochastic stacking of splats in the coating with regard for the topography of its surface, which varies dynamically at the spraying, as well as the formation of lamellar structure and porosity of the coating. Results of stochastic simulation of lamellar structure TBCs and its characteristics (microstructure, porosity, roughness) are given in conditions typical for atmospheric plasma spraying (APS) of YSZ top layer and for detonation (D-Gun Spraying) and supersonic flame spraying (HVOF) of bond layer. Comparison of characteristics of coatings for two model splats morphologies at their stochastic stacking on the sprayed surface is presented. The results of numerical experiments are presented through the examples of TBCs thermal spraying, illustrating the performance of the developed computational technology.


Introduction
Thermal spraying is widely utilized to produce thermal barrier coatings (TBCs) and chemical reactive coatings [1].
The TBCs are double-layered coatings consisting of MCrAlY (M=Ni, NiCo) dense metallic bond layer and porous low thermal conductivity 7wt.%Y 2 O 3 -ZrO 2 (YSZ) top layer. The bond layer has to provide the oxidation resistance and enough bond strength of the ceramic top layer. In turn, the top layer has to provide a sufficient thermal insulation and compatibility with bond coat and chemical stability at high temperatures up to ~ 1500 K [1].
Taking into account the significant differences in the functional purpose of the bond and top layers, as well as in the characteristics of the used powder materials, the spraying of TBCs is carried out in two stages. High-speed methods (HVOF, LPPS or D-gun spraying) and powders of the type MCrAlY with a particle size in the range of 10-40 μm are used for the spraying of a bond layer with a thickness of ~ 100-150 μm, while for the spraying of a top layer with a thickness of ~ 200-400 μm, as a rule, APS method and YSZ powders with a particle size in the range of 40-100 μm are used As is known [1], at the gas-thermal spraying (GTS), even with the maximum possible loading of the technological stream (high-temperature jet) with a powder material, the particles of which undergo a complete Since experimental study is difficult to obtain the optimum operating conditions, a computational experiment is supposed to be an effective method for optimizing plasma spraying process.
The reliability of the simulation of sprayed coating structure is largely determined by the adequacy of the physicomathematical model used and the accuracy of its numerical implementation in calculating the stacking of individual splats for given values of the KPPs at the impact of melt droplets with the substrate. It should be noted that the detailed numerical simulation of the processes in the formation of single splits, even on a smooth substrate, using a complete system of Navier-Stokes equations is associated with significant computer time costs commensurate, and in some cases exceeding the calculation time, actually, of a dust-laden technological flow. An alternative is to use rather simple, experimentally tested physical models, theoretical solutions or semiempirical dependencies that allow predicting the main parameters of splats, especially their thickness and diameter over a wide range of KPPs, with sufficient accuracy for practice, which significantly (by several orders of magnitude) accelerates the computational procedure of stochastic simulation of the formation of a layered structure of gas-thermal coatings.
The purpose of this study is the further development of exactly such a computational technology, proposed by us in [2,3,4], for the layered structure simulation of double-layered TBCs, notably, "metallic sublayer -ceramic upper layer".

Simulation of TBC formation
A major weakness of TBC systems is the interface between the metallic bond layer and the top ceramic thermal barrier layer. At this interface an in-service degradation is observed often leading to a macroscopic spallation of the ceramic thermal barrier layer [5]. Therefore, an important aspect in predicting the properties of simulated double-layered TBCs is the possibility of calculating not only the adhesion bonding strength of a bond layer with a substrate, but also the possibility of calculating the cohesive bonding strength of a bond layer with top layer. The software complex developed by the authors of this article has the ability to calculate the characteristics of adhesion and cohesive strength.

Summary of the software complex for the simulation of splats and coatings
The composition of the subsystems of the software package intended for simulation of gas-thermal coatings is shown in figure 1. Let's briefly consider the possibilities of two main subsystems ("SPLAT" and "COATING" subsystems). 3 Subsystem "SPLAT" with the use of the subsystem "DATA BASE" allows to calculate the geometric characteristics of splats for a given set of values of "key physical parameters (KPPs)" of particles, and also to determine the boundaries of the ranges of values of KPPs at which the "model" form of the splats remains integral (unbroken). This is an important aspect for tuning the technological spraying regime of coatings, which are formed when fixing splats on a sprayed surface [2,3,4]. The basis for calculating the characteristics of splat is experimentally tested theoretical self-similar solutions for the formation of a single splat that requires the calculation of a three-dimensional nonstationary boundary value problem with a free boundary for the Navier-Stokes equations in conjunction with the equations of conjugate convective-conductive heat exchange and phase transformations in a spreading particle, and sometimes -in the substrate [6]. In the process of splat formation on a sprayed surface [6], the spraying conditions (KPPs of the impacting melt particle with the surface and KPPs of the surface itself) are taken into account, which determine one of the possible 4 scenarios for the splat formation: (a) spreading and simultaneous solidification of the droplet on the solid base; (b) spreading and simultaneous solidification of the droplet, and local submelting of the base at the contact spot with the droplet; (c) spreading of the droplet over the solid base surface, and subsequent cooling and solidification of the spread layer; (d) spreading of the droplet accompanied with simultaneous local submelting of the base, followed by subsequent cooling and solidification of both.
The scenario corresponding to the spraying conditions makes it possible to calculate the thickness and diameter of the splat [6,7]. Some additional important aspects on the splat formation are set forth in [7].
The software package in the process of simulation the stacking of splats on the sprayed surface makes it possible to use one of two model "morphologies" of the splat: figure 2(а)in the form of a thin cylinder (disk); figure 2 (b) -in the form of a "smoothed" splat (the step of digitizing the X and Y coordinates, corresponding to the figure, was intentionally chosen at simulation not small, so the gradation of the shape profile of the smoothed splat is quite noticeable).  As can be seen from figure 3(b), the pores in the core zone of splat (marked on the right in the figure with a dark gray background, and the right-hand part of the splat periphery not shown) are not formed due to the considerable force of pressure head when the droplet collide with the substrate, but the pores increases in the direction of the splat periphery (on the left on figure 3(b)). The absence of pores in the contact zone of the splat core with a sprayed surface (radius of the splat core с0 ≅ 1.1 * ) is provided by the first algorithm for recognizing and filling depressions on the surface in the contact zone with the splat core. Further, another algorithm is used for recognizing and filling depressions behind obstacles on the surface in the contact zone with the peripheral annular part of the splat. Then, a wave algorithm is used to scan the sprayed surface with a variable topography to prepare an array of "reference vertices" with the consequential approximation of the bottom surface of the splat by means B-spline surfaces. In the procedure for simulation the stacking of a splat on the surface, the XY plane is discretizated with certain steps Δx and Δy (usually Δx = Δy) in the form of a twodimensional grid of nodal points whose indices coincide with the indices of the elements of the twodimensional array of "reference vertices". The values of these elements of the array are the zcoordinates of the "reference vertices" of the sprayed surface in "floating-point arithmetic" without any discretization along Z axis.
From the "residual" volume of the liquid melt of the spreaded drop, after the above-mentioned filling of the depressions, the thickness h s of the splat is corrected. Spline approximation based on the "reference vertices" array calculates the Z-coordinates of the bottom surface of the splat corresponding to the topography of the sprayed surface, taking into account the angle and trajectory of the melt flowing into the depressions of the sprayed surface [7]. The degree of flowing of the melt particle into the depressions of the surface is greater the higher the velocity of the particle's collision with the surface (head pressure of the particle is proportional to the square of the velocity). Consequently, as the normal component of the velocity of the spraying particles increases, the porosity should decrease, which is confirmed by the results of computational experiments and experimental data [7,8]. Taking into account the chosen morphology of the splat (figure 2) and the corrected thickness h s of the splat, the Z-coordinates of the top surface of the splat are calculated, the relief of which is similar to the bottom surface ( figure 3(b)). The Z-coordinates of the bottom and top surfaces of the splat corresponding to the X and Y coordinates of the nodal points of the grid in the XY plane, as a result of the process of stacking the splat on the sprayed surface are stored in the corresponding file on the external disk on the basis of which in Subsystem "Functional characteristics of coatings" the coating characteristics (such as porosity P%, roughness Ra, adhesion strength of bond of the coating with substrate, etc.) are calculated.
The angle and the flowing trajectory of the melt into the surface deepening is calculated according to the expression obtained from certain model concepts on the splat morphology and the physical foundations of the spreading and solidification of the melt particle on the sprayed surface [7]. To calculate the flowing trajectory, we use the empirical constant C 0 , which was determined in calibration computational experiments on the basis of "reference" data from [8], corresponding to the detonation spraying of aluminum oxide coatings from α-phase Al 2 O 3 powder (M28 corundum) with a porosity of 0.6%. At the calibration stage at the variation of a parameter C 0 , the use of the cylindrical splat model does not yield acceptable results (the porosity of the simulated coatings is on the order of 7-8%, which is much higher than the 0.6% value given in [8]). The agreement with the experimental value of 0.6% porosity is obtained only for the model of the smoothed splat (i.e., at the variation of a parameter C0 for a given set of KPPs, the calculated values of porosity changes from one percent to tenths of a percent). The value of the constant C 0 obtained by the calibration is used by the authors as a universal value of this constant for simulation not only ceramic, but also metal-ceramic and metal coatings.

Spreading factor of metal alloy splats
As noted in Subsection 3.1, subsystem "SPLAT" allows you to calculate and monitor KPPs for the stable formation of splats (ensuring their unbroken form). Therefore, computational experiments were performed, as a result of which the spreading factor ̅ = ⁄ was calculated, where and are the diameters of the splat and the initial melt particle. The KPPs values at which the spreading factor ̅ does not exceed 4.5 for metallic splats determine conditions for the stable splat formation and, in general, for the stable formation of the bond layer in the thermal barrier coating [9]. The results of the computational experiments for determining the spreading factor ̅ are given below ( figure 4). To simulate the simulation of the first layer (of the bond layer) from the NiAl metal powder, we used a set of KPPs: particle diameter D p0 = 20 μm, substrate temperature T b0 = 400 K, substratenickel superalloy Inc601-2, coating thickness H = 100 μm. Particle velocities: U p0 = 500, 600, 700, 800 m/s. Analyzing the results of computational experiments and taking into account the limiting value ̅ = 4.5 for the spreading factor of metallic particles, we can conclude that for NiAl particles forming the "bond layer" in TBCs, the particle temperature should not exceed from 1950 to 2000 K when sprayed in HVOF regime.

Spreading factor of metal oxide splats
Metal oxide splats were used to form the second layer (the "top YSZ layer") in the thermal barrier coating. The KPPs values, at which the spreading factor ̅ does not exceed 5 for metal oxide splats, determine the conditions for the stable splat formation and, in general, for the stable formation of the "top YSZ layer" in the thermal barrier coating [10]. In physical reference books and scientific literature there is practically no data on the thermophysical and other properties of the 7wt.% Y 2 O 3 -ZrO 2 (YSZ) particles, but data for zirconium dioxide ZrO 2 are known, which occupies 93% of YSZ material. Therefore, reference data corresponding to zirconium dioxide ZrO 2 is often used for YSZ particles in practical simulation.
The results of computational experiments to determine the spreading factor ̅ are given below ( figure 5). This is the dependences of spreading factor ̅ vs. the particle temperature T p0 and particle velocity U p0 . The diameter of the YSZ particles and the temperature of the substrate (the first NiAl layer) were fixed with the values D p0 = 60 μm and T b0 = 600 K.
Analyzing the results of computational experiments and taking into account the limiting value ̅ = 5 for the spreading factor for metal oxide particles, we can conclude that for YSZ particles forming the "top YSZ layer" of TBCs, the particle temperature should not exceed 3000 K when sprayed in APS regime.  Figure 5. Dependences of spreading factor ̅ of YSZ particles on the metal sublayer NiAl vs. the particle temperature Т р0 and particle velocity U р0 in the APS spraying regime (D p0 = 60 μmdiameter of YSZ particles, T b0 = 600 Ktemperature of metal sublayer NiAl): Analyzing the results of computational experiments and taking into account the limiting value ̅ = 5 for the spreading factor of metal oxide particles, we can conclude that for YSZ particles forming the "top YSZ layer" in TBCs, the particle temperature should not exceed 3000 K when sprayed in APS regime.

Porosity and roughness of simulated bond layer
It is well known that the microstructure and roughness of the bond layer are significantly affected by both the composition as well as powder size. In addition, the interface roughness (surface roughness) of the bond layer in turn is one of the most important factors that affect the lifetime of TBC systems [11].
To model the first layer (the "bond layer") from the NiAl metal powder, we used a set of KPPs: D p0 = 20 μm, T b0 = 400 K, the substrate was a nickel superalloy Inc601-2, the thickness H for the sprayed "bond layer" in the HVOF regime was set as H = 100 μm. For the HVOF spraying regime, the following particle velocities were used: U p0 = 500, 600, 700 and 800 m/s. The values of particle temperatures were chosen as follows: T p0 = 1950, 2050, 2150, 2250 and 2350 K. Based on the results of simulation of coating clusters (base 2х2 mm), the porosity and roughness of the bond layer were calculated, the values of which are shown in tables 1 and 2.
On the basis of the conclusion given at the end of subsection 2.2, the optimum spraying regime of the "bond layer" NiAl with the corresponding set of KPPs was determined: particle temperature T p0 = 1950 K, particle diameter D p0 = 20 μm, particle velocity U p0 = 500 m/s, temperature of the substrate T b0 = 400 K. The thickness of the layer is given by the size H = 100 μm, the size of the substrate is L c = 2 mm (the same in two dimensions). Using the cylindrical model of the splat in computational experiments, the porosity of the first layer turned out to be of the order of 6.1%, i.e. which is significantly higher than the experimental values of about 1%, and the roughness is of the order of 2.8 μm. Using the model of a smoothed splat in computational experiments, the porosity of the first layer turned out to be of the order of 1% (correct value), and the roughness was of the order of 2.2 μm.

Porosity and roughness of simulated top YSZ layer
As a basis for the application of the "top YSZ layer", a "bond NiAl layer" was selected with optimum sputtering conditions specified in the subsection 2.4.
As a "top YSZ layer", a thermal barrier layer of thickness H = 300-350 μm, obtained by plasma spraying (APS) on a sublayer of NiAl with a temperature T b0 = 600 K of YSZ particles with a diameter of D p0 = 60 μm was considered. During the simulation of the "top YSZ layer" the following parameters of the YSZ particles were used: a row of particle temperatures (in Kelvin) T p0 = (3000, 3100, 3300, 3400, 3500) and a row of particle velocities (m/s) of particles U p0 = (100, 150, 200, 250). Based on the simulation results of coating clusters of (substrate 2х2 mm), the porosity and roughness of the coatings were calculated and graphical dependencies were constructed (figures [6][7][8][9]. The graphs of the porosity dependences P% of the top YSZ layer of the thermal barrier coating NiAl-YSZ vs. the particle velocity U р0 in the APS regime for particle temperatures Т р0 above 3000 K are not given, since they are similar to the graphs corresponding Т р0 = 3000 K and reflected on figure  7. From the analysis of figures 6 and 7 it follows that the model of a smoothed splat, in contrast to the cylindrical model, correctly reflects the decrease in porosity with increasing of particle temperature Т р0 and particle velocity U р0 .  Similarly, the graphs of the roughness dependences R a of the top YSZ layer of the thermally barrier coating NiAl-YSZ vs. the particle temperature Т р0 and the particle velocity U р0 in the APS spraying regime are not all presented, since they are similar to the examples of the dependencies shown in figures 8 and 9.  From figures 8 and 9 it follows that the roughness Ra increases with an increase in the temperature T p0 and the particle velocity U p0 by no more than 1 μm, and reaches 7 μm. The structure of porosity and roughness of the simulated double-layered thermal barrier coatings NiAl-YSZ (figure 10) is consistent with the data reflected in figures 6 -9.

Conclusions
The results of the computer-aided design and prediction of the functional properties of double-layered thermal barrier coatings NiAl-YSZ using the smoothed splat model and presented in section 2 are in good agreement with the experimental data, for example, with the data of [12]. In paper [12] a powder (7)(8) wt.% Y2O3-ZrO2 with abbreviated designation YSZ and APS spraying regime (particle velocities from 150 to 250 m/s) was used. The nominal size of YSZ particles was varied from 44 μm to 66 μm [12]. With increasing particle size, the porosity of coatings linearly increased [12] from 10.1 to 12.7%, which agrees well with the authors results of this paper, reflected in figure 7(b).
The computational technology and the corresponding software complex developed by the authors of this article have made it possible to establish a relationship between the spraying conditions and the properties of the simulated coatings.
The results of numerical experiments are presented through the examples of simulated doublelayered thermal barrier coatings, illustrating the performance of the developed computational technology.