Application of Improved Dynamic Substructure Finite Element Model-Based State-Space Techniques in Mistuned Blisks

State Key Laboratory of Reliability and Intelligence of Electrical Equipment, Hebei University of Technology, Tianjin 300401, China School of Mechanical Engineering, Hebei University of Technology, Tianjin 300401, China National Engineering Research Center for Technological Innovation Method and Tool, Hebei University of Technology, Tianjin 300401, China School of Electrical Engineering and Automation, Tiangong University, Tianjin 300387, China


Introduction
The bladed disk is a complex rotational mechanical structure, which is one of the key functional transformation components of aeroengine. The complex mistuning reduces the reliability and durability of bladed disk. The influence of mistuning on the vibration localization gradually becomes one of the key problems to be solved urgently in engineering practices. Therefore, the mistuned bladed disk is extensively investigated by scholars.
Usually, the vibration characteristics of mistuned bladed disk are closely related to the natural frequency and modal shape. For example, the effect of bladed mistuning on the blade-disk coupling regions, modal localization, and amplitude magnification were investigated [1,2]. The modal frequency and modal shape of the mistuned bladed disk were studied [3,4]. Based on the investigation of the modal, many researchers studied the output response, including forced responses and transient responses [5][6][7]. Furthermore, Bonhage et al. [8,9] proposed a semianalytical solution to calculate the transient exceeding amplitudes. They established the lumped mass model; although this model has certain directive effect to theoretical research, it has little application in engineering.
In these studies, the models mainly include lumped mass model, continuous parameter model, and finite element model. An electromechanical lumped mass model was established to study the vibration suppression on a mistuned bladed disk using a biperiodic piezoelectric network [10]. A finite element model was built to initiate the sensitivity analysis; however, the computational efficiency of this model is very low. Thus, a new methodology called reduced-order model (ROM) is developed. For instance, Salas et al. [11] considered bladed disk sector as individual substructure with free-interface approach known as Craig-Chang. The ROM was extended and devoted to the development of small geometric mistuning, which required only sector-level calculations [12]. To prevent damages caused by high cycle fatigue, the ROM technique was used to investigate the nonlinear dynamic responses of bladed disk. The performances of bladed disk were examined via ROMs [13]. In addition, the stiffness and mass mistuning, pristine-rogue-interface modal expansion, displacements and deformations, and stress amplitude of bladed disk were investigated by establishing ROM [14][15][16][17]. These scholars have made positive contributions to the development of ROM, but the sector is taken as the research object. Although the computational efficiency is improved, energy cannot be transferred in the whole circumference direction due to the existence of mistuning; therefore, it is not reasonable to take the sector as the research object, and it is very necessary to take the integrally bladed disk (blisk) as the research object.
The design of the mistuned bladed disk is getting thinner and thinner of modern aeroengine, the coupling between blades and disk is getting stronger, the modal distributions are denser, and it is very sensitive to the small mistuning. In addition, the localization caused by mistuning will aggravate the vibration of sectors, which makes a few blades appear higher fatigue stress and even fracture. Therefore, the single blade or sector is regarded as the research object using Campbell diagram to avoid resonance is not suitable; the mistuned blisk is regarded as the research object is necessary to study its vibration response characteristics. Furthermore, it will take a lot of time, and the requirement of the computer configuration is very high if the high-fidelity finite element model method (HFFEMM) is directly adopted. In order to improve the computational efficiency in the condition of ensuring the computational accuracy, a methodology termed as improved dynamic substructure finite element model method based on state-space technique (IDSFEM-SST) is proposed to analyze the vibration characteristics of the mistuned blisk for aeroengine, which is of great significance to theoretical research and engineering application.

Forced Vibration Equation-Based State-Space Techniques
State-space techniques are not only used in modern control theory but also can be used in vibration system.

Modal Equation.
Assuming that the blisk with stiffness and mass mistuning rotates in the airflow under the action of simply harmonic excitation force, the steady-state response equation is expressed as where K∧ 0 and M∧ 0 are, respectively, the stiffness and mass matrix of tuned blisk; ΔK and ΔM are the changes of the stiffness and mass matrix caused by mistuning, which are called mistuned stiffness and mass matrices;Ĉ andĜ are, respectively, the damping and gyro matrix; u is the vibration amplitude; f̂e is the given excitation force; and f̂m is the aerodynamic force associated with vibration. In fact, Equation (1) is a nonlinear equation involving aerodynamic force and damping, which is very complicated. But the main purpose of this work is to verify the effectiveness of the proposed method; the aerodynamic and mechanical forces are considered periodically harmonic excitation force; and the damping is regarded as a constant in order to the convenience of research; thus, Equation (1) is simplified as a linear equation in the application in Sections 3 and 4.
The characteristic equation of the tuned blisk is written as where Φ∧ 0 is the modal shape matrix; Λ 0 is the eigenvalue matrix. They are described as where N is the number of degrees of freedom (DOFs) of tuned blisk; φ 0 i and λ 0 i are, respectively, the modal shape vector and eigenvalue of tuned blisk.
The corresponding modal subset is selected as the modal change matrix Φ∧ 0 , the amplitude vector u is expressed as and where n is the number of nominal modes of selected modal set; d is the number of low-order truncated modes; and α j is the corresponding modal coordinates. The aeroelastic load related to vibration is represented as International Journal of Aerospace Engineering where pðφ 0 j Þ is the unstable aeroelastic load vector caused by the vibration amplitude.
Substituting Equation (5) and Equation (8) into Equation (1), the reduced order modal equation of mistuned blisk can be obtained, which is rewritten as where K 0 and M 0 are, respectively, the modal stiffness and mass matrix of tuned blisk, and they are diagonal matrices; ΔK and ΔM are the changes of modal stiffness and mass matrix caused by mistuning, which are called mistuned modal stiffness and mass matrices; C and G are, respectively, the modal damping and gyro matrix; u is the vibration amplitude;f e is the modal excitation force; Z m is the reduced order modal aerodynamic impedance matrix. They are expressed as where K m , C m , and M m are, respectively, the modal aerodynamic stiffness, damping, and mass matrix; Φ 0H is the Hermitian of Φ 0 .

Forced Vibration Equation.
State variables can completely describe the motion of a system. The present state variables can completely determine the motion state of system at future time if the external input of system is known. The relationship between internal state variables and external input and output variables can be established by state variable description. Thus, the state space equation of Equation (10) can be described as where The state space equation of mistuned blisk, i.e., Equation (17), is described by the reduced order modal equation of tuned blisk. Although D and B are the full matrices, their orders are reduced and the DOFs of model can be declined, which simplifies the calculation process and significantly improves the computational efficiency.
The modal eigenvalue equation of Equation (17) can be obtained when q = 0, which is expressed as where r j is the right eigenvector of mistuned blisk state space; λ j is the jth order characteristic value, and it is expressed as Equation (23) is a complex number consisting of the real component and imaginary component. The jth order mistuning modal has positive damping, and the vibration is stable when real component λ R,j < 0; the jth order mistuning modal has negative damping, and the blisk generates flutter due to the transient response increasing with the rise of the time when real component λ R,j > 0.
Assuming that all λ R,j < 0, the forced vibration response of Equation (17) can be expressed as a linear combination of the state space feature vectors r j , which is written as where β j is the mistuned modal amplitude caused by harmonic vibration. Substituting Equation (24) into Equation (17), the β j is decoupled using the mode orthogonality and represented as where l j is the jth order left eigenvector of mistuned blisk state space of Equation (22), and 3 International Journal of Aerospace Engineering Then, the u in physical coordinates can be obtained combining Equation (5), Equation (17), Equation (21), and Equation (25), which is rewritten as where r dj is the displacement for the eigenvectors of state space r j . The state space techniques do not increase the complexity of system with the rise of state variables, input variables, and output variables. This method can reveal the relationship between the internal and external variables of system; thus, many important features of system, which have not been recognized in the past, can be found. It indicates that the description of state variables to the input and output is more comprehensive than that of transfer functions from the perspective of system structure.

Modal Analysis of Mistuned Blisk
3.1. Modeling. Mistuned blisk modeling is the foundation for investigating the modal and response. According to the actual situation, the accidents or failures of aeroengine commonly occur in blades and less in disk; therefore, the blade is regarded as mistuned substructure, and disk is regarded as tuned substructure. The model parameters are given at first, which are listed in Table 1. The variations of blade stiffness and density are shown in Figures 1 and 2.
The finite element model of mistuned blisk is established by IDSFEM-SST, which is shown in Figure 3.
As seen in Figure 3, the disk is regarded as the superelement, and its internal DOFs are reduced, the contact DOFs between blades and disk are regarded as constraint modal, and the internal DOFs of blades are regarded as the main modal. Figure 3 indicates there are only 128208 units and 254712 nodes in improved dynamic substructure FEM, but there are 2796984 units and 4218072 nodes in HFFEMM. The number of units and nodes of the proposed model is only 4.58% and 6.04% of HFFEMM. So, the computational time of IDSFEM-SST is cheaper than that of HFFEMM.

Verification of Effectiveness.
To verify the rationality of IDSFEM-SST, the first 40 modal frequencies of mistuned blisk with w = 1060 rad · s -1 , t = 1050°C are analyzed by this method compared with that of by classical dynamic substructure finite element model method (CDSFEMM) and HFFEMM in Figure 4. The first 10, 20, 30, and 40 modal frequencies are, respectively, calculated by IDSFEM-SST, CDSFEMM, and HFFEMM in Figure 5 and Table 2 to study the influence of different order frequencies on the computational efficiency.
As seen in Figure 4, the errors of modal frequencies of mistuned blisk are, respectively, 0.0035%~0.1146% and 0.0047%~0.1241% with IDSFEM-SST and CDSFEMM relative to that of HFFEMM which is regarded as a criterion. The modal frequencies obtained by IDSFEM-SST and CDSFEMM are almost consistent with that of HFFEMM. Thus, the computational accuracy of IDSFEM-SST and CDSFEMM is comparable to that of HFFEMM. Furthermore, the computational accuracy of IDSFEM-SST is higher than that of CDSFEMM. Figure 5 and Table 2 indicate that the computational time of mistuned blisk is, respectively, reduced by 18.78%~23.98% and 19.06%~30.74% via CDSFEMM and IDSFEM-SST compared with that of via HFFEMM. It is obvious that the computational efficiency of IDSFEM-SST is higher than that of CDSFEMM. Obviously, the computational efficiency and accuracy of IDSFEM-SST are superior to those of CDSFEMM. Therefore, the proposed IDSFEM-SST is reasonable.

Natural Frequency Analysis
3.3.1. Influence of Stiffness on Natural Frequency for Mistuned Blisk. The centrifugal force generated by rotational speed causes blades to subject to the tensile action, and this tensile force has a tendency to take it back to the equilibrium position, which may result in increasing of stiffness, while the high temperature may cause the decrease of elastic coefficient for the material, which results in decreasing of the stiffness. The first 40 order natural frequencies with stiffness mistuning of blades are shown in Figure 6.
It can be seen from Figure 6 that the stiffness mistuning has a greater influence on high-order natural frequencies than that of low-order natural frequencies. Compared with the modal of tuned blisk, the 1st modal decreases by 1.75 Hz, 1.19 Hz, and 0.7 Hz, and the 40th modal decreases by 29.9 Hz, 20.3 Hz, and 11.9 Hz when the stiffness mistuning is -5%, -3%, and -2%. However, the 1st modal increases by 1.68 Hz, 1.00 Hz, and 0.66 Hz, and the 40th modal increases den d is the disk density; e d is the disk elasticity modulus; pr d is the disk Poisson's ratio; α d is the disk thermal expansion coefficient; k d is the disk thermal conductivity. den b is the blade density; e b is the blade elasticity modulus; pr b is the blade Poisson's ratio; α b is the blade thermal expansion coefficient; k b is the blade thermal conductivity; w is the rotational speed of mistuned blisk; t is the temperature of mistuned blisk. 4 International Journal of Aerospace Engineering by 28.4 Hz, 16.9 Hz, and 11.2 Hz when the mistuning stiffness is +5%, +3%, and +2%. This manifests that the decrease of stiffness has a slightly greater impact on the natural frequencies than that of the increase of stiffness, which is unfavorable to the aeroengine.

Influence of Density on Natural Frequency for Mistuned
Blisk. The processing, manufacturing, installation, and wear may cause the uneven of quality for the mistuned blisk; furthermore, the mistuned blisk works in severe environment; thus, it is necessary to study the influence of the density on natural frequencies of mistuned blisk. The first 40 order natural frequencies with the density mistuning of the blade are shown in Figure 7. It can be seen from Figure 7 that the variation tendency of density mistuning is the same as that of stiffness mistuning, but the increase or decrease of density mistuning causes the modal frequency to change, which is contrary to that of stiffness mistuning. In addition, the mistuned blisk becomes lighter, which results in widening of frequency band when the wear of blade is severe, which increases vibration and is unfavorable for aeroengine.    International Journal of Aerospace Engineering

Modal Shape Amplitudes Analysis.
The blisk is very sensitive to mistuning which can trigger vibration localization, high fatigue stress, and fracture. Therefore, it is very necessary to investigate the modal shape under centrifugal force and thermal load. The difference from previous studies (localization factor as research object [39]) is that the maximum modal strain energy amplitudes (MMSEAs) is regarded as the research object so as to make sure the variation tendency of the most dangerous position and reasonably prevent blisk from destroying by mistuning in this work.

Maximum Modal Strain Energy Amplitudes under
Centrifugal Force. Aeroengine works in high rotational speed, and it is very necessary to investigate the effect of centrifugal force on the maximum modal strain energy amplitudes; the influence of centrifugal force on mistuned blisk is a nonnegligible factor. The first 40 order MMSEAs are calculated using IDSFEM-SST in different rotational speeds as shown in Figure 8(a); the deviations of MMSEAs are shown in Figure 8(b).
As seen in Figure 8, the fluctuation of MMSEAs of mistuned blisk in different rotational speeds from the 7th order is drastic, and their variation tendencies are jagged, and their values are larger than those of tuned blisk. However, the variation tendency of MMSEAs is different from the maximum modal displacement and stress amplitudes. The deviations gradually rise with the increase of rotational speeds, and the 36th order MMSEA achieves the maximum, which makes up the insufficient for the investigations of the maximum modal displacement and stress amplitudes [4]. This investigation is very useful to engineering practices because the centrifuge force is the main cause of blisk damaged when aeroengine works at high rotational speed.

Maximum Modal Strain Energy Amplitudes under
Thermal Load. The turbine blisk is also affected by high temperature. The first 40 order MMSEAs and their deviations are calculated using IDSFEM-SST in different temperatures as shown in Figure 9. Figure 9 indicates that the MMSEAs begin to mutate from the 7th order, and oscillation amplitudes gradually increase. Therefore, the sensitivity of high-order strain energy to the temperature significantly rises. In addition, the deviations of MMSEAs for mistuned blisk constantly oscillate and gradually increase. Meanwhile, the deviations of mistuned blisk change smoothly, and values are smaller in the first vibration period than that of in the second vibration period compared with tuned blisk. Therefore, the deviation changes of MMSEAs are different from those of MMDAs and MMSAs.

11
International Journal of Aerospace Engineering

Maximum Modal Strain Energy Amplitudes under
Centrifugal Force and Thermal Load. In fact, the mistuned blisk is subjected to common action of rotational speed and temperature; thus, it is more consistent with the actual situation to investigate the maximum modal strain energy amplitudes in temperature and rotational speed. The first 40 order MMSEAs and their deviations are calculated using IDSFEM-SST in different rotational speeds and temperatures as shown in Figure 10. Figure 10 manifests that the high-order MMSEAs are sensitive to the temperature, and the mistuned blisk is more likely to occur failure when the blisk enters into the second vibration period. Therefore, the sensitivity of high-order MMSEAs to the temperature will heighten when the temperature rises to a certain degree. In addition, the deviations of MMSEAs increase sharply in the first vibration period. Furthermore, the destructive of mistuned blisk increases when the rotational speed is invariant, and the temperature is enhanced to a certain degree. Meanwhile, the 9th order frequency occurs larger change when the temperature is invariant, and rotational speed is changed; thus, MMSEAs are impacted obviously by rotational speed in low-order frequency; large shocks of high order are observed.
It can be seen from Section 3 the investigation of the maximum modal strain energy amplitudes is a beneficial complement for that of the displacement and stress. In fact, the probability distribution of the mode shape can be further studied so as to find out the main factors affecting the mode shape to reasonably control these factors and improve the robustness of the system, as well as carry out prediction research by the inverse probability design.

Vibration Response Analysis of
Mistuned Blisk 4.1. Maximum Vibration Response Analysis. According to the actual situation, the displacement response is the biggest in the radial direction relative to that of in the axial and tangential directions; moreover, the radial deformation has the greatest influence on the colliding of blisk with casing. Therefore, the maximum radial displacement response is investigated. This work mainly studies the small mistuning (A int = 0%-5%), because when the mistuning exceeds 5%, the blisk has been damaged, and there is no need to study it. Thus, assuming that the damping ratio of the system is 0.18%, and the most dangerous point at the front tip of the blade edge is as the extraction point of response amplitude. Only the maximum displacement responses rather than all the displacement responses are investigated, which not only can improve the computational efficiency but also has practical engineering significance. This is different and novel from previous studies. The stiffness mistuning is, respectively, 0%, 2%, 3%, and 5% when w = 1046 rad s −1 , t = 1050°C; the radial maximum displacement responses, i.e., frequency response functions of mistuned blisk, are shown in Figure 11. Figure 11(a) indicates the peak of tuned blisk is caused by resonance when the excitation frequency is 834 Hz, i.e., its excitation frequency equaling or closing to a certain natural frequency. Figures 11(b)-11(d) indicate that the mistuned blisk appears multiple peaks, which are caused by resonance and mistuning, and the number of peaks rises with the increase of mistuning. The main reason is that when the blisk works in high rotational speed, the strain energy cannot be uniformly transferred to each blade, but part of the energy is reflected back and gathered on some blades due to the existence of mistuning so that multiple peaks appear in the blades, and the peak appears earlier than that of tuned blisk.

Verification of Computational Efficiency.
To verify the effectiveness of IDSFEM-SST, the computational time and the saving rate of maximum displacement responses for the mistuned blisk are analyzed compared with that of by HFFEMM and CDSFEMM. The calculation results are shown in Table 3.
It can be seen from Table 3 that the computational time of IDSFEM-SST and CDSFEMM is greatly shortened compared with that of HFFEMM; the computational efficiencies are improved by 43.82%-62.58% and 39.15%-58.32%, respectively. The computational efficiency of IDSFEM-SST is higher 1.74%-10.21% than that of CDSFEMM. Therefore, the IDSFEM-SST is more advantageous; the main reason is that the symmetry of mistuned blisk is broken, and the energy is concentrated in a few blades. However, the mistuned blisk has a lot of DOFs using HFFEMM, and the computational time is greatly extended, but the DOFs of mistuned blisk are significantly reduced with IDSFEM-SST, the computational time is greatly decreased, and the computational accuracy is comparable to that of HFFEMM. Therefore, the IDSFEM-SST is a very applicable and effective method for mistuned blisk analysis.
This investigation indicates the blisk is more likely to be destroyed with the increase of the mistuning level; thus, the research of the mistuning is very important as well as that of the resonance, and reasonable reduction of mistuning is very necessary. Usually, the mistuning level of the blisk can be reduced by controlling the distribution of input parameters or by introducing the intentional mistuning so as to improve the robustness of the blisk and reduce the destruction.

Conclusions and Outlooks
(1) A methodology named IDSFEM-SST is proposed, and the ROM is derived.  Figure 11: Radial frequency response functions of mistuned blisk. A int is the mistuning level; it is tuned blisk when A int = 0%. 13 International Journal of Aerospace Engineering computational time of maximum displacement responses is reduced by 43.82%-62.58% via IDS-FEM-SST, and the computational efficiency is improved 1.74%-10.21% compared with CDSFEMM when the mistuned level is 0%-5%. Thus, the computational efficiency and accuracy by IDSFEM-SST are superior to by CDSFEMM for mistuned blisk (3) It is only investigated that why there are multiple peaks after mistuning of blisk and the number of peaks increases with the rise of mistuning. In addition, it is pointed out the displacement responses are the largest in the radial direction compared with that of in the axial and tangential directions. However, there is not a reasonable method is proposed to control the radial deformation. This is the next step to study. Besides, the influence factors such as temperature, rotational speed, and geometry size are random in working for aeroengine. Therefore, the nondeterministic investigation is further carried out based on the dynamic analysis, i.e., the probability is very necessary. Thus, another work of next step is to study the reliability of the mistuned blisk. The probability analysis is the study of the sensitivity of output results to input variables so as to find out the main factors affecting the mode shape and vibration response to reasonably control these factors or by introducing the intentional mistuning, improve the robustness and reliability of the system, and reduce the destruction

Data Availability
The (data type) data used to support the findings of this study are currently under embargo while the research findings are commercialized. Requests for data (6/12 months after publication of this article) will be considered by the corresponding author.

Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.