Abstract
Numerical study on AHWR fuel bundle has been carried out to assess influence of circumferential and cross flow rewetting on the conduction heat transfer. The AHWR fuel bundle quenching under accident condition is designed primarily with radial jets at several axial locations. A 3D (, , ) transient conduction fuel pin model has been developed to carry out the study with a finite difference method (FDM) technique with alternating direction implicit (ADI) scheme. The single pin has been considered to study effect of circumferential conduction and multipins have been considered to study the influence of cross flow. Both analyses are carried out with the same fluid temperature and heat transfer coefficients as boundary conditions. It has been found from the analyses that, for radial jet, the circumferential conduction is significant and due to influence of overall cross flow the reductions in fuel temperature in the same quench plane in different rings are different with same initial surface temperature. Influence of cross flow on rewetting is found to be very significant. Outer fuel pins rewetting time is higher than inner.
1. Introduction
Rewetting of hot surface is a process in which a liquid wets a hot surface by displacing its own vapour that otherwise prevents the contact between the solid and liquid phases. This has generated immense interest in studying rewetting through both theoretical simulation by Yamanouchi [1] and Coney [2] and experimental study which was carriedout by Yamanouchi [1] and Duffey and Porthouse [3]. Falling film rewetting for several vertical geometries such as plates [2, 4], rods by Blair [5], and tubes by Satapathy and Sahoo [6] had been modeled by a number of researchers. In general, in all models, a moving rewetting front that divides the solid into two distinct regions is considered. Most of the models also consider a constant rewetting velocity that reduces the problem into a quasistatic one. Initial efforts were made to formulate one-dimensional conduction models by Yamanouchi [1] that are reasonably successful in correlating rewetting phenomena at low Peclet number. Tien and Yao [4] presented the asymptotic solutions of a two-dimensional conduction model which clearly demonstrates the different physical pictures for the cases of high and low coolant flow rates. A variety of techniques have been used for solving two-dimensional conduction models for falling film rewetting. Some of the important studies are elaborated. Because of mathematical difficulty, most two-dimensional analyses are either approximate or numerical ones. Duffey and Porthouse [3] first considered solving the rewetting problem by separation of variables. They retained only the first term in the series solution. However, Coney [2] reported that using a small number of terms in a series yields inaccurate results due to slow convergence of the series. An approximate solution to the same model for a cylindrical rod was presented by Blair [5]. Tien and Yao [4] first applied the Wiener-Hopf technique to a two-dimensional rewetting problem of a rectangular slab, while an exact solution to the same problem was presented by Castiglia et al. [7], employing the method of separation of variables. Numerical solutions of conduction controlled rewetting were provided by Satapathy and Sahoo [6], Thompson [8], and Raj and Date [9] by using the finite difference technique. Heat balance integral method (HBIM) is one of many semianalytical methods used to solve conduction problems by Eckert et al. [10, 11]. This is analogous to classical integral technique used for fluid flow and convective heat transfer analysis. A numerical study has been made to investigate the effect of internal heating and precursory cooling during quenching of an infinite tube and was studied by Satapathy and Sahoo [6]. An experimental assessment of rewetting of an advanced heavy water reactor (AHWR) specific nuclear fuel bundle with radial jets has been carried out by Patil et al. [12]. The experimentation involved cooling of a 54-rod bundle by in-bundle injection and demonstrated that quenching occurs for all fuel pins along the radial direction. But maximum temperature is limited to be 400°C and no heat generation has been considered during the quenching of experiment.
The numerical study presented in this paper aims to bring out a comparison in rewetting pattern for a radial jet rewetting versus axial rewetting for AHWR fuel bundle. The proposed reactor is a vertical pressure tube type, heavy water moderated boiling light water coolant reactor. Fuel bundle is housed in the pressure tube (PT) which in turn is housed in calandria tube (CT). The 3.5 m height fuel bundle is having 54 fuel pins arranged in three different concentric rings as shown in Figure 1. The inner, middle, and outer rings have 12, 18, and 24 fuel pins and centre of the fuel pin has a water rod which is used to inject emergency coolant in radial direction at different elevations of the fuel bundle under a pipe break scenario. The water rod is having eight holes of 1.5 mm dia at one elevation along the circumference of water rod. The radial injection is designed for 13 different axial locations along the 3.5 m length of water rod. More details of AHWR fuel bundle and its injection arrangement are furnished by Sinha and Kakodkar [13]. Bottom reflooding is also an option with the designer where water may be injected from bottom most holes instead on 13 axial elevations.
Under this study, a 3D (, , ) transient conduction fuel pin model has been developed based on finite difference method (FDM) technique with alternating direction implicit (ADI) scheme. The number of nodes in , , and -direction has been considered to be 50, 20, and 50, respectively. A 1.6 ms time step has been used as predicted from stability criteria. A total time of 4 hrs is required for individual rewetting study.
Single pin from first circle has been considered for which bottom and radial reflooding are studied. The 1st circle pin is selected as it experiences a strong radial jet and experiences a large circumferential temperature gradient. For bottom rewetting case all the pins will experience the same velocity front irrespective of fuel pin location. However, the developed computational model is equally applicable for 2nd and 3rd circle pins with different boundary conditions. Both analyses are carried out with same fluid temperature and heat transfer coefficients as boundary conditions. The single fuel pin has been considered for circumferential and axial rewetting. The three-dimensional partial differential equation for unsteady state conduction equation for cylindrical rod is as follows:
Alternating direction implicit (ADI) numerical method has been adopted. With an ADI method the heat diffusion equation is first solved implicitly in the -direction while leaving the other two directions explicit. The heat diffusion equation is then solved implicitly in a similar way in the - and -direction. This scheme reduces a three-dimensional problem to a series of one-dimensional implicit problems. Figure 2 shows the three different steps in ADI.
Many methods [14–17] are also available in ADI for solving the heat diffusion equation, such as Peaceman-Rachford pure-ADI method, Brain method, Douglas method, and the method based on superposition principle by Jules Thibault [18]. All above schemes having the same problem of conditionally stable criteria as right side of equations have a negative coefficient. This paper discusses conduction effect on axial, circumferential direction as well as influence of rewetting of fuel pins in different concentric circles due to cross flow model.
The tridiagonal matrix has been used to solve the equation which corresponds to given direction. Normally in ADI scheme, first solution is obtained in -direction in an implicit manner and leaves both directions explicit. But here all the possible six combinations, like --, --, --, --, --, and --, have been used to solve the heat diffusion equation and compared results.
2. 3D Model Development
2.1. Finite Difference Formulation
FDM is used to formulate the heat diffusion equation in cylindrical coordinate (see (1)) and further semi-implicit method used to formulate in ---direction. The differential equation has been divided into the three categories centre of fuel pin: between centre of fuel pin and surface, and surface of fuel pin. The length increment in , , -directions are , , , respectively. The semi-implicit scheme has been used to formulate the differential equation in , , and -directions and is represented by (2), (5), and (8) respectfully. As per ADI norms, first 1/3 time increment has been taken to solve the temperature at different grid points in -direction, next 1/3 time increment for solving the temperature in -direction, and last time step for solving the temperature in -direction. For stability criteria, the coefficient of , , in right side of the equation must be positive.
Finite difference formulation between centre and surface of fuel rod: -direction or Condition for stability -direction or Condition for stability -direction or Condition for stability
2.2. Centre and Surface of the Fuel Pin
The centre line temperature for fuel pin has been calculated by using the concept of energy balance and finite difference formulation is used for calculating the centre line temperature. It is shown by (11). The centre line grid point is represented by and next to centre line temperature is represented by and is shown in Figure 3. In a similar way, surface temperature of fuel pin has been formulated in (14): or Condition for stability At surface of fuel pin.
Implicit Scheme or
2.3. Boundary Condition and Solver
Once the formulation of differential equations is over then, with help of boundary, it is required to calculate the temperature at different nodes of fuel pin. There are two cases considered for validation of numerical code; first case is constant coolant temperature and heat transfer coefficient and second case is coolant temperature varies along the length of fuel pin when it flows along the length of pin. The governing equation which calculates the coolant average temperature between the corresponding nodes of fuel pin is represented by (17).
Energy balance at surface of fuel pin Equation (18) is used for calculating the outlet coolant temperature for a given volume. After applying the boundary condition, it is required to solve the discretised partial differential in each direction by ADI method. At each incremental time step, the differential equations form a tridiagonal matrix and are solved by Thomas algorithm (see (19)). In first time increment (i.e., 1/3 s) it solves temperature in -direction at different grids points and takes old values of temperature in , -direction. Similarly for higher time step (2/3 s), this algorithm solves in -direction and remaining directions take previous value: Different heat transfer correlations have been used to determine the heat transfer coefficient at the surface of fuel pin. Thom correlation (see (20)) is used for nucleate boiling and modified Bromley correlation (see (21)) is used for film boiling: where
2.4. Steady State Computation for Uniform and Nonuniform Heat Generation
Figure 4 shows the 3D temperature behavior in the fuel pin in different cases; the first case is uniform heat generation (1.16E + 8 W/m3) with a constant heat transfer coefficient and coolant temperature and second case is with sinusoidal heat generation with variation of coolant temperature along the length of the fuel pin.
(a)
(b)
Temperature in fuel pin has been predicted uniform across the fuel pin and at any cross section minimum and maximum temperatures are at the surface and centre of fuel pin, respectively. For sinusoidal heat generation, as shown in Figure 4, the fuel surface temperature is highest at the middle of the fuel pin since pin has a maximum heat flux. Temperature at any section of fuel pin along the length of fuel pin decreases in both directions from the middle of the fuel pin.
3. Benchmarking of 3D Model
A benchmarking exercise has been carried out for the 3D conduction model with the available 1D analytical solution for a steady state case. As analytical solution of 3D is not available, the benchmarking is carried out with 1D solution. The analytical solution for estimating the coolant temperature and fuel pin surface temperature along the length of fuel pin is represented by (23) [19] and used to calculate the fluid temperature and surface temperature for 1D analytical solution. The heat flux along the fuel pin has been considered as a cosine profile. The Fourier conduction equation and its analytical for calculating the temperatures along the radius of fuel pin solution are given in (24) and (26), respectively. Table 1 provides the input parameters for this exercise: 1D Fourier conduction equation is Apply boundary condition At outer surface.
Heat generation in body = heat transfer to coolant Two-step benchmarking exercise has been done. In the first step, boundary conditions used for calculating the fuel radial temperature distribution for 3D model and analytical solution are compared and in the subsequent step the radial temperature distribution is compared. Figure 5 shows the coolant temperature and fuel surface temperature variation along the length of the fuel pin.
The calculated coolant temperature and surface temperature using (17) and (18) act as boundary conditions for the 3D code. A good agreement is found between the numerical code result with analytical solution obtained using (23). Figure 6 shows the comparison of temperature distribution in -direction with 3D numerical model and analytical solution for fuel rod. Result shows a good agreement between the two solutions.
Further, this single pin model has been extended to multiplepins model. The multiplepin model considers that pins are arranged in three concentric rings with radial jet simulation. A validation exercise has been carried out with experimental results reported by Patil N. D. The experimental case with initial temperature of 168°C for all pins with an injection flow rate of 73 lpm has been considered. The validation exercise results are shown in Figure 7 at an elevation of 2.4 m.
The exercise shows that the multipin model is able to capture the rewetting pattern well. However, in this study circumferential rewetting pattern could not be produced as measurements across the circumference of a simulated fuel pin are not being done/reported.
4. Study on Influence of Radial Rewetting on Conduction Heat Transfer Assessment
As the circumferential conduction phenomenon is predicted by numerical technique, influence of ADI solution technique on the above mentioned findings is investigated. The schematic for radial reflooding is shown in Figure 12. For this study, boundary conditions like high rewetting heat transfer coefficient of 30000 W/m2 K (only water) are applied over half of the circumference (180°) of fuel pin where the water jet is likely to impinge. At the same time, a low heat transfer coefficient of 500 W/m2 K (only steam) is applied on the other half of the circumference (180°) where the jet is unlikely to reach. A linear variation of heat transfer coefficient from 500 W/m2 K (only steam) to 30000 W/m2 K (only water) for the water impinged surface is considered as described earlier. Variation of fluid temperature from 300°C (only steam) to 30°C (only water) is assumed over 20 s period at the node where rewetting front has reached. Figure 8 shows the transient cooling temperature for four circumferential locations covering the half of the circumference which experienced the water jet.
The result shows that for the chosen boundary conditions, influence of circumferential conduction is significant. Temperature of the unwetted portion dropped by 250°C with rewetting of the other half of the fuel over a minute of the beginning of rewetting of font section. The influence of one-sided rewetting has caused the maximum temperature of fuel pin to shift from centre to unrewetted portion of fuel pin, as shown in Figure 9.
(a)
(b)
5. Study on Influence of Cross Flow in Fuel Bundle
This work has been extended to multipins which are arranged in three concentric rings. The (1/4)th portion of fuel bundle, direction of radial injection of coolant flow, is shown in Figure 10.
To study the influence of the cross flow and axial flow on rewetting of fuel pin in coolant channel, the four and five control volumes have been modeled in radial and axial direction, respectively. The (1/5)th length of fuel pins is attached to one axial control volume. The innermost volume in radial direction is attached with only half of the total heat from 1st ring of fuel pin, whereas the second radial volume receives heat from remaining portion of the total heat from 1st ring of fuel and heat from front side of 2nd ring. It is a similar way of modeling for third and fourth volumes. In the axial direction, each innermost volume receives coolant from water tube. The coolant then goes into radial as well as axial direction. In top control volume coolant is heated by hot fuel pins and comes down due to gravitational force, mixes with cold coolant injected into second volume. In radial flow modeling, coolant first gets heated from first volume and then heated fluid enters into second volume and similar way to third and fourth volumes. In this flow distribution, top volume receives lowest coolant and next volume just below top volume gets higher, bottom most control volume receives highest amount of coolant. Two cases have been considered for modeling of power in the fuel pins; in the first case all pins in different rings have uniform power and in second case radial power variation (0.69 : 0.79 : 1.33) in different fuel rings has been considered.
Figure 11 shows the fuel surface temperature behavior in different rings. Initial temperature of fuel pins in all the rings is 666.7°C. The coolant flow rate for rewetting is considered to be 90 lpm [20]. Fuel surface temperature starts decreasing with injection of coolant into the coolant channel. Inner ring pins quenching takes place at s. on both sides of fuel pin. Due to the influence of circumferential conduction, fuel surface temperature in opposite sides of fuel pins of first ring falls slower than front side. Quenching time for second ring and third rings is 3.3 s and 6.8 s, respectively. Initially surface temperature of fuel pin falls slowly since heat transfer coefficient is poor due to the presence of strong steam film between heated surface and coolant. Once fuel surface temperature reaches below the Leidenfrost point, surface gets rewetted which results in rapid fall of fuel surface temperature. The axial fuel surface temperature of outer ring at different locations is shown in Figure 12. Rewetting time of all nodes is close to s but temperature fall rate at top node is a little slower than lower node. It is due to the lower coolant amount received by top node.
During quenching of inner ring fuel pins, the heat released from the fuel pins is more as compared to the outer ring fuel pins. Therefore, the rise in coolant temperature observed during the quenching period is more for the inner fuel pins. Figure 13 shows the coolant temperature in different volume. The average temperature of coolant rises when it moves from first volume to second volume and continues till the last volume in radial direction. Fourth volume gets maximum coolant temperature. Oscillation of coolant temperature is due to continuous supply of coolant to volume.
For the case of rewetting with radial power distribution, the initial temperatures obtained in 1st, 2nd, and 3rd rings are 605°C, 635°C, and 742°C [16], respectively, as shown in Figure 14. The same mass flow rate (90 lpm) has been considered for quenching of fuel pins. Similar to the previous case, during initial period of quenching, strong steam film is formed between the fuel surface and liquid so that the fuel surface temperature decreases slowly because of poor heat transfer. Once surface temperature reaches Leidenfrost point it rewets the surface. The rewetting time of heated fuel surface temperatures is 1.9 s, 3.1 s, and 7.7 s. Influence of circumferential conduction effect has been also observed. Front side clad surface temperature falls faster than opposite side.
6. Conclusion
Numerical studies on the influence of circumferential effect and cross flow for AHWR fuel bundle reflooding phenomena has been studied with the help of 3D conduction numerical model. A circumferential conduction is found to be significant for radial jet rewetting. The conservative assessment shows that radial jet is able to bring down the unrewetted portion temperature by 250°C. A radial direction jet is found to be effective to cool the hot surface. The numerical studies on influence of cross flow have been studied. It has been found that even for the clad surfaces having the same initial surface temperature, rewetting time of inner fuel pins is lesser than the outer pins and influence of circumferential conduction has been also observed. The axial fuel surface temperature behavior at different locations is observed to be nearly same.
Nomenclature
| : | Area (m2) |
| : | Specific heat capacity (J/kgK) |
| : | Local acceleration due to gravity (m/s2) |
| : | Gravitation constant (m/s2) |
| : | Heat transfer coefficient at outer surface (W/m2K) |
| : | Latent heat of vaporization (W/m) |
| : | Thermal conductivity of steam (W/mK) |
| : | Thermal conductivity of fuel (W/mK) |
| : | Length of fuel (m) |
| : | Coolant flow rate (kg/s) |
| : | Pressure (MPa) |
| : | Local heat generation in fuel (W/m3) |
| : | Fuel radius changing in -direction |
| : | Outer radius of fuel pin (m) |
| : | Temperature (K) |
| : | Fuel surface temperature (K) |
| : | Coolant outlet temperature (K) |
| : | Coolant inlet temperature (K) |
| : | Average coolant temperature |
| : | Temperature difference between fuel surface and saturation temperature of coolant (°C) |
| : | Time (s) |
| : | Time increment (s) |
| : | Volume (m3) |
| : | Distance increment in fuel (m) |
| : | Distance increment in -direction (m). |
| : | Density of steam (kg/m3) |
| : | Density fuel (kg/m3) |
| : | Surface tension (N/m) |
| : | Angle (radians) |
| : | Thermal diffusivity (m2/s) |
| : | Density difference between liquid and steam (kg/m3) |
| : | Angle increment in circumference direction (radians). |
| . |
| : | Fuel |
| : | Mesh point in -direction |
| : | Mesh point in -direction |
| : | Mesh point in -direction. |
| : | Time step. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Authors acknowledge Dr. R. K. Singh, Mr. H. G. Lele, and Dr. P. K. Vijayan giving their valuable suggestions during the development of code.