Investigation of the channel disassembly behaviour of Indian 200MWe PHWR – A numerical approach

Nuclear Engineering and Design 339 (2018) 137–149

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Investigation of the channel disassembly behaviour of Indian 200MWe PHWR – A numerical approach

T

⁎

Ankit R. Singha, Pradeep K. Sahoob, , Andallib Tariqa a b

Department of Mechanical and Industrial Engineering, Indian Institute of Technology Roorkee, Roorkee, Haridwar, India Department of Mechanical, Energy and Industrial Engineering, Botswana International University of Science and Technology, Palapye, Botswana

A R T I C LE I N FO

A B S T R A C T

Keywords: Disassembly Coupled field Channel deformation

In absence or non-functioning of a failsafe system, the postulated accident in Pressurised heavy water nuclear reactor may lead to severe degradation of structural integrity of channels. It is utmost important to discover the nature of channel disassembly or sagging behaviour following the loss of channel integrity. An experimental facility is presently being developed at Indian Institute of Technology Roorkee (IIT R), India to study such a severe event. The experiments are grouped under two categories, one being study under inert environment to check channel deformation due to material flow at the higher temperature and other one is the investigation on the effect of an metal oxidation in presence of a steam surrounding. The Latter is the case similar to an actual condition exists in a reactor core under the postulated accident scenario. Initially, experiments will be carried out on a single channel, and later it will be extended up to three channels arranged one below another. This paper discusses various aspects considered towards the facility development. A detailed numerical analysis of a single channel sagging in an inert environment is carried out. It provides information beforehand about the deformation pattern which has helped in designing the experimental facility. The numerical model was simulated using direct coupled fields approach inside finite element domain to get a preliminary insight into structural and thermal coupled response on the channel deformation.

1. Introduction

Calandria tube sheets of end shield at both ends while PT is placed such a way that its one end is fixed, and other end is free to move in axial direction to accommodate the thermal expansion. CO2 gas flows through the annulus of a PT and CT, and it acts as thermal barrier to restrict lateral heat transfer to moderator during normal operation. Primary Heat Transport System (PHTS) is schematically shown in Fig. 1 of a 220MWe IPHWR reactor. It has two passes through the reactor; in one pass it connects 153 channels from south to north, and in other pass the flow is from north to south direction (Mukhopadhyay et al., 2002). An accident scenario pertaining to PHWR can be imagined in two ways; one way is in which moderator system is functioning, and thus availability of moderator surrounding Calandria tubes removes decay heat during Loss of Coolant Accident (LOCA) events. Such scenario is termed as Limited Core Damage Accident (LCDA). Second condition is imagined when the moderator system malfunctions, and along with there is a loss of Emergency Core Cooling System (ECCS). Such conditions lead to continuous decrease in moderator level with the passing transient and this result in a reduction in the decay heat removal from fuel bundles and from coolant channels. This scenario is termed as a Severe Core Damage Accidents (SCDA) events (IAEA,

A typical 220MWe core of Indian Pressurised Heavy Water Reactor (IPHWR) consists of 306 horizontal channels and each channel has two concentrically placed tubes. The inner tube is called a Pressure Tube (PT) while the outer tube is called a Calandria Tube (CT). Twelve fuel bundles are placed inside the PT along its length. Out of these, ten bundles lay inside the active region of the core (Bajaj and Gore, 2006). The channels assembly are placed inside a large cylindrical vessel termed as Calandria. The Calandria is filled with heavy water (termed as moderator) at the atmosphere pressure. The moderator acts as a heat sink under an accident conditions (Thomson and Kohn, 1983). The Calandria is kept inside a large shield tank called Calandria vault that stores a large volume of light water at atmospheric pressure. The Calandria vault acts as the ultimate heat sink under severe accident condition and it avoids or delays the melting of the Calandria vessel (Nitheanandan and Brown, 2013). Pressurized heavy water coolant flows through PT passage, and it takes heat generated from fuel bundles during controlled fission reaction. Then, exchange heat to secondary side inside steam generators. The CT is rolled joint and fixed to

⁎

Corresponding author. E-mail addresses: [email protected] (A.R. Singh), [email protected] (P.K. Sahoo), [email protected] (A. Tariq).

https://doi.org/10.1016/j.nucengdes.2018.09.008 Received 23 March 2018; Received in revised form 7 September 2018; Accepted 10 September 2018 0029-5493/ © 2018 Published by Elsevier B.V.

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

[Ct] [Kt] [Kut] [Ctu] {u} {T} {F} {Fth} {Q} {Qted}

Nomenclature a, b, c ε̇ σ ε T [M] [C] [K]

Creep coefficients Creep strain rate (/s) Stress (MPa) Equivalent creep strain Temperature (K) Structural mass matrix Structural damping matrix Structural stiffness matrix

Thermal specific heat matrix Thermal conductivity matrix Thermo-elastic stiffness matrix Thermo-elastic damping matrix Displacement vector Temperature vector Force vector Thermal strain force vector Heat flow rate vector Heat generation rate vector for thermo-elastic damping

Fig. 1. Schematic of PHTS of a typical IPHWR.

Semenov. In the postulated event of LOCA with the loss of ECCS and loss of moderator cooling system, Pressure inside the PHTS remains high for quite certain duration after reactor trip. This lead to ballooning of PTs in the top most rows of fuel channels, and PTs make contact circumferentially to respective CTs. Still pressurised channel might fail due to crack or burst which in turn reduces the pressure in whole PHTS. Thus, remaining fuel channels experience lower inside pressure. Due to this, PTs start to sag as explained earlier, and it lead to sagging mode of failure for channels. Few experiments are reported in literature (Gupta et al., 1996; Nandan et al., 2010; Negi et al., 2017) for study of deformation of PT under LCDA event for IPHWR. However, no experimental investigation has been carried out to develop the understanding for progression of channel sagging in case of SCDA event. To address this gap, it is planned to study the deformations aspects under sagging event through the experiments with a scaled-down (1:3) fuel channels of 220MWe IPHWR. The experimental facility is presently being developed at IIT-Roorkee, India. In order to perform scaled-down experiments, there are certain criteria related to geometry, mechanical properties, and applied load on system that are needed to be considered a priori with their implication. This paper takes into account of these important criteria and tries to explain in detail the effect of these parameters on system design. The sagging of a single channel in an inert environment is numerically analysed to acquire information beforehand performing the actual experiments. To imitate the postulated severe accident events in lab, one needs to know the power requirement to achieve such temperature transients. The presented simulation study helped to select the power

2008). This kind of events can also be possible in the case of prolonged Station Black-Out (SBO) scenario. With continuous decrease of moderator, upper channels will be exposed to dry-out condition leading to the deformation of channel under fuel bundle weight and its self-weight due to overheat. Upper channel contacts the channel below and this process continues as more and more channel gets uncovered. Afterward, channels may fail at the region of high stress level (Mathew et al., 2008) and get collected at the bottom of Calandria as the terminal debris. Channel disassembly following a postulated SCDA events was investigated experimentally (Mathew et al., 2001) for CANDU6 reactor channel. The experiment was conducted with scaled down channels having a scale of one fifth. The composite channel was assumed in sense that PT ballooned and contacted completely to CT peripherally. Experiments showed that above 800 °C, channel sag was significant, and bundle geometry did not add rigidity to PT deformation. Which can be anticipated as the presence of a large gap between bundles and PT. Channel break-up temperature was observed to be less than 1400 °C. A simulation using beam element was also performed to validate the experimental findings (Mathew et al., 2003). A simulation model was developed to check criteria for fuel channel break under the load of sagged contacting channels (Singh et al., 2011). In total, five channels were simulated placing one above another. The analyses were performed considering uniform channel temperature, and channel disassembly was assumed to happen when total strain at particular location along the channel exceeded failure strain for a given temperature. Reference of failure strain as developed for PT material (Zr-2.5%Nb), in the simulation was taken from Rodchenkov and 138

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

due to thermally induced creep of channel material, it is warranted to validate the creep equations sought for the simulation work. The methodology adopted for the validation of available creep correlation is also discussed in the paper.

range required to heat up the channel, which in turn helped in selection of DC rectifier output range. Under influence of high temperature transients, channel sags enormously under inert media owing to ductility of Zirconium materials. Simulation under inert condition helped to select the test section dimensions. Numerical model was simulated using direct coupled field approach inside finite element domain to get a preliminary insight into structural and thermal coupled response on the channel deformation. Since, this type of deformation are produced

2. Experimental facility and setup layout Experimental facility is schematically detailed in Fig. 2(a).

Fig. 2. Experimental facility. 139

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

εcṙ = f (σ , ε , t , T )

Trapezoidal parts are provided on front-side and top-side of the rectangular tank. The front section facilitates the mounting of digital and thermal imaging camera while the top hood facilitates the mounting of laser light generator. This type of configuration provides the required distance between the channel and camera lens to get the correct field of view (FOV) for cameras. Circular openings are provided on each part as an observation window for camera and laser light generator, and these opening are sealed with toughened quartz glass. End enclosures are designed to house the channel ends and corresponding connection, and to create an inert media inside PT passage and the annulus of PT and CT. As illustrated in Fig. 2(b), four cameras, two on back-section and two on front-section, will be used to measure sag of channels. The test section can be used to check the sagging of maximum three channels placed one above other. Channel end connections and test section walls are water cooled using a separate water cooled circuit. Test section walls will be flooded with the circulating water. An inert environment inside the PT and an annulus of PT-CT is maintained by purging Argon gas to mitigate the effect of oxidation by air. The channels are to be resistively heated using thyristor controlled DC rectifier (86 kW total capacity) with the help of Copper bus bars as shown in Fig. 2(a). Deformation of channels will be measured using edge detection technology with digital image processing. The images of transient deformation will be taken with four monochromatic digital cameras (Imperx, Inc. made). Two laser line generators operating in blue wavelength bands (Z-Laser, Inc. made) are mounted on top of channel for the generation of a straight line along top most periphery of the channel, and this laser line will be used as a tracer during image processing. Band filters will be installed over camera lenses to filter specific wavelength as of the laser line (400–450 nm). Use of band filter ensures accurate measurement of deformation at high temperature using digital imaging technique (Pan et al., 2011). The channel will also be thermally observed with the help of an infrared camera in near IR range to identify any possible failure (like crack initiation location) during the transient run. End load on the channel and more specific on CT end will be measured with a miniature type compressive load cell. S-type and Ktype thermocouples will be used for acquisition of temperature of PT and CT respectively. Wall temperature (likely to be below 400 °C) will be measured with the optical pyrometer. During oxidation (steam condition) experiments, exhaust steam from test setup will be passed through the cooling unit where steam will be condensed and simultaneously drained while Hydrogen will be forced through gas analyser to measure the amount of Hydrogen generated per volume basis. The condensation of steam will be carried out by re-circulating the chilled water that would be produced from chiller unit. All the data will be acquired continuously with DAQ system for further analysis. Flow chart for the experimental cases is illustrated in Fig. 3. The sagging experiments are broadly categorised into an inert media and steam media experiments. Single channel and multiple with maximum three channels will be performed under both the cases. The results from the inert media experiments will be used to validate the analytical or numerical findings; because, the effect of oxidation and hydrogen embrittlement are not considered in the numerical models.

(1)

Creep strain rate curve with time as abscissa has three distinct regions namely primary, secondary, and tertiary creep. The most engineering investigation are anchored to primary and secondary creep region. There are many empirical correlations are available which represent relation between the various dependent variables with creep strain rate. However, simplify the understanding procedure a generalised form of creep strain rate i.e. Norton form for the secondary creep is used, and it is expressed as Eq. (2). Temperature dependency on creep strain rate is given in the Arrhenius terme−c / T .

ε ̇ = a. σ b. e−c / T

(2)

A geometric scale of approximately 1/3rd is considered for the investigation of channel sagging. This includes scaling on length of the core, fuel bundles, and lattice pitch and diameters of PT, CT, and fuel bundles. For getting the desired material behaviour the scaled down PT and CT are fabricated with Zr-2.5%Nb and Zr-2 respectively, which are the same Zirconium alloys as use for the actual reactor channels. This ensures the same thermo-mechanical properties for both the tubes and by maintaining the exact composition in scaled down channel gives the same material flow characteristics during sagging. To maintain same mechanical stress level as in the prototype channel, tungsten fuel simulators are used to mimic weight of fuel bundles. Ten fuel simulators will be placed inside and along PT length. Keeping the same material compositions, stress intensity (σ ), and heatup rate as of the prototype reactor channel, Eq. (2) will result in the same creep strain rate; thus, scale down factor will only depends on geometric scaling. Herewith, a geometrical scaling of 1/3rd can be achieved on the deformation or sagging of a channel. The finalised dimensions for channels based on scaled down approach are listed in Table 1. 3. Coupled fields approach A numerical analysis was performed for a case related to a single channel heat up under inert media in order to gain an insight into deformation or sagging of the channel beforehand any experimentation. A coupled physics technique is used for numerical formulation and details of approach and formulation are discussed in the onward sections. There are two types of coupling exist in the scope of finite element physics, and these are direct coupling and sequential coupling (ANSYS, 2016a). Under sequential coupling, a result from one physics is used as load for the second physics. It allows the problem to be solved independently with different meshes for each one. However, if physics of a problem are strongly dependent on each other i.e. change in a parameter of one physics affects the strongly the parameters of second physics and vice-versa, then the sequential coupling is not the viable

2.1. Scale down channel approach It is observed that the sagging of channel would result as the channel temperature transient shoots up (IAEA, 2008). At temperature above the 30–40% of the melting point of Zirconium alloys, a thermally activated creep is a dominant cause of material flow. Henceforth, it is of a great importance to ensure that scaled down channel should replicate the creep behaviour of an actual size reactor channel. A creep strain rate generally depends upon variables such as stress, strain, time, and temperature, and it can mathematically be represented as Eq. (1).

Fig. 3. Tree structure of formulated experimental cases. 140

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

Table 1 Prototype and Scaled down dimensions of 220MWe IPHWR channel. Variables

Geometrical consideration Core length (m) PT inner diameter (mm) PT outer diameter (mm) CT inner diameter (mm) CT thickness (mm) vertical pitch (mm) Fuel simulator diameter (mm) Fuel simulator length (mm) Mechanical consideration Mid-Span Max. Stress of PT (MPa) Mid-Span Max. Stress of CT (MPa)

Prototype channel

Scaled-down channel

Table 2 Variation of k2 with temperature.

Scaling ratio

5.085 82.55 89.19 107 1.5 228.6 78.82

1.8 28.01* 30.03* 35.9* 0.5* 80 26

2.825 2.9472 2.97 3 3 2.86 3.03

495.3

175

2.83

21.597

21.5

1.0045

31.455

28.108

1.119076

k2 (MPa-1)

800 825 850 875 900

0.0119 0.0114 0.0167 0.0134 0.0100

ε 2.4 −26670 ε ̇ = 8 × 1010. σ . ⎛k2− ⎞ e ( T ) σ⎠ ⎝

(7)

k2 is a temperature dependent factor, and it is given in Table 2. Values of k2 can be assumed as zero for temperatures below 750 °C, as equation exhibits negligible strain rate for lower temperature region (Shewfelt et al., 1985).

(* As received dimensions).

3.2. Implementation of creep and its validation against analytical calculation

formulation to use. Such problems can be solved using direct coupling between the physics, and such formulation uses same mesh element for each physics. ANSYS v17.2 Mechanical solver is well capable to handle both the kinds of coupling. With usual notation, structural and thermal system equations in matrix form can be written as;

[M ]{u¨} + [C ]{u}̇ + [K ]{u} = {F }

(3)

[Ct ]{T ̇ } + [Kt]{T} = {Q}

(4)

Creep Eq. (7) is not in a standard form to be directly used in ANSYS solver (ANSYS, 2016b) environment. To incorporate such a model into ANSYS solver, a Usercreep subroutine, a FORTRAN program, has been developed and linked to solver path. Usercreep subroutine requires creep strain rate and its derivative with respect to stress and strain to calculate the global tangent matrix at each element integration points. Moreover, subroutine takes the solver variables like temperature and stress and calculates the incremental creep strain as output. The calculated strain is then used by solver to calculate the overall strain in model. ANSYS v17.2 uses implicit algorithm for creep calculation as it is being more robust, accurate, and unconditionally stable. It can easily handle a temperature dependent creep coefficients. However, it is necessary to check the feasibility of the subroutine in terms of its output for the creep strain. Hence, single mesh element model was developed to test the creep strain output of the subroutine with analytical calculation performed through Excel (MS Office16) spreadsheet. Details of the testing model are described in Fig. 4, it is a simple horizontal bar with a circular cross-section. One end of the bar is fixed, while other end is subjected to a pull force of 1000 N. Stress value is calculated analytically with the equation σ = F A , which gives a stress of 100 MPa. Using this stress value and Eq. (7), Creep strain was calculated analytically for a time period of 1 s and at three different temperatures of 800 °C, 850 °C, and 900 °C respectively. Temperature range were chosen keeping in mind the temperature dependency ofk2 . ANSYS Mechanical APDL model for same problem was developed with one element mesh. An axial 1D LINK180 element with one degree of freedom (DOF) per node was used for mesh and subsequent simulation. Creep Eq. (7) has been coded into Usercreep subroutine and linked with solver input file with “/upf ” command. Results of analytical approach and subroutine output is shown in Fig. 5. ANSYS creep strain values strongly agree with the analytical calculation at each temperature, thus subroutine is good to use for the subsequent analysis. A large strain jump from 800 °C to 850 °C can observed which is in accordance with observation made by Shewfelt et al.

Direct coupling between thermal and structural system is considered in two ways, and they are represented by Eq. (5) and (6) respectively. Equation (5) represents a strong coupling (e.g. Piezo-electricity problems) while Eq. (6) shows a vector coupling (e.g. Piezoresistivity problems). Direct thermal-structural problem can be solved using vector i.e. load coupling (ANSYS, 2016b). ut ⎡[M ] [0]⎤ ⎧ {u¨} ⎫ + ⎡ [C ] [0] ⎤ ⎧ {u}̇ ⎫ + ⎡[K ] [K ]⎤ ⎧ {u} ⎫ = ⎧ {F } ⎫ tu t ̇ ¨ ⎢ ⎥ {T } ⎬ ⎣ [0] [K t ] ⎥ ⎢ [0] [0]⎦ ⎥ ⎨ {T } ⎬ ⎢ ⎨ {Q} ⎬ ⎣ ⎭ ⎭ ⎩ ⎩ {T } ⎬ ⎦⎨ ⎩ ⎭ ⎣[C ] [C ]⎦ ⎨ ⎩ ⎭

(5)

⎡[M ] [0]⎤ ⎧ {u¨} ⎫ + ⎡[C ] [0] ⎤ ⎧ {u}̇ ⎫ + ⎡[K ] [O] ⎤ ⎧ {u} ⎫ t ⎥ t ⎥ {T } ⎬ ̇ ¨ ⎢ ⎢ ⎢ ⎣ [0] [0]⎥ ⎦⎨ ⎩ ⎭ ⎩ {T } ⎬ ⎭ ⎣ [0] [C ]⎦ ⎨ ⎩ {T } ⎬ ⎭ ⎣ [0] [K ]⎦ ⎨ {F } + {F th} ⎫ =⎧ ted ⎨ ⎩ {Q} + {Q } ⎬ ⎭

Temperature (°C)

(6)

In this problem, the contact between PT and CT is a transient event and the heat transfer due to thermal contact plays an important role for heat flow from PT to CT. The sequential coupling formulations (i.e. different mesh for different physics) cannot be employed for such transient problem and hence, direct coupling formulation (i.e. same mesh for all concerned physics) is used in this study. Moreover, load coupling is used in the direct formulation. 3.1. Material creep model As mentioned in the preceding section, thermally activated creep is the most dominating factor intruding the channel structural integrity. The postulated hypothesis for such severe accident leading to channel disassembly is the bending dominated event. This kind of deformation can be represented mathematically by a longitudinal strain rate equation and one such equation, as given in Eq. (7), is available in literature (Shewfelt et al., 1985). This equation is extensively used by large number of researchers (Gillespie et al., 1987; Mathew et al., 2003; Majumdar et al., 2014; Zhou et al., 2018) for the study of sagging deformation of PT.

Fig. 4. 1D Problem for creep subroutine. 141

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

locations. The PT was electrically heated to simulate resistive heating with the help of DC rectifier (42 kW). Experiments were performed till the contact between PT and CT was established; thereafter, heating was stopped. Multiple experiments were performed with heating power of 32.5 kW, 27.8 kW, and 15 kW respectively. To validate applicability of creep equation, an objective was set to find out time of first contact with ANSYS simulation and compare result with experimental time. If simulation is able to predict the time of first contact correctly then it can be assumed that the creep equation is sufficient to model high temperature sagging deformation in severe accident scenario. The simulation was carried out in two step; in first step temperature profile was developed and in the second step temperature profile was used to measure the deformation. ANSYS model consisted of the PT and CT as shown in Fig. 6(a). Heat generation as equivalent to heating power was applied to the PT. Radiosity solver was used to model the radiation between the PT and CT. No convective heat transfer was included inside annulus. Extra portion (250 mm on each side) outside tank was modelled with convective heat loss (h = 10 W/ m2.K) to the environment. Contact between the PT and CT was modelled using Conta174 and Targe170 surface elements and are laid out on each surface respectively. Creep strain rate as predicted by Eq. (7) is very small at lower temperature. As the temperature increases, it increases rapidly. A better agreement with experimental findings at higher temperature is expected (Zhou et al., 2018). From Fig. 7(c), it can be seen that default equation gives negligible sagging up to 800 °C, as anticipated. Hence, the parameter k2 for temperature up to 800 °C is set to be 0.019 (as used by Mathew et al., 2003 for their study) to alter the default creep initiation temperature. Even with this adjustment, it can be seen from Fig. 7(c) that the creep equation lags behind to predict the PT-CT contact as in the experiment. Hence, it is concluded that though creep equation is valid for predicting such deformation mechanism, it underpredicts displacement at lower temperature. One reason for this can be attributed towards development of different stress level because of difference in the methods of fabrication of PT and CT for the Indian reactors. Keeping this in the mind, it was then assumed to enhance the creep strain at lower temperature range by using a magnification factor to check whether an integer multiplication on creep equation can be sufficient to match the deformation of PT. A hit and trial method with

Fig. 5. Comparison of creep strain calculated by Mechanical APDL and Analytical approach.

3.3. Evaluation of creep model against experimental results The creep Eq. (7) was developed for CANDU reactor material. The Indian PT and CT are fabricated with a different manufacturing method, and it is observed that the longitudinal creep strain rate equation might under-predicts the deformation of the channel (Majumdar et al., 2014). Hence, it is very important to evaluate its effectiveness before it is being used in disassembly simulation analysis. The results from the experiments (Nandan et al., 2010) were used to evaluate the applicability of creep equation. These experiments were performed with the Indian PT and CT to study the events of LCDA at IIT R, India under the aegis of Bhabha Atomic Research Centre (BARC), India. Schematic of the experimental setup as well as the points for displacement and temperature measurements are represented in Fig. 6(a) while dimensions of critical parts are shown in Fig. 6(b). During sagging experiment, the CT was wrapped with ceramic wool to avoid any heat loss to environment by convection and radiation mode of heat transfer and air oxidation. Linear variable displacement transducer (LVDT) was used to measure the displacement of PT in vertical direction at different

Fig. 6. Details of experimental setup: (a) Schematic, (b) Important dimensions. 142

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

From the above validation for creep equation, it can be surmised that magnified creep equation can be used for channel disassembly simulation. It might over-predict the deformation between its operating temperature ranges. Designing the setup with this information is safe and is conservative. Henceforth, Creep Equation (7) with magnification factor was used in the analysis for an insight into possible progression of the channel deformation. 4. Coupled field modelling To investigate the channel disassembly in scope of understanding its deformation, a methodology is developed to study the thermo-structural response of the channel. Simulation model description along with the various aspects are detailed in the following section. Though, it is difficult to anticipate the experimental results beforehand, the results from such modelling can be viewed as a blind calculation for the prediction of probable deformation under the assumed accident hypothesis. The following assumptions are considered during the formulation of simulation model; 1. Depending on pressure inside the channel, the PT will come in contact with the CT, either as ballooned or sagged. In this paper, sagging mode of contact was considered. 2. Failure/Crack of the PT and CT was not considered in the present analysis. This simulation was done to investigate the deformation of channel at very high temperature. The corresponding situation in experiments would be the one that will be carried out in inert media. 3. Bundles were assumed to be remained intact during accident condition progression. 4.1. Model formulation Schematic of the model used for the analysis is shown in Fig. 8(a). Model setup consists of a PT and CT which represent the reactor channel. Dimensions of tubes are taken as that of the scaled down dimensions from Table 1. Only half portion was used in simulation to take the benefits of symmetry about y-z plane. Core length in Fig. 8(a) represents the channel length which resides inside the reactor core and this portion of channel is in direct contact with the moderator during normal operation of reactor. Fig. 9(a) to (c) show the material properties used in this simulation. Specific heat was kept same for both the PT and CT. Density and Poisson’s ratio are taken as 6520 kg/m3 and 0.35 for both the PT and CT respectively. Yield strength (Fig. 9(d)) are taken from the literature (Rodchenkov and Semenov, 2005). Direct-coupled field element Solid226, was used for meshing both PT and CT (shown in Fig. 8(b) for thermo-structure coupled effect. The Solid226 is a 3D 20-node hex (Brick) element, and it can handle the geometric non-linearity (i.e. large deformation), the material non-linearity (i.e. creep, plasticity etc.), and the surface-to-surface contact problems (ANSYS, 2016a). It has 4 degree of freedom (three translational dofs in x, y, z directions and one thermal dof) per node. Conta174 and corresponding Targe170 is used as a contact and target elements

Fig. 7. Comparison of Temperature and displacement measured by Mechanical APDL with Nandan Experiment (at 32.5 kW).

different integer value was run. It was found that if creep equation is multiplied by 6 (six) and temperature of the PT is kept same as in experiments, the simulation contact time matches with the experimental results within the acceptable limit. Simulated and experimental results (for 32.5 kW heating power) are shown in Fig. 7. It observed that though simulated deformation profile vary significantly from the experiment, the magnification factor do help to predict the first contact time between PT and CT correctly. To have a confidence over this magnification value, Creep equation with same integer multiplication was applied to other two sets of experiments (27.8 kW and 15 kW), and results are summarized in Table 3 for time of first contact between the PT and CT.

Table 3 Comparison of time of first contact for various heating powers. Cases

Ceramic wrapped at 32.5 kW (Nandan et al., 2010) Ceramic wrapped at 27.8 kW (Nandan, 2010) Ceramic wrapped at 15 kW (Nandan, 2010)

143

Contact time between PT/CT

% difference

Experimental

Simulation

280

274.43

1.989

170

165

2.94

400

392.5

1.875

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

Fig. 8. Description of simulation model: (a) Schematic, (b) Meshing of PT and CT.

Fig. 9. Thermal structural Properties of Channel: (a) Young’s modulus, (b) Thermal conductivity, (c) Specific heat, (d) Yield strength. Fig. 10. Thermal boundary conditions.

144

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

Fig. 11. Temperature distribution over the PT and CT at different longitudinal positions.

4.1.2. Thermal load boundary conditions Heat generation corresponding to effective power of 18 kW was applied to the PT, which could imitate the electrical heating. The thermal contact conductance depends on the texture of contacting surfaces, the contact pressure, surface roughness, and other parameters. Effect of these parameters can be included into ANSYS solver as a tabular data or as a user subroutine for contact elements. However, a constant value of 11,000 W/m2.K was assumed (Majumdar et al., 2014). Other boundary conditions, for possible modes of heat transfer that could be observed in the experiments, are illustrated in Fig. 10. Outer surfaces of extruded portion of PT was modelled with convective mode of heat transfer with a constant value as mentioned in the Fig. 10. During experiments, the CT ends will be water cooled, hence that cooling effect was modelled with convective mode of heat transfer with a considerably high value of convection coefficient. The outer surface of CT, inside the core, was modelled with both convective and radiative mode of heat transfer. Since, CT would be surrounded by Argon gas, and it is expected that there would not be much variation in convection heat transfer coefficient circumferentially over CT. Moreover, at higher temperature, the contribution of convective heat transfer will be less insofar as the radiative heat transfer. Therefore, a constant coefficient value was used over the outer surface of CT. Open enclosure condition was considered to model the surface-to-surface radiation heat transfer between outer surface of PT and inner surface of CT. Radiative mode of heat transfer was solved using radiosity solver which is an independent solver, and it is more robust than AUX-12 solver. It also facilitates the use of temperature dependent emissivity values. Due to non-linear nature of problem, direct sparse algorithm was used as Newton-Raphson solver for equilibrium iteration calculation while Gauss – Siedel iterative method was used as a radiosity solver. Total 30 circumferential divisions was used to mesh PT and CT. In order to reduce hourglass effect, through thicknesses two mesh layers were

respectively to model surface-to-surface interaction between the PT and CT. The formulation assumes the frictional contact (coefficient of friction, µ = 0.2) between the surfaces. Thermal-structural coupling was activated by setting keyopt 1 to 11 for Solid226 element. The keyopt settings in ANSYS are the mechanism to change the behaviour of elements. Load vector coupling was used between structural and thermal dofs through the proper keyopts setting. Time integration effect was not considered for the structural dofs, but it was used for the thermal dof.

4.1.1. Structural load boundary conditions Structural boundary conditions are shown in Fig. 8(a). One end of CT (z = 0) is fixed by constraining all three degree of freedom while other end’s movement is controlled by a non-linear spring. The spring had Stiffness of 2.E06 N/mm in the tension direction while it is fixed at 1.E03 N/mm in compression. Both ends of PT were fixed. Extended portions of the PT and CT on either sides of core length are imitating the experimental conditions that will be imposed (i.e. support for the CT and PT). Fuel bundle weight was assumed to be a distributed uniformly along inside bottom surface of PT. Force due to bundle weight was taken as 87.916 N considering symmetry assumption. Modified creep equation was used as per explanation in creep model evaluation section till temperature of 950 °C while beyond it, a transverse creep correlation (Shewfelt et al., 1984) as given in Eq. (8) was used. At higher temperature, creep in both longitudinal and transverse direction are assumed to be the same (Mathew et al., 2003).

ε ̇ = 10.4σ 3.4 exp(−19600/ T )

(8)

Effect of plasticity is expected to be less as compared to creep rate for such thermally driven deformation, hence to reduce the complexity in the model, it was not included in the this study. 145

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

Fig. 12. Circumferential temperature distribution at mid-section: (a) of PT, (b) of CT. Fig. 13. Temperature distribution along length: (a) of PT, (b) of CT.

used. Total elements and nodes stood to 27,600 and 1,55,970 respectively.

this region. Overall circumferential gradient due to heat conduction in circumferential direction is observed up to 65° and 42° for the PT and CT respectively at the end time of analysis. Thus, much lesser circumferential temperature gradient can be seen for such postulated accident event. Fig. 13 shows temperature distribution for PT and CT along respective lengths (detailed in each inset). Uniform temperature is seen over most of the core length at each timestamp. Temperature increases to 100 s, then it drops successively over the time. The drop in temperature after 100 s can be expected due to increase in surface contact area between the PT and CT which results in better heat leakage due to thermal contact conductance between mating surfaces. It may also be due to increase in radiative heat loss from the outer surface of CT. Sagging or vertical deformation of the channel is detailed in Fig. 14 at the point PT1 for multiple locations (370 mm, 570 mm, 770 mm, and 1070 mm respectively from the fixed end). Maximum deformation is observed at location Z:1070 mm which can be inferred based on assumed boundary conditions. The inset image in Fig. 14(a) shows deformation in more detail for initial 40 s of the analysis. A 3 mm (Initial gap between the PT and CT) deflection of point PT1 is achieved at the location Z:370 mm which is close to supports in 20 s from the start of simulation. It indicates that around 60% of PT core length got deformed to 3 mm within 20 s. Fig. 14(b) shows channel deformation of CT at mid-span. Initial deformation of CT is insignificant as compared to PT owing to no additional weight except its self-weight. It is observed that the difference between deformation of PT and CT at Z:1070 mm is about 3 mm (Initial gap). Large deformation of CT is produced after 25 s due to successive transfer of weight of PT and fuel bundles as they

4.2. Results and discussion Circumferential temperature distributions of the PT and CT at different axial locations are shown in Fig. 11. The locations represented in inset are only for representation and are not to the scale. A linear temperature rise at the rate of 30 °C/s is observed for the point PT3 at all locations for the first 25 s of simulation. Temperature difference between the points PT2 and PT3 and similarly between CT2 and CT3 was observed to be very small. This can be due to the assumptions of a uniform convective heat transfer coefficient over the outer surface of CT, and heat conduction in the circumferential direction can only be taken place through the contact area between the PT and CT. The departure in the temperature trend of point CT1 from that of points CT2 and CT3 at location Z:370 mm is seen at 25 s from the start of simulation. Thus, within 25 s contact has progressed up to the length of 700 mm on either side of mid-section (i.e. the location Z:1070 mm). Temperature at points PT1 and CT1 is seen to converge at approximately same value after 100 s on all locations. The circumferential transient temperature distribution at Z:1070 mm axial position is shown in Fig. 12 for the both PT and CT. The angles 0° and 180° represent the bottom and top point respectively for both tubes. In upper region of tubes, circumferential temperature gradient is observed to be less, again could be due to simplified approach for the convection heat transfer boundary condition. A soot up temperature up to angle of 24° in circumferential direction is seen at 50 s for CT and correspondingly drop-in is observed for PT. It can be inferred that surface contact between PT and CT could take place over 146

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

Fig. 15. Equivalent (Von-Mises) Stress along length: (a) of PT, (b) of CT.

Fig. 14. Sagging deformation of channel: (a) of PT1 at different locations, (b) of CT at mid-span.

can be deduced that the failure length can correspond to lengths of L1, L2, L3, and L4 with the transient pass as shown in Fig. 16(b). Length L4 could be the most likely the length of disassembly as the CT thickness reduction is nearly 10% of its initial thickness. Increase in the thickness towards very ends of channel may be resulted because of the imposed frictionless boundary conditions which allowed the movement in all direction but normal to applied surface. Load transfer due to contact between the PT and CT is reflected in Fig. 17. It is seen that load of approximately 83 N i.e. 94.41% of total fuel bundle weight, was transferred to the CT within 40 s of the simulation. However, after that it decreases. Since, contacting load on the contact element is calculated as a reaction force in ANSYS. The Plot can be envisaged as the CT supported the loads until 40 s due to its elastic stiffness giving load output which can also be seen through equivalent stress plot in Fig. 15(b). After that the CT deforms with PT therefore results in a relaxation of reaction force. The power in reactor core varies among the channels with maximum power channels are concentrated in centre region with surrounding rest channels are at lower powers (Bajaj and Gore, 2006). Hence, disassembly rate will differ throughout the core during accident transient. Fig. 18 shows the effect of heat-up rate variation on the channel sagging and on the load transfer from PT to CT. The same trend in the channel sagging is observed for each heatup case but maximum deformation is observed to be less for lower heat-up power; thus, elucidate the nature of thermally dependent creep. Large deformation due to sagging is seen to start from the more or less at same time at around 40 s from start of transient and is consistent with contacting load from

deform. The Gradual deformation of PT due to creep and fuel bundles weight, lead to progressive deformation of CT as weight gets transferred from PT. It can be concluded that the bulk deformation of the channel will only be predicted once most length of PT rests on CT. Fig. 15 shows an equivalent stress (Von-Mises) plot for the both PT and CT along respective length of contact. A stress of 22 MPa at the very start of simulation and at the mid-section of PT, represents the initial bending stress while a spike in stress (to 12 MPa) for CT mainly arises due to its contact with PT. Large stress value (at 240 mm and 1920 mm) on either end of PT, may result from end constraint effect. Stress level throughout the analysis is observed to be less than yield stress for the both PT and CT material (as referenced from Fig. 9(d)); hence, our assumptions for precluding the plasticity model from the simulation is here justified. Moreover, it establishes the creep as dominating model for disassembly or sagging of the channel under assumed accident hypothesis. The change in thicknesses of PT and CT due to creep dominated material flow in elucidated in Fig. 16 for the end time (i.e. 500 s). Largest reduction in thickness is observed toward the ends of each tube and can be contributed due to elongation of channel as it sags. Disassembly will result as CT loses structural integrity as its thickness reduces significantly. More gradual reduction in thickness is seen for PT while for CT step change in thickness is prevailed. To avoid more complexities in the numerical solver, failure model was not included to anticipate the crack location. However, the possible crack initiation can be guessed through changes in channel thickness. Based on this note, it 147

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

Fig. 16. Channel thickness along length at end time: (a) of PT, (b) of CT.

Fig. 18. Effect of heat-up rate: (a) On channel sagging, (b) On contact load transfer.

m2.K) are shown in Fig. 19. It is clear that the TCC has no or negligible effect either on the channel displacement or on the contact load transfer. Possible reason for such observation can be elaborated as reason of the lesser contact angle between the PT and CT (24° as seen through Fig. 12 for base case) thus producing lesser contacting area. And, heat flow in circumferential direction is observed to be far less as explained in base case results. Hence, effect of heat transfer due to TCC is restricted to very small sector of the PT and CT while through the rest of the region radiative mode is the dominated mode of heat transfer.

5. Conclusion A postulated accident, under low pressure event and completely voided condition, progresses to severe phase, a very first anticipated event is failure of channel integrity. It becomes very important to know the nature of event in order to provide any mitigating countermeasures. But, it is nearly impossible to perform a full scale channel disassembly experiments in laboratory owing to many inherent difficulties. A scale down methodology is adopted to process the experiments. This paper points out certain criteria related to geometry, mechanical properties, and applied load on system that are needed to be scoured for scaled down experiments. This reduced scale sagging experiment will help to understand more about the accident event of disassembly progression for the channels. Furthermore, a numerical model was performed using direct coupled field approach with the help of ANSYS workbench mechanical solver. The results from this numerical investigation helped to

Fig. 17. Load transfer through contact.

Fig. 18(b). Thus, load transfer from PT to CT also play an important role in sagging of channel as whole, and time for load getting transfer is seen to be roughly same. Thermal contact conductance (TCC) plays a significant role in dissipating heat rate from PT to CT after contact. To check the effect of TCC as a sensitivity parameter on the channel sagging and correspondingly on the contact load transfer, simulations were performed by changing the TCC to a random lowest and largest value. The results of simulations and its comparison with base case (case with 11,000 W/

148

Nuclear Engineering and Design 339 (2018) 137–149

A.R. Singh et al.

dedicated lab setup and also thankful to Institute computer Centre (IIC, IIT Roorkee) for giving an access to the simulation facility. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.nucengdes.2018.09.008. References Ansys ANSYS Mechanical Application User ’s 2016 Guide. ANSYS Inc. Ansys Workbench Mechanical APDL 2016 ANSYS Inc. Bajaj, S.S., Gore, A.R., 2006. The Indian PHWR. Nucl. Eng. Des. 236, 701–722. https:// doi.org/10.1016/j.nucengdes.2005.09.028. Gillespie, G.E., Moyer, R.G., Hadaller, G.I., 1987. An experimental investigation of the creep sag of pressure tubes under LOCA conditions. In: Presented at the 8th Canadian Nuclear Society, pp. 68–72. Gupta, S.K., Dutt, B.K., Raj, V.V., Kakokar, A., 1996. A Study of The Indian PHWR Reactor Channel Under Prolonged Deteriorated Flow Conditions. Presented at the Proc. IAEA TCM on Advances in Heavy Water Reactors. IAEA, 2008. Analysis of Severe Accidents in Pressurized Heavy Water Reactors (No. June). International Atomic Energy Agency, Vienna, Austria. Majumdar, P., Chatterjee, B., Lele, H.G., Guillard, G., Fichot, F., 2014. ASTEC adaptation for PHWR limited core damage accident analysis. Nucl. Eng. Des. 272, 273–286. https://doi.org/10.1016/j.nucengdes.2013.10.011. Mathew, P.M., Kupferschmidt, W.C.H., Snell, V.G., Bonechi, M., 2001. CANDU-Specific Severe Core Damage Accident Experiments in Support of Level 2 PSA. In: Presented at the 16th International Conference on Structural Mechanics in Reactor Technology (SMiRT-16), Washington DC, USA, pp. 1–8. Mathew, P.M., Nitheanandan, T., Bushby, S.J., 2008. Severe Core Damage Accident Progression within a CANDU 6 Calandria Vessel. Presented at the European Review Meeting on Severe Accident Research (ERMSAR-2008). AECL Chalk River Laboratories Reactor Safety Division, Nesseber, Bulgaria. Mathew, P.M., White, A.J., Snell, V.G., Bonechi, M., 2003. Severe Core Damage Experiments and Analysis for CANDU Applications. In: Presented at the 17th International Conference on Structural Mechanics in Reactor Technology (SMiRT17). Czech Republic, Prague, pp. 1–8. Mukhopadhyay, D., Majumdar, P., Behera, G., Gupta, S.K., Raj, V.V., 2002. Thermal analysis of severe channel damage caused by a stagnation channel break in a PHWR. J. Press. Vessel Technol. 124, 161. https://doi.org/10.1115/1.1463036. Nandan, G., 2010. Study of sagging and ballooning behaviour of Indian PHWR Pressure tube during LOCA. PhD Thesis. http://hodhbhagirathi.iitr.ac.in:8081/jspui/handle/ 123456789/6689. Nandan, G., Sahoo, P.K., Kumar, R., Chatterjee, B., Mukhopadhyay, D., Lele, H.G., 2010. Experimental investigation of sagging of a completely voided pressure tube of Indian PHWR under heatup condition. Nucl. Eng. Des. 240, 3504–3512. https://doi.org/10. 1016/j.nucengdes.2010.05.042. Negi, S., Kumar, R., Majumdar, P., Mukopadhyay, D., 2017. Full length channel pressure tube sagging study under postulated LOCA with un-availability of ECCS in an Indian PHWR. Nucl. Eng. Des. 320, 361–373. https://doi.org/10.1016/j.nucengdes.2017. 06.017. Nitheanandan, T., Brown, M.J., 2013. Backup and ultimate heat sinks in CANDU reactors for prolonged SBO accidents. Nucl. Eng. Technol. 45, 589–596. https://doi.org/10. 5516/NET.03.2013.711. Pan, B., Wu, D., Wang, Z., Xia, Y., 2011. High-temperature digital image correlation method for full-field deformation measurement at 1200 °C. Meas. Sci. Technol. 22. https://doi.org/10.1088/0957-0233/22/1/015701. Rodchenkov, B.S., Semenov, A.N., 2005. High temperature mechanical behavior of Zr2.5% Nb alloy. Nucl. Eng. Des. 235, 2009–2018. https://doi.org/10.1016/j. nucengdes.2005.05.032. Shewfelt, R.S.W., Lyall, L.W., Godin, D.P., 1984. A high-temperature creep model for Zr2.5 wt% Nb pressure tubes. J. Nucl. Mater. 125, 228–235. https://doi.org/10.1016/ 0022-3115(84)90548-8. Shewfelt, R.S.W., Nal, L., Lyall, L.W., 1985. A High Temp Longitudinal strain rate equation for Zr-2.5wt% Nb Pressure Tubes. J. Nucl. Mater. 132, 41–46. https://doi. org/10.1016/0022-3115(85)90391-5. Singh, R.J., Ravi, K., Gupta, S.K., 2011. Methodology for developing channel disassembly criteria under severe accident conditions for PHWRs. Ann. Nucl. Energy 38, 1884–1890. https://doi.org/10.1016/j.anucene.2011.05.011. Thomson, P.D., Kohn, E., 1983. Fuel And Fuel Channel Behaviour In Loss Of Coolant without availability of ECCS. Presented at the Specialist Metting on Water Reactor Fuel Safety and Fission Product Release in Off-Normal and Accident Conditions, IWGFPT/16. Riso National Laboratory, Denmark. Zhou, F., Novog, D.R., Siefken, L.J., Allison, C.M., 2018. Development and Benchmarking of Mechanistic Channel Deformation Models in RELAP/SCDAPSIM/MOD3.6 for CANDU Severe Accident Analysis. Nucl. Sci. Eng. 190, 209–237. https://doi.org/10. 1080/00295639.2018.1442060.

Fig. 19. Effect of Thermal contact conduction (TCC): (a) On channel sagging, (b) On contact load transfer.

understand the accident phenomena to a certain extent. Presently, the available longitudinal creep equation is not sufficient to predict the deformation at lower temperature in context to Indian PHWR reactor channel; thus, it was modified using a multiplying factor. It was then validated with the available experimental results to estimate the first contact time for PT with the CT. Conclusions from the analysis performed under inert media/environment is summarised below: 1. Much lesser circumferential temperature gradient can be seen for such postulated accident event. 2. Large deflection or sag of the channel can be observed once more than a half portion of the PT contacts on the CT. 3. Large sag of channels with different heat-up power can be seen to start from more or less at the same time. 4. Creep is observed as a dominant mechanism for the channel sagging. Standalone plasticity model might not be sufficient to model the channel sag. 5. Thermal contact conductance has negligible effect on the channel sag and contact load transfer. Acknowledgments Authors are grateful to IIT Roorkee, India for providing the space for

149