引言

航空发动机外部管路是航空发动机的重要组成部分，如同血管一般为航空发动机输送燃料、液压油等工作介质.由于航空发动机外部管路具有结构的复杂性、激励的多源性以及耦合的关联性，在其工作过程中极易发生管路系统共振，从而降低整个系统的安全可靠性^[1，2].目前航空发动机管路的振动特性分析多聚焦于流固耦合方面，缺少对管道模型结构本身的研究.为此，针对航空发动机外部基本管道建立精确的有限元模型是十分必要的.

在工程设计阶段，基于已有参数建立的有限元模型由于受到局部网格划分、材料参数不确定以及模型固有的阶次误差、离散误差等因素的影响^[3]，有限元模型预测结果与实验结果之间会存在一定误差，从而使得初始模型很难反映实际几何结构.为保证模型的状态预测、分析与评估具有实际意义，需要基于模态试验数据对初始有限元模型加以修正.有限元模型修正的本质在于通过修正有限元模型中的相关参数，缩小有限元模拟结果与试验实测结果之间的误差，从而建立出更高精度的有限元模型，以替代实际结构进行分析与监测.

传统的模型修正方法主要分为矩阵法和参数法^[4].参数法通过调整设计参数以实现修正目标，同矩阵法相比更容易在工程中实现应用.代理模型方法是参数法中应用十分广泛的一种.所谓代理模型，通常是指在分析和优化设计过程中可替代那些比较复杂和费时的数值分析的近似数学模型，也称为响应面模型、近似模型或元模型^[5，6].使用代理模型方法可以极大地提高计算效率，还能同其他优化算法相结合.常用的代理模型有多项式响应面模型^[7，8]、径向基函数模型^[9，10]和Kriging模型^[11，12]等.Kriging模型是一种基于随机过程的代理模型，被称为样本点的无偏线性估计函数^[13].Kriging模型在结构优化与可靠性分析领域中应用广泛，其优势在于针对非线性结构进行模型修正具有较强的适用性，在较少样本点的情况下依然能保障较高的修正精度.该模型不仅在地质、水文、气象、环境科学等自然科学领域得到应用，也在航空航天、汽车等工程科学领域得到研究、发展和应用^[14].例如，王文竹等^[15]引入Kriging模型，对鼓式制动器稳定性进行了优化设计，抑制了制动噪声；张文鑫等^[16]使用自适应Kriging结合Monte Carlo模拟以及所提停止准则计算了涡轮盘盘心低周疲劳寿命可靠性，验证了该方法的可靠性；Cui等^[17]将Kriging模型引入空间仪器的热模型修正中，使模型区域的温度值与热平衡实验结果相接近；Yin等^[18]构建了Kriging模型、二次多项式模型和径向基函数模型三种代理模型，针对平面桁架有限元模型进行了有限元模型修正，并对比了三种方法的异同.

本文以航空发动机管道模型为研究对象，通过模态试验和有限元建模得出管道模态试验值和有限元分析值，然后选取合适的待修正参数和试验设计（Design of Experiment，DOE）方法，构建出Kriging代理模型.利用构建好的代理模型代替有限元模型，基于多目标遗传算法（multi objective genetic algorithm，MOGA）实现对管道模型的有限元模型修正.通过将采用Kriging模型的修正结果与灵敏度分析方法的修正结果进行对比，验证了所用模型修正方法的可行性.

1 管道模态试验

1.1 模态试验方案

选取某航空发动机管道模型作为研究对象.管道部件由导管、堵头、锁紧螺母、支座、垫圈、直通管接头、扩口式平管嘴和外套螺母构成，管道模型结构如图1所示.

模型中导管为两端带扩口的直径6mm，壁厚0.6mm的500mm不锈钢直管，密度为7930kg/m³，弹性模量为194GPa，泊松比为0.3.导管与支座之间采用螺栓连接方式，外套螺母和管嘴起到导管与直通管接头之间的密封连接作用.整体模型通过沉头螺栓紧固安装在铸铁试验台上，最大程度上降低了环境的影响.由于本文主要分析管道系统中硬管干模态下的振动特性，故忽略模型两端软管及其连接件，模型各组件材料参数如表1所示.

为分析系统模态特性，选用 LMS 数据采集仪、PCB冲击锤和PCB单轴加速度传感器对管道模型进行模态试验，试验模态测试系统如图2所示.

在试验模态测试过程中，利用冲击力锤施加脉冲激励，采用单点激励多点响应的形式^[19]，将力传感器信号和拾振点处单轴加速度传感器信号反馈给数据采集系统，采集输入输出信号，基于多参考点最小二乘复频域方法（Polyreference least-squares complex frequency-domain method，PolyMAX）^[20]对输入输出信号进行分析处理，进而获得系统模态参数.在导管表面从端部起每隔50mm布置一个测点，共布置11个测点，如图3所示.为确保试验结果真实有效，每次敲击均由同一名试验人员依次对测点进行横向激励并进行5次有效平均.

1.2 试验结果分析

基于上述试验测试方案，得出管道模型前5阶模态频率为128.76Hz、218.62Hz、409.37Hz、651.56Hz、1049.762Hz.试验模型的前5阶模态振型如图4所示.

管道模型整体在-z方向力锤激励下发生横向振动，随着模态阶次的增加，振型中波数和节点数随之增加.由第二阶振型、第四阶振型和第五阶振型可以明显看出，管道端部存在较大的位移.因此，建模时不能仅将导管端部视为固定端约束，需要考虑到管道端部连接结构对试验模型整体模态的影响.

2 管道有限元建模及分析

针对直管管道试验模型，采用SOLIDWORKS三维建模软件和ANSYS有限元软件进行了有限元建模、网格划分和动力学特性分析.在ANSYS中建立的有限元模型如图5所示.导管部分采用六面体单元（HEXA）划分网格，两端连接组件及支座部分采用四边形单元（QUAD4）划分网格，共25309个单元.根据接触设置分析算例^[21]并结合工程实际，将与直通管接头外圆柱面接触的外套螺母、锁紧螺母、堵头等组件与直通管接头之间的接触方式设置为绑定接触，其余接触均设置为粗糙接触，支座底面设置为固定端约束面.

使用ANSYS对有限元模型进行模态分析，提取前5阶模态参数结果，模态频率的对比结果如表2所示.有限元分析结果与实测模态结果中，模态频率误差在20%以内，各阶振型无较大差异，具有一定的相关性.同时，实测模态频率与有限元分析结果相比最大误差达到8.60%，各阶误差均超过5%，超出工程实际中可接受的范围，有必要进行有限元模型修正.

为探究支座及管接头组件对有限元模型的影响，建立了去除全部管接头组件的有限元模型A和仅去除支座的有限元模型B并分别进行模态分析，模态频率结果如表2所示.模型A的最大误差达到76.58%，模型B的最大误差达到93.75%，误差远超过工程实际中可接受的范围.模型A和模型B与实际结构之间存在振型不匹配现象.相比于模型A和模型B，完整模型的模态参数更加接近于实测模态参数.可见支座和管接头连接组件对试验模型整体模态参数具有较大影响，建模时应予以保留.

3 管道有限元模型修正

基于Kriging模型的有限元模型修正流程如图6所示.可大致概括为3个主要流程：①将模态试验实测数据与有限元模型模拟结果进行自由度匹配与相关性分析，确保有限元模型与实际结构具有一定的相关性.其中，若有限元分析频率与实测模态频率的误差在20%以内且有限元分析振型与实测模态振型之间无较大差异，则视为满足相关性条件.②Kriging模型的生成.首先选取合适的待修正参数并进行试验设计，使得在一定的样本空间内最大化获取修正空间内信息.常用的试验设计方法有正交设计、均匀设计、D最优设计、中心复合设计抽样、拉丁超立方抽样等.之后基于训练样本拟合Kriging模型并检验模型精度直至符合精度要求.③针对目标函数的修正.在所拟合的Kriging模型的基础上，构建目标函数，之后引入如遗传算法、退火模拟算法、粒子群算法等优化算法对目标函数进行迭代处理，获取参数修正值.

3.1 参数筛选与试验设计

根据工程经验，初选结构中部分组件的材料参数和几何参数作为待修正参数.管道模型中选择的修正对象有导管密度、导管弹性模量、支座密度、支座弹性模量、外套螺母密度、外套螺母弹性模量等6个材料参数及导管内径、直通管接头长度、垫圈厚度、外套螺母长度等4个几何参数共10个参数.基于斯皮尔曼（Spearman）等级相关系数^[22]，分析各初选几何参数和材料参数对前5阶模态频率的灵敏度.选取显著性水平为0.05，对各待修正参数设置上下2%的微小变化，如表3所示.

计算前五阶频率对各参数的灵敏度，如图7所示.结合工程经验，剔除灵敏度相对较小的参数，选取导管内径d、导管密度ρ₁、导管弹性模量E₁和支座弹性模量E₂四个参数作为待修正参数.

试验设计方法采用拉丁超立方抽样方法，共选取了25个样本点，如表4所示.拉丁超立方抽样是一种基于分层抽样思想的方法，将变量的概率分布函数分成数个互不重叠的等概率子区间，然后在每个子区间内进行独立随机采样，确保能在每个子区间内精确采样^[23].

3.2 Kriging模型构造

Kriging模型是一种基于随机过程的半参数化模型，由线性回归部分和非参数随机部分组成.对于一组给定的输入参数和输出响应$\mathit{X}={\left[{\mathit{x}}_1，{\mathit{x}}_2，\cdots ，{\mathit{x}}_m\right]}^{\mathrm{T}}$和输出响应$\mathit{Y}={\left[y_1，y_2，\cdots ，y_m\right]}^{\mathrm{T}}$，Kriging模型的数学形式可表示为

式中，f_l（x_i）是基函数，β_l是基函数系数，z（x_i）为拟合偏差函数，是对模型局部偏差的近似，$\mathit{f}\left({\mathit{x}}_i\right)={\left[f_1\left({\mathit{x}}_i\right)，f_2\left({\mathit{x}}_i\right)，\cdots ，f_p\left({\mathit{x}}_i\right)\right]}^{\mathrm{T}}$是关于样本x的多项式向量，$\mathit{\beta }={\left[{\beta }_1，{\beta }_2，\cdots ，{\beta }_p\right]}^{\mathrm{T}}$是线性回归系数向量，z（x）服从正态分布N（0，${\sigma }_z^2$）并具有非零协方差，其协方差表达式为

式中，z（x_i），z（x_j）为任意两个样本点x_i和x_j的局部偏差；R是由R_ij（x_i，x_j）组成的相关函数矩阵.R_ij（x_i，x_j）是一种可选定的相关函数，表示任意两个样本点间的空间相关性.本文选取工程中常用的高斯函数作为相关函数，其形式如下

式中，$x_i^k$和$x_j^k$分别为样本点x_i和x_j的第k个分量；n为设计变量的数目；θ_k为未知的相关系数，控制不同维度上的相关衰减率.

在训练Kriging模型中的参数时，通常使用最大似然估计.样本点的似然函数可表示为：

式中|R|是相关函数矩阵R的行列式，F为样本点处f（x）向量值组成的矩阵.β和${\sigma }_z^2$的最小二乘估计可以表示为

将β和${\sigma }_z^2$最小二乘估计代入样本点的似然函数，忽略常数项并取对数，即可得到最大似然函数的对数形式为

${\sigma }_z^2$和R均为θ_k的函数，通过优化算法求得式（7）的最大值，即可确定θ_k的值.

对于任意待测点x₀，最终响应的最佳线性无偏预测结果可表示为

式中

${\mathit{r}}^{\mathrm{T}}\left({\mathit{x}}_0\right)=\left[R\left(x_0,x_1\right),R\left(x_0,x_2\right),\cdots ,R\left(x_0,x_n\right)\right]$

基于上述Kriging模型理论，采用表4所示的拉丁超立方设计中的样本，拟合管道模型模态频率与四个设计变量之间的函数关系，绘制Kriging响应曲面，部分结果如图8所示.

3.3 Kriging模型拟合精度分析

为确定该模型能否代替有限元模型进行后续结构响应值的计算分析，对Kriging模型进行精度校验.常用的精度校验评价指标有：决定系数R²、修正可决系数R²_adj、均方根误差RMSE等^[24].本文选用决定系数和均方根误差RMSE两种检验标准，如式（9）和式（10）所示.

式中，y_i为设计点处有限元计算值；${\widehat{y}}_i$为响应面计算值；$\stackrel{-}{y}$为有限元计算结果的平均值，K为样本点总数.R²→1.表明代理模型与原模型相似度高，RMSE→0表示代理模型误差小.Kriging模型精度检验结果如表5所示.结果显示前五阶决定系数均为1且RMSE都接近于0，这表明所拟合的Kriging代理模型误差小、拟合精度高，可以作为原有限元模型的替代模型.

3.4 模型修正及结果分析

本文以修正后的管道模型有限元模型与实际管道模型响应尽可能接近为总体目标，利用模态试验所得模态频率构建的目标函数为

$minF\left(x\right)=\sum _{i=1}^m{\left(\frac{{\tilde{f}}_i-f_i}{{\tilde{f}}_i}\right)}^2$

式中，${\tilde{f}}_i$为第i阶实测模态频率；f_i为第i阶计算模态频率；x_j为第j阶修正参数，$x_{\text{low}}^j$和$x_{up}^j$为修正参数的最小值和最大值；m为参与修正的样本个数；n为修正参数的个数.

基于拟合的Kriging模型，本文选用MOGA遗传算法对目标函数进行求解优化.MOGA遗传算法是基于受控精英主义策略的快速非支配排序遗传算法（NSGA-Ⅱ）^[25]的变体，它可以允许多个目标和约束，能获得分布均匀且兼具多样性的Pareto解集.基于MOGA遗传算法进行优化求解，最初生成4000个样本，每次迭代生成800个样本，并在最多20次迭代中筛选出3个数据较好的样本点，在6937次评估后收敛，得到最终的Kriging模型修正结果.为验证本文模型修正方法的修正精度和修正效率，将Kriging模型修正方法与传统模型修正方法灵敏度方法进行对比，两种方法获得的管道模型修正前后的结构参数和模态频率的结果分别如表6和表7所示.

从表中可以发现：（1）通过灵敏度方法和Kriging模型修正方法修正后，管道各参数修正率较为接近，可以相互证实所提修正方法的可行性.（2）采用Kriging模型修正方法，修正前有限元模型与试验值之间各阶误差均超过工程允许的5%的误差范围，修正后将有限元模型误差同修正前原始误差相比，最大误差由8.60%变为1.43%，各阶相对误差均在5%以内，已满足工程实际的计算需求.（3）灵敏度分析方法得到的最大相对误差达到2.56%，修正精度略逊于Kriging模型.同时，在同一系统上进行计算时，采用灵敏度分析方法进行修正需耗时2h，而采用Kriging模型进行修正只需5min，可见Kriging模型在计算效率方面具有显著的优越性.

图9给出了重点修正的模态阶次所对应的模态振型.从结果中可以看出，通过Kriging模型修正后的有限元计算振型与实测振型具有较好的一致性，进一步说明了采用Kriging模型修正方法对管道结构进行有限元模型修正的适用性.

4 结论

本文以某航空发动机管道模型为研究对象，采用力锤法进行了结构的模态试验并得出模态参数，通过ANSYS有限元分析软件建立有限元模型并进行模态分析，对比实测模态数据与有限元分析结果以确定有限元模型的精确度.结合工程经验与灵敏度分析方法筛选出待修正参数，采用拉丁超立方抽样设计进行试验设计，拟合出模态频率与结构参数之间的Kriging模型，利用MOGA遗传算法完成了有限元模型修正.根据模态试验与有限元模型修正的结果，可得出以下结论：

（1）通过对管道模型的模态试验结果进行分析，得到了管道模型的模态参数.其中，管道振型中端部有较大位移，管接头连接组件和支座对模态参数有较大影响，建模时应予以保留.

（2）利用Kriging模型代替复杂的管道有限元模型进行有限元模型修正可以明显提高计算精度和效率.同修正前的有限元模型相比，所得各阶模态频率最大误差由8.60%降为1.63%，各阶相对误差均在5%以内，满足工程实际的计算需求.

（3）对于管道结构，Kriging模型修正方法同灵敏度分析方法相比具有更高的修正精度和修正效率.灵敏度分析方法最大误差为2.56%，Kriging模型修正方法最大误差为1.43%，Kriging模型修正方法修正精度优于灵敏度分析方法.修正后的有限元模型能起到更好的预测效果，可为此类管系结构模型修正提供参考.

参考文献