# R&D on a Nonlinear Dynamics Analysis Code for the Drop Time of the Control Rod

Authors: Daogang Lu; Yuanpeng Wang; Qingyu Xie; Huimin Zhang; Muhammed Ali
Journal: Science and Technology of Nuclear Installations
 (2017)

---

## Abstract

Whether the control rod can drop down in time is one of the important guarantees for the safe operation of the nuclear power plant. The drop-down process of the control rod is very complicated. For a long time, the researchers have done a lot of work on that, but it is hard to consider all the nonlinear factors. This paper considers the main factors together. Based on the theoretical analysis, we developed the nonlinear dynamics response analysis software for the nuclear power plant, which can be used to calculate the rod’s drop-down time. Compared with the results of the experiments, the software we developed proves to be applicable and reliable.


---

## Body

## 1. Introduction

Reactivity control is one of the three main factors to ensure the safety of the reactor. The drop-down time of control rod is the key to the reactivity control. The time in which the control rod drops from its highest position by gravity fall into the buffer section of the guide tube is called the drop-down time of the control rod (Figure1). The drop-down time of control rod, which is a very important factor for safety, decides the rate of the reactor shutdown. The analysis of the control rod’s drop-down is very complex which is affected by temperature, load, flow field, pressure, and so forth. Now researchers mainly rely on experiments to determine the control rod’s drop-down time. They are also committed to the development of the software of the control rod’s drop-down analysis. Scholars have done a lot of work in this area.Figure 1
Simple model description for the control rod and the guide tube.Yu et al. [1] analyzed the falling process of control rod assembly, considering the resistance from fluid comprehensively. It researched on the free fall of control rod assembly in the absence of external excitation, without considering the impact under accident conditions. It also ignores the collision between control rod assembly and control rod guide tubes. Dou et al. [2] analyzed the entire process of the control rod’s drop-down. The author developed software for the control rod’s drop-down basing on his theory. What the author did has provided a strong reference for the safety analysis of reactor control rod drive line components. While dealing with the contacts collision between control rod assembly and the control rod guide tube, the author considered that the impact force is proportional to the size of the gap between them. The relationship between the impact force and the gap is complex actually. Sun and Wei [3] made a further analysis on the control rod’s drop-down time basing on the existing research. The author developed a finite element program of control rod’s drop-down, built a collision model, and could calculate fluid resistance not relying on the experiment, while the author only simplified the force of the fluid on the control rod as added mass, having made a limited analysis on the fluid-structural model, without taking into account the impact from the gap between control rod guide tube and control rod assembly.In summary, there is a certain discrepancy between the present calculation values and the experimental results. The reason is that it is hard to consider all the nonlinear factors. Therefore, this paper will make further improvement by considering all the main effects together in order to develop software to calculate the control rod’s drop-down time more accurately.


## 2. Theory and Model

### 2.1. Equations of Transverse Vibration

The control rod assembly has transverse vibration under flow excitation. We consider the single control rod and the guide tube as elastomers and establish an equation of transverse vibration [4], as shown in the following equation:(1)ρiVi+mai∂2uit∂t2+ci+cai∂uit∂t+kiuit=pit+FCit,where i is used to make a difference between control rod and control rod guide tube; ρ is the density; V is the volume; ma is the added mass; t is the time; u(t) is the transverse deflection of the structure at time t; c is the damping coefficient; ca is the added damping coefficient; k is the stiffness coefficient; p(t) is the external excitation at time t; F(t) is the impact force at time t. The transverse motion equations of the control rod and the control rod guide tube are coupled through the collision force between them.


### 2.2. Lateral Fluid Forces Analysis

Suppose what surrounds the control rod assembly is two-dimensional incompressible viscous fluid. The entire structure does forced transverse vibration, and no fluid flows along the axial direction of the structure.Fluid-structure interaction force includes fluid resistance and the inertia force. The vibration of the structure in fluid generates the fluid resistance. The fluid resistance can be divided into shape resistance and frictional resistance according to the direction of the force [5].The forward resistance caused by the pressure difference which is caused by the relative velocity of the particle and fluid is called shape resistance. Frictional resistance is a lateral shear stress due to edge effects arising from the roughness of the structure.Inertial force is caused by the acceleration of fluid particle. Inertial force includes Froude-Krylov force and added mass force. Froude-Krylov force can be considered as the force generated by the motion of the particle in the flow field when the structure does not exist. The added mass force may be regarded as the force that the structure needs to promote the fluid flow.Due to different boundary conditions, the force that the structure suffers is divided into two parts. First part (Fiff) represents that the structure is completely immersed in the infinite flow field and the second part (Fngf) represents that a portion of the structure boundary suffers the force of fluid between narrow gaps, as shown in the following [6]:(2)F=λ·Fiff+Fngf,where λ is the correction factor for the force Fiff. When the narrow gap does not exist, λ=1; and when the narrow gap exists, the correction coefficient λ is slightly less than 1, and the more the boundary surfaces suffer the force of fluid between narrow gaps, the smaller the value of λ is.Since the control rod aspect ratio is greater than 40, we only need to consider the fluid-structure interaction in the radial direction of the control rod assembly. The added damping of the control rods in the narrow gap model is as follows:(3)Ca=μBh2Bh+32L+λCd.The added mass is as follows: (4)Ma=ρB2B10h+1140L+λρ∀Ka,where μ is the fluid viscosity coefficient; B is the control rod diameter; h is the size of the gap between the control rod guide tubes and the control rod; L is the control rod length; λ is the dimensionless correction factor; Cd is the equivalent linear damping coefficient; ρ is the density of the liquid; ∀ is the control rod volume; Ka is the added mass coefficient. Because the narrow gap is too small compared with the control rod diameter, it plays an important role in the total force. We can even ignore the force Fiff and the dimensionless correction coefficient λ is close to zero.


### 2.3. Vertical Motion

Control rod’s drop-down experiments show that [2] we can introduce a delay coefficient of rod’s drop-down during the process of the control rod assembly falling from the highest position to the buffer section of the guide tube. The delay coefficient reflects the resistance from the fluid. It can generally be determined by experiment. It is related to the shape of the flow channel, the fluid pressure within the flow channel, temperature, motion, and other factors. Vertical movement of the control rod assembly drop satisfies the following relationship:(5)md2ztdt2+mαdztdt=mg-μn2∑1n1Fct-ρwVwg,where m is the total mass of the control rod assembly; z(t) is the control rod displacement in the vertical direction at time t; t is time; α is the delay coefficient; g is gravitational acceleration; μ is the friction coefficient between the control rod guide and the control rod; n1 is the number of the control rod nodes; n2 is number of the control rods; Fc(t) is the collision force at time t; ρw is the density of the fluid; Vw is volume of the control rod immersed in the fluid.


### 2.4. Collision Model

During earthquake conditions, neighboring structures are prone to have a collision, such as the collision between the control rod and the guide tube. This will cause huge damage to the structure. We need to establish a reasonable contact-collision model for analysis and research.Generic processing algorithm about contact problems includes nodes constraint method, parameter assignment method, and symmetric penalty function method. We used the symmetric penalty function method in this paper.Firstly, we should check whether each slave node goes through the active surface of the nodes. We do not deal with the nodes that do not go through the active surface. If there is a penetration, we bring a contact force called penalty function value between the slave node and the active surface. It is equivalent to placing a spring in one direction between the slave node and the active surface physically, to limit the penetration. This method is a relatively new method that can better handle the collision problems. It is widely used.We use the model Hertzdamp [7] to describe the collision between control rod assembly and the control rod guide tubes.Hertzdamp model (Figure2) is actually an improved model of Hertz model. It is, namely, the generalized Hertz model, and this model can take into account the energy loss during the contact-collision between two objects. It can better simulate the entire contact-collision dynamic process and has been widely used [8]. The model is the Hertz model in parallel with a nonlinear spring damper, to describe the energy loss during the contact-collision between two objects.(1)
During the stage of compression, the normal collision force and the normal plunge depth along the contact surface have the following relationship [9–12]:(6)Fc=kh·δn·1+34·1-e2·1-2·kh·δn+1n+1·δ0˙2·m1m2/m1+m2,δ≥1.(2)
During the stage of recovery, the normal collision force and the normal relative plunge depth along the contact surface have the following relationship:(7)Fc=kh·δn-ch·δ˙=kh·δn·1-34·1-e2·1-2·kh·δn+1n+1·δ0˙2·m1m2/m1+m2,δ≥1,where Fc is the collision force; kh is the Hertz model stiffness parameters; δ is the normal relative plunge depth along the contact surface between the control rod guide tube and the control rod; n is Hertz model coefficient, with n usually taken as 3/2; λ is the damping coefficient; e is the recovery coefficient; δ0˙ is the normal relative velocity along the contact surface when the collision occurs; m1 is the quality of the guide tube part of the collision that occurred; m2 is the quality of the control rod part of the collision that occurred.Figure 2
Hertzdamp model.


## 2.1. Equations of Transverse Vibration

The control rod assembly has transverse vibration under flow excitation. We consider the single control rod and the guide tube as elastomers and establish an equation of transverse vibration [4], as shown in the following equation:(1)ρiVi+mai∂2uit∂t2+ci+cai∂uit∂t+kiuit=pit+FCit,where i is used to make a difference between control rod and control rod guide tube; ρ is the density; V is the volume; ma is the added mass; t is the time; u(t) is the transverse deflection of the structure at time t; c is the damping coefficient; ca is the added damping coefficient; k is the stiffness coefficient; p(t) is the external excitation at time t; F(t) is the impact force at time t. The transverse motion equations of the control rod and the control rod guide tube are coupled through the collision force between them.


## 2.2. Lateral Fluid Forces Analysis

Suppose what surrounds the control rod assembly is two-dimensional incompressible viscous fluid. The entire structure does forced transverse vibration, and no fluid flows along the axial direction of the structure.Fluid-structure interaction force includes fluid resistance and the inertia force. The vibration of the structure in fluid generates the fluid resistance. The fluid resistance can be divided into shape resistance and frictional resistance according to the direction of the force [5].The forward resistance caused by the pressure difference which is caused by the relative velocity of the particle and fluid is called shape resistance. Frictional resistance is a lateral shear stress due to edge effects arising from the roughness of the structure.Inertial force is caused by the acceleration of fluid particle. Inertial force includes Froude-Krylov force and added mass force. Froude-Krylov force can be considered as the force generated by the motion of the particle in the flow field when the structure does not exist. The added mass force may be regarded as the force that the structure needs to promote the fluid flow.Due to different boundary conditions, the force that the structure suffers is divided into two parts. First part (Fiff) represents that the structure is completely immersed in the infinite flow field and the second part (Fngf) represents that a portion of the structure boundary suffers the force of fluid between narrow gaps, as shown in the following [6]:(2)F=λ·Fiff+Fngf,where λ is the correction factor for the force Fiff. When the narrow gap does not exist, λ=1; and when the narrow gap exists, the correction coefficient λ is slightly less than 1, and the more the boundary surfaces suffer the force of fluid between narrow gaps, the smaller the value of λ is.Since the control rod aspect ratio is greater than 40, we only need to consider the fluid-structure interaction in the radial direction of the control rod assembly. The added damping of the control rods in the narrow gap model is as follows:(3)Ca=μBh2Bh+32L+λCd.The added mass is as follows: (4)Ma=ρB2B10h+1140L+λρ∀Ka,where μ is the fluid viscosity coefficient; B is the control rod diameter; h is the size of the gap between the control rod guide tubes and the control rod; L is the control rod length; λ is the dimensionless correction factor; Cd is the equivalent linear damping coefficient; ρ is the density of the liquid; ∀ is the control rod volume; Ka is the added mass coefficient. Because the narrow gap is too small compared with the control rod diameter, it plays an important role in the total force. We can even ignore the force Fiff and the dimensionless correction coefficient λ is close to zero.


## 2.3. Vertical Motion

Control rod’s drop-down experiments show that [2] we can introduce a delay coefficient of rod’s drop-down during the process of the control rod assembly falling from the highest position to the buffer section of the guide tube. The delay coefficient reflects the resistance from the fluid. It can generally be determined by experiment. It is related to the shape of the flow channel, the fluid pressure within the flow channel, temperature, motion, and other factors. Vertical movement of the control rod assembly drop satisfies the following relationship:(5)md2ztdt2+mαdztdt=mg-μn2∑1n1Fct-ρwVwg,where m is the total mass of the control rod assembly; z(t) is the control rod displacement in the vertical direction at time t; t is time; α is the delay coefficient; g is gravitational acceleration; μ is the friction coefficient between the control rod guide and the control rod; n1 is the number of the control rod nodes; n2 is number of the control rods; Fc(t) is the collision force at time t; ρw is the density of the fluid; Vw is volume of the control rod immersed in the fluid.


## 2.4. Collision Model

During earthquake conditions, neighboring structures are prone to have a collision, such as the collision between the control rod and the guide tube. This will cause huge damage to the structure. We need to establish a reasonable contact-collision model for analysis and research.Generic processing algorithm about contact problems includes nodes constraint method, parameter assignment method, and symmetric penalty function method. We used the symmetric penalty function method in this paper.Firstly, we should check whether each slave node goes through the active surface of the nodes. We do not deal with the nodes that do not go through the active surface. If there is a penetration, we bring a contact force called penalty function value between the slave node and the active surface. It is equivalent to placing a spring in one direction between the slave node and the active surface physically, to limit the penetration. This method is a relatively new method that can better handle the collision problems. It is widely used.We use the model Hertzdamp [7] to describe the collision between control rod assembly and the control rod guide tubes.Hertzdamp model (Figure2) is actually an improved model of Hertz model. It is, namely, the generalized Hertz model, and this model can take into account the energy loss during the contact-collision between two objects. It can better simulate the entire contact-collision dynamic process and has been widely used [8]. The model is the Hertz model in parallel with a nonlinear spring damper, to describe the energy loss during the contact-collision between two objects.(1)
During the stage of compression, the normal collision force and the normal plunge depth along the contact surface have the following relationship [9–12]:(6)Fc=kh·δn·1+34·1-e2·1-2·kh·δn+1n+1·δ0˙2·m1m2/m1+m2,δ≥1.(2)
During the stage of recovery, the normal collision force and the normal relative plunge depth along the contact surface have the following relationship:(7)Fc=kh·δn-ch·δ˙=kh·δn·1-34·1-e2·1-2·kh·δn+1n+1·δ0˙2·m1m2/m1+m2,δ≥1,where Fc is the collision force; kh is the Hertz model stiffness parameters; δ is the normal relative plunge depth along the contact surface between the control rod guide tube and the control rod; n is Hertz model coefficient, with n usually taken as 3/2; λ is the damping coefficient; e is the recovery coefficient; δ0˙ is the normal relative velocity along the contact surface when the collision occurs; m1 is the quality of the guide tube part of the collision that occurred; m2 is the quality of the control rod part of the collision that occurred.Figure 2
Hertzdamp model.


## 3. Solution of the Equation

In this paper, we will use the step method [13] proposed by Newmark to perform numerical analysis. Newmark assumed the final velocity and displacement were as follows: (8)u˙1=u˙0+1-γΔtu¨0+γΔtu¨1(9)u1=u0+Δtu˙0+12-βΔt2u¨0+βΔt2u¨1.Coefficientγ indicates how the initial and the final acceleration change the speed. Coefficient β indicates how the initial and final acceleration change the displacement. When γ equal 1/2 and β equals 1/4, the method is called the Newmark-β average acceleration method. This can satisfy the engineering requirements, with second-order accuracy; when γ equal 1/2 and β equals 1/6, the method is called the linear acceleration method. Compared with the Newmark-β average acceleration method, the linear acceleration is only conditionally stable.In this paper, we will use Newmark-β average acceleration method and apply the basic assumption to the vibration equation. The second-order differential equation could be reduced to a first-order vibration equation. Then (1), (5), (3), (4), (6), and (7) can be solved simultaneously to get the collision force, rod’s drop-down time, and the displacement, velocity, and acceleration of the control rod from the results. Because the control rods are connected by a spider, we can assume the lateral movements of the control rods are the same to simplify the calculation. We can calculate the lateral vibration equations of one control rod to simulate the lateral vibration of the control rod assembly.


## 4. Development and Verification of the Code

We discretize the control rod and guide tube and use the finite element model. By developing a FORTRAN program based on the above analysis (Figure3 shows the program flow chart) we can do the calculation. When the control rod assembly is falling from the place of 2.8-meter height into the buffer section of the guide tube (about 0.48 m) under an external load excitation, it takes about 1.4540 s at the speed of 1.97 m/s. The displacement, velocity, acceleration, and friction of the control rod assembly are shown in Table 1.Table 1
Displacement, velocity, acceleration, and friction of the control rod at each time.
Time/s
Rod displacement/m
Velocity/(m⋅s−1)
Acceleration/(m⋅s−2)
Friction/N
0
0
0.02
8.64
0
0.2
0.12
0.99
3.51
26.73
0.4
0.38
1.55
0.57
42.03
0.6
0.71
1.76
0.93
19.85
0.8
1.07
1.83
0.14
29.79
1
1.45
1.93
0.03
24.52
1.2
1.83
1.92
1.36
0.09
1.4
2.22
1.91
−0.22
30.41
1.43
2.32
1.92
−0.06
27.93Figure 3
The program flow chart.Where Figure4 shows the relationship between the displacement of the control rod assembly and the rod’s drop-down time and Figure 5 shows the relationship between the velocity of the control rod assembly and the rod’s drop-down time, Figure 6 shows the relationship between the acceleration of the control rod assembly and the rod’s drop-down time. After research it is known that some scholars calculate the control rod’s drop-down time without considering the collision effect.Figure 4
Displacement time history of the control rod.Figure 5
Velocity time history of the control rod.Figure 6
Acceleration time history of the control rod.In addition, we also calculated without considering the collision. It was found that the rod’s drop-down time was 1.2820 s at the speed of 2.25 m/s. It means that the collision effect has a great impact on the calculation and thus cannot be ignored. The rod’s drop-down experiment shows that control rod assembly falls into the buffer section in 1.24 s and at the speed of 2.67 m/s. The results calculated by the program are basically consistent with the data from rod’s drop-down experiment [14, 15], with a deviation less than 10%. Therefore the computing method of this paper is reliable and conservative.


## 5. Conclusions

This paper analyzes the transverse vibration equation of elastomer under flow excitation and uses Newmark-β method to do numerical analysis, in addition to the fluid resistance that the elastomer suffers, the calculation method, and the effect of the gap on the fluid resistance. Hertzdamp improved model is applied to the control rod assembly and control rod guide tubes. Based on the theoretical analysis, we developed the nonlinear dynamics response analysis software for the nuclear power plant, which can be used to study rod’s drop-down time calculations. By calculating control rod’s drop-down instance, we can obtain the following conclusions:(1)
During the drop-down of the control rod in the power plant, there are several factors that affect the drop time. Different from front researches, we have considered the main factors and established a model. We put them into the program developed by ourselves.(2)
With this code, we can calculate the displacement, velocity, acceleration, and friction force of the control rod during its drop-down process.(3)
We found that the collision has large effect on the drop time of the control rod.



---

*Source: 1494385-2017-06-19.xml*