CN104850759A - 一种风洞强迫振动动稳定性导数试验数据处理方法 - Google Patents

Abstract

本发明公开了一种强迫振动动导数试验数据数字处理方法。所述动导数数据处理方法能从低信噪比（信噪比接近0dB）的原始动态数据中分离出风洞试验仿真模拟的强迫振动力矩和试验模型角位移时间序列数据，并利用Hilbert等幅90°移相变换和三角函数系在周期内正交性特点，提取隐含在原始动态数据中的强迫振动力矩幅值、模型角位移幅值，及两者间的相位差，相位差分辨精度能达到0.1°。

Description

一种风洞强迫振动动稳定性导数试验数据处理方法

技术领域

本发明涉实验流体力学领域，具体的说是一种风洞强迫振动动稳定性导数试验数据处理方法。

背景技术

动稳定性导数(简称动导数)是飞行器气动力系数和气动力矩系数对飞行器无量纲旋转角速度或姿态角变化率的导数，是飞行器导引系统和控制系统设计以及动态品质分析的原始气动参数，主要用于描述飞行器的动态稳定性。强迫振动法是风洞试验获取飞行器动导数的主要方法之一。强迫振动法采用驱动机构带动模型在某一自由度作简谐振动，应变天平测量气动力(矩)随运动的动态响应来计算动导数。

强迫振动动导数试验可以分为俯仰、偏航和滚转三种类型的动导数试验。以俯仰动导数试验为例说明强迫振动动导数试验数据处理原理。模型作单自由度俯仰强迫振动运动时，其振动微分方程为：

$I_{\mathrm{zm}}\stackrel{\·\·}{\mathrm{\θ}}+(C-M_{\stackrel{\·}{\mathrm{\θ}}})\stackrel{\·}{\mathrm{\θ}}+(k-M_{\mathrm{\θ}})\mathrm{\θ}=M_y{e}^{i{\mathrm{\ω}}_{y}t}---\left(1\right)$

上式中：I_zm为绕穿过质心的俯仰轴的转动惯量，θ为模型的角位移，通过角位移传感器动态采集输出，C为振动系统的机械阻尼，为气动俯仰阻尼力矩导数，k振动系统的弹性铰链常数，M_θ为气动俯仰恢复力矩导数，ω_y为振动的角速度，为外加强迫振动力矩，通过应变天平采集输出，其中和M_θ为动导数试验最终需要获取的变量。

方程(1)是一个常系数、线性、二阶非齐次微分方程。方程的解由两部分组成，一部分是方程所对应的齐次形式的通解，另一部分是非齐次方程的特解。齐次方程的通解随时间很快就衰减掉了，因此关心的是其特解，即：

式中：θ₀—最大俯仰角位移，rad；—外加力矩和振动角位移之间的相位角，rad。

代入方程(1)，可得：

由公式(3)和(4)可见，强迫振动动导数试验数据处理的关键是获取施加的外加力矩振幅M_y、外加力矩与振动角位移之间的相位角但对于高速风洞动导数试验，风洞中存在的各种形式的扰动(如气流速度或压强的脉动、噪声等)所产生的气动力矩可能与被测量的气动力矩数量级相同，并且应变天平反映的总力矩信号中还混杂有大的谐波和噪声干扰信号，导致风洞动态气动力数据信噪比低，典型状态下接近0dB。另一方面受载荷变化影响，模型振动频率不是单一恒定的。故，实际采集到的角位移和力矩信号为：

$\mathrm{\θ}\left(t\right)={\mathrm{\θ}}_0\sin \left(\mathrm{\ωt}\right)+\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}{\mathrm{\θ}}_0\sin \left(\mathrm{n\ωt}\right)+{\mathrm{\θ}}_c\left(t\right)+{{\mathrm{\θ}}_0}^{\′}\sin \left[(\mathrm{\ω}+\mathrm{\Δ\ω})t\right]---\left(6\right)$

式中：

M₀—静力矩信号；

M_y θ₀sin(ωt)—基波信号；

—各次谐波信号；

M_c(t)、θ_c(t)—随机噪声信号。

θ₀′sin[(ω+△ω)t]—变频干扰信号

正是上面所述特点给动导数数字处理带来了较大困难。为了解决此问题，动导数团队对主流的三类方法进行了测试：

1)傅里叶分析法

通过直接对位移和力矩信号进行离散傅里叶变换获得其幅值和相位关系。数值仿真表明，该方法可在无噪声的情况下应用，但在位移和气动力(矩)数据分别加入1％幅值的白噪声时，傅里叶分析法相位辨识处理误差达到0.2°～1°，误差量值与有效信号量值(在典型状态下有效量值为1°)接近，而动导数试验中噪声幅值水平接近10％。

2)采用延时移相的相关法

该方法需要对位移信号移相90°，因模型名义运动频率已知，可直接对信号进行延时获得信号移相90°的数据。该方法的主要缺点是不具备抗频率漂移能力，因其使用的基本假定是模型运动和气动力(矩)响应都是单频信号，但实际动导数试验信号由于驱动电机随机失步影响，并不满足此前提，由此会造成移相误差。数值仿真也表明：模型运动频率偏移1％引起的偏差就与有效信号同量级。

3)采用差分法移相的相关法

该方法通过对信号进行差分并除以角频率的方法获得信号移相90°的数据，能克服频率偏移或漂移影响，但移相处理时会对数据信号幅值进行系数加权，且系数权值与信号自身的频率成正比，而模型强迫振动频率相对噪声频率较低，导致噪声信号被放大，降低了数据处理结果的信噪比。

发明内容

本发明的目的是在上述三种方法受均不能满足高速动导数试验数据处理需求下，通过数值仿真和理论分析，提出一种能克服噪声和振动频率漂移影响动导数数据处理方法，从存在噪声干扰和/或频率偏移的强迫振动动导数试验原始数据获取动导数。

为实现上述目的，本发明采用以下技术方案：

一种风洞强迫振动动稳定性导数试验数据处理方法，包括以下步骤：

步骤一：将采集到的角位移信号θ(t)和力矩信号M_fe(t)通过同一带通滤波器，削弱信号中的直流信号、谐波信号和背景噪声，得到滤波后的位移信号θ′(t)和力矩信号M′(t)；

步骤二：Hilbert变换是一个理想宽带相移全通网络，能使正信号滞后90°，但不改变信号幅值。结合所述特点，将所述位移信号θ′(t)在通带频率范围内进行等幅值移相90°，即将所述位移信号θ′(t)中的θ₀sin(ωt)变成θ₀′sin[(ω+△ω)t]变成在这里将移相得到的信号记为过渡信号θ₁(t)；

步骤三：将所述位移信号θ′(t)、力矩信号M_fe′(t)、过渡信号θ₁(t)进行组合得到四组新的数据序列：[θ′(t)]²、[M′(t)]²、θ′(t)M′(t)、θ₁(t)M′(t)。结合利用三角函数系在2π长周期内具有正交性的特点，以模型振动周期2π/ω为积分区间对所述四组新的数据系列进行滑动平均处理，去掉数据系列中的交流分量，将包含模型角位移幅值、力矩幅值及两者的相位差的物理量信息转移到数据序列中的直流分量内。这里将滑动平均处理的四组数据序列分别记为A、B、C、D；

步骤四：分别对四组数据序列A、B、C、D进行低通滤波，滤波上限均取为0.1Hz，以进一步去掉高频噪声信号，然后对滤波处理后的四组数据序列在整个数据长度内进行平均，就得到了以模型角位移幅值、力矩幅值、角位移与力矩信号的相位差及变频干扰信号幅值为自变量的四元二次静定方程组；

步骤五：求解四元二次静定方程组，并使用已知的模型振动角频率ω_y可得到组合量进而得到动俯仰阻尼力矩导数和气动俯仰恢复力矩导数M_θ；

在上述技术方案中，所述滤波器以振动角频率ω_y为中心选取，最小带宽(-3dB处)1Hz、谐波抑制能力大于50dB。

综上所述，由于采用了上述技术方案，本发明的有益效果是：所述动导数数据处理方法能从低信噪比(信噪比接近0dB)的原始动态数据中提取强迫振动力矩幅值、模型角位移幅值，及两者间的相位差，且相位差分辨精度能达到0.1°。

具体实施方式

本发明的实施通过以下步骤完成。

步骤一：将采集到的角位移信号θ(t)和力矩信号M_fe(t)通过同一带通滤波器，削弱信号中的直流信号、谐波信号和背景噪声，得到滤波后的位移信号θ′(t)和力矩信号M′(t)。其中：θ₀′为变频干扰噪声幅值。带通滤波器以振动角频率ω_y为中心选取，最小带宽(-3dB处)1Hz、谐波抑制能力要大于50dB。由于谐波噪声、背景噪声和直流噪声信号经滤波后对后期数字信号处理影响微弱，故这里在表达式中直接忽略

M′(t)＝M_ysin(ωt+φ_f)    (7)

θ′(t)＝θ₀sin(ωt)+θ₀′sin[(ω+△ω)t]    (8)

步骤二：Hilbert变换是一个理想宽带相移全通网络，能使正信号滞后90°，但不改变信号幅值。结合所述特点，将所述位移信号θ′(t)在通带频率范围内进行等幅值移相90°，即将所述位移信号θ′(t)中的θ₀sin(ωt)变成θ₀′sin[(ω+△ω)t]变成在这里将移相得到的信号记为过渡信号θ₁(t)：

${\mathrm{\θ}}_1\left(t\right)={\mathrm{\θ}}_0\cos \mathrm{\ωt}+{{\mathrm{\θ}}_0}^{\′}\sin [(\mathrm{\ω}+\mathrm{\Δ\ω})t+\frac{\mathrm{\π}}{2}]---\left(9\right)$

步骤三：将所述位移信号θ′(t)、力矩信号M_fe′(t)、过渡信号θ₁(t)进行组合得到四组新的数据序列：[θ′(t)]²、[M′(t)]²、θ′(t)M′(t)、θ₁(t)M′(t)。结合三角函数系在2π长周期内具有正交性的特点，以模型振动周期2π/ω为积分区间对所述四组新的数据系列进行滑动平均处理，去掉数据系列中的交流分量，将包含模型角位移幅值、力矩幅值及两者的相位差的物理量信息转移到数据序列中的直流分量内。这里将滑动平均处理的四组数据序列分别记为A、B、C、D：

$\begin{array}{c}A=\underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}{\left[{\mathrm{\θ}}^{\′}\left(t\right)\right]}^2\mathrm{dt}\\ =\underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}{({\mathrm{\θ}}_0\sin \left(\mathrm{\ωt}\right)+{{\mathrm{\θ}}_0}^{\′}\sin [(\mathrm{\ω}+\mathrm{\Δ\ω})t\left]\right)}^2\mathrm{dt}\\ =\frac{1}{2}{\mathrm{\θ}}_0^2{{\mathrm{\θ}}_0}^{\′2}\end{array}---\left(10\right)$

$\begin{array}{c}B=\underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}{\left[M^{\′}\left(t\right)\right]}^2\mathrm{dt}\\ \underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}{\left(M_y\sin \left(\mathrm{\ω}{t+\mathrm{\φ}}_f\right)\right)}^2\mathrm{dt}\\ =\frac{1}{2}{M_y}^2\end{array}$

$\begin{array}{c}C=\underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}{\mathrm{\θ}}^{\′}\left(t\right)\·M^{\′}\left(t\right)\mathrm{dt}\\ =\underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}\{{\mathrm{\θ}}_0\sin \left(\mathrm{\ωt}\right)+{{\mathrm{\θ}}_0}^{\′}\sin [(\mathrm{\ω}+\mathrm{\Δ\ω})t\left]\right\}\·M_y\sin \left(\mathrm{\ω}{t+\mathrm{\φ}}_f\right)\mathrm{dt}\\ =\frac{1}{2}{\mathrm{\θ}}_0{M}_y\cos {\mathrm{\φ}}_f\end{array}$

$\begin{array}{c}D=\underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}{\mathrm{\θ}}_1\left(t\right)\·M^{\′}\left(t\right)\mathrm{dt}\\ =\underset{-\mathrm{\π}/\mathrm{\ω}}{\overset{\mathrm{\π}/\mathrm{\ω}}{\∫}}\{{\mathrm{\θ}}_0\cos \mathrm{\ωt}+{{\mathrm{\θ}}_0}^{\′}\sin [(\mathrm{\ω}+\mathrm{\Δ\ω})t+\frac{\mathrm{\π}}{2}\left]\right\}\·M_y\sin \left(\mathrm{\ω}{t+\mathrm{\φ}}_f\right)\mathrm{dt}{\mathrm{\θ}}_0\\ =\frac{1}{2}{\mathrm{\θ}}_0{M}_y\sin {\mathrm{\φ}}_f\end{array}$

步骤四：分别对四组数据序列A、B、C、D进行低通滤波，滤波上限均取为0.1Hz，以进一步去掉高频噪声信号，然后对滤波处理后的四组数据序列在整个数据长度内进行平均，就得到了以模型角位移幅值、力矩幅值、角位移与力矩信号的相位差及变频干扰信号幅值为自变量的四元二次静定方程组：

${\left[A\right]}_{\mathrm{DC}}=({\mathrm{\θ}}_0^2+\frac{1}{2}{{\mathrm{\θ}}_0}^{\′2})/2$

[B]_DC＝M_y ²/2    (11)

[C]_DC＝θ₀M_ycosφ_f/2

[D]_DC＝θ₀M_ysinφ_f/2

步骤五：求解四元二次静定方程组，求解即可得到模型振动幅值θ₀、力矩幅值M_y及相位差φ_f，但实际上为得到气动俯仰阻尼力矩导数和气动俯仰恢复力矩导数M_θ，只需要计算组合参数和即可，这样能简化求解过程。最后，将得到的所述组合参数代入公式(3)和(4)即得到了最终所需要的所述气动俯仰阻尼力矩导数和气动俯仰恢复力矩导数M_θ；

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims (2)

1.一种风洞强迫振动动稳定性导数试验数据处理方法，其特征在于包括以下步骤：

步骤一：将风洞强迫振动动导数试验采集到的模型振动角位移和模型所受力矩动态信号序列通过同一带通滤波器，削弱两组动态序列中的直流信号、谐波信号和背景噪声，使所述两组信号从信噪比接近于0，提高到50dB以上；

步骤二：利用Hibert变换将滤波后的模型振动角位移动态信号序列在通带频率范围内进行等幅值移相90°，生成一组新的包含模型振动角位移的三角函数数据序列；

步骤三：利用三角函数集在一个周期长度内具有正交性的特点，将前两步得到的三组数据进行组合生成四组新的数据序列：

滤波处理后的模型振动角位移信号的自乘积生成数据序列A；

滤波处理后的模型力矩动态信号序列的自乘积生成数据序列B；

滤波处理后的模型振动角位移与模型力矩信号的乘积生成数据序列C；

步骤二中的包含模型振动角位移的三角函数数据序列与滤波处理后的模型力矩动态信号的乘积生成数据序列D;

对四组新的数据序列以一个单位强迫振动周期作为积分长度进行滑动平均处理，去掉信号的交流分量，将包含模型角位移幅值、力矩幅值及两者的相位差的物理量信息转移到数据序列中的直流分量内；

步骤四：对滑动平均处理的四组信号进行低通滤波，然后对整个数据序列进行平均构成四元二次静定方程组，求解该方程组就可以得到最核心的三个物理量：角位移幅值、力矩幅值、角位移与力矩信号的位移差；通过计算可以直接得到气动俯仰阻尼力矩导数和气动俯仰恢复力矩导数。

2.根据权利要求1所述的一种风洞强迫振动动稳定性导数试验数据处理方法，其特征在于所述滤波器以振动角频率w_y为中心选取，最小带宽（-3dB处）1Hz、谐波抑制能力大于50dB。