CET Volume 86 DOI: 10.3303/CET2186241 Paper Received: 1 November 2020; Revised: 11 March 2021; Accepted: 15 April 2021 Please cite this article as: Maarawi A., Anxionnaz-Minvielle Z., Coste P., Di Miceli Raimondi N., Cabassud M., 2021, 1d/3d Numerical Approach for Modelling a Milli-structured Heat Exchanger/reactor, Chemical Engineering Transactions, 86, 1441-1446 DOI:10.3303/CET2186241 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor Antoinette Maarawia*, Zoé Anxionnaz-Minviellea, Pierre Costea, Nathalie Di Miceli Raimondib, Michel Cabassudb a Univ. Grenoble Alpes, CEA, LITEN, DTBH, Laboratoire Echangeurs et Réacteurs, F-38000 Grenoble b Laboratoire de Génie Chimique, Université de Toulouse, CNRS, INPT, UPS, Toulouse, France antoinette.maarawi@cea.fr Milli-structured heat exchangers/reactors are one of the intensified technologies, which have been successfully used in the case of exo- or endo-thermal chemical syntheses. Enhanced performances of such technologies compared to conventional apparatuses have been observed. This work is focused on the numerical investigation of the “DeanHex” heat exchanger/reactor with a 2D-meandering square cross-section channel with a hydraulic diameter of 2 mm. Throughout the years, only a part of the corrugated channel had been modelled due to the high ressource consumption of CFD simulations. Therefore, in order to model the reactor in a less resource demanding way than a full 3D approach, a new numerical approach is developed. The 2×2 mm2 meandering flow channel, which is the core of the heat exchanger/reactor design, is depicted by a 1D homogeneous model coupled to a 3D model for the surrounding solid in which the channel is embedded. This offers a compromise between the simulation at small scale and the simulation at the reactor scale, giving a great advantage in computational efficiency over meshing and computing 3D channels. In this paper, the 1D/3D approach is validated with heat transfer calculations between the 1D monophasic flow milli-channel and its 3D surrounding using an experimental work previously conducted on the “DeanHex” heat exchanger/reactor. 1. Introduction Process intensification (PI) has been defined in the literature (Stankiewicz and Moulijn, 2000) as a multidisciplinary engineering approach with the aim of developing innovative technologies, such as compact heat exchangers, that offer more energy-efficient and environmental-friendly processes, with higher productivity and safer operation. Compact heat exchangers with milli- and micro- channels have experienced some developments, like the implementation of corrugated geometries, in the purpose of offering new possibilities of heat transfer enhancement (Kapustenko P.O. et al., 2015). Tan Wei Q. et al. (2018) investigated the effects of corrugated plate geometries on the improvement of the thermo-hydraulic performances of a plate heat exchanger. Moreover, these compact heat exchangers can turn into heat exchangers/reactors when the reacting fluids are introduced at the inlet on one side. The removal or addition of heat for the reaction is thus enhanced because of the intensified heat transfer performances of such technologies. The plate HEX reactor (DeanHex) built by the Laboratoire de Génie Chimique (LGC, France) and the Commissariat à l’Energie Atomique et aux Energies Alternatives (CEA, France) is a fine example of such intensified devices. It consists of a succession of process plates (where the reaction takes place) and utility plates (for cooling or heating), where millimetric channels are engraved by laser machining. In this work, a square process channel of a hydraulic diameter dh equal to 2 mm with periodic zigzag units is considered. The geometry of the channel is represented in Figure 1. The radius of curvature of the bends Rc, the length of the straight section Ls and the angle between two straight sections θ are equal to 1.5 mm, 7 mm and 90 °, respectively. The interest of this particular corrugated geometry lies in the generation of secondary flow (Dean vortices) which helps improving mixing and heat transfer. In fact, the transition between laminar and turbulent flow in these DeanHex zigzag channels is around a Reynolds 1441 number equal to 200 (Théron et al., 2014). This allows reaching a quasi-plug flow and heat and mass transfer intensification while sufficient residence time is achieved to complete chemical conversion. Figure 1: Schematic of the zigzag channel with geometric parameters. However, from a numerical point of view, the complexity of the geometry, along with the high channel length to hydraulic diameter ratio, make the modelling of such devices with a full 3D approach difficult and very time- consuming regarding the high mesh refinement required in such geometries notably with transitional flow. The full-scale DeanHex reactor has a channel of total developped length of about 6.6 m. In the previous years, only a part of the channel could have been modelled using CFD approaches. Moreau et al. (2017) built a 3D model for the DeanHex channels, using COMSOL Multiphysics 4.3b, allowing the estimation of axial dispersion coefficients from calculated velocity and concentration fields in laminar flow. The simulated channel had a total length of about 0.27 m. Tetrahedral and cubic meshes were used in order to compromise between the simulation accuracy and the computing time. Shi et al. (2019) characterised the flow and heat transfer in the DeanHex square zigzag channel in laminar flow conditions by conducting 3D simulations using the finite volume CFD code, ANSYS CFX 16. Wavy channels with different geometries (Ls=2-12 mm) and a total length of 0.1 m were modelled. The total number of elements needed to mesh the wavy channels of Ls=3.5 mm and Ls=12 mm was about 4.3 million and 7.6 million, respectively. This shows that the 3D modelling of a portion of the wavy channel is highly resource demanding. Therefore, in order to overcome the limitations of CFD simulations of such corrugated channels, this work proposes a new 1D/3D numerical approach that makes the modelling of the whole reactor possible. The full description of this approach is detailed in this paper as well as its validation by a previous experimental work conducted in the DeanHex millimetric zigzag channel. 2. 1D/3D numerical approach description A full 3D numerical approach allows the direct modelling of local transfer mechanisms in the HEX reactor. However, given the complex geometry of the zigzag channel and the high L/dh ratio, the 3D simulations are very costly in resources, especially in computing time. Therefore, in order to model the HEX reactor in a less resource demanding way, a new 1D/3D approach has been developed, using COMSOL Multiphysics V5.5. A 1D model is built for the meandering flow channel, using COMSOL’s Pipe Flow Module. The latter solves for homogeneous cross sections, which means it performs averaged calculations in the channel’s cross-sections. The modelled variables vary only along the length of the channel. Correlations, which were obtained during previous experimental works, are integrated in this 1D model. Consequently, this approach does not only help to avoid meshing the cross section of the channel with a full 3D mesh, but also offers a simulation method that is capable of predicting the performance of the reactor (reaction conversion, selectivity, etc.). The 3D channel is represented by a 1D curve "C" following its axial direction, embedded in a 3D surrounding solid "S". However, as seen from the solid matrix, the boundary condition between the channel and “S” is imposed by the channel on its axis “C” which is a line without any cross-section. In reality, this boundary condition is imposed on the walls of the channel with a 2×2 mm² cross-section, not on the axis. Therefore, in the present work, the cross- section of the millimetric zigzag channel, being non-negligible compared to the whole HEX reactor, must be taken into account in the model in order to correctly represent the DeanHex reactor. Consequently, the proposed solution is to define, around the curve “C”, the three-dimensional channel "S2" as shown in Figure 2. Then a thermal chain “C-S2-S” is established to transmit the boundary condition from "C" to the walls of "S2". A sensitivity study was carried out in order to determine the material properties of “S2” (density ρ (kg‧ m-3), specific heat capacity Cp (J‧ kg-1‧ K-1) and thermal conductivity λ (W‧ m-1‧ K-1)). The variations of the density and heat capacity of "S2" have showed no influence on the simulated thermal profiles. However, the effect of the radial thermal conductivity is clearly observed. As expected, by increasing λ to about 1000 W‧m-1‧K-1, the propagation of the temperature boundary condition from the channel “C” to the wall of “S2” is improved. Above this value, no more effect is observed on the thermal profiles. Consequently, to stay on the safe side, a radial thermal conductivity of 105 W‧m-1‧K-1 is applied to “S2” to transmit the “thermal information” from "C" to the walls of "S2", regardless of the system conditions. 1442 Figure 2: Illustration of the simulated geometry using the 1D/3D approach. Since the geometry of the studied millimetric channel is composed of periodic zigzag units, imposing a radial infinite conductivity to "S2” requires the definition of a new curvilinear coordinate system (cc) for “S2”. It allows the generation of a new base of vectors (ccx, ccy, ccz) at any point of “S2”, as shown in Figure 3. The thermal conductivity of "S2" is thus defined according to this curvilinear coordinate system: infinite (equal to 10 5 W‧m- 1‧K-1) along ccy and ccz (green and blue vectors respectively) and zero along ccx (red vector). Figure 3: The new curvilinear coordinate system applied to the 3D zigzag channel. 3. Model equations For incompressible fluid at steady state, the 1D pipe flow module calculates the pressure, velocity and temperature by solving the following momentum (Eq(1)), continuity (Eq(2)) and energy balance (Eq(3)) equations inside the channel referred to as the “process channel ” (subscript p): ρ𝑝𝑢𝑝‧∇𝑢𝑝 = −∇𝑝𝑝 − 𝑓𝐷 𝜌𝑝 2𝑑ℎ,𝑝 𝑢𝑝|𝑢𝑝| (1) 𝛻‧(𝐴𝑝𝜌𝑝𝑢𝑝) = 0 (2) 𝜌𝑝 𝐴𝑝 𝐶𝑝,𝑝𝑢𝑝‧𝛻𝑇𝑝 = 𝛻‧𝐴𝑝𝜆𝑝𝛻𝑇𝑝 + 𝑓𝐷 𝜌𝑝𝐴𝑝 2𝑑ℎ,𝑝 |𝑢𝑝| 3 + 𝑄𝑤𝑎𝑙𝑙 (3) 𝑄𝑤𝑎𝑙𝑙 = ℎ𝑝Z𝑝(𝑇𝑆 − 𝑇𝑝) (4) 𝑢𝑝 is the cross-section averaged velocity (m∙s-1), 𝑝𝑝 the pressure (Pa), 𝑇𝑝 the process fluid temperature (K), 𝑓𝐷 the Darcy friction factor (dimensionless). The process fluid (water) physico-chemical properties such as ρ𝑝 the density (kg‧m-3), 𝐶𝑝,𝑝 the heat capacity at constant pressure (J‧kg-1‧K-1) and 𝜆𝑝 the thermal conductivity (W‧m- 1‧K-1) are assumed variable in function of 𝑇𝑝. 𝐴𝑝 is the theoretical pipe cross-sectional area (m2) and 𝑑ℎ,𝑝 the mean hydraulic diameter (m) of the process channel. 𝑓𝐷 is correlated to Reynolds number 𝑅𝑒𝑝 for the millimetric zigzag channel as shown in Eq(5) (Théron et al., 2014). The Reynlods number in the channel is calcuated using Eq(6), where μ is the fluid viscosity (Pa‧s). 𝑓𝐷 = 24.3𝑅𝑒𝑝 −0.71 𝑓𝑜𝑟 20 < 𝑅𝑒𝑝 < 200 𝑓𝐷 = 6𝑅𝑒𝑝 −0.43 𝑓𝑜𝑟 𝑅𝑒𝑝 > 200 (5) 𝑅𝑒 = 𝜌‧𝑢‧𝑑ℎ 𝜇 (6) 1443 Furthermore, 𝑄𝑤𝑎𝑙𝑙 is the heat exchanged (W∙m-1) across the channel’s wall (Eq(4)), ℎ𝑝 the convective heat transfer coefficient inside the channel (W∙m-2∙K-1), Z𝑝 the wetted perimeter of the channel (m) and 𝑇𝑆 the temperature (K) of the 3D surrounding solid. In this work, ℎ𝑝 is calculated from a Nusselt number correlation Nu𝑝 (dimensionless) established for the DeanHex millimetric zigzag channel (Anxionnaz et al., 2011): Where 𝑅𝑒𝑝 and 𝑃𝑟𝑝 are the Reynolds and Prandtl dimensionless numbers in the process channel, respectively. The Prandtl number is calculated as follows: 𝑃𝑟 = 𝐶𝑝‧𝜇 𝜆 (8) In addition, the resolved heat balance equation for the 3D surrounding solid “S” is given by: 𝛻‧(−𝜆𝑆𝛻𝑇𝑆) = 𝛻‧𝑞𝑢 (9) Where 𝜆𝑠 is the thermal conductivity (W‧m-1‧K-1) of “S” and 𝑞𝑢 the heat flux boundary condition (W∙m-2) imposed on the walls of “S”. The latter wil be further described in the following, where the subscript u refers to the “utility channel ”. 4. Experiments vs 1D/3D simulations 4.1 Description The validation of the proposed 1D/3D approach is discussed by comparing the simulated results to experimental ones. The experimental work of Théron et al. (2014), in which thermo-hydraulic performances of the DeanHex millimetric zigzag channel had been characterised, is used. During this study, the HEX reactor consisted of a succession of three process plates sandwiched between four utility plates. 2 mm square cross-section channels, with periodic zigzag units, were engraved on both plates by laser machining. The process and utility streams are separated by a 2 mm metal thickness (closing plates). In the present work, only the first process plate is modelled (Figure 4a). The computational domains consist of a three-dimensional zigzag channel of total length equal to 2.2 m “S2” with a radially infinite thermal conductivity, centered by a one-dimensional zigzag curve “C”. “S2” is surrounded, from both sides, by a metallic solid matrix “S” of 2 mm thickness, representing the closing plates. The thermal chain “C-S2-S” is thus established as shown closely in Figure 4b. Nonetheless, the utility channels do not take part in the computational domains but are replaced by heat flux boundary conditions 𝑞𝑢 at the upper and lower walls of “S”: 𝑞𝑢 = ℎ𝑢(𝑇𝑆 − 𝑇𝑢) (10) Where 𝑇𝑢 is the utility fluid temperature (K). The convective heat transfer coefficient ℎ𝑢 (W∙m-2∙K-1) in the zigzag utility channels is deduced from the Nusselt correlation (Eq(11)) proposed by Shi (2019): 𝑁𝑢𝑢 = ℎ𝑢𝑑ℎ,𝑢 𝜆𝑢 = 0.2(𝑅𝑒𝑢 0.67 + 8.9)𝑃𝑟𝑢 0.3 ( 𝑑ℎ,𝑢 𝐿𝑠,𝑢 ) 0.4 (11) Where 𝑅𝑒𝑢 and 𝑃𝑟𝑢 are the Reynolds and Prandtl dimensionless numbers for the utility side, repectively. 𝑑ℎ,𝑢 and 𝐿𝑠,𝑢 are, correspondingly, the hydraulic diameter (m) and the straight length (m) between the bends in the utility zigzag channel. 𝜆𝑢 is the thermal conductivity of the utility fluid (water) estimated at 𝑇𝑢 as for all the other physicochemical properties of the utility fluid. A correction factor 𝐴𝑢 𝐴𝑆 is applied to 𝑞𝑢 to consider the difference of heat transfer areas between the utility channels (𝐴𝑢, m2) and the walls of “S” (𝐴𝑠, m2). 𝐴𝑢 is the heat transfer area of the utility channels which can be calculated using Eq(12): 𝐴𝑢 = 𝑃𝑢. 𝐿𝑢. 𝑁𝑢 (12) Where 𝑃𝑢 and 𝐿𝑢 are correspondingly the perimeter (m) and the total developed length (m) of the utility channel. 𝑁𝑢 is the number of utility channels engraved on one utility plate. Note that, the utility plate that is below the first process plate of the pilot, exchanges heat with the first and second process plates. Therefore, 1 2 𝑃𝑢 is used when imposing the heat flux boundary condition on the bottom side of “S”. A free tetrahedral meshing is applied to the 3D domains where only thermal conduction occurs. Consequently, there is no need for a refined finely constructed meshing. Furthermore, the experimental conditions in Théron et al. (2014) are applied to the model and are presented in Table 1. Nu𝑝 = ℎ𝑝𝑑ℎ,𝑝 𝜆𝑝 = 0.16𝑅𝑒𝑝 0.66𝑃𝑟𝑝 0.33 (7) 1444 Table 1: Operating conditions applied for the simulations Utility Side Process Side Qu (kg‧h-1) Tu (°C) Qp (kg‧h-1) Tp (°C) ≈152 ≈15.6 2.4 ≈76 4 5.5 6.9 Figure 4: Illustrations from the 1D/3D model of: (a) the process plate with the inlet and outlet of the flow (in red: the position of thermocouples in the experiments) and (b) a zoom-in on the process inlet where the C-S2-S chain is showed. 4.2 Results Since for the experiments the heat transfer process is not limited by the utility side, the comparison of the 1D/3D simulation results with experimental results will reveal any error that could be made on the process side and thus the proposed approach will be appropriately evaluated. In Figure 5, the temperature evolution is represented in function of the axial position of the process channel, for different “process” flowrates. In the experiments, 5 K-type thermocouples were located along the first process channel, thus a part of the thermal profile is observed (Théron et al., 2014). As shown in Figure 5, the simulated process thermal profiles are in good agreement with the experimental measurements. One aberrant point is observed and corresponds to the temperature measured at 0.28 m from the inlet. This discrepancy, observed for all the curves, can be explained by a misplacement of the thermocouple during the experiments. All the simulated temperatures at this point are higher than that measured by this thermocouple. In fact, the fluid flowing through the process channel is normally hotter than the plate material. Therefore, the probe maybe in contact with the material of the plate instead of being inserted in the flowing fluid. Regardless this point, the maximum relative error is less than 8 %, which is assumed acceptable. Consequently, the 1D/3D approach with the radially infinite thermal conductivity imposed for "S2" is validated considering these experimental data. Figure 5: Evolution of experimental and simulated process fluid temperature along the process channel for different process flowrates. 1445 5. Conclusions A new 1D/3D numerical approach is proposed in this work in order to model the DeanHex zigzag channel in a less resource demanding way than CFD approaches. The latters allow the simulation of only few zigzag periodic units with a high CPU cost. However, avoiding the 3D meshing of the channel, by using the proposed 1D/3D approach, reduces significantly the computational time to about few minutes and allows the modelling of the full- scale heat exchanger/reactor at a reasonable CPU cost. By comparing the 1D/3D simulated temperature results with experimental ones for one process plate of a 2.2 m long zigzag channel, a maximum relative error less than 8 % is observed. Thus, the 1D/3D approach seems promising and shows a great potential for the purpose of predicting the performance of the whole reactor. This leads to the the next step which consists of completing the heat transfer model by a mass transfer model and to validate it with experimental results. The present work is a first step towards addressing the scale-up issue of “DeanHex” reactors. Expecting a future industrialisation of such equipment, a generalised scaled-up model will be used to predict the performance of the reactor in terms of conversion rate and energy efficiency at industrial scale. References Anxionnaz, Z., Theron, F., Tochon, P., Couturier, R., Bucci, P., Vidotto, F., Gourdon, C., Cabassud, M., Lomel, S., Bergin, G., Bouti, S., Peerhossaini, H., Lemenand, T., Habchi, C., 2011. RAPIC project: Toward competitive heat-exchanger/reactors. The 3rd European Process Intensification Conference (EPIC). Kapustenko P.O., Kukulka D.J., Arsenyeva O.P., 2015. Intensification of heat transfer processes. Chemical Engineering Transactions 45, 1729–1734. https://doi.org/10.3303/CET1545289 Moreau, M., Di Miceli Raimondi, N., Le Sauze, N., Gourdon, C., Cabassud, M., 2017. A new numerical method for axial dispersion characterization in microreactors. Chemical Engineering Science 168, 178–188. https://doi.org/10.1016/j.ces.2017.04.040 Shi, H., 2019. Etude du transfert de chaleur en canaux millimétriques de type zigzag pour le développement et l’extrapolation de réacteurs-échangeurs compacts. Université Toulouse 3 - Paul Sabatier. Shi, H., Di Miceli Raimondi, N., Fletcher, D.F., Cabassud, M., Gourdon, C., 2019. Numerical study of heat transfer in square millimetric zigzag channels in the laminar flow regime. Chemical Engineering and Processing - Process Intensification 144, 107624. https://doi.org/10.1016/j.cep.2019.107624 Stankiewicz, A., Moulijn, J., 2000. Process Intensification: Transforming Chemical Engineering. Chemical Engineering Progress 22–34. Tan Wei Q., Sun Li, Jia Zhigang, Li Jie, 2018. Thermal-hydraulic performance analysis of corrugated plate heat exchanger. Chemical Engineering Transactions 70, 1153–1158. https://doi.org/10.3303/CET1870193 Théron, F., Anxionnaz-Minvielle, Z., Cabassud, M., Gourdon, C., Tochon, P., 2014. Characterization of the performances of an innovative heat-exchanger/reactor. Chemical Engineering and Processing: Process Intensification 82, 30–41. https://doi.org/10.1016/j.cep.2014.04.005 1446