A new finite element model for vibration analysis of a coated blisk with variable coating thickness

Abstract

In the damping design of blisks with hard coating, changing coating thickness is one of the frequently used ways for achieving better damping performance. Nevertheless, for the vibration analysis with variable coating thickness by using commercial finite element software, it is necessary to regenerate the finite element mesh, which results in inconvenience in computation. In this paper, a new laminated element is developed to solve this problem, which needs solely the node coordinates corresponding to the maximum coating thickness. The additional virtual layers to simulate the variations of coating thickness are introduced for the presented element to avoid regenerating new node coordinates. By comparing the analysis results obtained by SOLID185 in ANSYS software, experiment and the developed element, this new element is verified that it remains effective for variable coating thickness parameter input. Furthermore, the effects of the coating thickness on the natural frequency and the forced response of an actual blisk are analysed and due to using the same node coordinates, the computational efficiency is improved significantly. The analysis results indicate that with the increase of coating thickness, the resonance response is reduced significantly and even the thinner coating could also achieve better damping performance.

Introduction

Because of the blades and the disk being made of one piece, blisks can improve the thrust-to-weight ratio significantly and has been used widely in aeroengine. ¹ Due to lack of the friction damping like traditional bladed disks,^2,3 the blisk is easily resonated under the broadband excitation from aerodynamic loads, mechanical loads and temperature field, resulting in fatigue damage in aeroengine applications.^4,5 Therefore, it is necessary to minimise the excessive vibration of blisk in the design process.

Many studies have reported the approaches to enhance the damping of blades, such as friction dampers, ⁶ piezoelectric materials,^7,8 passive surface treatments ⁹ and active magnetic bearings. ¹⁰ These methods are effective in increasing the resistance of blisks against vibrations. At present, a new passive surface treatment technique using hard-coating damping has been presented due to its perfect performance to improve surface wear resistance by increasing surface hardness^11–13 and reduce the vibration stress of structures in high temperature and corrosion environment.^14,15 Even the coating spraying technology could also have an effect on the coating damping property, and the smart adhesion ¹⁶ could be used to improve its performance when microscopic characteristics of coating are studied. Thus, hard-coating damping technology is applied widely.

Over the past decades, the modelling methods for blisks mainly include lumped-parameter model,^17–21 continuous parameter model,^22,23 and high-fidelity finite element model (HFEM).^24,25 In general, the lumped-parameter model and the continuous parameter model (containing the FEM and analytical model simultaneously) for simplified blisks could not capture the real characteristic of actual blisks ⁹ and are only used for researching the mechanism. In contrast, the HFEM is usually used to obtain the reasonable result in industrial applications. However, an HFEM could contain thousands of degree of freedoms (DOFs), which requires expensive calculation. Therefore, many scholars use model reduction techniques such as SNM (Subset of Nominal Modes), ²⁶ FMM (Fundamental Mistuning Method), ²⁷ CMM (Component Mistuning Method), ²⁸ AMM (Asymptotic Mistuning Model) ²⁹ and N-PRIME (New Rogue Interface Modal Expansion) ³⁰ to obtain the analysis model to study the vibration characteristics of blisks. These Reduced Order Models can improve the efficiency for vibration response analysis significantly.

Due to the high-damping performance of hard coating, many researchers have studied the applications of hard coating on simple structures recently, such as beams, plates and shells.^31–34 A few of studies^35–37 researched the HFEM and the vibration characteristics of hard-coating blisks. These tentative researches are of great significance for reducing the vibration level of blisks and improving their working life. However, in the above-mentioned researches including the modelling method of uncoated blisk, the stiffness and mass matrices were usually obtained from the commercial finite element software such as ANSYS. It should be noted that in the coating damping design of blisks, it is especially essential to vary the coating thickness to improve damping performance. When the coating thickness is changed in commercial finite element software to solve the new stiffness and mass matrices, the finite element mesh of coated blisk needs to be restructured. It is exceedingly inconvenient in the damping design process and will be specifically described in Disadvantages of vibration analysis for coated blisk based on ANSYS section. Therefore, to study the modelling method for blisks with variable coating thicknesses is of great significance.

The aim of this paper is to develop a new FEM of hard-coating blisks. By introducing the additional virtual layers, this new FEM can be used to incorporate variable coating thickness in the damping design process. Some modelling methods for changing the coating thickness on blades are given by using ANSYS, and the disadvantages of using the commercial software to perform the damping design of blisks are also presented in the next section. Then, the FEM of a hard-coating blisk is established based on the presented three-dimensional eight-node incompatible laminated element (named as LI-SOLID element). This is followed by a section in which the dynamic analysis method of the coated blisk is given. In Model verification and application section, the rationality of the created model is verified by ANSYS software and experiment. The advantage of the LI-SOLID element is also demonstrated. Finally, some important conclusions are listed.

Disadvantages of vibration analysis for coated blisk based on ANSYS

The commercial finite element software ANSYS is a powerful tool to dynamics analysis of many structures. However, it could be very inefficient to solve stiffness and mass matrices when the coating thickness is changed in the damping design process of hard-coating blisks. In this section, we only introduce two methods to simulate the coated blisk with variable coating thickness.

For the first method, the laminated element SOLID185 is used to simulate the coated blades as shown in Figure 1. Note that the element number of a single blade is 165, with the nodes number being 384. When the coating thickness is changed, the finite element mesh needs to be manually regenerated. In this way, the nodes coordinates of the coated blades and those of the disk connected to the coated blades must be changed simultaneously as shown in Figure 1.

Figure 1.

Laminated FEM with SOLID185.

For the second method^35,38 shown in Figure 2, the coated blade is divided into three layers elements of SOLID185. Then, the element number of a single blade is 495, with the nodes number being 768. The problem similar to the first method still exists when coating thickness is changed.

Figure 2.

FEM with three layers elements of SOLID185.

Although the element number of the first method is less than the second one, the two methods based on ANSYS software are also computationally expensive for vibration damping design of blisks with coating. There may be other modelling methods we have not introduced, but the problem could be same as the above-mentioned methods when the commercial software is used. So it is essential to develop a composite element for coated blisks to satisfy the requirement of variable coating thickness, which can be conducive to vibration reduction techniques using hard coating.

FEM of coated blisk

In this section, we present a new laminated incompatible solid element, named as LI-SOLID element, to derive the stiffness and mass matrices for the blisk with hard coating on blades as shown in Figure 3(a). This element is similar to SOLID185 in ANSYS, with eight nodes and three DOFs of each node as shown in Figure 3(b), and Figure 3(c) is the element sector of coated blades. However, it does not need to restructure the finite element mesh when coating thickness is changed, which could improve the computational efficiency in the design process with coating.

Figure 3.

FEM of the coated blisk. (a) Coated blisk, (b) LI-SOLID and (c) element sector of coated blade.

Stiffness and mass matrices

The LI-SOLID element introduces nine additional DOFs, with being eliminated by static condensation while calculating the stiffness matrix. Therefore, it could avoid the shear-locking effect of compatible element and not increase the DOFs of the system. The corresponding displacements for an arbitrary point in the LI-SOLID element can be described as

$$\begin{array}{c}u={\displaystyle \sum _{i=1}^8{N}_i{u}_i+}{N}_9{\alpha }_1\text{+}{N}_{10}{\alpha }_2\text{+}{N}_{11}{\alpha }_3\\ v={\displaystyle \sum _{i=1}^8{N}_i{v}_i+}{N}_9{\alpha }_4\text{+}{N}_{10}{\alpha }_5\text{+}{N}_{11}{\alpha }_6\\ w={\displaystyle \sum _{i=1}^8{N}_i{w}_i+}{N}_9{\alpha }_7\text{+}{N}_{10}{\alpha }_8\text{+}{N}_{11}{\alpha }_9\end{array}$$

(1)

where ${\alpha }_j\left(j=1,\dots ,9\right)$ are the introduced nine additional DOFs, and the above various shape functions can be expressed as

$$N_i=\frac{1}{8}\left(1+{\xi }_i\xi \right)\left(1+{\eta }_i\eta \right)\left(1+{\varsigma }_i\varsigma \right),(i=1,2,\dots ,8),N_9=1-{\xi }^2,N_{10}=1-{\eta }^2,N_{11}=1-{\varsigma }^2$$

(2)

For a composite structure with n layers, the block static equilibrium equation can be expressed as

$$\left[\begin{array}{l}{\mathit{K}}_{11}{\mathit{K}}_{12}\hfill \\ {\mathit{K}}_{21}{\mathit{K}}_{22}\hfill \end{array}\right]\left(\begin{array}{c}{\mathit{\delta }}^{\text{e}}\\ \mathit{\alpha }\end{array}\right)=\left(\begin{array}{c}{\mathit{F}}_1\\ \mathbf{0}\end{array}\right)$$

(3)

where ${\mathit{\delta }}^{\text{e}}$ represents the displacement vector corresponding to the main DOFs, with $\mathit{\alpha }$ denoting the displacement vector of the additional DOFs, and ${\mathit{F}}_1$ is the force vector. The stiffness matrices in equation (3) can be computed by

$${\mathit{K}}_{11}^{\hspace{0.17em}}\text{=}{\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle \underset{S}{\iint }{\displaystyle {\int }_{z_m}^{z_{m+1}}{\mathit{B}}_{\hspace{0.17em}}^{\text{T}}{\mathit{D}}_m{\mathit{B}}_{\hspace{0.17em}}\text{ d}z\text{d}S}}}$$

(4a)

$${\mathit{K}}_{12}^{\hspace{0.17em}}\text{=}{\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle {\iint }_S{\displaystyle {\int }_{z_m}^{z_{m+1}}{\mathit{B}}_{\hspace{0.17em}}^{\text{T}}{\mathit{D}}_m{\overline{\mathit{B}}}_{\hspace{0.17em}}\text{ d}z\text{d}S}}}$$

(4b)

$${\mathit{K}}_{22}\text{=}{\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle {\iint }_S{\displaystyle {\int }_{z_m}^{z_{m+1}}{\overline{\mathit{B}}}_{\hspace{0.17em}}^{\text{T}}{\mathit{D}}_m{\overline{\mathit{B}}}_{\hspace{0.17em}}\text{ d}z\text{d}S}}}$$

(4c)

where the elastic matrix ${\mathit{D}}_m$ and the strain matrix ${\mathit{B}}_{\hspace{0.17em}}$ are consistent with the compatible element and $\overline{\mathit{B}}$ is an enhanced strain matrix generated by introducing additional DOFs. The specific expressions of each parameter are

$$\mathit{B}=\left[\begin{array}{cccccccc}{\mathit{B}}_1& {\mathit{B}}_2& {\mathit{B}}_3& {\mathit{B}}_4& {\mathit{B}}_5& {\mathit{B}}_6& {\mathit{B}}_7& {\mathit{B}}_8\end{array}\right],\overline{\mathit{B}}=\left[\begin{array}{ccc}{\mathit{B}}_9& {\mathit{B}}_{10}& {\mathit{B}}_{11}\end{array}\right]$$

(5)

$${\mathit{B}}_i={\left[\begin{array}{cccccc}\hfill \frac{\partial N_i}{\partial x}\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial y}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial z}\hfill \\ \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial y}\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial x}\hfill & \hfill \frac{\partial N_i}{\partial z}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial z}\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial y}\hfill & \hfill \frac{\partial N_i}{\partial x}\hfill \end{array}\right]}^{\text{T}}\begin{array}{cc},& (i=1,\dots ,8)\end{array}$$

(6)

$${\overline{\mathit{B}}}_i={\left[\begin{array}{cccccc}\hfill \frac{\partial N_i}{\partial x}\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial y}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial z}\hfill \\ \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial y}\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial x}\hfill & \hfill \frac{\partial N_i}{\partial z}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial z}\hfill & \hfill 0\hfill & \hfill \frac{\partial N_i}{\partial y}\hfill & \hfill \frac{\partial N_i}{\partial x}\hfill \end{array}\right]}^{\text{T}}\begin{array}{cc},& (i=9,\dots ,11)\end{array}$$

(7)

$${\mathit{D}}_m=\frac{E_m^{\hspace{0.17em}}}{\left(1+{\mu }_m\right)\left(1-2{\mu }_m\right)}\left[\begin{array}{cccccc}\hfill 1-{\mu }_m\hfill & \hfill {\mu }_m\hfill & \hfill {\mu }_m\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill {\mu }_m\hfill & \hfill 1-{\mu }_m\hfill & \hfill {\mu }_m\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill {\mu }_m\hfill & \hfill {\mu }_m\hfill & \hfill 1-{\mu }_m\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \frac{1-2{\mu }_m}{2}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \frac{1-2{\mu }_m}{2}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \frac{1-2{\mu }_m}{2}\hfill \end{array}\right]$$

(8)

where $E_m^{\hspace{0.17em}}$ and ${\mu }_m$ are the Young’s modulus and the Poisson’s ratio for the mth layer, respectively.

Then, with static condensation to eliminate the additional DOFs $\mathit{\alpha }$ , the stiffness matrix of the LI-SOLID element can be given by

$${\mathit{K}}_e={\mathit{K}}_{11}-{\mathit{K}}_{12}{\left({\mathit{K}}_{22}\right)}^{-1}{\mathit{K}}_{21}$$

(9)

It should be noted that the mass matrix of the LI-SOLID element is obtained by the shape function of the compatible one, which is expressed as

$${\mathit{M}}_{\text{e}}={\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle {\iint }_S{\displaystyle {\int }_{z_m}^{z_{m+1}}{\rho }_{\text{m}}{\mathit{N}}^{\text{T}}\mathit{N}\text{ d}z\text{d}S}}}$$

(10)

where the shape function matrix $\mathit{N}$ is

$$\mathit{N}\text{=}\left[\begin{array}{ccccccc}{N}_1& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& N_8& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& N_1& \hspace{0.17em}& \dots \text{.}& \hspace{0.17em}& N_8& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& N_1& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& N_8\end{array}\right]$$

(11)

In order to use the Gauss rule for each layer conveniently in equation (4) and equation (10), the coordinate transformation of the mth layer is employed and is given by

$$\varsigma \text{=}\frac{{\varsigma }_m-{\varsigma }_{m\text{+}1}}{2}{\varsigma }^{\text{v}}+\frac{{\varsigma }_m\text{+}{\varsigma }_{m\text{+}1}}{2}$$

(12)

where ${\varsigma }_m$ and ${\varsigma }_{m\text{+}1}$ are the thickness-direction coordinates for the mth layer in local coordinate $\xi \eta \varsigma$ . Furthermore, it is easily to obtain that

$${\mathit{K}}_{11}^{\hspace{0.17em}}={\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\mathit{B}}_{\hspace{0.17em}}^{\text{T}}{\mathit{D}}_m^{\hspace{0.17em}}\mathit{B}\left|\mathit{J}\right|\frac{{\varsigma }_{m\text{+}1}-{\varsigma }_m}{2}\text{d}\xi \text{d}\eta \text{d}{\varsigma }^{\text{v}}}}}}$$

(13a)

$${\mathit{K}}_{12}^{\hspace{0.17em}}={\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\mathit{B}}_{\hspace{0.17em}}^{\text{T}}{\mathit{D}}_m^{\hspace{0.17em}}\overline{\mathit{B}}\left|{\mathit{J}}_0\right|\frac{{\varsigma }_{m\text{+}1}-{\varsigma }_m}{2}\text{d}\xi \text{d}\eta \text{d}{\varsigma }^{\text{v}}}}}}$$

(13b)

$${\mathit{K}}_{22}^{\hspace{0.17em}}={\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\overline{\mathit{B}}}_{\hspace{0.17em}}^{\text{T}}{\mathit{D}}_m^{\hspace{0.17em}}\overline{\mathit{B}}\left|\mathit{J}\right|\frac{{\varsigma }_{m\text{+}1}-{\varsigma }_m}{2}\text{d}\xi \text{d}\eta \text{d}{\varsigma }^{\text{v}}}}}}$$

(13c)

$${\mathit{M}}_{\text{e}}={\displaystyle \sum _{m\text{=1}}^{\text{n}}{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_{-1}^1{\rho }_{\text{m}}{\mathit{N}}^{\text{T}}\mathit{N}\left|\mathit{J}\right|\frac{{\varsigma }_{m\text{+}1}-{\varsigma }_m}{2}\text{d}\xi \text{d}\eta \text{d}{\varsigma }^{\text{v}}}}}}$$

(13d)

where $\mathit{J}$ is the Jacobian matrix and is given by

$$\mathit{J}=\left[\begin{array}{ccc}\hfill \frac{\partial x}{\partial \xi }\hfill & \hfill \frac{\partial y}{\partial \xi }\hfill & \hfill \frac{\partial z}{\partial \xi }\hfill \\ \hfill \frac{\partial x}{\partial \eta }\hfill & \hfill \frac{\partial y}{\partial \eta }\hfill & \hfill \frac{\partial z}{\partial \eta }\hfill \\ \hfill \frac{\partial x}{\partial \varsigma }\hfill & \hfill \frac{\partial y}{\partial \varsigma }\hfill & \hfill \frac{\partial z}{\partial \varsigma }\hfill \end{array}\right]$$

(14)

Note that the Jacobian matrix should use ${\mathit{J}}_0$ while calculating ${\mathit{K}}_{12}$ , i.e. ${\mathit{J}}_0\text{\hspace{0.17em}=\hspace{0.17em}}\mathit{J}|{\hspace{0.17em}}_{\left(\xi \text{=}0,\eta \text{=}0,\varsigma \text{=}0\right)}$ , just as existence of the one-layer incompatible element. In that way, the above-mentioned element mass and stiffness matrices can be efficiently calculated by using the two-points Gauss rule.

Application method of the LI-SOLID

The application method of the LI-SOLID element would be introduced in detail for coated blades and the uncoated disk of blisks. Note that the developed modelling method can be used to any actual coated blisk where the blades are uncoupled. It needs only the element node coordinates corresponding to the maximum coating thickness, and even the thickness of the blades is smaller, it only needs to increase the number of elements.

In this work, Figure 4(a) shows the coated element with the thickest coating. With the thickness of the coating being changed, it only needs to divide the above-mentioned lamination element into five layers as shown in Figure 4(b). Here, the third layer represents the substrate of blade, and the second and fourth layers represent the hard coating, with the fifth and first layers being considered as the virtual layers. By setting the Young’s modulus, the density and other material parameters of the virtual layers to 0, the new stiffness and mass matrices can be recalculated quickly when coating thickness is changed. Even the coating thickness is 0, the LI-SOLID element can also be used here by setting the material parameters of the first, second, fourth and fifth layers to 0. Note that for the variable coating thickness, it only needs to change the thickness values of each layer in equation (12) under local coordinate, which do not need to regenerate the finite element mesh manually. In this way, if the coating thickness is random, the method is also suitable, and the efficiency of the solution can be significantly improved due to avoiding a re-meshing process.

Figure 4.

Sector of LI-SOLID element. (a) With the thickest coating and (b) with the variable coating thickness.

For the uncoated disk, the element without hard coating is just a special case of the three-layered LI-SOLID element, as it only needs to set the Young’s modulus, the density and other material parameters of all layers to the same one.

Therefore, the developed element is comparatively more convenient to simulate the coated blisk with the variable coating thickness and could improve the computational efficiency significantly.

Analysis method of coated blisk

Equations of motion

The forced vibration equation in frequency domain of a blisk with hard coating on blades can be expressed as

$$\left[\mathit{K}\text{+i}\left({\eta }_{\text{b}}{\mathit{K}}_{\text{b}}+{\eta }_{\text{h}}{\mathit{K}}_{\text{h}}\right)-{\omega }^2\mathit{M}\right]\mathit{X}=\mathit{F}$$

(15)

where $\mathit{K}$ and $\mathit{M}$ denote the stiffness and mass matrices of the coated blisk with hard coating on both sides of each blade, respectively. ${\eta }_{\text{h}}$ and ${\eta }_{\text{b}}$ represent the loss factor of the hard coating and the blisk. ${\mathit{K}}_{\text{b}}$ and ${\mathit{K}}_{\text{h}}$ are the contribution of the substrate and the hard coating to the stiffness matrix. $\mathit{X}$ is the displacement amplitude vector, and $\mathit{F}$ is a harmonic external force. Note that the system stiffness and mass matrices can be solved in the following process.

For a blisk, the system stiffness and mass matrices can be represented by those of N sectors, given by

$$\mathit{K}={\mathit{G}}^T\left[\begin{array}{c}\hfill {\mathit{k}}_1\hfill \\ \hspace{0.17em}& \hfill \ddots \hfill \\ \hspace{0.17em}& \hspace{0.17em}& \hfill {\mathit{k}}_N\hfill \end{array}\right]\mathit{G}={\mathit{G}}^T\overline{\mathit{K}}\mathit{G}$$

(16a)

$$\mathit{M}={\mathit{G}}^T\left[\begin{array}{c}\hfill {\mathit{m}}_1\hfill \\ \hspace{0.17em}& \hfill \ddots \hfill \\ \hspace{0.17em}& \hspace{0.17em}& \hfill {\mathit{m}}_N\hfill \end{array}\right]\mathit{G}={\mathit{G}}^T\overline{\mathit{M}}\mathit{G}$$

(16b)

where $\mathit{G}$ is the structural system matrix and $\overline{\mathit{K}}$ and $\overline{\mathit{M}}$ are the stiffness and mass matrices of the uncoupled system. Symbols ${\mathit{k}}_{\text{1}}\dots {\mathit{k}}_N$ and ${\mathit{m}}_{\text{1}}\dots {\mathit{m}}_N$ mean the stiffness matrices and mass matrices of N sectors, where for each sector it is sorted according to the DOFs of the interface t, the internal nodes g and the interface p as shown in Figure 5, specifically expressed as

$${\mathit{k}}_j=\left[\begin{array}{ccc}\hfill {\mathit{k}}_j^{\text{tt}}\hfill & \hfill {\mathit{k}}_j^{\text{tg}}\hfill & \hfill {\mathit{k}}_j^{\text{tp}}\hfill \\ \hfill {\mathit{k}}_j^{\text{gt}}\hfill & \hfill {\mathit{k}}_j^{\text{gg}}\hfill & \hfill {\mathit{k}}_j^{\text{gp}}\hfill \\ \hfill {\mathit{k}}_j^{\text{pt}}\hfill & \hfill {\mathit{k}}_j^{\text{pg}}\hfill & \hfill {\mathit{k}}_j^{\text{pp}}\hfill \end{array}\right],{\mathit{m}}_j=\left[\begin{array}{ccc}\hfill {\mathit{m}}_j^{\text{tt}}\hfill & \hfill {\mathit{m}}_j^{\text{tg}}\hfill & \hfill {\mathit{m}}_j^{\text{tp}}\hfill \\ \hfill {\mathit{m}}_j^{\text{gt}}\hfill & \hfill {\mathit{m}}_j^{\text{gg}}\hfill & \hfill {\mathit{m}}_j^{\text{gp}}\hfill \\ \hfill {\mathit{m}}_j^{\text{pt}}\hfill & \hfill {\mathit{m}}_j^{\text{pg}}\hfill & \hfill {\mathit{m}}_j^{\text{pp}}\hfill \end{array}\right],j=1,\dots ,N$$

(17)

Figure 5.

Local Cartesian coordinate systems of the blisk.

To obtain the structural system matrix $\mathit{G}$ , the derived process is given as follows. The displacement vector corresponding to $\overline{\mathit{K}}$ and $\overline{\mathit{M}}$ in equation (16), named as the expanded displacement vector, can be written as

$${\mathit{x}}_{\mathrm{expand}}=\left[{\mathit{x}}_{1t}^T\hspace{1em}{\mathit{x}}_{1g}^T\hspace{1em}{\mathit{x}}_{1p}^T\hspace{1em}{\mathit{x}}_{2g}^T\hspace{1em}{\mathit{x}}_{2p}^T\dots \hspace{1em}{\mathit{x}}_{\mathit{Nt}}^T\hspace{1em}{\mathit{x}}_{\mathit{Ng}}^T\hspace{1em}{\mathit{x}}_{\mathit{Np}}^T\right]$$

(18)

where ${\mathit{x}}_{j\text{t}}^{\text{T}}$ , ${\mathit{x}}_{j\text{g}}^{\text{T}}$ and ${\mathit{x}}_{j\text{p}}^{\text{T}}$ denote the displacement vectors of the interface t, internal nodes g and interface p for the jth sector. In fact, the interface t of the previous sector is consistent with the interface p of the next sector, i.e. the displacement vector in equation (18) is repeated. So it is necessary to eliminate the repeated displacement vector of the interfaces. Then, the real displacement vector for a blisk becomes

$${\mathit{x}}_{\mathrm{basic}}={\left[{\mathit{x}}_{1t}^T{\mathit{x}}_{\lg }^T{\mathit{x}}_{2t}^T{\mathit{x}}_{2g}^T\hspace{1em}\cdots \hspace{1em}{\mathit{x}}_{\mathit{Nt}}^T{\mathit{x}}_{\mathit{Ng}}^T\right]}^{\text{T}}$$

(19)

The displacement coordination relationship between the expanded displacement vector and the real displacement vector can be expressed as

$${\mathit{x}}_{\mathrm{expand}}=\left[\begin{array}{ccc}\hfill \mathit{H}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \mathit{H}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \ddots \hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \mathit{H}\hfill \end{array}\right]\left[\begin{array}{ccccc}\hfill {\mathit{E}}_1\hfill & \hfill {\mathit{E}}_2\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathit{E}}_1\hfill & \hfill {\mathit{E}}_2\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill {\mathit{E}}_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill {\mathit{E}}_1\hfill \end{array}\right]={\mathit{x}}_{\mathrm{basic}}=\mathit{G}{x}_{\mathrm{basic}}$$

(20)

where,

$$\mathit{H}=\left[\begin{array}{ccc}{\mathit{I}}_{\text{t}}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& {\mathit{I}}_{\text{g}}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \mathit{R}\end{array}\right],\hspace{1em}{\mathit{E}}_1=\left[\begin{array}{cc}{\mathit{I}}_{\text{t}}& 0\\ 0& {\mathit{I}}_{\text{g}}\\ 0& 0\end{array}\right],\hspace{1em}{\mathit{E}}_2=\left[\begin{array}{cc}0& 0\\ 0& 0\\ {\mathit{I}}_{\text{t}}& 0\end{array}\right]$$

(21)

where ${\mathit{I}}_{\text{t}}$ and ${\mathit{I}}_{\text{g}}$ are identity matrices with the dimension equal to the number of DOFs for interface t and interface g, respectively. $\mathit{R}$ is the linear transformation matrix, which represents the displacement vector transformation of the interface t for the (j + 1)th sector to the coordinate system of the jth sector, as the FEM is created in Cartesian coordinate system. The linear transformation matrix $\mathit{R}$ can be expressed as

$$\mathit{R}={\mathit{I}}_{\text{t}}\otimes \left[\begin{array}{ccc}\cos \theta & -\sin \theta & \hspace{0.17em}\\ \sin \theta & \cos \theta & \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& 1\end{array}\right],\theta =\frac{2\pi }{N}$$

(22)

where $\theta$ indicates the sector angle and $\otimes$ represents Kronecker product.

Then, the specific expression for the assembly matrices $\mathit{G}$ can be given by

$$G=\left[\begin{array}{ccc}\hfill \mathit{H}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \mathit{H}\hfill & \hfill \hspace{0.17em}\hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \ddots \hfill \\ \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \hspace{0.17em}\hfill & \hfill \mathit{H}\hfill \end{array}\right]\left[\begin{array}{ccccc}\hfill {\mathit{E}}_1\hfill & \hfill {\mathit{E}}_2\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathit{E}}_1\hfill & \hfill {\mathit{E}}_2\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill {\mathit{E}}_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill {\mathit{E}}_1\hfill \end{array}\right]=\left[\begin{array}{ccccc}\hfill \mathit{H}{\mathit{E}}_1\hfill & \hfill \mathit{H}{\mathit{E}}_2\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathit{H}{\mathit{E}}_1\hfill & \hfill \mathit{H}{\mathit{E}}_2\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill \mathit{H}{\mathit{E}}_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill \mathit{H}{\mathit{E}}_1\hfill \end{array}\right]$$

(23)

Forced response analysis

For a single sector with variable coating thickness, the stiffness matrix ${\mathit{k}}_j$ and the mass matrix ${\mathit{m}}_j$ can be solved by the LI-SOLID element. Then, equation (15) can be obtained by the method described in Equations of motion section. Solving the characteristic equation, the eigenvalue ${\lambda }_r^{*}$ and eigenvector ${\mathit{\phi }}_r^{*}$ can be computed by

$$\left[-{\lambda }_r^{*}\mathit{M}+\mathit{K}\right]{\mathit{\phi }}_r^{*}=0r=1,2,\dots n$$

(24)

Furthermore, the truncated mode matrix ${\mathit{\phi }}^{*}$ can be obtained by forming the first rth order eigenvectors.

By multiplying ${\mathit{\phi }}^{*}$ , equation (15) can be rewritten as

$${\mathit{\Lambda }}^{\text{o}}\mathit{p}-{\omega }^2\mathit{p}+\text{i}\left({\eta }_{\text{b}}{\widehat{\mathit{K}}}_{\text{b}}+{\eta }_{\text{h}}{\widehat{\mathit{K}}}_{\text{h}}\right)\mathit{p}=\widehat{\mathit{f}}$$

(25)

where ${\mathit{\Lambda }}^{\text{o}}$ is the truncated eigenvalue matrix corresponding to the truncated mode matrix ${\mathit{\phi }}^{*}$ and others are

$$\{\begin{array}{l}{\widehat{\mathit{K}}}_{\text{h}}={\mathit{\phi }}^{*\text{H}}{\mathit{K}}_{\text{h}}{\mathit{\phi }}^{*}\\ {\widehat{\mathit{K}}}_{\text{b}}={\mathit{\phi }}^{*\text{H}}{\mathit{K}}_b{\mathit{\phi }}^{*}\\ \widehat{\mathit{f}}={\mathit{\phi }}^{*\text{H}}\mathit{F}\end{array}$$

(26)

Furthermore, the forced response in frequency domain of the blisk subjected to engine order (EO) excitation can be obtained by the classic mode superposition method. Note that the above-mentioned matrices are all sparse matrices.

Model verification and application

In this work, an actual blisk attached with hard coating on both sides of each blade is used to verify the developed modelling method. The geometric parameters of the vibration model used in this section are listed in Table 1, and the boundary constrains are shown in Figure 6. In order to avoid excessive additional mass and ensure better damping performance, a range of 0–0.3 mm for coating thickness is used, i.e. the maximum thickness of coating is 0.3 mm (10% of blade thickness).

Table 1.

Geometrical parameters of the blisk and hard coating.

Disk		Hard-coating blades
Inner radius of hub	50 mm	Height of blade	50 mm
Outer radius of hub	80 mm	Width of blade	25 mm
Thickness of hub	19 mm	Thickness of blade	3 mm
Inner radius of rim	182 mm	Number of blade	18
Outer radius of rim	200 mm	Coating thickness	0.3 mm
Thickness of rim	20 mm	Thickness of sternum	7 mm

Figure 6.

Vibration model of the coated blisk.

Description of the experimental test

Here, an experimental test for a steel blisk with 0.3 mm coating on both sides of each blade is used to verify the modelling method based on ANSYS in literature^35,38 (the method was also described in Disadvantages of vibration analysis for coated blisk based on ANSYS section) and the developed method. It should be noted that the steel blisk is usually used in reported researches.

The Young’s modulus, the density, the loss factor and the Poisson’s ratio of the Q235-A(F) steel blisk are 208 GPa, 7860 kg/m3, 0.0004 and 0.28, respectively, and those of the coating material NiCoCrAlY + YSZ are 54.45 GPa, 5600 kg/m3, 0.0212 and 0.3, respectively. Figure 7 shows the hard-coating blisk. Then, the natural frequencies obtained by experimental test, the modelling method based on ANSYS and the developed method are listed in Table 2.

Figure 7.

Hard-coating blisk.

Table 2.

Natural frequencies and relative differences.

Experiment (A)	ANSYS (B)	Developed method (C)	Difference |A−B|/A	Difference |A−C|/A
702.5	697.2	695.98	0.75	0.92
713.42	721.94	722.81	1.19	1.31
739.5	729.76	726.00	1.32	1.82
746	781.11	781.08	4.71	4.70
750.25	792.15	793.51	5.58	5.76
761.5	795.52	797.37	4.47	4.71
767.7	796.92	798.98	3.81	4.07
776.75	797.58	799.75	2.68	2.96
786.55	797.89	800.11	1.44	1.72
801.25	797.98	800.22	0.41	0.13

The results show that the maximum difference between the results obtained by experiment and ANSYS for the coated blisk is 5.58%. And the one between the results obtained by experiment and the developed modelling method is 5.76%. In that way, it can be concluded that the modelling method in literature^35,38 and the developed method could be used to calculate the vibration characteristics of the coated blisk.

Verification of the developed modelling method

In practice, the blisk made of titanium material is widely used in aeroengine and its modelling method needs to be focus on. Since the blisk made of titanium material has not yet been obtained in the study, so only a titanium blisk attached with hard coating on both sides of each blade is used to verify the developed LI-SOLID element.

For titanium material, the Young’s modulus is 103.213 GPa with the density being 4200 kg/m³, the loss factor being 0.0007 and Poisson’s ratio being 0.33. In addition, the coating material is also NiCoCrAlY+YSZ and its thickness is 0.3 mm.

Two verification cases are analysed. It should be noted that the element and node number are 32778 and 46800, respectively, for the FEM by applying LI-SOLID element, and those for the FEM by applying SOLID185 in ANSYS software are 38728 and 53712, respectively.

For case 1, the LI-SOLID element is used to solve the natural frequencies and the forced responses subjected to EO excitation for the tuned coated blisk. Figure 8 shows the natural frequencies sorted with nodal diameters (NDs) obtained by the developed LI-SOLID element model solved in MATLAB and the HFEM using three layers SOLID185 in ANSYS. Here, only the first three orders frequencies of ND 3 and 6 are listed in Table 3. Note that the EO excitation is applied to the z direction of the tip nodes of each blade, and the pick-up points are same to the excitation points. The forced responses subjected to three EO and six EO excitations are also obtained based on the above two elements, and the relevant results are shown in Figure 9(a) and (b), with the resonance responses being listed in Table 4. Note that for the tuned blisk, the forced responses of each pick-up point are exactly the same.

Figure 8.

Natural frequencies for the coated blisk.

Table 3.

First three orders frequencies of ND 3 and ND 6 for the coated blisk (Hz).

Order	1	2	3
ND 3
Three layers SOLID185 (S)	730.14	2276.7	2640.8
LI-SOLID (L)	727.74	2319.54	2621.58
Difference $\left(\left|S-L\right|/S\right)$	0.33%	1.88%	0.73%
ND 6
Three layers SOLID185 (S)	745.5	2595.9	4499.7
LI-SOLID (L)	741.69	2568.89	4498.60
Difference $\left(\left|S-L\right|/S\right)$	0.51%	1.04%	0.025%

ND: nodal diameter.

Figure 9.

Forced responses with coating. (a) Forced responses under three EO excitations and (b) forced responses under six EO excitations.

Table 4.

Resonance response under three EO and six EO excitations for the coated blisk (m).

Order	1	2	3
ND 3
Three layers SOLID185 (S)	9.33 × 10−4	7.23 × 10−5	8.39 × 10−5
LI-SOLID (L)	9.82 × 10−4	7.51 × 10−5	8.61 × 10−4
Difference $\left(\left|S-L\right|/S\right)$	5.25%	3.87%	2.62%
ND 6
Three layers SOLID185 (S)	9.37 × 10−4	8.54 × 10−5
LI-SOLID (L)	9.85 × 10−4	8.91 × 105
Difference $\left(\left|S-L\right|/S\right)$	5.12%	4.02%

ND: nodal diameter.

It can be seen that the maximum difference of natural frequencies is 1.88% for the modelling methods using the LI-SOLID and SOLID185 in ANSYS software. For the resonance responses, the maximum difference is 5.25%. So the presented element can be used to solve the vibration characteristics of coated blisks.

For case 2, the natural frequencies and the forced responses for the uncoated blisk are solved by the developed modelling method and the above-mentioned methods based on ANSYS software. The relevant results are shown in Figures 10 and 11 and Tables 5 and 6. Note that the pick-up points and the excitation points are same to case 1.

Figure 10.

Natural frequencies without coating.

Figure 11.

Forced responses for uncoated blisk. (a) Forced responses under three EO excitations and (b) forced responses under six EO excitations.

Table 5.

First three orders frequencies of ND 3 and ND 6 for the uncoated blisk (Hz).

Order	1	2	3
ND 3
SOLID-185 (S)	707.5	2390.6	2598.9
LI-SOLID (L)	706.72	2389.92	2598.23
Difference $\left(\left|S-L\right|/S\right)$	0.11%	0.026%	0.027%
ND 6
SOLID-185 (S)	717.6	2484.6	4403.5
LI-SOLID (L)	716.71	2483.82	4400.01
Difference $\left(\left|S-L\right|/S\right)$	0.13%	0.031%	0.078%

ND: nodal diameter.

Table 6.

Resonance responses under three EO and six EO excitations for the uncoated blisk (m).

Order	1	2	3
ND 3
SOLID-185 (S)	1.704 × 10−2	1.072 × 10−3	8.601 × 10−4
LI-SOLID (L)	1.706 × 10−2	1.071 × 10−3	8.598 × 10−4
Difference $\left(\left|S-L\right|/S\right)$	0.12%	0.091%	0.035%
ND 6
SOLID-185 (S)	1.612 × 10−2	1.851 × 10−3
LI-SOLID (L)	1.614 × 10−2	1.852 × 10−3
Difference $\left(\left|S-L\right|/S\right)$	0.12%	0.054%

ND: nodal diameter.

The results indicate that both the natural frequencies and the forced responses for the modelling method using the LI-SOLID element are in good consistency with those using SOLID185. Note that using the LI-SOLID element to obtain the vibration characteristics for the uncoated blisk is equivalent to change the maximum coating thickness from 0.3 mm to 0 mm, which can be considered as a special case of changing the coating thickness.

Therefore, it can be concluded that the LI-SOLID element can be used to simulate the variable coating thickness.

Model application

An exhaustive research is conducted to obtain the effect of different thicknesses of coating on the forced response. First, the stiffness and mass matrices for the tuned coated blisk with 0.3 mm coating on blades can be solved by the presented LI-SOLID element. When the coating thickness is changed, the new stiffness and mass matrices of blisk can be resolved by the LI-SOLID element. Note that here by using the LI-SOLID with the virtual layers as shown in Figure 4(b), the finite element mesh does not need to regenerate, and we only need to change the thickness values of each layer in equation (12) under local coordinate. Then, the method described in Analysis method of coated blisk section can be used to obtain the forced response. But with the coating thickness being changed in ANSYS software, it needs to regenerate the finite element mesh, which can be very inconvenient and result in an inefficient computation in in the process of damping design. Herein, only the uniform change of the coating thickness for each blade is analysed.

The variable coating thicknesses are 0 mm, 0.05 mm, 0.1 mm, 0.15 mm, 0.2 mm, 0.25 mm and 0.3 mm. Tables 7 and 8 list the obtained natural frequencies of ND 3 and ND 6. Figures 12 and 13 show the obtained forced responses under three EO and six EO excitations, with the resonance responses being listed in Tables 9 and 10.

Table 7.

First three orders frequencies of ND 3 with the variable coating thicknesses (Hz).

Order	1	2	3
0 mm	706.72	2389.92	2598.23
0.05 mm	708.73	2380.10	2590.96
0.1 mm	710.64	2368.94	2587.54
0.15 mm	713.21	2356.83	2588.51
0.2 mm	716.37	2343.79	2593.53
0.25 mm	720.026	2329.99	2602.17
0.3 mm	727.74	2319.54	2621.58

Table 8.

First three orders frequencies of ND 6 with the variable coating thicknesses (Hz).

Order	1	2	3
0 mm	716.71	2483.82	4400.01
0.05 mm	719.25	2491.29	4408.57
0.1 mm	721.72	2500.67	4417.76
0.15 mm	724.90	2512.48	4430.32
0.2 mm	728.70	2526.47	4445.55
0.25 mm	733.07	2542.42	4462.72
0.3 mm	741.69	2568.89	4498.60

Figure 12.

Forced responses for coated blisk under three EO excitations (m). (a) Forced responses containing the first order and (b) forced responses containing the second and third order.

Figure 13.

Forced responses for coated blisk under six EO excitations. (a) Forced responses containing the first order and (b) forced responses containing the second order.

Table 9.

Resonance response of three EO excitations with the variable coating thicknesses (m).

Order	1	2	3
0 mm	1.706 × 10−2	1.071 × 10−3	8.598 × 10−4
0.05 mm	4.998 × 10−3	4.134 × 10−4	3.861 × 10−4
0.1 mm	2.876 × 10−3	2.474 × 10−4	2.456 × 10−4
0.15 mm	1.995 × 10−3	1.716 × 10−4	1.776 × 10−4
0.2 mm	1.513 × 10−3	1.284 × 10−4	1.376 × 10−4
0.25 mm	1.211 × 10−3	1.006 × 10−4	1.113 × 10−4
0.3 mm	9.82 × 10−4	7.51 × 10−5	8.61 × 10−5

Table 10.

Resonance responses of six EO excitations with the variable coating thicknesses (m).

Order	1	2
0 mm	1.614 × 10−2	1.852 × 10−3
0.05 mm	4.613 × 10−3	5.175 × 10−4
0.1 mm	2.633 × 10−3	2.934 × 10−4
0.15 mm	1.815 × 10−3	2.013 × 10−4
0.2 mm	1.368 × 10−3	1.513 × 10−4
0.25 mm	1.090 × 10−3	1.199 × 10−4
0.3 mm	9.85 × 10−4	8.91 × 10−5

The results listed in Tables 7 and 8 indicate that the natural frequencies for one certain order changed monotonously with the increased coating thickness. For the forced response, with the increase of coasting thickness, the resonance response is reduced. Especially for coated 0.05 mm, the resonance response for the uncoated blisk is reduced sharply. But when the coating thickness exceeds 0.1 mm, the resonance response is reduced a little. In short, the damping properties and the additional mass of the coating should be considered together in the process of damping design.

Conclusion

This paper highlights a new FEM, which can be used to incorporate variable coating thickness in the damping design of coated blisk by introducing the additional virtual layers of the LI-SOLID element. Then, based on the developed element, the effect of different coating thickness on forced response was analysed. Some important conclusions are listed as follows.

To simulate the blisk with the variable coating thickness in the damping design, the virtual layers are introduced to the LI-SOLID element and the thicknesses of each layer can be changed by the coordinate transformation under local coordinate. With the material parameters of the virtual layers being setting to 0, the new stiffness and mass matrices can be solved without regenerating the finite element mesh.

Compared with the natural frequencies and forced responses obtained by the HFEM in ANSYS, the presented LI-SOLID element is verified that it can be used to improve the computational efficiency with variable coating thickness. For the blisk with the maximum coating thickness, the results between the developed method and the one in ANSYS software show that the maximum natural frequencies difference is 1.88% and the maximum difference of the resonance responses 5.25%.

With the increase of coating thickness, the natural frequencies for one certain order changed monotonously and the forced response is reduced significantly. Even the thinner coating could also achieve better damping performance. The damping properties and the additional mass introduced by the coating should be considered together in the damping design process with coating.