# Fretting Wear and Fatigue Life Analysis of Fuel Bundles Subjected to Turbulent Axial Flow in CEFR

Authors: Yafeng Shu; Jianjun Wu; Yongwei Yang; Zelong Zhao
Journal: Science and Technology of Nuclear Installations
 (2019)

---

## Abstract

In a fast spectrum reactor, the fuel rod bundle is mainly positioned radially by the wire which can make contact with the adjacent fuel rods, and then it is inevitable that flow-induced vibration (FIV) will cause fretting wear and vibration fatigue of the fuel cladding at the contact position. Therefore, the computational model of fretting wear and fatigue life about the fuel rod bundle caused by FIV will be studied in this paper. Based on the random vibration model of the fuel rod bundle, Hertz contact theory, and Archard wear theory, the fretting wear life computational model and the fatigue life computational model of the wire-to-adjacent fuel rod (WAFR) contact have been established. Finally, taking CEFR design parameters as an example, the fretting wear life and vibration fatigue life of the cladding are calculated, and it is found that fatigue affects the service life of the fuel rod more seriously than fretting wear. The calculation model and method lay a foundation for further accurate prediction and analysis of the fuel rod life.


---

## Body

## 1. Introduction

In a pressurized water reactor (PWR), the FIV causes grid-to-rod fretting wear, and the fretting wear is responsible for 70% of fuel cladding leakage accidents [1]. Similarly, there are some typical accidents caused by flow-induced vibration in fast spectrum reactors [2, 3]. In general, the fuel rods in fast reactors such as CEFR in China, FFTF in the United States, JOYO in Japan, and China initial Accelerator Driven System (CiADS) use the wire wrapping for radial positioning. With the change of fuel burnup under the reactor operation, pellet swelling and thermal expansion of the cladding will lead to the wire-to-adjacent fuel rod (WAFR) contact. Turbulent flow-induced vibration occurs when the coolant flows through the fuel rod bundle. Therefore, the turbulence-induced vibration is an important factor that causes the fretting wear and vibration fatigue of the cladding [4, 5]. In order to study the wear life and fatigue life of the contact between the wire and adjacent fuel rods, a theoretical computational model must be established. Bakewell, Clinch, and Chen and Wambsganss obtain pressure data and semiempirical pressure expressions [6–8]. According to the empirical formula of the pressure spectrum, the displacement power spectrum of the rod bundle can be obtained and the stress power spectrum can be calculated, which lays a foundation for the establishment of fretting wear life calculation and random vibration fatigue life calculation models. At present, there is little research on the complex life analysis of the contact between the wire and the fuel rod. In order to ensure the safe operation of the reactor, the theoretical computational model of the WAFR contact fretting wear and vibration fatigue under the random excitation of turbulence pulsation in the China Experimental Fast Reactor (CEFR) is studied in this paper.Firstly, the FIV dynamic model of the fuel rod bundle is established. According to the relationship between the pressure power spectrum and the displacement power spectrum, the displacement power spectrum of the fuel rod is obtained. Based on the Hertz contact theory, the WAFR contact stress model is established. Furthermore, the WAFR contact wear model is established by using Archard wear theory, and the fatigue life model of the WAFR contact is established, which is based on the random vibration fatigue theory. Finally, the theoretical computational model is used to predict the life of the contact danger points of fuel rods in CEFR.


## 2. Problem Description

In normal fast reactor component design, the initial installation gap is reserved between one fuel rod and another, as shown in Figure1. During reactor operation, the fuel cladding will undergo radial deformation because of thermal expansion and other factors, such as the error of installation clearance, irradiation swelling, and irradiation creep. About the effect of the error of installation clearance, we will use different equivalent spring stiffness to analyze the effect in Section 3.2. In general, the natural deformation which is caused by other factors for the cladding is mainly considered plastic after some operation time [9], and its contribution to the contact stress is very little to be considered so that the point contact between the wire and adjacent fuel rods occurs. Because of FIV of the turbulent pulsating pressure, the fretting wear and vibration fatigue will appear in the WAFR contact position. For the safe operation of the reactor, the fretting wear and vibration fatigue at the contact position must be calculated and studied. Table 1 shows the main parameters of the cold state design of fuel rods for the China Experimental Fast Reactor (CEFR). The following sections will take these parameters as examples for analysis.Figure 1
Illustration of contact between the wire wrap and the adjacent fuel rod. (a) Initial stateδ0. (b) Contact shape.
(a)
(b)Table 1
Design parameters of the CEFR [10, 11].
Rod diameter (mm)
Cladding thickness (mm)
Pitch (mm)
Wire wrap diameter
Wire pitch (mm)
Rod length (mm)
Cladding average temperature (°C)
Average velocity (m/s)
6.0
0.3
7.0
0.95
100
1350
500
4.51
Source:http://www.chinafra.com/newsitem/271381546.In the fast reactor, 316 stainless steel is regarded as the fuel cladding material (316SS). The thermal expansion coefficient is [12](1)αT=1.789×10−5+2.398×10−9T+3.269×10−13T2,where T is the cladding average temperature (K).For the sodium-cooled fast reactor CEFR, the cladding average temperature is about 500°C and the installation clearanceδ0 is 0–0.05 mm. Utilizing equation (1) to compute the radial thermal deformation of the cladding, Δ, it is found that Δ = 0.06 mm > δ0. As shown in Figure 1, the installation WAFR gap will close and make contact. Then, the stress concentration will occur at the contact point, which is the design danger point. The fretting wear and vibration fatigue of random vibration at dangerous points will be analyzed below.


## 3. Mathematical Model and Analysis Method

In general, the displacement response of the fuel rod is caused by turbulence excitation at the subcritical flow velocity (defined as the flow velocity before buckling and divergence). Therefore, in order to calculate the fretting model and vibration fatigue of the WAFR contact, the random vibration response model of the rod bundle caused by the turbulence of the fuel rod is firstly determined, and then the stress response spectrum is calculated through the contact stress and displacement response spectrum so as to establish the fretting wear life model to predict the failure life of the fuel rod. Finally, the fatigue life computational model of the WAFR contact is established through the stress response spectrum.


### 3.1. Vibration Response Model of Bundle

The WAFR contact stress is composed of static stress and dynamic stress, and turbulence pressure pulsation is the main factor causing dynamic stress change. In order to obtain the dynamic stress power spectral density, the dynamic equation governing the vibration motion of the fuel rod bundle can be expressed in the matrix form as [13, 14](2)Mq¨+Kq+Cq˙=F,where M, C, and K are, respectively, the mass, damping, and stiffness matrices and F is a generalized force matrix caused by turbulent pressure fluctuation. For a given system, the elements of the matrices are constant and the dimension is 2KN, where K is the number of cylinders and N is the Galerkin discrete order. Through the Fourier transform and frequency response function, the power spectral density (PSD) of fuel rod bundle displacement in equation (2) is expressed as [14](3)Sηiηjξ,ξ′,ω=RL3EI2Sppω∑l=1N∑n=1Nϕlξϕnξ′×∑p=12KN∑q=12KNH∗βl,i,pω×Hβn,j,qω×∫01dξ1dξ2ϕapϕaqSfbpfbqξ1,ξ2,ω,where Sηiηj(ξ, ξ′, ω) is defined as the cross-spectral density (CSD), H(ω) is the frequency response function and equal to [−ωM + iωC + K]−1, and ϕ(ξ) is the beam comparison function; here, the simply supported condition is applied. N is the discrete order, K is the number of rods, EI is the flexural rigidity, R is the rod radius, L is the length, a(p) = {largest integer ≤ (2K + p − 1)/2K}, b(p) = p − 2K[a(p) − 1], and β(n, j) = 2K(n − 1) + j. It is difficult and expensive to measure the pressure field of a wire-wrapped rod bundle. Therefore, numerical simulation can be used instead of the experiment to obtain the turbulent pressure spectrum. Abbasian et al. have done 43 fuel rods of CANDU reactor wall pressure test experiments and CFD simulation based on LES, DES, and RSM turbulence models. The comparative analysis shows that the numerical results are in good agreement with the experimental pressure spectrum, but the numerical simulation is time expensive [15]. And thus in this paper, the more accurate pressure spectrum formula, which is based on Curling and Païdoussis measured data of the bundle pressure field, is applied, and it is written as [16, 17](4)ΦPPΩ=2K2Rρ2U3e−2K1RΩ/U,RΩU≥0.2,2K3Rρ2U3RΩU,0≤RΩU≤0.2,where Ω is the circular frequency, U is the average axial flow velocity, and parameters K1, K2, and K3 are 0.25, 2.0 × 10−6, and 7.94 × 10−6. According to the relation Spp(ω) = [EI/(m + ρA)/L^4]0.5ΦPP(Ω), the displacement power spectrum of the arbitrary rod i can be calculated by substituting equation 4 into equation 3. In the fast reactor, the rods of the fuel assembly are arranged in a triangular arrangement, and for the convenience of calculation, only a three-cylinder system is selected for analysis, as shown in Figure 2. And the in-plane direction means the selected motion direction is in the same plane. The out-plane direction means the selected motion direction is in the parallel plane and perpendicular to the in-plane direction.Figure 2
Sketch of the three-cylinder system.Because of the very small amplitude, it is reasonable to assume that the change of contact depth between adjacent rods is caused by the displacement of one of the rods. The stress power spectrum can be obtained from the displacement power spectrumSx(ω). As shown in Figure 3, the theoretical PSD is in agreement with the experimental PSD, where the system parameters are the same as those in the experiment [18].Figure 3
Comparison of theoretical and experimental PSD of midcylinder displacements in a four-cylinder system in the radial direction. The scales on the left are for the first mode group, while those on the right are for the second group.


### 3.2. Contact Stress Model of WAFR Contact

The WAFR contact will inevitably produce contact stress. In order to obtain the stress power spectrum, it is necessary to establish the contact stress model. According to Hertz contact theory [19], for upper and lower surfaces, the contact area is generally an elliptical region, as shown in Figure 4.Figure 4
2D elliptical contact area.The relation between contact pressure and stress along the pressure direction is expressed as(5)σyp0=−be2a1−T2T,where a is the long half-axis of the ellipse, b is the short half-axis of the ellipse, eccentricity of the ellipse e = (1 − (b/a)2)1/2, and T = ((b2 + y2)/(a2 + y2))1/2. p0 is the maximum contact pressure (MPa). At y = 0, the maximum contact pressure is the maximum contact stress. The maximum contact pressure and maximum indentation depth, δ, are expressed as follows [19]:(6)p0=6PLE∗2π3Re21/3F1e−2/3,(7)δ=9PL216E∗2Re1/3F2e,from Figure 4.4 in Reference [19], we can get the correction factors F1(e) and F2(e). E∗ and Re are the equivalent Young’s modulus and the equivalent radius, respectively. They satisfy the relations(8)1E∗=1−ν12E1+1−ν22E2,(9)Re=12A⋅B−1/2,where E1, E2 and ν1, ν2 are Young’s modulus and Poisson’s ratio of the wire wrap and cladding. The relationship between positive real constants A and B is(10)A+B=121R1+1R1′+1R2+1R2′,B−A=121R1−1R1′2+1R2−1R2′2+21R1−1R1′1R2−1R2′cos2α1/2,where R1, R1′ and R2, R2′ are the maximum and minimum curvature radius of the wire wrap and fuel rods, respectively, and α is the slope angle between the wire wrap and the fuel rod axis (Figure 3). Curvature radius is expressed as(11)R1=2πR2+H24πR2,R1′=rw,R2=∞,R2′=R,where H is the helical pitch, R is the rod radius, and rw is the wire radius. The wire and cladding materials are 316L stainless steel. Table 2 shows the mechanical parameters of 316L. E∗ = 109.89 GPa is obtained from equation (8).Table 2
Main mechanical properties of 316L [20].
Mechanical properties
E (GPa)
ν
σ
0.2 (MPa)
σ
b (MPa)
HV
316L
206
0.3
285
555
135The equivalent spring stiffness at the contact pointKstiff = 1.42 N/mm [21], initial displacement h∼10 μm, and then the normal initial load is expressed as(12)PL=Kstiff∗h.Considering turbulent excitation, the normal load is a function of rod displacement; then,(13)PL=Kstiffh−xt,where x(t) is the rod displacement caused by turbulent excitation.Substituting equation13 into equation (6), and carrying out the first-order Taylor expansion at x(t) = 0, the formula is written as(14)p0=6E∗2π3Re21/3F1e−2/3Kstiff1/3h1/3+13h−2/3xt.Therefore, the stress power spectral density corresponding to equation (14) is(15)Sσω=6E∗2π3Re22/3F1e−4/3Kstiff219h−4/3Sxω.The ratio of the principal relative radii of the curvature and the ratio of the long half-axis and short half-axis are defined as [19](16)R′R″=BA,(17)ba≈R′R″−2/3.According to the relationship between load and pressure distribution on the semiellipsoid, the maximum pressure can be expressed as(18)p0=32PLπab.Combining equations (17) and (18), the long half-axis a and the short half-axis b can be determined. When the maximum static stress σymax exceeds the yield strength of the material, as shown in Table 2, plastic deformation occurs at the contact position until the contact area can resist external forces.In order to verify the correctness of the contact stress model in this paper, we use the contact stress model to compute the contact problem as in Reference [22], as shown in Figure 5, and it is found that both results are in good agreement. This example verifies the correctness of the model in this paper.Figure 5
Relation between the contact pressure and the contact deformation.


### 3.3. Wear Life Computational Model

In the previous section, the calculation of the WAFR contact stress has been analyzed. However, turbulence-induced vibration will cause the WAFR contact wear. In engineering, the maximum wear depth of the fuel rod cladding has to be less than 10% cladding thickness [23]. Therefore, the fretting wear analysis of the cladding is very important for reactor safety. Based on the Archard wear formula, written as equation (19), the fretting wear life calculation model of the cladding will be established in this section:(19)V=KAA0PLl,where A is the WAFR wear contact area, A0 is the reference contact area, and l is the total sliding distance. According to the experiment, the wear coefficient of 316L stainless steel for a sliding speed of 0.06 m/s at a frequency of 2 Hz is K (mm3/N/m), i.e., the slope of the fitting straight line, as shown in Figure 6. FN is the normal load on the contact surface, Rd is the roughness of the abrasive paper, and the diameter of the reference contact surface D0 = 8 mm in Reference [24]. The WAFR wear contact area A is approximately a semiellipsoid surface, and the formula is written as(20)A=23πab+bδ+aδ,where a, b, and c are the half-axis lengths of the contact semiellipsoid, respectively.Figure 6
Wear rate as a function of normal load.For the turbulent random vibration, the full-contact sliding friction distance of the wire on the cladding surface per second,s, is defined. According to the displacement power spectral density, the mean square root drrms of the wire displacement amplitude in the compression direction can be obtained.As shown in Figure7, from equation (13), the normal load of the WAFR contact is expressed as(21)P′=KstiffR2−srms2sin2πfavet2−R−drms.Figure 7
Sliding wear.The contact trajectory is written as(22)s=Rarcsinsrmssin2πfavetR.By deriving formula (22) with respect to time, the relative friction rate between the wire and the cladding can be obtained as follows:(23)s˙=2πRsrmsfavecos2πfavetR2−srms2sin22πfavet.From equation (19), we can obtain the wear volume of T0/4 time, and it is written as(24)VT0/4=∫0T0/4KAA0P′s˙dt,where the average period T0 = 1/fave. The average frequency fave, that is, the average number of times the positive slope crosses the mean value in unit time, is consistent with the definition in Section 3.4:(25)fave=∫−∞+∞f2Sfdf∫−∞+∞Sfdf1/2.AfterNc cycles, as shown in Figure 7, the wear morphology can be determined according to the sliding path of the wire on the cladding surface. The wear volume can be divided into three parts: two quarter ellipsoids and a body whose cross section is a semiellipse rotating at a certain angle 2arcsin (srms/R). And thus, the total wear volume is written as(26)V=4NcVT0/4=2V1+V2=23πabc+12πacR⋅2arcsinsrmsR.Then, the maximum wear depthc is expressed as(27)c=V2/3πab+πaRarcsinsrms/R≤10%h.The longest allowable life of cladding wear can be calculated by equation (27), as follows:(28)t≤10%h2/3πab+πaRarcsinsrms/RT0/4VT0/43600∗24.Finally, according to equation (28), the wear life of the cladding can be determined in units of hours.In order to verify the correctness of the wear model in this paper, we use Wang et al.’s example [9] which calculates the contact wear of the fuel rod to the grid. As shown in Figure 8, the selected parameters are in Reference [9]; from the comparison of both results, there is good agreement between the established model in this paper and the Wang computational model.Figure 8
Validity of the fretting wear life computational model.


### 3.4. Random Vibration Fatigue Life Computational Model

The turbulence excitation will generate random contact stress at the contact danger points, which will lead to the cladding fatigue. The following will establish the vibration fatigue model caused by random contact stress based on the statistical average and probability function method. The random stress PSDW(f) can be calculated through equation (15). The peak probability density function of the stationary stochastic process can be written as [25](29)Pp=1σ1−α2πe−1/21−αx/σ2+xα2σ21+erfxσα21−αe−1/2x/σ2,where erf(x) is the error function and σ is the standard deviation, which means the root mean square value of stress σrms, and the random signal irregular factor α is expressed as(30)α=∫0∞f2Wfdf2σ2∫0∞f4Wfdf.For narrowband random signals,α = 1, the probability density function of stress amplitude obeys Rayleigh distribution. However, for broadband random signals, the amplitude probability density function is the combination of Rayleigh distribution and Gauss distribution, and especially α = 0; it obeys Gauss distribution.The mean square values of the continuous spectrum and the corresponding discrete spectrum are expressed as follows:(31)σ2=∫−∞+∞Wfdf,(32)σ2=∑−∞+∞WfiΔfi.For the discrete spectrum, the zero-crossing expectation of stressSσ(t) is expressed as equation (33), that is, the average number of times during unit time passing through the mean value at a positive slope:(33)f0=∑−∞+∞fi2WfiΔfi∑−∞+∞WfiΔfi1/2.For the discrete spectrum, the expected peak value of stressSσ(t), i.e., the average number of peak occurrences per unit time, is written as(34)n0=∑−∞+∞fi4WfiΔfi∑−∞+∞fi2WfiΔfi1/2.Then, according to equation (30), the irregular factor of random signals is(35)α=f0n0,and the bandwidth coefficient is defined as ε = (1 − α2)0.5. Two parameters, α and ε, characterize the irregularity and bandwidth of stochastic processes, respectively. When α = 1, i.e., ε = 0, it is narrowband. It is generally believed in engineering that when ε ≤ 0.35, the random process can be approximated as a narrowband treatment [26].According to Miner’s linear cumulative damage theory, the cumulative damage of the material can be expressed as(36)Di=niNi.From equation (36), the following structural failure criteria can be established:(37)∫0∞Dids=1;the discrete formD=∑i=1kniNi=1,where the subscript i represents the ith operating condition, ni denotes the number of true cycles when the peak stress value is s, and Ni denotes the number of cycles for material failure when the peak stress is s.The S-N curve describing the material is generally modeled by power functions as follows:(38)simNi=C,where C and m are parameters related to the material, stress ratio, loading, and so on.The Basquin model of the stress control conditions can be written as [27](39)σa=σf′2Nib.In order to consider the average stress effect, the SWT parameter method introduces equivalent stress as follows [28]:(40)σaeq=σmaxσa.Combining equations (39) and (40), the fatigue life under any type i working condition is written as(41)Ni=12σaeqσf′1/b.Letm = −1/b, C = 1/2(σf′)−1/b, and si = σaeq. That is, equation (41) is written in the form of equation (39), where the strength coefficient σf′ of 316 stainless steel is 2100 MPa and the fatigue strength index b is −0.1 [29], and then substituting equation (41) into equation (37) leads to(42)∑i=1kniCsim=1.AssumingNt is the total number of cycles, the number of cycles in the s to s + ds interval is(43)dni=NtPSads.By combining equation (42) with equation (43), the total cycles of material failure are written as(44)Nt=1∫0∞PSasm/Cds.For the narrowband process, from equation (19), the peak stress probability density function is(45)Ps=s2σ2e−s2/2σ2,s≥0,0,s<0.Substituting equation (45) into equation (44), the gross number of stress cycles is obtained as(46)Nt=12mCΓ2+m2σ−m,where Γ(x) is the gamma function.In practical engineering, most of them are broadband random vibration. For broadband stochastic processes, Wirsching and Light introduce a correction factorλ as [30](47)λ=1−am1−εbm+am.The values ofa and b in equation (47) are related to the material parameter m, which can be expressed as(48)am=0.926−0.033m,bm=1.587m−2.323.Thus, a life estimation formula of the broadband stochastic process is obtained as(49)Nt=1λ12mCΓ2+m2σ−m.According to the value of the bandwidth coefficientε, the random vibration fatigue life of the structure can be determined by using equation (46) or (49).In order to verify the correctness of the random vibration fatigue life model, we use the model to compare with the experimental data [31], as shown in Figure 9, and it is found that most of the data are in the double error band and very few data are within the quadruple error band. The life prediction errors are acceptable in engineering, and thus, the applicability of the model is verified.Figure 9
Comparison between the fatigue life model and the experimental data [31].


## 3.1. Vibration Response Model of Bundle

The WAFR contact stress is composed of static stress and dynamic stress, and turbulence pressure pulsation is the main factor causing dynamic stress change. In order to obtain the dynamic stress power spectral density, the dynamic equation governing the vibration motion of the fuel rod bundle can be expressed in the matrix form as [13, 14](2)Mq¨+Kq+Cq˙=F,where M, C, and K are, respectively, the mass, damping, and stiffness matrices and F is a generalized force matrix caused by turbulent pressure fluctuation. For a given system, the elements of the matrices are constant and the dimension is 2KN, where K is the number of cylinders and N is the Galerkin discrete order. Through the Fourier transform and frequency response function, the power spectral density (PSD) of fuel rod bundle displacement in equation (2) is expressed as [14](3)Sηiηjξ,ξ′,ω=RL3EI2Sppω∑l=1N∑n=1Nϕlξϕnξ′×∑p=12KN∑q=12KNH∗βl,i,pω×Hβn,j,qω×∫01dξ1dξ2ϕapϕaqSfbpfbqξ1,ξ2,ω,where Sηiηj(ξ, ξ′, ω) is defined as the cross-spectral density (CSD), H(ω) is the frequency response function and equal to [−ωM + iωC + K]−1, and ϕ(ξ) is the beam comparison function; here, the simply supported condition is applied. N is the discrete order, K is the number of rods, EI is the flexural rigidity, R is the rod radius, L is the length, a(p) = {largest integer ≤ (2K + p − 1)/2K}, b(p) = p − 2K[a(p) − 1], and β(n, j) = 2K(n − 1) + j. It is difficult and expensive to measure the pressure field of a wire-wrapped rod bundle. Therefore, numerical simulation can be used instead of the experiment to obtain the turbulent pressure spectrum. Abbasian et al. have done 43 fuel rods of CANDU reactor wall pressure test experiments and CFD simulation based on LES, DES, and RSM turbulence models. The comparative analysis shows that the numerical results are in good agreement with the experimental pressure spectrum, but the numerical simulation is time expensive [15]. And thus in this paper, the more accurate pressure spectrum formula, which is based on Curling and Païdoussis measured data of the bundle pressure field, is applied, and it is written as [16, 17](4)ΦPPΩ=2K2Rρ2U3e−2K1RΩ/U,RΩU≥0.2,2K3Rρ2U3RΩU,0≤RΩU≤0.2,where Ω is the circular frequency, U is the average axial flow velocity, and parameters K1, K2, and K3 are 0.25, 2.0 × 10−6, and 7.94 × 10−6. According to the relation Spp(ω) = [EI/(m + ρA)/L^4]0.5ΦPP(Ω), the displacement power spectrum of the arbitrary rod i can be calculated by substituting equation 4 into equation 3. In the fast reactor, the rods of the fuel assembly are arranged in a triangular arrangement, and for the convenience of calculation, only a three-cylinder system is selected for analysis, as shown in Figure 2. And the in-plane direction means the selected motion direction is in the same plane. The out-plane direction means the selected motion direction is in the parallel plane and perpendicular to the in-plane direction.Figure 2
Sketch of the three-cylinder system.Because of the very small amplitude, it is reasonable to assume that the change of contact depth between adjacent rods is caused by the displacement of one of the rods. The stress power spectrum can be obtained from the displacement power spectrumSx(ω). As shown in Figure 3, the theoretical PSD is in agreement with the experimental PSD, where the system parameters are the same as those in the experiment [18].Figure 3
Comparison of theoretical and experimental PSD of midcylinder displacements in a four-cylinder system in the radial direction. The scales on the left are for the first mode group, while those on the right are for the second group.


## 3.2. Contact Stress Model of WAFR Contact

The WAFR contact will inevitably produce contact stress. In order to obtain the stress power spectrum, it is necessary to establish the contact stress model. According to Hertz contact theory [19], for upper and lower surfaces, the contact area is generally an elliptical region, as shown in Figure 4.Figure 4
2D elliptical contact area.The relation between contact pressure and stress along the pressure direction is expressed as(5)σyp0=−be2a1−T2T,where a is the long half-axis of the ellipse, b is the short half-axis of the ellipse, eccentricity of the ellipse e = (1 − (b/a)2)1/2, and T = ((b2 + y2)/(a2 + y2))1/2. p0 is the maximum contact pressure (MPa). At y = 0, the maximum contact pressure is the maximum contact stress. The maximum contact pressure and maximum indentation depth, δ, are expressed as follows [19]:(6)p0=6PLE∗2π3Re21/3F1e−2/3,(7)δ=9PL216E∗2Re1/3F2e,from Figure 4.4 in Reference [19], we can get the correction factors F1(e) and F2(e). E∗ and Re are the equivalent Young’s modulus and the equivalent radius, respectively. They satisfy the relations(8)1E∗=1−ν12E1+1−ν22E2,(9)Re=12A⋅B−1/2,where E1, E2 and ν1, ν2 are Young’s modulus and Poisson’s ratio of the wire wrap and cladding. The relationship between positive real constants A and B is(10)A+B=121R1+1R1′+1R2+1R2′,B−A=121R1−1R1′2+1R2−1R2′2+21R1−1R1′1R2−1R2′cos2α1/2,where R1, R1′ and R2, R2′ are the maximum and minimum curvature radius of the wire wrap and fuel rods, respectively, and α is the slope angle between the wire wrap and the fuel rod axis (Figure 3). Curvature radius is expressed as(11)R1=2πR2+H24πR2,R1′=rw,R2=∞,R2′=R,where H is the helical pitch, R is the rod radius, and rw is the wire radius. The wire and cladding materials are 316L stainless steel. Table 2 shows the mechanical parameters of 316L. E∗ = 109.89 GPa is obtained from equation (8).Table 2
Main mechanical properties of 316L [20].
Mechanical properties
E (GPa)
ν
σ
0.2 (MPa)
σ
b (MPa)
HV
316L
206
0.3
285
555
135The equivalent spring stiffness at the contact pointKstiff = 1.42 N/mm [21], initial displacement h∼10 μm, and then the normal initial load is expressed as(12)PL=Kstiff∗h.Considering turbulent excitation, the normal load is a function of rod displacement; then,(13)PL=Kstiffh−xt,where x(t) is the rod displacement caused by turbulent excitation.Substituting equation13 into equation (6), and carrying out the first-order Taylor expansion at x(t) = 0, the formula is written as(14)p0=6E∗2π3Re21/3F1e−2/3Kstiff1/3h1/3+13h−2/3xt.Therefore, the stress power spectral density corresponding to equation (14) is(15)Sσω=6E∗2π3Re22/3F1e−4/3Kstiff219h−4/3Sxω.The ratio of the principal relative radii of the curvature and the ratio of the long half-axis and short half-axis are defined as [19](16)R′R″=BA,(17)ba≈R′R″−2/3.According to the relationship between load and pressure distribution on the semiellipsoid, the maximum pressure can be expressed as(18)p0=32PLπab.Combining equations (17) and (18), the long half-axis a and the short half-axis b can be determined. When the maximum static stress σymax exceeds the yield strength of the material, as shown in Table 2, plastic deformation occurs at the contact position until the contact area can resist external forces.In order to verify the correctness of the contact stress model in this paper, we use the contact stress model to compute the contact problem as in Reference [22], as shown in Figure 5, and it is found that both results are in good agreement. This example verifies the correctness of the model in this paper.Figure 5
Relation between the contact pressure and the contact deformation.


## 3.3. Wear Life Computational Model

In the previous section, the calculation of the WAFR contact stress has been analyzed. However, turbulence-induced vibration will cause the WAFR contact wear. In engineering, the maximum wear depth of the fuel rod cladding has to be less than 10% cladding thickness [23]. Therefore, the fretting wear analysis of the cladding is very important for reactor safety. Based on the Archard wear formula, written as equation (19), the fretting wear life calculation model of the cladding will be established in this section:(19)V=KAA0PLl,where A is the WAFR wear contact area, A0 is the reference contact area, and l is the total sliding distance. According to the experiment, the wear coefficient of 316L stainless steel for a sliding speed of 0.06 m/s at a frequency of 2 Hz is K (mm3/N/m), i.e., the slope of the fitting straight line, as shown in Figure 6. FN is the normal load on the contact surface, Rd is the roughness of the abrasive paper, and the diameter of the reference contact surface D0 = 8 mm in Reference [24]. The WAFR wear contact area A is approximately a semiellipsoid surface, and the formula is written as(20)A=23πab+bδ+aδ,where a, b, and c are the half-axis lengths of the contact semiellipsoid, respectively.Figure 6
Wear rate as a function of normal load.For the turbulent random vibration, the full-contact sliding friction distance of the wire on the cladding surface per second,s, is defined. According to the displacement power spectral density, the mean square root drrms of the wire displacement amplitude in the compression direction can be obtained.As shown in Figure7, from equation (13), the normal load of the WAFR contact is expressed as(21)P′=KstiffR2−srms2sin2πfavet2−R−drms.Figure 7
Sliding wear.The contact trajectory is written as(22)s=Rarcsinsrmssin2πfavetR.By deriving formula (22) with respect to time, the relative friction rate between the wire and the cladding can be obtained as follows:(23)s˙=2πRsrmsfavecos2πfavetR2−srms2sin22πfavet.From equation (19), we can obtain the wear volume of T0/4 time, and it is written as(24)VT0/4=∫0T0/4KAA0P′s˙dt,where the average period T0 = 1/fave. The average frequency fave, that is, the average number of times the positive slope crosses the mean value in unit time, is consistent with the definition in Section 3.4:(25)fave=∫−∞+∞f2Sfdf∫−∞+∞Sfdf1/2.AfterNc cycles, as shown in Figure 7, the wear morphology can be determined according to the sliding path of the wire on the cladding surface. The wear volume can be divided into three parts: two quarter ellipsoids and a body whose cross section is a semiellipse rotating at a certain angle 2arcsin (srms/R). And thus, the total wear volume is written as(26)V=4NcVT0/4=2V1+V2=23πabc+12πacR⋅2arcsinsrmsR.Then, the maximum wear depthc is expressed as(27)c=V2/3πab+πaRarcsinsrms/R≤10%h.The longest allowable life of cladding wear can be calculated by equation (27), as follows:(28)t≤10%h2/3πab+πaRarcsinsrms/RT0/4VT0/43600∗24.Finally, according to equation (28), the wear life of the cladding can be determined in units of hours.In order to verify the correctness of the wear model in this paper, we use Wang et al.’s example [9] which calculates the contact wear of the fuel rod to the grid. As shown in Figure 8, the selected parameters are in Reference [9]; from the comparison of both results, there is good agreement between the established model in this paper and the Wang computational model.Figure 8
Validity of the fretting wear life computational model.


## 3.4. Random Vibration Fatigue Life Computational Model

The turbulence excitation will generate random contact stress at the contact danger points, which will lead to the cladding fatigue. The following will establish the vibration fatigue model caused by random contact stress based on the statistical average and probability function method. The random stress PSDW(f) can be calculated through equation (15). The peak probability density function of the stationary stochastic process can be written as [25](29)Pp=1σ1−α2πe−1/21−αx/σ2+xα2σ21+erfxσα21−αe−1/2x/σ2,where erf(x) is the error function and σ is the standard deviation, which means the root mean square value of stress σrms, and the random signal irregular factor α is expressed as(30)α=∫0∞f2Wfdf2σ2∫0∞f4Wfdf.For narrowband random signals,α = 1, the probability density function of stress amplitude obeys Rayleigh distribution. However, for broadband random signals, the amplitude probability density function is the combination of Rayleigh distribution and Gauss distribution, and especially α = 0; it obeys Gauss distribution.The mean square values of the continuous spectrum and the corresponding discrete spectrum are expressed as follows:(31)σ2=∫−∞+∞Wfdf,(32)σ2=∑−∞+∞WfiΔfi.For the discrete spectrum, the zero-crossing expectation of stressSσ(t) is expressed as equation (33), that is, the average number of times during unit time passing through the mean value at a positive slope:(33)f0=∑−∞+∞fi2WfiΔfi∑−∞+∞WfiΔfi1/2.For the discrete spectrum, the expected peak value of stressSσ(t), i.e., the average number of peak occurrences per unit time, is written as(34)n0=∑−∞+∞fi4WfiΔfi∑−∞+∞fi2WfiΔfi1/2.Then, according to equation (30), the irregular factor of random signals is(35)α=f0n0,and the bandwidth coefficient is defined as ε = (1 − α2)0.5. Two parameters, α and ε, characterize the irregularity and bandwidth of stochastic processes, respectively. When α = 1, i.e., ε = 0, it is narrowband. It is generally believed in engineering that when ε ≤ 0.35, the random process can be approximated as a narrowband treatment [26].According to Miner’s linear cumulative damage theory, the cumulative damage of the material can be expressed as(36)Di=niNi.From equation (36), the following structural failure criteria can be established:(37)∫0∞Dids=1;the discrete formD=∑i=1kniNi=1,where the subscript i represents the ith operating condition, ni denotes the number of true cycles when the peak stress value is s, and Ni denotes the number of cycles for material failure when the peak stress is s.The S-N curve describing the material is generally modeled by power functions as follows:(38)simNi=C,where C and m are parameters related to the material, stress ratio, loading, and so on.The Basquin model of the stress control conditions can be written as [27](39)σa=σf′2Nib.In order to consider the average stress effect, the SWT parameter method introduces equivalent stress as follows [28]:(40)σaeq=σmaxσa.Combining equations (39) and (40), the fatigue life under any type i working condition is written as(41)Ni=12σaeqσf′1/b.Letm = −1/b, C = 1/2(σf′)−1/b, and si = σaeq. That is, equation (41) is written in the form of equation (39), where the strength coefficient σf′ of 316 stainless steel is 2100 MPa and the fatigue strength index b is −0.1 [29], and then substituting equation (41) into equation (37) leads to(42)∑i=1kniCsim=1.AssumingNt is the total number of cycles, the number of cycles in the s to s + ds interval is(43)dni=NtPSads.By combining equation (42) with equation (43), the total cycles of material failure are written as(44)Nt=1∫0∞PSasm/Cds.For the narrowband process, from equation (19), the peak stress probability density function is(45)Ps=s2σ2e−s2/2σ2,s≥0,0,s<0.Substituting equation (45) into equation (44), the gross number of stress cycles is obtained as(46)Nt=12mCΓ2+m2σ−m,where Γ(x) is the gamma function.In practical engineering, most of them are broadband random vibration. For broadband stochastic processes, Wirsching and Light introduce a correction factorλ as [30](47)λ=1−am1−εbm+am.The values ofa and b in equation (47) are related to the material parameter m, which can be expressed as(48)am=0.926−0.033m,bm=1.587m−2.323.Thus, a life estimation formula of the broadband stochastic process is obtained as(49)Nt=1λ12mCΓ2+m2σ−m.According to the value of the bandwidth coefficientε, the random vibration fatigue life of the structure can be determined by using equation (46) or (49).In order to verify the correctness of the random vibration fatigue life model, we use the model to compare with the experimental data [31], as shown in Figure 9, and it is found that most of the data are in the double error band and very few data are within the quadruple error band. The life prediction errors are acceptable in engineering, and thus, the applicability of the model is verified.Figure 9
Comparison between the fatigue life model and the experimental data [31].


## 4. Results and Discussion

In this section, the results of the fretting wear model and vibration fatigue model are analyzed, taking CEFR design parameters as an example. The physical parameters of the CEFR areρUO2 = 10 g/cm3, ρNa = 0.8443 g/cm3 at 500°C, the equivalent mass m of the fuel rod per unit length of 0.2934 kg, and the average flow velocity of 4.51 m/s, and other parameters are shown in Table 1. When the parameters are dimensionless, the corresponding dimensionless parameters in the matrix of equation (2) are as follows: β = 0.0752, ε = 225, γ = −0.7505, cf = 0.016, u0 = 0.2984, α = c = cb = Π = 0, and Gc = Gw = 0.17. If the fuel rod is regarded as the Euler–Bernoulli beam model, the natural frequencies of the vibration of the pinned-pinned beam in vacuum are as follows: fv = π/2/L2(EI/m)1/2 = 4.37 Hz.According to equation (3), PSD of rod No. 1 in the in-plane direction and out-plane direction can be obtained, as shown in Figure 10.Figure 10
PSD of the displacement at the fuel rod midpoint.Figure10 shows that the first-order frequency is 3.49 Hz, which is about 20% less than that in vacuum because of the effect of the added mass of the coolant. The corresponding frequency points of peak values in both directions are basically the same. The root mean square value of displacement is about 70 μm. The value is used as amplitude drrms. As shown in Figure 7, Srms = (R2 − (R − drrms)2)0.5. From equations (5)–(28), the calculation parameters of contact wear can be obtained, as shown in Table 3.Table 3
Contact wear zone parameters.
Long half-axis (mm)
Short half-axis (mm)
σ
y
max (MPa)
Indentation depth (mm)
Wear life (days)
0.1063
0.001
56.22
0.337∗ 10−6
1517The equivalent spring stiffnessKstiff has a great influence on the wear and also varies with the operation of the reactor, as shown in Figure 11, and the maximum stress varies with the depth under different equivalent spring stiffness. The maximum depth is about 50 μm, and the stress is 0. According to equation (28), the relationship between the wear life and the equivalent spring stiffness is calculated. As shown in Figure 12, it is found that, with the increase of equivalent spring stiffness, wear life decreases, maximum stress increases, and the effect on wear life becomes smaller.Figure 11
Stress variation along the depth direction under different equivalent spring stiffness.Figure 12
Wear life (left) and the maximum contact pressure (right) with the equivalent spring stiffness change.Table4 shows that ε ≤ 0.35, and then the stress spectrum corresponding to the pressure spectrum, i.e., equation (4), is a narrowband stochastic process. By substituting the equivalent stress amplitude of 214 MPa into equation (46), it is calculated that, after 81 h, the fatigue crack initiation of the cladding at the contact position begins. As shown in Figure 13, the S-N curve calculated from the random load model in this paper is compared with the S-N curve calculated from the SWT model. It is found that the random load has a more severe effect on fatigue life. By comparing the wear life with the fatigue life, it can be judged that fatigue failure is the main cause of fuel rod failure. But there are no corresponding theoretical and experimental results about the WAFR contact.Figure 13
S-N curve.Table 4
Computational parameters of random vibration fatigue.
e
f
ave (Hz)
σ
rms (MPa)
Fatigue life (h)
0.0196
3.68
214
81


## 5. Conclusion

Considering the dynamic behavior of rod bundles under turbulent axial flow and the environment in which the fuel rods are mainly subjected to high temperature and strong irradiation effect, it is found that the wire wrap and adjacent fuel rods will make contact by estimation. A theoretical model for life analysis of the WAFR contact is established. Taking the parameters of the CEFR as an example, the conclusion can be summarized as follows.Based on the study of Païdoussis’ team, the dynamic equation of flow-induced vibration of rod bundles is established and the displacement power spectrum is obtained. The displacement amplitude is about 70μm; because of the effect of the additional mass of the fluid, the natural frequency in the coolant is about 20% less than that in vacuum.If the contact equivalent stiffness is 1.42 N/mm, the wear life of the cladding is 1517 days, and it is found that the effect of the contact equivalent spring stiffness on the wear life decreases with the increase of the contact equivalent spring stiffness.The random load irregularity factorα = 0.98 is calculated; that is, the bandwidth coefficient ε = 0.0196. Therefore, the random process is treated as a narrowband process, and the random vibration fatigue life of the contact position is calculated to be 81 h according to the random load model.The results show that fretting wear and vibration fatigue are important problems in WAFR contact life analysis. The work preliminarily established the mechanical model of cladding failure about the WAFR contact for the first time. In the follow-up, the elastic-plastic problems in the WAFR contact will be further analyzed so as to establish a more accurate model to predict the life of fuel rods.



---

*Source: 5613737-2019-09-17.xml*