一种卡箍-管路系统频响函数的不确定性分析方法与流程

本发明属于机械动力学技术领域，尤其涉及一种卡箍-管路系统频响函数的不确定性分析方法。

背景技术：

在航空发动机中，管路系统是由很多复杂管道构成的，这些管道通过卡箍固定在航空发动机上。卡箍-管路系统是润滑油、燃料油、液压油等流体输送的重要部件，在航空航天工业中得到了广泛的应用。管路系统的故障已成为引起航空发动机故障的主要问题之一。如果能根据外部条件估计出卡箍-管路系统振动的不确定区间，就能减少这种故障。

不确定性方法大致可分为三大类：概率不确定性方法、模糊不确定性方法和区间不确定性方法。

概率不确定方法的研究对象是服从一定分布特征的不确定参数，其中多项式混沌展开法是最常用的概率不确定方法之一。xiu等提出了广义多项式混沌方法，这种方法可用于分析服从不同分布状态的随机输入参数。wan等在多项式混沌展开法的基础上提出了一种考虑多因素影响的方法，该方法适用于具有大量随机参数的问题。xu等人将多项式混沌展开法与稀疏网格方法相结合，在保证计算精度的前提下提高计算效率。manan和cooper利用非侵入的多项式混沌展开法预测了各参数的不确定性回归模型，这些参数构成了频响函数曲线。

模糊不确定方法的目标是边界或区域模糊的不确定参数，需要用模糊值来描述。klimke等采用稀疏网格插值方法得到模糊值。puig等人提出了一种不随时间变化的模糊参数线性系统的求解方法。

区间不确定法的目标是边界清晰的不确定参数，但难以描述其分布规律。常用的方法有泰勒级数展开法、切比雪夫多项式逼近法等等。邱等人利用泰勒级数展开法得到了桁架结构动力响应的不确定区间。wu等提出了基于切比雪夫多项式逼近法的chebyshev包含函数，并将该方法应用于非线性系统的不确定性区间分析中。wu等基于稀疏网格方法，提高了切比雪夫包含函数的效率和精度。fu等分别将切比雪夫多项式逼近法和勒让德多项式逼近法应用于裂纹转子的不确定区间估计。fu等人提出了一种基于切比雪夫多项式逼近法的区间估计代理模型。muscolino等考虑杨氏模量和杆长的不确定性，利用有理级数展开法得到了24杆桁架结构极限承载力的不确定区间。

综上，不确定性分析方法在许多研究领域都有应用，但对卡箍-管路系统的不确定性的研究相对较少。

技术实现要素：

(一)要解决的技术问题

针对现有存在的技术问题，本发明提供一种卡箍-管路系统频响函数的不确定性分析方法，分析结果准确且效率高。

(二)技术方案

本发明提供一种卡箍-管路系统频响函数的不确定性分析方法，包括如下步骤：

a1、基于梁单元模型和弹簧模型，建立卡箍-管路系统的有限元模型；

a2、将预先获得的卡箍刚度的不确定区间作为区间变量，应用于所述卡箍-管路系统中的有限元模型中，利用切比雪夫多项式逼近法得到卡箍-管路系统的频响函数的不确定区间。

进一步地，所述步骤a1包括：

a11、基于梁单元模型，得到管路的刚度矩阵kp和质量矩阵m；

其中，梁单元模型中包括梁单元的刚度矩阵和梁单元的质量矩阵

a12、基于弹簧模型，得到卡箍的刚度矩阵kh；

a13、根据所述管路的刚度矩阵kp和卡箍的刚度矩阵kh，得到卡箍-管路系统的刚度矩阵k；

a14、根据所述管路的质量矩阵m和卡箍-管路系统的刚度矩阵k，得到卡箍-管路系统的有限元模型。

进一步地，所述步骤a13中，卡箍-管路系统的刚度矩阵k＝管路的刚度矩阵kp+卡箍的刚度矩阵kh。

进一步地，所述管路的刚度矩阵kp满足如下公式：

式中：v为梁单元的数量；tw为转换矩阵：

tw＝diag[tptptptp]

式中：τ为全局坐标轴和局部坐标轴的余弦值，τxxb＝τzzb＝1,τzxb＝τxzb＝0；

所述管路的质量矩阵m满足如下公式：

进一步地，所述梁单元模型中包括梁单元的刚度矩阵和梁单元的质量矩阵

梁单元的刚度矩阵满足如下公式：

式中：ei为抗弯刚度；l为梁单元长度；

梁单元的质量矩阵满足如下公式：

式中：ρ为梁单元密度；a为梁单元横截面面积。

进一步地，所述卡箍-管路系统的有限元模型满足如下公式：

式中：c为卡箍-管路系统的阻尼矩阵；f为外力矢量；u分别为卡箍-管路系统的加速度、速度和位移矢量。

进一步地，所述卡箍刚度包括卡箍在z方向的平动刚度kz和卡箍在x方向的扭转刚度kθx。

进一步地，所述卡箍在z方向的平动刚度kz满足如下公式：

kz＝δfz/δz；

所述卡箍在x方向的扭转刚度kθx满足如下公式：

kθx＝δtx/δθx；

式中：δfz和δz分别为卡箍在z方向上的加载力和位移；δtx和δθx分别为卡箍在z方向上的扭矩和角位移。

进一步地，所述步骤a2包括：

a21、通过实验预先获得卡箍刚度的不确定区间：对卡箍的灵敏度进行分析，确定有效刚度，根据所述有效刚度获取刚度不确定区间；

a22、将所述刚度不确定区间作为不确定区间变量，利用切比雪夫多项式逼近法得到卡箍-管路系统的频响函数的不确定区间。

(三)有益效果

本发明提供的方法采用切比雪夫逼近法计算卡箍-管路系统频响函数的不确定性区间，计算结果准确性高、效率高，同时引入容差百分比量化实验和仿真结果的差异。

附图说明

图1为本发明提供方法的流程图；

图2为本发明中卡箍-管路系统的物理尺寸图；

图3为本发明中欧拉梁模型的示意图；

图4为本发明中卡箍-管路系统的有限元模型示意图；

图5为本发明中测量卡箍刚度的试验台示意图；

图6为本发明中卡箍刚度的实验结果与仿真结果对比图；

图7为本发明中平动刚度测试仪示意图；

图8为本发明中扭转刚度测试仪示意图；

图9为本发明实施例中卡箍刚度区间的实验测试流程图；

图10为本发明实施例中卡箍1在第一、二、三次装配的平动刚度示意图；

图11为本发明实施例中卡箍1在第一、二、三次装配的扭转刚度示意图；

图12为本发明实施例中卡箍2在第一、二、三次装配的平动刚度示意图；

图13为本发明实施例中卡箍2在第一、二、三次装配的扭转刚度示意图；

图14为本发明中卡箍-管路系统的频响函数不确定性区间结果示意图；

图15为本发明中两侧螺栓相同拧紧力矩均为3、5、7、9、11n·m时仿真和实验的频响不确定性区间对比图；

图16为本发明中两侧螺栓不同拧紧力矩的情况下，仿真和实验的频响不确定性区间对比图。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

本发明提供一种卡箍-管路系统频响函数的不确定性分析方法，如图1所示，包括如下步骤：

a1、基于梁单元模型和弹簧模型，建立卡箍-管路系统的有限元模型。

卡箍-管路系统的物理尺寸图如图2所示，参数如表1所示。每个节点有4个自由度，分别为x方向和z方向的平移位移和角位移，采用如图3所示的欧拉梁模型建立管路，得到如图4所示的卡箍-管路系统的有限元模型。

表1卡箍-管路系统的参数

梁单元的自由度d如公式(1)所示：

d＝[xa，,za,θxa,θza,xb，,zb,θxb,θzb](1)

式中：xa、za分别为梁单元节点a在x、z方向的平动位移；θxa、θza分别为梁单元节点a在x、z方向的角位移；xb、zb分别为梁单元节点b在x、z方向的平动位移；θxb、θzb分别为梁单元节点b在x、z方向的角位移。

梁单元的刚度矩阵如公式(2)所示：

式中：ei为抗弯刚度；l为梁单元长度。

梁单元的质量矩阵如公式(3)所示：

式中：ρ为梁单元密度；a为梁单元横截面面积。

将每个卡箍等效为4个平动弹簧和4个扭转弹簧，每个弹簧的刚度是卡箍的一半，卡箍的刚度矩阵kh如公式(4)所示：

式中：kx为卡箍在x方向的平动刚度；kz为卡箍在z方向的平动刚度；kxθ为卡箍在x方向的扭转刚度；kzθ为卡箍在z方向的扭转刚度。

在本实施例中，利用如图5所示的实验台，采用光传感器采集振动信号，敲击点节点坐标为(0,0.075,0)，加速度传感器节点坐标为(0,0.225,0)，对卡箍刚度进行多次测量，对试验结果进行最小二乘拟合得到拟合刚度值。然后对拟合刚度值进行线性化，得到等效平均刚度。等效刚度为：kx＝3.53mn/m、kz＝4.98mn/m、kxθ＝58.84n·m/rad、kzθ＝27.37n·m/rad。仿真结果与实验结果对比如图6所示：除第四阶固有频率外，仿真结果与实验结果吻合较好，第四阶固有频率误差为9.1％。

管路的刚度矩阵kp和质量矩阵m分别如公式(5)、(6)所示：

式中：v为梁单元的数量；tw为转换矩阵，如公式(7)所示：

式中：τ为全局坐标轴和局部坐标轴的余弦值，由于梁单元的全局坐标与局部坐标相一致，因此τxxb＝τzzb＝1,τzxb＝τxzb＝0。

卡箍-管路系统的刚度矩阵k如公式(8)所示：

k＝kp+kh(8)

卡箍-管路系统的有限元模型满足公式(9)：

式中：c为卡箍-管路系统的阻尼矩阵；f为外力矢量；u分别为卡箍-管路系统的加速度、速度和位移矢量。

其中，卡箍-管路系统的阻尼矩阵c由瑞利阻尼确定，满足公式(10)：

式中：α、β分别为管路质量矩阵m和卡箍-管路系统刚度矩阵k的比例系数；ω1、ω2分别为卡箍-管路系统的一阶和二阶圆频率，ω1＝2πf1，ω2＝2πf2，f1、f2分别为卡箍-管路系统的一阶和二阶固有频率；ξ1、ξ2分别为卡箍-管路系统的一阶和二阶模态阻尼比，ξ1＝0.02，ξ2＝0.03。

卡箍-管路系统的位移传递函数hd(ω)由拉普拉斯变换确定，如公式(11)所示：

hd(ω)＝(-mω²+jωc+k)^-1(11)

式中：j为虚数；ω为固有频率。

根据卡箍-管路系统振型的正交性，卡箍-管路系统的位移传递函数hd(ω)还满足公式(12)：

式中：ai、分别为第i阶的特征向量和其转置；ξi为第i阶的特征值；ωi为第i阶的固有频率；为模态截断数。

根据拉普拉斯变换的性质，卡箍-管路系统的加速度传递函数ha(ω)满足公式(13):

激振点和测量点的加速度频响函数满足公式(14)：

a2、将预先获得的卡箍刚度的不确定区间作为区间变量，应用于所述卡箍-管路系统中的有限元模型中，利用切比雪夫多项式逼近法得到卡箍-管路系统的频响函数的不确定区间。

a21、通过实验预先获得卡箍刚度的不确定区间。

首先，对卡箍的灵敏度进行分析，确定有效刚度。

当参数r在m、c、k时，参数r对卡箍-管路系统的一阶固有频率的灵敏度可表示为如公式(15)所示：

计算结果如公式(16)-(18)所示：

式中：mpq、cpq和kpq分别为m、c和k中第p行第q列的元素。

经实验，kx、kz、kxθ和kzθ的敏感度分别为0、0.15、0.04和0。可以看出，kz对卡箍-管路系统的一阶固有频率影响最大，其次是kxθ，kx和kzθ对卡箍-管路系统的一阶固有频率没有影响。这是因为kx和kzθ的作用方向与z方向垂直，所以对kx和kzθ与z方向上的频响函数，因此，当分析z方向上的频响函数时，kx和kzθ可以忽略。

其次，根据有效刚度获取刚度不确定区间。

如图7、8所示，分别为平动刚度测试仪和扭转刚度测试仪。卡箍在z方向的平动刚度kz和卡箍在x方向的扭转刚度kxθ分别满足公式(19)、(20)：

kz＝δfz/δz(19)

kθx＝δtx/δθx(20)

式中：δfz和δz分别为卡箍在z方向上的加载力和位移；δtx和δθx分别为卡箍在z方向上的扭矩和角位移。

在本实施例中，测试夹具刚度区间(包括平动刚度kz和扭转刚度kxθ)时，为考虑不同螺栓拧紧力矩对卡箍刚度区间的影响，选择五种螺栓拧紧力矩(3n·m、5n·m、7n·m、9n·m、11n·m)进行测试。由于卡箍重新装配后刚度区间会发生变化，因此在每一螺栓拧紧力矩下需要进行7次重复试验。当一组实验完成后，取下卡箍，重新装配卡箍-管路系统，收集三组刚度区间的数据，得到最终的卡箍刚度区间。实验测试过程如图9所示。

本实施例采用统计中常用的箱形图来描述卡箍的刚度区间。如图10、11所示，从左至右分别为卡箍1在第一、二、三次装配的平动刚度示意图和扭转刚度示意图；如图12、13所示，从左至右分别为卡箍2在第一、二、三次装配的平动刚度示意图和扭转刚度示意图。由图可知，随着螺栓拧紧力矩的增加，平动刚度kz和扭转刚度kθx的平均值近似线性增加，扭转刚度kθx平均值变化率大于平动刚度kz的平均值变化率。

如表2所示，为卡箍1、卡箍2的刚度测试结果。

表2卡箍1、卡箍2的刚度测试结果

由表2可以看出，刚度测试结果具有较大的分散性。本实施例中，最大平动刚度是最小平动刚度的1.6倍以上(卡箍1的拧紧力矩为7n·m时)，最大扭转刚度是最小扭转刚度的1.4倍以上(卡箍2的拧紧力矩为3n·m)。对比两种刚度数据的标准差可以看出：平动刚度kz的分散性要大于的扭转刚度kθx的分散性，而分散性与拧紧力矩之间的关系不明显。

a22、以表2中刚度(平动刚度和扭转刚度)不确定区间作为不确定区间变量，利用切比雪夫多项式逼近法得到卡箍-管路系统的频响函数的不确定区间。

首先，对切比雪夫多项式逼近法进行描述，如公式(21)所示：

b∈[b^l,b^h](21)

式中：b为切比雪夫多项式的区间变量；b^l和b^h是区间变量b的变化范围，当有多个区间变量时，可表示为公式(22)：

切比雪夫多项式hj(x)如公式(23)所示：

式中：n为切比雪夫多项式的展开阶数；

将切比雪夫多项式hj(x)与权函数正交，得到公式(24)：

式中：hy(x)、hz(x)分别为y阶切比雪夫多项式和z阶切比雪夫多项式。

未知函数g(x)的一维n阶切比雪夫多项式i(x)可近似表示为如公式(25)所示：

式中：b0和bj为多项式的展开系数，bj如公式(26)所示；n为多项式的展开阶数；展开系数和展开阶数由未知函数g(x)确定。

由于未知函数g(x)的表达式未知，用高斯-切比雪夫插值积分表示公式(26)，如公式(27)所示：

式中：xp为插值点，q为插值点个数；bp为高斯积分系数，g(xp)为未知函数g(x)在插值点xp处的值；hj(xp)为切比雪夫多项式hj(x)在插值点xp处的值。结合上述公式，得到bj表达式，如公式(28)所示：

根据卡箍-管路系统的有限元模型中插值点xp处的频响函数值和切比雪夫多项式i(x)，即可确定bj。

由于切比雪夫正交逼近的区间为[-1,1]，因此需要将实际区间转化为[-1，1]。在计算频响函数的变化范围时，需要注意极值点和边界点，极值点和边界点可以通过偏导运算得到。由于单个区间变量互不干扰，将各个变量的最大值和最小值代入运动方程[25，28]，即可得到频响函数的不确定区间。

本实施例中，在两侧螺栓相同拧紧力矩的情况下，以平动刚度kz和扭转刚度kθx的区间(如表2所示)作为卡箍刚度不确定区间变量，采用切比雪夫多项式拟合法计算卡箍-管路系统的频响函数不确定性区间。

以卡箍-管路系统的第一、二阶固有频率为感兴趣的频率范围进行分析，卡箍-管路系统的频响函数不确定性区间结果如图14所示。由图可知，随着螺栓拧紧力矩的增大，卡箍-管路系统的一阶、二阶固有频率的区间值逐渐增大，不确定性区间的宽度减小，这是一种频移现象。

放大3n·m的拧紧力矩的卡箍-管路系统的一阶、二阶固有频率的不确定性区间，可见频响振幅随频率的增加而逐渐增大。当振幅最大时，卡箍1和卡箍2的平动刚度kz分别为5.0621×10⁶n/m和4.8761×10⁶n/m。卡箍1和卡箍2的扭转刚度kθx分别为48.8876n·m/rad和49.8705n·m/rad。

为了验证切比雪夫多项式逼近法的效率和有效性，采用扫描法作为比较方法(扫描法是不确定度分析中常用的一种方法，通过将变化范围划分为等步长采样点来计算区间不确定性)，比较145hz和457hz的频响振幅范围，结果如表3所示。

表3扫描法和切比雪夫多项式逼近法的对比结果

由表3可知，扫描法和切比雪夫多项式逼近法计算的频响振幅范围基本相同，但切比雪夫多项式逼近法耗时更短。

如图15所示，为两侧螺栓相同拧紧力矩均为3、5、7、9、11n·m时仿真和实验的频响不确定性区间对比图，对比数据如表4所示。

表4两侧螺栓为相同拧紧力矩时的实验和仿真结果对比

测量卡箍-管路系统的频响时，在输入脉冲信号之前，加速度传感器存在微弱信号，这是测试系统的硬件问题。因此，低频区(0hz-50hz)将出现弱峰。为了消除这种干扰，对比仿真结果与实验结果进行时，截取50hz-800hz频段进行分析。

由图可知，仿真的频响不确定性区间基本涵盖了频响不确定性区间。当螺栓拧紧力矩较小时，卡箍-管路系统的频响分布较分散，当螺栓拧紧力矩较大时，卡箍-管路系统的频响分布较聚集。这是因为当螺栓拧紧力矩小(或大)时，夹紧力小(或大)。

由表4可知，一阶固有频率下边界的容差百分比大于上边界的，而二阶固有频率的容差百分比则相反。说明，实验结果的一阶固有频率接近仿真的上边界，实验结果的二阶固有频率接近仿真的下边界。区间宽度随拧紧力矩的增大而减小。中心频率的容差百分比不超过5％，频率区间的容差百分比小于10％，验证了仿真结果的正确性。

在两侧螺栓不同拧紧力矩的情况下，采用切比雪夫多项式逼近法计算卡箍-管路系统的频响函数不确定性区间，其中将左右侧螺栓的拧紧力矩分别定义为tl和tr，结果如图16所示，对比数据如表5所示。tl＝9,tr＝5和tl＝11,tr＝3的情况下，左右两侧螺栓的平均拧紧力矩为7n·m。将两种情况与两侧螺栓相同的拧紧力矩7n·m的不确定性区间情况进行比较，可以看出，以上三种情况的不确定性区间基本相同。因此，当螺栓两侧的拧紧力矩不同时，螺栓两侧的拧紧力矩可以近似的平均。当两侧螺栓的不同拧紧力矩等效为两侧螺栓的相同拧紧力矩时，固有频率的不确定性区间存在左偏现象。这可能是因为较小拧紧力矩的一侧比较大拧紧力矩的一侧对系统频率影响要大。因此，这种情况相对于两侧螺栓均为7n·m的拧紧力矩的频率要降低，会出现频率向左侧偏移的现象。

表5两侧螺栓为不同拧紧力矩时的实验和仿真结果对比

由表5可知，实验结果的一阶固有频率接近仿真结果上边界，实验结果的二阶固有频率接近仿真结果下边界，且容差百分比大于两侧螺栓相同的拧紧力矩的情况。

综合上述实验，得到以下结论：

(1)随着螺栓的拧紧力矩的增加，平动刚度kz和扭转刚度kθx的平均值都近似线性的增加。扭转刚度kθx平均值的变化率大于平动刚度kz平均值的变化率。实验得到的刚度数据具有较大的分散性，就本发明实施例中，最大平动刚度是最小平动刚度的1.6倍以上，最大扭转刚度是最小扭转刚度的1.4倍以上，平动刚度kz的分散性大于扭转刚度kθx。

(2)当螺栓的拧紧力矩较小时，卡箍-管路系统的频响函数分布较分散；当螺栓的拧紧力矩较大时，卡箍-管路系统的频响函数分布较为集中。实验测试的频率区间表明，实验结果的一阶固有频率接近仿真的上边界，实验结果的二阶固有频率接近仿真的下边界。

(3)分析两侧螺栓在不同拧紧力矩下的频响函数的不确定性时，左侧螺栓的拧紧力矩(tl)和右侧螺栓的拧紧力矩(tr)可以被近似的平均，这种情况的固有频率的不确定区间会存在左偏现象。当tl和tr的差异大时，这种现象更为明显。随着螺栓拧紧力矩的增大，卡箍-管路系统的一阶固有频率和二阶固有频率的区间值逐渐增大，不确定区间的宽度逐渐减小。在卡箍刚度不确定的情况下，系统的固有频率会发生频移现象。

以上结合具体实施例描述了本发明的技术原理，这些描述只是为了解释本发明的原理，不能以任何方式解释为对本发明保护范围的限制。基于此处解释，本领域的技术人员不需要付出创造性的劳动即可联想到本发明的其它具体实施方式，这些方式都将落入本发明的保护范围之内。