Study of Different Alternatives for Dynamic Simulation of a Steam Generator Using MATLAB

Abstract

This work presents the simulation of a steam generator or water-tube boiler through the implementation in MATLAB^® for a proposed mathematical model. Mass and energy balances for the three main components of the boiler—the drum, the riser and down-comer tubes—are presented. Three alternative solutions to the ordinary differential equation (ODE) were studied, based on Runge–Kutta 4th order method, Heun’s method, and MATLAB function Ode45. The best results were obtained using MATLAB^® function Ode45 based on the Runge–Kutta 4th Order Method. The error was less than 5% for the simulation of the steam pressure in the drum, the total volume of water in the boiler, and the mixture quality in relation to what was reported.

1. Introduction

The world should prepare for the demand for energy to skyrocket in the next 20 years. By 2040, it will rise by 30% with an annual growth rate of 1.4% [1]. This growth will be like adding another China and India to global demand, warns the annual report of the International Energy Agency (IEA). The IEA indicates that the global economy is growing at an average rate of 3.4% per year. Additionally, the IEA estimates that the population will expand from 7.4 billion to 9 billion people by 2040, and there will be an urbanization process that will add the equivalent of a city the size of Shanghai to the world’s urban population every four months. The energy sector will experience profound changes, with new production powers and a shift in energy sources that will supply energy to humanity.

For these reasons, the energy supply on the planet can be considered a major concern for all countries. In addition, most energy production is carried out with fossil fuels, along with the environmental problems associated with their use, such as global warming, air pollution, and acid rain [2].

One of the most feasible strategies to mitigate the impact generated by the use of fossil fuels in power generation is to improve the efficiency of power plants. The efficient operation of steam generators or boilers is indispensable for the correct functioning of this type of plant known as thermal power plants. The boiler is responsible for transforming the chemical energy of a fuel such as coal, oil, gas, or nuclear energy into thermal energy, using this heat to convert water into steam [3].

The steam generator is used in the oil [4], pharmaceutical [5], and thermoelectric industries, among others. It is necessary to know and understand the functioning of these systems. This has been possible through the formulation of mathematical models that represent their behavior with excellent results [6,7,8,9,10,11,12].

Recent studies have used the mathematical model defined by ordinary differential equations (ODEs) [6]. They used this model to simulate the behavior of thermal power plants with heat recovery systems. This is the case of the works that reported the dynamic behavior of key parameters involved in the steam generation process, such as the appropriate rate of water and fuel consumption according to the rate of heat required in the boiler risers [13]. Other studies reported the response of water consumption and steam generation in the boiler, based on the variables of water level, pressure, and mixture quality in the boiler drum [14].

In relation to the works shown in the literature where numerical models generally include functions that have already been implemented, this work aims to show a comparison between conventional numerical models (Runge–Kutta of order 4 and Heun’s method) and the Ode45 function implemented in MATLAB^®(MathWorks Inc., Natick, MA, USA),implementing the mathematical model and validating the results shown in [6] and, in this way, simulating the operation of a steam generator, since this model adequately represents the dynamics of the system and allows a practical analysis of moderate complexity for its development [14]. A methodology is proposed for the implementation of the mathematical model in MATLAB^®, which can be used in real practical applications. The algorithm for the execution of the model is programmed and three options are studied for the solution of the ordinary differential equation (ODE).

The model is programmed in MATLAB^® and validated with the data reported by [6]. Three numerical methods were analyzed to solve the implicit ODE for the implementation of the model. The results obtained were analyzed using the Runge–Kutta 4th order method used by [13], MATLAB^® function Ode45 implemented by [14], and Heun’s method as a solution alternative. The error is determined for the three methods studied with reference to the results reported in [6].

2. Methodology

This work was developed following the methodology proposed in Figure 1. The objective is to present a systematic process for the implementation of a mathematical model that allows the dynamic simulation of the operation of a water-tube boiler. Figure 1 represents the stages in which the project was carried out, with the (*) each programmed numerical method is indicated, as well as each process variable calculated with the application of the studied methods.

2.1. Model Analysis

The study carried out by [6] has become a reference in recent years for carrying out new research on the behavior of steam generators. This work opened the way for other research studies [13,14,15]. The mathematical model proposed by [6] was studied and analyzed. The main characteristic of this model is that, despite its moderate complexity, it can easily be used to represent any drum boiler with high accuracy.

In Figure 2, the main element of a recovery boiler is detailed, Where, for its study, three main components can be highlighted: the drum, which is where the water–steam mixture is located and, in turn, has the input of the feed water ${(\dot{\mathrm{m}}}_{\mathrm{f}})$ and the output of generated steam ${(\dot{\mathrm{m}}}_{\mathrm{s}})$, and two sections of pipeline: one that transports the water known as the down-comer ${(\dot{\mathrm{m}}}_{\mathrm{dc}})$, and the risers through which the water rises and begins to evaporate thanks to a heat input rate $\left(\dot{\mathrm{Q}}\right)$.

Equation (1) represents the overall mass balance in the drum and Equation (2), the global energy balance.

$$\frac{\mathrm{d}}{\mathrm{dt}}\left[{ \mathsf{\rho }}_{\mathrm{w}}{\forall }_{\mathrm{wt}}{+\mathsf{\rho }}_{\mathrm{s}}{\forall }_{\mathrm{st}}\right]{=\dot{\mathrm{m}}}_{ \mathrm{f}}{-\dot{\mathrm{m}}}_{\mathrm{s}}$$

(1)

$${\dot{\mathrm{Q}}+\dot{\mathrm{m}}}_{\mathrm{f}}{\mathrm{h} }_{\mathrm{f}}{ -\dot{\mathrm{m}}}_{\mathrm{s}}{\mathrm{h}}_{\mathrm{s}}=\frac{\mathrm{d}}{\mathrm{dt}}\left[{\mathsf{\rho }}_{\mathrm{s}}{\forall }_{\mathrm{st}}{\mathrm{h}}_{\mathrm{s}}{+\mathsf{\rho }}_{\mathrm{w}}{\forall }_{\mathrm{wt}}{\mathrm{h}}_{\mathrm{w}}-\mathrm{P}{\forall }_{\mathrm{t}}{+\mathrm{m}}_{\mathrm{m}}{\mathrm{C}}_{\mathrm{p}}{\mathrm{T}}_{\mathrm{m}}\right]$$

(2)

By solving the above balances, the ordinary differential Equations (3) and (4) are obtained, which describe the mathematical model in a simple way for a steam generator.

$${\mathrm{e}}_{11}\frac{\mathrm{d}{\forall }_{ \mathrm{wt}}}{\mathrm{dt}}{+\mathrm{e}}_{12}\frac{\mathrm{dP}}{\mathrm{dt}}{=\dot{\mathrm{m}}}_{ \mathrm{f}}{-\dot{\mathrm{m}}}_{\mathrm{s}}$$

(3)

$${\mathrm{e}}_{21}\frac{\mathrm{d}{\forall }_{ \mathrm{wt}}}{\mathrm{dt}}{+\mathrm{e}}_{22}\frac{\mathrm{dP}}{\mathrm{dt}}{=\dot{\mathrm{Q}}+\dot{\mathrm{m}}}_{ \mathrm{f}}{\mathrm{h}}_{\mathrm{f}}{-\dot{\mathrm{m}}}_{\mathrm{s}}{\mathrm{h}}_{\mathrm{s}}$$

(4)

where

$$\begin{array}{l}{\mathrm{e}}_{11}{=\mathsf{\rho }}_{\mathrm{w}}{-\mathsf{\rho }}_{\mathrm{s}}\\ {\mathrm{e}}_{12}={\forall }_{\mathrm{st}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{s}}}{\mathrm{dP}}+{\forall }_{\mathrm{wt}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{w}}}{\mathrm{dP}}\\ {\mathrm{e}}_{21}{=\mathrm{h}}_{\mathrm{w}}{\mathsf{\rho }}_{\mathrm{w}}{-\mathrm{h}}_{\mathrm{s}}{\mathsf{\rho }}_{\mathrm{s}}\\ {\mathrm{e}}_{22}={\forall }_{\mathrm{st}}{(\mathrm{h}}_{\mathrm{s}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{s}}}{\mathrm{dP}}{+\mathsf{\rho }}_{\mathrm{s}}\frac{{\mathrm{dh}}_{\mathrm{s}}}{\mathrm{dP}})+{\forall }_{\mathrm{wt}}{(\mathrm{h}}_{\mathrm{w}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{w}}}{\mathrm{dP}}{+\mathsf{\rho }}_{\mathrm{w}}\frac{{\mathrm{dh}}_{\mathrm{w}}}{\mathrm{dP}})-{\forall }_{\mathrm{t}}{+\mathrm{m}}_{\mathrm{m}}{\mathrm{C}}_{\mathrm{p}}\frac{{\mathrm{dT}}_{\mathrm{m}}}{\mathrm{dP}}\end{array}$$

Equation (1) can be multiplied by ${\mathrm{h}}_{\mathrm{w}}$ and then subtracted from Equation (2). The result is shown in Equation (4) and this process results in Equation (5).

$$\begin{array}{l}\frac{\mathrm{dP}}{\mathrm{dt}}[{\forall }_{\mathrm{st}}{\mathsf{\rho }}_{\mathrm{s}}\frac{{\mathrm{dh}}_{\mathrm{s}}}{\mathrm{dP}}+{\forall }_{\mathrm{st}}{\mathrm{h}}_{\mathrm{c}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{s}}}{\mathrm{dP}}+{\forall }_{\mathrm{wt}}{\mathsf{\rho }}_{\mathrm{w}}\frac{{\mathrm{dh}}_{\mathrm{w}}}{\mathrm{dP}}-{\forall }_{\mathrm{t}}{+\mathrm{m}}_{\mathrm{m}}{\mathrm{C}}_{\mathrm{p}}\frac{{\mathrm{dT}}_{\mathrm{m}}}{\mathrm{dP}}]+\frac{\mathrm{d}{\forall }_{\mathrm{wt}}}{\mathrm{dt}}{(- \mathsf{\rho }}_{\mathrm{s}}{\mathrm{h}}_{\mathrm{c}})=\\ {\dot{\mathrm{Q}}-\dot{\mathrm{m}}}_{\mathrm{f}}{(\mathrm{h}}_{\mathrm{w}}{-\mathrm{h}}_{\mathrm{f}}{)-\dot{\mathrm{m}}}_{\mathrm{s}}{\mathrm{h}}_{\mathrm{c}}\end{array}$$

(5)

Considering that ${\mathrm{h}}_{\mathrm{s}}{-\mathrm{h}}_{\mathrm{w}}$ represents the enthalpy of condensation ${\mathrm{h}}_{\mathrm{c}}$, Equation (5) is reduced to:

$${\mathrm{e}}_1\frac{\mathrm{dP}}{\mathrm{dt}}{=\dot{\mathrm{Q}} -\dot{\mathrm{m}}}_{\mathrm{f}}{(\mathrm{h}}_{\mathrm{w}}{-\mathrm{h}}_{\mathrm{f}}{)-\dot{\mathrm{m}}}_{\mathrm{s}}{\mathrm{h}}_{\mathrm{c}}$$

(6)

The dominant terms in the coefficient ${\mathrm{e}}_1$ in Equation (6) are $\frac{{\mathrm{dh} }_{\mathrm{w}}}{\mathrm{dP}}$ and $\frac{{\mathrm{dT}}_{\mathrm{m}}}{\mathrm{dP}}$. This is because the energy concentrated in the mass of water and metal are the physical phenomena of the system that determine the dynamic behavior of the pressure in the drum. As a result, ${\mathrm{e}}_1$ can be operated as shown in Equation (7).

$${\mathrm{e}}_1={\forall }_{\mathrm{wt}}{ \mathsf{\rho }}_{\mathrm{w}}\frac{{\mathrm{dh}}_{\mathrm{w}}}{\mathrm{dP}}{+\mathrm{m}}_{\mathrm{m}}{\mathrm{C}}_{\mathrm{p}}\frac{{\mathrm{dT}}_{\mathrm{m}}}{\mathrm{dP}}$$

(7)

Riser and Down-Comer Mathematical Model

Equation (8) below allows the determination of the quality of mixture ${\mathsf{\alpha }}_{\mathrm{r}}$ as a result of the mass and energy balance in the riser and down-comer of the steam generator.

$${\mathrm{e}}_{31}\frac{\mathrm{dP}}{\mathrm{dt}}{+\mathrm{e}}_{32}\frac{{\mathrm{d}\mathsf{\alpha }}_{\mathrm{r}}}{\mathrm{dt}}{=\dot{\mathrm{Q}}- \mathsf{\alpha }}_{\mathrm{r}}{\mathrm{h}}_{\mathrm{c}}{\dot{\mathrm{m}}}_{\mathrm{dc}}$$

(8)

where

$$\begin{array}{l}{\mathrm{e}}_{31}=({\mathsf{\rho }}_{\mathrm{w}}\frac{{\mathrm{dh}}_{\mathrm{w}}}{\mathrm{dP}}\hspace{1em}{\mathsf{\alpha }}_{\mathrm{r}}{\mathrm{h}}_{\mathrm{c}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{w}}}{\mathrm{dP}}\left)\right(1-{\overline{\mathsf{\alpha }}}_{\mathrm{v}}){\forall }_{\mathrm{r}}{+\mathrm{m}}_{\mathrm{t}}{\mathrm{C}}_{\mathrm{p}}\frac{{\mathrm{dT}}_{\mathrm{s}}}{\mathrm{dP}} - {\forall }_{\mathrm{r}}+\\ \left(\right((1 - {\mathsf{\alpha }}_{\mathrm{r}}) {\mathrm{h}}_{\mathrm{c}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{s}}}{\mathrm{dP}}{+\mathsf{\rho }}_{\mathrm{s}}\frac{{\mathrm{dh}}_{\mathrm{s}}}{\mathrm{dP}})){\overline{\mathsf{\alpha }}}_{\mathrm{v}}{\forall }_{\mathrm{r}}+((1 - {\mathsf{\alpha }}_{\mathrm{r}}) {\mathsf{\rho }}_{\mathrm{s}}- {\mathsf{\alpha }}_{\mathrm{r}}{\mathsf{\rho }}_{\mathrm{w}}) {\mathrm{h}}_{\mathrm{c}}{\forall }_{\mathrm{r}}\frac{{\mathrm{d}\overline{\mathsf{\alpha }}}_{\mathrm{v}}}{\mathrm{dP}}\end{array}$$

$${\mathrm{e}}_{32}=\left((1-{\mathsf{\alpha }}_{\mathrm{r}}) {\mathsf{\rho }}_{\mathrm{s}} - {\mathsf{\alpha }}_{\mathrm{r}}{\mathsf{\rho }}_{\mathrm{w}}\right){ \mathrm{h}}_{\mathrm{c}}{\forall }_{\mathrm{r}}\frac{{\mathrm{d}\overline{\mathsf{\alpha }}}_{\mathrm{v}}}{{\mathrm{d}\mathsf{\alpha }}_{\mathrm{r}}}$$

The average steam volume fraction ${\overline{\mathsf{\alpha }}}_{ \mathrm{v}}$ in the down-comer tube is calculated using Equation (9) determined by [6].

$${\overline{\mathsf{\alpha }}}_{\mathrm{v}}=\left(\frac{{\mathsf{\rho }}_{\mathrm{w}}}{{\mathsf{\rho }}_{\mathrm{w}}{- \mathsf{\rho }}_{\mathrm{s}}}\right)(1- \frac{{\mathsf{\rho }}_{\mathrm{s}}}{{(\mathsf{\rho }}_{\mathrm{w}}\hspace{1em}{\mathsf{\rho }}_{\mathrm{s}}{) \mathsf{\alpha }}_{\mathrm{r}}}\mathrm{Ln}1+\frac{{(\mathsf{\rho }}_{\mathrm{w}}{- \mathsf{\rho }}_{\mathrm{s}}{) \mathsf{\alpha }}_{\mathrm{r}}}{{\mathsf{\rho }}_{\mathrm{s}}})$$

(9)

The volume fractions are derived depending on the drum pressure, and the quality values of the mixture are as follows:

$$\frac{{\mathrm{d}\overline{\mathsf{\alpha }}}_{\mathrm{v}}}{\mathrm{dP}}=\frac{1}{{{(\mathsf{\rho }}_{\mathrm{w}}{- \mathsf{\rho }}_{\mathrm{s}})}^2}{(\mathsf{\rho }}_{\mathrm{w}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{s}}}{\mathrm{dP}}{ - \mathsf{\rho }}_{\mathrm{s}}\frac{{\mathrm{d}\mathsf{\rho }}_{\mathrm{w}}}{\mathrm{dP}}\left)\right(1+\frac{{\mathsf{\rho }}_{\mathrm{w}}}{{\mathsf{\rho }}_{\mathrm{s}}}\frac{1}{1+\mathsf{\eta }} - \frac{{\mathsf{\rho }}_{\mathrm{w}}{+\mathsf{\rho }}_{\mathrm{s}}}{{\mathsf{\eta }\mathsf{\rho }}_{\mathrm{s}}}\mathrm{Ln} (1+\mathsf{\eta }))$$

$$\frac{{\mathrm{d}\overline{\mathsf{\alpha }}}_{\mathrm{v}}}{{\mathrm{d}\mathsf{\alpha }}_{\mathrm{r}}}=\frac{{\mathsf{\rho }}_{\mathrm{w}}}{{\mathsf{\rho }}_{\mathrm{s}}\mathsf{\eta }}\left(\frac{1}{\mathsf{\eta }}\mathrm{Ln}(1+\mathsf{\eta }) - \frac{1}{1+\mathsf{\eta }}\right)$$

where $\mathsf{\eta }=\frac{{\mathsf{\alpha } }_{\mathrm{r}}}{{\mathsf{\rho }}_{\mathrm{s}}}{(\mathsf{\rho }}_{\mathrm{w}}{ - \mathsf{\rho }}_{\mathrm{s}})$.

2.2. Determination of Terms Using Polynomials

The terms of Equations (3)–(8) were determined from the water steam tables as a function of saturation pressure ${\mathrm{P}}_{\mathrm{sat}}$. Engineering Equation Solver (EES^®) software was used to plot enthalpy, density, and temperature of the metal in the drum wall. The latter was assumed to be equal to the steam saturation temperature ${\mathrm{T}}_{\mathrm{m}}{=\mathrm{T}}_{\mathrm{sat}}$. Then, the curve was adjusted using linear regression to each graph, thus obtaining a polynomial of degree n as a function of the saturation pressure ${\mathrm{P}}_{\mathrm{sat}}$.

2.3. Methods for Ordinary Differential Equation Solution

To simulate the proposed model [6] and validate the results, it is necessary to solve the ordinary differential equations through a numerical method that allows for obtaining results adjusted to the dynamic behavior of the boiler. Considering the characteristics of each method, such as complexity for the programming, precision in the results, steps for the solution, and computational cost [16,17,18], it was decided to use the Runge–Kutta method of order 4 and Heun’s method and compare the results with those obtained using the MATLAB^® Ode45 function that solves this type of equation based on an explicit formula of RK 4 and RK5 [18].

2.3.1. Heun’s Method

This method is used to improve the estimation of the slope and involves the determination of two derivatives for the interval, one at the start point and one at the end point [16]. The two derivatives are averaged to obtain an improved estimate of the slope for the whole interval. Equations (10) and (11) shown below are used to solve differential equations with this method.

$${\mathrm{H}}_1=\mathrm{h} \mathrm{f}({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}})$$

(10)

$${\mathrm{H}}_2{=\mathrm{h} \mathrm{f}(\mathrm{x}}_{\mathrm{i}}+\frac{2}{3}{\mathrm{h},\mathrm{y}}_{\mathrm{i}}+\frac{2}{3}{\mathrm{H}}_1)$$

(11)

After having the values of ${\mathrm{H}}_1$ and ${\mathrm{H}}_2$, each of the points in the plane (x, y) are determined using the following Equations (12) and (13).

$${\mathrm{x}}_{\mathrm{i}+1}{=\mathrm{x}}_{\mathrm{i}}+\mathrm{h}$$

(12)

$${\mathrm{y}}_{\mathrm{i}+1}{=\mathrm{y}}_{\mathrm{i}}+\frac{1}{4}{(\mathrm{H}}_1{+3\mathrm{H}}_2)$$

(13)

Figure 3 shows the flow chart representing the code programmed for the solution of the ODE using Heun’s method.

2.3.2. Runge–Kutta 4th Order Method

The solution of this method was programmed in MATLAB^®, which develops multiple slope estimates to obtain an improved average slope for the interval. Each k represents a slope [16]. Equations (14)–(17) used by the Runge–Kutta fourth order method are shown below, where h is the step width.

$${\mathrm{k}}_1{=\mathrm{h} \mathrm{f}(\mathrm{x}}_{\mathrm{i}}{,\mathrm{y}}_{\mathrm{i}})$$

(14)

$${\mathrm{k}}_2{=\mathrm{h} \mathrm{f}(\mathrm{x}}_{\mathrm{i}}+\frac{1}{2}{\mathrm{h},\mathrm{y}}_{\mathrm{i}}+\frac{1}{2}{\mathrm{k}}_1)$$

(15)

$${\mathrm{k}}_3{=\mathrm{h} \mathrm{f}(\mathrm{x}}_{\mathrm{i}}+\frac{1}{2}{\mathrm{h},\mathrm{y}}_{\mathrm{i}}^{}+\frac{1}{2}{\mathrm{k}}_2)$$

(16)

$${\mathrm{k}}_4{=\mathrm{h} \mathrm{f}(\mathrm{x}}_{\mathrm{i}}{ + \mathrm{h},\mathrm{y}}_{\mathrm{i}}{ + \mathrm{k}}_3)$$

(17)

When each k is obtained, the points on the (x,y) axes of the plane are determined using Equations (18) and (19), as follows:

$${\mathrm{x}}_{\mathrm{i}+1}{=\mathrm{x}}_{\mathrm{i}}+\mathrm{h}$$

(18)

$${\mathrm{y}}_{\mathrm{i}+1}{=\mathrm{y}}_{\mathrm{i}} + \frac{1}{6}\left({\mathrm{k}}_1{ + 2\mathrm{k}}_2{ + 2\mathrm{k}}_3{ + \mathrm{k}}_4\right)$$

(19)

Figure 4 shows the flow chart representing the code programmed for the solution of the ODE using the Runge–Kutta 4th order method.

2.3.3. MATLAB^® Function Ode45

This method is based on an explicit formula of the Runge–Kutta (4th and 5th) order method, which is used to solve ODEs. The syntax shown in MATLAB^® is:

$$[\mathrm{t},\mathrm{x}]=\mathrm{Ode}45(\mathrm{odefun},\mathrm{tspan},\mathrm{x}0)$$

where

• x is a matrix where each column corresponds to the dependent variables, and t is the time vector,
• odefun is the name of the function to be evaluated,
• ${\mathrm{t}}_{\mathrm{span}}$ specifies the time interval, a vector of two numbers ${\mathrm{t}}_{\mathrm{span}}=[\mathrm{ti},\mathrm{tf}]$, start and end time,
• x0 is a vector that has the initial conditions of the variables to be evaluated.

2.4. Programming and Simulation of the Model in MATLAB^®

A code was developed in MATLAB^® for the implementation of the mathematical model that represents by simulation the behavior of the drum pressure P, the total water volume ${\forall }_{\mathrm{wt}}$, and mixture quality ${\mathsf{\alpha }}_{\mathrm{r}}$ in the steam generator. Figure 5 shows the algorithm that represents the code development in MATLAB^®.

2.5. Validation of Results

The data used for the simulation were taken from the P16–G16 unit operating in a power plant in Malmö, Sweden, which was reported by [6].

Table 1. Data used to perform the dynamic simulation of the steam generator.

Operation Data
Variable	Units	Value
Heat flux	$\left[\mathrm{kW}\right]$	86,000
Mass flow of water	${[\mathrm{kg}·\mathrm{s}}^{-1}]$	50
Mass flow of steam	$[{\mathrm{kg}·\mathrm{s}}^{-1}]$	50
Specificheat	$[\mathrm{kJ}·\left(\mathrm{kg}\right)·{\mathrm{K}}^{-1}]$	500
Enthalpy of water	${[\mathrm{kJ}·\mathrm{kg}}^{-1}]$	1080.9
Raiser volume	${[\mathrm{m}}^3]$	37
Down-comer area	${[\mathrm{m}}^2]$	0.1123
Total volume	${[\mathrm{m}}^3]$	88
Total mass	$\left[\mathrm{kg}\right]$	300,000
Gravity	${[\mathrm{m}·\mathrm{s}}^{-2}]$	9.81

The table was recreated according to [6].

shows the data used.

A time interval of 200 s was considered for the simulation of the dynamic behavior of the boiler. The percentage error was determined using Equation (20). To calculate the relative error, the reported reference value ${\mathrm{V}}_{\mathrm{Ref}}$ and the value obtained in the simulation ${\mathrm{V}}_{\mathrm{Sim}}$ are taken into account. An approximation error below 5% is acceptable when performing this type of study [19].

$$\%\mathrm{Error}=\left|\frac{{\mathrm{V}}_{\mathrm{Ref}}{-\mathrm{V}}_{\mathrm{Sim}}}{{\mathrm{V}}_{\mathrm{Ref}}}\right|\times 100\%$$

(20)

3. Analysis of Results

3.1. Polynomial Solution Using EES^®

Figure 6 shows the behavior of water and steam densities ${(\mathsf{\rho }}_{\mathrm{w}}{,\mathsf{\rho }}_{\mathrm{s}})$ as a function of saturation pressure ${\mathrm{P}}_{\mathrm{sat}}$. It is observed that when the saturation pressure increases, the saturated liquid reaches the saturated vapor point. For this reason, the density of the liquid ${\mathsf{\rho }}_{\mathrm{w}}$ decreases and the density of the steam ${\mathsf{\rho }}_{\mathrm{s}}$ increases until reaching the critical point.

Knowing this phenomenon, the curve fit is performed, resulting in the polynomials shown in the figure.

Figure 7 shows the behavior of the enthalpy of water and the enthalpy of steam as a function of the saturation pressure. It is necessary that as the saturation pressure increases, the saturated liquid reaches the saturated vapor point. For this reason, the enthalpy of water increases and the enthalpy of steam decreases until reaching the critical point.

Understanding what was done in the graph, the curve fit is performed, resulting in the polynomials shown in the figure.

Figure 8 shows that the temperature is increasing since, as the saturation pressure increases, the temperature approaches the critical point where the critical temperature of the water is 374 °C. Knowingthis, it is assumed that the saturation temperature of the liquid is equal to the saturation temperature of the metal ${(\mathrm{T}}_{\mathrm{sat}}{=\mathrm{T}}_{\mathrm{m}})$.

Thus, by knowingthe temperature behavior, an adjustment of the cross is made and thus we obtain the polynomial shown in the figure.

The polynomials shown in the previous figures are used to solve the differential equations that represent the behavior of the steam generator.

These polynomials apply only to the drum, since this is where there are two-phase liquid–steam conditions (mixture).

3.2. Analysis of Simulation Results

Figure 9 shows the behavior of the pressure, which is obtained using the mentioned numerical methods. These results are compared with the data shown in [6] and it is observed that the Runge–Kutta method and the Ode45 function are much closer to the data of the reported pressure graph, while the Heun’s method result is less exact.

Figure 10 shows the behavior of the total volume of water obtained using the numerical methods of Runge–Kutta of the fourth order, Heun, and the Ode45 function of MATLAB^®. In the figure, it can be seen that the Runge–Kutta method and the Ode45 function are the ones that would produce results close to the results reported in [6].

Figure 11 shows the behavior of the quality of the obtained mixture using the numerical methods of Runge–Kutta of the fourth order, Heun, and the Ode45 function of MATLAB^®. The figure shows that the results of the Runge–Kutta method and Heun’s method are further from the results reported in [6], while the Ode45 function of MATLAB^® is more accurate to the reported data.

3.3. Validation Using Approximation Error of the Results

To understand the dynamic operation of the boiler, it is important to select variables to achieve the physical interpretation of the process, which describe the storage of mass, energy, and momentum in the system. The three predominant variables in the process are: the total volume of water represents the accumulation of water in the system, the pressure P in the drum represents the energy received by the system, and, finally, the quality of the mixture that represents the distribution of water–steam that represents the fraction of steam in the risers [14]. Figure 12 shows the difference between the pressure data obtained, using the different numerical methods, and the pressure reports shown in [6]. This difference in results is obtained using Equation (20) and thus determines the approximation error with each of the methods and the report [6].

When conducting this study, it can be seen that the three solution goals are in an error range between 2.2% and 2.5%.

In Figure 12, it is observed that the three methods used present errors below 2.5% in the calculation of the pressure, and this indicates that they can be valid to determine this variable, however, with the RK4 method, a lower error is obtained in the results, stabilizing its value around 2%. For the pressure calculation, the highest error percentages are presented when using the Ode45 function with a maximum value of 2.5% for the simulated time. Figure 13 shows the difference between the total water volume data obtained, using the different numerical methods, and the volume reports shown in [6]. This difference in results is obtained using Equation (20) and thus determines the approximation error with each of the methods and the report [6].

Analyzing Figure 13, it can be seen that the three methods present a good behavior to determine the total volume of water with practically equal percentages of error and minimum values that are below 0.6% during the simulated operating time.

Figure 14 shows the difference between the quality data of the mixture obtained, using the different numerical methods and the quality reports shown in [6]. This difference in results is obtained using Equation (20) and thus determines the approximation error with each of the methods and the report in [6].

Figure 14 shows that for Heun’smethod and the RK4 function in the time range between 75 and 100 s, the error obtained in the calculation of the quality of the mixture tends to stabilize at a value of around 7% after an accelerated increase during the first seconds of the simulation. Undoubtedly, the Ode45 function allows for obtaining the lowest percentage of error in the calculation of the quality of the mixture, stabilizing its value around a maximum of 4.5%, much lower than those reached by RK4 and Heun’smethod of 7.5% and 8.5%, respectively.

4. Discussion

In the present work, the Runge–Kutta numerical method of order 4, Heun’s numerical method, and the MATLAB^® Ode 45 function were used to reproduce experimental values of the reference [6], which simulates the behavior of pressure, total volume of water, and quality of the mix in the boiler drum. Based on the results obtained, it is valid to affirm that the three methods can be used solve the ordinary differential equation and represent in general terms the dynamic behavior of the three most important variables in the operation of the boiler, however, when analyzing the percentage of error in the results, it is evidenced that the RK4 method is the best option to represent the behavior of the pressure in the drum and the Ode4 function to determine the quality of the mix. The three methods present very good results to determine the total volume of water with minimum percentages of error that do not exceed 0.6%.