Modeling of Creep Deformation Behavior of DZ411 and Finite Element Simulation of Turbine Blade

Abstract

Creep tests were conducted on DZ411 material at 930 °C and 850 °C, and creep curves were recorded and employed in normalization creep model building. A yield function suitable for directional solidification nickel-based materials was proposed in an ascending-order approach. Combined with the normalized creep model and the proposed function, a creep subroutine was compiled to simulate the creep deformation behavior of a turbine blade. The typical boundary conditions of the blade were determined and used for finite element analysis. According to the analysis results, the assessment positions for the actual application of a turbine blade were determined and checked for endurance intensity. The phenomenon of deviation angle between crystal axis and blade height direction in actual casting was further analyzed. Multiple angles and directional deviation angles were simulated for 10,000 h creep deformation. Considering the difficulties and challenges of the complex geometric structures of blades, it is necessary to conduct creep tests of DZ411 material and a simulation analysis of a real blade. Based on the above analysis and discussion, the present work sheds light on finite element analysis and has great potential for structural analyses in the engineering applications of complex high-temperature structures.

1. Introduction

With the continuous increase in the turbine inlet temperature of aero-engines and gas turbines, the creep problem is becoming increasingly prominent in high-temperature structural components [1]. In long-term high-temperature and high-load environments, creep deformation inevitably occurs in turbine blades [2,3]. There is a complex internal structure in rotor blades, including elements such as the film-cooling hole, exhauster window, inner cavity, and transition fillet. The structural design should not only satisfy the requirements of weight reduction and cooling, but also meet security requirements. Therefore, the need for the creep behavior simulation of turbine blades in service environments is necessary, together with the analysis of stress relaxation phenomena.

In the second half of the 20th century, the performance of high-temperature alloys significantly improved, propelled by rapid advances in both alloy materials and casting processes. In 1953, Rutter and Chalmers [4] put forward the theory of composition supercooling. Later that same year, the conditions for constitutional supercooling were established on a quantitative basis by Tiller, Rutter, Jackson, and Chalmers [5]. The foundation of directional solidification technology was laid by the research results. Through directional solidification technology, the solidification process of castings could be controlled to obtain a parallel columnar crystal structure and eliminate the transverse grain boundary perpendicular to the stress axis [6]. Since the 1980s, directional solidification high-temperature materials have rapidly developed along their unique path. Compared with cast alloy K403 with the same composition, directional solidification alloy DZ3 has 2–3 times higher tensile plastic properties, 60–130 MPa higher endurance strength, and nearly 100 times higher thermal fatigue life [7]. Due to the application and optimization of directional solidification nickel-based alloys, the creep property of turbo blades has been greatly improved.

In practical engineering applications, the creep behavior simulation of nickel-based high-temperature alloys is mainly analyzed using phenomenological macroscopic models [8]. Macroscopic models are realized using anisotropic tensors, describing the inelastic anisotropic creep deformation behavior. This research on creep models is mainly divided into three sections. (1) On the basis of creep activation energy, a creep model is established with creep strain rate as the main measurement index [9,10,11]. A model in this category has the advantage of clear physical meaning and the disadvantage of low accuracy and lack of efficiency. (2) With the essential variables of temperature, stress, and time, a creep model is put forward to accurately describe creep strain curve with time [12,13,14,15]. The advantage is that the model establishment and parameter acquisition are relatively simple, with high accuracy and feasibility in engineering applications. The disadvantage is that the physical meaning is unclear. (3) Applying the principle of damage mechanics and introducing the evolution law of damage variables, the creep damage model is established [16,17,18]. This type of model has the advantage of definite physical significance but suffers from the disadvantage of cumbersome expression.

Many scholars dabble in this issue, and quite a number of works propose various models to describe creep behavior. At the same time, the finite element method is used to numerically simulate the creep deformation of structural components such as turbine blades. Despite the significant amount of research being conducted, there are still some shortcomings in engineering applications. The commonly used creep models cannot fully describe the creep deformation of the three stages, especially the tertiary creep stage. The general finite element software ANSYS 19.0 (Canonsburg, PA, United States) provides 13 types of creep models [19], all of which can only describe the first and second stages of creep deformation. The method used to calculate the three stages of creep, via the introduction of damage parameters into partial classic creep models or viscoplastic constitutive models, is relatively complex in calculation and difficult to apply in actual structural analysis. Therefore, there is an urgent need to research creep models and numerical simulation in order to meet design needs in the creep analysis of actual engineering structures.

Aiming to understand the deformation behavior of nickel-based superalloy materials, many domestic and foreign scholars have researched its constitutive equation, with the yield criterion as one of the core contents. Most of the existing models have the disadvantages of complex form, difficult parameter acquisition, low calculation accuracy, and poor efficiency. Viscoplastic constitutive theory is usually divided into two types [20]: (1) associating the inelastic flow rule with a certain yield condition [21,22], and (2) introducing internal variables into the inelastic strain rate tensor equation directly to represent the resistance of the material to inelastic flow [23,24]. Cazacu et al. [25] proposed a yield criterion containing six undetermined coefficients on the basis of the generalized second and third deviator stress invariant. Monotonic stretching data from at least six orientations were required to solve a series of nonlinear equations. Darrieulat and Piot [26] believed that the slip system was activated above the threshold shear stress and proposed the corresponding yield function. As the parameter n tended to infinity, the associative yield surface approximated the Schmid hyperplane infinitely so as to avoid the helplessness in derivation. Ramaglia [27] proposed a yield function for the interaction between the octahedral main slip system and the cubic slip system. The parameters in the yield function were related to the material and temperature characteristics, but the part of cubic slip system was taken as zero at lower temperatures.

It should be pointed out that virtually all major hollow turbine blades have complex cavity structures, resulting in the mesh number reaching into the hundred-thousands and even millions. Simultaneously, the stress distribution of hollow turbine rotor blades is complex, and the temperature range is wide. For the simulation of the elastic–plastic deformation and viscoplastic deformation of nickel-based superalloys’ turbine rotor blades, the constitutive model should have high computational efficiency and robustness, but it is often difficult for the models in the existing literature to meet the engineering requirements.

In this paper, the creep tests of DZ411 at 850 °C and 930 °C had been designed and carried out under multiple stress levels. With the creep curve fitting to the experiment data, normalization model parameters were obtained. In view of the shortages of the existing yield function, this article proposed a reasonable and lower-order yield function form by increasing the order of the yield function polynomial gradually, and compiled a corresponding creep subroutine. Combined with the model of a practical engineering turbine blade built from nickel-based superalloy DZ411, a creep deformation simulation for 10,000 h was accomplished, which met a certain computational accuracy and had high engineering efficiency. The stress relaxation phenomenon throughout the entire creep process was discussed, with further analysis of the simulation results of creep deformation behavior. Meanwhile, with no guarantee of the absolute accuracy of crystal orientation in the directional solidification process, the influence of deviation angle on the creep resistance of nickel-based superalloy DZ411 was analyzed to simulate the corresponding situation in practical engineering applications more accurately.

2. Materials and Methods

2.1. Composition and Specimen

The nominal compositions of DZ411 directional solidification superalloy [28] are shown in

Table 1. Nominal compositions of DZ411, data from [28].

Element	C	Cr	Co	Ni	W	Mo	Al	Ti	Ta
Composition/wt.%	0.07~0.12	13.5~14.3	9.0~10.0	Bal.	3.5~4.1	1.3~1.7	2.8~3.4	4.6~5.2	2.5~3.1
Element	B	Si	P	S	Pb	Bi	As	S	Sb
Composition/wt.%	0.007~0.02	≤0.2	0.005	0.01	0.0005	0.0001	0.005	0.002	0.001
Comment: $\omega \left(\mathrm{Al}+\mathrm{Ti}\right)\ge 7.80$

. A vacuum melting casting process was adopted. The standard heat treatment process of casting is: solution treatment (1225 ± 10 °C, 2 h, air cooling), primary aging (1120 ± 10 °C, 2 h, air cooling), and secondary aging (850 ± 10 °C, 24 h, air cooling).

In order to complete the creep tests of directional solidification superalloy DZ411 under different temperature and stress combinations and obtain the creep strain versus time curve, a round bar specimen was designed, as shown in Figure 1. The gauge length of the specimen used for measuring creep strain was 25 mm, with a diameter of 5 mm. The processed specimen is shown in Figure 2.

2.2. Method and Equipment

The crystallization direction of directional solidification alloy DZ411 was selected to carry out the experiment. First, the temperature of the specimen was raised to 850 °C or 930 °C. Then, the temperature was maintained for one hour. After temperature homogenization, the load was increased to target and a constant tensile load was maintained until rupture. The constant tensile load was a constant tensile force without considering the changes in the cross-sectional area of the specimens. Finally, the endurance life was recorded, and the creep curve was drawn.

The QBR-100 microcomputer control electron creep testing machine (Changchun Qianbang, Changchun, China) was used for the test, as shown in Figure 3. Constant force control mode was used for creep testing. A three-zone fired furnace was adopted for creep testing, with upper, middle, and lower sections equipped with an S-type thermocouple for temperature measurement. The temperature in the high-temperature furnace was evenly distributed, ensuring the uniform heating of the specimen. During the experiment, a high-temperature extensometer was used to measure the strain in the gauge range to obtain the complete creep curve.

2.3. Endurance Life and Fracture Elongation

The endurance life and fracture elongation obtained from the experiments under different temperature and stress conditions are shown in

Table 2. Endurance life and fracture elongation of DZ411 specimens.

Serial Number	Temperature /°C	Load /MPa	Life /h	Elongation /%	Serial Number	Temperature /°C	Load /MPa	Life /h	Elongation /%
A-1	930	360	18.1	19.21	B-1	850	540	17.1	14.56
A-2	930	340	21.6	14.77	B-2	850	510	36.9	15.17
A-3	930	320	38.0	10.56	B-3	850	480	67.4	10.72
A-4	930	320	28.1	9.85	B-4	850	450	128.0	13.41
A-5	930	290	57.3	10.26	B-5	850	450	226.3	11.83

. According to the test, the average fracture elongation of DZ411 directional solidification superalloy is 12.9%.

2.4. Creep Curve

The creep strain curves of the specimens against time are drawn in Figure 4, classified according to temperature. The experiments were carried out under constant tension, without considering changes in the cross-sectional area of the specimens. The obvious three stages of creep are presented in all specimens, namely the periods of attenuation, steady increase, and acceleration.

2.5. Microstructures after Testing

Optical micrographs (Figure 5) taken from the fracture region of creep rupture specimens to present the effect of accumulated strain on the microstructure.

Figure 5a shows the fracture feature of a dimple under 540 MPa tensile stress at 850 °C, and Figure 5b shows intergranular fracture under 290 MPa tensile stress at 930 °C. Intergranular fracture characteristics are observed in all specimens at 930 °C. Dimples, a typical feature of high-temperature metal ductile fracture, were observed on the fracture surface of all specimens at 850 °C. The residual stress after mechanical processing of directionally solidified alloy caused surface recrystallization of the sample under high temperatures and long-term creep, resulting in intergranular fracture morphology of the fracture surface. All tests were conducted in atmospheric furnaces, which may be a reason for the surface recrystallization. Reference [29] suggested that surface recrystallization of materials led to a decrease in endurance life.

3. Theoretical Model of Directional Solidification Superalloy

3.1. Yield Criterion of Directional Solidification Superalloy

Aiming at the directional solidification material, the following three requirements need to be met:

(1) Each grain of the directional solidification material is an orthotropic cubic symmetric material, and the material properties of three mutually orthogonal crystal axes are identical, that is, the yield functions possess symmetry in x, y and z directions.

(2) The proposed yield function has the capacity to degrade into the Mises yield function applicable to isotropic materials.

The yield functions corresponding to the Mises criterion and the Hill’s criterion are

$$f=\frac{1}{2}[{({\sigma }_{yy}-{\sigma }_{zz})}^2+{({\sigma }_{zz}-{\sigma }_{xx})}^2+{({\sigma }_{xx}-{\sigma }_{yy})}^2]+3({\sigma }_{yz}{}^2+{\sigma }_{zx}{}^2+{\sigma }_{xy}{}^2)$$

(1)

$$f=F{({\sigma }_{yy}-{\sigma }_{zz})}^2+G{({\sigma }_{zz}-{\sigma }_{xx})}^2+H{({\sigma }_{xx}-{\sigma }_{yy})}^2+2L{\sigma }_{yz}{}^2+2M{\sigma }_{zx}{}^2+2N{\sigma }_{xy}{}^2$$

(2)

The quadratic polynomial form was generally adopted by classical yield functions, thus assuming that the yield function for the directional solidification superalloy was in polynomial form. For the convenience of expression, the polynomial was expressed in the form of deviatoric stress. Before simplification, considering that there were 5 independent terms of the first-order polynomial and 25 independent terms of the second-order polynomial, the number of terms of the third-order polynomial would increase exponentially. In order to balance computational efficiency, only the first- and second-order polynomials were considered.

The coordinate system composed of three principal axis directions was defined as the material coordinate system oxyz. Each grain of the directional solidification superalloy had the same material property in the three crystal axis directions. This resulted in the symmetry of the x, y and z directions. The operation of combining like terms was performed on the parameters with symmetry term. The preliminary yield function form of directional solidification superalloy was expressed as

$$\begin{array}{ll}f& =A(s_{xx}+s_{yy}+s_{zz})+B(s_{xy}+s_{yz}+s_{xz})+C(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)+D(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)\\ & +E(s_{xx}{s}_{yy}+s_{yy}{s}_{zz}+s_{xx}{s}_{zz})+F(s_{xy}{s}_{yz}+s_{xy}{s}_{xz}+s_{xz}{s}_{yz})\\ & +G[s_{xx}(s_{xy}+s_{yz}+s_{xz})+s_{yy}(s_{xy}+s_{yz}+s_{xz})+s_{zz}(s_{xy}+s_{yz}+s_{xz})]\end{array}$$

(3)

where s_xx, s_yy, s_zz, s_xy, s_yz, and s_zx were all deviatoric stress components, satisfying $s_{xx}+s_{yy}+s_{zz}=0$. And according to the formula derivation, the components met the requirement $\frac{1}{2}(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)=-\left(s_{xx}{s}_{yy}+s_{yy}{s}_{zz}+s_{xx}{s}_{zz}\right)$. The directional solidification superalloy was considered to have a lack of tension-compression asymmetry, leading to $B=0$. The yield function was given normalized treatment by taking parameter $C=1$. The yield function was further simplified to

$$f=(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)+D(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)+F(s_{xy}{s}_{yz}+s_{xy}{s}_{xz}+s_{xz}{s}_{yz})$$

(4)

When parameter $D=2$ and $F=0$, the above equation corresponded to the Mises yield criteria. When parameter $F=0$, the above equation corresponded to the Hill’s yield criteria, taking into account the symmetry of x, y, and z. The yield strengths corresponding to uniaxial loading in the $[111]$ and $[\overline{1}11]$ direction of the grain were identical. If the $[111]$ and $[\overline{1}11]$ direction was selected as the loading axis, ${\sigma }_{[111]}={\sigma }_{[\overline{1}11]}$, the deviator stress tensors in the crystal axis coordinate system were written as

$$\left[\begin{array}{ccc}0& \frac{\sigma }{3}& \frac{\sigma }{3}\\ \frac{\sigma }{3}& 0& \frac{\sigma }{3}\\ \frac{\sigma }{3}& \frac{\sigma }{3}& 0\end{array}\right] \mathrm{and} \left[\begin{array}{ccc}0& -\frac{\sigma }{3}& -\frac{\sigma }{3}\\ -\frac{\sigma }{3}& 0& \frac{\sigma }{3}\\ -\frac{\sigma }{3}& \frac{\sigma }{3}& 0\end{array}\right]$$

The above equation was substituted into the yield function, and the parameter $F=0$ was obtained. Therefore, the yield function only contained parameter D, corresponding to Hill’s yield criterion. At this time, the proposed yield function had insufficient ability to describe the yield strength of the directional solidification superalloy in different orientations. By increasing the order of the yield function, the yield strength of the directional solidification superalloy could be described more accurately in different orientations and under different load conditions.

The cubic polynomial was obtained by multiplying the simplified quadratic polynomial and a linear polynomial, as shown in Equation (5)

$$\begin{array}{ll}f& =\left[C(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)+D(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)+F(s_{xy}{s}_{yz}+s_{xy}{s}_{xz}+s_{xz}{s}_{yz})\right]\\ & \times \left[A(s_{xx}+s_{yy}+s_{zz})+B(s_{xy}+s_{yz}+s_{xz})\right]\end{array}$$

(5)

From the above derivation process, it was easy to draw the conclusion that the coefficient of the cubic term of any single term is 0. However, the combination of specific cubic terms satisfied the characteristics of the yield function of nickel-based directional solidification material, such as the third invariant form of the deviator stress tensor

$$J_3=s_{xx}{s}_{yy}{s}_{zz}+2s_{xy}{s}_{xz}{s}_{yz}-\left(s_{xx}{s}_{yz}{}^2+s_{yy}{s}_{xz}{}^2+s_{zz}{s}_{xy}{}^2\right)$$

(6)

Considering that the proposed yield function should be able to degenerate into the Mises criterion form, the cubic polynomial form was not considered for the time being. Furthermore, the forms of quartic and quadratic polynomials were selected. Using the above conclusion and for the convenience of being reduced to the Mises criterion form, the function was preliminarily simplified as follows

$$\begin{array}{ll}f& =H{(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)}^2+I{(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)}^2+J(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)\\ & +K(s_{xy}{}^2{s}_{yz}{}^2+s_{xy}{}^2{s}_{xz}{}^2+s_{xz}{}^2{s}_{yz}{}^2)+C(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)+D(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)\end{array}$$

(7)

To ensure that the proposed yield function had the ability to degrade into isotropic yield function, we set parameter $K=4$. To ensure that the yield function was dimensionless, we set parameters $C=D=0$. The function was normalized by defining parameter $H=1$. The yield function with the capability of characterizing different orientations was obtained as

$$f={\left(s_{xx}^2+s_{yy}^2+s_{zz}^2\right)}^2+I{\left(s_{xy}^2+s_{yz}^2+s_{zx}^2\right)}^2+J\left(s_{xx}^2+s_{yy}^2+s_{zz}^2\right)\left(s_{xy}^2+s_{yz}^2+s_{zx}^2\right)$$

(8)

where parameter $I=J=4$. The proposed yield function was capable of degeneration into the Mises yield criterion form, thereby transforming from a nickel-based superalloy yield function to a yield function suitable for isotropic materials. From the above formula, it was evident that the proposed yield function had symmetry in x, y, and z directions, meaning that the performance of the material was identical in the three crystal axis directions. Due to the form of deviatoric stress, the proposed yield function was independent of hydrostatic stress. However, the asymmetry of tension and compression was not considered in the proposed yield function. The yield function was applicable to anisotropic materials, such as DD6 [30]. The defined yield function was rewritten as a stress tensor expression, as shown in Equation (9)

$$\begin{array}{ll}f& =\frac{1}{9}{[{({\sigma }_{yy}-{\sigma }_{zz})}^2+({\sigma }_{zz}-{\sigma }_{xx}{)}^2+({\sigma }_{xx}-{\sigma }_{yy}{)}^2]}^2+I{({\sigma }_{yz}{}^2+{\sigma }_{zx}{}^2+{\sigma }_{xy}{}^2)}^2\\ & +\frac{J}{3}[{({\sigma }_{yy}-{\sigma }_{zz})}^2+{({\sigma }_{zz}-{\sigma }_{xx})}^2+{({\sigma }_{xx}-{\sigma }_{yy})}^2]({\sigma }_{yz}{}^2+{\sigma }_{zx}{}^2+{\sigma }_{xy}{}^2)\end{array}$$

(9)

For directional columnar materials, the parameters in the yield function were presented as

$$\{\begin{array}{l}k=\sqrt{\frac{2}{3}}{\sigma }_{\mathrm{YV}}\hfill \\ I=\frac{4}{9}\frac{{\sigma }_{\mathrm{YV}}^4}{{\tau }_{\mathrm{YV}}^4}\hfill \\ J=\frac{32}{3}\frac{{\sigma }_{\mathrm{YV}}^4}{{\sigma }_{\mathrm{YH}}^4}-\frac{2}{3}-\frac{2}{3}\frac{{\sigma }_{\mathrm{YV}}^4}{{\tau }_{\mathrm{YV}}^4}\hfill \end{array}$$

(10)

where ${\sigma }_{\mathrm{YV}}$ and ${\sigma }_{\mathrm{YH}}$ were longitudinal and transverse tensile yield strength, and ${\tau }_{\mathrm{YV}}$ was longitudinal torsional yield strength.

Based on the above derivation, the equivalent stress of the directional solidification superalloy was rewritten as

$$\sigma =\sqrt{\frac{3}{2}}\sqrt[4]{\begin{array}{l}\frac{1}{9}{\left[{\left(s_{xx}-s_{yy}\right)}^2+{\left(s_{yy}-s_{zz}\right)}^2+{\left(s_{zz}-s_{xx}\right)}^2\right]}^2+I{\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)}^2\\ +\frac{J}{3}\left[{\left(s_{xx}-s_{yy}\right)}^2+{\left(s_{yy}-s_{zz}\right)}^2+{\left(s_{zz}-s_{xx}\right)}^2\right]\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)\end{array}}$$

(11)

3.2. Normalized Creep Model

The model developed in this article, based on normalized parameters, expressed creep strain as

$${\epsilon }_{\mathrm{c}}={\eta }_1(1-e^{-{\eta }_4\zeta })+{\eta }_2\zeta +{\eta }_3{\zeta }^{{\eta }_5}$$

(12)

where $\zeta =t/t_{\mathrm{c}}$ was dimensionless time, $t_{\mathrm{c}}$ was the endurance life at a given temperature and stress, and $\zeta \in \left[0,1\right]$. ${\eta }_i\left(i=1,2,3,4,5\right)$ was a function of dimensionless temperature. If the creep deformation test was conducted until fracture, the fracture time of the creep curve was the endurance life $t_{\mathrm{c}}$. Otherwise, the endurance life at this temperature and stress was obtained based on the endurance test or material manual data. The endurance life at a given temperature and stress was obtained through the endurance equation.

From the data provided in the material manual, it was acknowledged that the Larson–Miller (L-M) equation was adopted for most materials, such as GH4169 [31] and DZ125 [32]. For other materials, the Manson–Succop (M-S) equation was employed, such as K417 [33] and DD6 [30].

After normalization of experimental time by endurance life, the parameters ${\eta }_1$, ${\eta }_2$, and ${\eta }_3$ in the model had a clear physical meaning, as shown in Figure 6. The vertical axis in the figure represented creep strain (excluding instantaneous strain), with the specific meaning as follows:

(1)

At the time of , ξ = 1, ${\epsilon }_{\mathrm{r}}={\eta }_1+{\eta }_2+{\eta }_3$, ${\epsilon }_{\mathrm{r}}$ was the creep strain (fracture elongation) at rupture time.

(2)

The three terms in the model described the three stages of creep, respectively, and the parameters ${\eta }_1$, ${\eta }_2$, and ${\eta }_3$ represented the creep strain of the three stages. The steady-state creep strain rate ${\eta }_2$ after time normalization was also the slope of the straight line in the second stage of creep.

(3)

${\eta }_1$ was equal to the intercept of the second stage straight line on the vertical axis. ${\eta }_3$ was equal to the distance between the intersection of the second stage straight on the line of $\zeta =1$ and the point of rupture.

(4)

The curvature changes of the curves in the first and third stages of creep were controlled by ${\eta }_4$ and ${\eta }_5$, reflecting the speed of creep change in each stage.

After normalization, the variation ranges of temperature and stress were not large, but the five parameters (especially ${\eta }_4$) in the model were relatively large. The five parameters, which were considered to depend on temperature and stress, were taken as a natural logarithm, to express as the following form

$$\ln {\eta }_i=a_i+b_i\frac{T}{T_{\mathrm{m}}}+c_i\frac{\sigma }{{\sigma }_{0.2}}+d_i\frac{T}{T_{\mathrm{m}}}\frac{\sigma }{{\sigma }_{0.2}}$$

(13)

where ${\eta }_i\left(i=1,2,3,4,5\right)$ represented five parameters, $T/T_{\mathrm{m}}$ and $\sigma /{\sigma }_{0.2}$ represented normalized temperature and normalized stress, $\sigma$ was equivalent stress of directional solidification superalloy, and $a_i$, $b_i$, $c_i$, and $d_i$ were four sets of undetermined constant coefficients. Creep parameters at a given temperature and stress were obtained through Equation (13), with an ability to interpolate and extrapolate.

The work flow chart of the creep model is shown in Figure 7.

4. Creep Deformation Simulation of Turbine Blade

4.1. Parameters of Directional Solidification Superalloy DZ411

DZ411 is a nickel-based precipitation hardening directional solidification columnar superalloy with a melting point of about 1300 °C. The elastic properties of DZ411 [28] are listed in

Table 3. Elastic properties of DZ411, data from [28].

T/°C	${\mathit{E}}_{11}, {\mathit{E}}_{22}, {\mathit{E}}_{33}$/GPa	${\mathit{G}}_{12}, {\mathit{G}}_{23}, {\mathit{G}}_{13}$/GPa	${\mathit{\nu }}_{11}, {\mathit{\nu }}_{22}, {\mathit{\nu }}_{33}$
25	130	60.01	0.361
100	128	57.49	0.363
200	126	54.98	0.365
300	123	52.91	0.368
400	118	50.62	0.372
500	114	48.91	0.376
600	110	47.06	0.381
700	106	44.55	0.385
800	101	41.08	0.389
900	95	37.99	0.407
1000	86	34.48	0.418
1100	80	33.00	0.430
1200	75	32.00	0.450

. The yield strength data at different temperatures were obtained using linear interpolation fitting.

Endurance life curve equation (L-M equation)

$$\lg t=g_1+g_2\frac{1}{T}+g_3\frac{x}{T}+g_4\frac{x^2}{T}+g_5\frac{x^3}{T}$$

(14)

where $x=\lg \sigma$, $T=(9/5\theta +32)+460$ and the unit of $\theta$ is degrees Celsius. The fitting coefficients of equation are shown in

Table 4. Fitting coefficients (median) of DZ411 endurance equation.

${\mathit{g}}_1$	${\mathit{g}}_2$	${\mathit{g}}_3$	${\mathit{g}}_4$	${\mathit{g}}_5$
−20.3546	$2.3813\times {10}^6$		$1.0767\times {10}^6$	$-1.4160\times {10}^5$

.

Analyzing the creep fracture elongation in the tensile direction, the significant dispersion of creep fracture elongation was found, with little correlation to the levels of stress and temperature. Referring to Section 2.3, the average fracture elongation of DZ411 material was about 12.9%.

The five parameters, expressed in Equation (13) and listed in

Table 5. Fitting parameters of creep curves.

Serial Number	Temp/°C	Stress/MPa	${\mathit{\eta }}_1$	${\mathit{\eta }}_2$	${\mathit{\eta }}_3$	${\mathit{\eta }}_4$	${\mathit{\eta }}_5$
A-1	930	360	0.29	5.54	9.18	322.28	7.57
A-2	930	340	0.29	3.93	7.78	192.84	5.75
A-3	930	320	0.35	3.35	6.30	146.93	5.78
A-4	930	320	0.24	2.89	5.87	270.54	6.10
A-5	930	290	0.29	3.64	6.07	199.87	6.76
B-1	850	540	0.34	4.09	8.57	193.70	5.26
B-2	850	510	0.32	3.75	7.94	353.53	5.25
B-3	850	480	0.36	3.94	6.70	355.88	5.54
B-4	850	450	0.29	3.99	6.72	378.16	5.35
B-5	850	450	0.18	3.37	5.45	390.00	5.29

, in the creep model were taken as natural logarithm, which was considered to depend on temperature and pressure. The fitting results are shown in

Table 6. Fitting parameters of creep model.

Parameter	${\mathit{a}}_{\mathit{i}}$	${\mathit{b}}_{\mathit{i}}$	${\mathit{c}}_{\mathit{i}}$	${\mathit{d}}_{\mathit{i}}$
${\eta }_1$	−34.2135	43.1135	41.0394	−53.6568
${\eta }_2$	22.6810	−30.6188	−27.6279	39.6027
${\eta }_3$	5.7368	−07.7139	−05.8624	11.2015
${\eta }_4$	73.2063	−91.1465	−78.1460	105.3585
${\eta }_5$	5.2302	−04.9121	−07.7107	10.7151

.

In the second section, a yield criterion for directional solidification superalloys was proposed. The yield function without the tension-compression asymmetry was written in the following form

$${(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)}^2+I{(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)}^2+J(s_{xx}{}^2+s_{yy}{}^2+s_{zz}{}^2)(s_{xy}{}^2+s_{yz}{}^2+s_{xz}{}^2)=k^4$$

(15)

According to the uniaxial tensile test data of directional solidification superalloy, the parameters in the yield function were expressed as

$$\{\begin{array}{l}{k}^4=\frac{4}{9}{\sigma }_{0.2}{}^4\hfill \\ I=0.0019T+1.6884\hfill \\ J=-0.0373T+45.84\hfill \end{array}$$

(16)

The yield surface was considered to be similar to the potential energy surface. Equation (15) was applied to the creep deformation process of directional solidification materials. Finally, the relevant parameters of DZ411 material in this section were applicable for the endurance life prediction.

4.2. Programming and Verification of Subroutine for Normalized Parameter Model

A usermat subroutine based on ANSYS was adopted and compiled. The involved usermat subroutine mainly completed two functions: generating an elastic matrix and updating stress. The iterative process for solving the internal force displacement relationship [34] is shown in Figure 8.

The iterative method was used to update the stress. The invocation process and the iteration process of subroutines are shown in Figure 9. The variable directly related to the nickel-based superalloy model during the iteration process was the equivalent creep strain rate ${\dot{\epsilon }}_c$.

4.3. Simulation of Creep Deformation Behavior of Turbine Blade

4.3.1. Turbine Blade Grid and Boundary Conditions

The real turbine blade grid, as shown in Figure 10, was used for creep deformation simulation.

The boundary conditions of the turbine blade were set as follows:

(1)

The first pair of tenon teeth were constrained on the normal displacement, respectively.

(2)

The second, third, and fourth pairs of tenon teeth were subjected to normal pressures of 180, 160, and 140 MPa, respectively.

(3)

The groove displacement was constrained at the baffle.

(4)

The given speed was 9500 r/min.

(5)

The typical temperature characteristics [35] at the tenon’s root, blade body root, and blade body top were given as 550, 650, and 1000 degrees Celsius, as shown in Figure 10b.

It is believed that the centrifugal force borne by the first pair of tenon teeth is the largest among the four pairs, so the surface displacement of the first pair was constrained in the normal direction. The other three pairs of teeth had a decreasing bearing area and capacity, with a decreasing applied load.

4.3.2. Stress Analysis of Turbine Blade

The crystal axis <100> direction was selected as the blade height direction in turbine blade stress analysis, without casting error. The finite element software was used to calculate the stress distribution, and the calculation was set as example A. The finite element analysis results showed that the maximum stress point of the blade (assessment point L) was located at the bottom of the turbine blade tenon, and the equivalent stress nephogram is shown in Figure 11a. The maximum stress point of the blade body (assessment point M) was located at the turbulence column near the root of blade, and the equivalent stress nephogram is shown in Figure 11b.

4.3.3. Creep Analysis of Turbine Blade

The subroutine was used to simulate the creep behavior of example A for 10,000 h, with a focus on stress concentration areas such as tenon and spoiler columns [36]. The creep strain and equivalent stress curves versus time at the stress concentration areas of the spoiler column are detailed in Figure 12. It was noted that the maximum creep strain area of the overall blade was not located at the maximum stress area, but at the root of the uppermost spoiler column (assessment point N). The equivalent stress in this area was low, but the environmental temperature was high and the creep strain was large. The creep strain and equivalent stress nephograms at this point are plotted in Figure 13, and the curves of the equivalent stress and creep strain at this point are plotted in Figure 12.

It was observed that the stress concentration area at the bottom of the tenon had lower temperature and higher stress, with a strain of about 0.06% after a 10,000 h creep process. The stress concentration area of the spoiler column had a higher temperature and lower stress, and the strain at the root area of the top spoiler column was also larger after creep for 10,000 h. The roots of spoiler column should be paid attention to in actual turbine blades. With the continuous progress of creep, there was a significant stress relaxation phenomenon at the areas with obvious creep strain. The equivalent stress at the upper spoiler column decreased from 73 MPa to 51 MPa, with a stress relaxation of about 30%. The creep strain at the lower spoiler column was lower, and the equivalent stress decreased from 205 MPa to 201 MPa, with a stress relaxation of 2%. The creep strain at the bottom of the tenon was too low, and there was no significant change in equivalent stress. As the creep time increased, the stress further relaxed.

4.3.4. Endurance Analysis of Turbine Blade

According to Equation (12) and Table 4, the stress levels before and after 10,000 h creep at assessment points L, M, and N were evaluated for their endurance life. The predicted life results are shown in

Table 7. Predicted endurance life of assessment points.

Assessment Point	Temperature /°C	Without Creep		After 10,000 h Creep
		Equivalent Stress /MPa	Predicted Endurance Life/h	Equivalent Stress /MPa	Predicted Endurance Life/h
L	550.0	354.5	$6.97\times {10}^{12}$	354.5	$6.97\times {10}^{12}$
M	731.5	205.0	$1.71\times {10}^8$	200.7	$2.00\times {10}^8$
N	937.8	072.7	$4.53\times {10}^5$	051.2	$9.00\times {10}^6$

.

From the table, it is evident that the endurance life of assessment points M and N markedly increased. The endurance life of assessment point N increased from about $4.53\times {10}^5$ h to about $9.00\times {10}^6$ h, with the latter being about 20 times that of the former. If the stress relaxation effect caused by creep was not considered, the predicted structural life would be too conservative. However, due to the existence of stress relaxation phenomenon, the results obtained by adopting the stress after a certain period of creep for durability life prediction were closer to the actual results. The stress relaxation of assessment point L was only 2%, with an addition of 17% predicted endurance life. For the assessment point N, the creep strain caused by the stress/temperature combination was very small, without obvious stress relaxation phenomena. The prediction of the endurance life for point N remained unchanged.

Based on the above analysis, it could be concluded that the stress relaxation effect existed in turbine blades used for long-term service. The stress level was reduced in the stress concentration area, accompanied by the increase in endurance life. The more obvious the stress relaxation phenomenon was, the more significant the improvement in endurance life was. In the designing and analyzing of actual structures, creep should be taken into account to reduce unnecessary margins and ensure the accuracy of prediction results.

5. Analysis of Crystal Axis Deviation Angle of Turbine Blade

It is difficult to precisely control the crystallization direction of turbine blades in the casting process. There is a certain deviation angle between the crystallization direction of the directional solidification alloy and the blade height direction. In engineering, it is generally believed that a deviation angle within 10° [37] or 15° [38] is an allowable range.

5.1. Definition of Crystal Axis Deviation Angle

In the actual casting process of turbine blade production, it is hard to ensure that the blade height direction fully coincides with the crystal axis direction. With a deviation angle between the crystal axis direction and the blade height direction, the performance of the turbine blade is changed. The larger the deviation angle, the greater the impact on the creep performance of turbine blades [39]. In light of the asymmetry of turbine blades, the deviation directions have different effects on the performance of the blade. The deviation angles in different directions were analyzed. The deviation angles of 5° and 10° for <100> crystal axis were analyzed for the creep performance, along with four orthogonal deviation directions. The schematic diagram of crystal axis deviation angle is shown in Figure 14. The red arrows in the figure represent the deviation directions related to the <100> crystal axis, while B, C, D, and E represent four orthogonal directions. The actual deviation angles between the crystal axis direction and the blade height direction are all 5° and 10°.

5.2. Comparison of Finite Element Analysis Results for Crystal Axis Deviation Angle

The results of multiple examples are compared in Table 7, and the equivalent stress and creep strain at the maximum creep strain positions of different examples were extracted for comparison. As described in

Table 8. Comparison of results of stress concentration position after 10,000 h creep.

Example	Deviation Angle /°	Equivalent Stress /MPa	Creep Strain /%
A	00	51.23	0.0592
B1	05	52.14	0.0572
B2	10	52.11	0.0592
C1	05	59.89	0.0594
C2	10	59.05	0.0614
D1	05	56.66	0.0642
D2	10	47.60	0.0701
E1	05	49.72	0.0650
E2	10	47.87	0.0708

, the maximum creep strain of Example A after 10,000 h of creep was 0.0592%. After the existence of a deviation angle, the creep strain of all other examples except Example B1 increased, reaching a maximum of 0.0701%, with an increase of 18.4%. The deviation angle of Example B1 was 5°, and the creep strain was 0.572%, which was 3.4% lower than Example A.

As the chart shows, the material properties varied with angle and orientation, affecting the performance of the entire blade. The crystal axis <100> orientation was the direction of centrifugal load. Without a deviation angle in the crystal axis, the creep resistance was excellent, and the creep strain was small. As the orientation of the crystal axis changed, the creep strain increased or decreased within a certain range of deviation angle, but the overall trend was up. In other words, when the crystal axis orientation coincided with the leaf height direction, the creep performance of the material may not necessarily be optimal, but it could be considered relatively superior overall.

The deviation angles of examples B2, C2, D2, and E2 were all 10°, which was the maximum value within the engineering allowable range of 10°. The maximum creep strain change amplitude did not exceed 20%. The maximum creep strain of Example B1 was reduced by 3% compared to Example A. Considering the complexity of the geometric structure of the turbine blade, the coincidence of the crystal axis direction with the blade height direction may not be the optimal crystal axis position. The crystal axis <100> position of Example A was only near the optimal orientation.

In practical engineering applications, the acceptable allowable range of deviation angle is selected as 10°, which can not only meet the material performance requirements of engineering, but also reduce the defect rate of turbine blade products and improve economic benefits.

Furthermore, the equivalent stress results, deviating from different angles along the deviation directions of Examples B and D, were calculated, with 10,000 h creep. After the deviation angle exceeded 10°, as the deviation angle between the crystal axis and the blade height direction increased, the local creep strain increased, and the creep performance of the material decreased. The maximum creep strain value with the deviation angle is plotted in Figure 15, and the maximum local creep strain value of the turbine blade is shown in

Table 9. Creep strain at different deviation angles.

Direction	Deviation Angle/°	Creep Strain/%	Direction	Deviation Angle/°	Creep Strain/%
B	00	0.0594	D	00	0.0594
	05	0.0572		05	0.0642
	10	0.0642		10	0.0701
	15	0.0642		15	0.0750
	20	0.0696		20	0.0777
	25	0.0729		25	0.0775
	30	0.0731		30	0.0749

.

As the deviation angle exceeded 10°, the creep performance of the material deteriorated severely, and the creep strain rapidly increased with further increase in the angle. At a deviation angle of 20°, the maximum creep strain increased by 31%. In the actual casting process, it is necessary to control the deviation angle within 10° to avoid a serious decrease in the creep resistance of the turbine blade.

The increase in creep strain was caused by the deterioration of material properties. For a directional solidification superalloy, the creep performance in the crystalline direction is generally the best, and the endurance life is the longest. With a deviation angle over 10°, the creep performance decreases, as does the endurance life. For the same creep time, the creep strain increases.

The fundamental reason is the presence of the shear stress component in the crystal axis coordinate system. In the directional solidification equivalent stress defined in this article, under multiaxial stress states, the contribution of the shear stress component is greater than that of the tensile stress component in most cases (only a small portion of cases close to pure shear are exceptions). As a consequence, the equivalent stress increased and the endurance life decreased.

6. Conclusions

A three-stage creep model was adopted in this article for simulating the creep deformation of actual turbine blades based on finite element analysis. The creep tests were conducted to complete modeling and the subroutine was compiled for subsequent simulation. According to the creep subroutine, the phenomenon of a deviation angle between the crystal axis and blade height in engineering structures was analyzed. Although there are some differences in the boundary conditions employed in this article and the actual turbine blade, the results are strict and have generality. The following conclusions can be drawn on the preceding analysis:

(1)

The yield function suitable for nickel-based high-temperature alloys was proposed by order upgrading. The creep subroutine was compiled combined with the three-stage normalized creep model. The yield function and creep subroutine had high computational efficiency and robustness, meaning they are suitable for the creep analysis of complex turbine blades.

(2)

In practical engineering structures, stress concentration areas are inevitable. For turbine blades, the important locations, such as the roots of the spoiler columns or the corners of the cooling channels, deserve special attention. At the service temperature, the stress concentration areas will experience stress relaxation with the increase in creep strain, causing stress redistribution.

(3)

The stress redistribution in the vicinity of the assessment points causes changes in the material strength, which is beneficial for engineering structures and the endurance life of assessment points. Therefore, creep analysis should be conducted first when analyzing the endurance life of turbine blade. The simulation results after creep analysis should be used for strength verification, instead of the static strength analysis results.

(4)

There is a certain deviation angle between the crystal axis direction and the blade height direction of the turbine blade material. The finite element simulation results proved that a 10° deviation angle within the engineering allowable range is appropriate. This allowable range not only meets the material performance requirements of engineering application, but also enhances the economic benefits of the product.