UDK 539.37:519.61/.64
Original scientific article/Izvirni znanstveni članek
ISSN 1580-2949 MTAEC9, 46(3)197(2012)
COMPUTER SIMULATION OF FATIGUE, CREEP AND THERMAL-FATIGUE CRACKS PROPAGATION IN GAS-TURBINE BLADES
RAČUNALNIŠKA SIMULACIJA NAPREDOVANJA UTRUJENOSTNIH RAZPOK, RAZPOK PRI LEZENJU IN TERMIČNO-UTRUJENOSTNIH RAZPOK V LOPATICAH PLINSKIH TURBIN
Artem Semenov1, Sergey Semenov1, Anatoly Nazarenko1, Leonid Getsov2
1St. Petersburg State Polytechnical University, St. Petersburg, Russia 2NPO CKTI, St. Petersburg, Russia guetsov@yahoo.com
Prejem rokopisa - received: 2011-07-06; sprejem za objavo - acce^^ed for publication: 2012-03-01
Methods and computational algorithms for determining the growth rate of fatigue creep and thermal-fatigue cracks are considered. The rate of crack growth is dependent on the stress-intensity factor (or J-integral) for fatigue, on the C*-integral for creep and on the stress-intensity factor (or J-integral) and C*-integral for thermal fatigue. Simulations of the crack propagation under fatigue, creep and thermal fatigue at the edge of the blade of a gas turbine are carried out and discussed.
Keywords: fatigue, creep, thermal fatigue, crack, C*-integral, turbine blades, finite element simulation
Članek obravnava metode in računske algoritme za določanje hitrosti rasti utrujenostnih in toplotno utrujenostnih razpok pri lezenju. Hitrost rasti razpoke je odvisna od intenzitete napetostnega faktorja (ali J- integrala) pri utrujenosti, od C*-integrala pri lezenju in od intenzitete napetostnega faktorja (ali J-integrala) in od C*-integrala pri toplotni utrujenosti. Predstavljene in komentirane so simulacije širjenja razpoke po robovih lopatic plinske turbine pri utrujenosti, lezenju in toplotnem utrujanju. Ključne besede: utrujenost, lezenje, toplotna utrujenost, razpoka, C*-integral, turbinske lopatice, metoda končnih elementov
1 INTRODUCTION
Blades of a gas-turbine engine (GTE) are subjected to extreme non-steady operating conditions which can give rise to cracking. With the non-destructive testing of the turbine blades after operation, cracks of various sizes, configuration and locations could be detected. Fracto-graphic studies allow us to identify the nature of the detected cracks. Usually, the cracks are1 because of:
•	fatigue (Figure 1a),
•	creep fracture (Figure 1b),
•	thermal fatigue (Figure 1c).
For an accurate assessment of the durability and life prediction of blades, it is expected during calculations to take into account not only the stage of the crack initiation, but also the stage of the crack propagation and consider the differentiation of specific types of cracks and their growth patterns. To estimate the remaining life of individual blades with a known configuration of crack, the most reliable prediction is based on calculations in the context of a three-dimensional analysis, which includes the kinetics of the crack growth and takes into account changes of the shape of the crack front and of its growth direction.
Figure 1: Cracks of different nature in gas-turbine blades: a) fatigue crack; b) creep crack; c) thermal fatigue cracks
Slika 1: Razpoke različne narave v turbinskih lopaticah: a) utrujenostna razpoka; b) razpoka zaradi lezenja, c) razpoke zaradi termične utrujenosti
The accumulation of experimental data on the rate of the propagation of cracks of different natures, the development of improved phenomenological models of inelastic deformation, improvement of methods of computational mechanics and computer technology advances make it possible to implement solutions to complex, non-linear, boundary value problems arising with the modelling of crack propagation in turbine blades. In some cases, a direct mathematical modelling of crack growth eliminates costly and time-consuming experiments to confirm the state of the blades.
The aim of this work was to develop and implement methods for the resource calculation of the blades in which cracks were detected during the control. The approach is based on a direct-step simulation of the crack propagation using the finite-element method (FEM) and experimental data on the dependence of the crack growth rate from the stress-intensity factor (SIF) amplitude ^K and the C*-integral obtained for the blade material. The main steps of the practical implementation techniques are discussed below:
•	determination of the size and the crack location in the blade based on the results of an inspection (or analysis of statistical data on failures) and the identification of the nature of the defects identified with fractographic methods;
•	identification of the alleged operation modes GTE, which caused the formation of the cracks;
•	solving problems of heat conduction and aerodynamics in order to determine the distribution of the temperature fields in the bulk of the blade and the distribution of the gas pressure on the surface of the blade, as well as their dependence on time;
•	solution to the problems of thermo-elasticity, thermo-elasto-plasticity and creep to determine the stress-strain state of the blades in the presence of growing cracks;
•	computational-experimental determination of the rate of growth of cracks in different modes (based on the diagrams AK - da/dN, C* - da/dt at different operating temperatures);
•	identification of the critical state of the blade that causes its destruction (or maximum values of permissible stress that could be reached according to the strength standards and material characteristics);
•	wording of the requirements for the minimum time before the next audit.
2 COMPUTATIONAL METHOD FOR THE CRACK GROWTH PREDICTION IN TURBINE BLADES
This methodology is used to determine the kinetics of crack propagation, estimate the number of cycles (or time) to reach a critical crack length (or to determine its length for a given number of load cycles or duration of operation). The initial distribution of defects is accepted
as surface cracks of a given length and the calculations are made according to actual and (or) predictive models of operation (modes of loading, the load levels).
2.1 Models for the c^ack propagation rate
To determine the growth rate of fatigue cracks the following Paris power-type approximation was used:
da
- = B( AK^f) -
dN
(1)
where B and m are material constants, AKeff — Kmax — Kop < AK — Kmax — Kmin are the effective scope of the SIF, which in the simplest model that takes into account only the stage of steady growth and simply reflects the effect of the crack closure, defined by the relations:
AKff =
Kmax - K^„, R = ^Km^ ^ 0,
Kmax
Kmax, R= ^Km^ < 0, K
(2)
In the presence of additional experimental information on the form of the kinetic diagrams of fatigue fracture (KDFF), more complicated equations than (1) may be used. The effect of cycle asymmetry and the presence of transient regions in the KDFF can be taken into account, for example, by using the equations: Forman2
_da _ B(AK)-dN " (1-R)K^-AK
Walker3
Aa_ dN
_ b[(1-r) " ak]"
(3)
(4)
or other, more complex, equations4-8.
While analyzing the growth of short cracks in modes of loading, corresponding to the threshold region SIF, it is necessary to take into account the specific nature of KDFF and, for instance, explicitly use the terms da/dN — 0 when AK < AKth and da/dN ^ 0 under AK > AKth — the threshold range of SIF, which depends on the material, on the cycle asymmetry and on the aggressive environment.
To determine the creep crack growth rate, the following expression is used:
da dt
_ A(C *) ^
(5)
where A and q are the material characteristics depending, in general, on the temperature and C*-integral9, which is invariant when considering the steady creep stage. In general, the three-dimensional case for an arbitrarily oriented crack with a curved edge, uses a vector-integral defined by:
c; -j
dil J
^ *	^ J
dsS
(6)
where S is the surface, covering the front of the crack and nk - kth-vector component normal to the surface.
The parameters of growth of a thermal fatigue crack (low-cycle fatigue) for an arbitrary cycle form were defined using the expression:
da_ dN
- B(AKff)" +_ A(C *(r)) ^ dr
(7)
where the values of the material parameters B, m, A and q are the same as in equations (1) and (5). The integration is performed within one cycle (from 0 to rc). The first term in (7) characterizes the growth of thermal fatigue cracks due to thermoelastic stresses during starting and stopping of GTE and the second, the growth of cracks in operating conditions between starts. It should be noted, however, the following features of AKeff value in low-cycle (thermal) fatigue: during thermocycling the cycle of stress tends to a symmetry1 and during the half cycle of compression the crack closes. For irregular regimes of thermal cycling under conditions of frequent and abrupt changes of level and duration of exposure, instead of (7) the following expression should be used:
da - B(AKff)m dN + A(C * (r))q dr
(8)
2.2 Calculation of fatigue c^ack growth kinetics
the original direction of growth based on the criterion of maximum tensile stress Ki sinAÖ + Kii (3 cosAÖ - 1) = 0. Hence, we have:
Tl^1+8( KII / K i)2
The increment of crack length is determined by the relations:
Aa=
B( AKeff)m AN, AKeff > AK th
0
AKeff ^AK th
(9)
The problem is solved in a linear-elastic formulation under the assumption of small strains. The external exposure is considered as the action of centrifugal forces, gas pressure and vibration loads. If necessary, the influence of the temperature fields can be accounted, also. In the FE model of the blade fracture, the crack is defined geometrically by introducing at different banks nodes with identical coordinates but with different numbers. When meshing the region with finite elements, it is desirable to use a focused mesh around the crack's front line and use hexahedral finite elements.
The whole operation period is divided into intervals (with a given number of cycles AN), each interval for a typical cycle and two elastic problems are solved in the presence of cracks. First, the maximum value of SIF Kmax is determined, corresponding to the positive direction of the application of vibration loads. Then the minimum value of SIF Kmin is determined, corresponding to the opposite direction of application of vibration loads. The values of Kmax and Kmin are defined for each node at the front of the crack.
During the determination of the maximum SIF in the case of the conditions KI >> KII and KI >> KIII, the direction of crack growth is preserved (crack of normal separation, I mode). When these inequalities violate the correction of the crack-growth direction AO should be accounted for by determining the angle of deviation from
Based on the values of crack increment, the crack length for the next iteration,
ai+1 - ai +Aa	(10)
is determined, the FE mesh is modified and the previous steps of the calculation are repeated until the critical crack length is reached.
When using the Paris law in order to minimize the computational cost, the increment of crack length is determined for the most loaded point of the front using the expression:
Aa - A(C ;)q At	(11)
The calculation of the SIF K,, K,,, K,,, is based on an analysis of the distribution of crack displacement fields in the vicinity of its tip. Except for the extreme points of the crack front, the asymptotic behaviour of stresses in the region near the crack tip is assumed to be flat deformable. When using the FE software ANSYS version 12 and above, the SIF can be considered automatically.
2.3 Calculation of creep crack growth kinetics
In solving the boundary problems in a FE analysis a nonlinear viscoelastic material model with a steady-state creep law for stage II is used (the Norton's law).
The whole operation period is divided into time steps At and the stress-strain state of the blade is determined at each interval. Based on the obtained values, the C;-inte-gral is calculated for each node at the crack front. The increment of crack length is determined by the integration of equation (5). Using the explicit Euler method, we obtain the expression:
Aa - Aa„
K
V KImax 7
(12)
Based on the values of the crack increment, the crack length for the next iteration is determined with (10), for which the FE mesh is modified and the previous steps of the calculation are repeated until the critical crack length is reached.
The order of the calculation is as follows:
•	Creation of a FE model with a crack.
•	Setting the properties of the material.
m
•	Evaluation of C*-integral.
•	Determination of crack-length increment.
The calculations of C*-integral can be made automatically using the FE software package ABAQUS version 6.10 and above.
2.4 Calculation of thermal fatigue c^ack growth kinetics
When the loading conditions consist of relatively short start/stop periods versus the exposure time at elevated temperature, the solution may be simplified to a formulation based on a separate consideration of the exposure time using the creep model and the elasto-plastic model for periods of start/stop. A further simplification is possible while analyzing the start/stop period when the stresses, remote from the crack tip, do not exceed the yield stress and the elastic material model may be used for the calculations.
The whole operation period is divided into intervals (with a given number of cycles AN, the duration of each cycle rc). At each interval and for a typical cycle in the presence of a crack in the blade, two boundary problems are solved: the analysis of the creep processes within the cycle and the analysis of the fatigue in the thermo-elastic formulation to solve the problem when the engine is stopped (cooling). During start-up (heating) the surface cracks tend to be closed. The increment of crack length is determined by integrating (7) and for one typical cycle (block of cycles with similar values of AK^jf and C*) we have:
Aa =
r
B(. AK eff)'

A(C *(r))
1
'dr AN
(13)
Based on the values of crack increment, the crack length for the next iteration is determined with (10), for which the FE mesh is modified and the previous steps of the calculation are repeated until the critical crack length is reached. The calculation of the parameters characterizing the fatigue-crack growth may be performed using the FE software packages ANSYS or ABAQUS, and the parameters related to creep, with the help of ABAQUS.
3 SIMULATION OF CRACK GROWTH IN A TURBINE BLADE
Let us consider the results of calculations made for the blade in Figure 2. In these calculations we assumed that the crack (idealized defect) is located at the output edge in the plane orthogonal to the blade and has an initial length of 1 mm. The fatigue, creep and thermal fatigue crack are located, respectively, at a height of (15, 50 and 80) mm (1/3 height of the blade) from the root section of the blade.
The blade was fixed in the direction of its axis over all nodes on the lower grounds and also three degrees of freedom in the plane of ground are fixed for the elimi-
Figure 2: FE model of the blade with the crack (1 mm initial length) located at a height of 15 mm from the root section Slika 2: FE-model lopatice z razpoko (1 mm začetna dolžina) na razdalji 15 mm od korenskega dela
nation of rigid body translations in this plane and rotation around its axis. The calculations used a model of linear elastic isotropic material (for the fatigue cracks), Norton's model (for the crack creep) and the model thermo-visco-elastic-plastic (for the thermal fatigue cracks). The problems were solved under the assumption of small strains. The material parameters were determined on the basis of references (for operating temperature at 850 "C)10,11.
3.1 Simulation of the fatigue crack growth
The FE model of the blade with the crack of initial length is shown in Figure 3. The FE mesh in the cracks plane of blades consists of quadratic isoparametric 20-node elements. The parameters of the FE models with an initial fatigue crack are given in Table 1.
Table 1: Parameters of FE model for the blade with a crack Tabela 1: Parametri FE-modela za lopatico z razpoko
For the initial crack length		For a crack length of 1.7 mm	
Number of nodes	73 197	Number of nodes	65 021
Number of elements	16 382	Number of elements	14 476
Number of degrees of freedom	219 591	Number of degrees of freedom	195 063
In this problem the load was been taken as follows: action of centrifugal forces (m = 3000 rpm), gas pressure on the lateral surface (p (x, y, z),), vibration loads (± Fx, ± Fy). The distribution of the stress intensity fields for the case of a positive direction of the application of vibration loads is shown in Figure 4.
The SIF values, computed for the case of positive direction of application of vibration loads (the definition of Kmax), are presented in Table 2 for all the corner nodes lying on the front of the crack.
Figure 3: Edge of the blade and fatigue crack front (in the plane of propagation): a) initial configuration, b) crack after 107 cycles Slika 3: Rob lopatice in ~elo razpoke (v ravnini propagacije): a) za~etna konfiguracija in b) razpoka po 107 ciklih
Figure 4: Field distribution of von Mises stress intensity in the blade with: a) fatigue-crack length of 1 mm and b) 1.7 mm after 107 cycles Slika 4: Polje razdelitve von Misesove intenzitete napetosti v lopatici z utrujenostno razpoko z: a) dolžino 1 mm in b) 1,7 mm po 107 ciklih
Table 2: SIF values for the case of a positive direction of the application of vibration loads (definition of ^max) for a fatigue-crack length of 1 mm. The distance is measured along the crack front (from the back through).
Tabela 2: SIF-vrednosti za primer pozitivne smeri obremenitve zaradi vibracij (definicija Kmax) za dolžino razpoke 1 mm. Dolžina je merjena vzdolž ~ela razpoke (od hrbta skozi).
Distance (mm)	Ki/ (MPa/m05)	kii/ (MPa/m05)	kiii/ (MPa/m0.5)
0	4.20	0.27	0.02
0.56	8.92	0.51	0.05
1.11	10.76	0.66	0.04
1.67	12.17	0.76	0.03
2.22	12.92	0.82	0.06
2.78	13.12	0.85	0.07
3.33	12.54	0.84	0.06
3.89	11.15	0.83	0.06
4.44	5.66	0.60	0.05
Figure 5: Evolution of fatigue-crack front with an increasing number of cycles
Slika 5: Evolucija ~ela utrujnostne razpoke pri pove~anju števila ciklov
the FE solution of the nonlinear boundary value problem for the stress, a step-incremental iterative procedure was used. The value of the initial size of the time step was assumed to be equal to 10-10 s. In the process of the boundary solution, the step size was adaptively changed and increased, gradually. The precision of the creep strain was assumed to be of 10-6, which corresponds to an error in the calculation of stresses 0.12 MPa and the corresponding integrals of the asymptotic values of C*-integrals were determined based on the resulting FE solutions of the nonlinear boundary problem for the stress fields and displacements (and velocities). The time dependence of changes for the five characteristic points on the front of the crack, as shown in Figure 6. The crack growth during the year is shown in Figure 7.
The change of position of the crack front with the increase in the number of cycles is shown in Figure 5. Originally, a straight edge is close to semi-elliptical for N = 107 and with a further increase to N = 2 ■ 107 a progressive development of cracks in the region adjacent to the back is observed.
3.2 Simulation of creep c^ack growth
As a model problem, a blade with a crack 1 mm long, located on the edge of the output at 80 mm from the root section was chosen. The crack was identified as described above.
The creep-crack growth was investigated for a time of 100 000 h, which roughly corresponds to 11 years. In
Figure 6: Dependence of CfJj-integral on the time for the characteristic points on the front of the crack
Slika 6: Odvisnost C(t)-integrala za karakteristi~ne to~ke ~ela razpoke
Figure 7: Evolution of the crack front Slika 7: Evolucija ~ela razpoke
3.3 Simulation of thermal fatigue crack growth
Schematic representation of the effect of inhomoge-neity of temperature field, observed at the start and stop of GTE, was set in accordance with Figure 8.
The determined values of the scale SIF at the first stage of the modelling process of the growth for all the nodes lying on the front of the crack are shown in Table 3. It should be noted that the condition of fatigue-crack propagation AK > AKth is satisfied for all nodes on the crack front. Table 3 shows the crack-growth rate calculated in accordance with (7), as well as separate components of the crack rate caused by fatigue and creep.
It should be noted that in this case the component due to creep dominates. In Table 3, also the values of the increment of crack length after 1.5 ■ 104 cycles, calculated using equation (5), are shown. The increment of the crack length Aa (based on 1.5 ■ 104 cycles) from the arc coordinate, is calculated in this case along the crack front (for different front sites). Based on the values of crack increment, the crack length for the next iteration is determined with (10), for which the FE mesh is modified and the previous steps of the calculation are repeated until the critical crack length is reached.
Figure 8: Schematic representation of the effect of the inhomogeneity of the temperature field in the blade during starting up and shutting down the turbine
Slika 8: Shemati~na predstavitev vpliva nehomogenosti temperaturnega polja v lopatici pri zagonu in zaustavitvi turbine
4 CONCLUSIONS
1.	The methods and computational algorithms for the simulation of propagation process of fatigue, creep and thermal fatigue cracks in the blades were developed and verified for an uncooled stationary GTE turbine blade. Finite-element calculations were performed using the ANSYS and ABAQUS FE software.
2.	It was established that the shape of front of cracks of different nature varies significantly with the accumulation of the number of cycles and the operating time.
3.	It is shown that for the considered blade geometry and the assumed initial crack configurations, the values of KI are dominating on KII and KIII.
4.	During the thermal cycling, the maximum SIF scale determining the components of fatigue-crack growth are observed at the stop of GTE.
Table 3: SIF AKi, contributions from fatigue, creep, and the total value of thermal fatigue crack growth rate Aa/AN, as well as the increment of crack length Aa for 1.5 ■ 104 cycles in the first step of modelling the growth (for crack length of 1 mm)
Tabela 3: SIF AKj, deleži utrujenosti, lezenja in skupna velikost hitrosti rasti utrujenostne razpoke Aa/AN in prirastek pove~anja dolžine razpoke Aa za 1,5 ■ 104 cikle v prvi stopnji modeliranja rasti razpoke (za razpoko z dolžino 1 mm)
Distance along the crack front from the back through (mm)	AKi/ (MPa m0 5)	Aafat/AN/ mm per cycle	Aacreep/AN/ mm per cycle	Aa/AN/ mm per cycle	Aa (based on 1.5 ■ 104 cycle), mm
0	12.22	5.53 ■ 10-8	0	5.53 ■ 10-8	0.0008
0.56	26.22	1.10 ■ 10-6	3.50 ■ 10-6	4.60 ■ 10-6	0.0690
1.11	30.41	1.97 ■ 10-6	1.79 ■ 10-5	1.99 ■ 10-5	0.2980
1.67	32.93	2.69 ■ 10-6	2.97 ■ 10-5	3.24 ■ 10-5	0.4860
2.22	33.61	2.92 ■ 10-6	3.07 ■ 10-5	3.36 ■ 10-5	0.5040
2.78	32.78	2.65 ■ 10-6	2.38 ■ 10-5	2.65 ■ 10-5	0.3970
3.33	30.17	1.91 ■ 10-6	1.04 ■ 10-5	1.23 ■ 10-5	0.1850
3.89	25.85	1.04 ■ 10-6	1.58 ■ 10-7	1.20 ■ 10-6	0.0180
4.44	11.97	5.11 ■ 10-8	0	5.11 ■ 10-8	0.0008
5.	The areas of possible practical applications of the developed techniques are suggested.
6.	The practical implementation of the calculations of the blade crack requires the systematic accumulation of experimental data for the blade material for the construction of diagrams ^K - da/dN, C* - da/dt under different operating temperatures.
5 REFERENCES
1L. B. Getsov, Materials and durability of parts for gas turbines. Fourth edition in two volumes. Rybinsk. Ed. House gas turbine technology, 2010
2R. G. Forman, V. E. Kearney, R. M. Engle, Numerical Analysis of Crack Propagation in a Cyclic-Loaded Structure, Journal of Basic Engineering, Transactions of the ASME, 1967
3	K. Walker, The Effect of Stress Ratio During Crack Propagation and Fatigue for 2024-T3 and 7075-T6 Aluminum. In: Effects of Environment and Complex Load History on Fatigue Life, ASTM STP 462 1970, 1-14
4	F. Erdogan, M. Ratwani, Fatigue and fracture of cylindrical shells containing a circumferential crack, International Journal of Fracture Mechanics, 6 (1970), 379-392
5	Jr. Newman, A crack closure model for predicting fatigue crack growth under aircraft spectrum loading. In Methods and Models for predicting Fatigue Crack Growth under Random Loading. ASTM STP 748 American Society for Testing and Materials, Philadelphia, PA., 1981, 53-84
6	V. S. Balin, G. G. Madyakshas, Strength, durability and crack resistance of structures with long-term cyclic loading. St. Petersburg: Polytechnics, 1994, 204C
7	A. V. Prokopenko, V. N. Trade, L. B. Getsov et al., Influence of the loading regime on the growth rate of fatigue cracks in stainless steels in air and sea salt solution, Strength of Materials, (1983) 12, 41-45
8	A. V. Prokopenko, Method of determining the stress intensity factor in spades GTE, Strength of Materials, (1984) 4, 21-24
9	J. D. Landes, J. A. Begley, A fracture mechanics approach to creep crack growth. In: Mechanics of Crack Growth, ASTM STP 590. Am. Soc. Testing Mat. 1976, 128-148
10	Y. A. Nozhnitsky, E. R. Golubovsky, On the Strength Reliability of single-crystal turbine blades of high prospective GTE. In the book. Strength of materials and resource elements of power equipment, Proceedings of CKTI vyp.296, St. Petersburg, 2009, 74-82
11	M. Tabuchi, K. Kubo, K. Yagi, A. T. Yokobori, A. Fuji, Results of Japanese round robin on creep crack growth evaluation methods for Ni-base superalloys, Engineering Fracture Mechanics, 62 (1999), 47-60