Acta Polytechnica CTU Proceedings doi:10.14311/AP.2016.4.0073 Acta Polytechnica CTU Proceedings 4:73–79, 2016 © Czech Technical University in Prague, 2016 available online at http://ojs.cvut.cz/ojs/index.php/app CFD SIMULATION OF UPWARD SUBCOOLED BOILING FLOW OF FREON R12 Tomas Romsy∗, Pavel Zacha Faculty of Mechanical Engineering, Czech Technical University in Prague, Technická 4, Prague 6, 166 07, Czech Republic ∗ corresponding author: tomas.romsy@fs.cvut.com Abstract. Subcooled flow boiling under forced convection occurs in many industrial applications to maximize heat removal from the heat source by the very high heat transfer coefficient. This work deals with CFD simulations of the subcooled flow boiling of refrigerant R12 solved by code ANSYS FLUENT r16. The main objective of this paper is verification of used numerical settings on relevant experiments performed at DEBORA test facility. Also comparisons with previously provided simulation on NRI Rez are presented. Data outputs from this work are basis to subsequent calculations of steam-water mixture cooling of Pb-Li eutectic. Keywords: CFD simulation, subcooled boiling, freon R12. 1. Introduction Subcooled flow boiling can solve a wide range of high power thermal challenges in the following areas: Cen- tral Processing Unit (CPU) and computational ap- plication Specific Integrated Circuits (ASICs) [1, 2], ultra-high brightness Light Emitting Diodes (LEDs) and lasers sources [3], automotive power electronics [4], or avionics [5], where research findings have made pos- sible the development of a new generation of cooling hardware, which promises order of magnitude increase in heat dissipation compared to today’s cutting edge cooling schemes. Subcooled flow boiling is expected to realize the high heat flux cooling for electronic devices [6, 7], or for high heat flux cooling in mi- crogravity [8]. Subcooled flow boiling also plays a significant role in case of cooling of plasma facing component of thermonuclear reactors. The subcooled flow boiling is a high complex form of two-phase liquid flow in which the single-phase flow merges into the two-phase flow and back to the single- phase flow. A local appearance and disappearance of the two-phase flow allow the use of the latent heat of vaporization for a heat removal from the surface into the liquid without causing the boiling crisis. In comparison with a single-phase fluid it represents the efficient way to improve heat transfer from the wall to the flowing liquid. Therefore, it provides an ability of the sufficient heat transfer at high heat fluxes. Detailed investigation of this phenomenon can help for its better understanding to make wider improvement. One of the perspective application of subcooled flow boiling is so-called cold trap facility. This device serves for corrosion impurity separation from the metal eutec- tic (Pb, Pb-Li), whereas the separation is carried out by the cooling of the eutectic down to the required temperature. This facility considered in Centre of Advanced Nuclear Technology project (CANUT) [10] provides sufficient heat removal from eutectic across the wall to the subcooled steam-water mixture (cool- ing loop) and subsequently to the tertiary water loop providing ultimate heat removal. Part of the cold trap design is cooling concept which strongly depends on steam-water mixture behaviour. Description of cooling loop behaviour is provided with the aid of Computational Fluid Dynamics (CFD) code ANSYS FLUENT. Verification of suitable numerical settings and computational methodology should be determined in the first step. Selected experiments on the DEB- ORA facility were used for this purpose. A detailed description of this device can be found in [11] and [9]. DEBORA project at CEA Grenoble (French Alter- native Energies and Atomic Energy Commission) was focused on determination of subcooled boiling charac- ter by using refrigerant dichlorodifluoromethane (R12). Provided experiment outputs contain especially void fraction, gas velocity, temperature of the fluid and the bubble size. Because both freon R12 and water have relative physical similarity, we can use this freon for his lower operating parameters in comparison with water. This leads to a much simpler and safer mea- surements. The geometry of the DEBORA facility is shown in Figure 1. It is the vertical heated tube with inner diameter of 19.2 mm. Freon is heated along the 3.5 m long part of pipe. The rest of 5 m long tube is unheated. It corresponds to 1 m of inlet section and 0.5 m of outlet section. A radial profile of the vapour volume fraction and its speed at the end of the heated part were measured by optical probe. Ad- ditionally, the bubble size profiles were available for this area. Axial wall temperatures were measured using thermocouples, but radial temperatures were not available [9]. CFD model was performed based on the DEBORA facility parameters and the boundary conditions and standard numerical approaches were set according to the experiment conditions. The results from simula- 73 http://dx.doi.org/10.14311/AP.2016.4.0073 http://ojs.cvut.cz/ojs/index.php/app Tomas Romsy, Pavel Zacha Acta Polytechnica CTU Proceedings Figure 1. a) Sketch of the DEBORA test geometry [9] b) Computational model sketch with detail of used grid. Measured quantity p [MPa] v [m/s] qw [W/m2] Tin [◦C] Tsat [◦C] Case 1 3.008 0.86 58 260 67.89 94.24 Case 2 3.006 0.891 58 260 73.7 94.21 Table 1. Specified conditions of the chosen experimental cases. tions provided at Nuclear Research Institute Rez (NRI Rez) [12] are included in Section 4 for better compar- ison. Main purpose of this comparison are different approaches in calculation which lead to opportunity to compare them. While the previous simulation from NRI Rez was performed at ANSYS FLUENT r13, presented work was processed at ANSYS FLUENT r16. Earlier release 13 missed necessary subcooling numerical models and these calculations were imple- mented by individual User Defined Functions (UDF). However, starting from release 15, numerical models describing this phenomenon are implemented inter- nally. The main output of this work is to determine selected parameters (pressure, velocity, temperatures) in the specified geometry position and compare them against the measured experimental data provided from DEBORA test facility and against the previously cal- culated values from NRI Rez simulations. 2. Model description Simulation model is based on geometry of the DEB- ORA experimental device. The 10° V-cutout, shown in Figure 1, was considered due to the symmetry and to reduce the computational time. Slightly rough hexahedral grid was tested in the first steps of the cal- culation. Further enhancements were applied after ob- taining better knowledge of the calculation behaviour on this grid. Also much finer grid was tested. This grid was only the 10° V-cutout of tested geometry for computational time reduction due to a greater number of cells closer to the wall. However, there was a too small angle of the tip in the axis (centre of the tube) and therefore a smooth auxiliary wall had to be applied with tip neglecting. As the result, flow cross-section has been little changed, but hexahedral cells could be used for the whole geometry. The simu- lation results for this grid were similar to case with 45° V-cutout geometry. The 45° V-cutout grid has been selected for the final evaluation of the results, be- cause no geometry adjustments are used there. Also, it is not appropriate to use fine meshes with respect to numerical approach for the used boiling computa- tional model, like accretion of bubbles near the wall etc. This final grid contains 90 000 hexahedral cells with their higher concentration in the wall vicinity, shown in Figure 1. Cells quality was evaluated for Equi-size skew, which does not exceed value of 0.5 and by the Aspect-ratio with the maximum of 10.3. Two cases with different boundary conditions based on the DEBORA experiments are considered. The first one – Case 1 – has a much larger subcooled inlet state then the Case 2, but the saturation conditions and heat flux are almost the same for both of them. The main specified conditions of these experiments are given by Table 1 and setting of the boundary conditions in the simulation model corresponds with them. Usage of these two experimental cases brings the opportunity to try if the simulation model can solve various ratio of boiling states; from subcooled up to almost boiling crisis – Departure from Nucle- ate Boiling (DNB). The physical values setting for the freon R12 liquid and vapour is based on these conditions. Therefore linear functions depending on the temperature were used for density, specific heat capacity, viscosity and thermal conductivity. These quantities were obtained from NIST database [13]. 74 vol. 4/2016 CFD Simulation of Upward Subcooled Boiling Flow of Freon R12 Setting of the model boundary conditions: • Inlet: Velocity inlet – the same speed for both phases Turbulence intensity – estimated at 3 % • Outlet: Pressure outlet – 0 backflow volume fraction • Heated wall: Heat flux The measuring line was created on one symmetry wall directly at the interface between the heated and outlet non-heated section. This line contains 11 radi- ally distributed evaluation points which were created to monitor calculated data. They are situated in the centre of each appropriate cell and their distribution is shown in the Figure 2. radial distance in [m] . ........ . . Figure 2. Distribution of the evaluation points. 3. Solver setting The solver setting is mainly supported by the docu- mentation [12, 14–17]. Time dependent solution has to be considered with time step set to 0.01 s for each calculated case. The Realizable k−� turbulence model with the Non-Equilibrium wall function has been cho- sen. This turbulence model is solved for each phase, but the wall function is calculated only for single- phase. The Multiphase coupled model was used for the sequential solution of velocity and pressure fields and the Least square Cell-based method was set for gradients. All discretization schemes were set to 1st order, otherwise the simulation is very susceptible to the formation and disintegration of bubbles along the simulation domain. Calculations were carried out on the sixteen Intel Xeon E5-2660 processors. The complete final calcu- lation took about 8 hours of computer time for each case. Setting of boiling in ANSYS FLUENT 16: Multiphase Eulerian – RPI Boiling Model This boiling model is appropriate to use for subcooled flow, where condition Tsat − Tbulk > 3 K is fulfilled. Model does not calculate the vapour temperature, but it is fixed at the saturation temperature instead. As an alternative, Non-Equilibrium Boiling Model can be used with criterion Tsat − Tbulk ≤ 3 K, where the vapour temperature is included in the solution process. This condition could occur at the higher model levels, where the liquid can be already sufficiently heated close to saturation conditions. However, the RPI Boiling Model is considered for the next steps of this work. Figure 3 shows scheme of the selected sub- models for RPI Boiling Model. Figure 3. Scheme of the selected sub-models for RPI boiling model. The Wall Lubrication Model (effect of virtual mass) was neglected, because of big impact to the solution grid fineness, which is not very desirable for the model. 4. Calculation results The results of performed calculations are compared with data calculated by NRI Rez (ANSYS FLUENT r13) and measured values from the DEBORA experi- ment given by [12]. The radial profiles of important monitored parameters for both calculated cases are described below. These profiles are obtained from measuring line, for contained evaluation points respec- tively. Dispositional (orthographic) views are also given for contours of selected parameters in the com- putational area, see Figure 12. If they had not been applied, the side view would have showed just the long and thin tube with nothing obvious to see in the 75 Tomas Romsy, Pavel Zacha Acta Polytechnica CTU Proceedings computational domain. Given the similarity of the boiling behaviour of both cases, only results for Case 1 are listed. 4.1. Radial temperature profiles Unfortunately due to unavailable experimental data for radial temperature profiles, only the comparison with the results provided by NRI Rez is mentioned. For both calculated cases the evaluated radial profile has a similar course, as shown in Figure 4 for Case 1 and in Figure 5 for Case 2. . . . . . . . . . . . . . . . . . . . . . . Figure 4. Radial temperature profile, Case 1. . . . . . . . . . . . . . . Figure 5. Radial temperature profile, Case 2. Only slight differences between the temperature values obtained by NRI Rez and calculation up to 0.78 ◦C in the fluid midstream and up to 1.36 ◦C near the wall for Case 1 are seen. As regards the Case 2, here is obvious difference between temperature profiles in the midstream. The NRI Rez results for Case 2 are very close to saturation temperature almost across the whole radial profile, except the wall. The character of CTU calculation temperature profile has a decreasing character up to the axis – tube centre (little bit subcooled state). But near the axis the temperature difference does not exceed 0.58 ◦C. 4.2. Radial void fraction profile Radial profiles for void fraction with measured data and also NRI Rez results have been compared already. This is shown in Figure 6 (Case 1) and Figure 7 (Case 2). In Case 1 is seen that all the data are in good agreement near the wall, but the curves have slightly different courses closer to axis – tube centre. In case of this work, such disagreement is probably caused by collecting of non-condensed secondary phase in the midstream and efforts of liquid go closer toward the wall. This can be seen from void fraction rendered along the calculated area in Section 4.5, where the results evaluation along the model is described. Figure 6. Radial void fraction profile, Case 1. Figure 7. Radial void fraction profile, Case 2. 76 vol. 4/2016 CFD Simulation of Upward Subcooled Boiling Flow of Freon R12 The greater disagreement with other data is obvious for Case 2. The CTU simulation model estimates al- most abrupt increase of void fraction near the heated wall pretty well. But in the midstream the void frac- tion profile is lower than the others. This probably correlates with the results for temperature quite well. For the treatment of such disagreements it would be possible to try a different approach in model setting, or using of Non-Equilibrium Boiling Model for a future simulations respectively. 4.3. Radial vapour velocity profile Vapour velocity radial profile comparison for Case 1 is illustrated in Figure 8. It can be seen that the CTU and NRI Rez results have similar character, but compared with the measurement they have a higher drop of vapour velocity in the fluid midstream. This is probably caused by using of non-ideal single phase wall function. This will be the subject of further investigation in a future simulations. Figure 8. Radial vapour velocity profile, Case 1. In Case 2, see Figure 9, it is well visible that curve of CTU calculation follows the measured data near the wall pretty well. On the other hand, the same problem as in the previous case occurs in the fluid midstream. Figure 9. Radial vapour velocity profile, Case 2. 4.4. Radial bubble diameter profile Very important parameter in terms of subcooled boil- ing behaviour is diameter of vapour bubbles. For experiment measured data; the Sauter mean diameter is defined as the diameter of a sphere, which has the same ratio of volume and surface area as the base particle (bubble) and the Average diameter Dg is the diameter of the bubble, which would have the equiva- lent two-phase flow with the same density of bubbles amount, the same density of interfacial area and the same diameter for all bubbles, see [12]. For Case 1, the results are depicted in Figure 10. Here is evident that the bubble diameters along the radius have similar size and character as the measured data. It means that the bubbles are breaking up near the wall and their diameter is increasing in the centre of the stream. Slighter difference in the fluid midstream probably corresponds with the above listed diagram for void fraction. Also bubble diameter along the tube for this case is shown in Figure 12. Figure 10. Radial bubble diameter profile, Case 1. Figure 11. Radial bubble diameter profile, Case 2. The results for Case 2, shown in Figure 11, are in very good agreement with Experiment – Dg and NRI Rez near the wall. Whereas the bubble diameters in the midstream have little bit lower values compared with other data, they are still in very good agreement 77 Tomas Romsy, Pavel Zacha Acta Polytechnica CTU Proceedings Figure 12. Axial liquid temperature, void fraction and bubble diameter, Case 1. with Experiment – Dg, because it is seen, that the character of radial profile have nearly the same course. Smaller values for bubbles are probably caused by lower temperature of fluid in the midstream, shown in Figure 5. Therefore, the question of using the Non-Equilibrium Boiling Model arises here. 4.5. Axial results profiles Figure 12 (Case 1) due to similarity of the both cases of boiling behaviour mentioned above is presented for a better understanding, how the fluid behaves along the tube and where boiling occurs. There are orthographic views for the fluid temperature, void fraction and bubble diameters. There is visible, that all of these quantities correlate with each other. For example; the fluid temperature approximately in the middle of the model reaches the saturation condition near the wall and the void fraction begins to occur there. The higher level the fluid stream in tube reaches, the higher temperature in the midstream occurs. It also means, that the higher concentration of void fraction is there. On the profile of bubble diameters, it is shown how to bubbles reach the greater diameters in the higher levels of the tube and how they are trying to get into the centre of fluid stream. Considering this view, it should be observed, that for setting of bubbles modelled by Yao-Morel Model, the smallest possible value of the bubbles have to be set. For both cases this value was set to −105 m. 5. Conclusion The main goal of this work was to try the possibility to utilize the newly implemented boiling models in CFD solver ANSYS FLUENT r16. The DEBORA experi- ments were chosen for comparison with data received from these simulations. Some of the earlier works, such as NRI Rez (FLUENT r13), were also tried to simulate boiling experiments by using older version of Fluent CFD solver. However, suitable boiling models had to be supplied by User defined functions (UDF) in this earlier version. For comparison, how the solver with newly implemented boiling models (FLUENT r16) stand up against the older version with UDF’s, the evaluation of these data compared with CTU is described. Two cases (Case 1 and Case 2) of DEBORA ex- periments were chosen. Only inlet conditions for fluid temperature are different in these two cases. It means that the first one is much more subcooled case then the second one. The comparison between measured data, calculations from NRI Rez and results of this work are in the high level of agreement. Only few small differences occur. The vapour velocity disagree- ment of both calculations can be primarily caused 78 vol. 4/2016 CFD Simulation of Upward Subcooled Boiling Flow of Freon R12 by numerical approach, or by using of non-ideal sin- gle phase wall function respectively. Also, for radial void fraction closer to axis the curve have slightly different course compared to experiment and NRI Rez. This disagreement is probably caused by collecting of non-condensed secondary phase in the midstream and effort of liquid go closer toward the wall. And also in Case 2 the void fraction has the lower values in the fluid midstream. Therefore the question of using the Non-Equilibrium Boiling Model arises here. It will be subject of further research, because it can treat the fact, that for Case 2 the midstream fluid temperature could reache the values close to saturation conditions in the measured section. After evaluation of all the results we can say, that the CFD solver with newly implemented boiling mod- els (FLUENT r16) can be successfully used for prob- lems where the subcooled boiling occurs. But some minor deficiencies are still present. Unfortunately due to very complicated and sensitive numerical model used for description of this boiling phenomenon, there are many degrees of freedom at the model settings. Also variables that are tracked interact with each other and individual numerical approaches have limited va- lidity. This leads to time-consuming study to solve these minor deficiencies and verifying all of the used hypotheses (and their interaction) by experiments. All the knowledge and experiences described above are currently being used for simulations of more com- plicated situation of subcooled boiling of water, which occurs for example in Pb-Li cold trap cooling system, and they are continuously being improved. List of symbols p Pressure [MPa] v Inlet liquid velocity [m/s] qw Heat flux [W/m2] T Temperature [°C] Tin Inlet liquid temperature [°C] Tbulk Bulk liquid temperature [°C] Tsat Saturation temperature [°C] % Liquid density [kg/m3] cp Specific heat capacity [kJ/kg K] ν Dynamic viscosity [Pa s] λ Thermal conductivity [W/m K] References [1] S. G. Kandlikar. High flux heat removal with microchannels – A roadmap of challenges and opportunities. Heat Transfer Engineering 26(8):5–14, 2005. doi:10.1080/01457630591003655. [2] N. Khan, K. C. Toh, D. Pinjala. Boiling heat transfer enhancement using micro-machined porous channels for electronics cooling. Heat Transfer Engineering 29(4):366–374, 2008. doi:10.1080/01457630701825481. [3] A. Bar-Cohen, P. Wang, E. Rahim. Thermal management of high heat flux nanoelectronic chips. Microgravity Science and Technology 19(3):48–52, 2007. doi:10.1007/BF02915748. [4] F. Ramstorfer, H. Steiner, G. Brenn, et al. Subcooled boiling flow heat transfer from plain and enhanced surfaces in automotive applications. ASME J Heat Transfer 130(1), 2008. doi:10.1115/1.2780178. [5] I. Mudawar. Assessment of high-heat-flux thermal management schemes. IEEE Transactions on Components and Packaging Technologies 24(2):122–141, 2001. doi:10.1109/6144.926375. [6] J. L. Parker, M. S. El-Genk. Enhanced saturation and subcooled boiling of FC-72 dielectric liquid. International Journal of Heat and Mass Transfer 48(18):3736–3752, 2005. doi:10.1016/j.ijheatmasstransfer.2005.03.011. [7] Y. Madhour, J. Olivier, E. Costa-Patry, et al. Flow boiling of R134a in a multi-microchannel heat sink with hotspot heaters for energy-efficient microelectronic CPU cooling applications. IEEE Transactions on Components, Packaging and Manufacturing Technology 1(6):873–883, 2011. doi:10.1109/TCPMT.2011.2123895. [8] K. Suzuki, H. Kawamura. Microgravity experiments on boiling and applications: Research activity of advanced high heat flux cooling technology for electronic devices in Japan. Annals of the New York Academy of Sciences 1027(1):182–195, 2004. doi:10.1196/annals.1324.017. [9] E. Krepper, R. Rzehak. CFD for subcooled flow boiling: Simulation of DEBORA experiments. Nuclear Engineering and Design 241(9):3851–3866, 2011. doi:10.1016/j.nucengdes.2011.07.003. [10] O. Frybort. Background report – Pb-Li17 cold trap. Nuclear Research Centre Rez, 2014. [11] J. Garnier, E. Manon, G. Cubizolles. Local measurements on flow boiling of refrigerant 12 in a vertical tube. Multiphase Science and Technology 13(1&2), 2001. doi:10.1615/MultScienTechn.v13.i1-2.10. [12] L. Vyskočil. Simulation of the subcooled boiling in Fluent 13. Nuclear Research Institute Rez, 2011. [13] E.W. Lemmon, M.O. McLinden, D.G. Friend. Thermophysical Properties of Fluid Systems. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69, P. Linstrom, W. Mallard (Eds.). National Institute of Standards and Technology. http://webbook.nist.gov. [14] T. Romsy, P. Zacha. CFD simulation of subcooled boiling, 2016. ERIN, The 8th International Conference for Young Researchers and PhD Students, April 23–25, 2014, Blansko, Czech Republic. [15] ANSYS, Inc. Fluent User’s Guide, 2013. Release 15.0. [16] ANSYS, Inc. Fluent Theory Guide, 2013. Release 15.0. [17] ANSYS, Inc. FLUENT Multiphase 15.0 Optional Lecture 03 – Wall Boiling Models: Training materials, 2014. 79 http://dx.doi.org/10.1080/01457630591003655 http://dx.doi.org/10.1080/01457630701825481 http://dx.doi.org/10.1007/BF02915748 http://dx.doi.org/10.1115/1.2780178 http://dx.doi.org/10.1109/6144.926375 http://dx.doi.org/10.1016/j.ijheatmasstransfer.2005.03.011 http://dx.doi.org/10.1109/TCPMT.2011.2123895 http://dx.doi.org/10.1196/annals.1324.017 http://dx.doi.org/10.1016/j.nucengdes.2011.07.003 http://dx.doi.org/10.1615/MultScienTechn.v13.i1-2.10 http://webbook.nist.gov Acta Polytechnica CTU Proceedings 4:73–79, 2016 1 Introduction 2 Model description 3 Solver setting 4 Calculation results 4.1 Radial temperature profiles 4.2 Radial void fraction profile 4.3 Radial vapour velocity profile 4.4 Radial bubble diameter profile 4.5 Axial results profiles 5 Conclusion List of symbols References