# NUCLEAR SAFETY DESIGN AND INNOVATION

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

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-663-1 DOI 10.3389/978-2-88963-663-1

### 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 SAFETY DESIGN AND INNOVATION

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

Citation: Wang, J., Shi, K., Meng, Z., Revankar, S. T., eds. (2020). Nuclear Safety Design and Innovation. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-663-1

# Table of Contents


# Editorial: Nuclear Safety Design and Innovation

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

*<sup>1</sup> Department of Engineering Physics, University of Wisconsin-Madison, Madison, WI, United States, <sup>2</sup> College of Nuclear Science and Technology, Harbin Engineering University, Harbin, China, <sup>3</sup> Department of Nuclear Engineering, Purdue University, West Lafayette, IN, United States*

Keywords: nuclear safety, core coolant, thermal hydraulics, simulation, nuclear power plant

**Editorial on the Research Topic**

### **Nuclear Safety Design and Innovation**

The electrical capacity of nuclear energy has consistently increased from the 1950s, first with development prototype reactor and ultimately installation of commercial nuclear power plants (NPPs) (BP, 2019; Budley, 2019). Since then we have had three major severe accidents despite the advances in technology and regulatory framework on the safety and security of NPPs (Budley, 2019). The Fukushima Daiichi nuclear accident has lead Japan to scale back on operations on existing NPPs. Germany has decided to phase out the NPP's operation. However, there are new constructions of NPPs in other countries and China with the fast-developing nuclear market due to the demand for clean energy. Thus nuclear safety is an important and continuing issue with NPPs. As there are new developments and innovations in nuclear safety of NPPs, many experts from different organizations are invited to share their advances in technology and innovations on nuclear reactor safety (Zhou and Zhang, 2010; Zhou et al., 2011). In this topic, many experts from different organizations are invited to share their most advanced progress and innovations.

The first part of this topic focuses on the heat transfer and bubble formation on the surface of fuel bundles. This research includes: "Diameter effect on the wall temperature behaviors during supercritical water heat transfer deterioration in circular tubes and annular channels" (Cheng et al.); "Preliminary experimental investigation on the filtration performance of submicron insoluble aerosol in a bubble column" (Li, Shi et al.); "Bubble formation characteristic of submerged singlehole orifice in aerosol suspension" (Sun et al.); "Visualized experiment of bubble behaviors in a vertical narrow rectangular channel under natural circulation condition" (Yan et al.).

The second part of this topic focuses on the heat transfer in the primary loop of a nuclear power plant and includes: "Numerical analysis of turbulent flow and heat transfer in internally finned tubes" (Liu et al.); "Startup thermal analysis of a supercritical-pressure light water-cooled reactor CSR1000" (Yuan et al.); "Effect of liquid injection arrangements on injection flow rate of a laboratory-scale venturi scrubber" (Zheng et al.); "Design of the container for the sampling and detection monitoring system of N-13 in pressurized water reactor primary loop water leakage based on the coincidence method" (Zhao et al.).

The third part of this topic talks about safety simulations and includes: "Research on timedependent failure modeling method of integrating discrete dynamic event tree with fault tree" (Xu et al.); "A numerical investigation on gaseous stratification break-up phenomenon of air fountain experiment by fountain experiment by code\_sature" (Zhang C. et al.); "Sequential failure modeling and analyzing for standby redundant system based on FTA method" (Zhang M. et al.); "Analysis of categorical subgroup method for resonance self-shielding treatment" (Li S. et al.); "Study on calculation method of soluble aerosol removal efficiency under high humidity conditions" (Li, Ma et al.).

#### Edited and reviewed by:

*Shoaib Usman, Missouri University of Science and Technology, 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: *25 January 2020* Accepted: *24 February 2020* Published: *10 March 2020*

#### Citation:

*Wang J, Shi K, Meng Z and Revankar ST (2020) Editorial: Nuclear Safety Design and Innovation. Front. Energy Res. 8:36. doi: 10.3389/fenrg.2020.00036*

**4**

After a few months' preparation, over 10 manuscripts were collected, to show the ongoing safety research from three different perspectives. Due to limits of time and effort, a mass of excellent work on other important aspects was not included. However, these topic are expected to provide an example to display nuclear safety research at Frontiers in Nuclear Energy. More high-quality topics

### REFERENCES

BP (2019). BP Statistical Review of World Energy. BP World Energy.


and articles are welcomed at this open platform in the future.

### AUTHOR CONTRIBUTIONS

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

**Conflict of Interest:** 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 © 2020 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.

# Sequential Failure Modeling and Analyzing for Standby Redundant System Based on FTA Method

Min Zhang<sup>1</sup> \*, Zhijian Zhang<sup>1</sup> \* and Gangyang Zheng1,2

<sup>1</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, College of Nuclear Science and Technology, Harbin Engineering University, Harbin, China, <sup>2</sup> NeUtron eXploration Team, Beijing, China

Fault Tree Analysis (FTA) has been a well-established and widely used method to deduct system failure scenarios for large complex systems like Nuclear Power Plants (NPPs). Redundant design is usually adopted in NPPs to improve system reliability, including parallel design and standby design. Sequential failures exist among the modules in a standby redundant system, which have not been detailed considered in FTA in industry, leading to an overestimation of system failure probability. Then if FTA is used to compare the reliability of the two designs, it will be found that parallel design is more reliable than standby, which is just the opposite of the conclusion from Reliability Block Diagram (RBD) analysis. To solve this problem, an improved Fault Tree methodology is proposed in this paper, using Priority-AND (PAND) gate and Condition-AND (CAND) gate to model the sequential failures. And the Boolean laws of logic is extended correspondingly for qualitative analysis, as well as the mathematic formulas for quantitative analysis. A case study is also presented to demonstrate the process and benefits for using the proposed approach.

### Edited by:

Jun Wang, University of Wisconsin-Madison, United States

### Reviewed by:

Carlos Antonio Sartin, Universidade de São Paulo, Brazil Keyou Mao, Purdue University, United States

#### \*Correspondence:

Zhijian Zhang zhangzhijian\_heu@hrbeu.edu.cn Min Zhang zhangmin@hrbeu.edu.cn

#### Specialty section:

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

Received: 06 May 2018 Accepted: 06 June 2018 Published: 26 June 2018

#### Citation:

Zhang M, Zhang Z and Zheng G (2018) Sequential Failure Modeling and Analyzing for Standby Redundant System Based on FTA Method. Front. Energy Res. 6:60. doi: 10.3389/fenrg.2018.00060 Keywords: fault tree analysis, standby redundant system, sequential failure, Priority-AND gate, Condition-AND gate

## INTRODUCTION

Fault Tree Analysis (FTA) is a widely used method in system reliability analysis (Lee et al., 1985; Guimarees ˜ and Ebecken, 1999), and has been applied to nuclear power industry since early 1970s (NUREG, 1975).

It is a deductive methodology, which connects Top Event (system failure) with a set of Intermediate Events and Basic Events (component failures, human errors, etc.) by logic units like AND/OR gates. System failure paths and contributions of components/events to system failure can be located by FTA results (Vesely et al., 2002).

Redundancy is an effective measure to improve system reliability, including parallel redundant and standby redundant.

**Figure 1** shows a typical two-redundant system design, in which the blue parts are unique to standby redundant design while the others for both. In parallel design, A and B are activated simultaneously. While in standby, A will be activated firstly, and then B by switching unit S if A fails.

The Minimal Cut Sets (MCSs) of the two systems are:

Parallel: A T B Standby:A T B, A T S

With "no failure in standby" assumption, the standby system failure probability from FTA is bigger than the other, and the conclusion is opposite to that from Reliability Block Diagram analyzing (Bilintion and Allan, 1992). Detailed analysis about the MCSs of the two redundant systems is carried out, and it is found that the FTA has overestimated the failure probability of standby system, by involving two failure paths which should not be considered. They are: (1) B fails before A, which should not happen because B is in standby before A fails; (2) S fails after A, which may happen but should not lead to system failure as B should have been activated before S fails.

These are so called sequential failure problems in this paper, and are defined as:


Sequential failure problems are very common for process systems in which system/components may be placed in service orderly. And several solutions have been proposed to solve these problems in recent years. Representative works including:


The first two kinds of methodology focus on only one sequential issue. DFT is designed with house event tables to take system configuration changes into consideration, but not to analyze

the sequential issues discussed here. ESD is comprehensive in theory, but too complicated to construct a detailed model for large complex systems.

Hence, this research is to develop an FTA method to analyze both of the two sequential issues. The method suggested should be convenient to be implemented in model construction and computer aided analyzing.

### SEQUENCE-DEPENDENT FAILURE MODELING AND ANALYZING

Priority-AND (PAND) gate is adopted to model the SDFs, which has been defined in Fault Tree Handbook (Vesely et al., 2002) but without qualitative and quantitative analysis laws.

The graphic symbol of PAND gate is shown in **Figure 2**, and the mathematic expression is written asQ = A∀B, which means that event Q will occur if and only if:

1) B has occurred, following the occurrence of A; or

2) Both A and B occur simultaneously.

Most of the basic Boolean logic laws (Verma et al., 2010) are still available for Minimal Cut Sets Analysis for FTs involving PAND gates. For the convenience of reading, they are listed as follows:

Distributive Law:

(A S B)∀(C S D) = (A∀C) S (A∀D) S (B∀C) S (B∀D)

Idempotent Laws: A∀A = A Absorption Law: A S A∀B = A S (A∀B) = A

But it should be noticed that the Exchange Law is no longer available for PAND gate, that isA∀B 6= B∀A.

And for multiple SDF, the extended Boolean Laws are:

(A∀B)∀C = A∀B∀C;

A∀(B∀C) = (A T B)∀C, which means:

1) All of A, B and C should occur;


To calculate the probability of MCSs with SDFs, use

$$P\{A\forall B\forall C,t\} = \int\_0^t f\_C(t\_3) \int\_0^{t\_3} f\_B(t\_2) \int\_0^{t\_2} f\_A(t\_1)dt\_1dt\_2dt\_3$$

If B is an event with constant probability of QB, then

$$P\{A\forall B\forall C,t\} = Q\_B \int\_0^t f\_C(t\_3) \int\_0^{t\_3} f\_A(t\_1)dt\_1dt\_3$$

Where,f(t) is the probability density function of an event—the probability that an event may occur at time t.

Any A, B, or C may be a combination of more than one events. IfA = A<sup>1</sup> T A2, then

$$f\_A(t) = \frac{dF(A\_1 \bigcup A\_2)}{dt} = f\_{A\_1}(t)F\_{A\_2}(t) + F\_{A\_1}(t)f\_{A\_2}(t)$$

Where, F(t)is the distribution function of an event—the probability that an event may occur in time duration of t.

Then for the system in **Figure 1**, the MCS A T S becomesS∀A, and its probability is

$$P\{S\forall A,t\} = \int\_0^t f\_A(t\_1) F\_S(t\_1) dt\_1$$

### CONDITION-DEPENDENT FAILURE MODELING AND ANALYZING

Condition-AND (CAND) gate is constructed to model CDFs. The graphic symbol is shown in **Figure 3**, and the mathematic expression is written asQ = hA |Bi , which means that:

1) Event Q will occur if and only if both A and B occur; and 2) B never occurs before A.

The basic Boolean logic laws for CAND gate are as follows:

Distributive Law:

$$\left< (A \bigcup B) \left| (C \bigcup D) \right> = \left< A \left| C \right> \bigcup \left< A \left| D \right> \bigcup \left< B \left| C \right> \bigcup \left< B \left| D \right> \right> \right| $$

Idempotent Laws: hA |Ai = A Absorption Law: A S hA |Bi = A S (hA |Bi) = A

Exchange Law is not available for CAND gate either. And for multiple CDF, use the extended Boolean Law:

$$
\langle\langle A \mid B \rangle \mid \mathcal{C} \rangle = \langle A \mid \langle B \mid \mathcal{C} \rangle\rangle = \langle A \mid B \mid \mathcal{C} \rangle
$$

To calculate the probability of MCSs with CDFs, use

$$P\{\left\right>,t\right>\}=\int\_{0}^{t}f\_{A}(t\_{1})\int\_{t\_{1}}^{t}f\_{B}(t\_{2}-t\_{1})\int\_{t\_{2}}^{t}f\_{C}(t\_{3}-t\_{2})dt\_{3}dt\_{2}dt\_{1}$$

or

$$P\{\left\langle A\left| B\left| C\right\rangle, t \right\rangle = Q\_B \int\_0^t f\_A(t\_1) \int\_{t\_1}^t f\_C(t\_3 - t\_1) dt\_3 dt\_1\right\}$$

Then for the system in **Figure 1**, the MCS A T B becomes hA |Bi, and the probability is

$$P\{\left\langle A\left| B \right\rangle, t \right\} = \int\_0^t f\_A(t\_1) \int\_{t\_1}^t f\_B(t\_2 - t\_1) dt\_2 dt\_1$$

### HYBRID LOGIC ANALYSIS

The two sequential logic gates discussed separately above usually exist simultaneously in one MCS for high-redundancy system, and the logic is hybrid.

**Figure 4** shows a simplified typical composition of a standby triple-redundant system with function modules of A, B, and C. The switching unit is consisted by three sensors (S1–S3) and a processing modular (P). Sensors S1, S2, and S3 are used to monitoring the parameters out from A, B, and C correspondingly. P processes the parameters from sensors and then generates a signal to activate the standby unit orderly if necessary. There are three operation modes for this system:


These modes are completely symmetrical. We assume that the system is operated in Mode1.

The possible hybrid logics and laws to do qualitative and quantitative analysis are:

1) (S1∀A) (S2∀B) , and

$$\begin{aligned} &P\left\{ \left< (\mathcal{S}\_1 \forall A) \left| (\mathcal{S}\_2 \forall B) \right> \right\}, t \right\} \\ &= \int\_0^t f\_A(t\_1) F\_{\mathcal{S}\_1}(t\_1) \int\_{t\_1}^t f\_B(t\_2 - t\_1) F\_{\mathcal{S}\_2}(t\_2) dt\_2 dt\_1 \end{aligned}$$

2) (S1∀A) (S2∀B) = -(S<sup>1</sup> T S2)∀A |B , if B is a failure on demand with a constant probability, and

$$\begin{aligned} &P\left\{ \left< \left( S\_1 \forall A \right) \left| \left( S\_2 \forall B \right) \right> , \, t \right\} \right. \\ &= Q\_B \int\_0^t f\_A(t\_1) F\_{S\_1}(t\_1) F\_{S\_2}(t\_1) dt\_1. \end{aligned}$$

3) A (P∀B) , with P as the unit shared by A and B. Then

$$\begin{aligned} &P\left\{ \left< A \left| \left( P \forall B \right) \right> , t \right\} \right. \\ &= \int\_0^t f\_A(t\_1) \int\_{t\_2}^t f\_B(t\_2 - t\_1) \int\_{t\_1}^{t\_2} f\_P(t\_3) dt\_3 dt\_2 dt\_1. \end{aligned}$$

4) For any A and B, A (B∀A) = B∀A.

### CASE STUDY

Taking the system shown in **Figure 4** as an example to demonstrate the FTA process using the proposed methodology. And define events as follows:


5. SW-component name: The component is not activated because of the failure of relative switching unit.

The FT is shown in **Figure 5**.

S3 does not appear in the FT. It is because that A and B have failed when C is activated and there would be no standby remained once C is also failed. Then S3 is only used to monitor system parameters and find whether system has failed or not, with no contribution to system function failure.

Apply the extended Boolean logic rules, then

SYS = hA |B∗ |C∗i = hA |B ∗ |Ci + hA |B∗ |C − SW i = hA |B ∗ |Ci + A |B∗ (SW <sup>−</sup> <sup>C</sup>∀B∗) = hA |B ∗ |Ci + A (SW <sup>−</sup> <sup>C</sup>∀B∗) B∗ = B + B − SW = B + (SW − B∀A) = B + - (S1 + P)∀A = B + (S1∀A) + (P ∀A) hA |B ∗ |Ci = A [<sup>B</sup> <sup>+</sup> (S1∀A) <sup>+</sup> (P <sup>∀</sup>A) ]|<sup>C</sup> = hA |B |Ci + A (S1∀A) <sup>|</sup><sup>C</sup> + A (P <sup>∀</sup>A) <sup>|</sup><sup>C</sup> = hA |B |Ci + (S1∀A)|C + (P∀A)|C 

$$\begin{aligned} & \left< A \left| \{SW - C\forall B \* \} \right> \\ &= \left< A \left| \{SW - C\forall \left[ B + (S1\forall A) + (P\forall A) \right] \right> \right> \\ &= \left< A \left| \{(S2 + P)\forall \left[ B + (S1\forall A) + (P\forall A) \right] \} \right> \\ &= \left< A \left| (S2\forall B) \right> + \left< A \left| \{S2\forall (S1\forall A) \} \right> + \left< A \left| \{S2\forall (P\forall A) \right> \right> \right> \right> \right> \end{aligned} $$

$$\begin{aligned} &+ \left< A \, \left[ \left( P \forall B \right) \right) + \left< A \, \left[ \left[ P \forall \left( S1 \forall A \right) \right] \right> + \left< A \, \left[ \left[ P \forall \left( P \forall A \right) \right] \right> \right] \\ &= \left< A \, \left[ \left( S2 \forall B \right) \right) + \left< S2 \bigcap S1 \right> \forall A + \left< S2 \bigcap P \right> \forall A \\ &+ \left< A \, \left[ \left( P \forall B \right) \right) + \left< S1 \bigcap P \right> \forall A + P \forall A \\ &= \left< A \, \left[ \left( S2 \forall B \right) \right) + \left< S2 \bigcap S1 \right> \forall A + \left< A \, \left[ \left( P \forall B \right) \right) \right> + P \forall A \end{aligned} \right. \end{aligned}$$

$$\begin{aligned} \text{SYS} &= \langle A \, \big|B \, \big|C \rangle + \left\langle \text{(S1\forall A)} \, \big|C \right\rangle + \left\langle A \, \big|\left(\text{S2\forall B}\right) \right\rangle \\ &+ \left\langle \text{S2} \bigcap \text{S1} \right\rangle \forall A + \left\langle A \, \big|\left(P\forall B\right) \right\rangle + P \forall A \end{aligned}$$

Ultimately, there are six MCSs for system failure:


The MCSs of a classical FT are:


Using the same failure rate of 1E-04 h−<sup>1</sup> for all the components in system, then the system failure probability in 24 h evaluated by the proposed FT method is 2.8892e-06, while the value of classical FT is 5.8013e-06, which is almost twice as much as the former.

From the case study, it is very clear that:


And also it should be reasonable to suspect the accuracy of importance analysis result of the classic FT.

The proposed FT method provides a new perspective to avoid the issues above. It has been applied to an NPP system, under a

### REFERENCES

Azaron, A., Katagiri, H., Kato, K., and Sakawa, M. (2006). Reliability evaluation of multi-component cold-standby redundant systems. Appl. Math. Comput. 173, 137–149. doi: 10.1016/j.amc.2005.02.051

Project on NPP risk monitor. The qualitative analysis to generate MCSs is implemented manually, and MATLAB program is used to do the quantitative calculation. And a software tool to do a complete analysis is under development.

### CONCLUSION

It has been recognized for many years that the order in which the components may fail will affect system behavior in a standby redundant system. But, this issue has not been comprehensively considered yet in industry. The study of this paper aims to find a way to take account of this issue based on FTA methodology.

Firstly, the issue is divided into two categories:


Accordingly, two gates, Condition-AND gate and Priority-AND gate, are adopted in this study. The former gate is new constructed, while the later one has been defined in Fault Tree Handbook but not applied because of the lack of qualitative and quantitative analysis rules in the Handbook.

Then, the extended Boolean logic laws are proposed to implement the qualitative analysis of a FT, generating the MCSs with sequential characteristics. The mathematical formulas to calculate MCS probability (quantitative analysis) are also developed, as well as the variations in application.

Finally, a case study is presented to demonstrate the modeling and analyzing process using the improved FTA method. Compared with the results from classical FTA, it is found that the proposed method can lead to more realistic failure scenarios and reduce the conversation of classical FTA.

The future work will focus on how to improve the efficiency of the method for applying to large complex system (e.g., NPP system) when combined with other problems like common cause failure. Also, it is necessary to develop a software to do analysis automatically.

### AUTHOR CONTRIBUTIONS

ZZ found the engineering problem and proposed to do this research. MZ was in charge of and implemented the research to find the reason and solution of the problem. GZ helped to verify the method in a industrial system.

### ACKNOWLEDGMENTS

This study was supported by the National Science and Technology Major Project Research on Living PSA and Online Risk Monitoring and Management of NPPs [2014ZX06004-003].

Bilintion, R., and Allan, R. (1992). Reliability Evaluation of Engineering System. New York, NY: Springer.

Cepin, M., and Mavko, B. (2002). A dynamic fault tree. Reliab. Eng. Syst. Safety 75, 83–91. doi: 10.1016/S0951-8320(01)00 121-1


**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 and Zheng. 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.

# Visualized Experiment of Bubble Behaviors in a Vertical Narrow Rectangular Channel Under Natural Circulation Condition

Meiyue Yan<sup>1</sup> \*, Tingting Ren<sup>1</sup> , Kailun Chen<sup>1</sup> , Changqi Yan<sup>1</sup> \*, Yongyong Yang<sup>2</sup> , Chunping Tian<sup>1</sup> , Kuan Yang<sup>1</sup> and Dechun Ai <sup>3</sup>

*<sup>1</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China, <sup>2</sup> China Nuclear Power Operation Technology Corporation, LTD, Wuhan, China, <sup>3</sup> School of Mining and Civil Engineering, Liupanshui Normal University, Liupanshui, China*

#### Edited by:

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

#### Reviewed by:

*Yacine Addad, Khalifa University, United Arab Emirates Claudio Tenreiro, University of Talca, Chile*

#### \*Correspondence:

*Meiyue Yan yanmeiyue@hrbeu.edu.cn Changqi Yan changqi\_yan@163.com*

#### Specialty section:

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

Received: *09 August 2018* Accepted: *21 September 2018* Published: *23 October 2018*

#### Citation:

*Yan M, Ren T, Chen K, Yan C, Yang Y, Tian C, Yang K and Ai D (2018) Visualized Experiment of Bubble Behaviors in a Vertical Narrow Rectangular Channel Under Natural Circulation Condition. Front. Energy Res. 6:105. doi: 10.3389/fenrg.2018.00105* The characteristics of bubble behavior have been in particular interest for decades due to its significant contribution to understanding the mechanism of heat transfer. In the present work, visualized experiment is conducted to study the bubble characteristics in subcooled flow boiling of a narrow rectangular channel under natural circulation. The experiments were performed at pressures of 0.2 MPa, with inlet subcooling ranging from 20 to 60 K and heat flux ranging from 100 to 300 kW/m<sup>2</sup> . A high-speed digital camera is used to capture the pictures of bubble behaviors. A sequence of image processing algorithms is used deal with the original bubble images to get relevant bubble parameters. We observe the whole process of a single sliding bubble lifetime and found most of bubbles slide along the heating surface after detaching the nucleation sites. Five typical sliding bubble growth paths are observed in the present experimental conditions. According to the analysis of the experimental data, it can be found that the liquid subcooling and wall superheat are the main factors that affect the bubble size during sliding in narrow rectangular channel under natural circulation condition. Due to the difference of driving force, the sliding velocity of bubble in forced circulation is always larger than that in natural circulation. At the same time, the bubble velocity changes significantly at different heat flux and shooting location.

Keywords: narrow rectangular channel, natural circulation, subcooled boiling, sliding bubble, bubble behavior

## INTRODUCTION

Flow boiling under natural circulation has been found in many industry applications including nuclear reactor, External Reactor Vessel Cooling (ERVC) (Ha et al., 2004), Passive containment cooling System (PCCS) (Park et al., 2008). In the studies concerning about narrow rectangular channel, more attention has been attracted in sliding bubbles for their great improvements on the capacity of heat transfer. Therefore, study on sliding bubble behaviors is helpful to understand the mechanism of heat transfer in a narrow rectangular channel under natural circulation.

Many researchers have explored sliding bubble characteristics. In the experiment of subcooled flow boiling of water by Prodanovic et al. (Prodanovic et al., 2002), bubble slid for the length of several bubble diameters and migrated toward the bulk fluid. And it was found that bubble sliding is the only way for bubble departure from the nucleation site in microchannel in the study carried out by Yin et al. (Yin et al., 2015). Okawa et al. (Okawa et al., 2005a,b) investigated the bubble rise after the departure from a nucleation site. They observed two types of sliding bubble path in vertical upflow boiling: (1) bubbles slide along the wall for long distance without detachment; (2) bubbles detach from the wall after sliding for a certain distance and then migrate toward the bulk liquid. The two types of sliding bubble were also observed by Li et al. (Li et al., 2013) in subcooled upward forced convective flow boiling in narrow channel. There are also many studies on sliding bubble characteristics parameters. From a visualized experiment in a narrow rectangular channel conducted by Xu et al. (Xu et al., 2010, 2013), we found that the maximum bubble dimension in the direction normal to the heating surface is >1 mm in the flow direction, while most bubble diameters are in the range of 0.09–0.4 mm, and the upstream and downstream contact angles were almost equal during sliding process. Chen et al. (Chen et al., 2011) performed a study of vapor bubble behavior under different system pressures in a narrow channel. The sliding bubble size and growth rate decrease with the increase of system pressure. The sliding bubble phenomenon also was found in gas-liquid flow. Zaruba et al. (Zaruba et al., 2007) experimentally and numerically studied the motion of single bubble in a vertical gas-liquid flow, and they found that small bubbles always slide along the wall. However, it is not enough to study the behavior of a single bubble. To understand the bubble population characteristics, Puli et al. (Puli and Kumar, 2012) analyzed the bubble size density distribution under different pressures, and the number of smaller bubbles increase with the increasing pressure.

The mechanism of sliding bubble to enhance surface heat transfer attracts many attentions. Thorncroft et al. (Thorncroft et al., 1998; Thorncroft and Klausner, 1999) conducted experiments in subcooled flow boiling using refrigerant FC-87 in both upward and downward flow. In upward flow, bubble slid along the wall and was not detached from the wall. While in downward flow, bubbles were detached directly or just slid for several millimeters from the nucleation sites, it was found that the heat transfer coefficient in upflow is significantly higher

than that in downflow. Thorncroft believed that the reason for this phenomenon is the existence of sliding bubble enhancing the heat transfer. Ozer et al. (Ozer et al., 2011) studied the subcooled nucleate boiling in horizontal narrow channel using R-11 and Novec 649. He found that the position where the bubble first nucleated was caused primarily by sliding bubbles. B. Donnelly et al. (Donnelly et al., 2009), Y. Yan et al. (Yan and Kenning, 1997), and D. K. Hollingsworth (Hollingsworth et al., 2009) also observed bubbles sliding phenomenon on the inclined heated wall. They concluded that the wave behind the sliding bubbles was responsible for the heat transfer enhancement.

Natural circulation is driven by the buoyancy force due to the density gradient in the loop, does not need any moving part. The mass flow rate depends on the operating parameters, and the mass flow rate is lower comparing with forced circulation under same conditions. The bubble behaviors are various and complicated under natural circulation. At the same time, natural circulation itself is related to bubble behavior. Most of the researches were carried out under forced convective flow, while the sliding bubbles characteristics under natural circulation was less explored. The mechanism of sliding bubble in natural circulation has not been thoroughly understood. Only Zhou et al. (Zhou et al., 2013) studied the single bubble characteristics under natural circulation. Compared with the experimental results under the forced circulation with same thermal parameters, the growth time was longer, the departure velocity was higher and the waiting time was shorter. The objective of the present work is to investigate sliding bubble characteristics in upward subcooled flow boiling in a vertical narrow rectangular channel under natural circulation. In order to describe the features.

### EXPERIMENTAL METHODS

In order to study bubble behaviors in a narrow rectangular channel under natural circulation condition, an experimental loop shown in **Figure 1** is designed and constructed. The experimental loop mainly consists of four major parts, including a pressurizer, a preheater, a test section and a condenser. The pressurizer is connected with nitrogen cylinder to keep the system pressure at 0.2 MPa. The 40 kW preheater is used to control the water temperature at the inlet of the test section. Two K type shield thermocouples with an accuracy of ±0.5 K are used to measure the inlet and outlet temperature of the preheater. After heated by the preheater, the water is delivered to the test section. After exiting the test section, the work liquid is sent to a condenser. The condenser is the only heat sink of the loop where the fluid temperature is reduced before returning to the preheater. To measure the temperature along the test section, 17 N-type shielded thermocouples are arranged along the flow channel, and the temperature distribution points are illustrated in **Figure 2**.

The setup of visualization system is shown in **Figure 3**. A high speed camera is used to record the sliding bubble process. The

camera is mounted on a 2-D traverse system, which can be moved vertically and horizontally. The camera is located at different height of channel to record the bubble behaviors under each operating condition. The height is the vertical distance between the center of camera lens and location of fluid inlet. Two 150 W fiber lights are placed before the test section in order to provide sufficient light for shooting.

**Figure 4** shows cross section of the test section. The total heating length is 550 mm. The test section is made up of stainless steel and glass. To enable visualization of the bubble behaviors, the channel is designed to be single-side heated. The stainless heater plate is heated by a DC power with capacity of 2,000 A/50 V.

Deionized water is used as working fluid. In order to avoid the effect of dissolved non-condensable gases in the water on experimental phenomena, the water kept boiling by preheater and test section for degassing before the experiment. The non-condensation is released from the valve on the top of experimental loop. After degassing, the secondary cooling was activated to create density difference between the riser and downcomer of primary loop. Then centrifugal pump is shut off in order to switch the working mode from forced circulation to natural circulation. Main parameter ranges are summarized in **Table 1**.

In the experiment, the filming frame rate and resolution are set to 5,000 fps and 1,240 × 1,240 pixels. The corresponding physical scale of the view window is about 20 × 20 mm<sup>2</sup> . The uncertainty for locating the position is within ±2 pixels, thus the maximum error of bubble location is limited to ±0.04 mm.

The uncertainties of main parameters are summarized in **Table 2**.


TABLE 2 | Estimated errors of main measurement parameters.


### RESULTS AND DISCUSSION

### Observation of Single Sliding Bubble Behaviors

The growth behavior of single bubble is a complex phenomenon, which has a coupling effect to the wall temperature field, fluid temperature field, fluid velocity, etc. **Figure 5** shows the whole process of a single sliding bubble lifetime. Because of the restriction of narrow channel, there will more sliding bubbles generated on the heating surface. The mass existence of sliding bubble will increase heat transfer coefficient strongly. A whole sliding bubble cycle is mainly divided into two stages: growth stage and waiting stage. The growth stage represents the time during a bubble still attached to its nucleation site. In general, a bubble generates at a nucleation site then grows in the radial direction by adsorbing heat from superheated liquid layer between the bubble and heating wall, and then departs from the nucleation site after its diameter reaches a critical value. The waiting stage describes the time after bubble departs from its site until the next bubble forms at the same nucleation site. After normalizing the bubble diameter, it is found that the in natural circulation the bubble growth time accounts for a larger proportion in bubble cycle than in the forced circulation (Zhou et al., 2013).

The bubble size is determined by the balance of the evaporation rate at the bubble interface in the superheated liquid layer and the condensation at the interface surrounded by the subcooled bulk liquid. When the evaporation rate not equal to the condensation rate, the bubble diameter will change. Consecutive images of bubbles obtained in the experimental condition of 1Tin,sub = 41 K, q = 141 kW/m<sup>2</sup> , five typical bubble growth are shown in **Figure 6**. The image sequences of bubbles 1–5 are shown on **Figures 7**–**11**, respectively. These bubbles are obtained from the same shooting window, but due to the differences of these bubbles in each picture, in order to highlight the bubble shapes and facilitate the acquisition of subsequent bubble parameters, the operations performed in the image processing software are different, thus the backgrounds seem to be different. In particular, the bubble in **Figure 10** disturb the surrounding liquid, which cause the background to be different greatly from other figures. **Figures 7**, **8** show the growth characteristics of two typical types of sliding bubbles in narrow rectangular channel. The diameter change rate of bubble 1 is rapidly, while the bubble sliding distance and survival time are short. By observing the bubble videos, most sliding bubbles were similar to bubble 3 shown in **Figure 9**. The diameter of these bubble keep nearly constant during its sliding stage, and the sliding velocity is constant during the sliding process. The bubble diameter shown in **Figure 10** gradually decrease, such bubbles will eventually disappear which means the condensation in advantage all the time. There is also a kind of bubble shown in **Figure 11**. After the detachment, the bubble diameter will decrease slowly due to the subcooled bulk fluid. However, this kind of bubbles does not collapse in the subcooled liquid but reattaches to the heating wall, then the bubble diameters begin to increase again. The reason for this type bubble existing is that the bubble departing from the heating wall for a period of time then reattaching to the wall due

to the limitation of the gap size of narrow rectangular channel. This kind of bubble is called bouncing bubble, and this kind of bubbles were clearly observed as a typical bubble behaviors in the present experiments.

### Observation of Boiling Phenomenon Under Natural Circulation

Massive amounts of bubble image data are obtained in the experiment. Therefore, a batch-processing method is developed to deal with the bubble images. The pre-processing of bubble images is conducted by Image Pro Plus 6.0, by which a sequence of digital image processing algorithms is applied to deal with the original video in order to extract bubble location and shape information. Based on the image processing results, the bubble diameter of each frame existed in the filming window can be obtained. For the average diameter of a large number of bubbles, the formulas can be calculated as followed

$$D\_{\text{ave}} = \frac{1}{\sum\_{i=1}^{m} n\_i} \left( \sum\_{i=1}^{m} \sum\_{j=1}^{n\_i} D\_{e,j} \right) \tag{1}$$

in the formula, the Dave is the average diameter of all the bubble existing in the filming window. m represents the number of frames, and n<sup>i</sup> represents the number of bubbles at each frame.

The mean sliding diameters under different heights and heat flux are depicted in **Figure 12**. Each point value is the mean value of all captured bubble diameters in the filming. It can be seen from the figure that when the heat flux is less than a certain value, the mean bubble diameter increases with increasing heat flux. This is because an increase in the heat flux will cause an increase in evaporation at the bottom of the bubble, resulting an increase in the bubble diameter. When the heat flux exceed a certain value, the mean diameters change randomly. It is mainly attributed to two reasons: the bubble size prefer to oscillate when the heat flux is high. At the same time, the number of nucleation sites on the heating surface is increased, which are favorable to the generation of both big bubbles and small bubble. This results in a random change in the mean bubble diameter as heat flux increases.

The height H in the **Figure 12** represents the distance between the shooting position of high-speed camera lens and the lower side of the rectangular channel. For a certain operating condition, the mass flow rate and the heat flux are the same at different heights in a channel, thus the wall superheat and liquid subcooling are the main factors that affect sliding diameter. It can be seen that the wall superheat increases and liquid subcooling decreases with the increasing height. The liquid temperature is calculated from the temperatures, which is measured at the inlet and outlet of the test section assuming a linear profile between the two temperature measurement points. Zhang et al. (Zhang et al., 2015) had reported that the bulk flow temperature of subcooled boiling changes almost linearly along the flow direction with the outlet void fraction about 50%. From the bubble filming videos, we can find the void fraction is much <50% in the experiment, thus the assumption is reasonable. The liquid layer evaporation rate under the bubble increases with the increase of wall superheat, and the bulk flow condensation rate decreases with the decease of liquid subcooling. Both of the two factors above are benefits to the sliding diameter increase in the channel.

### Comparison Between Natural Circulation and Forced Circulation

It can be seen from the video that the bubble movement in horizontal direction can be negligible compared with that in vertical direction. So the bubble velocity is calculated according the successive images in the vertical direction. The velocity of the centroid of a bubble represents its velocity can be given by

$$V\_{\varepsilon} = \frac{\chi\_{t\_2} - \chi\_{t\_1}}{t\_2 - t\_1}C\tag{2}$$

where y is the bubble coordinate in vertical direction; t<sup>2</sup> – t<sup>1</sup> denotes the time internal of two images; C is the calibration of bubble image.

$$V\_{\text{ave}} = \frac{1}{\sum\_{i=1}^{m} n\_i} \left( \sum\_{i=1}^{m} \sum\_{j=1}^{n\_i} V\_{e,j} \right) \tag{3}$$

To figure out the differences of sliding bubble under natural circulation and forced circulation, the comparison of mean sliding bubble velocity is shown in **Figure 13**. It can be seen that there are obvious differences of sliding bubble diameter between natural and forced circulation while the operating parameters are the same in the test section. We can found the mean velocity of sliding bubble in forced circulation is larger than in natural circulation. This is mainly because the circulation of the working medium in the forced circulation depends on the drive of the main pump in the experimental loop; while the driving force in natural circulation is the density difference between riser and down-comer in experimental loop. Under the same operating condition, the drag force of the liquid flow on the sliding bubble in forced circulation is larger than that in natural circulation, thus the bubble velocity in the forced circulation is generally larger than the natural circulation. The Lx/Ltot in **Figure 13** represents the ratio of the height of the camera to the length of the entire channel. Since the bubble is subjected to buoyancy force and drag force, it is in acceleration state all the time, so the bubble velocity will increase as the sliding distance increases. For the same reason, the location of ONB will moves down when the heat flux increases. In other words, the bubble sliding distance becomes longer when the shooting position is unchanged. The

other reason is the heat flux (<200 kW/m<sup>2</sup> ) will make the bubble diameters larger, and the increased bubble diameter can result in a larger drag force acting on the bubble, which can cause an acceleration of bubble velocity.

### REFERENCES


### CONCLUSION

The sliding bubble behaviors in a narrow vertical channel under natural circulation subcooled flow boiling condition was studied. In the experiment, water is used as working medium. The test pressure is 0.2 MPa, the inlet liquid subcooling ranged from 20 to 60 K, and the heat flux varied from 100 to 300 kW/m<sup>2</sup> . A highspeed camera was used to capture the behaviors of sliding bubble. The main conclusions of the present work can be expressed as:


### AUTHOR CONTRIBUTIONS

MY carried out experiments and analyzed the experimental results, as well as the writing of paper. TR and KC assisted in experimental work and dealt with the experimental data. CY and DA gave the guideline of this research and modified the language of manuscript. YY developed tools to analyze the bubble behaviors. CT and KY designed the experimental loop and guided the experimental work.

### ACKNOWLEDGMENTS

The authors greatly appreciate support of the Natural Science Foundation of China (Grant No. 11675045). The author thanks to the 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). Additionally, the authors are thankful for the support of the Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, China.

annular gap between a reactor vessel and insulation system. Int. Comm. Heat Mass Transfer 31, 43–52. doi: 10.1016/S0735-1933(03) 00200-8


**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 Yan, Ren, Chen, Yan, Yang, Tian, Yang and Ai. 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(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.

## NOMENCLATURE

Nu Nusselt number


Frontiers in Energy Research | www.frontiersin.org 10

# Startup Thermal Analysis of a Supercritical-Pressure Light Water-Cooled Reactor CSR1000

Yuan Yuan<sup>1</sup> \*, Jianqiang Shan<sup>2</sup> , Li Wang<sup>1</sup> , Dongqing Wang<sup>1</sup> and Xiaoying Zhang<sup>1</sup>

*<sup>1</sup> Sino-French Institute of Nuclear Engineering and Technology, Sun Yat-sen University, Zhuhai, China, <sup>2</sup> School of Nuclear Science and Technology, Xi'an Jiaotong University, Xi'an, China*

Supercritical-pressure light water-cooled reactors (SCWR) are the only water cooled reactor types in Generation IV nuclear energy systems. Startup systems, and their associated startup characteristic analyses, are important components of the SCWR design. To analyze the entire startup system, we propose a wall heat transfer model in a paper, based on the results from a supercritical transient analysis code named SCTRAN developed by Xi'an Jiao tong University. In this work, we propose a new heat transfer mode selection process. Additionally, the most appropriate heat transfer coefficient selection method is chosen from existing state-of-the-art methods. Within the model development section of the work, we solve the problem of discontinuous heat transfer coefficients in the logic transformation step. When the pressure is greater than 19 Mpa, a look-up table method is used to obtain the heat transfer coefficients with the best prediction accuracy across the critical region. Then, we describe a control strategy for the startup process that includes a description of the control objects for coolant flow rate, heat-exchange outlet temperature, system pressure, core thermal power, steam drum water-level and the once-through direct cycle loop inlet temperature. Different control schemes are set-up according to different control objectives of the startup phases. Based on CSR1000 reactor, an analytical model, which includes a circulation loop and once-through direct cycle loop is established, and four startup processes, with control systems, are proposed. The calculation results show that the thermal parameters of the circulation loop and the once-through direct cycle meets all expectations. The maximum cladding surface temperature remains below the limit temperature of 650◦C. The feasibility of the startup scheme and the security of the startup process are verified.

Keywords: Supercritical Water Reactor, SCTRAN, heat transfer coefficient, control system, startup

### INTRODUCTION

Supercritical Water Reactor (SCWR) is the only water cooled reactor type in Generation IV nuclear energy systems. SCWR is an innovative system which is aimed for high thermal efficiency and economy. SCWR works at a high pressure, 25 MPa, with a core outlet temperature up to 500◦C. Moreover, Canada's pressure vessel-type reactor outlet temperature even up to 625◦C. So the cladding temperature can reach 650◦C, far beyond the current pressurized water reactor. Supercritical water properties rapid changes at the trans-critical region, and neutron moderating ability would be weakened. The University of Tokyo was the first to study the startup procedure of

#### Edited by:

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

#### Reviewed by:

*Yacine Addad, Khalifa University, United Arab Emirates Claudio Tenreiro, University of Talca, Chile Hui Cheng, City University of Hong Kong, Hong Kong*

\*Correspondence: *Yuan Yuan yuanyuan5@mail.sysu.edu.cn*

#### Specialty section:

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

Received: *08 September 2018* Accepted: *07 November 2018* Published: *26 November 2018*

#### Citation:

*Yuan Y, Shan J, Wang L, Wang D and Zhang X (2018) Startup Thermal Analysis of a Supercritical-Pressure Light Water-Cooled Reactor CSR1000. Front. Energy Res. 6:127. doi: 10.3389/fenrg.2018.00127*

SCWR. Since the startup procedure involves the process of cooling a reactor from a subcritical to supercritical state, analyzing SCWR-based startup thermo-hydraulic characteristics becomes an important consideration (Oka and Koshizuka, 2001). There exist two main startup strategies for SCWRs. One is constant pressure startup, and the other one is sliding pressure startup (Yi et al., 2004). During constant pressure startup procedure, the reactor starts at a supercritical pressure. This procedure is that there is no need to concern about critical heat flux or two-phase instability in two-phase flow. The disadvantage is that large pumping power is needed to make the system operate. During sliding pressure startup procedure, the reactor starts at a subcritical pressure. Using the traditional variable pressure startup system, it must be assured that the maximum cladding surface temperature is below the criterion and flow instability is avoided during startup at subcritical pressures. To avoid the above issues, Ishiwatari et al proposed e recirculation system, instead of a bypass system (Oka, 2013). The recirculation system was separated from the main line of once-through directcycle. Compared with the constant pressure startup procedure, the advantage of the sliding pressure startup system is that the pumping power consumption is decreased, which results in a higher efficiency during low load operations. Therefore, the recirculation system startup procedure is applied in this paper. The past investigation of startup was carried out by several researchers. Yi et al analyzed detailed thermal behavior for a high-temperature supercritical-pressure light water-cooled thermal reactor Super LWR (SCLWR-H) (Yi et al., 2004) and a Super FR with downward-upward two-pass coolant flow (Cai et al., 2009). A sliding pressure startup scheme of recirculation system for SCWRs was proposed to provide a pressurization system which was raised independently from the power (Yamada and Ishiwatari, 2009). With the modified ATHLET-SC code, which realizes the simulation of trans-critical transients using two-phase model the thermal-hydraulic dynamic behavior of the mixed SCWR core during the startup process is simulated (Fu et al., 2011). It can be concluded that the focus of previous work are mainly on the thermal hydraulics behavior of SCWR during startup to validate the limiting criteria, but there are few work on the control system design to achieve startup procedure.

The SCWR coolant loop is a once-through direct cycle that is very sensitive to disturbances. The startup procedure needs to be operated by the control system. Yuki et al. designed a control system of the SCLWR-H under steady state. The turbine inlet pressure is controlled by the turbine control valves. The main steam temperature is controlled by the feedwater flow rate. The core power is controlled by the control rods (Ishiwatari et al., 2003). Wang et al. sing FORTRAN program to make code, when perturbation happens such as pressure setpoint changes, feed-water temperature drop, reactor parameters changes with time were researched at the control rod, steam turbine control valves and there actor coolant pumps as the plant control system (Wang et al., 2012). To analyze the control characteristics of SCWRs, Ishiwatari et al. (2003); (Ishwatari et al., 2010) optimized the SCR control parameters, modified the feed water control system of Super LWRs and improved the delay problem of the temperature control system. Sun and Jiang (2012) studied the control system of Canadian SCWRs and obtained the control system linearization equations via the solution of a complete nonlinear dynamic equations. Dong et al. (2016) proposed a Moving Boundary Method to control the steam temperature. To meet the startup requirements of the entire system, and based on the three typical PID control systems designed by Nakatsuka et al. (1997), control systems are proposed that have controls for coolant flow rate, heat exchanger outlet temperature, system pressure, core thermal power, steam drum water level and the once-through direct cycle loop inlet temperature. By adjusting the control parameters, the startup procedure of the SCWR can be operated according to the starting curve.

The SCWR startup procedure must ensure that the Maximum fuel cladding surface temperature (MCST) does not exceed the limit 650◦C. Accordingly, the heat transfer coefficient value is very important for the calculation of the cladding temperature. A complete set of heat transfer coefficients is needed to meet the wall heat transfer requirements for a smooth transition between the subcritical and supercritical heat transfer coefficients. Ishiwatari et al. (2005) used the "oka-koshizuka correlation" (or "Kitoh correlation") to calculate the heat transfer deterioration under supercritical region of SCWR. oka-koshizuka correlation was developed on the basis of the typical singlephase turbulence model, which is in good agreement with the experimental data of Yamagata (Yamagata et al., 1972; Koshizuka et al., 1995). However, the calculated values of heat transfer coefficient of oka-koshizuka correlation is larger than Watts correlation and Bishop correlation in the deteriorating area (Ishiwatari et al., 2005). SPRAT-DOWN code which is developed by Ishiwatrai adopts oka-koshizuka correlation in rated, transient and accident conditions. Thermophysical properties and transport properties of coolant change greatly during subcritical pressure to supercritical pressure under SCWR startup procedure. It needs a wall heat transfer model with wide-scope parameter which can satisfy the heat transfer coefficient smooth transition requirements from subcritical and supercritical region. By summarizing a large number of supercritical regional databases and literature, Zahlan et al. (2015) proposed a look-up table to calculate trans-critical heat transfer in water-cooled tubes. The SCTRAN code is chosen to analyze the SCWR startup control characteristics. This choice is made because SCTRAN (Wu et al., 2013) can smoothly calculate the physical properties of water in the crosscritical region. The wall heat transfer model of the existing SCTRAN code has a heat transfer coefficient that changes drastically near the critical point and the heat transfer coefficient jump in supercritical regions. The paper puts forward the new wall heat transfer model. The new wall heat transfer model proposes a new logic for the selection of the heat transfer modes and the most appropriate model from the literature is chosen to determine the heat transfer coefficient. Moreover, the discontinuous heat transfer coefficients in the logic transformation are removed. The purpose of this paper is to design a control system for the SCWR's circulation sliding pressure startup process based on an improved SCTRAN code. This work is the foundation of a complete and reliable startup analysis system.

#### TABLE 1 | The heat transfer correlations used in SCTRAN.


### IMPROVEMENT AND VERIFICATION SCTRAN

SCTRAN was developed in Xi′ an Jiao Tong University. It is a one-dimensional safety analysis code for SCWRs (Wu et al., 2013). SCTRAN code has been verified by relap5 and other codes. Typical accidents such as LOCA, LOFA and so on of CANDU-SCWR, CGNPC-SCWR and CSR1000 are calculated (Gou et al., 2012). A homogenous equilibrium mixture model with an optional phasic slip formulation was applied in code SCTRAN.



#### TABLE 3 | SCWR control system.


#### TABLE 4 | Experimental conditions.



The homogenous model made the following assumptions: (1) the velocities of gas and liquid were the same; (2) gas and liquid were in thermodynamic equilibrium. Here the basic equations of homogenous equilibrium mixture model are summarized. Mass conservation equation

$$\frac{\partial}{\partial t}\rho A + \frac{\partial}{\partial z}W = 0\tag{1}$$

where:

$$
\rho = \alpha\_{\mathfrak{g}} \rho\_{\mathfrak{g}} + \alpha\_{l} \rho\_{l} \tag{2}
$$

$$W = W\_{\mathcal{S}} + W\_{l} = \left(\alpha\_{\mathcal{S}} \rho\_{\mathcal{S}} \nu\_{\mathcal{S}} + \alpha\_{l} \rho\_{l} \nu\_{l}\right) \tag{3}$$

Momentum conservation equation

$$\frac{\partial}{\partial t}W + \frac{\partial}{\partial z}W\nu = -A\frac{\partial p}{\partial z} - \frac{2A\rho\nu\,|\nu|}{D\_h}f\_{lp} + \rho A g\_z \tag{4}$$

Energy conservation equation

$$\begin{split} \frac{d}{dt}U &= -\frac{1}{2}\frac{L}{A}\frac{d}{dt}\left(\frac{W^2}{\rho}\right) \\ &- \sum\_{j} \left[ \left(W\_{\mathcal{S}}h\_{\mathcal{S}} + W\_{l}h\_{l}\right) + \frac{1}{2}\left(W\_{\mathcal{S}}\nu\_{\mathcal{S}}\nu\_{\mathcal{S}} + W\_{l}\nu\_{l}\nu\_{l}\right) \right. \\ &\left. + W\_{\mathcal{S}}\left(z - z\_{j}\right) \right] + Q \end{split} \tag{5}$$

In numerical calculation, SCTRAN code adopts the staggered grid in fluid space discretization and adopts the control volume balance method to discrete fluid basic equations.

The SCTRAN module call diagram is shown in **Figure 1**. SCTRAN code mainly includes the following several modules: input module, output module, conservation equation module, pressure calculation module, thermal power calculation module, structure temperature calculation module, heat transfer coefficient and heat flux calculation module, subcritical region and supercritical region fluid physical parameters calculation module. The calculation flow chart of SCTRAN is shown in **Figure 2**. The iterative calculation mainly determines judges whether the change speed of internal energy 1U/U is convergent and whether it reaches the end time of calculation.

### Wall Heat Transfer Model With a Wide-Scope Parameter

The wall heat transfer model mainly includes three parts: heat transfer modes, selection procedures of heat transfer modes (logic) and heat transfer correlations. To analyze the wall heat transfer characteristics of the SCTRAN code, a simulation of a pipe was made. The pipe fluid inlet temperature is 100◦C, the outlet temperature is 500◦C, the fixed flow rate per unit area is 500 kg·m−<sup>2</sup> ·s −1 , and the pipe is uniformly heated along the axial direction; the pressure varies from 1 to 28 MPa. The curves of the heat transfer coefficient are associated with pressure and enthalpy; they are shown in **Figures 3**, **4**. There are some problems with the SCTRAN calculations of the heat transfer coefficient. They include: (1) SCTRAN adopts Jackson heat transfer correlation for supercritical regions, which leads to a

large error. (2) The wall heat transfer model is not comprehensive for all modes. (3) The heat transfer coefficient drastically changes near the critical point and the mode transition region. To solve these problems, the project introduces a new wall heat transfer model that can be calculated from the subcritical to supercritical region. The model selects the appropriate heat transfer correlation and a lookup table (Zahlan et al., 2015) to improve the accuracy of the calculation (**Table 1**). In the sub-critical region, since the heat transfer coefficient of the convection term of the Chen correlation is calculated by using the DB correlation for hmac, after the heat transfer mode is converted, the DB relational expression can be smoothly transferred to the Chen correlation. Then the values of qCHF and TCHF are

the results of the iterative calculation of the table and Chen correlation, so the nucleate boiling correlation can be smoothly transitioned to transition boiling. While the rest of the relations are transformed by linear interpolation. In the supercritical region, only the look-up table for trans-critical region is used.

### Control System Design

Per the SCWR startup sequence (**Table 2**), the startup procedure is divided into four phases: raising of feed water temperature, pressurization, switching from recirculation line to once-through line and power-raising. The control of the SCWR circulation startup system should be able to meet the control requirements

needed for the pressure, thermal power, temperature, steam drum water level and coolant flow rate. Therefore, we designed six different control systems for the SCWR startup procedure (**Table 3**). The basic control methods and equations can be introduced by reference (Nakatsuka et al., 1997; Ishiwatari et al., 2003). Based on the HPLWR thermal cycle, SCWR circulation loop, CSR1000 design parameters and the control method proposed in this paper, a complete CSR1000 control system (**Figure 5**) is designed.

the first process flow rate. (H) Transient performance of the second process flow rate.

## Verification of the Wall Heat Transfer Model

Bailey experiment and Subbotin experiment (Wu et al., 2016) were selected to validate the calculation accuracy of the wall heat transfer model with wide-scope parameter. The inner diameter and wall thickness of Bailey and Subbotin experiment are 12.7762 mm and 2.3368 mm, respectively. The heated length of the test section of Bailey and Subbotin experiment are 3.048 and 6.000 m, respectively. The experimental conditions are shown in **Table 4**. By comparing the wall temperature obtained by the codes with the experiment data, we can verify the modified model.

The CHF calculation correlation used in the new wall heat transfer model is the Groeneveld CHF look-up table developed in 2006 (Groeneveld et al., 2007). From the validation results of Bailey experiment (**Figure 6**) and Subbotin experiment (**Figure 7**), we can see that the modified model better predicted the wall temperature during film boiling under all cases, and it got the relatively good results of the locations where CHF occurred in most cases. In the prediction of the CHF occurrence location, the SCTRAN new calculation results were early, which was caused by the inherent error of the CHF look-up table. The low heat transfer coefficient value obtained from the look-up table leads to a higher wall temperature and makes the CHF occurrence location in advance.

### Verification of the Control System

Based on the various disturbances, this part demonstrates the steady ability of control system. The introduced linear perturbations include main steam temperature, system pressure, core power, coolant flow, and so on (**Table 5**).

When subjected to various disturbances, under the control of the control system, each control system can effectively controls the control target quickly and steadily (**Figure 8**). It shows that the control system designed in this thesis is fast and effective. Response of the main steam temperature control system.

### ANALYSIS OF RECIRCULATION SLIDING PRESSURE STARTUP SYSTEM

Typical startup schemes use a sliding startup system that includes a circulation loop for startup and a once-through direct loop (**Figure 5**). The whole startup procedure needs to be completed in approximately 61.59 h. The variation of thermodynamic parameters is shown in **Figure 9**. During the startup procedure, with the thermal power control system and the flow control system controlling, the core power and the flow rate operates according to the requirements (**Figures 9A,B**). The doppler reactivity, void fraction reactivity and reactivity inserted by the control rod can be seen in **Figure 9B**. At the beginning of the fourth stage, the flow rate of the water rod in the second process is negative, and the coolant is counter-current. At about 41.2 h, the flow rate of the water rod in the second process becomes positive, and the coolant flows normally. At this time, the power and flow rate is about 59%.With the rise of the coolant flow, the buoyancy resistance is not enough to resist the pressure drop increase with the driving force, caused the second process water rod flow rate from the larger negative value changes into larger positive value. It makes the value of void fraction reactivity dramatic changes in this place. The water level can be controlled,

### REFERENCES


up to a height of about 3 m (**Figure 9C**), via the steam drum water level control system. The steam drum water level control system is deactivated at 20 MPa (15.6 h). Under the control of the once-through direct-cycle loop-inlet control system, the coreinlet coolant temperature can be maintained at 280◦C. The outlet coolant temperature gradually increases with time until the rated operating temperature is 500◦C (**Figure 9D**). The MCST does not exceed the limit temperature of 650◦C throughout the startup procedure (**Figure 9E**). The lower plenum will not generate steam during the entire startup procedure and the supercritical steam occurs only in the lower plenum after the supercritical pressure has been achieved. Thus, the coolant remains in a singlephase state and there is no boiling crisis or flow instability problem within the core (**Figure 9F**). Because the buoyancy is greater than the driving force, during the lowering of the water in the second process, the coolant reverses from the subcritical region to the supercritical region (**Figures 9G,H**). The mass flow rate greatly affects the surface temperature of the cladding and the value of the reactivity feedback. The abrupt flow generally causes avoid reactivity great changes (**Figure 9A**).

### CONCLUSION

A new wall heat transfer model is developed for SCTRAN applications to analyze the SCWR startup characteristics. In the model, drastic changes to the heat transfer coefficient (calculated by the SCTRAN) near the critical region is resolved and the heat transfer coefficient from a subcritical to supercritical pressure is forecast precisely and smoothly.

Because the thermo-physical properties and transport properties of coolant change significantly from a subcritical to supercritical pressure, a control system is required to adjust the parameter changes during the startup procedure. Under the control strategy of the startup procedure, the system pressure, temperature, thermal power and flow rate can be regulated according to the startup objectives. The calculation results show that the thermal parameters of the circulation loop and the once-through direct cycle meet the requirement and the MCST remains below the limit temperature of 650◦C.

### AUTHOR CONTRIBUTIONS

YY wrote the manuscript. JS guided this research. LW, DW, and XZ critical revised the article.


**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 Yuan, Shan, Wang, Wang and Zhang. 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.

## NOMENCLATURE

*A*: flow area/m<sup>2</sup> *e(t)*: main steam temperature deviation from set point, % *Dh* : Hydraulic equivalent diameter, m *ftp*: friction coefficient *g*: gravitational acceleration, m·s –2 *h*: enthalpy, J·kg–1 *K1*: lead-lag coefficient *K2*: integral coefficient *KP*: proportional gain *KI* : integral gain *L*: length, m *T1*: lead time, s *T2*: lag time, s *Tmean*: measured temperature by the thermometer, ◦C *Tset*: main steam temperature set point, ◦C *Tstream*: main steam real temperature, ◦C *t*: time, t *U*: internal energy, J·kg–1 *u*(*t*): feed water flow rate signal, kg/s 1*V*: valve opening relative to initial time *v*: fluid velocity, m·s –1 *W*: flow rate, kg·s –1 *z*: distance, m *v*ρ: insertion reactivity speed, cent/s *v*ρmax: the maximum insertion reactivity speed, cent/s <sup>ρ</sup>*mix*: mixture density, kg/m<sup>3</sup>

# Study on Calculation Method of Soluble Aerosol Removal Efficiency Under High Humidity Condition

Yingzhi Li <sup>1</sup> \*, Qianchao Ma<sup>1</sup> , Zhongning Sun<sup>1</sup> , Haifeng Gu<sup>1</sup> , Yanmin Zhou<sup>1</sup> \* and Kaiyi Shi <sup>2</sup>

<sup>1</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China, <sup>2</sup> School of Mining and Civil Engineering, Liupanshui Normal University, Liupanshui, China

Pool scrubbing is a potential method to remove aerosol particles under accident conditions of nuclear power plants. The relative humidity of aerosol laden gas will increase significantly when passing through the water pool, which will most likely induce hygroscopic growth of soluble aerosol. The hygroscopic growth of soluble aerosol can lead to the deviation of the size distribution of aerosol at the outlet of the water pool, resulting in a large evaluation error of the removal efficiency. In order to solve this problem, the size distributions of sodium chloride at the upstream and downstream of bubble column were measured at a gas flow rate of 4 lpm and a liquid height of 80 cm, respectively. Two methods of calculating aerosol actual size-dependency removal efficiency were proposed and compared. One method is directly calculating the removal efficiency by adding a diffusion drying tube installed at the upstream of Scanning Mobility Particle Sizer (SMPS) to reduce the relative humidity of sample gas below the efflorescence point. Another method is modifying the aerosol size distribution and concentration curve by using the hygroscopic growth theory of soluble aerosols. The experiment results obtained by the two methods are in good agreement.

### Edited by:

Muhammad Zubair, University of Sharjah, United Arab Emirates

#### Reviewed by:

Anca Melintescu, Horia Hulubei National Institute for R&D in Physics and Nuclear Engineering (IFIN-HH), Romania Ivo Kljenak, JoŽef Stefan Institute (IJS), Slovenia

#### \*Correspondence:

Yingzhi Li lyz1993@hrbeu.edu.cn Yanmin Zhou zhouyanmin@hrbeu.edu.cn

#### Specialty section:

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

Received: 27 November 2018 Accepted: 18 February 2019 Published: 12 March 2019

#### Citation:

Li Y, Ma Q, Sun Z, Gu H, Zhou Y and Shi K (2019) Study on Calculation Method of Soluble Aerosol Removal Efficiency Under High Humidity Condition. Front. Energy Res. 7:26. doi: 10.3389/fenrg.2019.00026 Keywords: soluble aerosol, hygroscopic growth, removal efficiency, pool scrubbing, FCVS

## INTRODUCTION

Containment is the last shield for protecting the safety of nuclear power plants (NPPs). It is of significance to guarantee its integrity. However, once a long-term loss of core coolant accident (LOCA) occurred, the core molten materials would melt through the pressure vessel. Then, large amounts of coolant would flow into the containment and experience flash evaporation, accompanied by large amounts of non-condensable gas generated by the interaction of a core meltdown and concrete in the pitch. It would cause the pressure inside the containment to rise persistently. Once the containment pressure rises beyond the threshold pressure, rupture of containment would appear and cause leakage of radioactive material into the environment. Filter Containment Venting System (FCVS) is designed to reduce the pressure in the containment through actively venting and reducing the radioactive release through multi-scrubbers. Because of the advantages of a larger specific transfer area and long bubble resident time in liquid pool, bubble column possesses a potential capability for aerosol retention. So far, the bubble column has been used in FCVS, such as the FCVS system composed of wet scrubber and metal fiber designed by Control Components Inc. (CCI) which has been installed in Swiss Nuclear Power Plants (Basu, 2014).

Aerosol as the main component of fission products behaves strong diffusion and migration properties, thus, it is important to prevent aerosol leak into the environment. According to its solubility, aerosol can be divided into soluble aerosol (such as CsI, CsOH, etc) and insoluble aerosol (such as <sup>90</sup>Sr, SnO2, etc) (Jokiniemi, 2007; Tamotsu et al., 2007; Miwa et al., 2015). Soluble aerosol has hygroscopic properties. It means soluble aerosol can absorb water when the relative humidity (RH) of the surrounding environment increases. In the process of increasing relative humidity, significant particle growth occurs only once a critical relative humidity is reached (the exact value of which depends on the temperature and the initial size of the dry particle), called the deliquescence point (DRH). At the same time, a droplet of saturated solution is instantaneously formed. However, the dissolved aerosol will lose water when the relative humidity keeps decreasing. Similarly, it will transform into a crystalline particulate once a critical relative humidity is attained which is called the efflorescence point (ERH). Besides, the hysteresis phenomenon (DRH>ERH) also exists (Hämeri et al., 2001; Gysel et al., 2002). The relative humidity of aerosol laden gas will increase significantly when passing through the liquid pool, which will cause the hygroscopic growth of soluble aerosol. For the aerosol removal efficiency evaluation method based on the mass concentration or the count concentration, this hygroscopic growth phenomenon will lead to a significant deviation in the calculation results of aerosol removal efficiency (Peyres and Herranz, 1996).

Similar scenes of accidents also include the suppression pool in the boiling water reactor and the rupture of the steam generator tube in the pressurized water reactor. All of the above accidents involve the process of aerosol laden gas passing through a liquid pool, which will cause the relative humidity of carried gas to increase (Herranz et al., 1997). The research on the hygroscopic properties of soluble aerosol mainly focuses on the atmospheric environment field. In terms of experimental research, Gysel studied the hygroscopic properties of 100 nm sodium chloride (NaCl), ammonium sulfate ((NH4)2SO4) and sodium nitrate (NaNO3) under different humidity conditions through H-TDMA experimental facility (Gysel et al., 2002). The experiment results showed that NaCl and (NH4)2SO<sup>4</sup> had an obvious deliquescence point and efflorescence point, while NaNO<sup>3</sup> did not show an obvious deliquescence and efflorescence phenomenon. Biskos measured the hygroscopic growth characteristics of 6–60 nm NaCl, and found that the hygroscopic growth factor of NaCl decreased with the decrease of particle size, yet the deliquescence point increased with the decrease of particle size (Biskos et al., 2006). In terms of theoretical research, the hygroscopic growth model of aerosol is mainly based on the Kohler theory (Kohler, 1936), which is based on the spherical hypothesis and takes into account the Kelvin effect of particle size, and proposes a theoretical relationship between particle size and the relative humidity when aerosol particle reaches the hygroscopic equilibrium. Subsequently, Biskos et al. (2006) considered the shape correction factor related to particle size and modified the hygroscopic growth model. In this study, a Scanning Mobility Particle Sizer (SMPS, TSI3082), equipped with a condensation particle counter (CPC) were used to analyze the aerosol size distribution and concentration at the inlet and outlet of a bubble column. It is expected to modify the size distribution and concentration curve of aerosol at the outlet of the bubble column with the help of the hygroscopic growth theory to obtain the actual size-dependency of removal efficiency.

### EXPERIMENTAL SETUP

**Figure 1** shows the schematic diagram of the aerosol test facility. The experiments were conducted in a rectangular bubble column, which is made of polymethyl methacrylate with a cross section of 150 × 150 mm and a variable height of up to 1.5 m. As model soluble particles, sodium chloride was used because the hygroscopicity of NaCl has been deeply investigated in the previous research. The sodium chloride particles were generated by atomization of 3 wt% salt solutions. The salt aerosol generator(TSI Model 8118A)uses compressed air, supplied to five atomizer jets within the liquid reservoir producing an aerosol mist. The atomized solution droplets were dried in a heating tube and further dried by clean and dry air in the mixed chamber. The geometric mean diameter of the number size distribution is 66 nm with a geometric standard deviation of 1.6, which was satisfied with lognormal distribution (**Figure 2**). Additionally, it is possible to guarantee the aerosol concentration constant by setting the mass flow controller FM1 and FM2.

Bubbles are generated in a stagnant liquid by two steel nozzles with an inner diameter of 1.4 mm. Deionized water was selected as the retention liquid to prevent formation of secondary aerosol particles when the bubbles burst at the liquid surface. By setting the valves V3 and V4, it allows to adjust the gas flow injected into the column. The mass flow controller FM3 allows to know the bypass flow and thus to determine the flow injected into the column by balance. This gas flow can also be measured with the mass flow controller FM4 (V7 closed).

The characterization of the aerosol size distribution is analyzed upstream and downstream of the bubble column with a TSI Scanning Mobility Particle Sizer (SMPS), which associates a Condensation Particle Counter (CPC) with a Differential Mobility Analyzer (DMA). To determine the aerosol particles concentrate upstream of the bubble column, the value V5 is opened so that the SMPS can pump a suitable flow. Similarly, to determine of the aerosol particles concentrate downstream of the bubble column, the value V6 is closed and the value V7 is opened so that the SMPS can collect the aerosol laden gas flow. For each particle size, the removal efficiency is determined from these concentration measurements upstream and downstream of the bubble column.

The experiment was conducted at the gas flow rate of 4 lpm and liquid height of 80 cm. The gas flow rate was measured by GFM37 mass flow meter with measurement accuracy of 1%. The temperature and RH of outlet aerosol laden gas were determined by measurement of a temperature and humidity sensor which inserted into the top of the bubble column.

**Figure 3** shows the schematic diagram of the long DMA which is the key component of SMPS. It is a cylindrical configuration. A filtered flow loop recirculates the sheath flow. Sheath air enters

the Sheath Flow+ port on the top of the DMA and passes to an annular chamber at the top of the DMA. The flow then goes through a Dacron mesh to straighten the flow. The air flows downward axially through the classifier region. The polydisperse aerosol was drawn into the inlet at the top of the DMA. To reduce aerosol loss due to Brownian diffusion, the length of the annular aerosol passage was shortened while maintaining a laminar and steady flow. This thin annular flow is introduced into the classifier region and smoothly merged with the laminar sheath air flow. The aerosol flow rate of 0.3 lpm and sheath air flow rate of 3 lpm were preset. Because aerosol of different sizes has different electric charges, the polydisperse aerosol was classified in the electric field which is provided by the highvoltage rod at the center of the DMA. The size-classified aerosol continues to the external CPC to be counted. Besides, the RH of the sheath air was measured by a humidity transducer installed in circulation loop.

### HYGROSCOPIC GROWTH THEORY OF SOLUBLE AEROSOL

Aerosol particles that are soluble in water exhibit hygroscopic properties such that they can absorb moisture from an atmosphere with relative humidity less than 100%. This effect will lead to a growth of the particle size as water vapor condenses into the soluble particle. Maxwell's classical differential equations, which describe particle size and mass changes that result from water vapor diffusion to or form a single spherical droplet, was used to track growth and evaporation of atmospheric condensation nuclei (Hiller, 1991; Pruppacher and Klett, 1997). This approach was also applied to study the effect of hygroscopic growth of soluble aerosol. Previous study shows that the equilibrium time of nanoscale soluble

aerosol is very short, usually less than 1 s under RH < 90% (Broday and Georgopoulos, 2001).

The growth factor (GF) is used to indicate the hygroscopic ability of soluble aerosol. It is defined as the ratio between the particle equilibrium diameter dRH in humid air and the volume equivalent diameter d<sup>0</sup> of the solid particle. Considering that a spherical salt droplet will be formed due to the surface tension of droplet if deliquescence of sodium chloride occurs. Theoretical growth factors can be calculated from the solution concentration csol, the solution density ρsol and the salt density ρ<sup>s</sup> by:

$$\lg\left(RH\right) = \frac{d\_{RH}}{d\_0} = \left(\frac{100\rho\_s}{c\_{sol}\rho\_{sol}}\right)^{1/3} \tag{1}$$

The Kohler theory gives the equilibrium relationship between particle diameter and relative humidity (Kohler, 1936). In the equilibrium state, the RH is equal to the ratio between the vapor pressure over the solution droplet P<sup>d</sup> and the saturation vapor pressure of water P<sup>s</sup> , which is also called the saturation ratio Sr. The inclusion of the effect of the surface tension and aqueous solution lead to

$$RH[\%] = S\_r = \frac{P\_d}{P\_s} = a\_w S\_{Kelvin} = a\_w \exp\left(\frac{4M\_w \sigma\_{sol}}{RT \rho\_w d\_m}\right) \tag{2}$$

Where a<sup>w</sup> denotes the water activity and SKelvin is the Kelvin correction factor. M<sup>w</sup> is the molecular weight of water, and σsol

is the surface tension of the solution. R is the ideal gas constant, and T is the temperature at the particle surface. The empirical values of solution of density, surface tension and water activity as a function of concentration and temperature are used for the theoretical calculations (Pitzer, 1991; Gysel et al., 2002).

It is well-known that dry sodium chloride particles are not spherical but instead have a cubic shape (Hämeri et al., 2001). Therefore, the measured mobility diameter of sodium chloride must be shape corrected to get the volume equivalent diameter by adding a shape correction factor to ensure the aerodynamic drag force unchanged. The drag force on non-spherical particles are described with a modified Stokes' law (Willeke and Baron, 2011)

$$F\_{drag} = \frac{3\pi\,\eta\nu\chi\,a\_0}{C\_c\,(a\_0)}\tag{3}$$

Where

$$\mathcal{C}\_{\rm c} \left( d \right) = 1 + \frac{2\lambda}{d} (1.165 + 0.483 \exp(-\frac{0.997}{2\lambda/d})) \tag{4}$$

$$\chi = \begin{cases} 1.08 & 2\lambda/d\_0 < 0.1 \\ 1.08 \frac{C\_\text{\textdegree C} \{d\_0\}}{C\_\text{\textdegree C} \{1.15d\_0\}} & 0.1 < 2\lambda/d\_0 < 10 \\ 1.24 & 2\lambda/d\_0 > 10 \end{cases} \tag{5}$$

Here η is the gas viscosity, v is the velocity, and λ is the dynamic shape factor, which is concerned with the size of the aerosol particles (Decarlo et al., 2004). C<sup>c</sup> is the Cunningham slip correction factor, a0is the dried sodium chloride size, λ is the mean free path of gas molecules.

Therefore, the modified coefficient of hygroscopic growth factor is:

$$f\_{cub\varepsilon} = \frac{d\_o}{d\_{\nu\varepsilon}} = \frac{1}{\chi} \frac{C\_\varepsilon \left(d\_0\right)}{C\_\varepsilon \left(d\_{\nu\varepsilon}\right)}\tag{6}$$

Growth factor of 65 nm NaCl as a function of RH for T = 20◦C calculated by Equations (2)–(8) was plotted in **Figure 4**. It is obvious that the hygroscopic growth factor of NaCl particles is closely related to the relative humidity of the environment. Furthermore, the growth factor of NaCl particles increases faster with the increase of the relative humidity. Besides, growth factor of NaCl particles as a function of particles size for T = 20◦C and RH = 72.4% calculated by Equations (2)–(8) was plotted in **Figure 5**. It can be seen that the hygroscopic growth factor of NaCl particles increases with the increase of aerosol size.

### RESULTS AND DISCUSSION

In order to clarify whether entrainment droplets caused by the bubble bursting on the surface of the liquid pool will have noticeable impact on the measurement results, the aerosol concentration at the inlet and outlet of the bubble column were measured with an aerosol generator turned on and turned off, respectively. The gas flow rate was 4 lpm and liquid level was 80 cm. The experimental results are shown in **Figure 6**. It can be seen that when the aerosol generator was turned on, the aerosol particles number concentration in the inlet and outlet of the bubble column are much higher, almost four orders of magnitude higher than that without the aerosol supplied. This means that the background aerosol fraction produced by the bubble bursting is negligible. In addition, the stability of aerosol concentration at the inlet and outlet also indicated the reliability of the experimental results.

The aerosol laden gas enters the bubble column through stainless nozzles, the final gas phase will exist as discrete bubbles in the liquid pool. During the rising period of bubbles, the relative humidity of the gas in the bubbles will increase gradually due to the strong heat and mass transfer between gas and liquid phases. According to our experiment's results, the relative humidity of the gas above the liquid surface exceeded 85%, which is much greater than the deliquescence point of NaCl (76%). Therefore, the NaCl particles will absorb water and grow up which will result in the difference of aerosol size distribution between the inlet and outlet of the bubble column. A 0.3 lpm of the sample gas will enter the DMA and mix with a 3.0 lpm of sheath air, so the sample gas will experience a process of decreasing RH. Meanwhile, the relative humidity of the sheath air will increase, and then become stable after waiting enough time.

**Figure 7** shows the sample measurement results of aerosol particles size distribution and concentration curve at the inlet and outlet of bubble column. The RH of inlet laden gas was 30% which indicated the NaCl particles were dry crystal substance. The relative humidity of outlet gas represented the relative humidity of sheath air in SMPS. The RH of the sheath air exceeded the efflorescence point (42%) which indicated the NaCl

FIGURE 6 | Inlet and outlet aerosol concentration with or without aerosol supplied.

FIGURE 7 | Aerosol size distributions under different relative humidity of sheath air.

particles were metastable salt droplet because the outlet sample gas lies in reducing humidity process. It is obvious that the outlet aerosol size distributions shifted to right compared with the inlet aerosol size distribution. Besides, the larger the relative humidity of sheath air, the more apparent aerosol size distributions shifted to the right. It can be seen that the aerosol number concentration at the outlet is significantly higher than the inlet aerosol number concentration of the same aerosol particle size when the aerosol particle size is greater than 100 nm. This experiment's results indicate that during the process of the aerosol laden gas passing through the bubble column, the hygroscopic growth of soluble aerosol will lead to the increase of the aerosol particles size. So, it is obviously unreasonable to directly calculate the removal efficiency by using the data results in **Figure 6**.

To get the removal efficiency of aerosol particles with different sizes, the outlet aerosol size distribution and concentration under same RH with the inlet condition should also be achieved. One approach is adding a diffusion drying tube before the sample gas enters the SMPS to reduce the RH of sample gas

FIGURE 8 | Aerosol size distributions at inlet and outlet of bubble column with diffusion drying tube.

below the ERH. The circle line in **Figure 8** represents the outlet aerosol size distribution and concentration curve under dried condition. It can be seen that the aerosol size distribution and concentration curve in the inlet and outlet under dried condition have a similar peak particle size. Furthermore, there is no case that the outlet aerosol concentration is higher than the inlet aerosol concentration of the same aerosol particles size, which is reasonable and in agreement with the actual situation.

In view of the partial loss of aerosol particles when the aerosol laden gas passes through the diffusion drying tube, it is necessary to verify the performance of the diffusion drying tube. The performance test of the diffusion drying tube was also conducted under a gas flow rate of 4.0 lpm, and the result was shown in **Figure 9**. It can be seen that the aerosol size distribution and concentration curve at the inlet and outlet of the diffusion drying tube is close. When the aerosol number concentration is greater than 5000 per cubic centimeter, the partial loss fraction of aerosol particles is almost 5%. Therefore, to gain more accurate aerosol removal efficiency, it is necessary to modify the outlet aerosol concentration with a partial loss fraction.

According to the hygroscopic growth theory of soluble aerosol, growth factors of aerosol particles with different sizes can be calculated at given temperature and RH. Because the equilibrium time of nanometer aerosol particles is very short under RH<90%, it's reasonable to consider aerosol particles reach hygroscopic equilibrium in the DMA. Another approach to achieve the aerosol removal efficiency under high humidity condition is obtained by using the hygroscopic growth theory of soluble aerosol to modify the outlet aerosol particles size distribution and concentration curve. The hygroscopic growth theory of soluble aerosol described by equations (2)–(8) is used to modify the aerosol particles size distribution and concentration curve at the outlet under stable humidity environment. And the modified curve is compared with the curve measured by the diffusion drying tube (**Figure 10**). It can be seen that there is a subtle difference between the modified aerosol size distribution and concentration curve and the aerosol size distribution and concentration curve obtained by adding a diffusion drying tube when the RH of sheath air is 60.7%. It may be caused by the unstable humidity in DMA. While, the modified aerosol size distribution and concentration curve basically coincides with the aerosol size distribution and concentration curve obtained by adding a diffusion drying tube when the RH of sheath air exceeded 70% and maintain constant for more than 2 min, and all particles size was located at the inlet aerosol size distribution range. It shows that the aerosol particle size distribution and concentration curve modified by the aerosol hygroscopic growth theory can accurately represent the actual outlet aerosol situation under dry condition.

In order to further verify the feasibility of using a hygroscopic growth model to calculate soluble aerosol removal efficiency, aerosol removal efficiency calculated by the modified aerosol size distribution and concentration curve was compared with that directly, measured with the diffusion drying tube, and the results were shown in **Figure 11**. As can be seen, aerosol removal efficiencies obtained by the two methods are in good agreement with the relative error within 25% when the RH maintains constant. The large error at both sides of the curve is mainly because the low aerosol number concentration results in the large uncertainty of measurement results.

### CONCLUSIONS

Pool scrubbing is a key measure to capture soluble aerosol in case of severe accident at nuclear power plants. Under high humidity conditions such as pool scrubbing, the hygroscopic growth of soluble aerosol will lead to the increase of aerosol size, which further causes a large deviation in calculating the aerosol size-dependency removal efficiency. Kohler theory coupling shape correction factor related to the aerosol particles size can be used to accurately describe the hygroscopic growth of soluble aerosols in the equilibrium state. To account for the aerosol's actual size-dependency removal efficiency, two methods of calculating the removal efficiency were proposed. One method is directly calculating the removal efficiency by adding a diffusion drying tube to reduce the RH of sample gas below the ERH. Another method is modifying the aerosol size distribution and concentration curve by using the hygroscopic growth theory of soluble aerosol. The experiment results obtained by the two methods are in good agreement.

### AUTHOR CONTRIBUTIONS

YL and QM carried out experiments, analyzed the experimental results, and wrote the paper. YZ and HG assisted in experimental work and dealt with the experimental data. ZS and KS gave the guideline of this research and modified the language of the manuscript. YL and YZ developed the hygroscopic growth theory of soluble aerosol. YZ and KS designed the experimental loop and guided the experimental work.

### ACKNOWLEDGMENTS

The authors greatly appreciate the support of the Natural Science Foundation of China (Grant No. 11675045). The author thanks the Key Supported Discipline of Guizhou Province (Qian Xuewei He Zi ZDXK[2016]24), 2011 Collaborative Innovation Center of Guizhou Province (Qian Jiao he xietongchuangxin zi [2016]02). Additionally, the authors are thankful for the support of the Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, China.

### 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 © 2019 Li, Ma, Sun, Gu, Zhou and Shi. 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.

# Analysis of Categorical Subgroup Method for Resonance Self-Shielding Treatment

### Song Li, Zhijian Zhang\*, Qian Zhang, Chen Hao and Qiang Zhao

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

To meet the requirements of one step lattice calculation on resonance effect, a self-developed design and construction of a resonance treatment code are composed based on subgroup method and HELIOS-1.11 library. Subgroup fixed source equations are solved by method of characteristics to get subgroup fluxes, which are subsequently used to deduce effective resonance cross sections combined with subgroup weights and subgroup levels. Bondarenko method is employed to handle resonance interference effect and a resonance category scheme and resonance geometry simplification method are introduced to improve efficiency. Benchmarks of single pin cells and assemblies of light water reactor are adopted for numerical validation and the calculating results indicate that this method can treat resonance effect both precisely and effectively.

#### Edited by:

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

#### Reviewed by:

*Zhitong Bai, University of Michigan, United States Limin Liu, University of California, Berkeley, United States*

> \*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: *13 January 2019* Accepted: *30 April 2019* Published: *15 May 2019*

#### Citation:

*Li S, Zhang Z, Zhang Q, Hao C and Zhao Q (2019) Analysis of Categorical Subgroup Method for Resonance Self-Shielding Treatment. Front. Energy Res. 7:48. doi: 10.3389/fenrg.2019.00048* Keywords: resonance self-shielding, subgroup method, resonance interference effect, reactor physics, one-step method

### INTRODUCTION

In traditional lattice physics codes, the three-step method and pin-by-pin method are widely used in light water reactor calculation (Cacuci, 2010). Three-step method firstly adopts detailed numerical analysis on fuel cell or lattice scale, the regional homogenization is carried out subsequently and the last step is the diffusion or transport calculation for the whole core. By contrast, pin-by-pin method has two steps (Choi et al., 2017), namely assembly transport calculation and a reactor core calculation. However, three-step method and pin-by-pin method both consider the fuel cell or assembly as a unity, so nuclide reaction inside the fuel area cannot be taken into consideration in detail. To handle this problem, one-step method without any approximation for fuel cell and assembly has become the researching focus currently (Downar et al., 2016). One of the most significant parts of one-step method is treating resonance selfshielding problem effectively and precisely since it provides material cross section for the whole calculating process. Traditional methods for resonance treatment are equivalent method (Askew et al., 1965; Zhang et al., 2015), which is based on the equivalence between heterogeneous and homogenous problems and the ultra-fine group method with extremely detailed division of groups (Ishiguro and Takano, 1971). The former method has good efficiency but shows obvious drawbacks in accuracy and geometry adaptability, while although the latter one has a satisfactory accuracy, it's hard to be applied to lattice scale problems since the calculating burden is unacceptable. Other methods such as embedded self-shielding method and pseudo-isotope-method have been further developed in recent years, but problems for detailed region in fuel cell is still up in the air (Liu et al., 2018; Zhang et al., 2018).

In recent years, subgroup method is proved to have advantages in solving complex geometry problems with high accuracy and also has a good performance for sub-pin scale (Nikolaev et al., 1971; Li et al., 2018), so it has been one of the research hotspots of resonance treatment. Different from traditional method, subgroup method divides the subgroups by cross section value and the effective resonance cross section is calculated by subgroup parameters. By combining with transport module capable of handling problems with arbitrary geometry, subgroup method can be applied to any kind of lattice geometry. In this work, subgroup method combing with method of characteristics is adopted to treat resonance self-shielding effect.

Subgroup fixed source problems (SGFSP) are established to obtain subgroup fluxes (Jung et al., 2013), which are used subsequently to deduce effective resonance cross sections. For problems with various number of resonant isotopes, resonance interference effect is taken account by Bondaronko iteration method (Casal, 1991). In this procedure, while treating current resonant isotope, cross sections of other resonant isotopes are firstly considered as hypothetical values and then be updated by the values of the previous iteration step. It can be seen that SGFSP and resonance interference effect treatment occupy the most calculating time so optimizing method has been taken in this work to improve efficiency. In this case, we use a smaller number of subgroups in SGFSP for simplification while a larger number of subgroups in deducing effective cross section to ensure accuracy, and an extrapolation of relevant data obtained by SGFSP is used between these two set of subgroups. In addition, a resonance category is introduced to divide resonant isotopes into different categories and each category has been assigned with a representative isotope. SGFSPs are carried out according to representative isotopes and only interference of isotopes in the same category are taken into consideration. What's more, geometry setting of transport and resonance calculation are not necessarily to be the same and resonance treatment could adopt a simpler division with fewer sub-regions to predigest calculation burden.

In this work, a resonance calculating code is developed based on the method above and HELIOS-1.11 library. To verify the calculating results of this work, a series of benchmark including single cell and two-dimensional lattice released by Japan Atomic Energy Research Institute (JAERI) is employed for numerical verification and analysis (Yamamoto et al., 2002), and the calculating results indicate a good performance of method introduced in this work. Since this work provides effective material cross sections for the whole reactor physics calculation precisely and efficiently, it will serve as an important component for project of research on key technology of numerical reactor engineering in Harbin Engineering University. This work will also be coupled with transient, burn-up, and thermal hydraulic feedback module developed in this project for the further researches, which will not be illustrated here since we only concentrate on resonance calculation in this paper. This work is a redevelopment of HELIOS-1.11 program as they share the same neutronic functions. However, Different from HELIOS-1.11 code, which is written by outdated FORTRAN77 language, program of this work adapts C language for redesigning since it is more compatible to different operation systems and has the direct access to physical addresses as well as hardware operating. Besides, the resonance treatment procedure is coupled with an inhouse developed MOC transport code rather than CCCP module

used in HELIOS1-11. The independent design and construction of a completed lattice physics analyze code is realized.

### THEROTICAL MODEL

### Subgroup Method

Subgroup parameters consist of subgroup level and subgroup weight. As shown in **Figure 1**. The resonance group is divided into 4 subgroups and the subgroup levels indicate the 4 discrete cross sections from σ<sup>1</sup> to σ4, while subgroup weights w<sup>1</sup> to w<sup>4</sup> stand for the probability of the neutrons from this resonance group locating in each subgroup respectively. Take the absorption cross section for an example, if there are I subgroups in resonance groupg, effective cross section can be shown as follows:

$$\sigma\_{a,\emptyset} = \frac{\int\_{\Delta E\_i} \sigma\_{a,i}\left(E\right)\varphi\_i\left(E\right)dE}{\int\_{\Delta E\_i} \varphi\_i\left(E\right)dE} = \frac{\sum\_{i}^{I} \sigma\_{a,i}\varphi\_i \frac{\Delta E\_i}{E\_{\emptyset}}}{\sum\_{i}^{I} \varphi\_i \frac{\Delta E\_i}{E\_{\emptyset}}} = \frac{\sum\_{i}^{I} \sigma\_{a,i}\varphi\_i\omega\_i}{\sum\_{i}^{I} \varphi\_i\omega\_i} \tag{1}$$

In Equation (1), ϕ<sup>i</sup> is the subgroup flux. To obtain this value, subgroup fixed source equations which have the same format as Boltzmann transport equation are solved. SGFSP can be solved by any type of transport method, such as method of characteristics. The deduction process of this procedure is introduced as follows.

In light water reactors, energy of neutrons produced by fission reaction is relatively high, and nearly up to 99.5% of these neutrons are among the fast group range (Hebert, 2009). Therefore, it can be assumed that the fission source of resonance groups is negligible. In addition, up-scattering effect is also insignificant in resonance range, so the source term in this energy range only consists of down-scattering reaction from groups with higher energy. In this condition, for subgroup i in resonance group g, a fixed source equation can be written as follows:

$$\begin{split} \mathfrak{Q} \cdot \nabla \varphi\_{\mathfrak{F},i} \left( r, \mathfrak{Q} \right) &+ \Sigma\_{t, \mathfrak{F},i} \left( r \right) \varphi\_{\mathfrak{F},i} \left( r, \mathfrak{Q} \right) \\ &= \frac{1}{4\pi} \left[ \left( 1 - \lambda\_{\mathfrak{g}} \right) \Sigma\_{s, \mathfrak{g},i} \varphi\_{\mathfrak{F},i} \left( r \right) + \Sigma\_{b, \mathfrak{g}} \right] \end{split} \tag{2}$$

In Equation (2), λ<sup>g</sup> is the intermediate resonance factor and 6b,<sup>g</sup> is macroscopic background cross section and its definition is shown in Equation (3), in which M is the total number of nuclides types, Nmand 6p,g,mare the number density and potential scattering cross section of resonant nuclide m respectively.

$$
\Delta\Sigma\_{b,\emptyset} = \lambda\_{\emptyset}\,\Sigma\_{p,\emptyset} = \sum\_{m}^{M} \lambda\_{m,\emptyset} N\_m \sigma\_{p,\emptyset,m} \tag{3}
$$

In Equation (2), scattering cross section are made up by potential scattering cross section and resonance scattering cross section, while the latter is only non-zero for resonant isotopes. However, for simplification, resonance scattering cross section is also considered to be zero for resonant isotope in this work as resonance scattering integrals are not stored in the HELIOS library (Stammal'er, 2008). To compensate for this approximation, resonance absorption integrals in HELIOS library is adjusted during the producing process. Therefore, resonance scattering will not be taken into consideration and the final fixed source equation can be expressed as Equation (4). It can be seen that the source item of Equation (4) has no connection with flux, so it means that SGFSP can avoid the source iteration process of transport calculation, which will be much easier and more efficient for MOC module to solve. When subgroup fluxes are solved by transport module, effective

TABLE 1 | Absorption resonance integral of U-238 in 5.72–7.34 eV and its subgroup parameters.


#### (B) TWO SETS OF SUBGROUP PARAMETERS


resonance cross sections can be obtained by Equation (1).

$$\begin{aligned} \left[\boldsymbol{\Omega} \cdot \nabla \boldsymbol{\varphi}\_{\mathcal{G},i} \left(\boldsymbol{r}, \boldsymbol{\Omega}\right) \right. &+ \left[\boldsymbol{\Sigma}\_{a, \mathcal{G},i} \left(\boldsymbol{r}\right) + \lambda\_{\mathcal{G}} \boldsymbol{\Sigma}\_{\mathcal{P}, \mathcal{G},i} \left(\boldsymbol{r}\right)\right] \boldsymbol{\varphi}\_{\mathcal{G},i} \left(\boldsymbol{r}, \boldsymbol{\Omega}\right) \\ &= \frac{1}{4\pi} \lambda\_{\mathcal{G}} \boldsymbol{\Sigma}\_{\mathcal{P}, \mathcal{G},i} \left(\boldsymbol{r}\right) \end{aligned} \tag{4}$$

### Optimizations for SGFSP

Although solving SGFSP is much less time-consuming compared with Boltzmann transport equation, it still accounts for the vast majority of calculating time of subgroup method. Theoretically, subgroup method needs to carry out SGFSPs for each nuclide and each resonance group respectively, which is very time-consuming and accounts for the majority of the whole resonance calculating process. Therefore, the main purpose of the optimizations for SGFSP is to reduce the total number of subgroup fixed source equations, and the key part of this procedure is to reduce the total number of subgroups in all resonance groups.

Learned from the definition of subgroup method, each resonance group has a set of subgroup parameters and SGFSPs are carried out group by group. However, from Equation (4), we find that only subgroup levels are used in SGFSP, while subgroups weights are only used in the subsequent process of effective cross section calculation. For subgroups with the same level, no matter what values their subgroup weights are, they will have the same solution for SGFSP according to Equation (4). In this condition, the more subgroups have the same level, the less calculating burden it will be. If all resonance groups have the same set of subgroup levels, SGFSP would only need to be carried out in one resonance group.

In the procedure of generating traditional subgroup parameters (Cullen, 1977), the fuel is supposed to be composed of one resonant nuclide R and one moderator nuclideH. Therefore, the flux of the uniform condition can be written as follows

$$\varphi\_{\mathcal{S}} = \frac{\lambda\_R \Sigma\_{p,R} + \lambda\_H \Sigma\_{p,H}}{\left(\Sigma\_{a,R} + \lambda\_R \Sigma\_{\mathfrak{s},R} + \lambda\_H \Sigma\_{\mathfrak{p},H}\right)} = \frac{\sigma\_b}{\left[\sigma\_{a,R} + \lambda\left(\sigma\_{\mathfrak{s},R} - \sigma\_{\mathfrak{p},R}\right) + \sigma\_b\right]} \tag{5}$$

In Equation (2), 6<sup>a</sup> and σ<sup>a</sup> are absorption cross sections, 6<sup>p</sup> and σ<sup>p</sup> are potential scattering cross sections. In this condition, the unknown values in Equation (5) are only



subgroup levels and subgroup weights. As is discussed above, it should be arranged that every resonance group has the same set of subgroup levels but different subgroup weights. In this case, the subgroup level should be pre-specified by trial and error (Joo et al., 2009), and subgroup weights will be the only un-known value left. The multi-group libraries always store resonance integral data, which can be interpolated with background section and temperature. Take the absorption cross section for example, by selecting K resonance integrals of different background cross sections in the library, the subgroups weights with subgroup number of I can be obtained through optimum fitting method, which is shown in Equation (6)

$$\min F(\boldsymbol{w}) = \sum\_{k=1}^{K} \left[ \frac{RI\_a \left( \sigma\_{b,k} \right) - \sum\_{i=1}^{I} w\_{a,i} \sigma\_{a,i} \phi \left( \sigma\_{b,k} \right)}{RI\_a \left( \sigma\_{b,k} \right)} \right] \tag{6}$$

On the other hand, observed from **Figure 1**, the more subgroups number is adopted, the more precise to describe the resonance effect. However, it's not necessary to calculate SGFSP with

that many subgroups since subgroup fluxes are able to be interpolated (Park and Joo, 2018). Therefore, a set of subgroup parameters with fewer group number could be used in SGFSP to improve efficiency while another set of subgroup parameters with an extensional number of groups are used to calculate effective cross section to guarantee the accuracy. In this work, the number of subgroups in SGFSP is 4 and the latter one is set to be 7. To determine the subgroup levels, an iteration process is adopted. First, an initial set of subgroup levels is selected randomly, such as by geometric progression. Then the subgroup weighs can be generated by Equation (6). Through this set of subgroup parameters, resonance integral can be reproduced and its root-mean-square error (rmse) compared with the reference value is obtained. Subsequently, the values of subgroup levels are modified by a slight increase. If the new rmse also increase then the subgroup level should be adjusted to be diminished, vice versa, until the variations of rmse meet the converge criteria, which is set to be 0.05 and the maximum iteration number is 50. Take the absorption cross section of U-238 for an example, the variation of its resonance integral (RI) by background cross section (BCS) in energy range between 5.72 and 7.34 eV and temperature of 300 K is shown in **Table 1A**. By comparison, subgroup parameters of both 4 and 7 subgroups are given in **Table 2**.

Observed from **Table 1A**, absorption resonance integral of U-238 in this energy range shows a rapid rise by background cross section, which indicates a significant resonance peak occurs in this range. In **Table 1B**, M indicates for subgroup number of 4 and N indicates for subgroup number of 7. It is worthwhile to mention that if the sum of M or N subgroup weights are not equal to 1, another subgroup whose level is 0 is added and its weight is calculated by 1 − Pw<sup>i</sup> . What's more, it can be noticed that some subgroup parameters are not physical, such as 0 level, negative weights and the largest M and N level are both more than the actual largest cross section value in **Table 1A**. This phenomenon results from the procedure of Equation (6). Since it is not subgroup parameters values but the final effective cross section in Equation (1) that matters, the non-physical subgroup


values are acceptable as long as the accuracy of Equation (1) is ensured. After SGFSP is carried out by M levels, an extrapolation method is adopted to get relevant values for N subgroups. The detailed procedure is shown as follows.

$$\sigma\_a = \frac{\sum\_{i}^{I} \sigma\_{a,i} \omega\_i \frac{\sigma\_{b,i}}{\sigma\_{a,i} + \sigma\_{b,i}}}{\sum\_{i}^{I} \omega\_i \frac{\sigma\_{b,i}}{\sigma\_{a,i} + \sigma\_{b,i}}} \tag{7}$$

Combing Equations (1) and (5), the relationship between effective absorption cross section and background cross section can be deduced as Equation (7). After SGFSP is solved and subgroup flux is obtained, σb,<sup>i</sup> can be calculated respectively for 4 subgroups according to Equation (8). However, it is not σb,<sup>i</sup> that is used in interpolation. In practical heterogeneous system, there is an equivalence theorem which indicating σb,<sup>i</sup> should be augmented by an equivalence cross sectionσe,<sup>i</sup> . In this way, σb,<sup>i</sup> = λσ<sup>p</sup> + σe,<sup>i</sup> (Stammal'er, 2008). When σb,<sup>i</sup> is acquired, σe,<sup>i</sup> can be calculated just by subtracting λσ<sup>p</sup> from σb,<sup>i</sup> . Therefore, a set of σe,<sup>i</sup> with dependence of σa,<sup>i</sup> for 4 subgroups is established and the interpolation process is carried out between σe,<sup>i</sup> and ln σa,i .

$$
\sigma\_{b,i} = \frac{\sigma\_{a,i}\varphi\_i}{1 - \varphi\_i} \tag{8}
$$

### Resonance Interference Effect and Resonance Category

Method above is deduced on the basis that there is only one resonant isotope existing in the fuel area, which is not corresponding with actual condition. Traditional fuel used in light water reactors is composed of different kinds of uranium isotopes and in the process of burn-up, varieties of heavy isotopes will be generated and a considerable portion of them are resonant. Different resonance peaks of different isotopes will overlap with each other and this phenomenon is known as resonance interference effect.

Under these circumstances, while calculating the current resonant isotopes, the impact of resonance peaks of other resonant isotopes must be taken into consideration and Equation (7) is modified to capture this influence. Take the absorption cross section for an example, we define a pseudo cross section

σ<sup>x</sup> to represent the absorption contributions of the other resonant isotopes and the effective cross section can be shown as Equation (9).

$$\sigma\_{a} = \frac{\sum\_{i}^{I} \sigma\_{a,i} \omega\_{i} \frac{\sigma\_{b,i}}{\sigma\_{a,i} + \sigma\_{b,i} + \sigma\_{x}}}{\sum\_{i}^{I} \omega\_{i} \frac{\sigma\_{b,i}}{\sigma\_{a,i} + \sigma\_{b,i} + \sigma\_{x}}} \quad \text{with} \quad \sigma\_{x} = \frac{\sum\_{n \neq m} N\_{n} \sigma\_{a,n}}{N\_{m}} \tag{9}$$

In Equation (9), m represents the resonant isotope in concern while n means other resonant isotopes. However, since calculating σ<sup>x</sup> requires the absorption cross section of other resonant isotopes, whose value is unknown in the beginning, the Bondarenko iteration process (Cacuci, 2010) will be introduced in this process. At the very start, σ<sup>x</sup> is assumed to be zero in Equation (9), then the newly acquired σ<sup>a</sup> will be used in the second round of calculation and this process will be repeated to convergence. The convergence criterion is set that σ<sup>a</sup> should not be different more than 10−<sup>3</sup> compared with the previous iteration. The maximum number of iteration is set to be 30 and a relaxation factor of 0.5 is applied to accelerate the convergence process.

However, according to Bondarenko method above, each type of resonance isotopes should be taken into account and the number of iterations will sharply increase with the number of resonant isotopes. During the burn-up process, the calculating efficiency will become unacceptable. Therefore, there is a critical need for simplification of the resonance interference effect. In this work, according to the shape of resonance peak, different resonant isotopes can be divided into several categories (Stammal'er, 2008). Isotopes in the same category have similar property and the isotope with the most prominent absorption character is chosen as the representative isotope in each category. As U-238 is the most domain resonance isotope, it is classified as a single category containing only itself. Other Heavy mass isotopes and fission neutron productions are specified as a category due to its absorption or fission characters, and U-235 is assigned as its representative isotope. If Nature Zr is treated as resonance isotope, particular subgroup parameters can be used. However, in many cases resonance of Nature Zr does not matter much so only temperature interpolation is needed to get its cross section from pre-stored values. For other resonance isotopes such as control rod materials or strong absorbers like gadolinium isotopes, they are treated as another category and Hf-177 is assigned as representative isotope. The detailed resonance categories and representative isotopes are shown in **Table 2**. In SGFSP, only resonant isotopes of the same category will be treated and the influence of other categories is assumed to be

Frontiers in Energy Research | www.frontiersin.org

FIGURE 6 | Errors of U-235 absorption cross section of single cell. (A) UO2 pin, (B) MOX pin.

insignificant. In addition, isotopes of the same category all adopt subgroup parameters of representative isotope in SGFSP.

It can be noticed that cross sections in Equation (4) is group dependent, so there is a group-averaged calculation by weight of infinitely diluted resonance integrals. For regions with resonant isotopes, group-averaged calculation can be shown as the left side of Equation (10) while the right side for regions without resonant isotopes.

$$\Sigma\_{\mathcal{X}} = \frac{\sum\_{m \in k}^{M} \sum\_{\mathcal{S}} \Sigma\_{\text{xg}} N\_m R I\_{m\mathcal{S}, \infty} \Delta u\_{\mathcal{S}}}{\sum\_{m \in k}^{M} \sum\_{\mathcal{S}} N\_m R I\_{m\mathcal{S}, \infty} \Delta u\_{\mathcal{S}}}$$

$$or \Sigma\_{\mathcal{X}} = \frac{\sum\_{\mathcal{S}}^{G} \Sigma\_{\text{xg}} R I\_{r, \mathcal{S}, \infty} \Delta u\_{\mathcal{S}}}{\sum\_{\mathcal{S}}^{G} R I\_{r, \mathcal{S}, \infty} \Delta u\_{\mathcal{S}}} \tag{10}$$

In Equation (10), m ∈ k represents the isotopes belonging in the same category, RIm,g,∞is infinity diluted resonance integral of group for isotopem, 1u<sup>g</sup> is the lethargy width of resonance group g. On the right side, RIr,g,∞is the infinity diluted resonance integral for representative isotope of the category in calculation. In this condition, extrapolation variable ln σa,i is also adjusted by multiplying a modifying factor F toσa,<sup>i</sup> , which is shown in Equation (11).

$$F = \frac{\sum\_{\mathcal{G}}^{G} R\_{r, \mathcal{G}, \infty} \Delta u\_{\mathcal{G}}}{\sum\_{\mathcal{G}}^{G} R\_{m, \mathcal{G}, \infty} \Delta u\_{\mathcal{G}}} \tag{11}$$

### Geometry Setting and Program Design

In transport module, a fine mesh will be set to decompose the region and each mesh is considered to be a flat source region, in which the source is assumed to the same (Jung et al., 2013). As is shown in **Figure 2A**, each cell is divided into 4

TABLE 4 | Calculation burden of the 73rd path of UO2.


rings along the radial direction and each ring is also divided into 8 sectors. According to this geometry, there are 64 flat source regions containing resonant isotopes, where SGFSP and Bendarenko iteration will be carried out respectively. For regions only with non-resonant isotopes, no specific treatment is needed other than interpolation with temperature. However, on the beginning of life for a reactor, regions in the same fuel pin will have exactly the same materials. In this case, it can be assumed that resonance cross section of a certain isotope will only vary dramatically along radius due to space self-shielding, while regions with the same radius will have the same cross section. Therefore, geometry setting for resonance treatment can be simplified as shown in **Figure 2B**, which only reserve the region decomposition along the radius direction while sectors along the ring are neglected. Through this simplification, only 8 regions have resonant isotopes and numbers of iterations and memory usage will decrease dramatically. Calculating process of this part is shown as follows.

1) Based on the transport module, SGFSPs are solved according to geometry in **Figure 2A** to obtain fine mesh subgroup fluxes.


In this work, both resonant and non-resonant cross section treatment will be carried out after completing geometry treatment of transport module, which will provide all the geometry information needed in resonance treatment. Geometry information is defined in transport input file while the material information is defined in material input file separately. After SGFSP cross sections are set, fixed source solver of transport module will be called to obtain SGFSP fluxes, which will be used afterward in resonance interference and removal correction to obtain effective material cross sections. The final effective material cross section will be stored in variables defined by the transport module when the resonance treatment is completed. Transport module will continue to carry out the



eigenvalue calculation. By coupling with MOC and CMFD module, this work is capable to treat 2D or 3D scale problems. Moreover, as cross section treatment can be handled by each region independently, the parallel calculation based on domain decomposition is applied in this work. Flow chart of the whole calculating process is shown in **Figure 3**.

### NUMERICAL VALIDATION

This work adopts the 47 energy group structure of HELIOS-1.11 (Kim et al., 2015) based on the ENDF VI library, in which the resonance range is from 1.855 to 9118 eV. Benchmark used in this work is selected from light water reactors which contains UO<sup>2</sup> and MOX fuels with both fuel pin cell and PWR fuel assembly (Yamamoto, 2004) and subgroup method illustrated in this work will be used to calculate these problems. It is worthwhile to mention that HELIOS-1.11 program is applied as a code-to-code benchmark for the verification of the new developed code of this work. The reference values of pin cell and lattice problems are both provided by the report issued by Japan Atomic Energy Research Institute (Yamamoto et al., 2002).

### Single Cell Problems

Single cell problems contain two typical cases. One is a typical UO<sup>2</sup> pin and the other is a MOX pin, number densities of isotopes contained in these two problems are shown in **Table 3**. These two cases have the same geometrical configuration and dimension, which is shown in **Figure 4**. In addition, the burn-up processes of these two problems are also analyzed.

In the first problem, resonant isotopes include only U-238 and U-235 while the latter one also has plutonium and actinium isotopes. It has a great increment of isotopes accumulated during the burn up process and more than 100 isotopes would be generated, and more than 20 of them are resonant isotopes. Pin-averaged absorption cross sections of U-238 and U-235 of these two cases and Pu-239 and Pu-241 of MOX pin are shown in **Figures 5**–**7**, respectively. **Figure 8** gives the eigenvalue calculated during the whole burn-up process with 73 steps. The relative error of cross sections and error of eigenvalue are calculated by Equation (12) and Equation (13) respectively, where subscript x stands for values calculated by subgroup method while ref stands for reference value.

$$Relative\ Error = \frac{\sigma\_{\text{x}} - \sigma\_{ref}}{\sigma\_{ref}} \times 100\% \tag{12}$$

$$\text{Eigenvalue Error} = \left(k\_{\text{x}} - k\_{\text{ref}}\right) \times 10^5 \tag{13}$$

Learned from **Figure 5**, improved subgroup method (ISGM) introduced in this work has a good performance to describe the variations of U-238 absorption cross section both in UO<sup>2</sup> and MOX pin, of which the largest relative error is <0.6%. It should be noticed that cross sections with energy range from 5∼7 to 10∼100eV, especially in energy around 6.7 eV, have significant resonance peaks, so that its error shows a moderate rise compared with other energy range. In general, ISGM precisely captured the dramatic change of cross sections as well as resonance interference effect of U-238.

From **Figures 6**, **7**, it can be seen that ISGM has an obvious improvement in figuring up cross sections for U-235 and plutonium isotopes compared with U-238. It is reasonable since resonance peaks them are not as remarkable as those of U-238. The largest relative error of U-235 absorption cross section of UO<sup>2</sup> and MOX pin are both <0.3% and errors are no more than 0.1% in most of the energy range. For Pu-239 and Pu-240, errors slightly tend to be upward compared with U-235 since their resonance peaks overlap more intensely with U-238. However, the largest error of plutonium isotopes is still <0.5%, indicating a good performance of ISGM for resonance treatment.

Eigenvalue of the burn-up process is shown in **Figure 8**. In this process, as the number of resonant isotopes increases


dramatically, influence of resonance category is analyzed. In **Figure 8**, SGM indicates for traditional subgroup method without resonance category and ISGM has adopted resonance category shown in **Table 2**. It can be seen that error for k effective is <10 pcm for UO<sup>2</sup> pin in the beginning for both SGM and ISGM. As the burn-up goes on, errors have an increasing trend with a downwards in the end. The largest error in the whole burnup process is 29 pcm for ISGM while −98 pcm for SGM, which indicates both methods in this work has a good performance to treat UO<sup>2</sup> problems and its burn-up process. As for MOX pin, the error becomes a little bit of larger since resonance isotopes increase. The largest error occurs at the 19th burn-up step with 60 pcm for ISGM while −106 pcm for SGM at the last step, which is also acceptable. ISGM only takes resonant isotopes in the same category into consideration, so that influence of resonance peaks of other resonant isotopes is neglected, making values of cross section smaller than that of SGM and this accounts for the phenomenon that eigenvalue of ISGM is larger compared with SGM. **Table 4** compares the calculating burden of SGM and ISGM in the 73rd burn-step of UO2. It can be seen that ISGM can make a considerable reduction of the number of SGFSP equations and also shows an improvement for resonance iteration number. The total time for resonance treatment of ISGM is reduced by 35.5% compared with SGM while the accuracy is also guaranteed.

### Assembly Problem

Assembly problems consist of a PWR UO<sup>2</sup> fuel assembly and a MOX fuel assembly, which are the same geometrical configuration as a traditional 17 × 17 type PWR fuel design. The UO<sup>2</sup> assembly is composed of UO<sup>2</sup> and UO2-Gd2O<sup>3</sup> fuel rods. The geometrical description and the configuration of the assembly geometry are given in **Figure 9A**. The MOX fuel assembly is the same geometrical configuration as the UO<sup>2</sup> fuel assembly. Different from MOX single cell, the assembly is composed of low, middle, and high Pu content fuel rods, which is shown in **Figure 9B**. Materials of guide tubes and thimble tubes of these two assemblies are the same as clad. Besides, the clad and moderator material are the same as those of single cell problem. In addition, UO<sup>2</sup> in the first assembly is the same as that of single cell. Isotope number densities of UO2-Gd2O<sup>3</sup> and 3 types of MOX are shown in **Table 5**. The aim of these lattice problems is to examine the ability of the subgroup method used in this work to capture the resonance effect in strong absorbers and large scale problems. In the calculating procedure, the boundary condition is set to be reflective in all four sides.

Due to space self-shielding, fluxes would suffer a sharply decrease from the pin surface to the inner side, making resonance cross sections change dramatically along the radial direction even with the same number density. To describe the detailed cross section distribution, the Gadolinium pin is divided into 8 equal volume rings along the radial direction while the UO<sup>2</sup> pins have 3 rings. Each ring is also divided into 8 equal volume sectors in the transport module as fluxes would vary in different region of the same ring because of the asymmetric configuration of the lattice. However, as illustrated before, sectors in the same ring would have the same cross sections despite in reality they are different. To analyze the influence of this simplification, radial distribution of resonance cross sections from 5.72 to 7.34 eV of U-238 in Gadolinium pin and high Pu pin are shown in **Figure 10**. In this figure, SGM represents for subgroup method with the same geometry as transport module while ISGM represents for resonance effect is calculated in ring-wise scale. In addition, both SGM and ISGM adopts resonance category introduced in **Table 2**. The reference value is given by fine mesh identical to SGM, which is more precise theoretically.

For Gadolinium pin, resonance cross sections will change more acutely due to the existence of strong absorbers, especially for regions near the pin surface. It can be seen that even in the

TABLE 7 | Calculation burden of UO2 lattice.


same ring, cross sections of the outer-most side vary sharply. SGM can describe this phenomenon accurately and its cross section curve is almost entirely coincided. Errors of SGM are <0.3% in all regions. As for ISGM, it is more like an average value of cross section in the ring, cross sections in the sector-wise region will have a larger relative error compared to the fine-mesh reference value. Calculating results for high-Pu pin have the same tendency, in which errors of SGM are <0.3% while the largest error of ISGM will be more than 5%. The influence of these errors to eigenvalue and power distribution is shown in **Table 6** and **Figures 11**, **12**, respectively.

Compared with the reference eigenvalue, errors of SGM and ISGM are both <100 pcm, indicating these two calculating options can both give an accurate effective multiplication factor in spite of different cross sections. Particularly, difference between eigenvalues of SGM and ISGM is negligible as only 1 pcm for both UO<sup>2</sup> and MOX assembly. From **Figures 11**, **12**, errors of pin power distribution are also slight as the largest one is <1%. For UO<sup>2</sup> assembly, pins around Gadolinium pin have a relatively higher error as strong absorber is harder to deal with in referred to resonance effect. For MOX assembly, relative error for pins in the outer side tend to rise as absolute values of pin power are smaller in this region. In addition, differences of pin power distribution between SGM and ISGM are also insignificant with only 0.1% deviation. In general, accuracy for eigenvalue and pin power distribution of SGM and ISGM are the same satisfying for lattice physics calculation. **Table 7** compares the calculating burden of SGM and ISGM in terms of time-consuming and memory cost. It can be seen that ISGM has an obvious improvement in both two aspects as resonance interference effect and geometry setting are simplified. Only 7.8% of total time and 13.5% of total memory is needed for ISGM to provide effective region-wise macro cross sections while accuracy is also assured. In this case, ISGM is a more appropriate method to handle the resonance effect both accurately and effectively.

### REFERENCES


### CONCLUSIONS

Subgroup method based on HELIOS−1.11 library is adopted to develop the resonance calculating code in this work, which is capable to satisfy the requirement for resonance treatment in one-step reactor lattice physics calculation. This method could be combined with the transport module for arbitrary geometry to get subgroup fluxes for effective cross sections deduction. Through Bondarenko iteration method, resonance interference effect is taken into consideration. In addition, a simplified set of subgroup parameters is used subgroup fixed source problems to improve efficiency while a detailed subgroup flux is interpolated to secure accuracy. Moreover, a resonance category is introduced to simplify the process of resonance interference treatment and a resonance geometry setting is analyzed to reduce the calculating burden of time and memory usage. Series of simplified and original subgroup levels and the integrated resonance category are detailed illustrated in this work. The main innovation point of this work is the realization of independent research and redevelopment of a resonance treatment program, which is capable to incorporate with the in-house transport and burn-up code to compose a complete lattice physics analyzing program. The performance of resonance treatment program is revalidated by a series of publically released benchmarks. Numerical results show that methods used in this work have a good performance in calculating resonance cross sections, the eigenvalue and pin power distribution accurately and effectively.

### AUTHOR CONTRIBUTIONS

SL, QZhang and CH contributed conception and design of the study. SL organized the database, performed the statistical analysis, and wrote the first draft of the manuscript. ZZ and QZhao decide the research contents of this work and give hardware-software support. All authors contributed to manuscript revision, read and approved the submitted version.

### ACKNOWLEDGMENTS

This work is supported by the Research on Key Technology of Numerical Reactor Engineering [J121217001], the Fundamental Research Funds for the Central Universities [GK2150260154], the Heilongjiang Province Science Foundation for Youths [QC201803] and China Scholarship Council.

Choi, S., Smith, K. S., Kim, H., Tak, T. and Lee D. (2017). On the diffusion coefficient calculation in two-step light water reactor core analysis. J. Nuclear Sci. Technol. 54, 705–715. doi: 10.1080/00223131.2017.1299648


Stammal'er (2008). HELIOS Methods. Studsvik Scandpower


**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 © 2019 Li, Zhang, Zhang, Hao and Zhao. 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.

# Effect of Liquid Injection Arrangements on Injection Flow Rate of a Laboratory-Scale Venturi Scrubber

Guangzong Zheng1,2, Bo Chen<sup>3</sup> , Jishen Li <sup>1</sup> , Kaiyi Shi <sup>4</sup> and Haifeng Gu<sup>1</sup> \*

*<sup>1</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China, <sup>2</sup> China Institute of Atomic Energy, Beijing, China, <sup>3</sup> China Nuclear Power Engineering Co., Ltd., Beijing, China, <sup>4</sup> Department of Chemistry and Chemical Engineering, LiuPanshui Normal University, LiuPanshui, China*

#### Edited by:

*Muhammad Zubair, University of Sharjah, United Arab Emirates*

#### Reviewed by:

*Yacine Addad, Khalifa University, United Arab Emirates Keyou S. Mao, Purdue University, United States Francesco Di Natale, University of Naples Federico II, Italy Khurram Mehboob, King Abdulaziz University, Saudi Arabia*

> \*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: *27 November 2018* Accepted: *16 May 2019* Published: *04 June 2019*

#### Citation:

*Zheng G, Chen B, Li J, Shi K and Gu H (2019) Effect of Liquid Injection Arrangements on Injection Flow Rate of a Laboratory-Scale Venturi Scrubber. Front. Energy Res. 7:51. doi: 10.3389/fenrg.2019.00051* As the core device of a filtered containment venting system (FCVS), Venturi Scrubber is an efficient device to scrub the radioactive gases and aerosols before release into the atmosphere. The design concept of Multi-Venturi Scrubber System makes the laboratory-scale venturi scrubber research more valuable. This paper studied the injection flow rate of Venturi Scrubber in different injection arrangements. The liquid is injected horizontally and vertically to the throat at different radial position of (i) at the wall of throat, (ii) at the half radius of throat, and (iii) at the center of Venturi Scrubber throat with a nozzle of diameter 4 mm. Throat gas velocities range from 0 to 190 m/s. A constant level water tank was installed to keep water level constant during the injection process. The results showed that liquid injection modes significantly affect the injection performance. The arrangements of straight tube at center and elbow tube at center had larger injection flow rate among the others, and the injection flow rate increased as the throat gas velocity increased. The conventional wall opening (i.e., straight tube at the wall) injection method had the worst injection performance. This study provides a valuable reference for the liquid injection arrangement and structural design of the Venturi Scrubber.

Keywords: Venturi Scrubber, injection arrangement, injection flow rate, filtered containment venting system, FCVS

## INTRODUCTION

As the core device of a filtered containment venting system(FCVS), Venturi Scrubber is an efficient device to scrub the airborne source items before release into the atmosphere. The radioactive airborne source items mainly contain aerosol, iodine and methyl iodide (Eckardt and Losch, 2012). A Venturi Scrubber has convergent section, throat section, and divergent section. Due to small cross-sectional area of the throat, the air velocity reaches the maximum at the Venturi Scrubber throat. Under the action of the shear force of the high-speed airflow, the injection water is atomized into numerous small droplets, which provide a sufficient surface area for dedusting (Hills, 1995; Das and Biswas, 2006). The airborne sources contact with the droplet surface and settle inside the droplets. This process contains many mechanisms, including inertial collision, interception capture, diffusion capture, gravity sedimentation, and electrostatic adsorption (Pulley, 1997; Ali et al., 2013).

Venturi Scrubbers are widely used in industry because of their simple structure and high dust removal efficiency (Rudnick et al., 1986; Talaie et al., 1997; Ahari et al., 2008). Early studies of Venturi Scrubbers focus on internal two-phase flow, including jet characteristics (submergence depth, jet velocity, jet diameter, and jet trajectory, etc.) and droplet characteristics (initial droplet concentration, droplet size, and droplet velocity, etc.). In engineering, researches focus on the effect of operating parameters and structural parameters on injection flow rate and filtration efficiency. The operating parameters include gas flow rate, operating pressure and submergence depth, etc. The structural parameters contain throat length, divergent section angle and the placement of the ejector tube, etc. Many people studied the different structural parameters of the Venturi Scrubber including area ratio, projection ratio, divergent angle, throat length, and convergent section shape, to explore the performance of the Venturi Scrubber (Kroll, 1947; Panchal et al., 1991; Cramers and Beenackers, 2001; Sriveerakul et al., 2007; Yadav and Patwardhan, 2008). Das and Biswas (2006) proposes that the filtration efficiency of the Venturi Scrubber is closely related to its structure design including the suction chamber, the mixing throat, the divergent diffuser, the forcing nozzle, and divergence angle, etc. Zhou et al. (2015) explores the effect of operating conditions and structures on the removal of aerosols in the Venturi Scrubber. The study finds that the retention efficiency in a Venturi Scrubber increases with the increase of both gas velocity and injection flow rate, and the influence of gas velocity on efficiency is more effective at low injection flow rate. Moreover, the Venturi Scrubber with a long throat length or small diffuser angle performs excellent retention performance for small size aerosols.

Westinghouse offers the FILTRA-MVSS (Multi-Venturi Scrubber System) to effectively mitigate the consequences of a severe reactor accident by signifcantly reducing the level of radioactive release to the surrounding environment (Westinghouse Electric Company, 2012). The number of Venturi Scrubbers in operation is determined by the actual mass flow. This design method can maintain high throat gas velocity and injection flow rate of each Venturi Scrubber.

The aim of this paper is to present some experimental analysis of the effect of different liquid injection arrangements on injection flow rate of a self-priming Venturi Scrubber. The experiment is carried out using two types of tubes, straight tubes and elbow tubes. The liquid is injected into the throat of the Venturi Scrubber parallel or perpendicular to the airflow at different radial position of (i) at the wall of throat, (ii) at the half radius of throat, and (iii) at the center of Venturi Scrubber throat. Firstly, the injection flow rate of 6 different injection methods is measured under several operating conditions. Then, the range of working conditions are expanded for the configurations of straight and elbow tubes at center, which have advantages in injection performance. And the influence of throat gas velocity and submergence depth are investigated. As a comparison of the above two types of injection, the conventional wall opening configuration(straight tube at the wall of throat) is also studied at the same time. This study provides a valuable reference for the liquid injection arrangement and structural design of the Venturi Scrubber.

### EXPERIMENTAL DETAILES

### Experimental Loop

A schematic illustration of the experimental loop is shown in **Figure 1**. The system is mainly composed of screw air compressors, air tank, air filters, different types of valves, mass flowmeter, Venturi Scrubber, constant level tank, weigh, and lifter. The air compressors provide a stable gas source to the system with flow rate of 10∼1,200 kg/h. After the gas flow is buffered in the air tank and enters the experimental pipeline, it flows through air filters, shutoff valve, mass flowmeter, throttle valve, check valve, and finally into the Venturi Scrubber. The injection water from the constant level water tank is atomized into small droplets under the effect of high-speed air shear force in the Venturi Scrubber throat, which provides a sufficient surface area for dust removal.

The water tank designed in this experiment can maintain a constant water level while supplying water for the Venturi Scrubber. The injection flow rate can be obtained by measuring the water tank weight at different times. When the water tank provide liquid to the Venturi Scrubber, it calls for the valve➀ closed, the valve➁ and valve➂ open. At this time, the right chamber of the water tank vents to the atmosphere while the left chamber does not. The water level on the right side of the Venturi Scrubber maintains constant while the left side level drops gradually. The difference in height between the right chamber water level to the injection orifice of the Venturi Scrubber becomes submergence depth, supplying power to liquid injected. The constant level tank, connecting to the injection tube of the Venturi Scrubber through rubber hose, is placed on the weighing device with weight measurement accuracy to 0.1 g, and can realize submergence depth adjustable by manipulating the lifter.

### Test Section

The equipment used in the experimental tests was a Venturi Scrubber constructed of acrylic, with a circular cross-section and throat dimensions of diameter 24 mm. Liquid was injected into the throat through six different nozzle arrangements. The main dimensions and nozzle arrangements of the Venturi Scrubber are shown in **Figure 2**. The gas velocity at the Venturi Scrubber throat is evaluated in the present study as follows.

$$U\_{\mathfrak{g},th} = \frac{M\_{\mathfrak{g}}}{\rho\_{\mathfrak{g},th} A\_{th}}$$

M<sup>g</sup> is the gas mass flow rate. ρg,th is the throat gas density and Athrepresents the throat cross-sectional area. Contrary to conventional method with orifice in throat, the liquid was injected horizontally and vertically to the throat with nozzles of 4 mm diameter at different radial positions. The nozzle positions were located on the wall surface, 0.5 times radius and the central axis of the Venturi Scrubber throat.

### Experimental Procedure

Air was supplied by the air supply line and the flow rate could be measured. After the air compressor achieved desired pressure 0.7 Mpa, the shutoff valve and throttle valve were opened in turn. The gas flow rate was adjusted to a certain value by controlling the throttle valve. Gas flow rate values were measured by mass flowmeter and recorded by computer.

The constant level water supply system was capable of supplying water of a constant submergence depth to the Venturi Scrubber and measuring the injection flow rate. When offering water to the Venturi Scrubber, valve➀ was closed, and valve➁➂ were open. The liquid level in the right chamber was constant and parallel to the small holes. The bubble entered the left chamber through the small holes, and the liquid level in the left chamber gradually decreased. The injection flow rate can be calculated by measuring the weight difference of the water tank and time. Operating parameter ranges are shown in **Table 1**.

### RESULTS AND DISCUSSION

The effect of liquid injection arrangements was carried out using two types of tubes (straight tubes and elbow tubes),

#### TABLE 1 | Experimental operating parameters.


and three radial positions (i.e., wall, 0.5r and center of the Venturi Scrubber throat). Section Effect of Radial Positions and Type of Injection Tubes involves the effect of radial positions and type of injection tubes. After that, three particular nozzle arrangements (i.e., straight tube at the wall, straight tube at the center, and elbow tube at the center) concerned were further studied, considering the influence of throat gas velocity and submergence depth. Injection characteristics of the three nozzle arrangements were studied in section The Influence of Throat Gas Velocity under different throat gas velocity, which, respectively, took place at three submergence depths. Section The Influence of Submergence Depth focuses on the influence of varied submergence depths.

### Effect of Radial Positions and Type of Injection Tubes

**Figures 3**, **4** show the effect of radial positions and type of injection tubes on injection flow rate under particular throat gas velocity and submergence depth. It is found that the injection flow rate, with the same nozzle radial position but two different types of tubes, straight and elbow tubes, shows an increasing trend as nozzle radial positions vary from wall, 0.5r and the center. This is due to the highest gas velocity at the center of the throat, and the lowest velocity at the wall due to the presence of the boundary layer. Correspondingly, the center of the throat produces a greater negative pressure than the wall surface, providing more power for injection water. Although the nozzles locating at the center always prevail in injection flow rate, it still remains a question which tube, straight or elbow, acts better under certain operating parameters. In **Figures 3B,C**, **4C**, the elbow tubes, compared to straight tubes, prevail in injection flow rate at the same position, while in **Figure 4A**, the elbow tubes are at a disadvantage in injection flow rate. Particularly, the injection flow rate curves of the two types of tubes cross in **Figures 3A**, **4B**. It can be clearly seen that whether straight tubes or elbow tubes prevail in injection flow rate is a complicated process, which is an effect of the coupling of throat gas velocity and submergence depth. Still, many implicit rules can be obtained after analysis. It reveals that, as throat gas velocity increases, the injection flow rate of elbow tubes at the same radial position does not dominate at first and gradually exceeds the straight tube, as is shown in **Figures 4A–C**. **Figures 3A–C** shows the same law as well. Comparing injection flow rate (nozzle at the throat center) in **Figures 3B**, **4B**, it can be seen that with the increase of the submergence depth from 2 to 42 cm, the injection flow rate of straight tube does not dominate in the beginning and finally exceeds the elbow tube. It indicates that, although the effect of injection flow rate is a complicated process of coupling of throat gas velocity and submergence depth, the increase of the throat gas velocity is more conducive to the elbow tubes, while the increase of submergence depth is more conducive to the straight tubes.

### The Influence of Throat Gas Velocity

Due to the study in section Effect of Radial Positions and Type of Injection Tubes on six different nozzle arrangements, in which the injection flow rate of two types of tubes at the center always dominates, it is necessary to do further research on the central straight tubes and elbows. The configuration of the straight tube at the wall is also taken into consideration as a reference, which is mainly due to the fact that the straight tube at the wall is a structure commonly adopted by the traditional Venturi Scrubber. **Figures 5A–C** show the injection flow rate curves of three nozzle arrangements (i.e., straight tube at the center, elbow tube at the center, and straight tube at the wall) as throat gas velocity changes, which are developed at submergence depths of 2, 12, and 42 cm, respectively. It can be clearly found that, with the increase of throat gas velocity, the injection flow rate of the straight tube and elbow tube at the center is significantly higher than that of the straight tube at the wall, and the greater the throat gas velocity is, the more obvious the advantages will be. When the throat gas velocity is 190 m/s and the submergence depths are 2, 12, and 42 cm, the elbow tube at the center, compared with the traditionally adopted straight tube at the wall, can increase the injection flow rate by 2.6, 2.3, and 2 times, respectively. Meanwhile, it can be found that under the condition of global throat gas velocity, the injection flow rate of straight tube at the center is always larger than the straight tube at the wall. This meets and reinforces the discovery in section Effect of Radial Positions and Type of Injection Tubes that for the same type of tube, as nozzle position moves closely to the throat center, the injection flow rate increases.

The detail that cannot be ignored is that when the throat gas velocity is very low, compared with the straight tube at the wall, the injection flow rate of the elbow tube at the center is not dominant. As the throat gas velocity increases, the elbow tube at the center begins to reverse the straight tube at the wall, thus resulting in an intersection of the curves. And the intersections move backward with the increase of submergence depth, that is, as the submergence depth increases from 2, 12–42 cm, the intersection occurs at throat gas velocity of 14, 31, and 59.5 m/s, respectively. This is due to the shape resistance of the elbow, so that when the throat gas velocity is low, the injection flow rate of the elbow tube at the center is lower than that of the straight tube at the wall under the same operating parameters. It can be seen from **Figures 5A–C** where the throat gas velocity is 0, that the difference in the injection flow rate between the two types of nozzle arrangements is more pronounced as the depth of submergence increases. When the submergence depth is 42 cm, a difference of 960 g/min between the two arrangements is generated, which is due to the increase of submergence depth, making the effect of the elbow tube shape resistance more

prominent. Accordingly, greater throat gas velocity is need as a counteracting of elbow tube shape resistance, to make the elbow tube at the center overtake the straight tube at the wall, which also explains why the previously mentioned increase of submergence depth makes the backward moving of intersection.

Comparing the straight tube and elbow tube both at the center, it is found that when the throat gas velocity is low, the elbow tube is not dominant. As the throat gas velocity gradually increases, the injection flow rate of the elbow tube gradually exceeds that of the straight tube. The intersections in **Figure 5A** wherein the submergence depth is 2 cm are more complicated, especially when the throat gas velocity is in the range from 24 to 71 m/s. In **Figures 5B,C**, however, intersections occur at throat gas velocity of 59.5 and 124 m/s, respectively. This is also the result of the coupling effect of the elbow tube shape resistance, submergence depth and throat gas velocity.

It is worth noting that by observing the injection flow rate curve of the straight tube at the wall in **Figures 5A–C**, it is found that the injection flow rate shows different trends at different submergence depths. When the submergence depth is low, i.e., 2 cm, the injection flow rate generally increases with the increase of throat gas velocity, while the injection flow rate decreases when the submergence depth is high, i.e., 42 cm. The overall injection flow rate change is small when the submergence depth is 12 cm. This signals that for a straight tube at the wall, increasing the throat gas velocity does not necessarily contribute to the increase of injection flow rate. Nevertheless, from the point of view of the injection flow rate under different submergence depths, the injection flow rate under the high submergence depth always dominates, despite the different trends under different submergence depths. This newly discovered phenomenon, wherein the underlying causes are related to internal two-phase flows and flow fields, needs further exploration.

### The Influence of Submergence Depth

The above studies were conducted at several submergence depths, and the injection flow rate often showed different trends under different submergence depths. From the above studies, it is also found that, with different throat gas velocity, the submergence depth seems to always promote the increase of injection flow rate. Since the above studies are merely conducted at several

submergence depths, it is necessary to study the conditions at more depths of submergence. **Figures 6A–C** show the injection flow rate trend with the submergence depth in three nozzle arrangements (i.e., straight tube at the center, elbow tube at the center, and straight tube at the wall), at an throat gas velocity of 48, 119, and 167 m/s. It can be clearly seen that the injection flow rate of different nozzle arrangements increases as the submergence depth increases. However, the injection flow rate curves in **Figure 6A**, compared with those in **Figures 6B,C**, are not significantly different from each other, while the injection flow rate of straight tube and elbow tube both at the center in **Figures 6B,C**, compared to that of straight tube at the wall, has obvious advantages.

In **Figure 6A**, the injection flow rate of the straight tube at the center, despite the slight advantage, is generally higher than that of the elbow tube at the center and the straight tube at the wall. While the injection flow rate of elbow tube at the center, compared to that of straight tube at the wall, prevails at first and then is overtaken by that of straight tube at the wall when submergence depth reaching 22 cm. The above phenomena are caused by the radial pressure distribution generated by gas flow in the Venturi Scrubber throat and elbow tube shape resistance. The shape resistance of the elbow tube hinders the flow of water in the elbow tube, wherein the process is affected by the submergence depth and the throat gas velocity. That is, the higher the submergence depth is, the more obviously the effect of shape resistance behaves, and the greater the throat gas velocity is, the less the effect becomes. In **Figure 6A**, the throat gas velocity is relatively low (48 m/s), the radial pressure distribution at the Venturi Scrubber throat is not obvious enough. Nevertheless, the injection flow rate of the straight tube at the center, due to the radial pressure distribution, is still higher than that of straight tube at the wall.

As for the elbow tube at the center, the injection flow rate of which is lower than that of the elbow tube at the center, the shape resistance emerges in it, so that the elbow tube does not dominate over the straight tube under the same operating parameters. The

existence of the elbow tube shape resistance also explains the phenomenon that after the submergence depth exceeds 22 cm, the injection flow rate of the elbow tube at the center is overtaken by that of the straight tube at the wall. Before the submergence depth reaches 22 cm, the radial pressure distribution of the throat make the elbow tube prevail in injection flow rate, despite the relatively low throat gas velocity. As the submergence depth grows, shape resistance of the elbow tube start to merge, resulting in the injection flow rate of central elbow tube being overtaken by straight tube at the wall. It should be noted that the throat gas velocity at this time is still small, equaling to 48 m/s. When throat gas velocity is 119 or 167 m/s, as is shown is **Figures 6B,C**, raising the submergence depth does not make the injection flow rate of the straight tube at the wall catch up with that of elbow tube at the center. This is due to the large throat gas velocity hinders the effect of elbow tube shape resistance, and increasing the submergence depth within the experimental range (2–42 cm) is insufficient to make the elbow tube shape resistance effect appear.

It should be pointed out that the above comparison has a nonnegligible disadvantage for the straight tube at the wall. That is, the negative pressure at the wall caused by the radial pressure distribution is significantly lower than that at the center, so that the injection power at the wall is lower than that at the center. When comparing the straight tube and elbow tube both at the center, as is shown in **Figure 6B** where the throat gas velocity is 119 m/s, the injection flow rate of central straight tube reverses that of elbow tube as the submergence depth rises to 37 cm. The same phenomenon occurs in **Figure 6C** where the throat gas velocity is 167 m/s. When raising the submergence depth to 47 cm, the injection flow rate of the straight tube at the center surpasses that of the elbow tube. Thus, the comparison of straight tube and elbow tube both at the center also proves the existence and role of the elbow tube shape resistance.

### CONCLUSIONS

In this study, effect of six nozzle arrangements on injection flow rate were studied. Particularly, three nozzle arrangements, i.e., straight tube at the center, elbow tube at the center and straight tube at the wall were further examined, considering the influence of varied operating parameters of throat gas velocity and submergence depth. Conclusions are summarized as follows.


### AUTHOR CONTRIBUTIONS

GZ, as the main executor of the project, participated in the facility construction, experiment, analysis, and thesis writing. BC, as a number of my lab and now graduated, help me analysis the experimental results, including different nozzle arrangements' effect on injection. JL, as a member of my lab, assisted me in the whole process, including the reading of experimental data and the construction of the test bench. KS, help me build data acquisition system. Many of his suggestions helped me improve the design of the Venturi scrubber. HG, my teacher, also expert in the field of filtered containment venting system, give me some key guidance on the experimental process.

### FUNDING

These authors are profoundly grateful to the financial supports of the Fundamental Research Funds for the Central Universities (NO. HEUCFM181203). 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).

### REFERENCES


### NOMENCLATURE


**Conflict of Interest Statement:** BC was employed by company China Nuclear Power Engineering Co., Ltd.

The remaining 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 © 2019 Zheng, Chen, Li, Shi 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(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.

# A Numerical Investigation on Gaseous Stratification Break-Up Phenomenon of Air Fountain Experiment by Code\_Saturne

Chunhui Zhang\*, Jiyang Yu, Xianping Zhong and Muhammad Saeed

*Department of Engineering Physics, Tsinghua University, Beijing, China*

Since the Fukushima nuclear power plant accident, the nuclear safety experts devote a number of efforts to explore the potential risk of hydrogen detonation and measures for the mitigation of hydrogen safety during a severe accident in the nuclear power plant. A series of experiments have been performed to investigate the behavior of hydrogen distribution during a severe accident in the nuclear power plant. This study aims to investigate the gaseous stratification break-up phenomenon in the presence of a vertical gas fountain on the top of containment using Code\_Saturne. The computed results not only analyze the flow field of the gas and the density distribution in the containment, but also measure the adequacy of Code\_Saturne concerning the multi-component gas flow calculation.

### Edited by:

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

### Reviewed by:

*Mohammad Alrwashdeh, Khalifa University, United Arab Emirates Weiqian Zhuo, Virginia Tech, United States*

> \*Correspondence: *Chunhui Zhang sarszch@163.com*

#### Specialty section:

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

Received: *20 May 2019* Accepted: *29 May 2019* Published: *19 June 2019*

#### Citation:

*Zhang C, Yu J, Zhong X and Saeed M (2019) A Numerical Investigation on Gaseous Stratification Break-Up Phenomenon of Air Fountain Experiment by Code\_Saturne. Front. Energy Res. 7:57. doi: 10.3389/fenrg.2019.00057* Keywords: numerical simulation, air fountain experiment, code\_saturne, hydrogen distribution, containment

## INTRODUCTION

A severe accident occurred as a result of huge earthquake at the Fukushima Daiichi Nuclear Power Plant, located in the East Pacific Ocean in Japan, on March 11, 2011 (IAEA, 2015). In this accident, the loss of backup power caused the core to lose continuous cooling of the coolant and finally led to the melting of the core. Then the hydrogen generated during the process of cladding oxidation. As the generated hydrogen reached to a certain concentration, an explosion occurred in the nuclear power plant. Since then, the research on hydrogen safety in the containment of nuclear power plants during a severe accident has become a hot topic for the nuclear safety experts.

Due to the existence of the buoyancy force, hydrogen will accumulate in the upper region of the containment (Rodi, 1982). When the local hydrogen and oxygen concentration reaches to a specific value, the explosion or combustion will occur, which results in worse consequences. It is worth noting that in this process, the mixed gas will be stratified due to the difference of density. As the hydrogen rises due to its lower density, the gas between the top of the containment and the source of hydrogen will be converted into hydrogen plus air. The mass fraction of hydrogen in the mixed gas increases with height, so that the density of the mixed gas decreases with height, and the stratification will take place in the mixed gas. The situation mentioned above is an ideal state and does not consider the influence and destruction of other phenomena. In the actual containment of nuclear power plant, there are several phenomena that may break the gas stratification during severe accident such as injection of a sprinkler system, hydrogen recombiner and gas injection, etc. Moreover, the steam will condensate on the wall, which leads to passive heat transfer and mass transfer. It will cause buoyancy-driven flow and exacerbate gas mixing. However, it does not

mean that the hydrogen stratification will not occur during a severe accident. Taking the LOCA accident of the pressurized water reactor as an example, for the pressurized water reactor which using vertical steam generator and higher containment, the main circuit is located in the lower position of the containment while the blasting valve of the pressure regulator relief box is located in the higher position. During the accident, the hydrogen and water vapor will be released if the main circuit break and a huge amount of substances and energy will be released in the lower portion of the containment. If the containment has internal and external parts, it will easily generate circulating flow, which exacerbates mixing. If the rupture is near the regulator release valve, the mass source and energy source will appear at a higher position in the containment. Then the gas stratification may occur. This phenomenon has already been observed in the PHDR(Projekt HDR-Sicherheitsprogramm) experiment (IAEA-TECDOC-1661, 2011).

In order to study these phenomena in depth, the researchers have performed a lot of experimental projects and also developed some code to tackle this problem. In 2006, Auban et al. performed an experiment with horizontal jet steam in the air using the PANDA(located at PSI in Switzerland) experimental facility (Auban et al., 2007). The ThAI (Thermal-hydraulics, Hydrogen, Aerosol, and Iodine) experimental facility was also used to perform horizontal steam jet in a layered helium mixture (Allelein et al., 2007). In 2008, Gupta et al. performed several tests focused on the possibility of the permanence of hydrogen stratification just beneath the containment dome under conditions determined by low velocity and a plume or jet of higher density (Gupta et al., 2008). In 2009, Mott and Woods investigated the mixing of a stratified fluid of finite volume by a turbulent buoyant plume (Mott and Woods, 2009). In 2007, Guyez et al. made an investigation on turbulent mixing at a stable density interface (Guyez et al., 2007).

### MODEL OF CODE\_SATURNE

Code\_Saturne is a general purpose computational fluid dynamics (CFD) free computer software package developed by EDF(Electricite De France). It solves the Navier-Stokes equations for 2D, 2D-axisymmetric and 3D flows, steady or unsteady, laminar or turbulent, incompressible or weakly dilatable, isothermal or not, with scalars transport if required (Archambeau et al., 2004).

The mass conservation equation is given by<sup>1</sup> :

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \,\overrightarrow{U}\right) = \Gamma \tag{1}$$

where ρ refers to the density of the fluid, U refers to the velocity of the fluid and Γ refers to the mass source term. This term is included in the mesh located near the air injection nozzle in this work.

The momentum conservation equation is given by

$$\frac{\partial}{\partial t} \left( \rho \overrightarrow{U} \right) + \nabla \cdot \left( \rho \overrightarrow{U} \otimes \bigotimes \overrightarrow{U} \right) = \nabla \cdot \left( \overset{\circ}{\sigma} \right) + \overrightarrow{ST} - \overset{\circ}{K} \overset{\circ}{U} \text{ (2)}$$

where −→ST – ⇀⇀ K −→<sup>U</sup> refers to additional momentum source term. The mesh located near the air injection nozzle products momentum source due to mass source in this work. σ refers to the Reynolds stress tensor. Here we take the k-ε model as an example, The expression of σ is:

$$\stackrel{\cong}{\sigma} = \stackrel{\cong}{\tau} - \left(P + \frac{2}{3}\rho k\right) \stackrel{\cong}{Id} \tag{3}$$

### EXPERIMENT INTRODUCTION AND MODEL BUILDING

### Experiment Introduction

In 2008, CEA(French Atomic Energy Commission) and LEEF(Experimental Fluid Mechanics Laboratory) performed an experiment of vertical fountain in the erosion of gaseous stratification on small scale. Helium is used to simulate hydrogen released by cladding metal water reaction in the experiment for safety reasons. The purpose of this experiment was to study the effect of air fountain in the low-density gaseous stratification in the top part of the containment (INSAG, 1996; Deri et al., 2010).

The facility of this experiment is a transparent box whose height is 1.29 m and the bottom is a square of 0.92 m side length as shown in **Figure 1**. There is a gap of 0.06 m height on the wall near the bottom of the box which is the outlet to the environment. The environmental condition is 1 atm and room temperature. There is a vertical air inlet of 0.02 m diameter at the center of the bottom. Six thermoelectric sensors were installed at different heights on one wall to measure the concentration of helium at these positions. The six

<sup>1</sup>Theoretical manual of Code\_Saturne. Available online at: http://code-saturne.org

positions are 0.49 m, 0.69 m, 0.83 m, 0.96 m, 1.09 m, and 1.29 m high, respectively.

At the beginning of the experiment, the box was filled with air. Then a 9.1 g helium was injected horizontally at 0.99 m high during 300 s and the overflowed air was exhausted through the outlet. The gaseous stratification was formed and the gas free diffused for 60 s after the helium injection. The air was injected vertically from the inlet in the center of the bottom during 360 s to 660 s. The experiment tested six different velocities of air injection to observe the influence on the erosion of the gaseous stratification. The gas free diffused after 660 s and the experiment continued till 2000 s. The six different velocities tested in the experiment are 0.411 m/s, 0.849 m/s, 1.286 m/s, 1.929 m/s, 2.803 m/s, 5.143 m/s, respectively.

According to the experimental results, the erosion effect caused by different air inject velocities are divided into three regimes. When the velocity is small, the molecular diffusion is dominant in the mix of gas. As the velocity increases, the buoyancy become more dominant in the mix of gas and the stratification zone was push upward. If the velocity continuously increases, the injected air will reach to the top and will break the stratification zone directly. Under this circumstance, the momentum is dominant in the mix of gas.

In order to distinguish these three regimes, a dimensionless number called Froude number was introduced as following:

$$Fr = \frac{U\_{ref}}{NL} \tag{4}$$

where L refers to local inject diameter; Urefrefers to the velocity of air injection in the stratification zone; N is the characterized parameter of stratification zone defined as:

$$N = \left(\frac{2g}{\rho\_1 + \rho\_2} \frac{\rho\_1 - \rho\_2}{H}\right)^{0.5} \approx 2s^{-1} \tag{5}$$

where ρ1and ρ<sup>2</sup> are, respectively, the density of the bottom and the upper zone.

According to the experimental data, six Froude numbers(0.16, 0.33, 0.50, 0.75, 1.09, 2.00) are chosen to observe the three regimes mentioned above. When Froude number = 2, an immediate mixing occurs. On the other hand, when Froude number less than 1, the buoyancy overwhelms the momentum and the air injection is thus unable to penetrate the stratification. So the Froude number 1.09 is known as the demarcation of the buoyancy-dominant regime and the momentum-dominant regime (Deri et al., 2010). In order to verify the ability to simulate the three regimes(free diffusion, buoyancy and momentum) by Code\_Saturne, we choose Froude number = 1.09 for the final simulation.

### Geometry Model

The whole geometry was simplified as a three-dimensional quarter size to reduce the simulation time and simplify the computation. Symmetry condition was imposed for the whole domain and simulated the same problem in the remaining experimental domain. The diameter of air inlet is 0.02 m, so we

use a square inlet with 0.01 m on each side which is nearly the size of one mesh cell as shown in **Figure 2**. We add a mass source in the mesh above the inlet instead of an inlet boundary condition for the concision of the boundary condition.

### User-Defined Model

Because the gas in the experiment facility remains constant (room temperature), the energy conservation equation is not considered in this case. This case contains two different kinds of gas (air and helium), therefore, a user-defined diffusion equation for helium mass fraction calculation is added to the total mass conservation equation of Code\_Saturne:

$$\frac{\partial \rho X\_{He}}{\partial t} + \nabla \cdot \left(\rho \,\overrightarrow{U} \, X\_{He}\right) - D\_{He} \nabla^2 \left(\rho X\_{He}\right) = 0 \tag{6}$$

where DHerefers to the diffusion coefficient of helium in air, which is achieved by experiment (7.35<sup>∗</sup> 10−5m<sup>2</sup> s −1 ) (Deri et al., 2010).

### RESULTS AND DISCUSSION

Before the calculation for the final results, some pre-work need to be done. First, the calculation of the initial density distribution was made to verify the accuracy of this simulation. The simulations starts at 300 s (after the helium injection) with no air injection till the end. Then the mesh sensitivity analysis was performed to ensure that the final results are mesh-independent solutions. At last, the final calculation included three stages: the free diffusion stage(300 s−360 s), the air injection stage(360 s−660 s), the free mixing stage(660 s−2000 s).

### Initial Condition

The experiment focuses on the break regime caused by the air fountain after the gas stratified, so we ignore the process of helium injection. Therefore, we choose the moment which stopped the injection of helium at 300 s as the initial time for our simulation. At this phase of the simulation the gaseous stratification formed and became stable. The total mass of helium will no longer increase. This initial time is also recommended by the experiment researchers and the article provides the detailed initial density distribution in the facility at 300 s.

In order to prove the accuracy of the initial density distribution, we use Code\_Saturne to simulate a case without air inlet (i.e., there is only molecular diffusion) and compare with the experimental data before the case of Fr = 1.09.

For the purpose of reducing the influence introduced by the difference of air density, the recorded experimental data adopted a dimensionless form:

$$
\rho' = \frac{\rho\_{\rm ref} - \rho}{\rho\_{\rm ref}} \tag{7}
$$

where ρref refers to the reference density i.e., the density of pure air in the experimental condition.

For the convenience of comparison, the results calculated by Code\_Saturne are also transformed by the formula above. The comparison is shown in **Figure 3**.

The calculation starts at 300 s of the physical time. It has been observed that the results of the three positions near the top are in better agreement with the experimental data well. At the beginning, the density of the mixture gas is very small which means the mass fraction of helium is large, so the diffusion process is very fast. But as the diffusion continues, the gradient of the mass fraction becomes smaller and finally stabilized. The processes of the three positions near the bottom are opposite to the above positions and the density decrease. The tendencies of the three curves are the same with the experimental data but the deviation is larger than the curves above. The calculation results are under estimated compared to the experimental data in late stage which mean that the total helium mass is a little larger than the experimental data. That's probably because there is no gas leak in numerical calculation.

### Calculation With Air Injection

### Test on Sensitivity of Mesh

After the verification of initial density distribution, some analysis tests on sensitivity of mesh has been performed. If there is no air injection, the velocity of the gas in the experimental facility is very low. Therefore, the size of the mesh has no effect on the calculation results. When the air injection velocity was set to 2.803 m/s, the size of the mesh will have a significant effect on the results.

This paper tested three kinds of mesh:


The three different curves represent the calculation results of three different fineness of mesh as shown in **Figure 4**. For the sake of simplicity, only the time-varying graph of density at Z = 0.2 m is shown. Fixed time steps i.e., 0.01 s is used in all three cases (which is also the best value obtained by time step sensitivity analysis). From **Figure 4**, we can see that there is a significant requirement for the fineness of the mesh if the time step is 0.01 s. If the mesh size is too coarse, it will create instability, resulting in concussions in the final results. And the rougher the mesh

size, the greater the oscillation amplitude. According to the test results, we have chosen the most refined mesh for the following calculation in this paper.

This test of sensitivity for the mesh was also applied to the initial condition calculation(without air injection). As the velocity of the gas in this case is very low, there is no difference for three kinds of mesh size, consequently, the coarse mesh is fine enough to have a mesh-independent solution.

### Calculation of the Case Using Fine Mesh

After the verification of initial density distribution and tests on sensitivity of mesh and time step, finally we got results as following figure:

In the experiment, the stage from 300 s to 360 s was free diffusion after helium injection. So we also considered this stage in our simulation. We can see the curves were smooth at the beginning according to **Figure 5** and **Supplementary Table 1**. Then the air injected vertically upward from the center of the bottom at speed of 2.803 m/s. The three curves in the upper region decline fast which means that the densities of three positions near the top increase quickly. This is not because the air went into the zone directly. As we see in **Figures 6**, **7**, at the first 2 s of air injection(360 s−362 s), the high-speed air directly flowed into the top part of the box. But after that, it will split into two eddies separately above and under the plane of 0.85 m height. The upper eddy entrained the gas of large density from the lower part to the top thus resulting in the increase of density in the upper zone. The lower eddy also entrained the gas of small density from the upper part so the density in the lower zone decreased. The three curves in the lower zone increase slowly during the stage of air injection(360 s−660 s). After the injection, the three lower curves in the stage of free mixing were similar to the case without air injection.

According to **Figure 5**, at the stage from 300 s to 660 s when the air injection stopped, the calculated results are in better agreement compared with the experimental data well. But at the next stage of free diffusion, the calculated results of density at six positions are under estimated as compared with the experimental data. We see that this deviation is close to the case without air injection. So we can inference that it is the deviation caused by the free diffusion that leads to the same problem. As it mentioned in section Initial Condition, that's probably because there is no gas leak in numerical calculation.

### CONCLUSIONS

This paper presents a numerical investigation on gaseous stratification break-up phenomenon of air fountain experiment by Code\_Saturne. The numerical investigation includes the verification of initial condition, the sensitivity test of time step, the sensitivity test of mesh size and the final numerical calculation of the air fountain experiment.

The numerical simulation for the initial condition shows the result of gaseous stratification phenomenon. This result verifies the model of buoyancy driven flow in Code\_Saturne and establishes a good foundation for subsequent calculation. After the sensitivity test of mesh size, the best parameters and mesh for the numerical simulation of the air fountain experiment were found.

The final calculation results show the process of gaseous stratification break-up phenomenon of air fountain. At the beginning, the air injected to the zone of gaseous stratification directly. After the injection the flow was driven by buoyancy force and molecular diffusion. The largest relative error of the density in all six positions is <3% which verifies the accuracy of the Code\_Saturne.

### DATA AVAILABILITY

All datasets generated for this study are included in the manuscript/**Supplementary Files**.

### REFERENCES


### AUTHOR CONTRIBUTIONS

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

### ACKNOWLEDGMENTS

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

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg. 2019.00057/full#supplementary-material

Supplementary Image 1 | Fine mesh.

Supplementary Table 1 | Data of calculation result.

IAEA (2015). The Fukushima Daiichi Accident. Vienna: IAEA.


Rodi, W. (1982). Turbulent Buoyant Jets and Plumes. Pergamon Press.

**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 © 2019 Zhang, Yu, Zhong and Saeed. 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.

# Numerical Analysis of Turbulent Flow and Heat Transfer in Internally Finned Tubes

Zhanwei Liu<sup>1</sup> , Yanwei Yue<sup>2</sup> , Luchao She<sup>1</sup> and Guangming Fan<sup>1</sup> \*

*<sup>1</sup> Fundamental Science on Nuclear and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China, <sup>2</sup> Research and Development Center, China Nuclear Power Engineering Co., Ltd, Beijing, China*

In this study, single-phase heat transfer enhancement in internally finned tubes is investigated numerically. The influence of fin number, helix angle, fin height, fin width, and shape on the flow and heat transfer characteristics is studied. The research results indicate that the resistance coefficient and Nusselt number both increase with the increment of these parameters, among which the helix angle has the largest impact on the heat transfer enhancement. In addition, the shape of fins also has a small effect on the flow and heat transfer, and the heat transfer effect of triangular fins is the best.

Keywords: internally finned tube, heat transfer enhancement, numerical simulation, single -phase flow, heat transfer coefficient

### Edited by:

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

#### Reviewed by:

*Han Bao, Idaho National Laboratory (DOE), United States Zhitong Bai, University of Michigan, United States*

> \*Correspondence: *Guangming Fan fanguangming007@hotmail.com*

#### Specialty section:

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

> Received: *30 April 2019* Accepted: *02 July 2019* Published: *18 July 2019*

#### Citation:

*Liu Z, Yue Y, She L and Fan G (2019) Numerical Analysis of Turbulent Flow and Heat Transfer in Internally Finned Tubes. Front. Energy Res. 7:64. doi: 10.3389/fenrg.2019.00064*

### INTRODUCTION

Single-phase convection heat transfer enhancement techniques are widely applied in industries such as petroleum, chemical engineering, etc. In recent years, its applications in the field of nuclear technology has been an engaging research area of heat transfer augmentation. In the advanced nuclear power plant like AP1000, HPR1000 and CAP1400, passive residual heat removal system (PRHR) (Schulz, 2006; Xing et al., 2016; Zheng et al., 2016) can discharge residual heat through single-phase convection heat transfer under normal condition, therefore enhancement of its heat transfer efficiency is of great significance to the safe operation of nuclear reactors.

In order to strengthen the single-phase convection heat transfer technology, many heat exchanger elements have been developed, the internally finned tube is one of them. Ji et al. (2015) found that their thermal-hydraulic performance is the best among all enhanced tubes through a comprehensive literature review. Much work has been done to investigate the heat transfer characteristics of internally finned tubes in condensation (Seo and Kim, 2000; Kim et al., 2002; Yu et al., 2002; Yun et al., 2002) and evaporation (Chamra and Webb, 1996; Chamra et al., 2004, 2005; Chamra and Mago, 2006), but the heat transfer characteristics of single-phase in these tubes are not thoroughly understood for lack of research. The internally finned tube has high requirements for water quality, and the working fluid is required to be clean enough with almost no impurities, otherwise they may affect the heat transfer performance. Compared with other heat exchanger elements, the tremendous advantage of internally finned tubes is that its heat transfer effect is greater than the increase of pressure drop. Generally, this effect can be seen in turbulent state, but neither the enhancement of heat transfer nor that of pressure drop is obvious in laminar state, thus they are recommended to be used in turbulent flow (Al-Fahed et al., 1999; Copetti et al., 2004).

Previous experimental studies have been carried out on single-phase flow of different working media to investigate the influence of different geometrical sizes (fin height, fin width, helix angle, diameter, fin number) on the heat transfer characteristics. Wang et al. (1996) experimentally tested seven internally finned tubes with different geometric sizes such as diameter, fin height, helix angle and fin numbers, and established heat transfer and pressure drop correlations in a wide range of Reynold number. Han and Lee (2005) studied the single-phase heat transfer and flow characteristics of water in four internally finned tubes with different diameters, fin heights, helix angles and fin pitchs experimentally, and developed correlations for friction factor and heat transfer coefficients. Jensen and Vlakancic (1999) carried out experimental studies on turbulent flow in eight internally finned tubes with different geometrical parameters of helix angles, fin heights, fin numbers and fin widths, they obtained the correlations between friction coefficient and Nusselt number. Later Zdaniuk et al. (2008) determined and compared friction factors as well as heat transfer coefficients of eight internally finned tubes experimentally. Empirical correlations for one particular internally finned tubes

was developed by Siddique and Alhazmy (2008) through experimental research. And a critical Reynolds number Recr was discovered by Li et al. (2007) through experimental study. Ji et al. (2011) measured the turbulent heat transfer and friction factors of sixteen internally finned tubes with various geometrical parameters including fin number, fin height, diameter, helix angle, and developed an equation for the heat transfer in the fully developed flow as an extension of the Gnielinski equation. Friction formula was also developed by Celen et al. (2013). Besides, Experimental studies of R134, R12 (Eiamsa-ard and Wongcharee, 2012) and the water-CuO nanofluid (Eiamsa-ard and Wongcharee, 2013) were also conducted.

Due to the limitations of manufacturing and experimental conditions, only a few tubes with specific geometrical parameters have been studied in previous experiment research, so it is difficult to obtain the detailed and comprehensive heat transfer characteristics in different tubes.

There are only a few numerical studies on the heat transfer characteristics of internally finned tubes. Agra et al. (2011) studied the heat transfer characteristics and pressure drop of five internally finned tubes in the turbulent state numerically, they

compared their simulation results with the experimental data of Zdaniuk et al. (2008), and discovered that CFD predicted the experimental data more accurately than the Blasius equation. Their results also indicated that the internally finned tube has a higher heat transfer coefficient than the corrugated tube. Celen et al. (2014) studied the pressure drop in one internally finned tube with TiO2-water nanofluids numerically, and described the average and local values of temperature, pressure, and velocity distributions of the tube. Dastmalchi et al. (2017) investigated the characteristics of laminar oil flow in internally finned tubes with different fin heights and helical angles numerically. They

angles. (A) Nu. (B) f.

TABLE 1 | Increment of heat exchange area under different fin number and helix angle.


concluded that there is an optimal height and an optimal helix angle which leads to the maximum performance of heat transfer.

The foregoing literatures indicate that a comprehensive study on the flow and heat transfer characteristics considering a wide range of parameters including fin number, helix angle, fin height, fin width and shape of fins remains to be conducted. Influence of geometric parameters on characteristics of heat transfer could be implemented with modern visualization techniques and computational fluid dynamics (CFD) tools. Therefore, the purpose of the study is to conduct a comprehensive numerical study on the factors affecting the heat transfer efficiency of the internally finned tubes.

### NUMERICAL METHOD AND PROCEDURE

### Geometric Model

The fins are uniformly distributed in circumferential direction and along the axis with the same helix angle, so they are periodic in geometric structure. The rectangular fins studied have the diameter D = 20 mm, length L = 20 mm, fin number N = 30, the fin height e = 0.5 mm, the width s = 0.5 mm, the helix angle of 30◦ . The cross section is shown in **Figure 1**.

Using the periodicity of the internally finned tube and taking one pitch of the fully developed segment as the computational domain, the difficulty of grid generation is greatly reduced. The established model is shown in **Figure 1**.

### Grid Generation and Boundary Conditions

In order to determine a proper grid for the numerical simulations, a grid independence study is carried out for the smooth and internally finned tubes. A section of circular tube with a diameter of 20 mm and a length of 20 mm is taken as the computing domain, and meshing is used to generate grids.

The mesh number chosen in this paper is the minimum of the mesh number that does not affect the simulation results. When the number of grids is lower than a certain value, the change of the number of grids will have an impact on the simulation results. The larger the number of grids, the closer to the stable result. When the number of grids is higher than a certain value, the encrypted grid will have no impact on the simulation results.

As shown in **Figures 2A,B**, the grid of smooth tube is made by independence verification. Grid quality measurement method is skewness, the value <0.7 is acceptable, so the grid generation is credible as most values are below 0.6.

Realizable k-epsilon model satisfies the constraint conditions of Reynolds stress, so it can be consistent with the real turbulence at Reynolds stress. This is something neither the standard kepsilon model nor the RNG k-epsilon model can do. The advantage of this feature in the calculation is that it can more accurately simulate the diffusion velocity of plane and circular jets. Meanwhile, in the calculation of rotating flow, boundary layer with directional pressure gradient and separation flow, the calculation results are more in line with the real situation.

TABLE 2 | Increment of heat exchange area under different fin height and helix angle.


The Realizable k-epsilon model is a newly emerging k-epsilon model, and although its performance has not been proven to be superior to the RNG k-epsilon model, studies on separated flow calculations and complex flow calculations with secondary flows show that the Realizable k-epsilon model is the most excellent turbulence model among all the models. Given the advantages of the Realizable k-epsilon model, this paper chooses this model for calculation.

We used the standard k − ε model in fluent to solve the turbulent kinetic energy equation and turbulent dissipation rate

equation, then calculated the turbulent viscosity with the value of the sum, and finally obtained the solution.

The enhanced wall treatment is selected to obtain the detailed flow condition near the tube wall. In the two-layer zone model, the near-wall flow can only be divided into two regions, namely the region affected by viscosity and complete turbulence. The two regions are distinguished by Re<sup>y</sup> based on the distance y to the wall.

$$\mathrm{Re}\_{\mathcal{V}} = \frac{\rho \sqrt{k} \mathcal{Y}}{\mu} \tag{1}$$

The optimal grid division requires that the first grid is at the position y+ =1. It is acceptable as long as it is in the bottom viscous layer even a little bit bigger, like y+ =4∼5.

After the grid is read into FLUENT, the parameters and boundary conditions are set as follows: single phase water with constant physical properties, the constant heat flow boundary condition and the given heat flux is q = 1,00,000 W/m<sup>2</sup> , the countercurrent volume temperature is set as 300 K, the coupling solution of pressure and velocity is SIMPLEC algorithm, and the momentum and energy equations are solved by the second-order upwind scheme.

### Verification of the Numerical Procedure

Based on the pressure drop value per unit length obtained, the average f can be calculated using Darcy formula:

$$f = \frac{2\,\Delta PD}{\rho u^2} \tag{2}$$

Based on the difference between the average wall temperature and the average fluid temperature, the average Nusselt number can be calculated:

$$Nu = \frac{qD}{\Delta T \lambda} \tag{3}$$

For forced convection heat transfer in smooth circular tubes, the resistance coefficient can be calculated using the Filonenko formula (Filonenko, 1954):

$$f\_{\mathcal{F}} = (1.82 \lg Re - 1.64)^{-2} \tag{4}$$

The nussel number can be calculated according to the Gnielinski (Gnielinski, 1975) formula:

$$Nu = 0.012 Pr^{0.4} (Re^{0.87} - 280) \tag{5}$$

The selected grid scheme is used for verification. The comparison results with the classical experimental formula are shown in **Figures 3A,B**. The calculated values of f and Nu are compared with the values of Filonenko formula (Filonenko, 1954) and Gnielinski formula (Gnielinski, 1975), respectively. It can be seen that the calculated value of the resistance coefficient and the Nusselt number are in good agreement with the theoretical value. The maximum error between the calculated value and the theoretical value of the resistance coefficient is 6.25%, most of which are within 3%.The maximum error between the calculated value and the theoretical value of nusserl number is 11.2%, most of which is about 5%. The maximum error occurs in the lower Reynolds number region. The results show that the application of periodic boundary conditions and the established calculation model are correct and can accurately simulate the turbulent flow and heat transfer in the tube.

RESULTS AND DISCUSSION

In this section, the effect of fin number, fin height, fin width, helix angle and the shape of fins on the heat transfer characteristics of internally finned tube are studied.

**Figure 4** shows the effect of the fin number on the Nu and f at different helix angles. The fin height and fin width are fixed at 0.5 mm, and the fin number are 15, 30, and 45, respectively. It can be seen that the Nu and f both increase with the augment of the fin number and helix angle.

**Table 1** shows the variation of the increment of heat exchange area with fin number and the helix angle. Here the increment of heat exchange area is defined as:

$$\delta = \frac{A\_{fin} - A\_{smooth}}{A\_{smooth}} \tag{6}$$

$$A\_{smooth} = 2\pi DL\tag{7}$$

$$A\_{fin} = \frac{(2\pi D + 2Ne) \, L}{\cos \chi} \tag{8}$$

It can be seen that the rate at which Nu increases is bigger than the increase speed of heat exchange area and f, which indicates that the enhancement effect of heat transfer is greater than the increase of resistance and heat transfer area.

**Figure 5** shows the velocity cloud diagram and temperature cloud diagram of internally finned tubes with fin number of 15, 30, and 45. It can be seen that at the sharp corners of the windward side, the wall is the most strongly impacted by the fluid, the boundary layer is the thinnest, and the velocity gradient is the largest, so a great shear stress is generated which means that there is a great loss of energy. Meanwhile, the corresponding fluid temperature is the lowest, and the heat transfer convection is the strongest. As the fin number increases, the number of sharp corners increase, resulting in an increase in both f and Nu. Along with the increase of fin number, the intercostal region decreases accordingly, the fluid viscosity of the intercostal channel increases, which reduces the fluid velocity of fin root and surface, thickens boundary layer, and inhibits the heat transfer.

**76**

(A) Nu. (B) f.

### Effect of Fin Height

**Figure 6** shows the changes of Nu and f with the fin height under different helix angles. The fin number is 30, the fixed fin width is 0.5 mm, and the fin height is 0.3, 0.5, 0.7, and 0.9 mm, respectively (the dimensionless fin height H changes from 0.03 mm to 0.09 mm). It can be seen that the Nu and f both increase with the augment of the fin height, as it not only increases the heat transfer area on both fin sides, but also promotes the interactions between the fin top. When the helix angle is >20◦ , Nu and f grow significantly compared with small helix angle such as 10 and 20◦ with the increase of fin height, which indicates that the helix angle plays an important role in the heat transfer.

**Table 2** shows the changes of heat exchange area with the fin height and helix angle. Nu and f both increase with fin height, because the higher the fin height, the greater the disturbance to the fluid near the wall surface. Therefore, the heat transfer coefficient is enhanced and the resistance coefficient is increased. Increasing the spiral angle will further enhance the heat transfer and the resistance.

The increment of the fin height increases the interaction between fins and fluid, resulting in rapid increase in f and Nu.

### Effect of Fin Width

Generally, the fin width has a complex effect on the flow and heat transfer in internally finned tubes. The increment of fin width does not increase the heat transfer area, but it will cause the change of the size of the intercostal passage. Due to the limitations of time and conditions, only the effect of fin width within a small range between 0.5 and 0.9 mm is studied under the condition that the number of fins is 30.

**Figure 7** shows the changes of Nu and f with the fin width. As we can see, with the increment of the fin width, Nu and f both increase slowly.

**Figure 8** shows the velocity cloud and temperature cloud when the fin width is 0.5, 0.7, and 0.9 mm, respectively. It can be seen that the increase of fin width means that the area of the fin top increases, equivalently the area in contact with the fluid is larger, and fins are more impacted by the fluid, so heat transfer and friction resistance are both enhanced.

Similar to the increase of fin number, as fin width increases, the intercostal region decreases accordingly, the viscous effect of the fin side and the inner wall against the intercostal fluid also increases, which reduces the flow rate in the intercostal region, thickens the boundary layer and suppresses the heat transfer between the fin root and the wall, thus causes the decrease of Nu and f. The opposite effect of the two factors makes Nu and the f changed slowly.

### Effect of Helix Angle

**Figure 9** shows the changes of f and Nu with the helix angle under different fin numbers to study the effect of helix angle. The fixed fin width is 0.5 mm, and the fixed fin height is 0.5 mm. The reason that we chose the spiral angle of 10 to 40◦ is that the spiral Angle is usually 10 to 45◦ in industrial applications, so we chose the most representative parameters for numerical simulation.

As we can see, Nu and f both increase with helix angle dramatically. The bigger the helix angle is, not only the area of the fin exposed to the incoming flow windward increases, but also the change of fluid flow direction increases, that is, the greater the obstruction to the incoming flow. Therefore, the incoming flow has a stronger impact on the windward side of the fin, and f and Nu both increase.

**Figure 10** shows the velocity vector of intercostal area. The arrows in the bottom parts of **Figures 10A–C** point to left bottom, while the arrows in the bottom part of **Figure 10D** point to left, this is because the larger the spiral angle is, the more the fluid near the wall is affected by the spiral angle, and the greater the included Angle with the mainstream is. So it tends to the horizontal.

TABLE 3 | Comparison of parameters of three kinds of fin tubes.


### Effect of Fin Shape

**Figure 11** shows the changes of f and Nu with three fin shapes, as shown in the picture, the average f in rectangular fins is the highest and that of triangular fins is the lowest, and that of circular fins is between them. The maximum difference between rectangular and triangular fins is 9.85%, and with the increase of Re, the difference between f gradually decreases. As can be seen from **Figure 11**, the Nu of rectangular and triangular fins is basically the same, and the triangular fin is slightly larger. With the increase of Re, the difference of Nu between the circular fin and the other two kinds of fins gradually increases, and the maximum difference is 5.9%.

**Table 3** compares the increment in heat transfer area, Nu and f caused by the three types of fins in turbulent state. It can be seen from the table that the heat transfer area of rectangular fins increases the most, while that of triangular fins increases the least, and the difference of the three types of fins lies in the shape of fin top. For rectangular fins and triangular fins, the top has a sharp angle, while the top of circular fins is a smooth surface, which lead to different flow and heat transfer characteristics.

**Figure 12** shows the velocity field cloud and temperature field cloud of rectangular fin, triangular fin and circular fin respectively when the Re is 70000 and ρ is 100000. As we can see, the sharp angles of rectangular fins and triangular fins are impacted by the incoming flow, which greatly enhance the heat transfer in the region. Therefore, the average Nu is higher than that of circular fins, and as the Re increases, the interaction between the fluid and the sharp angle of rectangular fins and triangular fins becomes more intense. It is noted that the Nu of the triangular fins is slightly higher than that of the rectangular fins, because compared with rectangular fins, the intercostal passage between adjacent fins of the tube in the triangular fins is larger, and the larger intercostal area could pump more turbulence into the root area of the fin, thinning the boundary layer, thus enhancing heat transfer. The heat transfer area of the triangular fins increases the least, f increases the least and Nu increases the most, so the heat transfer performance is the best.

### CONCLUSION

In this study, FLUENT software and periodic boundary conditions are used to simulate the flow and heat transfer in internally finned tubes comprehensively, and to study the influence of various geometric parameters. The parameters are fin number (15–45), fin height (0.3– 0.9 mm), fin width (0.5–0.9 mm), helix angle (10–40◦ ) and the operating condition is Re = 44430. Meanwhile, the influence of the fin shapes on the heat transfer and the influencing mechanism are also explored. The conclusions are as follows:


### REFERENCES


### DATA AVAILABILITY

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

### AUTHOR CONTRIBUTIONS

In the course of the completion of this thesis, all authors have substantial contributions to design of the work, and the contribution of YY is equivalent to that of the first author. GF provided the research direction of the subject, all authors worked together to define the methodology and procedures. LS built the geometric model. ZL carried out numerical simulation and obtained, collated, and analyzed the data. GF and LS also provided guidance and assistance in the process. ZL wrote the first draft of the paper and corrected it by GF and LS. YY provided important guidance in the process, and his contribution is equivalent to that of the first author. All authors approved the final version to be published. All authors agreed to be responsible for all aspects of the work.

### FUNDING

This work was supported by National Natural Science Foundation of China (11875117).

### ACKNOWLEDGMENTS

ZL would like to extend his sincere gratitude to his supervisor GF for his instructive advice and guidance on his thesis. ZL would also like to thank YY for his important guidance in this process, whose contribution is equivalent to that of the first author.


**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 © 2019 Liu, Yue, 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(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.

## NOMENCLATURE


·K)-1


### Greek symbols


# Diameter Effect on the Wall Temperature Behaviors During Supercritical Water Heat Transfer Deterioration in Circular Tubes and Annular Channels

#### Hui Cheng<sup>1</sup> , Aaron Bai Yanlin<sup>2</sup> , Jun Wang<sup>3</sup> and Jiyun Zhao<sup>1</sup> \*

<sup>1</sup> Department of Mechanical Engineering, City University of Hong Kong, Hong Kong, China, <sup>2</sup> Shanghai Chuangyi Fluid Machinery Co., Ltd., Shanghai, China, <sup>3</sup> College of Engineering, University of Wisconsin-Madison, Madison, WI, United States

Diameter effect on heat transfer deterioration (HTD) for supercritical water upward flow in circular tubes and annular channels at low mass flux and high heat flux was studied numerically based on validated turbulence model. When the same boundary conditions were applied, i.e., mass flux, heat flux, and inlet temperature, it was observed that for circular tubes the first peak of wall temperature moves upstream and its magnitude increases at first then decreases as tube diameter increases, the second peak moves downstream as tube diameter increases. For annular channels with a fixed inner diameter, HTD is suppressed with small outer diameter, while HTD occurs gradually as outer diameter increases. The phenomenon agrees with previous experiment studies. The mechanism was analyzed based on fully developed turbulent velocity profiles at the heated sections inlet. Increasing inner diameter for circular tubes or outer diameter for annular channels will result in the decrease of maximum velocity and velocity gradient in the near heated wall region, which makes velocity profile in this region much easier to be flattened by the buoyancy. Then an M-shaped axial velocity profile is formed, which will significantly decrease the Reynolds shear stress and the turbulent kinetic energy and hence impair the heat transfer and cause HTD. At the same flow conditions, HTD is much easier to occur in circular tubes than in annular channels with the same hydraulic diameters.

Keywords: supercritical water, heat transfer deterioration, diameter effect, maximum velocity, velocity gradient, turbulent kinetic energy

### INTRODUCTION

The supercritical water reactor (SCWR) is one of the six generation IV nuclear systems. One advantage of SCWR is to increase the nuclear power plants thermal efficiency, which is now about 36–38%, to about 45–50% (Pioro and Mokry, 2011). Because no phase change occurs for supercritical water, SCWR can reduce the capital and operational costs through eliminating steam related facilities (Cheng et al., 2003; Pioro and Duffey, 2005).

Pioro and Mokry (2011) provided the pressure-temperature diagram of water (**Figure 1**). The critical temperature and critical pressure of water are 373.95 and 22.06 MPa, respectively. **Figure 2**

### Edited by:

Bruno Panella, Polytechnic University of Turin, Italy

#### Reviewed by:

Walter Ambrosini, University of Pisa, Italy Ivo Kljenak, JoŽef Stefan Institute (IJS), Slovenia

> \*Correspondence: Jiyun Zhao jiyuzhao@cityu.edu.hk

#### Specialty section:

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

Received: 27 May 2019 Accepted: 19 July 2019 Published: 02 August 2019

#### Citation:

Cheng H, Yanlin AB, Wang J and Zhao J (2019) Diameter Effect on the Wall Temperature Behaviors During Supercritical Water Heat Transfer Deterioration in Circular Tubes and Annular Channels. Front. Energy Res. 7:73. doi: 10.3389/fenrg.2019.00073 shows the thermophysical properties of water at 25.3 MPa. Because of the strong variation of the thermophysical properties in the area near the pseudocritical point, the characteristics of supercritical water heat transfer are quite different from

FIGURE 1 | Pressure-temperature diagram of water, reprinted from Pioro and Mokry (2011) with permission of IntechOpen.

subcritical water (Shitsman, 1963; Cheng et al., 2003, 2007; Hu and Gu, 2018; Xiao et al., 2018).

In the design of SCWR, one important safety issue due to the drastic property variation is the heat transfer deterioration (HTD). Shitsman was one of those who early pointed out the phenomenon (Shitsman, 1963). HTD usually occurs in a situation when fluid at supercritical pressure flows upward in a circular tube or an annular channel while it goes through the pseudocritical point at low mass fluxes and high heat fluxes conditions. When HTD occurs, the heated wall temperature increases drastically and then recovers, which may cause the fuel rod failure in SCWR. This phenomenon has been studied for decades and it is concluded that buoyancy force leads to the HTD (Jackson and Evans-Lutterodt, 1968; Hall and Jackson, 1969; Hall, 1971; Jackson et al., 1989; Kurganov and Kaptil'ny, 1992; Wang et al., 2004; Bae et al., 2005, 2008; Bae and Kim, 2009; Jackson, 2013; Cheng et al., 2017).

This paper aims to investigate the diameter effect on the behavior of HTD and to reveal the mechanisms behind these phenomena, which could provide some valuable guidance and support for the risk assessment in the design of SCWR. In the following, studies of diameter effect in heat transfer experiments of supercritical fluids in the past 50 years are reviewed (Yildiz and Groeneveld, 2014).

Yamagata et al. (1972) experimentally studied heat transfer characteristics of supercritical water flowing vertically upward in tubes. They reported the relation between the limit heat flux and mass flux for the flow in 10 mm ID tube (**Figure 3A**).

For a given mass flux, when heat flux is larger than the limit value, HTD occurs. However, this relation was achieved at a fixed tube diameter (10 mm ID) and may not be applicable to other diameter tubes. For a given mass flux, heat flux and inlet temperature, if we know the diameter effect on the heat transfer, then we can estimate whether HTD is more or less severe in a different diameter tube based on previous experimental data and get a deeper understanding on HTD phenomenon.

pressure with different heat fluxes, adapted from Zahlan et al. (2015) with permission of Elsevier.

Ackerman (1970) conducted experiments on supercritical water heat transfer at extensive flow conditions in tubes of several different diameters. For a given pressure and tube diameter, the limit curve between the heat flux which will cause HTD and mass flux is presented in **Figure 9** in their paper (Ackerman, 1970). The limit curve shows that operation in the lower right region of the curve for a given diameter will cause HTD. Operation above and to the left will result in normal heat transfer. For a given pressure and mass flux, increasing the diameter will move the curve to the left and HTD will occur at a lower heat flux. They found that the HTD onset heat flux decreased about 28.6% for the 24.4 mm tube when compared with the 9.4 mm tube.

Shiralkar and Griffith (1970) studied HTD of supercritical CO<sup>2</sup> flowing upward in two tubes with 3.18 and 6.35 mm diameter, respectively. It was observed that the peaks of wall temperature during HTD of the 6.35 mm ID tube were much sharper than the 3.18 mm ID tube.

Yamashita et al. (2003) conducted experiments on heat transfer of a supercritical pressure fluid flowing upward in uniformly heated tubes with different diameters, using HCFC22 as a test fluid. They reported the relation between the HTD onset heat flux and bulk enthalpy considering the tube diameter effect. It was found that the limit heat flux increases with smaller tube diameters, which is similar to Ackerman et al.'s results.

Song et al. (2008) experimentally investigated the heat transfer of supercritical CO<sup>2</sup> flowing upward in tubes at two different diameters (4.4 mm and 9.0 mm ID) for the same mass flux and heat flux. They found that in larger diameter tubes the flow was more susceptible to HTD due to buoyancy.

Kim et al. (2008) conducted lots of experiments on supercritical CO<sup>2</sup> heat transfer for various geometries including tubes with 4.4 and 9.0 mm ID, respectively, and an annular channel with 10 mm OD and 8.0 mm ID. They found that under the same flow conditions with low mass fluxes, HTD occurred in a 9.0 mm tube, while HTD was suppressed in the 4.4 mm tube and in the annular channel.

Bae et al. (2010) studied heat transfer to supercritical CO<sup>2</sup> flowing upward and downward in tubes with 6.32 mm ID. The results were compared with the data provided by Song et al. (2008) and Kim et al. (2008) of 4.4 mm and 9.0 mm tubes and they found that under the relatively low heat flux of 30 kW/m<sup>2</sup> , HTD occurred in the 9.0 mm tube while normal heat transfer persisted in the tubes of 4.4 and 6.32 mm.

Zahlan et al. (2015) performed extensive thermal hydraulic tests for CO<sup>2</sup> upward flow in two heated tubes with 8 mm and 22 mm ID at supercritical pressure conditions. The heat transfer characteristics of supercritical CO<sup>2</sup> in the two tubes are shown in **Figure 3B**. When q equals 50 kW/m<sup>2</sup> , normal heat transfer occurs in both tubes. While when q equals 70 kW/m<sup>2</sup> , some HTD is observed in the 22 mm ID tube but not in the 8 mm ID tube. When q equals 125 kW/m<sup>2</sup> , obvious HTD occurs in both tubes and it can be seen the wall temperature peak moves upstream and its value is larger in 22 mm ID tube.

It can be concluded from these previous experimental results that the HTD onset heat flux is decreased with larger diameter tubes, HTD is much easier to occur in larger diameter tubes for the same flow conditions, i.e., mass flux, heat flux, and inlet temperature. However, these conclusions are all drawn directly from experimental data and the mechanism behind these phenomena is still not clear. Compared with experiments, numerical method has got its advantage to extract more information from the fluid field and has already been successfully applied in the researches of supercritical fluid heat transfer (Koshizuka et al., 1995; Yang et al., 2007; Palko and Anglart, 2008; Cheng et al., 2017). In the present work, numerical simulations were conducted for supercritical water flowing vertically upward in circular tubes and annular channels at low mass flux and high heat flux based on validated turbulence model to study the diameter effect on HTD and the corresponding mechanisms.

### NUMERICAL METHODOLOGY

### Governing Equations and Turbulence Model

Because the geometries of circular tubes and annular channels are axisymmetric, the heat flux is uniformly applied on the heated wall, 2D axisymmetric model is used in the current simulation, which aims to greatly reduce the computational cost. The continuity, momentum and energy governing equations in cylindrical coordinates are given by:

Equation of continuity:

$$\frac{1}{r}\left\{\frac{\partial}{\partial \boldsymbol{x}}\left(r\rho U\right) + \frac{\partial}{\partial r}\left(r\rho V\right)\right\} = 0\tag{1}$$

Equation of U-momentum:

$$\begin{split} \frac{1}{r} \left\{ \frac{\partial}{\partial \mathbf{x}} \left( r\rho U^2 \right) + \frac{\partial}{\partial r} \left( r\rho V U \right) \right\} &= -\frac{\partial p}{\partial \mathbf{x}} + \rho \mathbf{g} + \frac{1}{r} \left\{ 2\frac{\partial}{\partial \mathbf{x}} \left[ r\mu\_{\varepsilon} \left( \frac{\partial U}{\partial \mathbf{x}} \right) \right] \right\} \\ &+ \frac{\partial}{\partial r} \left[ r\mu\_{\varepsilon} \left( \frac{\partial U}{\partial r} + \frac{\partial V}{\partial \mathbf{x}} \right) \right] \end{split} \tag{2}$$

Equation of V-momentum:

$$\begin{split} \frac{1}{r} \left\{ \frac{\partial}{\partial \mathbf{x}} \left( r\rho UV \right) + \frac{\partial}{\partial r} \left( r\rho V^2 \right) \right\} &= -\frac{\partial p}{\partial r} + \frac{1}{r} \left\{ \frac{\partial}{\partial \mathbf{x}} \left[ r\mu\_{\varepsilon} \left( \frac{\partial V}{\partial \mathbf{x}} + \frac{\partial U}{\partial r} \right) \right] \right. \\ &\left. + \ 2 \frac{\partial}{\partial r} \left[ r\mu\_{\varepsilon} \left( \frac{\partial V}{\partial r} \right) \right] \right\} - 2 \frac{\mu\_{\varepsilon} V}{r^2} \end{split} \tag{3}$$

where U is the x-direction velocity, V is the r-direction velocity, ρ is the density. µ<sup>e</sup> is the effective viscosity given by µ<sup>e</sup> = µ+µ<sup>t</sup> , µt is the turbulent viscosity defined by:

$$
\mu\_t = \rho \mathcal{C}\_{\mu} f\_{\mu} \frac{k^2}{\varepsilon} \tag{4}
$$

where fµ is the damping function considering near-wall effect, Cµ is a constant.

Equation of energy:

$$\begin{split} \frac{1}{r} \left\{ \frac{\partial}{\partial \boldsymbol{x}} \left( r \rho \boldsymbol{U} \boldsymbol{h} \right) + \frac{\partial}{\partial r} \left( r \rho \boldsymbol{V} \boldsymbol{h} \right) \right\} &= \frac{1}{r} \left\{ \frac{\partial}{\partial \boldsymbol{x}} \left[ r \left( \frac{\mu}{Pr} + \frac{\mu\_t}{\sigma\_t} \right) \frac{\partial \boldsymbol{h}}{\partial \boldsymbol{x}} \right] \\ &+ \frac{\partial}{\partial r} \left[ r \left( \frac{\mu}{Pr} + \frac{\mu\_t}{\sigma\_t} \right) \frac{\partial \boldsymbol{h}}{\partial r} \right] \right\} \tag{5} \end{split}$$

where Pr is molecular Prandtl number, h is enthalpy. σ<sup>t</sup> is turbulent Prandtl number, which is 0.9 in all simulations.

In the past two decades, various turbulence models based on Reynolds Averaged N-S equations were used to investigate supercritical fluid heat transfer by a lot of researchers. However, these models are not universally applicable and have their own limitations, especially in predicting heat transfer deteriorations. It is widely accepted that k-ε models are usually found to overestimate the heat transfer deterioration while k-ω models usually underestimate it (Pucciarelli et al., 2015). Hence, it should be very careful when choosing these models to quantitatively study supercritical fluid heat transfer.

The Shear Stress Transport (SST) k-ω turbulence model has the advantage to combine the accuracy of k-ω model in near wall region with the good performance of k-ε model in the main flow region. Several researchers (Palko and Anglart, 2008; Wen and Gu, 2010, 2011; Jaromin and Anglart, 2013) claimed that SST k-ω model has better performance in predicting HTD in their applications. To qualitatively study the diameter effect on HTD, the SST k-ω turbulence model is used here. Details of SST k-ω turbulence model is referred to the help documents of Fluent 15.0.

### Computational Domains

To investigate diameter effect on HTD, 7 circular tubes with different IDs ranging from 7 to 32 mm and 6 annular channels with fixed 8 mm inner diameter and different outer diameters ranging from 12 to 32 mm are used as computational domains. The geometric parameters for these circular tubes and annular channels are shown in **Table 1**. Schematic of the boundary conditions of circular tubes and annular channels are shown in **Figure 4**. Among all these cases, geometric parameters of Case C2 is taken from experimental data of Shitsman (Shitsman, 1963). In Shitsman's experiment, heat transfer of the supercritical water flowing upward in 8 mm ID tubes was investigated.

The computational domains of both circular tubes and annular channels are divided into three connecting sections. The inlet extension is used to achieve fully developed turbulence flow without heat flux. In all cases, the length of the unheated inlet extension is set to be 60 times of the hydraulic diameter. In the heated section, heat flux is uniformly applied at its inner wall. The outlet extension is to ensure no backflow occurs on the outlet surface.

In all cases as shown in **Table 1**, the grids along axial xdirection are uniform and are set to 2 mm. While the grids along the radial r-direction are non-uniform, to increase the accuracy (Cheng et al., 2007; Yang et al., 2007), the grids near the wall are fine enough to make sure that y<sup>+</sup> at the first element to be smaller than 1. In all cases, the size of the first grid is set to be 1e-3 mm, the size ratio is 1.08. The r-direction node numbers are also shown in **Table 1**. The values of y<sup>+</sup> in all cases are <0.3.

### Boundary Conditions

The thermophysical properties are already embed into Fluent, which are based on the database of The National Institute of Standards and Technology (NIST). The boundary conditions in **Table 2** come from Shitsman's experiment. The operating pressure is 25.3 MPa. Inlet is set as mass flow inlet and the


corresponding turbulence intensity and turbulence viscosity ratio are 5% and 10, respectively. Outlet is set as pressure outlet. Heat flux is uniformly applied on the inner wall for both circular tubes and annular channels. Bulk temperature at inlet is 324.5 at subcritical condition. Buoyancy force is taken into account, the acceleration of gravity is 9.8 m/s<sup>2</sup> . Steady-state simulations are conducted for all cases.

TABLE 2 | Boundary and operating conditions.


### Numerical Method Validation

In our previous work, Case C2 was simulated with the same computational domain, flow conditions and turbulence model to study HTD phenomena. Grid independence tests were conducted and the calculated results agreed with Shitsman's experiment results. In this paper, details of the

### RESULTS AND DISCUSSIONS

### Bulk Enthalpies and Temperatures

Using the following equation, the bulk enthalpies can be calculated:

$$h\_b(\mathbf{x}) = h\_{in} + \frac{q \pi ID \mathbf{x}}{GA\_{flow}} \tag{6}$$

where x is the axial position, hin is the inlet bulk enthalpy corresponding to Tin. q is heat flux, ID is inner diameter of circular tubes or annular channels, G is mass flux, and Aflow is cross-sectional area. With the bulk enthalpy and the operating pressure, the bulk temperature can be obtained.

FIGURE 7 | Wall temperature peak positions for different diameter tubes.

### Diameter Effect in Circular Tubes

In Cheng et al.'s recent work (Cheng et al., 2017), two peaks of wall temperature were found in the simulation results of mixed convective heat transfer at the same flow condition corresponding to Case C2 in the present study when HTD occurred, which was consistent with Shitsman's experimental results. The wall temperature distribution in the axial direction is shown in **Figure 5** and the mechanism was explained in detail with axial velocity profiles and turbulent kinetic energy (TKE) distributions in radial direction.

In the present study, the simulated wall temperature distributions in the axial direction in circular tubes with different IDs at the same flow conditions are shown in **Figure 6**. The bulk temperature in the axial direction of 8 mm ID circular tube is presented in this figure. It shows that two wall temperature peaks always occur at present flow conditions, no matter tube diameter increases or decreases. The positions of wall temperature peaks

in the axial direction for different diameter tubes are shown in **Figure 7**. From these figures, it can be seen that with the increase of tube diameter, the first peak of wall temperature moves upstream gradually and its magnitude increases first when inner diameter is ≤16 mm and then decreases when inner diameter is larger than 16 mm. This phenomenon is consistent well with the experiment result of Shiralkar and Griffith (1970) and Zahlan et al. (2015) as presented in **Figure 3B**. For the second peak, it moves downstream gradually as tube diameter increases.

To explore the mechanism behind these phenomena, the fully developed turbulent axial velocity profiles at the heated sections inlet with different inner diameters for circular tubes are extracted and shown in **Figure 8**. In order to facilitate comparison of these fully developed turbulent velocity profiles, the radial positions of the heated walls are aligned to the right. It shows clearly that for the same mass flux and inlet temperature, the boundary layer velocity profiles are different. As tube diameter increases, the maximum velocity in the main flow region gradually decreases and the velocity gradient in the boundary layer also gradually decreases. These conclusions can also be derived by theoretical analysis. In fact, there are numerous empirical fully developed velocity profiles for turbulent circular tube flows. The most famous one is the power-law velocity profile given by:

$$\frac{U}{U\_{\text{max}}} = \left(1 - \frac{r}{R}\right)^{1/n} \tag{7}$$

where Umax is the maximum axial velocity at the centerline. n is a constant, which increases as Reynolds number increases. The average velocity can be derived by integrating Equation (7):

$$\overline{U} \cdot \pi R^2 = \left[ \int\_0^R \left( 1 - \frac{r}{R} \right)^{1/n} \cdot 2\pi r dr \right] \cdot U\_{\max} \tag{8}$$

$$\frac{U}{U\_{\max}} = \frac{2n^2}{(n+1)\left(2n+1\right)}\tag{9}$$

$$U\_{\max} = \frac{\overline{U}}{2} \left( 1 + \frac{1}{n} \right) \left( 2 + \frac{1}{n} \right) \tag{10}$$

$$U = \frac{\overline{U}}{2} \left( 1 + \frac{1}{n} \right) \left( 2 + \frac{1}{n} \right) \left( \frac{R - r}{R} \right)^{1/n} \tag{11}$$

Since the average velocity U is the same for all circular tubes due to the same mass flux and inlet temperature, the exponent n increases with the increase of Reynolds number and diameter, we can conclude from Equation (10) that the maximum velocity in the main flow region decreases with the increase of tube diameter. Based on velocity profile in Equation (11), we can see that at the same position (R-r has the same value in **Figure 8**) away from the heated wall, U decreases monotonously as tube diameter increases. Hence, there is no intersection point for these fully developed turbulent velocity profile when aligned to the right. We can conclude from theoretical analysis that the gradient of axial velocity in the boundary layer gradually decreases with the increase of tube diameter.

According to Kurganov et al.'s experimental study (Kurganov and Kaptil'ny, 1992), direct numerical simulations (Bae et al., 2005, 2008; Bae and Kim, 2009) and Cheng et al.'s analysis (Cheng et al., 2017), HTD occurs in upward flow when the axial velocity profile in the near heated wall region is flattened due to buoyancy

tube, first wall temperature peak at x = 364 mm.

effect and an M-shaped axial velocity profile is formed. At the same time, Reynolds shear stress and TKE in the near heated wall region are reduced significantly.

In the present study, when tube diameter increases, the maximum axial velocity in the main flow region and the gradient of axial velocity in the near heated wall region are both reduced at the heated section inlet. Hence, for larger diameter tubes, the axial velocity profile in the near heated wall region is much easier to be flattened due to buoyancy, an M-shaped axial velocity profile is formed earlier and hence HTD occurs earlier along the axial position, which moves the first peak upstream. The axial and radial velocity, Reynolds shear stress and TKE distributions at four different axial positions (inlet of heated section, 100 mm before, at and 100 mm after the first wall temperature peak) for Case C2, C4, and C6 are shown in **Figure 9**. We can see that at the axial position of the first wall temperature peak, an M-shaped axial velocity profile is formed because of buoyancy. The radial velocity has the maximum value, which means that

FIGURE 10 | Wall temperature distributions in the axial direction in annular channels.

the natural convection is strongest at this position. While the Reynolds shear stress and TKE have the minimum values at the peak position in the near heated wall region as shown in the blue circles, which means that the heat transfer due to turbulence is strongly suppressed. We can also find from these figures that the M-shaped axial velocity profile is formed earlier with the increase of tube diameter.

Regarding the first peak magnitude, it firstly increases, then decreases as tube diameter increases as shown in **Figure 6**. When the diameters are small, the maximum axial velocity decreases significantly as tube diameter increases as shown in **Figure 8**. The decrease of velocity gradient and TKE is relatively large, hence HTD is much more severe with the increase of tube diameter. However, when tube diameters are large enough, the maximum axial velocity and near wall velocity profile changes very little as tube diameter increases. In this situation, the effect of mass flow rate becomes larger than velocity gradient effect, more water is used to cool the heated surface, which results in the HTD less severe and the decrease of the first peak magnitude.

The formation process of the second wall temperature peak is described in detail in Cheng et al.'s recent work (Cheng et al., 2017). In phase III of **Figure 5**, the shear stress makes the axial velocity profile flat in the main flow region and leads to the second peak. Different from phase I, the wall temperature response speed in phase III depends on the tube diameter. With larger tube diameter, the wall temperature response speed is slower and then the second peak moves downstream.

### Diameter Effect in Annular Channels

The simulated inner wall temperature distributions in the axial direction in annular channels with fixed inner diameter and various outer diameters are shown in **Figure 10**. The bulk temperature in the axial direction of the annular channel with 12 mm outer diameter is also shown here. The characteristics of wall temperature in annular channels are quite different from the circular tubes. It is interesting to find that when OD is ≤24 mm, even though the flow conditions and the diameter of heated wall are the same as in the circular tubes, HTD is suppressed and wall temperature increases monotonically. The result agrees with the experimental result of Kim (Kim et al., 2008). More interestingly, it is found that when OD is larger than 24 mm, HTD occurs and wall temperature peaks becomes more and more noticeable with the increase of OD. However, compared with the first peaks in circular tubes, the peaks in annular channels are much smaller and smoother.

The fully developed turbulent axial velocity profiles at the heated sections inlet with different outer diameters for annular channels are extracted and shown in **Figure 11**. Similar to circular tubes, the maximum axial velocity and the velocity gradient in the near heated wall region decreases gradually as outer diameter increases. When OD is small, the velocity gradient is so large that the velocity profile in the near heated wall region cannot be flattened by the buoyancy effect. The Reynolds shear stress and TKE in the near wall region keeps large enough to suppress HTD. With the increase of OD, the maximum axial velocity and the velocity gradient decreases gradually and HTD occurs gradually as shown in **Figure 10**. The axial and radial velocity, Reynolds shear stress and TKE distributions at four different axial positions (inlet of heated section, 100 mm before, at and 100 mm after the wall temperature peak) for Case A6 is shown in **Figure 12**. We can see that a half M-shaped axial velocity profile is formed near the heated wall when HTD occurs in the annular channel. Similar to circular tubes, the Reynolds shear stress and TKE have the minimum values at the peak position in the near heated wall region as shown in the blue circles, which means that the heat transfer due to turbulence is suppressed. However, these minimum values in the near heated wall region in annular channels are much larger than those in the circular tubes, which explains why the peaks in annular channels are much smaller and smoother.

### Comparison Between Circular Tubes and Annular Channels

Even though the hydraulic diameters of Case C2 and Case A2 are both 8 mm, it is interesting to find in **Figure 13A** that the heat transfer behaviors are quite different at the same flow conditions. HTD occurs in the circular tube but not in the annular channel, where the heated wall temperature just increases monotonously. **Figure 13B** shows the fully developed turbulent axial velocity profiles at the heated section inlet for Case C2 and Case A2, where the heated walls are aligned to the right. We can see that the differences of maximum velocity and velocity gradient between circular tube and annular channel are very small. Hence, the reason that results in the different behaviors of heat transfer does not just depend on the initial velocity profile at the heated section inlet. To explain this phenomenon, let's equally divide the flow regions along radial direction in circular tube and annular channel into pieces as shown in **Figure 13C**, where each piece is a small annular flow region with the same thickness. For the yellow and green flow region with the same distance away from the inner wall in **Figure 13C**, the area of the green flow region in annular channel is larger than the yellow one in the circular tube. And the difference becomes more pronounced along the direction of heat flow. Since the heat flux is the same, with larger area of flow region to absorb heat in annular channel, the water temperature in the near heated wall region increases significantly slower than in the circular tube, then buoyancy effect is much weaker and hence HTD is much easier to be suppressed.

## CONCLUSIONS

In this paper, diameter effect on HTD for supercritical water upward flow in circular tubes and annular channels at low mass flux and high heat flux was studied numerically. The simulations

were carried out using SST k-ω turbulence model with Fluent 15.0 software. When the same boundary conditions were applied, i.e., mass flux, heat flux, and inlet temperature, it was observed that for circular tubes the first peak of wall temperature moves upstream while the second peak moves downstream as tube diameter increases. The magnitude of the first peak of wall temperature increases first when inner diameter is ≤16 mm and then decreases when inner diameter is larger than 16 mm. While at the same boundary conditions in annular channels with fixed inner diameters, we found that HTD is suppressed in the situation when the outer diameter is ≤24 mm, and when outer diameter is larger than 24 mm, it was interesting to find that HTD occurs. These phenomena are consistent with previous experiments. The mechanism of these phenomena was analyzed in detail based on fully developed turbulent velocity profiles in the near wall region at the inlet of the heated sections. Increasing inner diameter for circular tubes or outer diameter for annular channels will result in the decrease of maximum velocity and velocity gradient in the near heated wall region, which makes velocity profile in this region much easier to be flattened by the buoyancy. Then an M-shaped axial velocity profile is formed, which will significantly decrease the Reynolds

### REFERENCES


shear stress and TKE and hence impair the heat transfer and cause HTD. At the same flow conditions, HTD is much easier to occur in circular tubes than in annular channels with the same hydraulic diameters.

### DATA AVAILABILITY

All datasets generated for this study are included in the manuscript/supplementary files.

### AUTHOR CONTRIBUTIONS

HC conducted the research and wrote the manuscript. AY worked together with the first author in the research. JW provided lots of valuable suggestions during the research. JZ guided the research.

### FUNDING

This study was funded by Shenzhen Science Technology and Innovation Committee with the project number: JCYJ20150601102053068.


**Conflict of Interest Statement:** AY was employed by company Shanghai Chuangyi Fluid Machinery Co., Ltd.

The remaining 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 © 2019 Cheng, Yanlin, Wang and Zhao. 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.

## NOMENCLATURE


# Design of the Container for the Sampling and Detection Monitoring System of N-13 in Pressurized Water Reactor Primary Loop Water Leakage Based on the Coincidence Method

Yue Zhao\*, Guo-Pu Qu and Jian-Liang Zhou

School of Nuclear Science and Technology, University of South China, Hengyang, China

The N-13 coincidence method is an effective approach for monitoring the primary loop leakage in pressurized water reactors. The high coincidence efficiency and low background are crucial to achieving a lower limit of measurement for such a monitoring system. In this paper, we proposed four types of geometric designs of the sampling composing NaI(TI) detector. For varying container volumes (V), the detection efficiency (ε) of these containers was investigated through the Geant4 simulation and experimental measurement. The value of ε•V, which is a key combined parameter, was obtained accordingly. Finally, we obtained the optimal size of the suitable sampling and detection container for the coincidence method.

#### Edited by:

Muhammad Zubair, University of Sharjah, United Arab Emirates

#### Reviewed by:

Ivo Kljenak, Jožef Stefan Institute (IJS), Slovenia Ming Tang, Los Alamos National Laboratory (DOE), United States

> \*Correspondence: Yue Zhao 430000089377@usc.edu.cn

#### Specialty section:

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

Received: 28 December 2018 Accepted: 08 July 2019 Published: 07 August 2019

#### Citation:

Zhao Y, Qu G-P and Zhou J-L (2019) Design of the Container for the Sampling and Detection Monitoring System of N-13 in Pressurized Water Reactor Primary Loop Water Leakage Based on the Coincidence Method. Front. Energy Res. 7:67. doi: 10.3389/fenrg.2019.00067 Keywords: N-13, coincidence method, sampling container, coincidence efficiency, simulation, Geant4

### INTRODUCTION

During the operation of nuclear power plants, a series of uncertainties, such as the stress relaxation of material and weld joint pressure corrosion, can cause the leakage of the primary loop radioactive water (D'Auria et al., 2017). Such leakage may affect the safety of operation of the nuclear power station, and so it should be accurately monitored. Methods to detect the release of reactor coolant to the containment include the indication and/or monitoring of changes in (a) airborne particulate radioactivity, (b) airborne gaseous radioactivity, (c) containment atmosphere humidity, (d) containment atmosphere pressure and temperature, and (e) condensate flow rate from the air coolers (Dissing and Svansson, 1980; Aoki, 1991; Zheng et al., 2016). Among these methods, a reliable leakage monitor method is to measure the concentration of <sup>13</sup>N in the containment vessel, e.g., by using the γ-ray energy spectrum method. The use of <sup>13</sup>N as tracer isotopes for a reactor coolant pressure boundary leakage detection system has been experimentally demonstrated at the R2 research reactor at Studsvik, Sweden (Dissing and Svansson, 1980). Ling-Qiu et al. applied this method to study the leakage of <sup>13</sup>N in the primary loop of nuclear power stations. The lower limit of detection was 11 L/h, which occurred at the pressure vessel for the 600-MW nuclear power plant in Qin Shan, China. The lower limit of detection of the <sup>13</sup>N energy spectrum method is high, and the small-amount leakage of the first loop is difficult to find. For the energy spectrum method, the detector system requires good stability, and it is challenging to achieve long-term stability for such system. The coincidence method is a common method to measure the radionuclide activity of simultaneously emitting several rays, which can effectively reduce the background, improve the signal-to-noise ratio and reduce the stability requirements of the measurement system (Ember et al., 2004; Vidmar et al., 2007; Antovic and Svrkota, 2009; Volkovitsky and Naudus, 2009; Tillett et al., 2017; Zhang et al., 2018). However, the disadvantages are the low coincidence efficiency and long measurement time. Because <sup>13</sup>N is a short-lived positron emitter, its positron annihilation can produce two gamma photons with equal energy of 0.511 MeV in opposite directions. Hence, we propose a method to measure the concentration of <sup>13</sup>N in the containment using the coincidence method. In this paper, we focus on designing an appropriate sampling and detection container with a high coincidence efficiency. Geant4 simulation and experimental measurements were performed to optimize the detection efficiency of the container for several geometric schemes.

### METHODS

### Theoretical Descriptions

The elastic collisions between a neutron from the nuclear fuel fission and a hydrogen nucleus in water produce recoil protons:

$$
\mathfrak{l}\_I^I H + \mathfrak{n} \to \mathfrak{n} + \mathfrak{p} \tag{1}
$$

When the produced protons have energy higher than 5.5 MeV, they may interact with <sup>16</sup>O and produce <sup>13</sup>N via:

$$\rm{H}\_{8}^{16}O + p \rightarrow {^{13}\_{7}N} + \alpha \tag{2}$$

Here, <sup>13</sup>N is a positron emitter, which has a half-life of 9.96 min. The emitted positrons lose their energy when they pass through matter, annihilate with abundant electrons surrounding them, and decay into two γ-ray protons with equal energies of 0.511 MeV.

The concentration of <sup>13</sup>N (defined as N1) in primary loop water is proportional to reactor power P:

$$N\_I = K\_I P \tag{3}$$

where K<sup>1</sup> is the coefficient in the unit of L−<sup>1</sup> (MW)−<sup>1</sup> .

If the water in the primary loop leaks, it can be gasified due to its high temperature and pressure, and <sup>13</sup>N exists in a gaseous form in the containment vessel. The gas in the vessel was pumped into a sampling container and measured with the coincidence system. The leakage rate V<sup>L</sup> is obtained from the net coincidence counting rates (n), transmission coefficient (K2), and reactor power (P) as follows:

$$\text{V}L = \frac{n}{\text{K}\_2\text{N}\_I} = \frac{n}{\text{K}\_I\text{K}\_2\text{P}} \tag{4}$$

Here, the units of the true coincidence rate and leakage rate are cph and L/h, respectively.

Transmission coefficient K<sup>2</sup> is the true coincidence rate measured by the coincidence system when the leakage rate of the primary loop radioactive water is 1 L/h in the containment. K<sup>2</sup> depends on many factors such as coincidence efficiency ε of the gamma ray with an energy of 0.511 MeV, gasification of water leaked from the primary loop, uniform mixing time t1, gas transmission time of gas in the sampling pipe (t2) and sampling container (t3) as follows:

$$K\_2 = \frac{\varepsilon \mathcal{Q} e^{-\lambda t\_I} \ e^{-\lambda t\_2} (I - e^{-\lambda t\_3})}{\lambda \, V\_I} \tag{5}$$

where λ is the decay constant of <sup>13</sup>N in units of h−<sup>1</sup> , Q is the sampling flux in units of L/h, and V<sup>1</sup> is the volume of hood containment of the reactor driving mechanism ventilation in units of L. As a result, leakage rate V<sup>L</sup> can be expressed as:

$$V\_L = \frac{m\lambda V\_I}{\varepsilon Q e^{-\lambda t\_1} e^{-\lambda t\_2} (1 - e^{-\lambda t\_3}) K\_I P} \tag{6}$$

When the reactor operates at full power, by substituting Equation (3) into Equation (6), the leakage rate becomes:

$$V\_L = \frac{n\lambda \, V\_I}{\varepsilon \, \mathbf{Q}e^{-\lambda t\_I} e^{-\lambda t\_2} (1 - e^{-\lambda t\_3}) N\_I} \tag{7}$$

The lower limit is an important parameter for the detection instrument system. At 95% of confidence probability, if only statistical fluctuations are considered, lower limit of net count N<sup>D</sup> can be expressed as (Knoll, 2000):

$$N\_D = 4.65\sqrt{N\_b} \tag{8}$$

Where N<sup>b</sup> is background counts. To convert N<sup>D</sup> in counts to minimum detectable activity A, the additional factors of radiation yield per disintegration (f), measurement time of background t<sup>b</sup> and absolute detection efficiency (ε) must be taken into account (Knoll, 2000):

$$A = \frac{N\_D}{f\varepsilon t\_b} = \frac{4.65\sqrt{N\_b}}{f\varepsilon t\_b} \tag{9}$$

The activity of the gas in the sampling container (A1) is determined by concentration of <sup>13</sup>N (N) in the gas and volume V of the sampling container. It takes a form:

$$A\_1 = \lambda N V = \lambda K V\_L V \tag{10}$$

Here K is a constant.

When A<sup>1</sup> equal to A, the lower limit of leakage rate VLD can be expressed as:

$$V\_{LD} = \frac{4.65\sqrt{\frac{N\_b}{t\_b^2}}}{Kf\lambda\_\varepsilon\varepsilon V} \tag{11}$$

During the measurement, Kfλ can be considered a constant. From Equation (8), we find the designing sampling and detection container, which can enhance the coincidence efficiency and ε•V, which is a key issue to reduce the lower limit of detection of the 13N leakage monitoring system.

### Preliminary Design of the Sampling and Detection Container

The Marinelli beaker (MB) as a sampling container has been used for the <sup>13</sup>N leakage monitoring system based on the γ-ray energy spectrum analysis (Dissing and Svansson, 1980; Erbeszkorn et al., 1996; Ahmed et al., 2009). The design parameters of the MB container enable one to detect gamma ray with energy of 0.511 MeV. We design the current sampling and detection container based on the MB container. Since the positron annihilation can produce two gamma photons in opposite directions, the sampling container is selected as a sealed stainless-steel (1Cr18Ni9Ti) cylindrical container. The stainless-steel container is 1 mm thick. The inner and total heights are 78 and 80 mm, respectively. Five Φ3"×3"NaI(TI) detectors were used as coincidence detectors. Among them, one was installed at the center of the container, and the others were symmetrically distributed around the container. To obtain an optimal coincidence detection efficiency of the container, four types of geometric designs of the container were selected: four detectors were installed outside the container (NO.1), a quarter of the detector circumference was located in the container (NO.2), the center of the detector in the circumference of the container (NO.3) and all detectors in container (NO.4). A schematic of these sampling and detection container arrangements is shown in **Figure 1**.

### Detection Efficiency Simulation

The Geant4 toolkit (Agostinelli et al., 2003) presents an objectoriented programming that enables one to select among a wide range of physical processes or implement them according to

TABLE 1 | Technical dimensions of detector provided by the manufacturer.


experiment requirement. In this work, Geant4 was used to simulate the entire sodium iodide detector assembly, materials, and sampling containers with different volumes, as shown in **Figure 2**. The dimensions of the sodium iodide crystal, thickness of reflective-layer magnesium oxide, and thickness of cladding aluminum were provided by the manufacturer, as shown in **Table 1**; the thickness of the stainless steel was physically measured using a Vernier caliper. The physical processes activated in the detector model of this work are: Compton scattering, pair production, photoelectric effect for photons, ionization processes, multiple scattering and Bremsstrahlung for electrons, and ionization processes, multiple scattering and annihilation for positrons. The <sup>13</sup>N nuclei were uniformly and randomly generated in the sampling container. To ensure that the simulation was most similar to the experimental measurement, a useful signal was taken when the energy deposition of gamma ray in the NaI(TI) crystal exceeded 150 keV, and the coincidence only formed only when any two detectors simultaneously generated useful signals. Sampling numbers of 10<sup>5</sup> , 10<sup>6</sup> , 10<sup>7</sup> , and 10<sup>8</sup> were used to test and ensure statistical uncertainties.

### Experiment Measurement Gaseous <sup>13</sup>N Production

Gaseous <sup>13</sup>N was produced with cyclotron protons (whose energy was 16.5 MeV) by bombarding deionized water in the target chamber. The bombardment lasted approximately 10 min. The produced <sup>13</sup>N existed in the form of liquid <sup>13</sup>N-NO<sup>−</sup> 2 and <sup>13</sup>N-NO<sup>−</sup> 3 . After pressured by dense argon, they were delivered to a collection container with Devarda's alloy and NaOH. The collection container was also shielded by a lead chamber. After a few minutes, the produced <sup>13</sup>N-NH<sup>3</sup> in the collection container was transferred into a 5 ml vacuum bottle. The activity of <sup>13</sup>N was measured using an activity meter CRC-25R produced by CAPINTEC.IN, and activity A<sup>0</sup> was recorded. Furthermore, <sup>13</sup>N-NH<sup>3</sup> was pushed under the effect of CCl<sup>4</sup> as a booster and sampling container for the following measurement of the coincidence efficiency, which will be later discussed. A detailed flow chart of the <sup>13</sup>N production is shown in **Figure 3**.

### Sampling Container and Measurement System

To verify the accuracy of the simulation detection efficiency, the detection efficiency of <sup>13</sup>N was experimentally studied in the laboratory. Four arrangements of containers with the capacity of two liters were produced. All parameters are identical to those in the preliminary design, as discussed. The radii of four types of containers (NO.1–NO.4) were 100, 110, 120, and 130 ± 1.0 mm. The sampling containers and detector are shown in **Figure 4**. The measurement system consists of five Φ3" × 3"NaI(TI) detectors, which were purchased from Bei Jing ZhongGuang Detector CO.LTD CHINA, and a data acquisition system, which

was purchased from Beijing Nuclear Instrument CO.LTD. The detector and sampling container are shown in **Figure 4**.

Five Φ3" × 3"NaI(TI) detectors were installed into locations of the sampling container. The background of the coincidence detection system and coincidence efficiency were measured without the lead shield. In our experiments, the threshold of a single channel depends on the pulse amplitude of gamma rays of 0.511 MeV. The lower threshold of a single channel was set to 20% of the amplitude of the pulse, and the upper threshold was 120%. The input pulse width in the coincidence circuit was set to 0.5 µs, which satisfies the requirement of time resolution for the discussed coincidence measurement. According to the arrangement of five detectors and characteristic of the positron annihilations, we see that even if there is a radioactive source, it has very low probability to form true coincidence between detector 1 and detector 3 or between detector 2 and detector 4. To reduce the background of the coincidence system, these two types of coincidence events are considered negligible.

## RESULTS AND DISCUSSION Data-Based Geant4 Simulation for the Efficiency

Due to the limitation of space in the containment, the sampling container volume of the primary loop leakage monitoring system is small, so the max radius of the container is only 16 cm in the simulation calculation. For different container radii, coincidence efficiency ε and ε•V are shown in **Figures 5**, **6**. The coincidence efficiency decreases with the increase in container radius due to the decreased solid angle of the NaI(TI) detector. Unlike the coincidence efficiency, ε•V first increases and subsequently decreases with the container radius. Maximal coincidence efficiency and ε•V are obtained to be 11.55% and 0.103, respectively. According to the relationship between detection lower limit VLD and ε•V, we can find that when the fourth container radius is R = 130 mm, the detection limit is minimal in all simulated sampling containers, which is theoretically considered the optimal sampling container size.

### Background of Single Detector and Coincidence System

The backgrounds of the single detector and coincidence detection system were measured without the lead shield within 10,000 s. The results for three repeated measurements and their average value are shown in **Table 2**. The coincidence measurement system rapidly suppresses the background count rate for a single TABLE 2 | Background results.


detector to the order of 0.1 Hz, which is only 1/660 times of that from a single detector.

### Coincidence Efficiency

Under the same condition of measuring the background coincidence, we performed the coincidence measurements. The measurement integral is selected as 300 s. For the short-lived <sup>13</sup>N isotope, its activity visibly decreases during the measurements. Both numbers of decayed nuclei in the measurement period and true coincidence events should be considered to calculate the coincidence efficiency, which takes a form:

$$\varepsilon = \frac{n\lambda}{A\_0(e^{-\lambda(t\_1 - t\_0)} - e^{-\lambda(t\_2 - t\_0)})} \times 100\% \tag{12}$$

Here, n is the true coincidence events, λ is the decay constant, A<sup>0</sup> is the activity of <sup>13</sup>N at time t0, and t<sup>1</sup> and t<sup>2</sup> are the initial and end time of the coincidence measurement, respectively. The results for three repeated measurements, which correspond to three types of <sup>13</sup>N activities are shown in **Table 3**. For different <sup>13</sup>N activities, the coincidence efficiencies obtained with Equation (9) are consistent with one another within systematic and statistical uncertainties.

Coincidence efficiency ε and ε•V of the four containers are summarized in **Table 4**. The comparison between simulation and experimental results shows a slight discrepancy of ≤5%, which further validates the Geant4 predictions. This difference may be caused by the effect of <sup>18</sup>F and <sup>15</sup>O produced from the cyclotron accelerator, since these short-lived isotopes with tiny quantity are present in the liquid of <sup>13</sup>N-NO<sup>−</sup> 2 and <sup>13</sup>N-NO<sup>−</sup> 3 due to the lack of separation.

In our work, four arrangements of containers with the capacity of two liters were used in the experiment, and the results are in good agreement with simulation. However, due to the lack of experimental or theoretical verification, the results of other volumetric containers are uncertain. If the method was applied to engineering, further experiments are needed to verify the reliability of the results.

### CONCLUSION

Based on the <sup>13</sup>N decay characteristic, we proposed the gammagamma coincidence method to measure the concentration of



TABLE 4 | Comparison of the simulated and measured results.


<sup>13</sup>N. This coincidence method can be applied to monitor the leakage of <sup>13</sup>N from the primary loop of nuclear power plants. To improve the coincidence efficiency and ε•V value, four types of sampling and detection containers were designed. The coincidence efficiencies in four arrangements of sampling and detection containers were investigated through Geant4 simulations and experimental measurements. The results show that when the NO.4 sampling and detection container (five detectors in the container) has radius R = 130 mm, the coincidence detection system has a higher detection efficiency and maximum ε•V value for <sup>13</sup>N. Thus, the NO.4 container of this size has optimal arrangement sampling and detection. This work was performed in the laboratory, and the effects

### REFERENCES


of <sup>18</sup>F and <sup>15</sup>O produced by the accelerator were ignored in the experiment.

### AUTHOR CONTRIBUTIONS

YZ completed simulation, experiment research, and wrote the manuscript. G-PQ and J-LZ completed experiment research and revised the article.

### FUNDING

This work was supported by the National Science Foundation of China under Grant No. 11175083.


**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 © 2019 Zhao, Qu and Zhou. 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.

# Research on Time-Dependent Failure Modeling Method of Integrating Discrete Dynamic Event Tree With Fault Tree

#### Anqi Xu<sup>1</sup> , Zhijian Zhang<sup>1</sup> \*, Min Zhang<sup>2</sup> , He Wang<sup>1</sup> , Huazhi Zhang<sup>1</sup> and Sijuan Chen<sup>1</sup>

*<sup>1</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China, <sup>2</sup> China Nuclear Power Engineering Co., Ltd., Beijing, China*

#### Edited by:

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

#### Reviewed by:

*Weiwen Peng, Sun Yat-sen University, China Tangtao Feng, University of Wisconsin-Madison, United States Han Bao, Idaho National Laboratory (DOE), United States*

\*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: *29 May 2019* Accepted: *19 July 2019* Published: *13 August 2019*

#### Citation:

*Xu A, Zhang Z, Zhang M, Wang H, Zhang H and Chen S (2019) Research on Time-Dependent Failure Modeling Method of Integrating Discrete Dynamic Event Tree With Fault Tree. Front. Energy Res. 7:74. doi: 10.3389/fenrg.2019.00074* Classical PRA methods such as Fault Tree Analysis (FTA) and Event Tree Analysis (ETA) are characterized as static methods due to predetermined event sequences and success criteria of frontline systems. They are widely accepted for risk analysis of nuclear power plants. Unlike classical PRA, Dynamic PRA (DPRA) couples the stochastic random failures of system with deterministic analysis (by simulation) to determine the risk level of complex systems. It considers the safety significance of the timing and order of events on accident progression and consequences. However, it is time-consuming to establish a complicated full-scope system simulation model. Meanwhile, thousands of accident scenarios are generated due to randomness of state transition, uncertainty of model and parameters. An overload of modeling, calculating, and post-processing will arise. So, it is a prospective and challenging idea to integrate the classical PRA method with the dynamic PRA method. The objective of this paper is to address an integrated method of risk quantification of accident scenarios. It points out how to treat time-dependent interactions of accident dynamics including random failures, temporal events, configuration changes, and physical process parameters explicitly. Possible dependencies and configuration consistency issues accounting for Discrete Event Tree (DET) branch probabilities are discussed. For DET simulation, some of non-safe-related components to be analyzed could be modeled by FTs for conditional branching probability, instead of a computationally expensive simulation model. A method of integrating FT into DET is introduced which emphasizes on computing the conditional branch probability with FTs online, as well as developing a DET model in case of temporal relations of failure. Finally, a simple case of a Low Pressure Injection System in Large Break Loss-of-Coolant Accident (LBLOCA) scenario is provided as a demonstration.

Keywords: dynamic event tree, fault tree, hybrid PRA, time-dependent, performance-based

### INTRODUCTION

### Overview of Classical PRA and Dynamic PRA Method

The classical PRA methods, such as Event Tree Analysis (ETA), Fault Tree Analysis (FTA), Reliability block Diagram (RBD) are widely accepted throughout nuclear industry. ET/FT are generally based on static Boolean logic structure. ET is used to inductively model the accident progression to dictate all possible accident sequences based on engineering judgment and thermo-hydraulic analysis. FT is used to deductively model the system failure by a top-down, hierarchical tree to analyze all the possible combinations of failure events. But PRA methodology faces challenges including the treatment of time dependent interactions (accident dynamics) and the propagation of physical process uncertainties to risk.

Dynamic PRA (DPRA) have been developed since 1980s. Under the framework of DPRA, many methods can lead to a more realistic risk assessment for nuclear power plant, originating from Dynamic Event Tree Analysis Method (DETAM; Acosta and Siu, 1992; Siu, 1994), Dynamic Logical Analytical Methodology (DYLAM; Cojazzi, 1996), Dynamic Event Tree (DET; Acosta and Siu, 1993). DPRA evaluates the timing and sequencing of events in accident progression and identifies the failure paths under all possible accident scenarios. DPRA is also treated as simulation-based PRA (Mosleh, 2014), or Integrated Deterministic and Probabilistic Safety Assessment (IDPSA), and related reviews and literatures can be found in the references (Aldemir, 2013; Zio, 2014). Among DPRA methods, DET can partially solve the problem of timing and ordering by coupling the stochastic analysis (reliability) with accident simulation. Along with theoretical research, mature modeling, and computational tools of DPRA have been developed for risk quantification and uncertainty analysis, like Accident Dynamics Simulator (ADS; Hsueh and Mosleh, 1996), Analysis of Dynamic Accident Progression Trees (ADAPT; Hakobyan et al., 2008), Simulation-based PRA (SimPRA; Mosleh et al., 2004), Risk Analysis Virtual Environment (Alfonsi et al., 2017), so as to assist PRA practitioners in improving the modeling and computing efficiency. However, DPRA still suffers from several difficulties:


The first challenge has been partly overcome by time discretization, process discretization, appropriate parameter sampling strategies, branching rules, pruning rules, or leveraging the similar branches. The investigation of the third issue is still in progress, such as the traditional scenario clustering algorithm (Mandelli et al., 2013) and topological clustering algorithm (Maljovec et al., 2016). In addition, the extensive computing problem can be addressed by parallel processing and configuring high-performance computing resource. But for the second issue, a refined simulation model is not required for some components of non-safety systems or supporting systems in DPRA which could reduce the work of modeling and debugging.

Currently, DET is one of the most prospective methods in discrete-time DPRA. DET is similar to traditional ET as it evolves the accident sequences by spawning branch events. However, the difference is that the timing and order of ET accident sequences following an initiating event is predetermined by PRA analysts without any interaction of system responses, while that of DETs are determined "online" by a time-dependent system evolution model and a set of branching rules. In general, DETs are classified into two kinds depending on the treatment of continuous stochastic variables: (1) discretization, such as Discrete Dynamic Event Tree (DDET); (2) sampling-based, such as Monte Carlo Dynamic Event Tree (MCDET; Kloos and Peschke, 2006) DDET discretizes continuous aleatory variables (e.g., failure time or recovery time of equipment, operator response time, etc.) into various branches, while MCDET adopts Monte Carlo sampling with DET simulation (with each sample represented by a DDET). This distinction could impact the accuracy of final results when many continuous variables are concerned in the accident sequences.

In recent years, the concept of "hybrid PRA" (Mandelli et al., 2018) has been proposed as a new methodology which refers to the integration of classical PRA into dynamic PRA. The integrated model of DET and FT has been discussed in previous dynamic PRA research. For DET, the conditional probability of header events/frontline system is regarded as a DET branching probability, which accounts for supporting system dependencies and other dependencies among the safety functions in DET simulations, as mentioned in section Treatment of Dependencies in DET Modeling During Accident Analysis. In most of current DET research, the conditional probabilities are estimated offline (Chao and Chang, 2000) before simulation, in order to simplify the DET simulations. An established FT or a known probability distribution is used to represent the DET branch probability of a certain system state. The limitations of conditional probabilities computed offline is that it ignores the dependencies of conditional branch probability on process variables or operating conditions. Thus, it is not numerically accurate enough.

**Abbreviations:** ACC, Accumulator; ADS, Accident Dynamics Simulator; CCF, Common Cause Failure; CCS, Component Cooling Water; CDF, Cumulative Distribution Function; DDET, Discrete Dynamic Event Tree; DET, Dynamic Event Tree; DETAM, Dynamic Event Tree Analysis Method; DYLAM, Dynamic Logical Analytical Methodology; DPRA, Dynamic Probabilistic Risk Assessment; ECCS, Emergency Core Cooling System; ET, Event Tree; ETA, Event Tree Analysis; FT, Fault Tree; FTA, Fault Tree Analysis; HPI, High Pressure Injection; HX, Heat Exchanger; IDPSA, Integrated Deterministic and Probabilistic Safety Assessment; INL, Idaho National Laboratory; LBLOCA, Large Break Loss-of-Coolant Accident; LPI, Low Pressure Injection; LPR, Low Pressure Recirculation; MCDET, Monte Carlo Dynamic Event Tree; MCS, Minimal Cut Set; MCSQ, Minimal Cut Sequence; MGL, Multiple Greek Letter; NPP, Nuclear Power Plant; PAND, Priority-AND; PBRI, Performance-based reliability indicator; PRA, Probabilistic Risk Assessment; RAVEN, Risk Analysis Virtual Environment; RBD, Reliability block Diagram; RCSS, Reactor Containment Spray System; SCAIS, Simulation Code System for Integrated Safety Assessment; SEQ, Sequence Enforcing Gate; SF, Sequential Failure; SimPRA, Simulation-based PRA.

As an alternative, the conditional branch probabilities of DET can be calculated online. It intends to take advantages of classical PRA for treating dependencies and simplifying simulation modeling work (Karanki and Dang, 2016). Even though basic ideas are comprehensible, it is complex to couple the Boolean algorithms and probabilistic calculation with dynamic accident progression. Simulation Code System for Integrated Safety Assessment (SCAIS; Ibánez et al., 2015; Izquierdo et al., 2017) is one of such DET computational tools, but the consideration of how to calculate the conditional probability online in DET simulation is not thoroughly discussed, especially when the branch probability is time-dependent or effected by process variables. Also, this issue is not clearly clarified in other DET software tools.

### Goal of This Paper

The objective of this paper is to investigate how to incorporate classical PRA models into simulation-based PRA in a consistent manner. It addresses a new framework of risk quantification of accident scenarios. The integration method points out how to treat time-dependent interactions of accident dynamics including random failures, temporal failure events, system configuration changes, and physical process parameters explicitly. Unlike most of DET previous literatures, this paper does not estimate the time-dependent conditional branch probabilities offline before DET simulation but intends to update FT online in order to obtain DET simulation results more realistically. So, the issue of time-dependent failure modeling is illustrated including:


The paper is organized into the following five sections. The mathematical explanation of integration FT into DET with time-dependent and condition-dependent characteristics is illustrated in section Mathematical Basis of Integration Method. The treatment of potential dependencies in DET modeling and their effect of the DET branch probability is presented in section Treatment of Dependencies in DET Modeling During Accident Analysis. Section Integration Method of FT into DET focuses on technical solutions of time-dependent updating, system configuration updating and temporal relations of failure. Finally, a simplified case is provided for demonstration in section Case Study.

### MATHEMATICAL BASIS OF INTEGRATION METHOD

The PRA analysis which describes the failure events of a system or equipment is characterized by specific state, process of state transition, while the mechanics simulation of system provides a time-dependent system response of operating parameters during accident evolution. The corresponding control logic is required in Technical Specification, Accident Procedures, etc. From an integrated point of view, an integrated model is developed which couples probabilistic risk analysis with plant condition analysis.

The system state at time t in phase space is characterized by a set of continuous variables and a set of discrete variables as shown in **Figure 1**. Generally, the continuous variables refer to process variables that have evolved with time, such as pressure and temperature. The discrete variables mainly represent stochastic behaviors such as the random failure of components, uncertainty in the initial and boundary conditions, uncertainty of the system model, etc.

The state of system SE(t)=[θE(t), Ec(t)] can be described by the following equations (Rabiti et al., 2013),

$$\begin{cases} \frac{\partial \vec{\theta}(t)}{\partial t} = \Xi(\vec{\theta}(t), \vec{c}(t), t) \\\\ \frac{\partial \vec{c}(t)}{\partial t} = \Gamma(\vec{\theta}(t), \vec{c}(t), t) \\\\ \vec{\theta}(t\_0) = \vec{\theta}\_0 \\\\ \vec{c}(t\_0) = \vec{c}\_0 \end{cases} \tag{1}$$

Where θE(t) describes the NPP status vector of operating parameters, such as the primary pressure and temperature. Ec(t) is a vector of discrete states for all components, such as operating or failure state, valve open or closed. It is also named as system configuration, i.e., possible combinations of component state. 4 represents the simulation model of system, which is used to calculate continuous variable, especially operating parameters response. Generally, 4 is a multi-physics and multi-scale model of thermal-hydraulic, neutronics, material aging, etc. Γ is a timedependent function of equipment state transitions. It describes the randomness of component states as a probability function. Meanwhile, it is regulated by the control logic (e.g., setpoint values) of a system in operating procedures. For instance, if the water level of pressurizer exceeds an opening threshold, the relief valve is required to action from closed to open.

For a continuously operating device, the time-dependent failure probability can be described by Equation (2).

$$\begin{aligned} F(t) &= P(\mu < t) \\ &= 1 - \exp(-\int\_0^t \lambda(\mu) d\mu) \end{aligned} \tag{2}$$

Where λ(t) is the failure rate of component, t is working time.

The reliability method used in classical PRA is "failurebased reliability" modeling. In this method, the accumulation of failure time and other information are required. Furthermore, a priori life probability function should be determined, and the parameters of probability function should be quantified based on the collected failure data. As a result, the failure-based reliability model is developed to evaluate the remaining lifetime and reliability level of components. Basically, all the failure modes can be classified as traumatic failures and degradation failures. The basic assumption of failure-based reliability modeling is that all failures occur at an instant or in a short period of time. According to this rule, all the failure modes can be regarded as traumatic failures. However, this rule is not always correct because 70–80% of failure modes of mechanical components belongs to degradation failures.

Concerning traumatic and degradation failures, the research on failure probability related to operating environment of system is called "performance-based reliability." It is based on performance-based reliability indicators (PBRI), that is, some of the operating variables which are strongly dependent on reliability.

In fact, reliability parameters are usually both conditiondependent and time-dependent, where **X**t=[X1, X2, .., Xm] T is the vector of PBRI. The failure rate function equipment λ(t,**X**t) is defined as

$$\lambda(t, \mathbf{X}\_t) = \frac{f(t, \mathbf{X}\_t)}{R(t, \mathbf{X}\_t)} \tag{3}$$

Where f(t,**X**t) is the probability density function of equipment, R(t,**X**t) is the reliability function of equipment.

Under different performance levels **X**t1 and **X**t2, it is assumed that the ratio of failure rate is a constant related to performance level at time t. Then the relationship between failure rate and performance level is in accordance with the Proportional Hazard Model (John, 2008)

$$\frac{\lambda(t, \mathbf{X}\_{t1})}{\lambda(t, \mathbf{X}\_{t2})} = \mathbf{C} \tag{4}$$

Where C is a constant which does not vary with time, only influenced by performance level.

Then the traumatic failure rate λ(t,**X**t) can be split into two parts as

$$
\lambda(t, \mathbf{X}\_t) = \lambda\_{r0}(t)\mathbf{g}(\boldsymbol{\mathfrak{g}}, \mathbf{X}\_t) \tag{5}
$$

Where λr0(t) is the baseline failure rate which is a function of time alone, independent of performance variables. β is a vector of coefficients, each element of which indicates the importance or weight of each variable of **X**<sup>t</sup> . g(β,**X**t) is a function of performance variables **X**<sup>t</sup> and coefficient β.

To be simplified, a general form considered for g(β,**X**t) is

$$\log(\mathfrak{g}, \mathbf{X}\_t) = \exp\{\sum\_{i=1}^m \beta\_i X\_i\} \tag{6}$$

Where X<sup>i</sup> is an element of **X**t=[X1, X2, .., Xm] T , i = 1, 2, . . . , m.

Therefore, the unavailability function is described as

$$\begin{aligned} F(t, \mathbf{X}\_{\ell}) &= P(s < t \, | \mathbf{X}\_{\ell}) \\ &= 1 - \exp(-\exp(\sum\_{i=1}^{n} \beta\_{i} \mathbf{X}\_{i}) \int\_{0}^{t} \lambda\_{r0}(s) ds) \end{aligned} \tag{7}$$

Equation (7) points out the FT model of continuous running components is not only related with the system configuration and state duration, but also dependent on the variation of process variables. However, it is quite difficult to quantify the value of β<sup>i</sup> corresponding to X<sup>i</sup> in engineering practice because of the coupling relationship between failure mechanism and performance level, especially for mechanical components.

### TREATMENT OF DEPENDENCIES IN DET MODELING DURING ACCIDENT ANALYSIS

In order to simplify the DET modeling process, one of the generic modeling strategies is to select failures of pivotal functional events/frontline systems as the headers of DET branch. DET branch probabilities are supposed to account for dependencies and configuration consistency among the safety functions in DET simulations. The dependencies here include the functional dependency, configuration dependency, component failure dependency, and human error dependency. Other dependencies are explicitly accounted for as the event tree headers, such as the correlations between functional events, correlations among initiating event and each functional events.

The functional dependency in PSA practice mainly refers to the supporting system dependencies among functions which result from shared supply and shared components. The safety systems in NPP have shared electrical power, water supply, and cooling systems. In terms of the "Small DET- Large FT" strategy, the unavailability of supporting systems are modeled as common modules of FT, and treated by quantitatively determining the Minimal Cut Sets (MCS) of each sequence. Also, it is not recommended to be explicitly modeled in DET headers, so as to reduce the number of branches. In addition, the operating conditions and state duration of subsequent functioning systems depend on the continuous operating time and performance level of previous safety systems, the response time of required human actions. These timing and conditions are not considered in the classical ET method.

The configuration consistency means that the unavailable components during the accident remains irreparable until the end of simulation time except for offsite power recovery and EDG recovery. The component failure dependency includes the common cause failure (CCF) of components due to similar design defects, process of fabrication or installation, procedures of operation. CCF events have been employed in the PRA model to represent all possible dependent reasons in ET/FT logical models. Some of CCF methods (NUREG/CR-6268, NUREG, 2007) have been applied in PRA modeling, such as β Factor Model, Multiple Greek Letter (MGL) model and Alpha Factor Model. Another type of component failure dependency is sequential failure (SF) of components, which are normally not explicitly considered in classical PRA, except for employment of some dynamic gates, such as Priority-AND Gate (PAND) Gate, Sequence Enforcing Gate (SEQ). Accordingly, the CCF modeling and SF modeling also need to be considered in DET.

As for the human error dependency, in order to simplify the analysis and calculation process, the correlation between human tasks is usually divided by specifying correlation parameters. In addition, the probability of human error in the task should be modified according to the level of correlation.

## INTEGRATION METHOD OF FT INTO DET

The INL report (Mandelli et al., 2018) explains how the four main classical PRA methods (Markov, ET, FT, and RBD) is extended

TABLE 1 | Mapping rules of functional failure simulation for CCS pump.


to time domain. It gives a case study of the LBLOCA accident to demonstrate how FTs can be linked to RELAP5-3D PWR model and RAVEN Ensemble Model, and compares the classical PRA results with dynamic PRA. The considered FTs of LBLOCA are of four main frontline systems, i.e., Accumulator (ACC), Low Pressure Injection (LPI) System, High Pressure Injection (HPI) System, and Low Pressure Recirculation (LPR) System. For each frontline system, basic events indicate the states of trains or components, such as failure on demand, unavailable due to test or maintenance, fail to open, etc. Monte-Carlo sampling method is employed to sample the set of stochastic parameters from the calculated probability distribution of basic events. As a result, a series of simulations are generated for the calculations of safety parameters. However, the probability distributions of basic events have been determined before RELAP5-3D simulation. Although the "mission time" of safety system is not illustrated, it is usually set as 24 h in classical PRA practice conservatively. Thus, the "mission time" of safety system is not equivalent to the commissioning time of a frontline system for mitigating accident successfully. Besides, all the branch conditions and configuration changes are embedded into the simulation. The states and state delay time of ACC, HPI, LPI, and LPR are mapped into RELAP5-3D input file. Therefore, INL integration method is only applicable when the probability distribution of stochastic parameters would not change with process variables, that is, the failure probability of a component/train is "failure-based," not "performance-based."

To further extend the integration method of INL, a combination of DDET and FT is proposed with online calculation of conditional branch probability. The basic assumptions of events' behaviors are:


### Key Technical Issues of Integration Method System Configuration Changes Due to Branching Rules

Unlike ET, DET is more flexible to allow for multiple branches, if necessary. The branching rules include:

a. The branching occurs whenever the control logic is fulfilled, setpoint values of system are exceeded or operator actions are required. These requirements mostly exist in NPP design, Technical Specifications, and Emergency Operating Procedures. The system parameters evolve with time for each branch. Based on the possible outcome of system response, new demands of frontline systems or human actions are generated to perform safety functions. Thus, descendent branching occurs and generates more scenarios because of different events at different timings.


The state transition of components determines the system configuration of each branch. For probabilistic analysis, it is required to update logic values of events. For the simulation model of the system response, a set of mapping rules of the functional failure simulation are proposed to update the states of spawned branch nodes according to system configuration changes. The mapping rules transfer the state of components into related operating parameters or control variables. The implementation of such rules is introduced as follows. Firstly, the main functional requirements of the system should be identified to obtain a set of system parameters that characterize its functional state. Secondly, the function(s) of the system is decomposed according to a hierarchical level of system-subsystem-component. Thirdly, FMEA analysis for each of the main functional components is carried out to determine the characteristic parameters of components. Equivalently, the system simulation model can be updated using these characteristic parameters or control variables in case of configuration changes. For example, Component Cooling Water (CCS) System is designed to provide adequate cooling water for certain systems with nominal temperature and flow rate. It is composed of two redundant trains A and B. Each train has two redundant 100% CCS pumps and two 50% heat exchangers (HX), one tank and related valves. In normal conditions, CCS operates with one train. To ensure the safe operation, the normal states of the components/system are guaranteed by a set of characteristic parameters in **Figure 2**, thereby one of the ways to update the system configuration in simulation model is shown in **Table 1**.

### Time-Dependency on Conditional Branch Probability

Unlike ET, the functional demand of a header event/frontline system in DET is determined by a simulation model, not by predetermined success criteria of system. Furthermore, the actual operating time between DET nodes is determined according to actual accident progression, which is noted as "online." In this view, the conditional branch probability can be understood as "the probability of a specific system configuration, given a set of history events under an accident scenario." For discrete variables of safety functions, such as valve open or closed, the probability distribution of demand failure of a component is a Bernoulli distribution. For continuous variables of safety functions, such as operating failures, or human action failures, the probability distribution should be updated including logic values and state duration of events in an accident scenario. The duration of a certain state is from its demand moment to the terminated moment, where "termination" occurs because of exhaustion of water, power, or random failures of components.

The demand failures occur when a component state transfers from a non-functional state to functional state which include:


Thus, the demand failure probability FSD consists of two parts:


*FB, standby failure FD; failure on demand; FW, operating failure; RC, refuse to close; RO, refuse to open; KP, failure to keep the position; RU, rupture; BK: block; EL, external leakage; IL, internal leakage.*

*TOP is the running time interval from the last period test to the moment of IE.*

*t<sup>1</sup> is the time when the demand of component is generated since IE occurs.*

*t<sup>2</sup> is the time when functional demand of component ends since IE occurs. It is determined by DET branching conditions and simulation results.*

λ*<sup>F</sup>* (*t*) *is the failure rate of corresponding failure events.*

*Qs*(*t*) *is the unavailability of a switch-type component. Qs*(*t*) *is regarded as a constant value.*


For online-monitored components, the failure state can be monitored by Digital I&C system of NPP in a quite short period after it occurs. It is guaranteed that the state of component is available when IE occurs. So FSD can be described by Equation (8) in which t = 0 means the end time of the last periodic test.

$$\begin{aligned} F\_{\text{SD}}(t, T\_s \left| F\_{\text{SD}}(0, T\_s) = 0 \right\rangle &= \left. F\_{\text{S}}(t, T\_s \left| F\_{\text{SD}}(0, T\_s) = 0) \right. \\ & \left. + \left. Q\_0 \cdot \left[ 1 - F\_{\text{S}}(t, T\_s \left| F\_{\text{SD}}(0, T\_s) = 0) \right] \right] \right. \right. \\ & \left. = 1 - \exp \left( - \int\_{t\_1}^{T\_s + t\_1} \lambda\_s \left'(u) du \right) \right. \\ & \left. + \left. Q\_0 \cdot \exp \left( - \int\_{t\_1}^{T\_s + t\_1} \lambda\_s \left'(u) du \right) \right. \right. \right. \end{aligned}$$

Where λs(t) is the standby failure rate. T<sup>s</sup> is the standby time interval from the last period test to the instant when IE occurs. t<sup>1</sup> is the time when generating the demand of component after IE occurs.

For other non-monitored components, the failure probability is

$$\begin{aligned} F\_{\text{SD}}(t, T\_s) &= F\_{\text{S}}(t, T\_s) + Q\_0 \cdot [1 - F\_{\text{S}}(t, T\_s)] \\ &= 1 - \exp\left(-\int\_0^{T\_s + t\_1} \lambda\_s^{'}(u) du\right) \\ &+ Q\_0 \cdot \exp\left(-\int\_0^{T\_s + t\_1} \lambda\_s^{'}(u) du\right) \end{aligned} \tag{9}$$

To be conservative, T<sup>s</sup> is chosen to be test interval (TI). So Equations (8, 9) can be written as

$$\begin{aligned} F\_{SD}\{t, T\_s \left| F\_{SD}(0, T\_s) = 0 \} \right) &= 1 - \exp\left(-\int\_{T\_l}^{T\_l + t\_l} \lambda\_s \, ^{'}(u) du\right) \\ &+ Q\_0 \cdot \exp\left(-\int\_{T\_l}^{T\_l + t\_l} \lambda\_s \, ^{'}(u) du\right) \end{aligned} \tag{12}$$

(10)

$$F\_{SD}(t, T\_s) = 1 - \exp\left(-\int\_0^{T l + t\_1} \lambda\_s^{'}(u) du\right)$$

$$+ Q\_0 \cdot \exp\left(-\int\_0^{T l + t\_1} \lambda\_s^{'}(u) du\right) \tag{11}$$

The operational failure means the failure of a component whose functions should be maintained. It includes:


For an irreparable component, the time-dependent probability model related to a certain DET branch is listed in **Table 2**. It illustrates the updating process of logic value and probability when state transition occurs in DET branches.

### The Temporal Failure

The temporal failure relationship is dependent on the timing and order of events. It is assumed that all the failures occur instantaneously. If X and Y are both events which lead to system failure, the temporal failure relationship of two events X and Y can be summarized as:


Obviously, the static fault tree cannot directly describe all of temporal relationship with AND or OR Gates. The issue of temporal failure is focused on standby redundant systems, basically sequential failures. There are two common kinds of sequential failures (Zhang et al., 2018):

• When one of the redundant units/trains fails, standby unit/train actuates and continues running. The sequencedependent failure usually occurs between control/monitoring units and redundant units in a standby train. The failure order of two units determines whether the system fails or not, represented by PAND Gate.


*It is assumed that the failures of components are independent of each other.*

*tX*, *t<sup>Y</sup> are the failure moments of X and Y, tSX , tSY are the demand moments of X and Y. x<sup>i</sup> is the sequence value of the ith event. f<sup>X</sup>* (*t*)*, f<sup>Y</sup>* (*t*) *are the probability distribution functions of X and Y. F*(·) *is the function of failure probability. T is the required operating time of system, also named as mission time.*

• For some standby redundant units/trains, failure of the first operating unit is a prerequisite of the second operating unit to fail. It is called condition-dependent failure, represented by CAND Gate.

In accident analysis of NPP, the temporal failure relations are characterized by supporting systems, not by safety systems, because the actuating conditions of different trains in a safety system are the same, as well as the demanding requirements. In general, all the available trains of safety systems are demanded to operate after actuation. So, switching among redundant trains is considered for supporting systems. For instance, the sequential FT of CCS is shown in **Figure 3**.

Under some circumstances, the order of events is more important than the exact moments when they occur. The sequence value, like in Pandora (Walker and Papadopoulos, 2009), is used to establish a set of temporal laws, and analyze qualitatively to obtain Minimal Cut Sequences (MCSQ). This process can identify some sequential contradictions which reduces the complexity of hybrid model with dynamic and static gates. **Table 3** represents sequence values, temporal relations, and failure probability for logical gates of dynamic fault tree and static fault tree.

### Procedures of Integration Method

In static FT, events are represented by Boolean variables with logical values (True, False, and Normal) and probability value. FT logic gates are treated as operators. In order to ensure that the timing relationship of different variables is not hindered by FT model in accident scenarios, this paper proposes "missionbased DDET" as a framework in consideration of timing and probability characteristics. "Mission phase" is determined according to different functional requirements of safety functions and triggered by branching conditions. In other words, a mission phase will not change until a new event occurs or ending conditions is reached.

To illustrate the modeling and updating process of missionbased DDET, it is necessary to define the state of each branch. Each branch corresponds to a specific state of safety systems, or human actions. The response of safety functions could be represented by either discrete or continuous variables.

1) For the discrete stochastic variables, the number of possible system states is finite, so the DDET simulation runs to implement is finite. For example, the demand states might be success on demand and failure on demand. Besides, it is common practice to spawn branches with the consideration

of how many available trains successfully start-up on a demand, such as the number of ECCS available trains (0/3, 1/3, 2/3, 3/3). Thus, the DDET allows non-binary branching nodes and a specific FT is constructed for a branch state.

2) For the continuous stochastic variables, it is characterized with a continuous probability distribution, such as the response time of operator actions, failure time of components. In DDET method, it is assumed that the response time of continuous safety function is described as a known cumulative distribution function (CDF(t)). Then according to a user-defined percentile/probability threshold, CDF(t) could be discretized into several intervals. The occurrence probability of a representative point for each interval is regarded as branch probability.

TABLE 4 | Nominal parameters for power operation mode.


TABLE 5 | Accident sequence of SB-LOCA in FSAR report (Cd = 0.4).


The progression of Mission-based DDET is shown in **Figure 4**. The skeleton framework of system failures is developed by DDET method while the failures of subsystem and components are represented by FTs. Each node of DDET is labeled by a specific branching event. Therefore, given previous nodes, the conditional branch probability could be calculated by the analysis of FTs. However, this is different from the traditional PRA method which predetermines the success criteria of systems offline. The mission-based FTs are used to calculate the probability of a specific state in a time interval, under the specific performance level. They are mostly used to model operating failures and obtain "the conditional probability of specific state or specific action from the last DET node to this DET node," such as the conditional branch probability of S2 system between Node c and Node d, given only one available train. Meanwhile, the timing and order of standby redundant trains switching is considered in FT modeling process on the basis of the temporal failure relationship as described in section The Temporal Failure. DDET evolves with time according to the branching rules described in section System configuration changes due to branching rules. As a result, the consequences of accident scenarios are determined by DDET simulation results such as safety parameters of interest and the probability of accident scenario is determined by the conditional branch probability.

The procedures of DDET simulation coupled FT is shown in **Figure 5**. After the initialization of DDET simulation, the calculation of simulation model is performed at each discrete time step. The response of process variables, safety components and operator actions are examined, in order to judge whether

branch conditions occurs. If one of the branch conditions is satisfied, the process parameters and system configuration of this node are recorded. After that, time-dependent probability updating is performed as described in section Time-Dependency on Conditional Branch Probability and new branches are triggered. For the sake of simplicity, only binary branch is generated in this paper. Then the system configuration updating of failure branches is carried out. On the one hand, the FT of failure branch is updated as described in section Time-Dependency on Conditional Branch Probability. On the other hand, process variables and control variables are updated in the simulation model, as required by mapping rules of functional failure simulation in section System Configuration Changes Due to Branching Rules. But configuration updating is not required for successful branches. After this updating process is completed, each new branch will restart the simulation process, and continues the above procedures until any terminating conditions of simulation is reached. Finally, simulation results are output which include time response of process parameters and conditional branch probability.

In addition, it is recommended to regard the timing and state of each event as another attribute of logical gates, so that the state of upper level can be estimated in advance under some circumstances. To a large extent, it can prejudge some of system responses in this way, thereby simplifying the system simulation model.

### CASE STUDY

A simple case of a typical PWR with two loops is used for demonstration in this section. The initiating event is Large Break of Loss-of-Coolant Accident (LB-LOCA) on cold leg. The range of break size is from 400 cm<sup>2</sup> to double-ended guillotine break. In the earlier period, the reactor suffers a sub-cooled blowdown due to the sudden break on cold leg. The pressure of primary system drops dramatically from 15.5 to 12.74 MPa, which triggers the reactor to trip.

When the primary pressure falls to 13.0 MPa, ECCS system is actuated. To complement the large amount of coolant inventory, ECCS system including HPI, LPI, and ACC. ACC automatically actuates at the beginning of the transient accident. When the reactor coolant system pressure drops below 4.91 MPa, the N<sup>2</sup> in ACC automatically injects boron water into the reactor coolant system under the pressure drop of N2 inside. LPI keeps injecting the cooling water with boron from the Refueling Tank to the reactor vessel when the primary system pressure drops below 0.98 MPa, until the low level (2.26 m) of Refueling Tank is reached. The pressure and temperature of containment goes up after the transient accident occurs. When the pressure of containment increases to the setpoint value of 0.149 MPa, then Reactor Containment Spray (RCSS) System actuates to cool down the containment. The large amount of coolant leaked from the primary pipe and the cooling water injected are all collected in the bottom of containment. In order to continue cooling down the reactor, the water source of LPI and RCSS changes to the sump, noted as LPI recirculation mode and RCSS recirculation mode.

**Table 4** is the steady state parameters of NPP (Ouyang, 2000). **Table 5** lists the timeline of accident progression in Final Safety Analysis Report. Some of the possible scenarios are illustrated in the ET model, as shown in **Figure 6**, but it cannot cover all the possible scenarios. The progression of accident was modeled with RELAP5/Mod3.1 program. In order to demonstrate the proposed integration method discussed in section Integration Method of FT into DET, a set of phase-based FTs are coupled with the RELAP5 PWR model.

The LPI is taken as an example of conditional branch probability. LPI system is composed of two parallel trains Train A, Train B. Each train consists of one LPI pump and related valves, sensors, etc. In order to compensate the coolant inventory and remove the decay heat from the reactor, the main function of LPI in direct injection mode is to pump the cooling water from the Refueling Tank to the reactor vessel during the first injection period. So, the performance indicator of LPI is mass flow of cooling water.

Assume that:


For the direct injection phase of LPI, the nominal mass flow of each LPI pump is M<sup>f</sup> kg/h, and then the LPI state of direct injection phase is divided by the availability state of Train A, B. The combinations of available trains are various, as listed in **Table 6**. In fact, the earlier timing of LPI failure would lead to a decreased total mass flow of LPI. That may result in variation of state duration, which will influence the subsequent branching timing and conditions of LPR, RCSS.

For State 1, it requires Train A and B both work at nominal state, so FTs of State 1 branch (**Figure 7**) is easily to understand. Assume that the time when actuating LPI is actuated at time t1, and continuously operates until time t2, then the conditional branch probability of functional state 1 is defined as "the LPI system maintains to cool down the reactor from Refueling Tank during the period of t<sup>1</sup> to t2, with the required mass flow of 2 <sup>∗</sup>M<sup>f</sup> ." Thus, the FTs of State 1 is constructed with LHIA0100 and LHIA0200 via OR gate. The conditional probability of State 1 branch is

$$\begin{aligned} P\{M\_{\rm 1}(t)\} &= 1 - F(LHIA0000) \\ &= P\{M\_{\rm 1}(t) = M\_{\rm f} \bigcup M\_{\rm B1}(t) = M\_{\rm f}\} \end{aligned} \tag{12}$$

Where MA(t), MB(t), Ms(t) are the characteristic parameters of Train A, Train B of LPI system.

P(Ms1(t)) is the conditional probability of State 1 branch. To calculate the conditional probability, a set of reliability parameters (**Table 7**) is adopted in FTs and others like LHIA 2000, LHIA 2100, LHIA 2200, and LHIA 2300 is treated with transfer gates, whose probabilities are treated with the consideration of CCF failures.

The FT of State 2 branch can be automatically reconstructed with LHIA0100, LHIA0200 via NOR gate after time-dependent updating. Similarly, the FT of State 3 branch is reconstructed with the same two gates (LHIA0100, LHIA0200) via AND gate. For the consideration of a conservative simulation output, the system configuration of State 2 is modeled that only one train successfully operation, while the other one is not actuated. Similarly, the system configuration of State 3 is modeled as no mass flow.

$$M\_{A2}(t) = M\_f, \; M\_{R2}(t) \; = \begin{cases} M\_f, t\_1 & < t \le \iota \\ 0, \tau < t < t\_3 \end{cases} \text{ or }$$

$$M\_{R2}(t) = M\_f, \; M\_{A2}(t) \; = \begin{cases} M\_f, t\_1 & < t \le \iota \\ 0, \tau < t < t\_3 \end{cases} \tag{13}$$

$$P(M\_{\ell 2}(t)) = P\{M\_{A2}(t) \bigcup M\_{B2}(t)\}\tag{14}$$

$$M\_{A3}(t) = 0, \ M\_{B3}(t) = 0, \ or$$


#### TABLE 7 | Parameters adopted in FTs for LPI system.


#### TABLE 8 | Conditional probability of LPI branches.


$$M\_{A3}(t) = \begin{cases} \quad M\_f, t\_1 < t \le \iota \\ \quad 0, \mathfrak{r} < t < t\_4 \end{cases},$$

$$M\_{R3}(t) = \begin{cases} \quad M\_f, t\_1 < t \le \iota \\ \quad 0, \mathfrak{r} < t < t\_4 \end{cases} \tag{15}$$

$$P(M\_{\mathbb{S}^3}(t)) = P\{M\_{\mathbb{A}^3}(t) \bigcup M\_{\mathbb{B}^3}(t)\}\tag{16}$$

The final results of each LPI branch are listed as **Table 8**, where T is the state duration obtained from the simulation outputs of a certain accident scenario. Compared to the classical PRA results, the mission time is conservatively chosen to be 24 h of all time-dependent events of LPI in accident scenarios, but the state duration of the integration method are based on simulation results, which varies from each accident scenarios.

### CONCLUSION

In Dynamic PRA, the timing and order of events are explicitly considered. But classical FTA/ETA are based on Boolean logic structures with predetermined event sequences and success criteria of systems. This paper aims to incorporate the classical PRA models into simulation-based PRA in a consistent manner. The mathematical basis of integration proves that a timedependent and condition-dependent model is necessary to describe the relationship among system configurations and state duration, influenced by system response of process variables. To better couple the classical PRA with DPRA, an integration method of FT into DDET is investigated. It analyzes about the treatment of dependencies accounted for DET branch probabilities. To ensure the safety level of NPP, the integration method proposes a mission-based DDET framework to describe the time-dependent interactions between physical phenomena, equipment failures, control logic, and operator actions. Mission-based DDET considers both time discretization and discretization of state transition process based on functional demand. It spawns different branches triggered by branching rules and each FT is modeled for a specific system state. As the operating conditions and state duration are determined by the output of simulation model, the conditional branch probability of DDET in the integration method is calculated online by the automatically reconstructed and updated FTs, according to the time-dependent updating rules and system configuration updating rules. A case study of LPI system in LBLOCA accident is taken to demonstrate the feasibility of this integration method.

For the benefits of this paper, the integrated method can provide a means to identify and characterize a priori unknown vulnerable accident scenarios in the safety analysis of NPP, instead of only binary and predetermined logic by analysts. Besides, it reduces the reliance on expert judgment and simplifying (or overly conservative) assumptions about interdependencies. Furthermore, it provides a way to partially reduce the difficulties of complicated modeling and calculation in DET simulation, which gives support for decision-making in safety margin analysis, modifications of operational procedures, changes of system design.

To improve this study, the future work will focus on how to maintain the accuracy and efficiency of integrated framework, how to balance the time step of simulation process and updating frequency of branch probability, how to enhance the capability of constructing FT automatically in a computational tool. In addition, it is prospective to incorporate the performance-based reliability analysis into this integration method.

### DATA AVAILABILITY

All datasets generated for this study are included in the manuscript.

### AUTHOR CONTRIBUTIONS

ZZ instructed and pointed out the problem and proposed to do this research. AX was in charge of and implemented the research to find the reason and solution of the problem. MZ assisted in the case study. HW, HZ, and SC helped to improve and verify the method.

### FUNDING

This study was supported by National Key R&D Program of China on Risk-informed Safety Margin Characterization Technology (2018YFB1900301).

### ACKNOWLEDGMENTS

The authors also appreciate the dedication from Dr. Xiaomeng Dong.

### REFERENCES


**Conflict of Interest Statement:** MZ was employed by China Nuclear Power Engineering Co., Ltd.

The remaining 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 © 2019 Xu, Zhang, Zhang, Wang, Zhang 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.

# Bubble Formation Characteristic of Submerged Single-Hole Orifice in Aerosol Suspension

### Qingyang Sun, Haifeng Gu\*, Xiang Yu and Weikai Yin

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

Bubble columns are regarded as an effective containment filtration method. A large amount of radioactive aerosol will be retained in a liquid pool and generate decay heat. The surface tension, viscosity and other parameters of the liquid phase are affected by increasing temperature of the liquid pool and the presence of aerosol in the pool, which in turn affect the formation characteristics of bubbles. Another important factor affecting bubble generation characteristics is the proportion of steam in the gas phase, that is, the exchange of mass, and energy of gas phase. This paper will investigate the effects of liquid temperature, aerosol and gas phase composition on bubble formation characteristics by analyzing the variation of bubble detachment volume of the orifice in liquid pool.

#### Edited by:

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

#### Reviewed by:

*Yen-Shu Chen, Institute of Nuclear Energy Research (INER), Taiwan Zeyong Wang, Converge Science, United States*

> \*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: *31 May 2019* Accepted: *19 August 2019* Published: *03 September 2019*

#### Citation:

*Sun Q, Gu H, Yu X and Yin W (2019) Bubble Formation Characteristic of Submerged Single-Hole Orifice in Aerosol Suspension. Front. Energy Res. 7:92. doi: 10.3389/fenrg.2019.00092* Keywords: bubble column, containment filtration method, single-hole orifice, aerosol suspension, formation characteristics

### INTRODUCTION

The Containment Filtration and Exhaust System is one important measure to guarantee the integrity of the containment. Bubble columns are regarded as a kind of wet filtration method. The gas in the containment is bubbled through the orifice into the bath, and the radioactive aerosol, or inorganic salt is retained in the bath. Radioactive inorganic salts such as iodine can be reacted with sodium thiosulfate in the liquid pool (Wen et al., 2017), and the aerosol is directly retained in the pool by means of mass transfer sedimentation. Bubbles are the form of filtering gas into the liquid pool, and the bubble detachment volume will directly determine the area of gas-liquid contact and affect the efficiency of filtration during the filtering process. The bubble generation characteristics determine the detachment volume of bubble, so the exploration of that is of great significance in evaluation of the filtration system effect.

The bubbling of the orifice plate is a relatively complicated process, and various forces act on the bubble during formation, including its own momentum, surface tension, drag during rising and buoyancy (Gaddis and Vogelpohl, 1986). One of the leading factors affecting the force of the bubble is the flow of the bubble. Most researchers (Davidson and Schueler, 1960; Davidson and Harrison, 1963) have suggested that the decisive factors affecting the size of bubble formation are surface tension and liquid phase density at low flow rates. At high flow rates, the flow will be the main factor. A variety of bubble-off volume prediction models have also been given by different researchers (Davidson and Schueler, 1960; Davidson and Harrison, 1963; Kumar and Kuloor, 1970) which use flow as a parameter of interest, with a simple bubble volume formula:

$$V\_b = K \frac{Q\_G^{6/5}}{\mathcal{g}^{3/5}} \tag{1}$$

K is the empirical coefficient obtained by the experiment. Davidson proposed two different coefficients of 1.378 (Davidson and Schueler, 1960) and 1.138 (Davidson and Harrison, 1963), and K = 0.976 was obtained in the experiment of Kumar (Kumar and Kuloor, 1970). It is the ignoring of the influence of surface tension, viscosity and orifice diameter in formula (1) that leads to diverse empirical results in different experiments. Another bubble diameter calculation model that combines more factors was proposed by Gaddis (Bhuiyan et al., 2015):

$$d\_b = \left[ \left( \frac{6d\_o \sigma}{\rho\_l g} \right)^{\frac{4}{3}} + \left( \frac{81\nu Q}{\pi g} \right) + \left( \frac{135Q^2}{4\pi^2 g} \right)^{\frac{4}{5}} \right]^{\frac{1}{4}} \tag{2}$$

The influence of liquid phase properties on the bubble generation characteristics has always been the focus of scholars. The addition of inorganic or solid particles in the liquid pool changes the liquid phase properties. For example, when the liquid phase is a sodium thiosulfate solution, as the concentration of the solution increases, the dimension of generated bubbles decreases (Wen et al., 2017), and the presence of inorganic salts also inhibits bubble fusion (Kracht and Finch, 2009). The addition of solid particles in the liquid phase is believed to increase the load of bubble rise and reduce the rate of rise, and the solid particles also promote the fusion of bubbles (Vazirizadeh et al., 2016). There have been many studies on the nature of the addition of inorganic salts in the liquid phase, but studies on aerosol-based solid particles have been fewer as more have focused on the characteristics of the bubble rise process, and there is not enough investigation into bubble formation characteristics. Therefore, there is even less attention paid to liquid phase temperature and gas phase properties during bubbling. Based on this, the present paper focuses on the effects of liquid phase addition of aerosol particles, high temperature liquid phase, gas phase temperature, and composition on the formation characteristics of single bubbles.

### EXPERIMENTAL METHOD

The schematic of the test facility is shown in **Figure 1**. The overall experimental circuit consists of four parts for air supply, steam generation, tracheal transport, bubble generation and imaging, and is also equipped with a parameter measurement, and acquisition device. A rectangular open water tank with a size of 300 × 200 × 700 is provided at the bottom of the bubble generating portion of the experimental table. The facility is designed to run at atmosphere pressure and temperature ranging from 25 to 98◦C. The experimental gas flow rate was measured by a mass flow meter, and the flow rate was controlled from 0.05 to 0.20 l/min.

The bubble image is captured by a high-speed camera, and the bubble is processed by Matlab software to extract its characteristic parameters. Pre-processing of the captured image is required because of the large amount of noise in the image taken directly from the experiment. A clearer bubble profile will be obtained after image background removal, median filtering and binarization.

The program calculation relies on the processed image contour to obtain the most suitable center point and symmetry axis of the bubble, and then adopts the processing method of the slice to cut the bubble into a plurality of small slices in a direction perpendicular to the axis of symmetry, each of which can be regarded as a cylinder. Under the idea of differentiation and integration, the volume of the entire bubble is obtained by accumulating the volume of the cylinder. The calculation formula is as shown in equation (3).

$$V\_{\mathbf{b}} = \sum\_{\mathbf{i}=1}^{n} \frac{\pi \,\mathrm{d}\_{\mathbf{i}} \,^2 \mathrm{h}\_{\mathbf{i}}}{4} \tag{3}$$

The result obtained is the pixel volume of the bubble in the picture, an actual scale k will be given for each image in the

experimental shooting and finally volume represented by the pixel can be converted into the actual bubble volume:

$$V\_{\mathbf{b}}{}^{\prime} = \mathbf{k}^{3} V\_{\mathbf{b}} \tag{4}$$

### RESULTS AND DISCUSSION

### Effect of Orifice Diameter, Gas Phase Flow on Bubble Formation

The influence of orifice diameter, gas phase flow and liquid phase properties during bubbling has been a concern. In previous studies, it was found that the momentum of the gas phase and other factors change with the increase of pore size and gas phase flow, which in turn makes the size of the bubble increase (Badam et al., 2007). In this paper, as shown in the results of the deionized water experiment in **Figure 2**, the bubble increased significantly as the orifice diameter increased from 1 to 2 mm and the gas flow rate increased from 0.05l/min to 0.2l/min. On the other hand, for the study of liquid phase properties, the researchers found that bubbling is not affected by liquid level changes in the range of liquid level studied (600–2,100 mm) (Jamialahmadi et al., 2001). In this paper, the aerosol in the liquid phase as well as the liquid phase and gas phase temperatures were taken as the research object. And to ensure the proper experimental variables, diameters of 1 mm for the orifice plate and of 700 mm for the liquid level were uniformly used in the study.

### Effect of Aerosol Particles on Bubble Formation

Polydisperse aerosol particles TiO<sup>2</sup> with a median particle size of 300 nm were chosen in this experiment. To ensure the visibility of the experimental suspension, two concentrations of 0.004 g/l and 0.008 g/l were selected for the experiment. The variation of bubble detachment volume with flow rate in different concentrations of TiO<sup>2</sup> suspension at 25◦C is shown in **Figure 2**. As a result, the addition of TiO<sup>2</sup> particles in the liquid causes a larger detachment volume of the bubbles compared to the

deionized water, and the size increases as the concentration of the TiO<sup>2</sup> particles increases. Nanofluids such as TiO<sup>2</sup> have been found in previous studies to increase surface tension (Bhuiyan et al., 2015). Therefore, the surface tension was measured by the platinum plate method in the experiment, and the surface tension of the TiO<sup>2</sup> suspension of 0.008 g/l at 25◦C (72.2 mN/m) was slightly larger than that of the deionized water (71.9 mN/m). When the bubble is detached and the force is balanced, the surface tension is the resistance, which is directly affected by the liquid phase property when the bubble contact angle is zero (Gaddis and Vogelpohl, 1986). And surface tension will affect the buoyancy and detachment size of the bubble when other external conditions such as flow rate are constant.

### Effect of Gas-Liquid Phase Temperature on Bubble Formation

The water tank and the vent pipe were insulated in the experiment to maintain the temperature of the gas phase, as shown in **Figure 3**, the bubble size at different liquid temperatures when the gas was 25◦C and flow rate was 0.2 l/min. The results show that the bubble size increases after the liquid temperature rises, and its increase is more severe as the temperature approaches saturation.

The nature of the bubble must be affected by the nature of the liquid or gas phase. On the one hand, in the research and analysis of the changes in liquid phase properties, previous researchers have shown that the increasing of liquid temperature will greatly decrease the surface tension at the gas-liquid interface, which in turn will reduce the size of the bubble. However, it can be seen from the experimental results that the high-temperature liquid phase increases the size of the bubble. Therefore, it is suggested that most scholars only pay attention to the change in liquid phase properties but ignore gas phase. If the gas phase expands or the liquid phase evaporates into the gas phase during bubble generation, the result can be explained. Among them, the expansion will be derived from the temperature difference between the gas and the liquid, so the effect on the gas phase temperature was verified in the experiment. The bubble volume comparison at different gas phase temperatures at a flow rate of 0.2 l/min and a liquid temperature of 60◦C is shown in **Figure 3B**. It can be seen from the figure that the bubble generation size does not change significantly under different gas-liquid temperature differences. Therefore, it is considered that the gas phase does not undergo a significant expansion phenomenon, and it is judged that the increase in the bubble size is derived from the evaporation of the liquid phase.

### Effect of Gas Phase Composition on Bubble Formation

Since the gas source for bubbling in the experiment comes from the gas storage tank and the boiler, respectively, the bubbling characteristics when the gas phase in the untreated air or the steam-containing air can be analyzed in the experiment. As

**Figure 4B** shows the shape of the bubble when eliminating the influence of liquid phase evaporation and the effect of liquid

liquid temperature.

phase evaporation is eliminated, and **Figure 4C** shows the form case of the liquid phase for low temperature bubbles formed when the liquid phase is cold. It can be seen from the comparison that the difference in the morphology of the bubble formed under these two conditions is not obvious. Therefore, it can be considered that the influence of the liquid phase temperature rise is mainly manifested by the evaporation of the liquid phase. When the influence of the liquid phase evaporation is removed, the influence of the liquid phase temperature change is relatively weak.

### Establishment of Bubble Generation Volume Model

Some prediction models for the volume of orifice bubbling are given in the preface, and the results obtained in this experiment are compared with the existing formula model, as shown in **Figure 5**.

In the model comparison, the experimental values at the liquid temperature of 25◦C are better fitted with some models, but when the temperature is higher, the results have a larger deviation since the model does not consider the evaporation factor of the liquid phase.

In view of the result fitting, the calculation formula of Gaddis and Vogelpohl (1986) with better degree of fitting is selected for correction. Inspired by the experimental results of the introduction of steam-containing gas phase, the flow of untreated air can be converted into a wet-saturated gas flow equivalent to the air partial flow through the gas state equation, taking into account the liquid phase portion of the evaporation. The modified model can be obtained by substituting the converted flow value into the original model:

$$Q' = \frac{PQT'}{TP'} \tag{5}$$

$$V\_b = \frac{\pi}{6} \ast \left[ \left( \frac{6d\_o \sigma}{\rho\_l g} \right)^{\frac{4}{3}} + \left( \frac{81\nu}{\pi g} \frac{PQT'}{TP'} \right) + \left( \frac{135}{4\pi^2 g} \left( \frac{PQT'}{TP'} \right)^2 \right)^{\frac{4}{3}} \right]^{\frac{3}{4}} (6)$$

Where Q, T, and P are the flow rate, temperature and pressure of untreated air, Q′ , T′ , and P′ are the flow rate, temperature and pressure of wet-saturated air at the same temperature, and d<sup>o</sup> is the diameter of the orifice. Equation 6 is obtained by the cubic equation of the bubble size after forcing analysis by Gaddis, and each set of variables in the formula represents a set of coefficients.

The experimental results of different working conditions were fitted with the corrected bubble volume calculation model—the fitting error was within ±30%, which is shown in **Figure 6**.

### CONCLUSION


### DATA AVAILABILITY

The datasets generated for this study are available on request to the corresponding author.

### AUTHOR CONTRIBUTIONS

QS conducted the experiments, data processing, and subsequent analysis. HG guided the direction of the experiment and the writing of the paper. XY is an experimental collaborator and coauthors of the article. WY is an experimental collaborator and collects the experiment data.

### 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 © 2019 Sun, Gu, Yu and Yin. 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.

# Preliminary Experimental Investigation on the Filtration Performance of Submicron Insoluble Aerosol in a Bubble Column

Yingzhi Li <sup>1</sup> \*, Kaiyi Shi <sup>2</sup> , Zhongning Sun<sup>1</sup> , Haifeng Gu<sup>1</sup> and Yanmin Zhou<sup>1</sup> \*

*<sup>1</sup> Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin, China, <sup>2</sup> School of Mining and Civil Engineering, Liupanshui Normal University, Liupanshui, China*

Pool scrubbing is a representative process for aerosol retention involved in severe LWR accidents and has proved to be a potential method for the removal of micron-sized aerosols. Bubble scrubbing, as the key component of pool scrubbing, has the advantages of a large gas-liquid contact area and long residence time. However, the filtration mechanism of submicron aerosol by bubble scrubbing is still not clear. Therefore, the goal of this research is to evaluate the filtration performance of submicron insoluble aerosol in a bubble column. Some preliminary experiments were carried out to investigate the influence of operational parameters on aerosol filtration efficiency, such as inlet aerosol concentration, gas flow rate and pool temperature. Mass method was used to analyze the aerosol filtration efficiency. The results show that the aerosol filtration efficiency is independent of inlet aerosol concentration when inlet aerosol concentration is in the range of 0.2–0.5 g/m<sup>3</sup> . Besides the aerosol filtration, efficiency is influenced by the superficial gas velocity. With the increase of superficial gas velocity, homogeneous bubbling, transition and a heterogeneous bubbling regime will appear in turn. The aerosol filtration efficiency decreases in the transition regime due to the appearance of large bubbles but increases in the heterogeneous bubbling regime account for intense gas-liquid interaction. Moreover, the aerosol filtration efficiency significantly decreases with pool temperature because the evaporation rate at the gas-liquid interface will significantly increase as the boiling point is approached, which gives retroaction to aerosol filtration efficiency.

Keywords: insoluble aerosol, operation parameter, filtration efficiency, bubble column, pool scrubbing

## INTRODUCTION

Since the Fukushima nuclear accident occurred, more attention has been paid to the issue of severe accidents. To prevent radioactive material from leaking into the environment, some dry or wet filtration methods are used. Pool scrubbing is a representative process for aerosol retention involved in severe LWR accidents. Examples include the core reflooding process in PWR, gas injection into suppression pool in BWR and also the wet scrubbing process in Filter Containment Venting Systems (FCVS) (Basu, 2014). The behavior of pool scrubbing can be divided into 4 parts, which include gas inject region, globule-breakup region, bubble swarm

#### Edited by:

*Wenxi Tian, Xi'an Jiaotong University, China*

#### Reviewed by:

*Jiejin Cai, South China University of Technology, China Shimpei Saito, Asahi Glass, Japan*

#### \*Correspondence:

*Yingzhi Li lyz1993@hrbeu.edu.cn Yanmin Zhou zhouyanmin@hrbeu.edu.cn*

#### Specialty section:

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

Received: *03 July 2019* Accepted: *26 August 2019* Published: *11 September 2019*

#### Citation:

*Li Y, Shi K, Sun Z, Gu H and Zhou Y (2019) Preliminary Experimental Investigation on the Filtration Performance of Submicron Insoluble Aerosol in a Bubble Column. Front. Energy Res. 7:96. doi: 10.3389/fenrg.2019.00096* region and entrainment droplet region (Owczarski and Burk, 1991). The corresponding filtration mechanism of each region is different. The bubble scrubbing process is especially considered as the migration process of particles inside the bubble toward the liquid phase.

The deposition characteristics of aerosol in the rising bubble involve various mechanisms including gravitational settlement, Brownian diffusion, inertial impaction, and they also include thermophoresis and diffusiophoresis in the case of thermal nonequilibrium conditions. The calculation formula of deposition efficiency of particles in spherical bubbles was first proposed by Fuch (1964), which takes into account the gravitational settlement, Brownian diffusion, and inertial impaction. Additionally, the influence of each deposition mechanism is considered independent. The change of particle concentration in the bubble is equal to the algebraic superposition of concentration variation caused by each mechanism. From the model of Fuchs, the removal efficiency of particles due to Brownian diffusion decreases with the increase of particle diameter, while the removal efficiency due to gravitational sedimentation and inertial impaction increases with the increase of particle diameter. Because each removal mechanism has an effective size range, the particle size that has a minimum efficiency is in the nanometer size range, which is also called the most penetrating particle size (MPPS).

Experimental investigations on aerosol filtration performance by pool scrubbing have been conducted since the 1980s, and typical experiments include the Advanced Containment Experiments (ACE), the Battel Columbus Laboratories experiment (BCL), the Pool Atomic Energy Research Institute programme (JAERI), and the Pool Scrubbing Effect on Iodine Decontamination programme (POSEIDON), among others (Hakii et al., 1990; Kaneko et al., 1992; Escudero Berzal et al., 1995; Peyres and Herranz, 1996). Later, the POSEIDON-II experiment was carried out by Dehbi et al. (2001), which included 17 groups of different operating conditions with the aim of revealing the impact of various operation parameters on the decontamination process. The results showed that DF increased with the increase of submergence depth, steam fraction, and aerosol particle size. The author also studied the influence of several removal mechanisms by controlling the single change of entrance parameters, for example the effect of an inertial impaction mechanism can be verified by a single change of particle size. Kanai et al. (2016) studied the relationship between aerosol concentration and submerged depth by using an aerosol spectrometer, and a simple formula was proposed to calculate the decontamination factor (DF). Charvet et al. (2011) and Cadavid-Rodriguez et al. (2014) studied the influence of different operation parameters on nanoparticle removal efficiency, but the experimental gas flow rate was too small, which limits the generalization of experimental results. Some investigations indicate that the particle removal efficiency increases with the inlet particle concentration due to the coagulation of particles (Meikap and Biswas, 2004). However, Sun et al. (2018) experimentally investigated the relationship between DF and aerosol concentration with three measurement methods and drew the conclusion that DF increases with the decrease of aerosol concentration. The same experimental phenomenon also appeared in the study of Bandyopadhyay and Biswas (2006). Recently, Abe et al. (2018) visually studied the bubble dynamics with aerosol and flow structure during pool scrubbing and found that the interface oscillation of globules leads to the amount of aerosol collected by liquid phase.

Previous investigations confirmed pool scrubbing as a potential method to remove the micron-sized aerosol. Nonetheless, the influence of some experimental factors has not been clearly studied—for example, inlet aerosol concentration, gas flow rate, and subcooling degree. Besides, the aerosol removal mechanism is not clear and there is a lack of reliable experiment data to valid predictive models. Therefore, the goal of this paper is to study the effect of some operating parameters on the filtration efficiency of submicron aerosol by bubble scrubbing. Furthermore, this research is also expected to provide theoretical support for the design of pool scrubbing and source term analysis for severe nuclear accidents.

### EXPERIMENTAL FACILITY

**Figure 1** shows the schematic diagram of the aerosol experimental facility. The main test facility can be used for high-temperature and -pressure experiments, and the main test facility performance is listed in **Table 1**. The test section is made of stainless steel with an inner diameter of 260 mm. Additionally, the test vessel is 2.5 m in height. An air compressor and electric boiler provide air flow rate and steam flow rate, which can be adjusted to a desirable value. The air mass flowrate is measured by two mass flowmeters with operating ranges of 0–100 kg/h and a measurement accuracy of 0.2%. A uniform gas flow distribution is generated in the bubble column by a perforated stainless steel plate with 500 orifices, and the orifice diameter is 2 mm. Two heaters with different electric powers were used to heat up the water in the bubble column and maintain a constant pool temperature during the experiment process. Average gas hold-up was measured by a differential pressure transducer with an operating range of 0–20 KPa and a measurement accuracy of 0.055%. Pressure probes were installed at different positions away from the orifice plate and the pool surface to avoid fringe effect.

Insoluble barium sulfate (BaSO4) powder was used as an aerosol simulant material. **Figure 2** gives the aerosol size distribution, which was measured by a welas 2000 laser particle size analyzer. The median diameter of the size distribution is 0.5µm. Aerosol particles with a particle size <0.5µm account for almost 45% of the total volume. The inlet aerosol concentration can be adjusted by switching the supplied gas flow rate through the aerosol generation. A mixed chamber was used to ensure the aerosol concentration was uniform and steady. Additionally, a bypass flow was set to guarantee the inlet aerosol concentration constant if the inlet laden gas flow rate was changed.

The mass method was used to analyze the removal efficiency of aerosol particles by passing sample gas through a glass fiber filter paper mounted on a filter hold. The mass of filter paper was measured by a high precision microbalance with an accuracy of 0.1 mg. No less than two measurements were conducted to

TABLE 1 | The main test facility performance.


reduce experimental error under each design condition. In order to prevent the measurement error of aerosol sampling caused by vapor condensation, the sampling pipeline needs to be heated and insulated, and isokinetic sampling should be guaranteed. The mass difference of the filter paper gave the total mass of aerosol particles collected, and the sample air volume was measured by two mass flowmeters, with operation ranges of 2–50 L/min and a relative deviation of 1%. Steam was removed in the snake condenser before entering the flowmeter when the sample gas contained water vapor, and the condensate water could be collected. Thus, the total aerosol mass removal efficiency was determined as follows:

$$\eta = 1 - \frac{c\_{out}}{c\_{in}} = 1 - \frac{m\_{out}/\nu\_{\text{gout}}}{m\_{in}/\nu\_{\text{gin}}} \tag{1}$$

Where is the mass of unremoved aerosol in outlet air retained on filter paper, g; is the mass of particles in inlet air retained on filter paper, g; is the outlet sample air volume, m<sup>3</sup> ; and is the inlet sample air volume, m<sup>3</sup> .

### RESULTS AND DISCUSSION

### Influence of Inlet Aerosol Concentration on Filtration Efficiency

The inlet aerosol concentration increases with the air flow rate supplied for aerosol generation, which is in the range of 0.006– 0.5 g/m<sup>3</sup> . The test air flow rate is 10 kg/h and the submergence depth is 1.2 m. **Figure 3A** presents the aerosol filtration efficiency against inlet aerosol concentration. It can be seen that the aerosol filtration efficiency is almost equal to 100% when the inlet aerosol concentration is lower than 0.1 g/m3, but the aerosol filtration efficiency is independent of inlet aerosol concentration when inlet aerosol concentration is in the range of 0.2–0.5 g/m<sup>3</sup> . It can be explained that the inlet aerosol concentration is so low that the aerosol coagulation is limited. Because the outlet aerosol concentration is so low, there exists some measurement uncertainty by the mass method when inlet aerosol concentration is lower than 0.1 g/m<sup>3</sup> .

### Influence of Superficial Gas Velocity on Filtration Efficiency

In order to verify the effect of superficial gas velocity on aerosol filtration efficiency, control experiments were conducted at different superficial gas velocities (8–100 kg/h) and at constant submerged depths of 0.9 m and 1.2 m. **Figure 3B** shows the aerosol filtration efficiency achieved at different air flow rates. It can be seen that when superficial gas velocity is >0.12 m/s, the aerosol filtration efficiency increases with the increase of superficial gas velocity and tends to be stable when superficial

temperature on evaporation rate.

gas velocity exceeds 0.2 m/s. Furthermore, the aerosol filtration efficiency with submerged depth of 1.2 m is higher than the aerosol filtration efficiency with submerged depth of 0.9 m at the same gas flow rate, which also indicates that the aerosol removal efficiency increases with the submerged depth.

**Figure 3B** also shows the relationship of fraction gas hold-up and superficial gas velocity. With the increase of gas flow rate, homogeneous bubbling, transition, and heterogeneous bubbling regimes will appear in turn. It shows that the regime is a homogeneous bubbling regime at superficial gas velocities lower than 0.04 m/s. There is a linear relationship between gas holdup and apparent gas velocity. The bubble size distribution is narrow and the gas hold-up is radially uniform. But the heterogeneous bubbling regime occurred when superficial gas velocity was higher than 0.12 m/s. The gas-liquid interaction is strong and the bubble size distribution is wide. These two flow patterns are separated by the transition regime, which is caused by local liquid circulation in the bubble column (Vial et al., 2000).

The two-phase flow pattern is a heterogeneous bubbling regime when the superficial gas velocity exceeds 0.12 m/s. When the superficial gas velocity increases, the fractional gasholdup increases as well as the intensity of interface oscillation and bubble breakup and coalescence. The inertial impaction effect would become a dominant mechanism, which contributed positively to higher particle collection efficiency (Abe et al., 2018). However, when the superficial gas velocity is too large, the coalescence effect of the bubbles was enhanced, which would reduce the interfacial area between gas phase and liquid phase. Thus, the aerosol filtration efficiency tends to stabilize finally. When the superficial gas velocity is lower than 0.12 m/s, the two-phase flow pattern is the transition regime in this experimental study, and the aerosol filtration efficiency decreases with superficial gas velocity. The increase of superficial gas velocity causes increased fractioning of large bubbles, which further leads to the decrease of aerosol filtration efficiency. In this experiment, the two-phase flow pattern did not include a homogeneous bubbling regime, because an orifice plate with 2 mm orifice diameter was used.

### Influence of Pool Temperature on Filtration Efficiency

In order to verify the effect of pool temperature on aerosol filtration efficiency, pure air experiments were conducted at different pool temperatures, a constant submerged depth of 0.7 m, and a gas flow rate of 10 kg/s. **Figure 3C** shows the aerosol filtration efficiency achieved at different pool temperatures. It can be seen that the aerosol filtration efficiency decreases with the increase of pool temperature. The pool temperature increases from 35 to 95◦C, and the aerosol filtration efficiency significantly decreases from 90 to 40%. The aerosol filtration efficiency is somewhat anomalous between 75 and 90◦C, which may require further repetitions of the experiment in order to clarify.

The increasing temperature leads to a decrease in the surface tension coefficient and kinematic viscosity. With the decrease of liquid surface tension with temperature, not only does the bubble detachment volume decrease (Bhavaraju et al., 1978), but the film drainage speed decreases (Lin et al., 2010) and thus leads to a decrease of bubble size. With the decrease of liquid viscosity, the bubble rise velocity increases, which results in a reduction of bubble residence time. However, the interface evaporation of gas-liquid phases leads to an increase of bubble size. Additionally, the gas density decreases with increasing temperature, which also increases the bubble size and decreases gas holdup. Thus, the total influence of pool temperature on bubble size and gas holdup depends on the competitive effect. As reported by Pohorecki et al. (2001) and Feng et al. (2019), the water temperature has a slight influence on bubble size, gas holdup, and interfacial area.

When the pool temperature increases, the subcooling degree decreases while the saturation vapor pressure increases, which influences the evaporation rate at the gas-liquid interface. **Figure 3D** gives the relationship between the evaporation rate and pool temperature. The evaporation rate is defined as the ratio of the mass of water vapor in the outlet gas to the mass of air. The mass of water vapor can be obtained by condensate water and the air mass can be achieved by the mass flowmeter installed at the sample line. It can be seen that the evaporation rate at the gas-liquid interface will significantly increase as the boiling point is approached, which gives retroaction to aerosol filtration efficiency. From the experimental results, it can be obtained that the interfacial evaporation of water vapor, which is called Stefan flow, plays a dominant role in collecting aerosol by bubbly scrubbing under low-subcooling conditions.

## CONCLUSIONS

In this study, the filtration efficiency of insoluble aerosol in a bubble column was experimentally investigated under pure air condition, varying several operation conditions such as inlet aerosol concentration, gas flow rate, and pool temperature. The main conclusions are as follows:


## DATA AVAILABILITY

The datasets generated for this study are available on request to the corresponding author.

## AUTHOR CONTRIBUTIONS

YL and YZ carried out experiments and analyzed the experimental results as well as the writing of the paper. YL and KS assisted in experimental work and dealt with the experimental data. ZS and YZ set the guideline of this research and modified the language of the manuscript. HG and KS designed the experimental loop and guided the experimental work.

### ACKNOWLEDGMENTS

The authors greatly appreciate support of the Natural Science Foundation of China (Grant No. 11675045). The author thanks

### REFERENCES


the 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). Additionally, the authors are thankful for the support of the Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, China.

in 22nd DOE/NRC Nuclear Air Cleaning and Treatment Conference (Denver, CO).


**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 © 2019 Li, Shi, Sun, Gu and Zhou. 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.