## Abstract

Turbine bladed disks or blisks, which constitute critical components of most modern turbomachinery, are known for their complex vibratory behavior. The nonlinear dynamics observed in most operational regimes of blisk with contact interfaces are dominated by one of two typical contact behaviors. Frictional contacts are dominated by Coulomb friction forces, while intermittent contacts are characterized by multiple separation events. Other factors such as the dispersion in material or geometric properties across blades, known as mistuning, also affect the dynamics significantly. Presently, probabilistic analysis is the widely accepted design methodology to account for mistuning, which is unknown prior to manufacture. Thus, reduced order modeling of these blisks is essential as high fidelity models are prohibitively expensive for such simulations. This paper provides a technical discussion of dynamic modeling and reviews projection-based techniques used for creation of reduced models of blisks with contacts.

## 1 Introduction

Bladed disks, also alternatively called blisks or rotors, which constitute the compressor and turbine stages of all turbomachinery experience high thermal and cyclic structural stresses during their operation, which involves extremely high rotational speeds. As a consequence, a significant thrust of research in the field of turbomachinery has been dedicated to understanding, quantifying, and modeling the dynamic behavior of bladed disks with the aim of accounting for their inherently harsh operating conditions and producing reliable and robust designs. Early studies into the subject during the 1960s and 1970s involved quantifying the (linear) mass, stiffness, and damping characteristics of the blisks employing lumped mass representations [1,2]. With the improvement in modeling techniques and computational capabilities, these simplistic models were replaced by high fidelity finite element (FE) models, which became the standard for industry by the mid 1990s- early 2000s [3]. Although a single linear analysis for even a large FE model, which often contains millions of degree-of-freedom (DoFs), can be performed within a reasonable period of time typically on the order of hours or days with current computational speeds, the multiple analyses that need to be carried out during the design procedure make a design approach based only on such high fidelity simulations cumbersome. As a result, the development of reduced order models (ROMs) of blisks, which allow faster simulations for accurate response predictions, emerged as an important research area in the field.

The most obvious reduction one may apply to any blisk model is based on its property of cyclic symmetry. If any sector of the blisk is considered to be geometrically and materially identical to any another, the DoFs of the blisk dynamic may be reduced by a factor of the number of blades by the application of suitable constraints. However, it is well known that even small differences, which exist between different sectors in practice due to manufacturing tolerances, collectively called mistuning, can have significant effects on blisk response amplitudes [414]. Mistuning may be small when the variability in parameters between the sectors does not significantly affect the modes of the structure, or large when the response of the mistuned blisk in a certain frequency range cannot be captured by a set of nominally cyclic symmetric or tuned modes whose natural frequencies lie in the corresponding range. Thus, the phenomenon of mistuning limits the applicability of a cyclic symmetric model to approximately determining the natural frequencies of a system with small mistuning. The response amplification, which occurs due to mistuning, also necessitates the use of probabilistic analyses [15,16] in the design of blisks. These analyses account for uncertainties associated with the values of mistuning parameters such as material stiffness and density, which cannot be exactly specified in the design as they are subject to manufacturing uncertainties.

Reduced order models are indispensable for Monte Carlo (MC) type probabilistic analyses, which require multiple runs of different mistuning patterns applied to the same underlying tuned system. The effects of mistuning were studied widely, and a number of different ROMs were proposed for linear mistuned blisks. These include modifications of general linear reduction methods such as Craig-Bampton component mode synthesis (CB-CMS) [1722], as well as more specialized techniques developed especially for blisks such as subset of normal modes (SNM) [23], fundamental model of mistuning (FMM) [24], component mode mistuning (CMM) [2527], and asymptotic model of mistuning [2830]. Alternatives for linear systems with parametric variations include parametric reduced order models [31,32].

However, not all blisks are designed to behave linearly. In fact, nonlinearities are often introduced into blisk design to lower resonant responses and stresses. Traditionally, friction interfaces on shrouds, under-platform dampers or ring dampers have been employed to provide additional damping. Nonlinearities can also occur unintentionally such as cracks or defects in blades. Modeling these nonlinear effects is paramount to achieve accurate estimates of responses and stresses in blisks. The earliest attempts at modeling nonlinear contacts in blisks involved approximating the nonlinear forces using lumped parameters [33,34] at every sector. The research focus then shifted to modeling complex localized behavior of frictional forces [35,36]. Iwan models, which were initially developed for modeling hysteretic behavior under constant unidirectional normal load, have been employed extensively to model certain types of interfaces [3740]. Later, models which predict the local frictional forces under different regimes of stick, slip, and separation under varying periodic normal loads were proposed [4144]. To model nonlinear Coulomb friction accurately, a number of these models must be used at every contact interface to calculate the local values of normal and tangential frictional forces from the local relative displacements. Moreover, to find steady-state solutions, the friction forces over the entire period must be calculated iteratively until convergence is achieved. These additional calculations greatly increase the computational effort required to predict forced responses for blisks with frictional interfaces. The computational cost associated also increases sharply with the number of nonlinear degree-of-freedom. Practically, simulation times for calculation of the steady-state forced response of a high fidelity blisk model with nonlinearities at single frequency using time marching or transient dynamic analysis are in the order of days or weeks, even with modern computational hardware. Hence, reduction in the dynamic simulation times is essential to make computational predictions feasible for nonlinear blisks.

Reductions in calculations were initially obtained in the time domain by using well-known general approaches to nonlinear differential equation solving such as the harmonic balance method (HBM) [4547]. HBM converts the nonlinear differential equations in the time domain associated with the blisk dynamics into algebraic equations in the frequency domain associated with periodic harmonics of the response. Reductions in calculations are obtained by retaining a finite number of harmonics for calculating the solution. The conversion to the frequency domain entails a complication due to the nonlinear forcing functions in the time domain equations, which appear as functions of displacements at the interfaces. In general, it is not possible to find a closed form expression for the equivalent harmonic nonlinear forcing functions, which must appear in the frequency domain HBM equations. To address this issue, an alternating frequency-time technique (AFT) [48,49] was proposed, to numerically obtain the harmonic nonlinear forces by calculating them in the time domain as functions of displacements and converting them back into the frequency domain by using Fourier and inverse Fourier transforms. Calculation of the forced frequency responses for nonlinear blisk using HBM and AFT has become standard practice in the field [5060].

Due to the local nature of the forces, which act only at contact interfaces, it is possible to obtain further reductions. Researchers have used linear reduction techniques to obtain spatial reductions at the noncontact linear DoFs only [46,52,61,62]. However, these techniques retain all the nonlinear DoFs from the original models and must perform associated HBM and AFT simulations, which contribute to the majority of the calculation costs. Another perspective on understanding the response of nonlinear systems is the formulation of nonlinear normal modes (NNMs) [6367], which are nonlinear synchronous motions of the system where the motion of a single DoF of a system can describe the motion along all other DoFs. Unlike linear modes of a system, these modes are energy dependent and change with the level of activation and consequently with response frequency [68,69]. The calculation of these modes is no easier than calculating the nonlinear response itself. In fact, HBM is commonly used to calculate these modes for large and complex systems such as nonlinear blisks [68,69]. Since the concepts of superposition and (linear) orthogonality do not generally apply for NNMs [68,70], model reduction by using a reduced order basis of NNMs is not an option. Instead, researchers have sought to find invariant manifold descriptions [66,7073] of these NNMs to reduce the nonlinear dynamics. In the past, such methods have been applied to simple models with relatively few DoFs [70]. Reduction of the dynamics of large models such as blisks with complex friction nonlinearities using NNMs is quite challenging due to high computational costs and remains a field of active research [74,75].

More commonly for such nonlinear systems, researchers have used other methods to calculate a linear basis, which approximately spans the nonlinear motion space and projected both the linear and nonlinear dynamics onto that basis to obtain ROMs. One method to obtain such a basis is to augment a basis obtained from a linear system with other basis vectors. These other vectors might be obtained using derivatives of the original basis vectors with respect to parameters [31,32,76]. Such system parameters may affect the system dynamics nonlinearly. Commonly, this nonlinear parameter is the input amplitude of the system, in which case the additional basis vectors perform a function similar to describing functions in control theory [77] where a nonlinear system is approximated by a linear one which changes with amplitude. A similar approach is an equivalent linearization [78] of the nonlinear system where nonlinear terms in the dynamics may be estimated by equivalent linear terms. Equivalent linearization type techniques may also be used to find additional vectors to augment projection bases [79]. Extending these ideas, many approaches for the formulation of projection-based ROMs discussed here involve the combination of basis vectors or modes of linear systems, which closely represent the physics of the nonlinear system in some specific/selected operational regime. Note that projection-based ROMs thus formulated may be numerically ill-conditioned due to linear dependency between vectors. Consequently, intelligent truncation and conditioning of the basis is often carried out using methods such as singular value decomposition [80,81].

Many popular ROMs use proper orthogonal modes (POMs) [8284] as a projection basis to estimate the span of the nonlinear response. The main drawback of such POM based methods and other, more advanced methods referred to as hyper-reduction [8589] is the computational power required to calculate the nonlinear responses from which the POMs are calculated. Even with the associated computational cost, POMs may be feasible for creating ROMs of specific, well-understood systems, which are not expected to change their dynamic behavior between instances and applications. However, they remain a weak choice for systems like blisks where fundamental uncertainties of mistuning and associated probabilistic calculations would require that nonlinear simulations be performed every time the model is changed slightly. There also exist other projection-based methods, which obtain reductions by capturing the physics of the dynamic problem at hand. While such ROMs may not be as broadly applicable to different structures and nonlinearities (as POMs are), they serve well for the specific application they are designed for.

One of the early methods pertinent to blisks with intermittent contact or cracks is based on ideas of bilinear modes (BLMs) representing the dynamics of localized piecewise linear systems [9092]. BLMs are linear normal modes for the system with special boundary conditions (BCs) at the surface where the intermittent contact takes place. In Ref. [91], it was shown that BLMs are able to capture the nonlinear dynamics of a cracked plate with intermittent contact by approximating the dominant POMs calculated from the nonlinear response. The concept of BLMs was extended in Ref. [93] to directly calculate amplitudes of the periodic nonlinear steady-state response at resonant frequencies using modes similar to BLMs. This method is referred to as bilinear amplitude approximation (BAA), and is considerably faster than the BLM ROMs used in Ref. [91]. BAA was used to construct ROMs for blisks in Ref. [94] also. More recently, in Ref. [92], BLMs were used to reduce a cracked plate model where the effects of frictional forces and contact prestresses were also considered. In Ref. [95], ROMs were developed using piecewise linear modes (PLMs), which are normal modes of the piecewise linear systems, which approximate the instantaneous structural dynamics of the nonlinear system.

Another important operational regime for nonlinear blisks, especially with friction damping mechanisms, is microslip [35,36,44,46,80,9699]. A structure in microslip is dominated by frictional effects and exhibits complex stick-slip behavior at contact interfaces both spatially and temporally. ROMs for blisks with ring dampers were presented in Ref. [62] where a reduction basis comprised of linear modes of the damper in full stick and gross slip were used. The same reduction basis was used in Ref. [79], which achieved further simplification of the dynamics by approximating equivalent single harmonic modal damping and stiffness parameters for the nonlinear forces using an energy equivalence calculation. The dynamics of a shrouded blisk is significantly different from the afore-mentioned structures with cracks or from blisks with ring dampers. Unlike in cracked structures, friction plays a more dominant role in shrouded blisks. Moreover, due to the contacts at the shrouds being in proximity of the blade tip, there is a significant effect of the frictional forces on the blisk response to such an extent, that the nonlinear displacements no longer lie on the subspace spanned by the linear modes corresponding to stick and slip. A ROM for such shrouded blisks was developed in Ref. [80] using an adaptive microslip projection (AMP) basis, which can span the required subspace by approximating the spatial correlations exhibited by the system. The AMPs are obtained from modes of intermediate linear systems corresponding to special BCs at the contact interfaces by using concepts similar to those used for PLMs. An evolution of the AMP method, developed especially for applications to nonlinear systems with dominant multiharmonic dynamics, called the Jacobian projection (JP) method, was presented recently also [100].

The focus of this paper is to review the techniques and the concepts of projection-based reduction methods for blisks with frictional or intermittent contacts. To facilitate the discussion of projection-based ROMs, a technical and theoretical background in the dynamics of blisks relevant to the topic is provided first. The nuances of the dynamics of nominally cyclic symmetric structures and the effects of mistuning are discussed in greater detail. A brief mathematical overview of HBM and AFT is used to impart a more concrete understanding of the nonlinear blisk dynamics and provide a rigorous framework for formulating the reduction problem using projection bases. This is followed by a discussion of contact nonlinearities and their behavior in different regimes of operation, supplemented by a brief description of popularly used contact models from literature and their applicability in capturing the physics of these nonlinear effects. Finally, an in-depth technical discussion and review of projection-based reduced order models for blisks are provided.

## 2 Nominal Cyclic Symmetry in Linear Blisk Dynamics

### 2.1 Cyclic Coordinates.

A cyclic symmetric blisk or tuned blisk is comprised of nmax materially and geometrically identical contiguous sectors $Sn (n=1,2,…,nmax)$. Each sector contains interior DoFs, which are not shared with other sectors, as well as high and low interface DoFs, which are shared with the adjacent sectors as shown in Fig. 1.

Fig. 1
Fig. 1
Close modal
The high and the low nodes are usually matched in a cyclic symmetric model such that one set may be obtained from the other by rotation about the axis of the blisk by an amount equal to the sector angle $β=2π/nmax$. The stiffness matrix for any sector Sn with free BCs at high and low nodes is in cylindrical coordinates for all n, and may be represented as
$KS=[KLLKLIKLHKILKIIKIHKHLKHIKHH]$
(1)

where subscript L refers to low DoFs, subscript I refers to interior DoFs, and subscript H refers to high DoFs. Since all the operations applied to the system matrices in this discussion are similar, only the operations applied to the stiffness matrix are shown. The matrices corresponding to mass and damping will have similar forms. The full blisk is formed by merging the interface DoFs of adjacent sectors as shown in Fig. 2. This is equivalent to the addition of stiffnesses at interfaces.

Fig. 2
Fig. 2
Close modal
The corresponding full system stiffness matrix K is given by
$K=[KIIKIH000⋯0KILKHIKHH+KLLKLIKLH0⋯0KHL0KILKIIKIH0⋯⋯00KHLKHIKHH+KLLKLI⋯0000⋱⋱⋱⋱⋱⋮⋮⋮⋱⋱⋱⋱⋱⋮00000⋯KIIKIHKLIKLH000⋯KHIKHH+KLL]$
(2)
where 0's represent a matrix of zero elements of the appropriate size, which depends on the number of nodes and the DoFs per node. It may be observed that K may also be written as
$K=[Kgen,1Kgen,2Kgen,3⋯Kgen,nmaxKgen,nmaxKgen,1Kgen,2⋯Kgen,nmax−1⋮⋮⋮⋱⋮Kgen,2Kgen,3Kgen,4⋯Kgen,1]$
(3)
$Kgen,n$ are called generating matrices of K and are given by
$Kgen,1=[KIIKIHKHIKHH+KLL];Kgen,2=[00KLIKLH]Kgen,j=[0000];(j=3,4…,nmax−1);Kgen,nmax=[0KIL0KHL]$
(4)
A matrix with this special structure is called a block circulant matrix. The symmetries associated with the geometric structure of the blisk, which are reflected in the mathematical interpretation of its block circulant system matrices, make it conducive to mathematical treatment using group theory. A comprehensive mathematical theory for circulant matrices was first developed in Ref. [101]. A recent summary of this theory as applied to blisk dynamics may be found in Ref. [102]. An especially convenient special property of these block circulant matrices is that they may be block diagonalized using the transformation matrix $(E⊗I)$ as follows:
$(E⊗I)HK(E⊗I)=[K̃10⋯00K̃2⋮⋮⋮⋮⋱⋮0⋯⋯K̃nmax]=BDiag(K̃p+1); p=0,1,…,nmax−1$
(5)
where BDiag(.) represents an operator, which constructs a block diagonal matrix from its arguments of smaller matrices, I is the identity matrix of the appropriate size, superscript H represents the complex conjugate transpose (or Hermitian) of the matrix, and the symbol ⊗ refers to the direct (or Kronecker) product of two matrices, $K̃p$ are the blisk stiffness matrices corresponding to the pth spatial harmonic in cyclic coordinates, and E is the complex Fourier matrix defined as
$E=1nmax[111⋯11eiβe2iβ⋯e(nmax−1)iβ1e2iβe4iβ⋯e2(nmax−1)iβ⋮⋮⋮⋱⋮1e(nmax−1)iβe2(nmax−1)iβ⋯e(nmax−1)(nmax−1)iβ]$
(6)
where i is the square root of −1. The transformation in Eq. (5) is referred to as the transformation to cyclic coordinates. The term $eiβ$ seen in Eq. (6) is the $nmaxth$ primitive root of unity and the complex Fourier matrix is composed of its integral powers. The effect of the transformation in Eq. (5) is a Fourier decomposition in the spatial domain for each cyclic symmetric nodal DoF across the sectors. Submatrices $K̃p$ along the diagonal in Eq. (5) are given by
$K̃p=∑n=1nmaxKgen,nei(n−1)(p−1)β=[KIIKIH+KILe−i(p−1)βKHI+KLIei(p−1)βKHH+KLL+KLHei(p−1)β+KHLe−i(p−1)β]$
(7)
A more physical understanding of what these cyclic matrices $K̃p$ represent may be obtained by studying an alternative derivation for them. One may apply constraints between the nodal DoFs of matching nodes on low and high surfaces of the free sector [53] described by
$q̃Lp=q̃Hpe−ipβ$
(8)
where q represents a displacement vector, and $̃$ represents its constrained nature. Eq. (8) also illustrates why is sometimes referred to as the interblade phase angle [19], as it is similar to imposing a phase difference between the complex displacements of the high and low nodal DoFs in a traveling wave. The transformation between the free and constrained nodal DoFs may be expressed by a matrix Tp as follows:
$[qLqIqH]=Tp[q̃Ipq̃Hp];Tp=[0Ie−ipβI00I]$
(9)
Then, the stiffness matrix in cyclic coordinates $K̃p$ of Eq. (7) may be recovered by applying the transformation matrix Tp in Eq. (9) to the free sector stiffness matrix KS in Eq. (1) as follows:
$K̃p=TpHKSTp$
(10)

### 2.2 Natural Frequencies and Modes of a Tuned Blisk.

The block diagonalization in Eq. (5) decouples the tuned full blisk system into smaller independent systems whose sizes are of the order of a single sector model and are obtained by applying constraints in Eq. (8) on the sector model with free high and low nodes. As the transformation in Eq. (5) is linear, the generalized eigenvalues of the stiffness and mass matrices in cyclic coordinates $K̃p$ and $M̃p$ are the same as the generalized eigenvalues of the full system matrices K and M. Note that these smaller systems are not approximations and together comprise an exact description of the full linear tuned blisk since the spatial Fourier basis used in the projections of these transformations spans the full space. Hence, the natural frequencies of the full system may be obtained by nmax smaller sector-level eigenvalue problems as opposed to solving it for the coupled matrices of the full blisk which is much more expensive computationally.

Due to these specific properties of a cyclic symmetric system, the natural frequencies and modes also exhibit certain unique characteristics. The modes of the full blisk may be obtained from the modes calculated from the cyclic matrices by using the transformation EI. As the modes are subject to constraints of Eq. (8), the corresponding full blisk expansions display regular patterns corresponding to the number of spatial harmonics corresponding to p. Each mode has a specific number of diameters along which the nodal displacements are zero known as nodal diameters (ND) of that mode. Sometimes, the number of NDs exhibited by a mode is also referred to as the harmonic index.

One may also observe that the powers of the $nmaxth$ root of unity also have the property $ei(nmax−p)β=e−ipβ$. An examination of Eqs. (7) and (8) in light of this property shows that the cyclic matrices $K̃p$ and $K̃(nmax−p)$ may be formulated by applying the same constraints except that the high and low nodes are reversed in definition. However, the eigenvalues yielded by both these matrices will be the same. Consequently, tuned blisks have repeated natural frequencies corresponding to all $K̃p$ except p =0 (and p = nmax/2 if nmax is even). The corresponding modes for repeated eigenvalues contain the same number of NDs and have similar spatial patterns, differing only in circumferential position. These mode pairs are also called forward and backward traveling modes [23,24,102]. The backward traveling modes correspond to values of p higher than the maximum number of NDs (nmax/2 or $(nmax−1)/2$ for blisks with even and odd nmax, respectively). Hence, as per the theory of Fourier transforms, an alternative but consistent nomenclature may be used, where these backward traveling modes may be ascribed negative spatial harmonic number based on the equivalent negative phase value yielded by the powers of $nmaxth$ root of unity higher than the maximum number of NDs.

The natural frequencies for tuned blisks are often arranged in plots such as the one shown in Fig. 3 where the natural frequencies are plotted against the number of NDs exhibited by the corresponding modes. Figure 3 shows such a plot for a blisk with 24 sectors, which has modes with up to a maximum of 12 NDs. The natural frequencies are represented by circle markers. There are two modes corresponding to every natural frequency except for modes with ND 0 or 12. The horizontal lines connecting these frequencies are known as modal families with the family of lowest frequencies being the first family, the next one the second family, and so on. Figures 4 and 5 show typical mode pairs corresponding to ND 1 and ND 8, respectively. Within a family, there exist some flat or nearly horizontal regions of high modal density where the natural frequencies are very close to each other. For conventional blisk designs, these usually occur for modes with higher NDs, which are dominated by displacements at the blade, such as the ND 8 pair in Fig. 5. Hence, these are known as blade-dominant modes. The modes in these flat regions of each family may thus be approximated by applying the cyclic transformation to a single mode of the cantilevered blade. Modes with lower number of NDs such as the ND 1 pair in Fig. 4, where the families are not flat, have significant displacement at the disk DoFs. The regions in the frequency versus nodal diameter plot where two families veer away from each other, such as the one shown in Fig. 3, are characterized by disk-blade interactions as well as contributions from multiple blade modes.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal

### 2.3 Linear Dynamics of a Tuned Blisk.

The NDs of the modes are important also in the context of the typical excitation patterns encountered by blisks in operation. The internal structures of axial compressors or turbines comprise of alternating stators and rotors. Since the fluid encounters a stator prior to the rotating blisk, the pressure wave encountered by the blades of the blisk has a dominant spatial harmonic pattern related to the number of blades in the stator. While studying the dynamics of the blisk computationally, an equivalent dynamic problem is considered for convenience. The problem is analyzed from the perspective of the rotating frame of reference, where the blisk is stationary and the pressure wave is rotating [4,103,104]. Such a rotating wavefront, called a traveling wave excitation, is conceptually illustrated in Fig. 6. When a traveling wave has a single spatial harmonic component, it is referred to as an engine order (EO) type excitation where the EO is equal to the dominant spatial harmonic of the traveling wave. To capture and study the effects of such excitation in simulation, it is often sufficient to apply discrete time varying forces at certain nodal locations (commonly the blade tips) as shown in Fig. 6. This excitation fE at sector n at time t may be expressed as
$fE(t,n)=|F| cos(ωt+2πEO(n−1)nmax)$
(11)

where $|F|$ is the excitation force amplitude. Note that ω = ΩEO, where Ω is the rotational speed of the blisk.

Fig. 6
Fig. 6
Close modal
Designers are interested in the studying the dynamics of the blisk for long-term steady-state operation which leads to high cycle fatigue. Hence, the steady-state frequency responses of blisk are calculated using the dynamic equation of motion (EOM) in the frequency domain, which is expressed as
$Dq¯=f¯E ; D=[−ω2M+iωC+K]$
(12)
where $¯$ represents complex quantities in the frequency domain, ω is the frequency of excitation, D is the dynamic stiffness matrix, and C is the linear damping matrix commonly obtained by assuming proportional or structural damping. In the frequency domain, the excitation vector $f¯E$ is a spatial harmonic pattern, and is populated at the excitation DoFs by complex numbers with amplitudes $|F|$. Applying the cyclic transformation, one may decompose the system into smaller systems as before and the excitation vector is transformed as follows:
$[−ω2M̃p+iωC̃p+K̃p]q¯̃p=f¯̃Epq¯̃p=(ep⊗I)q¯ ; f¯̃Ep=(ep⊗I)Hf¯E$
(13)

where ep is the pth column of the complex Fourier matrix E. For EO type excitation, only the spatial excitation component $f¯̃Ep$ corresponding to the number of harmonics in the traveling wave will be nonzero. Hence, only the reduced cyclic subsystem corresponding to those harmonics needs to be simulated to obtain the response. In general, any mode-pair of the blisk may only be excited by a forcing with the same EO as its ND and by higher EO = jnmax ± ND where j is any integer. Hence, not only are significant computational savings obtained in the simulation of linear tuned blisks by using the cyclic reduction, but in most cases, the dynamics may be captured by simulating only one or two of these subsystems corresponding to the expected EO excitation.

## 3 Mistuning and Its Consequences

### 3.1 Mistuned Modes and Responses.

Mistuning occurs when there is a variability across the sectors of the blisk in a parameter, which affects its dynamics. Variations in different material properties such as mass, stiffness, and damping as well as in geometry across sectors have been studied in literature [4,5,814,27,105107]. The consequence is that the sector-level system matrices are nonidentical across sectors. The mistuned stiffness matrices for the free sector and full blisk may be represented as
$KnS=KoS+ΔKnS ; K=Ko+ΔK ; ΔK=BDiag(ΔKnS)$
(14)
where the subscript o represents a tuned matrix or mode, $ΔKnS$ is the mistuning stiffness component at sector n, and ΔK is the change in stiffness due to mistuning for the entire blisk. The transformation to cyclic coordinates does not yield an exact block-diagonal structure as in Eq. (5). Instead, one obtains
$(E⊗I)HK(E⊗I)=BDiag(K̃op)+(E⊗I)HΔK(E⊗I)$
(15)

The projection of mistuning onto cyclic coordinates is a full matrix in general. Hence, the mistuned systems may not be decomposed perfectly. The modes of a mistuned blisk also do not exhibit perfect ND patterns in general, and multiple modes will respond to an EO excitation even in a relatively narrow frequency range spanning an isolated mode family. Figure 7 shows such typical frequency responses of a tuned blisk and a mistuned blisk with small stiffness mistuning. The maximum normalized amplitude response across all the blades is plotted. In the tuned blisk response, only modes with the corresponding ND respond to the EO excitation. The response amplitudes of all the blades are identical, differing only in phase. Hence, the maximum response amplitude across all blades is the same as the response of any given blade. The frequency of the single peak is close to the natural frequency of the responding mode-pair. In the mistuned blisk response, multiple modes respond, each with a different natural frequency. The modal displacements at each blade are also nonidentical for these mistuned modes. Hence, different blades respond differently as per the mistuned modal contributions and multiple peak frequencies occur corresponding to when each blade responds the most.

Fig. 7
Fig. 7
Close modal

In most cases, the maximum response amplitude over all blades across frequencies also increases when mistuning is added. The ratio of this maximum response of a mistuned blisk to that of a tuned blisk for identical excitations is called the amplification factor (AF). This amplification in response introduces an element of uncertainty in the analysis because the level (the variance across blades or sectors), pattern and source (material or geometric) of mistuning are unknown at the design stage. The amplification may be understood as a strain energy localization phenomenon [1,4], and is dependent on the energy transfer between sectors due to interconnectivity between the blades through the disk. It is due to this reason that the veering region dominated by modes with blade-disk interaction is most sensitive to mistuning. Note that sensitivity does not necessarily imply amplification of the response when a specific mistuning pattern is present, but there is a high likelihood of encountering some amplification. Different blade stiffness patterns with the same variance may exhibit drastically different AFs. In fact, sometimes, intentional mistuning patterns are introduced into design to ensure small AFs [3,14,106,108,109]. Several studies have been conducted to calculate the theoretical upper bounds for AF when mistuning is introduced in a linear blisk [9,28,110]. However, these limits are conservative, do not convey information regarding the specific patterns, which lead to the highest AFs, and are of limited use for operational designs. Currently, the industry standard when designing for mistuning is to investigate thousands of randomly generated patterns applied to a specific nominal tuned design, and conduct probabilistic analyses at different levels of mistuning [4,15,16,111,112].

However, any prediction of the AF requires solving the mistuned blisk dynamics. Industrial FE models, which typically contain millions of DoFs per sector, are cumbersome for such analyses. Hence, the use of ROMs becomes essential. Mistuned ROMs may also be useful in other scenarios where repeated simulations of blisks are necessary, such as mistuning identification from blisk response measurements [27,113116], optimizing intentional mistuning patterns in blades [3,14,109], deciding mistuned damper arrangements [117,118], or detecting abnormalities in the structure such as cracks [119,120].

### 3.2 Dynamic Reduction and Probabilistic Analyses.

Several ROMs of gradually increasing sophistication have been developed for predicting linear dynamics of ROMs with small mistuning in the blades. One of the earliest methods is based on the observation that the mistuned modes in a given frequency range may be a linear combination of the tuned modes of the families in that region $Φo$. Hence, the dynamics of any mistuned blisk may be projected onto this subset of nominal modes (SNM) to obtain the model reduction [23]. The projection basis is formulated by collating tuned modes, which can be calculated in cyclic coordinates and expanded based on their corresponding spatial harmonic p using the cyclic transformation as follows:
$Φo=[⋯(ep⊗I)ϕ̃op⋯]$
(16)
where $ϕ̃op$ is a tuned mode in cyclic coordinates. It then becomes possible to project the additional mistuned matrices onto the modal coordinates. For instance, the stiffness matrix undergoes the transformation as follows:
$ΦoHKΦo=ΦoH(Ko+ΔK)Φo=Λo+AΛo=ΦoHKoΦo; A=ΦoHΔKΦo$
(17)

where $Λo$ is a diagonal vector of eigenvalues corresponding to the selected tuned modes.

The FMM, developed in Ref. [24], is an extension of the SNM projection shown in Eq. (17). In the FMM procedure, the choice of nominal modes is limited to a single family. The projected matrix of mistuned parameters A is not diagonal in general. However, its values may be determined using sector level calculation exploiting the tuned cyclic modes in Eq. (16) and sector level mistuning components shown in Eq. (14). FMM further uses the approximation that the tuned modes in flat blade-dominant regions may be generated from cantilevered blade modes. FMM calculates the projected mistuning components as functions of blade alone frequencies and the average of system frequencies in the flat region for fast ROM formulation directly in the reduced space. However, FMM ROMs are only accurate for analysis near flat isolated families. The asymptotic model of mistuning [2830] generalizes the concept and extends the applicability FMM by relaxing assumptions regarding dominance of blade motion in the tuned system modes, seeking to instead approximate which tuned modes actively contribute the most to the mistuning components assuming small mistuning and damping in the modes.

Other studies have used variants of the CB-CMS [17,18] substructuring approach to generate mistuned blisk ROMs [16,19,20,22,25,121123]. In CB-CMS, a system is first partitioned into substructures, which are individually reduced and then coupled to each other to obtain a reduced formulation of the system. The whole system may be treated as a single structure in which case the second coupling step is not necessary. For each substructure, the DoFs of the system and consequently the system matrices may be partitioned as follows:
$q=[qMqS]; K=[KMMKMSKSMKSS]$
(18)
where subscript M refers to the master DoFs; subscript S refers to the slave DoFs; and subscripts MM, MS, SM, or SS represent a matrix partition based on the master–slave partition of the displacement vector. In CB-CMS reduction [17], the master nodes are retained unchanged while applying a reduction to the slave DoFs to reduce them to some generalized DoFs $ηS$ corresponding to a set of slave modes. The reduction is represented as follows:
$qCB=[qMηS]=[I0ΨCBΦCB][qMqS]$
(19)
where the subscript CB represents quantities related to CB-CMS reduction. qCB are the set of reduced CB-CMS coordinates, qM are displacements of master DoFs retained in the CB-CMS reduction, $ηS$ are reduced modal coordinates corresponding to a set of slave modes. $ΨCB$ is a matrix whose columns are referred to as constraint or attachment modes and $ΦCB$ is a matrix whose columns are called the normal or slave modes of the system. The constraint modes are obtained by applying a unit displacement to the master DoFs individually. Hence, they must satisfy the relation
$[KMMKMSKSMKSS][IΨCB]=[fM0]$
(20)
where fM is the force applied at the master DoFs to enforce the constraints. Simplifying the second row of the matrix Eq. (20), it follows that the constraint modes may be determined using the equation:
$ΨCB=−KSS−1KSM$
(21)
The normal modes are the linear modes of the system obtained by completely constraining all the master DoFs, and can be determined as the solution to the eigenvalue problem
$KSSΦCB=MSSΦCBΛSS$
(22)

where MSS represents the mass matrix, and $ΛSS$ represents a diagonal matrix containing eigenvalues of the constrained system. Further details regarding coupling of sub-structures may be found in Ref. [17]. Representations of the CB-CMS reduction modes and reduced matrices may also be obtained directly in cyclic coordinates [19,22,112,124] offering computational savings.

Craig-Bampton component mode synthesis, which uses constrained interface modes, is only one of the popular methods among the many different flavors of available substructuring methods [125127]. Other methods use free interface modes or a combination of both free and constrained interface modes [124].

One of the more general techniques developed for mistuned blisk reduction based on CB-CMS is known as CMM [2527]. CMM employs a substructuring approach where the change in the tuned system due to mistuning is treated as an additional component to the tuned system. The dynamics of the system are then formulated using DoFs corresponding to a subset of normal modes of the tuned system and constraint modes of the mistuning component. In its most general formulation, CMM can be applied to both small and large mistuning. However, since it is known from SNM studies that small mistuning at the blades may be captured using tuned normal modes only, a more specialized formulation of CMM ignores the additional reduced DoFs corresponding to the attachment modes for such cases. The second major assumption in CMM is that small mistuning at the blades may be captured by projection onto a set of nominal cantilevered blade modes. The CMM ROM is formulated by projecting the blade mistuning onto tuned cantilevered blade modes and further projecting onto the set of system modes (which may or may not include attachment modes) by employing modal participation factors, which can be obtained using sector level calculations and cyclic expansions. Mistuning parameters are expressed directly in the reduced space, and the participation factors need only be calculated once for a given nominal system. This allows for fast generation of CMM ROMs for new mistuning patterns without any involved calculations with larger models.

Another alternative perspective to ROM generation for linear mistuned blisk was provided in Ref. [128], which focuses on the calculation of the inverse of the mistuned dynamic stiffness matrix (D in Eq. (12)) using a series expansion representation. Specific methods also exist for creating ROMs for large mistuning effects where the mistuned modes may no longer be captured by the tuned modes [5,10,34,129132]. These cases usually arise due to large differences in geometry across structures for example due to missing material. Such ROMs require the development of specific modal bases, which can capture the dynamics.

Any of these ROMs may be used to carry out probabilistic analyses by generating mistuned models randomly and using MC simulations. The AFs for each randomly generated pattern may be obtained statistically. Worst case AFs are especially of interest and usually a high percentile such as the 95th percentile in the distribution of AF values for any given level of mistuning is used as a benchmark for design. MC simulations can only provide estimates of the percentile values of the actual distribution. As AFs are obtained for more cases with randomly generated mistuning patterns, the distribution of simulated AFs approaches the actual distribution, and the accuracy of the percentile estimates calculated from them increases. However, for reasonable estimates of high percentiles such as the 95th percentile, simulation of many cases (about 1000) may be required. It was shown in Ref. [19] that since the forced vibration response of mistuned linear blisks is bounded [9,28,110], the distribution of the AFs approaches a 3-parameter Weibull distribution. Hence, another estimate of the distribution may be obtained by fitting AFs from a relatively few cases (around 50–100) to a Weibull distribution. Figure 8 shows a comparison between the cumulative distribution function (CDF) of the MC distribution obtained from 1000 cases and Weibull fits to different sets of 50 cases. It is seen that these estimated CDFs are quite accurate and can be used to generate estimates of the percentile values for different levels of mistuning requiring fewer simulations for each mistuning level. Figure 9 shows such an analysis, where the circle markers represent the AFs for each individual simulation of response calculation of a blisk with randomly generated stiffness mistuning. For many blisk designs, the AF percentiles first increase and then decrease as the standard deviation (mistuning level) increases, giving rise to a local maxima [19]. This plot allows the determination of worst-case mistuning standard deviations where the blisk responses are likely to show highest amplification over the nominal case, and stresses in the structure are likely to be highest. This information may then aid engineers to design with adequate safety margins to ensure operational robustness.

Fig. 8
Fig. 8
Close modal
Fig. 9
Fig. 9
Close modal

## 4 Simulation Methods for Blisks With Nonlinearities

### 4.1 Contact Modeling.

Contact is a very complex phenomenon, and its modeling has developed into an entire branch of mechanics. Interaction between contacting surfaces may even require scrutiny at nanoscale resolution [133]. These interactions are highly localized and depend on a number of factors such as the size and shape of local surface features and asperities, local material properties, and ambient conditions [134]. Contacts in FE models are usually modeled as time-varying constraints, which may be imposed using different techniques. The penalty method uses an explicit penalty stiffness to penalize penetration between the contacting surfaces [135]. Another approach is the Lagrange multiplier method, where terms containing extra DoFs representing the contact forces are added to the dynamic equations and simultaneous solution of the dynamics and constraints ensures near-zero penetration [136,137]. An augmented Lagrange method contains concepts from both these methods and uses a penalty stiffness as well as a Lagrange multiplier [138,139]. For contact modeling in blisks, it is desirable to strike a balance between accurately capturing dominant effects of contact interactions, which are pertinent to blisk dynamics and maintaining feasible simulation times. The modeling method currently preferred for nonlinear blisk simulations [11,46,50,52,61,80,124,140148] is a penalty-based method, which involves using arrays of local contacts on every surface and approximating the contact effects on small regions around each location. For convenience, the local contact may be described between nodes at the contact surfaces of FE models, as it is easier to obtain the responses of these nodes during calculations.

Two types of such node-to-node contact models, which have been commonly used in blisk simulations, are shown in Fig. 10. These contact models calculate the normal and tangential contact forces as a function of relative displacements between the nodes at any instance of time.

Fig. 10
Fig. 10
Close modal
The two-dimensional (2D) contact model [41,149152] is so named as it calculates the force vector in the tangential direction as components along two orthogonal directions. Sometimes, this may also be referred to as a three-dimensional (3D) model in the literature referring to the consideration of all three orthonormal direction [41,146]. While this is the more accurate model [46,146], many studies use two separate one-dimensional (1D) contact models in orthogonal directions to independently estimate the tangential forces in those directions [42,43,50,80,140,153]. It has been shown that this strategy can obtain accurate estimates of the contact forces (while offering significant computational time savings) if the relative motion at the interface is aligned with the tangential direction of one of the independent 1D models during most of the cycle [46,146]. For both models, penetration is allowed, and the normal contact force is proportional to the magnitude of penetration when there is contact (i.e., in the stick or slip conditions) between the surfaces. The normal force is zero when there is separation. The normal force Ni at the ith at node pair, along the local normal direction zl as shown in Fig. 10, may be obtained as follows:
$Ni=max(kz,ivi,0)$
(23)

where kz,i is the normal contact stiffness, and vi is the normal relative displacement between the nodes, assumed negative when there is no penetration between the two surfaces. Ni may vary with time, and includes any prestresses due to constant normal loads. Prestress loads lead to a constant component of $vi$. An initial gap may also be modeled where the initial value of vi is negative.

The tangential contact forces $Tx,i$ at any point in time may be calculated based on the current contact state at the node pair, which can be described by one of three distinct cases: stick, slip or separation. For the 1D model, these relationships in the tangential direction $xl$ in Fig. 10 are given by
$Tx,i={kx,i(ux,i−wx,i)sticksign(w˙x,i)μNislip0separation$
(24)

where μ is the coefficient of friction, $kx,i$ is the appropriate tangential contact stiffness, and $ux,i$ is the tangential relative displacement between the two contact nodes. $wx,i$ is known as the slip in the tangential direction and describes the position of the free end of the tangential contact stiffness relative to one of the nodes, as shown in Fig. 10. The model is said to be in stick when $wx,i$ remains constant and the spring is allowed to change length, and in slip otherwise. $wx,i$ is an internal variable whose value is unknown at the beginning of simulation. The values of $wx,i$ and consequently all tangential contact forces are usually determined by running the model multiple times and obtaining convergence of the values over a predetermined time period. When not in separation, the state of the model is determined by the rate of change of $wx,i$. The Coulomb force friction limit $μNi$ may not be exceeded by the tangential force $Tx,i$ for any tangential spring deformation, thereby leading to motion of the free end of the tangential stiffness and slip. During slip, $wx,i$ changes so that the tangential force limit is maintained. Similar equations for the independent 1D model along the orthogonal tangential direction $yl$ in Fig. 10 may be obtained by replacing subscripts x with y.

The 2D model calculations differ from Eq. (24) only during the slip state [146]. The difference between 1D and 2D models during the stick state arises from the manner in which the Coulomb friction limit of $μNi$ is enforced on the tangential force. In the 2D model, the limit applies to the vector sum of the tangential forces $Ti=Tx,i2+Ty,i2$, which is constrained to be parallel to the tangential slip velocity vector $w˙i=w˙x,i2+w˙y,i2$. For 1D models, the same force limit is applied independently to the tangential forces along two orthogonal directions, which results in good accuracy when the tangential motion is aligned with either one of the 1D models, and low accuracy otherwise [146,154]. The tangential force components in the 2D model may be given by
$Tx,i={kx,i(ux,i−wx,i)stickw˙x,iw˙x,i2+w˙y,i2μNislip0separation$
(25)

$Ty,i$ may be calculated similarly. Although subtle, the differences in between Eqs. (24) and (25) have major implications in terms of computational time, especially when applied to blisk models [146]. Instead of the independent convergence of a single slip variable, the 2D model requires simultaneous convergence of two dependent slip variables under the constraint established by the tangential force direction. There is only one contact state shared across the two orthogonal directions at any given point in time for the 2D model, whereas they may be different in the 1D model. The transition times between contact states may also be different for the two models depending on the dynamics being analyzed.

Different aspects of contact modeling are emphasized for different dynamic behaviors. For intermittent contacts, the small magnitudes of tangential forces, which are encountered during motion, are often completely ignored [91,93,124]. In most operational regimes, contact surfaces in blisks are designed with high preloads, and the system generally does not encounter separation at any of the node pairs [46,80]. Both the 1D and 2D models represent approximations of reality, and as such both experimental methods [155157] and modeling techniques [158] have been developed to fit these models for contact behavior between surfaces with known materials, sizes, and finishes by estimating contact parameters (contact stiffnesses and coefficient of friction).

### 4.2 Relative Coordinates.

For a blisk with contact nonlinearities, the vector of displacements in absolute (global) physical coordinates is given by
$q=[qS1,aqS2,aqO]$
(26)

where subscripts S1 and S2 correspond to the two contacting surfaces such as the ones shown in Fig. 11. Subscript a refers to the representation of a quantity in physical absolute coordinates (such as the Cartesian or cylindrical coordinate system of the FE model). The nonlinear contact forces will be nonzero along DoFs $qS1,a$ and $qS2,a$, while they will have zero values along the other DoFs corresponding to $qO$. The DoFs along which the nonlinear forces act are known as nonlinear DoFs and the computational cost of solving the dynamics of the blisk increases greatly as their number grows.

Fig. 11
Fig. 11
Close modal
A critical observation, which may be made here, is that the forces at the contact are all functions of relative displacement at the contact [159,160] and not the absolute displacements of the contacting DoFs. Hence, one may transform the vector of absolute global displacements at contacts in q in Eq. (26) to the local relative coordinates [129,159] such that
$q=[qS1,aqS2,aqO]=[0R0RR000I]qr; qr=[qS2,l−qS1,lqS1,lqO]$
(27)
where l refers to the local coordinates corresponding to each contact interface as shown in Fig. 11, and R is the rotation matrix required for transformation from local to global coordinates. In these relative coordinates, the nonlinear forces will act only along the DoFs $qN=qS2,l−qS1,l$. Hence, the number of nonlinear DoFs is halved. The rest of the DoFs $qL=[qS1,l qO]T$ are the linear DoFs, where $T$ represents a matrix or vector transpose. The transformation may be applied to corresponding physical system matrices. The dynamic EOM for a blisk with contact in the time domain in relative coordinates may be expressed as
$Mfq¨r+Cfq˙r+Kfqr=fE+fC(qr)$
(28)

where $Mf,Cf,Kf$ are the mass, damping, and stiffness matrices of the (free) system where no contact conditions have been enforced, expressed in relative coordinates. fE is the vector of excitation forces usually an EO excitation. fC is the vector of contact forces.

### 4.3 Harmonic Balance.

Equation (28) may be simulated in the time domain using time marching techniques [147,161,162]. However, such calculations require multiple expensive simulations at different excitation frequencies. When steady-state solutions are desired, much of the computational effort is wasted is simulating transients due to initial conditions in such time marching procedures. Consequently, the problem is recast by making the assumption that the desired steady-state solution is periodic and may be expressed as a sum of finite number of temporal harmonics h as follows:
$qr=∑h=0hmaxℜ(q¯rheihωt)fE=∑h=0hmaxℜ(f¯Eheihωt)fC=∑h=0hmaxℜ(f¯Cheihωt)$
(29)
where $ℜ(.)$ represents the real part of a complex quantity, and hmax is the highest harmonic number retained in the approximation. Substituting Eq. (29) into Eq. (28) and equating the coefficients of the linearly independent harmonic functions $eihωt$, one may obtain
$Dhq¯rh=f¯Eh+f¯ChDh=[−(hω)2Mf+(ihω)Cf+Kf]$
(30)

This procedure is called the HBM [46,50,52] and yields a set of coupled nonlinear algebraic equations, which describe the nonlinear blisk dynamics in the frequency domain. The coupling among harmonics occurs due to the nonlinear contact force vector $f¯Ch$, which is a function of the displacements in time, and hence contributes to multiple harmonics in the nonlinear EOM. The excitation usually only has nonzero components corresponding to prestress $f¯E0$, which does not vary with time, and a periodic excitation $f¯E1$ with a fixed frequency ω. When the excitation is applied at a higher harmonic h (instead of h =1), the HBM-based solution can capture integral subharmonic responses.

The nonlinear algebraic equations in Eq. (30) may be solved using optimization algorithms, which minimize its residual. Most commonly, iterative gradient-based optimization techniques such as trust-region based algorithms are employed. Although the gradients of the residual with respect to the harmonics of the displacements may be approximated numerically using finite differences, such calculations increase the solution times by orders of magnitude and are infeasible for solutions of large systems. Hence, it is often necessary to provide the analytical gradients of the residuals to the algorithm [46,50,163]. The contact forces fC are described as functions of displacements in time by most contact models. Since it is impossible to predict the switching times between multiple contact states in most scenarios, obtaining a general analytical expression for the nonlinear contact force harmonics $f¯Ch$ as a function of the harmonic displacements $q¯rh$ is also not possible. Hence, most simulations employ an AFT or hybrid frequency time procedure [48,49] as shown in Fig. 12. During each iteration of the solver, the relative local displacements in the frequency domain are extracted from the estimated displacement vector $q¯rh$ at that step. These are then converted into local displacements in the time domain using an inverse fast Fourier transform (IFFT), and they are used to calculate the local temporal contact forces. A fast Fourier transform (FFT) is then used to recover the Fourier coefficients of these forces, which are then used to construct the nonlinear force vector in the frequency domain $f¯Ch$ for various harmonics. A similar procedure is also used to obtain the derivatives of the Fourier coefficients of the contact forces with respect to the Fourier coefficients of the displacements at each iteration step, which are then used to formulate the gradients of the residuals analytically.

Fig. 12
Fig. 12
Close modal

In some practical applications, continuation methods [46] are employed, where the dynamic EOMs are augmented with the solution frequency ω. A predictor-corrector method is used to solve the augmented system while following the response frequency curve. In this case, similar gradient-based optimization is still required for the corrector step. Continuation methods possess the advantage of having variable frequency steps, which do not need to be prespecified by the user, possessing an in-built logic for nonconvergence situations and not requiring multiple restarts when the system dynamics has multiple solutions at the same frequency.

### 4.4 Nonlinear Responses.

Figure 13 shows different types of contact interfaces in a blisk, which may lead to nonlinear behavior. While fir-tree joints (which join the blades to the disk) with improper fits or cracked blades [119,120,124,159,164] result in intermittent contacts, frictional damping mechanisms such as shrouds [11,50,74,80,143,146,165,166], under-platform [96,141,154,167170], or ring dampers [5,62,68,79,145] operate in the microslip regime (spatially and temporally varying stick and slip conditions in a contact region).

Fig. 13
Fig. 13
Close modal

Intermittent contact and microslip result in markedly distinct dynamic behavior of the blisk. Their effects are most dominant close to different frequency ranges corresponding to different linear conditions. Assuming that all contact node pairs at all interfaces have the same contact condition at all times in the period, three distinct linear blisk conditions may be obtained. One of these is full stick, which corresponds to all the contact node pairs being in the stick condition. Full stick is therefore equivalent to enforcing linear stiffnesses in the local normal and tangential directions at all contact node pairs. Another condition, called gross slip, which is characterized by a constant friction force, is modeled by enforcing linear contact stiffnesses only in the local normal directions and leaving the tangential direction free. The third condition is that of the free linear blisk with no contacts being enforced, and the interfaces being free to separate from or penetrate into each other. In a frequency range dominated by an isolated mode, the frequency responses of a tuned blisk in these three distinct conditions resemble the schematic in Fig. 14.

Fig. 14
Fig. 14
Close modal

Intermittent contacts represent nonlinear conditions, usually encountered during low preloading or high excitations, where local interface regions come in and out of contact and hence exhibit responses in frequencies between those corresponding to the free and gross slip conditions. Typical frequency responses for a dynamical system with intermittent contacts are shown in Fig. 15. At excitation frequencies far from the peak linear frequencies, the nonlinearities are not activated significantly, and the contacts are either free or in gross slip throughout the cycle. Consequently, the structure exhibits a near linear response in these ranges. Near (linear) resonant frequencies, the contact exhibits complex behavior with localized regions of alternating slip and separation conditions, and the frequency response may deviate significantly from the linear behavior. For large excitation amplitudes, multiple response amplitudes may also occur for the same excitation conditions corresponding to multiple equilibria of the nonlinear dynamics [50,142,171], indicated by bends in the frequency response, as shown in Fig. 15.

Fig. 15
Fig. 15
Close modal

Frequency responses for a structure with microslip contacts are shown in Fig. 16. Most commonly, these contacts are designed to operate under normal preloads, which prevent separation, and hence peak nonlinear response frequencies usually occur between the full stick and gross slip resonances. Ignoring separation, the dynamics in microslip depends on the coefficient of friction μ, the excitation amplitude, and the normal preload. Consequently, at the system level, nonlinear microslip behavior may be approximately predicted by the nondimensional parameter $ρ=μ|N0|/|F|$ where $|N0|=∑iN0,i=∑ikz,ivi0$ is the sum of normal contact preloads at all the interfaces, and $|F|$ is the amplitude of excitation. As ρ decreases, the conditions at the interfaces become more conducive to slip, leading to more localized contact regions entering slip condition for longer fractions of the cycle period. This increase in the microslip level is reflected in the frequency response. As the energy dissipation at the contact increases, there is a decrease in the normalized response amplitude. There is also a reduction in the peak frequency. Intuitively this is due to an effective decrease in the interface stiffness caused by the increased time spent in the slip state by a larger part of the contact surface. Energy dissipation due to friction at full stick is theoretically zero. The energy dissipated per cycle at large gross slip is proportional to the amplitude of motion. In contrast, the energy dissipated per cycle in viscous damping is proportional to the square of the amplitude of motion, i.e., it increases much faster with the amplitude of motion. Thus, both the case of stick and the case of large gross slip do not lead to the maximum energy dissipation per cycle. The maximum energy dissipation and minimum normalized response occur in a microslip state between the two linear states [140,153,154]. Eventually as the microslip level increases, the optimum level of frictional damping is achieved where the normalized response is minimum at the desired DoF. If the value of ρ is further decreased, the amount of energy dissipated per cycle decreases and the normalized response starts increasing until the gross slip condition is reached.

Fig. 16
Fig. 16
Close modal

## 5 Projection-Based Reduced Order Models

### 5.1 Spatial Reduction by Projection.

Projection-based ROMs seek a reduction basis $ΦROM$, which captures the dominant spatial correlations in the dynamics and can therefore be used to project them into the reduced generalized DoFs p such that
$qr≈ΦROMp$
(31)
Applying the reduction in Eq. (31) to the frequency domain dynamics obtained using HBM in Eq. (30), one may obtain
$DROMhp¯h=ΦROMTf¯Eh+ΦROMTf¯Ch DROMh=ΦROMTDhΦROM$
(32)
Many ROMs use the formulation in Eq. (32) where both linear and nonlinear DoFs are reduced simultaneously by projection after the formulation of the nonlinear tions. Alternatively, it is possible to employ separate reductions for the linear and nonlinear DoFs. The matrix equation in Eq. (30) may be separated into two sets of equations, corresponding to the linear and nonlinear DoFs from $qr$ in Eq. (27) as follows:
$[DNNhDNLhDLNhDLLh][q¯Nhq¯Lh]=[f¯E,Nhf¯E,Lh]+[f¯C,Nh0]; q¯rh=[q¯Nhq¯Lh]$
(33)
Individual reductions may be applied to the linear and nonlinear DoFs to project them onto the reduced DoFs sL and sN, respectively, as follows:
$qN≈ΦROM,NsN; qL≈ΦROM,LsL$
(34)
where $ΦROM,N$ and $ΦROM,L$ are the nonlinear and linear reduction bases, respectively. Substituting Eq. (34) into Eq. (33) and simplifying, one may observe that the reduced linear DoFs are functions of the reduced nonlinear DoFs and excitation force as follows:
$DLL,ROMhs¯Lh=−DLN,ROMhs¯Nh+ΦROM,LTf¯E,Lh ; DAB,ROMh=ΦROM,ATDABhΦROM,B; A={L,N} , B={L,N}$
(35)
Using Eqs. (33)(35), it is possible to describe the dynamics of the system only in terms of the nonlinear reduced DoFs as follows:
$HROMhs¯Nh=g¯E,Lh+ΦROM,NTf¯E,Nh+ΦROM,NTf¯C,Nh HROMh=DNN,ROMhDNL,ROMh(DLL,ROMh)−1DLN,ROMhg¯E,Lh=−DNL,ROMh(DLL,ROMh)−1ΦROM,LTf¯E,Lh$
(36)

Any relevant linear reduction technique, such as those discussed in Sec. 3.2, may be applied to obtain the linear reduction basis $ΦROM,L$ [125127]. Before the advent of methods for reduction of the nonlinear DoFs, CB-CMS was commonly used by many researchers as a method of reducing only the linear DoFs in the nonlinear blisk model. For such reductions, the nonlinear DoFs are retained as master DoFs, and linear DoFs are reduced as slave DoFs [46,50]. Sometimes, reductions may also be applied sequentially, by carrying out only the reduction of the linear coordinates first and then applying reduction again to nonlinear and reduced linear DoFs [80,92]. Alternatively, some methods apply reductions to the frequency response function matrices, which describe the relationship between nonlinear force and response [148]. The rest of Sec. 5 is dedicated to the discussion of methods to formulate a nonlinear reduction basis where $ΦROM,N$ in Eq. (34) and the nonlinear portion of $ΦROM$ in Eq. (31) are not identity matrices. Such reduction bases may be obtained by augmenting linear reduction bases with modal derivatives to capture the effect of nonlinearity on the linear modes [173175]. However, this strategy can be computationally expensive and the number of reduction modes yielded might be quite large, and consequently not offer sufficient reduction in DoFs. Hence, specific techniques have been developed to estimate these nonlinear modal dependencies and generate the appropriate reduction basis in a computationally efficient manner.

### 5.2 Proper Orthogonal Modes.

The proper orthogonal decomposition (POD), also sometimes called the Karhunen–Loeve decomposition, was first proposed independently by a number of sources [176180]. POD has since gained popularity for nonlinear analysis in a number of different fields where it may be variably referred to as principal component analysis, empirical orthogonal function, or factor analysis [82]. As applied to structural dynamic analysis, it is used to synthesize a set of orthogonal vectors or modes, which capture the dominant (most energetic) motions of the system. This is achieved by carrying out the singular value decomposition (SVD) of the matrix $Q=[qr(t1) qr(t2)…]$, which is comprised of columns containing snapshots of the time domain responses of the system at times $t1,t2,…$ at all the DoFs (linear and nonlinear). The left singular vectors $Θ=[θ1θ2…]$ of Q are known as POMs. Each of the POMs is associated with a corresponding singular value, whose magnitude represents the relative dominance of POMs in the motion described in Q. The POMs may alternatively obtained as the eigenvectors of $QQT$ and the relative dominance indicated by the corresponding eigenvectors [82]. POMs may also be interpreted as the best least-squares fit of a linear representation of a single synchronous NNM describing the motion [181].

Proper orthogonal modes may be used as basis vectors in the reduction basis $ΦROM$ in Eq. (31). As more POMs are included, the reduction is able to approximate the dynamics with sufficient accuracy. However, the main drawback of this method lies in the determination of the matrix Q, which requires extensive experimental observations or computationally expensive full order simulations. In the case of blisks, the practical applications of POD for ROM formulation is further limited by the fact that uncertainty is fundamentally tied into the reduction problem due to mistuning and variability in system behavior between operating conditions. Hence, researchers focus on obtaining reductions by addressing the dominant physical phenomenon underlying the blisk dynamics such as the behavior of contacts, which remains largely invariant to these uncertainties. However, POMs do serve as a convenient validation tool for these physics-based ROMs, under the specific conditions they were generated [80]. For instance, POMs should not be used for validation of responses involving modes, which were not excited in the snapshots. Thus, a POM ROM developed from the response of a mistuned blisk to a given EO excitation, generally cannot be used to accurately predict responses to a different EO excitation. Nevertheless, POMs may be used for model reduction in other systems, which do not present large uncertainties in their operational envelope [87].

### 5.3 Reduced Order Models for Intermittent Contacts.

One of the concepts that led to the development of modern ROMs for systems with intermittent contacts such as cracked blades or blisks [124,182,183] was the identification of the piecewise linear behavior of these systems [184,185]. While such a structure may exhibit highly nonlinear characteristics, at any instant of time during its motion, it may be described by as a linear system with particular constraints at the interfaces. In particular, researchers postulated that the nonlinear system behavior may be captured by linear modes of these systems. In Ref. [90], it was shown that the free and gross slip linear system frequencies could be used to formulate an effective natural frequency of the nonlinear structure called the bilinear frequency when the difference between the linear regions was small. This was validated with NNMs obtained using a two degree-of-freedom model representing a cracked beam. Order reduction for mass-spring systems with local nonlinearities was also obtained in Ref. [186] by using projection bases composed by augmentation of linear normal modes with basis functions with appropriate discontinuities at the locations of nonlinearity.

These ideas were combined and modified to develop BLMs as a feasible projection basis for cracked structures. BLMs were treated as approximations of POMs of the nonlinear system [91]. As shown in Fig. 17, BLMs $ΦBLM=[ΦfΦg]$ are comprised of selectively chosen normal modes of the linear systems corresponding to the free ($Φf$) and gross slip ($Φg$) conditions. Two selection criteria for these modes are provided in Ref. [91]. In one of the methods, no POM calculations are required and the BLM basis is first populated by linear modes from both the free and gross slip systems lying within a frequency range of interest (usually the frequency range of excitation) since they most resemble the nonlinear motion. However, these modes are usually not sufficient to span the nonlinear motion space and hence must be augmented with more linear modes. The additional modes are augmented in the increasing order of some metric, which determines the closeness of the linear modes, which are potential BLMs to the already selected BLMs. This metric is either the Euclidean angle between the vectors describing the modes, or another angle metric which uses a stiffness-based norm [91]. In the other method described in Ref. [91], POMs $Θcm$ obtained from simulations of another FE model of the structure with a coarser mesh are used. These POMs are less expensive to calculate for the coarser mesh model than the fine mesh model for which the BLMs are developed. The POMs in $Θcm$ can then be spatially interpolated to obtain approximate POMs for the finer mesh model $Θfm$. The projection basis vectors in $ΦBLM$ may then be chosen from among the linear modes, to minimize the error between $Θfm$ and its projection on $ΦBLM$, which is given by $||Θfm−ΦBLM(ΦBLMTΦBLM)−1ΦBLMTΘfm||$. Thus, the BLM basis is optimized to describe as much of the space described by the approximate POMs as possible. Figure 18 shows the comparison between BLM ROM response near the first bending mode of a cracked plate to the FE model and other CMS-based ROMs [91]. It is seen that the response may be captured with just a few BLMs.

Fig. 17
Fig. 17
Close modal
Fig. 18
Fig. 18
Close modal

Another ROM formulation technique for intermittent contacts called the BAA was proposed in Ref. [93]. The BAA method aims to directly approximate the amplitude of the nonlinear periodic motion and assumes that the structure enters the free and gross slip stage exactly once every period as shown in Fig. 19. The motion is linear at all times. During the free and gross slip stages, the motion may be completely described by $Φf$ and $Φg$. It is assumed that the transition between the stages occurs very fast, and therefore any motion corresponding to partially closed contacts may be ignored. Consequently, the motion during the transition is spanned by a the columns of the matrix $Φt$, which is a subset of the spaces spanned by $ΦBLM$. Numerically, $Φt$ is obtained as the left singular vectors corresponding to the few largest singular values of $ΦBLM$. The piecewise linear periodic response is then solved in reduced coordinates of the basis corresponding to the applicable linear condition, while enforcing compatibility criteria at the transition times. This allows the determination of the amplitude of the nonlinear motion and calculation of the response at each frequency.

Fig. 19
Fig. 19
Close modal
Additional effects introduced into the dynamics by presence of preloads at the contact interfaces require careful consideration. In general, the problem of determining the static deflection based on a known preload can be overdetermined and hence can have multiple solutions [92]. In most simulations, the static deflection is usually estimated by a quasi-static application of preloads at the interfaces, which admits a unique solution $qPS0$ [80,92]. This is a good model of the preload effects when the mean deflection during motion is not affected by the dynamics of the structure [50,80], i.e., there is minimal coupling between the dynamic EOMs (corresponding to h >0 in Eq. (30)) and the static or 0th harmonic term. In this case, the BLMs, along with the quasi-statically determined deflection due to prestress, form an adequate projection basis for the dynamics. However, this assumption regarding the decoupling of the statics and dynamics is not always true, especially for prestressed structures with intermittent contacts where frictional effects are significant. Such a case was studied in Ref. [92], where it was determined that in order to formulate accurate ROMs, which captured such dynamics, more vectors were required in the projection basis in addition to the BLMs and the quasi-static deflection. The role of these additional vectors is to capture the effect of the dynamics on the statics. This was achieved, by analyzing the fictitious normal static frictional forces, which would be generated at the contacts if the deflection of the structure resembled the BLMs. At the ith node pair, these fictitious normal force $NBLM,i$ may be expressed as
$NBLM,i=kz,ivBLM,i$
(37)
where $vBLM,i$ is the component of a BLM mode expressed in relative coordinates along the normal direction at ith contact node pair. A vector of these fictitious forces $fBLM,N$ may then be formulated for all nonlinear DoFs, whose entries correspond to $NBLM,i$ at the appropriate locations and zeros everywhere else. A static deflection $qBLM0$ for the entire structure may be calculated based on the application of these fictitious forces for each BLM as follows:
$qBLM0=Kf−1[fBLM,N0]$
(38)

When considered together, the quasi-statically determined displacement $qPS0$, the BLMs themselves $ΦBLM$, and $qBLM0$'s calculated for all BLMs may be used to formulate the projection ROM. This combined static deflection and BLM ROM or SD + BLM ROM, as it is referred to in Ref. [92], performed well in replicating the nonlinear dynamics for a blade-like cracked plate when regular BLMs were insufficient. Note that the static solution in Eq. (38) may be carried out in reduced coordinates if a reduction (such as CB-CMS) has been pre-applied to the linear DoFs as discussed in Sec. 5.1.

BLMs are based on information generated from only two limiting linear cases of a structure with contact nonlinearity, which is piecewise linear. Though sufficient in many cases, complex structural motion and preload distributions at the interfaces may render this bilinear assumption invalid if various distinct spatial subregions of the contact interface open and close single or multiple times during the motion, with different subregions possibly undergoing the transition at different times during the cycle. Under such circumstances, the modes of other linear conditions must be considered for inclusion into the ROM projection basis. Aptly called PLMs, these are the modes of all piecewise linear systems (observed at different time instances of nonlinear motion), which contribute to the nonlinear dynamics. The creation of ROMs using PLMs was discussed in Ref. [95]. The challenging aspect of determining the PLMs is the identification of the specific linear systems corresponding to time varying BCs at the interfaces. If the structure only has localized contact nonlinearities, even its nonlinear motion is dominated by a few linear modes at any specific excitation frequency. Hence, the BCs may be estimated by examining the kinematics of the contact assuming motion along these dominant modes as illustrated in Fig. 20.

Fig. 20
Fig. 20
Close modal
Given a static normal relative displacement due to prestress $vPS,i0$ at the ith contact node pair, the total estimated normal displacement $v̂i$ along the same node pair due to the prescribed kinematics is
$v̂i=vPS,i0+α vDOM,i$
(39)
where $vDOM,i$ is the component of the dominant mode $ϕDOM$ expressed in relative coordinates along the ith contact node pair in the normal direction and $α∈[αminαmax]$ is a scalar value representing the assumed modal amplitude. If $v̂PLM,i>0$, the ith node pair is estimated to be in contact (slip), or otherwise it is assumed open (free). Consequently, for each value of α, a set of linear BCs is obtained over the entire interface, with linear normal contact stiffnesses being enforced only for nodes in contact. αmin and αmax are determined such that no new BCs are predicted beyond these values. By enforcing the BCs obtained for the entire range of α, one may obtain the piecewise linear systems and consequently the PLMs, which likely contribute to the response. Note that the BCs appear as stiffness constraints only since the mass of the structure remains unchanged during contact. The PLM $ϕPLM$ for the BC corresponding to a given α and dominant mode may be obtained as normal modes of the constrained system by solving the generalized eigenvalue problem
$KPLMϕPLM=λPLMMfϕPLM$
(40)

where KPLM is the constrained stiffness matrix corresponding to a specific BC and λPLM is the eigenvalue associated with $ϕPLM$. The reduction basis is then formulated by collating all such $ϕPLM$ for which λPLM lie in the frequency range of interest. As these $ϕPLM$ are not linearly independent, conditioning (such as SVD-based conditioning) might be required to prevent numerical problems during simulations. In Ref. [95], the PLM ROM was validated for various structures with intermittent contacts such as a cracked plate and contacting coaxial cylinders considering cases with and without preloads.

### 5.4 Reduced Order Models for Microslip Using Linear Modes.

Projection ROMs for microslip contacts share many of the ideas related to ROMs for intermittent contacts discussed in Sec. 5.3. The complexity of the ROMs is often predicated on the complexity of the full order nonlinear dynamics, specifically the number of dominant linear modes in the nonlinear motion. For instance, in blisks with ring dampers, the frequency separation between full stick and gross slip modes is relatively small. Consequently, the overall motion of the blisk itself does not vary much in the microslip regime between full stick and gross slip [5,62,79]. This is due to the simple geometry of the ring damper and also its location near the root of the blade, where the motion is generally small compared to the blade tips where the blisk motion is usually the highest. Hence, using an argument similar to that for BLMs, ROMs may be obtained for blisks with ring dampers by projecting the dynamics onto $ΦROM=[ΦsΦg]$, where $Φs$ and $Φg$ represent a subset of the stick and gross slip modes in the frequency range of interest, respectively [5,62]. Moreover, for tuned blisks, only the modes with the same ND as the EO of excitation need be considered.

It was shown in Ref. [79] that due to the simplicity of the reduction basis, it is possible to obtain an equivalent linear representation of the nonlinear dynamics for given modal amplitudes. A ROM may be formulated by approximating the nonlinear friction forces projected into the reduced domain in terms of equivalent stiffness and damping terms as follows:
$ΦROMTf¯C,eq1≈g¯C,eq1=[iΓC,eq°KC,eq+KC,eq]p¯eq1 ΓC,eq=ΓC,eq(p¯eq1) ; KC,eq=KC,eq(p¯eq1)$
(41)

where $f¯C,eq1$ and $g¯C,eq1$ are the equivalent first harmonic frictional force vectors in the full order and reduced order domains. $p¯eq1$ is the equivalent modal displacement vector in the reduced domain. $ΓC,eq$ is the equivalent damping matrix. $KC,eq$ is the equivalent stiffness matrix, and $°$ represents the elementwise (Hadamard) product of two matrices. After substituting $g¯C,eq1$ in Eq. (32), the reduced system may be solved for the first harmonic to obtain $p¯eq1$ whose amplitude is approximately equal to that of the nonlinear system for a given frequency. Prestress effects can also be addressed using this method through inclusion of an additional stiffness term in the reduced dynamic equation, details regarding which may be found in Ref. [79].

$ΓC,eq$ and $KC,eq$ are functions of $p¯eq1$, and are full matrices in general. General individual entries of these matrices $γC,eq$ and $kC,eq$, respectively, represent the equivalent damping and stiffness due to the nonlinear force encountered along one reduction mode due to motion along another mode. The procedure for the calculation of these elements is shown in Fig. 21. First, a quasi-static harmonic motion is applied to the blisk along a reduction mode with amplitude α. This is similar to the modal amplitude based estimates used in the PLM procedure discussed in Sec. 5.3. However, unlike the static approximation of modal amplitude used for PLM, the motion is applied for an entire harmonic period. Based on this motion, a hysteresis cycle is constructed in the modal reduced coordinates, and the average elastic energy stored and the total energy dissipated during a period are calculated. A work-energy equivalence is then applied, where these dissipated and stored energies are equated with equivalent quantities for a linear system, which are expressed in terms of $γC,eq$ and $kC,eq$, respectively. From this equivalence, $γC,eq$ and $kC,eq$ may then be determined as functions of the modal amplitude α. Note that although $ΓC,eq$ and $KC,eq$ are functions of the motion along multiple modes $p¯eq1$, their individual elements $γC,eq$ and $kC,eq$ may be precalculated as single dimensional functions of α.

Fig. 21
Fig. 21
Close modal

This procedure enables the creation of very fast ROMs, where the nonlinearity is captured by equivalent single harmonic calculations. The need for AFT is also completely eliminated as the equivalent nonlinear forces may be calculated directly in the reduced coordinates in the frequency domain. Another advantage of this ROM is that full knowledge of the contact model and its parameters is not required. The contact may be treated as a black box, which outputs frictional forces for enforced displacements. Results from the successful application of this ROM to a tuned blisk with a frictional ring damper are shown in Fig. 22, where the nonlinear frequency domain responses of the ROM are compared to the results from the full order FE model obtained using time marching [79].

Fig. 22
Fig. 22
Close modal

Unlike blisks with ring dampers, the motion of blisks with other types of microslip contacts may be much more complex. For instance, shrouded blisks are designed to be frictionally damped by contacts between shrouds attached to the blades near the tips. The full stick and gross slip modes of shrouded blisks are very different from each other and are well separated from each other in terms of their natural frequencies [11]. Consequently, the change between the full stick and gross slip motion with changes in microslip level is much more gradual than in the case of shrouded blisks than it is for blisks with ring dampers. Hence, coherences in the nonlinear dynamics of the shrouded blisk in microslip are not spanned by the union of $Φs$ and $Φg$ and other projection bases are required for reduction. A method for creating ROMs for such shrouded blisks and other structures, which show such gradual transition during microslip, using an AMP basis, was proposed and validated in Ref. [80]. The procedure for obtaining the AMPs is illustrated in Fig. 23.

Fig. 23
Fig. 23
Close modal
The AMP procedure is very similar to the one for PLMs as both were developed somewhat concurrently for different applications. The AMP method focuses on estimating the tangential BCs at the contact interfaces (which were ignored for PLMs). These BCs may then be used to obtain the intermediate linear systems between the full stick and gross slip whose normal modes can approximate the nonlinear dynamics. Similar to PLMs, due to the localized nature of contact nonlinearities some dominant bulk motion $ϕDOM$ (usually a mode belonging to $Φs$ or $Φg$) is assumed and applied to the structure with some modal amplitude $α∈[αminαmax]$ in addition to the deflection due to prestress $qPS0$. Then, the statically estimated forces at contact interfaces due to this assumed deflection may be calculated at each node pair i in the normal and tangential directions as follows:
$N̂i=kz,i(vPS,i0+α vDOM,i); T̂x,i=kx,i(α uDOM,x,i)$
(42)
where $N̂i$ and $T̂x,i$ are the estimated normal and tangential forces, respectively, expressed in relative coordinates at contact node pair i. $vDOM,i$ and $uDOM,x,i$ are the respective normal and tangential components of $ϕDOM$ at the same location. $vPS,i0$ is the normal component of $qPS0$ at node pair i, which is obtained assuming contact at gross slip, and hence does not have tangential components. The predicted BC may then be obtained by comparing $μN̂i$ and $T̂x,i$ at each node pair. Node pair i is predicted to be in stick when $T̂x,i<μN̂i$ or slip otherwise. Based on these BCs, one may then mask the diagonal stiffness matrix $Ks$, which is the difference between the linear stiffness of the fully stuck ($Ks$) and free ($Kf$) structure and may be expressed in relative coordinates as
$ΔKs=Ks−Kf=[⋱0⋯⋯⋯0kx,i00⋯⋮0ky,i0⋯⋮00kz,i⋯⋮⋮⋮⋮⋱]$
(43)
As illustrated in Fig. 23, the tangential stiffnesses may be removed or retained from $ΔKs$ based on the estimated BC at the corresponding node pair, to formulate the matrix $ΔKAMP$. The stiffness matrix corresponding to the intermediate linear system may be obtained as $KAMP=Kf+ΔKAMP$ and the AMP reduction mode $ϕAMP$ may be obtained by solving the following eigenvalue problem in the frequency range of interest:
$KAMPϕAMP=λAMPMfϕAMP$
(44)

where λAMP is the eigenvalue associated with $ϕAMP$. The modes $ϕAMP$ obtained by repeating this procedure for multiple $ϕDOM$ and various values of α may be collated into a reduction basis. When used directly without any validating full order or experimental results available for comparison, convergence studies might be required to determine the ideal number of AMP modes that should be retained in the reduction basis. However, for many cases, it is sufficient to calculate more AMP vectors than strictly necessary (since the linear calculations are fast) and use SVD conditioning (which is also recommended for numerical conditioning) to obtain a smaller basis, which can capture the dynamics [80].

Results for the application of the AMP procedure to a mistuned shrouded blisk with 27 sectors are shown in Figs. 24 and 25.

Fig. 24
Fig. 24
Close modal
Fig. 25
Fig. 25
Close modal

The frequency response of the full order baseline blisk model under full stick (linear) and a high level of microslip to an EO 1 excitation in a high modal density region where natural frequencies of the system are clustered close to each other may be seen in Fig. 24. Since the blisk is mistuned, each plotted line representing the frequency response amplitude of each blade is different. Figure 24 also shows the linear modes at the contact DoFs corresponding to the full stick frequency peaks and the spatial FFT of a vector of modal amplitudes for identical nodes at each sector. For a tuned blisk, the spatial FFT would only have a single nonzero component corresponding to its ND. However, for this blisk, which has small stiffness mistuning, multiple ND components are seen in the modes. There are also three modes with dominant ND1 patterns in a narrow frequency range, which can only occur for a mistuned blisk as the ND1 nodes for tuned blisks appear in pairs. At high microslip levels, all these modal ND components also interact with each other at the nonlinear interfaces to give rise to very complex dynamics. However, as shown in Fig. 25, the AMP ROM is able to reproduce the nonlinear response with only 238 reduction modes and offers significant computational time savings compared to the original FE model with more than 50,000 DoFs.

For the AMP procedure, the same reduction basis is used to reduce EOMs of all harmonics. Details regarding dominant mode selection and the generation of the basis for cases where higher harmonics (h >1) have significant dynamic contribution are provided in Ref. [80]. An alternative reduction technique called the JP is described in Ref. [100]. The method for generating the BCs is the same for AMP and JP. However, in the JP method, the BCs are enforced not on the linear system (as it is for AMP) but instead on the system represented by the multiharmonic Jacobian ($BDiag(Dh)$ in Eq. (30)) after a static condensation. The JP method was validated for a cantilevered shrouded blade in Ref. [100], which also provides an error metric for assessing the accuracy of the ROM without a convergence study.

## 6 Summary and Discussion

This main goal of this paper is to discuss recent advancements in creating nonlinear reduced order models for bladed disks (blisks) with contact interfaces, which are widely used in industrial turbomachinery such as jet engines or power generating gas and steam turbines. The harsh operating conditions and high cyclic stresses encountered by blisks during operation render their design against failure especially challenging. Hence, extensive research into the dynamics of these structures has been carried out over the past half-century. There are inherent uncertainties associated with blisk manufacturing and its operating conditions, which may have significant unwanted effects on their dynamic characteristics. Presently, the probabilistic approach to blisk design is the preferred method, which requires simulations of many blisk models to characterize the effects of these uncertainties. Limits on available computational power have placed a premium on fast, adaptable, and accurate reduced order models for such studies.

More recently, with improvement in simulation tools and methods, researchers have also begun to realize the importance of accounting for nonlinearities such as contact damping mechanisms or cracks in these blisks. The state of the art in the modeling of such intermittent and frictional contact nonlinearities using both high fidelity and reduced order models was discussed in this paper, and a summary of the technical background and important developments in the field were provided. The cyclic symmetric nature of a nominal blisk and the characteristic linear dynamics of such structures were discussed. The effect of uncertainties and the linear reduced order models employed to study them were described also.

All the nonlinear reduction techniques discussed here focus on obtaining an accurate basis of reduction vectors onto which the nonlinear dynamics may be projected. It has been known that if the nonlinear motion is captured by the projection basis to a sufficient extent, the reduced dynamics can achieve the desired accuracy. The challenge lies in predicting these motions without any a priori nonlinear full order simulations, which would defeat the purpose of such reduced order model generation. To address this problem, the piecewise linear nature of these systems may be leveraged, where at different instants of time during the nonlinear motion, the system is in different linear conditions. When change in the localized contact conditions occurs fast, and the system spends most of the time in one of two limiting linear conditions (e.g., fully open or closed contact interface for a crack), it is sufficient to use these linear modes for reduction. In other cases, due to the geometric properties of the structure, prestress, and other factors, the nonlinear behavior might be more complex. Modes from multiple intermediate systems corresponding to intermediate contact boundary conditions (e.g., a partially open crack) may interact at the interfaces to give rise to these complex dynamics. It becomes difficult to predict which intermediate linear systems are responsible for the dynamics. However, due to the localized nature of contact nonlinearities, dominant motions observed in the overall nonlinear response of the system are similar to the linear modes of the systems under the limiting conditions. Although these modes approximate the bulk motion, using them alone for projection leads to inaccuracies as they do not accurately capture small motions at the interfaces, which dictate the nonlinear forces and dynamics. Instead, most methods use static or quasi-static calculations based on this dominant bulk motion, to predict the most likely boundary conditions at the interface and the corresponding linear systems, whose modes are then used in a basis to achieve the reduction.

These reduction methods do not require a priori nonlinear full order calculations and are often orders of magnitude faster than high fidelity models. The models are also physics-based rather than data-driven and are hence reasonably robust in terms of accuracy when relatively small changes are made to the properties or structure of the blisk, making them excellent candidates for use in probabilistic analyses.

Reduced order models are also based on certain assumptions and approximations. For instance, the presence of dominant nonlinear spatial correlations is assumed. The models lose accuracy when these conditions are not satisfied and are, therefore, not broadly applicable to every type of nonlinearity. In addition, the accuracy of the models depends on the quality of the contact models, which are employed. These contact models often have significant assumptions, which limit their range of applicability.

There is scope for future research in establishing mathematical bounds on reduced order model accuracy and convergence and establish guarantees of performance. Combination and evolution of ideas from various available reduction techniques discussed here are also possible for broadening their application to other types of systems and nonlinearities. Currently, reduced order models are almost always developed based on some high-fidelity computational model of a blisk. A lot of time and effort is required to accurately characterize these large models. Determining the contact parameters at the interfaces is especially challenging. Often, the accuracy of a reduction method is rendered moot by the inaccuracies in the high fidelity model that is being reduced. An example is the need for advanced contact models, either microscopic or at the mesoscale. Hence, future research in the area may be trending toward creating sufficiently flexible and capable reduced models, which may be tuned directly using experimental results. It is also often difficult to identify all the sources of nonlinearities in a complex structure like a blisk, based on limited experimental measurements. Newer system identification methods and advancements in neural network technology might hold the key to creating black-box or gray-box models for these systems wherein not all model parameters might be easily accessible or interpretable. Correspondingly, the paradigm for model reduction techniques may need to evolve to address such models in the future.

## Acknowledgment

The authors would like to acknowledge Mr. Adegbenga Odofin, Dr. Seunghun Baek, and Dr. Weihan Tang for kindly providing explanatory materials for their work which is discussed in this paper.

## Nomenclature

### Symbols

Symbols
A =

C =

damping matrix

D =

dynamic stiffness matrix

e =

column of complex Fourier matrix

E =

complex Fourier matrix

EO =

engine order

f =

force vector

$|F|$ =

amplitude of excitation force

g =

vector of generalized forces in reduced domain

h =

harmonic index number

H =

receptance matrix for nonlinear degrees of freedom

i =

square root of −1

k, k =

nodal contact stiffness value and vector of contact stiffnesses

K =

stiffness matrix

l =

vector of generalized forces in reduced domain

M =

mass matrix

n =

N, N =

scalar and vector of normal contact force(s) at single or multiple node pairs

ND =

number of nodal diameters

p =

spatial harmonic

p =

vector of displacements in generalized reduced coordinates

q =

vector of displacements in physical or higher order coordinates

Q =

matrix of collated displacements in time domain

R =

rotation matrix

s =

vector of displacements in generalized reduced coordinates

S =

sector

t =

time

T, T =

scalar and vector of tangential contact force(s) at single or multiple node pairs

u =

relative tangential displacement of node pair

v =

relative normal displacement of node pair

w =

slip at node pair

α =

estimated scalar modal amplitude

β =

sector angle

$γ,Γ$ =

structural damping value and matrix

Δ =

change in quantity

$η$ =

vector of displacements in generalized reduced coordinates

$θ,Θ$ =

individual modes and matrix of proper orthogonal modes

$λ,Λ$ =

eigenvalues and diagonal matrix of eigenvalues

μ =

coefficient of friction

Σ =

sum

$ϕ, Φ$ =

single mode and matrix of modes

$ψ,Ψ$ =

single mode and matrix of modes

ω =

circular frequency of excitation

Ω =

rotational speed

≈ =

approximated as

$°$ =

Hadamard/Schur product for entrywise multiplication of matrices

⊗ =

Kronecker product of matrices

$ℜ$ =

real part of a complex number

### Subscripts

Subscripts
a =

physical global coordinate system

A =

variable for degrees of freedom in matrix partition

AMP =

B =

variable for degrees of freedom in matrix partition

BLM =

reduction using bilinear modes

C =

contact

cm =

coarse mesh

CB =

Craig-Bampton component mode synthesis

DOM =

dominant modes

E =

excitation

eq =

equivalent single harmonic quantity

f =

free system

fm =

fine mesh

g =

gross slip system

gen =

generating matrix

H =

high surface sector degrees of freedom

L =

low surface sector degrees of freedom

i =

index for contact node-pair number

I =

interior sector degrees of freedom

j =

general index or integer

l =

local contact coordinate system

L =

linear degrees of freedom

M =

master degrees of freedom

max =

maximum value

min =

minimum value

N =

nonlinear degrees of freedom

o =

tuned system quantity

O =

noncontact degrees of freedom

PLM =

reduction using piecewise linear modes

PS =

prestress

r =

relative coordinate system

ROM =

reduced order model

s =

full stick system

S =

slave degrees of freedom

S1 =

contacting surface one degree-of-freedom

S2 =

contacting surface two degrees-of-freedom

t =

transition between systems

x, y, z =

local tangential and normal degrees of freedom

0 =

static quantity

### Superscripts

Superscripts
h =

harmonic index number

H =

complex conjugate transpose of a matrix

S =

sector level matrix

T =

transpose of a matrix

### Abbreviations

Abbreviations
AF(s) =

amplification factor(s)

AFT =

alternating frequency-time

AMP(s) =

BAA =

bilinear amplitude approximation

BC(s) =

boundary condition(s)

BLM(s) =

bilinear mode(s)

CB-CMS =

Craig-Bampton component mode synthesis

CDF =

cumulative distribution function

CMM =

component mode mistuning

DoF(s) =

degree(s) of freedom

EO =

engine order

EOM(s) =

equation(s) of motion

FE =

finite element

FMM =

fundamental model of mistuning

HBM =

harmonic balance method

JP =

Jacobian projection

ND(s) =

nodal diameter(s)

NNM(s) =

nonlinear normal mode(s)

PLM(s) =

piecewise linear mode(s)

POD =

proper orthogonal decomposition

POM(s) =

proper orthogonal mode(s)

ROM(s) =

reduced order model(s)

SNM =

subset of nominal modes

SVD =

singular value decomposition

### Accents

Accents
$¯$ =

quantity in frequency domain

$˙$ =

single time derivative

$¨$ =

double time derivative

$̂$ =

estimate

$̃$ =

constrained quantity

## References

1.
Srinivasan
,
A. V.
,
1997
, “
Flutter and Resonant Vibration Characteristics of Engine Blades
,”
ASME J. Eng. Gas Turbines Power
,
119
(
4
), pp.
742
775
.10.1115/1.2817053
2.
Ewins
,
D. J.
,
2010
, “
Control of Vibration and Resonance in Aero Engines and Rotating Machinery—An Overview
,”
Int. J. Pressure Vessels Piping
,
87
(
9
), pp.
504
510
.10.1016/j.ijpvp.2010.07.001
3.
Castanier
,
M. P.
, and
Pierre
,
C.
,
2002
, “
Using Intentional Mistuning in the Design of Turbomachinery Rotors
,”
AIAA J.
,
40
(
10
), pp.
2077
2086
.10.2514/2.1542
4.
Castanier
,
M.
, and
Pierre
,
C.
,
2006
, “
Modeling and Analysis of Mistuned Bladed Disk Vibration: Current Status and Emerging Directions
,”
J. Propul. Power
,
22
(
2
), pp.
384
396
.10.2514/1.16345
5.
Tang
,
W.
,
Baek
,
S.
, and
Epureanu
,
B.
,
2017
, “
Reduced-Order Models for Blisks With Small and Large Mistuning and Friction Dampers
,”
ASME J. Eng. Gas Turbines Power
,
139
(
1
), p.
012507
.10.1115/1.4034212
6.
Baek
,
S.
, and
Epureanu
,
B. I.
,
2017
, “
Reduced-Order Models of Blisks With Small Geometric Mistuning
,”
ASME J. Vib. Acoust.
,
139
(
4
), p.
041003
.10.1115/1.4036105
7.
D'Souza
,
K.
,
Jung
,
C.
, and
Epureanu
,
B. I.
,
2013
, “
Analyzing Mistuned Multi-Stage Turbomachinery Rotors With Aerodynamic Effects
,”
J. Fluids Struct.
,
42
, pp.
388
400
.10.1016/j.jfluidstructs.2013.07.007
8.
Fitzner
,
C.
,
Epureanu
,
B. I.
, and
Filippi
,
S.
,
2014
, “
Nodal Energy Weighted Transformation: A Mistuning Projection and Its Application to Flade (TM) Turbines
,”
Mech. Syst. Signal Process.
,
42
(
1–2
), pp.
167
180
.10.1016/j.ymssp.2013.08.027
9.
,
D. S.
,
1998
, “
The Maximum Factor by Which Forced Vibration of Blades Can Increase Due to Mistuning
,”
ASME J. Eng. Gas Turbines Power
,
120
(
1
), pp.
115
119
.10.1115/1.2818061
10.
Avalos
,
J.
,
Murthy
,
R.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2014
, “
On the Effects of Blade-Disk Interface Mistuning on the Response of Integrated Bladed Rotors
,”
ASME
Paper No. DETC2013-13352.10.1115/DETC2013-13352
11.
Mitra
,
M.
,
Zucca
,
S.
, and
Epureanu
,
B. I.
, 2016, “
Effects of Contact Mistuning on Shrouded Blisk Dynamics
,”
ASME
Paper No. GT2016-5781210.1115/GT2016-57812.
12.
Besern
,
F. M.
,
Kielb
,
R. E.
, and
Key
,
N. L.
,
2015
, “
Forced Response Sensitivity of a Mistuned Rotor From an Embedded Compressor Stage
,”
ASME
Paper No. GT2015-43032.10.1115/GT2015-43032
13.
Bhartiya
,
Y.
, and
Sinha
,
A.
,
2012
, “
Reduced Order Model of a Multistage Bladed Rotor With Geometric Mistuning Via Modal Analyses of Finite Element Sectors
,”
ASME J. Turbomach.
,
134
(
4
), p. 041001.10.1115/1.4003224
14.
Han
,
Y.
,
Murthy
,
R.
,
Mignolet
,
M. P.
, and
Lentz
,
J.
,
2014
, “
Optimization of Intentional Mistuning Patterns for the Mitigation of the Effects of Random Mistuning
,”
ASME J. Eng. Gas Turbines Power
,
136
(
6
), p. 062505.10.1115/1.4026141
15.
Kielb
,
R. E.
,
Hall
,
K. C.
,
Hong
,
E.
, and
Pai
,
S. S.
,
2006
, “
Probabilistic Flutter Analysis of a Mistuned Bladed Disks
,”
ASME
Paper No. GT2006-90847.10.1115/GT2006-90847
16.
Beck
,
J. A.
,
Brown
,
J. M.
,
Slater
,
J. C.
, and
Cross
,
C. J.
,
2013
, “
Probabilistic Mistuning Assessment Using Nominal and Geometry Based Mistuning Methods
,”
ASME J. Turbomach.
,
135
(
5
), p.
051004
.10.1115/1.4023103
17.
Craig
,
R. R.
, and
Bampton
,
M. C. C.
,
1968
, “
Coupling of Substructures for Dynamic Analyses
,”
AIAA J.
,
6
(
7
), pp.
1313
1319
.https://arc.aiaa.org/doi/10.2514/3.4741
18.
Hurty
,
W. C.
,
1965
, “
Dynamic Analysis of Structural Systems Using Component Modes
,”
AIAA J.
,
3
(
4
), pp.
678
685
.10.2514/3.2947
19.
,
R.
,
Castanier
,
M. P.
, and
Pierre
,
C.
,
2001
, “
Component-Mode-Based Reduced Order Modeling Techniques for Mistuned Bladed Disks—Part I: Theoretical Models
,”
ASME J. Eng. Gas Turbines Power
,
123
(
1
), pp.
89
99
.10.1115/1.1338947
20.
,
R.
,
Castanier
,
M. P.
, and
Pierre
,
C.
,
2001
, “
Component-Mode-Based Reduced Order Modeling Techniques for Mistuned Bladed Disks—Part II: Application
,”
ASME J. Eng. Gas Turbines Power
,
123
(
1
), pp.
100
108
.10.1115/1.1338948
21.
Bonney
,
M. S.
,
Kammer
,
D. C.
, and
Brake
,
M. R. W.
,
2015
, “
Fully Parameterized Reduced Order Models Using Hyper-Dual Numbers and Component Mode Synthesis
,”
ASME
Paper No. DETC2015-46029.10.1115/DETC2015-46029
22.
Beck
,
J. A.
,
Brown
,
J. M.
,
Cross
,
C. J.
, and
Slater
,
J. C.
,
2014
, “
Component-Mode Reduced-Order Models for Geometric Mistuning of Integrally Bladed Rotors
,”
AIAA J.
,
52
(
7
), pp.
1345
1356
.10.2514/1.J052420
23.
Yang
,
M. T.
, and
Griffin
,
J. H.
,
2001
, “
A Reduced-Order Model of Mistuning Using a Subset of Nominal System Modes
,”
ASME J. Eng. Gas Turbines Power
,
123
(
4
), pp.
893
900
.10.1115/1.1385197
24.
Feiner
,
D. M.
, and
Griffin
,
J. H.
,
2002
, “
A Fundamental Model of Mistuning for a Single Family of Modes
,”
ASME J. Turbomach.
,
124
(
4
), pp.
597
605
.10.1115/1.1508384
25.
Lim
,
S. H.
,
,
R.
,
Castanier
,
M. P.
, and
Pierre
,
C.
,
2007
, “
Compact, Generalized Component Mode Mistuning Representation for Modeling Bladed Disk Vibration
,”
AIAA J.
,
45
(
9
), pp.
2285
2298
.10.2514/1.13172
26.
Joshi
,
A. G. S.
, and
Epureanu
,
B. I.
,
2012
, “
,”
ASME J. Vib. Acoust.
,
134
(
5
), p. 051015.10.1115/1.4006880
27.
Nyssen
,
F.
, and
Golinval
,
J. C.
,
2016
, “
Identification of Mistuning and Model Updating of an Academic Blisk Based on Geometry and Vibration Measurements
,”
Mech. Syst. Signal Process.
,
68–69
, pp.
252
264
.10.1016/j.ymssp.2015.08.006
28.
Martel
,
C.
, and
Corral
,
R.
,
2008
, “
Asymptotic Description of Maximum Mistuning Amplification of Bladed Disk Forced Response
,”
ASME J. Eng. Gas Turbines Power
,
131
(
2
), p.
022506
.
29.
Khemiri
,
O.
,
Martel
,
C.
, and
Corral
,
R.
,
2014
, “
Forced Response of Mistuned Bladed Disks: Quantitative Validation of the Asymptotic Description
,”
J. Propul. Power
,
30
(
2
), pp.
397
406
.10.2514/1.B34631
30.
Khemiri
,
O.
,
Martel
,
C.
, and
Corral
,
R.
,
2013
, “
Asymptotic Description of Damping Mistuning Effects on the Forced Response of Turbomachinery Bladed Disks
,”
J. Sound Vib.
,
332
(
20
), pp.
4998
5013
.10.1016/j.jsv.2013.04.030
31.
Lu
,
J.
,
D'Souza
,
K.
,
Castanier
,
M. P.
, and
Epureanu
,
B. I.
,
2017
, “
Nonlinear Parametric Reduced-Order Model for the Structural Dynamics of Hybrid Electric Vehicle Batteries
,”
ASME J. Vib. Acoust.
,
140
(
2
), p.
021018
.10.1115/1.4038302
32.
Hong
,
S. K.
,
Epureanu
,
B. I.
, and
Castanier
,
M. P.
,
2013
, “
Next-Generation Parametric Reduced-Order Models
,”
Mech. Syst. Signal Process.
,
37
(
1–2
), pp.
403
421
.10.1016/j.ymssp.2012.12.012
33.
Griffin
,
J. H.
, and
Sinha
,
A.
,
1985
, “
The Interaction Between Mistuning and Friction in the Forced Response of Bladed Disk Assemblies
,”
ASME J. Eng. Gas Turbines Power
,
107
(
1
), pp.
205
211
.10.1115/1.3239684
34.
Sinha
,
A.
, and
Griffin
,
J. H.
,
1985
, “
Effects of Friction Dampers on Aerodynamically Unstable Rotor Stages
,”
AIAA J.
,
23
(
2
), pp.
262
270
.10.2514/3.8904
35.
Menq
,
C. H.
,
Bielak
,
J.
, and
Griffin
,
J.
,
1986
, “
The Influence of Microslip on Vibratory Response, Part I: A New Microslip Model
,”
J. Sound Vib.
,
107
(
2
), pp.
279
293
.10.1016/0022-460X(86)90238-5
36.
Menq
,
C. H.
,
Griffin
,
J. H.
, and
Bielak
,
J.
,
1986
, “
The Influence of Microslip on Vibratory Response, Part II: A Comparison With Experimental Results
,”
J. Sound Vib.
,
107
(
2
), pp.
295
307
.10.1016/0022-460X(86)90239-7
37.
Iwan
,
W. D.
,
1966
, “
A Distributed-Element Model for Hysteresis and Its Steady-State Dynamic Response
,”
ASME J. Appl. Mech.
,
33
(
4
), pp.
893
900
.10.1115/1.3625199
38.
Segalman
,
D. J.
, and
Starr
,
M. J.
,
2012
, “
Iwan Models and Their Provenance
,”
ASME
Paper No. DETC2012-71534.10.1115/DETC2012-71534
39.
Bonney
,
M. S.
,
Robertson
,
B. A.
,
Mignolet
,
M.
,
Schempp
,
F.
, and
Brake
,
M. R.
, 2016, “
Experimental Determination of Frictional Interface Models
,”
Dynamics of Coupled Structures
(Conference Proceedings of the Society for Experimental Mechanics Series, Vol. 4), M. Allen, R. Mayes, and D. Rixen, eds., Springer, Cham, Switzerland.
40.
Koh
,
K. H.
,
Griffin
,
J. H.
,
Filippi
,
S.
, and
Akay
,
A.
,
2005
, “
Characterization of Turbine Blade Friction Dampers
,”
ASME J. Eng. Gas Turbines Power
,
127
(
4
), pp.
856
862
.10.1115/1.1926312
41.
Yang
,
B. D.
, and
Menq
,
C. H.
,
1998
, “
Characterization of 3D Contact Kinematics and Prediction of Resonant Response of Structures Having 3D Frictional Constraint
,”
J. Sound Vib.
,
217
(
5
), pp.
909
925
.10.1006/jsvi.1998.1802
42.
Yang
,
B. D.
, and
Menq
,
C. H.
,
1998
, “
Characterization of Contact Kinematics and Application to the Design of Wedge Dampers in Turbomachinery Blading—Part 1: Stick-Slip Contact Kinematics
,”
ASME J. Eng. Gas Turbines Power
,
120
(
2
), pp.
410
417
.10.1115/1.2818138
43.
Yang
,
B. D.
, and
Menq
,
C. H.
,
1998
, “
Characterization of Contact Kinematics and Application to the Design of Wedge Dampers in Turbomachinery Blading—Part 2: Prediction of Forced Response and Experimental Verification
,”
ASME J. Eng. Gas Turbines Power
,
120
(
2
), pp.
418
423
.10.1115/1.2818139
44.
Cigeroglu
,
E.
,
An
,
N.
, and
Menq
,
C. H.
,
2007
, “
A Microslip Friction Model With Normal Load Variation Induced by Normal Motion
,”
Nonlinear Dyn.
,
50
(
3
), pp.
609
626
.10.1007/s11071-006-9171-4
45.
Cardona
,
A.
,
Coune
,
T.
,
Lerusse
,
A.
, and
,
M.
,
1994
, “
A Multiharmonic Method for Nonlinear Vibration Analysis
,”
Int. J. Numer. Methods Eng.
,
37
(
9
), pp.
1593
1608
.10.1002/nme.1620370911
46.
Firrone
,
C. M.
, and
Zucca
,
S.
,
2011
,
Modelling Friction Contacts in Structural Dynamics and Its Application to Turbine Bladed Disks, Numerical Analysis—Theory and Application
,
J.
Awrejcewicz
, ed.,
Intech
,
Rijeka
, Croatia, pp.
301
334
.
47.
Gastaldi
,
C.
, and
Berruti
,
T. M.
,
2017
, “
A Method to Solve the Efficiency-Accuracy Trade-Off of Multi-Harmonic Balance Calculation of Structures With Friction Contacts
,”
Int. J. Non-Linear Mech.
,
92
, pp.
25
40
.10.1016/j.ijnonlinmec.2017.03.010
48.
Cameron
,
T. M.
, and
Griffin
,
J. H.
,
1989
, “
An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems
,”
ASME J. Appl. Mech.
,
56
(
1
), pp.
149
154
.10.1115/1.3176036
49.
Poudou
,
O.
, and
Pierre
,
C.
,
2003
, “
Hybrid Frequency-Time Domain Methods for the Analysis of Complex Structural Systems With Dry Friction Damping
,”
AIAA
Paper No. 2003-1411.10.2514/6.2003-1411
50.
Siewert
,
C.
,
Panning
,
L.
,
Wallaschek
,
J.
, and
Richter
,
C.
,
2010
, “
Multiharmonic Forced Response Analysis of a Turbine Blading Coupled by Nonlinear Contact Forces
,”
ASME J. Eng. Gas Turbines Power
132
(
8
), p. 082501.10.1115/1.4000266
51.
Zucca
,
S.
, and
Firrone
,
C. M.
,
2014
, “
Nonlinear Dynamics of Mechanical Systems With Friction Contacts: Coupled Static and Dynamic Multi-Harmonic Balance Method and Multiple Solutions
,”
J. Sound Vib.
,
333
(
3
), pp.
916
926
.10.1016/j.jsv.2013.09.032
52.
Petrov
,
E. P.
, and
Ewins
,
D. J.
,
2005
, “
Method for Analysis of Nonlinear Multiharmonic Vibrations of Mistuned Bladed Disks With Scatter of Contact Interface Characteristics
,”
ASME J. Turbomach.
,
127
(
1
), pp.
128
136
.10.1115/1.1812781
53.
Petrov
,
E. P.
,
2004
, “
A Method for Use of Cyclic Symmetry Properties on Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks
,”
ASME J. Turbomach.
,
126
(
1
), pp.
175
183
.10.1115/1.1644558
54.
Frey
,
C.
,
Ashcroft
,
G.
,
Kersken
,
H.-P.
, and
Voigt
,
C.
,
2014
, “
A Harmonic Balance Technique for Multistage Turbomachinery Applications
,”
ASME
Paper No. GT2014-25230.10.1115/GT2014-25230
55.
Forster
,
A.
, and
Krack
,
M.
,
2016
, “
An Efficient Method for Approximating Resonance Curves of Weakly-Damped Nonlinear Mechanical Systems
,”
Comput. Struct.
,
169
, pp.
81
90
.10.1016/j.compstruc.2016.03.003
56.
Cigeroglu
,
E.
,
An
,
N.
, and
Menq
,
C. H.
,
2007
, “
Wedge Damper Modeling and Forced Response Prediction of Frictionally Constrained Blades
,”
ASME
Paper No. GT2007-27963.10.1115/GT2007-27963
57.
Cigeroglu
,
E.
,
An
,
N.
, and
Menq
,
C. H.
,
2009
, “
Forced Response Prediction of Constrained and Unconstrained Structures Coupled Through Frictional Contacts
,”
ASME J. Eng. Gas Turbines Power
,
131
(
2
), p. 022505.10.1115/1.2940356
58.
Cigeroglu
,
E.
, and
Ozguven
,
H. N.
,
2006
, “
Nonlinear Vibration Analysis of Bladed Disks With Dry Friction Dampers
,”
J. Sound Vib.
,
295
(
3–5
), pp.
1028
1043
.10.1016/j.jsv.2006.02.009
59.
Grolet
,
A.
, and
Thouverez
,
F.
,
2012
, “
On a New Harmonic Selection Technique for Harmonic Balance Method
,”
Mech. Syst. Signal Process.
,
30
, pp.
43
60
.10.1016/j.ymssp.2012.01.024
60.
Grolet
,
A.
, and
Thouverez
,
F.
,
2015
, “
Computing Multiple Periodic Solutions of Nonlinear Vibration Problems Using the Harmonic Balance Method and Groebner Bases
,”
Mech. Syst. Signal Process.
,
52–53
, pp.
529
547
.10.1016/j.ymssp.2014.07.015
61.
Petrov
,
E. P.
, and
Ewins
,
D. J.
,
2005
, “
Effects of Mistuning on the Forced Response of Bladed Discs With Friction Dampers
,”
Evaluation, Control and Prevention of High Cycle Fatigue in Gas Turbine Engines for Land, Sea and Air Vehicles
, RTO, Neuilly-sur-Seine, France, p. 38.
62.
Tang
,
W. H.
, and
Epureanu
,
B. I.
,
2017
, “
Nonlinear Dynamics of Mistuned Bladed Disks With Ring Dampers
,”
Int. J. Non-Linear Mech.
,
97
, pp.
30
40
.10.1016/j.ijnonlinmec.2017.08.001
63.
Rosenberg
,
R. M.
,
1962
, “
The Normal Modes of Nonlinear n-Degree-of-Freedom Systems
,”
ASME J. Appl. Mech.
,
29
(
1
), pp.
7
14
.10.1115/1.3636501
64.
Grolet
,
A.
, and
Thouverez
,
F.
,
2012
, “
Free and Forced Vibration Analysis of a Nonlinear System With Cyclic Symmetry: Application to a Simplified Model
,”
J. Sound Vib.
,
331
(
12
), pp.
2911
2928
.10.1016/j.jsv.2012.02.008
65.
Shaw
,
S.
, and
Pierre
,
C.
,
1993
, “
Normal Modes for Non-Linear Vibratory Systems
,”
J. Sound Vib.
,
164
(
1
), pp.
85
124
.10.1006/jsvi.1993.1198
66.
Vakakis
,
A.
,
1997
, “
Non-Linear Normal Modes (NNMs) and Their Applications in Vibration Theory: An Overview
,”
Mech. Syst. Signal Process.
,
11
(
1
), pp.
3
22
.10.1006/mssp.1996.9999
67.
Kerschen
,
G.
,
Peeters
,
M.
,
Golinval
,
J. C.
, and
Vakakis
,
A. F.
,
2009
, “
Nonlinear Normal Modes—Part I: A Useful Framework for the Structural Dynamicist
,”
Mech. Syst. Signal Process.
,
23
(
1
), pp.
170
194
.10.1016/j.ymssp.2008.04.002
68.
Laxalde
,
D.
,
Thouverez
,
F.
,
Sinou
,
J.
, and
Lombard
,
J.
,
2007
, “
Qualitative Analysis of Forced Response of Blisks With Friction Ring Dampers
,”
Eur. J. Mech. A/Solids
,
26
(
4
), pp.
676
687
.10.1016/j.euromechsol.2006.10.002
69.
Joannin
,
C.
,
Chouvion
,
B.
,
Thouverez
,
F.
,
Mbaye
,
M.
, and
Ousty
,
J. P.
,
2016
, “
Nonlinear Modal Analysis of Mistuned Periodic Structures Subjected to Dry Friction
,”
ASME J. Eng. Gas Turbines Power
,
138
(
7
), p. 072504.10.1115/1.4031886
70.
Pesheck
,
E.
,
Pierre
,
C.
, and
Shaw
,
S.
,
2002
, “
A New Galerkin-Based Approach for Accurate Non-Linear Normal Modes Through Invariant Manifolds
,”
J. Sound Vib.
,
249
(
5
), pp.
971
993
.10.1006/jsvi.2001.3914
71.
Jiang
,
D.
,
Pierre
,
C.
, and
Shaw
,
S.
,
2005
, “
Nonlinear Normal Modes for Vibratory Systems Under Harmonic Excitation
,”
J. Sound Vib.
,
288
(
4–5
), pp.
791
812
.10.1016/j.jsv.2005.01.009
72.
Shaw
,
S.
, and
Pierre
,
C.
,
1991
, “
Non-Linear Normal Modes and Invariant Manifolds
,”
J. Sound Vib.
,
150
(
1
), pp.
170
173
.10.1016/0022-460X(91)90412-D
73.
Pierre
,
C.
,
Jiang
,
D.
, and
Shaw
,
S.
,
2006
, “
Nonlinear Normal Modes and Their Application in Structural Dynamics
,”
Math. Probl. Eng.
,
2006
, p. 10847.10.1155/MPE/2006/10847
74.
Krack
,
M.
,
Panning-von Scheidt
,
L.
,
Wallaschek
,
J.
,
Siewert
,
C.
, and
Hartung
,
A.
,
2013
, “
Reduced Order Modeling Based on Complex Nonlinear Modal Analysis and Its Application to Bladed Disks With Shroud Contact
,”
ASME J. Eng. Gas Turbines Power
,
135
(
10
), p.
102502
.10.1115/1.4025002
75.
Peeters
,
M.
,
Viguie
,
R.
,
Serandour
,
G.
,
Kerschen
,
G.
, and
Golinval
,
J. C.
,
2009
, “
Nonlinear Normal Modes—Part II: Toward a Practical Computation Using Numerical Continuation Techniques
,”
Mech. Syst. Signal Process.
,
23
(
1
), pp.
195
216
.10.1016/j.ymssp.2008.04.003
76.
Pichler
,
F.
, and
Witteveen
,
W.
,
2014
, “
Efficient Reduction Method Based on Trial Vector Derivatives for Nonlinear Transient Analysis of Multilayer Sheet Structures
,”
PAMM
,
14
(
1
), pp.
37
38
.10.1002/pamm.201410011
77.
Krylov
,
N.
, and
Bogolyubov
,
N.
,
1947
,
Introduction to Nonlinear Mechanics
,
Princeton University Press
, Princeton, NJ.
78.
Caughey
,
T. K.
,
1963
, “
Equivalent Linearization Techniques
,”
J. Acoust. Soc. Am.
,
35
(
11
), pp.
1706
1711
.10.1121/1.1918794
79.
Baek
,
S.
, and
Epureanu
,
B. I.
,
2017
, “
Reduced Order Modeling of Bladed Disks With Friction Ring Dampers
,”
ASME J. Vib. Acoust.
,
139
(6), p. 061011.10.1115/1.4036952
80.
Mitra
,
M.
,
Zucca
,
S.
, and
Epureanu
,
B. I.
,
2016
, “
Adaptive Microslip Projection for Reduction of Frictional and Contact Nonlinearities in Shrouded Blisks
,”
ASME J. Comput. Nonlinear Dyn.
,
11
(
4
), p.
041016
.10.1115/1.4033003
81.
Lupini
,
A.
, and
Epureanu
,
B. I.
,
2019
, “
A Conditioning Technique for Projection-Based Reduced Order Models
,”
Comput. Methods Appl. Mech. Eng.
,
349
, pp. 251–265.10.1016/j.cma.2019.02.020
82.
Kerschen
,
G.
,
Golinval
,
J. C.
,
Vakakis
,
A. F.
, and
Bergman
,
L. A.
,
2005
, “
The Method of Proper Orthogonal Decomposition for Dynamical Characterization and Order Reduction of Mechanical Systems: An Overview
,”
Nonlinear Dyn.
,
41
(
1–3
), pp.
147
169
.10.1007/s11071-005-2803-2
83.
Bui-Thanh
,
T.
, and
Willcox
,
K.
,
2005
, “
Model Reduction for Large-Scale CFD Applications Using Balanced Proper Orthogonal Decomposition
,”
AIAA
Paper No. 2005-4617.10.2514/6.2005-4617
84.
Willcox
,
K.
, and
Peraire
,
J.
,
2002
, “
Balanced Model Reduction Via the Proper Orthogonal Decomposition
,”
AIAA J.
,
40
(
11
), pp.
2323
2330
.10.2514/2.1570
85.
Amsallem
,
D.
,
Zahr
,
M. J.
, and
Farhat
,
C.
,
2012
, “
Nonlinear Model Order Reduction Based on Local Reduced-Order Bases
,”
Int. J. Numer. Methods Eng.
,
92
(
10
), pp.
891
916
.10.1002/nme.4371
86.
Carlberg
,
K.
,
Farhat
,
C.
,
Cortial
,
J.
, and
Amsallem
,
D.
,
2013
, “
The Gnat Method for Nonlinear Model Reduction: Effective Implementation and Application to Computational Fluid Dynamics and Turbulent Flows
,”
J. Comput. Phys.
,
242
, pp.
623
647
.10.1016/j.jcp.2013.02.028
87.
Farhat
,
C.
,
Avery
,
P.
,
Chapman
,
T.
, and
Cortial
,
J.
,
2014
, “
Dimensional Reduction of Nonlinear Finite Element Dynamic Models With Finite Rotations and Energy-Based Mesh Sampling and Weighting for Computational Efficiency
,”
Int. J. Numer. Methods Eng.
,
98
(
9
), pp.
625
662
.10.1002/nme.4668
88.
Farhat
,
C.
,
Chapman
,
T.
, and
Avery
,
P.
,
2015
, “
Structure-Preserving, Stability, and Accuracy Properties of the Energy-Conserving Sampling and Weighting Method for the Hyper Reduction of Nonlinear Finite Element Dynamic Models
,”
Int. J. Numer. Methods Eng.
,
102
(
5
), pp.
1077
1110
.10.1002/nme.4820
89.
Ryckelynck
,
D.
,
2005
, “
A Priori Hyperreduction Method: An Adaptive Approach
,”
J. Comput. Phys.
,
202
(
1
), pp.
346
366
.10.1016/j.jcp.2004.07.015
90.
Chati
,
M.
,
Rand
,
R.
, and
Mukherjee
,
S.
,
1997
, “
Modal Analysis of a Cracked Beam
,”
J. Sound Vib.
,
207
(
2
), pp.
249
270
.10.1006/jsvi.1997.1099
91.
Saito
,
A.
, and
Epureanu
,
B. I.
,
2011
, “
Bilinear Modal Representations for Reduced-Order Modeling of Localized Piecewise-Linear Oscillators
,”
J. Sound Vib.
,
330
(
14
), pp.
3442
3457
.10.1016/j.jsv.2011.02.018
92.
Zucca
,
S.
, and
Epureanu
,
B. I.
,
2014
, “
Bi-Linear Reduced-Order Models of Structures With Friction Intermittent Contacts
,”
Nonlinear Dyn.
,
77
(
3
), pp.
1055
1067
.10.1007/s11071-014-1363-8
93.
Jung
,
C.
,
D'Souza
,
K.
, and
Epureanu
,
B.
,
2012
, “
Bilinear Amplitude Approximation for Piecewise-Linear Oscillators
,”
AIAA
Paper No. 2012-1793.10.2514/6.2012-1793
94.
Jung
,
C.
,
D'Souza
,
K.
, and
Epureanu
,
B. I.
,
2014
, “
Nonlinear Amplitude Approximation for Bilinear Systems
,”
J. Sound Vib.
,
333
(
13
), pp.
2909
2919
.10.1016/j.jsv.2014.01.029
95.
Zucca
,
S.
, and
Epureanu
,
B. I.
,
2017
, “
Reduced Order Models for Nonlinear Dynamic Analysis of Structures With Intermittent Contacts
,”
J. Vib. Control
,
24
(12), 2591–2604.10.1177/1077546316689214
96.
Gastaldi
,
C.
, and
Gola
,
M. M.
,
2016
, “
On the Relevance of a Microslip Contact Model for Under-Platform Dampers
,”
Int. J. Mech. Sci.
,
115–116
, pp.
145
156
.10.1016/j.ijmecsci.2016.06.015
97.
Botto
,
D.
,
,
C.
,
Gola
,
M. M.
, and
Umer
,
M.
,
2017
, “
An Experimental Investigation of the Dynamics of a Blade With Two Under-Platform Dampers
,”
ASME J. Eng. Gas Turbines Power
,
140
(
3
), p.
032504
.10.1115/1.4037865
98.
Battiato
,
G.
,
Firrone
,
C. M.
,
Berruti
,
T. M.
, and
Epureanu
,
B. I.
,
2017
, “
Reduced Order Modeling for Multi-Stage Bladed Disks With Friction Contacts at the Flange Joint
,”
ASME
Paper No. GT2017-64814.10.1115/GT2017-64814
99.
Hudson
,
R.
, and
Sinha
,
A.
,
2016
, “
Frictional Damping of Flutter: Microslip Versus Macroslip
,”
ASME J. Vib. Acoust.
,
138
(
6
), p.
061010
.10.1115/1.4034078
100.
Gastaldi
,
C.
,
Zucca
,
S.
, and
Epureanu
,
B. I.
,
2018
, “
Jacobian Projection Reduced-Order Models for Dynamic Systems With Contact Nonlinearities
,”
Mech. Syst. Signal Process.
,
100
, pp.
550
569
.10.1016/j.ymssp.2017.07.049
101.
Davis
,
P. J.
,
1979
,
Circulant Matrices
, 2nd ed,
Wiley
,
New York
.
102.
Olson
,
B. J.
,
Shaw
,
S. W.
,
Shi
,
C. Z.
,
Pierre
,
C.
, and
Parker
,
R. G.
,
2014
, “
Circulant Matrices and Their Application to Vibration Analysis
,”
ASME Appl. Mech. Rev.
,
66
(
4
), p. 040803.10.1115/1.4027722
103.
Petrov
,
E. P.
,
2010
, “
A Method for Forced Response Analysis of Mistuned Bladed Disks With Aerodynamic Effects Included
,”
ASME J. Eng. Gas Turbines Power
,
132
(
6
), p.
062502
.10.1115/1.4000117
104.
Petrov
,
E. P.
,
,
Z. I.
,
Beretta
,
A.
, and
Elliott
,
R.
,
2013
, “
A Study of Nonlinear Vibrations in a Frictionally Damped Turbine Bladed Disk With Comprehensive Modeling of Aerodynamic Effects
,”
ASME J. Eng. Gas Turbines Power
,
135
(
3
), p.
032504
.10.1115/1.4007871
105.
Wang
,
Y. R.
,
Fu
,
Z. Z.
,
Jiang
,
X. H.
, and
Tian
,
A. M.
,
2015
, “
Mistuning Effects on Aero-Elastic Stability of Axial Compressor Rotor Blades
,”
ASME J. Eng. Gas Turbines Power
,
137
(
10
), p. 102504.10.1115/1.4030280
106.
Nikolic
,
M.
,
Petrov
,
E. P.
, and
Ewins
,
D. J.
,
2008
, “
Robust Strategies for Forced Response Reduction of Bladed Disks Based on Large Mistuning Concept
,”
ASME J. Eng. Gas Turbines Power
,
130
(
2
), p.
022501
.10.1115/1.2799524
107.
Lupini
,
A.
, and
Epureanu
,
B.
,
2018
, “
On the Use of Mesh Morphing Techniques in Reduced Order Models for the Structural Dynamics of Geometrically Mistuned Blisks
,”
15th International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines
(
ISUAAAT-15
),
Oxford, UK
,
Sept. 24–27
, pp. 191–199.
108.
Petrov
,
E. P.
,
2011
, “
Reduction of Forced Response Levels for Bladed Disks by Mistuning: Overview of the Phenomenon
,”
ASME J. Eng. Gas Turbines Power
,
133
(
7
), p.
072501
.10.1115/1.4002619
109.
Murthy
,
R.
, and
Mignolet
,
M. P.
,
2013
, “
On the Benefits of Intentional Mistuning of Friction Dampers to Reduce the Response of Tuned and Mistuned Bladed Disks
,”
ASME
Paper No. GT2013-95792.10.1115/GT2013-95792
110.
Han
,
Y.
,
Xiao
,
B.
, and
Mignolet
,
M. P.
,
2007
, “
Expedient Estimation of the Maximum Amplification Factor in Damped Mistuned Bladed Disks
,”
ASME
Paper No. GT2007-27353.10.1115/GT2007-27353
111.
Lee
,
S.-Y.
,
Castanier
,
M. P.
, and
Pierre
,
C.
, “
Assessment of Probabilistic Methods for Mistuned Bladed Disk Vibration
,”
AIAA
Paper No. 2005-1990.10.2514/6.2005-1990
112.
,
R.
,
Pierre
,
C.
,
Castanier
,
M. P.
, and
Kruse
,
M. J.
,
2002
, “
Dynamic Response Predictions for a Mistuned Industrial Turbomachinery Rotor Using Reduced-Order Modeling
,”
ASME J. Eng. Gas Turbines Power
,
124
(
2
), pp.
311
324
.10.1115/1.1447236
113.
Bhartiya
,
Y.
, and
Sinha
,
A.
,
2014
, “
Geometric Mistuning Identification of Integrally Bladed Rotors Using Modified Modal Domain Analysis
,”
ASME J. Eng. Gas Turbines Power
,
136
(
12
), p. 122504.10.1115/1.4027762
114.
,
A. C.
,
Castanier
,
M. P.
, and
Epureanu
,
B. I.
,
2011
, “
Mistuning Identification of Blisks at Higher Frequencies
,”
AIAA J.
,
49
(
6
), pp.
1299
1302
.10.2514/1.J050427
115.
,
A. C.
,
Castanier
,
M. P.
, and
Epureanu
,
B. I.
,
2008
, “
Reduced-Order Model Construction Procedure for Robust Mistuning Identification of Blisks
,”
AIAA J.
,
46
(
11
), pp.
2890
2898
.10.2514/1.37314
116.
Pichot
,
F.
,
Laxalde
,
D.
,
Sinou
,
J. J.
,
Thouverez
,
F.
, and
Lombard
,
J. P.
,
2006
, “
Mistuning Identification for Industrial Blisks Based on the Best Achievable Eigenvector
,”
Comput. Struct.
,
84
(
29–30
), pp.
2033
2049
.10.1016/j.compstruc.2006.08.022
117.
Murthy
,
R.
, and
Mignolet
,
M. P.
,
2012
, “
Decreasing Bladed Disk Response With Dampers on a Few Blades—Part I: Optimization Algorithms and Blade-Only Dampers Applications
,”
ASME
Paper No. GT2012-69789.10.1115/GT2012-69789
118.
Murthy
,
R.
, and
Mignolet
,
M. P.
,
2012
, “
,”
ASME
Paper No. GT2012-69797.10.1115/GT2012-69797
119.
Zucca
,
S.
, and
Epureanu
,
B. I.
,
2014
, “
,”
International Conference on Noise and Vibration Engineering
(
ISMA
) and International Conference on Uncertainty in Structural Dynamics (USD),
Leuven, Belgium
,
Sept. 15–17
, pp.
2851
2865
120.
Jung
,
C.
,
Saito
,
A.
, and
Epureanu
,
B. I.
,
2012
, “
Detection of Cracks in Mistuned Bladed Disks Using Reduced-Order Models and Vibration Data
,”
ASME J. Vib. Acoust.
,
134
(
6
), p. 061010.10.1115/1.4007244
121.
Beck
,
J. A.
,
Brown
,
J. M.
,
Scott-Emuakpor
,
O. E.
,
Cross
,
C. J.
, and
Slater
,
J. C.
,
2015
, “
Dynamic Response Characteristics of Dual Flow-Path Integrally Bladed Rotors
,”
J. Sound Vib.
,
336
, pp.
150
163
.10.1016/j.jsv.2014.10.011
122.
Hohl
,
A.
,
Siewert
,
C.
,
Panning
,
L.
, and
Wallaschek
,
J.
,
2010
, “
A Substructure Based Reduced Order Model for Mistuned Bladed Disks
,”
ASME
Paper No. DETC2009-87459.10.1115/DETC2009-87459
123.
Salas
,
M. G.
,
,
R.
,
Martensson
,
H.
,
Petrie-Repar
,
P.
,
Fransson
,
T.
, and
Vogt
,
D. M.
,
2017
, “
Forced Response Analysis of a Mistuned, Compressor Blisk Comparing Three Different Reduced Order Model Approaches
,”
ASME J. Eng. Gas Turbines Power
,
139
(
6
), p.
062501
.10.1115/1.4035209
124.
Saito
,
A.
,
Castanier
,
M. P.
, and
Pierre
,
C.
,
2009
, “
Effects of a Cracked Blade on Mistuned Turbine Engine Rotor Vibration
,”
ASME J. Vib. Acoust.
,
131
(
6
), p. 061006.10.1115/1.4000458
125.
Friswell
,
M. I.
,
Penny
,
J. E. T.
, and
Garvey
,
S. D.
,
1995
, “
Using Linear-Model Reduction to Investigate the Dynamics of Structures With Local Nonlinearities
,”
Mech. Syst. Signal Process.
,
9
(
3
), pp.
317
328
.10.1006/mssp.1995.0026
126.
Friswell
,
M. I.
,
Penny
,
J. E. T.
, and
Garvey
,
S. D.
,
1996
, “
The Application of the Irs and Balanced Realization Methods to Obtain Reduced Models of Structures With Local Non-Linearities
,”
J. Sound Vib.
,
196
(
4
), pp.
453
468
.10.1006/jsvi.1996.0495
127.
Irons
,
B.
,
1965
, “
Structural Eigenvalue Problems—Elimination of Unwanted Variables
,”
AIAA J.
,
3
(
5
), pp. 961–962.10.2514/3.3027
128.
Han
,
Y.
, and
Mignolet
,
M. P.
,
2015
, “
A Novel Perturbation-Based Approach for the Prediction of the Forced Response of Damped Mistuned Bladed Disks
,”
ASME J. Vib. Acoust.
,
137
(
4
), p. 041008.10.1115/1.4029946
129.
Gan
,
Y. Q.
,
Mayer
,
J. L.
,
D'Souza
,
K. X.
, and
Epureanu
,
B. I.
,
2017
, “
A Mode-Accelerated Xxr (Max) Method for Complex Structures With Large Blends
,”
Mech. Syst. Signal Process.
,
93
, pp.
1
15
.10.1016/j.ymssp.2017.01.042
130.
Tang
,
W.
,
Epureanu
,
B.
, and
Filippi
,
S.
,
2017
, “
Models for Blisks With Large Blends and Small Mistuning
,”
Mech. Syst. Signal Process.
,
87
(
3
), pp.
161
179
.10.1016/j.ymssp.2016.10.019
131.
,
A.
,
Epureanu
,
B. I.
, and
Filippi
,
S.
,
2012
, “
Reduced-Order Modeling Approach for Blisks With Large Mass, Stiffness, and Geometric Mistuning
,”
AIAA J.
,
50
(
2
), pp.
366
374
.10.2514/1.J051140
132.
Bhartiya
,
Y.
, and
Sinha
,
A.
,
2013
, “
Reduced Order Modeling of a Bladed Rotor With Geometric Mistuning Via Estimated Deviations in Mass and Stiffness Matrices
,”
ASME J. Eng. Gas Turbines Power
,
135
(
5
), p. 052501.10.1115/1.4007783
133.
Jacobs
,
T. D. B.
, and
Martini
,
A.
,
2017
, “
Measuring and Understanding Contact Area at the Nanoscale: A Review
,”
ASME Appl. Mech. Rev.
,
69
(
6
), p.
060802
.10.1115/1.4038130
134.
Ibrahim
,
R. A.
,
1994
, “
Friction-Induced Vibration, Chatter, Squeal, and Chaos Part i: Mechanics of Contact and Friction
,”
ASME Appl. Mech. Rev.
,
47
(
7
), pp.
209
226
.10.1115/1.3111079
135.
Bayram
,
Y.
, and
Nied
,
H.
,
2000
, “
Enriched Finite Element-Penalty Function Method for Modeling Interface Cracks With Contact
,”
Eng. Fract. Mech.
,
65
(
5
), pp.
541
557
.10.1016/S0013-7944(99)00134-4
136.
Wriggers
,
P.
,
2006
,
Computational Contact Mechanics
,
Springer
,
Berlin
.
137.
Belytschko
,
T.
,
Liu
,
W.
,
Moran
,
B.
, and
Elkhodary
,
K.
,
2013
,
Nonlinear Finite Elements for Continua and Structures
,
Wiley
, Chichester, UK.
138.
Elguedj
,
T.
,
Gravouil
,
A.
, and
Combescure
,
A.
,
2007
, “
A Mixed Augmented Lagrangian-Extended Finite Element Method for Modelling Elasticplastic Fatigue Crack Growth With Unilateral Contact
,”
Int. J. Numer. Methods Eng.
,
71
(
13
), pp.
1569
1597
.10.1002/nme.2002
139.
Simo
,
J.
, and
Laursen
,
T.
,
1992
, “
An Augmented Lagrangian Treatment of Contact Problems Involving Friction
,”
Comput. Struct.
,
42
(
1
), pp.
97
116
.10.1016/0045-7949(92)90540-G
140.
Petrov
,
E. P.
, and
Ewins
,
D. J.
,
2003
, “
Analytical Formulation of Friction Interface Elements for Analysis of Nonlinear Multi-Harmonic Vibrations of Bladed Disks
,”
ASME J. Turbomach.
,
125
(
2
), pp.
364
371
.10.1115/1.1539868
141.
Petrov
,
E. P.
, and
Ewins
,
D. J.
,
2007
, “
Advanced Modeling of Underplatform Friction Dampers for Analysis of Bladed Disk Vibration
,”
ASME J. Turbomach.
,
129
(
1
), pp.
143
150
.10.1115/1.2372775
142.
Petrov
,
E. P.
,
2016
, “
Analysis of Bifurcations in Multiharmonic Analysis of Nonlinear Forced Vibrations of Gas Turbine Engine Structures With Friction and Gaps
,”
ASME J. Eng. Gas Turbines Power
,
138
(
10
), p.
102502
.10.1115/1.4032906
143.
Herzog
,
A.
,
Krack
,
M.
,
Panning-von Scheidt
,
L.
, and
Wallaschek
,
J.
,
2014
, “
Comparison of Two Widely-Used Frequency-Time Domain Contact Models for the Vibration Simulation of Shrouded Turbine Blades
,”
ASME
Paper No. GT2014-26226.10.1115/GT2014-26226
144.
Laxalde
,
D.
, and
Thouverez
,
F.
,
2009
, “
Complex Non-Linear Modal Analysis for Mechanical Systems: Application to Turbomachinery Bladings With Friction Interfaces
,”
J. Sound Vib.
,
322
(
4–5
), pp.
1009
1025
.10.1016/j.jsv.2008.11.044
145.
Laxalde
,
D.
, and
Thouverez
,
F.
,
2008
, “
Non-Linear Vibrations of Multi-Stage Bladed Disks Systems With Friction Ring Dampers
,”
ASME
Paper No. DETC2007-34473.10.1115/DETC2007-34473
146.
Afzal
,
M.
,
Arteaga
,
I. L.
, and
Kari
,
L.
,
2016
, “
An Analytical Calculation of the Jacobian Matrix for 3D Friction Contact Model Applied to Turbine Blade Shroud Contact
,”
Comput. Struct.
,
177
, pp.
204
217
.10.1016/j.compstruc.2016.08.014
147.
Saito
,
A.
,
Epureanu
,
B. I.
,
Castanier
,
M. P.
, and
Pierre
,
C.
,
2010
, “
Node Sampling for Nonlinear Vibration Analysis of Structures With Intermittent Contact
,”
AIAA J.
,
48
(
9
), pp.
1903
1915
.10.2514/1.J050061
148.
Petrov
,
E. P.
,
2011
, “
A High-Accuracy Model Reduction for Analysis of Nonlinear Vibrations in Structures With Contact Interfaces
,”
ASME J. Eng. Gas Turbines Power
,
133
(
10
), p.
102503
.10.1115/1.4002810
149.
Griffin
,
J. H.
, and
Menq
,
C. H.
,
1991
, “
Friction Damping of Circular Motion and Its Implications to Vibration Control
,”
ASME J. Vib. Acoust.
,
113
(
2
), pp.
225
229
.10.1115/1.2930173
150.
Menq
,
C. H.
,
Chidamparam
,
P.
, and
Griffin
,
J. H.
,
1991
, “
Friction Damping of 2-Dimensional Motion and Its Application in Vibration Control
,”
J. Sound Vib.
,
144
(
3
), pp.
427
447
.10.1016/0022-460X(91)90562-X
151.
Menq
,
C. H.
, and
Yang
,
B. D.
,
1998
, “
Non-Linear Spring Resistance and Friction Damping of Frictional Constraint Having Two-Dimensional Motion
,”
J. Sound Vib.
,
217
(
1
), pp.
127
143
.10.1006/jsvi.1998.1739
152.
Sanliturk
,
K. Y.
, and
Ewins
,
D. J.
,
1996
, “
Modelling Two-Dimensional Friction Contact and Its Application Using Harmonic Balance Method
,”
J. Sound Vib.
,
193
(
2
), pp.
511
523
.10.1006/jsvi.1996.0299
153.
Yang
,
B. D.
,
Chu
,
M. L.
, and
Menq
,
C. H.
,
1998
, “
Stick-Slip-Separation Analysis and Non-Linear Stiffness and Damping Characterization of Friction Contacts Having Variable Normal Load
,”
J. Sound Vib.
,
210
(
4
), pp.
461
481
.10.1006/jsvi.1997.1305
154.
Firrone
,
C. M.
,
Zucca
,
S.
, and
Gola
,
M. M.
,
2011
, “
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.
,
46
(
2
), pp.
363
375
.10.1016/j.ijnonlinmec.2010.10.001
155.
Botto
,
D.
,
Lavella
,
M.
, and
Gola
,
M. M.
,
2012
, “
Measurement of Contact Parameters of Flat on Flat Contact Surfaces at High Temperature
,”
ASME
Paper No. GT2012-69677.10.1115/GT2012-69677
156.
Gondhalekar
,
A. C.
,
Petrov
,
E. P.
, and
Imregun
,
M.
,
2009
, “
Parameters Identification for Nonlinear Dynamic Systems Via Genetic Algorithm Optimization
,”
ASME J. Comput. Nonlinear Dyn.
,
4
(
4
), p.
041002
.10.1115/1.3187213
157.
Battiato
,
G.
,
Firrone
,
C. M.
, and
Berruti
,
T. M.
,
2017
, “
,”
Mech. Syst. Signal Process.
,
85
, pp.
912
926
.10.1016/j.ymssp.2016.09.019
158.
Allara
,
M.
,
2009
, “
A Model for the Characterization of Friction Contacts in Turbine Blades
,”
J. Sound Vib.
,
320
(
3
), pp.
527
544
.10.1016/j.jsv.2008.08.016
159.
Marinescu
,
O.
,
Epureanu
,
B. I.
, and
Banu
,
M.
,
2011
, “
Reduced Order Models of Mistuned Cracked Bladed Disks
,”
ASME J. Vib. Acoust.
,
133
(
5
), p. 051014.10.1115/1.4003940
160.
Nacivet
,
S.
,
Pierre
,
C.
,
Thouverez
,
F.
, and
Jezequel
,
L.
,
2003
, “
A Dynamic Lagrangian Frequency time Method for the Vibration of Dry-Friction-Damped Systems
,”
J. Sound Vib.
,
265
(
1
), pp.
201
219
.10.1016/S0022-460X(02)01447-5
161.
Krack
,
M.
,
Panning
,
L.
, and
Wallaschek
,
J.
,
2016
, “
On the Interaction of Multiple Traveling Wave Modes in the Flutter Vibrations of Friction-Damped Tuned Bladed Disks
,”
ASME J. Eng. Gas Turbines Power
,
139
(
4
), p.
042501
.10.1115/1.4034650
162.
Siewert
,
C.
, and
Stuer
,
H.
,
2016
, “
Transient Forced Response Analysis of Mistuned Steam Turbine Blades During Startup and Coastdown
,”
ASME J. Eng. Gas Turbines Power
,
139
(
1
), p.
012501
.10.1115/1.4034151
163.
Borrajo
,
J. M.
,
Zucca
,
S.
, and
Gola
,
M. M.
,
2006
, “
Analytical Formulation of the Jacobian Matrix for Non-Linear Calculation of the Forced Response of Turbine Blade Assemblies With Wedge Friction Dampers
,”
Int. J. Non-Linear Mech.
,
41
(
10
), pp.
1118
1127
.10.1016/j.ijnonlinmec.2006.11.003
164.
Saito
,
A.
,
Castanier
,
M. P.
,
Pierre
,
C.
, and
Poudou
,
O.
,
2009
, “
Efficient Nonlinear Vibration Analysis of the Forced Response of Rotating Cracked Blades
,”
ASME J. Comput. Nonlinear Dyn.
,
4
(
1
), p. 011005.10.1115/1.3007908
165.
Hohl
,
A.
,
Siewert
,
C.
,
Panning
,
L.
, and
Kayser
,
A.
,
2008
, “
Nonlinear Vibration Analysis of Gas Turbine Bladings With Shroud Coupling
,”
ASME
Paper No. GT2008-50787.10.1115/GT2008-50787
166.
Krack
,
M.
,
Herzog
,
A.
,
Panning-von Scheidt
,
L.
,
Wallaschek
,
J.
,
Siewert
,
C.
, and
Hartung
,
A.
,
2012
, “
Multiharmonic Analysis and Design of Shroud Friction Joints of Bladed Disks Subject to Microslip
,”
ASME
Paper No. DETC2012-70184.10.1115/DETC2012-70184
167.
Gastaldi
,
C.
, and
Gola
,
M. M.
,
2017
, “
Pre-Optimization of Asymmetrical Underplatform Dampers
,”
ASME J. Eng. Gas Turbines Power
,
139
(
1
), p. 012504.10.1115/1.4034191
168.
Gola
,
M. M.
, and
Gastaldi
,
C.
,
2014
, “
Understanding Complexities in Underplatform Damper Mechanics
,”
ASME
Paper No. GT2014-25240.10.1115/GT2014-25240
169.
Berruti
,
T.
,
Firrone
,
C. M.
, and
Gola
,
M. M.
,
2011
, “
A Test Rig for Noncontact Traveling Wave Excitation of a Bladed Disk With Underplatform Dampers
,”
ASME J. Eng. Gas Turbines Power
,
133
(
3
), p.
032502
.10.1115/1.4002100
170.
Schwingshackl
,
C. W.
,
Petrov
,
E. P.
, and
Ewins
,
D. J.
,
2012
, “
Effects of Contact Interface Parameters on Vibration of Turbine Bladed Disks With Underplatform Dampers
,”
ASME J. Eng. Gas Turbines Power
,
134
(
3
), p.
032507
.10.1115/1.4004721
171.
Petrov
,
E. P.
,
2012
, “
Analysis of Flutter-Induced Limit Cycle Oscillations in Gas-Turbine Structures With Friction, Gap, and Other Nonlinear Contact Interfaces
,”
ASME J. Turbomach.
,
134
(
6
), p.
061018
.10.1115/1.4006292
172.
Odofin
,
A.
, and
Epureanu
,
B. I.
,
2019
, “
Frequency-Adaptive Reduced Order Models for Structures With Intermittent Contacts
,”
Nonlinear Dyn.
(in press).
173.
Wu
,
L.
, and
Tiso
,
P.
,
2016
, “
Nonlinear Model Order Reduction for Flexible Multibody Dynamics: A Modal Derivatives Approach
,”
Multibody Syst. Dyn.
,
36
(
4
), pp.
405
425
.10.1007/s11044-015-9476-5
174.
Witteveen
,
W.
, and
Pichler
,
F.
,
2014
, “
Efficient Model Order Reduction for the Dynamics of Nonlinear Multilayer Sheet Structures With Trial Vector Derivatives
,”
Shock Vib.
,
2014
, p.
1
.10.1155/2014/913136
175.
Pichler
,
F.
,
Witteveen
,
W.
, and
Fischer
,
P.
,
2017
, “
A Complete Strategy for Efficient and Accurate Multibody Dynamics of Flexible Structures With Large Lap Joints Considering Contact and Friction
,”
Multibody Syst. Dyn.
,
40
(
4
), pp.
407
436
.10.1007/s11044-016-9555-2
176.
Loeve
,
M.
,
1948
, “
Fonctions Aleatoires du Second Ordre
,”
Processus Stochastiques et Mouvements Browniens
,
P.
Levy
ed.,
Gauthier-Villars
,
Paris, France
.
177.
Karhunen
,
K.
,
1947
, “
Uber lineare methoden in der Wahrscheinlichkeitsrechnung
,” Annales Academiae Scientiarum Fennicae. Series A1, Mathematica-Physica,
37
, pp. 1–79.
178.
Obukhov
,
A. M.
,
1954
, “
Statistical Description of Continuous Fields
,”
Trans. Geophys. Int. Acad. Nauk USSR
,
24
, pp.
3
42
.
179.
Kosambi
,
D.
,
1943
, “
Statistics in Function Space
,”
J. Indian Math. Soc.
,
7
, pp.
76
88
.10.1007/978-81-322-3676-4_15
180.
Pougachev
,
V. S.
,
1953
, “
General Theory of the Correlations of Random Functions
,”
,
17
, pp.
1401
1402
.
181.
Feeny
,
B. F.
, and
Kappagantu
,
R.
,
1998
, “
On the Physical Interpretation of Proper Orthogonal Modes in Vibrations
,”
J. Sound Vib.
,
211
(
4
), pp.
607
616
.10.1006/jsvi.1997.1386
182.
Dimarogonas
,
A. D.
,
1996
, “
Vibration of Cracked Structures: A State of the Art Review
,”
Eng. Fract. Mech.
,
55
(
5
), pp.
831
857
.10.1016/0013-7944(94)00175-8
183.
Saito
,
A.
,
Castanier
,
M. P.
, and
Pierre
,
C.
,
2009
, “
Estimation and Veering Analysis of Nonlinear Resonant Frequencies of Cracked Plates
,”
J. Sound Vib.
,
326
(
3–5
), pp.
725
739
.10.1016/j.jsv.2009.05.009
184.
Chen
,
S. L.
, and
Shaw
,
S. W.
,
1996
, “
Normal Modes for Piecewise Linear Vibratory Systems
,”
Nonlinear Dyn.
,
10
(
2
), pp.
135
164
.10.1007/BF00045454
185.
Butcher
,
E. A.
, and
Lu
,
R. D.
,
2007
, “
Order Reduction of Structural Dynamic Systems With Static Piecewise Linear Nonlinearities
,”
Nonlinear Dyn.
,
49
(
3
), pp.
375
399
.10.1007/s11071-006-9129-6
186.
Segalman
,
D. J.
,
2007
, “
Model Reduction of Systems With Localized Nonlinearities
,”
ASME J. Comput. Nonlinear Dyn.
,
2
(
3
), pp.
249
266
.10.1115/1.2727495