0   引言

脉冲燃烧风洞的建造和运行成本较低，能够以不同尺度、方式进行基本流动和工程问题的研究，是目前国内开展大尺度超燃冲压发动机和机体推进一体化高超声速飞行器试验研究的主要设备之一^[1]。然而，目前脉冲燃烧风洞的工作时间仅能达到300ms左右^[2]，模型振动在试验时间内来不及衰减，测力天平得到的信号含有较大的振动成分，传统方法难以获得良好的结果。如果能先获得试验系统的动态输入输出特性，就可以用天平的输出信号反算出模型所受的载荷历程而不需要模型达到稳定，这种载荷辨识方法已经在激波风洞内得到应用^[3-4]。该方法能够得到载荷的时变过程，这对于飞行器带动力试验尤为重要。

载荷辨识的前提是已知系统的动态输入输出特性，动态标定试验是获得这一特性的可靠途径。昆士兰大学的Abdel-Jawad等对超高速膨胀管中的多分量应力波天平进行了动态标定，通过加载一系列标定载荷，得到了模型对分布载荷的脉冲响应矩阵^[7]。牛津大学的L.J.Doherty等在Stalker管风洞中对一带动力飞行器进行测力试验时，用力锤敲击方法对一台三分量应力波天平进行了动态标定。中国空气动力研究与发展中心的王锋等采用瞬时卸载方式对系统进行激励，分别用阶跃函数和斜坡阶跃函数描述卸载过程，对某单分量天平进行了动态标定。本文针对一个六分量测力试验系统，给出了一种动态标定方法并对这种方法的准确性进行了验证，为后续开展载荷辨识工作提供必要条件^[8]。

1   测力试验系统介绍

试验系统包括蒙皮框架结构的试验模型、天平以及固定底座3部分，如图 1所示。模型包含5个面，受载荷冲击时通过框架将振动传递至天平浮动框，天平固定框与底座相连，底座固定于地面。

图  1  试验系统示意图

Fig.  1  Schematic of the testing system

下载: 全尺寸图片 幻灯片

六分量天平由沿轴向布置的前后2个环形结构组合而成，如图 2所示，上部为浮动框，下部为固定框。浮动框与固定框之间是测量元件，其上贴有应变片，可测量应变并通过后续电路将应变信号输出。

图  2  六分量天平

Fig.  2  Six component balance

下载: 全尺寸图片 幻灯片

2   测力试验系统模态分析

模态分析是结构动力学分析的基础，系统在载荷作用下发生振动，可视为各阶模态振动的叠加^[9]。通过模态分析可近似了解系统的基本动力学特征，为开展动态标定提供重要的参考信息。利用有限元软件计算得到系统的前6阶模态，如表 1所示。

表  1  试验系统的前6阶模态

Table  1  The first six order modes of the testing system

Order	1	2	3	4	5	6
f/Hz	5.91	9.55	12.48	16.13	17.03	23.30
Mode
Description	around x axis	around y axis	around z axis	along z axis	along x axis	along y axis

下载: 导出CSV 

| 显示表格

3   测力试验系统动态标定方法

3.1   动态标定目的

动态标定的目的是获得单位脉冲载荷作用下系统的响应函数(UIRF)。这个单位脉冲载荷是指由模型表面等效到天平中心点(或模型质心)后的单位脉冲载荷。试验模型本身尽管整体刚度较大，但也是一个弹性体，在模型不同位置处加载对应的UIRF不同。这里近似认为模型同一个面上各位置对应的UIRF相同，因此可以在1个面上施加一系列集中载荷，通过输入载荷和输出应变求得该面对应的UIRF。在5个面上分别进行上述操作，然后将得到的5个UIRF加权即可得到整个试验系统的UIRF。加权系数与试验状态下对应面上的平均压力大小成正比^[10]。各面上的平均压力变化，加权系数要相应调整。仿真计算中，平均压力由事先的CFD计算近似得到；风洞试验时，在各面上布置测压传感器，每个时刻都可以测出1组压力值进而算得1组加权系数，实际上同一车试验中各面上平均压力的比例关系随时间变化不大，近似计算也可将不同时刻的加权系数取平均用于后续计算。

定义模型的5个表面为e₁~e₅，在模型的面e(e=1, 2, 3, 4, 5) 上施加一时变载荷，其相对于天平中心有6个分量，即x、y、z方向的力和这3个方向的力矩，分别记为F_x(t)，F_y(t)，F_z(t)，M_x(t)，M_y(t)，M_z(t)。天平输出Y有6个通道，将第k(k=1, 2, …, 6) 个通道的输出记为y_k(t)。测力系统的输入输出关系为

其中，“*”表示卷积。第k行表示第k通道的输入输出关系：在面e上施加载荷，载荷相对天平中心分别为x、y、z方向单位脉冲力或力矩时，第k通道输出分别为I_xek、I_yek、I_zek、I_mxek、I_myek、I_mzek。令

则(1) 式可写为

I_e即为面e对应的单位脉冲响应函数，是一个6×6的矩阵，整个系统的单位脉冲响应函数I为，

其中，α_e为面e对应的加权系数。对于某个时间点，I是一个6×6的矩阵，动态标定的目的即是获得I。

3.2   动态标定过程

本文仅针对模型的1个表面，以天平的1个通道为例介绍动态标定方法，即只求出I_e的第k行：I_ek。其他各面、各通道的动态标定方法是相同的。因此只需要考虑(1) 式中的第k行

本文中仅针对e₅面(e=5)，第1通道(k=1) 求解对应的I₅₁，为了简洁，下文中略去下标e和k，将(5) 式、(8) 式分别简写为:

下文中I均指代I₅₁，对于某个时间点为一1×6向量，下面介绍I的求法。

动态标定采用力锤敲击法，在e₅面上安装加载部件(图 1中有放大显示)，为力锤敲击提供5个不同方向的加载面。加载部件的加载面面积远大于力锤锤头的面积，这样便于操作并有利于保证敲击方向的准确性。定义加载方向为d₁~d₅，如图 3所示。

图  3  模型表面加载点与方向示意图

Fig.  3  Schematic of load points and directions

下载: 全尺寸图片 幻灯片

在点p(p=1, 2, 3, …)处沿方向d(d=1, 2, 3, 4, 5) 敲击，敲击力要足够大以保证天平有明显输出但又不能超出量程，力锤可记录其输入载荷大小的时间历程f(t)，结合加载位置和方向计算出其相对于天平中心的载荷向量F

其中

c_dp为单位载荷相对于天平中心的载荷矢量，其6个分量不独立，由加载位置和方向决定。给定d、p后，可由简单的几何关系求得c_dp，是已知的。

设输出信号为y_dp(t)，由(10) 式

定义

则

这里，I_dp即为模型表面加载点对应的脉冲响应函数，将在3.3节中给出其确定算法，此处暂认为已知。于是I_dp、c_dp已知，只有I待求。

由(9) 式，I有6个未知分量，需要6个系数不相关的方程组成非齐次线性方程组才可解。因此，需要改变d、p进行6次不相关的敲击。记敲击次数为h(h=1, 2, 3, 4, 5, 6)，并令:

由(14)，(16) 和(17) 式得到

因此，只需C可逆，即c_{dp_1}, …, c_{dp_6}线性无关，便可解得I

3.3   I_dp的求法

3.3.1   参数化I_dp

设有n个时间点t₁…t_j…t_n，那么就存在I_dp(t₁)…I_dp(t_j)…I_dp(t_n)共n个未知数，若直接将f和y_dp(t)代入(15) 式解卷积，问题的病态特征可能使结果误差很大。

考虑结构是线弹性的，可将脉冲响应函数I_dp(t)参数化

参数化将I_dp写成m个不同频率振动叠加的形式。其中，m为选取的模态数量，A_i是第i个模态的幅值，ω_i、γ_i分别是第i个模态的频率和阻尼比。参数化后，未知数变为3m个，通常，n > > m，未知数的个数大大减少。

设实际测量得到的输出信号为y_odp(t)，令R(t)=y_odp(t)-y_dp(t)，定义残差r

求I_dp(t)的问题转化为用I_dp*f拟合y_odp(t)的问题，即确定能使r取到最小值的m及其对应的A_i、ω_i、γ_i这一优化问题。

3.3.2   优化算法

分析3.3.1节中的这一优化问题，待定参数是m阶模态参数(ω₁, A₁, γ₁)…(ω_m, A_m, γ_m)，优化目标是使残差r取最小值。优化过程中，r逐渐减小，给定一个值r_l，以r＜r_l作为终止优化的判据。

阻尼比γ一般在10^-3量级或更小，在振动频率不太高、测量时长较短的情况下对优化结果影响很小，故先不考虑阻尼比，对每个模态仅考虑幅值和频率2个参数，I_dp(t)简化为

求解前m未知，可预估m值，然后对3m个参数同时优化。但这对预估的要求很高，m取小了导致I_dp(t)中必然缺少高阶振动成分；m取大了，则优化参数过多，优化难度较大。

这里设计了一个算法，将整个优化过程分若干次进行，每次应用一定的优化算法只辨识出1个模态对应的成分。第1次优化时，待拟合曲线为y_odp(t)。第i次(i=1, 2…)优化终止时，计算R_i(t)并输出辨识出的第i阶频率ω_i和幅值A_i。第i+1次优化时，以R_i(t)作为待拟合曲线，重复上一步的优化过程，得到ω_i+1和A_i+1。因为每次优化的目标都是使残差尽可能减小，被识别出来的都是当次待拟合曲线中最主要的振动成分，所以各振动成分是按照幅值由大到小的顺序被辨识出来的，即A_i+1≥A_i。当进行到第M次优化后，A_M < < A₁，主要的模态成分已经被辨识出来，若达到了r＜r_l的条件，就可以终止优化，M即为m，输出(ω₁, A₁)…(ω_m, A_m)，并由(20) 式计算I_dp(t)，供后续使用。上述过程如图 4所示。

图  4  优化算法流程示意图

Fig.  4  Flow diagram of the optimization algorithm

下载: 全尺寸图片 幻灯片

在第i次优化中，如图 4虚线框中部分，未知数为ω_i和A_i，待拟合曲线为R_i-1(t)(i=1时为y_odp(t))，目标函数为r_i(ω_i, A_i)。考虑r(ω_i, A_i)是一个多峰函数，使用传统的搜索方法极易陷入局部最优解。遗传算法是一种概率搜索算法，可以从一个种群开始并行搜索，并且不依赖于目标函数的梯度信息，具有很强的全局优化性能^[11-12]，适合提取当前幅值最大的振动成分。这里采用Matlab的遗传算法工具箱(GA)进行优化。收敛准则为本代残差r与上一代残差r_lg的相对差值小于ε，即

此处ε取10^-4。求解前给定初代种群的取值范围，求解中如果某一代满足了收敛准则，就输出结果；如果不满足收敛准则，就通过选择性复制、交叉、变异生成下一代种群，直到满足收敛准则为止。

在上述过程中，若r_l取得较小则会出现ω_i≈ω_j, A_i > > A_j，(i＜j)的情况。这种情况可以通过减小ε和提高y_odp(t)精度而避免，然而实际中，y_odp(t)是输出应变，其中必然会含有噪声，ε取得过小则会降低求解效率。实际上，这里第i次和第j次得到的就是同一个频率的振动成分，可以将对应的参数加权为1组参数作为初值继续优化

前面已经得到了(ω₁, A₁)…(ω_m, A_m)和I_dp(t)。现加入对阻尼比γ的考虑，若继续基于上述的遗传算法进行优化，减小ε和r_l，是可以得到预期的I_dp(t)的，但是求解效率很低。进一步优化时采用Matlab优化工具箱的fminsearch函数，它是基于单纯形算法的优化工具。将前面得到的(ω₁, A₁)…(ω_m, A_m)作为频率和幅值的初值，0作为阻尼比γ的初值。由于这些初值已经很接近最优解，fminsearch函数只需在最优解附近搜索即可，这就避免了其易于陷入局部最优解的缺点，而充分发挥了其收敛速度更快的优势^[15]。经过若干次迭代，残差迅速下降并达到稳定，此时终止优化，输出(ω₁, A₁, γ₁)…(ω_m, A_m, γ_m)，并代入(20) 式计算出I_dp(t)。

4   测力试验系统动态标定方法验证

第2节中介绍了动态标定试验的操作方法和求系统单位脉冲响应函数的算法。本节将利用有限元软件模拟试验过程，利用上述算法求得I，并检验I的准确性。

4.1   模拟动态标定试验

在Ansys Workbench中进行瞬态动力学分析，模拟上述的动态标定过程。分别在点p₁(d₁, d₂, d₄方向)、点p₂(d₁, d₂方向)、点p₃(d₁方向)模拟力锤加载，各加载点位置如图 3所示。由几何关系分别算得这6次加载的载荷矢量：c₁₁、c₂₁、c₄₁、c₁₂、c₂₂、c₁₃，代入(17) 式得到矩阵C

C的秩r(C)=6，载荷矢量满足线性无关的要求。

载荷历程为f₁(t)

测得对应的输出，记为y_o11(t)、y_o21(t)、y_o41(t)、y_o12(t)、y_o22(t)、y_o13(t)。

4.2   求I(t)

4.2.1   求6次敲击分别对应的I_dp(t)

利用3.3中的算法求这6次敲击对应的I_dp(t)，即I₁₁(t)、I₂₁(t)、I₄₁(t)、I₁₂(t)、I₂₂(t)、I₁₃(t)。

以I₁₁(t)为例，利用遗传算法得到前10个振动成分，如表 2所示。

表  2  前10个振动成分

Table  2  The first ten frequencies and amplitudes

f/Hz	A/με
12.29	15.41
5.77	6.87
16.92	6.03
9.39	-4.50
23.01	3.17
23.20	0.12
16.76	0.12
12.42	0.11
22.85	0.11
5.83	-0.08

下载: 导出CSV 

| 显示表格

可见，幅值递减，前5个幅值明显大于后续幅值，符合预期。前10个频率都在5.8Hz、9.4Hz、12.3Hz、16.9Hz、23.0Hz附近，这也与第2节中模态分析所得到的系统前几阶固有频率相符。进一步利用fminsearch函数优化时只需考虑这5个振动成分，按照(24)、(25) 式算得优化初值，如表 3所示。

表  3  优化初值

Table  3  Initial values of optimization

f/Hz	A/με	γ
5.77	6.80	0
9.40	-4.50	0
12.29	15.52	0
16.92	6.14	0
23.01	3.40	0

下载: 导出CSV 

| 显示表格

由于初值已经比较准确，残差r在前1000步迭代中便快速下降，经过约4000步稳定在10^-9量级，得到模态参数，如表 4所示。

表  4  优化结果

Table  4  Results of optimization

f/Hz	A/με	γ
5.77	6.89	2×10-4
9.39	-4.50	4×10-4
12.29	15.48	5×10-5
16.92	6.19	2×10-4
23.01	3.46	4×10-4

下载: 导出CSV 

| 显示表格

将表 4中的模态参数代入(20) 式，即可得到I₁₁(t)。同理，可求得I₂₁(t)、I₄₁(t)、I₁₂(t)、I₂₂(t)、I₁₃(t)。再按照(16) 中的定义即可得到

4.2.2   求I(t)

将(26)、(28) 式代入(19) 式，即可求得I(t)

由上文，此处I(t)省略了下标e和k，其中

(29) 式即为模型e₅面上天平1通道对应的UIRF，其他各面、各通道对应的UIRF求法相同。

得到其他4个面上天平1通道对应的UIRF，均与(29) 式相近，其中表征频率和阻尼比的S相同。表征幅值的系数矩阵略有不同，但其中较大的系数相对变化很小，说明在不同面上敲击激发的起主导作用的模态是相同的。由CFD计算算得5个面上的平均压力进而算出各面对应的加权系数分别为：0.3865，0.1725，0.0062，0.1150，0.3197。由(7) 式可得整个模型对天平1通道的UIRF,

测力系统在天平1通道的动态标定完成。

4.3   验证I(t)的准确性

在p₄点施加d₅方向的载荷，载荷历程为f₂(t)

测得输出应变y_o54(t)，并利用4.2中得到的I(t)计算出理论输出y₅₄(t)。将y_o54(t)、y₅₄(t)绘制在图 5中。两条曲线几乎完全重合，说明动态标定得到的I(t)是准确的。

图  5  计算应变y₅₄(t)与测量应变y_o54(t)对比

Fig.  5  Comparison between computed strain y₅₄(t) and measured strain y_o54(t)

下载: 全尺寸图片 幻灯片

4.4   验证标定方法在测得应变含噪声时的准确性

实际试验中，输出应变中含有噪声，这里验证上述动态标定方法在y_dp(t)含较大噪声时的准确性。

以p₁点处沿d₁方向的敲击为例，在y_o11(t)中加入以该信号标准差15%为幅值的随机信号，记为y_n11(t)，模拟含有较大噪声的输出信号，局部图如图 6所示。

图  6  输出信号y_o11(t)与加入噪声后的输出信号y_n11(t)

Fig.  6  Output y_o11(t) and output with noise y_n11(t)

下载: 全尺寸图片 幻灯片

利用3.3中的算法由y_n11(t)计算出I_n11(t)，对比I_n11(t)和I₁₁(t)。先利用遗传算法得到I_n11(t)的前10阶振动成分，发现与表 2中的结果十分接近。然后利用fminsearch函数继续优化，残差r稳定在约1.77，得到I_n11(t)的各振动成分，如表 5所示。

表  5  优化结果

Table  5  Results of optimization

f/Hz	A/με	γ
5.77	6.86	1×10-4
9.34	-4.50	5×10-4
12.29	15.47	4×10-5
16.92	6.17	2×10-4
23.01	3.46	4×10-4

下载: 导出CSV 

| 显示表格

表 5与表 4结果非常接近，只有低频成分的阻尼比稍有差异。这是因为本节为了验证算法，加入的噪声很大，低频阻尼比的影响远远小于噪声的影响。实际试验中噪声并不会这么大，低频阻尼比的误差会更小。绘制I₁₁(t)、I_n11(t)，局部图如图 7所示。两条曲线几乎重合，说明即使y_n11(t)中含有明显的噪声，利用3.3中的算法，得到的I_n11(t)仍然准确。

图  7  单位脉冲响应函数对比

Fig.  7  Comparison between UIRFS

下载: 全尺寸图片 幻灯片

由I_n11(t)、f₁(t)卷积计算出y′_n11(t)，绘制y′_n11(t)、y_n11(t)、y_o11(t)，局部图如图 8所示。

图  8  应变-时间曲线对比

Fig.  8  Comparison between strain responses

下载: 全尺寸图片 幻灯片

可见，y_n11(t)中的噪声基本没有被拟合，而是保留在了残差r中。

5   结论

本文给出了一种针对脉冲燃烧风洞测力系统的动态标定方法并利用ANSYS仿真对其进行了验证，得到了以下结论：

(1) 求解UIRF时通过参数化将解卷积问题转化为参数优化问题，可以有效回避问题的病态特性。

(2) 求解参数优化问题时先利用遗传算法搜索到全局最优解的近似值，再以其作为单纯形方法的初值继续优化，可以迅速求得全局最优解。

(3) 即使实际试验测得的应变信号含有较大的噪声，这种动态标定方法仍具有很高的精度。

(4) 本动态标定方法在实际应用中可能存在的误差主要由力锤敲击方向的偏差引起，试验人员在操作前需要进行适当的练习以尽量减小敲击方向的偏差。