A coupled approach to model the effect of wear on the dynamics of the shrouded bladed disk
1 A coupled approach to model the effect of wear on the dynamics of the shrouded bladed disk
Lakshminarayana Reddy Tamatam *, Daniele Botto and Stefano Zucca Abstract
This paper deals with modelling the effect of wear on the dynamics of the shrouded bladed disk with frictional contacts at the shrouds and the contact interface evolution. Prediction of fretting wear commonly occurring at the contacts of turbomachinery components, and its impact on the dynamics is increasingly researched due to the components subjected to their structural limits for performance and operating at high loading conditions. Over a lifetime, the fretting wear at these contacts could alter the global dynamic response of these bladed disks from the designed operating point and could lead to high vibration amplitudes. This study implements a coupled static/dynamic harmonic balance method (HBM) with wear energy approach and an adaptive wear logic to study the impact on the steady-state nonlinear dynamic response. Firstly, the methodology is applied to a cantilever beam with a contact patch and then to a shrouded bladed disk with shroud contacts using cyclic symmetry boundary conditions. The novelty of the paper is to show the effect of wear on the dynamics with the changing contact pre-load using a coupled approach. A wear acceleration parameter ๐ฏ ๐ค,๐๐๐ฅ is defined, and the impact of the choice of this parameter is discussed in detail on the computation time and the accuracy of results. The test cases demonstrate the impact of wear on the dynamic nonlinear response curves and the contact interface evolution for various scenarios. Keywords: non-linear structural dynamics, contact evolution, contact pre-load, fretting wear, harmonic balance method
1. Introduction
Mechanical assemblies made of two or more components have contact interfaces between them. Many a time relative sliding occurs during operation, sometimes intentional for design and functional purposes and at times unintentional. In turbomachinery components such as bladed disks, the relative sliding at blade roots, shroud contacts and under platform dampers are carefully designed to provide friction damping and to reduce the vibration amplitudes. These components have very high modal density, and resonances are unavoidable within the operating range. The contacts are precisely designed to undergo small relative motion and provide friction damping and reduce the peak vibration amplitudes. Consequently, this high-frequency small-amplitude motion leads to fretting wear during operation. This wear affects the global dynamics over a period of time and causes an undesirable shift in resonant frequency and amplitude. The turbomachinery components are designed to perform at its structural limits to obtain the maximum efficiency. Hence, they experience very high static and dynamic loads due to centrifugal forces, thermal strains and fluctuating gas forces, leading to high cycle fatigue (HCF). HCF is one of the common failure modes of the bladed disks. This leads to degraded performance, sometimes even catastrophic failure. Wear is a multi-physics and multi-scale phenomenon. The micro-slip or gross tangential slip occurring at the contact interfaces leads to contact hysteresis and energy dissipation and cause wear [1], [2]. There is ample literature available on characterizing and predicting fretting wear occurring at these contacts for various contact Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy * e-mail: [email protected] conditions ([3]โ[9]). This microscopic level surface changes can alter the dynamic system response in a significant way even for a relatively simple system, as shown experimentally in Fantetti et al. [10]. The friction contacts introduce nonlinearity to the system. To perform dynamic analysis considering the nonlinearity, one needs robust iterative solvers and accurate contact models to solve the differential equations of motion with friction nonlinearity. The two predominant methods are time integration methods ([11], [12]) and frequency-based methods ([11], [13]โ[19]). Harmonic balance method (HBM) is the most widely used frequency-based method to solve the non-linear differential equations due to its suitability for large systems, computational efficiency and directly obtaining steady-state response ([11], [13]โ[18]). Todayโs need is embedding the accurate wear prediction tool into these non-linear dynamic solvers to study the impact of wear and contact interface evolution. There is little research available on this joint topic. Brakeโs book [20] provides an up-to-date comprehensive overview of the developments of mechanics of jointed structures so far. Some of the other notable researches are [10], [21]โ[28]. Sfantos et al. ([3], [4]) have proposed the Boundary Element Method (BEM) techniques for wear simulation and optimization methods using incremental sliding BEM. Salles et al. ([24]โ[26]) proposed DLFT and HBM techniques to compute the effect of wear on the dynamics. They also further proposed a multi-scale method to predict the effect of fretting wear as the amount of wear in each fretting cycle is too small relative to the change in dynamics. A multi-scale approach would work on two different time scales for wear and the dynamics. Armand et al. ([21], [22]) extended the previous work. They proposed a multi-scale approach in time and space to predict the effect of fretting wear on the dynamic response for bladed disk under platform dampers application. Later, they also modelled the contact considering surface roughness and showed the change in the dynamic response of the system [23]. Albeit all these studies are only numerical so far. To the authorโs knowledge, the only research that studied the effect of wear on the evolution of dynamics of friction contact experimentally and characterized by a constitutive numerical model was by Fantetti et al. [10]. Fantetti et al. showed the impact of wear on the dynamics of structures experimentally and then validated using a constitutive numerical model - modified Bouc-Wen model. They studied the evolution of fretting wear and the dynamics during the running-in up to steady-state. They captured the evolution of fretting wear for a very short number of cycles up to a few millions of cycles. The study clearly showed the impact of the microscopic wear generated at the contact has on the resonant frequencies and damping and the corresponding change in dynamics in the running-in and steady-state conditions. The authors consider this research as the gateway to understand the effect of wear on the dynamics and for future predictive modelling. Overall, there is a need for a reliable and high-fidelity predictive model to include wear and to predict the impact of wear on the global dynamics of the systems. The current research is a numerical work and aims at developing and implementing a method to predict the effect of wear on the dynamics of structures with contact interfaces and also the evolution of contact interface using a proven static/dynamic coupled approach with harmonic balance method (HBM) [27]. Concerning the previous works existing in the literature ([29], [30]), the current approach allows the automatic update of the static pre-load distribution over the contact area during the non-linear dynamic analysis, without any need for a separate static analysis. The coupled approach is more relevant in shrouded blades, where the static pre-load at the contacts is strongly affected by the worn-out material at the shrouds, hence the vibration amplitudes. This is different from what happens at the under-platform dampers and blade root joints, where wear mostly affects the static load distribution. The wear is computed using wear energy approach [31], and an accelerated adaptive wear logic is defined. Two test cases are chosen to demonstrate the method โ a cantilever beam with a contact patch and a shrouded bladed disk with shroud contacts and cyclic symmetry boundary conditions. The novelty of this paper is two-fold: 1) Modelling the effect of wear on the dynamic forced response and the contact interface evolution with โchanging contact pre-loadโ. The impact of the choice of user-defined wear acceleration parameter ๐ฏ ๐ค,๐๐๐ฅ that is described in the paper.
2. Methodology
The contact interface introduces nonlinearity to the system of equations. Two ways to solve the non-linear differential equation is using time domain and frequency domain methods. A time-domain method such as Direct Time Integration (DTI) provides the transient as well as steady-state response, but is computationally demanding and is not a feasible solution for practical scenarios and systems with large degrees of freedom (DOF). The state-of-the-art frequency-domain method to compute the non-linear response is the Harmonic Balance Method (HBM). This method assumes a periodic response under periodic excitation. Krack and Grossโs book [32] provides a detailed description of the method and examples of non-linear vibration problems. This method offers speedy solution times assuming steady-state response and much less computationally demanding relative to the time domain methods. In our case, this is better as we are interested only in the steady-state response of the system subjected to a periodic excitation. The current study uses a more recent formulation of the HBM, which allows performing a coupled static and dynamic analysis of the system. This coupled method is proven to provide accurate results in comparison to the classical uncoupled approach [29]. It is comparable to the accuracy of direct time integration results with a sufficient number of harmonics.
The governing equation of motion of a generic system with contact interfaces is written as:
๐๐ฬ (๐ก) + ๐๐ฬ (๐ก) + ๐๐(๐ก) = ๐ (๐ก) + ๐ ๐ (๐, ๐ฬ , ๐ก) (1) where ๐ โ โ
NรN is the mass matrix;
๐ โ โ
NรN is the viscous damping matrix;
๐ โ โ
NรN is the stiffness matrices;
๐(๐ก) โ โ
Nร1 is a nodal displacement vector; ๐ฬ is the time derivative of ๐ ; ๐ (๐ก) โ โ
Nร1 is the excitation force vector and ๐ ๐ (๐, ๐ฬ , ๐ก) โ โ Nร1 is the nonlinear contact force interaction vector. The components of the displacement vector ๐ is given by: ๐ = [๐ช ๐๐ ๐ช ๐๐ ๐ช ๐๐ ] ๐ (2) where ๐ช ๐ , ๐ช ๐ and ๐ช are the displacements of contact node DOFs, accessory node DOFs, and other node DOFs of the system respectively. ๐ช ๐ node DOFs include input and output node DOFs. Since the contact DOF ๐ช ๐ and the accessory DOF ๐ช ๐ are usually much smaller in comparison to the other DOF ๐ช ๐ , a reduction of the size of the nonlinear model is typically performed before computing the nonlinear solution. In this paper, the classical Craig-Bampton reduction method [33] is applied, retaining as physical coordinates the contact DOFs ๐ช ๐ and accessory DOFs ๐ช ๐ โ cumulatively called as master nodes, and adding a certain number of fixed-interface modes ๐ผ . ๐ช = [๐ช ๐๐ ๐ช ๐๐ ๐ผ ๐ ] ๐ (3) where ๐ช โ โ nร1 and n โช N . The equation (1) can be rewritten as: ๐ฆ๐ชฬ (๐ก) + ๐๐ชฬ (๐ก) + ๐ค๐ช(๐ก) = ๐(๐ก) + ๐ ๐ (๐ช, ๐ชฬ , ๐ก) (4) where ๐ฆ โ โ nรn is the mass matrix, ๐ โ โ nรn is the viscous damping matrix and ๐ค โ โ nรn is the stiffness matrix, ๐ช(๐ก) โ โ nร1 is a nodal displacement vector; ๐(๐ก) โ โ nร1 is the excitation force vector; ๐ ๐ (๐ช, ๐ชฬ , ๐ก) โ โ nร1 is nonlinear contact force interaction vector obtained after performing the Craig-Bampton reduction method. To solve the equation (4) for periodic excitation using HBM ([32]), the periodic quantities (i.e. displacements and forces) with an angular frequency of ๐ are expressed as truncated series of harmonic terms: ๐ช(๐ก) = โ ( โ ๐ชฬ (โ) ๐ ๐โ๐๐ก๐ปโ = 0 ) ; ๐(๐ก) = โ ( โ ๐ฬ (โ) ๐ ๐โ๐๐ก๐ปโ = 0 ) ; ๐ ๐ (๐ช, ๐ชฬ , ๐ก) = โ ( โ ๐ฬ ๐(โ) (๐ชฬ)๐ ๐โ๐๐ก๐ปโ = 0 ) (5) where โ(โ) represents the real part of the quantity. By applying Galerkinโs procedure with the multi-harmonic approximation, the time domain nonlinear differential equation (4) is transformed into a nonlinear algebraic equation with Fourier coefficients as defined in equation (5), and the balance equation is written as: ๐ (โ) ๐ชฬ (โ) = ๐ฬ (โ) + ๐ฬ ๐(โ) ๐ค๐๐กโ โ = 0. . ๐ป (6) where ๐ (โ) = (โ(โ๐) ๐ฆ + ๐โ๐๐ + ๐ค) is the โ ๐กโ dynamic stiffness matrix of the system. The equation (6) consists the static ( โ = 0 ) and dynamic ( โ = 1. . ๐ป ) equations of the system coupled to each other by Fourier coefficients of the nonlinear contact force ๐ฬ ๐ depending on the Fourier coefficients of the displacement ๐ชฬ . The number of harmonics ๐ป is chosen in a way to approximate the dynamics of the structure with sufficient accuracy. The residual equation is formulated by rewriting equation (6) as: ๐๐๐ (โ) = ๐ (โ) ๐ชฬ (โ) โ ๐ฬ (โ) โ ๐ฬ ๐(โ) ๐ค๐๐กโ โ = 0. . ๐ป (7) Let us consider a shrouded bladed disk with shroud contact interface, as shown in Figure 1. The bladed disk is a tuned system because it is made of identical fundamental sectors arranged in a cyclical manner. Since the bladed disks are cyclic in nature, one can exploit the analysis by reducing the full bladed disk model to a single fundamental sector, by applying proper boundary conditions at the sector interfaces in the nonlinear dynamic analysis ([17], [34]). Figure 1: (a) Full bladed disk with shroud contacts (b) A fundamental blade-disk sector (c) Schematic view of cyclic symmetry for bladed disk model The equation of motion for a tuned bladed disk with contact interfaces with ๐ ๐ sectors can be obtained by modifying equation (1) as: ๐ ๐ฬ (๐) (๐ก) + ๐ ๐ฬ (๐) (๐ก) + ๐ ๐ (๐) (๐ก) = ๐ (๐) (๐ก) + ๐ ๐,๐(๐) (โ๐ช ๐,๐ , โ๐ชฬ ๐,๐ , ๐ก) + ๐ ๐,๐(๐) (โ๐ช ๐,๐ , โ๐ชฬ ๐,๐ , ๐ก) (8) where ๐ = 1. . ๐ ๐ is a sector number; ๐ โ โ NรN , ๐ โ โ
NรN , and
๐ โ โ
NรN denote the mass matrix, viscous damping matrix and the stiffness matrix respectively of ๐ ๐กโ isolated bladed disk segment without cyclic boundary conditions; ๐ (๐) (๐ก) โ โ Nร1 is a nodal displacement vector; ๐ (๐) (๐ก) โ โ Nร1 is the excitation force vector acting on ๐ ๐กโ blade, and ๐ (๐) ๐,๐ , ๐ (๐) ๐,๐ โ โ Nร1 are nonlinear contact force interaction vectors at the shrouds of the left and the right adjacent sectors respectively depending on the relative displacements and velocities. These contact forces depend on the relative displacements in a nonlinear fashion at the shroud contact. Figure 2: A schematic view of the isolated bladed disk segment The displacement vector at ๐ ๐กโ isolated bladed disk segment is written as: ๐ (๐) = [ ๐ช ๐,๐(๐) ๐ ๐ช ๐,๐(๐) ๐ ๐ช ๐(๐) ๐ ๐ช ๐,๐(๐) ๐ ๐ช ๐,๐(๐) ๐ ๐ช ๐(๐) ๐ ] ๐ (9) where ๐ช ๐,๐(๐) and ๐ช ๐,๐(๐) are the nodal displacements of the left boundary segment interface and right boundary segment interface, respectively, as shown in Figure 2. ๐ช ๐,๐(๐) and ๐ช ๐,๐(๐) are the displacements of left and right shroud contact nodes, respectively. ๐ช ๐(๐) are the accessory nodes which contains the displacements of loading and response nodes, and ๐ช ๐(๐) contains the displacements of other interior nodes of the system of the ๐ ๐กโ segment. Similar to the generic system, to solve the equation (8) for periodic excitation using HBM ([32]), the periodic quantities (i.e. displacements and forces) with an angular frequency of ๐ are expressed as truncated series of harmonic terms: ๐ (๐) (๐ก) = โ ( โ ๐ฬ (โ)(๐) ๐ ๐โ๐๐ก๐ปโ = 0 ) ; ๐ (๐) (๐ก) = โ ( โ ๐ ฬ (โ)(๐) ๐ ๐โ๐๐ก๐ปโ = 0 ) ; (10) ๐ ๐,๐(๐) (โ๐ช ๐,๐ , โ๐ชฬ ๐,๐ , ๐ก) = โ ( โ ๐ ฬ ๐,๐(โ)(๐) ๐ ๐โ๐๐ก๐ปโ = 0 ) ; ๐ ๐,๐(๐) (โ๐ช ๐,๐ , โ๐ชฬ ๐,๐ , ๐ก) = โ ( โ ๐ ฬ ๐,๐(โ)(๐) ๐ ๐โ๐๐ก๐ปโ = 0 ) By applying Galerkinโs procedure with the multi-harmonic approximation, the time domain nonlinear differential equation (8) is transformed into a nonlinear algebraic equation with Fourier coefficients as defined in equation (10), and the balance equation is written as: ๐ (โ) ๐ฬ (โ)(๐) = ๐ ฬ (โ)(๐) + ๐ ฬ ๐,๐(โ)(๐) + ๐ ฬ ๐,๐(โ)(๐) ๐ค๐๐กโ โ = 0. . ๐ป (11) where ๐ (โ) = (โ(โ๐) ๐ + ๐โ๐๐ + ๐) is the โ ๐กโ dynamic stiffness matrix of the system. The equation (11) consists of the static ( โ = 0 ) and dynamic ( โ = 1. . ๐ป ) equations of the system coupled to each other by Fourier coefficients of the nonlinear contact force ๐ ฬ ๐,๐(โ)(๐) and ๐ ฬ ๐,๐(โ)(๐) depending on the Fourier coefficients of the relative displacements and velocities. The number of harmonics ๐ป is chosen in a way to approximate the steady state dynamics of the structure with sufficient accuracy. The equation (11) has to be solved iteratively for the unknown displacement vector of ๐ ๐กโ segment because of the nonlinear contact forces. Note the size of the dynamic stiffness matrix ๐ โ โ (๐(๐ป+1)ร๐(๐ป+1)) and the force vectors ๐ ฬ (๐) , ๐ ฬ ๐,๐(๐) , ๐ ฬ ๐,๐(๐) โ โ (๐(๐ป+1)ร1) . The system of equations still refer to the single fundamental bladed disk sector. The cyclic symmetry boundary conditions will be implemented in the next steps. The turbine bladed disk is assumed to be subjected to a travelling wave type excitation, see ref. [17], [34] for detailed explanation. The inter-blade phase angle (IBPA) is defined as: โ = . ๐ธ๐, where ๐ is the number of blades and ๐ธ๐ is the engine order that, in the Campbell diagram, crosses the resonance under investigation. Reducing the number of DOFs required for the nonlinear forced response analysis is done in two steps. First, the disk interface segment boundary DOFs are reduced using cyclic symmetry constraints and then similar constraints are applied to the shroud contact node DOFs. Reducing the disk segment boundary DOFs using cyclic symmetry constraints:
Since the bladed disk is subjected to the travelling wave type excitation, it is assumed that the steady-state vibration response also represents a travelling wave. This results in a relationship between the right segment interface nodal DOF, and the left segment interface nodal DOF. Owing to the cyclic symmetry conditions based on the assumption of the turbine bladed disk is made of the similar bladed-disk fundamental sector, one can write the following coupling relationship for the displacements at the segment boundary DOFs: ๐ชฬ ๐,๐(โ)(๐) = ๐ชฬ ๐,๐(โ)(๐) ๐ โ๐โโ (12) Redefining the reduced displacement vector of unknowns by considering the right-side boundary segment for โ ๐กโ harmonic component and ๐ ๐กโ sector as: ๐ฬ (โ)(๐) = ๐ฬ ๐ถ๐(โ) ๐ฬ ๐๐๐(๐) (โ) (13) Thereby making the new reduced displacement vector of unknowns as: ๐ฬ ๐๐๐(๐) (โ) = [ ๐ชฬ ๐,๐(โ)(๐) ๐ ๐ชฬ ๐(โ)(๐) ๐ ๐ชฬ ๐,๐(โ)(๐) ๐ ๐ชฬ ๐,๐(โ)(๐) ๐ ๐ชฬ ๐(โ)(๐) ๐ ] ๐ (14) ๐ฬ ๐ถ๐(โ) = [ ๐(๐ โ๐โโ ) 0 0 0 0๐ 0 0 0 00 ๐ 0 0 00 0 ๐ 0 00 0 0 ๐ 00 0 0 0 ๐ ] (15) Where ๐ฬ ๐ถ๐(โ) is the cyclic symmetry constraint matrix for โ ๐กโ harmonic component; ๐ is the identity matrix of an appropriate size corresponding to the size of DOFs. Substituting the vector ๐ฬ (โ)(๐) in equation (11) by the cyclic constraint relation as given in equation (13) for โ = 0. . ๐ป and left-multiplying by the corresponding complex conjugate transpose of ๐ฬ ๐ถ๐(โ) results in: ๐ฬ ๐๐๐ (โ) ๐ฬ ๐๐๐(๐) (โ) = ๐ ฬ (โ)๐๐๐(๐) + ๐ ฬ ๐,๐(โ)๐๐๐(๐) + ๐ ฬ ๐,๐(โ)๐๐๐(๐) ๐ค๐๐กโ โ = 0. . ๐ป (16) where ๐ฬ ๐๐๐ (โ) = (โ(โ๐) ๐ฬ (โ) + ๐โ๐๐ฬ (โ) + ๐ฬ (โ) ) is the โ ๐กโ dynamic stiffness matrix of the system. Note that the complex mass matrix ๐ฬ โ โ (๐(๐ป+1)ร๐(๐ป+1)) , viscous damping matrix
๐ฬ โ โ (๐(๐ป+1)ร๐(๐ป+1)) , and the stiffness matrix
๐ฬ โ โ (๐(๐ป+1)ร๐(๐ป+1)) are Hermitian because of the application of the cyclic symmetry constraints and ๐ ฬ ๐๐๐(๐) , ๐ ฬ ๐,๐๐๐๐(๐) , ๐ ฬ ๐,๐๐๐๐(๐) โ โ (๐(๐ป+1)ร1) where ๐(๐ป + 1) < ๐(๐ป + 1) because of the successful elimination of the left boundary segment interface DOFs. Reducing the shroud contact DOFs using cyclic symmetry constraints:
The displacement vector of unknowns can further be reduced by applying the cyclic symmetry constraints at the shroud contact DOFs and writing in relative displacement terms. Writing in relative displacements terms reduces the number of shroud contact DOFs by half as one can consider only the number of shroud contact DOFs of one side and reduce the size of the system required for nonlinear analysis. The relative displacement at the right-side shroud contact and the non-linear contact force can be written as: โ๐ชฬ ๐,๐(โ)(๐) = ๐ชฬ ๐,๐(โ)(๐) โ ๐ชฬ ๐,๐(โ)(๐โ1) ; โ๐ชฬ ๐,๐(โ)(๐) = ๐ชฬ ๐,๐(โ)(๐) โ ๐ชฬ ๐,๐(โ)(๐) ๐ ๐โโ ๐ฬ ๐,๐(โ)(๐) = โ ๐ฬ ๐,๐(โ)(๐) ๐ ๐โโ (17) The new displacement vector of unknowns in terms of relative displacements at the shroud contact becomes: ๐ฬ ๐๐๐(๐) (โ) = [ ๐ชฬ ๐,๐(โ)(๐) ๐ ๐ชฬ ๐(โ)(๐) ๐ โ๐ชฬ ๐,๐(โ)(๐) ๐ ๐ชฬ ๐(โ)(๐) ๐ ] ๐ (18) This leads to cutting down the shroud contact unknown DOFs by half by considering either right or left contact, hence reducing the size of the non-linear computation. Similarly, a relative displacement cyclic symmetry transformation matrix [35] can be defined, and the balance equation (16) can be rewritten as: ๐ฬ ๐๐๐ (โ) ๐ฬ ๐๐๐(๐) (โ) = ๐ ฬ (โ)๐๐๐(๐) + ๐ ฬ ๐,๐(โ)๐๐๐(๐) ๐ค๐๐กโ โ = 0. . ๐ป (19) where ๐ฬ , ๐ฬ, ๐ฬ โ โ (๐(๐ป+1)ร๐(๐ป+1)) and ๐ ฬ ๐๐๐(๐) , ๐ ฬ ๐,๐๐๐๐(๐) โ โ (๐(๐ป+1)ร1) where ๐(๐ป + 1) < ๐(๐ป + 1) . Since the number of shroud contact nodes, boundary nodes and accessory nodes DOFs are much smaller than the number of internal DOFs, as described for the generic system, the classical Craig-Bampton reduction method [33] is applied, retaining as physical coordinates the right segment interface nodes, right shroud contact nodes, accessory nodes DOFs โ collectively called as master nodes, and adding a certain number of fixed-interface modes ๐ผ . ๐ชฬ (๐) (โ) = [ ๐ชฬ ๐,๐(โ)(๐) ๐ ๐ชฬ ๐(โ)(๐) ๐ โ๐ชฬ ๐,๐(โ)(๐) ๐ ๐ผ (๐) ๐ ] ๐ (20) where ๐ชฬ (๐) โ โ ๐(๐ป+1)ร1 and ๐(๐ป + 1) โช ๐(๐ป + 1) . Finally, the equation (19) can be reformulated using equation (20) as: ๐ฬ (โ) ๐ชฬ (๐) (โ) = ๐ฬ (โ)(๐) + ๐ฬ ๐,๐(โ)(๐) ๐ค๐๐กโ โ = 0. . ๐ป (21) By rearranging equation (21), the residual equation is written as: ๐๐๐ (โ) = ๐ฬ (โ) ๐ชฬ (๐) (โ) โ ๐ฬ (โ)(๐) โ ๐ฬ ๐,๐(โ)(๐) ๐ค๐๐กโ โ = 0. . ๐ป (22) To use the nonlinear solvers, the residual has to be supplied in the real form. Hence, the complex form residual equation (7) and equation (22) is split into its real and imaginary components to real form to minimize the residual to an acceptable tolerance and find the unknown displacement vector as: ๐๐๐ = [๐๐๐ (0) , ๐ฝ(๐๐๐ (1) ), ๐ด(๐๐๐ (1) ) โฆ ๐ฝ(๐๐๐ (๐ป) ), ๐ด(๐๐๐ (๐ป) )] T (23) A contact model is necessary to compute the nonlinear contact forces mentioned in the previous section. There are various node-to-node and patch-to-patch contact models available in the literature ([36], [37], [46]โ[48], [38]โ[45]) such as Coulomb slider, Jenkins element โ 2D and 3D with constant and variable normal load, Iwan model and its variations, Bouc-Wen model, Valanis model, LuGre model etc. In our study, we chose the state-of-the-art node-to-node 2D Jenkins element with a variable normal load to compute the contact forces, as shown in Figure 3(a). A typical hysteresis loop generated for the given input force and displacements is shown in Figure 3(b). The solution to the balance equations requires the contact forces on the contact interface as input. The contact element is characterized by two linear springs in the tangential and normal direction with tangential ( ๐ ๐ก ) and normal ( ๐ ๐ ) contact stiffness, at each node pair over the contact interface. The contact model allows to characterize and simulate three possible contact states โ stick, slip and lift-off. Figure 3: (a) Jenkinโs element contact model with variable normal load, (b) A typical hysteresis loop At each element, the contact forces โ tangential force ๐(๐ก) and the normal force
๐(๐ก) depend on the relative tangential and normal displacements of the corresponding node pair, namely ๐ฎ(๐ก) and ๐ฏ(๐ก) respectively [29]. The friction limit value is defined by the Coulomb law as ๐๐(๐ก) , where ๐ is the coefficient of friction. To consider the friction phenomena occurring at the interface in the tangential direction, a slider is used to connect the two bodies. When the tangential force exceeds the limit value, the slider starts moving, and the amount of slip between the nodes is ๐ฐ(๐ก) . Since the contact is a unilateral constraint, the normal contact force ๐(๐ก) is defined at each time ๐ก as: ๐(๐ก) = max(๐ ๐ ๐ฏ(๐ก), 0) (24) When ๐ฏ(๐ก) is negative, no normal contact force is allowed, and separation occurs, hence referred to as lift-off state. For the tangential direction, the tangential contact force
๐(๐ก) is dependent on the contact state as defined by:
๐(๐ก) = {๐ ๐ก {๐ฎ(๐ก) โ ๐ฐ(๐ก)}๐๐(๐ก)๐ ๐๐๐(๐ฐฬ )0 ๐ ๐ก๐๐๐ ๐ ๐ก๐๐ก๐๐ ๐๐๐ ๐ ๐ก๐๐ก๐๐๐๐๐ก โ ๐๐๐ ๐ ๐ก๐๐ก๐ (25) The contact forces are computed using predictor-corrector logic [17]. At each time ๐ก a predictor step is performed assuming the stick contact conditions: ๐ ๐ (๐ก) = ๐ ๐ก [๐ฎ(๐ก) โ ๐ฐ(๐ก)] = ๐ ๐ก [๐ฎ(๐ก) โ ๐ฐ(๐ก โ โ๐ก)] (26) Where โ๐ก is the time step, then a corrector step is performed to compute the actual value of ๐(๐ก) given by:
๐(๐ก) = { ๐ ๐ (๐ก)๐๐(๐ก)๐ ๐๐๐(๐ ๐ (๐ก)) 0 ๐ ๐ก๐๐๐ ๐ ๐ก๐๐ก๐๐ ๐๐๐ ๐ ๐ก๐๐ก๐๐๐๐๐ก โ ๐๐๐ ๐ ๐ก๐๐ก๐ (27) The slider displacement ๐ฐ(๐ก) is computed as: ๐ฐ(๐ก) = { ๐ฐ(๐ก โ โ๐ก)๐ฎ(๐ก) โ ๐๐(๐ก)๐ ๐๐๐(๐(๐ก))/๐ ๐ก ๐ฎ(๐ก) ๐ ๐ก๐๐๐ ๐ ๐ก๐๐ก๐๐ ๐๐๐ ๐ ๐ก๐๐ก๐๐๐๐๐ก โ ๐๐๐ ๐ ๐ก๐๐ก๐ (28) The approximation introduced due to the predictor-corrector method is minimized by choosing a fine time step โ๐ก (for instance, 2e8 time steps). Figure 4: Alternating Frequency Time (AFT) method The balance equations are solved in the frequency domain using HBM, whereas the accurate contact forces are possible to compute only in the time domain. Hence an Alternate Frequency/Time (AFT) method ([49], [50]) is employed. The relative displacements are converted from the frequency domain to the time domain by applying Inverse Fast Fourier Transform and then run through the contact model. The contact forces are obtained in the time domain. Then applying Fast Fourier Transform to convert to the Fourier coefficient of contact forces in the time domain, as shown in Figure 4. Newton-Raphson logic is implemented for the iterative procedure. The Jacobian matrix needed to compute the residual for the Newton-Raphson method is computationally intensive if using MATLAB built-in finite difference method (FDM) for large systems. Hence an analytical Jacobian is implemented based on the works of [17], [51]. The contact forces and the partial derivatives of the contact forces are computed and assembled in the Jacobian matrix. This dramatically increases the speed of the solution by many folds. The contact forces and computation of the dynamic response are primarily based on the accurate input of the normal and tangential contact stiffness. The closed-form solutions needed to compute the accurate tangential ( ๐ ๐ก ) and normal ( ๐ ๐ ) contact stiffness to provide the contact model is provided based on the analytical equations as mentioned in Ref [52]โ[55]. Then the contact stiffnesses are distributed over the contact area according to the individual elemental contribution. To predict the effect of fretting wear on the dynamic forced response of the bladed disk, wear has to be evaluated. There are many wear models [56] available in the literature, such as Archardโs law [2], wear energy approach [31], thermodynamics approach, etc. Due to its simple nature in the formulation, good correlation with experiments, lower computation effort and easy implementation in the forced response solvers to evaluate the non-linear dynamics, a wear energy approach is used. Many research studies have shown the wear energy approach in good agreement with the experimental results ([6], [7], [57]). Initially, the dynamic forced response is computed with pristine contact surfaces without any wear. The nodal wear depth ๐ฏ ๐ค is defined as: ๐ฏ ๐ค = ๐ ๐ ๐ผ๐ธ๐ด (29) where ๐ด is the area associated with the node pair, ๐ผ is the wear coefficient for a particular contact pair, and ๐ธ is the energy dissipated over one cycle. It is worth mentioning that ๐ด depends on the finite element mesh over the contact surface, while ๐ผ is the experimentally obtained parameter for a specific material couple ([58]โ[60]) and ๐ธ is obtained by evaluating the area under the hysteresis loop. The control parameter of equation (29) is ๐ ๐ โ number of cycles. The process of updating the contact surface due to wear after one cycle is called a wear iteration. With materials typical of turbine blades, the wear depth obtained at each vibration cycle is so small that no significant effect on the forced response of the system is observed if a wear iteration is performed at each vibration cycle. Hence, a wear acceleration technique is used, where the wear depth ๐ฏ ๐ค associated at each node pair is considered for ๐ ๐ number of cycles. In this paper, an adaptive strategy is used, and the value of ๐ ๐ at each wear iteration is computed as: ๐ ๐ = ๐๐๐๐๐ ( ๐ฏ ๐ค,๐๐๐ฅ ๐๐๐ฅ(โโ ๐๐ )) (30) where, ๐ฏ ๐ค,๐๐๐ฅ is a user-defined parameter and โโ ๐๐ is the nodal wear depth at the contact patch for one vibration cycle. Choosing ๐ฏ ๐ค,๐๐๐ฅ is empirical at this stage, as there are no guidelines from the experiments what value is most suitable in terms of accuracy of results. However, we present a few thumb rules one can use to define the parameter ๐ฏ ๐ค,๐๐๐ฅ : โช the maximum wear depth allowed for each wear iteration, directly selected by the user For example, let us consider an assembly with a contact interface is functional up to a contact wear depth of 100 ยตm. In this case, the user can arbitrarily discretize this maximum allowable wear depth into 100 parts and define ๐ฏ ๐ค,๐๐๐ฅ as 1 ยตm. By doing this, the user can visualise the contact evolution and the effect on the non-linear response plots at sufficient intervals before the functional failure of the assembly. โช the percentage of the maximum static deflection for the given static loads acting at the contact This method is more quantitative than the previous arbitrary method. ๐ฏ ๐ค,๐๐๐ฅ is a parameter of the static loading condition. Irrespective of the geometry, size and loading scenario of the system, ๐ฏ ๐ค,๐๐๐ฅ can be defined as a percentage of the maximum static deflection at the contact. For example, ๐ฏ ๐ค,๐๐๐ฅ can be chosen as 0.1% to 10% of the maximum static deflection ( ๐ฟ ๐๐๐ฅ ). โช the maximum tolerable error in wear depth This is more of a special case. Let us consider that the user has a very strict requirement of the accuracy of wear depth needed in each wear iteration. Then ๐ฏ ๐ค,๐๐๐ฅ can be defined directly as the value of this wear depth. In this study, we define ๐ฏ ๐ค,๐๐๐ฅ as the percentage of the maximum static deflection at the contact patch for the given static loading condition. This percentage is varied between 0.5% to 10% of the maximum static deflection ( ๐ฟ ๐๐๐ฅ ). This parameter determines ๐ ๐ , in turn, is a trade-off between the accuracy of the results and computation time needed. A low value of ๐ ๐ will allow for a more accurate, but more time-consuming analysis, the opposite will happen if a high value of ๐ ๐ is selected. Once the term ๐ฏ ๐ค is computed using equation (29), the wear iteration is performed by updating the 0 th order Fourier coefficient of the normal relative displacement at the contact in the following way: ๐ฏ(๐ก) = โ ๐ฏฬ (โ) ๐ ๐โ๐๐ก๐ปโ = 0 โ ๐ฏ ๐ค (31) After the wear iteration, the forced response of the updated system is computed by solving the updated balance equations, taking into account the effect of the wear depth on the static pre-load and the non-linear dynamics of the system simultaneously. It is here worth mentioning that equation (31) allows for a direct update of the governing equations of the system since the static term (h = 0) is included in the set of non-linear equations to solve. This method allows computing the effect of wear on the forced response with โchanging pre-loadโ without the need for a separate non-linear static analysis routine to compute the contact pre-load distribution. The results are shown and discussed in the next section. Figure 5 gives an overview of the steps with pre-processing and solution phase to compute the effect of wear on the non-linear dynamic response. The pre-processing can be done using any standard commercial software such as ANSYS . In our case,
ANSYS V
MATLAB R
3. Test Case โ description, results and discussion
The main test case here is the tuned shrouded bladed disk with shroud contacts with cyclic symmetry boundary conditions. Before jumping on to the test case, the proposed method is first evaluated on a simple cantilever beam with a contact patch on one side and clamped at the other end. The two test cases demonstrate the proposed methodology with changing pre-load and the contact surface evolution during wear. The goal for the test cases is to choose the static load and the excitation in such a way that there is accelerated wear and total loss of contact due to wear with high energy dissipation at each cycle. This will aid validation in terms of the number of wear iterations, wear profile evolution and correlate the wear profile with the analytical free case. Also, the effect of the choice of parameter ๐ฏ ๐ค,๐๐๐ฅ is quantified in terms of accuracy and computational demand. In the following test cases, ๐ฏ ๐ค,๐๐๐ฅ is defined as a percentage of maximum static deflection under static loading conditions. Hence, a preliminary static analysis is performed to compute the maximum static deflection as if the contact was open without the presence of the second body. A steel cantilever beam with a contact patch and loading conditions as shown in Figure 6(a) and (b) and Table 1 is chosen as a test case to demonstrate the impact of wear on the dynamic response. The FE model mesh is shown in Figure 6(c). The contact elements at the contact patch are connected as node-to-node pair to the ground with the 2D Jenkins element contact model with the variable normal load, as shown in Figure 6(d). The beam is loaded at the mid-section so that as the wear progresses, the contact pre-load and the contact zone changes considerably. The contact patch is made up of 169 node pairs, and the grid distribution is as shown in Figure 7 with a response node as the centre node of the contact patch. Maximum static deflection ( ๐ฟ ๐๐๐ฅ ) in this case is considered as the maximum deflection at the contact patch for the applied static load, as shown in Table 1. Table 1: Parameters for the Cantilever beam and the contact patch Parameter Value
Material Steel Youngโs modulus ( ๐ธ ) 210 GPa Density ( ๐ ) 7860 kg/m Beam dimension 200mm L x 20 mm H x 20mm D Contact patch dimension 20mm x 20mm No. of contact elements 169 Tangential contact stiffness ( ๐ ๐ก ) 470 N/ยตm Normal contact stiffness ( ๐ ๐ ) 583 N/ยตm Friction coefficient ( ๐ ) 0.5 Wear energy coefficient ( ๐ผ ) 2e3 ยตm /J Static load ( ๐น ๐ ๐ก๐๐ก๐๐ ) 20 kN Static load/Excitation ratio ( ๐น ๐ ๐ก๐๐ก๐๐ /๐น ๐๐ฅ ) 2 Choice of ๐ฏ ๐ค,๐๐๐ฅ {0.5%, 1%, 2%, 5%, 10%} of ๐ฟ ๐๐๐ฅ Figure 6: (a) 3D Model of a cantilever beam, (b) a close-up view of the contact patch, (c) FE mesh of the beam and (d) contact elements connection Figure 7: Contact patch with 169 elements and the centre response node (red dot)
Static Load
Dynamic
Excitation Contact patch a)c)b)d) v u The beam is excited by a harmonic transverse force, as shown in Figure 8, while the contact pre-load is provided by a vertical static load. Figure 8 also shows the static deflection and the shape at the initial contact condition with the second body and at the free condition when the second body is absent. Figure 8: Excitation mode (first bending mode at 406 Hz) (top) and highlighting static deflection state with and without body 2 (bottom) The FE model is generated using
ANSYS , and the contact interface and boundary conditions are identified and defined in
ANSYS . A Craig-Bampton Component Mode Synthesis (CB-CMS) [33] is performed by defining contact patch nodes and loading nodes as master nodes and retaining first 15 fixed interface modes to reduce the size of the system needed for the non-linear dynamic analysis. Otherwise, the dynamic analysis with all DOFs would be large and prohibitively time-consuming. Reduced mass and stiffness matrices are extracted after the reduction process. The reduced system has 534 DOFs. The solution algorithm is followed, as shown in Figure 5. 0 th and 1 st order harmonics are used to solve the static and dynamic balance equations in the frequency domain using HBM. Analytically computed tangential and normal contact stiffness is computed, taking into account the dimensions of the contact patch. These stiffnesses are distributed over the contact patch according to the elemental area contribution. A preliminary linear analysis is performed to obtain free and stick states of the contact, and the responses at these states are assumed as reference. The computed non-linear response can then be compared with these two reference states. An initial gap/interference, if any, can also be defined using the current algorithm by offsetting in the direction normal to the contact. In the present case, the gap is set to 0 as the beam is resting on top of the body 2 before the static forces are applied to the system. Then both static and periodic forces are applied at the indicated excitation nodes. The contact forces, dissipated energy, contact status of each node and the shear/slip distributions are computed by solving equation (6). The wear energy approach, described in Section 2.5, is implemented to calculate the nodal wear depth at the contact patch for one vibration cycle. As the nodal wear depth is so minuscule at each vibration cycle, the dynamics is hardly affected. An adaptive wear acceleration logic is implemented by defining a parameter ๐ฏ ๐ค,๐๐๐ฅ as described in the previous section. The parameter ๐ฏ ๐ค,๐๐๐ฅ reduces the number of wear iterations, thus reducing the computation time, compared to running dynamic analysis for each vibration cycle. Since ๐ฏ ๐ค,๐๐๐ฅ is a user-defined value, the effect of the choice of ๐ฏ ๐ค,๐๐๐ฅ is also quantified. The iterative loop is performed until the required number of vibration cycles are obtained or until the complete loss of contact. Figure 9 shows the fineness of cumulative wear depth plots at the contact patch for different values of ๐ฏ ๐ค,๐๐๐ฅ . Figure 9: Cumulative wear depth plots at the contact patch for complete loss of contact for different values of ๐ฏ ๐,๐๐๐ Figure 10 shows the response curves of reference free and stick state along with the backbone curve of the maximum amplitude of response at each wear iteration for v ๐ค,๐๐๐ฅ = 5% of ๐ฟ ๐๐๐ฅ . It is evident and logical that, as the wear progresses, the contact interface loosens, and the contact tends towards the free state. The current formulation automatically takes into consideration the changing pre-load and is reflected in the response plots. Figure 11 provides an insight into the physical state of the system, cumulative wear depth and the contact status at different wear iterations. The wear iterations are chosen at the beginning, intermediate and at the end of the wear just before the complete loss of contact. (a) (b)(c) (d)(e) (a) (b) (c) (d) (e) v ๐ค,๐๐๐ฅ ๐ฟ ๐๐๐ฅ = Figure 10: Response plot showing the backbone of the non-linear response obtained as wear progresses around first bending mode with reference free and stick states for ๐ฏ ๐,๐๐๐ = 5% of ๐น ๐๐๐ Figure 11: A result matrix showing - starting, intermediate and ending wear iteration just before the full loss of contact for the given test case highlighting the beam physical state, cumulative wear depth until that particular wear iteration and the contact status at that wear iteration
Wear iteration 1(start)Wear iteration 164Wear iteration 410 (full loss of contact) ๐ ๐ก๐๐ก๐๐ ๐น ๐๐ฅ = ; v ๐ค,๐๐๐ฅ ๐ฟ ๐๐๐ฅ = 1 Contact status plot:
StaticLoad Dynamic
Excitation
Cumulativewear plot:
Dynamic Excitation
TOP VIEW
Static Load : s epa r a t i on - s li p - s t i ck : s epa r a t i on - s li p5 : s li p 2 : s li p - s t i ck : s epa r a t i on Figure 12 shows the plot of the number of vibration cycles at each wear iteration versus the number of wear iterations. For the shown plot, it took 183 wear iterations with the choice of ๐ฏ ๐ค,๐๐๐ฅ = 5% of ๐ฟ ๐๐๐ฅ . In other words, with the allowed wear depth of 5% of maximum static deflection at each wear iteration, it took 183 wear iterations to completely lose the contact with the second body. Depending on the loading configuration, the shape of this graph can be different. In this example, as shown in Figure 11 of wear iteration 1, the contact started with a line contact at the inner edge of the contact between the body 1 and body 2, and proceeded to establish an area contact with an increase in wear iterations. The close-up of the contact state at the beginning and at the end, is also shown in Figure 12. The criteria to complete one wear iteration is the amount of wear depth achieved, which is predetermined by the user input. At the start of the wear and towards the end, a relatively high number of cycles are needed to achieve this limit of wear depth to satisfy the wear iteration because of partial contact. However, when the contact is fully engaged, lower number of cycles can satisfy the wear depth criteria because of large contact, hence more energy dissipation. Figure 12: Wear iteration versus the number of vibration cycles plot showing the trend as the wear progresses from the start state to the complete loss of contact for ๐ฏ ๐,๐๐๐ = 5% of ๐น ๐๐๐ [Note: the number of cycles is adaptive at each wear iteration] Figure 13 shows the maximum response amplitude at a particular wear iteration versus the cumulative number of cycles. It is evident the effect of ๐ฏ ๐ค,๐๐๐ฅ on the number of cycles at which the maximum response starts increasing. The accuracy of the proposed method for different values of v ๐ค,๐๐๐ฅ is here investigated, by using as a reference the analysis performed with ๐ฏ ๐ค,๐๐๐ฅ /ฮด ๐๐๐ฅ = 0.5%, since no significant effects was observed by further reducing the value of v ๐ค,๐๐๐ฅ . Ideally, the best benchmark would be the response computed with ๐ ๐ = 1 (one wear iteration per vibration cycle), but it is prohibitively time-consuming. It is clear that the choice of ๐ฏ ๐ค,๐๐๐ฅ is the trade-off between the computational cost and the accuracy of results; the lower the value of ๐ฏ ๐ค,๐๐๐ฅ , the smoother the wear evolution and the contact patch wear profile is, due to the higher number of wear iterations. In general, a reasonable strategy should be based on a convergence analysis of the response as the value of ๐ ๐ is progressively reduced. In this particular case, for relatively high values of ๐ฏ ๐ค,๐๐๐ฅ (from 2% to 10% of ฮด ๐๐๐ฅ ) the analysis predicts the rapid increase in maximum response to occur earlier than the benchmark. Interestingly, when ๐ฏ ๐ค,๐๐๐ฅ /ฮด ๐๐๐ฅ = 1%, opposite results are obtained, showing that non-monotonic trend of the response during a convergence analysis can be observed. Loosing pre-load
Start state End state Figure 13: Maximum response at each wear iteration vs the cumulative number of cycles for different values of ๐ฏ ๐,๐๐๐ Table 2 summarizes the comparison of the effect of the choice of ๐ฏ ๐ค,๐๐๐ฅ on the computational time, the number of wear iterations necessary for the complete loss of contact, maximum wear depth and the error in the maximum wear depth taking ๐ฏ ๐ค,๐๐๐ฅ = 0.5% as the reference. The trade-off effect of ๐ฏ ๐ค,๐๐๐ฅ is evident by comparing computational time vs error in wear depth. Table 2: Comparison of the impact of wear analysis with non-linear dynamic response solver for different values of ๐ฏ ๐,๐๐๐ ๐ฏ ๐ค,๐๐๐ฅ ๐ฟ ๐๐๐ฅ Computational time
Total no. of wear iterations Max. total wear depth (mm) Error in wear depth (%) 0.5 % 19 hr 32 min 577 6.3919 - 1 % 16 hr 10 min 411 6.4230 0.49 2 % 10 hr 17 min 291 6.4849 1.45 5 % 6 hr 0 min 183 6.6260 3.66 10 % 3 hr 57 min 130 6.9668 8.99 A tuned shrouded bladed disk consisting of 40 blades with identical blades is considered, and the fundamental blade sector is as shown in Figure
14 and described in Table 3. This bladed disk is used as a test case to demonstrate the impact of wear on the dynamics with cyclic symmetry property. A fundamental sector consisting of only one blade is assumed with cyclic symmetry boundary conditions to reduce the size of the problem for non-linear analysis. The disk is assumed to be infinitely rigid and so fixed boundary conditions are applied at the blade root. In this way, the cyclic symmetry only operates through the blade-to-blade coupling at the shrouds. Maximum static deflection ( ๐ฟ ๐๐๐ฅ ) in this case is computed as the maximum deflection at the shroud contact patch for the applied static load in such a way to produce a twisting effect of the blade. Table 3: Parameters for the shrouded bladed disk and the contact patch Parameter Value
Material Steel Youngโs modulus ( ๐ธ ) 210 GPa Density ( ๐ ) 7860 kg/m Blade geometry As shown in Figure 14 Number of blades 40 Contact patch dimension 20mm x 30mm Contact elements 20 on the left and 20 on the right shroud Tangential contact stiffness ( ๐ ๐ก ) 93 N/ยตm Normal contact stiffness ( ๐ ๐ ) 113 N/ยตm Friction coefficient ( ๐ ) 0.5 Wear energy coefficient ( ๐ผ ) 2e3 ยตm /J Static load ( ๐น ๐ ๐ก๐๐ก๐๐ ) 60 kN Static load/Excitation ratio ( ๐น ๐ ๐ก๐๐ก๐๐ /๐น ๐๐ฅ ) 12 Choice of ๐ฏ ๐ค,๐๐๐ฅ {0.5%, 1%, 2%, 5%, 10%} of ๐ฟ ๐๐๐ฅ Figure 14: FE model of the blade sector, boundary and loading conditions Static torque is applied to the blade by two static forces ๐น ๐ ๐ก๐๐ก๐๐ , to simulate the twisting effect, actually due to the centrifugal force in rotating blades, which produces the static pre-load at the shrouds. A concentrated periodic excitation ๐น ๐๐ฅ is also applied at a mid-span node of the airfoil. Figure 15(a) highlights the shroud contact patch consisting of 20 contact elements each side with a 5x4 grid. Figure 15(b) shows the typical contact elements distribution between right contact patch with an adjacent blade left contact patch. The friction coefficient, contact stiffnesses and wear energy coefficient is defined as mentioned in Table 3. The non-linear dynamic analysis is performed by retaining the 0 th and the 1 st harmonics in the balance equation. First bending mode with first Engine Order (EO = 1) as shown in Figure 16 is chosen as the frequency of interest to analyse the effect of wear on the dynamics of the bladed disk. Figure 15: (a) Highlight of the right shroud contact showing 20 contact nodes (b) Pictorial representation of the applied node-to-node contact model To define the adaptive wear parameter ๐ฏ ๐ค,๐๐๐ฅ , the following logic is applied. The static deflection of the blade is computed under the action of the static forces. Then the average normal displacement of each contact surface is obtained. This value has been set as the maximum possible wear on each surface ๐ฏ ๐ค,๐๐๐ฅ , when such a depth is worn out the blades would vibrate freely without any contact with the adjacent blades. Finally, ๐ฏ ๐ค,๐๐๐ฅ is varied from 0.5% to 10% in steps indicated in Table 3 to investigate its effect on the accuracy of the model in predicting the evolution of the system dynamics. Figure 16 shows the nodal diameter versus frequency relationship up to 15 ND for the stuck case. Figure 17 shows the response for various engine order excitations along with the linear free state and stick state curves for the indicated loads. Figure 16: Nodal diameter diagram of the stuck contact shrouded blade (a)(b) Figure 17: Response to various EO excitations with reference linear free and stick state peaks Figure 18 shows the non-linear FRF for the given excitation forces with fixed static load and varying excitation loads. The figure also includes a reference linear free state and stick state for comparison. With the increasing excitation force, the resonance frequency moves towards the free state as the contact nodes experience a larger slip. Figure 18: FRF for various excitation loads with a fixed static load (EO = 1) Figure 19 shows the cumulative wear plot and the contact status plot at various wear iteration intervals. For the given loading conditions, the cumulative wear at the right shroud contact interface, the spatial wear distribution and the contact status at each node can be visualized. At the beginning of the wear, the contact is in the slip-stick state with uniform pressure distribution. As wear progresses, some of the contact nodes are in separation-slip-stick and eventually separation-slip before the complete loss of contact. For simplicity, only one loading condition is shown. The same procedure can be used to visualize various loading conditions and the effect of the adaptive logic parameter ๐ฏ ๐ค,๐๐๐ฅ with an accurate insight at the contact interface behaviour. Figure 19: Tabular representation of the cumulative wear plot and contact status plot of the right shroud contact at different wear iterations until the loss of contact for given loading conditions Figure 20 shows the response with the impact of wear and the backbone of the response curves with the progress of wear. The contact pre-load is continuously changing as the contact is worn and loosens. This loss of contact impacts the dynamic behaviour of the bladed disk. The coupled formulation used in this method is able to effectively update contact pre-load distribution in the non-linear analysis. Figure 21 shows the number of wear cycles ๐ ๐ at each wear iteration as wear progresses until the loss of contact. It is important to note that the number of cycles at each wear iteration is changing. It is adapted according to the geometry and loading conditions. It is worthy of mentioning that the shape of the curve shown in this figure is not unique. Rather it depends on the loading conditions and the pressure distribution at the contact as well as on the relative contact kinematics at the excited resonance.
1: slip-stick2: separation-slip-stick3: separation-slip
Wear iteration 1(start)Wear iteration 12Wear iteration 27Wear iteration 62(full loss of contact) Cumulative wear plot: Contact status plot: ๐ ๐ก๐๐ก๐๐ ๐น ๐๐ฅ = 1 v ๐ค,๐๐๐ฅ ๐ฟ ๐๐๐ฅ = Figure 20: Response graph showing the backbone of the non-linear response (EO = 1) obtained as wear progresses around first bending mode with free and stick states for reference for ๐ฏ ๐,๐๐๐ = 5% of ๐น ๐๐๐ Figure 21: Impact of wear on the number of vibration cycles at each wear iteration until full loss of contact for ๐ฏ ๐,๐๐๐ = 5% of ๐น ๐๐๐ [Note: the number of cycles is adaptive at each wear iteration] It is beneficial to track the maximum response with respect to the cumulative number of cycles as it provides a brief idea of the number of vibration cycles where the amplitude starts to rise considerably. Figure 22 shows the maximum response versus the cumulative number of cycles until the full loss of contact. Results show that for higher values of ๐ฏ ๐ค,๐๐๐ฅ , the rate of change of the response is under-estimated, leading to non-conservative predictions. Figure 22: Maximum response at each wear iteration vs the cumulative number of cycles for different values of ๐ฏ ๐,๐๐๐ Table ๐ฏ ๐ค,๐๐๐ฅ on the computational time, number of wear iterations necessary for the complete loss of contact, maximum wear depth and the error in the maximum wear depth. Table 4: Comparison of the impact of wear analysis with non-linear dynamic response solver for different values of ๐ฏ ๐,๐๐๐ ๐ฏ ๐ค,๐๐๐ฅ ๐ฟ ๐๐๐ฅ Computational time Total no. of wear iterations Max. total wear depth (mm) Error in wear depth (%) 0.5 % 5 hr 17 min 596 8.2566 - 1 % 2 hr 38 min 298 8.2572 0.01 2 % 1 hr 19 min 152 8.2916 0.4 5 % 0 hr 35 min 63 8.3808 1.5 10 % 0 hr 17 min 34 8.4697 2.6
The indicated computational times in Table 2 and Table 4 is the wall clock time obtained when the non-linear analysis was run on Intel Xeon 4 core processor @ 3.50 GHz and 32 GB RAM standalone workstation.
4. Conclusions
The current research work successfully studied the effect of wear on the dynamics of structures with frictional contacts using a multi-scale approach - the wear energy approach coupled with adaptive wear logic for accelerated wear and a coupled static/dynamic harmonic balance method formulation to capture the dynamics. The methodology is demonstrated on two test cases โ a cantilever beam with a contact patch and a tuned shrouded bladed disk with shroud contacts with cyclic symmetry properties. The proposed two-fold aim has been achieved: 1)
Modelling the effect of wear on the forced dynamic response and the interface evolution with โchanging contact pre-loadโ. This was achieved by loading the structure with a static pre-load distribution over the contact and automatically updated during the non-linear analysis as the contact loosens up due to wear. This method eliminates the need for separate non-linear static analysis routine. 2)
The effect of the choice of adaptive wear acceleration parameter ๐ฏ ๐ค,๐๐๐ฅ that is described in the paper. The effect of the choice of user-defined parameter ๐ฏ ๐ค,๐๐๐ฅ is studied on the dynamic response of the two test cases. With smaller values of ๐ฏ ๐ค,๐๐๐ฅ , the change in frequency and amplitude of response curves, wear evolution, and the contact profile is smoother. But with increased computation time and a higher number of wear iterations. In case of larger values of ๐ฏ ๐ค,๐๐๐ฅ , the results are coarse but with faster computation times. The choice of ๐ฏ ๐ค,๐๐๐ฅ is a trade-off between the computation time and the accuracy of results. Acknowledgements
This project has received funding from the European Unionโs Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 721865. References [1] R. B. Waterhouse, โFretting wear,โ
Wear , vol. 100, no. 1โ3, pp. 107โ118, 1984. [2] J. F. Archard, โContact and rubbing of flat surfaces,โ
J. Appl. Phys. , vol. 24, no. 8, pp. 981โ988, 1953. [3] G. K. Sfantos and M. H. Aliabadi, โApplication of BEM and optimization technique to wear problems,โ
Int. J. Solids Struct. , vol. 43, no. 11โ12, pp. 3626โ3642, 2006. [4] G. K. Sfantos and M. H. Aliabadi, โWear simulation using an incremental sliding Boundary Element Method,โ
Wear , vol. 260, no. 9โ10, pp. 1119โ1128, 2006. [5] P. Arnaud, S. Fouvry, and S. Garcin, โA numerical simulation of fretting wear profile taking account of the evolution of third body layer,โ
Wear , vol. 376โ377, pp. 1475โ1488, 2017. [6] S. Fouvry, P. Duรณ, and P. Perruchaut, โA quantitative approach of Ti-6Al-4V fretting damage: Friction, wear and crack nucleation,โ
Wear , vol. 257, no. 9โ10, pp. 916โ929, 2004. [7] S. Toumi, S. Fouvry, and M. Salvia, โPrediction of sliding speed and normal force effects on friction and wear rate evolution in a dry oscillating-fretting PTFE Ti-6Al-4V contact,โ
Wear , vol. 376โ377, pp. 1365โ1378, 2017. [8] A. I. Vakis et al. , โModeling and simulation in tribology across scales: An overview,โ
Tribol. Int. , vol. 125, pp. 169โ199, 2018. [9] L. Gallego, D. Nรฉlias, and C. Jacq, โA Comprehensive Method to Predict Wear and to Define the Optimum Geometry of Fretting Surfaces,โ
J. Tribol. , vol. 128, no. 3, p. 476, 2006. [10] A. Fantetti et al. , โThe Impact of Fretting Wear on Structural Dynamics: Experiment and Simulation,โ
Tribol. Int. , vol. 138, pp. 111โ124, 2019. [11] R. Lacayo et al. , โNonlinear modeling of structures with bolted joints: A comparison of two approaches based on a time-domain and frequency-domain solver,โ
Mech. Syst. Signal Process. , 2019. [12] N. M. Newmark, โNewmark A Method of Computation for Structural Dynamics,โ in
Journal of the Engineering Mechanics Division , 1959. [13] S. Zucca, โModelling Friction Contacts in Structural Dynamics and its Application to Turbine Bladed Disks.โ [14] S. Zucca, C. M. Firrone, and C. Gastaldi, โModels and Methods for the Dynamics of Mechanical Components With Contact Interfaces.โ [15] C. M. Firrone and S. Zucca, โEffect of static/dynamic coupling on the forced response of turbine bladed disks with underplatform dampers,โ
Power , pp. 1โ12, 2009. [16] L. Pesaresi, J. Armand, C. W. Schwingshackl, L. Salles, and C. Wong, โAn advanced underplatform damper modelling approach based on a microslip contact model,โ
J. Sound Vib. , vol. 436, pp. 327โ340, 2018. [17] C. Siewert, L. Panning, J. Wallaschek, and C. Richter, โMultiharmonic Forced Response Analysis of a Turbine Blading Coupled by Nonlinear Contact Forces,โ
J. Eng. Gas Turbines Power , vol. 132, no. 8, p. 082501, 2010. [18] C. Gastaldi, A. Fantetti, and T. Berruti, โForced Response Prediction of Turbine Blades with Flexible Dampers: The Impact of Engineering Modelling Choices,โ
Appl. Sci. , 2017. [19] S. Nacivet, C. Pierre, F. Thouverez, and L. Jezequel, โA dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems,โ
J. Sound Vib. , vol. 265, no. 1, pp. 201โ219, 2003. [20] M. R. W. Brake,
The Mechanics of Jointed Structures . 2018. [21] J. Armand, L. Pesaresi, L. Salles, and C. W. Schwingshackl, โA Multiscale Approach for Nonlinear Dynamic Response Predictions with Fretting Wear,โ
J. Eng. Gas Turbines Power , vol. 139, no. 2, pp. 1โ7, 2017. [22] J. Armand, โNonlinear dynamics of jointed structuresโฏ: a multi-scale approach to predict fretting wear and its effects on the dynamic response,โ 2017. [23] J. Armand, L. Salles, C. W. Schwingshackl, D. Sรผร, and K. Willner, โOn the effects of roughness on the nonlinear dynamics of a bolted joint: A multi-scale analysis,โ Eur. J. Mech. A/Solids , vol. 70, no. February, pp. 44โ57, 2018. [24] L. Salles, L. Blanc, A. Gouskov, P. Jean, and F. Thouverez, โDual Time Stepping Algorithms With the High Order Harmonic Balance Method for Contact Interfaces With Fretting-Wear,โ
J. Eng. Gas Turbines Power , vol. 134, no. 4, p. 032503, 2011. [25] L. Salles, A. M. Gouskov, L. Blanc, F. Thouverez, and P. Jean, โDynamic analysis of fretting-wear in joint interface by a multi-scale harmonic balance method coupled with explicit or implicit integration schemes,โ
ASME. Turbo Expo Power Land, Sea, Air, Vol. 6 Struct. Dyn. Parts A B , vol. 6, pp. 1003โ1013, 2010. [26] L. Salles, L. Blanc, F. Thouverez, A. M. Gouskov, and P. Jean, โDynamic Analysis of a Bladed Disk With Friction and Fretting-Wear in Blade Attachments,โ
Proc. ASME Turbo Expo 2009 Power Land, Sea Air , pp. 465โ476, 2009. [27] L. R. Tamatam, D. Botto, and S. Zucca, โEffect of wear on the dynamics of structures with contact interfaces by a coupled static/dynamic multi-harmonic balance method,โ in
First International Nonlinear Dynamics Conference. Book of abstracts , 2019, no. 1, pp. 129โ130. [28] L. R. Tamatam, D. Botto, and S. Zucca, โA coupled approach to model the effect of wear on the dynamics of bladed disks,โ in
Internation congress on sound and vibration (ICSV26) , 2019, pp. 1โ8. [29] C. M. Firrone, S. Zucca, and M. M. Gola, โThe effect of underplatform dampers on the forced response of bladed disks by a coupled static/dynamic harmonic balance method,โ
Int. J. Non. Linear. Mech. , vol. 46, no. 2, pp. 363โ375, 2011. [30] S. Zucca and C. M. Firrone, โNonlinear dynamics of mechanical systems with friction contacts: Coupled static and dynamic Multi-Harmonic Balance Method and multiple solutions,โ
J. Sound Vib. , vol. 333, no. 3, pp. 916โ926, 2014. [31] M. Z. Huq and J. P. Celis, โExpressing wear rate in sliding contacts based on dissipated energy,โ
Wear , 2002. [32] M. Krack and J. Gross,
Harmonic Balance for Nonlinear Vibration Problems . 2019. [33] M. C. C. BAMPTON and R. R. CRAIG, JR., โCoupling of substructures for dynamic analyses.,โ
AIAA J. , vol. 6, no. 7, pp. 1313โ1319, 1968. [34] E. P. Petrov, โA Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks,โ
J. Turbomach. , vol. 126, no. 1, p. 175, 2004. [35] S. Mehrdad Pourkiaee and S. Zucca, โA Reduced Order Model for Nonlinear Dynamics of Mistuned Bladed Disks with Shroud Friction Contacts,โ
J. Eng. Gas Turbines Power , vol. 141, no. 1, Jan. 2019. [36] G. Jenkins, โAnalysis of the stress-strain relationships in reactor grade graphite,โ
Br. J. Appl. Phys. , vol. 13, no. 1, pp. 30โ32, 1962. [37] W. D. Iwan, โA Distributed-Element Model for Hysteresis and Its Steady-State Dynamic Response,โ
J. Appl. Mech. , vol. 33, no. 4, pp. 893โ900, 1966. [38] K. Valanis,
Fundamental Consequences of a New Intrinsic Time Measure. Plasticity as a Limit of the Endochronic Theory , vol. 32. 1978. [39] J. H. Griffin, โFriction Damping of Resonant Stresses in Gas Turbine Engine Airfoils,โ
J. Eng. Power , vol. 102, no. 2, pp. 329โ333, 1980. [40] C.-H. Menq, J. Bielak, and J. H. Griffin, โThe influence of microslip on vibratory response, part I: A new microslip model,โ
J. Sound Vib. , vol. 107, no. 2, pp. 279โ293, Jun. 1986. [41] C. C. de Wit, H. Olsson, K. J. Astrom, and P. Lischinsky, โA new model for control of systems with friction,โ
IEEE Trans. Automat. Contr. , vol. 40, no. 3, pp. 419โ425, 1995. [42] B. D. Yang and C. H. Menq, โCharacterization of 3D Contact Kinematics and Prediction of Resonant Response of Structures Having 3D Frictional Constraint,โ
J. Sound Vib. , vol. 217, no. 5, pp. 909โ925, 1998. [43] D. J. Segalman, โA Four-Parameter Iwan Model for Lap-Type Joints,โ
J. Appl. Mech. , vol. 72, no. 5, pp. 752โ760, 2005. [44] M. Brake, โA reduced Iwan model that includes pinning for bolted joint mechanics,โ
Nonlinear Dyn. , vol. 87, no. 2, pp. 1335โ1349, 2017. [45] R. Bouc, โA Mathematical Model for Hysteresis,โ
Acustica , vol. 24, no. 1, pp. 16โ25, 1971. [46] Y.-K. Wen, โMethod for Random Vibration of Vibrating System,โ
Nanotechnology , vol. 8287, pp. 142โ145, 1976. [47] K. Y. Sanliturk and D. J. Ewins, โModelling Two-Dimensional Friction Contact and Its Application Using Harmonic Balance Method,โ
J. Sound Vib. , vol. 193, no. 2, pp. 511โ523, 1996. [48] E. P. Petrov and D. J. Ewins, โAnalytical Formulation of Friction Interface Elements for Analysis of Nonlinear Multi-Harmonic Vibrations of Bladed Disks,โ
J. Turbomach. , vol. 125, no. 2, p. 364, 2003. [49] T. M. Cameron and J. H. Griffin, โAn Alternating Frequency / Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems,โ
J. Appl. Mech. , vol. 56, p. 149, 1989. [50] O. Poudou and C. Pierre, โHybrid Frequency-Time Domain Methods for the Analysis of Complex Structural Systems with Dry Friction Damping,โ , no. April, pp. 1โ14, 2003. [51] A. Cardona, A. Lerusse, and M. Gรฉradin, โFast Fourier nonlinear vibration analysis,โ
Comput. Mech. , vol. 22, no. 2, pp. 128โ142, 1998. [52] M. Allara, โA model for the characterization of friction contacts in turbine blades,โ
J. Sound Vib. , vol. 320, no. 3, pp. 527โ544, 2009. [53] M. Ciavarella, D. A. Hills, and G. Monno, โThe influence of rounded edges on indentation by a flat punch,โ
Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. , vol. 212, no. 4, pp. 319โ327, 1998. [54] J. C. Jรคger, โA new principle in contact mechanics,โ
ASME J Tribol. , vol. 120, no. October 1998, pp. 677โ684, 1998. [55] K. L. Johnson,
Contact Mechanics . Cambridge: Cambridge University Press, 1985. [56] H. C. Meng and K. C. Ludema, โWear models and predictive equations: their form and content,โ pp. 443โ457, 1995. [57] S. Fouvry, T. Liskiewicz, P. Kapsa, S. Hannel, and E. Sauger, โAn energy description of wear mechanisms and its applications to oscillating sliding contacts,โ
Wear , 2003. [58] D. Botto and M. Lavella, โHigh temperature tribological study of cobalt-based coatings reinforced with different percentages of alumina,โ
Wear , vol. 318, no. 1โ2, pp. 89โ97, 2014. [59] M. Lavella and D. Botto, โFretting wear characterization by point contact of nickel superalloy interfaces,โ
Wear , vol. 271, no. 9โ10, pp. 1543โ1551, 2011. [60] D. Botto and M. Umer, โA novel test rig to investigate under-platform damper dynamics,โ