Modelling of the friction-vibration interactions in reactor core barrel couplings

The key-groove couplings between the lower part of core barrel and reactor pressure vessel in the nuclear VVER-type reactor show assembly side clearances. Vibration of reactor components, caused by coolant pressure pulsations generated by main circulation pumps, produces impulse and friction contact forces in the above mentioned couplings. These forces are used for calculation of power and work of friction forces in the slipping contact surfaces between the key and the groove of the particular couplings and for their fretting wear prediction in dependence on the clearance values. The paper presents an original approach to mathematical modelling and simulation analysis of reactor nonlinear vibrations respecting friction-vibration interactions at all the key-groove couplings with clearances. The method is applied to VVER 1000/320 type reactor operated in the nuclear power plant Temelı´n in the Czech Republic.


Introduction
The existence of clearances in the key-groove couplings between the lower part of core barrel and the reactor pressure vessel in the nuclear VVER-type reactors [12,13] is necessary and unavoidable, because of manufacturing and assembly tolerances of the reactor components.Clearances in couplings in a physical system generally lead to several complex dynamic phenomena such as surface contact, shock transmission and the development of different regimes of friction, the slider and free flight motion of connected parts and wear.A comprehensive survey is presented in [16].
Vibration of VVER-type nuclear reactors excited by seismic action was investigated in previous author's research works and papers [5,6,19] in co-operation with the Nuclear Research Institute R ˇez ˇin Czech Republic on linear clearance-free models with proportional damping.The original linearised spatial model of the VVER1000/320 type reactor (Fig. 1) was used for the calculation of the spatial motion of reactor components excited by pressure pulsations generated by main circulation pumps [21,22] and it was partially experimentally verified in [11].
The decomposition method [20] was used for development of the global reactor model with all couplings between its components.The reactor was decomposed into eight subsystems displayed in Fig. 1.The abbreviations of this subsystems, the list of generalized coordinates used in the reactor mathematical model and their sequence in the global vector of generalized

B
damping matrix of the reactor c i slip velocities in the central contact points of K-G couplings f (c i ) friction coefficient in K-G coupling i dependent on slip velocity F (k) P V amplitude of k-th harmonic component of the excitation force acting on PV f j vector of geometrical parameters corresponding to the excitation force generated by pump j (j = 1, . . ., 4) f (q, q) vector of nonlinear forces in all K-G couplings k stiffness of one key with cantilever in normal direction K stiffness matrix of the linearised reactor model with clearance-free and smooth K-G couplings M mass matrix of the reactor N i normal contact forces P i friction power in K-G coupling i q(t) vector of generalized coordinates dependent on time s i tangential shift of the groove with respect to key from central position of K-G coupling i t time T i,r radial components of friction forces T i,ax axial components of friction forces u i relative tangential displacement of the groove with regard to key of K-G coupling i W i work of friction forces in K-G coupling i δ j phase angle of the pump j Δ i clearances in K-G couplings Δm (H) i hourly fretting wear in K-G coupling i μ X (ω) experimentally obtained fretting wear parameter ω j angular frequency corresponding to the revolution of pump j  1.The components marked by grey were considered as rigid bodies with six degrees of freedom -three displacements in directions of axes collinear with global coordinate system x, y, z and by three angular rotation about these axes.A position of the particular components is described by relative displacements with respect to supporting bodies.The supporting body for CB, BPT, UB, DH and DA is reactor pressure vessel PV.Lower part CB3 of core barrel is supporting body for RC with fuel assemblies FA (see Fig. 1) and DH is supporting body for EM.Many components of the beam type (FAs, protection tubes PT1 and PT2 of BPT, tubes T and circular rods R of UB and components of DH and DA of linear stepping drives) are modelled as one-dimensional continua linked with rigid bodies by discrete couplings or ideal boundary conditions in longitudinal, torsional and transverse directions.Vibration of all uniform continua (e.g., 163 fuel assemblies) were supposed to be identical.The FA global Fig. 1.Scheme of the reactor assumed from [22] models were considered as clamped-clamped beams with concentrated masses in the chosen nodes on the FA axis on the basis of experimentally gained eigenfrequencies and eigenvectors coordinates corresponding to FA lateral displacements in the nodes [5].An influence of the coolant surrounding FAs was respected by additional masses resulting from decreasing of FA eigenfrequencies and conservation of FA eigenvectors.
Reactor mass matrix M was determined from the quadratic form of kinetic energy E K ( q) = 1 2 qT M q and reactor stiffness matrix K from the quadratic form of potential energy E P (q) = 1 2 q T Kq, where q is the vector of generalized coordinates (see Table 1).Substituting the kinetic and potential energy to Lagrange equations we obtain the conservative linear model of the reactor.The linear model of the reactor VVER1000 was used for the modal analysis and confronted with experimental frequency analysis continuously provided in Nuclear Institute R ˇez ˇin terms of project "Fuel cycle of NPP" [4,7,17].A detailed description of the reactor subsystems with couplings is presented, e.g., in paper [6], monograph [7] and full text on CD-ROM in [20].
The linear mathematical model of the reactor after completion of the damping approximated by a proportional damping matrix and an excitation by coolant pressure pulsations generated by main circulation pumps has the form where M , B, K are mass, damping and stiffness matrices of order 137 (see Table 1).Proportional damping matrix B satisfies the general condition in the form where V is the modal matrix of the reactor conservative linear model and Ω ν are its eigenfrequencies.The dimensionless damping factor D ν corresponding to particular mode shapes were estimated equal to 0.05 with the exception of the larger damping factors corresponding to the mode shapes, where the reactor core vibrates dominantly (an influence of coolant).
There are analytical methods to estimate the pressure pulsations [2] that are observed in the real reactors [8].The modelling of the reactor VVER1000 excitation, which will be further used to estimate vibration caused by pressure pulsations of the coolant, is based on the theory presented in work [10].The hydrodynamic forces CB cos(kω j t + δj), j = 1, . . ., 4 (2) in the space between the reactor PV and the CB walls (Fig. 1) generated by each main circulation pump j act on PV and CB in the middle between their bottom and PV-CB sealing in directions of the primary coolant loops (see Fig. 2).These forces must be transformed into the configuration space of reactor generalized coordinates [22].The excitation vector on the right side in (1) is expressed by means of amplitudes F (k) P V of k-th harmonic components of the hydrodynamic forces.Vectors f j of dimension 137 corresponding to particular pumps j (j = 1, 2, 3, 4), derived by Lagrange equations in [22], depend on ratio of exterior core barrel surface and interior reactor pressure vessel surface and on generalized coordinates structure (see Table 1).Each hydrodynamic force has the polyharmonic character having one another slightly different pump angular speeds ω j corresponding to the particular pumps and generally different phases δ j .As a result, the reactor vibration has a beat character.The fluctuation of the pump angular speeds is based on the measurement at NPP Temelı ´n blocks.Application of these forces in original clearance-free linear mathematical model (1) enables approximate calculation of friction forces in key-groove couplings (Fig. 3) during slip motion without the consideration of possible free flight motion followed by clearances.The aim of this paper is the investigation of reactor nonlinear vibration respecting the assembling side clearances and friction in eight key-groove couplings (Fig. 3) uniformly deployed to circumference between the lower part of CB (CB3) and PV in the nuclear VVER1000/320 type reactor operated in nuclear power plant (NPP) Temelı ´n in the Czech Republic.Other important phenomena associated with clearances in above-mentioned couplings, such as impact and friction forces, power and work of friction forces and fretting wear, are calculated and discussed in dependence on clearance values.

A mathematical model of the reactor with clearances in the core barrel couplings
Let us consider an arbitrary chosen key-grove (K-G) coupling i defined by contact plane angle α i (Fig. 4) from axis x parallel to basic reactor co-ordinate axis x (Fig. 1).The contact area is small with respect to distance r from central reactor axis y.Consequently, we can concentrate contact between the concrete key on reactor PV and the groove on CB3 in one central point.Clearance Δ i can be defined as the distance between the groove and key contact surfaces in a central position (Fig. 4).The nonlinear normal force characteristic (see Fig. 4b) was used to model the gap effect.
Relative normal displacement u i of the groove surface contact on CB3 from the starting position compared to the key can be written as where row vector d T i corresponding to coupling i in reduced form (only with non-zero generalized coordinates at the above marked positions according to vector q defined in Table 1 and depicted in Fig. 4) is Angle α i describes a position of K-G coupling, b 0 is the perpendicular distance of coupling plane from core barrel suspension on pressure vessel flange (level of point B), h 0 is the perpendicular distance of couplings plane below gravity centre CB3 level and r is the central contact radius of all couplings.Normal contact force N i transmitted by coupling i for general starting position of the CB3 on the level of K-G coupling plane relative to the reactor PV, defined by shift s in direction determined by angle α (see Fig. 4), can be expressed in the form where k is the stiffness in normal direction of one key connected with reactor PV by means of the cantilever (see Fig. 4) and Δ i is the clearance in the corresponding coupling i.Small shift s i of the groove with respect to key in direction of displacement u i from the ideal central position is Heaviside functions H in (5) are zero when the key contact with groove in coupling i is interrupted (−u i < Δ i + s i or u i < Δ i − s i ) and relative motion of connected part in normal direction has the free flight form.When the contact between the key and the groove takes place, normal impact force N i and friction force f N i with radial T i,r and axial T i,ax components (see Fig. 5) appear in the contact where c i,r and c i,ax are its components in radial and axial directions (see Fig. 5).Dependence f (c i ) can be approximated by smooth function [4,14] in the general form The friction characteristic is specified by means of four shaping parameters, where f 0 is the static friction coefficient, f d is the dynamic friction coefficient, ε 1 and d is the exponent in exponential function.A graphical map for applied shaping parameters f 0 = 1.5, f d = 1 [15], d = 100 and various values ε = 5×10 3 , 10 4 and 2×10 4 is presented in Fig. 6.In case of f 0 = f d , friction characteristic (8) transforms into the standard Coulomb friction model approximated by smooth function.
The radial and axial components of the relative velocities in particular K-G couplings can be written by means of row vectors of geometrical parameters corresponding to coupling i in the form The reduced form (only with non-zero coordinates at the marked positions according to vector of generalized speeds q) of vectors d T i,r and d T i,ax is The components of friction forces have opposite directions with respect to corresponding components of the slip velocities (see Fig. 5) and can be expressed in the form All contact forces transmitted by K-G couplings between the reactor PV and the CB including normal N i and friction forces T i,r , T i,ax for i = 1, . . ., 8 generate in the reactor nonlinear model force vector f (q, q) of dimension n = 137.This contact force vector can be derived using the virtual work principle.The virtual shift of the groove relative to key on the reactor PV in normal, radial and axial directions, according to (3) and ( 9), can be expressed as The virtual work in one K-G coupling corresponding to particular general coordinate q i from Table 1 is whereas the virtual shift only of one general coordinate is non-zero.Generalized forces Q i , i = 7, 8, 9, 16, . . ., 21 result from relations (see Fig. 5, where subscript CB was let-out from the CB generalized coordinates) Vector f (q, q) of contact forces in the mathematical model of the reactor includes frictionvibration interactions in all K-G couplings with clearances.Its reduced form is whereas order of non-zero coordinates is depicted on the right.The friction-vibration interactions in K-G couplings have an effect on linearised mathematical model (1) which is rearranged into the reactor nonlinear model The vector of elastic forces K C q(t) of all clearance-free and smooth couplings (for 1) is replaced by nonlinear force vector f (q, q).Softened global stiffness matrix K − K C corresponds to the linear reactor model without K-G couplings.

Modal analysis of the reactor conservative model with clearances in K-G couplings
To illustrate the modal properties of the reactor, chosen eigenfrequencies of the linear clearancefree model are presented in Table 2.An influence of the K-G couplings loosening was studied by harmonic balance method [18].In the first approximation, harmonic vibration with frequency Ω ν of the conservative nonlinear mathematical model of the reactor was supposed.The relative tangential displacements of the K with respect to G from central position can be expressed as u i,ν = a i,ν cos Ω ν t, i = 1, . . ., 8, ν = 1, . . ., 137.According to (3), the amplitude in K-G coupling i corresponding to mode shape ν can be expressed in the form a i,ν = d T i v ν , where v ν is the eigenvector in the configuration space of reactor generalized coordinates.Based on the harmonic balance method, equivalent stiffness of coupling i with clearance Δ i and amplitude a i,ν can be expressed for |a i,ν | ≥ Δ i as where k = 1.22 × 10 9 [N/m] is the stiffness of one key with cantilever (see Fig. 3) in contact phases.For |a i,ν | < Δ i is k i = 0.The reactor equivalent stiffness matrix for particular mode shape ν respecting concrete amplitudes a i,ν in couplings is where K i is the stiffness matrix corresponding to single K-G coupling i in linear reactor model (1) with unit stiffness.Eigenfrequency Ω ν and eigenvector v ν of the linearised conservative reactor model are calculated separately for chosen mode shape from the eigenvalue problem in several iteration steps in dependence on lossening parameter p = max i a i,ν /Δ i .The iterative process for estimation of amplitudes a i,ν starts from the clearance-free linear model of the reactor and is finished when relative error of the calculated eigenfrequency in following iteration steps is small.
As an illustration, eigenfrequencies Ω 34 to Ω 37 are plotted in Fig. 7 in dependence on loosening parameter p ∈ 2, 10 .Eigenfrequencies Ω 20 , Ω 34 and Ω 36 (see Table 2) in calculated frequency range up to 40 Hz are the most affected by clearances in K-G couplings.

Dynamic response of the reactor
The dynamic response of the reactor excited by coolant pressure pulsations generated by main circulation pumps is investigated by integration of nonlinear motion equations ( 14) in time domain.The vector of excitation and nonlinear contact forces on the right-hand side in model ( 14) can be formally rewritten into the global nonlinear vector Motion equations ( 14) are transformed into 2n nonlinear differential equations of the first order where Eqs. (19) are integrated using the MATLAB built in ode45 solver with adaptive step length under zero initial conditions u = 0.All the simulations are performed in large time interval 0; T , where T = 100 [s] is the large period of beat reactor vibration specified by the minimum difference in pumps angular speeds.

A prediction of the friction characteristic and the fretting wear of K-G contact surfaces
The vibration induced by pressure pulsations in the reactor produces relative motion between the key and groove in all couplings between the lower part of CB and reactor PV.Friction forces in the slipping contact surfaces may produce their fretting wear.The most of published papers about fretting wear deal with experimental investigation of the friction work per unit area which is converted into wear depth during a certain time interval [1] by means of the wear coefficient.Calculation of the friction work based on estimation of the average normal contact force [9] does not respect the time dependence of the normal contact forces and possible temporary loss of fretting contact in phases of possible free flight motion followed by clearances.Due to fretting wear of the contact surfaces, there is the potential for increasing of the clearances.
For an assessment of this phenomenon, friction coefficient f and fretting wear parameter μ [g/J] are obtained experimentally [13,15].The fretting wear parameter is defined as loss of mass in one contact surface generated by the work of friction force W =1 [J] at excited frequency ω [rad/s] of the pumps revolution [3].The main components of the experimental stand (see Fig. 8) create two hydraulic cylindersvertical (pressure) and horizontal (draft) -and two plates fixed lower representing groove and moving upper representing key.The experimental plates are produced from materials corresponding to K-G couplings in the VVER type reactor.Harmonic vibrations of the upper plate take place in the water reservoir [15] at measured the water temperature.The friction coefficient, after correction on inertia force and friction in joints of pressure cranks, is determined from the measurement of tractional force generated by the horizontal cylinder in dependence on slip velocity c.Harmonic motion frequency 16 Hz corresponds to main frequency of pressure pulsations.
On the basis of the experimental determination, the friction characteristic in the slipping contact surfaces in K-G couplings is modelled in the form (8), where ε, f 0 , f d , d are shaping parameters.The measured values of friction coefficients f 0 and f d have dominant influence.Fretting wear parameter μ has been calculated using relation μ = Δm/W , where Δm [g] is the measured loss of the mass [g] on one plate due to fretting wear caused by work of the friction force W = 1 [J].The methodology of μ parameter determination was published in [13,15].The upper plate is moved harmonically with frequency 16 [Hz] in the time that correspond to up to 5 [km] of plate's sliding.The mass loss of both plates in contact areas are precisely measured.
The criterion of the fretting wear of contact surfaces in coupling i can be expressed using the work of friction forces [4,17] during representative time interval t 1 , t 2 as where indicates the power of friction force in K-G coupling i.The fretting wear in grams in both opposite surfaces of the key (K) and the groove (G) can be expressed as where μ X (ω) is the experimentally obtained fretting wear parameter of the corresponding moving part (K or G) at harmonic excitation frequency ω.Time interval t 1 , t 2 for calculation of the friction work corresponds to beat period.The hourly fretting wear (the wear during one hour) of the coupling components is Fretting wear Δm i = Δm i,K + Δm i,G causes clearance increasing in K-G coupling i.

Application
The presented method is applied to the pressurized water reactor of the Russian construction marked with abbreviation VVER1000.Its nominal electric power is 1 000 [MW] and the thermal power is 3 000 [MW].The operational pressure of the coolant inside of pressure vessel is 15.7 [MPa] and maximal temperature of the coolant is 320 • C. The basic information about reactor proportions provides total mass about 9 ×10 5 [kg] with coolant inside of reactor pressure vessel.All applied geometrical, mass and stiffness parameters of reactor components correspond to design parameters of the VVER1000/320 type reactor [7,19].The main quantities which influence impact loading of the CB and the reactor PV and the fretting wear in K-G couplings are normal forces N i (5) transmitted by K-G couplings, slip velocities c i (7) and works W i (21) of the friction forces.Their values are calculated for coolant pressure pulsations [22] defined in Table 3 and zero phases δ j in (1).Harmonic components Basic revolution frequencies [Hz] of force Basic excitation frequencies ω j = 2πf j , j = 1, 2, 3, 4 correspond to the rotational speed of main circulation pumps.The friction characteristic (8), corresponding to different measured values of friction coefficients f 0 ≈ 1.5 ÷ 2.5 and f d ≈ 1 ÷ 1.5 [15] are approximated using shaping-parameters ε = 10    11 respect a sense of the normal force owing to information about of the contact side (N i > 0 for u i > Δ i − s i and N i < 0 for u i < −(Δ i + s i )).A dissipation is respected in (21) by means of the absolute value |N i |.Time interval 100 [s] of the vibration simulations includes the long period of the beating vibration caused by slightly different main circulation pumps revolutions [22].
The time behaviours of the state quantities in the maximal loaded coupling i = 3 for extreme friction determined by friction coefficients f 0 = 1.8, f d = 1.2 and maximal admissible design clearances Δ i = 75 [μm] in all eight couplings are shown in Fig. 12.The behaviours of normal contact force N 3 and relative velocity c 3 in the short time interval around maximal value of force N 3 are illustrated in Fig. 13.The areas with relative velocity c 3 approaching to zero values are given by low sliding velocities in phases of contact between key and groove corresponding to large impact force N 3 .On the contrary, the areas with normal contact force N 3 limiting to zero values correspond to free fligt motion of the key within clearance.It's evident that the normal force behaviour is almost periodical with period 0.06 [s] corresponding to main excitation frequency 16.64 [Hz], whereas the velocity behaviour is rather chaotic.The influence of the friction is illustrated in Table 4 for ratio of friction coefficients f 0 /f d = 1.5 on the basis of maximal values of the state quantities, friction work and fretting wear during large period T = 100 [s] of beat vibration.The coupling identification with maximal state quantity is depicted by means of index i.The design value of total clearances 2Δ i between the key and the groove range from 0.05 to 0.17 [mm].Due to supposed fretting wear, there is the potential for increasing the clearances.An influence of the clearance magnitude is illustrated in Table 5 for Coulomb friction (f 0 = f d ) in K-G couplings.Based on the results, increase of the friction has a negative influence on impact loading of the reactor pressure vessel and the core barrel and the fretting wear of contact surfaces in couplings.From the point of view of clearances, an optimal situation corresponds to minimal values of the total clearances 2Δ ≈ 50 [μm].The clearances increase as a result of fretting wear.The fretting wear can be positively influenced by the reduction of friction and an initial assembling gap between the K-G couplings.The noncentral position of the CB with regard to the PV has no practical effect on a steady-state vibration of the CB, impact forces and the fretting wear.

Conclusions
The main objective of the paper is to present the new basic method for estimation of impact and friction forces transmitted by K-G couplings between the lower part of CB and the reactor PV.These couplings centre the CB position carrying reactor core.A presence of the larger clearances leads to undesirable nonlinear dynamic phenomena such as impact loading, increasing of reactor vibration, FAs deformations and fretting wear of the contact surfaces in mentioned couplings.
The proposed method for the fretting wear prediction in K-G couplings is based on the numerical simulations of the reactor nonlinear vibration excited by the pressure pulsations of coolant and experimentally obtained fretting wear and friction parameters.The global nonlinear model of the reactor was developed from the original linear model whereof all clearance-free and smooth K-G couplings were replaced by nonlinear couplings with clearances and friction.All numerical simulations have been applied to pressurized water VVER 1000/320 type reactor with the friction and the fretting wear parameters evaluated experimental in Research Centre R ˇez ˇin Czech Republic.The influence of fretting wear parameters μ K and μ G on fretting wear values of the contact surfaces in K-G couplings is linear.The results, for more precise values of these parameters investigated experimentally in the conditions that are closer to real state in the reactor, could be modified linearly without new time demanding numerical simulations.The sensibility analyses focused on variation in the friction characteristic and clearances lead on as possible to reduce the contact forces transmitted by K-G couplings and their fretting wear.The method can be adapted for other mechanical systems with similar couplings and friction-vibration interactions.

Fig. 3 .
Fig. 3. Key-groove coupling between core barrel and reactor pressure vessel

Fig. 4 .
Fig. 4. Key-groove layout in contact plane between lower part of core barrel (CB3) and reactor pressure vessel (PV) (a) and stiffness characteristic of one coupling (b)

Fig. 5 .
Fig. 5. Contact forces transmitted by key (K)-groove (G) coupling i for u i > Δ i − s i

Fig. 8 .
Fig. 8. Experimental stand for assessment of the friction coefficient and the fretting wear parameter 4 and d = 100 [s/m].A prediction of the fretting wear of the K-G contact surfaces is calculated according to (21) and (23) for values μ K = 10 −9 and μ G = 5 × 10 −9 [g/J].

Fig. 12 .
Fig. 12.Time behaviour of normal contact force N 3 , relative velocity c 3 and friction power P 3 in the coupling 3

Fig. 13 .
Fig.13.Normal contact force N 3 and relative velocity c 3 in the short time interval (segment of Fig.12)

Table 2 .
Chosen eigenfrequencies of the clearance-free reactor model

Table 3 .
Parameter of pressure pulsations

Table 4 .
Maximal values of the state quantities-friction action (f 0

Table 5 .
Maximal values of the state quantities-clearance influence (f 0