Abstract
This paper reports two numerical simulation methods for modeling displacement instabilities around a surface groove in a metal substrate used in nuclear power plant. The amplitude change in the groove, the downward displacement at the base node, and the groove displacement at the periphery were simulated using ABAQUS to compare the results from two methods, as well as the tangential stress in the elements at the groove base and periphery. The comparison showed that for the tangential stress two methods were in close agreement for all thermal cycles. For the amplitude change, the downward displacement, the groove displacement, and the stress distribution, the two methods were in close agreement for the first 3 to 6 thermal cycles. After that, inconsistency increased with the number of thermal cycles. It is interesting that the thermal cycle at which the discrepancy between the two methods began to occur corresponded to a thermally grown oxide (TGO) thickness of 1 μm, which showed the accuracy of the present work over the classic method. It is concluded that the present work’s numerical simulation scheme worked better with a thinner TGO layer than the classic method and could overcome the limitation of TGO thickness by simulating any thickness.
1. Introduction
To advance the rapid development in nuclear energy and to meet the requirements for the future energy, ten countries across the world including Canada, France, and Japan have agreed on a framework to significantly enhance international cooperation in research for a future generation of nuclear energy systems, known as Generation IV. The Generation IV systems are intended to offer significant advances in the sustainability, safety and reliability, economics, proliferation resistance, and physical protection. The expected time for its wide application around the world sets about the year 2030 and its technology goal mainly involves sustainable energy generation with extension of nuclear fuel supply, the minimization and management of nuclear waste, the safe and reliable operation with improved accident management, competitive costs and financial risks of nuclear energy systems, and control and security of nuclear material and facilities [1]. The gas-cooled fast reactor (GFR) system, classed as a Generation IV reactor, maximizes the usefulness of uranium resources by breeding plutonium and can contribute to minimizing both the quantity and the radiotoxicity of nuclear waste by actinide transmutation in a closed fuel cycle. The helium cooled reactor operates at an outlet temperature of 850°C and uses a direct-cycle, helium turbine for electricity (42% efficiency at 850°C) and processes heat for thermochemical production of hydrogen. Protective coatings are visualized to protect various parts of the system and also protect the system in extreme cases where the functional temperature can increase up to 1250°C and there is depressurization from 70 bar to atmospheric pressure. Such coatings must withstand high temperature, depressurization, and specific conditions of wear linked to erosion by high-speed (about 280 m/s) helium gas flow. As an alternative solution, thermal barrier coatings (TBCs) are widely employed at ultra-high-temperature metallic components of gas turbines and high-temperature gas-cooled reactors (HTGRs) in nuclear power plants. Their purpose is to protect the metallic components from high temperatures and to significantly extend the life of the engines [2, 3] by their superior temperature, oxidation, and corrosion resistance.
Typically, a TBC system consists of four layers: a top coating (TC), thermally grown oxide (TGO), a bond coating (BC), and a metal substrate [4]. The coatings insulate components from large and prolonged heat loads by using thermally insulating materials that can sustain a large temperature difference between the load-bearing substrate alloys and the coating surface. However, three-layer TBCs without a TC [5], double-layer TBCs containing only TGO and a metal substrate [6], and double-TC layer TBCs [7] have also often been investigated to find the optimum high-temperature protection for an entire TBC system. The durability of any type of TBC is significantly influenced by the mechanical performance of the TGO layer. That layer forms as a result of a chemical reaction at high temperature between metal cations and oxygen anions located at the interface between the TC and BC layers. The TGO layer helps protect the underlying substrates from high-temperature corrosion. The main reason is that the formation of TGO (mainly of α-Al2O3) at elevated temperature always develops a large stress that can cause the TGO layer to be unstable against out-of-plane displacements and even induce this layer to break down, especially upon thermal cycling. This phenomenon is described as displacement instability due to TGO formation and has attracted much attention in gas turbine engine TBCs and in HTGRs, since it is often associated with spallation, buckling, and even debonding of the TC layer from the TGO interface, all common failure mechanisms for TBC systems.
The principal cause of displacement instability is TGO formation, which results in a large compressive stress (commonly 3 to 5 GPa) in the TGO layer at dwell time. Another cause is the thermal expansion misfit among constituent materials, which often develops plastic deformation in the metal substrate in the cooling and reheating phase. Such irrecoverable deformation makes tension stress at the surface groove base, finally bringing about TBC failure from spallation or buckling. Consequently, displacement instability in TGO is of great importance for investigating TBC system thermal performance because it can significantly exacerbate the compression already occurring in the TGO layer. Extensive investigations have been performed to broaden basic understanding of the displacement instabilities of TGO. A variety of research methods, including theoretical derivation, experimental observation, and numerical schemas, have been used.
Our review of the literature found that much work has been done on modeling displacement instability of nonplanar surfaces. Also researched was thermal residual stress distribution that directly results in displacement instability or failure between interfaces occurring in TBCs from multiple thermal cycles or from thermal-mechanical load cycles. Busso et al. [8] undertook a finite element study to investigate stress development within a TBC system in which particular attention was paid to the role of variables such as elastic anisotropy within the TC, interface roughness, variation in creep strength of the BC, and the volumetric strains associated with the formation of TGO. They concluded that another important factor determining TBC stress levels was the level of elastic anisotropy of the topcoat. Han et al. [9] employed 13 finite element models to simulate the influences of different TGO thicknesses on residual stress based on unified interface morphology. In their work, TGO dynamic oxidation and the relationship between TGO thickness and thermal cycling times were taken into account, and simulation results were close to the actual situation. Yang et al. [10] developed a finite element model for a turbine blade with TBCs to investigate its failure behavior under cyclic thermal loading and found dangerous regions in the ceramic coating, as determined by the maximum principal stress criterion. Ranjbar-Far et al. [11] used a viscoelastic model in their simulation to investigate the influence of the ceramic layer and high-temperature creep characteristics of the BC layer on coating stress. They found that ceramic layer and bond layer high-temperature creep constituted a relief mechanism for coating stress. Białas [12] investigated crack development induced by residual thermal stress in TBCs using a numerical simulation scheme in which cohesive zone elements at the oxide/ceramic interface were used to model the development of interfacial microcracks. TGO layer thickening and the creep deformation of all system constituents were also modeled. They concluded that the development of the interfacial crack allows for microcrack formation within air plasma spray TBCs. Subsequent TGO growth results in a tensional zone within the oxide layer. Freborg et al. [13] developed a finite element model to evaluate the stresses induced by thermal cycling and found that the oxidation of the BC had a strong effect on ceramic layer stresses and that oxidation-induced stresses were influenced by other factors such as BC creep, TC creep, and BC roughness. Sfar et al. [14] proposed a finite element model incorporating interface cracks to simulate residual stress fields to understand the failure mechanisms. They found that the modified crack closure integral method was found to be a very efficient tool to determine the mode-dependent energy release rates for modeling TBC failure. It is worth mentioning that Karlsson and her colleagues [15–18] have performed a large amount of original work in simulating displacement instability near a surface groove on an alumina forming alloy substrate. Their research results are widely accepted by worldwide peer researchers and can be considered as classic in the field. (Their work is referred to as the classic method in this paper to distinguish it from the present work.) In their work, to explore cyclic morphological instability incurred by TBCs, they proposed a spherically symmetric fundamental model that can be solved analytically. At the same time, they developed a novel numerical method to investigate displacement instabilities for a sinusoidal undulation embedded into the BC. They also investigated distortion occurring near a single groove subjected to thermal cycling using the numerical method. They performed a series of finite element simulations revealing that TGO displacement was reduced by increasing the high-temperature strength of the underlying BC, but displacement increased with the strength of the TGO and the curvature of the groove edge.
However, their numerical simulation works well only for thinner (usually less than 1 μm) TGO layers formed at dwell times between multiple thermal cycles. For thicker TGO layers of several microns, the simulation needs to apply a larger strain at the innermost layer of TGO elements to model TGO formation. This can result in an extremely severe deformation at the periphery of a groove, leading to a divergence solution problem in ABAQUS. Ding et al. [6, 19] developed their own procedure to simulate displacement instability in the TGO layer at the interface. They were able to model TGO growth at high temperature without the limitation of TGO thickness. Their procedure worked well even when TGO thickness attained 5 μm. However, when simulating TGO crack propagation, the procedure encountered problems due to the need to incorporate many more crack parameters. Consequently, it is necessary to develop another numerical scheme to model displacement instability for greater TGO thicknesses.
In this work, another numerical simulation procedure based on the ABAQUS UMAT (User-defined Material) subroutine was proposed to simulate TGO growth at high temperatures. This procedure would be used to investigate displacement instabilities at the periphery of a surface groove embedded in a metal substrate. To determine the correctness and effectiveness of this method in the present work, we compared two methods (the classic method and the present work). We used a series of finite element analyses for simulating the final amplitude change in the surface groove, the downward displacement at the base node of the groove, the groove displacement at the periphery, and the tangential stress distribution in the elements at the base and the periphery of the groove. The comparison showed that for tangential stress in the base elements the two methods were in close agreement for all thermal cycles. For the amplitude change, the downward displacement at the base node, the groove displacement, and the stress distribution in the elements near the groove, the two methods were more in agreement for the first 3 to 6 thermal cycles. After that, inconsistency increased with the number of thermal cycles. It is interesting that the thermal cycle at which the discrepancy between the two methods occurred corresponded to a thermally grown oxide (TGO) thickness of 1 μm, which showed the accuracy of the present work over the classic method. It is concluded that the present work’s numerical simulation scheme worked better with a thinner TGO layer than the classic method and could overcome the limitation of TGO thickness by simulating any thickness.
2. Materials and Methods
2.1. Specimen Description
The specimens were rectangles of commercially available heat-resistant alloy approximately 5 mm 50 mm with a nominal composition of Fe 72.8%, Cr 22%, and Al 5%. The previous studies [5] stated that nonplanar interfaces (due to imperfections, voids, and material heterogeneity) more easily incur displacement instability than planar interfaces. To aid in-depth exploration of the deformation mechanism, two groups of grooves, one horizontal and the other vertical, were carved in the metal substrate surface to highlight displacement instability in experimental measurements. Each group consisted of two parallel grooves 20 5 μm deep and 160 μm wide. The specimens were mounted on a microcreep tester for thermal cycling, as shown in Figure 1(a and b).
2.2. Experiment Observations
The specimens were held in equilibrium between upper and lower friction grips and surrounded by a quartz tube to shield them from turbulence. They experienced 24 thermal cycles and one typical cycle consisting of 600 sec cooling from 1,200°C to the ambient, 600 sec reheating, and 1,800 sec holding at high temperature. Each specimen was cut into foils to expose the cross sections of the vertical and horizontal grooves and then polished again and examined by scanning electron microscopy. In Figure 1(c), the dotted line is the shape of the groove prior to the experiment. It is evident that the larger upward deformation occurred at the periphery of the groove, while the smaller deformation occurred at the base. At the same time, notice that the specimens also expanded in length during the experiment. For the purpose of description, the outward deformation at the groove periphery is called groove displacement, , which is defined as the distance from the point at the curved section with the larger curvature of the final TGO deformation, as seen in Figure 1(d) [20].
2.3. The Classic Method to Simulate TGO Growth
2.3.1. TGO Growth Mechanism
Since the main objective in this work is to compare two different numerical simulation strategies for modeling TGO growth at high temperature, it is necessary to first understand TGO growth mechanism in this subsection. When oxidation occurs at high temperature in TBCs, Al ions coming from the substrate react chemically with O anions from air to form alumina upon thermal cycles. It follows 4Al + 3O2 → 2Al2O3. As illustrated in Figure 2, due to the peroxidation in air, alumina has already formed first at the interface between top ceramic layer and metal substrate while it typically exists as a columnar microstructural domain (shown in Figure 2). When TBCs were dwelled at the elevated temperature, due to the thermal expansion between grain boundaries for columnar TGO, the O anions come across the gap between grain boundaries and chemically react with the Al ions resided in metal substrate. This newly generated TGO forms at the interface between preexisted TGO and metal substrate which makes TBC system thicken upon thermal cycles. At the same time, another part of O anion meets Al ion at the gaps between the internal grain boundaries to form alumina which results in an increase in TGO length. Consequently, the TGO formation at elevated temperature was simulated by the two components, namely, lateral growth strain (increasing TGO length) and thickening growth strain (increasing TGO thickness). The subsequent sections focus on how to simulate such two growth strains in numerical simulation to better reproduce the actual TGO formation at high temperature.
2.3.2. The Numerical Simulation Procedure
This section summarizes the numerical scheme to simulate TGO growth at high temperatures for TBCs experiencing multiple thermal cycles, as appeared in [6, 15–18]. We call such a numerical simulation method the classic method. Figure 3 illustrates the numerical simulation scheme described in those studies. A surface groove in the BC layer was employed to model a flaw or imperfection at the interface between the TGO and the BC. Three layers, including the TGO, the BC, and the metal substrate, were arranged in top-down order in constructing the finite element method (FEM) model. Considering the geometric symmetry, a symmetric boundary condition was imposed at the central axis of the specimen configuration. A periodic boundary condition was imposed to model the random existence of imperfections in the TBC system. To keep the same displacement in direction for the nodes located at the right side, the keyword “” was used in the ABAQUS software to impose an displacement on the horizontal axis, and the displacements at the nodes at the bottom side of the metal substrate were constrained. Figure 3(b and c) show enlarged representations of the mesh configurations near the surface grooves. The initial TGO that was formed in air was first meshed into six layers of rectangular elements, with each layer having the same thickness. Consequently, the thickness of each element layer was 1/6 that of the initial TGO thickness. For simulating TGO growth at high temperatures, two stress-free strains, a lateral growth strain and a thickening strain, were used to represent TGO formation along the lengthwise direction and along the thickness direction, respectively, with the aid of user subroutine UEXPAN in the ABAQUS software. TGO growth in the thickness direction was simulated by applying a thickening strain at the sixth layer of TGO elements. (Those were closest to the BC layer and marked as green in Figure 3(c) at each thermal cycle, not related to the other five layers of TGO elements.)
2.3.3. TGO Growth Thickness in Simulation
The final TGO thickness after thermal cycles in the classic method is as follows [6]:where is the initial TGO thickness, taken as 0.5 μm in accordance with that observed in the TBC deposition, and thus equals 1/12 μm of TGO thickness. The value for the TGO thickening strain is obtained by the following equations:
So the final TGO thickness after 24 thermal cycles follows within the limitation
After a simple calculation, (3) shows that the final TGO after 24 thermal cycles is very thin. The newly increased TGO thickness seems far less than the initial TGO thickness of 0.5 μm. As a result, the final TGO thickness that can be simulated in the classic method cannot exceed 1 μm. However, previous work [20, 21] shows that the TGO thickness observed in experiments often attains approximately 4 to 5 μm. According to the theory of classic mechanics, when simulating TGO growth with a thickness of 4 to 5 μm, that is, if TGO begins from an initial thickness of 0.5 μm to 4 or 5 μm, it must undergo a very large strain at a value of 7 to 9. Such a severe large deformation occurring in the TGO layer cannot be smoothly simulated in the classic method. Figure 4 shows a finite element simulation result for TGO deformation using the classic method, in which a TGO thickness of approximately 4 μm that forms during 24 thermal cycles is what we want to simulate in a finite element analysis (FEA). Figure 4(b) shows the deformation shape for only the first four thermal cycles. It was observed that the layer of thickening elements was sharply deformed. Note that the thickening layer is also the sixth layer of the TGO element, exactly the same as that in the classic FEM. The thickening layer must experience a very large deformation in only the first four thermal cycles, many times as thick as that of the undeformed layer. Because of the large deformation, the deformed layer of TGO elements after four thermal cycles was extracted, and because of the deformed configuration of TBCs we remeshed again on the sixth layer of deformed TGO elements for the simulation for the upcoming 5 to 24 thermal cycles, as shown in Figure 4(a). Under the same boundary conditions as before, the remeshed FEM ran again in ABAQUS for another four thermal cycles, five to eight thermal cycles, where the TGO thickening strain was again applied at only the layer of TGO elements closest to the metal substrate layer. However, when ABAQUS ran the sixth thermal cycle, it stopped due to a very severe deformation at the periphery of the groove. The accumulated large deformation finally developed a mesh contact and even mesh penetration at that site (Figure 4(a)). Furthermore, the elements in the metal substrate next to the groove periphery were also distorted drastically. The trial and error procedure in the classic FEM was unavailable to model TGO thickening during thermal cycles for large final TGO thicknesses, although it worked very well for lesser thicknesses.
2.4. Numerical Simulation Procedure Based on UMAT
2.4.1. Implementation Based on UMAT
The ABAQUS UMAT software provides a robust simulation capability to model specified material constitutive relations [22]. As an effective and flexible tool for FEA, UMAT allows the implementation of particular functions to be typically written as FORTRAN code and must be included in a model when executing an analysis. UMAT relates directly to the constitutive stress response of a material, given prescribed conditions of deformation. The user specifies a range of information to the material module relating to both the beginning and end of a time increment. In particular, stress, strain, and a deformation gradient are provided at the beginning of the time increment. Strain and the deformation gradient are also provided at the end of the increment.
Figure 5 is a flowchart for programming FORTRAN code to simulate TGO formation and material property with a high-temperature substrate metal, using UMAT. First, the TGO material properties and the metal substrate, including the elastic modulus, Poisson’s ratio, and the yielding stress, are all incorporated into a UMAT subroutine. An initial calculation for strain increment and stress increment determines the sign of the two increments. Combined with the applied load condition, if the TBC is without load, UMAT directly calculates the stiffness matrix and then evaluates magnitudes of stress and strain in order to calculate the Jacobi matrix and store data into the solution state variable. However, once the load exists (in this case the thermal cycling load is considered as the applied load), it will automatically divide the load step into multiple increments and calculate the stiffness matrix, stress, and strain within each increment according to the input of constitutive relations for TGO and the substrate material. Generally, elastic deformation first develops when the material is under initial loading. As the load increases, once the stress generated in elastic deformation exceeds the yielding stress, plastic deformation occurs. Consequently, the user should prescribe the suitable yield criterion for the material. Since a TBC consists of metallic materials, the von Mises yield criterion supports the yielding behavior during multiple thermal cycles. An incorporated flow rule in the UMAT subroutine quantitatively describes the elastic-plastic deformation of the relationship. The solver can determine whether the results in each increment converge within the allowable tolerance. If so, verification and correction for the results can be validated and then used to update the strain increment iteratively over the cycle. Finally, ABAQUS calculates the Jacobi matrix and stores the data into the state variable for subsequent reads of the thermal cycling history.
2.4.2. Verification of Implementation
Based on the above-described procedures, the user subroutine for simulating TGO growth at high temperatures only at the dwell time of each thermal cycle is implemented using FORTRAN language based on the UMAT subroutine interface. Since numerical simulation for displacement instability near a surface groove for 24 thermal cycles was a difficult and complicated solving process, it was necessary for us to verify the correctness and validity of the programmed user subroutine before it was incorporated into ABAQUS for a full analysis. In this section, we focus on the verification for a single-element and a three-element simple FEA model and for elastic and elastic-plastic material properties, respectively.
Case 1 (a single element for three different material behaviors under three load steps). This case was directed at one single CAX4 element for an implicit static elastic deformation analysis using Complete ABAQUS Environment (ABAQUS/CAE) and using UMAT based on the procedure described in Section 2.4.1. Table 1 describes the material properties and the load boundary conditions used in the FEA. For the analysis using ABAQUS/CAE, the element was assigned three different material properties under a load of 100 N in three running analyses. For the analysis using UMAT, the elements were assigned the same material properties and load condition as those for ABAQUS/CAE.
Table 2 shows a comparison of displacement and stress between the results from running ABAQUS/CAE and running UMAT. It can be seen that the displacements from the two methods are identical, while the stress values show a small difference in magnitude between them. The biggest error developed in displacement and stress does not exceed 1%. As a result, the developed UMAT for this case worked very well for the latter FEA.
Case 2 (three elements, with each assigned one different material property). As a continuation to Case 1, in this case the verified mesh contained three elements, with each one assigned a different material property but with the same applied load of 100 N. Table 3 lists the details of the material properties and load boundary conditions. Two finite element analyses were performed using ABAQUS/CAE and UMAT. Table 4 shows a comparison for the calculated displacements and stress from the two calculations. It can be seen that the displacement and stress between the two calculations show larger errors than do those in Case 1. The displacement value had an error of 8% and the stress value 3%. This is mainly attributable to the element discretization, since three elements representing three different materials resulted in more serious nonlinearity than what occurred in Case 1.
Case 3 (three elements with a combination of various materials for elastic-plastic analysis). The previous two cases aimed mainly at a linear analysis. Case 3 involves material nonlinearity including yielding stress for verifying the subsequent complicated analysis. The FEM model for this case continued to use three elements to discretize the model. Table 6 describes three cases of material properties with a combination of three different material parameters. Each detail is listed in Table 5. Each case involved a different load and could be simulated by using ABAQUS/CAE and UMAT. That is, for number 1 analysis three elements were assigned the same material property with a load of 100 N. For number 2 analysis, elements 1 and 3 were assigned the same material property while element 2 was assigned material different from that of the other elements, with a load of 200 N. For number 3 analysis, elements 1 and 2 were assigned one material model, but element 3 had a different material model with a load of 300 N.
Table 7 shows a comparison of the displacement and stress values between the ABAQUS/CAE and UMAT running ways. The comparison shows that the error in displacement between the two running methods dramatically decreased as the material definition changed and the load increased, from 13.9% down to 2.7%, while the stress value indicates a similar trend in error from 4.9% to 0.54%. Also from the contour plot for von Mises stress, plastic deformation did not occur in number 1 analysis mainly due to the smaller applied load, while the von Mises stress for number 3 exceeded the yielding stress of 200 MPa. The error decreased as material yielding occurred.
In summary, through the comparison of calculated displacement and stress values between two different running ways, ABAQUS/CAE and UMAT, it can be seen that the obtained results from UMAT approached those from ABAQUS/CAE, since the CAE running way has commonly been considered the correct way. Although there was a small difference between them in Case 2, such a discrepancy can be reduced substantially by using more elements to discretize the models. In addition, accuracy improved substantially with the introduction of plastic deformation in Case 3. We concluded that using UMAT to simulate TGO growth at high temperatures could enable FEM results within a very acceptable tolerance limit.
3. Results and Discussions
3.1. The Change in Amplitude of Surface Groove
Based on the previously described work, the major difference between the two methods (the classic method and the present work) is that the numerical simulation scheme was used to model TGO growth at an elevated temperature through multiple thermal cycles. It is widely accepted that the numerical simulation scheme in the classic method can work very well for smaller TGO thicknesses, usually less than 1 μm, while we found that the scheme in the present work can work well also for larger thicknesses of TGO. In this section, by comparing the change in surface groove amplitude resulting from the two different methods we found that for a thinner TGO the methods coincided. However, for a thicker TGO the present work continued to calculate but the classic method stopped calculating due to some error. As for the amplitude of a surface groove, it can be defined in terms of the distance between the peripheries of the surface groove to the base of that, as illustrated in Figure 6.
Two finite element analyses, one based on the classic method and the other on the present work, were conducted to calculate the change in amplitude of a surface groove, where the FEMs were the same except for the numerical simulation scheme simulating TGO growth at high temperatures. Both schemes simulated morphological deformation for the vertical groove corresponding to the experimental observation, as illustrated in Figure 1(a). The FEM model consisted of 26,936 quadrilateral, first-order generalized plane strain elements, CPEG4 in ABAQUS, with the number of the elements determined after a mesh dependency check. Due to geometrical symmetry, half models were taken into account. The symmetric boundary condition was imposed on the side at . The periodic boundary condition was applied for those nodes located on the right side such that all those nodes had the same displacements in the direction by allowing the same movement in the direction. To avoid a rigid displacement of the overall system in the direction, any node in the model was chosen to fix the displacement in the direction. The material properties for this numerical simulation were taken similarly to our previous studies [6, 19, 20]. For the convenience of comparison, a TGO growth thickness of 1.2 μm formed at high temperature was assumed for the two methods. Thus, for the classic method, the initial TGO thickness when constructing the FEM model was taken to be 1 μm. As required, the TGO could be divided into six layers. The thickening strain, , was applied to only the sixth layer of TGO element (closest to the metal substrate elements). It was taken as a value of 0.05 in each thermal cycle in line with (1). For the present work, it was hypothesized that 1.2 μm of TGO needs to be uniformly formed within 24 thermal cycles. As a result, the thickness of TGO formed in each thermal cycle equaled 0.05 μm. That is to say, when constructing the FEM model of the present work, 24 layers of elements with an equal thickness of 0.05 μm should first be meshed, and the material property for those elements at the moment should be taken the same as the metal substrate. In the first thermal cycle, the outermost layer should be transferred in material property from the metal substrate into the TGO with the aid of the ABAQUS subroutine UMAT. The procedure was repeated 24 times, finally resulting in 1.2 μm thickness of TGO formed.
Figure 7 shows the variation in amplitude change of the surface groove between the two methods with the number of thermal cycles. It can be seen that when thermal cycles 1 to 9 were run, a higher coincidence in amplitude change between the two methods was attained, showing that the accuracy of the simulation scheme in the present work was comparable to that of the classic method. The ninth thermal cycle corresponds to a formed TGO thickness of 0.45 μm, further indicating that for a thinner TGO the two methods can work very well and that the accuracy of the present work is comparable with that of the classic one. However, after the first nine thermal cycles, a difference in amplitude change between the two methods appeared slight at first and then increased significantly as the thermal cycles continued. The value of amplitude change in the present work seemed smaller than that of the classic method. This was possibly because the stress-free strain used to model TGO formation exerted a slightly greater pressure on the metal substrate to relax the higher TGO growth stress at high temperature. This may have led to a more severely plastic deformation, which was considered a requirement for displacement instability around a surface groove in the BC of TBCs. As the number of thermal cycles increased, the thickness of the newly generated TGO became greater. At the 20th thermal cycle, the FEA using the classic method stopped running in ABAQUS due to a divergence problem, while the analysis using the present method continued running until the end of the 24 thermal cycles. The 20th thermal cycle coincided with 1.0 μm of TGO thickness, which further supported the deduction, as described previously, that TGO thickness should not exceed 1.0 μm for the classic method.
3.2. Tangential Stress in the Groove Base Element
The tangential stress in the groove base element contributed greatly to the displacement instability around the surface groove, leading to a morphological change of the BC surface [20] (the metal substrate surface in this work, due to the double-layer TBC system used). Consequently, it is important to study the tangential stress to better understand the displacement instability. Figure 6 shows the element of interest in this section, labeled number 1, which was at the base of the surface groove. Both the downward displacement at the base of the groove and the substantial upward displacement at the periphery of the groove were observed in the experiment and in the numerical simulation, which was attributed to tangential stress (called “in-plane” stress in some works [15–18]). That stress increased significantly at the dwell time of a thermal cycle where TGO grew at high temperature. It is worth mentioning that the classic method cannot work for TGOs thicker than 1.0 μm, which is verified in Figure 7. However, for comparison between the two methods, in this section the final TGO thickness upon 24 thermal cycles in all of the FEAs was just 1.0 μm. That assured that the two methods would work well for 24 thermal cycles and enabled a valid comparison between the two methods.
Figure 8 shows the tangential stress variation in number 1 element at the base of the surface groove, with thermal cycles. Most interesting was that when the final TGO thickness was taken as 1.0 μm, both FEAs worked well until the end of the 24 thermal cycles. The overall FEA results showed greater agreement between the two simulation methods, one based on a procedure using stress-free strain to model TGO formation and the other using a subroutine to calculate TGO formation cycle by cycle, for tangential stress in element number 1. This implies that the two methods were correct and maintained a higher consistency in the simulation scheme. To investigate tangential stress variation in detail, Figure 8(b) shows the amplified plot for the previous second thermal cycle (the dashed circle in Figure 8(a)). Even within each step in the first and second cycles, the tangential stresses from the two methods were still in close agreement.
(a)
(b)
Note in Figure 8(b) that the tangential stress first decreased during the thermal cycles because the temperature also decreased from a high temperature to the ambient temperature obtained during a compression of up to 3.0 Gpa. Such a great compressive stress forced the TGO material to move downward and the metal substrate to yield. When the high temperature returned, the compression in the TGO layer relaxed significantly to an almost zero stress state. But a slight tension stress at the end of the reheating period indicated that plastic deformation developed in the metal substrate. TGO grew within the third step. The tangential stress dropped to zero first and then decreased slightly to a steady compressive stress due to TGO plasticity at the high temperature (both of them were realized by using the temperature dependent material properties, and it can be referred to in [5, 15]).
3.3. Downward Displacement at the Base of the Groove
Except for tangential stress in elements at the base of the surface groove, the downward displacement for the node at the groove base was also investigated for comparison between the two simulation methods. The location of node 1 is illustrated in Figure 6. The main reason for choosing node 1 was that to investigate amplitude change in a surface groove the downward displacement of node 1 must be determined. Figure 9 shows the variation of downward displacement at node 1 with the number of thermal cycles. In comparison with the higher consistency between the two methods illustrated in Figures 8(a) and 8(b), there was a significant difference between the two methods from the third thermal cycle onward. Particularly from the 15th thermal cycle, they even showed a slight contrary trend in displacement between them. In the classic method, the displacement trended upward, not downward, even reaching a positive displacement. However, for the present work the displacement at node 1 was more downward as thermal cycles continued, showing an opposite displacement direction.
At the first few thermal cycles, the initial boundary condition and the applied temperature were same for both FEAs, although they used different simulation methods to simulate TGO formation at a high temperature. As also predicted by Figures 7, 8, and 9, for a thinner TGO, especially for a much thinner TGO corresponding to the first few thermal cycles, the simulation scheme was in essence similar between them. As a consequence, they showed higher consistency at the beginning of the thermal cycles. However, as thermal cycles continued, the newly generated TGO accumulated cycle by cycle and would influence mechanical performance in TBCs. In the classic method, TGO formation is simulated mainly by thickening and lengthening strains, with each thermal cycle being equal to a constant value. In that case the thickening strain was taken as 0.05 in each thermal cycle, while in the present method the newly increased TGO thickness was 0.042 μm at each thermal cycle. Furthermore, for the classic method the thickening strain was applied only at some layer of elements, usually with 1/6 μm of thickness. Obviously, the newly generated TGO thickness in each thermal cycle differed between the two methods although the final thickness was the same. TGO growth at each cycle in the present work was greater than that of the classic method. Consequently, more downward deformation was possibly produced in the present work than in the classic method due to the relatively thicker TGO thickness at the first few thermal cycles. Due to the cycle-by-cycle accumulation, a discrepancy was introduced in the downward displacement between the two methods, as shown in Figure 9.
3.4. Groove Displacement at the Periphery of the Surface Groove
The previous experiment observations illustrate that the specimen made of Fecralloy showed a significant change in groove profile after 24 thermal cycles as “pile-up” of the alloy at the periphery of the surface groove developed [5]. Our previous research [20] showed a similar deformation profile at the groove periphery, although the displacement at the base of the groove seemed smaller. Therefore, deformation occurring at the groove periphery is very important in studying depth displacement instability. To illustrate the extent of deformation around the surface groove, Figure 1(d) explicitly represents it by using groove displacement, which is defined as the distance measured from the displaced location of the most outward node (node 50 marked in Figure 6) to the original position of that node [20]. Accordingly, the elements surrounding that node and consisting of the groove periphery were considered of interest.
Figure 10 shows the final resultant displacement, equal to the groove displacement, against the number of thermal cycles. Figure 10(b) is the amplification plot for the area of Figure 10(a), surrounded by a dotted line, corresponding to the 7th to 12th thermal cycles. As shown in Figure 10, analogous to the previous analysis, in the first eight thermal cycles the groove displacement for node 50 for two methods shows good agreement, much earlier than Figure 9, of three thermal cycles. However, after eight thermal cycles, a greater and greater difference between them appears in the plot. Unlike what is mentioned in Section 3.3, the discrepancy was still small even at the end of the 24 thermal cycles. The detailed plot, as illustrated in Figure 10(b), shows that during 7 to 12 thermal cycles there remained a small difference between them. The result from the classic method seems slightly larger than that from the previous analysis.
(a)
(b)
Similarly, for the first few thermal cycles, due to the same boundary condition applied and the relatively thinner TGO, both FEM methods showed a minor difference in groove displacement. As described in Section 3.3, for the present work a thicker TGO than that in the previous method was employed to conduct an FEA. The occurrence of displacement instability around the surface groove could be attributed to plastic deformation in the metal substrate and the TGO layer following the conservation of overall volume. In Section 3.3, we stated that for the present work a greater downward displacement happened at the base of the surface groove. The TGO relaxed its higher compression by displacing the node at the base in a downward movement. Consequently, the elastic strain energy stored in the TGO for the present work seemed less than that for the classic method. As a result, in the area around the periphery of the groove there was less potential for the present work than the classic method to release the stored strain energy by displacing outward at the periphery, manifested as a groove displacement. That is the reason that the groove displacement for the present work was less than that for the classic method.
3.5. Stress Distribution at the Periphery of the Groove
The driving force to enhance the outward deformation around the periphery of the surface groove cycle by cycle was the stress generated in the TGO that originated from both a thermal expansion misfit between different constituent materials and TGO growth stress during TBC dwell time at the high temperature. The tangential stress generated in TGO elements around the surface groove contributed the most to the outward deformation at this site. Figure 11 shows the contour plot for tangential stress at TGO elements located at the periphery of the surface groove. Actually, tangential stress generated in these TGO elements was aligned with the circumferential direction of the groove, which dominated expansion of the circular section (as in Figure 11), further resulting in the occurrence of displacement instability of the surface groove.
(a) The tangential stress distribution in TGO layer from FEA using the classical method
(b) The tangential stress distribution in TGO layer from FEA in the present work
From the comparison, we can see that the tangential stress distribution for the present work appeared more uniform than for the classic method. It can be seen from Figure 11(a) that the greatest tangential stress was distributed primarily along the innermost layer of elements, which is just applied in the thickening strain at each thermal cycle. While at the middle part of the curved section, the biggest stress value (in tension) was at the outside of the TGO growth layer, corresponding to the projection of node 50 on the thickening layer of TGO elements. In the classic method, TGO growth was realized by applying a stress-free thickening strain at only the sixth layer of elements (the red layer elements in Figure 11(a)). When ABAQUS started running, the thickening strain was suddenly applied as an insert, which possibly caused the stress fluctuation and a nonuniform distribution around the periphery of the surface groove. However, as for the specific stress values, they were almost the same: 71.78 Mpa for the classic method and 68.69 Mpa for the present work, both in compression. The tolerance in stress value lay within acceptable limit. For the present work, as depicted in Figure 11(b), the stress values over the overall TGO elements showed a smooth transition, almost without a stress gradient, in that TGO growth was modeled by linearly transforming metal substrate material to TGO material property over the cycles. From the numerical simulation itself, the result for tangential stress for the present work seemed better than that from the present method.
Figure 12 shows the stress variation for element 690, located at the inner side of the curve section, with the number of thermal cycles. The stress state for the element upon 24 thermal cycles occurred in compression, with almost no stress in tension, satisfying the demand for the occurrence of displacement instability at the periphery of the surface groove. Figure 12 also shows a higher consistency in tangential stress for element 690 between the two methods in the first four thermal cycles. The possible reason for that should be the same as that previously described. However, after the fourth thermal cycle the comparison between them showed a difference in the maximum compressive stress value, and the discrepancy increased with the increasing number of thermal cycles. It is interesting that for a total of 24 thermal cycles the stress values in TGO growth (sometimes slightly in tension) were similar for the two methods. In one thermal cycle, the tangential stress always experienced three phases: compression, relaxation of compression, and constant values. These periods synchronized with the three steps of one thermal cycle. Although the two methods used different methods to simulate TGO growth at high temperatures in FEA, the consistency in stress state at TGO growth indicates that the methods were consistent, further implying the correctness of the present work. The greater and greater discrepancy in the compressive stress state, however, can be attributed primarily to the volume change induced by the newly generated TGO in each thermal cycle.
Figure 11 shows the volume change between them. The element at the sixth layer of the TGO element (counted from the outward to the inward direction) for the classic method is much thicker than that for the present work. This is because the sixth layer of TGO elements in the classic method expressed the overall TGO growth at high temperature, while for the present work the newly formed TGO upon 24 thermal cycles is not displaced in Figure 11. Actually, the final deformed TGO profile in the present work seems much thicker than that for the classic method. Notice that the final TGO thickness resulted not only from TGO growth at high temperature but also from other strains including elastic strain, thermal expansion, and creep strains. Even in a single FEA, final TGO thickness became greater and less uniform as the lengthening strain increased [5]. Accordingly, the different final TGO thicknesses contributed to the discrepancy in compression stress value between the two methods. For the present work, the variation in thickness of TGO formed between any two thermal cycles was greater than that for the classic method. The large volume change induced by TGO formation can affect a greater compressive stress state.
4. Conclusions
A comparison in this work between two different methods, the classic method and the present work, was made for the final amplitude change in a surface groove, the downward displacement at the base node of the groove, the groove displacement at the periphery, and the tangential stress distribution in the elements located at the base and at the periphery of the groove. The following conclusions can be drawn:(1)A comparison of the final amplitude change of a surface groove derived from the two methods showed greater consistency in the first nine thermal cycles. The discrepancy between the methods increased as the number of thermal cycles increased. The classic method stopped running at the 20th thermal cycle, at which TGO thickness was close to 1 μm, which indicated the correctness of the present work.(2)In a comparison made of the downward displacement at the base of the groove, the tangential stress in elements at the base and the periphery, and the groove displacement at the periphery, both methods had similar results in the first few thermal cycles. This demonstrated that the methods were similar in essence, although they each had a different numerical simulation scheme to model TGO growth at high temperature.(3)A comparison was made of different mechanical behaviors that have some effect on the displacement instability around a surface groove. The result showed that the numerical simulation scheme in the present work worked very well for a thinner TGO in comparison with the classic method but also overcame the limitation of TGO thickness and was able to simulate any TGO thickness.
Competing Interests
The authors declare no competing interests.
Acknowledgments
This work is financially supported by the Natural Science Foundation of China (11302272) and by NSAF (Grant no. U1530140).