Analysis of thermal and ﬂ ow phenomena in natural circulation boiler evaporator

Waterwall tubes


Introduction
Issues related to modelling transient thermal and flow phenomena occurring in a natural circulation boiler evaporator are attributed to their complex and strongly nonlinear characteristics in many studies. The complexity results from many factors, such as the high values of temperature and pressure, large surface areas of heating surfaces, the fact that the modelling process has to consider complex heat transfer processes occurring in the furnace chamber, local fouling of the waterwalls, configuration of burners, nonuniform supply of air and fuel to the burner system, and two-phase flow of the medium in the pipes. The nonlinearity is primarily caused by the dependence of the thermal and physical properties of water, steam, and flue gases on pressure and temperature or on temperature only.
The natural circulation of water in the boiler is possible due to the difference between the density of the wateresteam mixture in the risers and that of the water in the downcomers. This process is referred to as circulation. In large boilers, waterwall tubes (risers) located in heating zones, together with downcomers in nonheating areas with headers, comprise the circulation system or contour. Depending on the location, the circulation contours operate in different thermal load zones. Depending on the boiler design, their number varies from a few to approximately a dozen. The validity of the evaporator operation depends on the device design, quality of assembly, and the conditions in which the device operates.
The validity of water circulation in power boiler evaporators is the subject of numerous studies aiming to ensure the reliability and durability of the boiler components. In other words, their objective is to ensure a continuous and one-way water circulation in the evaporator tubes, ensuring that the evaporator is cooled appropriately.
A simplified lumped mass model of the boiler evaporator was defined previously [1,2]. However, the primary task of these studies is to determine the optimal change in the temperature and pressure of the medium in the boiler evaporator that will shorten the boiler start-up from the cold state.
It is noteworthy that water heaters (economisers), boiler evaporators, and steam superheaters/reheaters with a complex structure are still modelled using simple lumped parameter models [3,4] that are intended primarily to satisfy the needs of the boiler automatic systems of regulation. This is due to the intricacies of mathematically distributed parameter models. The advantage of models with lumped thermal parameters is their short computation time; meanwhile, the drawback is that excessive simplifications are made.
Many studies are related to modelling steam superheaters and feed water heaters. Some of them deserve special attention. For example, a distributed parameter model of a steam superheater and the feed water heater was proposed in Refs. [5,6]. The model in Ref. [5] was built using mass, momentum, and energy balance equations. Because the thermal and physical properties of steam, water, and flue gas depend on temperature and pressure, or temperature only, the equations are nonlinear. In Ref. [6], the dynamics of superheaters or economisers is modelled using simplified models based on the energy balance equation, omitting the equations of the balance of mass and momentum. It is noteworthy that the applied mathematical models were verified experimentally by comparing the results of measurements performed on a real facility to those obtained from numerical computations. The models proposed in the studies mentioned above are one dimensional (1D).
Many studies concern the simulations of thermal and flow phenomena occurring in the supercritical boiler furnace chamber waterwalls [7,8]. A distributed parameter mathematical model that enables such simulations is proposed in Ref. [7]. It is a 1D model based on solving the equations describing the principles of mass, momentum, and energy conservation. They are solved using appropriate difference schemes. The model is verified computationally by comparing its results with those obtained from the analytical solution for the transient state. Moreover, a comparison was also performed with the results obtained after solving the balance equations using the CrankeNicolson method. In both cases, excellent agreement was obtained between the results. The proposed model may also be used as an element of the power unit simulator. A conjugate heat transfer model was proposed in Ref. [8]. A 1D model solution is proposed on the side of the working medium. In the waterwall finned tube, 20 control volumes are separated in each analysed cross section and equations of the threedimensional (3D) transient-state heat transfer are written and solved for them. The results of the developed 1D/3D model are verified by the simulation calculations of the furnace chamber waterwalls installed in a supercritical boiler currently operating in a power station in Poland. Only the single-phase flow is considered in the studies mentioned above, and the results obtained from the simulation are not compared to the measurements.

Literature survey
The recent literature survey indicates that models of entire boilers or whole sub-and supercritical power plants are available. One example is [9], where Liu et al. developed a model of a 1000-MW once-through boiler using neural networks. The model results are compared to those obtained from measurements performed in a real facility. The comparison concerns the output electric power, live steam pressure, and fluid temperature at the separator outlet. It was proven that previously used linear models can well predict the selected working points. However, that effectiveness became worse if the models were used for a real installation in various operating conditions. The obtained results indicate the considerable limitations of linear models in terms of their ability to represent the highly nonlinear characteristics of the power plant dynamics.
A nonlinear approach to the problem of determining the medium enthalpy (temperature) between the water heater inlet and the once-through boiler evaporator outlet is presented in Ref. [10].
The control of steam parameters, such as pressure and temperature, is necessary in power plants to ensure their full reliability. Steam parameters downstream the evaporator and upstream the water heater should be controlled to obtain the required levels of loaddependent nominal values and to ensure that the boiler operates smoothly. This is achieved by controlling the feed water supply to the evaporator inlet and adjusting the mass flow according to the changes in loads. A lumped parameter model of the evaporator is used in Ref. [10]. The model considers the energy balance equation only.
The publications mentioned above concern simulations of the operation of once-through boilers performed using models based on considerable simplifications. One such simplification is by not considering the thermal load nonuniformity along the chamber height or using a lumped parameter model of the evaporator.
Many studies are devoted to the analysis of the waterwall tube operation. Dhanuskodi et al. [11] applied a model using artificial neural networks to determine the temperature of the waterwall tube walls. The model considers the steady-state distributions of pressure and temperature. The nonuniform distribution of loads is also considered. Zheng et al. [12] developed a heat-transfer model in spiral waterwall tubes. It is a distributed parameter model based on the fluid's three-dimensional temperature field, calculated temperature of the wall, mass flow rate (mass flux) distribution, and thermal load of the walls. The same model is proposed in Ref. [13]. Pan et al. [14] developed a model based on mass, momentum, and energy conservation equations. The equations were solved iteratively using the quasi-Newton method, which produced thermoehydraulic characteristics of the boiler heating surfaces.
Many studies are related to designing general models of entire power plants using commercial computer programs. A model of the dynamics of large power units was developed in Ref. [15]. The APROS software, developed by the VTT Technical Research Centre of Finland in collaboration with the Fortum company, is used for this purpose [16]. APROS includes two thermal hydraulic modules, i.e. the homogeneous and heterogeneous flow model. The homogeneous model is applied for single-phase flows. For two-phase flows (evaporator, condensers), the APROS heterogeneous slip model is used. The high accuracy of the mathematical model applied in the APROS software was confirmed in Refs. [17,18]. The developed model was verified against the power plant operating data.
An integrated model of a boiler in the steady state was developed in Ref. [19]. It is composed of a few submodels, and a modified NewtoneRaphson method was used to obtain the numerical solutions of the equations in each submodel. Each heat transfer surface was regarded as a lumped-mass one. The input data for each component of the model, such as temperature, pressure, and mass flow of the medium, were obtained from real installations. The results obtained from the model in steady-state conditions were compared to the real data of a power plant.
Two-phase flows in power boilers have been emphasised significantly in literature, for example in boiler evaporators. Boiler evaporators designed for subcritical parameters are characterised by two-phase flows in the entire range of the power unit loads. In this respect, their operation is different compared to supercritical boilers.
An open natural circulation cycle was proposed in Ref. [20] to design a passive containment cooling system. A numerical code was developed to simulate the process. It is used to analyse the flow characteristic in a small-scale model. The authors proposed a simplified scaling method that was verified through a comparative analysis. The scaling process was divided into heat-transfer scaling in the heat exchanger and pressure-loss scaling in risers and downcomers. The open natural circulation cycle was tested on a laboratory stand. The circulation cycle additionally included a heat exchanger. The natural circulation cycle was modelled using simplified mass, momentum, and energy balance equations. A homogeneous model was assumed to calculate the pressure losses in the two-phase flow.
A nonlinear dynamic model of a natural circulation evaporator including a single boiler drum was presented in Ref. [21]. It is a distributed parameter model based exclusively on the fundamental physical laws of mass, momentum, and energy conservation. The presented approach to modelling the boiler operation was based on an analysis of the physical phenomena occurring in the boiler drum, as well as in risers and waterwall tubes. It also considered the boiler's unique design features. The model effectiveness was verified by comparing its results with those obtained from a distributed control system. As a simplification, a lumped parameter model was used for outlet tubes and downcomers. An analysis was conducted to verify the model performance in irregularity conditions in the boiler operation. Special emphasis was placed on the impact of the irregularities on changes in the water level in the boiler drum. The simulation results were compared with the experimental results. Additionally, the simulation results were presented for other key variables of the boiler system that are typically difficult to measure.
A mathematical model based on mass and energy balance equations for the water heater, boiler drum, and superheater was presented in Ref. [22]. The results obtained from the model were compared to those measured in a power plant. The temperature and pressure values at the water heater and boiler drum outlet were analysed. The mean-square error (MSE) was calculated for the comparisons. Analysis was not performed in the studies of evaporator circulation cycle or two-phase flows. The focus was on the parameters of steam at the boiler drum outlet.
It is noteworthy to focus on [23] that presents a mixed model of a power plant operating in combination with a heat recovery steam generator. A simplified model of the boiler drum was used and the water natural circulation in the system was analysed using the SPORTS code developed by Chatoorgoon [24]. A simplified algorithm was applied in the code to analyse a two-phase flow. The results of the proposed model of the boiler dynamics were compared to the results obtained using two other models, i.e. the model developed by Astrom and Bell [25] and that proposed by Bhambare, Mitra, and Gaitonde [26]. The former is a lumped parameter model. The latter is based on the basic equations of mass and energy conservation, and on empirical correlations for twophase flows. Neither of the publications considers or verifies the boiler operation to a greater extent, e.g. during the boiler start-up. The model proposed in Ref. [23] uses distributed parameters, and the analysis concerns are one downcomer and one riser. The mass, momentum, and energy conservation equations are written such that time derivatives are used on their left side, whereas space derivatives are used on the right. The authors suggest that such a nonlinear system of equations should be solved using the available numerical codes. The model validation covers a wide range of boiler operation, i.e. start-ups from cold and the hot states. The results are compared to those obtained using the models described in Refs. [25,26].

Novelty of the present study
The performed literature survey indicates that no mathematical models can describe the complex transient-state thermal and flow processes occurring in power boiler evaporators with adequate accuracy.
Most studies do not pertain to distributed parameter models e the analysis is typically performed using lumped parameter models that are based on different forms of simplifications. Thus far, none have attempted to perform the experimental verification of the mathematical heat-transfer models concerning the evaporator tubes of natural circulation boilers. A detailed analysis of the evaporator (downcomers and risers) operation is still lacking. Therefore, this paper presents a new mathematical model describing the thermal and flow phenomena occurring in the boiler evaporator and the model is verified experimentally.
The proposed distributed parameter model is based on solving balance equations describing the principles of mass, momentum, and energy conservation, and enables an analysis of steady and transient conditions processes. A system of appropriate ordinary differential equations is written with space derivatives on the left side of the equation. The time derivatives on the right-hand side are replaced with backward difference quotients. The system is solved using the RungeeKutta method. Such an approach is characterised by high computational stability and accuracy. Owing to unique measurements, the developed model of the boiler evaporator enables the determination of the pressure distribution, the mass flow rate and the dryness factor of steam, as well as the water circulation multiplicity in the evaporator.
A flow meter was installed in the downcomer to enable the testing of natural circulation in the boiler evaporator. A series of measurements were performed in a natural circulation boiler with an output of 210 Â 10 3 kg/h to verify the proposed method of modelling transient-state processes occurring in the boiler evaporator. Considering the current requirements related to the steam power boiler operation flexibility, the start-up and shutdown processes and rapid changes in loads, e.g. by opening the turbine valves, should be performed while maintaining the necessary levels of natural circulation. Hence, the measuring points existing in the facility were used. In some cases, some new ones were installed. The measurements were performed during the boiler acceptance tests after a general overhaul. The results obtained during the boiler start-up, shutdown, and steady-state operation were compared to the numerical calculation results.
A novel contribution of the current study is the developed mathematical model, which is characterised by its high stability and accuracy in performing calculations during rapid changes in the boiler load. The results of the presented measurements are also unique.

Mathematical model
This section presents the method of determining the distribution of the fluid mass flow (mass flux), pressure (pressure drop), and enthalpy (temperature) along the furnace chamber waterwall tube. The relevant mass, momentum, and energy conservation equations are solved. The equations are derived in Ref. [27], and the expressions presented below are their final forms: Mass conservation equation: Energy balance equation: After some simplifications and transformations, the balance equations (1)e(3) can be written in a form where the space derivatives and time derivatives are obtained on the left and right sides of the equations, respectively [27]: The time derivatives on the right-hand side are replaced with backward difference quotients. This system of ordinary differential equations is solved using the RungeeKutta method.
After the changes above are introduced, the energy balance equation is of the following form: The fluid density is obtained as a function of enthalpy and pressure: Solving the mass and momentum conservation equations, the following are obtained, respectively: The changes in the fluid temperature are obtained as a function of enthalpy and pressure: The subscript 'j' in equations 7e11 denotes the subsequent number of analysed cross sections and varies from j ¼ 2 … M. All thermophysical properties of the fluid are calculated online. The segment of the fluid flow path from the boiler drum to the downcomers and back to the drum through risers is divided into nodes. The first and last point of the circulation cycle, i.e. 1 and M, are the downcomer and the riser nodes characterised by the same pressure, which is equal to the pressure value in the boiler drum. The equality of pressure in two extreme nodes is obtained iteratively using the RegulaeFalsi method.
Moreover, the presented method requires that the Courant condition should be satisfied [28]: The following data were assumed for the computations: This is a stability condition imposed on the time step to ensure the transposition of the numerical solution with a velocity Dz/Dt higher than the physical velocity w.
The flow occurring in a natural circulation boiler evaporator is a single-and a two-phase flow. The friction-related pressure decrease in formula (10) is defined based on the DarcyeWeisbach equation, which is valid for the fluid from the evaporator inlet to the point where the water reaches the saturation temperature: where: d in e inner diameter of a duct with a circular cross section or equivalent diameter, f e friction factor.
The friction factor is extremely difficult to determine in a mixture of water and steam. A number of mathematical models have been proposed to find the component being the effect of friction. The flow of the mixture can be considered as homogeneous, or a slip flow model can be applied. Different models of determination of pressure losses in two-phase flows are presented in Ref. [29], and the obtained results are compared to each other.
The use of the homogenous model involves assuming that the properties of the fluid are averaged for the entire mixture. Another assumption is that the velocity values of the gaseous and the liquid phase are the same, which implies that there is no slip between the phases. For the homogeneous model, the decrease in pressure due to friction is defined by the following equation: where: f TP e two-phase friction factor, Ge mass flux, kg/(m 2 s), r TP e two-phase density, kg/m 3 .
Classical relations, such as the Blasius formula, are used to determine the friction factor. In this case, the Reynolds number obtained for the averaged properties of the mixture is applied [30,31]: where: h TP , h G , h L e dynamic viscosity of the homogenous mixture, the gaseous phase and the liquid phase, respectively, x e dryness factor.
The two-phase mixture density is obtained using the following formula: where: r G , r L e density of the gaseous and the liquid phase, respectively.
Using slip models involves calculating the two-phase flow multiplier that enables the friction-related pressure decrease to be obtained. The concepts included in this group are the LockharteMartinelli model, Friedel model, Chisholm model, and MartinellieNelson graphical method [31].

The LockharteMartinelli model
To apply the LockharteMartinelli model the two-phase multipliers should be found. It makes it possible to determine the pressure gradient caused by frictional losses arising in the flow of the mixture [31]: The following formula is used to define the two-phase multiplier used in equation (18): The Martinelli parameter X appears in equation (19). It is defined as The determination of the LockharteMartinelli multiplier is related to establishing the flow character of individual phases. The parameter C is selected depending on whether the flow of the phases is laminar or turbulent (cf. Table 1).
The Friedel model was built using 25,000 measuring points. The model is applied for two-phase flows in vertical and horizontal tubes. The following formula is used to determine the two-phase multiplier [31]: Two-phase multiplier f 2 Lo is related to the flow of a mixture with properties related to the properties of the liquid phase. Coefficients E, F and H, as well as the Froude number (Fr) and the Weber number (We) in formula (21) are described in Ref. [31].
The Friedel model is recommended if the ratio between the liquid and gaseous phase dynamic viscosity is (h L /h G ) 1000. If the condition is satisfied, equation (21) may be used in the entire range of the steam dryness factor values.

The Chisholm Model
The Chisholm model is empirical. It can be used for a wide range of values of the dryness factor and steam pressure. The following relation is used to define the two-phase multiplier [31]: Bx 0:5ð2ÀnÞ ð1 À xÞ 0:5ð2ÀnÞ þ x 2Àn i ; (22) where: Y 2 e the ratio between pressure gradients of the gaseous and the liquid phase in single-phase flows B e mass flow rate (mass flux)-dependent coefficient n e exponent defining the coefficient of the flow resistance. For the Blasius formula, n ¼ 0.25.
The methods of finding parameters Y 2 and B are presented in Ref. [31].

The MartinellieNelson graphical method
The MartinellieNelson method was developed based on flows of the mixture of steam and water through horizontal pipes at pressure values ranging from 0.689 to 20.7 MPa. The two-phase multiplier needed to determine the losses arising due to friction is read from a chart developed using different steam pressures and different values of the dryness factor (x) (cf. Fig. 1). If pressure gets higher than the water critical point, the multiplier is f 2 Lo z1 [32]. The two-phase multiplier read from the chart makes it possible to determine the frictional decrease in pressure using the following relation: where. dp f ;LO dz e friction-related losses in the flow of a two-phase mixture characterised by properties related to the liquid phase properties, Pa/m.

Experimental verification
This section presents the experimental verification of the developed model. First, the values of frictional losses arising in a 210 Â 10 3 kg/h boiler waterwall tube and calculated using the homogenous and the phase slip models presented herein are compared. The experimental verification consists of comparing the calculated and measured values of the water mass flow in the downcomers and in obtaining the dryness factor and circulation multiplicity. The numerical computations were performed using in-house programs written in the Fortran language [34].
The water mass flow rate in the downcomer was measured in a boiler with the output of 210 Â 10 3 kg/h. It is a single-drum natural circulation boiler (cf. Fig. 2). The boiler generates high-pressure superheated steam fed into the power plant's primary duct. The boiler rated output totals 210 Â 10 3 kg/h at the steam outlet pressure of 9.8 MPa and temperature of 540 C. The platen superheater and the following stages of the steam superheaters are located in the combustion gas bridge: 1st stage (KPP-1), 2nd stage (KPP-2) and 3rd stage (KPP-3). The two stages of the economizer (ECO-1 and ECO-2), and a tubular-type air heater are placed in the convective section. For verification, the boiler evaporator operation was analysed during start-up, in the steady-state mode, and also during the boiler acceptance tests. The acceptance tests were made prior to the boiler re-commissioning after a general overhaul.
The steamewater system of the boiler evaporator comprises a drum, downcomers, risers, pipelines with fittings, and measuring instruments. The boiler water flows from the drum into 10 ∅273 Â 25 downcomers through 54 ∅108 Â 10 tubes and continues through four ∅130 Â 10 tubes to reach the waterwall tubes. Having passed through the waterwall tubes, the steamewater mixture flows into the drum through ∅76 Â 10 tubes. Fig. 3 presents a diagram of the analysed boiler circulation contour with a marked point of the water mass flow measurement in the downcomer. The meter was installed at a height of approximately 10 m. The method of meter fixing in the tube was described in Ref. [35]. The boiler water was transferred from the drum through five ∅108 Â 10 tubes to a ∅273 Â 25 downcomer; subsequently, from flowing through three ∅130 Â 10 tubes, it reaches into the header followed by 30 ∅60 Â 5 waterwall tubes. Twelve ∅76 Â 10 tubes carried the steamewater mixture into the drum.
The distribution of the thermal load along the waterwall tubes  was calculated using the zone method. Details of the calculation method are presented in Ref. [36]. The height-dependent load of the furnace chamber walls adopted for the calculations, for the boiler output of 210 Â 10 3 kg/h, is shown in Fig. 4. The saturation temperature is reached and the boiler water partially evaporates in the waterwall tubes. To compare the frictional losses, the calculation methodology is typically based on the evaporator division into a few zones [37]. In this study, the waterwall tubes are divided into two evaporation zones (cf. Fig. 3) [29]. At the first and second zone outlets, the steam dryness factors are x ¼ 0.04 and x ¼ 0.1, respectively. An increase in the dryness factor value involves an increase in frictional losses. The comparison between the resistances obtained using the homogenous and phase slip models is presented in Table 2.
The calculation results indicate that the friction-related pressure decreases are of similar values for most of the models. Some differences occur for the LockharteMartinelli model, which was developed based on the flow measurements of airewater and aireoil mixtures. The share of frictional pressure losses in the overall pressure drop is relatively small. Hence, the selection of model to determine friction-related pressure losses is irrelevant e the selection has little impact on the total pressure drop in the circulation contour of natural circulation boilers. Further, in equation (10), the frictional losses are calculated using the homogeneous model.
Figs. 5e7 present the measurement results and the results of the calculations performed during the boiler start-up. Fig. 5 shows a curve illustrating the measured changes in the steam pressure in the drum (saturation pressure).      (Table 3). It may therefore be stated that the agreement between the results is satisfactory, considering the fact that they relate to the boiler transient operation during start-up. The performance of the proposed mathematical model at such rapid changes is stable. Fig. 7 presents the determined circulation multiplicity and the dryness factor of the steamewater mixture flowing into the drum. As shown in Fig. 7, a certain instability occurs in the circulation multiplicity (j z 350 in time 3000e4000 s); this could be due to the small amount of generated steam at the beginning of the boiler operation. It is only after the boiler achieves the steady-state operation mode (after approximately 12,000e15,000 s) that the circulation multiplicity stabilises at 8.5e9.5.
Figs. 8e10 show the results of the measurements and calculations performed for the boiler operation under variable loads (acceptance tests). During the boiler acceptance tests, the boiler steam output was decreased from 180 Â 10 3 kg/h (50 kg/s) to approximately 150 Â 10 3 kg/h (41.7 kg/s) and subsequently increased to approximately 200 Â 10 3 kg/h (55.6 kg/s).
The measurements continued for a period of approximately 22 h. The changes in the drum pressure are illustrated in Fig. 8. Fig. 9 presents the results of the boiler steam output measurements (Curve 1) and the comparison between the measured and calculated values of the boiler water mass flow in the downcomer (Curve 2). In this case, the error between the results of the measurements and the calculations is MSE ¼ 1.862 (Table 3). The result is fully satisfactory. The computational algorithm remains stable during long periods of the boiler operation (t ¼ 80,000 s). Towards the end of the measurements (from the time of t~64,000 s), at a significant change in the steam mass flow at the boiler outlet from 160 Â 10 3 kg/h to 210 Â 10 3 kg/h, slight differences between the results can be observed. Fig. 10 presents the determined circulation multiplicity and the steamewater mixture dryness factor.
As shown in Fig. 10, regardless of the load, the flow in the downcomer was stable; only the circulation multiplicity varied within the range of j ¼ 7e9.5. In the time interval of (t ¼ 0e30,000 s) the boiler pressure totals 10.4 MPa (cf. Fig. 8), corresponding to a circulation multiplicity of approximately j ¼ 8. If the boiler drum pressure is decreased to the value of approximately p ¼ 10 MPa (t ¼ 30,000e63,000 s), the circulation multiplicity is approximately j ¼ 9.5. In the last interval of the boiler testing at t ¼ 63,000e80,000 s, the pressure measured in the boiler drum is  Table 3 Mean-square error.
Mass flow rate    Fig. 12 e Curve 1). Fig. 11 shows the curve illustrating changes in the pressure in the drum, and Fig. 12 presents a comparison between the results obtained from the water mass flow measurements and calculations (Curves 2). The error between the results of the measurements and the calculations is MSE ¼ 0.237 (Table 3). In this case, the result is excellent. During the boiler steady-state operation, the proposed mathematical model is highly stable and demonstrates high computational accuracy. Fig. 13 presents the determined circulation multiplicity and the dryness factor. As shown in Fig. 13, during the boiler steady-state operation, the obtained mean circulation multiplicity totals j ¼ 9.8 at the boiler drum mean pressure value of p ¼ 9.92 MPa. The circulation multiplicity value calculated for this pressure corresponds to the multiplicity occurring in real facilities.
To determine the effectiveness of the method, the mean-square    error is calculated. The MSEs for the performed calculations and measurements were determined using the following relation: where.
n e number of measurements; y i e measured value; y dt e expected (calculated) value. The results were obtained from the calculations in Table 3.
The proposed model can be a useful tool for online monitoring the thermal and flow processes occurring in the boiler pressure elements. It enables circulation control in the evaporator to prevent water stagnation in the boiler risers during fast changes in loads. Apart from monitoring flows through tubes, it is also possible to control the thermal processes occurring in the furnace chamber. Such a system, using the presented model and water flow measurements in the boiler downcomers, is installed in a power plant in Poland. The start-up process is performed and controlled such that the appropriate levels of circulation in the evaporator are ensured, thus enabling the levels to exceed the thermal stresses arising in the boiler thick-walled elements. The developed model can also be implemented in other steam-generating units in the range of sub-and supercritical pressure values.

Conclusions
The strict control of boiler evaporator operation is becoming increasingly important. Therefore, there has been strong motivation and growing interest towards the development of a mathematical model that can simulate the boiler evaporator performance not only at the boiler constant load, but also in a wide range of dynamic changes including start-ups and varying loads.
A mathematical model enabling the simulations of thermal and flow phenomena occurring in the evaporator of a natural circulation boiler was presented herein. It is a distributed parameter model based on solving the equations describing the principles of mass, momentum, and energy conservation. As proposed herein, the balance equations were solved using the RungeeKutta method.
The analysis of the calculations of friction-related losses performed in this study using different models enabled the differences between the obtained results to be observed. The differences were the effect of the different approaches to the two-phase flow multiplier determination. In the analysed conditions of the twophase flow in the natural circulation boiler evaporator, considering the evaporator pressure, the circulation multiplicity was high (j z 10); consequently, the dryness factor of the steamewater mixture at the drum inlet was small. This implied that in this type of flow, the friction-related pressure decrease was significantly smaller than the sum of the two other quantities comprising the total pressure decrease, i.e. the sum of the hydrostatic pressure drop and the decrease resulting from the change in the twoephase mixture momentum. In this situation, the accuracy of the method of the friction-related constituent determination is negligible. In this study, frictional pressure losses in the two-phase flow were determined using a homogenous model.
The developed method was verified experimentally. The natural circulation in the evaporator of a boiler with the steam output of 210 Â 10 3 kg/h was tested. The tests were performed during the boiler start-up and shutdown, and during the boiler steady-state operation. The calculation results were compared to the measurement results, and a fully satisfactory agreement was obtained. The developed mathematical model can be used to simulate thermal and flow phenomena occurring in both evaporators of natural circulation boilers and supercritical once-through boilers.