ML16053A094: Difference between revisions

From kanterella
Jump to navigation Jump to search
(Created page by program invented by StriderTol)
(Created page by program invented by StriderTol)
Line 35: Line 35:
Rayleigh (Ra) numbers. These values, combined with Nusselt (Nu) correlations, lead to theheat transfer coefficients which are used to find the energy transferred via convection.a. The Prandtl number is a measure of the fluid's kinematic diffusivity (v) to thermaldiffusivity (a) of the fluid:Pr =--,IaWhere thermal diffusivity of the fluid is calculated as:kcp
Rayleigh (Ra) numbers. These values, combined with Nusselt (Nu) correlations, lead to theheat transfer coefficients which are used to find the energy transferred via convection.a. The Prandtl number is a measure of the fluid's kinematic diffusivity (v) to thermaldiffusivity (a) of the fluid:Pr =--,IaWhere thermal diffusivity of the fluid is calculated as:kcp
* p(9S6)b. The Grashof number (and implicitly, values derived from the Grashof number) aredependent on both surface temperature and channel temperature. The Grashofnumber (for natural circulation heat transfer) is defined as:g./ (T- Tier) *(10)Where g is the acceleration due to gravity, /3 is the thermal expansioncoefficient of the coolant, s is the channel width, and/p is the dynamic viscosityof the coolant.c. The Rayleigh number is calculated [Kaminski4],:Ra =gIsr= fp2(Twaii -Tflu~d)s (1Ra Pr (1Where s is the total channel width.d. The Nusselt number relates the conductive and convective heat transfer effects ofthe fluid. The heat transfer is driven by the temperature difference and is foundthrough the Nusselt number. A correlation for natural circulation in verticalchannels for the Nusselt number4is:-1/2Nu3-.(12)e. The average heat transfer coefficient, h,, can be found from the average Nusseltnumber using:4D. A. Kaminski,, M. K. Jensen, "Introduction to Thermal and Fluids Engineering," John Wiley & Sons (2005)  
* p(9S6)b. The Grashof number (and implicitly, values derived from the Grashof number) aredependent on both surface temperature and channel temperature. The Grashofnumber (for natural circulation heat transfer) is defined as:g./ (T- Tier) *(10)Where g is the acceleration due to gravity, /3 is the thermal expansioncoefficient of the coolant, s is the channel width, and/p is the dynamic viscosityof the coolant.c. The Rayleigh number is calculated [Kaminski4],:Ra =gIsr= fp2(Twaii -Tflu~d)s (1Ra Pr (1Where s is the total channel width.d. The Nusselt number relates the conductive and convective heat transfer effects ofthe fluid. The heat transfer is driven by the temperature difference and is foundthrough the Nusselt number. A correlation for natural circulation in verticalchannels for the Nusselt number4is:-1/2Nu3-.(12)e. The average heat transfer coefficient, h,, can be found from the average Nusseltnumber using:4D. A. Kaminski,, M. K. Jensen, "Introduction to Thermal and Fluids Engineering," John Wiley & Sons (2005)  
--- 'k (13)4 The UT LOCA ModelThe University of Texas Loss of Coolant model is a combination of finite element analysis (FEA) forsteady state, and transient fuel conditions, as well as an air channel analysis sub-section to provideeffective estimation of air channel heating. This channel sub-model leads to proper parametricvariation analysis by giving a real world upper bounding temperature.4.1 Coolant Air TemperatureIn order to find the limiting values of the channel air temperature, a separate, one dimensionalvertical model was created independent of the FEA model and geometry used in finding the fueltemperatures. It utilizes an elemental, vertical, constant temperature surface interfacing withbuoyant air.The temperature entering the bottom of the channel surrounding the fuel element is thelimiting room air temperature following a loss of coolant, 20°C. The rise in temperature is foundby segmenting the pin vertically. Each iteration has a specific heat flux relative to thetemperature difference between the surface and the air, its specific dimensionless parameters,and a constant surface temperature that is user defined. The limiting conditions set the surfacetemperature at 950G. This was the value used later as a limiting factor in FEA parametricvariation.The change in air temperature across each segment is a function of the heat generated in thesegment and the heat transfer coefficient calculated from local non-dimensional parameters.Heat transfer characteristics in convection depend on intrinsic and extrinsic material propertiesand fluid temperature, with the heat transfer coefficient calculable though the use ofdimensionless numbers. The temperature rise of the fluid entering the region of heat transferfor subsequent segments is the exit temperature for the preceding segment, i.e., the rise intemperature across the previous segment added to the temperature of the coolant entering theprevious segment. The channel flow heat up model provided an order of magnitude estimationleading to proper parametric variation.4.1.1 Fluid flow, and thus the characteristic velocity, is driven by natural convection and isdominated by the buoyancy driven numbers in the Rayleigh number (Rag), the productof the Grashof and Prandtl numbers.a. The change in coolant temperature from fluid flow across a segment begins byfinding the appropriate Rayleigh number (eqn. 11) for the ith~ segment [Kaminski6],sT. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. DeWitt op cit & C. 0. Popiel and J. Wojtkowiak, "Simpleformulas for thermophysical properties of liquid water for heat transfer calculations (from 0 to 150 degrees C) (vol19, pg 87, 1998)," Heat Transf. Eng., vol. 19, no. 3, pp. 87-101, 1998.SD. A. Kaminski,, M. K. Jensen, "Introduction to Thermal and Fluids Engineering," John Wiley & Sons (2005)then the segment's Nusselt number (eqn. 12), then the heat transfer coefficient (eqn.13). With the heat transfer coefficient now found, the heat flux (q') is found using:q"= i (Ts -Tm~f, ) (15)Where TIs is the cladding surface temperature and Tinli is the heat sink temperature.b. Heat flux is used to find the Modified Rayleigh number. For uniform wall heating, themodified Rayleigh (Ra*)3 is:R =g "f3."q'" p2. c. s4Ra k2 (16)Where g is the acceleration due to gravity, fi is the thermal expansion coefficient,and/p is dynamic viscosity.c. The modified Rayleigh number leads to the characteristic channel velocity7 (Uz):Uz = (17)d. The change in temperature for fluid flow across a segment of the fuel element alongthe (axial) direction of flow can be calculated with:Q, = rh .cp. AT (18)e. Where the A~is calculated as:".AFE,AT = AlwU cP(19)Where Uz is calculated from eqn. (19). This AT is added to the segment's inlettemperature and becomes the inlet temperature for the next segment. The lastsegment's channel temperature represents the culmination of all the heating.4.1.2 As an independent calculation to determine limiting values of air temperature, thetemperature rise was found through standard gas laws.a. The characteristic velocity gives a stay time for the air (heated length divided bycharacteristic channel velocity). This allows the change in energy to be calculated asfollows:dE = q" " A5s* tstay (20)7K. Vafai, C. P. Desai, S. V. Iyer, and M. P. Dyko, "Buoyancy Induced Convection in a Narrow Open-Ended Annulus,"J. Heat Transfer, vol. 119, p. 483, 1997.
--- 'k (13)4 The UT LOCA ModelThe University of Texas Loss of Coolant model is a combination of finite element analysis (FEA) forsteady state, and transient fuel conditions, as well as an air channel analysis sub-section to provideeffective estimation of air channel heating. This channel sub-model leads to proper parametricvariation analysis by giving a real world upper bounding temperature.4.1 Coolant Air TemperatureIn order to find the limiting values of the channel air temperature, a separate, one dimensionalvertical model was created independent of the FEA model and geometry used in finding the fueltemperatures. It utilizes an elemental, vertical, constant temperature surface interfacing withbuoyant air.The temperature entering the bottom of the channel surrounding the fuel element is thelimiting room air temperature following a loss of coolant, 20°C. The rise in temperature is foundby segmenting the pin vertically. Each iteration has a specific heat flux relative to thetemperature difference between the surface and the air, its specific dimensionless parameters,and a constant surface temperature that is user defined. The limiting conditions set the surfacetemperature at 950G. This was the value used later as a limiting factor in FEA parametricvariation.The change in air temperature across each segment is a function of the heat generated in thesegment and the heat transfer coefficient calculated from local non-dimensional parameters.Heat transfer characteristics in convection depend on intrinsic and extrinsic material propertiesand fluid temperature, with the heat transfer coefficient calculable though the use ofdimensionless numbers. The temperature rise of the fluid entering the region of heat transferfor subsequent segments is the exit temperature for the preceding segment, i.e., the rise intemperature across the previous segment added to the temperature of the coolant entering theprevious segment. The channel flow heat up model provided an order of magnitude estimationleading to proper parametric variation.4.1.1 Fluid flow, and thus the characteristic velocity, is driven by natural convection and isdominated by the buoyancy driven numbers in the Rayleigh number (Rag), the productof the Grashof and Prandtl numbers.a. The change in coolant temperature from fluid flow across a segment begins byfinding the appropriate Rayleigh number (eqn. 11) for the ith~ segment [Kaminski6],sT. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. DeWitt op cit & C. 0. Popiel and J. Wojtkowiak, "Simpleformulas for thermophysical properties of liquid water for heat transfer calculations (from 0 to 150 degrees C) (vol19, pg 87, 1998)," Heat Transf. Eng., vol. 19, no. 3, pp. 87-101, 1998.SD. A. Kaminski,, M. K. Jensen, "Introduction to Thermal and Fluids Engineering," John Wiley & Sons (2005) then the segment's Nusselt number (eqn. 12), then the heat transfer coefficient (eqn.13). With the heat transfer coefficient now found, the heat flux (q') is found using:q"= i (Ts -Tm~f, ) (15)Where TIs is the cladding surface temperature and Tinli is the heat sink temperature.b. Heat flux is used to find the Modified Rayleigh number. For uniform wall heating, themodified Rayleigh (Ra*)3 is:R =g "f3."q'" p2. c. s4Ra k2 (16)Where g is the acceleration due to gravity, fi is the thermal expansion coefficient,and/p is dynamic viscosity.c. The modified Rayleigh number leads to the characteristic channel velocity7 (Uz):Uz = (17)d. The change in temperature for fluid flow across a segment of the fuel element alongthe (axial) direction of flow can be calculated with:Q, = rh .cp. AT (18)e. Where the A~is calculated as:".AFE,AT = AlwU cP(19)Where Uz is calculated from eqn. (19). This AT is added to the segment's inlettemperature and becomes the inlet temperature for the next segment. The lastsegment's channel temperature represents the culmination of all the heating.4.1.2 As an independent calculation to determine limiting values of air temperature, thetemperature rise was found through standard gas laws.a. The characteristic velocity gives a stay time for the air (heated length divided bycharacteristic channel velocity). This allows the change in energy to be calculated asfollows:dE = q" " A5s* tstay (20)7K. Vafai, C. P. Desai, S. V. Iyer, and M. P. Dyko, "Buoyancy Induced Convection in a Narrow Open-Ended Annulus,"J. Heat Transfer, vol. 119, p. 483, 1997.
Where tstay is the time the cooling air is in contact with the cladding surface.b. By using the density of air and the volume of the channel, the mass of the air in thespace at any given time can be found, by neglecting density changes. Using theequation below, the change in temperature can be found:q =mrcyAT (21)AT = q- =Tilit -+/-AT4.1.3 Results of calculations for limiting values of the channel air temperatureThese two methods routinely agreed across variations in surface temperatures, with thelimiting channel temperature, of 20°C inlet and 950°C surface temperature, being35.3°C.4.2 Finite Element Model Geometry and BasisThe calculation of temperature distribution in the fuel element is accomplished by using theprinciples of finite element analysis. The fuel element geometry is based on a cylindricalsegment. The axial height of the segment is the total heated length (0.381 m) divided by thenumber of segments (15). Radial dimensions are taken from General Atomics drawings.
Where tstay is the time the cooling air is in contact with the cladding surface.b. By using the density of air and the volume of the channel, the mass of the air in thespace at any given time can be found, by neglecting density changes. Using theequation below, the change in temperature can be found:q =mrcyAT (21)AT = q- =Tilit -+/-AT4.1.3 Results of calculations for limiting values of the channel air temperatureThese two methods routinely agreed across variations in surface temperatures, with thelimiting channel temperature, of 20°C inlet and 950°C surface temperature, being35.3°C.4.2 Finite Element Model Geometry and BasisThe calculation of temperature distribution in the fuel element is accomplished by using theprinciples of finite element analysis. The fuel element geometry is based on a cylindricalsegment. The axial height of the segment is the total heated length (0.381 m) divided by thenumber of segments (15). Radial dimensions are taken from General Atomics drawings.
dyINiFigure 1. Finite Element Radial GeometryThe Finite Element Model radii used in computation was selected based on both parametervalidation and computational power available. The limiting geometric figure of concern is theBiot number, which relates convective and conductive aspects of the element to its volume tosurface area ratio. It is determined using the equation below:Bi = h'Lc(2)Where, the characteristic length, 4c, is defined as the volume to surface area ratio:VL- (22b)Differential radii in the outer portions of themodel were chosen to most accurately subdividethe real geometry of the cladding and the gas gap. Internal fuel differential radii were chosen tominimize the Biot number. In addition to the Biot number, the Fourier number is a transientfigure of merit related to constants that determine time response and the geometry:a'tFo = -(23According to Bergman4, the Blot number must remain below 0.1, and the Fourier number mustremain below 0.5 for lumped parameter analysis to be valid. This was the merit to which thedifferential radii are chosen.
dyINiFigure 1. Finite Element Radial GeometryThe Finite Element Model radii used in computation was selected based on both parametervalidation and computational power available. The limiting geometric figure of concern is theBiot number, which relates convective and conductive aspects of the element to its volume tosurface area ratio. It is determined using the equation below:Bi = h'Lc(2)Where, the characteristic length, 4c, is defined as the volume to surface area ratio:VL- (22b)Differential radii in the outer portions of themodel were chosen to most accurately subdividethe real geometry of the cladding and the gas gap. Internal fuel differential radii were chosen tominimize the Biot number. In addition to the Biot number, the Fourier number is a transientfigure of merit related to constants that determine time response and the geometry:a'tFo = -(23According to Bergman4, the Blot number must remain below 0.1, and the Fourier number mustremain below 0.5 for lumped parameter analysis to be valid. This was the merit to which thedifferential radii are chosen.

Revision as of 02:39, 29 March 2018

University of Texas at Austin - Response to Request for Additional Information Regarding the License Renewal Request for the Nuclear Engineering Teaching Laboratory Triga Mark Ii Nuclear Research Reactor (TAC No. ME7694)
ML16053A094
Person / Time
Site: University of Texas at Austin
Issue date: 02/05/2016
From: Whaley P M
University of Texas at Austin
To: Balazik M F
Document Control Desk, Office of Nuclear Reactor Regulation
References
TAC ME7694
Download: ML16053A094 (25)


Text

Deparzmcnt of Me~chanical EngineeringTHE UNIVERSITY OF TEXAS AT AUSTINNuclear Engineering Teaching Laboratory: Austin, Texas 78758512-232-53 70 " FAX.512-471,-4589 herp :f/vww. me u texas, edu/- net, lFebruary 5, 2016ATTN: Document Control Desk,U.S. Nuclear Regulatory Commission,Washington, DC 20555-0001M. BalazikProject ManagerResearch and Test Reactors Licensing BranchSUBJECT: Docket No. 50-602, Request for Renewal of Facility Operating License R-229REF: UNIVERSITY OF TEXAS AT AUSTIN -REQUEST FOR ADDITIONAL INFORMATION REGARDING THE LICENSERENEWAL REQUEST FOR THE NUCLEAR ENGINEERING TEACHING LABORATORY TRIGA MARK II NUCLEARRESEARCH REACTOR (TAC NO. ME7694)Sir:Analysis has been completed on the two remaining technical issues, including a review of reactivity parametersbased on burnup calculations and a loss of coolant accident analysis. An updated validation of reactivity isattached along with a new LOCA thermal hydraulic analysis.Please contact me by phone at 512-232-5373 or email whaley(Wmail.utexas.edu if you require additionalinformation or there is a problem with this submittal.P. M. WhaleyAssociate DirectorNuclear Engineering Teaching LaboratoryThe University of Texas at AustinI declare under penalty of perjury that the foregoing is true and correct.Executed on February 5, 2016Steven R. BiegalskiNETL DirectorOg&

RQUEST FOR ADDITONAL INFORMATION (July 31, 2015: Adams Accession ML15211A362):The guidance in NUREG-1537 Section 4.5.2, "Reactor Core Physics Parameters," requests theapplicant provide calculations of certain core physics parameters and compare them withapplicable measurements. A request for additional information was sent to you in a letter datedJuly 25, 2012 (ADAMS Accession No. ML15211A638), regarding a comparison of calculated andmeasured values for reactivity parameters. In response, your submittal dated July 15, 2015(ADAMS Accession No. ML121500308), provided calculations for control rod worth and excessreactivity in Tables 11 and 12, respectively. In addition, you provided comparisons of these\calculations and measurements in Table 13.a. In Table 11, your calculations of excess reactivity and control rod worth on 3/16/ 1992,7/24/2007, 6/4/2008, 6/11/2008, 6/14/2010, 6/23/2010, 7/25/2011, and 8/2/2012indicate the shutdown margin for the UT TRIGA reactor would have a positive reactivity.A positive reactivity would not meet Technical Specification 3.2 Shutdown Margin.Please justify why the positive reactivity presents an acceptable level of safety regardingshutdown margin for the UT reactor.b. In Tables 11 and 12, your calculations for excess reactivity are consistently higher thanthe measured values. The average bias between the calculated and measured values is$2.56. Please justify why this bias presents an acceptable level of agreement inmodeling the UT reactor.c. Control rod worth in Table 13 appears to be inconsistent. For example, the RegulatingRod worth varies between +12.7% and -21.2%. Specifically, the values for two separateRegulating Rod data points on 7/13/2012, changes from -10.1% to +12.7%. Similarly,Shim-2 worth varies between -23.1% and +30.1%. Please justify why these variationsindicate an acceptable level of agreement in modeling the UT reactor.

RESPONSE:Calibrating model data to measured excess reactivity values indicates all shutdown margins arenegative by a wide margin.Most of the fuel elements in the current UT reactor had been partially burned at other facilitiesprior to use in this facility. The amount of uranium remaining in an element is based on totalcore burnup, distributed across the elements in the core. TRIGA peaking factors suggest asmuch as 40% variation in neutron flux across the core from center to peripheral elementsduring operation. Investigation of the sensitivity of reactivity calculations to uranium 235 massindicates a 10% change in mass has on the order of $5 impact on excess reactivity. Given thepotential uncertainty in the fuel mass in the partially burned elements, excess reactivity maynot be a reasonable function for validating modeling. When reactivity values are adjustedbased on calculations assume nominal uranium 235 and 90% of the nominal values and onmeasured excess reactivity values, integral control rods worth agree to measured values to areasonable level.Previous comparisons were based solely on burnup, and did not consider core configurations.As a result, reactivity values for different core configurations were identified by burnup datethat did not correspond in all cases to actual core configuration. In considering coreconfigurations and core burnup values, there is general agreement.ANALYSISAs previously noted, the initial UT TRIGA core was principally composed of previously (lightlyburned) irradiated fuel elements. Special Nuclear Material records were used as a basis for theoriginal uranium composition in material calculations. However, burnup calculations for theNuclear Materials Management and Safeguards System (NM MSS) track total facility inventorybased on core burnup, as opposed to individual fuel elements. Burnup in individual fuelelements varies significantly from average core burnup, and exacerbate potential uncertainty incalculating element-specific burnup.To evaluate the sensitivity of reactivity values to varying uranium 235 content, calculationswere performed for fuel composition using the nominal values from special nuclear materialrecords, and then with 90% of the uranium 235 specifications in the records. The differencebetween the two values was calculated. The current core contains 114 elements, andadditional data is provided for two burn intervals. The deviations in excess reactivity are muchhigher than the difference in individual integral control rods worth.

Table 1, Reactivity Based on 100% SNM 235 Material Composition and Difference at 90%EXCESS A RR A SHi A SH2 A TR A90 INITIAL $7.21 -$5.31 $4.33 $0.26 $4.30 -$0.16 $2.53 -$0.03 $2.37 $0.08FINAL $7.09 -$5.18 $4.41 $0.42 $3.94 $0.22 $2.38 $0.39 $2.38 $0.1389 INITIAL $8.12 -$5.05 $3.63 $0.56 $3.07 $0.74 $2.66 $0.35 $2.66 $0.26FINAL $7.69 -$5.07 $3.65 $0.46 $2.96 $0.24 $2.61 $0.61 $2.53 $0.1492 INITIAL $6.91 -$5.16 $4.55 $0.50 $3.08 $0.90 $2.35 $0.45 $2.05 $0.28FINAL $6.54 -$4.79 $4.07 $0.72 $3.09 $0.69 $2.06 $0.84 $1.67 $0.3995 INITIAL $7.44 -$5.09 $4.27 $0.04 $3.16 $0.40 $2.69 $0.21 $2.41 $0.05FINAL $6.32 -$5.36 $4.04 $0.28 $3.18 $0.00 $2.72 -$0.07 $2.24 $0.04103 INITIAL $8.04 -$5.69 $3.89 $0.31 $2.98 $0.35 $2.48 $0.20 $2.04 $0.27FINAL $7.71 -$3.59 $4.05 $0.07 $3.24 $0.12 $2.59 -$0.12 $2.19 -$0.06102 INITIAL $5.67 -$5.58 $2.69 $0.07 $3.28 $0.33 $2.26 $0.26 $53.21 -$0.08FINAL $4.87 -$4.97 $2.54 $0.81 $2.13 $0.46 $2.43 $0.39 $2.96 $0.30104 INITIAL $7.41 -$6.71 $3.32 $0.20 $2.77 $0.30 $3.00 $0.00 $2.76 $0.05FINAL $6.19 -$5.49 $3.43 $0.09 $3.20 -$0.13 $2.86 $0.14 $2.71 $0.06108 INITIAL $6.85 -$6.29 $3.09 -$0.22 $2.36 -$0.13 $2.68 -$0.16 $2.63 -$0.25FINAL $6.96 -$5.45 $3.07 $0.42 $2.51 $0.34 $2.98 $0.16 $2.66 -$0.03110 INITIAL $6.96 -$5.45 $4.24 -$0.09 $2.35 -$0.07 $3.03 -$0.07 $1.66 -$0.02FINAL $6.89 -$5.08 $3.99 $0.77 $2.91 $0.49 $2.91 $0.49 $1.55 $0.35114 INITIAL $4.90 -$3.42 $2.15 $0.38 $2.05 -$0.01 $0.48 $0.55 $0.72 $1.47114 11.09 MWD $7.60 -$5.06 $2.61 $0.45 $2.26 $0.48 $2.69 $0.29 $2.74 $0.03114 35.62 MWD $7.17 -$5.12 $2.96 $0.45 $2.24 $0.76 $2.70 $0.09 $2.50 $0.12114 59.88 MWD $6.91 -$5.18 $2.97 $0.32 $2.65 $0.14 $2.65 $0.14 $2.87 $0.09Since the excess reactivity (a measured value) appears to be extremely sensitive to uranium 235content, surveillance data for excess reactivity was used to standardize integral control rodworth values. Integral control rod worth values were linearly interpolated to measured excessreactivity values between the calculated excess and integral control rods worth for 100% and90% uranium 235 values. Shutdown margin was calculated as excess reactivity less the sum ofthe integral worth of all control rods. The worth of the most reactive control rod was added tothe shutdown margin to simulate the most reactive rod fully withdrawn.

,RQUEST FOR ADDITONAL INFORMATION (July 31, 2015: Adams Accession ML15211A362):The guidance in NUREG-1537 Section 4.5.2, "Reactor Core Physics Parameters," requests theapplicant provide calculations of certain core physics parameters and compare them withapplicable measurements. A request for additional information was sent to you in a letter datedJuly 25, 2012 (ADAMS Accession No. ML15211A638), regarding a comparison of calculated andmeasured values for reactivity parameters. In response, your submittal dated July 15, 2015(ADAMS Accession No. ML121500308), provided calculations for control rod worth and excessreactivity in Tables 11 and 12, respectively. In addition, you provided comparisons of these\calculations and measurements in Table 13.a. In Table 11, your calculations of excess reactivity and control rod worth on 3/16/1992,7/24/2007, 6/4/2008, 6/11/2008, 6/14/2010, 6/23/2010, 7/25/2011, and 8/2/2012indicate the shutdown margin for the UT TRIGA reactor would have a positive reactivity.A positive reactivity would not meet Technical Specification 3.2 Shutdown Margin.Please justify why the positive reactivity presents an acceptable level of safety regardingshutdown margin for the UT reactor.b. In Tables 11 and 12, your calculations for excess reactivity are consistently higher thanthe measured values. The average bias between the calculated and measured values is$2.56. Please justify why this bias presents an acceptable level of agreement inmodeling the UT reactor.c. Control rod worth in Table 13 appears to be inconsistent. For example, the RegulatingRod worth varies between +12.7% and -21.2%. Specifically, the values for two separateRegulating Rod data points on 7/13/2012, changes from -10.1% to +12.7%. Similarly,Shim-2 worth varies between -23.1% and +30.1%. Please justify why these variationsindicate an acceptable level of agreement in modeling the UT reactor.

RESPONSE:Calibrating model data to measured excess reactivity values indicates all shutdown margins arenegative by a wide margin.Most of the fuel elements in the current UT reactor had been partially burned at other facilitiesprior to use in this facility. The amount of uranium remaining in an element is based on totalcore burnup, distributed across the elements in the core. TRIGA peaking factors suggest asmuch as 40% variation in neutron flux across the core from center to peripheral elementsduring operation. Investigation of the sensitivity of reactivity calculations to uranium 235 massindicates a 10% change in mass has on the order of $5 impact on excess reactivity. Given thepotential uncertainty in the fuel mass in the partially burned elements, excess reactivity maynot be a reasonable function for validating modeling. When reactivity values are adjustedbased on calculations assume nominal uranium 235 and 90% of the nominal values and onmeasured excess reactivity values, integral control rods worth agree to measured values to areasonable level.Previous comparisons were based solely on burnup, and did not consider core configurations.As a result, reactivity values for different core configurations were identified by burnup datethat did not correspond in all cases to actual core configuration. In considering coreconfigurations and core burnup values, there is general agreement.ANALYSISAs previously noted, the initial UT TRIGA core was principally composed of previously (lightlyburned) irradiated fuel elements. Special Nuclear Material records were used as a basis for theoriginal uranium composition in material calculations. However, burnup calculations for theNuclear Materials Management and Safeguards System (NMMSS) track total facility inventorybased on core burnup, as opposed to individual fuel elements. Burnup in individual fuelelements varies significantly from average core burnup, and exacerbate potential uncertainty incalculating element-specific burnup.To evaluate the sensitivity of reactivity values to varying uranium 235 content, calculationswere performed for fuel composition using the nominal values from special nuclear materialrecords, and then with 90% of the uranium 235 specifications in the records. The differencebetween the two values was calculated. The current core contains 114 elements, andadditional data is provided for two burn intervals. The deviations in excess reactivity are muchhigher than the difference in individual integral control rods worth.

Table 1, Reactivity Based on 100% SNM 235 Material Composition and Difference at 90%EXCESS A RR A SH1 A SH2 A TR A90 INITIAL $7.21 -$5.31 $4.33 $0.26 $4.30 -$0.16 $2.53 -$0.03 $2.37 $0.08FINAL $7.09 -$5.18 $4.41 $0.42 $3.94 $0.22 $2.38 $0.39 $2.38 $0.1389 INITIAL $8.12 -$5.05 $3.63 $0.56 $3.07 $0.74 $2.66 $0.35 $2.66 $0.26FINAL $7.69 -$5.07 $3.65 $0.46 $2.96 $0.24 $2.61 $0.61 $2.53 $0.1492 INITIAL $6.91 -$5.16 $4.55 $0.50 $3.08 $0.90 $2.35 $0.45 $2.05 $0.28FINAL $6.54 -$4.79 $4.07 $0.72 $3.09 $0.69 $2.06 $0.84 $1.67 $0.3995 INITIAL $7.44 -$5.09 $4.27 $0.04 $3.16 $0.40 $2.69 $0.21 $2.41 $0.05FINAL $6.32 -$5.36 $4.04 $0.28 $3.18 $0.00 $2.72 -$0.07 $2.24 $0.04103 INITIAL $8.04 -$5.69 $3.89 $0.31 $2.98 $0.35 $2.48 $0.20 $2.04 $0.27FINAL $7.71 -$3.59 $4.05 $0.07 $3.24 $0.12 $2.59 -$0.12 $2.19 -$0.06102 INITIAL $5.67 -$5.58 $2.69 $0.07 '$3.28 $0.33 $2.26 $0.26 $3.21 -$0.08FINAL $4.87 -$4.97 $2.54 $0.81 $2.13 $0.46 $2.43 $0.39 $2.96 $0.30104 iNITIAL $7.41 -$6.71 $3.32 $0.20 $2.77 $0.30 $3.00 $0.00 $2.76 $0.05FINAL $6.19 -$5.49 $3.43 $0.09 $3.20 -$0.13 $2.86 $0.14 $2.71 $0.06108 INITIAL $6.85 -$6.29 $3.09 -$0.22 $2.36 -$0.13 $2.68 -$0.16 $2.63 -$0.25FINAL $6.96 -$5.45 $3.07 $0.42 $2.51 $0.34 $2.98 $0.16 $2.66 -$0.03110 INITIAL $6.96 -$5.45 $4.24 -$0.09 $2.35 -$0.07 $3.03 -$0.07 $1.66 -$0.02FINAL $6.89 -$5.08 $3.99 $0.77 $2.91 $0.49 $2.91 $0.49 $1.55 $0.35114 INITIAL $4.90 -$3.42 $2.15 $0.38 $2.05 -$0.01 $0.48 $0.55 $0.72 $1.47114 11.09 MWD $7.60 -$5.06 $2.61 $0.45 $2.26 $0.48 $2.69 $0.29 $2.74 $0.03114 35.62 MWD $7.17 -$5.12 $2.96 $0.45 $2.24 $0.76 $2.70 $0.09 $2.50 $0.12114 59.88 MWD $6.91 -$5.18 $2.97 $0.32 $2.65 $0.14 $2.65 $0.14 $2.87 $0.09Since the excess reactivity (a measured value) appears to be extremely sensitive to uranium 235content, surveillance data for excess reactivity was used to standardize integral control rodworth values. Integral control rod worth values were linearly interpolated to measured excessreactivity values between the calculated excess and integral control rods worth for 100% and90% uranium 235 values. Shutdown margin was calculated as excess reactivity less the sum ofthe integral worth of all control rods. The worth of the most reactive control rod was added tothe shutdown margin to simulate the most reactive rod fully withdrawn.

Table 2: Calculated Reactivity Values, Adjusted for Measured Excess ReactivityARO RR SHi SH2 TR SDM TS SDM90 0.00 $5.53 $4.41 $4.25 $2.52 $2.43 -$8.08 -$3.6731.07 $5.53 $4.54 $4.01 $2.49 $2.47 -$7.98 -$3.4489 31.07 $5.53 $3.92 $3.45 $2.83 $3.02 -$7.69 -$3.7731.31 $5.50 $3.85 $3.06 $2.87 $2.69 -$6.97 -$3.1292 31.31 $5.50 $4.69 $3.33 $2.47 $2.21 -$7.20 -$2.5136.82 $4.59 $4.37 $3.38 $2.40 $5.16 -$10.72 -$5.5695 36.82 $4.59 $4.29 $3.39 $2.81 $2.48 -$8.37 -$4.0845.72 $5.69 $4.07 $3.18 $2.71 $2.24 -$6.52 -$2.44103 45.72 $5.69 $4.02 $3.13 $2.56 $2.27 -$6.29 -$2.2781.19 $5.77 $4.09 $3.30 $2.53 $2.12 -$6.27 -$2.18102 81.19 $5.77 $2.69 $3.27 $2.25 $3.22 -$5.66 -$2.39106.12 $5.55 $2.43 $2.07 $2.37 $2.84 -$4.16 -$1.32104 106.12 $5.55 $3.38 $2.85 $3.00 $2.79 -$6.47 -$3.10121.82 $5.04 $3.45 $3.17 $2.89 $2.75 -$7.22 -$3.77108 121.82 $5.04 $3.13 $2.64 $2.79 $2.69 -$6.22 -$3.09186.53 $4.45 $3.27 $2.66 $3.06 $2.63 -$7.17 -$3.90110 186.53 $4.45 $4.20 $2.32 $3.00 $1.64 -$6.71 -$2.51204.88 $5.79 $4.16 $3.01 $3.01 $1.66 -$6.06 -$1.90114-2 226.17 $5.56 $2.79 $2.45 $2.80 $2.78 -$5.27 -$2.47114-3 226.17 $5.56 $3.10 $2.48 $2.73 $2.60 -$5.35 -$2.25Control rod worth data is measured periodically to verify that the minimum shutdown marginrequirements of Technical Specifications are met. For various reasons the burnup atsurveillances does not always correspond well to burnup assumed in analysis (used todetermine material compositions for the initiation and termination of core configurations).Previous work did not recognize this difference, comparing only calculated and measuredreactivity at the closest applicable burnup.This effort included ensuring that calculated and measured reactivity values are compared forsimilar burnup values and core configurations. Data indicates the model is consistent withoperating data. The comparison of calculated to measured reactivity data (Table 3) iscalculated as:SkM -SkcSkMWhere D is the deviation from measured values, t6kuis the reactivity from measured data, ande6kMis the reactivity based on the model.

Table 3, Comparison Measured and Calculated DataSURVIELANCE CALCULATION REACTIVITY COMPARISONRODDATE MWD CORE MWD RR SH1 SH2 TR SUM SDM07/01/92 0.00 90i 0.00 -8.02% -40.28% 20.55% 25.37% -0.51% -1.37%04/27/00 31.31 90f 31.07 -0.80% -15.16% 8.67% -4.72% -3.35% -5.81%89i 31.07 12.95% 0.96% -3.81% -27.96% -1.13% -1.96%89f 31.31 14.34% 12.11% -5.08% -13.91% 4.59% 7.56%92i 31.31 -4.17% 4.31% 9.62% 6.31% 2.86% 4.56%07/30/01 45.81 95f 45.72 2.76% 1.77% 7.92% 6.87% 4.47% 20.41%103i 45.72 4.13% 3.51% 12.86% 5.68% 6.28% 23.22%11/14/02 81.29 103f 81.19 4.93% 1.12% 8.16% 15.56% 6.70% 13.10%104i 106.12 -1.36% -2.59% 7.63% 16.08% 5.25% 6.45%07/18/05 121.93 104f 121.82 -12.45% -7.85% 7.95% 16.27% 1.37% -4.93%07/25/07 186.65 108f 186.53 -14.98% 3.19% 7.36% 20.68% 4.85% 0.02%06/29/10 226.30 1141 215.97 3.77% 3.36% 9.95% 11.46% 7.39% 10.75%__________1142 240.50 -6.82% 2.24% 12.31% 17.26% 6.71% 9.39%

LOSS OF COOLANT ACCIDENT ANALYSIS FOR THE UNIVERSITY OF TEXAS AT AUSTIN TRIGA REACTOR1. IntroductionThe loss of coolant accident (LOCA) analysis assumes steady state reactor operation at equilibrium(limiting core configuration conditions) followed by a reactor scram with the water coolingsimultaneously replaced with air cooling. The analysis models radial heat transfer from the center ofthe element outward to the air at the axial location/segment of the hot channel fuel element withthe maximum specific power.This LOCA analysis includes (1) an overview of the analysis, (2) specific characteristics of UT TRIGAsystem, (3) the basis of thermodynamic analysis, (4) development of the UT finite element analysismodel, (5) validation of the model against independent analytical method and against measureddata, and (6) analysis of the thermodynamic characteristics following a LOCA with initial conditionsestablished by the limiting core configuration.2. UT TRIGA CharacteristicsHeat generation following shutdown is a product of decay heat from fission products generatedduring operation, and has the same spatial distribution as power generation during operation.Analysis requires calculation of decay heat as a function of time. Thermodynamic properties ofTRIGA fuel are taken from reference material. A set of derived thermodynamic properties (i.e.,dimensionless numbers) is calculated. The decay heat, fuel geometry, and derived thermodynamicproperties are incorporated in model to simulate time dependent thermal dynamic response loss ofwater coolant.2.1. Decay HeatCalculations with TRACE indicate the maximum power for a fuel element with an acceptablecritical heat flux ratio of 2.0 is slightly less than 24 kW; the assumed initial condition for themaximum power in a fuel element is therefore 23 kW. Neutronic analysis with the fuel elementdivided into 15 equal axial segments shows the maximum power generation in a single axialsegment is 1.2 times the average segment or 1.84 kW for the initial conditions of the powergeneration in the maximum segment of the "hot channel." The decay heat is simulated as aheat source within the fuel element geometry.The decay power fraction remaining after an abrupt shutdown is found by equation1:0.04856 + 0.1189 .loglo t -0.103 * (log10 t)2 +F 0.000228 * (log10 t)3R(t) --(11 + 2.5481

  • log10 t -0.19632 * (log10 t)2 + 0.05417. (log10 t)3()The fuel temperature of the element producing the maximum power level in the core (hotchannel) is the most severe condition for heat transfer from the core during operation. For thelimiting case, the maximum specific power and the decay power fraction in the fuel element iscalculated from the maximum axial peaking factor for the fuel element:qgoen,j(t, r) = 1.2. qgen(r)" R(t) (21 Kansas State, "Kansas State University Safety and Analysis Report '06." KSU, Manhatten, 2006.

The radial distribution of power in each element remains constant, while the magnitudedecreases with time after shutdown according to eqn. (1).2.2. Fuel Element GeometryThe fuel element model in this analysis is a set of concentric cylinders representing a zirconiumrod at the center, the fuel matrix, a gas-gap between the fuel and cladding, and cladding. Thedimensions are taken from the GA drawings and UT Technical Specifications. The Zirconium fillrod diameter is 0.25 in (0.6125 cm) in diameter. The fuel matrix outer diameter is 1.47 in(3.6015 cm) diameter. The gas gap is approximately 0.005 in (1.97E-3 cm). Cladding is 0.020 in(0.0489 cm) thick. The total heated length of the fuel (section with Zr-U fuel matrix) is 15 in,segmented for thermal hydraulic analysis into 15 equal lengths. In this analysis only the verticalsegment with the highest heat generation rate is considered.2.3. Fuel Element Thermodynamic propertiesSimnad2 provides a number of mechanical characteristics and equations for fuel quantities. Thethermal conductivity (k) is given, density is calculated from a given equation for a specific Zr:Hratio of 1.6. Density is based off of an equation for the 8.5 wt% U:1P U wt% + (1-pzUWt%/)(3Where Ut,% is uranium weight per cent, pu is the density of uranium, and pu is the density ofzirconium. Simnad provides the temperature (T) dependent volumetric heat capacity (cp, vot):cp,1,o1 {3 = 2.04 + 4.17e -3

  • T (4a)Specific heat capacity (cp,ffiei) is calculated as the ratio of eqn. (4) to eqn. (3).c' et --- PF4bCp~fuelkg
  • KJ Cp,vo1(4b3. Basis of Thermodynamic AnalysisThe general thermodynamic basis in this analysis is based on an energy balance:Est =Egen + Ein, -- Eout (5a)2 M. T. Simnad, "The U-ZrHx Alloy: Its Properties and Use in TRIGA Fuel," Nuci. Eng. Des., vol. 64, pp. 403-422,1981.

Where, Est is the stored energy in the structure, Ege is energy generated within the structure,Ein1 is energy transferred into the structure, and [tout is the energy transferred out of thestructure. This model translates into:dTp.*V. Cp

  • j=qefl+ qcofld+ qconv (5b)Stored energy (and the associated temperature change) is a function of material density (p),specific heat (Cp), volume (1/), and the conduction, convection, and generation terms (qen qcond,and qcn respectively).3.1 Stored Energy (op V .c*T)Energy storage is related to material properties and temperature, an important factor incalculating the temperature transient analysis.3.2 Energy Generation (qgen)Energy generation in the core is a result of fission inside the element.3.4 Conduction Heat Transfer (qcond)Heat transfer through conduction within the radius of the fuel element and cladding ismodelled with Fourier's law of conduction using radial geometry:dTqc~ k s dr (6)Where k is thermal conductivity, A, is the surface area through which heat transfer occurs,and rris the rate of temperature change with respect to radial displacement. Asrecommended by Fenech3, the gas gap is approximated as thermal conductivity,calculated by the gas gap coefficient divided by the radial thickness of the gap.3.5 Convection Heat Transfer (qcon v)Convection applies to the surface element where heat is transferred from the fuel elementto the surrounding air. Convection heat transfer is modelled using Newton's law of cooling:Where the wall surface area is A5, the wall temperature is T5, and the bulk coolanttemperature is Ti[ and the heat transfer coefficient is h. The convection heat transfercoefficient is calculated from dimensionless numbers. For natural convection, thesignificant dimensionless numbers are the Prandtl (Pr), Grashof (Gr), and modifiedH. Fenech, "Heat Transfer and Fluid Flow in Nuclear Systems," Pergamon Press (1981)

Rayleigh (Ra) numbers. These values, combined with Nusselt (Nu) correlations, lead to theheat transfer coefficients which are used to find the energy transferred via convection.a. The Prandtl number is a measure of the fluid's kinematic diffusivity (v) to thermaldiffusivity (a) of the fluid:Pr =--,IaWhere thermal diffusivity of the fluid is calculated as:kcp

  • p(9S6)b. The Grashof number (and implicitly, values derived from the Grashof number) aredependent on both surface temperature and channel temperature. The Grashofnumber (for natural circulation heat transfer) is defined as:g./ (T- Tier) *(10)Where g is the acceleration due to gravity, /3 is the thermal expansioncoefficient of the coolant, s is the channel width, and/p is the dynamic viscosityof the coolant.c. The Rayleigh number is calculated [Kaminski4],:Ra =gIsr= fp2(Twaii -Tflu~d)s (1Ra Pr (1Where s is the total channel width.d. The Nusselt number relates the conductive and convective heat transfer effects ofthe fluid. The heat transfer is driven by the temperature difference and is foundthrough the Nusselt number. A correlation for natural circulation in verticalchannels for the Nusselt number4is:-1/2Nu3-.(12)e. The average heat transfer coefficient, h,, can be found from the average Nusseltnumber using:4D. A. Kaminski,, M. K. Jensen, "Introduction to Thermal and Fluids Engineering," John Wiley & Sons (2005)

--- 'k (13)4 The UT LOCA ModelThe University of Texas Loss of Coolant model is a combination of finite element analysis (FEA) forsteady state, and transient fuel conditions, as well as an air channel analysis sub-section to provideeffective estimation of air channel heating. This channel sub-model leads to proper parametricvariation analysis by giving a real world upper bounding temperature.4.1 Coolant Air TemperatureIn order to find the limiting values of the channel air temperature, a separate, one dimensionalvertical model was created independent of the FEA model and geometry used in finding the fueltemperatures. It utilizes an elemental, vertical, constant temperature surface interfacing withbuoyant air.The temperature entering the bottom of the channel surrounding the fuel element is thelimiting room air temperature following a loss of coolant, 20°C. The rise in temperature is foundby segmenting the pin vertically. Each iteration has a specific heat flux relative to thetemperature difference between the surface and the air, its specific dimensionless parameters,and a constant surface temperature that is user defined. The limiting conditions set the surfacetemperature at 950G. This was the value used later as a limiting factor in FEA parametricvariation.The change in air temperature across each segment is a function of the heat generated in thesegment and the heat transfer coefficient calculated from local non-dimensional parameters.Heat transfer characteristics in convection depend on intrinsic and extrinsic material propertiesand fluid temperature, with the heat transfer coefficient calculable though the use ofdimensionless numbers. The temperature rise of the fluid entering the region of heat transferfor subsequent segments is the exit temperature for the preceding segment, i.e., the rise intemperature across the previous segment added to the temperature of the coolant entering theprevious segment. The channel flow heat up model provided an order of magnitude estimationleading to proper parametric variation.4.1.1 Fluid flow, and thus the characteristic velocity, is driven by natural convection and isdominated by the buoyancy driven numbers in the Rayleigh number (Rag), the productof the Grashof and Prandtl numbers.a. The change in coolant temperature from fluid flow across a segment begins byfinding the appropriate Rayleigh number (eqn. 11) for the ith~ segment [Kaminski6],sT. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. DeWitt op cit & C. 0. Popiel and J. Wojtkowiak, "Simpleformulas for thermophysical properties of liquid water for heat transfer calculations (from 0 to 150 degrees C) (vol19, pg 87, 1998)," Heat Transf. Eng., vol. 19, no. 3, pp.87-101, 1998.SD. A. Kaminski,, M. K. Jensen, "Introduction to Thermal and Fluids Engineering," John Wiley & Sons (2005) then the segment's Nusselt number (eqn. 12), then the heat transfer coefficient (eqn.13). With the heat transfer coefficient now found, the heat flux (q') is found using:q"= i (Ts -Tm~f, ) (15)Where TIs is the cladding surface temperature and Tinli is the heat sink temperature.b. Heat flux is used to find the Modified Rayleigh number. For uniform wall heating, themodified Rayleigh (Ra*)3 is:R =g "f3."q'" p2. c. s4Ra k2 (16)Where g is the acceleration due to gravity, fi is the thermal expansion coefficient,and/p is dynamic viscosity.c. The modified Rayleigh number leads to the characteristic channel velocity7 (Uz):Uz = (17)d. The change in temperature for fluid flow across a segment of the fuel element alongthe (axial) direction of flow can be calculated with:Q, = rh .cp. AT (18)e. Where the A~is calculated as:".AFE,AT = AlwU cP(19)Where Uz is calculated from eqn. (19). This AT is added to the segment's inlettemperature and becomes the inlet temperature for the next segment. The lastsegment's channel temperature represents the culmination of all the heating.4.1.2 As an independent calculation to determine limiting values of air temperature, thetemperature rise was found through standard gas laws.a. The characteristic velocity gives a stay time for the air (heated length divided bycharacteristic channel velocity). This allows the change in energy to be calculated asfollows:dE = q" " A5s* tstay (20)7K. Vafai, C. P. Desai, S. V. Iyer, and M. P. Dyko, "Buoyancy Induced Convection in a Narrow Open-Ended Annulus,"J. Heat Transfer, vol. 119, p. 483, 1997.

Where tstay is the time the cooling air is in contact with the cladding surface.b. By using the density of air and the volume of the channel, the mass of the air in thespace at any given time can be found, by neglecting density changes. Using theequation below, the change in temperature can be found:q =mrcyAT (21)AT = q- =Tilit -+/-AT4.1.3 Results of calculations for limiting values of the channel air temperatureThese two methods routinely agreed across variations in surface temperatures, with thelimiting channel temperature, of 20°C inlet and 950°C surface temperature, being35.3°C.4.2 Finite Element Model Geometry and BasisThe calculation of temperature distribution in the fuel element is accomplished by using theprinciples of finite element analysis. The fuel element geometry is based on a cylindricalsegment. The axial height of the segment is the total heated length (0.381 m) divided by thenumber of segments (15). Radial dimensions are taken from General Atomics drawings.

dyINiFigure 1. Finite Element Radial GeometryThe Finite Element Model radii used in computation was selected based on both parametervalidation and computational power available. The limiting geometric figure of concern is theBiot number, which relates convective and conductive aspects of the element to its volume tosurface area ratio. It is determined using the equation below:Bi = h'Lc(2)Where, the characteristic length, 4c, is defined as the volume to surface area ratio:VL- (22b)Differential radii in the outer portions of themodel were chosen to most accurately subdividethe real geometry of the cladding and the gas gap. Internal fuel differential radii were chosen tominimize the Biot number. In addition to the Biot number, the Fourier number is a transientfigure of merit related to constants that determine time response and the geometry:a'tFo = -(23According to Bergman4, the Blot number must remain below 0.1, and the Fourier number mustremain below 0.5 for lumped parameter analysis to be valid. This was the merit to which thedifferential radii are chosen.

4.3 Steady State Finite Element AnalysisTo create a valid transient condition, a valid steady state initial condition must be found. Tofacilitate this, each element is assessed using an energy balance equation across the element.Since the steady state model is not time dependent, the energy balance is reduced to:Ef+/- +Egen = o0ut; oQUt=0 E=i+EeIn this analysis, energy flow is considered into the element. Fig. 2 illustrates element energybalance and temperature relationship.(24)qss 8enFigure 2. Finite Element Energy BalanceA matrix form of this energy balance is developed to solve for the temperature profile.=bWhere, is a vector representing the radial temperature profile, and b is a vector representingthe energy generation and non-temperature dependent terms. Below is the development ofthe steady state finite element equations. The cladding end element is the only elementcontaining a convection term, while fuel elements are the only ones containing generationterms. The following relationships are incorporated in the elements of the matrix equations:(25)

Conduction and Convection Terms:qgen,ss, r =qmax "q(r) " t. dy "(ri -ri_2qconv, ss = hwater" 11" rmax " dy" (Ts -Tinyf)2. l .dy .kfuel (Ti+_1-Ti)(26a)(26b)(26c)qcond,SS =In (2rargeGeneration and Temp~erature Independent Termsb1 = -qmax "ir. dy .(z -ri2)bend~l,4 0 ; (No heat generation in cladding/gas)bend = hwater *it rendt dy "T1ny[ a11a2Matrix Elements for A =LaendJ(26d)(26e)(26f)2.n dy .kfuel (T_1 -Ti)a1= --(26g)-2" i. dy. kfuel (T1+1 -Ti)in (26h)ai .=r2

  • t. dy" kfuet(gas, clad) *(T1.+1-Ti)I... -I in (r2~lI \r2~3)2 it dy kfuet(gasclad) (T1~1 -T1)+'~ ~r~i~i~)_(2. " n dy. kfuel(gas,cladl). (Til -Ti)In (r2i-1)2" r. dy. kfuel(gas,cladt) *(T/1~ -- Ti)' in (r21+1)

(26i)[2. *. dy"klad" (Tend-1l Tend)aend = .... in ( 'nd"ren )+nwadr-r*/ Y)krend-lMatrix Formula= /11. The energy generation term in the element is a function of both its axial and radial position. Thehighest axial peaking factor (1.2) was used to represent the cylindrical segment generating themost power. The radial peaking factor, q(r), is found through a curve fit to neutronic codeoutput, with the highest axial peaking factor of 1.2.MATLAB was utilized to build and solve the equation using native commands that maximize theefficiency and accuracy of the matrix inversion method.4.4 Transient Finite Element AnalysisHeat transfer analysis for systems that have time variation can be analyzed using lumpedparameter analysis where the internal resistance to heat transfer is small compared toconvection resistance, and the characteristic time constants are similarly related. The transientportion of the model takes the initial steady state temperature profile and systematically walksit forward with time. The basic concept of an energy balance as used in the steady state analysisis maintained, with the time dependent components now considered in addition to the otherterms. In the UT LOCA model the loss of coolant accident is considered to be instantaneous, andthus the cooling properties switch from water to air.Est=Etn Eout O 0 -(27a)dT= qcond + qlconv + qe (27b)+1- Ti)_PVZCP A -- qcond + qconv + qgn(27C)This leads to the transient analysis equation set which is related to the steady state equations asfollows:AtTPI= 1] + TP (28)The differential time element is selected based on the merit of the Fourier number previouslymentioned. Additionally, the code calculates a number of output values including a two-dimensional matrix B whose horizontal dimension represents the radial temperature distribution and whose vertical axis represents time. This allows three essential modelparameters to be extracted. First, the cladding surface temperature versus time is extracted andused to find peak cladding temperature. Second, the temperature profile across the pin at t, canbe found. Third, the maximum temperature both radially and through time can be found.5 Model Validation5.1 Comparison of TRACE and the UT MATLAB model Steady State Temperature ProfileThe core configuration contains 114 fuel elements, with a core radial peaking factorderived from SCALE physics calculation for the core (prior to January 2016) of 1.6, and amaximum axial peaking factor of 1.2. The current normal operating power is 950 kW.The power generated in the maximum segment of the hot channel for comparison usingdata prior to January 2016 is therefore 12.5 kW.The steady state solution using water coolant was developed for the maximum powerlevel in a fuel element operating at 12.5 kW and compared to the TRACE calculations(Fig. 3). The TRACE and FEA calculations are in substantial agreement with experimentaldata.4504003500J2..ES200150FUEL ELEMENT TEMPERATURE PROFILE(INITIAL CONDITION)i ~',I..-.-MatLAB ICTemp (C) e"1 0 0 .................................0.002 0.007 0.012 0.017 0.022Radial Displacement (in)Figure 3. TRACE and UT LOCA model steady state temperature profiles5.2 Comparison of FT2 Observations and Calculations (TRACE, UT MATLAB Model) SteadyState Temperature Response to Power OperationThe MATLAB finite element analysis was applied at power generation in an elementfrom 200 W to the 12.5 kW, and the maximum element temperature compared to theTRACE and FT2 measurements (taken prior to January 2016) across the range.

TRACE and the MATLAB based steady state temperature calculations in radial locationsassociated with thermocouples are essential the same. There is good agreementbetween calculated and observed values with some deviation at higher power levelswhere the heat transfer is presumably affected by the development of bubbles thatenhance heat transfer and reduce fuel temperature (Fig. 4).COMPARISON TRACE, MATLAB AND FT2 DATA450400350t3 300* 250200E 150wJI--S100IL5000 2000 4000 6000 800010000 12000 14000Fuel Element Power (W)--h--TRACE (Max) -a- MATLAB (Max) -.-FT2Figure 4. Comparison of Temperatures from Calculations and Observations at Varying Power Levels5.3 Comparison of FT2 Observations and Calculations (TRACE, UT MATLAB Model) TransientTemperature Response to Shutdown from Normal OperationsTransient fuel temperature was observed following a shutdown from power operationsat 950 kW (Fig. 5, FT2 Data). Calculations were performed to simulate the transientusing TRACE (Fig. 5, TRACE Calc) and MATLAB based model (Fig. 5, UT MATLAB). Thetemperature data is in good agreement.

L.EI-LL450400350300250200150100500FUEL TEMPERATURE FOLLOWING SHUTDOWNFROM STEADY STATE 950 KWI IJ ' I I_ _ _ _ i! i " " i !{! !ill ... ................. .................. '. ,.._ _.........._ iI110 100Time After Shutdown (s)1000-. UT MATLAB --TRACE Caic -aFT2 DataFigure 5. Fuel Temperature, Measuring Channel & Calculations Following Reactor Scram5.4 SummaryComparison of fuel temperature measuring channel data to calculated fueltemperatures during steady state and transient conditions is in good agreement. Theagreement between observations and calculations during steady state operationssuggests the method is fundamentally correct. The agreement between observationsand calculations during transient operations suggests the method will providereasonably accurate time-dependent calculations.6 ResultsThe UT MATLAB model calculation was performed for various values of both air channeltemperature and pin power. The radial temperature profile of the fuel element segmentgenerating the highest power in the core is provided in Fig. 6 following a shutdown from normalfull power operation at 950 kW operations with air cooling at inlet air temperature equal to UTrector bay nominal temperature.

LOCA CLADDING TEMPERATURE VS TIME FORINITIAL FUEL ELEMENT POWER OF 23 KW1000700500 [ .......................... ...!...........M axim um CladdingTernperature ......."o 40 !for 16°C Air: 935.40°C at 4701s400 After LOCA300 !Maximum CladdingTemperaturefor 20°C Air: 939.04°C at 4731s0 5000 10000 15000Time (s)Figure 6. LOCA Cladding Temperature vs TimeFuel element power level and inlet air temperature were varied to provide an indication tosensitivity to the parameters (Fig. 7). Line labels used in or significant to this analysis areprovided with label markers.PEAK LOCA FUEL TEMERATURE VS COOLINGTEMPERATURE1100lO5O : i-- -------------/1000 --.---........800 .750 i... .... .. ..~70 < --...... 0 ........65 ---. ... .600 .e. --15 25 35 45 55 65 75 85 95 105 115 125 135 145CoolingTemperature °(C)Element Power-&1.--1 17--1 19 2*--3--4Levels (RW) -2. ...14 ...17 -.- 8 ..... 9 ........2 --2 ...2Figure 7. Peak Fuel Temperature during Loss of Coolant AccidentFor reactor bay air at 16°C, the maximum fuel element power prior to LOCA initiation that couldachieve 950°C fuel temperature with air cooling is 23.6 kW. At 23 kW generated in the fuelelement during operation prior to the LOCA initiation (the maximum power generated in a fuel element in the limiting core configuration), air inlet temperature inlet less than 35°C iscalculated not to exceed 950°C fuel temperature. Therefore a LOCA following normal steadyoperation with a fuel element operating at 23 kW will not exceed the fuel temperature safetylimit.This analysis is extremely conservative in several important respects, including neglecting axialconduction, assuming an instantaneous loss of cooling water, assuming a complete loss ofwater, and assuming dry air.a. The UT LOCA model takes place at the point of highest axial power production and onlytransmits energy radially, while in reality the axial conduction effects would work to reducethe maximum fuel temperature prior to and during the transient.b. As shown in Fig. 5, water cooling immediately following shutdown reduces fuel temperaturesignificantly, with the measuring channel indicating l00°C decrease over about 13 seconds.A smaller quantity of stored heat reduces fuel temperature at the initiation of the LOCA.c. The most probable flow path for a LOCA is via failure in beam port casing. The beam portsare located slightly below core center, and a substantial fraction of the core structure andfuel elements will be in contact with pool water even if drained to the bottom of the beamports.d. The specific heat capacity of dry air is 1 Id/kg-K, but the reactor bay ventilation system isdesigned to control relative humidity for comfort. Specific heat capacity for moist air iscalculated8:cv= 1.005 + (8.

  • 0- T2 -+ 2.5
  • 107 T + 1.86)
  • HThe specific heat capacity of moist air increases with relative humidity, so that calculationswith dry air result in lower heat transfer and higher fuel temperatures. In addition, thenature of the event assures moist air in the cooling supply.8 http://www.engineeringtoolbox.com/