1 引言

在航空发动机风扇、压气机、涡轮上采用整体叶盘结构能够实现结构简化、降低重量的目标，有效提高发动机推重比。为此，众多航空发动机企业竞相研发整体叶盘(Integrally bladed disks，Blisk)技术^{[1, 2]}。在理想情况下，整体叶盘各个扇区的物理性质与几何参数是完全相同的，称为谐调性质；在实际情况下，制造误差，材料偏差和磨损等多种因素会导致扇区参数的改变，谐调性质被破坏，称为失谐。整体叶盘的盘体更薄，盘-叶刚度相差不多，因而振动模式多为盘-片耦合振动，而失谐会引起模态局部化现象(振动能量集中到少数叶片上)，从而加剧叶盘结构的局部化振动^[3]。过大的共振应力会使整体叶盘发生振动疲劳失效，严重威胁飞行安全。因此，能够利用随机分析(Monte Carlo仿真，摄动法和谱方法等^[4~6]准确预测失谐整体叶盘(失谐的形式、位置不确定)的最大动力学响应是非常重要的。

为了获得整体叶盘较准确的振动特性，一般采用高保真有限元模型进行动力学分析。但是，全尺寸有限元模型一般具有上百万甚至上千万的自由度，更重要的是计算过程中需要重复求解刚度矩阵的逆矩阵，这会使得结构分析变得非常困难。因此，学者们相继提出了许多模型降阶技术，主要包括CMS法，SNM法，CMM法和IMM法等^[7]；此外，许多改进的随机仿真技术也逐渐被发展出来，主要包括重要采样法，加速Monte Carlo仿真法和子集仿真法等^[8]。然而，降阶的有限元模型依然具有较大数量的模态和界面自由度，仍然需要通过重复求解刚度矩阵的逆矩阵来得到失谐整体叶盘的频域响应，并没有彻底解决有限元模型计算效率低这一关键难题。而在Monte Carlo仿真分析中，通常需要上千次的计算才能得到比较可信的统计规律，计算量仍然很大。

对于LU分解法，一般具有比较高的计算精度，而且其刚度矩阵表示为两个带状三角矩阵之积，在求解较小规模的分析模型时，可以通过两个子矩阵的逆变来有效提高求解动力学响应的效率，但是，在分析规模较大的失谐整体叶盘模型时效率非常低下。Neumann展开式是一种包含代数方程分级线性解的迭代过程，主要由确定矩阵和随机矩阵共同组成，是一种间接求解不确定动刚度矩阵的方法。因为Neumann级数可以用来对摄动矩阵进行系数拓展，所以刚度矩阵的逆矩阵也可以用一系列级数来表示。此外，由于Neumann展开式法可以较容易地拓展到高阶次项，所以这种方法所容许的随机波动范围比其他摄动法的波动范围更大^[9]。但是，当迭代矩阵的谱半径过大时，Neumann展开式法的计算性能将会失效^{[10, 11]}。Yamazaki等^[12]提出了一种可以在较大摄动情况下保证收敛的修正方法；Yuan等^[13]研究了该方法在失谐叶盘随机分析中的可行性，结果表明，该方法的计算效率是矩阵直接分解法的两倍，而最大误差却小于2%。但是，当激励频率接近共振频率时(尤其是系统具有较小阻尼或较大摄动时)，该方法并不能够保证收敛。此外，在航空发动机叶片的动力学特性分析中，Yuan等^[14]提出了一种参数化减缩模型，该模型具有较高的计算精度和灵敏度，但这是以牺牲计算效率为代价的。

基于以上调研，本文提出一种改进的Neumann-LU混合法来求解失谐整体叶盘的受迫响应。该方法可以实现Neumann级数法效率高和LU分解法精度高的优势互补：一方面，在非共振频率附近的大区域采用Neumann级数频域展开式，克服了LU分解法在大区域计算效率低下的缺陷，极大地提高了计算效率；另一方面，在共振频率附近小区域采用LU分解法以提高计算精度，克服了Neumann级数在共振区域失效的缺陷，有效提高了计算精度；同时，通过Neumann级数展开项数定义了一种自适应阀值，得到两种矩阵求逆方法在不同求解域的自动转换。首先，根据整体叶盘的动力学方程推导了失谐整体叶盘的频域响应公式；接着，分别详细阐述了Neumann展开式法，LU分解法的原理和Neumann-LU混合法求解过程；最后，以失谐整体叶盘作为研究实例，分别运用Neumann-LU混合法和单一LU分解法求解叶片在阶次激励下的频域响应，并从计算效率和计算精度方面进行比较。结果表明，在保证计算精度的基础上，改进的Neumann-LU混合法可以明显提高失谐整体叶盘随机分析的计算效率。

2 失谐整体叶盘受迫响应

在理想状态下，谐调整体叶盘在受外力作用时的动力学响应频域方程可以表示为

$\left\{ \begin{array}{l} \left[{{\mathit{\boldsymbol{K}}_{\rm{0}}}-{\mathit{\omega }^{\rm{2}}}{\mathit{\boldsymbol{M}}_{\rm{0}}} + i\mathit{\boldsymbol{D}}} \right]\mathit{\boldsymbol{q = F}}\\ \mathit{\boldsymbol{F}} = {F_0}{\left[{{{\rm{e}}^{{\rm{j}}{\phi _1}}}{\rm{\;\;\;\;}}{{\rm{e}}^{\mathit{j}{\phi _2}}}{\rm{\;\;\;\;}}{{\rm{e}}^{\mathit{j}{\phi _3}}} \cdots } \right]^{\rm{T}}}\\ {\phi _i} = \frac{{2\pi C\left( {i - 1} \right)}}{{{N_{\rm{b}}}}} \end{array} \right. $

(1)

式中K₀, M₀和D分别表示谐调的刚度矩阵、质量矩阵和阻尼矩阵, q和F分别表示位移向量(受迫响应)和外部激励力向量(本文采用谐波激励), ϕ_i表示第i个叶片的相位角, C表示激励阶次, N_b表示叶片的数量。

其中，当系统考虑结构阻尼(如滞后阻尼)时，阻尼矩阵D (ηK)与刚度矩阵成比例关系，η表示阻尼系数。

假设ΔK和ΔM分别表示由失谐造成的摄动刚度矩阵和摄动质量矩阵, K_m和M_m分别表示失谐刚度矩阵和失谐质量矩阵, 那么, 失谐整体叶盘在受外力作用时的动力学响应频域方程可以表示为

$\left\{ \begin{array}{l} \left[{{\mathit{\boldsymbol{K}}_{\rm{m}}}-{\mathit{\omega }_{\rm{2}}}{\mathit{\boldsymbol{M}}_{\rm{m}}} + i\mathit{\boldsymbol{D}}} \right]\mathit{\boldsymbol{q = F}}\\ {\mathit{\boldsymbol{K}}_{\rm{m}}} = {\mathit{\boldsymbol{K}}_{\rm{0}}} + \Delta \mathit{\boldsymbol{K}}\\ {\mathit{\boldsymbol{M}}_{\rm{m}}} = {\mathit{\boldsymbol{M}}_{\rm{0}}} + \Delta \mathit{\boldsymbol{M}} \end{array} \right. $

(2)

式(2)可以表示为谐调系统Z₀和摄动系统ΔZ两部分之和，经整理可得到

$\left\{ \begin{array}{l} \left[{{\mathit{\boldsymbol{Z}}_{\rm{0}}} + \Delta \mathit{\boldsymbol{Z}}} \right]\mathit{\boldsymbol{q = F}}\\ {\mathit{\boldsymbol{Z}}_{\rm{0}}} = {\mathit{\boldsymbol{K}}_{\rm{0}}} - {\mathit{\omega }^{\rm{2}}}{\mathit{\boldsymbol{M}}_{\rm{0}}} + i\mathit{\boldsymbol{D}}\\ \Delta \mathit{\boldsymbol{Z = }}\Delta \mathit{\boldsymbol{K}} - {\mathit{\omega }^{\rm{2}}}\Delta \mathit{\boldsymbol{M}} \end{array} \right. $

(3)

则失谐整体叶盘的受迫响应可以表示为

$\mathit{\boldsymbol{q = }}{\left[{{\mathit{\boldsymbol{Z}}_{\rm{0}}} + \Delta {\rm{Z}}} \right]^{ - 1}}\mathit{\boldsymbol{F = }}{\mathit{\boldsymbol{Z}}^{\mathit{\boldsymbol{ - }}{\rm{1}}}}\mathit{\boldsymbol{F}} $

(4)

3 改进的Neumann-LU混合法

由于矩阵Z^-1不是带状矩阵，所以式(4)右侧的乘积运算需要耗费大量的计算资源。本文提出一种基于LU分解法和Neumann展开式法优势互补的Neumann-LU混合法，其基本原理是在非共振区域的较大区域采用Neumann展开式法，克服了矩阵求逆计算效率低下的缺点，在共振区域附近的较小区域采用LU分解法求解响应，克服Neumann展开式法在该区域失效的缺点，从而可以快速得到全频域的精确响应。此外，本文还特别定义了一个自适应阀值，从而可以根据外部力激励频率来自动选择合适的矩阵求逆方法。

3.1 Neumann展开式法

当外部力激励频率ω远离共振频率时(较大区域)，利用Neumann展开式法就不需要对动刚度矩阵Z求逆。将摄动矩阵ΔZ展开，则逆矩阵Z ^-1可以表示为

$\left\{ \begin{array}{l} {\mathit{\boldsymbol{Z}}^{{\rm{ - 1}}}} = {\left[{{\mathit{\boldsymbol{Z}}_{\rm{0}}} + \Delta \mathit{\boldsymbol{Z}}} \right]^{ - 1}} = \left[{I-\mathit{\boldsymbol{E + }}{\mathit{\boldsymbol{E}}^{\rm{2}}}-{\mathit{\boldsymbol{E}}^{\rm{3}}} \cdots } \right]\mathit{\boldsymbol{Z}}_0^{ - 1}\\ \mathit{\boldsymbol{E = Z}}_0^{ - 1}\Delta \mathit{\boldsymbol{Z}} \end{array} \right. $

(5)

将式(5)带入式(4)，则失谐整体叶盘的受迫响应可以用下列级数表示为

$\mathit{\boldsymbol{q = }}{\mathit{\boldsymbol{q}}_{\rm{0}}} - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{q}}_{\rm{0}}} + {\mathit{\boldsymbol{E}}^{\rm{2}}}{\mathit{\boldsymbol{q}}_{\rm{0}}} - {\mathit{\boldsymbol{E}}^{\rm{3}}}{\mathit{\boldsymbol{q}}_{\rm{0}}} + {\mathit{\boldsymbol{E}}^{\rm{4}}}{\mathit{\boldsymbol{q}}_{\rm{0}}} - {\mathit{\boldsymbol{E}}^{\rm{5}}}{\mathit{\boldsymbol{q}}_{\rm{0}}} + \cdots $

(6)

式中q₀表示谐调整体叶盘的响应。假设q_j=E^jq₀, 则由式(6)可以得到

${\mathit{\boldsymbol{q}}_\mathit{j}} = \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{q}}_{\mathit{j - }{\rm{1}}}} $

(7)

假设r_err是预定义的误差值，一般取0.01^[12]，那么当满足式(8)的条件时，即可判定Neumann级数在第j阶收敛，式(5)展开式可在第j阶截断。

$\frac{{\left\| {\sqrt {\mathit{\boldsymbol{q}}_j^T{\mathit{\boldsymbol{q}}_j}} } \right\|}}{{\left\| {\sqrt {{{\left( {\sum\limits_{n = 0}^j {{{\left( { - 1} \right)}^n}{\mathit{\boldsymbol{q}}_n}} } \right)}^{\rm{T}}}\left( {\sum\limits_{n = 0}^j {{{\left( { - 1} \right)}^n}{\mathit{\boldsymbol{q}}_n}} } \right)} } \right\|}} \le {r_{{\rm{err}}}} $

(8)

3.2 LU分解法

在求解方面，LU分解法比Cholesk分解法具有两方面的优势，其一是当外部力激励频率ω位于共振频率附近时(小区域)，LU分解法具有较高的计算精度和计算效率；其二是由于复阻尼矩阵的存在，不能保证矩阵中的对角元素都是实数。因此，Cholesky分解法不适用于刚度矩阵求逆运算^[15]，本文选择LU分解法来求解刚度逆矩阵。动刚度矩阵Z可表示为两个三角矩阵之积，即

$\mathit{\boldsymbol{LU = Z}} $

(9)

此外，式(9)可以进一步分解成有关未知向量X和q的方程组，如下

$\mathit{\boldsymbol{LX = F, Uq = X}} $

(10)

由于分解后的矩阵L和U拥有和矩阵Z相同的带宽，所以未知向量X和q可通过式(10)快速求得。

3.3 自适应阀值判定

在Neumann展开式法和LU分解法转换时，需要一个合适的频率作为自适应阀值。若选取E的最大特征值λ _max作为自适应阀值，由于需要求解E的特征方程，所以计算效率并不高。而且，当λ _max > 1时，Neumann展开式法不能收敛；当λ _max < 1时，由于Neu- mann展开式法需要计算至少50项才会收敛，所以不能保证Neumann展开式法比LU分解法更有效。有鉴于此，本文采用Neumann展开式法的项数NET (Number of the expansion terms)作为Neumann-LU混合法的自适应阀值。随着项数j不断增大，当式(8)成立时，判定

改进的Neumann-LU混合法的计算过程如下：

$NET = j $

(11)

(1) 在Neumann展开式法中，分别以初始频率f₁和步长Δf进行分析，并采用NET作为转换阀值。

(2) 当到达低阶频率阀值f_b1(f_b1小于第1阶共振频率)时，分别以转换频率f₂(等于第1、2阶共振频率平均值)和-Δf为初始频率和步长，仍然使用Neumann展开式法进行分析，并采用NET作为转换阀值。

(3) 当到达高阶频率阀值f_b2时，Neumann展开式法转换到LU分解法，可求解失谐整体叶盘在第1阶共振区域附近(频率范围为(f_b1, f_b2))的受迫频域响应。

(4) 再次分别以f₂和Δf为初始频率和步长，依据步骤(1)~(3)的过程重复计算，可得到失谐整体叶盘在第2阶共振区域附近。以后的计算过程依次类推。

假设失谐整体叶盘的共振频率表示为[f₀¹ … f₀^j … f₀^J …]，则在上述步骤(2)中提到的转换频率可以通过图 1流程求解得到。

4 研究实例

图 2所示为在商业软件ANSYS中建立的整体叶盘有限元全尺寸模型。该模型大约包含4224个SOLID95 8节点单元，53406个节点，160218个自由度。整体叶盘由钛金属制成(杨氏模量为115GPa，密度为4800kg/m³，泊松比为0.33，阻尼比为0.03%)，叶片约长120mm，宽12mm，厚1.5mm轮盘内径2.5mm，外径80mm，厚2.5mm。以高性能计算机为平台(酷睿四核CPU i7 - 6700 3.4G，RAM 16GB)，使用Block Lanczos求解器对整体叶盘的有限元模型进行模态分析。

图 3所示为谐调整体叶盘结构的固有频率与节径数的关系。随着节径数的增加，轮盘的固有频率变大，并逐渐接近独立叶片的固有频率。节径数较大的水平线表示整体叶盘的振动由叶片模态主导，而节径数较小的斜线表示整体叶盘的振动由轮盘模态主导。在频率转向区域，叶片和轮盘发生强烈的耦合振动，而这对失谐整体叶盘结构的受迫频响产生重要影响^[16]。可以发现，第2族和第3族叶片弯曲模态的固有频率在特定的节径间隔内发生很大的变化，这就意味着叶片间的耦合作用会促使振动能量传输到轮盘上。当发生失谐时，上述机理会造成整体叶盘的模态局部化和过大的受迫频响^[17]。

4.1 失谐整体叶盘

整体叶盘的失谐是重要和不可避免的，在理论分析中必须考虑失谐的重要影响。其中，集中质量法^[18]是一种用来模拟结构失谐的有效方法，而且集中质量多放置在模态位移较大的位置，以确保固有频率能够产生特定的偏离量，大约(-4%，4%)。在ANSYS软件中，可用MASS21单元模拟放置在叶尖的集中质量(模拟整体叶盘的质量变化量^[19])，从而可以获得第1族弯曲模态所需要的摄动量。图 4所示为1~5族模态的叶盘扇区质量摄动量与固有频率偏移量的关系。可以看到，由集中质量引起的叶盘扇区质量摄动量范围约为(-1%，1%)，由质量失谐引起的固有频率偏移量范围约为(-4%，4%)。

图 5(a)和图 5(b)分别为谐调和失谐整体叶盘的六个叶片在C = 0时(5g加速度激励)下的频域响应(叶片编号如图 1所示)。可以清晰地看出，谐调整体叶盘的全部叶片都都有完全相同的频域响应，而且外部激励力仅能在特定节径数激发出模态。相比之下，失谐整体叶盘六个叶片的频域响应均不相同：有些叶片的最大频域响应比谐调状态的最大频域响应要大，有些叶片的最大频域响应却比谐调状态的最大频域响应还要小，而且六个失谐状态下的叶片的响应峰值与谐调状态下的叶片的响应峰值之间的差值也是不尽相同的。此外，失谐整体叶盘的频域响应还包含依赖于节径的谐波，文献[20]对这一现象进行了详细的研究和讨论。

图 6所示为传统Neumann展开式法和单一LU分解法在5，10，20，30，40和50的展开项数下的计算时间比率。可以清晰地发现，为了满足式(8)中的展开式截断误差要求(即r_err ≤ 0.01)，传统Neumann展开式法在分析计算时大约需要15项展开式，此时它与单一LU分解法的计算效率比较相近。此后，随着展开项数的逐渐增加(这意味着谱半径变大)，Neumann展开式法的计算效率要逐渐低于LU分解法的计算效率。因此，在Neumann展开式法和单一LU分解法进行自动转换的过程中，非常有必要选定一个合适的Neumann展开项数作为改进的Neumann-LU混合法中的自适应阀值，这对改进的Neumann-LU混合法的计算效率是至关重要的。

4.2 改进的Neumann-LU混合法的性能

图 7所示为通过Neumann-LU混合法与单一LU分解法分别获得的失谐整体叶盘的三个叶片的频响(叶片编号如图 1所示，固有频率摄动值为2%，阻尼比为0.03%)。从整体趋势上可以看到，Neumann-LU混合法的频响曲线与单一LU分解法的频响曲线具有非常好的一致性。此外，局部放大图(a)清晰显示了两种方法在主共振频率(此处频域响应最大)附近的频域响应的一致性，可以发现，两种方法的最大误差出现在转换频率处(即自适应阀值处)。在自适应阀值处进行矩阵求逆方法的转换是非常必要的，否则计算误差会继续增大，精度无法得到保证。局部放大图(b)和(c)分别清晰显示了失谐整体叶盘在42Hz和85Hz(偏离主频率)的频域响应，可以发现，在高阶模态的自适应阀值处的误差虽然存在，但是误差相对值却小于0.1%。总的来说，不论是在整体趋势上还是在局部相对误差上，改进的Neumann-LU混合法都具有比较高的计算精度。

图 8所示为Neumann-LU混合法与单一LU分解法在不同阻尼系数(0.1%:0.1%:1%)下分别求解50次Monte Carlo仿真分析的计算时间比率(其值越小，则蓝色越深，表示计算效率越高)。可以看到，改进的Neumann-LU混合法的计算时间约占单一的LU分解法计算时间的45%左右，相对于单一的LU分解法，本文提出的改进Neumann-LU混合法的计算效率平均提高了55%(45%~100%)左右。而且可以发现，当整体叶盘的结构阻尼增大时(如纵向区域(a)所示)，公称动刚度Z会呈量级增大，计算规模也会呈量级增大；而当整体叶盘的失谐量减小时(如横向区域(b)所示)，迭代矩阵的谱半径会随之减小，这就使得迭代次数(计算规模)随之增加；在上述两种情况下，改进的Neumann-LU混合法的高效优势将得到进一步突出。

综合以上计算数据和详细分析可知，在保证较高计算精度的重要前提下，本文提出的改进Neu⁃ mann-LU混合法在分析失谐整体叶盘模型时具有非常高的计算效率，而且其高效性能在分析整体叶盘结构阻尼较大或失谐量较小时尤其明显。

5 结论

通过对失谐整体叶盘全尺寸有限元模型的Mon-te Carlo仿真分析发现，本文提出的改进Neumann-LU混合法充分实现了Neumann展开式法和LU分解法的优势互补，其优越的计算性能主要体现在以下两个方面：

(1) 改进的Neumann-LU混合法具有比较高的计算精度，其最大相对误差出现在自适应阀值处，但是不大于0.1%。

(2) 改进的Neumann-LU混合法具有非常高的计算效率，其计算时间仅占单一LU分解法计算时间的45%左右，计算效率大约提高了55%；而且这种高效性能在分析结构阻尼较大或失谐量较小的整体叶盘模型时尤其突出。