# NUCLEAR THERMAL HYDRAULIC AND TWO-PHASE FLOW

EDITED BY : Jun Wang, Kaiyi Shi, Zhaoming Meng and Shripad T. Revankar PUBLISHED IN : Frontiers in Energy Research

#### Frontiers Copyright Statement

© Copyright 2007-2018 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-612-3 DOI 10.3389/978-2-88945-612-3

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# NUCLEAR THERMAL HYDRAULIC AND TWO-PHASE FLOW

Topic Editors:

Jun Wang, University of Wisconsin-Madison, United States Kaiyi Shi, Liupanshui Normal University, China Zhaoming Meng, Harbin Engineering University, China Shripad T. Revankar, Purdue University, United States

Scene of nuclear power plant. Image: Liang Dunhuang.

Nuclear energy is one of the most important clear energy and contributes more than 10% electric power to human society in the past decades of years. The nuclear thermal hydraulic and two-phase flow is one of the basic branches of nuclear technology and provides structure design and safety analysis to the nuclear power reactors. In the new century, the basic theoretical research of thermal hydraulic and two-phase flow, and innovative design for the next generation nuclear power plants (especially for the small modular reactor and molten salt reactor), along with other nuclear branches, constantly support the development of nuclear technology.

Citation: Wang, J., Shi, K., Meng, Z., Revankar, S. T., eds (2018). Nuclear Thermal Hydraulic and Two-Phase Flow. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-612-3

# Table of Contents

*05 Editorial: Nuclear Thermal Hydraulic and Two-Phase Flow* Jun Wang, Kaiyi Shi, Zhaoming Meng and Shripad T. Revankar

### CHAPTER I

### GENERAL THERMAL HYDRAULIC


Bin Sun, Cheng Peng, Di Yang and Hongwei Li

# CHAPTER II

#### COMPUTATION FLUID DYNAMICS SIMULATION


# CHAPTER III

# CORE THERMAL HYDRAULIC


# CHAPTER IV

#### SEVERE ACCIDENT ANALYSIS

*93 Review on Core Degradation and Material Migration Research in Light-Water Reactors*

Jun Wang, Chen Wang, Kaiyi Shi and Guanghui Su


# Editorial: Nuclear Thermal Hydraulic and Two-Phase Flow

#### Jun Wang<sup>1</sup> \*, Kaiyi Shi <sup>2</sup> , Zhaoming Meng<sup>3</sup> and Shripad T. Revankar <sup>4</sup>

*<sup>1</sup> Nuclear Engineering and Engineering Physics Institute, University of Wisconsin-Madison, Madison, WI, United States, <sup>2</sup> School of Chemistry and Materials Engineering, Liupanshui Normal University, Liupanshui, China, <sup>3</sup> College of Nuclear Science and Technology, Harbin Engineering University, Harbin, China, <sup>4</sup> School of Nuclear Engineering, Purdue University, West Lafayette, IN, United States*

#### Keywords: nuclear, thermal hydraulic, two-phase flow, experiment, computation fluid dynamics, core, severe accident

#### **Editorial on the Research Topic**

#### **Nuclear Thermal Hydraulic and Two-Phase Flow**

Nuclear Thermal hydraulic and two-phase flow research is always a typical component in nuclear engineering research from the 1960s (D'Auria, 2017). In American Nuclear Society (ANS) annual meetings, thermal hydraulic and two-phase flow is always one of the largest branches. International meeting on nuclear reactor thermal hydraulic (NURETH) starts at the 1980s, provides the communication platform for all the scientists in this region. As the first research topic in Frontiers in Energy Research, Nuclear Section, this topic is with great expectation.

Along with the development of computer technology, more modern simulation tools are involved, together with experiments, make contributions to the thermal hydraulic and two-phase flow research. This technology is more about the highly geometry-dependent transient phenomena, provide structure design and safety optimization of nuclear power plants. The existing research regions mainly include: general thermal hydraulic, thermal hydraulic experiment, two-phase flow modeling, boiling and condensation, safety analysis of current reactors and Gen IV reactors, computation fluid dynamics simulation, system code development and validation, critical heat flux (CHF), core thermal hydraulic and subchannel analysis, severe accident and so on.

Experiment is still the most important part of basic thermal hydraulic. Mechanism study of steam condensation with air in an imaged reactor component vertical tube experimental facility is conducted by Wang et al.. Some basic thermal hydraulic is studied by self-developed code. A onedimension gas turbine system transient code is developed for compressible flow simulation based on hybrid semi-implicit method by another Wang et al.

Thermal hydraulic modeling can now be studied by computation fluid dynamics simulation. Subcooled boiling, which may appear in the vertical pipe of a reactor system, is studied by Zhang et al.. The uncertainties from the software Fluent is also considered in Zhang's work. Another simulation work to study the flow distribution in reactor system heat exchangers is done by Zhou et al., who aims to reduce this situation to avoid system risk. Luchao studies the single-phase nanomaterial flow of CuO-water by looking into the hydromechanics and heat transfer mechanics (She and Fan). In the reactor, grids will affect the coolant flow, and bring in turbulence which may affect heat transfer. Dong et al. builds a model to study this situation by predicting the critical heat flux.

Thermal hydraulic is always coupled with neutronics, to provide highly accurate simulation of the core performance. For example, a theoretical Ultra-high temperature fuel is designed by Song et al. based on neutronics and thermal hydraulic analysis, which could provide capacity for the nuclear electric propulsion system. Core thermal hydraulic improves the primary use efficiency of nuclear power, and also gives the core safety performance analysis. Sun et al. works on this region, and also provide a core thermal hydraulic design to for a micro nuclear system.

#### Edited and reviewed by:

*Hyun Gook Kang, Rensselaer Polytechnic Institute, United States*

#### \*Correspondence:

*Jun Wang jwang564@wisc.edu*

#### Specialty section:

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

Received: *04 July 2018* Accepted: *23 July 2018* Published: *10 August 2018*

#### Citation:

*Wang J, Shi K, Meng Z and Revankar ST (2018) Editorial: Nuclear Thermal Hydraulic and Two-Phase Flow. Front. Energy Res. 6:80. doi: 10.3389/fenrg.2018.00080*

**5**

Severe accident turns to be one of the popular topics after Fukushima Daiichi nuclear disaster. Wang et al. made a review work on severe accident history and core degradation mechanics during his time at Xi'an Jiaotong University supervised by Prof. Guanghui Su. Yin et al. uses MELCOR, which is developed by Sandia National Laboratory, conducting a whole small modular reactor severe accident analysis to study its response during a station blackout. Containment is the final safety protective screen in nuclear power plants. Wen tries to find a way to help keeping containment integrity by testing surface tension of a kind of coolant material (Wen et al.).

Future thermal hydraulic research will have various difference. Multi-physical coupled code including neutronic, fuel, mechanics, and thermal hydraulic development is taking more important position. Meanwhile, more core thermal hydraulic experimental data are still required, especially for the critical heat flux, and accident tolerant fuels. Computation fluid dynamics technology is still developing, and it should be better for two-phase flow simulation. Severe accident and accident tolerant fuels thermal hydraulic are still important parts.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

### REFERENCE

D'Auria, F. (2017). "A historical perspective of nuclear thermal-hydraulics," in Thermal-Hydraulics of Water Cooled Nuclear Reactors, ed F. D'Auria (Duxford: Elsevier Ltd), 41–87.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wang, Shi, Meng and Revankar. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Experimental Study on the Condensation of Steam With Air Out of the Vertical Tube Bundles

Lu Wang<sup>1</sup> , Ping Chen<sup>1</sup> , Yi Zhou<sup>1</sup> , Wei Li <sup>1</sup> , Changbing Tang<sup>1</sup> , Yifei Miao<sup>1</sup> and Zhaoming Meng<sup>2</sup> \*

*<sup>1</sup> Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, China, <sup>2</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Heilongjiang, China*

In this paper, experimental study was carried out in an open natural circulation loop, to figure out the condensation heat transfer characteristic on the outside of the vertical pipe tube bundle and single- tube. Condensation heat transfer coefficients have been obtained under the wall subcooling ranging from 5 to 22◦C, total pressure ranging from 0.3 to 0.6 MPa and air mass fraction ranging from 0.05 to 0.65. The influence of mixed gas pressure, condensate depression of walls, and content of non-condensable gas on the condensation heat transfer performance is analyzed. Under the condition of the same air mass fraction, the condensation heat transfer coefficient increases with the increase of pressure. Also, as the air mass fraction is more than 30%, the effect of pressure will be weakened. The heat transfer features of pipe bundle and single tube is compared and studied, and an empirical correlation of the pipe bundle for the heattransfer coefficient is developed, covered all data points within 10%. According to the research results, heat transfer coefficient for pipe bundle decreases with the increase in non-condensable gas quality and wall subcooling, but increases with the increase in pressure.

Keywords: natural circulation, tube bundles, air-containing steam, condensation heat transfer, experimental investigation

# INTRODUCTION

In recent years, with the public's increasing concerns over the safety of the nuclear power plant,full deterministic evidence of evidence of nuclear power plants (NPP) safety is required (Wang et al., 2014), more natural circulation systems were applied in the designs of the advanced nuclear plants, because of the inherent passive safety. For example, Hualong Pressurized water reactor, Chinese Generation III nuclear reactor, uses an open natural circulation as the passive containment cooling system (PCCS) (Lakshmanan and Pandey, 2009). PCCS, which is one of the most important part in protecting the integrity of the containment under serve accidental conditions, can remove the heat in the containment by the condensation of steam over the specific heat transfer surface (Yu et al., 2008; Ye et al., 2012). The operating pressure of the open natural circulation is atmosphere pressure. It can make the vaporizing process easily occur and has a faster startup time.

However, the flow rate of the cooling water is constantly changing in the natural circulation process which is different from the forced circulation cooling water loop. Therefore, in order to understand the steam condensation phenomenon with non-condensing gas under the condition of natural circulation, and to provide the reliable empirical data and analysis models for the practical engineering applications, a more detailed research of the effect of the steam condensation is necessary.

#### Edited by:

*Jun Wang, University of Wisconsin-Madison, United States*

#### Reviewed by:

*Limin Liu, University of California, Berkeley, United States Gangyang Zheng, Independent Researcher, Beijing, China Kaiyi Shi, Liupanshui Normal University, China*

\*Correspondence:

*Zhaoming Meng mengzhaoming2011@hotmail.com*

#### Specialty section:

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

Received: *14 February 2018* Accepted: *04 April 2018* Published: *28 May 2018*

#### Citation:

*Wang L, Chen P, Zhou Y, Li W, Tang C, Miao Y and Meng Z (2018) Experimental Study on the Condensation of Steam With Air Out of the Vertical Tube Bundles. Front. Energy Res. 6:32. doi: 10.3389/fenrg.2018.00032*

**7**

Dehbi et al. (Dehbi, 1991) investigated the condensation phenomenon of steam with no-condensing gas in a vertical tube, and the calculation models of forced convection, natural convection and mixed convection were then established, which can be used to evaluate the effect of the molecular weight, the non-condensing gas content, the mixed gas pressure, the inlet flow rate of mixed gas et al. on the steam condensation process; Liu et al. (2000) conducted significant experimental research on the steam condensation process including the effect of non-condensing gas for the PCCS system of nuclear power plant, and evaluated the thermal elimination ability of the passive containment cooling system; Dehbi used the experimental model to calculate the total heat transfer for uncondensed gas vapor condensing outside the vertical tube and establishing forced convection, natural convection and mixed convection conditions. By using this model, the influence of noncondensable gas molecular weight, non-condensable gas content, mixed gas pressure and mixed gas inlet flow rate on the steam condensing can be judged. A large number of experimental studies have been conducted on the vapor condensation of gases by Liu, and the heat elimination capabilities of passive containment cooling systems have been evaluated (Wang, 2014; Su, 2016); Although the experiments of Liu and Dehbi took into account many factors, they had opposite results on the effect of undercooling on the heat transfer coefficient of the condensation.

Othmer (1929) carried out research by placing a single brass tube in stationary vapor atmosphere, and the empirical equations for the heat transfer coefficient correlated with the air volume and the supercooling degree of the tube surface had been further obtained. Ivashchenco (1989) carried out an experimental study on the condensation heat transfer process outside the vertical single tube with vapor containing nitrogen under the pressure of 0.5 MPa, the results showed that the heat transfer coefficient decreased rapidly with the increase of noncondensing gas content. Also, the empirical equations obtained by Uchida (Uchida et al., 1964; Tagami, 1965) associated the decrement of the condensation heat transfer coefficient with the percentage of air weight together. Zhuang et al. (2000) investigated the condensation heat transfer characteristics of the horizontal tube bundles and further analyzed the effect of cooling water flow rate as well as the impact of the mass fraction of non-condensing gas on condensation heat transfer of the tube bundles.

Despite of many contributions have been made, the previous research shows that the studies of the condensation heat transfer process mainly focus on the single tube and horizontal tube bundles. However, it is quite shortage for the experimental studies of condensation outside vertical tube bundles, especially with an open natural circulation.

In order to investigate the impact of various parameters on the external condensation heat transfer of tube bundles and obtain the empirical equations with the corresponding parameter ranges, experimental studies for the natural convection condensation process of vertical tube bundles and a single tube are carried out in this paper.

#### THE EXPERIMENT SYSTEM

The experiment system mainly consists of steam supply equipment, experimental vessel as well as its inner heat transfer tubes, water tank, experimental pipes and measuring instruments, as shown in **Figure 1.**

The internal diameter and the height of the experimental vessel are 416 and 3,550 mm, respectively. In the tube bundle experiment, the heating section consists of three parallel uniform stainless steel tubes with the inside and outside diameters of 34 and 38 mm separately, whose tube spacing and length are 45.76 and 2,234 mm, respectively. The inside and outside diameters of the heat transfer tube in the single tube experiment are identical with that of the heat transfer tubes in the tube bundle experiment.

The internal pressure of the experimental vessel is measured by a pressure sensor with a precision of 0.075%. Cooling water flux is measured by an electromagnetic flowmeter with a precision of 0.5%, and measuring range is−30∼30 m<sup>3</sup> /h. Steam flow measurement is achieved by the use of two vortex flowmeters with a precision of 0.75%, and measuring range are 0∼120 m<sup>3</sup> /h and 0∼40 m<sup>3</sup> /h, respectively. At the inlet and outlet of the experimental section, the temperatures of cooling water are measured by installing one nickel chromium-nickel silicon thermocouple.

The distribution of thermocouples in the experimental vessel is shown in **Figure 2**. **Figure 2A** shows thermocouples distribution in the single tube experimental vessel, and the thermocouples are vertically welded on the 9 cross sections of the experiment tube section (A∼I). Two thermocouples are symmetrically installed on the wall of each section. Meanwhile, one of the thermocouples is installed on each section to measure the temperature of the main mixture stream. In the experimental section of tube bundle, the thermocouples are vertically welded

on the 5 cross sections (A∼E), as shown in **Figure 2B**. Two thermocouples are installed on each section of the bracket which is arranged on the left and right pipes. For measuring the main temperature, two thermocouples are welded on both side of each tube in each section to measure wall temperature. A thermocouple is welded in the middle of the inlet and outlet header to measure the temperature of inlet and outlet header.

# EXPERIMENTAL DATA PROCESSING

The steam heat transfer rate through the tube is equal to the cooling water heat transfer rate, which is determined by the cooling water temperature change and the mass flow rate in the tube. The bulk temperature and the wall temperature as well as cooling water temperature are measured by themocouples. The bulk temperature and the wall temperature as well as cooling water temperature are measured by thermocouples. The enthalpy of cooling water is obtained from steam properties table. The flowrate of the cooling water is also measured by turbine flowmeter. Hence,the average heat transfer coefficient can be expressed as:

$$h = \frac{M\left(h\_{\rm gi} - h\_l\right)}{A\left(T\_b - T\_W\right)}\tag{1}$$

where M is the mass flow rate of the steam, kg/s; hgi and h<sup>l</sup> are the saturated steam enthalpy value and the condensate enthalpy value, kj/kg; A is the surface area of the heat transfer tube, m<sup>2</sup> ; T<sup>b</sup> and T<sup>w</sup> are standing for the mainstream temperature of the mixed gas and the outer surface temperature of the experimental tube, ◦C.

According to the measurement error of each instrument in the experiment,the total measurement error of the average heat transfer coefficient is <15%.

The non-condensable air quality can be calculated as:

Firstly, calculate the gas composition according to the Dalton partial pressure law:

$$\frac{N\_a}{N\_s} = \frac{P\_a}{P\_s} \tag{2}$$

Then, the non-condensable air quality:

$$\begin{split} \omega\_{a} &= \frac{m\_{a}}{m\_{a} + m\_{s}} = \frac{29n\_{a}}{29n\_{a} + 18n\_{s}} = \frac{29p\_{a}}{29P\_{a} + 18P\_{s}} \\ &= \frac{29\left(P - P\_{a}\right)}{29\left(P - P\_{a}\right) + 18P\_{s}} = \frac{29P - 29P\_{s}}{29P - 11P\_{s}} \end{split} \tag{3}$$

m is the mass, kg; n is the number of molecules, mol; M is the molecular weight; a is air; s is steam. The total pressure P is accurately collected by a pressure transmitter that condenses the experimental section.

When the saturated steam enters the condensation experimental section, the saturated steam becomes superheated steam due to the decrease of the partial pressure of the steam, and the cooling water will cool the superheated steam during the flow along the path so that the superheat loss is lost. Therefore, there is a certain degree of overheating in the process. However, since the condensation heat transfer coefficient is less affected by the superheat of the steam, it can be assumed that the steam is saturated. The partial pressure of steam can be obtained by checking the table of saturated steam properties, the formula is:

$$P\_s = f\left(T\_s\right) \tag{4}$$

Therefore, the non-condensable air quality is expressed as

$$\omega\_a = \frac{29P - 29f\,(T\_s)}{29P - 11f\,(T\_s)}\tag{5}$$

# ANALYSIS OF EXPERIMENTAL RESULTS The Effect of Supercooling Degree of the Tube Surface

In the condensation process, the pressure and the air mass fraction will have a great influence on the supercooling degree of the tube surface. The single impact of supercooling degree on the condensation heat transfer coefficient cannot be investigated precisely as the above two parameters vary. Therefore, we choose two start-up conditions for studying the effect of supercooling degree of the tube surface. One set of parameters are that the tube surface pressure P is 0.6 MPa and the air mass fraction W<sup>a</sup> is 75%. The other set of parameters are that the tube surface pressure P is 0.5 MPa and the air mass fraction W<sup>a</sup> is 75%.

As is shown in **Figure 3**, it can be known that the condensation heat transfer coefficient of the tube bundle decreases with the increase of the supercooling degree of the tube surface, and increases with the increase of pressure, which is similar to the trend of Dehbi correlation. The thickness of the condensation film will increase gradually with the increment of supercooling degree of the tube surface, which increases the heat and mass transfer resistance and leads to a decrease of the condensation heat transfer coefficient.

# The Effect of Mixed Gas Pressure and the Mass Fraction of Air

The experimental study for the steam condensation heat transfer process is carried out under the conditions of different pressure and air mass fraction. As is shown in **Figure 4**, under the condition of constant pressure, the condensation heat transfer coefficient decreases with the increase of air mass fraction, and the rate of descent decreases gradually. However, under the condition of the same air mass fraction, the condensation heat transfer coefficient increases with the increase of pressure, which is consistent with the results in **Figure 3**. Also, as the air mass fraction is more than 30%, the effect of pressure will be weakened.

This trend is due to the competition between gas pressure impact and air content impact. The partial pressure of the steam rises and saturation temperature increases accordingly along with the increment of mixture gas pressure, which strengthens the heat transfer. But the non-condensing gas layer will thickened due to the higher air content, which further leads the deterioration of heat transfer.

# A Comparison of Heat Transfer Characteristics of the Tube Bundle and the Single Tube

**Figure 5** shows the curve of the condensation heat transfer coefficient of the tube bundle and the single tube. The result indicate that the heat transfer coefficient of the tube bundle is higher than that of the single tube while the air content is low, but the heat transfer coefficient of the tube bundle will decrease rapidly with the air mass fraction grows, whereas the descending rate of the heat transfer coefficient of the single tube will decrease more slowly. After the air mass fraction reaches 40%, the heat transfer coefficient of the tube bundle will be almost consistent with that of the single. There are two main factors that affect the difference of heat transfer characteristics of the single tube and tube bundle: The steam flow rate of tube bundle is higher than that of the single tube in the experiment; The surface of the tube bundle in the process of condensation is covered by a convergence of non-condensing gas layer in the gap, which will lead to a heat transfer deterioration. **Figure 5** shows that the heat transfer coefficient of tube bundle is better than that of single tube, especially when the air content is low. This could be due to the steam flow rate of the tube bundle is greater than that of the single tube in the experiment, the higher steam flow rate near the tube surface may destroy the non-condensing gas layer and reduce the thermal resistance, thus the heat transfer coefficient is then strengthened. However, with the increase of non-condensing gas content, the convergence phenomenon of

FIGURE 6 | Comparison of correlation against experimental data.

the non-condensing gas layer in the tube bundle gap will be more obvious, which is equivalent to add a layer of condensation heat transfer thermal resistance on the surface of the vertical tube, so the condensation heat transfer coefficient then decreases, the heat exchange amount is reduced, and the steam flow rate of the decreases accordingly, which leads to a further deterioration of the heat transfer.

# The Derivation for Calculation Relations

As can be seen in **Figure 6**, we have compared the experimental values of condensation heat transfer coefficient with Liu' and Dehbi's correlations under the pressure of 0.4 Mpa and the air mass fraction of 40%. It can be seen that the experimental data fluctuates greatly and are very different from the calculation values of Liu and Dehbi, which is mainly due to the instability of the flow in the natural circulation system and the effect of the tube bundle. **Figure 6** shows that when the heat transfer coefficient experimental curve represents as a wave crest, the prediction curve of Liu will be a wave trough by coincidence, and the prediction curve of Dehbi will be a wave crest, which is basically the same for the Dehbi's with a defference of one-half cycle. This is due to both the proposed experimental correation and the Liu's correlation are established by the linear fitting course with the mean value of data sets, whereas Liu's data is relatively fewer, and the effect of the index of supercooling degree expressed by equations, is not taken into account. Hence, there is a deviation in the trend of the transient curve obtained.

In this paper, by combining with the influence factors of the condensation heat transfer coefficient of the tube bundle, i.e., the effect of the index of supercooling degree has been taken into account, and using the multiple linear fitting method and the experimental data, the experimental correlation among the condensation heat transfer coefficient of the tube bundle, the mixed gas pressure, supercooling degree of the tube surface, and the non-condensing gas content can be given as follows:

$$h = 1047 \times \frac{\mathcal{P}^{0.24}}{(T\_b - T\_w)^{2.18} w\_a^{0.87}} \tag{6}$$

The applicable ranges of the given equations are 0.05≤ w<sup>a</sup> ≤0.65,0.3≤ P ≤0.6 MPa and 5≤ T<sup>b</sup> − T<sup>w</sup> ≤22◦C. Comparison results between the predicted values of the experimental correlation and the experimental data are shown in **Figure 7.** The result indicates that the experimental correlation obtained by the linear fitting course are able to well predict the experimental condensation heat transfer coefficient, and experimental data almost all fall within the error band of ±10%.

As can be seen in **Figure 7**, we have compared the experimental values of condensation heat transfer coefficient with Liu' and Dehbi's correlations under the pressure of 0.5 MPa. The experimental values close to the Liu's correlation due to the

#### REFERENCES


proposed experimental correlation and the Liu's correlation are both established by the linear fitting course with the mean value of data sets.

#### CONCLUSIONS


#### AUTHOR CONTRIBUTIONS

LW: First author and the writer of the manuscript. ZM and YZ: Provided guidance for the manuscript. WL and PC: Provided help with writing this manuscript and the data analysis. CT and YM: Provided help with building the experimental bench and performing the experiment.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wang, Chen, Zhou, Li, Tang, Miao and Meng. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Hybrid Semi-implicit Method of 1D Transient Compressible Flow for Thermal-Hydraulic Analysis of (V)HTR Gas Turbine Systems

#### Jie Wang1,2, Xiaxin Cao<sup>1</sup> \*, Zhaoming Meng<sup>1</sup> , Jie Wang<sup>3</sup> and Ming Ding<sup>1</sup> \*

<sup>1</sup> College of Nuclear Science and Technology, Harbin Engineering University, Harbin, China, <sup>2</sup> Shanghai Nuclear Engineering Research and Design Institute, Shanghai, China, <sup>3</sup> Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing, China

#### Edited by:

Muhammad Zubair, University of Sharjah, United Arab Emirates

#### Reviewed by:

Ivo Kljenak, Jožef Stefan Institute (IJS), Slovenia Qazi Muhammad Nouman Amjad, University of Engineering and Technology, Pakistan

#### \*Correspondence:

Xiaxin Cao caoxiaxin@hrbeu.edu.cn Ming Ding dingming@hrbeu.edu.cn

#### Specialty section:

This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research

Received: 12 November 2017 Accepted: 30 April 2018 Published: 25 May 2018

#### Citation:

Wang J, Cao X, Meng Z, Wang J and Ding M (2018) A Hybrid Semi-implicit Method of 1D Transient Compressible Flow for Thermal-Hydraulic Analysis of (V)HTR Gas Turbine Systems. Front. Energy Res. 6:44. doi: 10.3389/fenrg.2018.00044 Transient thermal-hydraulic analysis of (very) high temperature gas-cooled reactors gas turbine systems (HTRGTSs) needs system transient analysis codes. However, compared with the mature system transient thermal-hydraulic codes of pressurized water reactors (PWRs), the system analysis codes of HTRGTSs have not been so fully developed. In this paper, a new hybrid semi-implicit (HSI) method is proposed based on the semi-implicit method and nearly-implicit method. In the HIS method, a new calculation strategy is devised: the convective term is treated explicitly to solve pressure and velocity, while density and temperature are solved in an implicit manner to get a convergent, stable, and accurate solution in multiple transient scenarios. The HSI method was further validated via the shock-tube benchmark problem and verified via FLUENT simulations. In FLUENT simulations, outlet pressure transient, inlet mass flow transient and inlet temperature transient were studied. It was found that the HSI method is capable of capturing both the fast and slow compressible flow transients with good convergence and stability. Furthermore, an adaptive time step scheme is proposed for faster calculations, considering the maximum relative density difference and Courant–Friedrichs–Lewy condition.

Keywords: hybrid semi-implicit method, system thermal-hydraulic analysis, one-dimensional compressible transient solver, high temperature reactor gas turbine system, one-dimensional

# INTRODUCTION

High temperature gas-cooled reactors gas turbine system (HTRGTS) is inherently safe and highly efficient with a reactor outlet temperature of 700∼1,000◦C. HTRGTS uses a Brayton cycle, which includes a reactor, one or two compressors, turbine, recuperator, precooler, and intercooler. The working medium is usually helium in the HTRGTS because of its chemical inertia. The helium in the Brayton cycle is compressed to very high pressure (e.g., 7 MPa) in compressors and furtherly heated through the reactor, and then expands to very low pressure (e.g., 2.7 MPa) to generate kinetic energy in the turbine.

There is a very big pressure difference in different parts of HTRGTS during normal operation, ranging from 2.7 MPa to 7.0 MPa, and moreover the pressure changes quickly during accidents. Since helium is compressible, a compressible transient system analysis code is very important and necessary to the accident analysis of HTRGTS, and it has been studied since 1970s. However, because the HTRGTS is not so widely and commercially used as the pressurized water reactors (PWRs), transient system analysis codes for HTRGTS are not so well developed as those of PWRs. Different from PWRs, HTRGTS uses helium as its working fluid, whose density strongly couples with pressure and temperature. Consequently, compressibility should be considered in the transient calculations of HTRGTS for more accurate solutions, especially when significant changes of pressure and temperature are involved.

To date, the developed system analysis codes related to HTRGTS include GTSim for MGR-GT and GTMHR-300 program in Japan, FLOWNEX for PBMR program in South Africa, and PLAYGAS and CATHARE2 in Europe. In GTSim(Yan, 1990), a one-dimensional compressible flow model was used to simulate the transient conditions. The numerical method used in GTSim is a non-iterative method that provides the strong coupling between energy and continuity equations, and the strong coupling between energy and momentum equations. However, in this method the coupling between continuity and momentum equations is not very strong. As a result, this method is not very good in solving pressure transients, where continuity and momentum equations are strongly coupled.

FLOWNEX (Greyvenstein, 2002; Rousseau et al., 2006; van Ravenswaay et al., 2006) uses an implicit pressure correction method, which is very similar to the- SIMPLE method (Patankar, 1980; Tao, 2001). This method strongly couples the continuity and momentum equations, but it does not couple the energy equation. The energy equation is solved independently after continuity and momentum equations are solved. Besides, as a system solver for one dimensional transient compressible flow, this method has two weaknesses: (1) it is an iterative method that may require many iteration steps for each time step, especially in solving the transient compressible flow, such as pressure transients; and (2) it requires the setting of a pressure relaxation factor, which is an empirical figure depending on transient conditions and user experiences, and thus greatly affects computational time. These two weaknesses limit its use in the system analysis codes of HTRGTS.

The CATHARE series codes (Widlund et al., 2005; Saez et al., 2006; Bentivoglio et al., 2008) use a fully-implicit time integration scheme for zero-dimensional and one-dimensional modules. This method represents a relatively larger computational effort, which might be compensated by the increased numerical stability and the possibility to choose greater time steps for slow and long transients (Bestion, 2008). However, in solving the transient compressible flow where there is a very strong coupling of pressure, temperature, and density, the fully implicit integration scheme requires too many iteration steps for a converged solution, especially in fast pressure transients.

Therefore, it is meaningful to devise a non-iterative method that requires less computational time in solving fast pressure transients, and allows greater time step in solving slow temperature transients. There is a potential of using the methods in other fully developed codes of PWRs, such as the non-iterative semi-implicit and nearly-implicit methods (Division, 2001), for the transient solvers of HTRGTS.

The structure of this paper is as follows. Section Analysis of Semi-Implicit and Nearly-Implicit Methods in Compressible Flow analyzes the performances of the semi-implicit and nearlyimplicit methods in solving compressible flow transients. Section A New Hybrid Semi-Implicit Method Proposes an Improved Non-Iterative Semi-Implicit Method, that is, Hybrid Semi-Implicit (HSI) method, based on the discussion of the semiimplicit and nearly-implicit methods. Section Verification of Hybrid Semi-Implicit presents the verification of the HSI method via shock tube benchmark problem and the comparisons of the HSI method with FLUENT simulations, in which an outlet pressure transient, an inlet mass flow transient and an inlet temperature transient are studied. Section Adaptive Time Step proposes an adaptive time step scheme for the HIS method for faster calculation speed. Section Conclusions presents the final conclusions.

### ANALYSIS OF SEMI-IMPLICIT AND NEARLY-IMPLICIT METHODS IN COMPRESSIBLE FLOW

# Solution Strategies of Semi-implicit and Nearly-Implicit Methods

One Dimensional Conservative Equations

In a pipe with a constant cross section area, the one-dimensional continuity, momentum and energy equations are, respectively,

$$\frac{\partial \rho}{\partial x} + \frac{\partial \rho u}{\partial x} = 0 \tag{1}$$

$$
\rho \frac{\partial u}{\partial \mathbf{r}} + \frac{1}{2} \rho c\_m \frac{\partial u^2}{\partial \mathbf{x}} = -\frac{\partial p}{\partial \mathbf{x}} - f\_\mathbf{w} \rho u \tag{2}
$$

$$c\_\nu \frac{\partial \rho \, T}{\partial \pi} + c\_\nu \frac{\partial \rho \, uT}{\partial x} = -\rho \frac{\partial u}{\partial x} + \frac{hD\_e}{A\_c} (T\_\nu - T) \tag{3}$$

In the discritization process of the semi-implicit and nearlyimplicit methods, the control volume is shown in **Figure 1**, in which [x<sup>i</sup> , xi+1] is control volume of density, pressure and temperature, and [xK, xL] is control volume of velocity. A first order upwind scheme is adopted for cell face density and temperature, and a central difference scheme is used for cell center velocity.

**Nomenclature:** ρ, u, T, P: density, velocity, temperature and pressure, respectively; 1τ , 1x: time step and space step, respectively; cv: specific heat capacity of helium under constant volume; fw: frictional factor; cm: averaged factor when 3D flow is integrated and averaged to 1D flow; De, A<sup>c</sup> : hydraulic diameter and cross Sec. area of the pipe; E,c, d, r, l, A, B, C, D, ap, aw, ae: coefficients used to simplify calculation; Superscript "∼": intermediate time variable; Superscript ".": the donored value at cell face; Superscript "<sup>n</sup> ": value of last time step; Superscript "n+1 ": value of this time step; subscripts K, L, M: cell center i-1/2, i+1/2 and i+3/2, respectively.

#### Semi-implicit Method

Discretize continuity and energy Equations (1) and (3), respectively, yields

$$
\tilde{\rho}\_{L}^{n+1} - \rho\_{L}^{n} + (\dot{\rho}\_{i+1}^{n} u\_{i+1}^{n+1} - \dot{\rho}\_{i}^{n} u\_{i}^{n+1}) \frac{\Delta \tau}{\Delta x} = 0 \tag{4}
$$

$$
\left(\tilde{\rho}\_{L}^{n+1} - \rho\_{L}^{n}\right) T\_{L}^{n} + \rho\_{L}^{n} \left(\tilde{T}\_{L}^{n+1} - T\_{L}^{n}\right) - \left[\frac{hD\_{c}}{c\_{\nu}A\_{c}} (T\_{w} - T)\right]\_{L}^{n}
$$

$$
\Delta \tau = -\left[\left(\dot{\rho}\_{i+1}^{n} \dot{T}\_{i+1}^{n} + \frac{P\_{L}^{n}}{c\_{\nu}}\right) u\_{i+1}^{n+1} - \left(\dot{\rho}\_{i}^{n} \dot{T}\_{i}^{n} + \frac{P\_{L}^{n}}{c\_{\nu}}\right) u\_{i}^{n+1}\right] \frac{\Delta \tau}{\Delta x} \tag{5}
$$

According to Taylor expansion, we get the expression of density ρ n+1 in term of P n+1 and T n+1 :

$$\tilde{\rho}\_{L}^{n+1} = \rho\_{L}^{n} + \left(\frac{\partial \rho}{\partial P}\right)\_{L}^{n} \left(P\_{L}^{n+1} - P\_{L}^{n}\right) + \left(\frac{\partial \rho}{\partial T}\right)\_{L}^{n} \left(\tilde{T}\_{L}^{n+1} - T\_{L}^{n}\right) \tag{6}$$

Substitute Equation (6) into energy Equation (5) and eliminate the density term ρ n+1 yields expression of the next-time step temperature T n+1 in term of P n+1 and u n+1 :

$$\begin{split} \tilde{T}\_{L}^{n+1} &= \left. T\_{L}^{n} - \left( \frac{\partial T}{\partial \rho} \right)\_{L}^{n} \left( \frac{\partial \rho}{\partial P} \right)\_{L}^{n} \left( P\_{L}^{n+1} - P\_{L}^{n} \right) \\ &- \left( \frac{\partial T}{\partial \rho} \right)\_{L}^{n} \left( \dot{\rho}\_{i+1}^{n} u\_{i+1}^{n+1} - \dot{\rho}\_{i}^{n} u\_{i}^{n+1} \right) \frac{\Delta \tau}{\Delta \mathbf{x}} \end{split} \tag{7}$$

Combining Equations (5)∼(7) to eliminate the terms of ρ n+1 and T n+1 yields expression of the next time step pressure P n+1 in term of u n+1 :

$$P\_L^{n+1} = P\_L^n + l\_L^n u\_i^{n+1} - r\_L^n u\_{i+1}^{n+1} - E\_L^n \Delta \tau \cdot \left(\frac{\partial \rho}{\partial T}\right)\_L^n / \left(\frac{\partial \rho}{\partial P}\right)\_L^n \tag{8}$$

Then, discretize the momentum equation in Equation (2)

$$\begin{aligned} \left[\rho\_i^n \left(u\_i^{n+1} - u\_i^n\right) + \frac{1}{2}\rho\_i^n c\_{m,i} \left[\left(u\_L^n\right)^2 - \left(u\_K^n\right)^2 - u\_L^n (u\_{i+1}^n - u\_i^n)\right.\right. \\ \left. + u\_K^n (u\_i^n - u\_{i-1}^n)\right] \frac{\Delta x}{\Delta x} &= -\left(P\_L^{n+1} - P\_K^{n+1}\right) \frac{\Delta x}{\Delta x} - f\_{w,i}^n \rho\_i^n u\_i^{n+1} \Delta x \end{aligned} \tag{9}$$

The expression of the next-time step velocity u n+1 in term of P n+1 can be derived:

$$u\_i^{n+1} = c\_i^n \left( P\_K^{n+1} - P\_L^{n+1} \right) + d\_i^n \tag{10}$$

It should be noted that the convection term in the momentum equation is treated explicitly in the semi-implicit method. The calculation strategy of the semi-implicit method is as follows:

Substitute Equation (10) into (8) to eliminate u n+1 yields the expression of P n+1 :

$$\begin{aligned} &-\boldsymbol{\varepsilon}\_{i}^{n}l\_{L}^{n}P\_{K}^{n+1} + (1 + \boldsymbol{\varepsilon}\_{i}^{n}l\_{L}^{n} + \boldsymbol{\varepsilon}\_{i+1}^{n}r\_{L}^{n})P\_{L}^{n+1} - \boldsymbol{\varepsilon}\_{i+1}^{n}r\_{L}^{n}P\_{M}^{n+1} \\ &= P\_{L}^{n} + l\_{L}^{n}d\_{i}^{n} - r\_{L}^{n}d\_{i+1}^{n} + \rho\_{L}^{n}RE\_{L}^{n}\Delta\tau \end{aligned} \tag{11}$$

Then substitute P n+1 into Equation (10) yields u n+1 .

Then substitute P n+1 and u n+1 into Equation (7) yields T n+1 .

Then substitute P n+1 and T n+1 into Equation (6) yields density ρ1.

Finally, substitute P n+1 and T n+1 into state equation yields another density ρ2, if the difference between ρ<sup>1</sup> and ρ<sup>2</sup> is small enough, then a converged solution is obtained.

#### Nearly-Implicit Method

In the nearly-implicit method, the continuity and energy equations are treated and the relationship between pressure and velocity is deduced. The convective term in the momentum equation is treated implicitly, as shown in Equation (12):

$$\rho\_i^n \left( u\_i^{n+1} - u\_i^n \right) + \frac{1}{2} \dot{\rho}\_i^n c\_{m,i} \frac{\Delta \tau}{\Delta x} \left[ \left( u\_L^n \right)^2 + 2 \mu\_L^n \right.$$

$$\left( u\_L^{n+1} - u\_L^n \right) - \left( u\_K^n \right)^2 - 2 \mu\_K^n \left( u\_K^{n+1} - u\_K^n \right) \Big]$$

$$- \frac{1}{2} c\_{m,i} \dot{\rho}\_i^n \left[ \mu\_L^n (u\_{i+1}^n - u\_i^n) - u\_K^n (u\_i^n - u\_{i-1}^n) \right]$$

$$\frac{\Delta \tau}{\Delta x} = - \left( P\_L^{n+1} - P\_K^{n+1} \right) \frac{\Delta \tau}{\Delta x} - f\_{w,i}^n \rho\_i^n u\_i^{n+1} \Delta \tau \tag{12}$$

The relationship between velocity u n+1 and pressure P n+1 can be deduced from Equation (12):

$$A\_i u\_{i-1}^{n+1} + B\_i u\_i^{n+1} + C\_i u\_{i+1}^{n+1} + D\_i = P\_K^{n+1} - P\_L^{n+1} \tag{13}$$

Density and temperature in the nearly-implicit method are solved implicitly in continuity Equation (14) and energy Equation (15), respectively.

$$
\tilde{\rho}\_{L}^{n+1} - \rho\_{L}^{n} + (\dot{\rho}\_{i+1}^{n+1} u\_{i+1}^{n+1} - \dot{\rho}\_{i}^{n+1} u\_{i}^{n+1}) \frac{\Delta \pi}{\Delta x} = 0 \tag{14}
$$

$$
\left[ \left( \tilde{\rho}\_{L}^{n+1} - \rho\_{L}^{n} \right) T\_{L}^{n} + \rho\_{L}^{n} \left( \tilde{T}\_{L}^{n+1} - T\_{L}^{n} \right) \right] \frac{\Delta x}{\Delta \pi} - \left[ \frac{hD\_{\varepsilon}}{c\_{\nu} A\_{\varepsilon}} (T\_{w} - T) \right]\_{L}^{n}
$$

$$
= \left( \dot{\rho}\_{i}^{n+1} \dot{T}\_{i}^{n+1} + \frac{P\_{L}^{n+1}}{c\_{\nu}} \right) u\_{i}^{n+1} - \left( \dot{\rho}\_{i+1}^{n+1} \dot{T}\_{i+1}^{n+1} + \frac{P\_{L}^{n+1}}{c\_{\nu}} \right) u\_{i+1}^{n+1} \tag{15}
$$

The calculation strategy of nearly-implicit method is in the following steps:


Frontiers in Energy Research | www.frontiersin.org May 2018 | Volume 6 | Article 44

#### TABLE 1 | Parameter settings.

#### Common setting


(a)Outlet pressure transient parameter settings for semi-implicit and nearly-implicit methods; (b)Outlet pressure transient parameter settings for HSI method and FLUENT simulation.

5) Finally, substitute P n+1 and T n+1 into state equation, another density ρ<sup>2</sup> can be solved. If the difference between ρ<sup>1</sup> and ρ<sup>2</sup> is small enough, a convergent solution is obtained.

In step (3) and step (4), density and temperature are solved in the following form:

$$a\_{\mathcal{P}^i} \phi\_L^{n+1} = a\_{\le i} \phi\_K^{n+1} + a\_{\le i} \phi\_M^{n+1} + b\_i \tag{16}$$

φ <sup>n</sup>+<sup>1</sup> mean density ρ n+1 in continuity equation and temperature T n+1 in energy equation.

### Performances of the Semi-implicit and Nearly-Implicit Methods

Two FORTRAN codes were developed based on the semiimplicit and nearly-implicit methods in order to analyze their performances in solving transient compressible flow. Pressure transient involves the strong coupling of pressure, density, velocity, and temperature, and it is used to testify the performances of the two methods. The transient condition is a step decrease of outlet pressure of a pipe from 2.7 MPa to 2.6 MPa. Other parameters, such as boundary conditions, time step, and space step, are shown in **Table 1**. The calculated pressure and temperature responses by the two methods are shown in **Figures 2**, **3**, respectively.

It can be seen from **Figure 2A** that the semi-implicit method is capable of capturing the fast pressure responses with little numerical diffusion. However, it has poor stability in capturing the slow temperature responses, as seen in **Figure 2B** that temperature diverged at 0.02 ∼ 0.04 s. It can be seen from **Figure 3B** that the nearly-implicit method has good stability in capturing the slow temperature responses. However, it cannot capture the fast pressure responses very accurately, as seen in **Figure 3A** that the pressure oscillations are suppressed and vanish very quickly.

One possible reason for the great numerical diffusion of the nearly-implicit method is that the time step 1τ = 10−<sup>5</sup> s is too large to capture the fast pressure responses accurately. The variation of time step from 10−<sup>6</sup> s to 10−<sup>8</sup> s was studied, and

some of the results are shown in **Figure 4**. It was found that if time step 1τ≤10−<sup>7</sup> s, the numerical diffusion of pressure can be very small. However, there are two weaknesses: (1) spurious oscillations exist at the pressure step change points; and (2) the time step is too small, which means a great computational effort.

As shown in Figures 2∼4, neither the semi-implicit method nor nearly-implicit method provides a stable and accurate solution in solving transient compressible flow.

#### A NEW HYBRID SEMI-IMPLICIT METHOD

# Further Discussion of the Semi-implicit and Nearly-Implicit Methods

For the semi-implicit method, velocity is implicit while density is explicit in continuity equation. In energy equation, only velocity is implicit, the other parameters, such as density, pressure, and temperature, are all explicit. In momentum equation, the pressure and frictional terms are implicit, while the convective

term is explicit. So the expressions of density, temperature, velocity, and pressure can be derived directly. In calculation, firstly, pressure, density and temperature are solved explicitly from the derived expressions, and then velocity is solved based on the solved pressure.

In the nearly-implicit method, the expression of pressure is derived from continuity and energy equations in the same way as in the semi-implicit method. Different from the semiimplicit method, the convective term in the momentum equation is treated implicitly. In calculation, velocity is solved first, then pressure is solved subsequently based on the solved velocity, and then density and temperature are solved implicitly by updating the solved pressure, velocity and density in continuity and energy equations.

There are mainly two differences between the semi-implicit and nearly-implicit methods: (1) the explicit or implicit treatment of the convective term in the momentum equation, which determines the different expressions of velocity, and (2) the explicit or implicit calculation of density and temperature. According to the two differences and the results in **Figures 2**–**4**, it is assumed that (1) in order to capture a detailed fast pressure transient response, the convective term in momentum equation should be treated explicitly, and (2) in order to get a stable solution for the slow temperature transients, coefficients in the continuity, and energy equations should be updated immediately using the previously solved parameters, such as pressure, velocity and density.

To testify this hypothesis, the convective term in momentum equation is treated implicitly as in the nearly-implicit method, and density and temperature are solved explicitly as in the semiimplicit method. If this hypothesis was correct, it is predicted that neither an accurate pressure response nor a stable temperature solution could be obtained. The results are shown in **Figure 5**, in which the time step is 1τ = 10−<sup>5</sup> s. It can be seen that for the fast pressure response, the artificial diffusion is too big, while for the slow temperature response, the solution diverged. The results are the same as predicted, so the hypothesis is reasonable.

# Proposal of a New Hybrid Semi-implicit Method

According to the hypothesis discussed above, a new HSI method is proposed. In the HSI method, the convective term in momentum equation is treated explicitly, and density and temperature are solved implicitly by updating coefficients in the continuity and energy equations using the previously solved parameters. The calculation strategy of the HSI method is in the following steps:


Under the same outlet pressure transient in section Analysis of Semi-Implicit and Nearly-Implicit Methods in Compressible Flow, and using parameter settings in **Table 1**, the pressure and temperature responses of the HSI method are shown in **Figure 6**. It can be seen from **Figures 6A,B** that the fast transient responses caused by pressure oscillations are captured and the slow temperature responses are stable and converged. As a result, the newly proposed HSI method has the advantages of both the semi-implicit and nearly-implicit methods but none of their weaknesses.

#### VERIFICATION OF HYBRID SEMI-IMPLICIT

The newly proposed HSI method has shown good performance in solving transient compressible flow. However, the accuracy of the results needs to be further verified. The HSI method was furtherly verified via the shock tube benchmark problem, and compared with FLUENT simulations.

#### Shock Tube Benchmark Problem

The shock tube problem is a typical Riemann problem, and has an exact solution (Toro, 2009). In the shock tube problem, a closed tube with a length of 100 m was separated into two equal parts by a membrane film, as shown in **Figure 7.** The left part of the tube is helium with a pressure of 2.0 MPa and a temperature of 400 K, and the right part of the tube is helium with a pressure of 1.0 MPa and a temperature of 400 K. If the film suddenly breaks in the center, the helium in two parts mixes together, and shock

wave, rarefaction wave, and discontinuity wave are produced in this process.

There is no exact closed-form solution to the Riemann problem for the Euler equations of shock-tube problem. However, it is possible to devise iterative schemes whereby the solution can be computed numerically to high accuracy. The space step is set as 0.05 m in both the iterative scheme and the HSI method. Comparison of the exact solution with that of the HSI method is shown in **Figure 8**. Three time steps of 10−<sup>4</sup> s, 10−<sup>5</sup> s, and 10−<sup>6</sup> s are calculated using the HSI method. It was found that when 1τ = 10−<sup>4</sup> s the solution is divergent, and when 1τ = 10−<sup>5</sup> s and 1τ = 10−<sup>6</sup> s, the results of the HSI method are almost

FIGURE 8 | Comparisons between the HSI method and the exact solution.

identical, as seen in **Figure 8**. So the time step is set as 10−<sup>5</sup> s in this case to get a convergent solution with shorter computational time in the following discussion.

It can be seen from **Figure 8** that the HSI method can capture the pressure transient with good accuracy. The slope of each wave nearly equals to the exact solution and there are no spurious oscillations at pressure step change points.

#### Comparisons With Fluent Simulations

In this section, the HSI method was compared with FLUENT simulations in which a two-dimensional laminar flow model was adopted. In the FLUENT simulations, constant outlet pressure, inlet mass flow, and inlet temperature were adopted for the boundary conditions of compressible flow. In order to test the performances of the HSI method under different transient scenarios, three typical transient conditions are studied and

FIGURE 10 | Pressure comparisons of the HSI method and FLUENT under outlet pressure transient using inviscid flow model.

compared with FLUENT simulations, namely an outlet pressure transient, inlet mass flow transient, and inlet temperature transient. Among the three transients, the pressure transient is a fast transient with a very short characteristic time that involves significant changes of pressure, velocity, and temperature. The mass flow transient is also a fast transient but does not involve significant changes of pressure or temperature. The temperature transient is a slow transient with a long characteristic time, and causes mild changes of pressure, velocity, and temperature.

#### Outlet Pressure Transient

The transient condition for outlet pressure transient is a step decrease of outlet pressure from 2.7 MPa to 2.6 MPa. The settings of other parameters are shown in **Table 1**, and the pressure responses at axial position x = 0.2 m are shown in **Figure 9**. It

can be seen that the frequencies of pressure response in both FLUENT and the HSI method are the same, while the amplitude in FLUENT decreases faster than that of the HSI method. The

(C) Temperature responses.

difference in amplitudes is caused by the entrance effect. In the HSI method, a simplified one-dimensional model was employed and the frictional loss caused by viscosity is simplified as a frictional coefficient f = 64/Re. However, in FLUENT a twodimensional laminar flow model was employed and thus there is a greater velocity gradient at the entrance, which leads to a greater energy loss than the fully developed position. Therefore, the pressure amplitudes of FLUENT simulations decrease faster than that of the HSI method at the entrance. If not considering viscosity and uses an inviscid flow model in both FLUENT and the HSI method, the pressure responses of the two methods at x = 0.6 m under the same outlet pressure transient is shown in **Figure 10**. It can be seen that the amplitudes of the two methods are almost the same.

#### Inlet Mass Flow Transient

The transient condition for inlet mass flow transient is a step increase in inlet mass velocity from 10 kg/(m2· s) to 20 kg/(m2· s), and the input parameters are shown in **Table 1**. For pressure, and temperature responses at x = 0.4 m, the comparisons between FLUENT and the HSI method are shown in **Figure 11**. It can be seen that the two methods have the same response frequency, while the amplitudes of the HSI method are slightly greater than that of FLUENT. Compared with the pressure transient in **Figure 10**, the response frequencies of the two transients are the same. There are ∼6.5 cycles within 0.02 s for both the pressure transient and mass flow transient. So the mass flow transient is propagated in the same mechanism as the pressure transient at the speed of sound.

#### Inlet Temperature Transient

The transient condition for inlet temperature transient is a step increase in inlet temperature from 300 to 350◦C, and the input parameter settings are shown in **Table 1**. For the responses of pressure, the mass velocity, and temperature at x = 0.2 m, the results of FLUENT and the HSI method are shown in **Figure 12**. It can be seen from **Figures 12A** that a step increase of inlet temperature results in a mild pressure oscillation, whose amplitude is very small (around 1 kPa). The frequency of the pressure oscillation is the same for both FLUENT and the HSI method, while the amplitude of the HSI method is slightly larger than that of FLUENT. Moreover, because pressure and velocity are strongly coupled, the pressure oscillation also leads to the mass velocity oscillation, which has the same frequency as pressure, as seen in **Figures 12B**. These are fast transients.

It can also be seen from **Figure 12B** that there is a sudden decrease of mass velocity at around τ = 0.04 s, this is a slow transient caused by temperature distribution. At around τ = 0.04 s, the temperature difference signal was transferred to the position x = 0.2 m, as seen in **Figure 12C**. This causes an increase in local temperature and a decrease in local density. According to continuity Equation (1), the decrease in local density leads to a decrease in the local derivative ∂ρ/∂t, which in turn leads to an increase in local space derivative ∂(ρu)/∂x. For the local control volume, this means that the mass flows out of the volume is greater than that flows into it. As a result, the total local mass flow decreases.

Another phenomenon can be seen from **Figure 12B** is that the mass velocity of FLUENT first drops to an intermediate value (around 10.3 kg/(m2· s)) at a steeper rate than that of the HSI method, and then gradually stabilizes at around 10 kg/(m2· s) with a diminishing oscillation. The mass velocity of the HSI method drops directly to the stable value [10 kg/(m2· s)] at a slower rate, and then gradually stabilizes with a diminishing oscillation.

There are two reasons for this phenomenon. First, FLUENT uses a second upwind scheme while the HSI method used a first upwind scheme, so the solution of FLUENT is more accurate than the HSI method, and the temperature distribution is steeper, as seen in **Figure 12C**. This in turn contributes to a steeper decrease of mass velocity in FLUENT. Second, because FLUENT used a two-dimensional laminar flow model, the velocity is zero at the wall and maximum at radial center. This velocity profile leads to the temperature distribution in **Figure 13**, which illustrates the average axial temperature distribution of FLUENT and the HSI method at τ = 0.1 s.

As a result, local temperature changes very quickly in the beginning and then gradually slows down, as shown in **Figure 12C**. This results in the intermediate value of mass velocity in FLUENT simulation, and its gradual stabilization to the stable value. In the HSI method, an averaged onedimensional model was employed and the radial temperature profile is uniform. Therefore, the mass velocity of the HSI method directly drops to its stable value when temperature difference signal arrives.

#### ADAPTIVE TIME STEP

In order to save computational time, two conditions are considered for adaptive time step in the HSI method to speed up calculations, namely the maximum relative density difference between the two densities (one by continuity equation and

the other by state equation), and the Courant–Friedrichs–Lewy (CFL) condition.

#### Maximum Relative Density Difference

The maximum relative density difference is

$$\Delta \rho\_{\text{max}} = \max\_{i=1}^{N} \left( \frac{|\rho\_{i,1}^{n+1} - \rho\_{i,2}^{n+1}|}{\rho\_{i,2}^{n+1}} \right) \tag{17}$$

In which ρ<sup>1</sup> is solved implicitly using the continuity equation and ρ<sup>2</sup> is solved by state equation. 1ρmax is the maximum relative density difference between ρ<sup>1</sup> and ρ2. If 1ρmax ≤ 1ρ<sup>L</sup> (the lower limit of 1ρmax), then the errors under this time step are small enough to allow greater time step for faster calculation while still retaining high accuracy. On the other hand, if 1ρmax > 1ρ<sup>H</sup>

FIGURE 14 | Relationship between temperature and time without considering CFL condition.

(1ρ<sup>H</sup> is the upper limit of 1ρmax), then the relative density error may be too large to converge, and the time step should decrease to ensure convergence. Based on the a lot of testing, the recommended values of 1ρ<sup>L</sup> and 1ρ<sup>H</sup> here are set as 1ρL= 10−<sup>8</sup> and 1ρH= 10−<sup>3</sup> for the HIS method for the transients of the HTRGTS.

### CFL Condition

As mentioned in section Maximum Relative Density Difference, a suitable value of 1ρ<sup>L</sup> saves computational time, and that of 1ρ<sup>H</sup> ensures convergence. However, the temperature response at x = 0.75 m, as shown in **Figure 14**, indicates that the solution is unstable under the outlet pressure transient condition.

To ensure stability, the CFL condition must be considered to control the time step below a certain value, so that the Courant number is lower than its maximum value (cmax). Under a fixed mass velocity of 38.90 kg/(m<sup>2</sup> · s), the relationship between cmax and the number of nodes under different pipe lengths are shown in **Figure 15**. It can be seen that (1) the value of cmax

decreases with the increase of node number, and (2) cmax was not influenced by pipe length. If the node number is >300, the value of cmax ranges between 1.2 and 1.3. To ensure stability, the value of cmax is set as 1.1 for different node numbers. After the CFL condition was considered, the relationship between time step and time is shown in **Figure 16A**, and the temperature response at x=0.75 m is shown in **Figure 16B**. It can be seen from **Figure 16A** that time step first increased to a certain value and then stabilizes, and a stable solution of temperature is obtained, as seen in **Figure 16B**.

#### CONCLUSIONS

To devise a non-iterative transient solver for one-dimensional compressible flow for high temperature gas cooled reactors gas turbine systems (HTRGTSs), the performance of the semiimplicit and nearly-implicit methods in typical compressible flow transients were studied. The results show that the semiimplicit method can capture the fast pressure responses with little numerical diffusion, but it has poor stability in capturing the slow temperature responses. The nearly-implicit method has good stability in capturing temperature responses, but it has too great numerical diffusion in capturing pressure responses.

Based on the discussion of the two methods, a new HSI method that combines the advantages of both the semi-implicit and nearly implicit methods was proposed. In the HSI method, a new calculation strategy is devised: the convective term in momentum equation is treated explicitly to solve pressure and velocity; density and temperature are solved implicitly using continuity and energy equations.

# REFERENCES


Tao, W. (2001). Numerical Heat Transfer. Xi'an: Xi'an Jiaotong University Press.

Toro, E. F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd Edn. Dordrecht; Heidelberg; London; New York, NY: Springer. 115.

The HSI method was verified via the shock tube benchmark problem, and furtherly compared with FLUENT simulations. In verification with the shock tube benchmark problem, it was found that the HSI method can capture pressure transient accurately, and there are no spurious oscillations at the pressure step change points. In comparisons with FLUENT simulations, the outlet pressure transient, inlet mass flow transient, and inlet temperature transient were studied. The results show that the response frequencies of the HSI method are the same as those of the FLUENT simulations, while the amplitudes of the HSI method are slightly larger than those of FLUENT simulations. The HSI method is capable of capturing both the fast transients (such as pressure transient), and the slow transients (such as temperature transient) with good accuracy.

An adaptive time step scheme was furtherly proposed for faster calculations of the HIS method while still retaining a good convergence, stability, and accuracy. Two conditions are considered in this scheme, namely the maximum relative density difference, whose upper limit is set as 10−<sup>3</sup> and lower limit 10−<sup>9</sup> , and the CFL condition, in which the maximum Courant number is set as 1.1.

### AUTHOR CONTRIBUTIONS

XC and MD contributed to the conception of the study. ZM contributed significantly to analysis and manuscript preparation. JW performed the data analyses and wrote the manuscript. JW helped perform the analysis with constructive discussions.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wang, Cao, Meng, Wang and Ding. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Effect of the Wick and the Working Medium on the Thermal Resistance of FPHP

Bin Sun\*, Cheng Peng, Di Yang and Hongwei Li

School of Energy and Power Engineering, Northeast Electric Power University, Jilin City, China

This paper introduces a novel micro-deep, triangle-grooved wick, and compares the performance of the flat plate heat pipe (FPHP) with various types of wicks, including foam metal copper wick, sintered copper powder wick, copper mesh wick, micro-deep rectangle-grooved wick, and micro-deep triangle-grooved wick. Test results verified the feasibility of the nanofluids applied to FPHP. Tested working fluids included deionized water, Cu-water nanofluids, carbon-coated copper nanofluids and multi-walled carbon nanotube (MWCNT) nanofluids. These nanofluids were applied on foam metal FPHP and the results were compared. The results indicated that the metal foam wick is better than the other heat pipe wicks, in the contrast experiment, when the heating power is 100 w, the thermal resistance is 0.1238◦C/K. The thermal resistance of micro-deep triangle-grooved wick FPHP is 15% lower than the micro-deep rectangle-grooved wick FPHP; The optimal filled ratio of 40% can make flat heat pipe have good performance; With the increase of mass fraction of nanofluid, the thermal resistance decreased first and then increased, the best mass fraction is 0.5%, the thermal conductivity of carbon-coated copper nanofluids FPHP is optimal, the average thermal resistance is 0.1369◦C/K, and is 17 and 9% lower than the multi-walled carbon nanofluids and copper FPHP, and is 46% lower than the deionized water FPHP; The gravitational stability of the heat pipe is verified by the experiments.

Keywords: metal foam, nanofluids, capillary wick, flat heat pipe, thermal resistance

# INTRODUCTION

Along with the continued demand for a rapid increase in electronic device cooling, heat dissipation components have also gradually become more miniaturized. Heat pipes for these devices must accommodate a small volume of working fluid. These pipes need to be lightweight and have high heat efficiency (Zou et al., 2016). Most traditional heat pipes are tube types that rely on the capillary structure of the wick to transport heat against gravity, thereby creating a driving force. The heat is transferred through phase transformation, movement of working medium, and evaporation from the adiabatic section to the condenser (Ghanbarpour et al., 2015). The shape of the flat heat pipe (FPHP) makes it difficult for the FPHP to be in full contact with the heating element because the contact thermal resistance is larger and because the heat pipe bending shape has a significant effect on heat transfer performance. The working principle of FPHP is exactly the same as that of the tubular heat pipe, except that the FPHP overcomes the contact problem because of the steam and condensated liquid reverse flow that causes a carry limitation. Realizing long-distance high efficiency heat transfer is difficult. When the work period is longer, the heat pipe size must be larger,

#### Edited by:

Zhaoming Meng, Harbin Engineering University, China

#### Reviewed by:

Han Zhang, Karlsruher Institut für Technologie, Germany Kaiyi Shi, Liupanshui Normal University, China Yandong Hou, Xi'an Jiaotong University, China

> \*Correspondence: Bin Sun sunbin@neepu.edu.cn

#### Specialty section:

This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research

Received: 28 February 2018 Accepted: 19 April 2018 Published: 24 July 2018

#### Citation:

Sun B, Peng C, Yang D and Li H (2018) Effect of the Wick and the Working Medium on the Thermal Resistance of FPHP. Front. Energy Res. 6:37. doi: 10.3389/fenrg.2018.00037 causing it to lose the miniaturization advantage of heat transfer elements (Littwin and Mccurley, 1981). At present, FPHP has got more attention.

Sadeghinezhad et al. (2016) carried out to examine the thermal performance of a sintered wick heat pipe using aqueous graphene nanoplatelets (GNP) nanofluids. It is observed after the experiments that the deposition of GNP creates a coating on the sintered wick surfaces in the evaporator section, and this coating layer increases the surface wettability, thereby enhancing the thermal performance of the heat pipe, The maximum reduction in the thermal resistance of a sintered wick heat pipe filled with 0.1 wt% of GNP is determined to be 48.4% compared with distilled water.

Aly et al. (2016) reports the thermal performance of a helically-micro-grooved heat pipe working with water-based alumina nanofluid, water-based alumina nanofluid with the volume concentration of Al2O<sup>3</sup> nanoparticles of 3.0% are experimentally investigated. Results revealed that the thermal performance of the heat pipe enhances with nanofluid compared to DW, The evaporation and condensation heat transfer coefficients increase as filling ratio increases, but thermal resistance decreases when filling ratio increases.

Tharayil et al. (2016) investigate the heat transfer performance of miniature loop heat pipe with graphene–water nanoflui. An optimum filling ratio of 30% of the total volume of the heat pipe is used in all the experiments. The experimental results indicate that the nanofluids improve the thermal performance of the miniature loop heat pipe and lower the evaporator interface temperature compared to distilled water. And the results confirm suitability of miniature loop heat pipe filled with graphene–water nanofluid for cooling applications.

Mehrali et al. (2016) investigated grooved heat pipes using aqueous nitrogen-doped graphene (NDG) nanofluids. The results indicated that a 90◦ inclination angle corresponds to the best thermal performance, and that the thermal performance of the grooved heat pipe can be improved significantly by using NDG nanofluids.

Tsai et al. (2013) conducted an investigation of standard flat heat pipe spreader and concluded that the tilt angle should be 90◦For maximum thermal resistance.

Hongchuan et al. (2015) proposed a new FPHP that used sintered nano-scale copper powder with a conical capillary wick used on a disc flat plate heat pipe. The steam escapes from largescale pores, which is then sucked up as liquid by small pores. This novel wick improved the heat pipe heat transfer performance. Mizuta et al. (Mizuta et al., 2015) studied a new form of flat heat pipe using neutron radiography. This flat heat pipe has an internal structure with multiple etched copper sheets stacked in a specific way. Li et al. (2010) observed the working medium flow condition inside a heat pipe by using infrared thermal imaging method.

Lefèvre et al. (2012) investigated a double wire mesh screen wick FPHP and a micro grooved FPHP covered with screen mesh. The working medium was methanol. After testing the performance of these two kinds of FPHP by changing the experiment conditions, they determined that nucleate boiling caused the reduction in efficiency of the capillary structure with double wire meshes at a certain angle. This case was only suitable for low heat flux. The unfavorable tilted position is not possible in micro-grooved FPHP.

Heat pipes are not confined to a single metal material. Hsieh and Yang (2013) used a No. 250 double copper mesh for the wick of FPHP and a copper and silicone rubber mixture for a more flexible pipe wall to create a new kind of polymer-based flat heat pipe capable of delivering a higher power of 12.67 W. The increased power was attributed to the vertical plane bending angle of 15◦ .

Li et al. (2010) used infrared thermal imaging method to compare different fins according to their size and number, as well as the Reynolds number of FPHPs with aluminum heat sinks and vapor chamber heat sinks. The thermal resistance of vapor chamber heat sink is reduced with the increase of Reynolds number, and the impact strength decreases. The vapor chamber heat sink outperformed the aluminum heat sink because the base plate of the vapor chamber heat sink has a more uniform heat distribution. Lower Reynolds number with appropriate fin number can yield the best heat conduction in a vapor chamber heat sink, and when the Reynolds number is high, thermal conductivity improves as the number of fins increases.

Moraveji and Razvarz (2012) used Al2O<sup>3</sup> (35 nm) and water as working fluid to test the thermal efficiency of a straight copper tube and circular tube with a 1 mm sintered wick. They found a 90◦ Curve in the circular heat pipe between the evaporator and condenser sections. They also found that charging the nanofluid in the circular heat pipe can significantly reduce wall temperature, reduce thermal resistance, and improve thermal efficiency when compared to conducting the test using only pure water as the working fluid.

Kim et al. (2016) used alumina-acetone as a flat heat pipe working medium, and compared the effect of three different nanoparticle shapes (sphere, brick, and cylinder). The thermal resistance of pure acetone FPHPs as working medium failed to sustain the performance of FPHPs using the shaped nanoparticles. Their improved performance was indicated by a decreased thermal resistance of 33% for the sphere shaped nanoparticles, 29% for the brick shaped nanoparticles and 16% for the cylinder shaped nanoparticles as compared to the alumnaacetone FPHP.

Wu et al. (2016) used the concentrations of 0.1, 0.2, and 0.3% of C60-ethanol nanofluids used in a flat-plate closed loop pulsating heat pipe (FCLPHP). The results show that although heat pipe performance improved, the heat load was reduced, causing the FCLPHP to dry out. The final analysis indicated that the heat load had a greater effect on thermal performance than the C60 nanofluid concentration.

Saleh et al. (2013) used the synthesis of ZnO-ethylene glycol nanofluids as a heat pipe working medium and found that an increase in the concentration and crystallite particle size resulted in a decrease in heat pipe surface temperature distribution and thermal resistance.

Sarafraz and Hormozi (2014) used an ultrasonic homogenizer, electromagnetic stirrer, and pH control system with constant temperature hot bath in preparation for testing different particle sizes in Al2O3-water-ethylene glycol and Al2O3-waterdiethylene glycol nanofluids. Their experiments were carried out using these working media operating in straight tubes. They found that the heat transfer coefficient of the heat pipe increased significantly.

The above studies show that FPHP is a primary concern in the research of heat pipe, and that this study has significant research value. Most domestic and foreign studies have focused on the shape of the heat pipe and structure of the wick. Improvement of the wick has also been a continual focus of study in heat pipe technology.

The structure of the wick plays a decisive role in the performance of heat pipe (Hassan and Harmand, 2015; Bhullar et al., 2016a,b; Su et al., 2016; Zhou et al., 2016; Nazari et al., 2018; Ozsoy and Corumlu, 2018). The metal foam is a metal substrate containing a certain number of pores with a certain diameter size that makes it porous. Used as a wick, this metal foam can improve the heat transfer performance of spread and other equipment and provide a new area for the heat pipe heat transfer enhancement in the field of heat pipe machining (Su et al., 2016), some of the most common types of wicks as groove type, copper screen, and sintered copper powder wick are also studied (Zhou et al., 2016) the working medium is also important to the FPHP.

Through the experiments, this paper studies the influence of variety wire and medium of FPHP. Experiments under the feasibility given conditions in this study indicated the effects of the wire and working medium on the performance of the heat pipe, which can provide a reference for heat pipe industrial production and processing.

# NANOFLUID PREPARATION AND STABILITY ANALYSIS

#### Nanofluid Preparation

Experiments using Cu nanofluids, multi-walled carbon (MWCNTs) nanofluids and carbon-coated copper water nanofluids, deionized water as based, and sodium dodecyl benzene sulfonate (SDBS) as dispersant were performed (Wusiman et al., 2013).

The "two-step" method (Haddad et al., 2014) was used and with a proportion of 1:1, the nanoparticles and dispersant were combined. At this point in the experiment, the stability of the nanofluids was the best (Hajian et al., 2012). The mass fractions of the three kinds of nanofluids were configured to 0.1, 0.1, 0.2, 0.4, and 0.5%. Oscillation in the ultrasonic oscillation device was performed for 60 min to obtain stable nanofluids. **Table 1** is the parameters of the selected nanoparticles.

# Analysis of the Stability of Nanofluids

In the present experiment, transmittance was used to evaluate the stability of three nanofluids. The transmittance was measured with a visible spectrophotometer (Model 721-100, Shanghai Quintessence Technology Instrument, Co., Ltd). This process is repeated three to four times, and the measured results are shown in **Figure 1**.

**Figure 1** display that, as time increase, nanofluids gradually precipitate, because the more uniform distribution of nanofluids,the lower light transmittance we get. This nanofluids TABLE 1 | The parameters of the nanoparticles.


is difficult to deposit in the long term and is suitable for subsequent experimentation and application.

# EXPERIMENTAL DEVICE AND EXPERIMENTAL METHOD

#### Experimental System

The experimental system depicted in **Figure 2** consists of a heating control system, a data acquisition system, a heat radiation system, and an experimental section. The heating zone, which is 20 mm in diameter, is composed of a brass column with three heating pipes of 100 W inserted below. The upper surface is used to heat the center of the FPHP. K-type sheathed thermocouples, 1 mm in diameter, are used to measure and data acquisition, and is directly connected to the data collector. Above the heating zone are the pillars from the heating surface in distances of 5, 15, and 25 mm and two holes of 1 mm diameter, which will be used to measure the temperature of the center line of pillars. Two thermocouple measurements were averaged as the point value of the temperature. A thermocouple is placed at the center of the lower surface of the heat pipe to measure the maximum temperature of the heat pipe. Eleven thermocouples were placed within the radius of the heat pipe upper side with equal spacing, as shown in **Figure 3A** 120 × 120 × 45 mm aluminum heat sink was placed on top of the heat pipe. The top of the heating zone is equipped with 3 W fans for forced convection heat transfer.

**Figure 3** shows the structure of the flat heat pipe and the distribution of temperature measurement points. In this experiment, to measure the surface temperature, the thermocouple and surface of heat pipe are fixed by SatlonD-3 and 606 catalysts. The thermocouple is in contact with the temperature measurement point at the surface, thereby reducing measurement error. Brass heating temperature is measured at the point of contact with the thermocouple and temperature measurement surface. Thermal grease was applied to reduce thermocouple thermal contact resistance. Insulation asbestos was placed around the heating element to insulate the process and lessen heat loss to the point of being negligible. Surface heating and cooling surfaces were evenly coated with thermal paste and fixing and pressing to reduce thermal contact resistance. Experimental data were collected by the data logger and transferred to a computer for monitoring and recording.

In this study, the working fluids included deionized water, Cu-water nanofluids, MWCNTs-water nanofluids, and carboncoated copper-water nanofluids. Before the experiment, a vacuum pump with pressure of up to 6 × 10−<sup>2</sup> Pa was used to vacuum the flat heat pipe, then a syringe was used to pack the amounts of liquid. Five kinds of filling ratio with values of 0%, 8%, 20%, 40%, 60% and 100% were used. The flat heat pipe placed at angle θ were 0, 30, 60, and 90◦ , and the level of placement of the heat pipe was 0◦ . The heating power used is from 0-200 W. When the temperature on the surface of the heat pipe near the center in 10 min fluctuation was <0.5◦C, the heat transfer can be considered to have reached steady state.

The flat heat pipe working fluid flow was vapor-liquid twophase flow, thus, the steam that the working fluid vaporization generated flowed to the condensing end of the heat pipe under slight pressure differential, when the flat end of the heat pipe is heated and evaporated. During the release of latent heat of vaporization and condensation, the working fluid condensed at the wick capillary returned from the condensing end to the evaporating end and achieved transfer of heat from the evaporating end to the condensing end. The normal operation of the flat heat pipe consisted of evaporation of the liquid, flow of steam, condensate, and steam condensate reflux composition. The operation must satisfy the following formula:

$$
\Delta p\_{\rm c} \ge \Delta p\_{\rm l} + \Delta p\_{\rm v} + \Delta p\_{\rm g} \tag{1}
$$

Where 1p<sup>c</sup> is the driving force of the working fluid circulating inside the heat pipe, which was used to overcome the vapor from the evaporator side pressure drop to the condensing end in flow 1p1, the pressure drop of reflux liquid from the condensing end to evaporating end 1pv, and pressure drop because of gravity 1pg. The largest of the capillary force of the foam metal can be obtained as (Peng Y. et al., 2013).

$$
\Delta p\_{\rm c} = \frac{4\sigma \cos \alpha}{d\_{\rm eff}} \tag{2}
$$

where σ is surface tension, α is a wetting angle, and deff is the effective aperture capillary force of the metal foam. In general, the deff is affected by the minimum diameter. In this experiment, because of the additional shortening of the liquid return path, the liquid reflux flow rate remained small and coupled with smaller pore size. The liquid flow Reynolds number was also small. The pressure drop of the liquid flow is calculated using Darcy's available (Hsieh et al., 2008):

$$
\Delta p\_{\rm l} = \frac{\mu u\_{\rm l}}{K} \Delta x \tag{3}
$$

where µ is the viscosity of the liquid, 1x is the fluid return path length, and K is the metal foam permeability coefficient. The literature shows that

$$K = \frac{\varepsilon d\_p^2}{32} \tag{4}$$

In equilibrium, the mass flow rate of the liquid evaporates and the liquid supplements are the same

$$q' = \frac{1}{8} \frac{\sigma \rho\_l h\_{\text{fg}}}{\mu} \frac{A\_{\text{l}} \varepsilon d\_{\text{p}}^2 \cos \alpha}{A\_{\text{e}} \Delta x d\_{\text{eff}}} = \frac{1}{8} NG \tag{5}$$

FIGURE 2 | Schematic of the experimental system. 1. Power supply; 2. Automatic voltage regulator; 3. Power meter; 4. Heating pipe; 5. Heating pillars; 6. Angle adjustable test bed; 7. Fixtures; 8. Thermocouple; 9. Flat plate heat pipe; 10. Heat sink; 11. Cooling fan; 12. Data acquisition; 13. Computer. (A) Experimental system diagram. (B) Physical map of the test loop.

where N = σ ρlhfg µ is related with working fluid, called the working medium quality factor. G = Alεd 2 p cos α Ae1xdeff is related to the wick and heat pipe structure. A flat heat pipe filled with working fluid should have high surface tension and latent heat of vaporization, because the former can provide a large capillary force, whereas the latter causes a certain amount of heat transfer fluid flow rate that requires less liquid reflux to reduce resistance. Lesser viscous liquid means higher quality of the working fluid.

The performance of the heat pipes can be produced by thermal conductivity and heat resistance. The thermal conductivity and heat resistance of the improved tablet wick was compared to the copper metal substrate. The heat transfer performance of the flat heat pipe was also analyzed. Flat plate heat pipe thermal conductivity is axially thermal. This parameter can be characterized through the temperature difference between the upper and lower surfaces of the heat pipe. Equivalent thermal resistance of heat conducting element is expressed as:

$$R = \frac{T\_{\text{max}} - T\_{\text{avg}}}{Q} \tag{6}$$

where Tmax is the maximum surface temperature of the heat pipe at the center of the lower surface temperature, Tavg is the average

temperature on the higher surface of the heat pipe, and Q is the lower surface of the heat pipe with heated copper cylinder heat input. Heat flux density q (W/m<sup>2</sup> ) can be calculated using the heat input. The heat flux formula is as follows:

$$q = \frac{Q}{A} \tag{7}$$

The use of heat insulation asbestos pillars resulted in minimal heat, which can be ignored. However, to reduce the reading errors of power matter as well as avoid cooling errors, heat flux can be obtained by the Fourier equation (Ji et al., 2012):

$$q = k \times \frac{\delta T}{\delta x} \tag{8}$$

where, k is the thermal conductivity of copper. The copper used in this experiment had a thermal conductivity of 401 W/(m·K). Through four temperature measurement points on the higher surface, the distribution of temperature on the surface on the plate can be obtained. Uniformity of temperature distribution reflects the average heat of the heat pipes. When the radial temperature gradient of the heat pipe is reduced, the temperature difference of the point on the surface is smaller, to that the surface temperature distribution is very uniform, thereby improving thermal uniformity.

The thermal resistance of the heat pipe is composed of two parts, pure phase changes thermal resistance and thermal resistance. To study the thermal conductivity of the filling heat pipe, the pure thermal resistance of the heat pipe should be measured when it is not filled. For comparison, the heat pipe and aluminum plate used in the same test must have the same sizes. After clamping fixtures, the power input heating tube is then connected. The experiment begins with a heating power of 10- 200 W, and a power interval of 10 W. The temperature changes after each increase in heating power. When the respective measurement points are less than every 5 min, 0.1◦C is reached, and is regarded as the steady state.

To enhance the precision of the experimental results, the system errors were analyzed and measurement uncertainty evaluated (Ghanbarpour et al., 2015). The uncertainty of the system mainly comes from the uncertainty of machining error, heat flux and thermal resistance. The computation formula of the uncertainty is as follows:

$$\frac{\Delta q}{q} = \sqrt{\left(\frac{\Delta k\_{\text{copper}}}{k\_{\text{copper}}}\right)^2 + \left(\frac{\Delta \delta T}{\delta T}\right)^2 + \left(\frac{\Delta \delta x}{\delta x}\right)^2} \text{ (9)}$$

$$\frac{\Delta R\_{\text{equ}}}{R\_{\text{equ}}} = \sqrt{\left(\frac{\Delta \delta T}{\left(T\_{\text{max}} - \overline{T}\_{\text{top}}\right)\_{\text{min}}}\right)^2 + \left(\frac{\Delta q}{q}\right)^2 + \left(\frac{\Delta A}{A}\right)^2} \text{ (10)}$$

kcopper is the heat conductivity coefficient of copper, it is consult from the standard property sheet, the uncertainty can be ignored, The uncertainty of temperature is caused by the measurement error of thermocouple, Machining error was 2%, Equation (9) is used to calculate the uncertainty of heat flux; The thermal resistance calculation formula shows that the uncertainty of thermal resistance is depends on the temperature and heat flux, the uncertainty of thermal resistance can be calculated by Equation (10). The uncertainties and parameters of the instruments used in the experiment are listed in **Table 1** while the uncertainty of the variables is presented in **Table 2** and uncertainty of the experimental variables is shown in **Table 3**.

#### RESULTS AND ANALYSIS

# Heat Balance Test of the Experimental Period

Before the experiment, the heat balance of the experiment system should be discussed by calculated the deviation between input power Q and Q<sup>l</sup> , which Q is measured by power meter and Ql is calculated by the temperature nearby the surface of the heating copper column. Heat balance deviation can be expressed as follows (Zhao et al., 2016):

$$\varepsilon = \frac{|Q - Q\_1|}{0.5 \ (Q + Q\_1)} \tag{11}$$

The heat balance deviation is got by computing and it is <5%, thus the total power of the working medium can reach close to the set power of 100 W. This shows that the thermal insulation device performance is good, the heat loss of the experiment can be ignored.

# The Surface Temperature Distribution of the Heat Pipe

**Figure 4** shows the spread performance of four metal foam FPHPs with different working mediums. The experiment shows that when the heating power is set to 100 W, the heat pipe started sufficiently and achieve a lower temperature, so choose


TABLE 3 | Experimental variables and uncertainty.


the heating power of 100 W to contrast the surface temperature distribution of FPHP. the mass fraction of all nanofluids is 0.5%, and the charging rate is 40%. As shown in the figure, when the heating power is 100 W, the maximum temperature difference of carbon-coated copper-water nanofluids, Cu-water nanofluids, MWCNT-water nanofluids and deionized water heat pipe on the surface is respectively 1.2, 2.8, 3.3, and 4.1◦C, the top temperature of the copper plate is nonuniform distribution, the maximum temperature difference can reach 10.4◦C, the FPHP filled of nanofluids is better than the deionized water heat pipe, and the FPHP is better than the copper plate. The carbon-coated copper-water nanofluids shows the best performance and has the smoothest curve compared to others, the temperature of the measurement points is low, Cu-water nanofluids heat pipe had the second-best performance. Meanwhile, compared with the previous two kinds of nanofluids, the spread performance of the multi-walled carbon nanotube-water heat pipe was poor. And the Cu nanofluids had larger fluid viscosity, and the viscosity had larger effect on the flow resistance of heat pipe system. From the surface temperature, the upper surface average temperature of carbon-coated-water heat pipe is 70.3◦C, is 4.3◦C lower than copper plate, that means the overall temperature of the heat pipe is lower than copper plate of the same size, also account for the cooling effect is better. Multi-walled carbon nanotubes-water nanofluid had both advantages of two nanofluids, high thermal conductivity and smaller viscosity of carbon nanotubes, and thus, the carbon-coated copper nanofluids heat pipe had the lowest thermal resistance.

# Thermal Resistance Analysis

Many factors affect heat transfer, including liquid ratio, capillary suction fluid wick structure and materials, heat pipe placed angle, internal working medium, heating power, etc. (Chen et al., 2013). Heat transfer performance evaluation was based on thermal resistance of heat pipe. When the thermal resistance is smaller, a better heat transfer performance can be obtained, and the main

factors affecting it as well as the transient temperature can be analyzed from the following several aspects.

#### (1) Effect on thermal resistance of the working medium

**Figure 6** shows thermal resistance under different powers of the four working mediums at 40% charging rate. The mass fraction of the three kinds of nanofluids is 0.5%, and the wick of the heat pipe is foam metal. As shown in **Figure 5**, for the four different working mediums, thermal resistance was larger heating power is small. Thermal resistance of the flat heat pipe decreased with the increase in heating power. This result can be attributed to the working medium inside the flat heat pipe producing heat absorption through the phase change with the gradual increase in heating power. Relative to the conduction and convection heat transfer, phase change heat transfer has the advantages of smaller heat transfer temperature difference and better heat transfer property. In addition, under the same heating power, flat heat pipe filled with nanofluids has smaller thermal resistance because the working medium added nanoparticles. The nanoparticles increased the surface area and thermal capacity, which in turn increased the coefficient of thermal conductivity of the working medium (Kumaresan et al., 2014). The collision and interaction between particles, water, water vapor, and metal surface, enhanced the heat transfer of the heat pipe. The carbon-coated copper nanofluids were chosen as the working medium, and the charging ratio is 40%. When the heating power is at 150 W and the thermal resistance is 0.1◦C/W, the foam metal designs for the wick with nanofluids are chosen as working medium of flat pipe and had good axial heat conduction.

(2) Effect on thermal resistance of the mass fraction of the working medium

The heat transfer efficiency of carbon-coated copper is better than that of other heat pipe working mediums. In this paper, the mass fraction is studied by using carbon-coated copperwater nanofluids heat pipe. **Figure 6** shows the thermal resistance of different mass fractions of carbon-coated copper-wanter nanofluids filled in foam metal FPHP. The heating power was 100 W, and the charging rate was 40%. With the increase of the mass fraction, the thermal resistance of foam mental FPHF is decreased first and then increased. In addition, the increase rate of the thermal conducted performance was reduced with the gradual increase of mass fraction. From the figure, the thermal resistance of deionized water heat pipe is 0.1391◦C/W. The figure shows that the thermal resistance of carbon coated copper-water nanofluids heat pipe with the mass fraction of 0.5% is 0.1121◦C/W. The increase of the mass fraction of the Nanoparticles caused high thermal efficiency, which could significantly improve the performance of the heat pipe. However, with the increase in mass fraction of the nanofluids, the viscosity of the working medium inside the heat pipe increased significantly and was averse to the heat transfer performance. According to the stability experiments of nanofluids introduced in section Introduction, precipitation will occur with the increase in mass fraction of nanofluids. Hence, a specific working medium concentration of nanoparticles that optimizes the heat transfer property exists. In the experiments, the thermal resistance of 0.5% of mass fraction minimum, when the mass fraction is 0.6%, heat pipe thermal resistance increased, that indicated the best quality score is between 0.5 and 0.6%.

(3) Effect on thermal resistance of the charging rate of working medium

The last group of experimental results show that the carboncoated copper foam metal heat pipe has the lowest thermal resistance, and hence, carbon-coated copper-water nanofluids was chosen as working medium to discuss the effect of different charging rates of filling ratio on heat resistance. The mass fraction of carbon-coated copper nanofluids is 0.5%. **Figure 7**, under the low heating power, the flat plate heat pipe of 20% charging ratio has smaller thermal resistance than that of the other liquid ratio, and as the heating power increases gradually, the thermal

resistance of all charging ratio is reduced. That is because with the heating power increasing, the boiling liquid film thickness at phase interface inside the heat pipe will reduce gradually, and the phase change thermal resistance will decrease until it reaches a certain value. In this experiment, 40% liquid filling ratio of the heat pipe shows good heat transfer effect during the high heating power, and over a wide range of heating power, it has maintained stable thermal resistance. When the heating power is at 150 W, the heat pipe carried the heat flux density of around 47.77 kW/m<sup>2</sup> , and the heat pipe boil dry phenomenon did not happen. When the heating power was adjusted from 150 W to 160 W, the measured temperature from the three thermocouples inside the copper heating column reached more than 30◦C in a matter of seconds. This phenomenon shows that flat heat pipe internal dry phenomenon has occurred, and the heat transfer is deteriorating.

#### (4) Effect on thermal resistance of the tilt angle

**Figure 8** shows the effect of the tilt angle of the carbon-coated copper nanofluids FPHP with mass fraction of 0.5% and charging rate of 40%. When the heating power is small and the inclination angles are 60 and 90◦ , thermal resistance value slowly decreased along with the increase in heating power. When the heat pipe is in a horizontal position or the tilt angle is small, thermal resistance rapidly decreased with the increase of heating power. When the heating power is within the range of 40-100 W, the thermal resistance of the flat plate heat pipe is at a minimum, whereas when the FPHP is upright, the resistance increased. However, with the further increase of the heating power, the heat resistance of the four angle conditions of heat pipe gradually converged, and when the heating power is nearly 150 W, the heat transfer deterioration did not occur.

#### (5) Effect on thermal resistance of the wick structure

As a key part of the heat pipe, the wick has three main functions: as a working medium flow channel, as a working medium of the phase interface, and as provider of backflow capillary force for the working medium. The effective pore radius of the wick should be smaller to guarantee the backflow power of the working medium. The wick should also have sufficient permeability to reduce the pressure drop of the backflow (Li et al., 2016). Wick thermal resistance should be smaller to reduce radial thermal resistance.

The wick is an important part of heat pipe. The most common wicks include screen, sintering, and groove type wicks. The performance of the grooved structure applied to FPHP is stable and does not require sintering. It can also conduct heat dissipation through the fin in the steam chamber. The experiment perfected and proposed a micro deep triangle-grooved wick, which involved cutting the copper board through three directions at an angle of 60◦ . The porosity of this kind of wick is larger, leading to a tri-prism unit structure with smaller pore diameter. The steam chamber increases the condensation area, which is good for steam condensation. Porosity and capillary force are increased by reducing the width of the pore and increasing the quality of porosity. The edge of the micro prism of groove structure with smaller angle has better infiltration of liquid than

different mediums of FPHP.

that of Right Angle. **Figure 9** shows the improvement of micro deep triangle-grooved wick and micro deep rectangle-grooved wick structure.

The porosity of copper screen wick and sintered copper wick can be calculated:

$$\theta = \frac{V\_{\text{pore}}}{V\_{\text{total}}} \tag{12}$$

$$V\_{\text{pore}} = V\_{\text{total}} - V\_{\text{copper}} \tag{13}$$

θ is the porosity of the wick, Vpore is the velocity of the internal porosity of the wick, Vtotal is the whole velocity of the wick. Vpore = Vtotal− Vcopper, the Vtotal is calculated by Equation (14):

$$V\_{\text{copper}} = \frac{1}{2}\pi D^2 \delta + \frac{\pi}{4}l[\left(D - 2\xi\right)^2 - \left(D - 2\xi - 2\delta\right)^2]l \tag{14}$$

Among them, D is the external diameter of the FPHP, δ is the thickness of the wick, ξ is the thickness of the pipe, mcopper is the mass of the wick.

TABLE 4 | The parameters of the sintering wick.


#### TABLE 5 | The parameters of the micro-deep grooved wick.


In this experiment, all the heat pipe is of the same size, and the shape is disc, its diameter is 90 mm, the thickness of the plate is 9 mm, the thickness of the tube wall on the top and lateral is 2 mm, the thickness of the lower tube wall is 1.5 mm. The diameter of copper tube connected with the pumping air hole is 3 mm, since the cut-off point is very close to the heat pipe, its influence can be neglected. Several parameters of the wick as the following **Tables 4**, **5**:

**Figure 10** shows the heat pipe using carbon-coated copper nanofluids as working medium, which had a mass fraction of 0.5%, and filling rate of 40% as applied to copper screen wick, sintered copper wick, foam metal wick, and two kinds of micro grooved wick of flat heat pipe. Under low heat input, when the input power is 30 W, the micro deep rectanglegrooved wick heat pipe had the lowest thermal resistance of 0.1897◦C/W, and the micro deep rectangle-grooved wick heat pipe had the second lowest thermal resistance of 0.1913◦C/W. The thermal resistance of the sintered copper heat pipe and foam metal was lower than that of the copper mesh heat pipe. When the heat input was increased, the change of heat transfer performance of the micro deep rectangle-grooved wick heat pipe was unobviously, the minimum thermal resistance is 0.1674◦C/W. However, with the increase of heat flux density, the thermal resistance of micro deep triangle-grooved wick is only 0.1460◦C/W, and it is obviously lower than the traditional heat pipe. Relative to the structure of the three wicks with larger cavity, the variations in thermal resistance of the sintered copper wick heat pipe is smaller compared to others, and the heat transfer performance of the foam metal and copper network wick further increased, causing the thermal resistance of the metal foam FPHP to be at the minimum. This result can be attributed to the groove type heat pipe having upper and lower connections of miniature frame column in its chamber; the rib can conduct axial heat directly. Because of the smaller internal cavity, the heat pipe can be used under low heat flux, and can easily reach a stable state; however, it also brings forward the boiling limit (Peng H. et al., 2013). The micro deep triangle-grooved wick heat pipe increased the space inside the heat pipe, had better porosity, and reduced the volume of the frame column. Thus, the thermal performance under low heat flow density is better than the traditional type. With better porosity, the heat pipe could withstand higher heat flux, and under low heat flux density, the heat pipe also had lower thermal resistance. The nanofluids can achieve good results in the appropriate aperture of the deep groove of the heat pipe, which strengthened the heat exchange through the fins. Relative to the cavity structure, the heat pipe with micro deep grooved wick remained stable within a certain range of heat flux. The size of the aperture of the bronze lotion wick was larger, making the effective radius bigger and the capillary force smaller, and causing the start-up temperature of heat pipe demand to become higher. The heat transfer performance is poor under low heat input. The sintered copper has a smaller aperture, which can effectively provide better capillary force. The small aperture will cause a blocking phenomenon when the working medium used is colloid containing nanoparticles. The small aperture will also lead to reduced permeability of the working medium, which increased the backflow resistance. Under high heat input, it cannot further enhance the heat transfer effect. Nanofluids in metal foam wick with appropriate aperture will not easily cause a blocking phenomenon, and can provide good capillary force, which results in high permeability.

Frontiers in Energy Research | www.frontiersin.org

transfer resistance, the following conclusions can be reached:

Therefore, the FPHP of metal foam wick has good heat transfer

# REFERENCES

properties.


#### AUTHOR CONTRIBUTIONS

BS: Provided experimental design ideas and overall project planning. CP: Carried on experiment, data processing and analysis. DY: Carried on experiment. HL: Digital signal process.

#### ACKNOWLEDGMENTS

This study was supported by the Program for Science and technology development projects of Jilin province (20160101282JC, 20160520032JH, 20170101123JC) and National Natural Science Foundation of China (No. 51406031).


sintered wick heat pipe using CuO nanofluids. Int. J. Heat Mass Transf. 72, 507–516. doi: 10.1016/j.ijheatmasstransfer.2014.01.029


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Sun, Peng, Yang and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.


# Investigation on the Subcooled Boiling in Vertical Pipe With Uncertainties From Boundary Conditions by Using FLUENT

Xiang Zhang, Rui Zhang, Xindi Lv and Tenglong Cong\*

*Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China*

Subcooled boiling flow taking place in the reactor system plays a critical role in the safety of nuclear power plants. It has been studied by experiments and system codes in the past decades. Now subcooled boiling can be predicted with CFD code based on the Eulerian two-fluid model, with the development in the computational technology and the understanding in the mechanism of two-phase flow. The published works on the validation of CFD code for two-phase flow were carried out based on the deterministic analysis by comparing the calculated and experimental nominal inputs and outputs, which is not sufficient for code validation since it didn't consider the inevitable uncertainty in the experiment measurements. In the current work, subcooled boiling was predicted by using a CFD code, FLUENT, with consideration of the uncertainties of boundary conditions. The resultant parameters with uncertainties were compared to the experiment data for validation purpose. Confidence intervals of the two-phase parameters were predicted. Besides, correlations between the boundary conditions and the outputs were analyzed.

#### Edited by:

*Jun Wang, University of Wisconsin-Madison, United States*

#### Reviewed by:

*Zeyong Wang, Rensselaer Polytechnic Institute, United States Xianping Zhong, Tsinghua University, China*

> \*Correspondence: *Tenglong Cong tlcong@hrbeu.edu.cn*

#### Specialty section:

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

Received: *02 February 2018* Accepted: *16 March 2018* Published: *11 April 2018*

#### Citation:

*Zhang X, Zhang R, Lv X and Cong T (2018) Investigation on the Subcooled Boiling in Vertical Pipe With Uncertainties From Boundary Conditions by Using FLUENT. Front. Energy Res. 6:23. doi: 10.3389/fenrg.2018.00023*

Keywords: subcooled boiling, uncertainty analysis, boundary conditions, two-phase CFD, Eulerian two-fluid model

# INTRODUCTION

In the past two decades, CFD codes based on the Eulerian two-fluid model have been employed to simulate the two-phase flow in the nuclear reactor system. The predicted results were compared with the experiment data to validate and assess the performance of codes. Good agreements have been obtained by tuning the input parameters or by altering the interaction models (Krepper and Rzehak, 2012). However, tens of parameters have been employed in the two-phase CFD code to model the boundary conditions and the phasic interactions (Zhang et al., 2015a). The accuracy and confidence of the two-phase CFD codes are still questionable due to the complexity of the interactions at heated surface and between liquid and bubbles. Besides, just like the system code, there might be uncertainties in the prediction of two-phase flow by CFD codes, which should be taken into consideration during simulation (Bestion et al., 2016). Besides, the accuracy and reliability of the codes should be assessed with considering the uncertainty introduced into the codes.

The sources of uncertainties in CFD codes are similar to these in system codes, including the boundary conditions, the physical models, the model parameters, the numerical methods, the geometry simplifications, the model simplifications, and the physical phenomena that have not been considered in the simulation (Bestion et al., 2016). From the point of view of code validation,

**37**

the uncertainties introduced by the boundary conditions should be analyzed in the first place since the experiment data are obtained with inevitable uncertainties.

In the process of code validation with considering the inputs and outputs uncertainties, rather than the nominal values, multiple sets of data that represent the uncertainties, i.e., the statistical characteristics of the inputs, are employed as the inputs of the code. Besides, the statistic of interested results can be obtained by analyzing the statistical characteristics of the resultant scatter data. Then, the calculated results with uncertainties should be compared with the experimental data with error bands. In general, the principles of uncertainty analysis for two-phase CFD codes are exactly the same with these for system codes. However, the uncertainty analysis in the field of CFD codes is far from mature, which is not like that in the system codes. Knio and Le Maitre (2006) made a brief review on the uncertainty analysis in CFD problems with polynomial chaos methodology and claimed that the combination of polynomial chaos method and Navier-Stokes equations for CFD uncertainty analysis was not mature. Badillo et al. (2013, 2014) and Prošek et al. (2016) modeled the GEMIX experiment for fluid mixing with consideration of the uncertainties from inlet boundary conditions by using the traditional uncertainty propagation method with random sampling. Bestion et al. (2016) summarized the applications of uncertainty analysis to single phase CFD simulations in the field of nuclear thermohydraulics and proved the high applicability and robustness of uncertainty propagation methodology when coping with CFD problems.

In the current work, a preliminary work to estimate the uncertainties in the subcooled boiling modeling was presented. The uncertainties on the two-phase parameters caused by the input boundary conditions have been investigated by using an uncertainty propagation method coupled with an improved Monte-Carlo based sampling method, Latin Hypercube Sampling (LHS) (Helton and Davis, 2003). The subcooled boiling experiment in a vertical pipe carried out by Bartolomei and Chanturiya (1967) was chosen as the benchmark test. The sources of uncertainties are introduced from the test boundary conditions, including inlet mass flux, the inlet temperature, the system pressure, and the wall heat flux. Effects of input uncertainties on the pressure drop, outlet vapor fraction, averaged wall temperature, net vapor generation (NVG) location, and the localized two-phase parameters were analyzed.

#### EXPERIMENT CASE TO ANALYZE

The boundary conditions employed in the numerical simulation come from the experimental measurements, such as the pressure, mass flow rate, and inlet temperature. The measurements in the experiment will have some uncertainty inevitably. The uncertainty in the boundary conditions can be transferred to the whole computation domain by the governing equations and then introduce uncertainties into the predicted results. In traditional, the nominal value or the most expected values are chosen as the boundary condition values specified in the simulation. However, the specification of definite nominal values without considering the uncertainty of inputs cannot predict the results with expectation and probability and cannot estimate the model accuracy reliably.

In this work, the propagation of uncertainty introduced by boundary conditions was analyzed based on the subcooled boiling experiment proposed by Bartolomei and Chanturiya (1967). In this experiment, subcooled boiling parameters of water-vapor two-phase flow were measured in a vertical pipe, including the vapor volume fraction, liquid temperature, and wall temperature. Besides, the boundary conditions were also recorded, including the inlet temperature, inlet mass flow rate, system pressure, and wall heat flux. The error bands of boundary conditions for inlet mass flux, wall heat flux, system pressure, and inlet temperature are 2.0, 3.0, 1.0%, and 1 K, respectively. However, the distributions of all the parameters in the relative error bands were not presented. Given that most of the distributions follow the normal distribution in the natural and social sciences, as noted in statistics theory (Chiasson, 2013). Thus, the distributions can be assumed to be normal if the real ones are unknown and this assumption works good in statistics (Chiasson, 2013). Besides, it should be noted that, the distribution of input parameters will affect the statistical characteristics for the outputs, such as the standard deviation, the 95% confidence interval and the error band. However, when compared to the other distribution profiles, such as uniform distribution, the normal distribution can give more credible results since it can present the statistical characteristics of the input parameters. Along with the statistical characteristics of the outputs, this work also focuses on the relationships between the inputs and the outputs, that is the correlation coefficients between inputs and outputs, which are independent from the distribution profiles of the inputs. Based on the above considerations, in this work, the random errors for all the input parameters are assumed to follow the normal distributions, that is,

$$X \sim N(\mu, \sigma)^2 \tag{1}$$

The probability density function (PDF) is:

$$f(\mathbf{x}) = \frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{(\mathbf{x}-\mu)^2}{2\sigma^2}}\tag{2}$$

where µ , σ 2 , and σ are the expectation, variance and standard deviation, respectively. The probability that the random parameters locate in the range of +3σ is larger than 99.73%. Thus, we can assume that the error band of the measured parameter is equal to the interval of +3σ for each variable. Following this assumption, the statistical parameters of all the boundary parameters can be obtained, as shown in **Table 1**.

#### ANALYSIS METHODOLOGY

# Brief Introduction on the Uncertainty Analysis Methodologies

The uncertainty of a simulation can be estimated by uncertainty propagation method (Badillo et al., 2013), the accuracy extrapolation method (D'Auria et al., 2012), and comparison method (Oberkampf et al., 1998), among which the uncertainty propagation method is the most mature and common used one and can be implemented easily. The input uncertainty can be propagated through the analysis code and then be transferred to the predicted results. With this method, samples representing PDFs of input parameters were used as the inputs of the code to obtain the PDFs of results. As for the samples, they can be obtained by the random sampling method or the deterministic sampling methods. The random sampling method is more reliable but needs numerous samples; while the deterministic one has not been fully validated and is far from mature (Bestion et al., 2016).

In the past few decades, most uncertainty analyses associated to system codes were performed based on the random sampling methods (Espinosa-Paredes et al., 2010), or the random methods with surrogate models, such as response surface method (Prošek et al., 2016) and artificial neural network (Secchi et al., 2008). Issues related to random sampling method have been investigated thoroughly, such as the determining of sample size (Wilks,

TABLE 1 | Statistical parameters for the measured input variables.


1942) or the improvements on the sampling methods (Helton and Davis, 2003; Iman, 2008). When it comes to the CFD applications, the deterministic methods with fewer sample requirements were employed to obtain the input parameters due to the increase of computational cost compared to the system code applications. Hessling (2013) proposed an efficient deterministic method for sampling. However, the accuracy and reliability of deterministic method in CFD applications were not fully validated and far from mature (Thiem and Schäfer, 2014). Thus, in this work, we employed the random sampling method to get the input parameters for uncertainty analysis in spite of the huge computational efforts. The LHS method, an efficient sampling methodology developed based on the traditional random sampling and domain decomposition technology (Helton and Davis, 2003), was used instead of the traditional Monte Carlo method to improve the sampling efficiency.

#### Implementation of the Uncertainty Analysis

The CFD simulation on subcooled boiling flow was carried out using the commercial code FLUENT. Eulerian–Eulerian twofluid model along with the interphase actions and wall boiling models used in FLUENT has been validated in our previous work (Zhang et al., 2015c). Besides, the sensitivities of the turbulence model, the mesh size and the near wall turbulence treatment have been studied (Krepper and Rzehak, 2012; Zhang et al., 2015b). According to previous study on turbulence models (Zhang et al., 2015b), the subcooled boiling simulations were carried out with the Realizable k-ε turbulence model and Enhanced wall function based on the grid with near wall cell Y-plus ranging from 20.2

to 33.1. Besides, the inlet boundary was set at the location 10 times upstream of the test section inlet with a fully developed velocity profile to avoid the effects of inlet velocity profile, since there was a long inducer prior to the inlet of heated section to ensure a fully developed flow in the experiment (Bartolomei and Chanturiya, 1967). The geometry and boundary conditions of the computation domain are illustrated in **Figure 1**.

The DAKOTA code (Adams et al., 2017), an uncertainty analysis code developed by the Sandia National Laboratory, was employed to generate the samples with LHS method, to drive the FLUENT code and to post-process the results. An interface was developed to allow the connection between DAKOTA and FLUENT codes by using Python script and the Internet Inter-ORB Protocol (IIOP) interface. The data flow chart is shown in **Figure 2**. A python subroutine is called by the DAKOTA code as an interface between DAKOTA code and the IIOP interface, which is an addon to allow the usage of FLUENT code as a server session. The physical properties of water may be effected by the input parameters, such as system pressure. Thus, a User-Defined Function (UDF) is compiled into standard FLUENT solver to call the IAPWS dynamics link library for liquid and vapor properties.

As mentioned above, the uncertainty propagation with random sampling method always needs numerous samples and runs to capture the statistical characteristics of inputs and outputs. Then the performance of the code could be evaluated. However, how many samples are enough from the view of statistics? Wilks presented an equation to correlate the confidence level, the coverage rate and the size of samples (Wilks, 1942), which has been extensively used in literature. Nevertheless, some other works increase the size of samples until the statistical parameters converge (Bestion et al., 2016). In our work, the Wilks's correlation under the condition of two-sided limits was employed to estimate the preliminary sample size for the inputs and then the sample size was increased until the value of expectation and standard deviation converge.

#### RESULTS AND DISCUSSIONS

Multiple sets of output can be obtained in the uncertainty analysis with error propagation method. Three-dimensional variables should be analyzed to understand the statistical characteristics of the two-phase flow features. Six variables were chosen among the numerous data to present the averaged and local two-phase characteristics, including,


Prior to the simulations, the preliminary size of samples was determined by using the Wilks's correlation (Wilks, 1942) and its variant (Pal and Makai, 2003). For problems with two-sided tolerance region, Wilks's correlations with first-order accuracy and high-order accuracy can be written as Equations (3) and (4), respectively.

$$\beta = 1 - \gamma^N - N(1 - \gamma)\gamma^{(N-1)}\tag{3}$$

$$\beta = \sum\_{j=0}^{N-2n} \binom{N}{j} \gamma^j (1-\gamma)^{N-j} \tag{4}$$

where β , N, n, and γ are the confidence level, sample size, order of accuracy and coverage probability of the confidence interval, respectively. The sample size N can be obtained by solve the correlation with the requirement of 95/95% rule, that is, the coverage probability of confidence interval is no < 95% with a confidence level no < 95%. The minimum sample requirements for first, second and third order accuracy are 93, 153, and 208, respectively. Liu (2014) investigated the uncertainty of velocity profile predicted by a CFD code for forced convection in a cavity with Monte Carlo sample method and found that 153 samples can reach convergence results. However, as pointed out by Bestion et al. (2016) it may not reach the convergence values for the expectation and standard deviation when the samples size determined by this Wilks's equation was used. Thus, we increased the sample size until the statistical characteristics converged. The convergence history is shown in **Figure 3**. As can be seen, the convergence results can be obtained when the sample size is as large as 740, which is the sample size calculated by Wilks's equation with fourteenth order of accuracy.

TABLE 2 | Statistical characteristics for the outputs.


averaged fluid temperature.

After carrying out 740 runs of calculation, we can obtain the resultant two-phase parameters with consideration of the input uncertainty, including the averaged parameters and localized parameters as shown in **Table 2**. First of all, the statistical characteristics of the calculated vapor fraction, fluid temperature, and wall temperature were analyzed and compared with the experiment data (Bartolomei and Chanturiya, 1967; Wang et al., 2016) 1 . The averaged value, 95% confidence interval and the +3σ interval are presented in **Figure 4**. It can be noted that the predicted cross-section averaged VOF and fluid temperature and their error bands agree quite well with the experiment data. However, there are slight deviations between the calculated and measured wall temperature, which might because the wall heat partition models used to simulate the boiling process at wall. Besides, the distribution of wall temperature is rather concentrated and its uncertainty is small. The effects of input uncertainties on the vapor fraction increase with the vapor fraction while wall temperature shows fewer impacts caused by the input uncertainty with increasing the vapor fraction.

**Figure 5** presents the histograms of above mentioned six representative outputs. As can be seen, the statistical distributions of these parameters follow the rule of high probability density in the middle of the interval and low value at the two ends of the interval, which is consistent with the probability density of inputs. However, the distribution is diverse from the normal distribution to some degree. This means, the statistical characteristics of the inputs change during the propagation of uncertainty in the FLUENT code.

To quantify the PDF characteristics of these six variables, the statistical characteristics of these parameters were given in **Figure 5**, including expectations, standard deviations, 95% confidence interval, error bands, relative error bands, and the p-value for normality test. Here we also use the +3σ as the error band. As can be noted, the ranges of deviations for the outputs vary from each other. The uncertainties of the outlet VOF, NVG position, ratio of maximum and averaged VOFs at outlet and wall maximum VOF are much larger than the input uncertainties, while the uncertainties of the pressure drop and wall averaged temperature are less than the input uncertainties. That means, the outlet VOF, NVG position, ratio of maximum, and averaged VOFs at outlet and wall maximum VOF are much more sensitive to the input uncertainty than the pressure drop and wall averaged temperature under the calculated condition. Moreover, the p-values resulted from normal distribution test show that the pressure drop, averaged outlet VOF, averaged wall temperature, ratio of maximum and averaged VOFs at outlet and wall maximum VOF follow the normal distribution well while the position of NVG doesn't.

The scatter plots between the inputs and outputs are given in **Figure 6** to illustrate the relationships between the boundary conditions and the interested parameters. As can be found, the mass flux has a significant linear effect on the pressure drop, that is, the pressure drop increases with mass flux linearly. The averaged outlet VOF, averaged wall temperature and maximum VOF at wall show slight decrease tendency with

<sup>1</sup>The experiment was performed by Bartolomei and Chanturiya (1967). However, this work only gave the error band for the input parameters, including the inlet temperature, heat flux, pressure, and mass flux. The error bands of for the fluid temperature, vapor fraction and wall temperature were not presented in this work. The error bands for these parameters were obtained from publication of Wang et al. (2016).

increasing mass flux. The NVG position and averaged outlet VOF increase with the mass flux, obviously. The temperature of fluid at inlet has dramatic impacts on the averaged outlet VOF, NVG position, and VOF ratio, while has little impacts on the pressure drop, averaged wall temperature and maximum wall temperature, which is because of the dependency of inlet enthalpy on the inlet temperature. The wall heat flux shows strong correlations with all the outputs except the pressure drop. The averaged outlet VOF, averaged wall temperature and maximum wall VOF increase with increasing the wall heat flux, while position of NVG and VOF ratio decrease. This is caused by the effects of wall heat flux on the nucleation process and the generation rate of vapor phase. Besides, the averaged wall temperature depends on the system pressure, due

to the relationship between saturated temperature and system pressure.

In the field of statistical analysis, covariance can be used to descript the dependency or correlation between two sets of samples. However, the values of covariance are associated with the absolute values of samples and cannot be used for the cross-comparison of the dependency between different sets of samples. The conception of correlation coefficient with values varying from zero to one was proposed to quantify the dependency of samples. The Pearson correlation coefficient was defined to quantify the linear dependency among samples, while the Spearman correlation coefficient was used to descript the

FIGURE 7 | Correlation coefficients between boundary conditions and outputs. (A) Correlation coefficients between pressure drop and boundary conditions. (B) Correlation coefficients between averaged outlet VOF and boundary conditions. (C) Correlation coefficients between averaged wall temperature and boundary conditions. (D) Correlation coefficients between NVG position and boundary conditions. (E) Correlation coefficients between VOF ratio and boundary conditions. (F) Correlation coefficients between maximum wall VOF and boundary conditions.

more complicated dependency among samples, no matter linear or not (Hauke and Kossowski, 2011). Pearson and Spearman correlation coefficients are given in **Figure 7** to quantify the relationship between boundary conditions and outputs. As can be seen, the effects of mass flux and heat flux on the subcooled boiling flow are much large than these of inlet temperature and system pressure. The absolute values of correlation coefficients between heat flux and most outputs, including averaged outlet VOF, averaged wall temperature, NVG position, and ratio of maximum and averaged VOFs at outlet, are larger than 0.5, which means there are strong dependencies among these variables. Besides, the Spearman and Pearson correlation coefficients are almost equal, which implies the dependencies between these parameters are linear. Besides, there is no dependency between the system pressure and pressure drop or maximum wall VOF since the absolute values of correlation coefficients are < 0.1.

# CONCLUSIONS

In this work, the subcooled boiling flow was modeled by a CFD code (ANSYS FLUENT) with considering the uncertainty of boundary conditions. The uncertainty was analyzed by using the uncertainty propagation method with Latin Hypercube Sampling method. Seven hundred and forty sets of samples were employed as inputs to achieve the convergence results from the point of view of statistical analysis. Effects of uncertainties from boundary conditions on the two-phase characteristics were obtained. Following conclusions can be drawn from the results.


# AUTHOR CONTRIBUTIONS

XZ and RZ: carried out the baseline case study on subcooled boiling, the sampling of boundary conditions, the results analysis; XL: wrote the python scripts for the uncertainty analysis and carried out the uncertainty analysis; TC: analyzed the results and wrote the paper.

# ACKNOWLEDGMENTS

This work is supported by National Natural Science Foundation of China (No. 11705035).

# REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Zhang, Zhang, Lv and Cong. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Reducing the Flow Distribution in Central-Type Compact Parallel Flow Heat Exchangers Having Optimized Tube Structure

#### Jian Zhou\*, Ming Ding, Zhaoming Meng, Yinxing Zhang and Zhongning Sun\*

*Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Heilongjiang, China*

Compact parallel manifolds are widely applied in flow heat exchangers. However, the flow maldistribution exists in the parallel manifolds when the flow is distributed in the header, because of the pressure maldistribution. For a heat exchanger, the uniformity of the flow distribution greatly influence the performance of the heat exchanger. It is significantly important to modify the design of the heat exchanger for a uniform flow distribution. In present study, to reduce the maldistribution in the manifolds, the numerical analysis has been made and a method by varying the insert length of tubes inside the header has been pointed out to solve this problem of flow maldistribution. To illustrate the reliability of this method to reduce the maldistribution, three base cases with different header diameters have been applied to be solved with the method. The results indicates that, by using the method, the maldistribution for three base cases can be reduced by 82% for the case1, 72% for the case2 and 68% for the case3. However, the insert length of tubes inside the header increases more pressure loss. The pressure drop for base cases are respectively increased by 2.83, 4.83, and 6.46%. This method is proved to be effective under all flow rates.

#### Edited by:

*Jun Wang, University of Wisconsin-Madison, United States*

#### Reviewed by:

*Xianping Zhong, Tsinghua University, China Claudio Tenreiro, University of Talca, Chile*

#### \*Correspondence:

*Jian Zhou zhoujian@hrbeu.edu.cn Zhongning Sun sunzhongning@hrbeu.edu.cn*

#### Specialty section:

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

Received: *13 February 2018* Accepted: *13 April 2018* Published: *07 June 2018*

#### Citation:

*Zhou J, Ding M, Meng Z, Zhang Y and Sun Z (2018) Reducing the Flow Distribution in Central-Type Compact Parallel Flow Heat Exchangers Having Optimized Tube Structure. Front. Energy Res. 6:35. doi: 10.3389/fenrg.2018.00035* Keywords: flow maldistribution, pressure maldistribution, manifolds, central-type heat exchanger, uniformity

# INTRODUCTION

Central-type compact parallel flow heat exchangers are used for removing the heat and ensure the safety of the system. These manifolds are widely applied in cooling systems such as passive containment cooling system of the pressure water reactor, spargers, electronic cooling equipment, and passive core cooling system.

However, the mass flow maldistribution occurred in the heat exchanger greatly reduces the heat transfer performance. Besides, the variations of heat transfer performance caused by maldistribution result in deviations from desired design performance, which may cause that the heat exchanger could not meet the need of the design performance in the practical application.

The flow maldistribution existed in the heat exchanger is caused by the pressure distribution in the header, therefore, the maldistribution in parallel manifolds is an inherent characteristic

**Abbreviations:** Dh, Header diameter; D<sup>t</sup> , Tube diameter; Dpt, The pitch between the tubes at the center; SE, Evaluation parameter of flow maldistribution; k, turbulent kinetic energy; ε, turbulent energy dissipation rate; ρ, density of the working fluid; u, velocity; σ<sup>k</sup> , σ<sup>ε</sup> , turbulent constants; µ<sup>t</sup> , turbulent dynamic viscosity.

of the heat exchanger. The way to flatten the flow distribution could be flattening the pressure distribution.

For instance, when dealing with tubes, Said et al. (2014) successfully reduced the maldistribution in a central-type heat exchanger with the method of applying orifice approach and nozzle approach in tubes respectively. Huang and Wang (2013) have done the investigation of design of uniform tube flow rates in a Z-type heat exchanger. Levenberg-Marquardt Method (LMM) was applied to decide the estimated optimal tube diameters and entrance length of the inlet. And the results showed that the flow distribution was nearly uniform. Tong et al. (2007) have applied gate-valve-like obstructions in manifolds in order to tailor the resistance of an individual channel to achieve a uniform distribution. And results showed that this method worked well, for example, in one of the case studies, the original flow imbalance of the untailored manifold system exceeded 100% and was reduced to less than 10% of the flow imbalance in one cycle of the method. In another study of Tong et al. (2009), the cross-sectional areas of the outflow channels have been modified in linear tapering or non-linear tapering, and results showed that the flow distribution become uniform after solutions were carried out.

As for the promotion of the header design, in the study by Wang et al. (2011b), five modified headers have been investigated, the results showed that the baffle tube performed best and it was applicable for all the flow rates. Wang et al. (2014) have proposed a numerical model of plate-fin heat exchanger to investigate the hydrodynamic characteristics with the method of applying porous media approach in the header. The results showed that the flow distribution was promoted and a correlation among maldistribution, pressure drop, and Reynolds number was derived based on the simulation. Similarly, a uniformly perforated grid has been used in the header of a cross flow heat exchanger in the study of Lalota et al. (1999) They reported that the grid was helpful to improve the fluid distribution.

Except for the promotion of structure for heat exchanger configuration, the influence of the structure parameters on the flow distribution has also been investigated. Gandhi et al. (2012) have done a lot of investigations on the influence made by configuration parameters on the flow distribution. According to their study, the tube diameter, number of tubes and the location of the inlet and outlet pipe are the key configuration parameters affecting the flow distribution. They also have done the experiments to validate simulation results. Bajura and Jones (1976) suggested that the area ratio(AR) < 1 for the parameter design. And Tong et al. (2009) found that enlarging the header diameter for increasing the flow ratio will help improve the flow distribution. Wang et al. (2011a) have done the investigations on the compact parallel flow heat exchangers. According to their studies, the uneven flow distribution is related to following parameters: (1) inlet flow condition, (2) tube diameter, (3) header diamter, (4) area ratio, (5) flow directions. And the effect of gravity is insignificant.

There is also some analytical method to investigate the flow distribution in the heat exchanger. Wang and Wang (2015) have developed a discrete model for flow uniformity and pressure drop in U-type arrangement as well as Z-type arrangement. The results show that the flow distribution in U-type arrangement is generally better than that of the Z-type in cases for large mass flow rate and small pressure drop coefficients. Tong et al. (Sparrow et al., 2007) set forth a quasi-analytical method which is fully validated for determining the flow distribution in multi-inlet collection manifold. And this method can help design a manifold in which the flow distribution is nearly uniform.

In the past studies, researchers focus on the optimization of the flow distribution in compact parallel flow heat exchangers, and a lot of valuable results have been obtained. However, little work has been done on tubes. Besides, it could be more easier and convenient to even the flow distribution with the modification with tubes. In this study, the attention are paid on the modification of tubes. For a conventional configuration of heat exchanger, the header and heat transfer tubes are seamlessly connected as shown in **Figure 1A**. However, this conventional design does not perform well because of the pressure maldistribution in the header. The modification of the tube in this study is applied to decrease the pressure malsitribution at the inlet of the tubes as shown in **Figure 1B**. The tubes are not connected to the header surface, but inserted inside the header. Besides, the insert length of tubes are not equal thus evening the pressure distribution at the inlet of tubes for a better flow distribution.

# PROBLEM DESCRIPTION

The main objective of this study is to search for a better flow distribution in a central-type compact parallel flow heat exchanger. Therefore, three basic configurations of exchangers with different headers as well as different flow distributions are applied of analysis. The three-dimensional model and schematic diagram of the configuration for this study are shown in **Figure 2**. Configuration models consist of two headers called dividing header and combining header respectively and 14 C-tubes are connected to the headers. The working fluid flows into dividing header and then is divided into 14 branches of flow through 14 C-tubes combining in the combining header at last.

# SOLVING METHOD

Three optimization based on three basic configurations with different size of headers are solved respectively. Three header diameters are 60, 80, and 100mm respectively denoted case1, case2 and case3. The objective of the problem is to vary the insert length of the tube inside the header for a better flow distribution as shown in **Figure 3**. Therefore, three different basic configurations with different diameters of headers are solved to determine the basic flow distribution. Based on the basic configuration, the insert of tubes are introduced to reduce the flow maldistribution. The problem is iteratively solved after each adjustment of insert length to produce a more even flow distribution. The iterations are finished until the flow distribution is basically no longer better by changing the insert length of tubes.

# SOLVING PROCESS

The CAD in the Star-ccm+ is used to create the heat exchanger model. Two different kinds of mesh elements are applied for the proper meshing of the whole geometry. The polyhedral mesh is applied to mesh headers, inlet, outlet and part of the Ctubes, and the rest of the C-tubes is meshed with generalized cylinder mesh. Besides, the prism layer model is used for all of the geometry walls to solve the boundary conditions as shown in **Figure 4**.

Star-ccm+ 10.04 software is applied to complete the simulation and constraints described by the boundary of a given system have been applied to solve the control equations for mass and momentum. The resulting partial differential equations together with a high turbulence model are solved numerically.

For the boundary conditions, velocity-inlet is chosen for the inlet, pressure-outlet selected for the outlet is set as zero gauge pressure and the walls are set as no slip condition and rough. The k–ε turbulent model is chosen as the turbulence model.When all of the residuals are less than 1 × 10−<sup>3</sup> , solutions are considered to be completely convergent.

In this study, water has been selected as the working fluid. Three assumptions have been made for the system as following listed.


To evaluate the flow distribution, dimensionless parameter SE has been utilized.

$$SE = \frac{\sqrt{\frac{\sum\_{i=1}^{N} (m\_i - m\_{i\nu})^2}{N}}}{M} \tag{1}$$

Where the mi and mav represent the mass flow rate through the ith tube and the average mass flow rate respectively. And the M represents the total flow rate.

# GOVERNING EQUATIONS

The steady-state continuity equation is written as

$$\frac{\partial \mu\_i}{\partial \varkappa\_i} = 0 \tag{2}$$

The steady-state momentum conservation equation is written as

$$
\rho \mu\_{\dot{j}} \frac{\partial \mu\_{\dot{i}}}{\partial \varkappa\_{\dot{j}}} = -\frac{\partial p}{\partial \varkappa\_{\dot{i}}} + \frac{\partial}{\partial \varkappa\_{\dot{j}}} [\mu\_{\ell} (\frac{\partial \mu\_{\dot{i}}}{\partial \varkappa\_{\dot{j}}} + \frac{\partial \mu\_{\dot{j}}}{\partial \varkappa\_{\dot{i}}})] \tag{3}
$$

The steady-state transport equation for k is written as

$$
\rho \mu\_j \frac{\partial k}{\partial \mathbf{x}\_j} = \frac{\partial}{\partial \mathbf{x}\_j} (\frac{\mu\_t}{\sigma\_k} \frac{\partial k}{\partial \mathbf{x}\_j}) + \mu\_t (\frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i}) \frac{\partial u\_i}{\partial \mathbf{x}\_j} - \rho \varepsilon \tag{4}
$$

The steady-state transport equation for ε is written as

$$\rho \mu\_{j} \frac{\partial \varepsilon}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} (\frac{\mu}{\sigma\_{\varepsilon}} \frac{\partial \varepsilon}{\partial \mathbf{x}\_{j}}) + \mathbf{C}\_{1} \mu\_{t} \frac{\varepsilon}{k} (\frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial u\_{j}}{\partial \mathbf{x}\_{i}}) \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} - \mathbf{C}\_{2} \rho \frac{\varepsilon^{2}}{k} \tag{5}$$

Where the k and ε represent the turbulent kinetic energy and turbulent energy dissipation rate, respectively. ρ is the density of the working fluid, u is the velocity, turbulent constants σ<sup>k</sup> = 1.0 and σ<sup>ε</sup> = 1.3 are used and µ<sup>t</sup> means turbulent dynamic viscosity. Empirical constants C<sup>1</sup> = 1.44, C<sup>2</sup> = 1.92 are used.

#### GRID INDEPENDENCE

As for the grid independence, four different number of grids are applied on the geometry. The number of grids are 2,365,663; 2,941,239; 4,539,910; 7,925,274 and denoted respectively as mesh1, mesh2, mesh3, and mesh4. Because the two most important objectives in this study are flow distribution and pressure drop, the evaluation criteria of the performance of the grid is based on the value of SE and the pressure drop in simulation results. The mass flow rate through 14 tubes is shown in the **Figure 5**. It can be seen that, the significant change of results stops when the number of grid elements increases to the mesh2. Therefore, mesh2 is chosen as the independent mesh size.

# MODEL VALIDATION

The experimental data from Gandhi et al. (2012) was used to validate the numerical results in present study. Gandhi et al. (2012) used 10 tubes with a cylindrical header and an inlet nozzle whereas in the present study, 14 tubes, two headers, an inlet nozzle and a outlet nozzle are used. However, the rest of the geometry is similar. Besides, in both studies, the water was used as the working fluid. Except for experiments, Gandhi et al. (2012) also has done the numerical work based on commercial software Fluent (version 6.2.16) and three different inlet velocities have been applied. Both of the numerical and experimental results done by Gandhi et al. (2012) as well as simulation results in present study are shown in **Figure 6**. It can be seen that the simulation for the present study suits well, and the results of the further work made by the simulation in this study are credible.

# RESULTS AND DISCUSSION

### Numerical Analysis

Basic equations for header and tubes have been separately formulated. In equations, the velocity and static pressure all represent average values. The configuration of dividing header and connected tubes are shown in **Figure 7**. For the control volume of the header section at the inlet of the ith tube, it is enclosed up with a dotted line as shown in **Figure 7**.

The Bernoulli equation for the flow in the tube is written as

$$\frac{1}{2}[P\_{lh}(i) + P\_{hr}(i)] = P\_{out}(i) + [1 + C\_t + C\_b + C\_{td} + \lambda\_t(i)\frac{L\_t}{D\_t}]$$

$$\frac{1}{2}\rho\_l V\_t^2(i)i = 1,2...n\tag{6}$$

Where C<sup>t</sup> represents local friction coefficient of tubes, for example, that caused by vortex flow at the tube inlet. C<sup>b</sup> represents bend loss coefficient. And Ctd represents inlet friction coefficient of tubes.

The continuity equation is written as

$$
\rho\_h A\_h(\mathbf{i})[V\_{hl}(\mathbf{i}) - V\_{hr}(\mathbf{i})] = \rho\_l A\_l V\_l(\mathbf{i})\mathbf{i} = 1, 2...n \tag{7}
$$

Where Ah(i)and A<sup>t</sup> represents effective flow cross section area of the header in front of the ithtube and the tube cross section area respectively. Besides, ρ<sup>h</sup> equals to ρ<sup>t</sup> , considering the assumption of incompressible flow.

The Equations (6, 7) show the influential parameters for the velocity in tubes. The way to achieve a uniform distribution lays on the adjustment of these influential parameters. For a certain heat exchanger with fixed structure parameters. The most effective and direct method is to vary the Ah(i) for each tube as shown in Equation (7). In addition, the variation of Ah(i) will unequally change friction coefficients such as C<sup>t</sup> and Ctd as well as pressure distribution in the header, which means that the influence of the insert length for each tube is not equal.

For compact parallel flow heat exchangers, the dynamic pressure is higher near the center of header, while the static pressure is higher near the edge, because of the influence of pressure recovery. Therefore, the sensitivity of variable insert is higher at the header center. And this will help save time for the iteration work.

The insert length of each tube for three cases with the best solution is listed in **Table 1**.

# Mass Flow Rate Distribution

The mass flow rate for three base cases with different header diameters is shown in **Figure 8** respectively. The flow patterns

of flow distribution on condition of without solution, best case solution and iterative solutions can also be seen in the **Figure 8**, For the flow distribution of the case1, the mass flow rate is highest at center tubes, and continues to drop until at tube 3 and tube 12 and then increases slightly from tube 3 and tube 12 to the extreme tubes of the header as shown in **Figure 8A**. Because of the small size of the header, the velocity of the water in the header is high, leading to a high pressure recovery in tubes at the center due to the reduction in velocity in the direction of the mainstream

TABLE 1 | The insert length of each tube for three cases with the best solution.


because of the flow into tubes. Therefore, the key to achieve a better flow distribution is to decrease the pressure recovery at the inlet of tubes which high mass flow rate of water flows in. The original process is to increase the extension of the tubes near the inlet inside the header. The extension of the tube inside the header leads to the decrease of the effective flow area above the tube, thus decreasing velocity and the mass flow rate in the tube. Iterative adjustments of variation of the extension for tubes are laid on the analysis discussed above, increasing the extension of tubes where high mass flow is and decreasing the extension of

tubes where low mass flow enters. Iterative solutions performed successfully and the best solution are shown in **Figure 8A**. By comparing the best solution and the base case which is without solution, it can be seen that the mass flow rate in tubes ranging from tube 5 to tube 10 which was high becomes lower and the mass flow rate in the other tubes becomes higher, thus reducing the maldistribution in the header.

For the mass flow rate distribution in the tubes of the case2, the maldistribution is better than the case1, the reason is that the bigger cross section area of the header causes the lower velocity of the flow thus leading to a lower pressure recovery and less mass flow rate in center tubes, therefore, the mass flow rate distribution in the case2 is better than that in the case1. In this case, the key to solve this problem is to focus on the tubes at the center and tubes close to extreme ends of the header, in which the mass flow rate is much higher than the average value of mass flow rate in tubes. Based on the analysis discussed before, iterative adjustments of variation of the extension for tubes have been made and the best solution as well as iterative solutions performed successfully is shown in **Figure 8B**. It can be seen that the mass flow rate in center tubes and tubes close to extreme ends of the header has been successfully decreased making the flow distribution more even.

For the case3, the bigger cross section area of the header than that in case2 continues to decrease the pressure recovery in center tubes leading to even worse flow distribution in the header causing that large amount of mass flow rate goes through tubes which are near extreme edges of the header. For this case, the promotion should be applied to tubes at the extreme edge of the header. The base case without solution, the best solution and iterative solutions performed successfully are shown in the **Figure 8C**. It can be seen that the flow distribution has been successfully promoted, showing the effectiveness of the method applied in this study.

The method applied in this study all performed well when dealing with these three base cases. The dimensionless parameter SE has been applied to evaluate the maldistribution of the flow distribution in the header. The best solution reduces the SE by 81, 69, and 69% for case1, case2, and case3 respectively. The main concept utilized in the method of reducing the maldistribution is to adjust the effective flow area and the pressure at the tube inlet thus achieving the goal of balancing the mass flow rate through tubes.

#### Pressure Distribution

The **Figure 9** shows the contour diagrams for the static pressure within the header and tubes for all cases under investigation. For all the cases as it can be seen from **Figure 9**. There is almost a significant pressure drop from the header to the tube for the flow resistance at the interface of the header and tubes. The flow recovery happens at the tube inlet increasing the static pressure in the header. And it can be found that the static pressure increase from the center of the header along the way to the edge of the header. For the cases with the best solution, the pressure gradient in the header and tubes is both smaller than that in the cases without solution. It is because that the different insert length of the tube inside the header promote the pressure distribution. For instance, the pressure at the edge of the header is bigger than other position. The longer extension of the tubes near the extreme edge increase more pressure drop from the header to the tube, thus making the pressure distribution in tubes more uniform.

#### Velocity Distribution

The **Figure 10** shows the velocity distribution within the header and tubes for all cases under study. It can be seen in **Figure 10A** that the flow velocity is bigger inside the tubes near the center of the header while smaller in tubes near edges of the header. With the solution, the extension of tubes near the center decrease the velocity inside these tubes near the center, making the velocity gradient in the header and tubes more uniform as shown in **Figure 10B**. For the case2 without solution, the velocity in the header near the center and the edge is much higher than that in the middle as shown in **Figure 10C**. After the solution is applied, the velocity becomes more uniform as shown in **Figure 10D**. For the case3 with bigger header diameter, the mass flow concentrates in the tubes near the edge, because of low pressure recovery in tubes at the center, pushing more mass flow into tubes at the edge, as shown in **Figure 10E**. The way should be to increase the extension of tubes near the edge in order to decrease the Ah(i), making the velocity distribution more uniform as shown in **Figure 10F**.

for Case 3.

TABLE 2 | The pressure loss for cases from the heat exchanger.


# Pressure Drop Through Heat Exchangers

The overall pressure drop through three cases either without or with the best solution are listed in **Table 2**. Because of the insert length of tubes inside the header, the frictional coefficients increase, bringing more pressure loss. For the best solutions, the pressure drops have been increased by 2.83, 4.83, and 6.46% for case1, case2, and case3 respectively, comparing to the cases without

solutions. For the case with smaller header diameter, the velocity in the header is larger. Therefore, the sensitivity of variable insert length is higher and less length of insert length for tubes is needed, which decrease the frictional coefficients brought by the extension of tubes.

# The Application of the Method for Different Inlet Mass Flow Rate

Considering the effect of the mass flow rate on the method in this study, the flow distribution in the case1 without the solution and the case1 with the best solution are numerically calculated under different inlet mass flow rate as shown in **Figure 11**. As the inlet mass flow rate increase, the maldistribution slightly increase for all cases. Besides, the insert method significantly even the flow distribution under all inlet flow rate, showing that the method is effective and applicable for all inlet flow rates.

# REFERENCES


# CONCLUSION

A method by varying insert length of tubes has been applied to reduce the flow maldistribution in manifolds of centraltype compact parallel flow heat exchangers. The key point of this method is to vary the effective flow area for each tube, adjusting mass flow rate through tubes for a more uniform flow distribution. Three base cases with different header diameters and flow distribution patterns have been used to compare with cases with the best solutions in terms of maldistribution and pressure drop. The results show that the maldistribution in the heat exchanger could be reduced by 81% at most and 69% at least. However, the pressure drop in the heat exchanger is increased by 6.46% at most and 2.83% at least when the header diameter is the largest and smallest respectively.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Huang, C. H., and Wang, C. H. (2013). The design of uniform tube flow rates for Z -type compact parallel flow heat exchangers. Int. J. Heat Mass Transf. 57, 608–622. doi: 10.1016/j.ijheatmasstransfer.2012.10.058

Lalota, S., Florentb, P., Langc, S. K. A., and Berglesc, E. (1999). Flow maldistribution in heat exchangers. Appl. Therm. Eng. 19, 847–863.

Said, S. A. M., Mansour, R. B., Habib, M. A., and Siddiqui, M. U. (2014). Reducing the flow mal-distribution in a heat exchanger. Comput. Fluids 107, 1–10. doi: 10.1016/j.compfluid.2014.09.012


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Zhou, Ding, Meng, Zhang and Sun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Numerical Simulation of Flow and Heat Transfer Characteristics of CuO-Water Nanofluids in a Flat Tube

Luchao She\* and Guangming Fan\*

Fundamental Science on Nuclear and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China

In this paper, water and CuO-water nanofluids are selected as working fluids. The characteristics of fluid flow and heat transfer in a heat transfer tube are analyzed by means of numerical simulation. The effect of nanofluid concentration on the flow and heat transfer characteristics of the flat tube under the condition with a Reynolds number of 6,000∼10,000 is analyzed and compared with that of the circular tube. The results show that the addition of nanoparticles in the nanofluid has little effect on the velocity and temperature distribution in the flow channel. With the increase of nanofluid concentration, the heat transfer coefficient and flow resistance of the fluid increased. And the comprehensive evaluation factor was 2.51∼3.29, which had the effect of enhancing heat exchange. There is a similar trend in the same concentration with different Reynolds numbers. Compared with the circular tube, the nanofluid has a better ability to enhance heat transfer when it flows in the flat tube. Under conditions with the same fluid concentration, the effect of the application of the nanofluids on the

Germany

#### Reviewed by:

Edited by: Uwe Schröder,

Claudio Tenreiro, University of Talca, Chile Khalil Ur Rahman, Pakistan Nuclear Regulatory Authority (PNRA), Pakistan

Technische Universitat Braunschweig,

#### \*Correspondence:

Guangming Fan fanguangming007@hotmail.com Luchao She sheluchaowy@163.com

#### Specialty section:

This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research

Received: 31 December 2017 Accepted: 05 June 2018 Published: 21 June 2018

#### Citation:

She L and Fan G (2018) Numerical Simulation of Flow and Heat Transfer Characteristics of CuO-Water Nanofluids in a Flat Tube. Front. Energy Res. 6:57. doi: 10.3389/fenrg.2018.00057 flat tube is better.

Keywords: nanofluids, flat tube, numerical simulation, flow characteristics, heat transfer enhancement

# INTRODUCTION

Nanoparticles are particles in the 1–100 nm, between the microscopic system and macroscopic system, and consist of a group of few atoms or molecules (Choi, 1995). Conventional fluids have poor thermal conductivity. An effective way to improve the heat transfer performance of conventional fluids is to suspend fine solids, such as metals and non-metallic particles. Nanofluids are fluids that add nanoparticles with diameters <100 nm in a conventional fluid (Ganvir et al., 2016). The main objective of the application of nanofluids is to achieve the highest thermal conductivity at the lowest possible concentration of nanoparticles. Due to the molecular chain behavior of nanofluids, the suitable dispersed nanoparticles in the base liquid have the advantages of high thermal conductivity, microchannel cooling, no blockage, reducing erosion probability, and increasing pump power (Solangi et al., 2015). Due to the effect of their small size, the coefficient of resistance is not increased significantly as the coefficient of fluid heat transfer increases. And this makes the nanofluid heat transfer technology become a hot research topic in the field of thermal science in recent years. The application of nanofluid in shell-and-tube heat exchanger can increase the heat transfer performance of the tube, without causing large resistance loss, wear, blockage and other undesirable consequences (Eastman and Choi, 2001).This characteristic endows the nanofluids with a wide application prospect in the field of enhanced heat transfer and a

potential application in nuclear engineering. Heat exchanger elements are important parts of heat exchangers, and the main form of them is a variety of heat transfer tubes. In recent years, many new heat transfer tubes have been developed to meet the requirements of high efficiency and compact heat exchanger, and the flat tube is one of them. Compared with the circular tube, the flat tube has a higher ratio of surface area and flow area, which can enhance the heat transfer rate and improve the efficiency of the heat exchanger. Moreover, the volume of the flat tube is relatively small, which makes the compact of the heat exchanger greatly improved. Therefore, it is necessary to study the flow and heat transfer performance of the nanofluid in a flat tube.

In recent years, various aspects of nanotechnology have been described as emphasizing applications in energy systems, such as nuclear energy, geothermal energy, and solar energy (Banerjee, 2017). In previous studies, nanofluids have been used as working fluids in the plate heat exchangers (Kumar et al., 2015), micro-channel heat exchangers (Qiang and Yi-Min, 2004; Hwang et al., 2009; Liu and Yu, 2010), Serpentine heat exchanger (Khoshvaght-Aliabadi et al., 2016). However, the application of flat tube heat exchanger is mainly based on experimental research (Li and Xuan, 2002; Wen and Ding, 2004; Chen et al., 2013; Ji et al., 2014), while numerical simulations are relatively rare. It is concluded that under the condition with the same Reynolds number, the convective heat transfer coefficient of the nanofluid in the flat tube is larger than that of the pure liquid, and increases with the increase of the particle volume fraction. Interestingly, the resistance coefficient of the nanofluid is not obviously increased. The addition of nanoparticles under laminar flow conditions has improved the heat transfer and pressure loss of the base fluids in all the flat tubes under various Reynolds numbers and temperatures (Zhao et al., 2016). All these conclusions are obtained under the condition of small channel or laminar flow, and the results are relatively few in the turbulent conditions of ordinary flat tubes. In this paper, the flow and heat transfer characteristics of a flat tube under conditions with different Reynolds number are numerically simulated with CuO-water nanoparticles as working fluids.

#### COMPUTATIONAL MODEL

In order to analyze the comprehensive heat transfer capacity of flat tube, the geometrical model of flat tube and circular tube is established, and its cross-section is shown in **Figure 1**. The model includes two parts of tube wall and fluid. In order to ensure the accuracy of the final simulation result and prevent the effect of entrance on the calculation result, the stable section of the heat transfer tube is set up respectively in the model (**Figure 2**). The geometrical dimensions of the heat transfer tube are shown in **Table 1**.

In this paper, the Workbench software is used to divide the hexahedral structured grids by sweep division. The three-stage fluid parts have the same setting. The minimum mesh size of the source surface is 0.8 mm, and the boundary layer grid is refined, and the solid wall surface is set to 5 nodes 4-layer grid. The final grid form is shown in **Figure 3**.

Grid quality measurement method is skewness, in general, the value <0.7 is acceptable. As shown in **Figure 4**, the maximum skewness number of the flat tube meshes is 0.5, and the circular tube is 0.55, both within the allowable range.

In this paper, the correlation between the heat transfer coefficient and the friction coefficient is also verified. The heat transfer coefficient is derived from analog data, and the calculation formula of friction coefficient is as follows:

$$
\xi = \frac{\Delta P}{\frac{L}{d} \cdot \frac{\rho u^2}{2}} \tag{1}
$$

**Nomenculture**: C<sup>p</sup> Specific heat at constant pressure, J·(kg·K)−<sup>1</sup> ; d equivalent diameter of the heat tube, m; f resistance coefficient; h heat transfer coefficient, W·(m<sup>2</sup> ·K)−<sup>1</sup> ; L test section length, m; Nu Nusselt number; 1Ptubeline pressure drop, Pa; Re Reynolds number; u fluid velocity, m/s

**Greek symbols**: ρ fluid density, kg/m<sup>3</sup> ; β thermal expansion coefficient, K−<sup>1</sup> ; µdynamic viscosity, Pa·s.

#### TABLE 1 | Geometrical parameters of heat transfer tube.


Type: 1P is the tubeline pressure drop; L is the test section length; d is equivalent diameter of the heat transfer tube; ρ is fluid density; u is the fluid Velocity.

As shown in **Figure 5**, when the number of tube meshes reaches 370,000, the average fluid heat transfer coefficient and the tube friction coefficient have no significant change, that is, 370,000 is the minimum number of meshes that can guarantee the accuracy of the simulation calculation. Therefore, in this paper, the circular tubes are calculated numerically by using 370,000 mesh grids.

Same as the circular tube, the number of flat tube grids is 570,000.

In the setting of boundary conditions, the inlet of the heat transfer tube adopts the velocity boundary, and the size is determined by Reynolds number at the temperature of 293 K. The outlet is provided with a pressure exit border of 0 Pa, which is discharged into the vacuum. The tube wall is provided with no slip boundary condition, and the heating wall is provided with constant wall temperature boundary of 400 K. At the same time, the fluid surface of the stable section and the end surface of the tube are adiabatic conditions, and the SST model is chosen.

It is known from the literature (Rao et al., 2003; Wang et al., 2003) that the stability of the fluid is the highest when the particle size of nanoparticles is 20 nm. Therefore, the concentrations of CuO particles in water are 0.05, 0.1, and 0.15 with particle size of 20 nm. The thermal conductivity and other physical parameters of the nano-flow vary with the concentration of nanoparticles, as shown in **Table 2**.

#### SIMULATION RESULTS ANALYSIS

#### Heat Transfer Characteristics Analysis

The change of heat transfer coefficient with Reynolds number in the tube is shown in **Figure 6**. The greater the concentration of the nanoparticles in the fluid, the greater the average heat transfer coefficient of the fluid. In the case with Reynolds number of 6,000, the heat transfer coefficient of nanofluid relative to water increased by 5.24, 9.4, and 12.63%, respectively. This is because the thermal conductivity of metal particles is much higher than that of water. Because of the presence of the boundary layer at the wall, the addition of appropriate amount of metal nanoparticles in water can reduce the heat transfer heat resistance while increasing the thermal conductivity, thus it can significantly improve the heat transfer performance of the fluid near the wall and increase the coefficient of heat exchange. At the same time, the average heat transfer coefficient h of water and the concentration of nanometer fluid in the flat tube is 11.38 ∼ 13.3% larger than that of the circular tube. This is due to the heat transfer intensity of the main heating boundary layer thickness, under certain temperature conditions, the thinner the thickness of the flow boundary layer, the thinner the thickness of the hot boundary layer, and the greater the heat transfer intensity of the fluid. According to Yuan-Xun (2005), the formula for calculating the thickness of flow boundary layer is:

$$\delta = 64.2 \frac{d\_e}{\text{Re}^{\frac{\gamma}{\gamma\_s}}} \tag{2}$$

Therefore, under the same section circumference of the heat transfer tube, the equivalent diameter of the flat tube is smaller, that is, the thickness of the flow boundary layer is thinner and the heat transfer effect is stronger under the same Reynolds number.

**Figure 7** shows the change of heat transfer coefficient of each concentration of nanofluids in the two heat tubes in relation to the increase of water and the Reynolds number. It can be seen from the diagram that the ability of the nanofluid to enhance the heat transfer coefficient is stronger in the circular tube. This is because the boundary layer in the flat tube is thinner than that of the circular tube, and it is known from above that the nanofluid can significantly enhance the heat transfer capability of the boundary layer, so the heat transfer coefficient of the thicker circular tube in the boundary layer increases greatly. This phenomenon is most obvious at Re = 6,000, and the heat transfer coefficient of nanofluid in the flat tube is up to 12.63%, while the maximum of the circular tube increases by 14.57%.

In order to understand the temperature change in the heat transfer tube, the axial and radial temperature distributions have been analyzed, and the temperature distributions of the average temperature along the axial and radial directions of each section of the two heat transfer tubes are given. As shown in **Figure 8**, TABLE 2 | Water and concentration of nanofluids physical parameters (Sun et al., 2015).


the temperature along the axial direction of the two heat transfer tubes is basically the same. However, in the same condition, the average temperature of the flat tube along the axial section is higher than that of the corresponding circular tube, and reaches to the maximum value at the outlet section. The temperature growth rate of the flat tube is also higher than that of the circular tube, and at the outlet section, the fluid concentrations are 5.41, 5.62, 5.74, and 5.8%, respectively.

As shown in **Figure 9**, the trend of temperature variation along the radial direction in the circular tube is basically consistent with that of the flat tube. Though the temperature of the point on the cross-section of the circular tube is lower than that of the corresponding point in the flat tube, the temperature gradient in the near wall of the circular tube is greater. This is because the heat transfer coefficient of the flat tube is larger, its absorption of heat is more so the temperature is higher, but the circular tube boundary layer is thicker, the temperature grows faster, so the temperature gradient of the circular tube is larger.

#### Flow Characteristics Analysis

As shown in **Figure 10**, because the flow area of flat tube is relatively smaller compared with that of the circular tube, the medium velocity is higher and the pressure loss is larger. In addition, the viscosity of the nanofluid increases with concentration, which means the larger the concentration of the nanoparticles, the greater the flow resistance is. Therefore, the pressure drop of circular tube tubelines in each Reynolds number is lower than that of the corresponding flat tubes, and the pressure drop of each tubeline increases with the increase of Reynolds number. But the increase rate in the circular tube is lower than the flat tube. Compared with the circular tube, the pressure drop of water in flat tube increases by 55.25 ∼ 66.91%, while the CuO-water nanofluids pressure drop increases by 72.76 ∼ 79.09%.

Due to the addition of nanoparticles in traditional fluids, their heat transfer capacity is enhanced, and the flow resistance is

FIGURE 8 | Temperature comparison of two tubes along axial section.

increased. Therefore, according to the literature (Qing and Zou, 2007), the comprehensive evaluation factor was introduced to analyze the heat transfer performance of the fluid. Its expression is as follows:

$$\eta = \frac{N\mu/N\mu\_0}{\left(f/f\_0\right)^{0.29}} \left(\frac{d}{d\_0}\right)^{-0.13} \tag{3}$$

Where:Nu, f , d are the Nusselt number, the resistance coefficient and the equivalent diameter of the water and each concentration nanometer fluid in the flat tube respectively; Nu0, f0, d<sup>0</sup> is the Nusselt number, the resistance coefficient and the equivalent diameter of the water and the concentration nanometer fluid in the circular tube respectively.

$$Nu = \frac{hL}{k} \tag{4}$$

Where:h is the average heat transfer coefficient of the fluid; L is the length of the test section; k is the coefficient of thermal conductivity.

$$f = \frac{\Delta P}{\frac{L}{d\_i} \cdot \frac{\rho u^2}{2}}\tag{5}$$

Where: 1P is the pressure drop; d<sup>i</sup> is the equivalent diameter of the heat tube; ρ is the fluid density; u is the average velocity.

As shown in **Table 3**, under the same Reynolds number, the comprehensive evaluation factors decreases with the increase of the concentration of nanoparticles in the fluids, but all >1. This shows that with the increase of nanofluids concentration, the comprehensive heat transfer ability of the flat tube relative to the circular tube decreases gradually, but in the range of this data, any working quality is always better than the circular tube. For the same concentration of fluid, the evaluation factor increases with the increase of Reynolds number, that is, the more general

TABLE 3 | Comprehensive evaluation factors of fluid in flat tube and circular tube at different Reynolds number and concentration.


Reynolds number of nanofluids applied to the flat tube, the better the effect.

# CONCLUSION

In this paper, the numerical simulation method is used to compare the flow and heat transfer characteristics of the water and the nanofluids in the two heat transfer tubes, and the following conclusions are drawn:

(1) The greater the concentration of the nanoparticles in the fluid, the greater the average heat transfer coefficient, and the maximum is reached at 0.15. Compared with the circular tube, the nanofluid has a higher heat transfer coefficient when it flows in the flat tube, but it is better to improve the heat transfer coefficient when it is applied to the tube.

# REFERENCES


### AUTHOR CONTRIBUTIONS

In the course of the completion of this thesis, all authors have substantial contributions to design of the work, GF provided the research direction of the subject, LS determined the research method and specific processes. LS carried out numerical simulation and obtained, collated and analyzed the data, in the process, GF has been instructed and helped. LS wrote the first draft of the paper and corrected it by GF also gave some suggestions. All authors approve the final version to be published. All authors agree to be responsible for all aspects of the work.

# ACKNOWLEDGMENTS

I would like to extend my sincere gratitude to my supervisor for his instructive advice and useful suggestions on my thesis. I am deeply grateful for her help in the completion of this thesis. This work is supported by Key Supported Discipline of Guizhou Provence (Qian Xuewei He Zi ZDXK [2016]24), 2011 Collaborative Innovation Center of Guizhou Province (Qian Jiao he xietongchuangxin zi [2016]02).

wire-bonded flat heat pipe using nanofluids. Int. J. Heat Mass Transf. 67, 548–559. doi: 10.1016/j.ijheatmasstransfer.2013. 08.060


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 She and Fan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Numerical Investigation of the Effect of Grids and Turbulence Models on Critical Heat Flux in a Vertical Pipe

Xiaomeng Dong<sup>1</sup> , Zhijian Zhang<sup>1</sup> \*, Dong Liu<sup>2</sup> , Zhaofei Tian<sup>1</sup> and Guangliang Chen<sup>1</sup>

*<sup>1</sup> College of Nuclear Science and Technology, Harbin Engineering University, Harbin, China, <sup>2</sup> Nuclear Power Institute of China, Chengdu, China*

Numerical simulation has been widely used in nuclear reactor safety analyses to gain insight into key phenomena. The Critical Heat Flux (CHF) is one of the limiting criteria in the design and operation of nuclear reactors. It is a two-phase flow phenomenon, which rapidly decreases the heat transfer performance at the rod surface. This paper presents a numerical simulation of a steady state flow in a vertical pipe to predict the CHF phenomena. The detailed Computational Fluid Dynamic (CFD) modeling methodology was developed using FLUENT. Eulerian two-phase flow model is used to model the flow and heat transfer phenomena. In order to gain the peak wall temperature accurately and stably, the effect of different turbulence models and wall functions are investigated based on different grids. Results show that O type grid should be used for the simulation of CHF phenomenon. Grids with Y+ larger than 70 are recommended for the CHF simulation because of the acceptable results of all the turbulence models while Grids with Y+ lower than 50 should be avoided. To predict the dry-out position accurately in a fine grid, Realizable k-ε model with standard wall function is recommended. These conclusions have some reference significance to better predict the CHF phenomena of vertical pipe. It can also be expanded to rod bundle of Boiling Water Reactor (BWR) by using same pressure condition.

#### Edited by:

*Kaiyi Shi, Liupanshui Normal University, China*

#### Reviewed by:

*Claudio Tenreiro, University of Talca, Chile Khalil Ur Rahman, Pakistan Nuclear Regulatory Authority (PNRA), Pakistan*

#### \*Correspondence:

*Zhijian Zhang zhangzhijian\_heu@hrbeu.edu.cn*

#### Specialty section:

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

Received: *20 April 2018* Accepted: *05 June 2018* Published: *03 July 2018*

#### Citation:

*Dong X, Zhang Z, Liu D, Tian Z and Chen G (2018) Numerical Investigation of the Effect of Grids and Turbulence Models on Critical Heat Flux in a Vertical Pipe. Front. Energy Res. 6:58. doi: 10.3389/fenrg.2018.00058*

Keywords: numerical investigation, Critical heat flux, turbulence models, wall functions, grids distribution

# INTRODUCTION

Critical Heat Flux (CHF) is a two-phase flow phenomenon that is characterized by a heat transfer mechanism change which rapidly decreases the efficiency of the heat transfer performance at the heater surface. When it occurs, heated surfaces are no longer wetted by boiling liquid and the vapor phase start to occupy the heat surface. As a result, the energy is directly transferred from the heat surface to vapor. It results in rapid reduction of the heat removal ability and sharp rise of the vapor temperature, as well as the rod surface temperatures which are important for the nuclear safety. During the design and operation of the nuclear reactors, CHF value should be calculated in advance. But experience and thousands of data points have shown the complexity of CHF phenomenon for different conditions. In the past decades, both experiment and numerical simulation are widely used to predict the CHF value during the design and operation process. During the experimental process, the wall temperature is monitored by use of thermocouples. Along with the increasing of heat power, CHF phenomenon is detected when the temperature of one thermocouple has a sharp rise. Similar to the experiment, this sharp rise of the wall temperature is also a signal of CHF in the numerical simulation. So it is of great significance to predict the wall temperature accurately and stably.

Nowadays, experimental method is widely used in the CHF prediction. But most of them are focused on the vertical pipe. Large scale experiments of rod bundle are few reported or open-access around the world. In the reactor design area, CHF look-up table is the main method to determine the CHF value. But both these two methods will give a large safety margin, and then reduce the power level. Compared to experimental measurements and CHF look-up table, numerical simulation has its own advantages. Thermal-hydraulics system or subchannel codes can efficiently characterize bulk flow behavior, while computational fluid dynamics (CFD) analysis may provide relatively accurate results for the velocity and temperature profiles around the rods. Although the subchannel analysis is more widely used for CHF prediction, CFD method can work as a supplement to provide more accurate results. Furthermore, numerical simulation can reduce the cost significantly when we use different geometries, materials, and test conditions. More accurate simulation techniques using these combined analysis tools can be used to optimize experimental designs and improve the test-analysis design cycle for advanced fuel rods.

Grids and turbulence models are found to be very important in the Computational Fluid Dynamic (CFD) simulation. At the moment, turbulence models used in two phase flow are still same with that of single phase flow. Especially in Eulerian two phase model, the liquid and vapor phases are applied with same turbulence model separately. In single phase flow, the results of Chen et al. (2017a) show different turbulence models will deliver to different secondary flow status. This conclusion is more obvious for isotropic and anisotropic turbulence models. However, the conclusion about applicability of these turbulence models is not consistent for the simulation of two phase flow. Standard k-ε model has been chosen by some researchers for its stability on the two phase calculation. Shin and Chang (2009) has used this model to study the effect of mixing vane on CHF enhancement by use of two-fluid model. It is also used by Shirvan (2016) to study the boiling crisis at high pressures. Shirvan's Results show same trend with the Russian CHF data at the condition close to Departure of Nucleate Boiling (DNB). In addition, a new mechanistic model was developed by Mimouni et al. (2016) to predict CHF in boiling flow of water and R12 refrigerant by use of k-ε models in Neptune-CFD code. Mimouni's results are compared with 150 tests and the mean relative error is equal to 8.3%. Apart from standard k-ε model, Shear Stress Transport (SST) k-ω model is also used by some researchers. Krepper et al. (2007, 2013) and Krepper and Rzehak (2011) predicted the subcooled boiling process in vertical heated tubes with water and R12 based on the SST k-ω turbulence model. Subcooled boiling flow in Internal Combustion (IC) engine cooling passages was studied based on SST k-ω turbulence model (Hua et al., 2015). The simulation results showed that present two-fluid model could get an accurate temperature field for cylinder head.

Although different turbulence models are used for the calculation of two phase flow, few researchers has ever compared the effect of different grids and turbulence models. In general, a grid dependency study should always be performed in high quality CFD simulation. Even though it is hard to obtain fully grid independent results in some cases, it is necessary to discuss the effect of key grid parameters. Thakrar et al. (2016) has investigated the nucleate boiling in rectangular geometries at high pressure by use of CFD code STAR-CCM+. Standard kε model and Reynolds Stress Transport models were compared to investigate the influence of turbulence. Results show the most mechanistic configurations produced remarkably good agreement with measurements of the area-averaged void. But the grid-dependent wall treatments are still needed to develop. Zhang et al. (2015) has analyzed the effect of grids and turbulence models for subcooled boiling flow in a vertical pipe. Her results show that the turbulence models could be chosen appropriately according to the Y+. This paper shows a relatively overall investigation for the options of grids and turbulence models in the subcooled boiling flow.

However, as we can see above, few researchers have ever shown the sensitivity analysis on turbulence models and wall functions for CHF phenomenon. The grid dependency study is also ignored in the CFD simulation of boiling two-phase flow, especially on the CHF phenomenon. It is still not clear that if the turbulence models have the same performance on CHF phenomena as that of subcooled boiling flow, especially on the value and position of highest wall temperature. So based on Eulerian two-phase model, this paper presents numerical simulations of CHF phenomenon by the use of CFD method to investigate the effect of different grids, turbulence models and wall functions. Lack of open CHF experimental data of rod bundle, especially the wall temperature distribution, we use vertical pipe instead. Considering both the flows in rod bundle and vertical pipe are vertical upward flow, the flow and heat characteristics are similar to each other as long as we use same pressure condition. The type of boiling crisis is dryout which can be identified by the distribution of wall temperature and void fraction. Results are focused on the wall temperature distribution to give a reference for the detector of CHF phenomena.

### MATHEMATICAL MODELS

Eulerian two-phase models are used as the basic models to simulate the two phase flow in a vertical pipe. All the interfacial mass, momentum and energy transfer models are based on the interfacial area density model. Furthermore, CHF model will also be shown below.

# Conservation Equation

The conservation equation of Eulerian two phase model includes six equations which are mass, momentum and energy equations for two phases separately. These six equations will be solved for six parameters which include pressure, velocities of liquid and vapor, temperatures of liquid and vapor, volume fraction.

(1) Mass equations

$$\begin{cases} \frac{\partial}{\partial t} \left( \alpha\_l \rho\_l \right) + \nabla \cdot \left( \alpha\_l \rho\_l \nu\_l \right) = \mathbb{S}\_l + m\_{\mathbb{v}l} - m\_{l\mathbb{v}}\\ \frac{\partial}{\partial t} \left( \alpha\_\nu \rho\_\nu \right) + \nabla \cdot \left( \alpha\_\nu \rho\_\nu \nu\_\nu \right) = \mathbb{S}\_\nu + m\_{l\mathbb{v}} - m\_{\mathbb{v}l} \end{cases} \tag{1}$$

#### (2) Momentum equations

$$\begin{cases} \begin{aligned} \frac{\partial}{\partial t} \left( \alpha\_l \rho\_l \nu\_l \right) + \nabla \cdot \left( \alpha\_l \rho\_l \nu\_l \nu\_l \right) &= -\alpha\_l \nabla p + \nabla \cdot \pi\_l \\ &+ \alpha\_l \rho\_l \mathbf{g} + m\_{vl} \nu\_{vl} - m\_{lv} \nu\_{lv} + F\_l \\ \frac{\partial}{\partial t} \left( \alpha\_\nu \rho\_\nu \nu\_\nu \right) + \nabla \cdot \left( \alpha\_\nu \rho\_\nu \nu\_\nu \nu\_\nu \right) &= -\alpha\_\nu \nabla p + \nabla \cdot \pi\_\nu \\ &+ \alpha\_\nu \rho\_\nu \mathbf{g} + m\_{lv} \nu\_{lv} - m\_{vl} \nu\_{vl} + F\_\nu \end{aligned} \end{cases} \tag{2}$$

(3) Energy equations

$$\begin{cases} \begin{aligned} \frac{\partial}{\partial l} \left( \alpha\_{l} \rho\_{l} h\_{l} \right) + \nabla \cdot \left( \alpha\_{l} \rho\_{l} \nu\_{l} h\_{l} \right) &= \alpha\_{l} \frac{\partial \rho}{\partial l} - \nabla \cdot q\_{l} \\ \quad + \mathcal{S}\_{l} + Q\_{\nu l} + m\_{\nu l} h\_{\nu l} - m\_{l\nu} h\_{l\nu} \\ \frac{\partial}{\partial l} \left( \alpha\_{\nu} \rho\_{\nu} h\_{\nu} \right) + \nabla \cdot \left( \alpha\_{\nu} \rho\_{\nu} \nu\_{\nu} h\_{\nu} \right) &= \alpha\_{\nu} \frac{\partial p}{\partial l} - \nabla \cdot q\_{\nu} \\ \quad + \mathcal{S}\_{\nu} + Q\_{l\nu} + m\_{l\nu} h\_{l\nu} - m\_{\nu l} h\_{\nu l} \end{aligned} \end{cases} \tag{3}$$

Where α,ρ,v,S, p,τ , g, and h refer to volume fraction, density, velocity, source term, pressure, stress-strain tensor, gravity and specify enthalpy. m and Q are the interfacial mass and energy transfers from phase v to phase l. The force F includes five parts, FD, Flift, Fwl, Fvm, Ftd. They are the drag force, lift force, wall lubrication force, virtual mass force and turbulence dispersion force which will be introduced in detail respectively.

#### Interfacial Momentum Transfers

The interfacial momentum transfers between liquid and vapor phases are decided by the five forces which are the drag force, lift force, wall lubrication force, virtual mass force and turbulence dispersion force.

For fluid-fluid flow, each secondary phase is assumed to form droplets or bubbles. This assumption gives a method to derivate the drag force which is the most important force of these five parameters. Drag force refers to the force that droplets or bubbles suffer because of the different velocities. As we know, because of viscosity, the general form of the drag force is:

$$F = C\_D \cdot \frac{1}{2} \rho A V^2 \tag{4}$$

When it is used for a particle in the fluid, this form will be changed to

$$F = C\_D \cdot \frac{1}{2} A \rho\_l \left| V\_l - V\_v \right| \left( V\_l - V\_v \right) \tag{5}$$

Where A is the projected area of the particle, ρ<sup>l</sup> is the density of the continuous phase, V<sup>l</sup> and V<sup>v</sup> are the velocities of different phase. C<sup>D</sup> is a coefficient. For the droplets or bubbles, A can be expressed as A = 1 4 πd 2 , where d is the diameter.

Associated with Re number whose form is

$$Re = \frac{\rho\_l d}{\mu\_l} \left| V\_l - V\_\nu \right| \tag{6}$$

We can gain a new form of drag force of the droplets or bubbles which is

$$F = C\_D \cdot Re \cdot \frac{1}{8} \pi d \mu\_l \left(V\_l - V\_\nu\right) \tag{7}$$

So in per unit mixture volume, the drag force can be expressed as

$$F\_D = \frac{F}{\frac{1}{6}\pi d^3} = \frac{C\_D Re}{24} \cdot \frac{\rho\_v d}{6 \cdot \frac{\rho\_v d^2}{18\mu\_l}} \cdot \frac{\pi d^2}{\frac{1}{6}\pi d^3} \ (V\_l - V\_v)$$

$$= \frac{\rho\_v f}{6\pi\_l} dA\_i \ (V\_l - V\_v) \tag{8}$$

Where f = CDRe <sup>24</sup> is always called drag force coefficient which is estimated by Ishii model (Ishii, 1979). τ<sup>v</sup> = ρvd 2 18µ<sup>l</sup> is called particulate relaxation time. A<sup>i</sup> = 6 d is the interfacial area concentration which is an important value.

Lift force act on a particle mainly due to velocity gradients in the continuous phase flow field. It is calculated by the function from Drew and Lahey (1993).

$$F\_{\rm lift} = -C\_{\rm lift} \rho\_l \alpha\_\nu \left(\nu\_l - \nu\_\nu\right) \times \left(\nabla \times \nu\_l\right) \tag{9}$$

In this equation, Clift is the lift force coefficient provided by Moraga model (Moraga et al., 1999). This model is applicable mainly to the lift force on spherical solid particles, though it can be applied to liquid drops and bubbles.

Wall lubrication force which tends to push the bubbles away from the wall can be expressed as

$$F\_{\rm wl} = C\_{\rm wl} \rho\_l \alpha\_{\nu} \left| (\nu\_l - \nu\_{\nu}) \right|\_{x}^{2} n\_{\rm w} \tag{10}$$

Where (v<sup>l</sup> <sup>−</sup> <sup>v</sup>v) z is the phase relative velocity component tangential to the wall surface, and n<sup>w</sup> is the unit normal pointing away from the wall. The coefficient Cwl is calculated by Antal Model (Antal et al., 1991).

Virtual mass effect should be included when the secondary phase accelerates relative to the primary phase. It is defined as

$$F\_{\rm vm} = 0.5 \alpha\_{\nu} \rho\_{l} \left( \frac{\partial \nu\_{l}}{\partial t} + (\nu\_{l} \cdot \nabla) \, \nu\_{l} - \left( \frac{\partial \nu\_{\nu}}{\partial t} + (\nu\_{\nu} \cdot \nabla) \, \nu\_{\nu} \right) \right) \tag{11}$$

The turbulent dispersion force acts as a turbulent diffusion in dispersed flows. It plays an important role in driving the vapor away from the vicinity of the wall to the center of the channel. Bertodano (1991) proposed the following formulation:

$$F\_{td,l} = -F\_{td,\nu} = C\_{TD} \rho\_l k\_l \nabla \alpha\_{\nu} \tag{12}$$

Where CTD is a constant, it is equal to 1 by default in the following calculation.

#### Interfacial Energy Transfers

The interfacial energy transfers includes two parts which are the heat from liquid to vapor phase at the near wall region and the heat transfer between vapor and liquid phases in the subcooled bulk. As the bubbles depart from the wall and move to the subcooled region, there is heat transfer from the bubble to liquid that is defined as

$$q\_{lt} = h\_{sl} \left( T\_{sat} - T\_l \right) \tag{13}$$

The heat transfer coefficient can be computed using Ranz and Marshall (1952) model.

The interface to vapor heat transfer is calculated using the constant time scale return to saturation method which was provided by Lavieville et al. (2005). It is assumed that the vapor retains the saturation temperature by rapid evaporation and condensation. The formulation is

$$q\_{\nu t} = \frac{\alpha\_{\nu} \rho\_{\nu} C\_{p,\nu}}{8t} \left( T\_{sat} - T\_{\nu} \right) \tag{14}$$

In this equation, δt is the time scale.

#### Interfacial Mass Transfers

Depends on the theory, interfacial mass transfer can also be divided in to two parts: liquid evaporation near the wall and liquid evaporation or vapor condensation in the main stream. The former is calculated on the basis of evaporation heat flux which will be introduced below.

$$m\_E = \frac{q\_E}{h\_{f\nu} + C\_{p,l}T\_{sub}}\tag{15}$$

The latter can be calculated directly based on the interfacial energy transfer and latent heat hfv.

$$m = m\_{lt} + m\_{vt} = \frac{q\_{lt} + q\_{vt}}{h\_{fv}} \tag{16}$$

#### Wall Boiling Model and CHF Model

Different from the heat transfer of single phase flow on the near wall region, the heat transfer of two-phase flow includes three different types heat transfer. As we can see from the subcooled boiling flow, the energy is transferred directly from the wall to the liquid. Part of this energy named convective heat flux will cause the temperature of the liquid to increase and part which is called evaporative heat flux will generate vapor. Interphase heat transfer will also cause the average liquid temperature to increase; however, the saturated vapor will condense. In addition, there is a quenching heat flux which model the cyclic averaged transient energy transfer related to liquid filling the wall vicinity after bubble detachment. These basic mechanisms are the foundations of the Rensselaer Polytechnic Institute (RPI) models developed by Kurul and Podowski (1991).

To model boiling departing from the nucleate boiling regime, or to model it up to the critical heat flux and post dry-out condition, it is necessary to include the vapor temperature in the solution process. In such cases, some of the energy may be transferred directly from the wall to the vapor. So in total, the wall heat partition will be expressed as

$$q\_W = f\left(\alpha \mathfrak{l}\right) \cdot \left(q\_C + q\_Q + q\_E\right) + \left(1 - f\left(\alpha \mathfrak{l}\right)\right) \cdot q\_V \tag{17}$$

The four heat fluxes qC, qQ, qE, q<sup>V</sup> is defined as

$$qc = hc \left(Tw - T\_l\right) \left(1 - A\_b\right) \tag{18}$$

$$q\_Q = \frac{2k\_l}{\sqrt{\pi \lambda\_l T}} \left( T\_W - T\_l \right) \tag{19}$$

$$q\_E = V\_d N\_\text{w} \rho\_\text{v} h\_\text{fv} f \tag{20}$$

$$q\_V = h\_V (T\_W - T\_V) \tag{21}$$

Where TW, T<sup>l</sup> , T<sup>V</sup> are the temperature of the wall, liquid and vapor. h<sup>C</sup> and h<sup>V</sup> are the heat transfer coefficient. k<sup>l</sup> is the thermal conductivity. λ<sup>l</sup> = kl ρlCpl is the diffusivity. hfv is the latent heat of evaporation. ρ<sup>v</sup> is the vapor density. V<sup>d</sup> is the volume of the bubble based on the bubble departure diameter Dw. T is the periodic time. A<sup>b</sup> is the area of influence which is covered by nucleating bubbles. N<sup>w</sup> is the nucleate site density.

The nucleate site density N<sup>w</sup> is usually represented by a correlation based on the wall superheat.

$$N\_{\bowtie} = \mathcal{C}^{\prime} \left( T\_{\bowtie} - T\_{sat} \right)^{n} \tag{22}$$

The empirical parameters C = 210 and n = 1.805 are from Lemmert and Chawla (1977) model.

The bubble departure diameter D<sup>w</sup> means the maximum diameters of bubbles when bubbles leave the wall. The empirical correlation is given by Tolubinski and Kostanchuk (1970).

$$D\_{\rm w} = \min\left(0.0014, 0.0006e^{-\frac{T\_{\rm sub}}{45}}\right) \tag{23}$$

The area of influence A<sup>b</sup> is based on the departure diameter and the nucleate site density.

$$A\_b = \min\left(1, K \frac{N\_\text{w} \pi D\_\text{w}^2}{4}\right) \tag{24}$$

The coefficient K is given by Del Valle and Kenning (1985).

$$K = 4.8e^{\left(-\frac{Ia\_{\rm subb}}{80}\right)}\tag{25}$$

$$J\_{a\_{sub}} = \frac{\rho\_l C\_{pl} T\_{sub}}{\rho\_v h\_{fv}} \tag{26}$$

The last parameter is the frequency of bubble departure f = 1 T . In general, the bubble in the wall suffers three forces, which are buoyancy force, surface tension and drag force. When the heat flux is high and bubbles are large, the surface tension could be ignored. So according to the equilibrium of buoyancy force and drag force which was described in paper (Cole, 1960).

$$\lg\left(\rho\_l - \rho\_\nu\right) \frac{\pi D\_w^3}{6} = F\_D \frac{\pi D\_w^2}{4} \frac{\rho\_l \mu^2}{2} \tag{27}$$

Associated with u = fD<sup>w</sup> and F<sup>D</sup> = 1, f can be expressed as.

$$f = \sqrt{\frac{4\text{g } (\rho\_l - \rho\_\nu)}{\Im \rho\_l D\_w}}\tag{28}$$

In the model, the critical heat flux is defined to occur when the volume fraction of vapor reach given values which is provided by Tentner et al. (2006).

$$f\left(\alpha\_{\boldsymbol{\nu}}\right) = 1 - f\left(\alpha\_{l}\right) = \begin{cases} 0 & \alpha\_{\boldsymbol{\nu}} < \alpha\_{\boldsymbol{\nu},1} \\ \frac{1}{2} \left(1 - \cos\left(\pi \frac{\alpha\_{\boldsymbol{\nu}} - \alpha\_{\boldsymbol{\nu},1}}{\alpha\_{\boldsymbol{\nu},2} - \alpha\_{\boldsymbol{\nu},1}}\right)\right) \alpha\_{\boldsymbol{\nu},1} \le \alpha\_{\boldsymbol{\nu}} \le \alpha\_{\boldsymbol{\nu},2} \\ 1 & \alpha\_{\boldsymbol{\nu}} > \alpha\_{\boldsymbol{\nu},2} \end{cases} \tag{29}$$

$$\text{Where } \alpha\_{\nu,1} = 0.9 \text{ and } \alpha\_{\nu,2} = 0.95.$$

# NUMERICAL SCHEME

The numerical analysis is performed using the commercial CFD code FLUENT 16.0 (Ansys Fluent, 2015) which is part of ANSYS WORKBENCH 16.0. The modules we use are based on the finite volume method, and uses RANS (Reynolds averaged Navier-Stokes) equations for the solution of mass, momentum, energy and turbulence model.

# Geometry and Boundary Condition

CHF experiment data in an upward flow vertical pipe from Becker (1983) are used to validate the simulation method. Although there are much data in the Becker's experiment, one case with 7Mpa pressure which is condition of boiling water reactor (BWR) is investigated about the effect of grids and turbulence models. The facility has a 7m long heated pipe and the inner diameter is 10mm. Detail information about boundary condition are shown in **Figure 1**.

# Grids, Turbulence Models and Wall Functions

It has been tested that the grid type has limit influence on the temperature distribution in the single-phase flow (Chen et al., 2017b). But in two-phase flow, this conclusion is doubtful for the complexity of wall boiling status. So the effect of different grid type should be tested first before the investigation of turbulence models. As we can see from **Figure 2**, O type grid, arbitrary Hexaprism type grid and Triprism grid are calculated under the same boundary condition. Results are shown in section Effect of Grid Type.

All the turbulence models can be divided into high Reynolds and low Reynolds turbulence models according to its application. The standard k-ε model (Launder and Spalding, 1972) is typical high Reynolds turbulence model. It should be used for fully developed turbulence and isotropic flow. Because of its reduced computational complexity, it has been widely used for decades of years. An important weak point is that the standard k-ε model cannot deal with rotational flow or flow in curved pipe. To solve these problems, RNG (Re-Normalization Group) k-ε model (Yakhot and Orszag, 1986) and Realizable k-ε model (Shih et al., 1995) are developed in 1986 and 1995 separately. These two models expanded the application of k-ε model to a larger field. The secondary flow status of shear flow and rotational flow could be predicted to a certain extent. The k-ω models are low Reynolds turbulence model which have a larger application. The standard k-ω model (Wilcox, 1998) is more suitable for wall limited flow and free shear flow. In SST (Shear Stress Transport) k-ω model (Menter, 1994), the effect of turbulent shear force is included when solve the turbulent viscosity. It is also worth to note that SST k-ω model is same as standard k-ω model in near wall region and similar to standard k-ε model in the fully turbulent region.

The turbulence model should be used with suitable wall function or resolve the boundary layer with a fine-enough mesh without any wall functions. As we know, the near-wall region can be subdivided into three sublayers, i.e., viscous sublayer, buffer layer and fully turbulent region (Pope, 2000). For k-ε model with

standard wall function, it uses an empirical correlation to solve the equations in viscous sublayer and buffer layer. It also requires the first node should be located in fully turbulent region where viscosity is of less effect. So according to reference (Pope, 2000), Y+ larger than 30 is recommended for high Reynolds turbulence models to ensure accuracy and stability. However, the enhanced wall function maybe helpful to use high Reynolds turbulence model in the region where Y+ is smaller than 30. The results sometimes show good agreement with experiment data. For k-ω model, it is essential to ensure the viscous sublayer to be covered with several cells to maintain the near wall Y+ to 1. However, kω model could be Y+-independent if works with enhanced wall function. This is also how the k-ω model is used in FLUENT.

In order to validate the correlation of Y+, grids and turbulence models, five sets of grids are used under different turbulence models and wall functions. The detailed information about grid is listed in **Table 1**. The computational matrix of grid, turbulence model and wall function is shown in **Table 2**. StdWF, EnhWF, NeqWF, and ScaWF represent Standard Wall Function (Launder and Spalding, 1974), enhanced wall function, Non-Equilibrium Wall Function (Kim et al., 1999) and Scalable Wall Function, respectively. The detail information about range of application is shown in **Table 3**.

Not all the possible combinations are calculated in this paper for the numerous computational time and resource. Because of the same basic mechanism for the three k-ε models, the results of each can be used to the others. So the simulation of ICEM1 to ICEM5 on Realizable k-ε models and enhanced wall function are calculated to investigate the effect of different grids. This work is also done for the standard k-ω models on ICEM3∼ICEM5. The effect of turbulence models are compared based on ICEM3 and ICEM4. We believe the conclusion can be extended to the

TABLE 1 | Grids information.


TABLE 2 | Computational matrixes of grids, turbulence model and wall function.


other grids. For the wall function part, different wall functions can only be used on different k-ε models. So in the k-ω models, no wall functions are calculated except the enhanced wall function.

#### RESULTS AND DISCUSSION

The effect of grid type, grids number, turbulence and wall functions are discussed below. Results are compared among the highest temperature values and dry-out location of the inner surface of the pipe. The temperature distribution is the key saferelative parameter in the CHF experiment and should be paid more attention to discuss.

### Effect of Grid Type

As has been validated, the grid type has limit influence on the temperature distribution in single phase flow. But for two phase flow, a little difference is found among three grid types. The grid types are shown in **Figure 2**, and the numbers of nodes are 1344000, 1298500, and 1228500 separately. **Figure 3** shows the temperature distribution of the inner surface of the pipe. Results show all the three grids have the same dry-out position while the highest temperature values are different. More detail information can be found in **Figure 4** which is the radial temperature distribution at 4.5 m. The distribution of Hexaprism grid shows a little decentered phenomenon which is not in accordance with theory. Furthermore, the distribution in Triprism grid shows a zigzag in the center of the channel. Though this phenomenon can be avoided through refining the grid, it will cost more computing resource and time. So to predict the dry-out position and highest temperature value accurately, O type or axis-symmetric grid should be chosen for the calculation.

# Effect of Grids Number on Turbulence Models

Grid independent calculation is always an important requirement in the numerical simulation. It is necessary to estimate the quality


of gird through the calculations of different grids' level. However, it is not easy to gain fully grid-independent results for all the simulation, especially for two phase flow. But It has been tested that all the k-ε models and k-ω models with enhanced wall function can be used for a large scale of Y+ values from one to hundreds in single phase flow. Zhang et al. (2015) had also investigated the effect of grids and Y+ values in the subcooled boiling flow. It is concluded that results of k-ε models match well with experiment data when near-wall grid Y+ is larger than 11.25 and enhanced wall function cannot cope with the grids with wall adjacent grid Y+ near to 1. However, results of CHF calculation show a different conclusion which has a more strict limitation for Y+ value.

**Figure 5** show the comparison of the surface temperature under different grids on standard k-ε model and enhanced wall function. Results show little difference is detected among different grids except the highest temperature. The simulation results show a good match with experimental data. The dry-out positions of simulation results and experimental data are also in accord with each other. What's more, the difference of highest temperature decreased along with the refinement of the grids. So it may have the trend to gain grid-independent results.

The surface temperature results including ICEM4 are shown in **Figure 6**, and it is calculated based on RNG k-ε model and enhanced wall function. As we can see from **Figure 6**, the dryout position of ICEM4 falls behind that of the other three grids. The highest value of surface temperature of ICEM4 is also much lower than that of the other three grids and the experimental data. So the results start to be worse along with the refinement of grid and the decrease of the Y+ value.

For a further validation, the comparison results including the ICEM5 are shown in **Figure 7** which is calculated based on Realizable k-ε model and enhanced wall function. The temperature distribution of ICEM 4 shows same trend as that in **Figure 6**. The dry-out position is found behind the others and the highest temperature value is lowest. We can also see from **Figure 6** that when compared with good performance of the other grids, the temperature distribution of ICEM5 shows a large fluctuation in 0.7-2.5m. What's more, the highest temperature is lower than the other four. After applying ICEM3, ICEM4 and ICEM5 with Standard k-ω model, we also found same phenomena in the same position. The results are shown in **Figure 8**. The Y+ value of ICEM 5 is shown in **Figure 9**. After

comparing this two figure, we found that this unreasonable phenomena mainly appear near the wall when the Y+ value is lower than 50. The reason is that the transient behavior of bubbles has a great influence on the space averaged property when the near wall cells are too small.

#### Effect of Turbulence Models

Different turbulence models are applied to the CHF calculation with enhanced wall function based on ICEM3 and ICEM4. The results are shown in **Figures 10**, **11**. Some of the results show good agreement with experimental data while the others not.

**Figure 10** shows the inner surface temperature distribution of the pipe which is calculated under different turbulence models on ICEM3. Little difference is detected between Standard k-ω model and SST k-ω model. Only the highest value of temperature shows a difference about 2-5K which can be ignored. The results of Standard k-ε model and RNG k-ε model are nearly the same while both of these two are slightly lower than the results of Realizable k-ε model. The highest value of temperature of k-ω models is 30K higher than that of k-ε models which is a considerable difference. However, when compared with experimental data, all the simulation results are acceptable. The predicted dry-out position and temperature distribution are in accord with experimental data.

When ICEM4 is employed to the CHF calculation, large differences are found among the different turbulence models. The results are shown in **Figure 11**. The two k-ω models still show similar results in dry-out position and temperature distribution. Furthermore, these two results are also in agreement with experimental data except the post-dryout position. For k-ε models, the dry-out positions are predicted behind the experimental data and the highest temperatures are also lower

FIGURE 6 | Surface temperatures under different grids on RNG k-ε model and enhanced wall function.

than the measurement. It is worth to say that all the above calculations are under enhanced wall function. The effect of wall functions will be discussed in next section.

FIGURE 11 | Surface temperatures under different turbulence models on ICEM4.

#### Effect of Wall Functions on k-ε Models

Heat transfer and flow characters in near wall regions are important to the calculation of temperature and velocity in two-phase flow. The wall functions and wall boiling model are also influenced by each other. Because the wall functions use different empirical correlations to solve the equations in viscous sublayer and buffer layer, different wall functions will have different performance. In FLUENT, k-ω models are low Reynolds models and set to be used with wall functions unless the grid Y+ doesn't meet the requirements of the low Reynolds number modeling. If the requirements are met, then FLUENT doesn't apply wall function but resolves the near wall region. So there is no need to discuss the effect of different wall functions when use k-ω models.

Based on ICEM3, the CHF calculation results of k-ε models under different wall functions are shown in **Figure 12** separately. All these three pictures show same trend for different wall functions. The temperature of non-equilibrium wall functions is 20K higher than that of the others while the enhanced wall function has the lowest temperature. The temperature distributions of standard wall function and scalable wall function are the same. This is determined by the correlation of scalable wall function.

The purpose of scalable wall functions is to force the usage of the log law in conjunction with the standard wall functions approach. This is achieved by adding a limiter in the Y+ calculations Y+ is calculated by

$$Y + =MAX\left(Y + , \text{ 11.25}\right) \tag{30}$$

So when the Y+ is larger than 11.25, Y+ is the same with original value, and then the scalable wall function is the same as standard wall function. Otherwise, Y+ is equal to 11.25 which bring different results. For the cases in this paper, all the Y+ values are larger than 11.25, so there is no difference for the results between these two wall functions. When compared with experimental data, the standard wall function and scalable wall function show the best performance while results of the other two are also acceptable.

Based on ICEM4, the comparison is only investigated between standard wall function and enhanced wall function on two kε models. The results are shown in **Figure 13**. As mentioned above, the predicted dry-out positions of enhanced wall function

are behind the measurement when ICEM4 is employed to the simulation. Among these four cases, only Realizable k-ε model with standard wall function can give an accurate dry-out position and acceptable temperature distribution.

# CONCLUSION

CHF phenomenon in a vertical pipe is investigated by use of CFD simulation. Simulations are calculated among different grids type, grids number, different turbulence models and wall functions to detect the effect of these parameters. Several conclusions are obtained after comparing the simulation results and experimental data. The conclusion about the effect of grid

# REFERENCES


and turbulence models can help someone reduce their workload on mesh validation and set up solve algorithm quickly.


# AUTHOR CONTRIBUTIONS

XD is the main author of the paper, he finished the research and write. ZZ is the supervisor and give the guide for the other authors. DL is the co-author who works for the comparison of the experiment and simulation. ZT is the associated supervisor and give the reference too. GC works with the other authors for the discussion part.

#### ACKNOWLEDGMENTS

This work is supported by the Research on Key Technology of Numerical Reactor Engineering (J121217001), and the authors also appreciate deeply the suggestions from Dr. Tenglong Cong.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Dong, Zhang, Liu, Tian and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Neutronics and Thermal Hydraulics Analysis of a Conceptual Ultra-High Temperature MHD Cermet Fuel Core for Nuclear Electric Propulsion

Jian Song<sup>1</sup> \*, Weijian An<sup>2</sup> , Yingwei Wu<sup>1</sup> \* and Wenxi Tian<sup>1</sup>

<sup>1</sup> Department of Nuclear Science Technology, Xi'an Jiaotong University, Xi'an, China, <sup>2</sup> Department of Reactor Engineering, China Institute of Atomic Energy, Beijing, China

Nuclear electric propulsion (NEP) offers unique advantages for the interplanetary exploration. The extremely high conversion efficiency of magnetohydrodynamics (MHD) conversion nuclear reactor makes it a highly potential space power source in the future, especially for NEP systems. Research on ultra-high temperature reactor suitable for MHD power conversion is performed in this paper. Cermet is chosen as the reactor fuel after a detailed comparison with the (U,Zr)C graphite-based fuel and mixed carbide fuel. A reactor design is carried out as well as the analysis of the reactor physics and thermal-hydraulics. The specific design involves fuel element, reactor core, and radiation shield. Two coolant channel configurations of fuel elements are considered and both of them can meet the demands. The 91 channel configuration is chosen due to its greater heat transfer performance. Besides, preliminary calculation of nuclear criticality safety during launch crash accident is also presented. The calculation results show that the current design can meet the safety requirements well.

Keywords: space reactor, neutronics, thermal hydraulics, reactor design, cermet fuel, MHD generator, electric propulsion

# INTRODUCTION

Interest in space exploration has developed rapidly in recent years and brought extensive research activity in the area of space propulsion that goes beyond present-day chemical systems (Valentian and Snecma, 2003; Wesley et al., 2011; James and Bruce, 2012). Among the available energy alternatives nuclear power offers important advantages and in many cases is the only viable alternative given actual operation conditions in interplanetary missions for that is the most compact form of energy available (Viorel, 2009). Nuclear propulsion systems, including nuclear electric propulsion (NEP) and nuclear thermal propulsion (NTP), are indispensable for manned deep space exploration. NTP can provide high specific impulse which is more than 2 times as the best chemical propulsion system, while NEP can give even much larger specific impulse than NTP (Thomas et al., 2005). Thus maximum payload capability and lowest launch cost can be achieved with NEP system. However, due to the low thrust capability of current electric propulsion system, large electric power is needed to maximize the thrust and minimize the trip time which is critical to the crew.

Power conversion efficiency is essential to NEP system. Among various power conversion systems, e.g., dynamic conversion cycles such as Stirling, Brayton, Rankine cycles, and static conversion systems such as thermoelectric, thermionic, Alkali-metal Thermoelectric Converter (AMTEC), and Magnetohydrodynamics (MHD), MHD possesses the highest efficiency (>40%), especially suitable for high-power space applications. Minimum thermal power, radiator area, and

#### Edited by:

Jun Wang, University of Wisconsin-Madison, United States

#### Reviewed by:

Xin Li, Japan Atomic Energy Agency, Japan Yago Rivera, Instituto Ingeniería Energética, Universidad Politécnica de Valencia (UPV), Spain

#### \*Correspondence:

Jian Song sjian\_eagle@163.com Yingwei Wu wyw810@mail.xjtu.edu.cn

#### Specialty section:

This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research

Received: 14 February 2018 Accepted: 29 March 2018 Published: 16 April 2018

#### Citation:

Song J, An W, Wu Y and Tian W (2018) Neutronics and Thermal Hydraulics Analysis of a Conceptual Ultra-High Temperature MHD Cermet Fuel Core for Nuclear Electric Propulsion. Front. Energy Res. 6:29. doi: 10.3389/fenrg.2018.00029 shielding mass can be realized with MHD power conversion system. The basic principle of MHD generation is to heat the working fluid to ionized plasma using reactor thermal power, then fast flow through a channel surrounded by perpendicular high magnetic field to produce electrical current under the effect of Lorentz force, demonstrated in **Figure 1**. The hot working fluid out from the power generation channel preheats the cold working fluid coming into the reactor in the heat exchanger and will be further cooled in the heat radiator, and finally passed through the compressor and turbine to the next cycle. Due to its high temperature, it is also possible to add a Brayton heat engine in the cold leg to enhance the efficiency of power generation. However, long-term operation in ultra-high temperature condition (>2,000 K) is a significant drawback to MHD system, and only NTP fuels may meet the requirements.

Recently, a concept of manned Mars exploration mission with NEP using closed MHD generator is proposed by China's Lunar and Deep Space Exploration Center (CLEP). The goal is to complete the space nuclear power supporting the interplanetary spacecraft, to achieve high efficiency MHD nuclear power conversion technology. **Table 1** lists the main parameters of the overall spacecraft power system. This paper focuses on the

TABLE 1 | Main parameters of the whole system.


**Abbreviations:** B4C, boron carbide; BeO, beryllium oxide; βeff, effective delayed neutron fraction; Cermet, ceramic-metal fuel type; CFD, computational fluid dynamics; CLEP, China's lunar and deep space exploration center; keff, effective neutron multiplication factor; Gd2O3, gadolinium oxide; LiH, lithium hydride; MW, megawatts of energy; MHD, magnetohydrodynamics; NEP, nuclear electric propulsion; NTP, nuclear thermal propulsion; NASA, national aeronautics and space administration; ρ, excess reactivity; UO2, uranium oxide; W, tungsten; Re, rhenium.

reactor core design and optimization. A preliminary reactor concept using cermet fuel is carried out. The current modeling tool includes an iteration between reactor physics and thermal hydraulics calculations. In order to take the advantage of different software, MCNP5 and MVP are used for the nuclear reactor core modeling, neutronics, shielding, depletion analyses, and critical safety analysis. ANSYS-CFX is used for thermal hydraulic analysis and optimization.

# REACTOR DESIGN

A reasonable design is provided that appear desirable from the standpoint of size, mass, safety, and reliability. The design work mainly consists of neutronics, thermal hydraulic and radiation shield.

# Fuel Consideration

Apparently, the critical performances of the reactor lie in the ultra-high temperature and long term duration. Three kinds of ultra-high temperature fuel derived from mainstream NTP programs are considered, i.e., (U,Zr)C graphite-based fuel, mixed carbide fuel, and cermet fuel (Jahshan and Kammash, 2005; Jonathan et al., 2011).

Extensive research work was conducted of (U,Zr)C graphite-based fuel for thermal spectrum reactors during the Rover/NERVA program (1955–1972) that proved NTP to be a viable technology. However, there was still some issue, i.e., cracking of ZrC coating on the coolant channels, remained unresolved till the end of the program. It is a significant limiting factor to the reactor lifetime because it may lead to a catastrophic loss of graphite and fuel (Stewart and Bruce, 2012; Benensky, 2013). Since the lifetime of the reactor required by NEP is much longer than NTP, (U,Zr)C graphite-based fuel is considered not suitable for the present MHD reactor design.

Mixed carbide fuel forms, including di-carbides and tricarbides, were developed by Soviet Union in heterogeneous reactor program for NTP (1950s-early 1990s) which shows several advantages. A "twisted-ribbon" geometry was chosen for fuel form to maximize heat transfer while preserving fuel shape and integrity under working conditions. However, the tricarbide is very brittle and instability at high temperatures testing which makes the development of these fuel forms challenging (Benensky, 2013). Additionally, the manufacturing processes of this kind of fuel forms are very complicated and costly. So mixed carbide fuel is beyond the scope of consideration in the expectable future.

Cermet, another option of NTP fuels, were mainly pursued in parallel to basic graphite fuel development in GE-710 (by General Electric) and ANL (Argonne National Laboratory) fast spectrum propulsion reactor programs in 1960s. Current research on cermet fuel is being conducted as part of the Nuclear Cryogenic Propulsion Stage (NCPS) Advance Exploration System (AES) technology project at NASA's Marshall Space Flight Center (Gomez, 2014). A supercritical carbon dioxide cooled Mars surface power reactor with Brayton cycle designed by the Center for Space nuclear Research (CSNR) of Idaho National Laboratory also uses low-enriched uranium cermet fuel (Kevin et al., 2016). Cermet fuel forms are composed of a tungsten or molybdenum fuel element matrix material (Benensky, 2013). It is considered that cermet fuels may provide a significant design advantage over graphite fuel elements in the NERVA/Rover program, such as its potential for very long operating life and ability to restart (Stewart and Bruce, 2013). In particular, researches show that deformation but not cracking of the cermet fuel and cladding may occur under high stresses (Stewart and Bruce, 2013). Fast spectrum reactors also tend to be much more compact than thermal spectrum reactors and an inherent "neutronic spectral shift" safety feature that helps maintain reactor subcritical when submerged in water or in wet sand of the reentry accident. Therefore, cermet fuel element developed under the GE-710 program is chosen as a reference fuel form for the MHD conversion reactor design.

The appearance and design of fuel elements are influenced by several factors. Important geometry features includes the element size, the number of coolant channels, the channel diameter and the web thickness of the fuel matrix (Bruce and Stanley, 2012). In the application of W-UO<sup>2</sup> cermet based fuel elements for MHD reactor power system, the cermet fuel matrix has a composition of W-60%UO2-6%Gd2O<sup>3</sup> with W-25%Re as the cladding material both on the coolant walls and exterior. Specially, volumetric ratio 60% UO<sup>2</sup> and volumetric ratio 6% Gd2O<sup>3</sup> (as the oxygen stabilizer of UO2) are dispersed in the tungsten matrix. Each fuel element has a hexagonal cross section with a flat-to-flat distance of 2.316 cm. In the axis direction, the first 10.24 cm of fuel element inlet is filled with a beryllium oxide neutron reflector and the latter 60.96 cm is filled with the cermet fuel mixture. There is no neutron reflector at the outlet of fuel element because the core outlet temperature is too high. In order to get the optimal heat transfer performance, two kinds of coolant channel configuration designs were proposed. One configuration design based on the GE-710 reactor was modeled, where each fuel hex contained 61 coolant channels lined with 0.286 mm thick W-25%Re cladding layer and another was modeled to have 91 coolant channels each lined with a 0.204 mm same cladding layer. Design parameters of each fuel element are shown in **Table 2**. The surface-area tovolume-ratio of the 61 coolant channel hex is 0.57 compared to the 0.73 ratio with the 91 coolant channel hex. **Figure 2** illustrates


the hexagonal prismatic cermet fuel element that contains 91 coolant channels.

#### Reactor Description

The reactor core configuration is shown in **Figure 3** and be designed as helium-cooled fast reactor. The radius of active zone is 21 cm. There're 246 fuel elements and 66 tungsten filler elements in the core. The <sup>235</sup>U enrichment of central 30 fuel elements is 70%, while the outer ones possess 93% <sup>235</sup>U enrichment. The purpose of this design is to flatten the power

TABLE 3 | Major design parameters of the reactor core.


distribution and reduce the nuclear hot spot factor. The radial and axial reflectors are both made of BeO with the thickness of 10 cm. There are two ways to control the reactor reactivity. 16 absorbed drums in the radial reflector are mainly used to control the power level and can also be used to shut down the reactor independently. These control drums are divided into eight sets, with two control drums in each set sharing one drive mechanism and rotating at the same time. Each of drum has 120 degree surface absorber made of B4C and the thickness of absorber is 2 cm. A safety rod in the center of core with 4 cm diameter is used to increase the shutdown depth, which is particularly important for critical safety in launching accident. In normal operation, the safety rod is pulled out of the core. The overall height of the core is 71.2 cm and the total mass is 1616.9 kg. All the major core design results are summarized in **Table 3**.

#### Shielding Design

Shielding is an important part of the reactor system, mainly used for absorbing or weakening the large amount of neutrons and gamma rays released from the core. In this research, the main design goal of shielding design is to make sure the maximum

irradiation dose of the astronauts is limited to 50 mSv per year at 51 m away from the reactor core. **Figure 4** provided a 17◦ truncated cone-shaped preliminary shadow shielding design meets the demand, which can protect crew or equipment in specific areas. The total mass of radioactive shield is 9.98 t, with heavyweight shields on both sides and lightweight shield arranged in the middle of them. The lightweight shield is made of LiH, the best neutron shielding material known so far, which is very advantageous for reducing the weight of shielding. The heavyweight shield is made of W and is a commonly used material for shielding gamma rays in space reactors. The design of heavyweight shields on both sides can not only shield the gamma rays released by the core, but also shield the induced gamma rays and protect the lightweight shield. This shielding structure can be further optimized in the future to reduce mass.

#### CORE NEUTRONICS ANALYSIS

Neutron physics calculation shows that the keff, ρ and βeff is 1.0417, 4.001 1k/k% and 6.91 × 10−<sup>3</sup> for cold state when the control drums are rotated outward at the beginning of core-life, respectively. While the keff, ρ and βeff is 1.0395, 3.78 1k/k% and 6.05 × 10−<sup>3</sup> for hot state, respectively. It can be seen from the data that the reactor has a negative temperature reactivity coefficient because the excess reactivity of hot state is lower than that of cold state, which is a key assurance for the inherent safety of nuclear reactor. When all the control drums are rotated inward (ρ = 6.136 1k/k%) without the safety rod, the reactor can be shutdown with the keff equal to 0.9779. The reactivity worth of safety rod is 3.254 1k/k%.

The total thermal power of the reactor core is 25 MWth. Normalized thermal power distribution of the reactor core is shown in **Figure 5**. It has a normalized peak-to-average thermal power ratio of 1.253. The corresponding maximum thermal power of the fuel elements (hot channel) is 0.1136 MWth. The average thermal power per fuel element (average channel) is 0.09158 MWth. **Figure 6** shows the axial power distribution of the hot channel as well as the average channel. There is an obvious local power peak near the axial reflector due to reflected neutrons increasing fission reaction rate. In the calculation, it is assumed that the radial power distribution of each fuel element is uniform.

Depletion calculation by MVP-BURN indicated that the operating lifetime of reactor is about 2 years which has 33% margin as the real operating time is about 1.5 years.

Apart from reactor physics design, preliminary critical safety analysis was also carried out using MCNP5. The reactor has to satisfy the stuck rod criterion: During the shutdown process, the

FIGURE 5 | Normalized thermal power distribution of 1/4 core.

reactor can be kept in subcritical with other control rods when the most valuable control rod is stuck outside of the core. As **Figure 7** illustrates, when a set of control drums are stuck facing outwards of the reactor core and other control drums are turned to the core, the keff of the core is 0.9899. Therefore, the current reactor physics design can meet the requirement of stuck rod criterion.

Launch failure reentry accident is unique to space reactor because the usual using of highly enriched fuel elements, which can easily lead to prompt supercritical accident. The thermal

neutron absorption cross-section of cladding material tungstenrhenium alloy (especially rhenium) is very large while the fast neutron absorption cross-section of that is very small, which is very helpful in operating efficiency and critical safety. Two severe conditions during launch are calculated, and the results show very inspiring safety performance. When the reactor is submerged in water and the inner voids (i.e., coolant channels and other clearance) are all flooded with water, the relevant keff is 0.9652. In another condition, the keff is 0.9664 when the reactor immersed in wet sand and the inner voids become full of water (as shown in **Figure 8**), which is the most serious critical accident scenario for space nuclear reactors. In the above calculation, all of the drums and safety rod were assumed at the largest position of reactivity worth.

# CORE THERMAL HYDRAULICS ANALYSIS

The basic task of thermal hydraulics design is to ensure that the heat transfer capability of reactor core matches with the thermal power, and all the thermal parameters are maintained within the prescribed limits during operation with suitable safety margin. On the other hand, the core outlet temperature must be high enough (more than 2,200 K) in operating conditions to ensure the efficiency of MHD conversion. The core thermal power distribution is the main input parameter for thermal hydraulics analysis. In this section, the initial reactivity lifetime thermal power distribution calculated by MCNP5 was converted into the volumetric heat release rate and then put into ANSYS-CFX for thermal hydraulics calculation. Because the heat release peak power factor of fuel element in the end of reactivity lifetime is lower than that in the beginning, which will reduce the corresponding heat transfer requirements. Single channel steady state model was used for CFD calculations, contains average channel fuel element and hot channel fuel element.

The physical properties of core materials and coolant used in the thermal hydraulics analysis are listed in **Table 4**. These physical properties are constants at rated operating conditions and are suitable for steady state thermal hydraulics analysis. The density of gaseous helium depends on pressure and temperature. It also can be seen from **Table 4** that the heat transfer performance of tungsten based cermet fuel is very good.

The reactor is cooled by helium at a very high velocity, which belong to compressible flow since the flow velocity is very fast in the core. The helium enters the core channel at a temperature of 800 K and reference pressure is 1 MPa. The assumption in the calculation was that the flow rate of coolant in each channel is evenly distributed. The fluid domain was adopting the RANS technique which solve time averaged equations instead of solving the Navier-Stokes equations directly, which has been widely

TABLE 4 | Physical properties of reactor core materials and coolant.


used in the engineering applications. The solid domain was solved by heat conduction equations with the no-slip boundary condition. The heat conduction and heat convection between fluid domain and solid domain were considered in this model. Thermal radiation inside the fuel element was neglected in the preliminary study. The shear stress transport (SST) κ-ω model was adopted to simulate the turbulent phenomenon. The most significant advantage of this model is that it can effectively integrate the advantages of different turbulent model: robust and accurate formation of the κ-ω model in the near wall region and free-stream independence of the κ-ε model in the far field.

The semi-structured mesh was generated using the ICEM software, show as **Figure 9**. A 1/6 sector (60◦ symmetry) was used to simplify the model and save calculation time due to

a zoomed-in view of center part.

the symmetry of fuel element geometry. Symmetrical boundary condition was adopted to the corresponding split surfaces of the sector and adiabatic boundary condition was adopted to the rest of solid outer surface. We used thermal equilibrium boundary conditions at the interface between different parts in the fuel element. An unstructured grid was adopted to the fuel matrix part, whereas a structured grid was applied to the coolant channels and cladding. In this way, it can be ensured that all the meshes quality are >0.7 when the aspect ratio of the fuel element is large. The models of two fuel element design cases contains approximately 1.89 and 2.76 million finite volume cells respectively, allowing the thermal power profile to be calculated with an axial mesh spacing of 1.78 mm. Mesh sensitivity analysis has been conducted.

The calculation results proved that there is neither obstruction nor supersonic flow in the process of coolant flow in the core. The calculation results of average channel determined that the core heat can be taken away safely and meet the design requirements when the coolant mass flow rate is 3.0996 kg/s. **Figure 10** shows the axial temperature profile of the average channel for 91 coolant channel fuel element case. It can be seen that both coolant and fuel element temperature gradually increase in the axial direction, and the maximum value occurs at the coolant outlet. The temperature distribution between bocks in different axial regions is not smooth transition because of the segmented distribution of volumetric heat release rate given by neutron physics. There are two reasons that cause the fuel temperature achieve the highest at coolant channel outlet (Stewart and Bruce, 2008): Firstly, the large specific heat capacity and fast flow velocity of helium, and the large coolant channel length to diameter ratio (799:1) are conducive to heat transfer; secondly, the thermal conductivity of tungsten based dispersion fuel is excellent which can help to promote the temperature flattening in the axial direction.

As demonstrated in **Figure 11**, the average coolant outlet temperature of 61 and 91 coolant channel element cases are about 2,250 and 2,217 K respectively while the peak fuel temperature are approximately 2,343 and 2,287 K respectively. The total core pressure drop are 1.067 MPa and 1.225 MPa respectively. For hot channel (as shown in **Figure 12**), the average helium outlet temperature for 61 and 91 coolant channel cases was about 2,595 and 2,558 K respectively while the peak fuel temperature is approximately 2,699 and 2,648 K respectively. The total core pressure drop is 1.162 and 1.328 MPa respectively. It can be seen that the average fuel temperature is always slightly higher than the coolant temperature, which also proves the good heatconducting property of W-UO<sup>2</sup> cermet fuel. Obviously, the fuel temperature of 61 channel case is 50–150 K higher than 91 channel case in the second half of hot channel. Larger volume ratio on the fuel element edge means more heat drains into edge channels, hence temperatures are higher. Due to the higher surface area to volume ratio, the peak temperature is lower for the 91 channel case than that for the 61 coolant channel case. It is foreseeable that this difference will be more distinct as the reactor power increases or the coolant flow rate drops.

The radial temperature profile of 91 coolant channel case and 61 coolant channel case in **Figures 13**, **14** show the distribution of outlet temperature are substantially uniform distribution. The fuel zone near the edges has coolant channels on only one side, so the heat exchange capacity is relative low. The edge temperature is approximately 8 K higher than center temperature under the assumption of adiabatic boundary condition. Coolant velocity distribution is shown in **Figure 15**. The maximum rate is located in the center of coolant channel proved that the boundary layer of mesh can simulate flow characteristics very well.

In conclusion, both of the configurations meet the requirements well. And the hot channel analysis indicates that the maximum temperature of the fuel is far below the

TABLE 5 | Thermal-hydraulic parameters of the reactor core.


limiting temperature (3,100 K) of cermet fuel (Robert et al., 2012). The heat transfer performance of 91 channel case is better than 61 channel case while pressure drop result is on the contrary. As for reactor design, we would prefer to the lower hot-spot temperature because of the mechanical properties requirement for long term operation in high temperature condition. So the 91 coolant channel configuration is chosen for preliminary design of the reactor finally. **Table 5** shows the reactor core thermal-hydraulic design parameters.

In normal operating condition, the minimum mass flow rate of coolant is 2.51 kg/s. The hot-spot temperature of fuel and pressure drop is approximately 3,110 K and 1.081 MPa,

#### REFERENCES


respectively. The results show that the reactor has sufficient thermal safety margin of 23.5%.

# CONCLUSIONS

Manned Mars exploratory missions under CLEP research are progressively enabling the combination of space technology and nuclear technology for the theoretical breakthrough of high power space nuclear power and electric propulsion. A cermetfueled reactor with MHD power conversion for NEP is designed and analyzed, which has the advantage of high energy conversion efficiency and high power density. Computation results show the current core design can meet the requirements very well. Furthermore, extensive research is necessary for preliminary ground demonstration verification, including full core CFD simulation, fluid-solid coupled analysis, transient analysis, and necessary experimental works.

However, cermet fuel technology, which is still being under development at NASA, is the key issue for the whole reactor concept and application. Many issues, especially degradation of mechanical integrity and loss of fuel, remain to be resolved (Haertling and Hanrahan, 2007). Meanwhile, manufacturing a reliable high-speed helium turbine is also a challenging technology (Isaac et al., 1993).

### AUTHOR CONTRIBUTIONS

JS: completed the thermal hydraulic design of the reactor system and wrote this article; WA: completed the reactor physics design of the reactor system; YW: funded the research and provided a computing platform for the reactor system design; WT: contributed to the cooperation between Xi'an Jiaotong University and the China Institute of Atomic Energy to jointly accomplish this research.

#### ACKNOWLEDGMENTS

The authors express their thanks to staff at the Department of Reactor Engineering in CIAE who supported this work for many useful discussions. The research works of this subject is funded by the National 863 High Technology Research and Development Program (2012AA050908) and National Natural Science Foundation of China (Grant No. 115 75141).

uranium dioxide fueled cermet materials. J. Nucl. Mater. 366, 317–335. doi: 10.1016/j.jnucmat.2007.03.024


in 41th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit (Tucson, AZ).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Song, An, Wu and Tian. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Hao Sun, Chenglong Wang\*, Xiao Liu, Wenxi Tian, Suizheng Qiu and Guanghui Su*

*Department of Nuclear Science and Technology, Xi'an Jiaotong University, Xi'an, China*

Underwater vehicle is designed to ensure the security of country sea boundary, providing harsh requirements for its power system design. Conventional power sources, such as battery and Stirling engine, are featured with low power and short lifetime. Micronuclear reactor power source featured with higher power density and longer lifetime would strongly meet the demands of unmanned underwater vehicle power system. In this paper, a 2.4 MWt lithium heat pipe cooled reactor core is designed for micronuclear power source, which can be applied for underwater vehicles. The core features with small volume, high power density, long lifetime, and low noise level. Uranium nitride fuel with 70% enrichment and lithium heat pipes are adopted in the core. The reactivity is controlled by six control drums with B4C neutron absorber. Monte Carlo code MCNP is used for calculating the power distribution, characteristics of reactivity feedback, and core criticality safety. A code MCORE coupling MCNP and ORIGEN is used to analyze the burnup characteristics of the designed core. The results show that the core life is 14 years, and the core parameters satisfy the safety requirements. This work provides reference to the design and application of the micronuclear power source.

#### *Edited by:*

*Jun Wang, University of Wisconsin-Madison, United States*

#### *Reviewed by:*

*Keyou Mao, Purdue University, United States Claudio Tenreiro, University of Talca, Chile*

> *\*Correspondence: Chenglong Wang*

*chlwang@mail.xjtu.edu.cn*

#### *Specialty section:*

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

*Received: 04 December 2017 Accepted: 28 February 2018 Published: 22 March 2018*

#### *Citation:*

*Sun H, Wang C, Liu X, Tian W, Qiu S and Su G (2018) Reactor Core Design and Analysis for a Micronuclear Power Source. Front. Energy Res. 6:14. doi: 10.3389/fenrg.2018.00014*

Keywords: underwater, heat pipe cooled reactor, MCNP, ORIGEN, core design

# INTRODUCTION

Nuclear power source can be used in the power systems of underwater vehicle. The future unmanned underwater vehicle (UUV) is designed as a multipurpose equipment (Navy, 2004). Its power system requires more compact structure, higher energy density, longer lifetime, higher power output, and higher reliability (Li and Wang, 2003). Compared with the conventional power systems like storage battery, nuclear power source featured with higher power density, and core lifetime of more than 10 years much better meets the future demands of UUV power system.

There are two types of micronuclear reactor power source (Rowe, 1991): (1) isotope generator, which converts decay heat into electric energy; (2) nuclear reactor power source, which converts fission heat into electricity energy. Comprehensively considering the reactor mass and volume, criticality safety, reliability, and maneuverability, heat pipe cooled fast reactor power source featured with compact structure, less movable parts, high reliability, and low noise level could be widely used in the power system of unmanned underwater vehicle in the future (Rowe, 1991).

Heat pipe cooled reactor has already been widely researched. Various micro heat pipe cooled reactor power sources for space missions have been designed in the United States. For instance, HOMER (Poston, 2000) is a series of heat pipe reactors applied for the moon and mars missions. Potassium or sodium heat pipe are adopted for cooling. Stirling engine is used to generate electric energy, whose efficiency is more than 20%. MSR (Bushman et al., 2004) is another micro heat pipe cooled reactor power source featured with electric power for the scale of 100 kilowatts, lithium heat pipe, and thermoelectric generator, thermal power of 1.2 MWt, and efficiency of more than 10%. SAIRS (El-Genk, 2008) is a kind of sodium heat pipe-cooled fast reactor featured with electric alkali metal thermoelectric converter, with efficiency of more than 20%. LEGO-LRCs (Bess, 2008) is a sodium heat pipe cooled reactor with Stirling engine. SAFE-400 (Poston et al., 2002) is a heat pipe-cooled reactor system featured with Braydon cycle producing 400 kWt thermal power. China institute of atomic energy has put forward a variety of heat pipe reactor designs for space missions, such as the mars surface power plant, and the lunar surface power plant HPCMR (Chengzhi et al., 2015), etc.

Based on a literature review, a micro heat pipe cooled nuclear reactor power source applied for underwater vehicle featured with 2.4 MWt power output and a lifetime of more than 10 years is designed. Lithium heat pipe cooled core, six control drums, tungsten, and an open water loop shield are adopted in this power source. Monte Carlo program and ORIGEN are used to preliminarily calculate design parameters and analyze the criticality safety and burnup characteristics of the design scheme, etc.

#### SYSTEM DESIGN

The nuclear reactor power source consists of the following parts: core, control drums, shield, heat pipes, and thermoelectric generator. The working principle is shown in **Figure 1**. In reactor core, fission fuels generate heat in chain reaction controlled by control drums. The heat is transferred by heat pipes and is converted into electrical power by thermoelectric generator. Waste heat is taken away by water.

The core structure is shown in **Figure 2**, and the key parameters of reactor core are shown in **Table 1**. The core is made up of 90 fuel pins, 37 heat pipes, BeO reflector, and 6 control drums. **Figure 3** shows the design parameters of fuels and heat pipes. High-enriched

uranium nitride (with 70% concentration of 235U) fuel and high temperature lithium heat pipes are adopted and are put in a hexagon Nb–1Zr alloy matrix in triangle array. The matrix is put in a barrel on whose inner surface is a coating of Gd2O3 burnable poison. The barrel is surrounded by a BeO reflector and 6 control drums with B4C. Reactor core is stored in a cylinder Mo–14Re alloy barrel. Tungsten and an open water loop are used as shields in the power source system. Thermoelectric generator is used to convert heat conducted from heat pipes into electricity of 120 kWe.

### NEUTRONICS CALCULATION

#### Calculation Methods and Model

Monte Carlo code MCNP5 is a probabilistic code developed by Los Amos National Laboratory (LANL) in early 1970s. Monte Carlo program has a significant advantage of simulating geometries without much approximation. Effective multiplication factor (*k*eff) is one of the most important parameters in criticality calculation using MCNP. It is defined as the following equation (Team, 2003):

$$\begin{split} k\_{\text{eff}} &= \frac{\text{fission neutrons in generation } i + 1}{\text{fission neutrons in generation } i} \\ &= \frac{\rho\_t \int\_{\upsilon} \int\_{\alpha}^{\upsilon} \int\_{x} \int\_{\Omega} \upsilon \sigma\_t \, \Phi \text{d} \, \text{Vdd} \, \text{Ed} \, \Omega}{\int\_{\upsilon} \int\_{\alpha}^{\upsilon} \int\_{x} \int\_{\Omega} \nabla \bullet \, \text{fdV} \text{d} \, \text{Ed} \, \Omega + \rho\_t \int\_{\upsilon} \int\_{\alpha}^{\upsilon} \int\_{x} \int\_{\Omega} (\sigma\_{\epsilon} + \sigma\_{f} + \sigma\_{m}) \Phi \text{d} \, \text{Vdd} \, \text{Ed} \, \Omega} \end{split} (1)$$

where the phase-space variables are *t, E,* and Ω for time, energy, direction, and implicitly *r* for position with incremental volume

Table 1 | Core design parameters.


As **Figure 4** shows, a 1:1 model is set up to simulate the core neutronics characteristics. The active zone is axially divided into 17 layers and 127cells per layer. In criticality calculation, the KCODE card is used to build a critical source. The importance of all the particles, neutron and photon, in all cells is defined to be equal to 1.

The calculation requires an ACE format nuclear database. Thus, NJOY (a nuclear data processing program developed by LANL) and the ENDF-B-VII.0 nuclear data are used to build a database for MCNP5 in further calculation. According to preliminary thermal analysis, the fuel temperature is assumed 2,000 K, temperature of the heat pipes and matrix is 1,200 K,

and temperature of the control drums and reflector is 900 K for neutronics analysis. A code MCORE (Zheng et al., 2014) coupling MCNP and depletion code ORIGEN is used to analyze the core lifetime. Only 1/6 of the reactor core is modeled due to the symmetry. In MCORE code, most isotopes react with neutrons are considered, and the flow chart (Team, 2003) of MCORE is shown in **Figure 5**.

The neutronics calculation using MCNP5 takes several hours for each case, while the depletion calculation using MCORE code takes 5–7 days for each case on a PC with an Intel i5-3470 CPU.

#### Power Distribution

MCNP5 is used to calculate the power distribution of the core. Axial and radial relative power distribution is shown in **Figures 6** and **7**, where θ is the rotation angle of control drums.

It can be seen from the **Figure 7** that the power reaches a maximum value at the center of the active zone and decreases from the center to the periphery where control drums and reflector locate. When θ = 0°, due to absorption of B4C absorber, the power density near the six control drums is lower than the core center. Reflector has influence on radial power distribution, due to the neutron reflection effect. When θ = 180°, power density near the reflector is relatively higher due to the reflection effect of BeO. The radial power peak factor obtained is 1.22. Axial normalized power density of the hottest, the average, and the coldest channel are shown in **Figure 6**. The axial power peak factor is 1.16. The power peak factor of the core (1.42, the product of axial and radial power peak factor) meets the design standard.

#### Reactivity Coefficient

The reactivity feedback is considered an important factor for reactor safety. Thus the reactivity temperature coefficient (αT), the change in reactivity with respect to temperature, needs to be negative. The impact of Doppler broadening effect and the density change due to thermal expansion are considered in the calculation. The thermal expansion coefficient of UN fuel is considered by the equation as follows:

$$
\alpha = 7.096 \times 10^{-6} + 1.4093 \times 10^{-9} \text{T} \tag{2}
$$

The thermal expansion coefficient of B4C on six control drums is considered by the equation as follows:

$$\alpha = 3.016 \times 10^{-6} + 4.3 \times 10^{-9} \text{T} - 9.18 \times 10^{-13} \text{T}^2 \tag{3}$$

The thermal expansion coefficient of BeO reflector is considered by the equation and data given by Kozlovskii and Stankus (2014).

Figure 8 | Reactivity varied with temperature.

The reactivity change with the core temperature varying from cold condition (300 K) to thermal condition (1,200 K) is shown in **Figure 8**. Reactivity decreases as the temperature increases. Reactivity is mainly caused by neutron leakage for thermal expansion of the reflector. By linear fitting, the reactivity temperature coefficient of the designed reactor core is −1.513 pcm/K.

#### Core Lifetime

A code MCORE coupling MCNP and ORIGEN is used to calculate the core life of the designed core. The thermal power is defined 2.4 MWt in the calculation. 613 depletion zones are considered and the time step is defined 0.25a. The effective multiplication factor change during the core life is shown in **Figure 9**. The lifetime of the heat pipe cooled reactor core is 14 years. Due to the depletion of burnable poison and fast neutron breeding effect, the effective multiplication factor increases at the beginning of life and then deceases when burnable poison run out until the end of life.

# THERMAL CALCULATION

A one-dimensional FORTRAN code is used in thermal design and analysis of the micronuclear power source. As **Figure 10** shows, a channel in active zone is defined as an area including a heat pipe and part of the surrounding fuels, air gaps, and matrix. One-dimensional calculation model for each channel is shown in **Figure 10**.

The temperature of the active zone channels is shown in **Figure 11**. Temperature distribution is consistent with the power density distribution obtained in Section "Power distribution." Temperature in the air gap grows significantly from the matrix to the fuel surface. The maximum temperature of fuels is 2,203 K. The maximum power of heat pipes is 79.1 kWt.

The temperature of thermocouple unit is also obtained by the same code. As shown in **Figure 12**, temperature of thermocouple unit varies from 924 to 1,240 K.

# CONCLUSION

In this paper, a core of a heat pipe cooled reactor power source applied for underwater vehicle is designed. The core featured with 2.4 MWt thermal power, lifetime of more than 10 years, and low noise level basically fulfills the application demands. Six control drums with B4C neutron poison have the capability of controlling the reactivity and keeping the reactor safe. The main design parameters are as follows:


# REFERENCES


This work provides a reference to the design and application of the micronuclear power source.

# AUTHOR CONTRIBUTIONS

HS: neutronic calculation; CW: system design; XL: thermal analysis; WT: system design revise; SQ: calculation results review; GS: paper writing guidance.

# ACKNOWLEDGMENTS

This work is carried out under the financial support of the China National Postdoctoral Program for Innovative Talents (Grant No. BX201600124), China Postdoctoral Science Foundation (Grant No. 2016M600796), and Key Supported Discipline of Guizhou Provence [QianXuewei He ZiZDXK(2016)24], 2011 Collaborative Innovation Center of Guizhou Province [Qian Jiao he xietongchuangxinzi (2016)02].


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Sun, Wang, Liu, Tian, Qiu and Su. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Review on Core Degradation and Material Migration Research in Lightwater Reactors

*Jun Wang1 \*, Chen Wang2 , Kaiyi Shi1 and Guanghui Su3 \**

*1College of Engineering, The University of Wisconsin-Madison, Madison, WI, United States, 2Department of Engineering Physics, Tsinghua University, Beijing, China, 3 Institute of Nuclear Science and Technology, Xi'an Jiaotong University, Xi'an, China*

Core degradation and material migration research is one of the key areas in severe accident research. A review of core degradation and melting materials immigration research in light-water reactor is important to encourage relevant research. In this paper, both relevant experiments and numerical analyses are reviewed. Due to their high cost, there have only been a few experiments on severe accidents performed. They focus on different aspects of core degradation and material migration, including early stage and late stage accidents in both PWRs and BWRs. All the current experimental data should be fully utilized due to the limited data available. On the other hand, the data available from numerical analyses for severe accidents is very extensive. There are already many systemic severe accident codes developed by different organizations. These codes provide severe accident sequences analysis, severe accident prediction, and are also important in security policy formulation. At the end of this paper, relevant work at Xi'an Jiaotong University is introduced.

Keywords: review, core degradation, melting materials immigration, experiment, numerical analysis

# INTRODUCTION

Since the inception of industrial nuclear power production there have been several severe accidents that have brought drastic damage to local populations and slowed down the advancement of the nuclear industry all over the world. As a consequence, it is very useful to study severe reactor accidents in order to gain a greater understanding of the processes and mechanisms for accidents, assist with relevant nuclear security policy, ensure reactor security, and to protect the public from radiation damage (Bromet, 2012). Core degradation and material migration research is one of the key areas in severe accident research (Hofmann, 1999). It is a multi-phase, multi-component complex physical chemical process, which is hard to accurately predict (Wang et al., 2014a). Due to the huge complexity and uncertainty, severe accident analysis faces a lot of difficulties (Soffer et al., 1995). Thus, a review of core degradation and melting material migration research in light-water reactors can be very helpful to encourage relevant research. A sample of core melting is shown in **Figure 1** (Zhang et al., 2015a,b).

In this paper, three nuclear severe accidents are overviewed: Fukushima, Chernobyl, and Three Mile Island (Rees, 2009; Morino et al., 2011; Petryna, 2013). These accidents caused a large amount of damage to local populations and remind people to respect the severity of their occurrences (Matzke, 1982; Simmons, 2013). Details of several core degradation and material migration experiments are then described (Hofmann et al., 1997; Schwarz et al., 1999). Some of the experiments focus directly on core degradation and material migration phenomenon (Hagen et al., 1996; Hofmann

#### *Edited by:*

*Arunkumar Nayak, Bhabha Atomic Research Centre, India*

#### *Reviewed by:*

*Han Zhang, Karlsruher Institut für Technologie, Germany Zeyong Wang, Rensselaer Polytechnic Institute, United States*

*\*Correspondence:*

*Jun Wang jwang564@wisc.edu; Guanghui Su ghsu@mail.xjtu.edu.cn*

#### *Specialty section:*

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

*Received: 10 October 2017 Accepted: 24 January 2018 Published: 27 February 2018*

#### *Citation:*

*Wang J, Wang C, Shi K and Su G (2018) Review on Core Degradation and Material Migration Research in Light-Water Reactors. Front. Energy Res. 6:3. doi: 10.3389/fenrg.2018.00003*

et al., 1997; Repetto et al., 2003; Van Dorsselaere et al., 2006), while other experiments deal with separate effects testing of core degradation and material migration phenomenon (Thomsen, 1998; Steinbrück et al., 2010). Next, this paper comes to the numerical investigation and code development of core degradation and material migration (Vierow et al., 2004; Van Dorsselaere et al., 2009). The physical models are reviewed and compared with provide references for relevant research. Finally, in Section "Severe Accident Review at XJTU," the core degradation and material migration research at Xi'an Jiaotong University (XJTU) is overviewed (Sugiyama et al., 2005; Su et al., 2006; Zhang et al., 2011; Chen et al., 2013).

# CORE DEGRADATION ACCIDENT REVIEWS AND PRIMARY RESEARCH

In the history of the nuclear industry, several severe accidents have brought drastic damage to local populations and slowed down the development of nuclear industry all over the world. The accidents at TMI-2, Chernobyl, and Fukushima are the foundation of current severe accident research and provide a valuable insight into probable accident scenarios as well as effectiveness of currently implemented safety systems.

# TMI Accident

On March 28, 1979, a core degradation accident occured in TMI-2 reactor near Harrisburg, Pennsylvania where half of the core melted (Friedman, 2011). This is the first accident that let people realize the speed of a core melt accident and the effectiveness of the implemented safety systems.

The accident began with of the loss of feed water to the steam generators which caused drying out of the steam generator's secondary side for approximately 10–15 min (Sehgal, 2012). The in-vessel pressure kept increasing, while the turbine stopped running and the reactor shut down. However, the pilot operated relief valve, used to bring down the vessel's pressure, did not close as it was designed. This caused the coolant to continue discharging. The high-pressure emergency-core cooling system (ECCS) was also intentionally closed by the operators, leading to the uncovering of the core for 130 min. The cladding temperature kept rising and oxidizing with steam (Moore et al., 1989). Finally, about 50% of the core melted down and the blockage formed near the bottom. If the operators had not restarted the cavitating pumps and filled the vessel with water, the core would have probably melted through the vessel and reacted with concrete in the containment, making the TMI-2 accident even more severe (Toth et al., 1986).

# Chernobyl Disaster

The Chernobyl disaster occurred in one of the four Reaktor Bolshoy Moshchnosty Kanalny (RBMK) reactors in Ukraine, Soviet Union on April 26, 1986 (Howard, 2008). This accident originated from a technology experiment in Chernobyl nuclear power plant which tries to test the available length of a spinning turbine of the electrical power to certain nuclear power plant system (Howieson and Snell, 1987; Sehgal, 2012). Normally, the power would have been reduced to about 30% to conduct the test, but the operators made a mistake and let the power fell to 1%, which is too low to conduct the test. Meanwhile, the absorption effects from Xenon and liquid water caused the power to continue to decrease. In order to complete the test, operators pulled out almost all the control rods and finally raised the power back to 7%. When the test began, the turbine was disconnected, so the water in the core moved more slowly and began to boil. The reactor power rose sharply due to the positive reactivity feedback of the RBMK reactor. When the operator tried to drive the emergency rods and shut down the reactor, it was too late. The power of the reactor increased up to about 100 times of full-loaded power and greatly destroyed the reactor. The radioactive pollution directly went into atmosphere due to the lack of containment, which caused a negative impact on large areas (Howieson and Snell, 1987; Friedman, 2011).

### Fukushima Accident

On March 11, 2011, the Fukushima accident took place in Japan, which was caused by the unexpected natural disasters of the earthquake and tsunamis. The earthquake caused the loss of the reactors' offsite power and initiated the start of the backup diesel generators as designed. However, the subsequent tsunami disabled the diesel generators, and most of the alternating current power in units 1–4 was lost (Holt et al., 2012; Iida, 2013). Thus, the core cooling capabilities were lost and the water levels in the reactor vessel dropped. Large amounts of hydrogen were generated due to the reaction between steam and the zirconium fuel cladding. This caused a large hydrogen explosion in unit 1 on March 20, causing widespread damage. As reported by TEPCO, all fuel in unit 1 melted to the bottom and most of the fuel even went into the primary containment vessel (PCV) (Friedman, 2011; Holt et al., 2012; Aoki and Rothwell, 2013). Finally, the accident was brought under control as the seawater was injected and offsite power was restored.

# EXPERIMENTAL RESEARCH ON CORE DEGRADATION AND MELTING MATERIAL MIGRATION

Due to the damaging effects of core degradation and material migration, it is very important to perform experimental research in this area. However, these experiments are very expensive and funding limits the number of experiments that can be performed. Therefore, we need to take full advantage of all data collected in these experiments.

CORA core degradation early stage experiments started in the 1980s by Germany's KFK National Lab aiming at International Standard Problem 31. The CORA experiment focused on PWR, BWR, WER core degradation and quench situation research (Schanz et al., 1992; Hofmann et al., 1997). It includes 19 tests; CORA-13 is represented here as a reference case (Firnhaber et al., 1993a). CORA-13 started November 15, 1990 at KFK. The testing parameters of CORA-13 include experimental boundary conditions, fuel temperature, hydrogen generation rate, postexperiment fuel morphology, and so on. Fuel morphology is recorded by pictures, which are not convenient for quantitative analysis. The primary and boundary conditions are controlled by the experiment input, so there is no need for extra tests for primary and boundary conditions. This operation method allows for better data in each test (Firnhaber et al., 1993b). The excessive high temperature in severe accident conditions, normally higher than 1,500 K, will bring severe damages to the reactor core, such as the interaction between different chemical materials, melting of core materials, relocation, core flow channel blockage, hydrogen generation, etc. In the early stage of severe accidents, the temperature of core is not extremely high. However, to alleviate this situation before too late, it is very important to learn about the core degradation early stage knowledge and make the probabilistic forecasting of severe accident sequence. The CORA experiments are very useful in helping us understand severe accident early stage phenomena and it also promotes the development of other severe accident early stage research (Haste et al., 2015).

The power burst facility–severe fuel damage (PBF–SFD) tests were performed at the Idaho National Engineering Laboratory, during 1982 and 1985. The objectives of these four in-pile experiments were to investigate fuel rod behavior, hydrogen generation, and the behavior of fission products in severe accidents (Knipe et al., 1986; Martinson, 1989; Petti et al., 1989; Hobbins et al., 1991). The primary components of the PBF reactor are a fission drive core, a central flux trap and an independent pressurized water coolant loop (IAEA TECDOC, 2011). The fuel bundle consisted of 32, 0.9-m long, trace-irradiated fuel rods (Knipe et al., 1986). The test process includes power calibration measurements, high-power operation, shutdown, low-power operation, boildown of the coolants, a short period of high temperature, shut down of the power, bundle storage, as well as the examination of the post-test bundle (IAEA TECDOC, 2011). The PBF–SFD in-pile fuel damage experiments were the first tests performed and they have provided most of the analytical results about the degrading of major phenomena.

The PHEBUS experiment aims to solve International Standard Problem 46 which are the core degradation and melting materials migration problems (Gonnier et al., 1992; Schwarz et al., 1999; Clément and Zeyen, 2013; Haste et al., 2015). The research target of PHEBUS is to study the key in-vessel physical phenomenon, including core degradation, fossil fuel transportation, transformation, and relevant physical chemical phenomenon. The PHEBUS-FP test matrix is shown is **Table 1**. PHEBUS FPT-1 tested a loss of coolant accident (LOCA) in high-density steam condition. This experiment provided scientists a good opportunity to understand the mechanism of core degradation and melting materials migration. To achieve better experimental results, the design of PHEBUS was an exact copy of a light-water reactor. The PHEBUS facility includes a core, steam generator, security vessel,


Table 1 | PHEBUS-FP test matrix (Clément et al., 2003).

and so on. There are 21 bundles in the core, with 1 control bundle in the center and 20 surrounding fuel rods. The control rod is made up of Ag-In-Cd, while the fuel rods of uranium dioxide. These experiments generated results on fission product behavior, and provided quantifiable information about the progression of the core melts in vessel (Clément et al., 2003).

PNL conducted the NRU-full-length high temperature (FLHT) experiment at the NRU reactor of Atomic Energy of Canada Ltd. by the Coolant Boil-away and Damage Progression program (CBDP) (Lanning et al., 1988; Lombardo et al., 1988). The objective of the CBDP program was to understand the impacts of the boiling-away of coolants and get data from the damage progress of the core at LWR conditions, as well investigate the severe accidents of full-length bundles. The NRU-FLHT hardware for testing includes one NRU reactor and other four parts, namely the steam closure cave, effluent control module, test train assembly, and a data acquisition and control system (IAEA TECDOC, 2011). The FLHT tests made data available on SFD behavior.

To explore the core components of the melt progression behavior of boiling water reactors, the annular core research reactor–damaged fuel (DF) tests were conducted between 1982 and 1989 at Sandia National Laboratory (SNL) (Gauntt et al., 1989; Gasser et al., 1990). The DF test bundles were produced from 9 to 14 half-meter long fresh UO2 fuel rods (IAEA TECDOC, 2011). According to the experiment, some significant data were obtained, including the response of the test bundle components, the oxidation of the Zr cladding and canister, and the output of H2 from metal oxidation.

In addition to the facilities already mentioned, there are many other core degradation and melting material migration integral facilities, such as CODEX, which aims at investigating the core degradation within the light-water reactors (Hózer et al., 2000; Hózer, 2002; Hozer et al., 2003); LOFT, which aims at investigating the PWR core behavior during LOCA-type sequences (Jensen et al., 1989; Cronenberg, 1992); SANDIA-XR, which aims to determine the conditions under which steady lower core blockages are formed and those where they cannot be formed; SCARABEE, which aims to solve fast reactor security analysis and the behavior of fuel pool caused by a sub-assembly melting at full power; QUENCH, which is to explicitly investigate the effect of re-flooding on bundle degradation (Sepold et al., 2007, 2009; Stuckert et al., 2010, 2011; Stuckert and Steinbrück, 2014); with the tool FARO, researchers are able to conduct large-scale experiments in order to gain better knowledge of things including structure integration, coolant or molten core with less uncertainties considering relocation and melt progression (Hohmann et al., 1987; Magallon and Huhtiniemi, 2001); and KROTOS, aiming at exploring the problem of steam explosion (Magallon et al., 1996; Huhtiniemi et al., 1997; Annunziato et al., 1999; Huhtiniemi and Magallon, 2001).

### CODE DEVELOPMENT REVIEW AND RELEVANT APPLICATIONS

The development of severe accident analysis codes (including core degradation and melting materials migration phenomenon) and numerical simulation are very important parts of severe-core melt research. This work can significantly reduce experiment research costs, cut down research time, and provide a more guided approach to obtaining research results. Several major accident analysis codes are described below along with their capabilities.

The Electric Power Research Institute develops the Modular Accident Analysis Program (MAAP) (Gabor and Henry, 1983; Kenton and Henry, 1983; Plys et al., 1993; Henry et al., 1994). It can perform fast-running full simulations of severe accidents of the light-water and heavy-water reactor. In the aspect of core degradation and melting materials migration calculations, MAAP mixes outside film flow and inner turbulent flow to calculate the mass of melting materials. When there is no contact between rods the outside film flow model is used, otherwise the inner turbulent flow model is used. Regardless of the flow regime, there is crust formation between the outside wall and flow materials. The crust formed will remain in the computational nodes, but other melted material is still allowed to pass by it and run down stream. This is what will remain in the receiving node. Due to the reduced conservative equations and the simplified discretized systems, MAAP runs much faster than other codes, but still provides credible results, thus it can be used in probabilistic safety analysis for existing reactors as well as more advanced light-water reactors (Petoukhov et al., 2009; Rychkov and Kawahara, 2015).

Methods for estimation of leakages and consequences of release (MELCOR) is developed at SNLs. It is a computer code modeling light-water severe accident progress integrally (Gauntt et al., 1998; Vierow et al., 2004; Wang et al., 2015b). It is a secondgeneration means to assess plant risks and it is considered to be a substitute of the source term code package (Carbajo, 1993; Ashbaugh et al., 2008). MELCOR can also be used to analyze design-basis accidents for advanced plant application (Tills et al., 2009). To analyze melting materials, MELCOR uses a candling flow model. Candling flow model means melting materials move down and are then relocated in different positions like a candle. Candling models are based on thermal/flow and mechanical theory. In MELCOR's candling model is shown in **Figure 2** (Gauntt et al., 1998), the melting materials migrate from top to bottom *via* gravity, until they are condensed on a model component, or stopped by another blockage. Those condensing materials are called condensed debris and they become part of the model's solid components, which are different from particulate debris.

SCDAP/Reactor Excursion and Leak Analysis Program (RELAP) is designed by Idaho National Laboratories to foresee the behavior of reactor systems in normal or severe accident situations (Allison et al., 1983, 1992; Birchley and Stuckert, 2011). LIQuefaction-flow-SOLidification (LIQSOL) SCDAP/RELAP is used to calculate the clad deformation, the oxidation and heat transfer during liquid clad migration (Mladin et al., 2009). The melting of fuel rods has a great effect on core damage. In some cases, the melting cladding may fall into a low-temperature environment, making the oxidation process slower. In other cases, the blockage in the flow channels will also slow down the oxidation process. Three calculation steps are involved in LIQSOL's model, first, calculating the fuel melting rate around the clad; second, calculating the time it takes melting fuel and clad to leave the oxidation level; third, calculating the geometry of the melting materials migration and oxidation generation.

Xi'an Jiaotong University developed Modular In-vessel Degradation Analysis Code (MIDAC) to simulate the process of

in-vessel severe accidents (Wang et al., 2014b,c; Hu et al., 2015). MIDAC includes five modules for calculations: early-behavior, core degradation, debris bed, molten-materials-in-vessel retention (IVR), and connection modules. Those modules can be used not only to demonstrate the entire severe accident, but also to serve as a link for the conjunction of multiple severe accident codes. Furthermore, the primary and transient thermal-hydraulic system of MIDAC was used to compare with the similar calculations obtained from the module of SCDAP/ RELAP5 and showed good comparison results. MIDAC will play a large role in many Chinese nuclear reactors severe accident analyses.

KESS [Institut Fur Kernenergetik Und Energiesysteme (IKE) modular program system to simulate and analysis core melt accident] is developed by IKE, Germany (Schatz and Hocke, 1995). Its intention of designing this program is to easily simulate the physiochemical processes that occur in the reactor core, during heat up, and the degradation process. It can simulate large amounts of phenomena including: heating up of structures, oxidizing of the cladding, production of hydrogen, expanding and collapsing of mechanical rods, dissolving of UO2, melting and relocating of materials, and releasing and transporting of fission products.

Severe accident analysis code with mechanistic, parallelized simulations oriented toward nuclear fields is an integral analysis code targeting severe accidents included in the IMPACT project (Ujita et al., 1999, 2002; Naitoh et al., 2000). It includes a mechanistic model that indicates many different phenomena, such as the reactor scram and the PCV damage with elaborated mathematical characteristics, and can explain physical behaviors for analytic results.

ICARE/CATHARE is developed at the French Institute for Nuclear Protection and Safety (IPSN), to evaluate severe accidents comprehensively under primary systems (Chatelard et al., 2006; Seiler et al., 2008). The ICARE/CATHARE code is designed to calculate, in a mechanistic way, reactor core damage and primary circuit behavior in PWRs. ICARE/CATHARE features a comprehensive set of models for late degradation allowing the model to follow the materials from their early melting in the core region to their later relocation.

Analysis of thermal-hydraulics of leaks and transients with core degradation is developed by Gasellschaft fur Anlagen und Reaktorsicherheit (GRS) in cooperation with IKE in Germany (Trambauer and Austregesilo, 2003; Austregesilo et al., 2007; Hollands et al., 2007; Repetto et al., 2007). The ATHLET code can be used to predict a spectrum of beyond design-basis accidents for BWRs and PWRs, but without core degradation. The CD part adds functions describing core melting, fission product release, and transport processes in the primary system. The code structure is highly modular so that many models can be included and further developments and changes to the models can be made easier by the users.

Thermal-hydraulic analysis of loss-of-coolant, emergency-core cooling and severe-core damage designed by Japanese Atomic Energy Research Institute, aims to predict the progression of the core in severe accident conditions of LWRs (Abe et al., 1986; Abe, 1990; Hashimoto and Soda, 1991; Ishikawa et al., 2002). The code was developed for Level 2 PSA and can be used in many assumed accidents. It has been applied to serial analyses of experiments, severe accident sequences, and evaluation of designs.

The code ASTEC is a joint research product of GRS (Germany) and IRSN (France) to simulate severe accidents (Van Dorsselaere et al., 2009; Coindreau et al., 2010; Dorsselaere et al., 2012; Chatelard et al., 2014). It is a relatively new and integrated code designed to predict all the behaviors during severe accidents in LWRs, such as any event in the beginning and any event with possible radiation related to a certain fission product. ASTEC has also been widely used in many applications like European pressurized reactor simulating. **Table 2** summarizes the codes and their features.

### SEVERE ACCIDENT REVIEW AT XJTU

Xi'an Jiaotong University is doing a lot of research on core degradation and melting materials migration. The oxidation of the metal core was investigated by Dr. Guanghui Su (Sugiyama et al., 2005). During his research from 2004 to 2006, he proposed an important theory, the "hot spot" during core degradation process (Su et al., 2006). Hot spot is the point which has the highest temperature in the core, which affects the performance of reactor pressure vessel.

Under severe accident conditions, clad oxidation cannot be ignored (Beuzet et al., 2011) since it always weakens the cladding material (Schanz et al., 2004). With the oxidation of cladding material, the accident process can be accelerated, and the integrity of vessel may be threatened by hydrogen generation (Volchek et al., 2004). Dr. Jun Wang made a comparison between CORA and MELCOR steam-oxidation results and showed that the hydrogen generation rate is not predicted as well as expected (Wang et al., 2014a,b). The oxidation model of MELCOR seems to lack mechanical changes during the quenching process which are needed to accurately predict oxidation results. Another model

Table 2 | Core degradation simulation codes.


developed by Mr. Keyou Mao named Clad solid-phase Oxidation Analysis Code (COAC) has shown better results through verification with a test from the CORA experiments, specifically CORA-13. This code COAC is also being used in AP1000 for further clad-oxidation calculations (Mao et al., 2015).

Dr. Ronghua Chen developed a solidified TEXAS-VI code model and applied this code into FARO L14 analysis (Chen et al., 2013). In Chen's work, a model was put forward that studies molten fuel breakup and solidified impacts for the TEXAS-VI code (Chen et al., 2012). The model concentrates on the thermal stress and the solid crust layer effects on fuel particles or fragments. In this model of solidification, the Fourier heat equation was used to study the temporary temperatures and the thickness of fuel particles' crust layers in boundary and initial circumstances. TEXAS-VI was compared with the FARO L14 experiment (Nilsuwankosit et al., 1996) to validate quench and mixed fuel coolants data. Results show that the FARO L14 pressure history, the liquid-water-pool temperature, and the vapor temperature were accordance with the revised model simulation results. FCI explosion energetics were significantly influenced by the mixing performance (Grishchenko et al., 2014). Researchers are still furthering investigating the solidified impacts to study energetics.

It is important to predict the probability, behavior, and the influence of core degradation to calculate the risk and corresponding mitigation (Wang et al., 2014c). However, even until now, all the models of severe accidents are not sufficient and lack precision (Sehgal, 2012). Thus, it is critical to investigate the core degradation mechanisms and to develop solutions for such accidents. Dr. Jun Wang made a core degradation and melting materials migration model that analyzed PHEBUS FPT-1 in MELCOR (Wang et al., 2015a). Through this work, the core degradation parameters, such as pressure, temperature, hydrogen generation, and mass distribution were analyzed. In addition, several parameters are shown by color maps to visualize the core degradation and melting materials migration throughout the accident process.

Dr. Y. P. Zhang is studying a possible solution to critical severecore accidents, named the IVR of core melts (Zhang et al., 2011). Another creative management process named external reactor vessel cooling for IVR analysis was recommended to be used in predicting the reactor cavity flooding and the depressurization process related to severe accidents to determine safety margins of IVR in the AP600. IVRASA proved to be applicable and accurate according to the results produced by UCSB using the FIBS benchmarks. It was also found that the thermal reactions of a couple molten configurations could be predicted.

In addition to the current work of those above, many others are currently pursuing other avenues of severe accident analysis including Dr. Xiaoli Wu's research on loss of pool-cooling accidents in PWRs by applying MAAP5 (Wu et al., 2014). This research involved two cases where the initial water levels differed. Dr. Wei Li also used MAAP5 to analyze the lower head of PWR RPV small break loss of coolant accident scenarios based on high-pressure injection system failures (Li et al., 2014). Dr. Luteng Zhang has performed a MAAP5 estimation of the PWR severe accident setting the pressurizer safety valve stuck-open as the initiating event (Zhang et al., 2015a,b). Results provide a clarified explanation of the process and such events as coolant release, ECCS operation, core exposure, etc. Dr. Liang Hu has also investigated the severe accident scenario of PWR reaction to LOCA along with SBO (Hu et al., 2015).

Xi'an Jiaotong University developed MIDAC for analysis of IVSA to satisfy the domestic requirement for independent scheduling of software (Wang et al., 2014c). Compared with SCDAP/ RELAP5, MIDAC demonstrated a high accuracy in the primary system thermal-hydraulic transient analysis in the CPR1000 station blackout, and the comparison results proved the validity of MIDAC code (Wang et al., 2014c).

In this paper, we reviewed core degradation and melting materials migration research in light-water reactors. The content contains experimental and numerical analyses as well as an overview of current work in the field. Based on analysis and discussion, we can come to the following conclusions:


### REFERENCES


well as being useful for predicating the processes in core melt accidents.

3. The research at XJTU is abundant and can be used for further Chinese research.

# AUTHOR CONTRIBUTIONS

The work is mainly done by JW. CW and KS both contributed to the accomplishment of this work, and GS is the advisor of JW.

### ACKNOWLEDGMENTS

Thanks for Mr. Robert John Armstrong's writing assistance, technical editing, language editing, and proofreading. This work is also supported by Chinese "Program for Changjiang Scholars and Innovative Research Team in University" (no: IRT1280) and Chinese "National Major Projects" (no: 2011 ZX 06004 008 004). This work is also supported by Key Supported Discipline of Guizhou Provence (Qian Xuewei He Zi ZDXK[2016]24), 2011 Collaborative Innovation Center of Guizhou Province (Qian Jiao he xietongchuangxin zi [2016]02).


of the primary core degradation mechanism. *Ann. Nucl. Energy* 85, 193–204. doi:10.1016/j.anucene.2015.05.015


Zhang, Y. P., Qiu, S., Su, G., and Tian, W. (2011). A simple novel analysis procedure for IVR calculation in core-molten severe accident. *Nucl. Eng. Des.* 241, 4634–4642. doi:10.1016/j.nucengdes.2011.03.055

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Wang, Wang, Shi and Su. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Accident Process and Core Thermal Response During a Station Blackout Initiated Study for Small Modular Reactor

Shasha Yin<sup>1</sup> \*, Yajing Tian<sup>1</sup> , Suizheng Qiu<sup>2</sup> , Ye Tian<sup>1</sup> , Huawei Fang<sup>1</sup> , Wei Huang<sup>1</sup> and Zhihui Chen<sup>1</sup>

*<sup>1</sup> Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, China, <sup>2</sup> School of Nuclear Science and Technology, Xi'an Jiaotong University Xi'an, Shanxi, China*

Interest in evaluation of reactor response and off-site consequences following beyond design basis accidents has significantly increased after Fukushima. Considering the inherent system safety of nuclear plants, experts begin to refocus on the development of small and medium nuclear reactors (SMRs). Based on proven and widely utilized traditional pressurized water reactor (PWR) technology, there are many kinds of small modular reactor design concepts, which are actively being developed in the USA, China, Japan, Russia, Korea, and other countries. But due to the significant differences between the traditional distributed arrangement and integrated arrangement of SMRs, the variation trend of key parameters may be different during normal operation or accident process. Hence, in this paper, we simulated the small modular reactor severe accident by MELCOR. This paper summarizes the core thermal hydraulic response for a hypothetical severe accident caused by station blackout at a small modular reactor using MELCOR. Analytical results for temperature distribution of the fuel pellets, the fuel cladding, flow rate of the coolant, and hydrogen mass changing over time are presented. The analyses are focused on safety assessment of the reactor core for severe accidents and are a part of the overall evaluation of safety features of the small modular reactor for residual risk posed by severe accidents.

#### Keywords: SMR, core degradation, melting materials immigration, severe accident, MELCOR

#### INTRODUCTION

After several severe accidents, such as in Three Mile Island, Chernobyl, and Fukushima, the safety has become the public focus. After the Fukushima accident, experts begin to refocus on the development of small and medium nuclear reactors (SMRs) (Shropshire, 2011), considering the inherent system safety of nuclear plants. Based on the proven and widely utilized traditional pressurized water reactor (PWR) technology, there are many kinds of small modular reactor design concepts, which are actively being developed in the USA, China, Japan, Russia, Korea, and other countries. The USA has basically completed the design of mPower (Halfinger and Haggerty, 2012), NuScale (Johnson et al., 2009; Ingersoll et al., 2014) and W-SMR (Fetterman et al., 2011). Korea has finished the design of SMART (Chang et al., 2000; Bae et al., 2001; Kim et al., 2011) and analyzed

#### Edited by:

*Jun Wang, University of Wisconsin-Madison, United States*

#### Reviewed by:

*Hui Cheng, City University of Hong Kong, Hong Kong Claudio Tenreiro, University of Talca, Chile*

> \*Correspondence: *Shasha Yin ysszjkong@163.com*

#### Specialty section:

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

Received: *23 March 2018* Accepted: *27 April 2018* Published: *22 May 2018*

#### Citation:

*Yin S, Tian Y, Qiu S, Tian Y, Fang H, Huang W and Chen Z (2018) Accident Process and Core Thermal Response During a Station Blackout Initiated Study for Small Modular Reactor. Front. Energy Res. 6:43. doi: 10.3389/fenrg.2018.00043* the system safety. The international co-operations also presented the concept of IRIS (Ricotti et al., 2002; Carelli, 2003; Carelli et al., 2004); meanwhile, the structure designs in other countries are ongoing as well. China has also carried out the related research and development (R&D) work. Accordingly, in this paper, we choose the SMR as the research object.

Before focusing on SMRs, we should strictly evaluate the safety of SMRs first, before its commercialization (Yin et al., 2014). Especially after the Fukushima accident, attention has been paid to the severe accident prevention and mitigation measures. There are many differences during the severe accident process between the small modular reactor and the traditional large-scale PWR because of the differences in structure and design issues for both reactors. Therefore, focusing on severe accident of small modular reactor and understanding the severe accident process and mechanism are the necessary foundation for the development of small modular reactor severe accident prevention and mitigation measures.

Based on the traditional PWR technology, core degradation, melting materials migration and clad oxidation are all priority research questions (Wang et al., 2014). Currently, SMRs mainly adopt a set of passive safety devices for mitigating the effects of severe accidents (Yin et al., 2016). However, due to the significant differences between the traditionally distributed arrangement and integrated arrangement of SMRs, the variation trend of key parameters may be different during normal operation or accident process. Hence, this article specially focusses the station blackout (SBO) accident for SMRs. Simulated the cladding and fuel pellets failure mode for radial and axial sections in the process of core melt. We calculated the SBO accident; analytical results of temperature distribution of the fuel pellets, the fuel cladding, flow rate of the coolant, and hydrogen mass changing over time are presented. The analyses are focused on safety assessment of the reactor core for severe accidents and are a part of the overall evaluation of safety features of the small modular reactor for residual risk posed by severe accidents.

# ANALYSIS MODEL

#### Brief Description of MELCOR Code

MELCOR is a fully integrated, engineering-level computer code that models the progression of severe accident in light water reactor nuclear power plant (Gauntt et al., 1998). MELCOR is executed in two parts, MELGEN and MELCOR (Gauntt et al., 2000a,b). The relationships between MELGEN and MELCOR are shown in **Figure 1**. MELGEN is also composed of a number of different packages, such as COR, CVH, and FL., and each of the models represents a different portion of the accident phenomenology or program control. The major inputs are specified in MELGEN, which are also processed and checked for the parameters. MELCOR is mainly responsible for transient calculation based on the input of MELGEN.

# Nodalization of China SMR Primary System

The nodalization of the severe accident analysis model for China SMR is shown in **Figure 2**. The major equipment are modeled as control volumes, which contain energy and mass, such as steam generator, pressurizer, core makeup tank (CMT), accumulator, and so on. Control volumes are collected by flow path module, which delivers energy and mass from one control volume to another with no mass and energy on itself. For China SMR, the riser and downcomer are similar to the hot and cold leg of traditional PWR. Hence, we can use the experience of traditional PWR directly. There are four reactor coolant pumps and 16 stream generators in the vessel. Therefore, these stream generators (SG) should be divided into four groups in order to correspond to the four coolant pumps. For every primary loop, the coolant flow through their respective SGs and coolant pumps and then merge into the same downcomer. Two core makeup tanks are modeled as control volumes and divided into four pipelines as the core makeup water. The pressure balance pipeline is injected into the middle annular of the reactor pressure vessel (RPV). The trip valves are installed on the pipeline between the CMT, and the RPV is used to turn on the CMT. The accumulator is also divided into four pipelines and controlled by check valves. The accumulator shares the same direct injection pipeline with the CMT, which makes up water for primary loop. Relief valves of stage 1 and stage 2 automatic pressure relief system (ADS) are installed on the pipeline between the top of pressurizer and pressure relief valve. The stage-3ADS is installed on the pressurizer wave pipe. When the relief valves of the stage-3 ADS are opened, primary pressure is fast balanced to containment pressure. The stage-4 ADS is installed on the pressurizer surge pipeline. For this system model, the secondary side of SG is modeled as 3 control volume from top to the bottom of the SG heat exchanger tube. The secondary side of the SG's feedwater is modeled as a mass source which may be prescribed in using control functions. Containment environment is modeled as control volumes, which includes upper ring cavity, lower ring cavity, and several compartments.

### Nodalization of Core

In the Core Package (COR) of MELCOR code, the core region (including the reactor core and lower plenum, two control volume) can be divided into some more detailed cells, as shown in **Figure 3**. The core region is axially divided into 14 nodes and radially divided into 4 rings. The core activity area contains the node 4 to node 3, and node 1 to node 3 belongs to the lower plenum area. Node 14 is the core upper inactive area. The Core Package also has the special definition to the lower head. In this article, we divide the lower head into four looped areas.

**Abbreviations:** ACC, accumulators; ADS, automatic depressurization system; CMT, core makeup tank; IRWIST, in-containment refueling water storage tank; RCS, reactor coolant system; PXS, passive core cooling system; PRHRS, passive residual heat removal heat exchange; PWR, pressurized water reactor; R&D, research and development; SBO, station blackout; SMRs, small and medium nuclear reactors; SG, stream generators.

#### DESCRIPTION OF ACCIDENT AND ANALYSIS CONDITIONS

Before we calculate the SBO accident, the steady state analysis was performed for evaluating the SMR MELCOR model. The main parameters of SMR are shown in the following **Table 1**.

The SBO accident process of SMR is simulated by MELCOR code and the whole process from the accident happened to the corium drop into the cavity is analyzed, especially the

TABLE 1 | The main parameter of SMR.


TABLE 2 | Sequence of SBO initiated severe accident.


passive core cooling system for accident mitigation effect of SMR. Before the accident happened (0 s), steady state operating conditions are obtained by the calculation for a period of time. After the accident has happened, the primary pump runs down, immediately following the reactor scram, and loss of feedwater occurs. Due to the loss of the hot trap in secondary loop, the temperature and pressure of the primary loop increase. When the pressure exceeds the setting value of the pressure relief valve (17 MPa), the pressure relief valves opens. The coolant escapes through the relief valve gradually, so the core begins to expose at the same time, and the cladding temperature rises. When the fuel cladding temperature rises to 1273.15 K, the cladding fails, then radioactive material leaked out. It gradually began to appear the core melt and collapsed. The process of the SBO accident is analyzed, and the characteristics of passive core cooling system (PXS) are also analyzed during the accident process.

#### RESULTS AND DISCUSSIONS

In severe accident analysis, reactor core thermal hydraulic response directly affects the safety of the reactor, and it is critically important in formulating the corresponding severe accident management strategy. The coolant flows through the lower plenum, and then it is divided into two parts in the core. The main part flows through the active core directly, and the other through the core bypass. Once the reactor shuts down, the coolant pump starts running down, and the cooling of core only relies on the natural circulation. When the natural circulation flow rate decreases, the pressure and temperature of the core increases, and when they get the set value, the pressure relief valve will be triggered to open. **Table 2** lists the time sequence of SBO-initiated severe accident for small modular reactor.

As shown in **Figure 4**, the water level remains volatile at 4.3 m about 20,000 s. It means that the effective cooling is provided through the new nature circulation created by CMT and lasts a long time. With the continuous natural circulation, the temperature of CMT rises gradually, and the water level begins to drop. When the valves of ADS open, the coolant and vapor are released from the primary system. The water level of the core starts to fall. The active core region is flooded underwater until 38,540 s(as shown in **Figure 5**), and then the temperature of the core outlet begin to rise.

The coolant mass flow rate of core flow channel and the accumulative coolant quality in core inlet are shown in **Figure 6** (before the transient accident, there are 5,000 s calculation for steady state). Hence, the coolant cooling capacity for core is gradually losing with curve declining in the figure.

In order to accurate the analysis of the process and consequences of severe accident, the present study carried out a detailed research in the distribution of transient temperature for core fuel, especially the period that core began to melt. When studying the transient temperature distribution of core fuel, we also took full consideration of the fission products, hydrogen source term, related mechanical deformation, the reactor core structure, and so on.

When a SBO accident happens without any operator intervention, the transient distribution of core fuel temperature is as shown in **Figure 7**. The number represents the core cell which is according to the axial and radial unit number in MELCOR model, as shown in **Figure 3**. In the beginning, due to the reactor shutdown, the fuel temperature curve accord with the heat load decline from full power to decay heat. The loss of forced circulation of primary loop does not lead to core fuel temperature rise at the start because the rest of feedwater in the secondary side also can be used as a heat sink to maintain core temperature for a period of time. Then at 2,843 s, CMT is triggered (when the pressure of pressurizer is less than a set value), a newly formed natural cycle can maintain the core temperature that does not rise in nearly 10 h. When the water level drops down to the set value (67.5%), the corresponding ADS system is triggered, and it is at about 34,485 s as shown in **Figure 7**. With the opening of ADS system, the water vapor of coolant system is released into the pressurizer relief tank, and therefore, the primary system temperature also has a slight fall. However, as the coolant loses from ADS system, the core begins to uncover since 34,500 s, and the fuel temperature also rises gradually. With the temperature rising, the zirconium cladding interacts with water or steam, and causes strong exothermic reaction. Hence, the related fuel temperature is also surged. As shown in **Figure 7**, melting begins at the top of fuel element, and then with the collapse, the melt area expands unceasingly. Some temperature oscillation in the lower part of core is mainly

due to the oxidation of cladding and the dynamic migration of fuel.

The fuel temperature distribution shown as **Figure 8** is for the right half of the core for a number of discreet periods of interest. At the onset of degradation, the high temperature regions are concentrated on the central higher power region, and then move downwards as exothermic reaction abated. At 43,572 s, the temperature of core cell 112 gets to 1,300 K. At about 47,000 s, the high temperature area is at the central part of the reactor core bottom.

In the process of severe accident, it can produce large amounts of hydrogen. These large amounts of hydrogen will be a serious threat for the integrity of reactor containment, so the production rate of hydrogen in core and the spread of hydrogen in containment need to be focused.

With the core uncovering, the coolant in primary loop gradually evaporates. The heat transfer performance is poor between the fuel and the steam. Therefore, fuel temperature rises rapidly. When the temperature exceeds 1,300 K, the zirconium claddings begin to fail and release a lot of heat, which will aggravate the melt of reactor core. The quality of hydrogen produced by core is shown in **Figure 9**. Between 4,000 and 8,000 s, due to the cladding failure, the interaction of zirconium–water produces large amounts of hydrogen. The peak of hydrogen quality change rate also appears at about 40,000 s. At 80,270 s, due to the failure of the lower head, the hydrogen releases into the containment, and the hydrogen quality in primary loop levels off.

The hydrogen mass change in containment is shown in **Figure 10**. Due to the opening of ADS system, when the hydrogen mass production rate gets to the peak in primary loop, the hydrogen will be released into containment through the ADS pipeline, so the time when the hydrogen mass production rate gets to the peak in the containment and in the core is consistent. Before the failure of the lower head, the hydrogen mass production rate is very small. However, as shown in **Figure 11**, continuing with time, the hydrogen quality will continue to accumulate, and it will be a great threat to the integrity of the containment vessel.

# CONCLUSIONS

This paper comprehensively analyzed the SBO severe accident through modeling a complete small modular reactor. The calculation results show that with the gradual development of accident process, severe accident will happen eventually. The

# REFERENCES


RPV fail at about 80270.0 s, and then inevitably containment failure accident happened. The detailed thermal hydraulic model of core can accurately evaluate the reactor core melting process and predict the radioactive source term. The intense interaction between melt and water in lower plenum caused the pressure vessel failure. The production and accumulation of hydrogen may result in containment failure. Hence, in the subsequent research work, the study may extend to sensitivity analysis of hydrogen, which will be helpful to analyze the severe accident mitigation of small modular reactor.

# AUTHOR CONTRIBUTIONS

SY: Thermal analysis and paper writing; YaT, HF and YeT: Calculation results review; SQ, WH and ZC: Paper writing guidance.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Yin, Tian, Qiu, Tian, Fang, Huang and Chen. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Measurement for surface Tension of aqueous inorganic salt

#### *Jiming Wen1 , Kaiyi Shi <sup>2</sup> , Qiunan Sun1 , Zhongning Sun1 and Haifeng Gu1 \**

*<sup>1</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China, 2Department of Chemistry and Chemical Engineering, LiuPanshui Normal University, LiuPanshui, China*

Bubble columns are effective means of filtration in filtered containment venting systems. Here, the surface tension has a significant influence on bubble size distribution and bubble deformation, which have a strong impact on the behavior of the bubble column. The influence of aqueous inorganic compounds on the surface tension depends on the electrolytic activity, Debye length, entropy of ion hydration, and surface deficiencies or excess. In this work, the surface tensions of same specific aqueous solutions have been measured by different methods including platinum plate method, platinum ring method, and maximum bubble pressure method. The measured surface tensions of both sodium hydroxide and sodium thiosulfate are less than that of water. As solution temperature ranges from 20 to 75°C, the surface tension of 0.5 mol/L sodium hydroxide solution decreases from 71 to 55 mN/m while that of 1 mol/L solution decreases from 60 to 45 mN/m. Similarly during the same temperature range, the surface tension of 0.5 mol/L sodium thiosulfate decreases from 70 to 38 mN/m, and that of 1 mol/L sodium thiosulfate is between 68 and 36 mN/m. The analysis for the influence mechanism of aqueous inorganic on surface tension is provided. In addition, experimental results show that the surface tension of solid aerosol suspension liquid has no obvious difference from that of distilled water.

Keywords: surface tension, aqueous inorganic, solid aerosol suspension liquid, filtered contain venting system, bubble column

# INTRODUCTION

Containment, as the last shield, plays a crucial role in protecting nuclear power plant (NPP), whose integrity defines the success of the safe strategy of NPP in severe accident. However, once a serious accident such as the lost of coolant accident occurred, the pressure in containment would rise persistently. When the containment pressure rose beyond the threshold pressure, rupture or crevasse of some parts in containment could appear. One strategy for solving this problem is to discharge the gas in containment automatically when the pressure of containment is close to the threshold pressure (OECD/NEA/CSNI, 2014). The gas to be vented is composed of not only air and steam but also various kinds of high-dose radioactive material. The representative radioactive materials that can do serious harm to the people and environment near NPP are iodine, aerosol, and methyl iodide. Filter containment venting systems (FCVS) are designed to filter radioactive materials when the gas within containment needs to be vented. After Fukushima Daiichi nuclear disaster, many countries, such as Japan and Korea, have made a specialized schedule of installing FCVS. Owning to the ability of enlarging specific transfer area and prolonging bubble resident time in solution, bubble column

#### *Edited by:*

*Uwe Schröder, Technische Universitat Braunschweig, Germany*

#### *Reviewed by:*

*Daniell Tincher, Virginia Commonwealth University, United States Keyou Mao, Purdue University, United States Yago Rivera, Universidad Politécnica de Valencia (UPV), Spain*

*\*Correspondence: Haifeng Gu guhaifeng@hrbeu.edu.cn*

#### *Specialty section:*

*This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research*

*Received: 12 November 2017 Accepted: 26 February 2018 Published: 12 March 2018*

#### *Citation:*

*Wen J, Shi K, Sun Q, Sun Z and Gu H (2018) Measurement for Surface Tension of Aqueous Inorganic Salt. Front. Energy Res. 6:12. doi: 10.3389/fenrg.2018.00012*

possesses a potential capability of filtering methyl iodide (Smith, 1991; Zhang and Fan, 2003; Gumulya et al., 2016).

The work performance of bubble column has relation with the bubble size distribution on which the surface tension of washing solution has significant effect (Kazakis et al., 2008). Surface tension is one force suppressing bubble detachment, so bubble detachment volume increases with surface tension (Loimer et al., 2004). The interface area between gas and liquid phases is influenced by the bubble deformation, which becomes more serious as surface tension decreases (Celata et al., 2007). It has been proved by many researchers that bubble coalescence probability increases with inorganic ions concentration. When inorganic ions' concentration is above the critical value, bubble coalescence even disappears (Jamialahmadi and Müller-Steinhagen, 1992; Ruzicka et al., 2008; Cachaza et al., 2011). One reason for this phenomenon is that the surface tension of washing solution increases due to the existence of inorganic ions. In practical application, the washing solution in bubbling filtration installation is added with sodium hydroxide and sodium thiosulfate, which are used to consume methyl iodide and iodine persistently (Wen et al., 2017). In addition, an amount of aerosols is suspended in washing solution. These kinds of materials may have influence on the surface tension of solution. The existing study pays major attention on the surface tension of high concentration inorganic solution (Li and Lu, 2001; Dutcher et al., 2010). But the concentration of inorganic ions in the washing solution of filtered containment venting system is below 1 mol/L. In this work, the surface tensions of aqueous solutions containing different materials are measured, and some applicable analyses have been given out.

#### EXPERIMENTAL PROCEDURE

The inorganic reagents are the analytical pure, which are produced by Sigma-Aldrich. As shown in **Table 1**, because the concentrations of sodium hydroxide and sodium thiosulfate in practical application are 0.5 and 0.2 wt%, respectively, which both belong to low concentration inorganic solution, their concentrations range 0–1 mol/L in this work. To avoid the influence of evaporation on the solution concentration, the temperature range of solution chosen as 15–75°C. In this work, barium sulfate and titanium dioxide are employed to simulate the radioactive aerosol, and the concentration is 0–1 g/L (Zhou et al., 2014). During


the experimental process, all glassware are cleaned in chromic acid and then rinsed at least five times to eliminate the organic contamination.

The surface tension is measured by different measurement methods. Most experimental points are conducted by K100 of KRÜSS GmbH. The instrument can use both the platinum plate and platinum ring methods to complete the measurements. The measurement principle of platinum plate method is expressed by Eq. 1. Platinum plate method possesses the advantage of short measurement period. The measurement principle of platinum ring method is expressed as Eq. 2. The measurement period of platinum ring method is obviously longer than that of platinum plate method, but the measured results of platinum ring method possess more perfect repeatability. To avoid the possible instrument error, some experiment points of inorganic solution are conducted back to back by SFZL-U of Innuo. For a few special experimental points, the surface tension is measured by BP100 of KRÜSS GmbH, which employs maximum bubble pressure method. The measurement principle of maximum bubble pressure method is expressed as Eq. 3.

The measure precision of the instrument K100 is tested by measuring the surface tension of distilled water as shown in **Figure 1**. It shows that the measured results of distilled water are consistent basically to the values measured by Abramzon (1994), which are accepted by many researchers. Before each point is conducted, the surface tension of distilled water is measured and compared with the standard value to guarantee no impurities contamination. The temperature can be measured by the resistance temperature detector assembled in the measurement instruction. The measurement precision of thermal resistant is verified by the second order accuracy mercury thermometer as shown in **Figure 2**. At each measurement point, it is pledged that the temperature remains constant for 5 min,

$$
\sigma = \frac{\mathcal{W}\_{\rm til} - \mathcal{W}\_{\rm p1}}{2l},
\tag{1}
$$

where σ is the measured value of surface tension, *W*t1 is the total weight of platinum plate and the solution adhering to the plate, *W*p1 is the net weight of platinum plate, and 2*l* is the wetted perimeter of platinum plate.

$$
\sigma = \frac{\mathcal{W}\_{\rm v2} - \mathcal{W}\_{\rm v2}}{4\pi R},
\tag{2}
$$

where *W*t2 is the total weight of platinum ring and the solution adhering to the ring, *W*r2 is the net weight of platinum ring, and *R* is the diameter of the platinum ring.

$$
\sigma = \frac{r}{2} \Delta P\_{\text{max}},
\tag{3}
$$

where *r* is bubble radius, which is obtained by calculating bubble volume according to the total amount of gas in a bubble and Δ*P*max is the maximum additional pressure caused by surface tension.

#### RESULTS AND DISCUSSION

The influence mechanism of inorganic ions on surface tension is very complicated. According to Weissenborn and Pugh (1996), surface tension decreases with increasing ions concentration when surface excess concentration or adsorption of the solute at the interface between liquid and gas is positive. *Vice versa*, surface tension increases for negative surface excess concentration or adsorption of the solute at the interface. As shown in **Figure 3**, the surface tension of 1 mol/L aqueous sodium thiosulfate solution has been measured six times at 20°C. The time interval between two adjacent measurements is around 15 min. The all measured surface tensions of sodium thiosulfate solution are smaller than the one of distilled water at the same temperature. This result suggests that the cation and anion in aqueous sodium thiosulfate can

platinum plate method (sodium thiosulfate concentration is 1 mol/L; solution temperature is 20°C).

result in surface deficiency. However, the measured surface tensions of sodium thiosulfate solution for the latter three measurements are larger than those for the former three measurements. This is because platinum plate is rinsed and roasted during the former three measurements. The ions adsorbed on platinum result in the ion lost, but the ions' concentration retrieval rate is slow near the solution surface, so the ion concentration at the interface between gas and liquid phases decreases, and the surface tension of solution approaches the value of water. Then, the latter three measurements are conducted to verify this point. During these three measurements, platinum plate is not processed at all. After one measurement, another measurement is conducted with that. As shown in **Figure 3**, the surface tension of aqueous sodium thiosulfate recovers to the initial value gradually. This phenomenon suggests that the ions' concentration reaches the equilibrium value gradually. Due to the adsorption characteristic between sodium thiosulfate and platinum, platinum ring method is employed to measure surface tension, and the measured results are shown in **Figure 4**. Because the contact area between platinum ring and solution is smaller, the influence of adsorption on the interface ion concentration can be neglected. The surface tension measured by platinum ring has good repeatability. Therefore, the surface tension of other solutions is measured with platinum ring.

**Figure 5** shows the measured results for 1 mol/L sodium hydroxide solution at 35°C. Because the measurement period of platinum ring is long, the time interval between adjacent measurements is 30–40 min. The surface tension of 20°C sodium hydroxide solution is 59.36 mN/m, but the former two measured surface tensions of 35°C solution are larger than that of 20°C solution as first and second measured results in **Figure 5**. This phenomenon goes against the principle of entropy increase. To verify whether this anomaly phenomenon is caused by contamination, the surface of solution is scraped to eliminate the impurity resulting possibly from air. Then, three more

Figure 4 | The history of measured surface tension of sodium thiosulfate by platinum ring method (sodium thiosulfate concentration is 1 mol/L; solution temperature is 20°C).

platinum ring method (sodium hydroxide concentration is 1 mol/L; solution temperature is 35°C).

measurements (third–fifth in **Figure 5**) are conducted, but the surface tension increases drastically. Then, the aqueous sodium hydroxide is placed in clean environment for 8 h. Three more measurements (sixth–eighth in **Figure 5**) are conducted, and the measured results are consistent to the principle of entropy increase compared with the value at 20°C. Therefore, the surface tension variance is not caused by the contamination but the standing time of aqueous sodium hydroxide. To verify this judgment, another set of experiments are conducted. Maximum bubble pressure method is adapted because this method is less

affected by impurity. In **Figure 6**, the surface tension of 1 mol/L sodium hydroxide is measured at 20 and 35°C. The bubble formation time can be altered from 0.05 to 200 s by adjusting gas flow rate. It can be seen that the measured surface tension by maximum bubble pressure method maintains almost constant. The maximum bubble formation time 200 s is enough compared with the enrichment time of active impurity (Kretzschmar and Miller, 1991). So the phenomenon in **Figure 5** is not caused by the contamination of active agent. It can be concluded that the surface tension decreases resulted from sodium hydroxide is self characteristic of solution. According to Meng (2011), hydroxyl ion clusters has significant influence on the molecular structure of water. The binging energy of hydrogen bond clusters is higher than the one of hydroxyl ion clusters. Because of hydronium existence, more hydroxyl ion clusters form, and the surface tension decreases. As the solution temperature is room temperature, the thermal motion of molecule is relatively weak, so the formation time of hydroxyl ion clusters is long. When the equilibrium state is destroyed, it will spend long time to retrieval.

Some further experiments have been conducted to study the variance characteristic of surface tension of sodium hydroxide. **Figures 7** and **8** both reflect the influence of agitation on surface tension. In **Figure 7**, the surface tension of 1 mol/L aqueous sodium hydroxide is measured under the condition where aqueous sodium hydroxide has been in equilibrium state, and the measured point at 35°C is employed to study the influence of agitation on surface tension because this measured point is locate at middle of temperature range and the thermal motion of molecule is not too strong. As shown in **Figure 7**, the surface tension after agitation at 40°C is obviously larger than the surface tension at other temperatures, because the experiments are conducted as solution temperature increases. The agitation has influence on the measured surface tension at higher temperature so that the surface tension at 50°C is larger than the value before agitation at 40°C. Different from the situation at 20°C, the retrieval rate

hydroxide (sodium hydroxide concentration is 1 mol/L).

of hydroxyl ion clusters at 40–50°C is faster. Then, the influence of agitation on surface tension is studied at 73°C as shown in **Figure 8**. The dynamic variance of surface tension after agitation can be observed during one measurement period. At 73°C, the time for complete retrieval is about 3 min. It can be concluded that the retrieval rate of hydroxyl ion clusters increases with solution temperature.

Because sodium thiosulfate belongs to the salt of strong alkali weak acid, the thiosulfate radical hydrolysis can produce hydroxyl ion. However as shown in **Figure 9**, the influence of agitation on the surface tension of sodium thiosulfate is weak. This is because the hydroxyl ion concentration is low and the influence of hydroxyl ion on surface tension can be ignored.

Because the surface tensions of aqueous sodium hydroxide and sodium thiosulfate are sensitive to the measurement operation, the possible steps that affect surface tension are conducted carefully to make sure that the steady surface tension can be measured accurately. As shown in **Figure 10**, the surface tension/ electrolyte concentration gradients are both negative for aqueous sodium hydroxide and sodium thiosulfate. The surface tension/ electrolyte concentration gradients of aqueous sodium hydroxide are larger than the one of aqueous sodium thiosulfate. Besides aqueous inorganic, the surface tensions of different aerosol suspension liquid have been measured. The surface tensions of barium sulfate and titanium dioxide are shown in **Figure 11**. It is seen that the surface tension of solid suspension liquid is close to the one of water. Although it has been reported that aerosol

can reduce the interface energy of liquid (Vafaei et al., 2009), the measured surface tension is statical variable. The aerosol particle is much larger than the liquid molecule or ion cluster, so the surface tension is not altered obviously by solid aerosol. The only exception is 0.1 g/L titanium dioxide suspension liquid whose surface tension is significantly smaller than the value of water. This is because the particle size of titanium dioxide is smaller, and the hydrophobicity of aerosol particle makes the liquid molecule structure loose. However, as the mass portion of titanium dioxide increases, the aerosol particle agglomerates, and the influence of aerosol on surface tension disappears.

# REFERENCES


# CONCLUSION

The surface tensions of aqueous sodium hydroxide, aqueous sodium thiosulfate, and solid aerosol suspension liquid have been measured mainly by platinum ring method. The surface tension/ electrolyte concentration gradients of aqueous sodium thiosulfate are negative, which proves that surface excess concentration or adsorption of the solute at the interface between liquid and gas is positive. Because sodium thiosulfate is easily adsorbed on platinum material, it is applicable to employ platinum ring for measuring aqueous sodium thiosulfate and other similar solution. Due to the formation of hydroxyl ion clusters, the surface tension of aqueous sodium hydroxide is smaller than the value of water. The retrieval time of hydroxyl ion clusters is much smaller than the one of active agent and decreases with temperature increase. Due to large particle diameter and agglomeration phenomenon, the surface tension of solid aerosol suspension liquid has no obvious difference from the surface tension of water.

# AUTHOR CONTRIBUTIONS

HG provides guidance for this study. ZS and KS subsidize this research work.

# FUNDING

These authors are profoundly grateful to the financial supports of the Fundamental Research Funds for the Central Universities (Grant Nos. HEUCF171502). This work is also supported by Key Supported Discipline of Guizhou Provence (Qian Xuewei He Zi ZDXK[2016]24), 2011 Collaborative Innovation Center of Guizhou Province (Qian Jiao He Xietongchuangxin Zi [2016]02).


Wen, J., Gu, H., Sun, Z., and Zhou, Y. (2017). A theoretical model and experiment validation on filtration characteristics of methyl iodide in bubble column. *Int. J. Heat Mass Transfer* 114, 1263. doi:10.1016/j.ijheatmasstransfer.2017.07.023

Zhang, J., and Fan, L. S. (2003). On the rise velocity of an interactive bubble in liquids. *Chem. Eng. J.* 92, 169. doi:10.1016/S1385-8947(02)00189-4

Zhou, X., Gu, H., and Li, F. (2014). Determination of experimental aerosols for filtered containment venting system and selection of technical parameter. *Nucl. Power Eng.* 124. doi:10.13832/j.jnpe.2014.05.0124

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Wen, Shi, Sun, Sun and Gu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*