Effect of Flexible Operation on Residual Life of High-Temperature Components of Power Plants

Abstract

Electricity generation from renewable energy sources is emerging as a result of global carbon emission reduction policies. However, most renewable energy sources are non-dispatchable and cannot be adjusted to meet the fluctuating electricity demands of society. A flexible operation process has been proposed as an effective solution to compensate for the unstable nature of renewable energy sources. Thermal load fluctuations during flexible operation may cause creep–fatigue damage to the high-temperature components of thermal power plants, as they are designed with a focus on creep damage under a constant power level. This study investigated the residual life of high-temperature components, such as a superheater tube and a reheater header, to failure under flexible operation conditions using finite element analysis and empirical models. First, we determined an analytical solution for the straightened superheater tube under thermal conditions and compared it with the numerical solution to verify the numerical models. Through the verified finite element model, the creep–fatigue life of the reheater header was estimated by considering flexible operation factors and employing the Coffin–Manson and Larson–Miller models. Although fatigue damage increases with decreasing minimum load and ramp rate, we confirmed that creep damage significantly affects the residual life during flexible operation. In addition, a surrogate model was proposed to evaluate the residual life of the reheater as a function of the flexible operation factors using the machine learning methodology, based on the results of finite element methods. It can be used to predict its residual life without performing complex thermo-structural analysis and relying on empirical models for fatigue and creep life. We expect our findings to contribute to the efficient operation of thermal power plants by optimizing the flexible operation factors.

1. Introduction

Global CO₂ emissions have steadily increased over time, and Belbute [1] has suggested that the emissions will increase by 27.4% by 2030. Prior studies have addressed the adverse effects of fossil fuels on CO₂ emissions and the need to generate electricity using renewable energy sources (RES) for low-carbon growth [2,3]. Most RES have non-dispatchable characteristics that prevent them from adjusting the energy supply to meet the demand. Climatic or geographic conditions cause variability in electricity production and can render power plant systems more unstable. According to the average net electricity production per week in Germany in 2021, there was a power fluctuation at approximately 12:00 p.m. owing to the influence of solar power [4]. Τhe concern about power outages caused by over- and under-generation increases with the proportion of RES in the power generation system. As a result, it is critical to address the power uncertainty caused by the increased demand for RES [5]. Flexible operation is a method of changing the output level of a thermal power plant to maintain constant output power according to the change in the output power of the RES. This is the most economical way to handle the instability of RES energy generation [6].

Flexibility indicates the efficiency of a concept. Although higher flexibility addresses the system stability problem, it can be accompanied by thermal load fluctuations that place conventional power plants in harsh conditions. The thermo-structural behavior of these facilities under thermal load variation must be considered as conventional power plants were designed by considering the fatigue and creep of the thermal parts under constant power levels. This study entails an investigation of the thermo-structural behavior of high-temperature components in conventional power plants by using computational analysis, bypassing the time efficiency weakness of practical experiments.

Previous studies primarily focused on assessing the life of power plants. In particular, researchers investigated the unsteady temperature and pressure changes in the superheater header during operation. They concluded that shutdown would cause significant damage to the plant components [7,8]. These studies have helped resolve the structural stability problems caused by creep–fatigue damage during power plant operation. However, these analyses were performed at a general load, and it was impossible to predict the change in the residual life of the plant component under flexible operating conditions. Several researchers [9,10,11,12,13] have studied the cumulative damage of each boiler component through a numerical approach. In addition, prior research focused on the thermo-structural fatigue life under specific loading conditions. Still, it did not evaluate the fatigue and creep life of high-temperature components under fluctuating thermal loading conditions. They described the risk of failure by focusing on creep–fatigue damage accumulation in the boiler header part and performed a quantitative evaluation of the residual life under specific boundary conditions. Although they performed a quantitative assessment of the residual life for different boundary conditions, previous works disregarded the variation in critical life to failure due to changes in thermo-structural conditions.

Recently, studies on flexible operations have progressed. Avagianos [14] summarized several prior studies that determined the allowable limits of operational flexibility features using commercial analysis programs and in-house codes. They argued that existing thermal power plants that do not account for flexible operation in their design might be vulnerable to thermal fatigue in high-flexibility situations. However, individual studies have paid little attention to the quantitative life evaluation of structural vulnerabilities and the prediction of the variation in residual life for the operation process. Several studies have investigated methods that can quickly respond to fuel combustion and load fluctuations during flexible operation [15,16]. They suggested that steady-state operation can yield higher efficiency in fuel usage and CO₂ emissions and that rapid load alternation reduces this efficiency. However, little attention has been paid to estimating residual life with fluctuating thermal loads by flexible operation. Therefore, a life evaluation study considering fatigue and creep life is required to investigate the reliability of high-temperature components because of thermal load fluctuations under flexible operating conditions.

The flowchart for studying the remaining life assessment of high-temperature pipes under flexible driving conditions is as follows. First, we calculated the stresses and strains using the thermo-structural simulation program ANSYS Mechanical. The boundary conditions of FEA were determined by the flexible operation cycle. Second, by utilizing the solutions (i.e., strain and stress range), we numerically evaluated the fatigue and creep damages through the Coffin–Manson equation [17] and the Larson–Miller method [18]. We considered the creep–fatigue interaction, which reduces the service time of the structure. Finally, a response surface model using a machine learning method was developed, which performs life assessment without repetitive experiments or numerical analysis. Our study aimed to estimate the residual life to failure under flexible operating conditions by performing finite element analysis (FEA). First, we calculated the strain range in the high-temperature pipe using thermo-structural analysis. We determined an analytical solution for high-temperature pipes under thermal conditions (i.e., convection, conduction, and radiation) and compared it to the numerical solution to verify our FEA platform. In addition, we evaluated the fatigue life using the Coffin–Manson equation [17] and the life by creep damage through the Larson–Miller method [18] under high pressure and temperature. Considering the creep–fatigue interaction, we calculated the overall damage and residual life of the power plant components. We also analyzed the response to the fatigue life for each variable (i.e., minimum load changes and ramp rate). As a result, we confirmed that creep had a greater effect on the overall life than fatigue, and flexibility significantly affected the creep–fatigue damage and residual life. Furthermore, this study introduced a response surface model using a machine learning method, which performs life assessment without repetitive experiments or numerical analysis.

2. Methodology

2.1. Thermo-Structural Analysis for Evaluation of Strain Range during Flexible Operation

In flexible operation, conventional power from fossil fuels is adjusted to generate power based on daily consumption. Flexibility is a factor that determines how the plant quickly controls the load in response to the changing demand for electricity. To increase the proportion of renewable energy, conventional power plants are required to control the load freely. However, high flexibility (low minimum load, fast ramp rate) causes rapid load changes, resulting in thermal damage to facilities with extreme load variations [19,20]. Therefore, this study investigated the relationship between the flexibility and thermo-mechanical behavior of high-temperature components in power plants.

We considered two main factors: ramp rate and minimum load. The ramp rate is the loading rate during flexible operation. A plant with a high ramp rate can rapidly respond to operating conditions. The minimum load is the lowest power that can be stably applied. The generation characteristics of the plant are more flexible with a lower minimum load. It is challenging to thoroughly analyze an actual situation because complex load fluctuations caused by natural phenomena should be considered. This focused on the tendency of the components to fail, thus simplifying the analysis. Despite the deviation from reality, this model proved sufficient for identifying the failure tendency. Figure 1 shows the ideal flexible operation cycle, divided into the following four sections.

The analysis was conducted for the flexible operation cycle presented in Figure 1, and the temperature and pressure of the working fluid varied according to the power generation output that changed from (a) to (d). This high-temperature steam applies a thermo-structural load to the header components and pipes. The FEA was performed with these components under ideal flexible operating conditions. The structure of the entire power plant was too complex and large to perform numerical analysis. As a result, we substituted the model with one column of header and piping by sub-modeling [9,21]. Plant facilities are subjected to loads by high-temperature fluids. Therefore, they are influenced by the hydrodynamics of their flows. However, the flow patterns of the fluids were not important in this study because we aimed to evaluate the structural behavior of high-temperature components under thermal conditions. Thus, we assumed that the flow field around the analysis model was sufficiently large and that the local temperature and heat transfer coefficient were equal to the empirically obtained average values [22,23].

Figure 2 presents a flowchart of the analysis. First, we derived the boundary conditions for the heat and internal pressure of each structural element based on heat transfer theories and solid mechanics. Equation (1) expresses the differential heat conduction in cylindrical coordinates and Equation (2) expresses the triaxial strain formula with the thermal effects.

$$\frac{1}{r}\frac{\partial }{\partial r}(r\cdot k\frac{\partial T}{\partial r})+\frac{1}{r^2}\frac{\partial }{\partial \varphi }(k\frac{\partial T}{\partial \varphi })+\frac{\partial }{\partial z}(k\frac{\partial T}{\partial z})+{\dot{e}}_{gen}=\rho c\frac{\partial T}{\partial t}$$

(1)

$$\begin{array}{l}{\epsilon }_t=\frac{1}{E}\{{\sigma }_t-\nu ({\sigma }_r+{\sigma }_z)\}+\alpha T\\ {\epsilon }_r=\frac{1}{E}\{{\sigma }_r-\nu ({\sigma }_t+{\sigma }_z)\}+\alpha T\\ {\epsilon }_z=\frac{1}{E}\{{\sigma }_z-\nu ({\sigma }_t+{\sigma }_r)\}+\alpha T\end{array}$$

(2)

where $\rho ,\hspace{0.17em}c,\hspace{0.17em}T,\hspace{0.17em}E,\hspace{0.17em}{\dot{e}}_{gen}$ corresponds to the density, specific heat, temperature, elastic modulus, and heat generation, respectively. ${\epsilon }_t,\hspace{0.17em}{\epsilon }_r,\hspace{0.17em}{\epsilon }_z$ indicate the hoop, radial, and axial strains, respectively. Finally, ${\sigma }_t,\hspace{0.17em}{\sigma }_r,\hspace{0.17em}{\sigma }_z$ represent the hoop, radial, and axial stresses.

Then, we calculated the stresses using the thermo-structural simulation program ANSYS Mechanical. The coupling analysis resolved the strain range problem corresponding to each condition. Additionally, utilizing the solution, we numerically evaluated the creep and fatigue damages using MATLAB (See Section 3.1 and Section 3.2). Finally, we considered the creep–fatigue interaction, which reduces the service time of the structure. Our model applied the condition in the creep–fatigue envelope corresponding to the ASME boiler section. The residual life of an object was evaluated by determining whether the structure was broken when it exceeded the threshold.

2.2. Creep–Fatigue Damage Theory

Thermo-structural loads can damage power plant components in the form of creep or fatigue. The fatigue life of the components is primarily evaluated using the Coffin–Manson model and low-cycle fatigue life equation [17].

$$\frac{∆\epsilon }{2}=\frac{{{\sigma }^{\prime }}_f}{E}{(2N_f)}^b+{{\epsilon }^{\prime }}_f{(2N_f)}^c$$

(3)

In Equation (3), $∆\epsilon$ is the strain range and $N_f$ is the critical cycle for fatigue failure. ${\sigma }_f{}^{\prime }$= 807 MPa and $b$= − 0.1486 are the coefficients and exponents for elastic deformation, and ${{\epsilon }^{\prime }}_f$= 0.1125 and $c$= − 0.4355 the ones for plastic deformation, respectively. These are material constants derived from fatigue test results [24]. Additionally, in the high-stress range, the alternating load generates a high mean stress ${\sigma }_m$, and derives a mean stress effect that further reduces the fatigue life [25,26]. In this paper, the fatigue life was predicted using Morrow’s equation (Equation (4)), considering the mean stress effect in the strain-life approach [27].

$$\frac{∆\epsilon }{2}=\frac{{{\sigma }^{\prime }}_f}{E}(1-\frac{{\sigma }_m}{{{\sigma }^{\prime }}_f}){(2N_f)}^b+{{\epsilon }^{\prime }}_f{(2N_f)}^c$$

(4)

Palmgren [28] and Miner [29] defined pure fatigue damage as the sum of the ratio of the progressive cycle $N$ to the failure cycle $N_f$, as expressed by Equation (5). The fatigue damage is a time-independent variable because $N_f$ follows the Morrow equation and depends only on the load change.

$$D_f={\displaystyle \sum _i\frac{N}{N_f}}$$

(5)

The creep deformation should also be considered in the design of thermal power plants, as they are exposed to high temperatures for an extended period while operating. The creep rupture time depended on the stress and temperature applied to the model, and the relationship could be determined experimentally. The National Institute of Material Science (NIMS) conducted creep tests at various stresses and temperatures on the SUS304 steel [30]. Several studies have employed the Larson–Miller method to predict creep life by defining the Larson–Miller parameter (LMP) using Equation (6) [18,31]. According to the experimental data, the stress was linearly proportional to LMP, as expressed in Equation (6). The coefficients $A$= −64.3 and $B$= 1474 are obtained by solving a least-squares problem.

$$\begin{array}{l}{\sigma }_r=A\times (LMP)+B\\ where\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}LMP=T\cdot {\log }_{10}{t}_r+20\end{array}$$

(6)

Similar to fatigue damage, creep damage can be calculated using the cumulative damage model. Robinson defined creep damage as the ratio of the loading time $t$ to the rupture time $t_r$ at a specific stress and temperature, as expressed in Equation (7) [32]. In other words, the rupture time can be predicted from the stress and temperature conditions using Equation (6), and the creep damage was evaluated using the following formula:

$$D_c={\displaystyle \sum _i\frac{t}{t_r}}$$

(7)

Equation (8) summarizes the creep–fatigue damage interaction. When the cumulative damage $D$ exceeded the permissible damage value $D_a$, the material was considered to have failed [33]. The simplest rule for evaluating creep–fatigue damage is the linear damage summation rule, which sets $D_a$ to 1 [18]. The $D_a$ value can be defined empirically, depending on the load type and material properties [34].

$$\begin{array}{l}D=D_f+D_c\\ D_f+D_c<D_a\end{array}$$

(8)

A conservative safety factor should be considered during the design stage to avoid unintentional power plant destruction. Therefore, the damage model for objects should be applied conservatively, which involves a higher safety factor. Figure 3 shows the ASME standard design of SUS304 obtained through experiments, which has a joint at (0.3, 0.3). In the bilinear region of ASME Section III for boilers and pressure vessels, $D_c$ and $D_f$ follow the two relational expressions in Equation (9) [35]. The critical cycle for safe operation under flexible conditions was estimated by solving these two inequalities.

$$\left\{\begin{array}{l}{D}_c\le -\frac{7}{3}{D}_f+1\hspace{1em}for D_f\le 0.3\\ D_c\le -\frac{3}{7}(D_f-1) for D_f\ge 0.3\end{array}\right.$$

(9)

2.3. Machine Learning Techniques

We derived the residual life from a specified set of flexibility factors through experiments or computational analyses. However, obtaining all outputs for the entire flexibility range using only the aforementioned methods is time consuming and inefficient. Therefore, this study entailed the development of a life assessment method that employs machine learning to predict the life for operating conditions that have not been previously analyzed.

2.3.1. Feedforward Neural Network Model

The feedforward neural network is an artificial neural network (ANN) model that propagates information from the input layer to the output layer. In this network, each hidden layer uses a rectified linear unit (ReLU) function, $f(x)\triangleq \max (0,x)$, as an activation function to provide the learning data. The $j^{th}$ node value of the $k^{th}$ hidden layer was given by Equation (10), and the model was trained to reduce the mean-squared error between the data and actual output. The upper subscript $k$ corresponds to the label number of the layer, the lower subscript $j$ to the label of the recipient node, and $i$ to the label of the giver node. Let $n_{layer}$ be the number of hidden layers; then, $a_i^0$ is the input value of the $i^{th}$ node, and $a_i^{n_{layer}+1}$ is the output value of the $i^{th}$ node.

$$a_j^k=f(b_j^k+{\displaystyle \sum _{i=0}^{j_{k-1}}{w}_{ji}^k{a}_i^{k-1}})$$

(10)

The ANN model in Figure 4 is trained using a back-propagation algorithm called Bayesian Regularization. The letter updates the weights $w_{ji}^k$ and bias $b_j^k$ using the partial differentiation of the objective function $V(x)$. The algorithm computes $∆x=-{[J^T(x)J(x)+\mu I]}^{-1}{J}^T(x)e(x)$ to update the $V(x)$ to $V(x+∆x)$. This process continues until $V(x)$ reaches the optimal value [36]. In the Bayesian Regularization algorithm, $V(x)$ consists of a linear combination of the sum of square errors ${\displaystyle \sum e{(x)}^2}$ and the sum of square weights ${\displaystyle \sum w{(x)}^2}$, where $a_{\alpha }$ and $a_{\beta }$ are adjusted through Bayesian optimization [37]:

$$V_{BR}(x)=a_{\alpha }{\displaystyle \sum e{(x)}^2+a_{\beta }{\displaystyle \sum w{(x)}^2{V}_{BR}}}$$

(11.a)

$${\displaystyle \sum e{(x)}^2}={\displaystyle \sum {(t(x)-a^{n_{layer}+1}{}_i)}^2}$$

(11.b)

2.3.2. Hyperparameter Optimization Using Random Search

To improve the performance of the machine learning model, the hyperparameters (e.g., learning rate and number of hidden layers) should be appropriately selected through an empirical method, such as a grid search or random search. The grid search algorithm discretizes the range of hyperparameter values into a finite set, or grid of values and trains the ANN model multiple times. The number of training trials is $N={\displaystyle \prod _k{n}_k}$, where $n_k$ refers to the number of grid values of the $k^{th}$ hyperparameter. A grid search performs well in low-dimensional problems, but its exhaustive strategy necessitates a large number of trials when there are numerous hyperparameters. For example, in a test including hyperparameters with 10 discrete values, a grid search will require ${10}^4$ trials to identify the best performing value.

The random search algorithm overcomes this difficulty. It generates random combinations of values from a given hyperparameter range and selects the combination set with the best performance. When tuning multiple hyperparameters, it is more likely to find the optimal value than the grid search. Furthermore, more time is reserved when determining the optimal parameters to use for training the machine learning model [38]. Numerous studies have demonstrated the algorithm’s accuracy through statistical and mathematical approaches [38,39,40]. Figure 5 shows a schematic comparison of the two search methods. Under the same conditions, the random search method offers better accuracy if the performance responses are too complex or sensitive to specific hyperparameters.

3. Numerical Examples

3.1. Validation for Thermo-Structural FE Model

We used the ANSYS Workbench software 2022 R2 for thermo-structural analysis. Numerical results were compared to the analytical solution to validate the analytical method. Heat transfer from the combustion gas and work fluid occurs during the boiler operation. The combustion gas interacts with the outer walls of the pipe components through natural convection and radiation. In addition, the working fluid flowing inside the pipe transfers heat to the inner wall through forced convection. The header part, located outside the combustion chamber, was not affected by the combustion gas and was under standard air conditions. The overall boundary conditions of the header and tube components are shown in Figure 6. Power plants are vulnerable to creep and fatigue failure. To this end we used Creep Strength Enhanced Ferritic Steel (CSEF) in this environment, which exhibits good heat resistance, and high creep strength. In this study, we used SUS304 steel, a common type of CSEF steel.

Table 1. Composition of the SUS304 steel.

Elements	C	Mn	Si	Cr	Ni	N	Nb	P
Composition (wt%)	0.07–0.13	1.00	0.010	17.0–19.0	7.5–10.5	0.05–0.12	0.30–0.60	0.040

lists the composition of the SUS304 steel. The temperature conditions change as the load fluctuates, owing to flexible operation. Therefore, the material properties of SUS304 are temperature dependent, and the temperature variation due to the operating conditions affects the properties.

Table 2. Temperature-dependent properties of the SUS304 steel.

Material Property	Value at 300 °C	Value at 500 °C	Value at 700 °C
Density, ρ $(kg/m^3$)	7790	7700	7610
Thermal expansion coefficient, α $(\times {10}^{-6}/°C)$	9.7	10.05	10.3
Elastic modulus E $\left(GPa\right)$	164.78	148.93	132.38
Poisson’s ratio ν	0.2874	0.2946	0.3018
Thermal conductivity, k $(W/m °C)$	21.461	24.923	30.98
Specific heat, c $(J/kg °C)$	542.62	579.71	616.81

and

Table 3. Thermo-structural boundary conditions during the operation of the power plant.

Boundary Condition			Under 100% Condition	Under 30% Condition
Tube	$\mathrm{Steam} \mathrm{temperature}, T_{\infty ,\hspace{0.17em}in}$ $(°C)$	Inlet	501.71	474.22
		Outlet	502.29	475.32
	$\mathrm{Flue} \mathrm{gas} \mathrm{temperature}, T_{\infty ,\hspace{0.17em}ex}$ $(°C)$	Inlet	1057.89	1066.17
		Outlet	843.19	857.56
	Internal pressure p $\left(MPa\right)$	Inlet	25.303	9.787
		Outlet	25.298	9.786
	$\mathrm{Convective} \mathrm{film} \mathrm{coefficient}, h_{conv }$ ($W/m^2$)	$\mathrm{Flue} \mathrm{gas}, h_{conv,\hspace{0.17em}ex}$	8.582	8.216
		$\mathrm{Steam}, h_{conv,\hspace{0.17em}in}$	5436.24	1620.95
	Emissivity	$\mathrm{Tube}, {\epsilon }_{tube}$	0.8
		$\mathrm{Gas}, {\epsilon }_{gas}$	0.281
Header	$\mathrm{Steam} \mathrm{temperature}, T_{\infty ,\hspace{0.17em}in}$$(°C)$		596	572
	Convective film coefficient, ($W/m^2$)		2403.8	980.54
	Internal pressure, p ($MPa$)		4.599	1.483

list the material properties and thermo-structural boundary conditions with the load level. The values between these categories were interpolated linearly.

Equation (12.a) defines the heat transfer coefficient ($h$) as the product of the thermal conductivity $k$ and the Nusselt number ($Nu$). $Nu$ is related to the Reynolds ($\mathrm{Re}=\rho VD/\mu$), Rayleigh ($\mathrm{Ra}=\frac{g\alpha }{\nu \phi }(T_s-T_{\infty })D^3$), and Prandtl numbers ($\mathrm{Pr}=\lambda /\phi$), where $\mu$ is the viscosity, $\phi$ is the thermal diffusivity, $\lambda$ is the kinematic viscosity, $V$ is the flow velocity, $D$ is the diameter of the pipe, and $g$ is the acceleration due to gravity. Equation (12.b) shows the empirical formulas used to define $Nu$ [23,41]. Given the value of $Nu$, the heat transfer coefficient can be calculated from Equation (12.a) [41].

$$h=Nu\times k$$

(12.a)

$$\mathrm{Nu}=\left\{\begin{array}{l}0.023{\mathrm{Re}}^{0.5}{\mathrm{Pr}}^{0.4} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for Forced convection\\ {\left[0.6+0.387{\left\{\frac{\mathrm{Ra}}{{\left(1+{\left(\frac{0.559}{\mathrm{Pr}}\right)}^{9/16}\right)}^{16/9}}\right\}}^{1/6}\right]}^{2}for Natural convection\end{array}\right.$$

(12.b)

Stefan–Boltzmann’s law was used to calculate the radiant heat transfer coefficient $h_{rad}$ of the combustion gas (see Equation (13)), where $\sigma$ is the Stefan—Boltzmann constant, ${\epsilon }_{tot}$ is the total emissivity, $T_g$ is the gas temperature, and $T_w$ is the wall temperature adjacent to the gas [42]. The composition and properties of the combustion gas mixture are described in the literature [43,44].

$$\begin{array}{c}q=h_{rad}(T_g-T_w)=\sigma \cdot {\epsilon }_{tot}(T_g{}^4-T_w{}^4)\\ where \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{h}_{rad}=\sigma \cdot {\epsilon }_{tot}\frac{T_g{}^4-T_w{}^4}{T_g-T_w},\hspace{0.33em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_{tot}=\frac{(1+{\epsilon }_{tube})}{2}{\epsilon }_{gas}\end{array}$$

(13)

Equation (1) provides a solution for the heat conduction equation. It is difficult to verify this model for header components with complex geometries. Therefore, this study validated the model by comparing these equations with the numerical results of a simple straight-pipe model. The boundary conditions used in the analysis are shown in Figure 6. This can be summarized as follows:

(a)

Radiant and convective heat transfer from combustion gas at the outer wall;

(b)

Conduction in the tube wall;

(c)

Convection at the inner wall of the working fluid.

The tube temperature changed from (a) to (b), and thermal deformation occurred. The heat flux during the heat transfer process is explained in the literature [23]. As each heat flux $\dot{q}$ is equal throughout the entire section, the heat flux at any time step follows Equation (14), which is a quaternary equation for heat flux. Therefore, we used Newton’s method to solve the complicated equation and derive its approximate solution, which was compared to the numerical solution calculated from the FE model. The validation was performed in the load-descending section, as shown in Figure 1b.

$$\dot{q}=\frac{T_{\infty .ex}-T_{\infty .in}+{\epsilon }_{tot}\sigma \{T_{\infty .ex}-{(\dot{q}(R_{cond}-R_{cond.in})+T_{\infty .in})}^4\}\times R_{conv.ex}}{R_{cond}+R_{conv.in}+R_{conv.ex}}$$

(14)

$T_{\infty .ex}$ and $T_{\infty ,in}$ are the external and internal temperatures of the tube, respectively, $R_{cond}$ is the thermal conductivity $k$ divided by the thickness $D$, $R_{conv,ex}$ and $R_{conv,in}$ are the reciprocal of the convective film coefficients $h_{ex}$ and $h_{in}$, which follow Equations (12) and (13), respectively.

Figure 7 compares the analytical solutions using Equation (14) and numerical results. The right-hand side of Equation (1) can be solved by utilizing the approximate equality between the heat flux and temperature over time. Using a finite difference method, changes in temperature with location can also be computed, in agreement with the laws of physics [23]. Therefore, the thermal boundary condition of the numerical model is well posed.

Temperature variations cause thermal deformation and stress in the object. Their mechanical behavior follows Equation (2). The triaxial thermal stresses of the superheater tube in a cylindrical coordinate system follow Equation (15) [22].

$$\begin{array}{l}{\sigma }_r=\frac{E\alpha }{1-\nu }\frac{1}{r^2}\left[\frac{r^2-r_i{}^2}{r_o{}^2-r_i{}^2}{\displaystyle {\int }_{r_i}^{r_o}Tr dr-{\displaystyle {\int }_{r_i}^{r}Tr dr}}\right]\hspace{1em}\hspace{1em}\hspace{1em}Radial stress\\ {\sigma }_t=\frac{E\alpha }{1-\nu }\frac{1}{r^2}\left[\frac{r^2+r_i{}^2}{r_o{}^2-r_i{}^2}{\displaystyle {\int }_{r_i}^{r_o}Tr dr+{\displaystyle {\int }_{r_i}^{r}Tr dr-T{r}^2}}\right]\hspace{1em}Hoop stress\\ {\sigma }_z=\frac{E\alpha }{1-\nu }\frac{1}{r^2}\left[\frac{2}{r_o{}^2-r_i{}^2}{\displaystyle {\int }_{r_i}^{r_o}Tr dr-T}\right]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}Axial stress\end{array}$$

(15)

From the numerical model, the thermal stress was estimated and compared with the analytical solution, as shown in Figure 8. Figure 8a compares the thermal results, confirming that the heat transfer and thermal load conditions in the numerical model follow the governing equation. In addition, the power plant equipment is subjected to internal pressure. Dowling [45] proposed stresses in thick-walled pressure vessels by internal pressure $p$(see Equation (16)).

$$\begin{array}{l}{\sigma }_r=-\frac{p{r}_i{}^2}{r_o{}^2-r_i{}^2}\left(\frac{r_o{}^2}{r^2}-1\right)\hspace{1em}\hspace{1em}\hspace{1em} Radial stress\\ {\sigma }_t=\frac{p{r}_i{}^2}{r_o{}^2-r_i{}^2}\left(\frac{r_o{}^2}{r^2}+1\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} Hoop stress\\ {\sigma }_z=0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}Axial stress\end{array}$$

(16)

The comparison of stress in the superheater tube under the structural load by pressure is presented in Figure 8b. Equations (15) and (16) can be linearly combined using the superposition method, as the stress generated by the heat and pressure is within the elastic range. Finally, Figure 8c shows the results of the coupled thermo-structural stresses by thermal loading and internal pressure. We validated the boundary conditions of the numerical model by comparing the results to those of the analytical solution.

3.2. Estimation of Fatigue Life under Cyclic Thermal Loads

In this analysis, all power plant components were divided into several submodels for time efficiency. Figure 9 shows the numerical models for each part. The finite element model of the components consist of SOLID 185 3D elements, which are linear hexahedral elements used for thermo-structural analysis.

Table 4. The number of nodes and the elements of the boiler components.

Part	Num. of Nodes	Num. of Elements
Tube	31,648	24,021
Header	24,324	17,527

lists the number of elements and nodes [46].

The thermo-structural boundary conditions are changed, since the factors of flexible operation modify the power level of the power plant. Each feature depends on the time $t$ and the power generation output $\varphi (t)$ at any given moment. Any boundary condition $X$ (e.g., the temperature of the working fluid and internal pressure at that time) follows Equation (17), where $∆X$ and $∆\varphi$ represent the changes in the boundary condition and power generation over time in a specific operation section, respectively. Additionally, $X_o$ and ${\varphi }_o$ indicate the initial values.

$$X(t)=\frac{∆X}{∆\varphi }\left\{\varphi (t)-{\varphi }_o\right\}+X_o$$

(17)

The electricity power level $\varphi (t)$ affects the condition $X$. Furthermore, it varied depending on the position in the boiler tube. Therefore, they are path dependent, and $X$ follows the binary function of Equation (18). The distance $d$ is $d_{in}$ at the inlet and $d_{out}$ at the outlet of the path.

$$X(d,t)=\frac{X(d_{out},t)-X(d_{in},t)}{d_{out}-d_{in}}\left\{d-d_{in}\right\}+X(d_{in},t)$$

(18)

Figure 10 shows the numerical contours of the boiler parts, where the boundary conditions of Equation (18) were applied under a minimum load of 20% and a ramp rate of 3%/min.

Table 5. Analysis results at minimum load ($P_{\min }$ ) of 20% and ramp rate of $3\hspace{0.17em}\%/\min$.

Value	At Header Component	At Tube Component
Maximum temperature over time (°C)	595.52 °	527.70 °
Minimum temperature over time (°C)	566.20 °	513.34 °
Maximum von Mises stress over time (MPa)	228.09	113.61
Minimum von Mises stress over time (MPa)	13.169	74.761
Maximum strain range	$14.648\times {10}^{-4}$	$5.175\times {10}^{-4}$

summarizes the stress and temperature results of the reheater header and superheater tube. The creep–fatigue damage in Equation (6) depends on the data listed in Table 5. The maximum stress in the header part was 1.92 times higher than that in the tube part. Therefore, we expected the header component to be more damaged by the thermo-structural load than the other parts. Based on this result, we determined that this part is vulnerable to creep–fatigue damage and selected it for assessing the residual life.

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

3.3. Creep and Fatigue Life of the Header

The creep–fatigue damage of the header in Equation (7) depends on the change in flexibility during the operation. We assumed that the power plant had an ideal flexible operation cycle, as shown in Figure 1, and that the holding time of the minimum load section was $t_c=3h$. Under these assumptions, the operating time in each section is defined as $t_a=24-(t_b+t_c+t_d), t_b=t_d=\frac{P_{\max }-P_{\min }}{\nu }, t_c=3$, where $P_{\min }$ is the minimum load and $\nu =dP/dt$ is the ramp rate The creep–fatigue damage was estimated under operating conditions using the data presented in Table 5. Fatigue damage occurred only in sections (b) and (d) of the optimized cycle shown in Figure 1, where the load alternated and caused strain. In contrast, creep damage, which occurs based on stress and temperature levels, appeared throughout the operation cycle. Equation (19) gives the total creep damage $D_c$ and fatigue damage $D_f$, that the header receives after operating $N$ cycles, with one cycle being $24h$. We can rearrange Equation (7) to solve for the operation cycle $N$, to obtain Equation (20). If $N$ reaches the critical cycle $N_{cr}$, this inequality satisfies the equal sign, and the damaged object is considered to have failed.

$$D=D_f+D_c=N\times (d_f+d_c)$$

(19)

$$N\le N_{cr}=\left\{\begin{array}{l}\frac{3}{3d_c+7d_f}{ \mathrm{for} \mathrm{D}}_f\le 0.3\\ \frac{3}{7d_c+3d_f}{ \mathrm{for} \mathrm{D}}_f\ge 0.3\end{array}\right.$$

(20)

The boundary conditions change when the minimum load changes and affect creep and fatigue damage. The variation in the ramp rate did not affect the conditions, but it affected the operating time, as shown in Figure 1. Since fatigue damage is time independent, it is unaffected by the ramp rate. In contrast, the creep damage varies with the flexibility factor. We estimated the critical cycle for the given creep–fatigue damage using Equation (20), Figure 11a shows the fatigue damage as a function of ramp rate and minimum load change. The results indicate that the fatigue damage was only affected by the variation in the minimum load. The fatigue damage decreased as the minimum load increased and was unaffected by the ramp rate changes. In other words, it is independent of the ramp rate. In contrast, creep damage depends on two flexibility factors (i.e., ramp rate and minimum load), as shown in Figure 11b. We demonstrated that the creep damage was higher for increasing ramp rate and minimum load. Additionally, the creep damage change for the minimum load over a certain ramp rate was not significant. This results from the fact that sections (b) and (d) in Figure 1, which cause fatigue damage, are shortened at high ramp rates, and creep damage becomes dominant.

The creep damage per cycle $d_c$ is approximately 100 times higher than the fatigue damage $d_f$ in most operating conditions. As a result, the expectancy of the critical cycle $N_{cr}$ is dominated by creep effects (Figure 12). The boiler parts generally have lower residual life under operating conditions that cause high creep damage levels (i.e., higher minimum load and ramp rate). Figure 13 shows the damage to the plant header during 3000 cycles in the creep–fatigue envelope, suggesting that the plant element is vulnerable to creep under most flexibility cases. During the flexible operation of thermal power plants, engineers must prioritize creep damage over fatigue damage.

Although the fatigue damage is minor compared to the creep damage, the effect of fatigue damage may increase as the ramp rate and minimum load decrease. As a result, fatigue damage must be considered when evaluating the residual life. Furthermore, this study assumed that only two load fluctuations (sections (b) and (d) of Figure 1) occurred per day, considering an ideal flexible operation. However, the proportion of fatigue to residual life is expected to increase as more complex operations progress. For these reasons, the fatigue effect should also be considered when assessing the residual life of boiler components.

3.4. Response Surface Model

Τo investigate the residual life, we estimated the parameters (e.g., strain range, temperature, and stress level). We proposed an ANN regression model that predicts the residual life under flexible operating conditions based on the above results. The proposed model follows the neural network structure shown in Figure 4. To improve its accuracy, we selected three hyperparameters for adjustment: the number of hidden layers ($n_{layer}$), the number of nodes per layer ($n_{nodes}$), and the initial learning rate (${\zeta }_{inital}$).

Table 6. Range of the hyperparameters to optimize.

Value	${\mathit{n}}_{\mathit{l}\mathit{a}\mathit{y}\mathit{e}\mathit{r}}$	${\mathit{n}}_{\mathit{n}\mathit{o}\mathit{d}\mathit{e}\mathit{s}}$	${\mathit{\zeta }}_{\mathit{i}\mathit{n}\mathit{i}\mathit{t}\mathit{a}\mathit{l}}$
Max	2	50	0.1
Min	1	5	0.000001

lists the ranges of the selected parameters. Twenty parameter sets were randomly extracted using the random search method and the model was trained to find the optimal set of hyperparameters.

To evaluate the performance of the neural network model, it is essential to normalize the input and output data to share the same scale. We used min–max normalization, which places the smallest value of the data at zero and the largest at one. The other values are scaled between them. Our ANN model sets each flexibility (i.e., ramp rate and minimum load) as an input value and the residual life $N_{cr}$ as an output value. In the training process using Equations (10) and (11), the machine learning network divides the data into the training set (60%), test set (20%), and validation set (20%) to use for cross-validation. For every iteration $i$, the network is trained using the training set and evaluated using the test set. The training performance was determined by the mean square error (MSE) $e_{test}(i)$ of the test set. The weights and biases were updated accordingly to reduce the error. The error $e_{val}(i)$ of the trained model was also evaluated on the validation set to avoid overfitting the test set. Our study identified the optimizer with the lowest validation error $e_{val}(i_{best})$ at the epoch $i_{best}$ and the model performance as the test error $e_{test}(i_{best})$ at this epoch. Subsequently, the ANN model was trained for each hyperparameter set. Among them, $(n_{layer},\hspace{0.17em}{n}_{nodes},\hspace{0.17em}{\zeta }_{inital})=(2,\hspace{0.17em}28,\hspace{0.17em}0.0931)$ exhibited the best performance, $e_{test}(i_{best})=2.1798\times {10}^{-5}$, and was selected as the optimal set. Figure 14 presents the error and coefficient of correlation $R$ for each dataset. Through this process, we defined a well-constructed ANN model exhibiting a small error for most inputs. Figure 15 shows the optimized structure of the ANN model. We propose an ANN model that can effectively evaluate the residual fatigue life within a given minimum load and ramp rate range.

4. Conclusions

This study aimed to assess the residual life during flexible operation using numerical methods and empirical models. We considered two main factors of flexible operation, ramp rate and minimum load, and investigated their effects on residual life. FEA was performed under the given temperature and pressure conditions to estimate the strain range of the superheater tube and reheater header. The finite element results were verified by comparing them to the analytical solution of the straightened tube under thermal load conditions. We estimated the creep–fatigue life of the reheater header under flexible operation factors using finite element modeling and the Coffin–Manson and Larson–Miller models. We numerically demonstrated that fatigue damage occurred owing to thermal cycling and increased with decreasing minimum load and ramp rate. Nevertheless, creep damage was the dominant damage mechanism. This originates from the fact that the thermal cycle is a low cycle that repeats up to twice a day, and creep is the dominant mechanism in terms of overall damage. In addition, by utilizing the ANN model, we proposed a response surface model for evaluating the residual life of the reheater heater, which is the most vulnerable component under flexible operation. This model can predict the residual life of the reheater header according to the flexible operation factors without performing complex thermo-structural analysis and empirical models for fatigue and creep life.

This is the first systematic study to analyze the effect of the flexible operation on the residual life of the tube and header. To evaluate the residual life under flexible operation conditions using experimental methods, there are some difficulties in terms of time and research execution. In this study, to overcome these difficulties, computational mechanics and empirical equations were used as an approach. It additionally presents a model that can predict the residual life using the response surface model. We believe that our findings will aid the efficient operation of thermal power plants by optimizing flexible operation factors in the future. Furthermore, it is expected that this study will contribute to the cost reduction and determination of maintenance cycles in power plant maintenance under flexible driving conditions.