脉冲燃烧风洞的建造和运行成本较低，能够以不同尺度、方式进行基本流动和工程问题的研究，是目前国内开展大尺度超燃冲压发动机和机体推进一体化高超声速飞行器试验研究的主要设备之一[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 |
模态分析是结构动力学分析的基础，系统在载荷作用下发生振动，可视为各阶模态振动的叠加[9]。通过模态分析可近似了解系统的基本动力学特征，为开展动态标定提供重要的参考信息。利用有限元软件计算得到系统的前6阶模态，如表 1所示。
| 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 |
动态标定的目的是获得单位脉冲载荷作用下系统的响应函数(UIRF)。这个单位脉冲载荷是指由模型表面等效到天平中心点(或模型质心)后的单位脉冲载荷。试验模型本身尽管整体刚度较大，但也是一个弹性体，在模型不同位置处加载对应的UIRF不同。这里近似认为模型同一个面上各位置对应的UIRF相同，因此可以在1个面上施加一系列集中载荷，通过输入载荷和输出应变求得该面对应的UIRF。在5个面上分别进行上述操作，然后将得到的5个UIRF加权即可得到整个试验系统的UIRF。加权系数与试验状态下对应面上的平均压力大小成正比[10]。各面上的平均压力变化，加权系数要相应调整。仿真计算中，平均压力由事先的CFD计算近似得到；风洞试验时，在各面上布置测压传感器，每个时刻都可以测出1组压力值进而算得1组加权系数，实际上同一车试验中各面上平均压力的比例关系随时间变化不大，近似计算也可将不同时刻的加权系数取平均用于后续计算。
定义模型的5个表面为e1~e5，在模型的面e(e=1, 2, 3, 4, 5) 上施加一时变载荷，其相对于天平中心有6个分量，即x、y、z方向的力和这3个方向的力矩，分别记为Fx(t)，Fy(t)，Fz(t)，Mx(t)，My(t)，Mz(t)。天平输出Y有6个通道，将第k(k=1, 2, …, 6) 个通道的输出记为yk(t)。测力系统的输入输出关系为
(1)
其中，“*”表示卷积。第k行表示第k通道的输入输出关系：在面e上施加载荷，载荷相对天平中心分别为x、y、z方向单位脉冲力或力矩时，第k通道输出分别为Ixek、Iyek、Izek、Imxek、Imyek、Imzek。令
(2)
(3)
(4)
(5)
则(1) 式可写为
(6)
Ie即为面e对应的单位脉冲响应函数，是一个6×6的矩阵，整个系统的单位脉冲响应函数I为，
(7)
其中，αe为面e对应的加权系数。对于某个时间点，I是一个6×6的矩阵，动态标定的目的即是获得I。
3.2 动态标定过程本文仅针对模型的1个表面，以天平的1个通道为例介绍动态标定方法，即只求出Ie的第k行：Iek。其他各面、各通道的动态标定方法是相同的。因此只需要考虑(1) 式中的第k行
(8)
本文中仅针对e5面(e=5)，第1通道(k=1) 求解对应的I51，为了简洁，下文中略去下标e和k，将(5) 式、(8) 式分别简写为:
(9)
(10)
下文中I均指代I51，对于某个时间点为一1×6向量，下面介绍I的求法。
动态标定采用力锤敲击法，在e5面上安装加载部件(图 1中有放大显示)，为力锤敲击提供5个不同方向的加载面。加载部件的加载面面积远大于力锤锤头的面积，这样便于操作并有利于保证敲击方向的准确性。定义加载方向为d1~d5，如图 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
(11)
其中
(12)
cdp为单位载荷相对于天平中心的载荷矢量，其6个分量不独立，由加载位置和方向决定。给定d、p后，可由简单的几何关系求得cdp，是已知的。
设输出信号为ydp(t)，由(10) 式
(13)
定义
(14)
则
(15)
这里，Idp即为模型表面加载点对应的脉冲响应函数，将在3.3节中给出其确定算法，此处暂认为已知。于是Idp、cdp已知，只有I待求。
由(9) 式，I有6个未知分量，需要6个系数不相关的方程组成非齐次线性方程组才可解。因此，需要改变d、p进行6次不相关的敲击。记敲击次数为h(h=1, 2, 3, 4, 5, 6)，并令:
(16)
(17)
由(14)，(16) 和(17) 式得到
(18)
因此，只需C可逆，即cdp_1, …, cdp_6线性无关，便可解得I
(19)
设有n个时间点t1…tj…tn，那么就存在Idp(t1)…Idp(tj)…Idp(tn)共n个未知数，若直接将f和ydp(t)代入(15) 式解卷积，问题的病态特征可能使结果误差很大。
考虑结构是线弹性的，可将脉冲响应函数Idp(t)参数化
(20)
参数化将Idp写成m个不同频率振动叠加的形式。其中，m为选取的模态数量，Ai是第i个模态的幅值，ωi、γi分别是第i个模态的频率和阻尼比。参数化后，未知数变为3m个，通常，n > > m，未知数的个数大大减少。
设实际测量得到的输出信号为yodp(t)，令R(t)=yodp(t)-ydp(t)，定义残差r
(21)
求Idp(t)的问题转化为用Idp*f拟合yodp(t)的问题，即确定能使r取到最小值的m及其对应的Ai、ωi、γi这一优化问题。
3.3.2 优化算法分析3.3.1节中的这一优化问题，待定参数是m阶模态参数(ω1, A1, γ1)…(ωm, Am, γm)，优化目标是使残差r取最小值。优化过程中，r逐渐减小，给定一个值rl，以r＜rl作为终止优化的判据。
阻尼比γ一般在10-3量级或更小，在振动频率不太高、测量时长较短的情况下对优化结果影响很小，故先不考虑阻尼比，对每个模态仅考虑幅值和频率2个参数，Idp(t)简化为
(22)
求解前m未知，可预估m值，然后对3m个参数同时优化。但这对预估的要求很高，m取小了导致Idp(t)中必然缺少高阶振动成分；m取大了，则优化参数过多，优化难度较大。
这里设计了一个算法，将整个优化过程分若干次进行，每次应用一定的优化算法只辨识出1个模态对应的成分。第1次优化时，待拟合曲线为yodp(t)。第i次(i=1, 2…)优化终止时，计算Ri(t)并输出辨识出的第i阶频率ωi和幅值Ai。第i+1次优化时，以Ri(t)作为待拟合曲线，重复上一步的优化过程，得到ωi+1和Ai+1。因为每次优化的目标都是使残差尽可能减小，被识别出来的都是当次待拟合曲线中最主要的振动成分，所以各振动成分是按照幅值由大到小的顺序被辨识出来的，即Ai+1≥Ai。当进行到第M次优化后，AM < < A1，主要的模态成分已经被辨识出来，若达到了r＜rl的条件，就可以终止优化，M即为m，输出(ω1, A1)…(ωm, Am)，并由(20) 式计算Idp(t)，供后续使用。上述过程如图 4所示。
|
| 图 4 优化算法流程示意图 Fig.4 Flow diagram of the optimization algorithm |
在第i次优化中，如图 4虚线框中部分，未知数为ωi和Ai，待拟合曲线为Ri-1(t)(i=1时为yodp(t))，目标函数为ri(ωi, Ai)。考虑r(ωi, Ai)是一个多峰函数，使用传统的搜索方法极易陷入局部最优解。遗传算法是一种概率搜索算法，可以从一个种群开始并行搜索，并且不依赖于目标函数的梯度信息，具有很强的全局优化性能[11-12]，适合提取当前幅值最大的振动成分。这里采用Matlab的遗传算法工具箱(GA)进行优化。收敛准则为本代残差r与上一代残差rlg的相对差值小于ε，即
(23)
此处ε取10-4。求解前给定初代种群的取值范围，求解中如果某一代满足了收敛准则，就输出结果；如果不满足收敛准则，就通过选择性复制、交叉、变异生成下一代种群，直到满足收敛准则为止。
在上述过程中，若rl取得较小则会出现ωi≈ωj, Ai > > Aj，(i＜j)的情况。这种情况可以通过减小ε和提高yodp(t)精度而避免，然而实际中，yodp(t)是输出应变，其中必然会含有噪声，ε取得过小则会降低求解效率。实际上，这里第i次和第j次得到的就是同一个频率的振动成分，可以将对应的参数加权为1组参数作为初值继续优化
(24)
(25)
前面已经得到了(ω1, A1)…(ωm, Am)和Idp(t)。现加入对阻尼比γ的考虑，若继续基于上述的遗传算法进行优化，减小ε和rl，是可以得到预期的Idp(t)的，但是求解效率很低。进一步优化时采用Matlab优化工具箱的fminsearch函数，它是基于单纯形算法的优化工具。将前面得到的(ω1, A1)…(ωm, Am)作为频率和幅值的初值，0作为阻尼比γ的初值。由于这些初值已经很接近最优解，fminsearch函数只需在最优解附近搜索即可，这就避免了其易于陷入局部最优解的缺点，而充分发挥了其收敛速度更快的优势[15]。经过若干次迭代，残差迅速下降并达到稳定，此时终止优化，输出(ω1, A1, γ1)…(ωm, Am, γm)，并代入(20) 式计算出Idp(t)。
4 测力试验系统动态标定方法验证第2节中介绍了动态标定试验的操作方法和求系统单位脉冲响应函数的算法。本节将利用有限元软件模拟试验过程，利用上述算法求得I，并检验I的准确性。
4.1 模拟动态标定试验在Ansys Workbench中进行瞬态动力学分析，模拟上述的动态标定过程。分别在点p1(d1, d2, d4方向)、点p2(d1, d2方向)、点p3(d1方向)模拟力锤加载，各加载点位置如图 3所示。由几何关系分别算得这6次加载的载荷矢量：c11、c21、c41、c12、c22、c13，代入(17) 式得到矩阵C
(26)
C的秩r(C)=6，载荷矢量满足线性无关的要求。
载荷历程为f1(t)
(27)
测得对应的输出，记为yo11(t)、yo21(t)、yo41(t)、yo12(t)、yo22(t)、yo13(t)。
4.2 求I(t) 4.2.1 求6次敲击分别对应的Idp(t)利用3.3中的算法求这6次敲击对应的Idp(t)，即I11(t)、I21(t)、I41(t)、I12(t)、I22(t)、I13(t)。
以I11(t)为例，利用遗传算法得到前10个振动成分，如表 2所示。
| 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 |
可见，幅值递减，前5个幅值明显大于后续幅值，符合预期。前10个频率都在5.8Hz、9.4Hz、12.3Hz、16.9Hz、23.0Hz附近，这也与第2节中模态分析所得到的系统前几阶固有频率相符。进一步利用fminsearch函数优化时只需考虑这5个振动成分，按照(24)、(25) 式算得优化初值，如表 3所示。
| 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 |
由于初值已经比较准确，残差r在前1000步迭代中便快速下降，经过约4000步稳定在10-9量级，得到模态参数，如表 4所示。
| 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 |
将表 4中的模态参数代入(20) 式，即可得到I11(t)。同理，可求得I21(t)、I41(t)、I12(t)、I22(t)、I13(t)。再按照(16) 中的定义即可得到
(28)
将(26)、(28) 式代入(19) 式，即可求得I(t)
(29)
由上文，此处I(t)省略了下标e和k，其中
(30)
(29) 式即为模型e5面上天平1通道对应的UIRF，其他各面、各通道对应的UIRF求法相同。
得到其他4个面上天平1通道对应的UIRF，均与(29) 式相近，其中表征频率和阻尼比的S相同。表征幅值的系数矩阵略有不同，但其中较大的系数相对变化很小，说明在不同面上敲击激发的起主导作用的模态是相同的。由CFD计算算得5个面上的平均压力进而算出各面对应的加权系数分别为：0.3865，0.1725，0.0062，0.1150，0.3197。由(7) 式可得整个模型对天平1通道的UIRF,
(31)
测力系统在天平1通道的动态标定完成。
4.3 验证I(t)的准确性在p4点施加d5方向的载荷，载荷历程为f2(t)
(32)
测得输出应变yo54(t)，并利用4.2中得到的I(t)计算出理论输出y54(t)。将yo54(t)、y54(t)绘制在图 5中。两条曲线几乎完全重合，说明动态标定得到的I(t)是准确的。
|
| 图 5 计算应变y54(t)与测量应变yo54(t)对比 Fig.5 Comparison between computed strain y54(t) and measured strain yo54(t) |
实际试验中，输出应变中含有噪声，这里验证上述动态标定方法在ydp(t)含较大噪声时的准确性。
以p1点处沿d1方向的敲击为例，在yo11(t)中加入以该信号标准差15%为幅值的随机信号，记为yn11(t)，模拟含有较大噪声的输出信号，局部图如图 6所示。
|
| 图 6 输出信号yo11(t)与加入噪声后的输出信号yn11(t) Fig.6 Output yo11(t) and output with noise yn11(t) |
利用3.3中的算法由yn11(t)计算出In11(t)，对比In11(t)和I11(t)。先利用遗传算法得到In11(t)的前10阶振动成分，发现与表 2中的结果十分接近。然后利用fminsearch函数继续优化，残差r稳定在约1.77，得到In11(t)的各振动成分，如表 5所示。
| 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 |
表 5与表 4结果非常接近，只有低频成分的阻尼比稍有差异。这是因为本节为了验证算法，加入的噪声很大，低频阻尼比的影响远远小于噪声的影响。实际试验中噪声并不会这么大，低频阻尼比的误差会更小。绘制I11(t)、In11(t)，局部图如图 7所示。两条曲线几乎重合，说明即使yn11(t)中含有明显的噪声，利用3.3中的算法，得到的In11(t)仍然准确。
|
| 图 7 单位脉冲响应函数对比 Fig.7 Comparison between UIRFS |
由In11(t)、f1(t)卷积计算出y′n11(t)，绘制y′n11(t)、yn11(t)、yo11(t)，局部图如图 8所示。
|
| 图 8 应变-时间曲线对比 Fig.8 Comparison between strain responses |
可见，yn11(t)中的噪声基本没有被拟合，而是保留在了残差r中。
5 结论本文给出了一种针对脉冲燃烧风洞测力系统的动态标定方法并利用ANSYS仿真对其进行了验证，得到了以下结论：
(1) 求解UIRF时通过参数化将解卷积问题转化为参数优化问题，可以有效回避问题的病态特性。
(2) 求解参数优化问题时先利用遗传算法搜索到全局最优解的近似值，再以其作为单纯形方法的初值继续优化，可以迅速求得全局最优解。
(3) 即使实际试验测得的应变信号含有较大的噪声，这种动态标定方法仍具有很高的精度。
(4) 本动态标定方法在实际应用中可能存在的误差主要由力锤敲击方向的偏差引起，试验人员在操作前需要进行适当的练习以尽量减小敲击方向的偏差。
| [1] | 乐嘉陵, 刘伟雄, 贺伟, 等. 脉冲燃烧风洞及其在火箭和超燃发动机研究中的应用[J]. 实验流体力学, 2005, 19(1): 1–10. Le J L, Liu W X, He W, et al. Impulse combustion wind tunnel and its application in rocket and scramjet research[J]. Journal of Experiments in Fluid Mechanics, 2005, 19(1): 1–10. |
| [2] | 刘伟雄, 谭宇, 毛雄兵, 等. 一种新运行方式脉冲燃烧风洞研制及初步应用[J]. 实验流体力学, 2007, 21(4): 59–64. Liu W X, Tan Y, Mao X B, et al. The development and preliminary application of a pulse combustion wind tunnel with new running way[J]. Journal of Experiments in Fluid Mechanics, 2007, 21(4): 59–64. |
| [3] | Robinson M J, Mee D J, Tsai C Y, et al. Three-component force measurements on a large scramjet in a shock tunnel[J]. Journal of Spacecraft and Rockets, 2004, 41(3): 416–425. DOI:10.2514/1.10699 |
| [4] | Robinson M J, Hannemann K, Schramm J M. Design and implementation of an internal stress wave force balance in a shock tunnel[J]. CEAS Space Journal, 2011, 1(1): 45–57. |
| [5] | 贺伟, 于时恩, 李宏斌. 高超声速一体化飞行器推阻特性测量研究[J]. 实验流体力学, 2010, 24(2): 65–68. He W, Yu S E, Li H B. Experimental investigation on thrust drag performance of hypersonic integrative vehicle[J]. Journal of Experiments in Fluid Mechanics, 2010, 24(2): 65–68. |
| [6] | 贺伟, 童泽润, 李宏斌. 单模块超燃发动机推力测量天平研制[J]. 航空动力学报, 2010, 25(10): 2285–2289. He W, Tong Z R, Li H B. Investigation of thrust balance for the single module scramjet[J]. Journal of Aerospace Power, 2010, 25(10): 2285–2289. |
| [7] | Abdel-Jawad M M, Mee D J, Morgan R G. New calibration technique for multiple-component stress wave force balances[J]. Review of Scientific Instruments, 2007, 78(6): 065101–1. DOI:10.1063/1.2744235 |
| [8] | 王锋, 任虎, 周正, 等. 载荷辨识方法用于脉冲风洞模型阻力测量研究[J]. 振动与冲击, 2015, 34(24): 202–208. Wang F, Ren H, Zhou Z, et al. Drag force measurement in impulse facilities by using load identification method[J]. Journal of Vibration and Shock, 2015, 34(24): 202–208. |
| [9] | 李东旭. 高等结构动力学[M]. 北京: 科学出版社, 2010: 273-298. Li D X. Advanced structural dynamics[M]. Beijing: The Science Publishing Company, 2010: 273-298. |
| [10] | Doherty L J, Smart M K, Mee D J. Measurement of three-components of force on an airframe integrated scramjet at Mach 10[R]. AIAA-2015-3523, 2015. |
| [11] | 刘国春, 费强, 赵武云, 等. 基于Matlab遗传算法优化工具箱的应用[J]. 机械研究与应用, 2014, 27(2): 71–73. Liu G C, Fei Q, Zhao W Y, et al. Application of genetic algorithm optimization toolbox based on matlab[J]. Mechanical Research and Application, 2014, 27(2): 71–73. |
| [12] | 林鸿彬. 基于遗传算法的数据拟合在MATLAB环境中的实现[J]. 湖南农机, 2010, 37(3): 92–97. Lin H B. Data fitting based on genetic algorithm implementation in MATLAB environment[J]. Hunan Agricultural Machinery, 2010, 37(3): 92–97. |
| [13] | 罗述全. 传统优化算法与遗传算法的比较[J]. 湖北工业大学学报, 2007, 22(3): 32–35. Luo S Q. Comparison between traditional optimized algorithm and heredity algorithm[J]. Hubei University of Technology Journal, 2007, 22(3): 32–35. |
| [14] | 郭海双, 梁佳雯, 张劭昀. MATLAB遗传算法工具箱GADS优化及应用[J]. 电子设计工程, 2015, 23(10): 27–30. Guo H S, Liang J W, Zhang S Y. Optimization and examples in Matlab GA toolbox GADS[J]. Electronic Design Engineering, 2015, 23(10): 27–30. DOI:10.3969/j.issn.1674-6236.2015.10.009 |
| [15] | 杨改强, 霍丽娟, 杨国义, 等. 利用MATLAB拟合van Genuchten方程参数的研究[J]. 土壤, 2010, 42(2): 268–274. Yang G Q, Huo L J, Yang G Y, et al. Research on fitting van genuchten equation parameter with MATLAB software[J]. Soils, 2010, 42(2): 268–274. |