Numerical Modeling of Cyclic Creep-Fatigue Damage Development for Lifetime Assessment of Steam Turbine Components

The paper presents the results of numerical analysis of creep-fatigue damage development in a steam turbine rotor under cyclic duty. Investigations were performed for a full cycle representing the most accurate description of real operation and three simplified cycles representative of different levels of simplification in inelastic strains modeling. It was shown by numerical simulations that significant inaccuracies in creep-fatigue damage predictions, reaching an order of magnitude in deviation of inelastic strains, cannot be excluded when simplified models are used. It is thus concluded that full consideration of creep-fatigue damage at real operating conditions requires proper modeling of stress/ strain histories using visco-elastic-plastic material models.


Introduction
The increasing demand for higher thermal efficiency and higher operational flexibility of modern power generation plants results in severe mechanical and thermal loading conditions experienced by plant components in service [1]. In order to increase the thermal efficiency, it is often necessary to increase the power plant operating temperature, while the higher operational flexibility requires frequent and faster turbine start-ups and shutdowns. The extremely severe service conditions of the plant components are subjected to produce a progressive material deterioration due to the accumulation of creep and fatigue damage and can lead to component failure [2].
Steam turbine components operating at high and highly variable temperatures, e.g., rotors, casings and valve chests, are among the most critical components of power plant units, and their lifetime is already limited at the design stage [3]. Safe and reliable operation of steam turbines can only be ensured with proper monitoring systems, periodic inspections and comprehensive lifetime assessment studies [4][5][6].
Typical operation cycle of a steam turbine consists of a start-up, steady state, load change, shutdown and natural cooling. During transient phases, the low-cycle fatigue damage is generated, while steady-state operation brings about the creep wear [7]. In order to accurately determine the creep-fatigue life exhaustion, inelastic material models and a proper approach to multiple cycles consideration have to be applied [8].
Traditional methods of lifetime assessment basing on the use of representative operation data and elastic material models can result in significant inaccuracies in damage prediction and lead to non-optimal decisions regarding future operation [9]. That is why the use of inelastic material models and long-term operation data is becoming mandatory for reliable calculation of creep-fatigue damage and residual life prediction.
The paper presents the results of investigations performed for a steam turbine rotor showing possible inaccuracies and, in particular, differences in predicted strain accumulation and damage depending on the material model and cycles analysis method used.

Cyclic operation of steam turbines
The lifetime of steam turbine components depends on variable loading conditions resulting in fatigue and creep-fatigue damage. Typical service start-stop cycle includes start-up, full load, part load and shutdown phase, generating variable temperature and stress-strain distributions in turbine components. A schematic illustration of different operation phases and the resulting temperature and stress variations on a rotor surface is shown in Figure 1. During turbine runup, rotor rotational speed increases from turning gear to rated speed (0-1) and steam parameters start rising from the initial values. After unit synchronization (1), turbine loading starts (1-2) with continued increase of steam temperature and pressure, and after some time, all the parameters reach their nominal values. The start-up process is finished and steady-state phase (2-3) with constant process parameters is reached. During steady-state operation, load changes (3-4, 5-6) with different rate can take place, which are accompanied by steam temperature and pressure variations. The operating cycle is closed by the shutdown (7-8) and natural cooling phases when all the process parameters go down.
The corresponding stress variations at a rotor surface are presented in Figure 1b. During startup, the rotor surface temperature increases with a higher rate than the center temperature, which results in compressive stresses at the surface [10]. The stress attains maximum and then slowly decays to a steady-state stress. The following load changes generate transient thermal stresses of different sign, but usually the highest tensile stress is generated on the rotor surface during the shutdown phase. The major stress range Δσ is thus formed by the start-up and shutdown stresses, while the minor stress range is defined by the load change stresses.
Full consideration of creep-fatigue damage at real operating conditions requires modeling of long stress/strain histories using visco-elastic-plastic material models. Such a modeling approach is impractical when 2D or 3D models are used to represent the real component geometry and the balance equations solved by means of a finite element method. For the purpose of numerical investigation of creep-fatigue damage development in steam turbine rotors, simplified cycles are defined with different material models prescribed to transient and steady-state phases. It was done in order to investigate the effect of material model applied and number of cycles analyzed on the creep-fatigue damage accumulation. The four cycle types are schematically shown in Figure 2. They are as follows: 1. Full multiple cycle: fatigue simulated with elastic-plastic model and creep with viscous model for multiple start-stop cycles.

2.
Simplified multiple viscous cycle 1: fatigue simulated with elastic-plastic model for one cycle and creep with viscous model for multiple start-stop cycles.   The first cycle describes the most realistic representation of actual operating conditions, while the other three cycles represent simplified approaches in terms of material models and number of cycles, which are used in engineering practice.

Problem formulation
The problem under consideration is cyclic heating up and cooling down of a steam turbine rotor taking place during start-up and the subsequent shutdown. Rotor heating and cooling are caused by hot steam flow through the turbine with varying temperature, pressure and velocity. Non-stationary heat transfer takes place via forced convection, and the thermal load is a primary load of the rotor [11]. The rotor rotates with a constant rotational speed ω in steadystate operation. Steam pressure and rotational body forces due to rotation are also considered in the model. Non-uniform temperature distribution in the rotor induces thermal stresses due to thermal expansion.
The following assumptions are adopted in the model: • The rotor material is isotropic and its physical and mechanical properties are temperature dependent.
• Heating and cooling process is described by the linear theory of heat conduction and nonlinear convective boundary conditions.
• The resulting thermal stresses are determined from the solution of system of equations of uncoupled thermoelasticity and plasticity.
• Due to the geometrical and load symmetry of the rotor about its axis, the computational region is assumed axisymmetric.
• Mechanical loads are modeled as surface pressure and volumetric body force.
Based on the above assumptions, a boundary problem of heat conduction and thermoelasticity is formulated and solved by means of a finite element method [12].

Heat conduction problem
Heat conduction in a homogeneous isotropic solid is described by the Fourier-Kirchhoff differential equation [13] div kgradT where T is the metal temperature, k is the thermal conductivity, g is the heat source, ρ is the density and c p is the specific heat.
Neglecting the heat source g and assuming heat conduction in the radial and axial direction only, Eq. (1) is rewritten in cylindrical coordinate system as follows: A non-uniform temperature distribution T r; z ð Þ is assumed as initial condition at t ¼ 0. For rotor free surfaces, the boundary conditions are: where α z; t ð Þ is the heat transfer coefficient varying with time and axial co-ordiante, T s is the metal surface temperature and T f z; t ð Þis the fluid temperature depending on time and axial coordiante. The outer radius of the rotor r z ð Þ varies with the axial coordinate z defining its outer contour.
The material physical properties are temperature dependent and their variation is described by polynomial functions. The variation of heat transfer coefficient in time and space α z; t ð Þ is considered by the Nusselt number Nu changing in time as a function of the Reynolds number Re and Prandtl number Pr: The Nusselt number is defined as follows: where d denotes a characteristic diameter. A detailed form of Eq. (6) depends on the surface type and flow character [11]. The heat transfer model adopted in this study is based on the well-proven formulae for heat transfer coefficients in steam turbine rotors and can be employed in online calculations of temperatures and thermal stresses [14], which is not possible when more advanced thermal FSI (fluid-structure interaction) modeling is adopted [15][16][17][18].
Geometrical model of the rotor is shown in Figure 3.

Thermoelasticity problem
Knowing the axisymmetric temperature field in the rotor, a solution of the problem of stress state induced by it can be obtained. The equilibrium conditions written for stresses are transformed, taking into account the relations between stresses, strains and displacements, to obtain differential equations with unknown displacements u and w in the r and z direction, respectively [19]: where the following definitions were employed: In Eq. (7), ν is the Poisson's ratio and β is the thermal expansion coefficient.
The system of partial differential equations (7) is solved with the following boundary conditions: • At the outer cylindrical surfaces of the discs σ r ¼ p b z ð Þ-pressure due to centrifugal forces of blades • At the outer surfaces in contact with steam σ n ¼ Àp z ð Þ-steam pressure acting normal to the surface • At the left end face The relations between strains and displacements for an axisymmetric body in the cylindrical coordinate system are expressed as [19]: Stress components are obtained from the solution of the constitutive equation of linear thermoelasticity [13] and for an axisymmetric body in the cylindrical coordinate system they are given as [19]: Þ is the shear modulus and E is the Young's modulus.
Eq. (9) is known as the Duhamel-Neumann relations and has a fundamental importance in thermal stress analysis within the validity of Hooke's law [20].

Plasticity model
At stress concentration areas, thermal stresses induced by nonuniform temperature field can exceed the material yield stress and give rise to plastic deformation. Instead of estimating the plastic strains based on thermoelastic stresses and using analytical notch stress-strain correction methods [4,14], a more accurate approach with direct use of elastic-plastic material was adopted. Numerical calculations were performed using Abaqus [21], employing elastic-plastic material model with the Prager-Ziegler linear kinematic hardening. The plasticity surface is defined by the Huber-Mises-Hencky function [22]: where f σ ij À α ij À Á is the equivalent stress related to the backstress α ij defining translation of the yield surface. The yield function is traditionally defined as follows where s ij is a deviatoric part of the stress tensor σ ij , while α d ij is a deviatoric part of the backstress tensor α ij .
The linear kinematic hardening model assumes the associated plastic flow rule in the form [23]: where _ ε p ij is a plastic flow rate and _ λ denotes here a plastic work. In the linear kinematic hardening model, translation of the yield surface is described by the backstress tensor α ij whose evolution in time is determined by the Prager-Ziegler linear hardening law where μ is a positive scalar coefficient.

Creep model
For creep strain predictions, the characteristic strain model proposed by Bolton was employed [24,25]. The model is simple and effective in description of creep deformation at long times. Based on the analyses of creep test data, Bolton extended the classical Norton model [26] and proposed the following isochronous relation between the stress exponent and stress expressed as a fraction of creep rupture strength σ R [24]: (14) where σ R denotes the creep rupture strength for a given time at constant temperature. Integrating the above relationship, he obtained an equation for the creep strain ε C in the form: where ε χ is a characteristic creep strain, which is a material constant at a given time and temperature. The characteristic creep strain can be evaluated using Eq. (15) and by knowing two stresses, i.e. the creep rupture strength σ R1 at time t 1 and the stress σ D1 to produce datum creep strain ε D at time t 1 . With this assumption, the isochronous stress-strain relation of Eq. (15) can be written as: To close the model, a relationship for the rupture strength described by a simple power-law relationship is adopted: where m denotes the exponent in the power-law relationship. The exponent m can be evaluated from two values of rupture strengths at time t 1 and t 2 : Finally, combining Eqs. (16) and (17), the model relationship between creep strain and time at a constant stress assumes the form: The creep model is thus described by three constants: σ R1 -creep rupture strength at time t 1 ; σ R2 -creep rupture strength at time t 2 and σ D1 -stress of creep strain ε D at time t 1 , which clearly have physical significance and can be derived from readily available data from standard creep tests.

Creep-fatigue damage modeling
Creep damage is assessed on the basis of accumulated creep strain and is calculated using the ductility exhaustion rule [27,28]. For a given load cycle i of duration t h , creep damage d cr i is calculated from the formula [29] where ε R is the creep ductility depending on the creep strain rate. The total creep damage D cr due to n cycles is obtained by summing up the fractional creep damage values: Fatigue damage is assessed on the basis of accumulated plastic strain. For a given load cycle i, fatigue damage d f i is calculated from the formula [30]: where ε f is the plastic strain at failure. The total fatigue damage D f due to n cycles is obtained by summing up the fractional fatigue damage values: According to the linear damage accumulation rule, the total creep-fatigue damage D is a sum of that due to creep and fatigue [31]: The total damage D cannot exceed a critical value D c above which creep-fatigue crack initiation can be expected. The rules of calculating the critical damage are given in different national standards [32][33][34][35].
Creep-fatigue damage calculation using the linear damage accumulation rule, creep and plasticity models and the proposed cycle descriptions can be used in lifetime calculations of both new and old turbines that have been in service for a longer period of time. For new machines, the following input data are typically required in order to perform lifetime calculations: • Number of cold, warm and hot starts that the turbine has to perform.
• Requested operation time with given steam parameters and loads.
• Variation in time of the basic process parameters like steam temperature and pressure, rotor rotational speed and turbine load for cold, warm and hot starts and shutdown (these data are provided in the turbine design start-up diagrams).
Creep-fatigue damage calculations are performed for the cycles defined using the above data.
For old machines, lifetime calculations are carried out in order to determine the current damage and predict the residual life assuming unchanged operating conditions. The accuracy of lifetime calculations strongly depends on the range and quality of operation data provided by the turbine user. In this case, long-term data characterizing turbine operation in the whole operation time or at least representative operation period should be used for computations. The range of data required for such analyses is similar to that required for new machines with the difference that the data were recorded and characterize real operation instead of the ideal operation assumed at the design phase. Creep-fatigue damage calculations are performed for the cycles defined using the real operation data.
The difference between design and real conditions can be high and result in a huge deviation of damage and residual life predictions exceeding an order of magnitude.

Damage development in a steam turbine rotor
According to the strain-based damage accumulation rules, damage due to a load cycle depends on the inelastic strain accumulated in the cycle and the strain at failure. For a given material and loading conditions, the failure strain can be assumed as constant, thus the accumulated strain can be adopted as a direct measure of creep or fatigue damage. This property was utilized in the present study to examine creep-fatigue damage development in a steam turbine rotor. In particular, the effect of prior creep and plastic deformation on the creep and plastic strain evolution is investigated, and its impact on damage development is discussed in detail.

Thermal and elastic analysis
Temperature and elastic stress distributions in the examined rotor were obtained by solving the heat conduction and thermoelasticity equations given in Sections 3.1 and 3.2. Numerical solutions were obtained by means of the finite element method [12] using Abaqus code [21]. Transient and steady-state temperature distributions in the rotor are shown in Figure 4. During turbine start-up, highest temperature gradients, predominantly in the radial direction, are present in the rotor hottest region close to the steam inlet (Figure 4a). In steady-state, temperature distribution shown in Figure 4b, mainly axial temperature gradients are present which result from steam expansion in the turbine steam path and gland system. Transient temperature gradients produce high thermal stresses in the rotor leading to thermal fatigue cracking in stress concentration areas. High temperature in the rotor hottest sections in combination with centrifugal stresses lead to the creep damage which can also be localized in the stress concentration zones [36,37]. In the examined rotor, the area of highest stress and temperature is located in the gland section where several heat relief grooves are present ( Figure 5). High elastic stresses are developed at the bottom of the grooves and they frequently exceed the material yield stress. This results in local plastic deformation whose repetitive occurrence leads to thermal fatigue cracking often found in the heat grooves [38]. The heat relief grooves are thus the most critical locations of the rotor and their creep-fatigue damage development is analyzed in detail in the subsequent sections.

Effect of creep deformation on plastic strain accumulation
Cyclic operation of the turbine and the resulting plastic strain accumulation in the heat grooves were simulated by three start-stop cycles, each consisting of cold start-up, steadystate operation and shutdown followed by natural cooling to the subsequent cold start. Throughout the entire analysis, the elastic-plastic material model was enabled with plasticity model defined in Section 3.3. Two types of cycles were utilized in order to examine the effect of creep deformation on plastic strain accumulation, namely full multiple cycle ( Figure 2a) and simplified multiple elastic-plastic cycle (Figure 2d). The former accounts for creep during steady-state operation, whereas the latter cycle neglects the viscous effects and stress relaxation under constant load. The duration of one full cycle was 386 h, which results in the total time to complete three cycles equal 1158 h. The comparison of equivalent plastic strain distribution after each cycle for both types of analysis is shown in Figure 6. The zone of plastic deformation increases in size during cycling and its shape remains almost unchanged being close to semi-elliptical. Both models, i.e. elasto-plastic and visco-elasto-plastic, predict similar development of the plastic zone but with slightly higher equivalent plastic strains at the bottom of the groove observed in the elastic-plastic model predictions (without creep effects). The evolution in time of the maximum equivalent plastic strain at the groove is shown in Figure 7 together with the corresponding temperature cycles. The equivalent plastic strain gradually increases over time and so does the difference between the two models. It is thus found that creep deformation decelerates the plastic strain accumulation, which can be explained by the stress relaxation taking place during steady-state phase. As it is seen from Figure 7b, the first deviation in equivalent strains occurs after c.a. 110 h, which corresponds to   the beginning of shutdown phase that started with different stress distributions in the groove computed with and without creep. Creep relaxation has thus a positive effect on the fatigue damage measured by the plastic strain accumulation. By neglecting the viscous effects in the steady-state operation, the plastic strains are overpredicted, and consequently the calculated fatigue damage will be too high as compared with the results of the visco-elastic-plastic model. Figure 8 presents the evolution of plastic strain components in the groove bottom over 3 cycles. The strain state is very close to plane strain, with the radial and axial components having the same value and opposite sign, and the circumferential plastic strain being close to zero. It is interesting to see that initially the radial and axial strain components have different signs during operation and natural cooling but after the second cycle the radial strain becomes positive and the axial component remains negative.

Effect of plastic deformation on creep strain accumulation
Similar to the effect of creep deformation on plastic strain accumulation, the effect of plastic deformation on creep strain accumulation can be expected. In order to investigate this phenomenon, different simplified cycles were adopted for analysis, namely simplified multiple viscous cycle 1 (Figure 2b) and simplified multiple viscous cycle 2 (Figure 2c). The duration of one full cycle was the same as in previous analyses, i.e. 386 h, but creep calculations with viscous model enabled were performed for steady-state phase lasting 100 h. This resulted in the total time of creep strain accumulation in 3 cycles equal 300 h.
The comparison of equivalent creep strain distribution after each cycle for all types of analysis is shown in Figure 9. Similar to the behavior of the plastic zone, the zone of creep deformation increases in size during cycling and its shape remains almost unchanged. But this time, the applied models predict visibly different development rate of the creep zone. The largest zone size and equivalent creep strain are obtained for the full multiple cycle (Figure 9a), which takes into account both viscous and plastic effects in each individual start-stop cycle (Figure 2a). A smaller creep zone size and equivalent creep strain are obtained for the simplified multiple viscous cycle 1 (Figure 9b), which includes viscous effects in each individual start-stop cycle but neglects plastic strain development in cycles 2 and 3 ( Figure 2b). The creep strain zone is nearly invisible in case of the simplified multiple viscous cycle 2, which completely neglects plastic effects.
The differences in creep strain accumulation are better seen from Figure 10a, which shows development of the maximum equivalent creep strain at the groove bottom. After the first start-stop cycle, a higher difference in the creep strain reaching 150% is found between the elastic and elastic-plastic models (cycles with and without plasticity). After the second cycle, the deviation increases, but a significant difference in creep strain reaching 100% occurs between the elastic-plastic models, i.e. full multiple cycle and simplified multiple viscous cycle 1. This is due to neglecting the plastic strain development in cycles 2 and 3. It follows from the evolution of the equivalent creep strains that the secondary creep regime is reached already after 100 h (the first cycle) for all simplified modeling approaches, while for the full multiple cycle the heat grooves are still in the primary creep regime after 300 h (third cycle). This qualitative difference is very important from the point of view of creep damage, which on one hand depends on the accumulated creep strain and on the other hand is a function of the failure strain (creep ductility) primarily depending on the creep strain rate, which is very different in the three creep regimes.
The evolution in time of the creep strain components at the groove is shown in Figure 10b. The creep strain state at the bottom of the groove is three dimensional with all three components having comparable values. The radial and circumferential strain components are negative, and the axial component is positive and has the largest absolute value. In contrast to the plastic strain components, all the creep strain components increase continuously with continuous slope within the cycle and discontinuous slope change between two subsequent cycles.

Summary
Creep-fatigue damage development in a steam turbine rotor in cyclic operation was numerically investigated using different material models and cycle modeling approaches. Creep strain was adopted as a creep damage parameter, and plastic strain was adopted as a damage parameter for fatigue. Four different operation cycles were defined and used in numerical simulations to examine the effect of various simplifications on the creep fatigue-damage prediction. The analyzed multiple cycles included different combinations of material models and number of individual cycles used in damage calculations during steady-state and transient conditions. Based on the performed investigations, the following observations can be formulated: • Prior creep deformation affects the plastic strain accumulation by reducing the plastic strain as compared with no creep conditions. The difference in the accumulated plastic strains increases over subsequent cycles.
• Local plastic deformation in heat grooves accelerates the creep strain accumulation during the cycle and over subsequent cycles.
• Neglecting the creep overestimates the fatigue damage (plastic strain) generated during cycling, and neglecting the plasticity in cyclic duty underestimates the creep damage (creep strain) in stationary conditions.
It was shown by numerical simulations that significant inaccuracies in creep-fatigue damage predictions, reaching an order of magnitude in deviation of inelastic strains already after first few cycles, cannot be excluded when simplified models are used. It is thus concluded that full consideration of creep-fatigue damage at real operating conditions requires appropriate modeling of stress/strain histories using visco-elastic-plastic material models. Lifetime assessment studies employing such methodologies ensure the highest possible accuracy of life predictions and optimal decisions regarding units' future operation.

Mariusz Banaszkiewicz
Address all correspondence to: mbanaszkiewicz@imp.gda.pl The Szewalski Institute of Fluid-Flow Machinery Polish Academy of Sciences, Gdańsk, Poland