Experimental and Numerical Investigation of Contact Parameters in a Dovetail Type of Blade Root Joints

: this paper focuses on the contact characteristics of the blade root joints subjected to the dry friction damping under periodic excitation. The numerical method and experimental procedure are combined to trace the contact behavior in the nonlinear vibration conditions. In experimental procedure, a novel excitation method alongside the accurate measurements is used to determine the frequencies of the blade under different axial loads. In numerical simulations, local behavior of contact areas is investigated using the reduction method as a reliable and fast solver. Subsequently, by using both experimental measurements and numerical outcomes in a developed code, the global stiffness matrix is calculated. This leads to find the normal and tangential stiffness in the contact areas of a dovetail blade root joints. The results indicate that the proposed method can provide an accurate quantitative assessment for investigation the dynamic response of the joints with focusing the contact areas.


Introduction
Complex mechanical structures usually include several subassemblies, which may have contact surfaces in common. These interfaces affect the behavior of the entire system. It can be realized grossly by introducing two contact parameters, namely contact stiffness and equivalent damping.
The major sources of damping in systems, beneficial for controlling the amplitude of vibrations, are the material damping, aerodynamic damping, and friction at contact interfaces. Regarding the last type, the energy is dissipated via friction when contact surfaces undergo a relative displacement. Micro-sliding or "microslip", macro-sliding or "gross-slip", and separation or "lift-off" are the three common states in the surfaces contact behavior. Partial sliding of the regions at the edge of the contact generates microslip.
Regarding the micro-scale measurements, sliding contacts between spherical surfaces were first analyzed by Mindlin et al. [1], who measured the hysteresis loops for convex surfaces. Johnson [2] measured the static and dynamic hysteresis of steel balls pressed on flat surfaces. Goodman and Brown [3] examined steel balls oscillating between two parallel planes in compression. Filippi et al. measured the friction coefficient, cycles of hysteresis and the tangential contact stiffness at room temperature [4] and at high temperature [5] for spherical surfaces in contact with flat surfaces. The sliding contacts between compliant surfaces were examined by Lavella et al. [6] by using a test bench, and measured the hysteresis cycles, the friction coefficient and the tangential stiffness of contact by varying normal load and excitation frequency. Schwingshackl et al. [7] and Stanbridge et al. [8] obtained the coefficient of friction and tangential stiffness of a plane contact for different materials and load conditions. To compute the Jacobian matrix for 3D friction contact modeling using an alternate frequency time domain method was presented by Afzal et al. [9]. In a separate study [10], they showed the damping potential of multiple friction contacts in a bladed disk, and the possibility of increasing the damping effectiveness. Recently, a few studies [11][12][13] devoted their attention to analyze the nonlinear dynamic behavior of blade-like structures, namely disk-beam, rotating beam, and rotating disk-beam, with dry friction dampers.
Regarding the macro-scale measurements, the first experimental analyses of the sliding contacts effect on global dynamic behavior of aeronautical turbine rotors were carried out by Goodman and Klumpp [14], who measured the hysteresis cycles for discblade joints subjected to a uniform compression load. Allara et al. [15] designed a test bench to study the behavior of dovetail joints between disc and blade as the centrifugal load and the vibrations amplitude vary. Firstly, the oscillations' free decay of an axially loaded beam under electromagnetic excitation was measured. Then, the obtained signals were analyzed to derive the frequency and the lost factor of the nonlinear system, as a function of the traction load and amplitude of vibration. The test bench is used in the present work for conducting experimental investigations. The same test equipment was utilized by Firrone and Bertino [16] to measure the forced response in the frequency domain to validate the numerical results based on the Harmonic Balance Method. Asai et al. [17] analyzed the behavior of a cantilever beam dovetail root joints of different geometries and evaluated the corresponding friction damping in the system. Simmons and Iyengar [18] employed the same method to contiguous palettes to study the tension distribution and notch effect in the blade attachments. Umer and Botto [19] evaluated the global dynamic behavior of a vane in contact with under platform dampers, simultaneously investigated the local contact behavior to realize the performance curves of these dampers. Very recently, Li and coauthors [20][21][22] conducted several in-depth studies, contributing to better understanding of different contact models. For instance, in Ref. [22], Iwan density function was utilized to model sphere-on-sphere contact and flaton-flat contact, and the results were then validated with analytical and experimental data.
The present work carries out an experimental procedure to measure the dynamic responses of a typical dovetail joint used for connecting the blades to the disc. Extracting the contact parameters of a blade simulacrum as a function of centrifugal force leads to deepen understanding of joints' dynamic performance. Furthermore, it provides a valuable measured database during free decay of oscillations to validate newly proposed methods of nonlinear temporal integration. Eventually, the current contribution is an attempt to obtain normal and tangential contact stiffness from the macroscopic measurements to find needed parameters for contact models development.

Experimental Configuration
In this section, the response due to an excitation is analyzed in terms of centrifugal force and the amplitude of vibration. The dummy blade beam, which replaces the real blade is shown in Figure 1.

Test Bench
Experimental equipment utilized in the current study was in use as described in Ref. [15], and it consists of the following five subsystems:

1.
The support and traction system of the blade simulacrum, shown in Figure 2. This system includes two crosspieces: a fixed piece integrated with a support of dovetail type (c); and a mobile piece connected to a hydraulic actuator that enables the traction (b). 2. The dynamic excitation system, which includes a signal generator, an amplifier, and a shaker. To obtain the response associated with the mode shapes, a type of input signal, shown in Figure 3, is applied. Performing excitation by the shaker was utilized for the first time in the present study unlike other contributions [6,15,16]. 3. The vibration measurement system, which consists of a laser pointer, Figure 4, and a laser controller. As the specimen were excited at frequencies close to the first and second flexural modes (1B and 2B), the laser was pointed towards the nodes of the corresponding mode shapes, where the maximum amplitude occurs.

4.
The traction force measurement system, incorporating two independent subsystems that allow for double checking the measurements. The first consists of a pressure gauge, providing pressure inside the hydraulic actuator. The second is composed of a Wheatstone bridge strain gauges ( Figure 5) for measuring the axial deformation δ of the beam and a data acquisition system.

5.
The data collection and acquisition system, which collects the velocity signals measured by the laser and the strain from the gauges and the input currents of the shaker, later to be used for further analysis.

Technical Limitations
After obtaining the velocity signal for each applied tensile force, it is possible to realize the dependence of nonlinear modal parameters, i.e., natural frequency and damping from the vibration amplitude. To this purpose, the free decay of the envelope velocity signal from the excitation starting points until the detachment of the shaker tip should be analyzed. However, the signal is still the superposition of all modes of the beam vibration, as shown Figure 6. A filter has therefore been implemented to extract the main modal component without altering the phase of the signal, necessary for subsequent identification of the modal parameters. A Parks-McClellan Finite Impulse Response (FIR) filter with a ripple of 1•10-40dB in the passband and a ripple of −40dB in the remaining frequency band was adopted (see Figure 7).  To identify the nonlinear dependence of the natural frequency and the damping η on the amplitude of vibration, the theory of analytical signals and Hilbert transform were utilized following the "FREEVIB" method proposed by Feldman in [23] for nonlinear systems. It can be only applied to systems that have frequency-independent damping [15,20]. To obtain a stable solution, particular methods such as that proposed in Ref. [20] should be taken into account due to Kelvin type structural damping in the model.
The measurement points, (P1) for the first and (P2) for the second flexural mode is displayed in Figure 8. Concerning the first mode, the dependence of damping on the vibration amplitude for different value of axial loads is shown in Figure 9. For small vibration amplitudes, it observes a constant damping, i.e., independent of amplitude. This suggests the complete adhesion at the contact region, implying the linearity of the system in question. For amplitudes greater than a critical value, the damping shows a steep increase, indicating the beginning of possible sliding in the contact region, and thereby the establishment of the hysteresis loop due to the alternation between adhesion and sliding states during each oscillation cycle. Furthermore, the dissipation peaks at a critical value of vibration amplitude, as can be seen in Figure 9, and eventually decreases as amplitude exceeds that critical value. Moreover, the peak (maximum damping) drops as the traction force F increases. This is because not all four contact regions are in the state of macro-creep, and some regions are in the state of microslip, i.e., sliding only at edges of the contact regions. Indeed, various combinations of the macroslip and microslip states at the four contact zones could influence the results in different manners. This hypothesis is compatible with the previous supposition on the variation in force normal to the contact region. The variation of the frequency against the vibration amplitude is depicted in Figure  10. Based on the figure, rather flat curves for small oscillations suggest linearity of the system, the indication of complete adhesion. As the amplitude grows, the frequency decreases sharply, signaling the occurrence of possible sliding in the contact region. Eventually, the stiffening phenomenon (increase in frequency) is observed for further larger values of vibration amplitude, which can be attributed to the variation of normal force at the contact region caused by the bending of the blade simulacrum.  Regarding the second flexural mode, variations of damping and frequency with respect to the vibration amplitude are shown in Figure 11. According to this, an initial stretch of constant damping is observed for small vibration amplitudes. Then, a steep growth is evident beyond a critical amplitude at which the sliding begins. Considering the same traction force, the critical amplitude depends upon the effective I_RMS current sent to the shaker. This current defines the amplitude of oscillation of the exciting force before the detachment of the shaker tip from the beam, and therefore the amplitude of the vibration at the beginning of the free decay. Close examination of the figure suggests that the highest damping is achieved when a traction force of approximately 8 kN is exerted to the beam, and more intense loading results in a drop in damping. Besides, a decrease in the critical amplitude is observed as the traction force increases, the phenomenon also presents in the first flexural mode.
The dependence of the frequency on the vibration amplitude is also illustrated in Figure 11, where the frequency follows a constant trait for low amplitude values, revealing the complete adhesion at the contact regions. Stiffening effect, already identified for the first bending mode, can be recognized from the frequency trend against the traction force.

Contact Model
In this part the normal and tangential stiffness components, i.e., and , are investigated by considering the equivalent model of the dovetail joint (see Figure 12). The following procedure was followed to search for the values of and : 1. The model of the beam with dovetail supports at both ends as well as the corresponding supports was discretized with the finite element method implemented in commercial finite element software, paying particular attention to the coincidence of the nodes in the contact interfaces. 2. Nonlinear static analyses were carried out in commercial finite element software on the beam alone as the axial traction varied to discern the stiffening effect of the force on the stiffness matrices of the blade simulacrum. 3. Reducing the model DOF in commercial finite element software using Craig-Bampton's Component Mode Synthesis technique [24] and importing the reduced mass and stiffness matrices into a numeric computing environment. 4. The reduced models of the blade simulacrum and supports (slots) were assembled in a numeric computing environment by introducing linear contact elements.
A script was developed to search for and , which compares the frequencies of the first and second flexural modes from the FEM model with those obtained from the measurements within the range of axial traction forces in which complete adhesion of the contacts was recognized.

Finite Element Method
To discretize the model using the ANSYS ® 2020-R2, SOLID185 8-node hexahedral linear elements (ANSYS, Canonsburg, PA, USA) were utilized for meshing, the size of which was gradually reduced in the vicinity of the contact areas. To employ Node-to-Node contact elements in the subsequent assembly in a numeric computing environment, the nodes at the contact interfaces were set to coincide geometrically, as shown in Figure 13.

Reduction Technique
The search for the values of the and was carried out by comparing the frequencies of the first and second flexural modes of the FEM model with those obtained from the measurements. This procedure required the use of a global method varying and . In this regard, the degrees of freedom of the beam and support were reduced. In particular, the technique proposed by Craig-Bampton [24] was adopted in the present study, already implemented in commercial finite element software.
The method implements a change of coordinates of the equations of motion to obtain a description of the system with a reduced number of generalized coordinates. The first step is to break down the degrees of freedom of the model into master DOFs (m), which will be also present in the final reduced model, and slave DOFs (s).
In a general case, the system of equations of the motion is of the form: where the matrices [M] and [K] and the vector {u} can be considered in the following form: (2) 0 0 Defining the static modes as the displacements of the slave DOFs due to the displacements of master DOFs, one can extract: with (4) The internal modes are defined as: (5) When all the "master DOFs" are bound, it is now assumed that the movements of the "slave DOFs" can be approximated by a combination of a subset of the internal modes as: Defining a new generalized coordinate system consisting of the master DOFs and the coordinates as: Thus, the movements of the slave DOFs can be written according to the new coordinates using the static modes and internal modes in the following form: The above relation implies that the initial degrees of freedom and the new generalized coordinate system are related by the following transformation matrix: 0 ̅ Defining 0 and ̅ , the equations of motion of the structure in the new coordinates are therefore: where (11) To analyze the problem of interest, the master nodes chosen for the beam are:  33 nodes for each of the four contact interfaces  21 nodes along the longitudinal axis of the beam on each side, see Figure 14b, to be capable of viewing the mode shapes of the reduced model.
Regarding each slot, the following master nodes were chosen:  33 nodes for each of the two contact interfaces  8 nodes corresponding to the vertices of the support for displaying the modes.  117 nodes at the base of the support, see Figure 14a, to be subsequently constrained in a numeric computing environment.
Once the master nodes for the blade simulacrum and the supports were defined, the reduction of each was carried out and the respective reduced matrices and were exported to a numeric computing environment.

Assembly of Stiffness and Mass Matrices
The derived local matrices and from the reduced models of the beam and two supports are used to calculate the global matrices and of the assembly. As described in Figure 15, the degrees of freedom of the reduced models were ordered for both beam and slots. By dividing the degrees of freedom into two categories, namely free (F) and the constrained (C), it was possible to split the mass and stiffness matrices as: The modal analysis was then carried out using the matrices and , obtained from the global matrices and by eliminating the rows and columns corresponding to the degrees of freedom of the bases of the two slots. The structure of the resulting matrices is exhibited in Figure 16. To insert the contact elements to the interfaces, it is necessary to define a connectivity matrix to associate the degrees of freedom of the coinciding nodes of the blade simulacrum and supports. Each row of the connectivity matrix corresponds to a pair of coinciding nodes between which a contact element operates. Coincident nodes were identified by looking for nodes with the same identification number used in commercial finite element software, according to the procedure explained above. Linear contact elements were created, i.e., contact elements that do not allow sliding nor the separation at the contact region. The stiffness matrix for the contact elements in the local coordinate system is: The rotation matrix [Λ] was used to switch from the local to the global coordinate system. Indeed, since the total potential energy is invariant with respect to the reference system, the global stiffness matrix becomes: (14) where with c and s denoting the trigonometric functions cos and sin. Eventually, the corresponding stiffness matric of contact elements for each pair of coincident nodes was inserted into the global matrix of the assembly (see Figure 17).

Determination of Stiffness, and
The global matrix is a function of the parameters and . To compare the natural frequencies with experimental outcomes, they need to be measured during the complete adhesion of the contact regions. Specifically, by examining the concerning curves (see Figures 9 and 11), the bandwidth in which the frequency has a rather constant trend was selected, to indicate the complete contact adhesion.
The frequencies of FEM model (for each axial traction force) were obtained from the classical linear modal analysis. Considering the global matrices and , the governing equations of motion are: With the following linear harmonic solution: Excluding the trivial solution {q} = 0, other possible solutions (eigenvectors) are obtained for the eigenvalues of the generalized problem , for which the matrix is not invertible. This condition is verified when: from which , are the roots of the characteristic polynomial and the eigenvalues: The shape of each depends on the coefficients , which are nonlinear functions of the parameters and . Thus, the frequencies depend non-linearly on these parameters, which required an iterative local or global method to find the solution ( , ). By making a comparison between the measured frequencies and those from the FEM model, for the first and second flexural modes (1B and 2B): , 0 which can be written in the following vector form: The function does not vanish at any point, i.e., there are no ( , ) such as to satisfy both Equations (20); however, it can be minimized. Therefore, the system of equations becomes: This system was solved with a global search within the domain ( , ) for satisfying the following condition: where || || is the Euclidean norm of the vector and || || is the error, calculated with the theory of error propagation and defines the stress.
Following the procedure explained above, the parameters and are obtained for i-th axial force and reported in Table 1. Eventually the frequencies obtained from the experiment are compared with those from the presented numerical approach and the analytical method in Figure 18. The theoretical data are extracted from the beam model of Euler-Bernoulli (constrained at both ends), modified to consider the stiffening effect of the traction force. The discrepancies between the calculated frequencies from the theoretical model for low traction values can be attributed to the microslip at the contact regions. In fact, that is why the axial traction forces greater than 12 kN are taken into account, where the contacts behave linearly, in the numerical model.

Conclusions
The present work is devoted to identifying the corresponding parameters of a type of contact element, namely normal stiffness, and tangential stiffness. The contact element intends to simulate the vibration of a turbine blade connected to its disc via a dovetail joint. For that matter, an experimental procedure was carried out, with the aim of realizing the modal parameters (frequency and damping) of system that resembles the disc-blade configuration. The first and second flexural modes of the experimental setup were extracted from its free decay, when excited with a frequency close to that of the mode shapes in question. Then, the collected data along with a comprehensive numerical method were utilized for determining the normal and tangential stiffness of the system. Along this path, finite element model of the configuration was built in commercial finite element software, and the corresponding stiffness and mass matrices were imported in a numeric computing environment after applying the Craig-Bampton reduction technique. Once the global matrices of the system were obtained and assembled with the global stiffness matrix of the contact element, the concerning parameters were extracted through an iterative process. Acceptable agreement was found between the obtained values and those from the measurements and analytical models.