Experimental Investigation and CFD Analysis of Pressure Drop in an ORC Boiler for a WHRS Implementation

Waste heat dissipated in the exhaust system of a combustion engine represents a major source of energy to be recovered and converted into useful work. The Waste Heat Recovery System (WHRS) based in an Organic Rankine Cycle (ORC) is an approach for recovering energy from heat sources, achieving a significant reduction in fuel consumption and, as a result, exhaust emissions. This paper studies pressure drop in an ORC shell-and-tubes boiler for a WHRS implementation experimentally and with computational simulations based on a 1-dimensional heat transfer model coupled with 3D calculations. An experimental database is developed, using ethanol in a pressure range of 10–15 absolute bar as working fluid, with mass fluxes inside the tubes in the range of 349.31 kg/s-m2 and 523.97 kg/s-m2, and inlet temperatures in the range of 60 °C and 80 °C. Thus, the friction factor of different regions of the boiler were estimated using both CFD simulations, experimental data, and bibliographic correlations. Simulations of operating points and the results of the experimental test bench showed good agreement in pressure drop results, with a mean absolute error of 15.47%, without a significant increment in the computational cost.

The ORC based on vapor generation in a secondary circuit represents an indirect method of waste heat recovery (WHR). This technique has advantages compared with the so-called direct WHR techniques, such as electrical turbo compounding [52], or mechanical turbo compounding [53], which uses a power turbine fitted to the vehicle exhaust and has a much higher impact on the engine pumping losses. In addition, an ORC achieves a high waste energy conversion and it is cheaper than other WHR techniques such as thermoelectric generators [2]. a well-known issue that most CFD software face convergence difficulties in multi-phase calculations, such as numerical instability in non-structured meshes, and higher computational cost. This fact is mainly due to an increase in the number of equations involved in the calculations. This limitation can be tackled in a nontrivial way by creating a mesh with an approximate grid size of the physical processes involved [78,79]. Nevertheless, the geometry of the industrial components is often excessively complex to afford a refined mesh that considers multi-phase fluid dynamics. Consequently, a multiphase CFD simulation would lead to an enormous computational effort and an increase in the calculation time, compromising or even preventing fast responses to the engineering industry.
There are not many studies in the literature that study evaporation by means of CFD simulations in heat exchangers due to the difficulties and challenges previously described. Mali et al. [57] studied pressure drop in vertical tubes for a boiler at high pressure, using the Eulerian-Eulerian model and dividing the length of the tube into different sections regarding vapor fraction. Wang et al. [80] and Shi et al. [81,82] used de RPI wall boiling model proposed by Kurul and Podowski [83] to simulate a particular region of an evaporator, taking advantage of the periodic distribution of the tubes inside the shell to predict the thermal-hydraulic characteristics of the tube. Mohammed et al. [84] used the VOF model to simulate the evaporation and condensation of the acetone in a horizontal circular tube. Cappelli et al. [85] used the VOF model to investigate gas flow rate, gas temperature, and liquid flow rate on mass transfer in a bubble column evaporator. Höhne [86] used the homogeneous multiphase model to simulate in a CFD software the heat and mass transfer process and compared the results with experiments from the literature. The studies reviewed have in common that neither specifically evaluate WF pressure drop via CFD, and they simulate simple geometries consisting in one straight tube or a specific region of the heat exchanger, and they are focused on studying mass transfer or thermal properties.
For all the above, there is an increasing interest in the study of pressure drop in multiphase flow in heat exchangers and the challenges involved in CFD simulations of multiphase flow in complex geometries. This paper presents a methodology to predict multiphase pressure drop in shell and tube heat exchangers by means of CFD simulation and our own experimental data. The limitations of CFD in the field of multiphase flow are addressed using a well-contrasted tool to couple a 1-dimensional model to a CFD software, presented by the authors of the present paper previously in 2019 [87]. The correlations to predict pressure drop used in the 1-dimensional model are based on bibliographic references for single-phase flow and multiphase flow, the experimental results obtained in our own experimental test bench, and full CFD simulations are used to model the local pressure drop factor of the boiler studied.
The present paper aims to study one of the main components of the ORC, the boiler, both experimentally and numerically. This work continues the studies published in a previous paper by the authors [87] where a 1-dimensional model was coupled to a computational fluid dynamics (CFD) commercial software in a boiler of an ORC in a WHRS. This previous paper focused on understanding heat transfer experimentally in a test bench via infrared thermography and using CFD software. However, the present paper is focused on analyzing and modeling the pressure drop of the multiphase flow in the WF in a boiler for an ORC by means of CFD simulations and experimental results, and proposing a methodology to estimate the local pressure drop factor in single-phase and multiphase flow.
In Section 2 of this paper, the test bench setup and the boiler of the study are presented. In Section 3, the experimental results are shown. Section 4 describes the meshing process. Section 5 explains the numerical approach. Section 6 presents the numerical and experimental methodology to obtain the friction factor and local pressure drop factor in the single-phase and multiphase stages of the WF inside the boiler. The CFD results are discussed and contrasted with experimental data in Section 7. Finally, conclusions regarding the methodology presented in this study are recapitulated in Section 8. Figure 1 depicts the experimental test bench used in the present work, which allows analysis of the boilers and evaporators in a broad range of mass flows, temperatures, and pressures. The test bench has three main circuits: air, ethanol, and cooling water. The ethanol accumulated in the tank is pumped and preheated before its entrance to the evaporator, where it receives heat from hot air, previously heated with an electric resistor. Subsequently, ethanol is cooled down in the condenser by transferring its heat to a cooling water circuit.  Figure 1 depicts the experimental test bench used in the present work, which allows analysis of the boilers and evaporators in a broad range of mass flows, temperatures, and pressures. The test bench has three main circuits: air, ethanol, and cooling water. The ethanol accumulated in the tank is pumped and preheated before its entrance to the evaporator, where it receives heat from hot air, previously heated with an electric resistor. Subsequently, ethanol is cooled down in the condenser by transferring its heat to a cooling water circuit. Figure 1. Schematic diagram of the test bench. The gas, WF, and water refrigerant circuits are represented in red, blue, and green, respectively. (T) Denotes temperature sensors, (p) absolute pressure sensors, and (dp) differential pressure sensors.

Experimental Setup
There are multiple temperature and pressure sensors in several points of interest in the test bench. Thus, temperatures, pressures, and mass flows can be measured. These points of interest include the evaporator inlet and outlet, and the ethanol pre-heater inlet and outlet, among others. K-type thermocouples were disposed for the gas side temperature measurement with an uncertainty of 1.5 °C, while T-type thermocouples were disposed for the WF and the metal sides with an uncertainty of 0.5 °C. The sensors used for pressure drop measurements were absolute pressure transmitters, which featured a single crystal silicon resonant sensor, with an accuracy span of ±0.055%, and a differential pressure sensor, which consists of a piezo-resistance silicon-type sensor with an accuracy of 0.25%. Schematic diagram of the test bench. The gas, WF, and water refrigerant circuits are represented in red, blue, and green, respectively. (T) Denotes temperature sensors, (p) absolute pressure sensors, and (dp) differential pressure sensors.
There are multiple temperature and pressure sensors in several points of interest in the test bench. Thus, temperatures, pressures, and mass flows can be measured. These points of interest include the evaporator inlet and outlet, and the ethanol pre-heater inlet and outlet, among others. K-type thermocouples were disposed for the gas side temperature measurement with an uncertainty of 1.5 • C, while T-type thermocouples were disposed for the WF and the metal sides with an uncertainty of 0.5 • C. The sensors used for pressure drop measurements were absolute pressure transmitters, which featured a single crystal silicon resonant sensor, with an accuracy span of ±0.055%, and a differential pressure sensor, which consists of a piezo-resistance silicon-type sensor with an accuracy of 0.25%.
The boiler used to validate the model is a flow shell and corrugated-tube heat exchanger in crossed flow and counter flow, made entirely of steel AISI 316l, with the excep- tion of the tubes, made of steel AISI 304. Air flows perpendicularly through a bundle of staggered tubes, while the WF flows inside those tubes. After the WF intake, WF mass flow is divided in rows of 14 tubes that converge in a mixing chamber at opposite sides of the boiler. In that mixing chamber, the WF converges and is also divided in another row of 14 tubes. This pattern is repeated in a total of 50 steps of rows of tubes inside the shell of the boiler, accounting for a total of 700 corrugated tubes of 5 mm of diameter and 0.25 mm of wall thickness. After the 50th row of tubes, the WF mass flow converges in the WF outtake and leaves the boiler into the WF circuit of the test bench. A schematic representation of a mixing chamber is shown in Figure 2, where each of the tubes displayed represents all the 14 tubes belonging to the same row. The boiler used to validate the model is a flow shell and corrugated-tube heat exchanger in crossed flow and counter flow, made entirely of steel AISI 316l, with the exception of the tubes, made of steel AISI 304. Air flows perpendicularly through a bundle of staggered tubes, while the WF flows inside those tubes. After the WF intake, WF mass flow is divided in rows of 14 tubes that converge in a mixing chamber at opposite sides of the boiler. In that mixing chamber, the WF converges and is also divided in another row of 14 tubes. This pattern is repeated in a total of 50 steps of rows of tubes inside the shell of the boiler, accounting for a total of 700 corrugated tubes of 5 mm of diameter and 0.25 mm of wall thickness. After the 50th row of tubes, the WF mass flow converges in the WF outtake and leaves the boiler into the WF circuit of the test bench. A schematic representation of a mixing chamber is shown in Figure 2, where each of the tubes displayed represents all the 14 tubes belonging to the same row.

Experimental Results
Two sets of experiments were performed: one with very low heat exchange, with air below the saturation temperature to ensure that the WF remains as liquid and to analyze pressure drop in the monophasic stage. The other experiment tested the operating conditions of the boiler to analyze the pressure drop in the multiphase stage. In every case, the fluids involved were ethanol as WF flowing inside the tubes and dry air at atmospheric pressure flowing outside the tubes and inside the shell of the boiler in perpendicular direction to the tubes.
The first set of experiments used to study monophasic pressure drop are listed in Table 1. Note that all the inputs in the test bench are fixed, except mass flow of the WF.

Experimental Results
Two sets of experiments were performed: one with very low heat exchange, with air below the saturation temperature to ensure that the WF remains as liquid and to analyze pressure drop in the monophasic stage. The other experiment tested the operating conditions of the boiler to analyze the pressure drop in the multiphase stage. In every case, the fluids involved were ethanol as WF flowing inside the tubes and dry air at atmospheric pressure flowing outside the tubes and inside the shell of the boiler in perpendicular direction to the tubes.
The first set of experiments used to study monophasic pressure drop are listed in Table 1. Note that all the inputs in the test bench are fixed, except mass flow of the WF. The second set of experiments was performed in operating conditions, where the WF enters the boiler as subcooled liquid and exists as an overheated vapor. Multiphase experiments were studied in several rounds of experiments, where WF inlet conditions were fixed and dry air conditions started at 500 • C and 50 kg/h and were gradually increased 10 • C and 1 kg/s every 20 min and at every measurement of the pressure drop. Outlet temperature was monitored and, if after 20 min the outlet temperature did not surpass the saturation temperature and provided the experiment had already reached a steady stage, the experimental point was dismissed and the inlet air temperature and inlet air mass flow were increased again. Only when the outlet temperature sensor showed a measurement significantly higher than the saturation temperature was it assumed that the WF had reached the stage of overheated vapor. The sensors exported measurements each second using a data acquisition system. Pressure drop final results are the average of the pressure drop measurements in the last 5 min of each experiment. The list of experiments analyzed in multiphase operating conditions and their results are shown in Table 2, where the mean average error (MAE) denotes the error between each measurement in every time step and the average pressure drop result. Considering the MAE in WF pressure drop results in Table 2, a moderate error has been obtained, considering the intermittent nature of multiphase flow, enhanced by the complex geometry in the WF side of the boiler studied. Before addressing predictive tools, it is important to point out some of the anomalies that influence two-phase flow behavior and, therefore, pressure drop. As mentioned earlier, a relatively high-pressure drop in micro channels can result in significant property changes, particularly specific volume and enthalpy of the individual phases, which lead to uncertainties and instabilities in the experimental results obtained and displayed in Table 2. These results contrast with the single-phase experiments shown in Table 1, whose experimental results were straightforward to obtain in the test bench and showed quite a linear response. In multiphase experiments, sensors occasionally showed sudden offsets in the measurements of pressure drop. This phenomenon is attributed to the number of mixing chambers of the boiler, which enhances the formation of vapor pockets in these low-velocity regions. When the vapor pockets are detached from the walls of the boiler, it is possible to modify the pressure drop measurements. This issue will be tackled in future work by inducing vibrations in the boiler to prevent vapor pockets from accumulating in the mixing chamber.

Meshing
The surfaces of the geometry, created in CATIA v5, were discretized with triangular elements, with an average size of the edge of each element of 1 mm, whereas the volume of the elements were meshed using ANSYS Fluent. The volume of the solid tubes was created by sweeping the triangular surface elements into prisms in the normal direction of each surface. In the gas side, to better model the near wall region [80], the triangular elements in the surface were swept into 10 prisms in the normal direction of each element of the surface, with the first element having a thickness of 0.015 mm with a growth ratio of every subsequent element of 1.15, to model the boundary layer. The meshing scheme used to fill the rest of the gas volume was Hexcore, filling the gaps between hexahedral elements and prisms with tetrahedral elements. The resulting mesh for air and metal domains sums a total element of around 60,000,000 volume elements, with a maximum skewness of 0.90 and volume average skewness of 0.173. Mesh element size was limited by tube distribution, with a narrow gap of about 2.5 mm between one another. As the ANSYS Fluent manual [76] suggests, a boundary layer should be meshed in the wall to properly model the law-of-thewall. A mesh configuration of 0.5 mm triangles was shown to be adequate for the geometry of the boiler since a coarser mesh generated many high skewness elements that could lead to instabilities or even divergence of the mesh. This mesh size and quality were already validated in a previous work [87] and did not represent a huge computational effort. A detail of the CFD mesh is displayed in Figure 3. The 1-dimensional WF mesh was created with an external routine created with the Used Defined Functions of ANSYS Fluent, which will be further explained in the following section of the present paper. of each surface. In the gas side, to better model the near wall region [80], the triangular elements in the surface were swept into 10 prisms in the normal direction of each element of the surface, with the first element having a thickness of 0.015 mm with a growth ratio of every subsequent element of 1.15, to model the boundary layer. The meshing scheme used to fill the rest of the gas volume was Hexcore, filling the gaps between hexahedral elements and prisms with tetrahedral elements. The resulting mesh for air and metal domains sums a total element of around 60,000,000 volume elements, with a maximum skewness of 0.90 and volume average skewness of 0.173. Mesh element size was limited by tube distribution, with a narrow gap of about 2.5 mm between one another. As the ANSYS Fluent manual [76] suggests, a boundary layer should be meshed in the wall to properly model the law-of-the-wall. A mesh configuration of 0.5 mm triangles was shown to be adequate for the geometry of the boiler since a coarser mesh generated many high skewness elements that could lead to instabilities or even divergence of the mesh. This mesh size and quality were already validated in a previous work [87] and did not represent a huge computational effort. A detail of the CFD mesh is displayed in Figure 3. The 1-dimensional WF mesh was created with an external routine created with the Used Defined Functions of ANSYS Fluent, which will be further explained in the following section of the present paper.

Numerical Model
The CFD software used to run the simulations of the present paper was ANSYS FLU-ENT 2022R2. The 1-dimensional model was programmed in an ad hoc algorithm and implemented using User Defined Functions (UDF) [88], which allows modification of their default heat transfer laws, among many other functionalities. As previously mentioned, multiphase CFD solvers have some limitations: difficulty to model heat and mass transfer mechanisms locally in many flow regimes, great computational costs, and numerical instabilities. The methodology in which this paper is based [87] was conceived as an effort to these limitations while still benefiting from the advantages of CFD software. The general approach, summarized in Figure 4, consists of the coupling of two computational

Numerical Model
The CFD software used to run the simulations of the present paper was ANSYS FLUENT 2022R2. The 1-dimensional model was programmed in an ad hoc algorithm and implemented using User Defined Functions (UDF) [88], which allows modification of their default heat transfer laws, among many other functionalities. As previously mentioned, multiphase CFD solvers have some limitations: difficulty to model heat and mass transfer mechanisms locally in many flow regimes, great computational costs, and numerical instabilities. The methodology in which this paper is based [87] was conceived as an effort to these limitations while still benefiting from the advantages of CFD software. The general approach, summarized in Figure 4, consists of the coupling of two computational domains, namely, the gas and metal side. Three-dimensional meshes are generated using the ANSYS meshing functionalities, and the WF side, where a one-dimensional mesh is generated and coupled with its metal frontier by specifically developed external algorithms. domains, namely, the gas and metal side. Three-dimensional meshes are generated using the ANSYS meshing functionalities, and the WF side, where a one-dimensional mesh is generated and coupled with its metal frontier by specifically developed external algorithms. The first step of the methodology implementation consists in discretizing the flow along the tubes in ordered cells into a 1-dimensional mesh flow along the WF path. This is achieved by programming a UDF macro called DEFINE_ON_DEMAND, which is a general-purpose macro that you can use to specify a UDF that is executed "on demand" in ANSYS Fluent, rather than having ANSYS Fluent call it automatically during the calculation. A 1-dimensional cell index is assigned to each CFD cell adjacent to the WF, based in its coordinates, and this cell index is stored using a user-defined memory. The calculations regarding this 1-dimensional mesh are not calculated by ANSYS Fluent itself, but with another UDF macro that updates the values of heat transfer coefficient and bulk temperature in each iteration by programming a DEFINE_EXECUTE_AT_END macro, which contains the physical model described in a previous work [87] and in the correlations described in Section 6 of this paper. Another function of this macro is to compute the increment of the enthalpy of each 1-dimensional cell by accounting the total energy absorbed in it through the walls adjacent to the WF fluid domain.
The DEFINE_EXECUTE_AT_END macro stores the heat transfer coefficient and bulk in another user-defined memory, and this information is transferred to ANSYS Fluent Solver using a DEFINE_PROFILE macro, which is implemented in the boundary adjacent to the WF side. This process is updated in each iteration of ANSYS Fluent solver, coupling the 1-dimensional and 3-dimensional calculations, and is repeated until there is convergence in the solution. The thermodynamic properties of the WF were obtained using NIST REFPROP.
The total WF mesh size in this simulation has 35,000 elements, with a total of 50 elements in each one of the 700 tubes. In a preliminary grid independence study, this 1-dimensional mesh size was proven to be acceptable, since finer grid elements generated oscillations in the boundary between phases and in regions with high wall heat flux gradients. On the other hand, coarser elements led to uncertainties in the onset of boiling. In any case, the 35,000 elements of the 1-dimensional mesh meant very little computational effort, and a finer mesh was shown to be a bigger concern than a coarser one. The first step of the methodology implementation consists in discretizing the flow along the tubes in ordered cells into a 1-dimensional mesh flow along the WF path. This is achieved by programming a UDF macro called DEFINE_ON_DEMAND, which is a general-purpose macro that you can use to specify a UDF that is executed "on demand" in ANSYS Fluent, rather than having ANSYS Fluent call it automatically during the calculation. A 1-dimensional cell index is assigned to each CFD cell adjacent to the WF, based in its coordinates, and this cell index is stored using a user-defined memory. The calculations regarding this 1-dimensional mesh are not calculated by ANSYS Fluent itself, but with another UDF macro that updates the values of heat transfer coefficient and bulk temperature in each iteration by programming a DEFINE_EXECUTE_AT_END macro, which contains the physical model described in a previous work [87] and in the correlations described in Section 6 of this paper. Another function of this macro is to compute the increment of the enthalpy of each 1-dimensional cell by accounting the total energy absorbed in it through the walls adjacent to the WF fluid domain.
The DEFINE_EXECUTE_AT_END macro stores the heat transfer coefficient and bulk in another user-defined memory, and this information is transferred to ANSYS Fluent Solver using a DEFINE_PROFILE macro, which is implemented in the boundary adjacent to the WF side. This process is updated in each iteration of ANSYS Fluent solver, coupling the 1-dimensional and 3-dimensional calculations, and is repeated until there is convergence in the solution. The thermodynamic properties of the WF were obtained using NIST REFPROP.
The total WF mesh size in this simulation has 35,000 elements, with a total of 50 elements in each one of the 700 tubes. In a preliminary grid independence study, this 1-dimensional mesh size was proven to be acceptable, since finer grid elements generated oscillations in the boundary between phases and in regions with high wall heat flux gradients. On the other hand, coarser elements led to uncertainties in the onset of boiling. In any case, the 35,000 elements of the 1-dimensional mesh meant very little computational effort, and a finer mesh was shown to be a bigger concern than a coarser one.
Note that the numerical methodology described is independent from the physical model chosen. The only requirement for implementing this methodology is that the WF flow can be parametrized into a 1-dimensional path.
Regarding the particularities of this geometry, namely, those regions that cannot be easily parametrized, in the mixing chambers and the inlet and outlet manifolds, perfect Sensors 2022, 22, 9437 9 of 21 mixture was assumed, averaging the enthalpies of the paths of fluids involved. Inlet enthalpies in any row of tubes that departs from one mixing chamber and converges in the following mixing chamber (referred to this paper as a row of tubes), is the mass flow average outlet enthalpy of the immediately previous row of tubes.
To model the mass flow distribution in tubes belonging to the same row, the same pressure drop was assumed for each row. Since not every tube in the same row of tubes absorbs the same amount of heat from the air and, as consequence, the WF in each tube of the same row is thermodynamically different from one another, the WF mass flow in each individual tube is modified in order to adjust the pressure drop in each individual tube to achieve a common pressure drop in every tube of the same row.
The numerical approach to estimate mass flow distribution consists in calculating, for every tube in every single row of tubes, the pressure drop of those outcomes with every mass flow possibility within a reasonable range. Having the pressure drop per tube as a function of its mass flow, the target is to find a combination of pressure drop and mass flow per tube that accomplishes the conditions of Equations (1) and (2) at the same time: where . m i stands for the mass flow of each tube of a row of tubes and . m T is the total mass flow in a column of tubes and the total mass flow of WF in the boiler; ∆p . m i is the targeted solution that stands for the pressure drop in a single tube that must be the same in every tube in the same row for parallel tubes, regardless of the mass flow in each tube. This equation was solved using the secant method (Equation (3)): Note that this will always find one solution in the single-phase stage, since the Darcy-Weisbach correlation (later discussed in Equation (4)) is linear. However, multiphase pressure drop nature is not linear, hence the secant method might find multiple solutions depending on the initial conditions, particularly in the case where some tubes started the multiphase stage and the rest remain in the subcooled liquid stage. In that case, every solution is recorded, and the lowest pressure solution of all of them is implemented in the CFD calculations. Note that the existence of multiple solutions could be translated into instabilities in the physics of the evaporation process. This scenario was not found in this study, but the methodology considers this possibility.

Pressure Drop Modeling
The general approach to obtain a predictive tool for the pressure drop in the boiler studied is to combine bibliographic correlations, experimental results, and CFD simulations. While bibliographic correlations for pressure drop are usually designed for simple and 1-dimensional geometries, the boiler studied, as already described, has several mixing chambers, which will contribute locally in the total pressure drop of the boiler. Hence, implementing only 1-dimensional bibliographic correlations in the CFD simulations will lead to an underprediction of the total pressure drop. This local pressure drop is estimated by calculating the difference between experimental pressure drop and CFD simulations using only 1-dimensional correlations.
In this section of the present paper, the proposed methodology of pressure drop modeling is broken down in two main steps: Single-phase pressure drop modeling and multiphase pressure drop modeling. Bibliographic correlation, numerical considerations, and results for local pressure drop factor are described.

Single-Phase Pressure Drop
The experiments shown in Table 1 were simulated in a full CFD single-phase calculation. The mesh described in Section 4 was modified to include a CFD 3-dimensional mesh to compute the WF domain as well. The resulting mesh has about 95 million cells, with approximately 500,000 cells each tube. A detail of the CFD of each tube is shown in Figure 5. A mesh convergence study was performed showing that the simulations reached convergence with y+ lower than 0.15 which, in the case of the meshes studied, means that the first boundary layer has a height of 0.015 mm, with 18 boundary layers and a surface mesh of triangles with an edge length of 1 mm. In this part of the study, a standard CFD simulation was launched, without using external routines.
In this section of the present paper, the proposed methodology of pressure drop modeling is broken down in two main steps: Single-phase pressure drop modeling and multiphase pressure drop modeling. Bibliographic correlation, numerical considerations, and results for local pressure drop factor are described.

Single-Phase Pressure Drop
The experiments shown in Table 1 were simulated in a full CFD single-phase calculation. The mesh described in Section 4 was modified to include a CFD 3-dimensional mesh to compute the WF domain as well. The resulting mesh has about 95 million cells, with approximately 500,000 cells each tube. A detail of the CFD of each tube is shown in Figure 5. A mesh convergence study was performed showing that the simulations reached convergence with y+ lower than 0.15 which, in the case of the meshes studied, means that the first boundary layer has a height of 0.015 mm, with 18 boundary layers and a surface mesh of triangles with an edge length of 1 mm. In this part of the study, a standard CFD simulation was launched, without using external routines. In Table 3, the results obtained with these simulations are shown using the nomenclature of the experiments previously detailed in Table 1. Comparing these 5 experiments with the CFD results in standard simulations, a mean average error (MAE) of 9.27% was obtained. This relatively good agreement between CFD capabilities and experimental results is a good starting point to adjust the bibliographic correlations to the experiments. In addition, full CFD pressure drop results can be broken down in each geometrical component, thus the contribution on the total pressure drop of the tubes and the mixing chambers can be quantified. In Table 3, the results obtained with these simulations are shown using the nomenclature of the experiments previously detailed in Table 1. Comparing these 5 experiments with the CFD results in standard simulations, a mean average error (MAE) of 9.27% was obtained. This relatively good agreement between CFD capabilities and experimental results is a good starting point to adjust the bibliographic correlations to the experiments. In addition, full CFD pressure drop results can be broken down in each geometrical component, thus the contribution on the total pressure drop of the tubes and the mixing chambers can be quantified. The objective of the implementation in the single-phase stage of the CFD coupled with the 1-dimensional model is the Darcy-Weisbach equation [89], where f is the friction factor and the parameter to be determined in the present section of this study. In laminar flow, the friction factor is given by the well-established Equation (5), which is a consequence of the Navier-Stokes equations.
Note that, although Equation (5) is used for circular section tubes, the Nikuradse diagrams predict that, for a low Reynolds number (laminar flow), friction factor is not a function of relative roughness. This statement has been proven true for corrugated tubes in preliminary simulations of this study. As expected, according to the Nikuradse diagram, for turbulent flow the friction factor is higher in a corrugated tube than in the equivalent smooth tube of the same diameter and Reynolds number. As a consequence, the next step in this study is the evaluation of the friction factor in the corrugated tube studied for turbulent flow. This was obtained by accounting for the pressure drop contribution of the tubes in CFD simulations. To this regard, a series of 20 CFD simulations were launched using 20 bar liquid ethanol as WF at 80 • C, with a mass flux in the tubes in a range from 12.48 kg/m 2 s to 449.11 kg/m 2 s. Taking only the simulations where the WF reached a turbulent regimen, and using a linear regression analysis, the correlation for monophasic friction factor of the corrugated tubes inside the boiler for single-phase turbulent flow ( f turb,s f ) obtained is shown in Equation (6) Once a friction factor correlation is obtained for the tubes, the next step of this study is to model the local pressure drop factor in the single-phase stage, k s f . The analytical expression for pressure drop of the WF in the boiler, including local pressure drop considerations, is shown in Equation (7): For the sake of simplicity, it is be assumed that this local pressure drop factor is the same for inlets and outlets where, in the case of the boiler of study, it has a total of 50 inlets in the tubes and 50 outlets from the tubes. Then, the expression for the local pressure drop factor is shown in Equation (8): Observing the CFD results of the experiments from Table 3, and taking into account the laminar and turbulent friction factor in the tubes discussed in Equations (5) and (6), the single-phase local pressure drop factor k s f in the mixing chambers was estimated with a linear regression analysis using MATLAB. Thus, a correlation for the local pressure drop factor in single-phase flow was obtained and is shown in Equation (9), where the values of a, b, and c coefficients are shown in Table 4: The correlation described in Equation (9) and in Table 4 for single-phase local pressure drop factor shows good agreement with experimental data, as is shown in the graph of Figure 6.

Re < 4000
Re > 4000 a −394.7 The correlation described in Equation (9) and in Table 4 for single-phase local pressure drop factor shows good agreement with experimental data, as is shown in the graph of Figure 6.

Multiphase Pressure Drop
The multiphase stage starts when the subcooled WF absorbs enough heat to reach saturation temperature. The progressive evaporation causes axial acceleration of the flow which increases both the wall shear stress and pressure gradient along the tube. There are

Multiphase Pressure Drop
The multiphase stage starts when the subcooled WF absorbs enough heat to reach saturation temperature. The progressive evaporation causes axial acceleration of the flow which increases both the wall shear stress and pressure gradient along the tube. There are a number of authors that studied pressure drop in multiphase flow in mini/micro channels, as previously discussed in the introduction. In this paper, the bibliographic correlation chosen to model multiphase pressure drop was proposed by Kim and Mudawar [68]. This correlation has been obtained from a big database of different fluids in a wide range of pressures and mass fluxes, that includes the operating conditions of the boiler of study. The correlation of Kim and Mudawar for pressure drop evaporating mini/micro channels is summarized in Table 5.
For the correlation for the pressure drop, the authors suggest in their paper that in the case of dryout, the correlation is no longer valid. The same authors, Kim and Mudawar, also studied and modeled the onset of dryout in mini/micro channels [90]. They gathered a database of experiments and correlations and proposed a new one based on statistical analysis. This correlation, shown in Equation (10), was used in the present paper to determine the vapor quality for dryout incipience (χ di ).

1.09
A dryout pressure drop correlation has not been found in the literature, so the pressure drop in this stage is obtained by blending the pressure drop for the unit of length at the point of dryout incipience and the pressure drop for the unit of length at the point of saturated vapor. The mathematical expression of this blending is described in Equations (11) and (12), where γ dryout stands for the blending factor, χ is the vapor quality, and dp dx go and dp dx g stand for the pressure drop gradient of saturated vapor and the pressure drop gradient of the vapor phase, respectively. dp dx dryout = γ dryout dp dx go + 1 − γ dryout dp dx g (11) When the correlations for the pressure drop shown in Equations (1)- (12) and Table 5 are implemented in preliminary CFD simulations and compared with the experimental data displayed in Table 2, it is clear from the simulation results that the model described so far under predicts the pressure drop in the boiler in the WF side. This outcome is attributed to the dynamics of the multiphase flow in the mixing chambers, where turbulence is severely enhanced and, as consequence, the interfacial shear stress between the liquid and the vapor phases increases. Thus, the single-phase pressure drop calculated previously is not valid to predict the pressure drop as well as the local pressure drop in multiphase flow. In order to numerically model this multiphase local pressure drop, a factor k m f is added to the total pressure drop model (Equation (13)): where n stands for the total number of cells of the path of the flow through the boiler, and m stands for the total number of singularities of the path of the WF inside the boiler, with this last one accounting for all the inlets and outlets of the tubes, including mixing chambers and the inlets and outlets of the WF manifold.
To obtain the local pressure drop factor in the multiphase stage, k m f , a linear regression study was performed using the multiphase experiments in operating conditions shown in Table 2 and another 13 experiments carried out in the experimental test bench, where WF did not reach the stage of the overheated vapor. This set of 13 experiments are not further discussed in this paper since the WF did not complete evaporation, but the results were used to feed data to estimate k m f .
Performing a study of the linear regression with the experimental results, with a similar procedure as conducted in the single-phase flow, the local pressure drop factor in the multiphase stage k m f j shown in Equation (13) needs to be corrected with the expression displayed in Equation (14), where χ j is the average vapor quality in the mixing chamber, and Re v is the Reynolds number with the properties of saturated vapor.
In the single-phase stage, k m f is set to 0, while k s f is set to 0 during the multiphase flow. The correlation of Equation (14), in the experiments analyzed, retrieved a value of k m f j , a maximum value in the range of about 30 to 50 in each experimental point.

Simulation Results and Discussion
The overall pressure drop results of the experiments are shown in Table 2 and the physical model described in the previous section implemented in CFD simulations in ANSYS Fluent are shown in Table 6. Figure 7 graphically displays, in a dispersion graph, the 9 results of Table 6 (color red) and the 13 experimental results where evaporation was not completed, as previously mentioned. The average mean error of all 9 experimental conditions analyzed is 15.47%, with 8 of 9 multiphase experimental points inside the 30% error band, and 4 of those experiments inside the 15% error band. The experimental point outside the 30% error band is barely outside it. When the 13 data points where evaporation was not completed are included in the analysis, the average mean error is 14.57%, with 10 experimental points outside the 15% error band and 4 experiments outside the 30% error band. The methodology and correlations previously described allow to obtain an estimation of the contributions to the total pressure drop in the different features of the boiler, resulting in a distribution of pressure drop as shown in Table 7. This table suggests that the greatest impact on the total pressure drop of the boiler is mainly due to the effect of the mixing chambers on the multiphase stage, in spite of WF remaining in the multiphase stage for much shorter length than the single-phase stage (see Table 8 and Figure 8).
The distribution of the length of the path of the WF that undergoes single-phase and multiphase stages, detailed in Table 8, can be graphically depicted in the graphs in Figure  8. Those graphs represent the average pressure drop in a normalized position in each path of the WF, where position "0" stands for the WF inlet and position "1" stands for the outlet of the WF. Note that when the graphs present a high increment on the value of pressure drop gradient ( / ), at around the position 0.8-0.85 in every experiment displayed, this means that the WF has begun the phase transition until it is quite close to the outlet of the boiler, where / is then low and stable again for a short fraction of the boiler, which represents the stage of overheated vapor. Shortly before the last single-phase stage, the WF undergoes a sharp decrease in / , which stands for the dryout stage. Only just before exiting the boiler, the / slightly increases in the stage of overheated vapor.  Table 6, while the blue points correspond to experiments where evaporation was not completed.  Table 6, while the blue points correspond to experiments where evaporation was not completed. The methodology and correlations previously described allow to obtain an estimation of the contributions to the total pressure drop in the different features of the boiler, resulting in a distribution of pressure drop as shown in Table 7. This table suggests that the greatest impact on the total pressure drop of the boiler is mainly due to the effect of the mixing chambers on the multiphase stage, in spite of WF remaining in the multiphase stage for much shorter length than the single-phase stage (see Table 8 and Figure 8). The distribution of the length of the path of the WF that undergoes single-phase and multiphase stages, detailed in Table 8, can be graphically depicted in the graphs in Figure 8. Those graphs represent the average pressure drop in a normalized position in each path of the WF, where position "0" stands for the WF inlet and position "1" stands for the outlet of the WF. Note that when the graphs present a high increment on the value of pressure drop gradient (dp/dx), at around the position 0.8-0.85 in every experiment displayed, this means that the WF has begun the phase transition until it is quite close to the outlet of the boiler, where dp/dx is then low and stable again for a short fraction of the boiler, which represents the stage of overheated vapor. Shortly before the last single-phase stage, the WF undergoes a sharp decrease in dp/dx, which stands for the dryout stage. Only just before exiting the boiler, the dp/dx slightly increases in the stage of overheated vapor.

Conclusions and Future Work
This paper expands the capabilities of the work previously presented in 2019 by the same authors of the present paper, in which a methodology to calculate heat transfer in multiphase flow was developed. The contribution in the present study consists of a methodology to estimate the pressure drop in multiphase flow by means of CFD simulations, and implement it in a boiler of an ORC in a WHRS. This was accomplished by combining CFD simulations, bibliographic correlations for pressure drop, numerical analysis and experimental results, to obtain a correlation for the local pressure drop factor to the specific Note that the pressure drop in the multiphase stage is non-linear and does not follow the intuitive pattern of the Darcy-Weisbach correlation (Equation (4)), where higher mass flow rate means higher pressure drop. In general, when the multiphase flow is involved, this will not be necessarily true. As shown in Figure 8, the multiphase stage has a remarkable impact on the pressure drop of the boiler. The highest peaks of the pressure drop gradient occur in the multiphase flow and the multiphase local pressure drop factor is significantly higher than the local pressure drop factor in the single-phase stage, given by Table 6. Therefore, comparing the pressure drop of similar experiments of the boiler, it seems that the highest pressure drop is found in the experiments where the WF remains in the multiphase transition for a longer portion of the WF path, not necessarily the experiment of higher mass flow. As a consequence, and for a lesser pressure drop purpose, counter flow involves a better performance for this specific boiler. The highest temperature gradient between the fluids involved in the boiler (WF and air, in this case) is produced near the multiphase stage of the WF, meaning a higher heat exchange between fluid and, therefore, the shortest presence on the multiphase stage of the WF.

Conclusions and Future Work
This paper expands the capabilities of the work previously presented in 2019 by the same authors of the present paper, in which a methodology to calculate heat transfer in multiphase flow was developed. The contribution in the present study consists of a methodology to estimate the pressure drop in multiphase flow by means of CFD simulations, and implement it in a boiler of an ORC in a WHRS. This was accomplished by combining CFD simulations, bibliographic correlations for pressure drop, numerical analysis and experimental results, to obtain a correlation for the local pressure drop factor to the specific geometry of the boiler studied and create a tool to be implemented in the CFD software ANSYS Fluent. Besides overall pressure drop, the methodology developed also identifies the contribution of different features of the WF flow and concluded that the main contributor on the total pressure drop of the boiler is the local pressure drop on the multiphase stage. The proposed methodology also allows to predict the mass flow distribution in a bundle of parallel tubes in the single-phase stage and the multiphase stage.
Regarding the results in the pressure drop, the simulations showed an acceptable agreement with experimental results obtained in the operating conditions of the boiler, with a mean average error of 15.47%, with 8 of 9 multiphase experimental points inside the 30% error band, and 4 of those experiments inside the 15% error band.
Future work will mainly target the experimental procedure. The work will focus on expanding the database of experimental results and adjusting the proposed correlation if needed. Additionally, the methodology should be tested with different geometries with different challenges in local pressure drop assessments.