This chapter presents two practical methods of thermoelastic stress calculation suitable for application in online monitoring systems of steam turbines. Both methods are based on the Green function and Duhamel integral and consider the effect of variable heat transfer coefficient and material physical properties on thermal stresses. This effect is taken into account either by using an equivalent steam temperature determined with a constant heat transfer coefficient or by applying an equivalent Green’s function determined with variable heat transfer coefficient and physical properties. The effectiveness of both methods was shown by comparing their predictions with the results of exact three-dimensional (3D) calculations of a steam turbine valve.
- thermoelastic stress
- steam turbine
- online calculation
The increasing demand for a higher thermal efficiency and a higher operational flexibility of modern power generation plants results in severe mechanical and thermal loading. Modern energy markets put a requirement of a high operational flexibility of power plant units . The flexible operation generates elevated loads and stresses in turbine components leading to material damage due to thermo-mechanical fatigue. In order to control the fatigue damage, thermal stresses should be computed online and limited to permissible values .
In power generation industry, the approach based on the Duhamel integral and the Green functions has been widely used for thermal stress calculations [3–6]. Such an approach has also been adopted for monitoring power boiler operation  and can also be used for calculating stress intensity factors at transient thermal loads .
The major issue in using Green’s function and Duhamel’s integral method in modeling transient thermal stresses in steam turbine components is the time dependence of material physical properties and heat transfer coefficients affecting proper evaluation of Green’s functions and making the problem nonlinear. There are known approaches assuming the determination of the influence functions at constant values of these quantities . The inclusion of temperature-dependent physical properties proposed by Koo et al.  relies on determining the weight functions for steady-state and transient operating conditions. A more important variation of the heat transfer coefficient can be taken into account by calculating the surface temperature using a reduced heat transfer model and employing Green’s function to calculate a stress response to the step change of metal surface temperature . However, numerical solution of a multi-dimensional heat transfer model for complex geometries present in steam turbines is complicated and time-consuming, and due to this, it cannot be used in online calculations. A full inclusion of time variability of the physical properties and heat transfer coefficients proposed by Zhang et al.  relies on the solution of nonlinear heat conduction problem by using artificial parameter method and superposition rule and replacing the time-dependent heat transfer coefficient with a constant value together with a modified fluid temperature. The effectiveness of the method has been proved by an example of a three-dimensional model of a pressure vessel of nuclear reactor.
This chapter presents and compares two alternative methods for steam turbine components both employing the Green functions and Duhamel’s integral and considering the variations of material physical properties and heat transfer coefficient. The first method uses a numerical solution of one-dimensional (1D) heat conduction problem and determination of equivalent steam temperature for use with the Green function determined with a constant heat transfer coefficient. The second method employs the equivalent Green function determined with variable physical properties and heat transfer coefficient. The methods have been validated and compared for a steam turbine valve with the three-dimensional (3D) numerical solutions.
2. Thermoelasticity problem
The problem under consideration is heating up and cooling down of a steam turbine taking place during start-up and the subsequent shutdown. Heating and cooling of steam turbine components is caused by a 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 turbine. The rotor rotates with a constant rotational speed ω in a steady-state operation. Steam pressure and rotational body forces due to rotation are also considered in the model. Non-uniform temperature distribution in the component induces thermal stresses due to thermal expansion.
The material is assumed isotropic, and its physical and mechanical properties are temperature-dependent. The 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. Due to the geometrical and loading symmetry of the turbine components like rotors or valve casings, the computational region is assumed axisymmetric. Mechanical loads are modeled as surface pressure and volumetric body force.
Knowing the axisymmetric temperature field in the component, 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 radial r- and axial z-direction, respectively ,
where T is the temperature, T0 is the initial temperature, ν is Poisson’s ratio, β is the thermal expansion coefficient, ρ is the density and G is the shear modulus. In Eq. (1), the following definitions were used:
The system of partial differential equations (1) is solved with the following boundary conditions:
At the outer cylindrical surfaces of the rotor
At the outer surfaces of the rotor in contact with steam
At the inner surfaces of the casing in contact with steam
At the end face of the rotor or casing
The relations between strains and displacements for an axisymmetric body in the cylindrical coordinate system are expressed as 
Stress components are obtained from the solution of the constitutive equation of linear thermoelasticity , and for an axisymmetric body in the cylindrical coordinate system, they are given as 
is shear modulus and E is the Young modulus.
3. Green’s function method
It is assumed that the thermal stress model used in online calculations is based on Green’s function and Duhamel’s integral method. This method allows for fast calculation of thermal stresses at supervised areas for any changes of fluid temperature causing heating up or cooling down of an element. In steam turbine components, a reliable and accurate measurement of the surface temperature at the critical region is usually not possible; thus, for online monitoring purposes, the steam temperature measurement is used as a leading signal. Consequently, in heat transfer model, Fourier’s boundary condition is employed :
where λ(r) is the metal thermal conductivity, T(r, t) is the surface temperature, Tst(t) denotes steam temperature, n is the normal direction and α(r) is the heat transfer coefficient.
where X(r, t) is Green’s function for temperature.
Using Green’s function X(r, t), we can calculate the metal temperature for arbitrary variation of steam temperature employing Duhamel’s theorem :
Considering elastic stresses, it can be assumed, according to the thermoelasticity theory, that stress distribution in an elastic body is a unique function of temperature distribution, and in this connection, the thermal stresses can be calculated using Duhamel’s integral as
where Gij(r, t) is a Green function for thermal stress tensor component σij, ij = r,θ,z.
For simple geometries like a plate, cylinder or sphere, Green’s functions for temperature and stress can be determined analytically by solving one-dimensional transient heat conduction problem. In case of more complicated shapes, it is necessary to use numerical methods, for example, finite element method in order to obtain accurate Green’s functions.
In order to solve Duhamel’s integral given by Eq. (11), it is transformed from a continuous into a discrete form:
where td is a cutoff time, Gs,ij(r) is a value of the influence function in steady state, while Gt,ij(r) is a transient part of the influence function.
Green’s functions are determined individually for each supervised region, and in case of thermal stresses, these functions must be determined for each component of the stress tensor σij.
An example of Green’s function for temperature and stress components at a step increase of fluid temperature by 1°C heating the external surface of a long cylinder is shown in Figure 1. The calculations were performed with a constant heat transfer coefficient α = 1000 W/m2/K and constant material physical properties evaluated at temperature 350°C. The temperature (Figure 1a) and stress (Figure 1b) responses are presented for the cylinder axis and its outer surface heated by steam.
Green’s function for surface temperature X_surf initially rapidly increases and reaches nearly 1 after approx. 104 s. The temperature at cylinder axis X_cent shows some inertia, and a delay time of temperature response equal to 500 s is seen. After this time, the temperature increases at approximately constant rate, and after c.a. 4000 s, the rate of variation visibly starts to decrease. The cylinder surface response is typical for surface-type sensor, while Green’s function at cylinder axis is a response of middle-type sensor. The stress response at both areas is a result of temperature variations: axial and circumferential stresses at cylinder surface are equal to each other and rapidly grow reaching a minimum (compressive stress) and subsequently tend slowly to zero; the stresses at cylinder axis grow more slowly, their maxima are lower than at surface (axial stress is two times higher than circumferential stress) and are tensile stresses (sign +). All stress components within the cylinder tend to zero and a stress-free state is reached after approx. 1.4 × 104 s, which corresponds to the instant when temperature within the whole cylinder is uniform.
4. Green’s function for variable physical properties and heat transfer coefficient
The above-presented Green functions were determined with constant physical properties and constant heat transfer coefficient. The physical properties of steel influencing temperature (specific heat , thermal conductivity ) and stress (Young’s modulus , Poisson’s number , thermal expansion coefficient ) distribution strongly depend on temperature, while the heat transfer coefficient varies with temperature and flow velocity (via the Reynolds number Re). Among the considered physical properties, the most varying is the specific heat which, for a typical low-alloy heat-resistant steel, increases by more than 50% in the temperature range of 20–500°C. The largest drop is observed for the Young modulus which decreases by more than 30% in the same temperature range. Heat transfer coefficient depends mostly on the flow velocity and during turbine start-up can change even by two orders of magnitude.
A consequence of the described variations of physical properties and heat transfer coefficient is a dependence of the Green function on these parameters. The dependence is illustrated in Figure 2 showing the Green function variations for axial stress at cylinder surface calculated by the finite element method at different temperatures and heat transfer coefficients. Each individual curve corresponds to stress variation determined at constant temperature and constant heat transfer coefficient . As it is seen from the plots, the largest effect on the Green function is observed for the heat transfer coefficient, which significantly affects the stress maxima, time of their occurrence and decay rate. The influence of temperature variation is much smaller and most visible from the instant of maximum stress to the moment of their vanishing.
5. Practical methods for thermoelastic stress calculation
In order to consider the variation of the Green function in Eq. (11), two methods are proposed:
Use of an equivalent steam temperature determined with a constant heat transfer coefficient
Use of an equivalent Green’s function determined with variable heat transfer coefficient and physical properties
5.1 Equivalent steam temperature
A flow chart of the proposed stress calculation algorithm using an equivalent steam temperature is presented in Figure 3 . The calculation process is split into two steps: thermal analysis (step 1) and stress analysis (step 2). In step 1, thermal calculations are performed online, and one-dimensional (1D) heat conduction problem with time-dependent boundary conditions and temperature-dependent physical properties is solved using the finite difference method (FDM). Boundary conditions are given by the measured steam temperature and heat transfer coefficient computed on the basis of the steam mass flow rate. By solving the 1D heat conduction equation, one obtains the radial temperature distribution in the component model. The use of FDM enables full consideration of nonlinearities introduced by time-dependent boundary conditions and temperature-dependent physical properties. Knowing the thermal boundary conditions and radial temperature distribution , the equivalent steam temperature for the constant heat transfer coefficient is computed from the surface heat flux (see Eq. (22)). Calculations of the equivalent steam temperature are performed online at each time step for the same constant heat transfer coefficient .
In step 2, thermal stress calculations are performed online with the help of the Duhamel integral. Input for the integral is given by the equivalent steam temperature from step 1 and the Green functions determined offline by means of two/three dimensional (2D/3D) finite element method (FEM) calculations carried out for the real component geometry and constant heat transfer coefficient .
Heat conduction in a homogeneous isotropic solid is described by the Fourier-Kirchhoff differential equation 
where is the metal temperature, is the heat source and is the specific heat.
Neglecting the heat source and assuming heat conduction in the radial direction only, Eq. (13) is re-written in a cylindrical coordinate system as follows:
A uniform temperature distribution is assumed as initial condition at . Boundary conditions are defined depending on the component modeled. For rotor sections modeled as an externally heated solid cylinder, the boundary conditions are.
where is the metal surface temperature.
For a hollow cylinder model, condition (16) directly applies to the outer surface, and for the inner surface, the boundary condition assumes the form.
In case of turbine and valve casings, their thermal behavior is approximated by an internally heated hollow cylinder with the following boundary conditions:
The one-dimensional heat conduction Eq. (14) is numerically solved using the finite difference method . The material physical properties are temperature-dependent, and their variation is described by polynomial functions. The variation of heat transfer coefficient is considered by the Nusselt number Nu changing in time as a function of the Reynolds number Re and Prandt number Pr:
The Nusselt number is defined as follows:
where denotes the cylinder characteristic diameter. A detailed form of Eq. (20) depends on the component geometry and flow character . Practical methods of determining heat transfer coefficients are presented in . The heat transfer model adopted in this study is based on the well-proven formulae for heat transfer coefficients in steam turbine components and can be employed in online calculations of temperatures and thermal stresses which are not possible when more advanced thermal FSI (fluid–structure interaction) modeling is adopted [25–28].
The convective surface heat flux defined by Eq. (16) or (18) with time-dependent heat transfer coefficient can be alternatively calculated using an equivalent steam temperature and a constant heat transfer coefficient :
This relationship can be used to determine the equivalent steam temperature which together with the assumed constant value of heat transfer coefficient ensures the same surface heat flux. The equivalent steam temperature evaluated as
is transferred to step 2 (Figure 3) for stress calculations.
The algorithm basing on the equivalent steam temperature has a capability of modeling different heat transfer mechanisms. The conventional application of the Duhamel integral is based on predefined Green functions evaluated for a constant heat transfer coefficient assuming convection heat transfer mode, which is not the only heat exchange mechanism occurring in steam turbines operation. When the turbine starts from a cold state, steam condensation can take place on the rotor outer surfaces or casing inner surfaces, and condensation heat transfer should be taken into account. It is assumed that condensation occurs on the wall, if the following conditions are satisfied :
The steam temperature is below the saturation temperature at the local pressureE24
The heat flux from the water film to metal is greater than the heat flux from steam to the water filmE25
Condensation heat transfer takes place with very high heat transfer coefficients resulting in high heat fluxes on the component surface. This, in turn, brings about high surface temperature rates producing elevated thermal stresses on the surface or in the stress concentration areas. The most favorable conditions for the condensation to take place prevail at the initial phase of cold start-ups when the metal temperature is low, that is, below 100°C, and the steam saturation temperature high enough to produce high heat fluxes. This phenomenon is minimized by pre-warming steam turbine components and supplying steam of possibly highest superheating.
5.2 Equivalent Green’s function
The second method relies on using an equivalent Green’s function in Eq. (11) determined, assuming a variable heat transfer coefficient and material physical properties. It can be done if the variation of process parameters is known in advance and used to determine the shape of the equivalent Green’s function. This is the case in steam turbines operation as steam parameters are defined in design start-up diagrams for specific types of starts and should be followed within some margins in real operation. The method assumes little deviations of real operating conditions from those defined by the start-up diagrams.
The equivalent Green function is evaluated by performing finite element thermal analysis with the following boundary conditions :
Steam temperature step characteristic for a given start-up type.
Variable heat transfer coefficient resulting from the variation in time of the steam mass flow rate characteristic for a given start-up.
The flow chart of the thermoelastic stress calculation algorithm based on the equivalent Green function is shown in Figure 4. The equivalent Green function is pre-computed offline by performing finite element calculations in 2D or 3D model of the real component. The finite element model is nonlinear due to the temperature-dependent material properties and time-dependent heat transfer coefficient. The equivalent Green functions are evaluated for each stress component and for various transient events like a cold, warm and hot start-up and shutdown. These events are characterized by specific variations of the heat transfer coefficients and the range of temperature variations affecting the material physical properties. The appropriate Green function is selected at the beginning of the transient event and used to calculate the thermal stress online with the help of the Duhamel integral.
6. Validation of the methods for turbine valve
First validations of both methods for a simple cylinder heated externally by steam were presented in  for the equivalent steam temperature and in  for the equivalent Green’s function. A more complex example of a steam turbine valve operating at highly variable steam conditions, like temperature, pressure and flow rate, was used in this study to compare and validate the proposed algorithms. Numerical calculations were performed using a FEM model, and the presented algorithms with transient boundary conditions corresponding to a cold start-up are shown in Figure 5. The equivalent steam temperature was obtained for the calculated surface temperature and the assumed heat transfer coefficient W/m2 K−1. The equivalent Green functions were computed assuming the variation of heat transfer coefficient shown in Figure 5. As it is seen, the equivalent steam temperature exceeds the real steam temperature at a time point where the actual heat transfer coefficient rises above the assumed constant value of , which is a consequence of the equal surface heat fluxes assumed.
Figure 6 presents the variation of surface heat flux and surface temperature rate dTsurf/dt during the cold start-up. Both curves exhibit two local minima: the first one after 2340 s and the second one after 3840 s. They correspond to the instants of the change in slope of the heat transfer coefficient curve shown in Figure 5 and will affect the shape of the Green functions and the variation of the equivalent Huber-Mises stress.
The calculations were performed with variable material physical properties depending on metal temperature. Figure 7 presents their variation during the cold start-up at the point of highest thermoelastic stresses located at the bowl inner surface. The largest drop is observed for the Young modulus which changes by around 30% during the cold start, while the largest rise occurs for the specific heat varying by more than 50% at the same time.
Figures 8 and 9 show the temperature and Huber-Mises stress distributions in the valve casing at three instants of the cold start calculated with the FEM model. These stresses were used to validate the proposed algorithms by comparing them with the predictions of both methods. The temperature field in the valve casing is two-dimensional, which results from both the non-regular geometry and non-uniform heat transfer conditions prevailing on the heated surfaces. The most intensive heat transfer takes place in the bowl and diffuser seat area, and the inner surface temperature is highest in these regions. The highest stress occurs at the inner surface of the valve bowl where both the temperature and stress gradients are predominantly present in the radial direction. Thermal stress distribution in the bowl wall (Figure 9) is similar to the typical distribution in a cylindrical or a spherical wall with stress maxima occurring on the outer and inner surfaces and the minimum stress present approximately in the middle of the wall thickness. The FEM model was also used for determining the equivalent Green functions and Green functions for a constant heat transfer coefficient = 500 W/m2/K in the bowl region and variable physical properties corresponding to the temperature range 30–525°C. The Green functions for each stress component at the bowl inner surface for the constant heat transfer coefficient are shown in Figure 10 and for the varying heat transfer coefficient are presented in Figure 11. All stresses are compressive, and except for the radial component, rapidly increase to their maxima and then slowly decay. The largest stresses are generated in axial and circumferential directions (and , respectively), perpendicular to the direction of the highest temperature gradient. These Green functions implicitly include also the effect of 2D temperature distribution resulting from different heat transfer conditions prevailing on different valve surfaces heated by steam. This effect was accounted for by applying to these surfaces the heat transfer coefficient values different than = 500 W/m2/K (in the equivalent steam temperature method) or α = var (in the equivalent Green function method) which was achieved by applying the appropriate multiplication factors to the base heat transfer coefficients. There is a clear qualitative difference observed between the Green functions obtained with the constant and variable heat transfer coefficient. When the applied heat transfer coefficient is constant during the heating process, all the Green functions exhibit one maximum at the very beginning of the process, followed by continuous stress decay (Figure 10). The effect of non-monotonic variation of the surface heat flux is not included here in the Green function but is taken into account in the equivalent steam temperature rate exhibiting two local minima (Figure 6). When the applied heat transfer coefficient is varying during the heating process, all the Green functions exhibit two maxima: the first one at the beginning of the process and resulting from the initial temperature step and the second one corresponding to the increase of the surface heat flux (Figure 11).
The variation of Huber-Mises stress at the bowl obtained from the FEM analysis is compared in Figure 12 with the predictions of both simplified algorithms. Calculations were carried out using the Green functions shown in Figures 10 and 11, and the steam temperatures presented in Figure 5. As it is seen, both methods predict the occurrence of two stress maxima with high accuracy as compared with the FE results. A slightly better agreement is found for the method using the equivalent Green function which underpredicts the first maximum by 1.1% and the second maximum by 2.5%. The instants of both maxima correspond to those obtained from the FE model. The method using the equivalent steam temperature overpredicts both the first and second stress maximum with an error of 2.8 and 3.3%, respectively. A worse agreement with the FE results is also observed for the instants of stress maximum: the first maximum occurs later while the second one sooner as compared with the FE calculation results.
The occurrence of two subsequent stress maxima predicted in the finite element analysis is a result of variable surface heat flux exhibiting in time local maxima and minima. This effect was reproduced by both simplified methods but in a different way. In the equivalent Green function method, this function had two maxima and in combination with the applied steam temperature profile resulted in the two stress maxima. In the equivalent steam temperature method, this effect was obtained as a result of specific variation of the equivalent steam temperature applied together with the Green function of a typical shape.
Also, the stress decay period persisting after 4000 s of the cold start-up is accurately modeled by the two simplified methods. The relative error in this phase is generally below 10%, and initially, both methods underpredict the FE equivalent stress with a tendency to overpredict it at longer times.
The presented methods for online calculation of the thermoelastic stresses in steam turbine components are based on the Green function and Duhamel integral and consider the effect of variable heat transfer coefficient and material physical properties on thermal stresses. This effect is taken into account either by using the equivalent steam temperature determined with a constant heat transfer coefficient or by applying the equivalent Green’s function determined with variable heat transfer coefficient and physical properties.
The first method is based on a numerical solution of 1D heat conduction problem and the Duhamel integral for thermal stress calculation. The equivalent steam temperature is calculated in step 1 and then used as an input to thermal stress model employing the Green function evaluated with a constant heat transfer coefficient. This approach allows for a full consideration of nonlinearities introduced by variable physical properties and heat transfer coefficient in temperature calculations. The effect of variable heat transfer coefficient significantly affecting the thermal stress results was taken into account by introducing the equivalent steam temperature evaluated on the basis of surface heat flux. The minor influence of variable physical properties on the stress response was considered by determining the Green functions from a step change of steam temperature in the range of its possible variation.
The second method relies on using the equivalent Green’s function determined assuming a variable heat transfer coefficient and material physical properties. The equivalent Green function is evaluated by performing finite element thermal analysis for a step change of steam temperature and variable heat transfer coefficient resulting from the variation in time of the steam mass flow rate.
The effectiveness of both methods was shown by comparing their predictions with the results of finite element calculations of a steam turbine valve. The multi-dimensional geometry and non-uniform thermal boundary conditions result in two-dimensional temperature and stress fields calculated using the FEM model. The comparison of equivalent stresses in the casing bowl calculated using the FEM model and the proposed methods revealed their good accuracy. The relative errors of maximum stress prediction are in the range of 1–3% and can be considered as very low taking into account the complex geometry of the valve, nonlinear boundary conditions and a broad range of temperature variation inducing significant variation of material physical properties. The methods can be applied for online calculation of thermoelastic stress utilized in thermal stress control and lifetime monitoring systems of steam turbines.
|cp [J/kg/K]||specific heat|
|d [m]||characteristic diameter|
|e [−]||volumetric strain|
|E [MPa]||Young modulus|
|g [W/m3]||heat source|
|G [MPa]||shear modulus|
|Gij [MPa/°C]||Green’s function tensor|
|Geq [MPa/°C]||equivalent Green’s function|
|Gs [MPa/°C]||steady-state Green’s function|
|Gt [MPa/°C]||transient Green’s function|
|H [°C]||Heaviside’s function|
|n [—]||normal direction|
|Nu [—]||Nusselt number|
|pb [MPa]||pressure due to centrifugal forces of blades|
|Pr [—]||Prandtl number|
|Qsurf [W/m2]||surface heat flux|
|rin [m]||inner radius|
|rout [m]||outer radius|
|Re [—]||Reynolds number|
|td [s]||cutoff time|
|Teq [°C]||equivalent steam temperature|
|Tsat [°C]||saturation temperature|
|Tst [°C]||steam temperature|
|Tsurf [°C]||surface temperature|
|T0 [°C]||initial temperature|
|u [m]||radial displacement|
|w [m]||axial displacement|
|X [°C/°C]||Green’s function for temperature|
|z [m]||axial coordinate|
|α [W/m2/K]||heat transfer coefficient|
|αcond [W/m2/K]||condensation heat transfer coefficient|
|αconv [W/m2/K]||convection heat transfer coefficient|
|αc [W/m2/K]||constant heat transfer coefficient|
|β [1/K]||thermal expansion coefficient|
|γrz [—]||shear strain component|
|εr [—]||radial strain component|
|εφ [—]||circumferential strain component|
|εz [—]||axial strain component|
|λ [W/m/K]||thermal conductivity|
|ν [—]||Poisson’s ratio|
|σn [MPa]||normal stress|
|σr [MPa]||radial stress component|
|σφ [MPa]||circumferential stress component|
|σz [MPa]||axial stress component|
|σrz [MPa]||shear stress component|
|σij [MPa]||stress tensor|
|τ [s]||artificial time|
|ω [rad/s]||angular velocity|