Printer Friendly

Multiphysical Simulations for the IAEA/ISCP Benchmark Model on the Contact of Pressure Tube and Calandria Tube in the Moderator System of CANDU-6 PHWR.

1. Introduction

The CANDU (CANadian Deuterium Uranium) reactors have been adopted in Korea since late 1980s [1], operated with four CANDU-6 units generating 2.6 GW overall output in Wolseong nuclear plant. One of the most important safety issues in the CANDU reactors is the moderator subcooling margin [2] during a large loss of coolant accident to ensure the integrity of the fuel channels submerged in the heavy-water moderator.

International Atomic Energy Agency (IAEA) has launched a committee proposing an International Collaborative Standard Problem (ICSP) to provide a new contact boiling experimental data to assess the subcooling requirements for a heated pressure tube plastically deformed to contact with the calandria tube during a postulated large break LOCA condition. Since participating in the first Workshop of IAEA/ICSP on "HWR Moderator Subcooling Requirements to Demonstrate Backup Heat Sink Capabilities of Moderator during Accidents" held in Ottawa, Canada, 2012 [3], the authors have introduced a multiphysical code based on COMSOL [4] to show the solution of blind test, also compared with one-dimensional CATHENA code result [5]. After the workshop, the two-dimensional simulation is evolved to consider the boiling model at the outer calandria tube bounded emerged on the reservoir of coolant in the moderator system.

The plastic deformation ofa pressure tube is brought forth with the high temperature from heat transfer of conduction and radiation where the model of Shewfelt and Godin [6] is applied to the structural dynamics equation of COMSOL. In this paper, the strain rate and the contact time of two tubes are computed for various heating-up conditions of the graphite heater located at the center of apparatus to investigate the sensitivity of open experimental conditions.

2. Numerical Method

2.1. IAEA/ICSP Benchmark. Figure 1(a) shows the IAEA/ ICSP for the benchmark of the LOCA model that consists of a graphite heater and the surrounding pipes called pressure tube (PT) and calandria tube (CT) [7]. As the gap between PT and CT is filled with carbon dioxide for the insulation, the center axis of graphite heater is located eccentric 9.5 mm lower from the centroid of tubes. However, the convection effect is very small during the simulation, so this error is ignored in the numerical model of this study. The thicknesses of tubes are 4.4 mm (PT) and 1.42 mm (CT), and also the inner diameters mark 103.6 mm (PT) and 129 mm (CT), respectively. The electrical heater exerts the heat power measured for time (red line) in Figure 1(b) where the total time is 140 seconds, so this scenario is modeled similarly for COMSOL (ver. 5.3a) simulation (blue line).

2.2. Computational Model. The governing equations are summarized such as (1) heat transfer (conductive for each tube; radiative among heater and tubes) and (2) structural dynamics where the range of thermal stress lies beyond the linear elastic range [5]. Therefore, this problem becomes multiphysics involving the coupling of two kinds of phenomena.

-[nabla] x [??] = [F.sub.V] (1)

[mathematical expression not reproducible] (2)

[??] = 1/2 ([nabla][u.sup.T] + [nabla]u) (3)

[rho][c.sub.p][[partial derivative]T/[partial derivative]t] = [nabla] x (k[nabla]T) + [Q.sub.V] (4)

In the mechanical equilibrium at each time step, (1): the stress tensor ([??]) is conserved for the external volumetric force ([F.sub.V]), which should be zero for this problem except for boundaries; and the difference stress from the initial condition, denoted as subscript 0 in (2), is linearly proportional to the strain vector ([epsilon]) with addition of thermal expansions depending on temperature (T) with thermal expansion coefficient ([??]), which can be expressed as the function of displacement gradient ([nabla]u) in (3). The elastic coefficient for the isotropic materials in (2) is composed of a matrix ([??]) dependent on Young's modulus (E) and Poisson's ratio (v). This problem becomes multiphysical because (2) contains the temperature (T) that is the solution of (4) where [rho], [c.sub.p], and k are material density, heat capacity, and thermal conductivity, respectively. The volumetric heat ([Q.sub.V]) distribution should be prescribed at the graphite heater located at the central axis in Figure 1(a).

Equation (1) is solved for components of the stress tensor, and the external force is given as boundary conditions. Note that the stress tensor, later reduced to the plane stress, should be three-dimensional though the computational domain is two-dimensional, which is due to the Poisson's ratio in the linear elastic stiffness matrix of (2). For isotropic materials, the stiffness is expressed as a symmetric square matrix with the six degrees of freedom [7]:

[mathematical expression not reproducible] (5)

where E is the Young's modulus and v is the Poisson's ratio. Plane strain is hypothesized for the two-dimensional approximation as the tubes in Figure 1(a) are slender enough. The computational domain consists of Zr 2.5Nb (PT) mounted concentrically inside a section of Zircaloy 2 (CT): see Figure 2(a), which are given as a material properties in E and v, functions of local temperature.

The boundary conditions at the inner circles of tubes should be prescribed by the displacement ([u.sub.r]) of the radial expansion from the initial radius [r.sub.0] as

[mathematical expression not reproducible] (6)

where only one-dimensional expansion is applied for the radial direction. Note that the thermal expansion term is modified from (2) where a is not a constant but a function of T. Therefore, four-point Gaussian quadrature is used for the integration in (6) [5].

In Figure 2(a), the boundary is pressurized inside the PT, and the pressure of carbon dioxide gas is almost atmospheric one, between PT and CT filled as a thermal insulator. Since this gas acts as an adiabatic thin layer between tubes, the natural convection is totally ignored in this study, only considering the dominant radiation effect. The volumetric heat ([Q.sub.V]) in (4) is modeled as Figure 1(b), and the thermal radiation at each surface is modeled with Stefan-Boltzmann law such as

[??] x (k[nabla]T) = e(G - [kappa][T.sup.4]) (7)

J = (1 - e)G + e[kappa][T.sup.4] (8)

where [??] is the outward-normal unit vector and [kappa] is a constant of 5.67 x [10.sup.-8] W/([m.sup.2][K.sup.4]). Equation (7) is solved for the local temperature in the boundaries to calculate the surface radiation (J) in (8), and the incident radiation (G) is an explicit function of J and the view factor determined iteratively by the mutual geometric configuration of surfaceto-surface radiation boundaries where e is emissivity marked in Figure 2(a), and their values of the core graphite heater, PT, and CT are 1.0, 0.8, and 0.34, respectively. Recall that Figure 2(a) is symmetric in the circumference direction, but is not simply one-dimensional because we should calculate J and G in (7)-(8) from the two-dimensional configuration of geometry. The outer temperature of water reservoir is fixed as 70[degrees]C, or subcooling of about 30[degrees]C under the standard atmospheric condition. However, after PT contacts with CT, the heat flux abruptly increases to full up the outer local temperature even to the boiling condition, which will be explained in Section 2.4.

Figure 2(b) is the triangular grids for the computation consisting of 3,424 elements with 952 edge and 20 vertex ones.

2.3. Nonlinear Thermal Deformation Model. As shown in (6), the thermal expansion a(T) is not a simple constant at the range of plastic deformation, and this nonlinear behavior is simplified with Shewfelt and Godin creeping model [13]:

[mathematical expression not reproducible] (9)

where the mechanical stress is s = [DELTA]p x r/w ([DELTA]p: pressure difference; r: radius, u>: thickness of the tube) and [t.sub.1] and [t.sub.2] are times at temperature [T.sub.1] = 973K and [T.sub.2] = 1123K, respectively. In (9), the time integration has no analytic solution, so the numerical method of four-point Gaussian quadrature is applied to compute it [5].

The strain ([epsilon]) is related to the strain rate, [??] = d[epsilon]/dt in (9), and a function of temperature from (6):

[mathematical expression not reproducible] (10)


The thermal expansion, the unknown in (10), is plotted as a function of temperature in Figure 3, and it is remained constant at the linear elastic range, but increased exponentially at the nonlinear plastic range or T > [T.sub.1]. The expansion rate is abruptly decreased at T [approximately equal to] [T.sub.2], but increased again as heat is added: definitions of [T.sub.1] and [T.sub.2] are given in (9).

As the temperature rate (dT/di) is computed in the solution of (4), (9) is transformed simply as a function of temperature [5]. Four-point Gaussian quadrature is used for the integration in (9) and (10). The integrated model in (10) is substitute to the right-hand side of (2) to get the stress field data in (1).

Figures 4(a) and 4(b) are the simulation result before and after contact, respectively, and the contact time is 74.6 s. After the expanded PT contacts with CT (see Figure 4(b)), the heat flux at the outer boundary in Figure 2(b) starts to be increased with natural convection in the water reservoir. However, the heat transfer soon transits to nucleate boiling after the outer temperature of CT exceeding the boiling point. It can be hardly observed that the boiling reaches to the critical heat flux. Therefore, the heat resistance between the contacted PT and CT due to the heat conduction at a very thin boundary of thickness [delta], or the heat conduction gradient is initially assumed as k/[delta] [approximately equal to] 1 kW/([m.sup.2]K), which is estimated from the experimental data of IAEA/ISCP. Because this value was not sufficient to bring about the critical heat flux, we numerically controlled this unknown value within a very restricted time of 0.1 seconds, for example, during the contact. The physical fact relying on this numerical correction is that the contact is initially an incomplete one, and highly amplified as the gap thickness ([delta]) reaches to a very small value in a local "spots", emitting very high heat flux to superheat the heavy water.

2.4. Boiling Heat Flux Model. When the temperature arise near 100[degrees]C, the boiling heat flux model is applied for the wall temperature higher than [T.sub.ONB] = 100[degrees]C, which is the boiling point of water under the standard atmospheric condition [14]. As the CT outer wall temperature ([T.sub.w]) increases in Figure 5, the boundary condition models are switched such as liquid convective, nucleate boiling, transition boiling, and film boiling. In Table 1, the models are listed with references.

2.4.1. Churchill and Chu Correlation. If the wall temperature is less than the boiling temperature, or [T.sub.w] < [T.sub.ONB] = 100[degrees]C, the wall heat coefficient is just that of a cylinder under the single-phase liquid free convection. There are many classical models, and Churchill and Chu Correlation [8] is applied for the mean convection heat transfer coefficient, [h.sub.m]:

[h.sub.m] = [k/D][Nu.sub.m] (11)

where k = 0.6627 W/(m x K) is thermal conductivity; D = 0.131 m is the outer diameter of CT; and [Nu.sub.m] is the mean Nusselt number:

[Nu.sub.m] = [[0.6 + [0.387[Ra.sup.1/6.sub.D]/[{1 + [(0.559/Pr).sup.9/16]}.sup.8/27]]].sup.2] (12)

In (12), the Nusselt number is correlated as a function of Rayleigh number, [Ra.sub.D], and Prandtl number, Pr:

[Ra.sub.D] = [g[beta]([T.sub.w] - [T.sub.[infinity]])[D.sup.3]/[v.sup.2]]Pr (13)

where g = 9.81 m/[s.sup.2], gravitational acceleration; [beta] = 0.00291/K, thermal expansion; [T.sub.[infinity]] = 70[degrees] C, the temperature of reservoir; and v = 4.03 x [10.sup.-7] [m.sup.2]/s, kinematic viscosity. The Prandtl number of heavy water is approximated as Pr = 2.417.

2.4.2. Rohsenow and Zuber Correlation. In the range of [T.sub.ONB] [less than or equal to] [T.sub.w] [less than or equal to] [T.sub.CHF], Rohsenow correlation [9] is used for the heat convection coefficient:

[h.sub.NB] = C[([T.sub.w] - [T.sub.ONB]).sup.2] (14)

where the coefficient is calculated from fluid and steam values: C = 149.2 W/([m.sup.2][K.sup.3]), and the subscript NB represents nuclear boiling.

The critical heat flux temperature ([T.sub.CHF]) is calculated from the solution of a simple equation based on Zuber correlation [10] where the critical heat flux ([q".sub.CHF]) is specified to 2.417 x [10.sup.6] W/[m.sup.2] from properties of the saturated steam. From the relation in (14) on [h.sub.NB] that is a quadratic function on the differential temperature,

[q".sub.CHF] = [h.sub.NB]([T.sub.CHF] - [T.sub.ONB]) (15)

Then, solving (15), [T.sub.CHF] = 125[degrees]C is temperature at the critical heat flux.

2.4.3. Bjornard and Griffith Correlation. In the transition region between the critical heat flux and the film boiling, or [T.sub.CHF] < [T.sub.w] < [T.sub.WET], Bjornard and Griffith [11] suggest a simple model of interpolation with the temperature ratio, [lambda] = [([T.sub.WET] - [T.sub.w]).sup.2]/[([T.sub.WET] - [T.sub.CHF]).sup.2]:

[h.sub.TR] = [lambda]([[??].sub.CHF]/[T.sub.w] - ([T.sub.[infinity]] + [T.sub.ONB])/2) + (1 - [lambda])[h.sub.FILM] (16)

where the critical heat flux ([q".sub.CHF]) in (16) is calculated from Zuber correlation in (15).

All the models are plotted as the heat flux and the heat power per unit area versus wall temperature in Figure 6. The heat increases about 40 times at the critical heat flux more than that at the boiling temperature, decreasing in the transitional boiling, and increases slowly in the film boiling.

2.4.4. Bromley Correlation. For the range of [T.sub.w] [greater than or equal to] [T.sub.WET], Bromley correlation [12] is used for the expression of heat convection coefficient:

[mathematical expression not reproducible] (17)

where the subscripts f and g stand for fluid (or saturated liquid water) and gas (or saturated pure steam) phase, respectively: therefore, [h.sub.fg] = [h.sub.g] - [h.sub.f] under a given saturated temperature, [T.sup.sat] [approximately equal to] [T.sub.ONB]. In the vapor mixture, denoted by the subscript v, all the values are linearly interpolated with the temperature difference, [T.sub.w] - [T.sub.ONB] from the empirical function in the superheated steam table.

The film boiling temperature is specified as [T.sub.WET] = 289.1[degrees]C, which is obtained from the graphical intersection point with curves generated with (16)-(17) such as Figure 6. This can be computed only after just several iterations.

3. Result and Discussion

3.1. Contact of PT and CT. In Figure 4, the contact occurs at 74.7 s, and the distribution of stress is reversed before (a) and after (b) the event in the radial direction. The maximum von Mises stress lies at the outer surface of PT before the contact: see Figure 4(a), but, however, it is located at the inner surface of PT after the contact: see Figure 4(b). The reason originates from the discontinuity of Shewfelt and Godin creeping model in Figure 3. The thermal expansion rapidly drops near T = 1123 K = 850[degrees]C, and this shock effects on the reversal of stress in the radial direction such as Figure 4(b).

The contact conductivity gradient elucidated in Section 2.3 is defined here to artificially induce a critical heat flux.

[mathematical expression not reproducible] (18)

where [t.sub.c] and [DELTA]t are the contact time and the time duration of amplified heat conductivity gradient, [DELTA]t = 0.1 s in this research, respectively.

3.2. Temperature on PT and CT. In the counterpart experiment [3], the temperature is measured at five section rings in spanwise stations of the apparatus shown in Figure 1(a) where the selected data in Figures 7(a) and 7(b) are selected at the midspan. The computational model is a two-dimensional but eccentric one, so the difference in the circumference is not simulated. Instead, we change the amplification value (Amp) in (18). However, the slight change in Amp does not show any remarkable change. When Amp = 20, the contact time and the peak temperature on the inner surface of PT are well predicted with small errors in Figure 7(a), but the transient temperature is affected with the high temperature at the outer boundary of CT. The combination of PT and CT is cooled in the computation faster than that of the experiment. As shown in Figure 7(b), the outer coolant water is quickly transit to the film boiling over the critical heat flux, coming back to nuclear boiling within ten seconds.

The mode at the outer boundary sensitively depends on the thermal resistivity at the contact surface between PT and CT. During the first stage of contact, for 0.1 seconds, it never reaches to the critical heat flux if we reduce the amplification factor below 15 (Amp [less than or equal to] 15): see Figures 8(a) and 8(b). The effect of outer boundary temperature is reduced (Figure 8(a)), and the overshoot of outer temperature is far mitigated for the smaller amplification (Figure 8(b)). As the thermal expansion is not uniform in the three-dimension, the real experimental data lies between these two amplifications at each azimuthal location, and the thermal spots distribute from nucleate to film boiling.

4. Concluding Remark

The computational simulation has been performed for the IAEA-ISCP benchmark problem concerning a severe accident such as LOCA with the commercial code COMSOL multiphysics. The pressure tube (PT) expands with the radiative heat transfer from the graphite heater, soon contacting with the Calandria tube (CT), and the boundary condition at the outer CT is modeled with various physical correlations. The experimental measured data are compared with the result of numerical analysis.

The expansion of PT depends on the creeping of Zircaloy at high temperature, and the Shewfelt and Godin model has a genuine discontinuity at the maximum PT temperature. Using four-point Gaussian quadrature, the numerical integration for the thermal expansion shows a sharp discontinuous point to reverse the stress distribution in the radial direction neat the contact. The thermal resistance at the contact point is so important for the change of modes in the boiling process of the coolant water that the amount of change in the amplification of thermal conduction gradient at a thin contact layer, if it exceeds some threshold, brings out sensitively a considerable difference at the outer heat flux of CT.

The heat transfer observed in the parallel experiment only fits with a tuning of an artificial parameter of amplification factor for the heat conductivity gradient at the early stage of contact. The instant discharge of heat exceeding the critical heat flux can only activate the film boiling, or the reservoir water will remain in the nucleate boiling. However, the numerical data of transient surface temperature are compared with the experimental measurement, showing reasonable agreement overall in spite of much approximation of geometry and physics.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was supported by (1) the Human Resources Development Program (Grant no. 20174010201350) of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grants funded by the Korea government (Ministry of Trade, Industry, and Energy) and (2) the Space Core Technology Development Program (NRF-2017M1A3A3A03015448) of Korea.


[1] "Republic of Korea: Nuclear Power Reactors -Alphabetic," IAEA Power Reactor Information System (PRIS),

[2] H. Z. Fan, R. Aboud, P. Neal, and T. Nitheanandan, "Enhancement of the Moderator Subcooling Margin using Glass-peened Calandria Tubes in CANDU Reactors," in Proceedings of the 30th Annual Conference of the Canadian Nuclear Society, vol. 2009, Calgary Canada, 2009.

[3] T. Nitheanandan, "Proposal for the IAEA International Collaborative Standard Problem on HWR Moderator Subcooling Requirements to Demonstrate Backup Heat Sink Capabilities of Moderator during Accidents," International Atomic Energy Agency, 2012.

[4] COMSOL Multiphysics, "Structural Mechanics User's Guide," Theory for the Solid Mechanics Interface, pp. 160-170, 2012.

[5] H. T. Kim, S.-M. Chang, and J. H. Park, "Unsteady two-dimensional multiphysical simulation of a pressure tube model expanded to contact with the outer concentric tube," Journal of Nuclear Science and Technology, vol. 53, no. 4, pp. 580-591, 2016.

[6] R. Shewfelt, L. Lyall, and D. Godin, "A high-temperature creep model for Zr-2.5 wt% Nb pressure tubes," Journal of Nuclear Materials, vol. 125, no. 2, pp. 228-235, 1984.

[7] T. H. G. Megson, Aircraft Structures for Engineering Students, Edward Arnolds, Bristol, UK, 1985.

[8] S. W. Churchill and H. H. S. Chu, "Correlating equations for laminar and turbulent free convection from a vertical plate," International Journal of Heat and Mass Transfer, vol. 18, no. 11, pp. 1323-1329, 1975.

[9] W. M. Rohsenow, "A Modeling of Correlating Heat Transfer Data for Surface Boiling Liquids," Transactions of ASME, vol. 74, pp. 969-976, 1952.

[10] N. Zuber, "Hydrodynamic Aspects of Boiling Heat Transfer," AECU-4439, NSA-14-007508, 1959.

[11] A. T. Bjornard and P. Griffith, "PWR Blowdown Heat Transfer," in Proceedings of the ASME Symposium on Thermal and Hydraulic Aspects of Nuclear Reactor Safety, vol. 1, pp. 17-41, ASME, New York, NY, USA, 1977.

[12] L. A. Bromley, "Effect of Heat Capacity of Condensate," Industrial & Engineering Chemistry, vol. 44, no. 12, pp. 2966-2969, 1952.

[13] K. J. Geelhood, C. E. Beyer, and W. G. Luscher, "PNNL Stress/ Strain Correlation for Zircaloy," Pacific Northwest National Laboratory, PNNL-17700, 2008.

[14] T. G. Beuthe and B. N. Hanna, "CATHENA MOD-3.5c/Rev 0 Theoretical Manual," CANDU Owners Group Report, COG-99-007, 1999.

Hyoung Tae Kim, (1) Se-Myong Chang (ID), (2) Young Woo Son, (2) and Taegee Min (3)

(1) Thermal Hydraulic and Severe Accident Research Division, Korea Atomic Energy Research Institute, Republic of Korea

(2) Department of Mechanical Engineering, Kunsan National University, Republic of Korea

(3) R&D Center, S&H Co. Ltd., Republic of Korea

Correspondence should be addressed to Se-Myong Chang;

Received 16 June 2018; Revised 23 September 2018; Accepted 30 September 2018; Published 1 November 2018

Academic Editor: Arkady Serikov

Caption: Figure 1: IAEA/ICSP benchmark for this simulation: (a) schematic drawing and (b) heater power versus time.

Caption: Figure 2: COMSOL computation: (a) boundary condition and (b) grids system.

Caption: Figure 3: Thermal expansion coefficient as a function of temperature.

Caption: Figure 4: Von Mises stress: (a) 72 s, before contact, and (b) 74.7 s, after contact.

Caption: Figure 5: Schematic diagram of boiling heat flux models.

Caption: Figure 6: Boiling heat flux curve at the outer surface of Calandria tube.

Caption: Figure 7: Measured and computed temperature at (a) the inner surface of PT and (b) the outer surface of CT, Amp = 20.

Caption: Figure 8: Dependency on the amplification factor of thermal conductivity gradient for 0.1 seconds of contact at (a) the inner surface of PT and (b) the outer surface of CT.
Table 1: Boiling heat flux models.

Wall temperature             Phenomena

[T.sub.w] <              Liquid convection

[T.sub.ONB] [less         Nucleate boiling
than or equal to]
[T.sub.w] <

[T.sub.w] =              Critical heat flux

[T.sub.CHF] <           Transitional boiling
[T.sub.w] <

[T.sub.w] [greater          Film boiling
than or equal to]

Wall temperature             Correlation             Reference
                                                  (published year)

[T.sub.w] <               Churchill and Chu          [8] (1975)

[T.sub.ONB] [less             Rohsenow               [9] (1952)
than or equal to]
[T.sub.w] <

[T.sub.w] =                     Zuber               [10] (1959)

[T.sub.CHF] <           Bjornard and Griffith       [11] (1977)
[T.sub.w] <

[T.sub.w] [greater             Bromley              [12] (1952)
than or equal to]
COPYRIGHT 2018 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Kim, Hyoung Tae; Chang, Se-Myong; Son, Young Woo; Min, Taegee
Publication:Science and Technology of Nuclear Installations
Date:Jan 1, 2018
Previous Article:Pressure Waves due to Rapid Evaporation of Water Droplet in Liquid Lead Coolant.
Next Article:Analysis of Steam Explosion under Conditions of Partially Flooded Cavity and Submerged Reactor Vessel.

Terms of use | Privacy policy | Copyright © 2020 Farlex, Inc. | Feedback | For webmasters