Analysis of Fluid Flow and Heat Transfer inside a Batch Reactor for Hydrothermal Carbonization Process of a Biomass

: This work analyzes the heat transfer and ﬂuid ﬂow within a batch reactor for hydrothermal carbonization (HTC) of raw olive pomace (ROP). The autoclave is partially ﬁlled with a mixture of ROP and distilled water and hence it is considered as a dispersed medium. The reactor is heated through its lateral surface, whereas the bottom wall and the upper surface of the mixture are thermally insulated. Under the effect of heat and pressure, the ﬂuid moves inside the reactor, while particles are subject to other forces. Additionally, the biomass (ROP) is decomposed into very ﬁne particles to produce a solid product (hydrochar). COMSOL Multiphysics software is used for the analysis of heat transfer and ﬂuid dynamics. Chemical kinetics of the reactions are modeled by a basic kinetics model. Numerical results are validated using experimental data carried out in similar operating conditions. They are in good agreement since the deviation between them does not exceed 6%. Isotherms, velocity ﬁelds, and isobars are evaluated within the reactor as well as velocity and distribution of particles. These amounts are inﬂuenced by the imposed heat ﬂux at the lateral wall (q 0 ). Also, it has been shown that the temperature and pressure values reached are above those required by the HTC process and, consequently, a HTC reactor could be designed with optimal operating conditions.


Introduction
The conversion of lignocellulosic biomass into biofuel is quite complex because it occurs following several parallel pathways whose mechanisms are still unclear. This transformation is usually performed by two main processes: the thermochemical and biochemical processes [1][2][3][4]. In this work, particular attention is paid to one of the most promising thermochemical processes, namely hydrothermal carbonization (HTC) which is a physico-chemical process for converting organic compounds into a valuable carbon material in the form of a hydrochar (a solid product) and other products (gases and aqueous phase). The HTC process can be defined as the combined dehydration and decarboxylation of a biomass (fuel) to increase its carbon content and, therefore, its higher heating value (HHV). It involves the use of wet biomass under temperature values ranging from 180 to 280 • C and high autogenous pressure (up to 8 MPa) in a closed reactor for several hours [5].
Available literature, concerning the above process, is mainly focused on physicochemical characterization and thermal analysis of obtained hydrochar [6][7][8][9][10], reaction kinetics modeling based on Arrhenius's law [11][12][13][14], and slightly less for heat transfer and fluid-dynamics modeling [15][16][17][18][19]. Indeed, Benevente et al. [6] used HTC technology to obtain useful bioenergy feedstock from three by-products, including olive mill, canned artichoke and orange wastes. Their work was conducted in a laboratory-scale reactor under high pressure and temperatures ranging from 200 to 250 • C for durations of 2, 4, 8 and 24 h. Effects of the residence time and the operating temperature on the hydrochar properties were analyzed using various characterization techniques. For the three carbonized organic wastes, the carbon content and the higher heating value (HHV) are close to those of black carbon. They also showed that the HTC process would allow up to 50% of energy savings compared with the torrefaction process. In 2019, Lucian et al. [11] utilized olive trimmings to achieve an experimental study and an original HTC kinetics model. The biomass was carbonized at a temperature ranging from 180 and 250 • C for a residence time of up to 8 h. Their experimental data and the findings of two other agro-wastes (grape marc and spineless cactus) previously were compared with the obtained model. They concluded that the model predicts the carbon distribution correctly inside HTC products with respect to time, including the thermal unsteady heat period. In the same year (2019), the HTC process was applied by Micali et al. [14] to olive pomace at temperatures ranging from 180 to 305 • C and a residence time between 60 and 180 min. The kinetic of hydrothermal carbonization was modeled by estimating the activation energy of hydrochar and the pre-exponential factor of Arrhenius's law using the least-squares method. The developed model predicts the temporal evolution of the product's (hydrochar) mass yield at different temperatures of the process.
The modeling of a HTC reactor has to cover the main thermo-fluid aspects such as heat transfer phenomena and fluid-dynamics connected with operating conditions. To the best of our knowledge, studies included heat transfer and computational models [15][16][17], coupled with kinetics modeling of the HTC process of lignocellulosic biomass, are still lacking. Nevertheless, Baratieri et al. [15] simulated the transient thermal behavior of the HTC batch reactor containing grape marc at several operating conditions (imposed temperatures and residence times). They established a simplified thermal dynamic model, based on lumped capacitance method, which was combined with a kinetic model to evaluate the HTC process performance. These types of model are not overly abundant because of the lack in data regarding the thermodynamics of HTC and the change of physical properties of the mixture [20]. As for the study of Alvarez-Murillon et al. [16], it modeled kinetic of cellulose hydrothermal carbonization taking into account heat transfer within an autoclave. The heat balance equation including conduction, convection and radiation terms was solved by COMSOL Multiphysics software under appropriate boundary and initial conditions. The authors showed the influence of the reactor (autoclave) temperature on the biomass decomposition represented by the pre-exponential and exponential factors evolving the reaction rate.
In 2020, thermo-fluid dynamic performance of the HTC process of olive pomace in a batch reactor has been developed by means of ANSYS Fluent software [17]. In this work, the authors considered the multiphase model composed of three phases (air, water and dried olive pomace). The mixture model, which is applicable for fluid-solid flow with dispersed phases, has been chosen to simulate the flow within the reactor. The conservation equations have been solved to simulate dynamic, thermal and kinetics behavior inside a batch reactor. This experimentally validated work is undoubtedly an important source of information on temperature, velocity and the spatial distribution of the mixture inside the reactor during the HTC process.
More recently (2021), Ischia and Fiori [18] published a review paper which examine and analyze some predicting tools such as kinetics, statistical, heat transfer and computational models, highlighting their potentialities and limits. A mathematical model for HTC has been developed by Heidari et al. [19], in 2021, and comparison with experimental results was performed in the case of pine wood particles. In this study, heat transfer rate, reaction kinetics, and the porous structure of the biomass have been estimated by using COMSOL Multiphysics software. Then, the predicted results were validated and operating parameters such as biomass to water ratio, temperature, residence time and power consumption were recorded.
From the above brief state of the art, it is obvious that studies on modeling and simulation of the hydrothermal carbonization (HTC) process are lacking in the literature. Accordingly, this work is a modest contribution that could help remedy this shortfall. It deals with a mathematical modeling and a computational fluid dynamics (CFD) simulation of the HTC conversion of biomass. Thus, we consider an autoclave (batch reactor) partially filled with a mixture of raw olive pomace (ROP) and distilled water in well-defined proportions. The heat flux (q 0 ) is imposed at its lateral surface (wall). Among the objectives of this study is to ensure that the conditions of temperature and pressure required by the HTC process will be reached (180 • C ≤ T ≤ 250 • C and 10 bar ≤ P ≤ 30 bar). The influence of q 0 on fluid flow and heat transfer inside the autoclave will be illustrated too. The assessment of autogenous pressure and the tracking of solid particles within the reactor are also original aspects of this work. Besides, the developed model will be used in designing an industrial reactor for the HTC process of biomass.

Materials
A photo of a treated ROP sample is provided in Figure 1. This agro-food residue is a product obtained from the olive oil production process. This biomass, coming from a crushing unit located in Meknès region (Morocco), has an average particle size of 5 mm, and an initial humidity rate ranging between 30% and 40%. Also, it has been dried using an indirect forced convection solar dryer and then characterized [4,[21][22][23][24]. Table 1 summarizes ultimate and proximate analyses of the selected biomass.
transfer rate, reaction kinetics, and the porous structure of the biomass have be mated by using COMSOL Multiphysics software. Then, the predicted results we dated and operating parameters such as biomass to water ratio, temperature, re time and power consumption were recorded.
From the above brief state of the art, it is obvious that studies on model simulation of the hydrothermal carbonization (HTC) process are lacking in the lit Accordingly, this work is a modest contribution that could help remedy this sho deals with a mathematical modeling and a computational fluid dynamics (CFD) tion of the HTC conversion of biomass. Thus, we consider an autoclave (batch partially filled with a mixture of raw olive pomace (ROP) and distilled w well-defined proportions. The heat flux (q ) is imposed at its lateral surface Among the objectives of this study is to ensure that the conditions of temperat pressure required by the HTC process will be reached (180 °C ≤ T ≤ 250 °C and 10 ≤ 30 bar). The influence of q on fluid flow and heat transfer inside the autoclave illustrated too. The assessment of autogenous pressure and the tracking of solid p within the reactor are also original aspects of this work. Besides, the developed will be used in designing an industrial reactor for the HTC process of biomass.

Materials
A photo of a treated ROP sample is provided in Figure 1. This agro-food resi product obtained from the olive oil production process. This biomass, coming crushing unit located in Meknès region (Morocco), has an average particle size o and an initial humidity rate ranging between 30% and 40%. Also, it has been drie an indirect forced convection solar dryer and then characterized [4,[21][22][23][24]. Table  marizes ultimate and proximate analyses of the selected biomass.

Physical Model and Hypotheses
The hydrothermal carbonization is a physicochemical process that occurs at high autogenous pressure (10-30 bar) and temperature ranging from 180 to 250 • C. It allows the conversion of organic compounds (biomass), placed in a batch reactor (autoclave), into a solid product called hydrochar (HC). Besides, the liquid mixture and gases obtained from this operation will not be treated in this study.
This modeling study consists of a circular cylinder (autoclave) with a height of 10 cm and a radius, R, of 3 cm (Figure 2). It is partially filled with a mixture of ROP and distilled water whose height is H and a biomass-to-water ratio, η, of 1/3. The reactor is heated through its lateral surface whereas the bottom wall and the upper surface of the mixture are thermally insulated. Under the effect of heat and pressure, the fluid moves inside the reactor, while particles are subject to other forces (drag force, Basset force etc.). Therefore, the mixture forms a dispersed medium which is composed of distilled water and suspended particles. Upon heating, the ROP is decomposed into very fine particles to produce a solid product (hydrochar) which is richer in fixed carbon than the ROP material [8].

Physical Model and Hypotheses
The hydrothermal carbonization is a physicochemical process that occurs at high autogenous pressure (10-30 bar) and temperature ranging from 180 to 250 °C. It allows the conversion of organic compounds (biomass), placed in a batch reactor (autoclave), into a solid product called hydrochar (HC). Besides, the liquid mixture and gases obtained from this operation will not be treated in this study.
This modeling study consists of a circular cylinder (autoclave) with a height of 10 cm and a radius, R, of 3 cm (Figure 2). It is partially filled with a mixture of ROP and distilled water whose height is H and a biomass-to-water ratio, η, of 1/3. The reactor is heated through its lateral surface whereas the bottom wall and the upper surface of the mixture are thermally insulated. Under the effect of heat and pressure, the fluid moves inside the reactor, while particles are subject to other forces (drag force, Basset force etc.). Therefore, the mixture forms a dispersed medium which is composed of distilled water and suspended particles. Upon heating, the ROP is decomposed into very fine particles to produce a solid product (hydrochar) which is richer in fixed carbon than the ROP material [8]. The mathematical modeling of the HTC process is undertaken by considering the following simplifying assumptions: • the physical model is 2D and axisymmetric; • thermophysical properties of the fluid (distilled water) depend on temperature; • the fluid is assumed to be Newtonian and the flow is laminar (Rep < 1); • solid particles are supposed to be homogeneous and spherical; • the solid particles and the fluid are in local thermal equilibrium; • the mixture is considered as a dispersed phase medium and it is dealt with according to the Eulerian-Lagrangian approach; • there is no interaction between particles (electrostatic or otherwise) nor collision; • viscous dissipation and radiative heat transfer are negligible; • there is no gaseous phase on the top of the free-surface of the mixture.

Continuous Phase (Fluid)
This phase is governed by fundamental equations describing the conservation of mass, momentum and energy, respectively:  The mathematical modeling of the HTC process is undertaken by considering the following simplifying assumptions:

•
The physical model is 2D and axisymmetric; • Thermophysical properties of the fluid (distilled water) depend on temperature; • The fluid is assumed to be Newtonian and the flow is laminar (Rep < 1); • Solid particles are supposed to be homogeneous and spherical; • The solid particles and the fluid are in local thermal equilibrium; • The mixture is considered as a dispersed phase medium and it is dealt with according to the Eulerian-Lagrangian approach; • There is no interaction between particles (electrostatic or otherwise) nor collision; • Viscous dissipation and radiative heat transfer are negligible; • There is no gaseous phase on the top of the free-surface of the mixture.

Continuous Phase (Fluid)
This phase is governed by fundamental equations describing the conservation of mass, momentum and energy, respectively: Energies 2022, 15 I is the indentity tensor.
Note that physical quantities used in the above equations are reported in the Nomenclature section.

Dispersed Phase
The dispersed phase model (Eulerian-Lagrangian approach) takes into account the motion of each particle separately and the equation of motion is written for each particle [25]: → F ma and → F PG represent body force, drag force, Basset force, virtual mass force and pressure gradient force, respectively. Most of them depend on the fluid velocity or pressure and they are all defined in the Appendix A.

Initial Conditions
The initial conditions associated with Equations (1)-(4) are: The process is supposed to start under vacuum condition to avoid possible chemical reactions (combustion) which could occur between gases (see Equation (7) below) and air. Therefore, P(t = 0, r, z) = 10 −6 bar (5d)

Boundary Conditions
The symmetry at the axis r = 0 leads to: At the vertical wall which corresponds to r = R: At the bottom of the cylinder (z = 0) and at the upper surface of mixture (z = H): The conservation Equations (1)-(4) have been resolved using COMSOL Multiphysics software which is based on the finite element method. Indeed, the time derivative of the continuous phase (fluid) was discretized by an implicit scheme of second order using the backward differentiation formula. Then, the problem was solved using the PARDISO solver combined with the α-generalized solver by considering relative tolerances of 10 −3 . The heat equation was discretized by linear elements whereas fluid velocity and pressure were determined using quadratic (P 2 ) and linear (P 1 ) elements, respectively [26]. Finally, the trajectory of the particles was simulated based on the calculation of velocity and pressure in the continuous phase (water).

Mesh Sensitivity
To perform mesh sensitivity, four densities (4980; 8260; 10,004 and 18,000) were tested in triangular mesh elements and the results are presented in Table 2 below: By comparing the computation times and the relative errors (error between the simulated and the experimental temperature results), it appears that the optimal density in our case is 10,004 triangular elements ( Figure 3), because the error does not vary much when increasing the density, while the computation time becomes very important.

CFD Modeling
The conservation Equations (1)-(4) have been resolved using COMSOL Multiphy ics software which is based on the finite element method. Indeed, the time derivative the continuous phase (fluid) was discretized by an implicit scheme of second order usi the backward differentiation formula. Then, the problem was solved using the PARDIS solver combined with the α-generalized solver by considering relative tolerances of 10 The heat equation was discretized by linear elements whereas fluid velocity and pressu were determined using quadratic (P2) and linear (P1) elements, respectively [26]. Final the trajectory of the particles was simulated based on the calculation of velocity an pressure in the continuous phase (water).

Mesh Sensitivity
To perform mesh sensitivity, four densities (4980; 8260; 10,004 and 18,000) we tested in triangular mesh elements and the results are presented in Table 2 below: By comparing the computation times and the relative errors (error between t simulated and the experimental temperature results), it appears that the optimal dens in our case is 10,004 triangular elements (Figure 3), because the error does not vary mu when increasing the density, while the computation time becomes very important.

Kinetics Analysis
As previously defined, hydrothermal carbonization (HTC) is a pre-treatment process to convert various raw materials (e.g., biomass) into solid fuel with high energy density. Understanding the reaction kinetics taking place during the biomass conversion is needed to design and optimize a HTC reactor. Certainly, thermo-fluid dynamics, chemical and mechanical phenomena occur throughout the degradation of a composite placed in the reactor, and they manifestly interact with each other in a significant way. A kinetics model must, therefore, consider all these phenomena in order to be realistic. We assume that the decomposition of biomass occurs inside the batch reactor through a first-order reaction rate according to the following reaction mechanisms: The reaction rate for ROP and hydrochar (HC) are respectively defined as: F HC is the fraction of ROP that is converted to HC (see Table 3) at the end of the biomass (ROP) decomposition. The reaction rate constant k (see Equation (8c)) is given by the Arrhenius law: E a (13.91426 kJ/mol) represents the activation energy; A (0.45 min −1 ) is the frequency factor (pre-exponential factor) and R g is the gas constant (8.314 J mol −1 K −1 ). They were deduced, by the least square method, from experimental results of Missaoui et al. [8] as well as the fraction F HC (see Table 3).
Finally, the reaction rate k (Equation (8c)) was obtained using the average temperature values obtained from the thermo-fluid dynamics model detailed above.

Results and Discussion
As mentioned above the thermophysical properties of water are subject to the pressure and temperature. They were obtained from reference [27] while those of ROP were previously published [23].
Simulations were performed for a batch reactor having a diameter d of 6 cm and a height of 10 cm. It was filled with a mixture (distilled water + ROP) with a biomass-to-water ratio of 1/3 (η = Dry biomass weight/Water weight). The processing time (transient heat transfer duration + residence time) considered in this study is equal to 3000 s.

Model Validation
To predict temperature distribution inside the batch reactor, Equations (1)-(3) were simultaneously solved taking into account their initial and boundary conditions. The obtained temperatures allow us to estimate temporal variations of the reaction kinetics of each solid species (biomass, hydrochar) inside the reactor.
The proposed model must be validated to check its accuracy and appropriateness. Validation tests were conducted in a volume (V = 300 mL) of the batch reactor (Top Industrie, France) made of Inconel 718 allowing to reach a maximal temperature of 300 • C and a . In order to carbonize the ROP, the biomass-to-water ratio was taken equal to 1/3 and m ROP = 60 g. The mixture was locked inside the autoclave under vacuum conditions and stirred with a speed rate of 300 rpm. Then, it was heated by an electric oven which delivers a maximum power of 2300 W. For more details concerning the experimental procedure and measuring instruments, readers should refer to Missaoui et al. [8].
Comparison between CFD modeling and experimental results are given in two Figure 4a,b. The first one (Figure 4a) depicts temperature variations inside the reactor and it shows that the steady-state heat transfer, considered in the HTC process, is reached at about 40 min. The difference between measurement data and simulation temperatures does not exceed 6% and hence they are in good agreement. Nevertheless, the difference is significant (about 20%) in the case of autogenous pressure, especially during the unsteady-state regime. This difference could mainly be explained by the fact that vacuum conditions were imposed in the upper part of the batch reactor during the simulation period while various gas from HTC process (wet pyrolysis) penetrate this part throughout the experimental test. Also, this discrepancy may have been caused by the level of non-homogeneity in the mixture. Despite this, the proposed model is acceptable because the HTC process is performed in the steady-state regime for which the measured and simulated pressures are in close agreement (the difference is of 6% in the worst case). dustrie, France) made of Inconel 718 allowing to reach a maximal temperature of 300 °C and a maximal pressure of 200 bar. Experiments were carried out at the ICARE laboratory (see Abbreviation). In order to carbonize the ROP, the biomass-to-water ratio was taken equal to 1/3 and m ROP = 60g . The mixture was locked inside the autoclave under vacuum conditions and stirred with a speed rate of 300 rpm. Then, it was heated by an electric oven which delivers a maximum power of 2300 W. For more details concerning the experimental procedure and measuring instruments, readers should refer to Missaoui et al. [8].
Comparison between CFD modeling and experimental results are given in two Figure 4a,b. The first one (Figure 4a) depicts temperature variations inside the reactor and it shows that the steady-state heat transfer, considered in the HTC process, is reached at about 40 min. The difference between measurement data and simulation temperatures does not exceed 6% and hence they are in good agreement. Nevertheless, the difference is significant (about 20%) in the case of autogenous pressure, especially during the unsteady-state regime. This difference could mainly be explained by the fact that vacuum conditions were imposed in the upper part of the batch reactor during the simulation period while various gas from HTC process (wet pyrolysis) penetrate this part throughout the experimental test. Also, this discrepancy may have been caused by the level of non-homogeneity in the mixture. Despite this, the proposed model is acceptable because the HTC process is performed in the steady-state regime for which the measured and simulated pressures are in close agreement (the difference is of 6% in the worst case).

Fluid Flow and Heat Transfer Analysis
The required conditions for the HTC process have already been defined in the Section 2.2 above: (180 °C ≤ T ≤ 250 °C and 10 bar ≤ P ≤ 30 bar). Hence, one of the main objectives of this study is to ensure that these conditions are properly met during the process. Consequently, we need to analyze more precisely the obtained results, especially isotherms, velocity fields, and isobars for the continuous phase (water), as well as the distribution of particles for the solid phase. Without loss of generality, results are depicted, near the wall, on a vertical plane perpendicular to both bases of the reactor (cylinder) and passing through r = 2.7. The simulations were recorded at t = 600, 1800, and 3000 s, immediately after the onset of the HTC process, for an imposed thermal power of 2300 W (q = 2300 W), and biomass-to-water ratio (η = 1/3).
Isotherms, temperature profiles, velocity field and isobars within the batch reactor, are respectively exhibited by Figures 5a-c, 6, 7a-c and 8a-c. It is obviously clear from

Fluid Flow and Heat Transfer Analysis
The required conditions for the HTC process have already been defined in the Section 2.2 above: (180 • C ≤ T ≤ 250 • C and 10 bar ≤ P ≤ 30 bar). Hence, one of the main objectives of this study is to ensure that these conditions are properly met during the process. Consequently, we need to analyze more precisely the obtained results, especially isotherms, velocity fields, and isobars for the continuous phase (water), as well as the distribution of particles for the solid phase. Without loss of generality, results are depicted, near the wall, on a vertical plane perpendicular to both bases of the reactor (cylinder) and passing through r = 2.7. The simulations were recorded at t = 600, 1800, and 3000 s, immediately after the onset of the HTC process, for an imposed thermal power of 2300 W (q 0 = 2300 W), and biomass-to-water ratio (η = 1/3).
Isotherms, temperature profiles, velocity field and isobars within the batch reactor, are respectively exhibited by Figures 5a-c, 6, 7a-c and 8a-c. It is obviously clear from Figures 5a-c and 7a-c that profiles are substantially symmetrical relative to the plane of symmetry (r = 0).         Moreover, the reactor bottom is slightly hotter than its top, regardless of time. Moving away from the inner lateral surface, Figure 6 reveals that there is an average temperature difference of 5 °C between both base surfaces (bottom and top surfaces). This fact can be explained by suspension particles [28,29] and their migration to the reactor top (Figure 9a-c) owing to their low bulk density [2] as well as their thermal conductivity [23] compared to those of water. Figure 5a-c and also Figure 6 prove that there are wide variations of temperature gradients in the vicinity of the vertical wall (r = R), while they become small far enough from it ( Figure 6). In the same region, the fluid remains practically stagnant because of the imposed no-slip condition applied at the solid-fluid boundary (Figure 7a-c). Far beyond this zone (vertical wall) and also from the lower and upper surfaces of the mixture, natural convection due to the bulk movement of water molecules takes place, and hence it generates convection cells (Figure 7a-c). In the same way, another set of large cells is developed along the symmetry axis. Figure 8a-c display isobars at three different times (t =600, 1800 and 3000 s). We observe that, at each time, autogenous pressure value remains constant with an increase from 1.8 bar to 28.09 bar at 600 and 3000 s, respectively. The last pressure value and its corresponding temperature (around 230 °C) are in the required intervals to perform the HTC process: 180 °C ≤ T ≤ 250 °C and 10 bar ≤ P ≤ 30 bar. The solid phase behavior in the batch reactor during the HTC process is analyzed through the distribution and velocity of particles as plotted in Figure 9a-c. At 600s (Figure 9a), the particles are accelerated and randomly scattered within the domain. Next, Moreover, the reactor bottom is slightly hotter than its top, regardless of time. Moving away from the inner lateral surface, Figure 6 reveals that there is an average temperature difference of 5 • C between both base surfaces (bottom and top surfaces). This fact can be explained by suspension particles [28,29] and their migration to the reactor top (Figure 9a-c) owing to their low bulk density [2] as well as their thermal conductivity [23] compared to those of water. Figure 5a-c and also Figure 6 prove that there are wide variations of temperature gradients in the vicinity of the vertical wall (r = R), while they become small far enough from it ( Figure 6). In the same region, the fluid remains practically stagnant because of the imposed no-slip condition applied at the solid-fluid boundary (Figure 7a-c). Far beyond this zone (vertical wall) and also from the lower and upper surfaces of the mixture, natural convection due to the bulk movement of water molecules takes place, and hence it generates convection cells (Figure 7a-c). In the same way, another set of large cells is developed along the symmetry axis. Figure 8a-c display isobars at three different times (t = 600, 1800 and 3000 s). We observe that, at each time, autogenous pressure value remains constant with an increase from 1.8 bar to 28.09 bar at 600 and 3000 s, respectively. The last pressure value and its corresponding temperature (around 230 • C) are in the required intervals to perform the HTC process: 180 • C ≤ T ≤ 250 • C and 10 bar ≤ P ≤ 30 bar. Moreover, the reactor bottom is slightly hotter than its top, regardless of time. Moving away from the inner lateral surface, Figure 6 reveals that there is an average temperature difference of 5 °C between both base surfaces (bottom and top surfaces). This fact can be explained by suspension particles [28,29] and their migration to the reactor top (Figure 9a-c) owing to their low bulk density [2] as well as their thermal conductivity [23] compared to those of water. Figure 5a-c and also Figure 6 prove that there are wide variations of temperature gradients in the vicinity of the vertical wall (r = R), while they become small far enough from it ( Figure 6). In the same region, the fluid remains practically stagnant because of the imposed no-slip condition applied at the solid-fluid boundary (Figure 7a-c). Far beyond this zone (vertical wall) and also from the lower and upper surfaces of the mixture, natural convection due to the bulk movement of water molecules takes place, and hence it generates convection cells (Figure 7a-c). In the same way, another set of large cells is developed along the symmetry axis. Figure 8a-c display isobars at three different times (t =600, 1800 and 3000 s). We observe that, at each time, autogenous pressure value remains constant with an increase from 1.8 bar to 28.09 bar at 600 and 3000 s, respectively. The last pressure value and its corresponding temperature (around 230 °C) are in the required intervals to perform the HTC process: 180 °C ≤ T ≤ 250 °C and 10 bar ≤ P ≤ 30 bar. The solid phase behavior in the batch reactor during the HTC process is analyzed through the distribution and velocity of particles as plotted in Figure 9a-c. At 600s (Figure 9a), the particles are accelerated and randomly scattered within the domain. Next, The solid phase behavior in the batch reactor during the HTC process is analyzed through the distribution and velocity of particles as plotted in Figure 9a-c. At 600 s (Figure 9a), the particles are accelerated and randomly scattered within the domain. Next, due to their lower density in comparison to water, they move towards the top of the cylinder where they can be stratified (Figure 9b,c) with low velocities. Additionally, the drag force has also an exceptional contribution to this fact. These results confirmed those found by Mendecka et al. [17]. In all conditions, particles motion must be controlled for the smooth running of the HTC process.
To highlight the influence of heat flux imposed at the lateral surface area of reactor, three values of q 0 (2.0, 2.3 and 3.0 kW) are considered. The results presented in Figures 10a-c, 11a-c and 12a-c are specifically performed for the steady-state heat transfer regime (3000s after the onset of the HTC process), η = 1/3 and r = 2.7 cm. Isotherms (Figure 10a-c) and isobars (Figure 12a-c) show that temperature and pressure increase with the imposed thermal power to reach maximum values (about 260 • C and 28 bar). As a result, the temperature and pressure are greater than the lower limit conditions for the HTC process (180 • C and 10 bar) and hence the thermochemical conversion may be efficient. In Figure 11a-c, the velocity fields within the reactor are plotted in the same conditions. The plots demonstrate that the vortices formation and velocity values are slightly affected by the heating process.
Energies 2022, 15, x FOR PEER REVIEW 11 of 18 due to their lower density in comparison to water, they move towards the top of the cylinder where they can be stratified (Figure 9b,c) with low velocities. Additionally, the drag force has also an exceptional contribution to this fact. These results confirmed those found by Mendecka et al. [17]. In all conditions, particles motion must be controlled for the smooth running of the HTC process.
To highlight the influence of heat flux imposed at the lateral surface area of reactor, three values of q (2.0, 2.3 and 3.0 kW) are considered. The results presented in Figures  10a-c, 11a-c and 12a-c are specifically performed for the steady-state heat transfer regime (3000s after the onset of the HTC process), η = 1/3 and r = 2.7 cm. Isotherms (Figure  10a-c) and isobars (Figure 12a-c) show that temperature and pressure increase with the imposed thermal power to reach maximum values (about 260 °C and 28 bar). As a result, the temperature and pressure are greater than the lower limit conditions for the HTC process (180 °C and 10 bar) and hence the thermochemical conversion may be efficient. In Figure 11a-c, the velocity fields within the reactor are plotted in the same conditions. The plots demonstrate that the vortices formation and velocity values are slightly affected by the heating process.
Additionally, Figure 13a-c divulge the sensitivity particles' velocity to the thermal power (q ) variation. Increasing this last parameter accelerates motion of particles which could cause collision between them. This increase should provide, in principle, much faster conversion reactions.    Additionally, Figure 13a-c divulge the sensitivity particles' velocity to the thermal power (q 0 ) variation. Increasing this last parameter accelerates motion of particles which could cause collision between them. This increase should provide, in principle, much faster conversion reactions.

Kinetics Analysis
Kinetics models are often needed to design a HTC reactor with optimal operating conditions which allow the production of good quality hydrochar (HC). Roman et al. [30] reported that most studies on HTC process have been restricted to liquid phase, whereas the evolution of solid phase during this process did not receive sufficient attention. Within this context, Equations (8a)-(8c) can model kinetics behavior of two solids (ROP and HC) involved in the HTC reactions. These models could quantify the decomposed mass of ROP as well as that of the produced HC. As the HTC process happens in steady-state conditions, it was assumed that the mixture is at an average temperature T av for estimating the reaction rate k through Equation (8c). Furthermore, each temperature (T av = 200, 230 and 260 • C) corresponds to one q 0 (2.0, 2.3 and 3.0 kW) respectively. Accordingly, the solid mass conversion (m ROP ), produced mass of hydrochar (m HC ) and their respective kinetics rates are assessed by using Equations (8a) and (8b). At the steady-state regime, the average temperature T av is obtained from temperature results computed using the CFD simulation explained above. For an initial mass, (m 0 = 420 g), of raw olive pomace (ROP), Figure 14a-c depicts simultaneously graphical curves giving instantaneous variation of the solid mass (m ROP ) conversion, produced mass of hydrochar (m HC ) and their respective kinetics rates.
As expected, the quantity of product (HC) increases while biomass amount (raw olive pomace) and the reaction rates decrease regardless of the mean batch reactor temperature (resp. the imposed heat flux). We emphasize also that the increase of these last important factors (T av = 200, 230 and 260 • C and q 0 = 2.0, 2.3 and 3.0 kW) accelerates the conversion reaction and, therefore, reduces the residence time of the product (ROP) within the reactor (180, 160 and 130 min). Moreover, the model proves that about 150, 160 and 200 g of hydrochar was produced from m 0 = 420 g of biomass for 200 • C (respectively 2.0 kW), 230 • C (respectively 2.3 kW) and 260 • C (respectively 3.0 kW) all the more so their residence times are very different.

Kinetics Analysis
Kinetics models are often needed to design a HTC reactor with optimal operating conditions which allow the production of good quality hydrochar (HC). Roman et al. [30] reported that most studies on HTC process have been restricted to liquid phase, whereas the evolution of solid phase during this process did not receive sufficient attention. Within this context, Equations (8a)-(8c) can model kinetics behavior of two solids (ROP and HC) involved in the HTC reactions. These models could quantify the decomposed mass of ROP as well as that of the produced HC. As the HTC process happens in steady-state conditions, it was assumed that the mixture is at an average temperature for estimating the reaction rate k through Equation (8c). Furthermore, each temperature (T = 200, 230 and 260°C) corresponds to one q (2.0, 2.3 and 3.0 kW) respectively. Accordingly, the solid mass conversion (m ), produced mass of hydrochar (m ) and their respective kinetics rates are assessed by using Equations (8a) and (8b). At the steady-state regime, the average temperature T is obtained from temperature results computed using the CFD simulation explained above.
For an initial mass, (m = 420 g), of raw olive pomace (ROP), Figure 14a-c depicts simultaneously graphical curves giving instantaneous variation of the solid mass (m ) conversion, produced mass of hydrochar (m ) and their respective kinetics rates. As expected, the quantity of product (HC) increases while biomass amount (raw olive pomace) and the reaction rates decrease regardless of the mean batch reactor temperature (resp. the imposed heat flux). We emphasize also that the increase of these last important factors ( T = 200, 230 and 260 °C and q = 2.0, 2.3 and 3.0 kW) accelerates the conversion reaction and, therefore, reduces the residence time of the product (ROP) within the reactor (180, 160 and 130 min). Moreover, the model proves that about 150, 160 and 200 g of hydrochar was produced from m =

Conclusions
In this work, CFD modeling of the HTC process of raw olive pomace (ROP) was conducted to predict fluid flow and thermal behavior within a batch reactor as well as the conversion of this biomass into hydrochar. This relatively new approach considered the mixture (distilled water +ROP) as a dispersed medium, which is analyzed according to the Eulerian-Lagrangian approach. The biomass (ROP) degradation was modeled using a basic kinetics model.
After the model validation, the main results were presented at the position r = 2.7 for three different times after the onset of HTC process namely isotherms, velocity field and isobars within the batch reactor. For q 0 = 2.0 kW and η = 1/3, it has been shown that far from the symmetry axis, isotherms and velocity fields are symmetrical with respect to this axis. Also, the bottom of the batch reactor is hotter than the top, regardless of time. As for the solid phase, it has been observed that particles are accelerated and randomly scattered within the domain during heating up period (t = 600 s). Next, they move towards the top of the cylinder where they can be stratified, and under all conditions, particles motion must be controlled for the smooth running of the HTC process.
The influence of q 0 on the thermal, dynamic and particles behavior inside the reactor was highlighted for the steady-state heat transfer regime (3000 s after the onset of the process). Isotherms and isobars profiles showed that temperature and pressure increase with the imposed thermal power to reach maximum values (about 260 • C and 28 bar). Thus, the temperature and pressure are above the lower limit conditions for the HTC process (180 • C and 10 bar) and hence the thermochemical conversion may be efficient. Moreover, the sensitivity of the velocity of particles to the thermal power (q 0 ) variation was examined. Increasing this last parameter accelerates motion of particles which could cause collision between them. This increase should provide much faster conversion reactions.
From reaction kinetics point of view, curves of biomass degradation and produced hydrochar with respect to time have an exponential form. As expected, the hydrochar (HC) produced increases while biomass amount (raw olive pomace) and the reaction rates decrease regardless of the mean batch reactor temperature (with respect to the imposed heat flux). Moreover, increasing these important factors (T av = 200, 230 and 260 • C and q 0 = 2.0, 2.3 and 3.0 kW) promotes the conversion reaction and therefore reduces the residence time of the product (ROP) within the reactor (180, 160 and 130 min).
In conclusion, the findings obtained through this CFD modeling could be used to design and engineer a HTC reactor with optimal operating conditions.

Conflicts of Interest:
We declare that this manuscript is original and also no conflict of interest associated with it.

Appendix A
The motion of particles is described by fundamental principle of dynamics (see Equation (4)) where several forces are involved especially: → F vm and → F PG which represent body force, drag force, Basset force, virtual mass force and pressure gradient force respectively [31,32]. These forces are defined below.

Appendix A.1. Body Force
This is the gravitational force defined by: This force depends on the particle size and it is due to the pressure gradient. Therefore,

.2. Virtual Mass Force
It is an unsteady force which has an important role in the dynamics of the mixture (particles + water). It represents the mass of liquid carried by the particle (C M ρ w τ p ), times relative acceleration between phases [32].
C M = 0.5, τ p and D Dt (.) are the virtual mass coefficient, the particle volume and the material derivative, respectively.
The presence of the Lagrangian derivative of velocity further complicates motion equation solving. Therefore, it was proposed to use the following relation [32]: The Basset force, known as the 'history term', describes the force due to boundary layer (effect of the viscosity) development. It is expressed by [32]: In the case of a spherical particle, the drag force is expressed by: The standard drag coefficient is that developed by Schiller and Naumann [32]: The flow regime around a particle can be identified from the particle Reynolds number (Rep): It should be noted that Equation (A5) is valid for the flow around a single particle in an infinite fluid medium. However, the empirical correlation (A8) developed in [32] takes into account the interactions between particles: