载荷辨识方法用于脉冲风洞模型阻力测量研究

王　锋， 任　虎， 周　正， 贺　伟

(1.中国空气动力研究与发展中心 吸气式高超声速技术研究中心，四川 绵阳　621000;2.中国空气动力研究与发展中心 高超声速冲压发动机技术重点实验室，四川 绵阳　621000)

载荷辨识方法用于脉冲风洞模型阻力测量研究

王锋1,2， 任虎1,2， 周正1， 贺伟1

(1.中国空气动力研究与发展中心 吸气式高超声速技术研究中心，四川 绵阳621000;2.中国空气动力研究与发展中心 高超声速冲压发动机技术重点实验室，四川 绵阳621000)

脉冲设备是进行高超声速气动问题研究的常用试验设备，如激波风洞，其特点是有效试验时间很短，在毫秒量级，为了在如此短的时间内对试验模型的气动力进行测量，澳大利亚昆士兰大学的Sanderson等[1-4]设计了应力波天平，基于事先获得的试验系统脉冲响应函数和试验中天平输出信号，用反卷积方法辨识试验模型在极短时间内所受气动载荷，是最早将载荷辨识用于脉冲风洞模型测力的研究工作。本文所针对的风洞是一种特殊的脉冲式设备[5-6]，其有效试验时间远大于激波风洞，可达300 ms左右。此类设备由于试验成本较低、试验段口径较大(目前最大直径为2.4 m)，是目前国内开展超燃冲压发动机和高超声速飞行器性能研究的主力设备。对于质量不太大的试验模型，通过适当提高天平的刚度，在300 ms左右的时间内，可以获得有效的输出信号，经低通滤波后，按照常规连续式风洞静态测量的处理方法得到模型的气动力[7- 8]。但随着飞行器试验模型的尺度和重量越来越大，为了在很短时间内得到有效的输出信号，必须进一步提高天平力敏元件的刚度，以保证试验系统的响应速度，这会降低天平的灵敏度和精度。未来的全尺寸试验飞行器可能重达数吨，用现在的方法可能难以获得满意的测力结果，有必要研究将载荷辨识的方法用于此类风洞的模型测力。辨识的方法在理论上不需要系统响应达到稳态，因此可以适当降低天平力敏元件的刚度，提高其灵敏度和精度。不过，针对300 ms左右的时间，应力波天平是不适用的，而仍使用常规的应变天平，前者测量的是应力杆中的应力波引起的应变信号，后者测量的是系统结构振动在天平力敏元件上引起的应变。本文以简单的钝头锥模型为对象，仅考虑单分量——阻力的测量，开展初步的方法研究，并进行试验验证，探索载荷辨识方法应用于脉冲燃烧风洞中模型测力的可行性。

1试验系统设计

试验模型是铝制的钝头锥，半锥角22.5°，底部半径125 mm，头部半径10 mm。

试验所用天平是单分量测力天平，由超硬铝整体加工而成，位于浮动框与固定框两端的两个薄板即为天平的力敏元件，其上下边缘的中部粘贴4个应变片，构成惠斯通电桥。天平结构如图1所示，(a)为俯视图，(b)为对称剖面图。

图1　天平结构Fig.1 Structure of the balance

钝头锥模型通过一个钢制的转接件与天平的浮动框连接，天平的固定框与钢制支座相连，支座下部固定于风洞内的地面。当钝头锥受到轴向载荷后，通过转接件带动天平浮动框产生位移，使力敏元件变形，通过应变片构成的电桥产生电压输出，实现测量。为了避免气流作用于转接件和天平，在模型后面用等直径的圆筒作外罩来包覆转接件和天平，外罩通过肋条固定于支座，外罩前端与钝头锥之间留有3 mm间隙，避免钝头锥在气流作用下与外罩相碰。此外，为了使运动部分的质心位于天平的几何中心，在转接件后面加了配重块。整个结构如图2所示，(a)为侧视图，(b)为外罩内部结构。

图2　试验系统组成Fig.2 Configuration of the test system

2试验系统动力学模型

载荷辨识需要知道系统动力学特性，可以由微分方程、系统脉冲响应函数或传递函数来描述。

从系统结构可以看到，相对于天平的力敏元件(两个薄板)，其它部分基本可看作刚性体。而作为力敏元件的矩形薄板在法向(面内拉伸)和侧向(面内剪切)的刚度也远高于轴向(弯曲变形)的刚度，因此，在仅受轴向载荷作用的情况下，试验系统可近似建模为沿轴向的单自由度振动系统，其振动方程为

(1)

式中:x是模型沿轴向的位移，k是天平的等效刚度，f是轴向载荷，ξ为阻尼比，ωn为无阻尼圆频率。在试验中直接测量的系统响应是天平的输出电压v，其与模型的位移存在如下对应关系

v/Ks=kx

(2)

式中:Ks为天平的轴向灵敏度，单位mV/N，于是可将式(1)变为以天平输出为变量的方程

(3)

式中：确定系统特性的参数就是ωn、ξ和Ks，下面通过动态标定来得到。

3试验系统动态标定

动态标定是通过给系统施加已知的激励，测量其动态输出，据此辨识得到系统的动力学特性参数，实现系统的动力学建模，用于后续的载荷辨识。

系统的激励通常采用阶跃载荷或脉冲载荷[3]，在此采用前者。动态标定在脉冲燃烧风洞试验现场进行，用细绳悬挂砝码给系统施加轴向载荷，待稳定后突然剪断细绳卸载，相当于给系统施加了一个反向的阶跃载荷。试验现场如图3所示。

图3　标定试验现场Fig.3 Dynamic calibration setup

细绳的断开过程需要一定的时间，假设在断裂过程中绳子的拉力是线性变化的，则卸载过程近似为一个反向的斜坡阶跃载荷。单位斜坡阶跃函数如图4所示，它本身含有一个参数，即上升时间τ，对应绳索的断开过程。

图4　单位斜坡阶跃函数Fig.4 Unit ramp-step function

式(3)描述的系统在斜率为f/τ的斜坡载荷作用下的理论输出[9]为

(4)

(5)

式(5)是加载过程的响应，对于卸载过程，相当于给加载后已处于稳态的系统施加一个反向的斜坡阶跃载荷，其响应为

v(t)=vRstp(∞)-vRstp(t)

(6)

系统参数的辨识通过使理论输出与动态标定的测量输出之间的误差平方和最小化得到，即求如下最小化问题的解

(7)

标定试验共进行了3次，所加砝码质量分别是10 kg、20 kg和40 kg，用三次试验数据得到的参数辨识结果列于表1，系统参数ωn、ξ和Ks一致性很好，而第三次试验的绳子断开时间τ明显大于前两次，这是因为加载40 kg砝码时用的绳子比前两次的粗，所以其断开时间更长。

表1　系统参数辨识结果

用辨识得到的系统参数计算得到的理论响应与实测响应的对比如图5，二者吻合很好，说明：第一，参数辨识结果是可靠的；第二，单自由度模型可以比较准确的描述系统。

图5　理论计算与实测响应对比Fig.5 Comparison between computed and measured responses

4模型阻力辨识

4.1载荷辨识方法

时域内的载荷辨识方法，通常首先将系统的动力学微分方程或卷积方程进行离散，建立待辨识时段上输入输出之间的线性关系方程[10]来求解。该方程一般是病态的，必须采取正则化处理，但在测量噪声和模型误差的影响下，获得好的解并不容易。在此，我们采用载荷分解的方法，将问题转化为一个函数参数辨识问题和一个低维的线性最小二乘问题。

根据脉冲燃烧风洞的工作特点[5-6]，在有效试验时间内，模型所受载荷总体上应比较接近一个形状如图4所示斜坡阶跃载荷，因此，用如下形式来近似模型的气动载荷

(8)

即用一个斜坡阶跃函数作为主要分量，用傅里叶级数逼近剩余分量。式中Rstp(*)代表斜坡阶跃函数，其参数包括上升时间τ、幅值f和开始时刻t0，t0即气流开始作用于模型的时刻，在此之前载荷为零，式中T为有效试验时间。之所以不直接用傅里叶级数来近似整个载荷，是因为傅里叶级数逼近阶跃型函数收敛较慢，通过用斜坡阶跃函数将主分量分离出来，可以减少所需傅里叶级数的项数，避免用到高频的三角函数，减小测量噪声的影响。后面还会看到，斜坡阶跃分量代表了试验通常所关心的载荷均值。

根据线性系统的叠加原理，在式(8)所示F(t)作用下的输出等于各分量载荷单独作用下的输出之和，将输出记作vF(t)，则有

(9)

式(9)中τ、f、t0等参数对vF(t)的影响是非线性的，而ak的影响是线性的，在此采用两步法分别来确定它们。首先令a0和ak、bk(k=1, …,N)取0，即假设模型仅受斜坡阶跃载荷，通过求解如下优化问题

(10)

(11)

vr(ti)=a0vStep(ti-t0)+

(12)

对应n个采样时刻，得n个这样的方程，将未知系数分离出来，写成矩阵形式有

ηr=Yα

(13)

(14)

最小二乘解一般不可能使方程精确满足，因此依然存在逼近误差，将其记作δ

(15)

用‖δ‖来表征误差大小，其显然与所取傅里叶级数的项数N有关，N取得越大，‖δ‖就可以越小。但是，由于测量本身存在噪声，并且系统模型本身也存在因略去高阶模态而引起的建模误差，令‖δ‖越小所得到的载荷并不一定就越接近真实载荷，需要一个准则来确定最佳的N。

(16)

式中：频率pL通过观察输出信号的频谱来确定。

4.2仿真验证

考虑系统的前两阶模态建立虚拟试验系统，构造一个虚拟的载荷输入，用计算的系统响应加上一定的噪声模拟实测输出，用来辨识输入载荷，通过与已知输入进行对比，对上述方法进行验证。

虚拟载荷历程取为带随机起伏的阶跃状曲线，如图6中实线所示，其引起的输出如图7中的实线所示，其上叠加了方差为σ2=0.52的白噪声。取截止频率pL=30 Hz，相对误差ε=0.001，载荷辨识结果如图6中的点划线所示，其中的斜坡阶跃分量如其中的虚线所示。辨识载荷产生的输出如图7中的虚线所示。可见，斜坡阶跃分量的幅值确实反映了载荷的均值，而辨识结果基本能够复现真实载荷的大幅波动，只是对阶跃处的快速变化载荷辨识精度不够高，要提高该局部的辨识精度，可以通过对式(13)求加权最小二乘解来实现，具体不在此赘述。

现在考察测量噪声对辨识结果的影响。定义辨识结果的相对误差为

(17)

表2　测量噪声的影响

现在考察系统参数误差的影响。假设用来进行载荷辨识的系统数学模型的一阶频率和阻尼比与实际系统存在误差，其定义如下

(18)

式中：x代表系统模型的频率wn或阻尼比ξ，x0代表实际值。表3给出了频率和阻尼比的独立误差对辨识结果的影响，可见频率误差对结果的影响很大，而阻尼比的误差对辨识结果影响甚微。所以，通过动态标定来比较准确地获取系统的参数十分必要。

表3　系统参数误差的影响

图6 真实载荷与辨识结果对比Fig.6Comparisonbetweenrealandestimatedloads图7 真实载荷与辨识载荷引起的输出Fig.7Outputsbyrealloadandestimatedload图8 第一车的天平输出信号Fig.8Balanceoutputofthefirsttest

4.3试验情况

吹风试验共进行了两车，由于实际配气条件不完全相同以及每车所用膜片存在差异，两次试验的来流条件亦有所不同，测得的天平输出信号分别如图8、图9所示。大约在0.75 s后，输出开始异常，这是由于风洞扩压段的气流反射回来作用于模型所致，因此有效试验段时间是0.75 s之前。

钝头锥模型的底部位于外罩之内，不直接受外流作用，但由于不是完全密封的，试验时在内外压差作用下会有气流进入，导致模型底部压力不断升高，为了比较准确的计算模型的阻力，在外罩内部靠近模型尾部的地方放置了压力传感器来测量实际压力，两车的测量结果如图10所示，取0.36s～0.75s间的数据作平均，分别是3.0 kPa和3.1 kPa，后面用于对CFD计算值作修正。

为了与辨识结果作比较，用CFD计算了钝头锥模型迎风面的阻力系数，参考面积取1 m2时，阻力系数CD为0.018 48。于是，在试验稳定段模型的平均阻力计算如下

X=qCD-PbS

(19)

式中:q是来流动压，Pb为模型底部平均压力，S是模型底部面积。即使在有效试验时段上，来流动压也是波动的，并非期望的恒定值。但因为无法直接测量动压的动态过程，只能根据配气条件、燃烧室平均总压测量值和风洞喷管参数计算给出试验段的平均动压。两车的计算结果分别是50.118 kPa和51.534 kPa。这样，式(19)给出两车试验中稳定段模型的平均阻力分别是778.92 N和800.18 N，将其视为模型实际阻力，用来考核用载荷辨识方法得到的结果。

4.4试验载荷辨识结果

取0.3 s～0.75 s时段的数据来进行载荷辨识。首先对该时段上的天平输出进行傅里叶变换，其幅频曲线如图11所示，信号主要能量分布在零频和20 Hz左右，在此取截止频率pL为30 Hz，取ε=0.001。

用4.1节的方法进行求解，两车由式(16)所确定的傅里叶级数的项数N分别是12和13，计算得到的天平输出与实测输出的对比如图12、图13所示，二者的低频特征吻合很好。

图9 第二车的天平输出信号Fig.9Balanceoutputofthesecondtest图10 模型底部压力测量曲线Fig.10Bottompressurevariation图11 天平输出的幅频曲线Fig.11Amplitude-frequencycurveofthebalanceoutputs

图12 第一车的计算与实测输出对比Fig.12Computedandmeasuredoutputsofthe1sttest图13 第二车的计算与实测输出对比Fig.13Computedandmeasuredoutputsofthe2ndtest图14 第一车的载荷辨识结果Fig.14Identifiedforcehistoryof1sttest

两车的载荷辨识结果分别如图14、图15所示，图中虚线所示为式(8)中的斜坡阶跃分量，实线为总载荷，可见，前者确实是载荷的主分量，而总载荷基本围绕斜坡阶跃分量波动，斜坡阶跃分量的稳态值代表了模型载荷的平均值，因此用它与式(19)给出的预测值作比较，来评价载荷辨识的精度，结果列于表4，二者的差别在1%左右。从总载荷曲线看，第二车比第一车更平稳一些，这是因为第二车破膜过程快，来流动压上升时间短，更接近阶跃形式。另外发现，两车载荷的起始时刻都比天平曲线上响应开始大幅上升的时刻要早，特别是第二车，在初始时刻载荷有个负170N的尖峰，这是由于在风洞启动时有个较重的快速阀门开启过程，造成较大冲击，通过风洞基础传到支架引起天平输出，图12、图13中3.5s时刻左右天平信号中的小幅波动即源于此，这可能导致初始时刻较大的载荷辨识误差。

图15　第二车的载荷辨识结果Fig.15 Identified force history of the 2nd test

车次辨识值预测值误差/%第一车770.14778.92-1.13第二车811.07800.181.36

5结论

本文以简单的钝头锥模型为对象，用载荷辨识的方法取代常规的直接测力方法，以解决短时工作的脉冲型风洞中大质量模型的测力问题。用动态标定方法获得了试验系统结构动力学参数，建立了系统的简化动力学模型。用斜坡阶跃函数和傅里叶级数之和来近似模型的气动载荷历程，将载荷辨识问题转化为函数参数辨识问题，采用分步求解方法，降低了求解难度，并提出了确定傅里叶级数项数的准则。针对试验结果，利用此方法辨识得到了模型的阻力随时间的变化曲线，与基于CFD计算和压力测量给出的平均阻力预测结果吻合良好，二者偏差在1%左右。

限于风洞现有测量手段，没有获得吹风过程中来流动压的动态变化数据，因此未能对动态载荷的辨识结果进行直接的验证。将来需要进一步完善测量手段，在模型表面分布一定数量的测压孔，用高精度的脉动压力传感器直接测量模型表面动态压力，进而获得模型的动态气动载荷，而不仅仅是平均载荷，以便对辨识结果进行更细致的精度评估。

辨识方法得到的是随时间变化的载荷历程，这对于研究超燃发动机在点火启动时的推力变化过程非常有用，而常规测力方法仅能通过低通滤波得到点火前后的稳态推阻特性[7-8]，难以较准确地反映启动的瞬态过程。因此，新方法在吸气式飞行器带动力试验研究中也有较大的应用价值。

参 考 文 献

[1] Sanderson S R, Simmons J M. Drag balance for hypervelocity impulse facilities[J]. AIAA Journal,1991, 29(12): 2185-2191.

[2] Tuttle S L, Mee D J, Simmons J M. Drag measurements at Mach 5 using a stress wave force balance[J]. Experiments in Fluids,1995, 19: 336-341.

[3] Mee D J. Dynamic calibration of force balances for impulse hypersonic facilities[J]. Shock Waves,2003(12): 443-455.

[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): 45-57.

[5] 乐嘉陵，刘伟雄，贺伟，等. 脉冲燃烧风洞及其在火箭和超燃发动机研究中的应用[J]. 实验流体力学,2005, 19(1): 1-10.

LE Jia-ling, LIU Wei-xiong, HE Wei, 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.

[6] 刘伟雄，谭宇，毛雄兵，等. 一种新运行方式脉冲燃烧风洞研制及初步应用[J]. 实验流体力学, 2007, 21(4): 59-64.

LIU Wei-xiong, TAN Yu, MAO Xiong-bing, 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.

[7] 贺伟，童泽润，李宏斌. 单模块超燃发动机推力测量天平研制[J]. 航空动力学报,2010, 25(10): 2285-2289.

HE Wei, TONG Ze-run, LI Hong-bin. Investigation of thrust balance for the single module scramjet[J]. Journal of Aerospace Power, 2010, 25(10): 2285-2289.

[8] 贺伟，于时恩，李宏斌. 高超声速一体化飞行器推阻特性测量研究[J]. 实验流体力学,2010, 24(2): 65-68.

HE Wei, YU Shi-en, LI Hong-bin. Experimental investigation on thrust drag performance of hyper-sonic integrative vehicle[J]. Journal of Experiments in Fluid Mechanics, 2010, 24(2): 65-68.

[9] Thorby D. Structural dynamics and vibration in practice: an engineering handbook[M]. Burlington: Elsevier Ltd, 2008: 53-55.

[10] 韩旭，刘杰，李伟杰，等. 时域内多源动态载荷的一种计算反求技术[J]. 力学学报,2009, 41(4): 595-602.

第一作者 王锋 男，博士，副研究员，1976年6月生

摘要：对载荷辨识方法用于脉冲燃烧风洞模型阻力测量的可行性进行了研究。设计了钝头锥模型阻力测量试验系统，通过突然卸载的方式对系统进行了动态标定，辨识得到了系统结构动力学参数，建立了系统动力学模型。基于脉冲风洞的工作特性，采用斜坡阶跃函数与傅里叶级数之和描述模型的载荷历程，将载荷辨识问题转化为函数参数辨识问题，采取先辨识斜坡阶跃函数参数再用最小二乘方法求傅里叶级数系数的求解策略。通过使理论计算输出在低频段逼近实测输出来确定傅里叶级数的项数。首先通过数值仿真检验了方法的有效性，然后在脉冲燃烧风洞进行了两车试验，基于天平输出对模型阻力进行了辨识，结果与CFD预测值吻合良好。

关键词：脉冲风洞;阻力测量;动态标定;载荷辨识

Drag force measurement in impulse facilities by using load identification method

WANGFeng1,2,RENHu1,2,ZHOUZheng1,HEWei1(1. Air-breathing Hypersonic Technology Research Center, China Aerodynamics Research and Development Center, Mianyang 621000, China;2. Science and Technology on Scramjet Laboratory, China Aerodynamics Research and Development Center, Mianyang 621000, China)

Abstract:The feasibility of load identification method used in drag measurements on a model in impulse combustion facilities was discussed. A test system was designed to measure the drag force of a blunt cone model. The system was dynamically calibrated by sudden relief of a preload to excite the system. The dynamic parameters were estimated from the balance out, and a dynamic structural model of the system was established. A ramp-step function plus Fourier series were adopted to model the load history, and the force identification was transformed to the estimation of parameters of the function. A two-step solution skill was proposed, which estimates the parameters of the ramp-step function firstly, and then solves the Fourier coefficients by least square method. The term number of the Fourier series was determined by making the theoretical response approximating to the measured response at low frequency band. The method was validated by numerical simulations at first and then by two experiments in the impulse combustion tunnel. The drag forces of the cone model were identified from the measured balance output and the results agree well with the predictions by the CFD tool.

Key words:impulse facilities; drag measurement; dynamic calibration; load identification

中图分类号:V231;O32

文献标志码：A DOI:10.13465/j.cnki.jvs.2015.24.034

收稿日期：2014-09-09修改稿收到日期：2014-11-19

基金项目：国家自然科学基金项目(11372339)；高超声速冲压发动机技术重点实验室基金项目(STSKFKT2012001)