Predicting the deflection of electric
heater rods in a critical heat flux test
loop
Master thesis by: Patrik Andersson
TRITA-FYS 2014:34 ISSN 0280-316X ISRN KTH/FYS/--14:34—SE
Abstract
This paper concludes the simulations done on the bending of heater rods due to electromagnetic forces in the ODEN test loop at Westinghouse Electric fuel testing facility in Västerås, Sweden. The by calculating the Lorentz force on the rods and using Euler-Bernoulli FEM calculation the deflection of the rods were simulated. A model for the stiffness of the rods was created after bending tests on a rod.
The MATLAB simulations show that the current placement of support grids can be improved, resulting in a decrease in maximum bending by 40% and at the same time reducing the number of support grids. This will result in more representative results from the experiments as the rods are closer to their proper position and flow interference from the support grids is reduced. An added bonus is that less time needed to assemble/disassemble the assembly.
Table of Contents
1. Introduction ... 1 1.1. Literature review ... 1 2. Theory ... 2 2.1. Magnetic field ... 2 2.2. Force ... 2 2.3. Bending ... 2 2.4. Stiffness ... 2 2.5. Temperature ... 2 2.6. MATLAB ... 3 2.7. COMSOL ... 3 3. Model ... 3 3.1. Magnetic field ... 3 3.2. Forces ... 4 3.3. Bending ... 5 3.4. Stiffness ... 7 3.5. Temperature ... 8 4. Experiments ... 104.1. Determining stiffness by bending ... 10
4.1.1. Setup ... 10
4.1.2. Results ... 11
4.2. Benchmarking against COMSOL multiphysics ... 15
4.2.1. One material simulation ... 15
4.2.2. Two material simulations ... 16
5. Program ... 18
5.1. What the files do ... 18
5.2. Using the program ... 18
6. Results ... 20
6.1. Decreasing deflection ... 22
7. Conclusion ... 24
7.2. Future work ... 24
Bibliography ... 25
Appendix B: Bending data points ... 8 Appendix C: Manual ... 15 Appendix D: Code ... 32
Acknowledgements
I would like to thank my supervisors Rasmus Enlund and Professor Waclaw Gudowski for their help and support in answering my questions. I would also like to thank the BTA departments at
Westinghouse for providing me with the opportunity and resources to do this work. Lastly I would like to thank Mattias Stenlund and Bo Malmström for helping me use their tensile test machine.
Abbreviations
CHF Critical Heat Flux DAS Data Acquisition System
DNB Departure from Nucleate Boiling FEM Finite Element Method
MVG Mixing Vane Grid
IAPWS International Association for the Properties of Water and Steam IF97 Industrial Formulation 1997
1. Introduction
The main benefit of using nuclear power for energy production, high energy density, is also the reason for the large consequences of a large scale failure. In the best case only the plant is destroyed, at a huge cost for the utility. In worse cases release of radioactive material affect a large number of people both in the country where the plant is located and in neighboring countries, at huge
socioeconomic costs. These consequences are the reason that such strict safety standards imposed on the nuclear industry has on both from regulators and form the utilities themselves. Because of the strict safety standards in the nuclear industry it is important that new components go through rigorous testing before they are put into operation. The testing is done both through computer simulations using codes like: OpenFOAM and in experimental facilities designed to represent different parts of a plant.
One such facility is the ODEN Critical Heat Flux Test Facility in Västerås, Sweden. It is built to find the critical heat flux (CHF) for different pressurized water reactor (PWR) assembly designs. To simulate the nuclear reaction the rods in ODEN have built in resistive heaters that can supply a maximum of 12 MW of power. The rods also have built in thermocouples to register the temperature excursion that occurs at departure from nucleate boiling (DNB). (Waldemarsson, 2013)
A problem that arises from the use of electric heating is that the electric current through the rods will generate a magnetic field that will pull the rods together. This will reduce the mass flow between the rods and result in the test showing a lower CHF (Anglart, 2010).
The standard way to prevent this is to add extra support grids to the assembly. They will prevent a majority of the bending, but while they are smaller and weaker than the ordinary grids to minimize flow disturbance they do affect the flow and thus the result. Because of this support grids should not be used excessively or when the deflection of the rods is insignificant compared to other errors. To see the effect of the magnetic field on the rods a MATLAB program was written that uses FEM (finite element method) calculation to find the deflection of all rods in a test assembly from Lorentz forces. The benefit from having a fast way of determining the deflection is that the position of the support rods can be optimized to reduce deflection during assembly design process, thus leading to more accurate results from the tests without requiring a lot of extra time. Another benefit is that if support grids can be removed it will reduce the time needed to assemble the test assembly.
1.1. Literature review
The author could not find any reports on studies of heater element deflection during the literature review. This lack of reports can stem from two reasons. Either the problem is specific for nuclear CHF test as they require a combination of high current, high flow, small distances and high precision. While for example in a furnace the heater elements can be thicker, thus stiffer, and have more supports while their deflection matters little thus there is no reason to try to calculate the deflection if the heater elements. The other reason for the lack of reports is that the problem is not complex as all underlying phenomena have been thoroughly studied thus publishing a report on the subject might be seen as unnecessary.
rods. This program was used as an inspiration for the MATLAB simulations but it uses some strong simplifications that are not needed today because of increased computing power. Firstly it uses one of the analytical solutions of a clamped Euler-Bernoulli beam (Lundh, 2000) not a FEM solution, thus it only looks at one length between two supports and not the entire length so it cannot see how the other parts of the rod affect the investigated part. Another problem is that the program regards the rods as homogenous when their properties in fact vary along their length. The program also does not give any information on the distance between the rods in the assembly only the maximum deflection of the rod therefore it is hard to determine if the deflection is too large.
2. Theory
The theory needed to calculate the deflection of the heater rods in the test loop is not complicated but this section will show how it is connected and at the end the computer programs that were used to simulate and benchmark will be presented.
2.1. Magnetic field
The heat in the rods is generated by applying a strong current through them and the current will create a magnetic field around the rods. The magnetic field strength of this field, at any point in space, can easily be calculated by Biot-Savart law from the geometry and current (Encyclopædia Britannica Online, 2014). When this is done for all the rods the total magnetic field acting on one rod is found by adding the field strengths together.
2.2. Force
The Lorentz force is the force on a charged particle from the electric and magnetic fields that it inhabits. From the total magnetic field at the rod the Lorentz force can be calculated with the current passing through the rod. Thus the force and direction of the force in all points along the rods is found. A deeper explanation can be found in Chapter 3.2
2.3. Bending
When the magnetic force acts on the rods they will bend. How much they bend is calculated using the Euler-Bernoulli beam theory as the rod is a thin beam. Because of the varying geometry and temperature of the rod it is divided into 1mm sections and FEM is used to find the deflection in all points along the rod.
2.4. Stiffness
One important component of the bending calculation is the stiffness as it measures the rods resistance to bending. The stiffness is calculated using the second moment of inertia multiplied by the Young’s modulus. As the calculations are done on 1 mm sections the value for the geometry is constant across the section.
2.5. Temperature
To find the stiffness of the rods the Young’s modulus of their component materials is needed and to accurately get that their temperatures must be found. For this a simple thermo hydraulic model was used to find the water temperature and then another correlation is used with the water properties gained from the international association for the properties of water and steam, in their document on industrial formulations from 1997 (IAPWS-IF97) they recommend formulas for finding water
properties in industrial applications (IAPWS, 2007). With this the average temperature of the materials can be found using simple thermodynamics.
2.6. MATLAB
The program is written in MATLAB R2012b a high-level programing language optimized for mathematical operations created by The MathWorks Inc. It was chosen because it has a large amount of built in functions that relate to engineering and mathematical operations, has an
extensive help library both built in and online and because the author has extensive experience with it. (MathWorks Inc., 2014)
2.7. COMSOL
When the models in the program where tested to see how accurate their results where they were tested against models built into COMSOL Multiphysics version 4.4.0.150 made by COMSOL Inc. COMSOL is a general purpose program for simulation physics based phenomena using numerical methods. It has many different modules for simulating different phenomena but in this work the ones used where solid mechanics, beam and heat transfer in solids. The reasons for using COMSOL where that it is easy to use and the author had some prior experience with it. (COMSOL Inc., 2014)
3. Model
Ones the theory needed to find the deflection was defined then the unnecessary pieces had to be removed to get a simplified model that could be converted into computer code and solved for different scenarios.
3.1. Magnetic field
While the resistive heating elements are more current carrying cylinders than wires simulation in COMSOL have shown that the magnetic field produced is the same, outside of the cylinder. As the stronger field from the parts of the rod closer to the measuring point is canceled by the weaker field by the parts further away from the point.
The Biot-Savart law states that the magnetic field from a current carrying wire is
𝑑𝑩 =4𝜋𝜇0𝐼 𝑑𝒍 × 𝒓�𝑟2 ( 1)
Where 𝜇0 is 4𝜋 10−7 𝑁𝐴2 the permeability of free space, 𝑑𝒍 is a vector pointing along the wire, 𝑟 is the
distance between the wire and point and distance from the wire to the point where the magnetic field is measured and 𝒓� is the unit vector pointing from the wire to the measurement. Figure 1 below shows the geometry. (Encyclopædia Britannica Online, 2014)
Figure 1: Explenation of the geomery of Biot-Savarts law Figure 2: Explanation of the geometry for the magnetic field from a straight wire
When integrating equation (1) for a straight wire the solution becomes
𝑩 =4𝜋 𝑎𝜇0 𝐼 (cos 𝜃1+ cos 𝜃2) ( 2)
Where 𝜃 is the angle between the point and the ends of the wire and 𝑎 is the shortest distance from the point to the wire (Subashki, 2014). The direction of the magnetic field can easily be found using the right hand rule (National High Magnetic Field Laboratory , 2014). In Figure 2 the geometry of this can be seen.
To get the magnetic field in one of the elements in a rod the field contribution from all the other rods are added together.
3.2. Forces
The Lorenz force equation gives the electromagnetic force on a charged particle and can be seen below
𝑭 = 𝑞𝑬 + 𝑞𝒗 × 𝑩 ( 3)
Where 𝑞 and 𝒗 is the charge and velocity of the particle and 𝑬 is the electric field. When the charged particles are contained in a wire the force exerted on them is transferred into the wire and the equation changer to (Encyclopædia Britannica Online, 2014)
𝑑𝑭 = 𝐼𝑑𝒍 × 𝑩 ( 4)
As the rods are parallel to one another the force pulling one element in a rod towards another rod is described by the following equation
𝐹 =𝑙 𝜇4𝜋 𝑎0 𝐼𝐴 𝐼𝐵(cos 𝜃1+ cos 𝜃2) ( 5)
Because the rods end at the top adapter but continue past the bottom plate the force are stronger in the bottom resulting in a force distribution like the one in Figure 3.
𝑑𝑩
𝑟
𝒓�
𝑑𝒍
𝜃
1𝜃
2𝑎
Figure 3: Distribution of force on the rods
3.3. Bending
The bending calculations are based on the Euler-Bernoulli beam theory that assumes that when under deformation all sections of a non-deformed plane remain plane. For the FEM calculation the rod is divided into nodes where each node has two degrees of freedom, deflection in y (𝑤) and rotation (𝜃). Superscript 𝑒 indicates that a single element is being worked with. The deflection vector 𝒘 is a column vector that contains the displacements both in 𝑤 and 𝜃.
𝒘𝒆= � 𝑤1 𝜃1 𝑤2 𝜃2 � ( 6)
An element is the beam between two nodes and it is affected by a force in the y direction. The element is the basic building block in the FEM calculation and they are used to build the matrixes needed to solve the equation. In Figure 4 a representation of an element is shown, 𝑢 and 𝑤 are the deflection in x and y direction.
Figure 4: A rod element (Ferreira, 2009)
To get the values of the deflection and rotation between the nodes, interpolation with the Hermite shape function is done as they yield a smooth continues function.
To interpolate the deflection Hermite shape functions are used, these give the values of the deflection and rotation between the nodes (Miller, 2012).
𝒘 = 𝑵(𝜉)𝒘𝒆 ( 7) Top adapter Bottom plates 1 2 Force Full force
Where 𝜉 is the dimentionless distance from -1 to 1 and 𝑵 is the Hermite shape functions in a row vector. The superscript 𝑇 indicates a transposed vector.
𝑵𝑻(𝜉) =1 4 ⎣ ⎢ ⎢ ⎢ ⎡ 2 − 3𝜉 + 𝜉2 1 − 𝜉 − 𝜉2+ 𝜉3 2 + 3𝜉 − 𝜉2 1 − 𝜉 + 𝜉2+ 𝜉3⎦⎥ ⎥ ⎥ ⎤ ( 8)
The stain in the element and virtual external work done by the distributed force 𝑝 to the element are given by 𝑈 and 𝛿𝑊. 𝑵′′ is the second derivative of the Hermite function, 𝑑2𝑵
𝑑𝜉2. 𝑈𝑒 = 𝒘𝒆𝑇𝐸𝐼𝑧 2𝑎3 � 𝑵′′𝑇 1 −1 𝑵 ′′𝑑𝜉𝒘𝒆 ( 9) 𝛿𝑊𝑒= 𝛿𝒘𝒆𝑇 𝑎𝑝 � 𝑵1 𝑇 −1 𝑑𝜉 ( 10)
From these equations the stiffness matrix, 𝑲, and force vector, 𝒇, are defined and then simplified by 𝐿 =𝑎2. 𝑲𝒆= 𝐸𝐼𝑧 2𝑎3 � 𝑵′′𝑇 1 −1 𝑵 ′′𝑑𝜉 = 𝐸𝐼𝑧 2𝑎3 � 3 3𝑎 −3𝑎 3𝑎 3𝑎 4𝑎2 −3𝑎 2𝑎2 −3 −3𝑎 3 −3𝑎 3𝑎 2𝑎2 −3𝑎 4𝑎2 � =𝐸𝐼𝐿3𝑧 � 12 6𝐿 −12 6𝐿 6𝐿 4𝐿2 −6𝐿 2𝐿2 −12 −6𝐿 12 −6𝐿 6𝐿 2𝐿2 −6𝐿 4𝐿 � ( 11) 𝒇𝒆= 𝑎𝑝 � 𝑵1 𝑇 −1 𝑑𝜉 = 𝑎𝑝 3 � 3 𝑎 3 −𝑎 � = ⎣ ⎢ ⎢ ⎡ 𝑝𝐿𝑝𝐿/22/12 𝑝𝐿/2 −𝑝𝐿2/12⎦⎥ ⎥ ⎤ ( 12)
Because the strain comes from the work they are the same and this gives the final equation for the element
𝒇𝒆= 𝑲𝒆𝒘𝒆 ( 13)
To calculate the deflection for the entire beam the stiffness and force matrixes for every element are added together and the n solved giving the solution to the deflection and rotation of all the nodes in the beam. (Ferreira, 2009)
The beam setup used to represent a rod is that of a clamped beam supported by springs at the position of the grids. The bottom end of the beam is displaced towards the center because the O-rings centering them in the holes in the bottom plate are elastic. A schematic figure of the rod setup can be seen in Figure 5, the bottom displacement is not shown as it only appears when force is applied.
Figure 5: Beam setup
To add a spring to the beam model an element representing the spring is added it consists of the node the spring is attached to and a dummy node with only one degree of freedom. As the spring does not affect the rotation of the beam it only changes the 𝑤 part of the stiffness matrix, thus it can be simplified into the equation below.
𝑲𝒔= 𝐶
𝑠�−11𝑤𝑠𝑤𝑠 −1𝑤𝑠𝑤𝑑
𝑤𝑑𝑤𝑠 1𝑤𝑑𝑤𝑑 � ( 14)
𝑲𝒔 = 𝐶
𝑠�1𝑤𝑠𝑤𝑠� ( 15)
Where 𝐶𝑠 is the spring constant and the subscripts note which rod and column in the total stiffness
matrix the value gets added to, where 𝑤𝑠 is the beam node connected to the spring and 𝑤𝑑 is the
dummy node. As the dummy node will not affect the beam it can be disregarded and is thus simplified to the second form. (Ferreira, 2009)
To add prescribed displacements to the model change all the values in the stiffness matrix on the row corresponding to the node with the displacement to zero except for the one in the column representing the node, change that to one. Also in the force vector change the value for the node to the prescribed value.
3.4. Stiffness
The stiffness of the rod (𝐸𝐼𝑧), as can be seen in the preceding section, is needed to find the deflection
of the rod. It contains two components the Young’s modulus, a material property, and the second moment of area, a geometry property. 𝐼𝑧 for a pipe is given below (Lundh, 2000)
𝐼𝑧 =𝜋4 �𝑟𝑜4− 𝑟𝑖4� ( 16)
𝐸 for Alumina is given in the source (Sánchez-González, et al., 2007) at different temperatures and these data points are interpolated in Microsoft Excel to get the temperature dependent equations in Table 1. In the same table the 𝐸 for Inconel 625, copper C11000 and nickel 200 tubing is given but the values are taken from the COMSOL material database where they are given temperature dependent
Table 1: Young's modulus for the materials in the rod dependent on temperature in Kelvin
Material Young's modulus [Pa]
Alumina 3.5612 1011− 5.58 10−2 𝑇 + 2 10−4 𝑇2− 4 10−7 𝑇3+ 2 10−10 𝑇4 Inconel 625 2.157999 1011− 2.697218 107 𝑇 − 28023.67 𝑇2 Nickel 200 2.198604 1011− 4.976173 107 𝑇 − 6940.452 𝑇2 Top adapter Bottom plate
With this the stiffness of the individual components in the rod can be calculated, but because the beam model calculates the entire rod as one beam then a single stiffness must be calculated. According to Bo Alfredsson because of the axisymmetric geometry the total stiffness is the sum of the components in the cross section (Alferedsson, 2014). This was further verified by COMSOL simulations.
As the alumina part of the rod is divided into pieces, the loss in strength from this is modeled by having the FEM element at the position of the break have a stiffness as if it only contained Inconel. This will give larger stiffness loss than in reality as the element is longer than the real break but the difference is small.
3.5. Temperature
To get the temperature of the water at a given point the enthalpy and other properties of the inlet water is calculated with the IAPWSs data and the power of all the rods at that element is calculated. From this the energy added to the water is in the element is calculated and the new enthalpy is used to find water properties in the next element. When this is complete and water properties of for the entire simulated length of the rod exist rod wall temperatures are calculated with the equation seen below
𝑇𝑟𝑜𝑑 = 𝑇𝑤𝑎𝑡𝑒𝑟+𝜋 𝑟 𝑁𝑢 𝑘𝐷𝐻 𝑃 ( 17)
𝐷𝐻 =4𝐴𝜋 𝑟𝑅 𝑓𝑙𝑜𝑤 ( 18)
𝑁𝑢 = 0.023𝑅𝑒45 Pr0.4 ( 19)
Where 𝑃, 𝑟, 𝐷𝐻 and 𝐴𝑅 𝑓𝑙𝑜𝑤 are the power, outer radius, hydraulic diameter and water flow area for
the rod. 𝑁𝑢, 𝑘, 𝑅𝑒 and 𝑃𝑟 are the Nusselt number from the empirical Dittus-Boelter correlation for heating (Subramanian, u.d.), thermal conductivity, Reynolds number and Prandtl number for the water. (Wallenius, 2011)
Some assumptions are taken during the temperature calculations to save computational power as the change in material strength from any reasonable temperature interval is small. The first assumption is that the water has the same temperature across the entire cross section, i.e. perfect mixing. The second assumption is that all the rods hear the same amount of water, the flow area for one rod is the total flow area divided by the number of rods. The last assumption is that the rods only heat the water at the heated length, there is 0 resistance in the other parts of the rod and no heat diffusion in the rod.
The correlation is also of a simplified nature as it works for single phase flow and thus it will not take into account any phenomena that occur due to boiling, so the rod temperature will only spike when all the water is vapor and there will never be dryout or DNB. But even if this temperature spike would strongly affect the stiffness of the rod the test loop is built to find them and will lower the power when they are found so any extra bending at that point is irrelevant.
To find the average temperature in the rod metal and ceramic the equation, seen below, for the temperature difference between the outer and inner wall of cladding in a reactor is used.
∆𝑇𝑐𝑙𝑎𝑑=2𝜋𝜆(𝑇𝑃
𝑎𝑣𝑒) ln �
𝑟𝑜𝑢𝑡𝑒𝑟
𝑟𝑖𝑛𝑛𝑒𝑟� ( 20)
Where 𝑟𝑜𝑢𝑡𝑒𝑟, 𝑟𝑖𝑛𝑛𝑒𝑟 and 𝜆 are the outer and inner radius and the thermal conductivity at the
average temperature of the cladding (Wallenius, 2011). The thermal conductivity of Inconel 625 was taken from COMSOL an can be seen in Table 2
Table 2: Thermal conductivity of the cladding material dependent on temperature
Material Thermal conductivity �𝑊
𝑚𝐾�
Inconel 625 5.482405 + 0.01380594𝑇 + 1.678069 10−6𝑇2
The equation above need to be modified because it assumes that the heat is not generated in the cladding as it is in the heater rod. Because of this the temperature in the ceramic will be the same as at the inner wall of the cladding. Assuming that the heat transport through the cladding increases linearly from 0 at the inner wall, compared to the constant one used in the equation. The
temperature profile in the cladding with normalized parameters would be
𝑦 = −𝑥2
2 + 𝑥, ℎ𝑒𝑎𝑡𝑒𝑟 𝑟𝑜𝑑 ( 21)
𝑦 = 𝑥, 𝑛𝑢𝑐𝑙𝑒𝑎𝑟 𝑟𝑜𝑑 ( 22)
Where 𝑦 and 𝑥 are the normaliced ∆𝑇 and distance from the outer wall respectively. This results in an average metal temperature of 𝑇𝑟𝑜𝑑+ ∆𝑇𝑐𝑙𝑎𝑑∗ 0.33325 and a ceramic temperature of 𝑇𝑟𝑜𝑑+ ∆𝑇𝑐𝑙𝑎𝑑
2 . The simulations done in COMSOL gave the same temperature profile, which can be seen in
4. Experiments
Two kinds of experiments where done to see if the program would generate valid results. All the graphs generated during the benchmark simulations can be found in Appendix A
4.1. Determining stiffness by bending
The second type of experiment was to see if the calculated rod stiffness matched one from a real rod that had been used in the test loop.
4.1.1. Setup
To do this a tensile test machine, consisting of a Chatillon MT-DFS2-025-EU dynamometer mounted on a Chatillon LTCM-500 motorized test stand and Limess RTSS video extensometer with software to measure the change in distance between two lines was used. It was converted to serve as a bending machine for a double overhang beam with a point load in the center. As the dynamometer could measure both tension and compression, the conversion only required a new tip for the
dynamometer to be able to securely press on a rod and new holders to be able to place the samples in an optimal position for the camera to see.
As the rods weight is negligible compared to the force of the point load and that no load acts on the overhang it can be seen a simple supported beam with point load in the center. Thus the stiffness of the rod is proportional to the slope of the force-deflection plot. [American forest and paper
assosiation, Design aid No.6: Beam design formulas with shear and moment diagrams, 2007]
∆𝑦 =48𝐸𝐼 → �𝐹 = 𝐶 ∆𝑦,𝐹𝐿3 𝐶 =48𝐸𝐼𝐿3 � ( 23) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized distance from outer wall [~] Normalized temperature change [~]
Figure 6: Comparison between proposed equation and COMSOL simulation of the temperature in the clad, straight line is heat generation behind cladding and curved line is heat generation in the cladding
The bending was done on rod WH-600-H3s heated length after it had been cut into pieces 600-710 mm long for ease of handling and transportation. The pieces where named A-G and can be seen in Figure 7below.
Figure 7: The sections of WH-600-H3
In addition to the rod, grids were tested to find their spring constant. The tested grids where two mixing vane grids, VSH 10029E35G02 #03956 that was new and an older of the same type, and two support grids, GRDT-04710 and GRDT-04719 both new. As seen below the spring constant can be found in the same way as the stiffness.
𝑘 =∆𝑦 →𝐹 {𝐹 = 𝐶 ∆𝑦 , 𝐶 = 𝑘} ( 24)
Because the equipment had not been used in this way before, two copper and two steel pipes with known diameter also had their stiffness measured to serve as reference to calibrate the results after. The last type of test conducted was done using rod section A and getting steady state reading over longer time periods than before to find the noise in the dynamometer and camera readings. For the tests all the rods had a length between the supports of 500 mm and every test run had a target load of 0-75-0 N, starting at 0 N going up to 75 N and then going down to 0 N again. The 75 N limit was imposed to have margin to the 100 N max load of the dynamometer. Because of the loading and then unloading errors from slipping or settling of the equipment could be found because then the rod would not return to its starting point. Another advantage of this load program is that every test supplies two slopes.
The test matrix for the rod sections can be seen below in Table 3, every section was tested twice at multiple points along its length to see the effect of the small alumina pieces.
Table 3: Test matrix for the rod sections
Test points measured from the middle of the rod section
-25 -20 -15 -10 -5 0 5 10 15 20 25
The tests on the grids were done by inserting a small piece of rod in one of the spaces in the grid and pressing equally on both sides of the rod. The force on the grid was increase until it was significantly compressed; this was done because the grid springs could not take the full 75 N load. Each grid had four spaces tested and each space was tested three times.
Both the calibration and uncertainty test were done the same way as the rod section tests but with fewer test points.
4.1.2. Results
For the interpretation of the results MATLAB was used to plot the data and using built in functions calculate a trend line from the data points. The slope of the trend line together with equation from
A B C D E F G
Chapter 4.2.1 was then used to find the stiffness of the rod. Built in functions were used to calculate standard deviation from the stiffness values.
In Figure 8 the data points from section E can be seen. As can be seen in the figure at lower forces the slope is lower, this is because the measuring equipment settles into place. The phenomena could be seen in every one of the test. To remove the error this causes the data points at low forces were discarded. At the maximum force for a test the pistons direction is switched and because this takes time the data acquisition system (DAS) samples more data points there, to not get any errors from
weighting the top part heavier than the rest of the curve these points are also disregarded. The interpretation of the data from the grid tests was done using the same principles but the grids had a smaller linear region because at higher forces the springs lost strength. There was also a much larger spread of values for the stiffness when measuring at different points.
The data figures from the other sections and the grids can be found in Appendix B.
The result form the calibration test showed that the results should be corrected according to the following formula
𝐸𝐼𝑐𝑜𝑟𝑟𝑒𝑐𝑡𝑒𝑑= 1.148 𝐸𝐼 − 3.56 ( 25)
When this is done the results in Table 4 are gained and the noise in measurements for the equipment can be seen in Table 5
0 0.5 1 1.5 2 2.5 3 x 10-3 0 10 20 30 40 50 60 70 80 Force [N] Deflection [m]
Table 4: Results from the tests Part measured EI [Nm2] 𝜎 [Nm2] Section A 75.81 0.416 Section B 62.44 0.303 Section C 79.21 2.209 Section D 70.97 1.469 Section E 71.22 1.145 Section F 75.11 0.784 Section G 76.60 0.561 Spring constant �𝑚𝑁� 𝜎 �𝑚𝑁� MVG old 832800 371900 MVG new 879100 470900 Simple support 130900 88910
Table 5: Noise in measurements for the equipment
Force at the measurement [N] Noise dynamometer [N] Noise camera [mm]
0 ±0.12 ±0.015
50 ±0.22 ±0.005
78 ±0.30 ±0.006
Comparing the stiffness of the rod with the theoretical stiffness, se Figure 9a, it can be seen that the theory does not match the result of the experiment. It seem that the alumina does not contribute to the thickness as much as in the theory, in section A, B and G the rod has the same strength as if it was empty. To reflect this strength of alumina in the stiffness calculation was reduced by 48% and the Inconel by 2%, this was done to better fit the virtual rods to the real ones. The result can be seen in Figure 9b.
0.5 1 1.5 2 2.5 3 3.5 4 4.5 40 50 60 70 80 90 100 110 Stiffness [Nm2]
Distance from bottom plate [m]
0.5 1 1.5 2 2.5 3 3.5 4 4.5 40 50 60 70 80 90 100 110 Stiffness [Nm2]
Distance from bottom plate [m]
Figure 9: a) Theoretical stiffness and experimental results, b) Simulated rod stiffness and experimental rod stiffness. The red line is the stiffness found during the experiment and blue line is the calculated stiffness. The vertical blue lines are because in parts of the rod the stiffness is calculated without the alumina as explained in section 3.4.
B A
4.2. Benchmarking against COMSOL multiphysics
4.2.1. One material simulation
The second experiment, to find if the FEM model generates proper results, was to test it against COMSOLs beam and solid mechanics models to see if they gave similar results and also against the analytical solutions for simple beam problems.
The first batch of simulations were done with solid alumina beams 500 mm long with a diameter of 10mm and a uniformly distributed load of 50 N/m. They were simulated with the FEM model and both the beam model and solid mechanics model in COMSOL. The resulting deflection plots were compared and if the setup had an analytical solution it was added as well. Five different beam setups were tested and the results are discussed in the sections below and all the results can be in Table 6.
4.2.1.1. Simple clamped beam
Consisting of only a clamped beam it is simple enough to have an analytical solution. The FEM solution matches the analytical and beam solution but the solid mechanics solution deflects ~0.1 µm more and this is acceptable.
4.2.1.2. Clamped beam with support
A clamped beam with a support at 200 mm, the beam can rotate freely at the support but is looked in place. Again the FEM and beam calculations match and the solid mechanics are not far off.
4.2.1.3. Clamped beam with spring
The clamped beam has a spring attached at 250 mm, and the spring has a spring constant of 5E7 N/m. Once again the beam model matches with the solid mechanics one is a tiny bit off.
4.2.1.4. Clamped beam with prescribed deflection
A clamped beam but the lower end is shifted 0.01 mm. The results are consistent with earlier simulations.
4.2.1.5. Long combination beam
This 5 m long beam has a spring every 500 mm and the lower end is shifted 0.01 mm. The spring constants in the springs are 7E5 N/m. When simulating an entire rod larger differences can be noticed while still being acceptably small. It should be noted that this increase can in the solid mechanics case stem from the known problem that COMSOL has with simulating shapes where one dimension is much larger than the rest, i.e. a long rod. The setup can be seen in Figure 10 below.
Figure 10: Long beam experiment setup
5m
0.5
m
0.5
m
50 N/m
0.01 mm
offset
Table 6: Maximum deflections for the first batch of simulations, all values in meters
Case FEM model COMSOL beam COMSOL Solid
mechanics Analytical Bernoulli beam Simple clamped 5.5262 ∗ 10−5 5.5262 ∗ 10−5 5.5330 ∗ 10−5 5.5262 ∗ 10−5 Support 8.7930 ∗ 10−6 8.7930 ∗ 10−6 8.8401 ∗ 10−6 − Spring 3.5795 ∗ 10−6 3.5795 ∗ 10−6 3.6157 ∗ 10−6 − Prescribed deflection 6.0389 ∗ 10 −5 6.0389 ∗ 10−5 6.0458 ∗ 10−5 − Combination 9.5298 ∗ 10−5 9.5289 ∗ 10−5 10.379 ∗ 10−5 −
From the first batch of simulations it can be concluded that the FEM model can simulate simple beams with satisfactory accuracy
4.2.2. Two material simulations
The second set of simulations where done to see how well the FEM model handle more complicated beams, a composite rod was used with the inner part of alumina and an outer part of Inconel was used. The inner part has a diameter of 5 mm and the outer part 10 mm. As in the previous simulations the beams were all 500 mm long and clamped at both ends. The COMSOL solid
mechanics model simulates three different connection conditions between the alumina and Inconel, they can either be connected, have contact or have a frictionless contact. When calculating the stiffness of the rod for the FEM model both the equation from the theory and the one adjusted to fit the bending experiment in chapter 4.1 where used. The different cases and their results are
discussed below and maximum deflections can be seen in Table 7.
4.2.2.1. Solid Alumina
In the first simulation the inner beam was solid. The COMSOL values coincide well with the theory FEM values but not well with the adjusted FEM values.
4.2.2.2. One gap in the alumina
In this case the inner part consisted of two alumina parts with a 1 mm gap between them. The gap was centered at 250 mm from the end, the middle of the rod. Once again the theory FEM fits better with the COMSOL values but overall the COMSOL simulations have been more effected by the addition of the gap.
4.2.2.3. One cut in the alumina
This case differs from the one above that in the COMSOL simulations the two alumina pieces are next to each other as it would be in the real case. This cannot be done in the FEM model so the results will be identical to the previous case. In all the COMSOL simulation the deflection increased a small amount, less than 1%.
4.2.2.4. Four gaps in the alumina
Similar to the simulation above but with five alumina parts and the gap centered at 100, 200, 300 and 400 mm from the beginning. In this case COMSOL could not run the frictionless contact as the alumina pieces where not connected to anything. Same as before less change from adding three more gaps, this is because none of the gaps are in the middle or ends of the rod where the strain will be highest.
4.2.2.5. Five gaps in the alumina
Similar to the simulation above but with six alumina parts and the gap centered at 100, 200, 250, 300 and 400 mm from the beginning. As before COMSOL could not run the frictionless contact. Over all a small increase in deflection compared to the four gap case, the increase was of the same size as when going from solid to one gap. The exception is the COMSOL simulation with friction where the deflection decreases but this is most likely due to a computational error resulting in a shifting of the peak in x direction.
4.2.2.6. Five cuts in the alumina
The same situation as the one cut but the cuts are at 100, 200, 250, 300 and 400 mm. When comparing to the five gaps case the difference is small but in the simulation with friction it is larger than 1%.
4.2.2.7. Changing thickness
In the last simulation the diameter of the alumina increased until the middle of the rod and then decreased back to original size while the thickness of the Inconel changed to keep the same outer radius. The same thing that used to create variable power profiles in the ODEN rods. In the
simulation the alumina diameter increases linearly from 3 mm to 7 mm in the middle of the rod and the n back to 3 mm again. In this case COMSOL could only handle that the alumina and Inconel where connected all the other cases failed. In this case the difference between the values from the theory and adjusted FEM simulations is extra-large because the decrease in alumina strength from the adjustment makes the rods stiffness decrease as the alumina thickens, the opposite if what happens in the theory FEM and COMSOL simulations.
Table 7: Maximum deflections from the second batch of simulations, all values are in meters
Case FEM adjusted FEM theory COMSOL
connected COMSOL friction COMSOL frictionless Solid 8.3691 ∗ 10−5 7.8595 ∗ 10−5 7.8368 ∗ 10−5 7.8785 ∗ 10−5 7.8037 ∗ 10−5 One gap 8.3712 ∗ 10−5 7.8626 ∗ 10−5 7.7863 ∗ 10−5 7.8570 ∗ 10−5 7.8339 ∗ 10−5 One cut 8.3712 ∗ 10−5 7.8626 ∗ 10−5 7.8531 ∗ 10−5 7.9013 ∗ 10−5 7.9063 ∗ 10−5 Four gaps 8.3714 ∗ 10−5 7.8628 ∗ 10−5 7.7881 ∗ 10−5 7.8356 ∗ 10−5 − Five gaps 8.3735 ∗ 10−5 7.8659 ∗ 10−5 7.8567 ∗ 10−5 7.8069 ∗ 10−5 − Five cuts 8.3735 ∗ 10−5 7.8659 ∗ 10−5 7.8532 ∗ 10−5 7.9174 ∗ 10−5 − Changing 8.4114 ∗ 10−5 7.7895 ∗ 10−5 7.7881 ∗ 10−5 − −
From the second set of simulations it can be concluded that the FEM model can handle the more complex geometries of the rods compared to an ordinary beam. The difference in maximum
deflection compared the more advanced COMSOL calculations did not exceed 1%. The adjusted FEM model did not match the COMSOL simulation equally well, but this was expected as it was adjusted to fix a mismatch between the theory and real rods. When comparing against all the COMSOL simulations the standard model gave on average 0.18% more deflection with a standard deviation of ±0.52% and the one adjusted to fit the stiffness from the bending experiment gave 6.8% more deflection with ±0.64% deviation. The difference between the models will increase when a larger fraction of the cross section contains alumina because it is the part weakened in the adjustment. This is because COMSOL uses the stiffness calculations from theory.
5. Program
The program is built in MATLAB R2012b and it is comprised of six main files seen in Table 8. The program also reads data from two typed of files, .xlsx files that contain assembly information and .m files that contain geometrical information on a type of rod. The program is written so that the main program does not contain any test condition or geometrical data thus all the proprietary 2
information is kept in the data files.
Table 8: Files required for the ODEN rod bending program
Main program
ODEN_RodBending.m rodAssembly.m rod.m
FEM.m addCircle.m IAPWS_IF97.m
Data files
InputData.xlsx rodData_XX.m crookedRods.m or
uncertainRodPosition.m
5.1. What the files do
ODEN_RodBending.m is the file that starts the program and it contains the user interface and handles the other files. It creates the assembly and determines what calculations to do next. rodAssembly.m is a file that defines the class rodAssembly. This class stores the information about the assembly and it also creates the rods that are housed in the assembly. All calculations that are relevant for the entire assembly are done here.
rod.m is the file that defines the class rod. Same as the previous but handles the rods.
FEM.m is based on two file created by Antonio Ferreira (Ferreira, 2009)and it handles the all the FEM calculations needed to find the deflection of the rods.
addCircle.m adds a circle to the active graph; it is used to draw the rods.
IAPWS_IF97.m is a file that contains the IAPWSs formulas for calculating properties of water and steam. It was created by Mark Mikofski and is freely available at Matalb Central (The MathWorks inc., 2014). It is used to find the water properties during the temperature calculations.
InputData.xlsx contains all the properties of a certain test configuration such as the power, pressure and assembly geometry. The data is read by the program during the creation of the assembly. rodData_XX.m contains the geometrical data for the rod type XX and it is used every time a rod of this type is created.
crookedRods.m contains the information on how the rods deviate from their proper position when not under magnetic load.
uncertainRodPosition.m contains information like crookedRods.m but not as detailed.
5.2. Using the program
The first step in using the program is to create an .xlsx file that contains the assembly geometry and conditions in ODEN for the test that one wants to simulate. The file needs to follow the formatting of the existing InputData.xlsx file but can have any name. Secondly the user needs to make sure that
the types of rods that make up the bundle are represented by a rodData files otherwise one must be created.
When the user loads an input data file the program reads the test conditions and geometrical data from it and the rod data files that it needs. It saves all the data as a rodAssembly object that contains rod objects. When the data gathering step is done the initial calculations start. By using the rod resistance the current and power distribution is calculated and the power profile is set by reading if the rod name. From the power distribution and inlet temperature the water temperature and material temperature profile in the rod is calculated. From the material temperatures and the geometric data of the rod bending stiffness is found in each mm section of the rod.
At this point the user is presented with a graph of the cross section of the assembly and information about the rods and test parameters. If they chose to more information about the temperatures or the gap between the rods can be displayed.
When the user starts the deformation calculations the program creates the stiffness matrix and force vector then it solves for the deflection. This is done for each rod in turn and when completed the deflection for all rods is shown in a graph together with information on the change in gap between the rods. The user can again choose if more information about the gap, temperature or assembly shall be displayed.
For more information on how to use the program consult the manual in Appendix C and the code can be found in Appendix D.
0 1 2 3 4 5 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Deflection of rods, black lines show 10% and 25% of original min rod-rod dist. Vertical lines show BOHL and EOHL
Distance from the bottom [m]
D ef lec ti on [ m m ] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Deflection of rods, black lines show 10% and 25% of original min rod-rod dist. Vertical lines show BOHL and EOHL
D ef lec ti on [ m m ]
6. Results
Simulations with the program have shown that deflection can be reduced from current levels at the same time as the number of support grids is reduced. The tradeoff for a decreased maximum deflection is an increased average deflection. In the examined test, ODEN10.0, the reason for the large deflections was that distance between the top and bottom MVGs and their respective ends significantly larger than the distance between the grids. By optimizing the placement of the top and bottom support grids and removing the rest the change in minimum rod-rod distance could be reduced from 35.8% to 18.6%. The bending profile for the two cases can be seen in Figure 11.
B A
-20 -10 0 10 20 30
Cross-section with the smallets gap to neigbors in mm, and the change from the drawings gap length
2.899 -6.5% 2.898 -6.5% 2.899 -6.5% 2.917 -5.9% 2.942 -5.1% 2.917 -5.9% 2.917 -5.9% 2.939 -5.2% 2.917 -5.9% 2.898 -6.5% 2.941 -5.1% 2.898 -6.5% 2.897 -6.6% 2.897 -6.6% 2.917 -5.9% 2.943 -5.1% 2.917 -5.9% 2.917 -5.9% 2.938 -5.2% 2.917 -5.9% 2.899 -6.5% 2.940 -5.2% 2.899 -6.5% 2.898 -6.5% 2.898 -6.5% 2.918 -5.9% 2.942 -5.1% 2.918 -5.9% 2.917 -5.9% 2.939 -5.2% 2.917 -5.9% 2.897 -6.6% 2.940 -5.2% 2.897 -6.6% 2.897 -6.5% 2.897 -6.5% 2.917 -5.9% 2.940 -5.2% 2.917 -5.9% 2.918 -5.9% 2.942 -5.1% 2.898 -6.5% 2.918 -5.9% 2.943 -5.1% 2.942 -5.1% 2.943 -5.1% 2.799 -9.7% 2.799 -9.7% 2.939 -5.2% 2.799 -9.7% 2.799 -9.7% 1.411 -18.4% 2.941 -5.1% 2.943 -5.1% 2.799 -9.7% 2.799 -9.7% 2.938 -5.2% 2.799 -9.7% 2.801 -9.7% 1.408 -18.6% 2.940 -5.2% 2.942 -5.1% 2.801 -9.7% 2.801 -9.7% 2.939 -5.2% 2.801 -9.7% 2.801 -9.7% 1.408 -18.5% 2.940 -5.2% 2.940 -5.2% 2.801 -9.7% 2.801 -9.7% 2.942 -5.1% 2.799 -9.7% 2.801 -9.7% 1.408 -18.6% 1.411 -18.4% 1.408 -18.6% 1.408 -18.5% 1.408 -18.6% 2.435 18.8% 2.510 22.4% 2.529 23.4% 2.511 22.5% 2.437 18.9% 2.437 18.9% 2.510 22.4% 2.533 23.6% 2.511 22.5% 2.435 18.8% 2.532 2.435 18.8% 2.437 18.9% 2.511 22.5% 2.529 23.3% 2.509 22.4% 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 -40 -30 -20 -10 0 10 20 30 40 -30 -20 -10 0 10 20 30
Cross-section with the smallets gap to neigbors in mm, and the change from the drawings gap length
Position [mm] 2.711 -12.6% 2.710 -12.6% 2.711 -12.6% 2.746 -11.4% 2.792 -9.9% 2.746 -11.4% 2.745 -11.5% 2.785 -10.2% 2.745 -11.5% 2.708 -12.6% 2.792 -10.0% 2.708 -12.6% 2.707 -12.7% 2.707 -12.7% 2.746 -11.4% 2.795 -9.8% 2.746 -11.4% 2.745 -11.4% 2.784 -10.2% 2.745 -11.4% 2.711 -12.6% 2.788 -10.1% 2.711 -12.6% 2.709 -12.6% 2.709 -12.6% 2.747 -11.4% 2.793 -9.9% 2.747 -11.4% 2.746 -11.4% 2.786 -10.1% 2.746 -11.4% 2.707 -12.7% 2.790 -10.0% 2.707 -12.7% 2.708 -12.7% 2.708 -12.7% 2.745 -11.4% 2.788 -10.1% 2.745 -11.4% 2.747 -11.4% 2.792 -9.9% 2.710 -12.6% 2.747 -11.4% 2.795 -9.9% 2.792 -9.9% 2.795 -9.9% 2.519 -18.7% 2.519 -18.7% 2.785 -10.2% 2.519 -18.7% 2.518 -18.8% 1.115 -35.5% 2.792 -10.0% 2.795 -9.8% 2.518 -18.8% 2.518 -18.8% 2.784 -10.2% 2.518 -18.8% 2.522 -18.6% 1.109 -35.8% 2.788 -10.1% 2.793 -9.9% 2.522 -18.6% 2.522 -18.6% 2.786 -10.1% 2.522 -18.6% 2.522 -18.6% 1.110 -35.8% 2.790 -10.0% 2.788 -10.1% 2.522 -18.6% 2.522 -18.6% 2.792 -9.9% 2.519 -18.7% 2.522 -18.6% 1.109 -35.8% 1.115 -35.5% 1.109 -35.8% 1.110 -35.8% 1.109 -35.8% 2.796 36.4% 2.942 43.5% 2.978 45.3% 2.943 43.6% 2.798 36.5% 2.798 36.5% 2.940 43.4% 2.986 45.7% 2.942 43.5% 2.796 36.4% 2.796 36.4% 2.937 43.3% 2.984 45.5% 2.941 43.4% 2.798 36.5% 2.796 36.4% 2.798 36.5% 2.942 43.5% 2.978 45.2% 2.939 43.4% 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
The change in rod-rod distance given above is somewhat misleading in a DNB perspective as in Figure 12 it can be seen that this distance is calculated against the thimble rod. As the thimble rod does not generate heat being closer to it is less of a problem for determining CHF (Anglart, 2010). When not taking the thimble into account the change goes from 18.8% to 9.7%.
Another result gained from the simulations is that while the rods in the outer parts of the assembly, in Figure 12 rods 1-16, experience the most bending. It is the center rod that has the largest change in rod to rod distance. This is because when the outer rods bend inwards the rods they approach also move inwards and thus away from them. This is true for all rods except those approaching the center rod as the forces on it cancel thus it will not move away.
The difference in the change in minimum distance between rods is not is effected by the lowering of the stiffness to fit the experimental values but the difference between the cases vary depending on position along the rod. If the minimum distance is at the end of the rod, as in ODEN10.0, then the difference is small 1-2 percentage points but if the setup is such that the critical distance is in the middle then differences of 50 percentage points have been seen but this was when simulating unsupported rods. For a properly supported rod the difference is again small as when the deflection is minimized then all the deflection peaks are the same height, thus when switching to the stiffness model from mechanical theory the deflection in the middle will decrease substantially while at the ends it does not change by more than a couple of percent
6.1. Decreasing deflection
When all the deflection peaks are the same height the maximum deflection is at its lowest. The easiest way to do this is to have the main supports, ends of the test section and MVGs, at the same distance from each other as longer sections will both have a larger force on them and longer between supports. This problem is worsened by the free rotation at the grids resulting in a larger deflection of the longer section and smaller deflection at the shorter ones as seen in Figure 13from a COMSOL simulation.
If the MVG setup required results in a few large peaks, like in ODEN10.0, then adding secondary supports at the peaks is a good first step. The supports will likely have to be moved to get optimal effect.
If the deflection is too large when all the peaks are the same size then either more supports has to be added or the power has to be reduced.
Figure 13: Effect of different length between supports, blue line: no rotation, green line: free rotation, a) same length, b) longer, c) shorter B A C
7. Conclusion
A MATLAB program for simulating the deflection of the heater rods in the ODEN test loop has been constructed and according to the benchmarking against COMSOL the results is accurate within 1% with the model giving slightly more deflection than the benchmark. The main question of the reliability of the model comes from the big difference in the theoretical and real stiffness of the rods found during testing. While the difference should be further studied it does not affect the end results to a significant degree because of the positioning of the MVG.
The deflection of the rods during a CHF test is a manageable problem that can be further reduced by placing the support grids in the places where they have the greatest effect. When comparing to the setup of the most recent test in the ODEN loop more efficient placement of the supports could reduce the deflection from magnetic forces by half and at the same time cut the number of supports from 9 to 2.
7.2. Future work
Given the big differences between the theoretical and experimental bending stiffness of the rods the most important addition to the accuracy of the program is to do more experiments on rods to get a better certainty in that value. One possible answer to why the stiffness was lower is that the alumina pieces could have been damaged during use or handling thus reducing their strength. At the same time this is done more tests on the grids should be done to reduce the uncertainty of that value. Another improvement to the program could be to make the element length variable. This would make it possible to focus the computational power in more important regions by giving them more elements and in less critical regions fewer elements would speed up the calculations.
Bibliography
ABB Combustion Engineering Inc., 1992. Magnet, s.l.: s.n. Alferedsson, B., 2014. private communication. Stockholm: KTH.
Anglart, H., 2010. Thermal-Hydraulics in Nuclear Systems. Stockholm: KTH.
COMSOL Inc., 2014. Multiphysics Simulation Software - Platform for Physics-Based Modeling. [Online] Available at: http://www.comsol.com/products
[Accessed 12 May 2014].
Encyclopædia Britannica Online, 2014. magnetism (physics) -- Britannica Online Encyclopedia. [Online]
Available at: http://www.britannica.com.focus.lib.kth.se/EBchecked/topic/357334/magnetism [Accessed 2014 April 3].
Ferreira, A. M., 2009. MATLAB Codes for Finite Element Analysis. Waterloo: Springer.
IAPWS, I. A. f. t. P. o. W. a. S., 2007. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, Lucerne, Switzerland: s.n.
Lundh, H., 2000. Grundläggande hållfasthetslära. Stocholm: KTH.
MathWorks Inc., 2014. MATLAB - The Language of Technical Computing - MathWorks Nordic. [Online]
Available at: http://www.mathworks.se/products/matlab/ [Accessed 12 May 2014].
Miller, A. R., 2012. Lagrange shape Functions. [Online]
Available at: http://infohost.nmt.edu/~es421/ansys/shapefnt.htm [Accessed 26 May 2014].
National High Magnetic Field Laboratory , 2014. MagLab - Right and Left Hand Rules Tutorial. [Online]
Available at: http://www.magnet.fsu.edu/education/tutorials/java/handrules/ [Accessed 12 June 2014].
Sánchez-González, E. et al., 2007. Temperature dependence of mechanical properties of alumina up to the onset of creep. Journal of the European Ceramic Society, Issue 27, p. 3345–3349.
Subashki, G., 2014. IIT Academic Resource Center. [Online]
Available at: http://iit.edu/arc/workshops/pdfs/Georgi_Subashki_Workshop_Biot_Savart_Law.pdf [Accessed 26 May 2014].
Subramanian, S. R., n.d. Heat transfer in Flow Through Conduits, s.l.: Dep. of Chemical and Biomolecular Engineering, Clarkson University.
The MathWorks inc., 2014. IAPWS_IF97 functional form with no slip - File Exchange - MATLAB Central. [Online]
form-with-no-slip
[Accessed 4 February 2014].
Waldemarsson, F., 2013. Data Report EDF Cosine Critical Heat Flux Tests 1 and 2 at ODEN , s.l.: Westinghouse Electric.
One material simulation Simple clamped beam
Clamped beam with support
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6x 10 -5 Deflection [m] X coorinate [m] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 -1 0 1 2 3 4 5 6 7 8 9x 10 -6 Deflection [m] X coorinate [m]
Clamped beam with spring
Clamped beam with prescribed deflection
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.5 1 1.5 2 2.5 3 3.5 4x 10 -6 Deflection [m] X coorinate [m] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7x 10 -5 Deflection [m] X coorinate [m]
Long combination beam
Two material simulations Solid Alumina 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.2 0.4 0.6 0.8 1 1.2x 10 -4 Deflection [m] X coorinate [m] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 8 9x 10 -5 Deflection [m] X coorinate [m]
One gap in the alumina
One cut in the alumina
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 8 9x 10 -5 Deflection [m] X coorinate [m] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 8 9x 10 -5 Deflection [m] X coorinate [m]
Four gaps in the alumina
Five gaps in the alumina
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 8 9x 10 -5 Deflection [m] X coorinate [m] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 8 9x 10 -5 Deflection [m] X coorinate [m]
Five cuts in the alumina Changing thickness 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 8 9x 10 -5 Deflection [m] X coorinate [m] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 8 9x 10 -5 Deflection [m] X coorinate [m]
Section A Section B 0 0.5 1 1.5 2 2.5 3 x 10-3 0 10 20 30 40 50 60 70 80 Force [N] Deflection [m] 0 0.5 1 1.5 2 2.5 3 3.5 x 10-3 0 10 20 30 40 50 60 70 80 Force [N] Deflection [m]
Section C Section D 0 0.5 1 1.5 2 2.5 3 3.5 x 10-3 0 10 20 30 40 50 60 70 80 90 Force [N] Deflection [m] 0 0.5 1 1.5 2 2.5 3 x 10-3 0 10 20 30 40 50 60 70 80 Force [N] Deflection [m]
Section E Section F 0 0.5 1 1.5 2 2.5 3 x 10-3 0 10 20 30 40 50 60 70 80 Force [N] Deflection [m] 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10-3 0 10 20 30 40 50 60 70 80 90 100 110 Force [N] Deflection [m]
Section G MVG old 0 0.5 1 1.5 2 2.5 3 x 10-3 0 10 20 30 40 50 60 70 80 Force [N] Deflection [m] -5 -4 -3 -2 -1 0 x 10-4 0 10 20 30 40 50 60 70 80 90 100 Force [N] Deflection [m]
MVG new Simple support 10 -4 -3 -2 -1 0 x 10-4 0 10 20 30 40 50 60 70 80 90 100 Force [N] Deflection [m] -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 x 10-4 0 5 10 15 20 25 30 35 40 45 50 55 Force [N] Deflection [m]
Simple support 19 -7 -6 -5 -4 -3 -2 -1 0 1 x 10-4 0 5 10 15 20 25 30 35 40 45 50 55 Force [N] Deflection [m]
Manual for ODEN rod bending
This manual will give a detailed explanation of how to prepare, run and interpret the data from the ODEN rod bending program.
Preparing the Data file
The input data excel file is where all the data about the assembly and rods is set. When creating a new one start by editing the ProgramInputs_templeat.xlsx file and saving it under a different name to get the formatting correct.
It is important to fill in all the non-blue fields for the program to be able to run. The blue fields are optional as the data for these are loaded from the rodData.m files if they are empty. For a
description of the different fields see table 1 and 2.
In figure 4 an example data file is shown, this one does not overwrite the standard data from the rodData.m files as no blue fields are filled.
Care must also be taken when constructing the input data as the program will work even if the rods are placed inside each other or if the power would melt the rods.
Table 9: Assembly data fields
Field Description
Flow area Cross section area with water flow. [m2]
Length Length of the assembly from the topmost O-ring plates wetted side, known as the bottom plate, to the wetted side of the top adapter. [mm]
Pressure Pressure in the test loop during the test with maximum power. [bar] Power Maximum total power of the test in the test matrix. [kW]
Tin Temperature of the water at the inlet during the test with maximum power. [C]
Mass flow Mass flow of water through the test section inlet during the test with maximum power. [m/s]
Grid position … Distance from the top adapter to where the grids clamp the rods, se figure 1 for description. [mm]
More grids can be added by separating the positions with a space. Grid spring … Spring constant for grid radial movement. [N/m]
Standard value: 8.3E5 for MVG 1.3E5 for support grid
More grids can be added by separating the values with a space.
Wall Corner … X The position of the corners of the walls surrounding the rod bundle starting from the top left corner and going clockwise, in horizontal direction, se figure 2. [mm]
Separate values with spaces
Wall Corner … Y The position of the corners of the walls surrounding the rod bundle starting from the top left corner and going clockwise, in vertical direction, se figure 2. [mm]
Separate values with spaces
By selecting ‘No’ the rods will be calculated as straight.
When precise information is known selecting ‘Known’ will make the program read the crookedRods.m file to get the deflection in X and Y direction. If only the amount of bending and not in what direction the bend is pointing select ‘Uncertain’ and the program will read uncertainRodPosition.m and for all distains measurements will consider the rods “thicker” as if the bend pointed in all directions at the same time.
Table 10: Rod data fields, * marks required fields
Field Description
Name Rod ID, to easier keep track of the rods. Not used by the program. [~] Type * Type of rod, this is the same as the rod data file that all the automatic values
are loaded from. [~]
Select from drop down list of available rod types. Heated length Heated length. [mm]
X position * Position of the center of the rod in a cross section when viewed from the top locking down, in horizontal direction. [mm]
Y position * Position of the center of the rod in a cross section when viewed from the top locking down, in vertical direction. [mm]
Resistance … * Measured electrical resistance over the heated length. [Ω]
Ceramic section … Distance from the beginning of the heated length to where a gap in the ceramic that exist because two ceramic pieces meet. [mm]
More pieces can be added by separating the positions with a space. Outer diameter Outer diameter of the rod. [mm]
Inner diameter Inner diameter of the ceramic, the diameter of the space for the thermocouple wires. [mm]
Figure 14: The distance should be measured to the line as the
beginning of the heated length.
Running the program
To run the program first make sure that the required files are in the same folder, see table 3. Then Start MATLAB and open the ODEN_RodBending.m file. Click run to start the program, if MATLAB cannot find the file in the current folder then tell it to change folder, and wait a short while
Table 11: Required files
Main files Function
ODEN_RodBending.m Create and run the user interface
rodAssembly.m Keep data and do calculations on an assembly
rod.m Keep data and do calculations on a rod
FEM.m Solve the FEM calculations to find the deflection
IAPWS_IF97.m Calculate the water properties
addCircle.m Add a circle to a graphics area
The rodData files for the rod types used, ex.
rodData_14ft95mmCosineHot_WH600.m Contain geometrical data on the hot WH600 rod rodData_14ft95mmCosineCold_WH600.m Contain geometrical data on the cold WH600 rod rodData_122mmThimble.m Contain geometrical data in a thimble rod with
outer diameter of 12.2 mm
If the rods are crooked one of the following
crookedRods.m Contain information on how the rods in an
assembly are crooked, change from original position in X and Y direction
uncertainRodPosition.m Contain information on how the rods in an assembly are crooked, distance from original position
A window like the one in figure 5 will appear, consisting of a graph area (1), a table (2), a text area (3), display selection buttons (4) and control buttons (5). To continue press the Load data file button, circled in red in the figure. Then select the xlsx file containing your data and wait while the data loads and some preliminary calculations are performed.
Ones the calculations are complete the window graph area will show a cross section of the assembly geometry with rod numbers as viewed from the top looking down, check if the rods are positioned as expected. The table and text area will show rod and assembly data respectively, check if the data is correct. If the rods are crooked then the graph area will also show, in dashed circles, the cross section where the rod-rod/rod-wall distance is smallest. The next step is to check the temperature profile and to do that press the Temperature data button, it can be seen circled in red in figure 6 that also shows an example of the Assembly data window.
When the temperature data display is selected the graph area will show simulated the rod and water temperature depending on distance from the bottom plate. The text area will show the maximum temperature of the rods and water. Attention should be paid to that the rod temperature is not
above its melting point. The temperature calculations are accurate in the single-phase region and not as accurate in the two-phase region. Next step is to check the gaps between the rods and walls by pressing the Rod/Wall gap data button, seen circled in red in figure 7 that also shows an example of the Temperature data window.
When selecting the rod/wall gap data display the graph area will show the distance between the rods and their neighbors in the cross section with the smallest distance. The largest and smallest gaps are emphasized in bold text and unlike the assembly data view this one will only show where the rods actually are. An example of this window is shown in figure 8 and circled in red is the Calculate deflection button that continues to the next stage of the program.
When the calculations are complete the window will show a graph of the deflection of all the rods depending on position and the position of the grids are shown with dots. The vertical black lines show where the heated length begin and end and the horizontal black lines show the 10% and 25% of the drawing minimum rod-to-rod distance for reference. The table will show deflection data for all the rods and the text area will show the relevant extreme values. An example of this can be seen in figure 9.
If the view changed to assembly data or rod/wall gap data when the deflection has been calculated the graph will update to show situation during deflection. In assembly data view the table will now show the coordinates for the deflected rod instead of the original.
Analyzing the data
The main view for analyzing the results is the deflection data view as it shows the deflection graph and tells you the change in rod-rod distance. When calculating where the maximum change is the program will disregard everything below BOHL and the first grid because the rods are thicker there resulting in a large change in rod-rod distance from a small deflection and there is also low
probability of DNB in the lower part of the rod. The program will mark maximum deflections above the EOHL even though there will be no DNB there; this is to help find optimal place for a support grid. It should be noted that the deflection does not translate directly to rod-rod change because as the outer rods move inwards the next layer of rods will also move decreasing the change in rod gap. This can be seen in the assembly data view. The rod-rod change will also be higher when thimble rods are used as they are thicker and do not bend so the initial distance will be lower and the rod will not deflect away. But because the rod is not heated the rod-rod distance to its neighbors maters less for DNB purposes.
Reducing the deflection
Extreme bending peaks form because of differences in distance between the grids as the deflection correlates to the forth power of the distance between grids. This is worsened by the fact that the rods can rotate at the grids resulting in longer section gaining more deflection while neighboring sections lose deflection, and the reverse with shorter sections. An example of this can be seen in figure 5 where the last section is so long that is causes deflection in the opposite direction in its neighbor. In figures 3 the effect if different length between sections is shown with COMSOL beam
bending simulations. The blue line is when there is no rotation at the grid and the green is with free rotation the phenomena talked about above can be clearly seen.
Figure 3: Bending with/without rotation at the grids, blue line is without rotation green is with rotation. A) All sections equal length. B) Middle section longer. C) Middle section shorter.
A
B
The most effective way to reduce the deflection is to avoid any extreme peaks. This can be achieved by ether having the grids evenly spaced or add support grids at the point with the most deflection. Minor tweaking of the grid position because of peaks coming from that the grids are not rigid and of variations in rod stiffness will improve the deflections but general not by much.
When moving the grids around to improve the deflection one important note is that the position in the input data is measured from the top adapter to move the grids up this number should be decreased.
ur e 4: E xa mp le o f a d ata fi le
