*Article* **0D Dynamic Modeling and Experimental Characterization of a Biomass Boiler with Mass and Energy Balance**

**Fateh Mameri 1, Eric Delacourt 1,2, Céline Morin 1,2,\* and Jesse Schiffler 3**

**Abstract:** The paper presents an experimental study and a 0D dynamic modeling of a biomass boiler based on the Bond Graph formalism from mass and energy balance. The biomass boiler investigated in this study is an automatic pellet boiler with a nominal power of 30 kW with a fixed bed. The balances allow to model as time function the flue gas enthalpy flux variation and the thermal transfers between the flue gas and the walls of the boiler subsystems. The main objective is to build a model to represent the dynamic thermal behavior of the boiler. Indeed, small domestic boilers have discontinuous operating phases when the set temperature is reached. The global thermal transfer coefficients for the boiler subsystems are obtained according to an iterative calculation by inverse method. The boiler has an average efficiency of 67.5% under our operating conditions and the radiation is the dominant thermal transfer by reaching 97.6% of the total thermal transfers inside the combustion chamber. The understanding of the dynamic behavior of the boiler during the operating phases allows to evaluate its energy performances. The proposed model is both stimulated and validated using experimental results carried out on the boiler.

**Keywords:** energy balance; biomass boiler; heat exchanger; 0D modeling; Bond Graph; global thermal transfers; inverse method

### **1. Introduction**

Biomass plays a significant role in the development of clean and sustainable heat production processes with a large reduction of CO2 emissions [1]. There are multiple ways to exploit energy potential of biomass, e.g., by pyrolysis [2], gasification or other bio-chemical processes using bacteria to generate gaseous and liquid biofuels or by direct combustion to generate heat and electricity [3–5]. Even if biomass has a lower calorific value than other fuels, such as fossil fuels, this source of energy remains cleaner with some reserves [6]. The biomass can be valued for the simultaneous production of heat and electricity from CHP (Combined Heat and Power) plants [7,8].

In the thermal conversion of biomass, there are multiple physical and chemical processes that have an influence on the performances of industrial and domestic applications, such as furnaces, industrial burners and biomass boilers [9], the exergy analysis must be used in order to find the best way to recover the maximum of mechanical work in a CHP (combined heat and power) unit. Biomass boilers provide a direct conversion of biomass into energy by combustion. They are widely investigated in several configurations according to delivered power: biomass domestic boiler of 24 kW, 27 kW and 32 kW [10–12], industrial biomass boiler of 4 MW [13].

Dynamic modeling of energy systems can be used for the design, the optimization or the control of the studied process. Tognoli and Najafi [14] provided a detailed dynamic model of an industrial fire-tube boiler with five different geometrical configurations. The dynamic model developed consists of two main sections separated on the flue gas side and

**Citation:** Mameri, F.; Delacourt, E.; Morin, C.; Schiffler, J. 0D Dynamic Modeling and Experimental Characterization of a Biomass Boiler with Mass and Energy Balance. *Entropy* **2022**, *24*, 202. https:// doi.org/10.3390/e24020202

Academic Editor: T M Indra Mahlia

Received: 30 November 2021 Accepted: 28 December 2021 Published: 28 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the evaporating shell. Both sides are integrated employing an energy balance. Then, a PID tuning was implemented for each boiler to control the vapor pressure, while responding to a demand with variable mass flow rate. The operation of the boilers was simulated to meet four different steam demand profiles. A wood pellet micro-cogeneration system with steam engine was modeled by Bouvenot et al. [15] and implemented in the TRNSYS code. Both theoretical and experimental approaches have been adopted to develop the model. The authors presented the dynamic response of the installation and took into account the steady and transient states. A dynamic model applied to two biomass boilers with nominal power of 6 and 12 kW was presented by Carlon et al. [16]. The model developed with TRNSYS calculates the mass and energy balances of the boilers under time variable inputs. It describes the operation of the boiler under dynamic conditions and provides the chemical composition of the flue gases from the chemical composition of the wood pellets and the value of the excess air and by adopting the assumption of a complete conversion of the mass of fuel. The model has been tested for two modes of boiler operating conditions: full and variable load and steady and transient states. The results of the modeling showed a better agreement with the experimental data during steady operation as well as in dynamic mode.

The modeling of thermofluidic systems related to heat and power generation are also described in terms of mechanical work generation processes. We can note for example, a study on the modeling of ORC (Organic Rankine Cycle) systems investigated by Ziviani et al. [17]. The authors presented an overview of the problems related to ORC modeling and developed an efficient and powerful simulation for an ORC system adapted to the exploitation of low-grade thermal energy. Other physical systems are modeled like ECE (External Combustion Engine), for example Stirling or Ericsson engines [18–20]. Due to their promising future paths for energy cogeneration by coupling them with thermodynamic systems, small biomass boilers were also studied from this type of approach. However, they raised several concerns, ranging from design to dynamic control [21]. Inappropriate power requirement definition and inadequate control can affect the boiler performances and reduce its efficiency. To facilitate the design process and overcome upstream design failing, the modeling represents a very interesting approach.

The dynamic behavior of this kind of system is generally described by non-linear differential equations. A suitable method, as the Bond Graph formalism, is necessary to well understand physical interactions in a such thermofluidic system. Therefore, an appropriate model that represents a system involving energy transfers can be extracted in a structured way [22]. The Bond Graph method is based on a graphic structure representing the power exchanges between different physical entities considered in multidisciplinary dynamic systems. It was initiated by Paynter in 1961 [23] and then developed by Karnopp and Rosenberg [24]. This tool is adapted to the modeling of the physical processes involved in different energy fields (hydraulic, mechanical, electrical, chemical and thermal). Bond Graph formalism allows to develop a parametrized model with an unified language that interprets the power transfers within the system considered explicitly through its graphic structure. Bond Graph investigations have been carried out on energy systems such as hot air engines (Ericsson engine [19]), Heating, Ventilation and Air-Conditioning (HVAC) systems [25], industrial biomass boiler [26], endoreversible heat engine [27], thermo-hydraulic system [28] and in chemical engineering [29]. Ould-Bouamama et al. [30] have developed a dynamic model using Bond Graph methodology for an industrial chemical reactor. The purpose of this application is to design a monitoring and survival platform in case of failure.

The modeling of biomass boilers operation has been the subject of several studies. Mathematical models based on thermodynamic laws have been developed to represent the dynamic behavior of the boilers during operating phases such as start-ups and load changes [31]. Åström and Bell [32] developed a simple non-linear model based on the first law of thermodynamics and configured with the basic design data of the boiler. Sandberg et al. [33] presented a dynamic model based on the mass and energy balances of a biomass boiler to study the effect of fouling on boiler performances. Table 1 summarizes

**Table 1.** Review of energy system modeling.

furnace, reactor and engine).

Published studies on the dynamic modeling of boilers often refer to black box or grey box models. There are some studies describing white box models but the detail of the modeling is often incomplete (use of components of commercial tool libraries rather opaque or description of physical phenomena modeled without specifying the interactions between them).

the experimental and numerical studies of the literature about different systems (boiler,

In this work, a 0D dynamic model of a domestic biomass boiler is provided using the Bond Graph formalism to simulate its dynamic behavior and to understand all the heat transfers involved in the boiler. The 0D model is based on mass and energy balances. It characterizes all the heat exchanges between the flue gas and the walls of the subsystems constituting the boiler. This dynamic modeling makes sense with domestic boilers whose operation is typically discontinuous unlike larger industrial boilers. The thermal needs of the house are variable which results in intermittent operation of the boiler. The strength of dynamic zonal modeling is to be able to predict the time evolution of different state variables of a complex system by coupling some fields of physics (mechanics, thermodynamics, . . . ). Moreover, it is possible, during the simulation, to insert time boundary conditions from in-situ measurements. The local evolution of the state variables is much less detailed than with CFD modeling but the dependencies of one zone with another are better taken into account with a 0D dynamic modeling. Moreover, CFD simulations are generally performed in steady state (averaged) because of the high computational cost in unsteady state, contrary to the dynamic 0D model which is able to predict the impacts on coupled systems. CFD and dynamic 0D modeling are therefore to be implemented according to the targeted objectives and the simulations results can hardly be compared. However, they can be efficiently coupled in multi-scale approaches. The objective of this study is to model the dynamic behavior of the boiler during the operating phases in order to take into account the variability of the heat production with regard to the thermal load of the heating network

or of any system which could be connected to it (hot air machine for example in the case of a CHP plant).

Compared to other dynamic modeling, the interest of Bond Graph methodology by its explicit graphic structure is to make clearer the modeling process of coupled multi-physical phenomena with blocks linked together by power links where effort and flow variables as well as causality are explicit. This methodology is very well adapted to model a system with thermal and mass transfers described with linear or non-linear differential equations.

The Bond Graph formalism was not developed in the literature to study the thermal transfers between the different fluids in a low power biomass boiler during the transient operating phases but it is increasingly used for the modeling of thermofluidic systems in general. Thanks to this formalism and its clarity, it is then possible to highlight the boiler components where the thermal transfers must be optimized and to understand the physical interactions. Moreover, there is a lack of experimental data for low power biomass boiler in the literature, these data are essential to develop a dynamic model by considering the real operating cycle of the boiler. An innovative way is used in this study by coupling experimental values and 0D modeling at each time step of the calculation with an analysis of energy performances for a domestic biomass boiler.

In the paper, the biomass boiler is described with all sensors used for the measurement of temperatures and mass flow rates. Then, the methodology is explained and the dynamic model 0D of the boiler is presented. Experimental and numerical results are discussed.

### **2. Description of the Biomass Boiler**

### *2.1. Experimental Setup*

The study is focused on an automatic domestic wood pellet boiler with a power of 30 kW Figure 1), equipped with a water-flue gas heat exchanger whose main role is to recover a part of the heat energy in the flue gas and transfer it in the water. The water circulation in the hydraulic circuit is ensured by a pump. The introduction of the pellets into the burner of the boiler is done by a screw that operates cyclically as long as the temperature of the outlet water is lower than the setpoint temperature. When the setpoint temperature is reached, the pellet supply stops. The primary air arrives through trapdoors located in the lower part of the boiler and its circulation is ensured by an exhaust fan mounted on the top cover of the boiler which is controlled by a lambda probe located in the chimney of the boiler. To dissipate the heat of the working fluid, the hydraulic circuit is connected to two air heaters located to outside of the test cell. In order to carry out an experimental characterization of the boiler, several sensors are installed at different locations in the boiler (Figure 1). An electromagnetic flowmeter with an operating range of 20 to 500 dm3/h with an uncertainty of 0.5% measures the water mass flow rate ( . mexp w ) circulating in the boiler heat exchanger. The flue gas mass flow rate ( . mexp fg ) is calculated from pressure and temperature measurement in the chimney (Pitot wing system connected to a micromanometer (uncertainty 5% and K-type thermocouple (uncertainty 0.75%) on their measurement ranges respectively). The water temperature at the inlet (Texp w,in) and outlet (Texp w,out) of the water-flue gas heat exchanger are recorded by two platinum Pt100 probes (uncertainty 0.8%). A K-type thermocouple is inserted at the chimney (uncertainty 0.75%) for the measurement of the flue gas temperature (Texp fg,exh). Another type S thermocouple (uncertainty 0.25%) is placed in the central axis of the combustion chamber to measure the instantaneous evolution of the flue gas temperature (Texp fg,cc).

K-type thermocouples (uncertainty 0.75%) are placed in the burner (Texp fg,bur), on the top and bottom sides of the combustion chamber (Texp fg,top and Texp fg,bot). Two other K-type thermocouples (uncertainty 0.75%) are also welded to each side of the combustion chamber wall (Texp wall,outer) and (Texp wall,inner). A K-type thermocouple (uncertainty 0.75%) is placed at the outlet of the heat exchanger tubes (Texp fg,exit) (Figure 1). The flue gas temperature measurements in the burner and the combustion chamber have been corrected from radiative effects. Indeed, with such temperature levels, the radiative dissipation of the thermocou-

ples is significant. Several methods exist to take into account this phenomenon which underestimate the true value of the temperature. The method used is the extrapolation method [35,36] which consists in using two thermocouples with wires of different diameters and therefore with different hot welds diameter (here 0.95 mm and 0.64 mm) placed at the same position. The radiative flux exchanged is assumed to be proportional to the surface of the hot weld, resulting in a zero radiative flux when the surface of the weld is infinitely small. From the two measured temperatures, an extrapolation allows to obtain the temperature value for a zero-weld surface corresponding to an absence of radiation. In our case study (in the flame and its vicinity), as an example, for a 1000 ◦C temperature measurement, the corrective value to be applied reach 170 ◦C. In this paper, the superscript "exp" corresponds to experimental measurements. The quantities calculated by the model have no superscript.

**Figure 1.** Thermocouples positions and hydraulic circuit. (1) Burner. (2) Combustion chamber. (3) Flue gas extraction.  (5) Air inlet. (6) Water-flue gas heat exchanger.

The flue gas temperature evolution in the combustion chamber measured at radius of 90 mm and height of 330 mm obtained during the boiler operating cycle is plotted versus time and correlated with the pellets mass flow rate (Figure 2). It shows a strong dependence between the quantities of pellet supplied by a feed screw and the temperature increase of the flue gas in the boiler combustion chamber. The burnt gas temperature varies between 600 and 1100 ◦C. It increases with the arrival of pellets and decreases with their complete consumption. As mentioned above, the pellets are introduced into the boiler burner by a feed screw that rotates with a PWM duty cycle as long as the water temperature at the outlet of the heat exchanger is lower than the set temperature. This operating mode is controlled by a pulse-width modulation control. When the set temperature is reached, the pellet supply stops. Thus, a long stop of the pellet supply (12 min) can be observed in Figure 2. The pellet supply disruptions induce a drop in the burnt gas temperature in the combustion chamber with a temporary delay. The time lapse between the increase and drop of flue gas temperature defines a thermal cycle. The pellets mass flow rate is deduced from the calibration of the feed screw, according to the angular position of the screw. All sensors are connected to a data acquisition system with a dedicated code developed under Labview software.

**Figure 2.** Boiler operating cycle.

### *2.2. Energy Balance of the Boiler*

Although estimated from the rotation speed of the screw during a calibration phase without combustion, the mass flow of pellets is difficult to obtain accurately over short operating times of the screw. In fact, the quantity of pellets introduced by the screw by PWM method is not constant between each cycle because more or less large pellet clusters are detached from the screw. Moreover, the combustion is not instantaneous, it would be necessary to introduce a dynamic combustion model of solid biomass to calculate the heat release as a function of time. The boiler model presented here can be modified in the future by integrating this combustion dynamics. The pellet combustion is therefore not modelled, so the heat generated during the combustion of the pellets has been calculated using the experimental mass flow rate and the experimental temperature of the flue gases in the burner. Then, the heat flux provided by the combustion of pellets is estimated from the following equation:

$$\overset{\bullet}{\mathbf{H}}\_{\mathrm{fg,bur}}(\mathbf{t}) = \overset{\bullet}{\mathbf{m}}\_{\mathrm{fg}}^{\mathrm{exp}}(\mathbf{t}) \cdot \left[ \mathbf{c}\_{\mathrm{P}} \left( \mathrm{T}\_{\mathrm{fg,bur}}^{\mathrm{exp}} \right) \cdot \left( \mathrm{T}\_{\mathrm{fg,bur}}^{\mathrm{exp}}(\mathbf{t}) - \mathrm{T}\_{\mathrm{fg,bur}}^{\mathrm{ref}} \right) + \Delta \mathrm{H}\_{\mathrm{ref}}^{0} \right] \tag{1}$$

With:.

Hfg,bur: heat flux released from pellet combustion (W) .

mexp fg : experimental flue gas mass flow rate (kg.s−1)

Texp fg,bur: experimental flue gas temperature in the burner (K)

Tref fg,bur: reference temperature for flue gas in the burner (298 K)

ΔH0 ref: standard formation enthalpy of gas in the burner (J.kg−1).

Considering the majority presence of N2 and O2 (air excess close to 1) in the mixture and for a first approximation, we assume that the mixture is composed as a gas including only pure species. We can therefore assume that ΔH0 ref = 0.

The energy balance of the boiler is established at each time step. It represents the heat exchanges between the flue gas and the boiler structure, the heat flux recovered by the water in the water-flue gas heat exchanger and the losses at the boiler exhaust. The outer wall of the boiler is assumed to be adiabatic because the boiler is very well insulated and the losses with the environment are negligible compared to the other heat flux. The losses are more significant from the boiler outlet through the exhaust pipe but this part is not modeled here. The heat flux released from the pellets combustion is then given by:

$$
\tilde{\mathbf{H}}\_{\text{fg,bar}}(\mathbf{t}) = \tilde{\mathbf{H}}\_{\text{fg,exh}}(\mathbf{t}) + \Delta \tilde{\mathbf{H}}\_{\text{w}}(\mathbf{t}) + \tilde{\mathbf{Q}}\_{\text{wall}}(\mathbf{t}) \tag{2}
$$

With:.

Hfg,exh: exhaust heat flux (W).

Δ . Hw: heat flux transferred to the water (W). .

Qwall: heat flux stored in the boiler structure (W).

The flue gases resulting from the combustion of the pellets go through the boiler subsystems (Figure 3 dashed red line) and exchange heat with their walls. Due to the transient phases, the walls store or yield a quantity of heat flux from or to the flue gases: the heat flux stored in the combustion chamber walls . Qwall,cc, in the inner wall of the heat exchanger . Qwall,HEx and in the walls of the flue gas tubes . Qwall,tub. The walls of the subsystems store some heat flux, consisting of three parts:

$$
\stackrel{\bullet}{\mathbf{Q}}\_{\text{wall}} = \stackrel{\bullet}{\mathbf{Q}}\_{\text{wall}, \text{cc}} + \stackrel{\bullet}{\mathbf{Q}}\_{\text{wall}, \text{HEx}} + \stackrel{\bullet}{\mathbf{Q}}\_{\text{wall}, \text{tub}} \tag{3}
$$

**Figure 3.** Energy balance of the boiler.

Each of these heat fluxes is calculated in the subsystems from inlet and outlet flux (Figure 3):

$$
\stackrel{\bullet}{\mathbf{Q}}\_{\text{wall},\infty} = \stackrel{\bullet}{\mathbf{Q}}\_1 - \stackrel{\bullet}{\mathbf{Q}}\_2 \tag{4}
$$

$$
\stackrel{\bullet}{\mathbf{Q}}\_{\text{wall}\downarrow\text{HEx}} = \stackrel{\bullet}{\mathbf{Q}}\_3 - \stackrel{\bullet}{\mathbf{Q}}\_4 \tag{5}
$$

$$
\stackrel{\bullet}{\mathbf{Q}}\_{\text{wall},\text{td}} = \stackrel{\bullet}{\mathbf{Q}}\_{6} - \stackrel{\bullet}{\mathbf{Q}}\_{5} \tag{6}
$$

The heat flux transferred to the inner walls of the heat exchanger and the flue gas tubes is partially transferred to the water.

$$\begin{array}{c} \stackrel{\bullet}{\rm \rm \mkern-1.0pt} \Delta \stackrel{\bullet}{\rm H}\_{\rm w} = \stackrel{\bullet}{\rm Q}\_{4} + \stackrel{\bullet}{\rm Q}\_{5} - \stackrel{\bullet}{\rm Q}\_{\rm w,st} \\ \stackrel{\bullet}{\rm \Delta \rm H}\_{\rm W} = \stackrel{\bullet}{\rm m}\_{\rm w}^{\rm exp} \cdot \left( \mathrm{c}\_{\rm W} \left( \mathrm{T}\_{\rm w,out}^{\rm exp} \right) \cdot \mathrm{T}\_{\rm w,out}^{\rm exp} - \mathrm{c}\_{\rm w} \left( \mathrm{T}\_{\rm w,in}^{\rm exp} \right) \cdot \mathrm{T}\_{\rm w,in}^{\rm exp} \right) \\ \stackrel{\bullet}{\rm \bullet \, H}\_{\rm W} = \stackrel{\bullet}{\rm H}\_{\rm w,out} - \stackrel{\bullet}{\rm H}\_{\rm w,in} \end{array} \tag{7}$$

With . Qw,st: heat flux stored by the water in the heat exchanger (W).

### **3. 0D Bond Graph Modeling**

The modeling of the main components of the boiler system, such as the combustion chamber, the flue gas tubes and the heat exchanger is performed using Bond Graph formalism. The boxes in Figure 4 represent the subsystems of the studied boiler, where the half-arrows characterize the thermal and hydraulic Bond Graph links between the subsystems. The word Bond Graph model describes here the thermal and mass transfers between subsystems. Causalities (**I**) are also present in order to indicate the variables at the origin of the system dynamics. In Figure 4, the combustion chamber box is not detailed, it includes the flue gas path from burner to the bottom of the heat exchanger. The combustion chamber temperature noted Tfg,cc is located inside this box but not appears in the Inlet/Outlet Bond Graph links.

**Figure 4.** Word Bond Graph model.

### *3.1. 0D Model of the Boiler*

The 0D dynamic model of the boiler is shown in Figure 5 with all the boiler subsystems (burner, combustion chamber, heat exchanger). The time variation of mass flow rate and temperature of both water and flue gas is considered. As input conditions, the experimental flue gas temperature Texp fg,bur in the burner and the experimental mass flow rate of the flue gas . mexp fg are introduced as time files. The water-flue gas heat exchanger is also modelled by providing the experimental water mass flow rate . mexp w and the experimental inlet water temperature Texp w,in as input conditions. The model is therefore stimulated with real limit

conditions and then validated with other experimental measurement obtained at the same time than these limit conditions values. This method improves the validation of the model.

**Figure 5.** 0D model of the boiler.

The storage and/or removal of thermal energy in the walls in the different zones of the boiler (refractory concrete, walls of the combustion chamber and walls of the heat exchanger) are modeled by the following law:

$$\stackrel{\bullet}{Q}\_{\text{wall}, \text{i}} = \mathbf{m}\_{\text{wall}, \text{i}} \times \mathbf{c}\_{\text{wall}, \text{i}} \times \frac{\text{d}\mathbf{T}\_{\text{wall}, \text{i}}}{\text{d}\mathbf{t}} = \sum \mathcal{Q}\_{\text{diss}} \tag{8}$$

With:.

Qwall,i: heat flux stored in the wall of the system i (W).

∅diss : dissipative fluxes between wall and flue gas (W).

mwall,i: wall mass of the system i (kg).

cwall,i: wall specific heat of the system i (J/kg−1.K−1).

This expression is traduced to a 'C' element (C1, C7, C3, C10) (Figure 5) in the bond graph formalism because the flux is a function of the derivative of the effort:

$$\frac{r\_{out}}{Q\_{out}} \cdot C$$

The causality applied to these 'C' elements is always a flux causality because at the beginning of the simulation, it is the temperatures (efforts) that are known and then allow calculation of dissipative heat fluxes (conductive, convective and radiative ones). So, in these differential equations solving, the value of temperatures of the next time step are obtained from the integration of the flux balance at the current time.

$$\mathbf{T}\_{\text{wall},i}(\mathbf{t}) = \frac{1}{\mathbf{m}\_{\text{wall},i} \times \mathbf{C}\_{\text{wall},i}} \int\_0^\mathbf{t} \overset{\bullet}{\mathbf{Q}}\_{\text{wall},i} \mathbf{dt} \tag{9}$$

where . Qwall,i is calculated from a heat flow balance between all dissipative fluxes ∅diss (Equation (8)) and therefore obtained with a '0junction centered in the wall.

The dynamic behavior of the boiler depends on the interaction between the both hydraulic and thermal systems. The 'RS' elements (R11, R12, R16, R18 and R23) have been used to couple them in order to calculate enthalpy flux (10) from inlet temperature and mass flow rate. The power input of each RS elements is defined with an effort causality which means that the temperature value (effort) of thermal power input is knows at the start of each calculation step. RS elements then calculate the enthalpy flux, which is necessary for each flux balance carried out by the zero junctions on the path of the flue gas labelled by the dashed red line. The enthalpy flux of flue gas at the inlet and outlet of each 'RS' element . Hfg,in and . Hfg,out (Figure 6) is calculated with the following Equations (10) and (11) and the same hypothesis than the Equation (1):

$$\overset{\bullet}{\mathbf{H}}\_{\text{fg,in}}^{\bullet}(\mathbf{t}) = \overset{\bullet}{\mathbf{m}\_{\text{fg}}}^{\exp}(\mathbf{t}) \cdot \left[ \mathbf{c}\_{\text{P}} \left( \mathbf{T}\_{\text{fg,in}} \right) \cdot \left( \mathbf{T}\_{\text{fg,in}}(\mathbf{t}) - \mathbf{T}\_{\text{fg}}^{\text{ref}} \right) + \Delta \mathbf{H}\_{\text{ref}}^{0} \right] \tag{10}$$

$$\stackrel{\bullet}{\mathbf{H}}\_{\text{fg,out}}(t) = \stackrel{\bullet}{\mathbf{H}}\_{\text{fg,in}}(t) \tag{11}$$

With:

**Figure 6.** RS element example.

. Hfg,in(t): calculated flue gas enthalpy flux at the RS-element inlet (W), equal to . Hfg,bur .

Hfg,out(t): calculated flue gas enthalpy flux at the RS-element outlet (W), equal to . Hfg,bur

Texp fg,bur(t): experimental temperature in the burner (K)

Tfg,cc(t): calculated flue gas temperature (K) imported from the following '0 junction in the combustion chamber. .

mexp fg (t): experimental flue gas mass flow rate (kg.s−1).

Pfg: experimental pressure in the boiler (Pa) supposed constant because pressure losses are low and not easy to model in 0D due to the complexity of the geometry.

The Figure 6 illustrates one of the 'RS' elements. This one is located between the burner outlet and the combustion chamber.

The thermal transfers by conduction and convection are modeled by equations of the following form:

$$\stackrel{\bullet}{Q}\_{\rm diss}(\mathbf{t}) = \frac{1}{\mathcal{R}\_{\rm th}} . \,\mathrm{AT}(\mathbf{t}) \tag{12}$$

With:.

Qdiss(t): dissipative flux (W).

R*th*: thermal resistor (K.W−1).

ΔT(*t*): temperature difference (K).

Here, the flux . Qdiss(t) directly depends on the effort ΔT(*t*) (not derivative link) whether it's linear or not. In this case the 'R' element is used in the bond graph formalism:

$$\xrightarrow[\partial\_{\rm diss}]{\mathfrak{A}\mathcal{T}} \mathcal{R}$$

ΔT(*t*) is obtained with a '1' junction which consist in an effort balance and then can calculate the temperature difference.

R4, R5, R21, R22, R8, R9 and R25 quantify the conductive exchanges through the walls. In a cylinder, the thermal conductive resistance Rcd is given by:

$$\mathbf{R}\_{\rm cd,i} = \frac{\ln\left(\frac{\mathbf{r}\_{2,i}}{\mathbf{r}\_{1,i}}\right)}{2\pi\lambda\_{\rm i}\mathbf{H}\_{\rm i}}\tag{13}$$

With:

r2,i: outside radius of the system i (m).

r1,*i*: inside radius of the system i (m).

λi: thermal conductivity of the system i (W.m−1.K−1).

Hi: height for the system i (m).

R3, R6, R20, R10, and R2 deal with the convective transfers between the flue gas and the different walls of the boiler. They are calculated from the convective resistance Rcv:

$$\mathbf{R}\_{\rm cv,i} = \frac{1}{\mathbf{h}\_{\rm g,i} S\_{\rm i}} \tag{14}$$

With:

hg,i: global thermal transfer coefficient of the system i (W.m−2.K−1).

Si: exchange surface of the system i (m2).

The global thermal transfer coefficients hg,i, including convective and radiative effects, for the different geometrical configurations in the boiler are obtained according to a first stage of simulation by inverse method. Indeed, the radiative effects of the flame or the gases with the walls are complex to model in 0D. This method is based on the energy balances presented in Section 2.2 for each zone. Heat fluxes are calculated by using experimental wall and flue gas temperatures as well as experimental water temperatures (and calculated temperatures by the dynamic 0D model when the experimental measurement is not available). These experimental temperature values are introduced into the model at each calculation time step. They thus allow at each time step to calculate the value of hg,i as illustrated in the following relations (Equation (15)) in order to use it for the calculation of the parietal fluxes in the model. This method allows to obtain temporal evolutions of hg,i coefficients like the one presented in Figure 7 and was implemented only once as a prerequisite to the main simulation. This allowed to determine the global coefficients even in areas where we were unable to place thermocouples probes by using temperatures calculated as close to reality as possible since in places where temperatures were measured, the model took them into account at each time step. This combination of measured and

calculated quantities inside a behavioral model is similar to a HIL (Hardware in the Loop) process.

**Figure 7.** Global thermal transfer coefficient inside the combustion chamber.

By example, the variation of the flue gas enthalpy flux between the inlet and the outlet of the combustion chamber Δ . Hfg,cc(t) is calculated at each time step. The global heat flux exchanged between the flue gas and the combustion chamber wall . Qfg,cc(t) is calculated at each time step also. By performing a balance between the two heat fluxes, the value of global thermal transfer coefficient is deduced for each time step (Equation (15)). The time evolution of a global thermal transfer coefficient, including radiative and convective transfers, near the inside combustion chamber is shown in Figure 7. The peaks observed are due to the low temperature difference between the flue gas and the wall of the combustion chamber. With the induced errors on the global thermal transfer coefficients higher than 2000 W.m−2.K−1), they cause discrepancies in the calculation carried out by the dynamic model. This adds complexity to the choice of the computation scheme.

$$\begin{aligned} \boldsymbol{\Lambda} \stackrel{\bullet}{\mathbf{H}}\_{\mathrm{fg,cc}}(\mathbf{t}) &= \stackrel{\bullet}{\mathbf{m}}\_{\mathrm{fg}}^{\mathrm{exp}}(\mathbf{t}) \cdot \left( \underbrace{\mathbf{c}\_{\mathrm{p,fg}} \left( \mathbf{T}\_{\mathrm{fg,cc}}^{\mathrm{exp}}(\mathbf{t}) \right) \mathbf{T}\_{\mathrm{fg,cc}}^{\mathrm{exp}}(\mathbf{t})}\_{\mathrm{cc} \text{ outlet}} - \underbrace{\mathbf{c}\_{\mathrm{p,fg}} \left( \mathbf{T}\_{\mathrm{fg,bur}}^{\mathrm{exp}}(\mathbf{t}) \right) \mathbf{T}\_{\mathrm{fg,bur}}^{\mathrm{exp}}(\mathbf{t})}\_{\mathrm{cc} \text{ inlet}} \right) \\ & \stackrel{\bullet}{\mathbf{Q}}\_{\mathrm{fg,cc}}(\mathbf{t}) = \mathbf{h}\_{\mathrm{g,cc}} \mathrm{S}\_{\mathrm{cc}} \left( \mathbf{T}\_{\mathrm{fg,cc}}^{\mathrm{exp}}(\mathbf{t}) - \mathbf{T}\_{\mathrm{wall,cc}}^{\mathrm{exp}}(\mathbf{t}) \right) \\ & \stackrel{\bullet}{\mathbf{h}}\_{\mathrm{g,cc}}(\mathbf{t}) = \frac{\mathbf{A} \mathbf{H}\_{\mathrm{fg,cc}}(\mathbf{t})}{\mathbf{S}\_{\mathrm{cc}} \left( \mathbf{T}\_{\mathrm{fg,c}}^{\mathrm{exp}}(\mathbf{t}) - \mathbf{T}\_{\mathrm{wall,cc}}^{\mathrm{exp}}(\mathbf{t}) \right)} \end{aligned} \tag{15}$$

With:

Δ . Hfg,cc: variation of the flue gas enthalpy flux between the inlet and the outlet of the combustion chamber (W) .

Qfg,cc: global heat flux exchanged between the flue gas and the combustion chamber wall (W)

cpfg: flue gas specific heat at constant pressure (J.kg−1.K−1).

Scc: combustion chamber exchange surface (m2).

hg,cc: global thermal transfer coefficient for the inner wall of the combustion chamber (W.m−2.K−1).

Texp wall,cc: experimental temperature of the inner wall of the combustion chamber (K).

For this example, apart from the peaks mentioned above, an average value of hg = 200 W.m−2.K−1 has been selected. The results presented in Figures 8 and 9 show that this approximation leads to some errors in the calculated water and flue gas temperatures. Indeed, we could identify here two operating regimes:

**Figure 8.** Comparison of experimental and calculated water temperatures at the outlet of the waterflue gas heat exchanger.

**Figure 9.** Comparison of experimental and calculated flue gas temperature at the outlet of the flue gas tubes.

These two operating regimes can be identified in Figure 10, the combustion is continuous until *t* = 170 min and then an operating cycle is set up as presented in Figure 2. The choice of only one value for the global coefficient generates an under estimation of the transfers on the first phase and an over estimation on the second one as it can be noticed in Figures 8 and 9.

**Figure 10.** Comparison of experimental and calculated flue gas temperatures in the combustion chamber.

Nevertheless, the change of operating regime is difficult to take into account here in an automated way without modeling the pellet supply mechanism and their subsequent combustion.

In order to differentiate radiative and convective heat exchanges in the boiler, the convective heat fluxes are calculated using the Newton's law and the convective coefficients from Equation (16). Knowing Reynolds number as well as Prandtl, Nusselt numbers was determined from the semi-empirical correlations of Dittus-Boelter [37] and Gnielinski [38], adapted to the studied configurations. Nusselt number then allows to calculate the convective exchange coefficient within the geometric configurations remaining inside the boiler. Table 2 includes the semi-empirical correlations used to calculate Nusselt number.

$$\text{Nu} = \frac{\text{h.D}\_{\text{h}}}{\lambda\_{\text{fg}}} \tag{16}$$

With:

Dh: hydraulic diameter (m).

λfg: flue gas thermal conductivity (W.m−1.K−1).

h: convective coefficient (W.m−2.K−1).

**Table 2.** Semi-empirical correlations used for the calculation of Nusselt number.

Finally, the calculation of the water temperature at each instant is obtained by a flux balance represented by the area inside the blue dotted lines in Figure 5.

The flux balance consists of the algebraic sum of the enthalpy input/output fluxes calculated by the 'RS' elements (respectively 'R14' and 'R15') with the convective heat fluxes calculated with the thermal resistances 'R2', 'R7' and 'R19'. The water temperature is then calculated from the integration of the flux balance performed in the 'C2' element.

### *3.2. Flue Gas Thermodynamic Properties*

To take into account the variation of the thermodynamic properties of the flue gas in the boiler, correlations have been used for each property as a function of the temperature.

In this section, all the correlations used to calculate the thermodynamic properties are detailed. The properties are: density ρfg from the perfect gas law with R = 8.314 J/mol−1/K−1, thermal conductivity λfg, dynamic viscosity μfg, and specific heat cpfg of the flue gas resulting from the combustion of the pellets. The mass fraction used to calculate some properties is evaluated from the molar fraction deduced from Equation (17). The correlation versus temperature of cpfg is introduced in the 0D model to calculate enthalpy flux. The other thermodynamic properties correlation are used to calculate Reynolds and Prandtl numbers in order to deduce the convective transfer coefficients introduced above.

The thermodynamic properties of flue gas are obtained by adding the properties of each species multiplied by the corresponding molar or specific fractions. The mixture of these species is given by the combustion reaction.

From an elementary analysis, the chemical formulation of pellet is C36.725H71.6O30.475. Then, the combustion reaction of pellets in air is given as below:

$$\begin{array}{c} \text{C}\_{36.725}\text{H}\_{71.6}\text{O}\_{30.475} + 39.388(1 + \text{e}\_{\text{air}})(\text{O}\_{2} + 3.76\text{N}\_{2})\\ \rightarrow 36.725\text{CO}\_{2} + 35.8\text{H}\_{2}\text{O} + 39.388.\text{e}\_{\text{air}}\text{O}\_{2} + 148.097(1 + \text{e}\_{\text{air}})\text{N}\_{2} \end{array} \tag{17}$$

With eair the air excess.

As discussed before, it is assumed that the pressure in the boiler remains constant and equal to the atmospheric pressure and the specific fractions of the combustion products also remain constant during the boiler cycles and in the different zones. Knowing that the boiler operates with an air excess of 80% (eair = 0.8), the correlations used are presented in Table 3.

**Table 3.** Correlations for the calculation of the flue gas thermodynamic properties (i = CO2, H2O, O2, N2)**.** The constants A, B, C, D and E were fixed for each species and for each property.

*3.3. Solver Scheme*

The Bond Graph method uses a system of algebraic-differential equations to describe the dynamic of the modeled system. The accuracy of the dynamic model is based on the choice of the computation scheme used to efficiently solve these differential equations. The resolution scheme used in our model is the Runge-Kutta fourth order formula (RK4) [42,43] which is a particulate case of Runge-Kutta method. This method is recommended when the required accuracy is very high but it requires more CPU time than simpler methods (for this study about 5 min). This method is based on the iteration principle, i.e., an estimation of the solution is calculated from the previous solution. The principle is to approach the next value yn+1 at time tn+1 by the current value yn obtained at time tn combined with a function taking into account the iteration step (δ) and the estimated slope. The slope is obtained by the weighted average of four slopes (k1, k2, k3 and k4), where each slope is the product of the iteration step and an estimated slope. The slope is specified by the function F on the right side of the differential equation [44,45].

The following problem is then considered:

$$\mathbf{y} = \mathbf{F}(\mathbf{t}, \mathbf{y}) \text{ with } \mathbf{y}\_0 = \mathbf{f}(\mathbf{t}\_0) \to \mathbf{y} = \mathbf{f}(\mathbf{t}) \tag{18}$$

From a known initial condition, the RK4 method is given by the equation:

$$y\_{n+1} = y\_n + \frac{\delta}{6}(k\_1 + 2k\_2 + 2k\_3 + k\_4) + \left(\delta^5\right) \tag{19}$$

$$
\delta = \mathbf{t}\_{\mathbf{n}+1} - \mathbf{t}\_{\mathbf{n}\prime} \ 1 < \mathbf{n} < \mathbf{N} \tag{20}
$$

where

$$\mathbf{k}\_1 = \mathbf{F}(\mathbf{t}\_{\mathbf{n}}, \mathbf{y}\_{\mathbf{n}}) \tag{21}$$

$$\mathbf{k}\_2 = \mathbf{F}\left(\mathbf{t}\_n + \frac{\boldsymbol{\delta}}{2}, \mathbf{y}\_n + \frac{\mathbf{k}\_1}{2}\right) \tag{22}$$

$$\mathbf{k}\_3 = \mathbf{F}\left(\mathbf{t}\_\mathbf{n} + \frac{\delta}{2}, \mathbf{y}\_\mathbf{n} + \frac{\mathbf{k}\_2}{2}\right) \tag{23}$$

$$\mathbf{k}\_4 = \mathbf{F}(\mathbf{t}\_n + \delta, \mathbf{y}\_n + \mathbf{k}\_3) \tag{24}$$

The RK4 method is of order 4, this means the error committed at each step is of the order of δ5, whereas the total accumulated error is of the order of δ4.

### **4. Results and Discussion**

The experimental results are discussed to characterize the boiler operation. They are then compared to the simulations. In order to validate the model, measurements of the flue gas temperature profiles in the combustion chamber and at the heat exchanger outlet as well as measurements of water temperature at boiler outlet are carried out.

The dynamic model input data are the experimental flue gas temperature in the burner (Texp fg,bur), the experimental water temperature at the boiler inlet (Texp w,in), the experimental flue gas mass flow rate ( . mexp fg ) and the experimental water mass flow rate ( . mexp w ).

The thermal behaviors of the flue gas in the boiler and the water in the heat exchanger are investigated.

From the analysis of the flue gas temperature profiles in the combustion chamber (Figure 10), sudden and fast temperature changes occur during the boiler start-up due to the uncontrolled combustion of a large mass of pellets during this step. Before the combustion start, pellets are heated during several minutes (about 15 min) with an air heater. During the entire control phase (regulation phase between 30–270 min), the flue gas temperatures remain very high and display fluctuations. Then, they decrease progressively during the cooling phase (270–355 min). The fast fluctuations observed are due to the quantities of pellets supplied every 20 s. These fluctuations are also observed on the temperature profiles calculated from the 0D model as a consequence of the limit condition that is an experimental measurement of the flue gas temperature in the burner. We can note that a notable difference exists between the calculated temperature and the one measured during the beginning of the combustion phase (just after the start-up jump). This difference is undoubtedly linked to the fact that the quantity of pellets burning in this phase is very important (accumulation before combustion start), also the gasification is such that the combustion continues in the upper combustion chamber (above the burner zone). The model only integrates the combustion in the burner and therefore does not integrate this excess of heat release in the combustion chamber area.

This problem would have been the same using a combustion model of the pellets in the burner area. It would be necessary to separate the combustion in the 2 zones (burner and combustion chamber) and thus to find a key of distribution of the combustible gases in each zone. This key is not easy to find because the problem is related to unsteady 3D aero thermochemical phenomena.

The water temperature at the heat exchanger outlet is also examined. As noted for the flue gas temperature, a drastic increase of the water temperature can be observed during the start-up phase. During this phase (0–30 min), the water circulates in closed circuit until to reach a temperature of 325 K (Figure 8). This process is imposed by the mixing valve (3-way valve) resulting in a significant increase of the water temperature. After this step, the hot water is redirected to the cooling circuit. The fluctuations observed during the cooling phase are due to the intermittent operation of the water pump to maintain as long as possible the boiler body closed to the operating conditions if the boiler needs to be restarted.

Here, the differences between the calculated and measured values are not significant. Differences of 5 ◦C are nevertheless noted in the cyclic operation zone (170–250 min), this is perhaps linked to the overestimation of the global exchange coefficients in this operation mode as mentioned at the end of Section 3.1.

The instantaneous evolutions of the experimental and calculated flue gas temperature at the outlet of the flue gas tube of the water heat exchanger are plotted in Figure 10. The flue gas temperature at the flue gas tube outlet has the same evolution as in the combustion chamber. Nevertheless fluctuations are filtered by the thermal inertia of the different parts of the boiler along the flue gas path. The temperature of the flue gas remains relatively high at outlet of the tubes (~573 K).

As mentioned at the end of Section 3.1, here the under estimation of the wall global exchange coefficients on the stabilized phase and the under estimation on the cyclic phase is notable. In the stabilized phase, the underestimation of the thermal wall fluxes limits the thermal dissipation of the gases and thus also the reduction of their temperature. In the cyclic phase, the overestimation of the fluxes increases abnormally the wall transfers and reduces the flue gas temperature. Significant temperature differences remain during the cooling phase of the boiler.

The time evolution of the water enthalpy flux variation between inlet and outlet, calculated by the dynamic model, is plotted in Figure 11. The heat flux drops to zero when the pump stops.

**Figure 11.** Variation of water enthalpy flux.

The instantaneous thermal power delivered to the water can be calculated as:

$$
\Delta \overset{\bullet}{\mathbf{H}}\_{\text{w}}(\mathbf{t}) = \overset{\bullet}{\mathbf{m}}\_{\text{w}}^{\text{exp}} \left( \mathbf{c}\_{\text{w}} (\mathbf{T}\_{\text{w,out}}), \mathbf{T}\_{\text{w,out}} - \mathbf{c}\_{\text{w}} \left( \mathbf{T}\_{\text{w,in}}^{\text{exp}} \right), \mathbf{T}\_{\text{w,in}}^{\text{exp}} \right) \tag{25}
$$

With:

Δ . Hw(t): water enthalpy flux variation (W).

Texp w,in : experimental water temperature at the heat exchanger inlet (K).

Tw,out: water temperature at the heat exchanger outlet calculated by the dynamic model (K).

cw: water specific heat at an average temperature of 328 K (4183/kg−1.K−1).

The amount of heat flux transmitted to the water remains very low during the startup phase of the boiler and then increases drastically after the start of the combustion. After the start of combustion, an increase of the heat transmitted to the water can be observed between 30 and 50 min in Figure 11. This progression exists because some heat from combustion is first accumulated by the metal walls inside the boiler before being completely transferred to the water when the walls reach an established thermal regime. Considering only the regulation phase represented by the period from the combustion start time (30 min) to the boiler shutdown (270 min), the heat flux transmitted from the flue gas to the water of the heat exchanger is quite stable and close to 37 kW. The fluctuations observed on the flue gas temperature during the cyclic phase are well absorbed by the inertia of the walls and the water. After the boiler shutdown and during the cooling phase, intermittent operation of the pump is observed. When the pump is shut down, the walls of the heat exchanger transmit heat to the volume of water became motionless in the heat exchanger, which explains the peaks of enthalpy flux as soon as the pump is started up again.

From the dynamic model of the boiler, it is also possible to calculate its efficiency, lost power, heat flux stored by the walls and released by combustion in the burner.

According to the manufacturer, the boiler must have an average efficiency of 85% under nominal operating conditions (thermal power of 30 kW). In this study, the duty cycle (PWM) of the pellet supply screw was modified to increase the power of the boiler in order to saturate the downstream thermal load and thus create thermal control cycles suitable to unsteady operating conditions.

The efficiency can be defined as the ratio of the enthalpy flux variation of the water heat exchanger and the heat flux released from pellet combustion:

$$\eta(\mathbf{t}) = \frac{\Delta \dot{\mathbf{H}}\_{\mathbf{w}}(\mathbf{t})}{\dot{\mathbf{H}}\_{\mathbf{fg,bur}}(\mathbf{t})} \tag{26}$$

The instantaneous evolution of the efficiency calculated by the model at each moment t is presented in Figure 12. Its evolution is drastically affected by the fluctuations of flue gas temperature. The boiler has an average efficiency of 67.5%. This low efficiency value is not surprising because the overpower generated in our test case cannot be fully absorbed by the capacity of the gas-water exchanger of the boiler. The thermal power of the boiler is nevertheless increased (37 kW instead of 30 kW).

Several parameters impacting the response of the 0D model can be highlighted. For example, an influence on the thermal behavior of the flue gas with the mass flow rate can be distinguished at the tube outlet. As the experimental flue gas mass flow rate is used as input condition, during the cooling phase a significant discrepancy between the evolution of the calculated and experimental flue gas temperatures is recorded as shown in Figure 13a.

**Figure 12.** Calculated boiler efficiency.

**Figure 13.** Influence of the flue gas mass flow rate at the flue gas outlet on the 0D model response. (**a**): With experimental flue gas mass flow rate, (**b**): With a flue gas mass flow corrected for each boiler operating phase.

This discrepancy can be explained by a low value of flue gas flow rates during the cooling phase (Figure 13a) according to the accuracy of the measurement chain (Pitot tube associated with a micromanometer and a thermocouple), which involves a maximum error of 40%. By adjusting the mass flow rate value of the flue gas, staying within the uncertainty range of the flowmeter, a clear improvement of the model response is observed (Figure 13b).

A heat flux balance in four zones of the boiler is performed by using the equations introduced in Section 2.2 and allows to compare the radiative and convective heat flux. The convective heat flux is calculated by using the convective coefficients obtained from the semi-empirical correlations given in Table 2 and the total heat flux φtot is obtained by using the global thermal transfer coefficients calculated from the inverse method by carrying out flux balances inside the boiler using Equations (3)–(6), previously introduced in Section 2.2. The radiative flux can be then deduced from the total heat flux, knowing the convective heat flux.

Due to the presence of large temperature gradients in the boiler, combustion products such as water vapor (H2O), carbon dioxide (CO2) and soot particles radiate significantly. Radiation is the dominant thermal transfer in the boiler and must be compared to the total thermal transfers (Table 4).

**Table 4.** Heat flux balance.

### **5. Conclusions**

A 0D dynamic modeling of a domestic biomass boiler of low power was developed by using Bond Graph formalism that allows to represent the coupled multi-physical phenomena, to study the thermal transfers between the different fluids during the transient operating phases, to evaluate the energy performances of the boiler and to take into account the variability of the heat production. The local evolution of the state variables is much less detailed than with CFD modeling but the dependencies of one zone with another are better taken into account with a 0D dynamic modeling. A biomass combustion model was not developed in this study but the combustion reaction of pellets in air allowed to calculate the thermodynamic properties of the flue gas in the boiler used in the 0D model. This model based on mass and energy balances was validated with experimental results, in particular the flue gas temperature in several locations of the boiler and the water temperature at the heat exchanger outlet. Some experimental data and 0D modeling at each time step of the calculation were coupled. The thermal transfers between the flue gas and the water circulating inside the heat exchanger and between these two fluids and the boiler structures were simulated. The experimental results showed a dependence of the evolution of the flue gas temperature in the combustion chamber as a function of the quantity of pellets supplied, according to the thermal cycle of the boiler. This directly affects the operating conditions of the boiler and generates important temperature fluctuations in the combustion chamber, which could significantly affect the operation of a hot air machine in the case of a conversion into a micro cogeneration unit. Indeed in this case, the air-gas exchanger of such an installation would be located in the zone where the temperature is the highest and thus closest to the flame. A calculation of the global thermal transfer coefficients by inverse method was done in the subsystems of the boiler. A good agreement between the experimental measurements and the simulation has been found and the origins of the differences have been identified, such as the excess of heat release in the combustion chamber above the burner zone not integrated in the model. It has been shown that the boiler has an average efficiency of 67.5% and the radiation is the dominant thermal transfer in the boiler by reaching 97.6% of the total thermal transfers inside the combustion chamber. The

0D dynamic model of the boiler during the operating phases allows not only to evaluate its energy performances but also to highlight the boiler components where the thermal transfers must be optimized.

The modeling of pellet combustion using a heat release law adapted to solid biomass combustion associated with an efficient identification of the pellet mass flow rate will make it possible to improve this model and make it independent of experimental boundary conditions. The radiative transfers being preponderant but difficult to model in 0D for mutual exchanges between gases, particles and walls, a detection of the different combustion phases and this according to the presence or not of flame in each of the zones would allow to better parameterize the global exchange coefficients which moreover will be able to be identified with the help of the proposed inverse method. The modeling methodology developed will allow the study of a complex unit, such as a CHP plant by coupling the different models for each component.

**Author Contributions:** Conceptualization, E.D.; Formal analysis, F.M.; Investigation, J.S.; Writing— Original Draft, F.M., C.M. and E.D.; Writing—Review and Editing, C.M. and E.D.; Software, F.M. and E.D.; Methodology, E.D. and C.M.; Validation, all authors; Resources: J.S.; Data curation, J.S.; Visualization, F.M.; Supervision, C.M.; Project administration, C.M.; Funding acquisition, C.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** F.M. was supported by joint PhD grant from ADEME (French Environment and Energy Management Agency) and Hauts-de-France Region. Grant Number: TEZ15-27.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors acknowledge CCM (ULCO, Dunkerque France) for the elementary analysis of wood pellets.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Nomenclature**

### **References**


**1. Introduction**

Recent years have brought attention to ferroptosis, an iron-dependent form of regulated cell death implicated in a broad range of diseases [1,2], and to the selective peroxidation of eicosatetraenoyl-phosphatidylethanolamine (ETE-PE) by 15-lipoxygenase (15LOX) complexed with PE-binding protein 1 (PEBP1) [3–5] as a ferroptotic mechanism. ETE-PE

**Citation:** Mikulska-Ruminska, K.; Anthonymuthu, T.S.; Levkina, A.; Shrivastava, I.H.; Kapralov, A.A.; Bayır, H.; Kagan, V.E.; Bahar, I. NO• Represses the Oxygenation of Arachidonoyl PE by 15LOX/PEBP1: Mechanism and Role in Ferroptosis. *Int. J. Mol. Sci.* **2021**, *22*, 5253. https://doi.org/10.3390/ijms22105253 **\***

†

Academic Editor: Masoud Jelokhani-Niaraki

Received: 3 April 2021 Accepted: 12 May 2021 Published: 17 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

peroxidation produces hydroperoxy-ETE-PE (HpETE-PE) molecules that, in turn, serve as pro-ferroptotic signals. 15LOX catalytic action includes the Fe-driven abstraction of a hydrogen atom from the carbon in bis-allylic position, formation of a carbon-centered radical followed by addition of molecular oxygen, culminating in the production of 15-HpETE.

The role of 15LOX in the onset of ferroptosis is now well-established [3–6]. Under physiological conditions, arachidonic acid (AA), a polyunsaturated fatty acid (PUFA) also termed *cis*-5,8,11,14-ETE or shortly ETE, is the preferred substrate of these non-heme iron-containing enzymes; but under inflammatory conditions, such as those stimulated by interleukin-13/4, the formation of the complex 15LOX/PEBP1 shifts the substratespecificity of 15LOX from free AA to AA esterified into PE (or ETE-PE). The AA-esterification of PE mostly occurs at the sn-2 acyl chain leaving the sn-1 chain for saturated or monounsaturated acyls, predominantly palmitic (C16:0), stearic (S) (C18:0), and oleic (C18:1). Of interest below will be 1-SA-2-ETE-PE, i.e., PE with stearic acid (SA) and arachidonic acid (AA or ETE) at the respective sn-1 and sn-2 chains, which will be shortly called SAPE. 15LOX-catalyzed peroxidation of SAPE at its sn-2 chain results in the production of 15- HpHETE-PE [3]—one of the major players in the ensuing induction of ferroptosis. As such, the 15LOX-2/PEBP1 complex emerged as a master promoter of ferroptotic cell-death [3]. This involvement of 15LOX-2/PEBP1 has been supported by the ability of ferrostatin-1, the most common ferroptosis inhibitor, to bind the complex [4] and inhibit the formation of 15-HpHETE-PE.

As a multistage process, ferroptosis is regulated by several enzymatic events occurring independently or concertedly: (i) biosynthesis of PUFA-PE peroxidation substrates by acyl Co-A synthetase 4 (ACSL4) and lysophosphatidylcholine acyltransferase (LPCAT) [7]; (ii) prevention of the formation of 15-HpETE-PE by flavin-containing oxidoreductase, ferroptosis suppressor protein 1 (FSP1) [8,9]; and (iii) reduction in 15-HpETE-PE into an innocuous alcohol, 15-HETE-PE, by glutathione peroxidase 4 (GPX4) [10,11]. Lately, this list of ferroptosis regulators has been extended to include new members, such as the Ca2+-independent phospholipase A2β (iPLA2β) [12], which can hydrolyze peroxidized phospholipids and, therefore, can eliminate the ferroptotic cell death signal, and the inducible isoform of nitric oxide synthase, iNOS, that generates NO• [13]. Notably, NO•'s reactivity towards different radical intermediates of lipid peroxidation has been discovered more than three decades ago [14]. Furthermore, NO• inhibits lipid peroxidation by lipoxygenase and cyclooxygenase [15].

In the present study, we first performed a lipidomics analysis to observe that NO• indeed suppresses the production of HpETE-PE by 15LOX-2/PEBP1, while promoting the formation of nitrosylated ETE-PE. As the 15LOX-2/PEBP1 complex, rather than 15LOX-2 alone, catalyzes the oxidation of 15-ETE-PE (rather than free ETE), we focused here on the interactions of NO• with SAPE and15LOX-2/PEBP1, and on its competition with the lipoxygenase substrate O2. Notably, the precise delivery of molecular oxygen to form a peroxyl radical from the carbon-centered radical of ETE at the 15LOX-2 catalytic site occurs via a specialized channel, whose structure and accessibility is altered upon binding of ETE to the enzyme [16]. With this in mind, we studied the effects of NO• peroxidation of SAPE by 15LOX-2/PEBP1 and the structural dynamics of the oxygen delivery process to the catalytic Fe3+ in the 15LOX-2/PEBP1 complex and the competition with NO•, in the presence and absence of the ETE-PE substrate compared to free ETE in 15LOX-2. Our study provides a mechanistic description of the time-evolved interactions between 15LOX-2 (bound to PEBP1) and its ligands, the ETE-PE- and PEBP1-induced changes in 15LOX-2 conformation or catalytic site accessibility which alter the delivery of O2 and NO•, and the competition between O2 and NO• for binding the active site. Several key residues mediating the events prior to peroxidation and the interference of NO• to repress them are identified, the critical roles of which are consolidated by in silico saturation mutagenesis experiments and comparisons with recent work.


**2. Contributions**

The exceptional properties of various metallic polycrystalline systems, such as high strength, excessive hardness, great ductility at room temperature, superior energy absorption capacity, and good corrosion resistance, make them outstanding candidates for a wide variety of applications where one or another of the mentioned qualities, or the combination of several, is of crucial importance. The presented contributions [1–9] evidently demonstrate that the properties of metallic materials are microstructure-dependent and, therefore, the thermomechanical processing (TMP) of the polycrystalline aggregates should be strictly controlled to guarantee the attainment of the desired suite of qualities. Given this, the assessment of microstructure evolution in metallic systems is of extraordinary importance. Since the trial–error approach is a time-consuming and rather expensive methodology, the materials research community tends to employ a wide spectrum of computational approaches to simulate each chain of TMP and tune the processing variables to ensure the necessary microstructural state for achieving the desired performance in the final product.

In the most general case, the TMP of metals involves a sequence of deformation and annealing processes (or the combination of both into one technological step) that lead to morphological changes in the polycrystalline aggregate, the evolution of crystallographic texture, or phase transformation. All these aspects of microstructure development are partially discussed in the presented contributions [1–9]. In the first manuscript [1], the deformation flow in a single-phase alloy with a face-centered crystal structure is analyzed

**Citation:** Sidor, J.J. Modelling the Deformation, Recrystallization, and Microstructure-Related Properties in Metals. *Metals* **2021**, *11*, 1759. https://doi.org/10.3390/met11111759

Received: 20 October 2021 Accepted: 29 October 2021 Published: 2 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

by means of a finite element (FEM) approach and computationally efficient flow-line model (FLM), which enabled the effective simulation of cold rolling in terms of the evolution of crystallographic texture. It was shown that by using the flow-line model, coupled with a Taylor-type homogenization crystal plasticity model, it is possible to carry out a texture simulation that is close to the one performed with the deformation history obtained by means of computationally expensive FEM. In addition, the correlation between the FLM model parameters and the rolling process quantitative indicators was defined, implying that the employed analytical approach was efficient in performing simulations of cold rolling without fitting constraints.

The second contribution [2] likewise addresses the crystallographic aspect of deformation by employing a well-established viscoplastic, self-consistent (VPSC) crystal plasticity approach. The numerical investigation is related to the deformation mechanisms in a magnesium alloy with a hexagonal crystal structure, where the polycrystalline system was subjected to uniaxial tension and compression. The influence of slip systems on the plastic deformation, mechanical response, and texture evolution was analyzed by the VPSC via engaging various slip and twinning activities.

The annealing phenomena are addressed in the third paper [3], where the authors investigated the phenomenon of dynamic recrystallization in a FeMnSiCrNi shape memory alloy subjected to hot compression. The analysis is based on a simulation of recrystallization using a cellular automaton approach, and apart from the microstructural features, both flow stress and dislocation density were predicted.

The latest advances in the modeling of mechanical performance based on microstructural features of polycrystalline systems are described in five contributions [4–8], which reveal the behavior of different materials under diverse processing operations. The springback phenomenon was simulated for the hot-stamping process in [4]. The calculations rely on the analysis of microstructures that evolved in the investigated boron steel. The relationship developed in the frame of the contribution [4] enables the estimation of the springback angle based on the area fraction of microstructures that evolved during the quenching. The following paper [5] deals with the analysis of the hot workability of 300 M steel, where the relation between the flow stress and microstructure evolution was examined. The investigation was conducted via in situ and ex situ compression tests. In another manuscript [6] equally concerned with hot deformation, the double-open multidirectional forging process is analyzed. In this investigation, the authors present a constitutive equation for GH4169 superalloy and simulate the microstructure evolution during dynamic recrystallization using the finite element approach. Anisotropy of plastic yielding is described in [7], where the authors investigate the mechanical response of the cross-rolled aluminum sheets. The simulation of cup earing is based on the evolved crystallographic texture. The model employed for the simulation of anisotropic behavior takes texture intensities into account, which appear in the {h00} pole figures. The performance of monocrystalline Cu-13.7% Al-4.2% Ni alloy is characterized in [8]. In that contribution, the authors investigate the so-called reversible martensitic transformation and describe both physical and structural changes in the complex metallic system induced by the thermocyclic treatments under applied loads.

Phase transformation is studied in [9]. In this study, high-strength steel was subjected to a hot-stamping process, and phase transformation was investigated via simulation of the transformation during the pressure-holding quenching process. The kinetics of phase alteration, which leads to the formation of ferrite, pearlite, and bainite, was analyzed according to the Kirkaldy–Venugopalan model.

#### **3. Conclusions and Outlook**

In this Special Issue, we present a wide spectrum of articles dealing with the modeling of microstructural aspects involved in deformation and recrystallization as well as simulation of microstructure-based and texture-based properties in various metals. The latest advances in the theoretical interpretation of mesoscopic transformations based on

experimental observations are also partially discussed. The studies dealing with the modeling of structure–property relationships are likewise analyzed in the present collection of manuscripts.

Although, in the submitted works, many hidden facets of various technological processes and related microstructural changes were revealed by employing advanced computational approaches, the contributions collected in this issue nevertheless clearly show that further efforts are required in the field of modeling to understand the complexity of the material's world. The ultimate goal of modeling efforts is, arguably, the development of a comprehensive model that is capable of describing many aspects of microstructure evolution during thermomechanical processing.

**Funding:** This research was funded by the EFOP-3.6.1-16-2016-00018 project: "Improving the role of research + development + innovation in higher education through institutional developments assisting intelligent specialization in Sopron and Szombathely". The work was performed in the frame of the Széchenyi 2020 program: "Innovative processing technologies, applications of energy engineering, and implementation of wide-ranging techniques for microstructure investigation".

**Conflicts of Interest:** The author declares no conflict of interest.