Coupling detailed radiation model with process simulation in Aspen Plus: A case study on ﬂ uidized bed combustor

(cid:129) Detailed radiation model coupled process simulation of ﬂ uidized beds. (cid:129) Combining the advantages of zone method and Aspen Plus. (cid:129) A case study on a 0.3 MW ﬂ uidized bed. (cid:129) A modest computing demand. While providing a fast and accurate tool for simulating ﬂ uidized beds, the major limitations of classical zero- dimensional ideal reactor models used in process simulations become irreconcilable, such as models built into commercial software (e.g. Aspen Plus ® ). For example, the limitations of incorporating heat absorption by the water wall and super-heaters and inferring thermal reciprocity between each reactor model/module. This paper proposes a novel modelling approach to address these limitations by incorporating an external model that marries the advantages of the zone method and Aspen Plus to the greatest extent. A steady state operation of a 0.3 MW atmospheric bubbling ﬂ uidized-bed combustor test rig was simulated using the developed modelling approach and the results were compared with experimental data. The comparison showed that the predictions were in agreement with the measurements. Further improvement is to be expected through incorporating more realistic zoned geometry and more complex reaction mechanisms. In addition, the developed model has a relatively modest computing demand and hence demonstrates its potential to be incorporated into process si- mulations of a whole power plant.


Introduction
Circulating Fluidized Bed (CFB) technology has been developed and served as an effective means for burning solid fuels for many industrial applications (such as coal combustion) since the 1980s due to its inherent advantages in fuel flexibility, high combustion intensity and low emissions [1]. Although technical knowledge about the design and operation of CFB is widely available for pilot plant and large scale units [2], few researchers have conducted modelling of the whole CFB power plant [3][4][5][6][7] and hence little is known about its dynamic performance and operational flexibility, like oxy-fuel combustion [8] and co-firing [9]. This was due to the challenges faced in modelling the highly complex gas/solid chemical reactions coupled with fluid flow and heat and mass transfer occurring in a CFB boiler.
One of the earliest simulation studies of fluidized beds is the National Aeronautics and Space Administration (NASA)'s 'Simulation of fluidized bed coal combustors' project conducted in the 1970s [10,11]. They developed a comprehensive model using in-house code to simulate the aforementioned physical and chemical processes in fluidized beds. Although the model is imperfect from the current point of view, it framed a modelling approach for fluidized beds and inspired the followup studies. de Souza-Santos [12][13][14][15][16][17][18] also developed a comprehensive model for the simulation of fluidized bed equipment, and the latest version includes all sub-models related to combustion and gasification of solid fuels and allows detailed simulation of boilers and gasifiers. With the maturity of commercial software for process simulations [19], various process simulators such as Aspen Plus® are available and have been widely employed for process simulation purposes by industrial entities since the late 1990s. The commercial process simulators usually have a powerful physical and chemical properties database, which is they're greatest advantage over models developed using in-house code.
Despite the advances in Aspen Plus for time-dependent dynamic simulation and process control, no provision is currently available to couple detailed heat transfer models with a process simulation as part of a complete engineering solution. As a result, in fluidized bed modelling the heat absorption by the water wall and super-heaters is usually specified as model inputs explicitly rather than predicted directly by the models. Furthermore, the modular modelling strategy does not involve thermal reciprocity between each reactor module [20][21][22][23][24][25][26][27][28]. These aspects of modelling are particularly useful for the design and evaluation of the operational flexibility of a fluidized bed boiler in supercritical conditions [29]. This is due to the fact that, under supercritical conditions, the boiler uses once-through circulation for the water wall and if excessive heat flux on the furnace side occurs, this cannot be automatically compensated for with a larger water flow rate in the tubes, thus resulting in over-heating of the water membrane.
On the other hand, although Computational Particle Fluid Dynamics (CPFD) is widely used to simulate fluidized bed combustion systems and can predict the detailed local heat and mass transfer information for the entire computational domain, it often takes several days, if not weeks, to provide useful results for large industrial cases. Therefore, CPFD models [30,31] are not feasible for incorporation into process simulations of a whole power plant. Moreover, Selçuk et al. [32,33], by using in-house codes investigated the radiative heat transfer characteristics on a bubbling fluidized bed combustor test rig. However, their developed fluidized bed combustion models relied on empirical models, rather than on first principles, to calculate gas temperature profiles and therefore may not generalise correctly in unknown situations.
As the combustion efficiency of fluidized bed combustor depends largely on the heat recovered in the freeboard region, where the dominant component of heat transfer is radiation [34], accurate modelling of the radiative heat transfer in such systems is essential. Previous work by Hu et al. [35] has been devoted to further improving the classical zone method of radiation analysis (hereafter called zone model) using in-house codes which can be applied to a wide range of industrial furnaces. In this paper, the authors further extend the application of the zone model by integrating it in an efficient manner with process simulation in Aspen Plus, thus providing a more complete and flexible solution for optimising the thermal performance of fluidized- bed boilers. As far as the authors are aware, no other process simulations of the fluidized bed boiler with similar features have been published. The developed model is expected to be relatively modest in computing demand and hence can be feasibly incorporated into dynamic simulations of a whole CFB power plant.

Object of the study
To demonstrate the proposed modelling approach, a 0.3 MW Atmospheric Bubbling Fluidized-bed Combustor (ABFBC) test rig, shown in Fig. 1, was chosen as the object for this study. The main feature of the test rig is the modular combustor formed of five 1 m high modules of internal cross-sectional area of 0.45 m × 0.45 m. The inner walls of the modules are lined and insulated with alumina-based refractory bricks of 6 cm in thickness. The first and fifth modules refer to the bed (bottom) and cooler (top), respectively. The ones between them are the freeboard modules. There exist two cooling surfaces in the modular combustor, as shown in Fig. 1, providing 0.35 m 2 and 4.3 m 2 of cooling surfaces, respectively.
Previously, Selçuk et al. [32,33] have conducted experimental and numerical studies on this test rig to investigate its heat transfer characteristics. The experiments, consisting of two combustion tests with and without recycling of fine particles from the cyclone, were carried out by burning a typical lignite coal. For the purpose of demonstrating the developed models, only the combustion test without particle recycling was simulated in this study. The coal analysis is as shown in Table 1 and the steady state operating conditions are presented in Table 2.

Methodology of fluidized bed modelling
Due to the complexity of coal combustion in fluidized beds, only the major steps of coal combustion are considered in the model, with some simplifying assumptions. When the coal particles travel through a fluidized bed combustor, drying, devolatilization, volatile combustion and char combustion occur consecutively, as illustrated in Fig. 2. A hydrodynamics model is also adopted to describe the influence of the gas-phase superficial velocity on the char particle pore size profile (which further influences the reaction rate of char combustion) along   the height of the freeboard region.

Devolatilization
With regard to the chemical formula for coal, it is treated for simplicity as a homogeneous mixture of elemental substances, such as carbon (C), hydrogen (H), oxygen (O), nitrogen (N) and Sulfur (S) as given in Reaction (R1). In this step, coal is converted into those virtual elementary substances in the simulation. As shown in Fig. 3, RYield 1 (Aspen Plus reactor) is used to model this coal decomposition by specifying the yield distribution vector (Eq. (1)) via an In-Line Fortran calculator 1 (Aspen Plus calculator block) according to the coal ultimate analysis.

Volatile combustion
During devolatilization, the Volatile Matter (VM) in coal is released and combusted instantaneously. The VM consists of volatile C (CISOLID 2 ), H 2 , O 2 , N 2 , S and H 2 O. The volatile combustion is assumed to be completed in the above bed region. Volatile C (CISOLID) is only oxidized to CO due the O 2 -lean environment in this region and is entirely consumed during the volatile combustion [20]. The amount of volatile C can be calculated by the difference between the Fixed Carbon (FC) in the proximate analysis and the total C in the ultimate analysis (Table 1) as Eq. (2), RStoic 3 (Aspen Plus reactor) is used to calculate this by specifying the fractional conversion of key components via an In-Line Fortran calculator 2 (Aspen Plus calculator block). For demonstrational purpose, three reactions are considered in this process, but any other reaction mechanisms can be easily incorporated:

Char combustion
The char particles resulting from the devolatilization process consist of the FC (NCSOLID 4 ) and ash. These particles are then burned to produce a mixture of CO and CO 2 in the entire fluidized bed boiler. Two reactions are considered in this process: CO .   2 CISOLID, conventional inert solid, is used for homogeneous solids that have a defined molecular weight [37]. 3 RStoic reactor, performs the calculations based on the reaction stoichiometry [37]. 4 NCSOLID, nonconventional inert solid, is used for heterogeneous solids that have no defined molecular weight [37]. 5 RCSTR reactor, performs simulation of ideal reactors operated under specific conditions (pressure and temperature or heat duty), valid phases, and user defined kinetics [37].  Table 3. Again, any other emission related reactions occurring in the entire fluidized bed boiler can be easily incorporated, such as SO 2 absorption [20], fuel and thermal NO x formations [38], via the user defined kinetics in the same way. Since those emission related reactions are not the main exothermal reactions and do not affect thermal performance appreciably, they are not considered in this study.

Hydrodynamics model
A lumped hydrodynamics model [40,41] was chosen to predict the voidage profile in the free board region of the studied fluidized bed combustor. According to the model, the free board region is divided into a discrete number of intervals and consists of two zones: a lower region and an upper region (Fig. 4).
In the lower region, because of the turbulent fluidization condition, the voidage (ε) can be assumed to be a constant [42] for a given gas superficial velocity (U g ); in the upper region, the value of ε in a certain interval between heights Z i−1 and Z i can be calculated based on Eqs. (3) and (4) via an In-Line Fortran calculator 3 (Aspen Plus calculator block) at every iteration: where ε * denotes the voidage under saturated conditions in the dilute phase at the top of fluidized bed riser column; ε b the voidage of the lower region, ε b ≈ 0.55 [40]; a represents the decay constant, and experimental data shows that U g is inversely proportional to the decay constant with a proportional constant of 4.0 s −1 [41,43]; ϕ represents the slip factor, which is a function of Froude number (F r ) as given in Eqs. (5)-(7), 6 3.

Zone method of radiation analysis
Due to the relatively simple geometry of the ABFBC, the freeboard region is split into a finite number of isothermal volume and surface zones along the height, as shown in Fig. 5. For modelling purposes, the bottom and the top surfaces, which also include the cooling tubes, are approximated each by an equivalent grey surface. An assumption of plug flow is applied to the flow pattern, which is reasonable given that the flow-induced enthalpy transport is predominantly in the longitudinal direction of the fluidized bed. However, even more complicated zoning arrangement and flow pattern are also possible, see [35], which depend on the requirements of a specific problem. An energy balance is then formulated for each zone, taking into account radiation interchange between all volume and surface zones, the enthalpy transport, and the source term associated with combustion and heat release [44]. Fig. 6 shows the flow sheet of the coupled solutions in Aspen Plus for ABFBC modelling. In the simulation, the RCSTR block temperatures are updated continuously by the zone model until an appropriate convergent tolerance is satisfied; i.e. negligible temperature difference between successive iterations.

Models with internal zone energy balance
The radiation term in the energy balance equations is written in terms of exchange factors known as the directed flux areas (DFAs) i j for gas-gas, gas-surface, surface-gas, and surface-surface exchange, respectively, in Eqs. (8) and (9)). The energy balances on all zones yield a set of simultaneous non-linear equations which can be solved to determine the temperature and heat flux on each zone.
For a system of N volume zones and M surface zones, the following energy balances can be written.
For the i-th volume (gas) zone, Qġ ,i represents the net rate of heat transfer to the volume (gas) zone; Qċ onv,i is the convection to all surfaces in contact with the volume zone; for the plug flow pattern assumed in this study where the enthalpy of the combustion products leaving one zone is the enthalpy input to an adjacent zone, Qė nth,i can be expressed as: Likewise, for the i-th surface zone, Qṡ ,i represents the net rate of heat transfer to the surface zone. If the surface temperature is not known, it  must be expressed as an additional equation in T i relating, for example, to a wall heat loss term or to the heat transfer to a load. The DFAs in Eqs. (8) and (9) represent the radiant flux between each zone-pair and can be expressed as an a-weighted summation of the total exchange areas (TEAs) (denoted by G G G S S G S S , , , i j i j i j i j ) for each of the radiating grey gases. This is a well proven method, known as the Weighted Sum of mixed Grey Gases (WSGG) model [44], to approximate the radiative behaviour of real gases, defined as follows: Gas -Surface directed flux area, Surface -Surface directed flux area, For a grey gas system, the total exchange areas are independent of temperature and need only be calculated once for a given geometry, and for fixed surface emissivity (ε) and attenuation coefficient (K). For a given boundary condition, these M + N non-linear equations (M unknown surface temperatures; N unknown gas temperatures) can be solved using the Newton-Raphson method [45]. This method provides an approach of computing successive approximations of the variables which converge towards the solution. Furthermore, an updated Monte-Carlo based Ray-Tracing (MCRT) algorithm [46] was used to calculate the radiation exchange areas. This, together with the retrieved enthalpy flow data in Fig. 6, provides the input data that needs to be supplied to the zone model.

Radiation properties of gas media (non-grey gas with particles)
In fluidized beds, particulate matter in the form of soot and ash particles exist in the gas medium and their sizes are close to, or larger than, the wavelengths of emission and absorption. Hence the impact of scattering due to particles can no longer be considered negligible compared with emission and absorption. When a beam of radiation is incident on an opaque particle, some of it will be absorbed and the remaining scattered away from the direction of the incident beam. The spectral attenuation of the beam by absorption and scattering can be expressed by Eq. (13): where, k g,λ , k s,λ , and k p,λ are the absorption coefficients of gas, soot, and particles at the specific wave number n, respectively. In this study, Truelove's WSGG model (N g = 3, N s = 2, ) [47] was used to calculate the gas radiative properties, and the grey gas parameters used in this model are listed in Table 4.
For the evaluation of absorption and scattering of radiation by the particles cloud, both the extinction efficiency factor (Q abs ) and the scattering efficiency factor (Q sca ) were evaluated via the classical Mie theory [48]. These are given by: The details related to these equations can be found in Ref. [49]. An in-house code was developed to evaluate these efficiency factors and the code had been validated with benchmark data (Table 4.1 in [49]). Then the absorption and scattering coefficients (k p,λ and σ sp,λ ) of the suspended particles can be obtained by: Similar to the non-grey gas properties, calculation of the absorption and scattering coefficient of particles requires the evaluation of these properties as a function of wavelength, which means a high computational cost. For engineering applications, grey approximation for the suspended particles is a reasonable compromise. A mean wavelength (λ m ) was used to represent the spectral distribution of the incident radiation energy. This wavelength can be calculated from the displacement law for black-body radiation [50], given by Eq. (19).
here T m is the mean bulk temperature of the gas media in the enclosure, and assumes that suspended particles have the same temperature as the gas phase since previous Computational Fluid Dynamics (CFD) simulation showed that their temperature difference is always less than 50 K [51]. Radiation properties of the gas media in this study are summarized in Table 5.

Thermal boundary conditions at walls in the zone model
In most instances, heat losses from furnace walls and/or heat transfer to a water-cooled loads exist, which are generally not known prior to measurements. In this study, the freeboard region is insulated with refractory bricks 6 cm in thickness and cooling tubes are arranged at the bed and exit areas and the mean bed temperature was measured as 1148 K in the experiment in Ref. [33]. Based on these facts, an adiabatic boundary condition was assumed for the side walls. The two imaginary surfaces were assumed as having the same temperature as their respective adjacent upstream gas zones where cooling tubes are arranged. Due to the plug flow assumption, it can be seen that only the cooling system in the bed region has a significant effect on convection terms in the zone model, as indicated by the red-dash 7 line in Fig. 7; the gas flowing through the cooling system arranged at the top region has limited effect on its upstream processes. Therefore, the overall convective heat transfer coefficient (U in Eq. (20)) of the cooling system at the bed region was tuned to match the bottom gas zone temperature (T g,1 ) to the measured bed temperature. The convective heat transfer is given by,

Results and discussion
To demonstrate the proposed modelling approach, firstly the voidage profile was evaluated and compared with available experimental data; then the direct exchange areas (DEAs) required by the zone model was prepared according to the geometry of the studied fluidized bed combustor using the MCRT algorithm. The computational efficiency of the MCRT was also compared to a direct numerical integration algorithm. Next, based on the developed model, the thermal performance of the fluidized bed combustor and reaction kinetics characteristics was predicted and compared with experimental data. Finally, the possibility of applying this model to large-scale fluidized bed boilers is discussed. It should be noted that the accuracy of the zone model is only significantly affected by the number of divisions along the height if a large temperature gradient exists, which was not the case in this study as indicated by the experimental measurements. Therefore 10 divisions were chosen which gave adequate resolution and yet meaningful comparison with actual measurement points.

Voidage profile
In fluidized bed design, the decay constant is an important parameter, but their reported values are scattered and sketchy. Based on available experimental data, Kunii et al. [41] drew an important conclusion that the mechanism of decay for solid fraction in the upper region of fast fluid beds is basically similar to the decay above bubbling and turbulent beds. This conclusion makes Eq. (3) applicable to all three types of beds. In a specific case, only the voidage of the lower region or solid fraction needs to be determined carefully. The authors also reviewed the literature and concluded the values for voidage at the lower region (ε b ) 8 as follows, bubbling bed: ε b = 0.45-0.60; turbulent bed: ε b = 0.60-0.78; fast fluidization: ε b = 0.78-0.84. The calculated voidage profile of the bubbling bed in the process simulation was Table 4 Grey gas parameters used in the WSGG model [47]. plotted in normalized height and compared with that of an ultra-tall CFB riser [53] in Fig. 8. From the comparison, it can be observed that appreciable differences exist at the lower region. These may be due to the inherent differences in features of various beds. More precise data that reflect how the voidage may vary with the imposed operating conditions would need to be further investigated.

Direction exchange areas (DEAs)
Although the radiation exchange areas used in the zone model include DEAs, TEAs, and DFAs, the main computing time is due to the calculation of the DEAs, which are functions of the geometry and gas attenuation coefficient. Since all the radiation emitted by a surface must go somewhere, the sum of all the DEAs for a surface is its area and from the symmetry of the integral, it follows that: The accuracy of the MCRT algorithm highly depends on the ray density. A large ray density (> 10,000 ray m −2 ) is required to ensure low errors [44]. Three ray densities (10,000 ray m −2 , 20,000 ray m −2 , and 50,000 ray m −2 ) were applied in the DEAs calculations. An algorithm for smoothing the approximate exchange areas [54] was used in the updated MCRT algorithm to ensure that both the summation and symmetry rules of Eqs. (21) and (22) are satisfied within a convergent criteria of 1e−10. Even with such a strict convergent criteria the     Table 6. On the right hand side, the accuracy of the direct numerical integration algorithm depends on the number of volume elements and increases exponentially with it. Without using the smoothing algorithm, relatively larger errors in energy balances may be expected, especially in the case with a small number of volume zones. Despite significant improvement in accuracy (by almost eight orders of magnitude) using the updated MCRT algorithm, the required CPU time is still much less than the direct numerical integration algorithm.

Thermal performance
The simulation of a steady state operation of the studied fluidized bed showed that the calculation took 0.40 s to reach a converged solution after 3 iterations, and the maximum temperature difference between the values calculated by RCSTR blocks and zone model was less than 10 −4 K. The measured and predicted thermal performance, including gas temperature profile and mean incident radiative heat flux, is shown in Figs. 9 and 10, respectively. The gas temperature profiles are found to be in good agreement, but appreciable discrepancies are found in the incident radiative heat flux at the locations close to the bed and the cooler. These discrepancies might be due to the equivalent surfaces which do not really exist and this would interdict 'cold' beams from the cooling tubes, and the radiometer probes located at those heights are also affected by the cooling tubes. Further sensitivity study on the emissivity of the wall surfaces showed that the emissivity has almost has no effect on the gas temperature profiles but a significant affect the incident radiative heat flux. In general, the incident radiative heat flux increases with increasing surface emissivity. This is because the temperature profile of the system will be constant when thermal equilibrium is attained under steady-state operation. The net incident radiative heat flux on a surface depends on the fourth-power temperature difference between gas and the surface and their emissivities. Hence, increasing the water wall emissivity will help to improve the boiler thermal performance. For the studied fluidized bed, since the inner walls of the freeboard are refractory lined and insulated, heat losses are negligible. The wall temperatures of the freeboard are therefore close to the adjacent gas temperatures during the steady state operation. Fig. 11 further presents the contour plot of average temperature and incident heat flux on each walls' surface zones along the freeboard region which are not available from traditional process simulations using Aspen Plus. Due to the plug flow and adiabatic wall assumptions, the temperature profile among the surrounding wall surfaces of each gas zone are uniform; and the single surface zone arrangement on the XZ plane also makes the radial temperature profile not available. However, these details can be made easily available, if necessary, by extending the zone model to multi-dimensions [35].
Furthermore, Fig. 12 shows the heat duty of each reactor due to enthalpy transport in the process simulation. Since the first reactor (module) is adjacent to the cooling tubes located in the bed region, most of the heat released (about 120 kW) is taken away by the cooling water. The heat duty of the reactors dramatically decreased with the development of the reaction process. The last reactor (module) has the minimum heat duty as the combustion is close to completion there.
In order to inspect the system energy audit, the energy input (370 kW = 13.2 MJ kg −1 × 101 kg s −1 ) in the simulation was broken  down into the heat duty of bed cooling, heat duty of top cooling plus loss in exhaust 9 , and loss in unburned coal, as plotted in Fig. 13. The combustion efficiency can be evaluated based on the unburned coal, which was about 95% in this study. This result is within a reasonable range of the normal combustion efficiency (over 90%) of bubbling fluidized beds in the technical report [55].

Reaction kinetics
Not only the thermal performance, but also the reaction kinetics characteristics of each gas/volume zone are available in the process simulation. Fig. 14 shows the reaction rate of key components in the exothermic reactions. A positive value means production and a negative value means consumption. From this figure it can be seen that O 2 is always being consumed and CO 2 is always being produced, while the reaction rates drop asymptotically to zero, as expected. The reaction rate of CO highly depends on the combustion atmosphere. Due to the assumption that Volatile C (CISOLID) is only oxidized to CO in the O 2lean environment at bed level, as mentioned in Section 3.2, a relatively large amount of CO is therefore present at the lower region of the freeboard region. In practical operation, appropriate air staging can be applied to ensure low pollutant emissions and completed combustion [56]. Fig. 15 shows the comparison of the measured and predicted O 2 , CO 2 and CO concentrations in the process simulation. Despite significant deviations, the predicted trends of concentration are in     agreement with the measurements. These results prove the effectiveness of the proposed modelling approach to some extent. The deviation in the absolute values is most likely due to the over simplification of the combustion mechanisms applied in this study, and the fact that combustion reactions are influenced by many factors such as the type of fuel, and the physical state of reactants and their dispersion, which are not fully considered. In addition, the assumption of virtual elemental substances during devolatilization may also cause the deviations since this is not the case in the real process. In order to improve the prediction, more complicated combustion mechanisms [2] have to be involved to accurately calculate the O 2 consumption rate in its related reactions. Further, a devolatilization mechanism which is more similar to the real process is also helpful to improve the prediction.

Discussion on the real application of the proposed modelling approach
The results shown above have demonstrated the application of the proposed modelling approach on a bubbling fluidized bed combustor. For real-time simulations of a whole power plant, the prerequisite is that the model execution is faster than real-time. The modest computing demand implies that the developed model can be potentially incorporated into process simulations of a whole power plant. When applying the proposed modelling approach to a large-scale fluidized bed boiler, the main difference would lie in boundary conditions of heat transfer through water walls. In subcritical operations, heat is transferred to a boiling fluid (vapour-liquid phase) at an operating temperature, as indicated by the solid red line in Fig. 16. Therefore, the Dirichlet (or first-type) boundary conditions [57] are applicable to the water walls. Even though excessive heat flux on the furnace side may occur, it can be automatically compensated for with a larger water flow rate in the tubes. Latent heat of phase changes takes the role of a buffer and maintains the temperature of the water membrane within a safe limit. In supercritical or ultra-supercritical cases, the conditions are much more sophisticated if excessive heat flux on the furnace side occurs. This is because, under supercritical conditions, only one state of phase exists. The temperature of the water membrane is very sensitive to the change of heat flux on the furnace side. If excess heat flux on the furnace side occurs and cannot be compensated in time, it can easily lead to heat transfer deterioration of the whole system. That is the reason why involving heat transfer to water walls is of particular importance to the evaluation of operating conditions and flexibility of a fluidized bed boiler in supercritical conditions. In that case, the water wall should be considered as heating load (or internal obstacle) individually in the zone model and the heat conduction through the water-cooled wall to the supercritical fluid needs to be considered as well. By doing so, the local heat absorption by water walls is incorporated and available in the process simulation.
In addition, the implication of the current work highlighted the existing limitation of the sequential modular strategy in Aspen Plus, in that there is no thermal reciprocity within its modelling framework. By incorporating the zone model with Aspen Plus interface in this work, the energy balance takes account of radiation interchange between all volume and surface zones, the enthalpy transport, and the source term associated with combustion and heat release. Although the developed model is one-dimensional in this case, the proposed modelling approach is equally applicable to two-and three-dimensional cases. For a multidimensional model with more complex geometry, the plug flow assumption is no longer applicable. The enthalpy transport term needs to be expanded to include the flows in all relevant directions. The flow data may be derived from other means such as physical models and isothermal CFD simulations [36].

Conclusions
This paper successfully demonstrated a novel modelling approach for fluidized beds which incorporated the classical zone method with Aspen Plus interface. A main advantage of the proposed modelling approach stems from the thermal reciprocity within a combustor derived from the zone model which resolves the limitation of the sequential modular strategy in Aspen Plus. A steady state operation on a 0.3 MW atmospheric bubbling fluidized-bed combustor test rig was simulated and completed in a few seconds of CPU time. The predicted thermal performance and the reaction kinetics of the studied fluidized bed combustor (such as the temperature profile, incident radiative heat flux) are in agreement with the experimental data; and the predicted overall combustion efficiency is also within a reasonable range of the normal combustion efficiency of bubbling fluidized beds. Relatively modest computing demand and acceptable accuracy make it possible for the developed model to be incorporated into process simulations of a whole power plant and to be used to study the operational flexibility of the whole power plant.