Dynamic Modeling and Projection-Based Reduction Methods for Bladed Disks With Nonlinear Frictional and Intermittent Contact Interfaces

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 [4–14]. 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) [17–22], 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) [25–27], and asymptotic model of mistuning [28–30]. 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 [37–40]. Later, models which predict the local frictional forces under different regimes of stick, slip, and separation under varying periodic normal loads were proposed [41–44]. 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) [45–47]. 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 [50–60].

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) [63–67], 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,70–73] 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) [82–84] 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 [85–89] 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 [90–92]. 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,96–99]. 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.

A cyclic symmetric blisk or tuned blisk is comprised of n_max 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.

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.

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 n_max 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 E ⊗ I. 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 = n_max/2 if n_max 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 (n_max/2 or $(nmax−1)/2$ for blisks with even and odd n_max, 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.

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

where e^p 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 = jn_max ± 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.

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.

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,113–116], 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].

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 [28–30] 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.

where M_SS 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 [125–127]. 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 [25–27]. 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,129–132]. 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.

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,140–148] 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.

where k_z,i is the normal contact stiffness, and v_i is the normal relative displacement between the nodes, assumed negative when there is no penetration between the two surfaces. N_i 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 v_i is negative.

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.

$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 [155–157] 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).

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.

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. f_E is the vector of excitation forces usually an EO excitation. f_C is the vector of contact forces.

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 f_C 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.

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.

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,167–170], 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).

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.

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.

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.

Any relevant linear reduction technique, such as those discussed in Sec. 3.2, may be applied to obtain the linear reduction basis $ΦROM,L$ [125–127]. 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 [173–175]. 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.

The proper orthogonal decomposition (POD), also sometimes called the Karhunen–Loeve decomposition, was first proposed independently by a number of sources [176–180]. 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].

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.

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.

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.

where K_PLM 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.

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.

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 α.

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].

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.

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.

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.

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.

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.

Symbols

A =

additional projected stiffness mistuning matrix

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 =

sector or blade index

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

a =

physical global coordinate system

A =

variable for degrees of freedom in matrix partition

AMP =

reduction using adaptive microslip projection

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

h =

harmonic index number

H =

complex conjugate transpose of a matrix

S =

sector level matrix

T =

transpose of a matrix

Abbreviations

AF(s) =

amplification factor(s)

AFT =

alternating frequency-time

AMP(s) =

adaptive microslip projection (modes)

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

$¯$ =

quantity in frequency domain

$˙$ =

single time derivative

$¨$ =

double time derivative

$̂$ =

estimate

$̃$ =

constrained quantity