本发明涉及民用飞机、战斗机、无人机等飞行器飞行安全地面综合试验方法，特别涉及飞行器颤振分析网格模型拉盖尔建模方法，属于航空航天与信息技术领域。
背景技术：
颤振是弹性结构在均匀气流中受到空气动力、弹性力和惯性力的耦合作用而发生的一种大幅度振动现象。对于飞机而言，在飞行中受到不确定扰动后会发生振动。此时，由于气流的作用，飞机的弹性结构如机翼、尾翼或操纵面将会产生附加气动力；作为一种激振力，附加气动力将加剧结构的振动。同时空气对飞机结构的阻尼力又试图减弱振动；在低速飞行时，由于阻尼力占优，扰动后的振动逐渐消失；当达到某一飞行速度即颤振临界速度颤振边界后，激振力占优，平衡位置失稳，将产生大幅度振动，导致飞机在数秒内解体，酿成灾难性后果；可以说，从航空工业起步的那一天起，颤振就一直是航空界研究的热门问题。
为避免颤振事故发生，新机研制必须经历颤振试验环节，以确定不发生飞行颤振的稳定飞行包线；开展颤振问题研究主要有两类途径，一是数值计算：这需要对分析对象进行数学建模，此过程需要在结构、气动等方面引入一定的假设，难以考虑真实存在的各种非线性因素和建模误差的影响，分析结果具有一定的参考价值，但可能与实际情况有较大的偏差；二是试验手段：与颤振有关的试验主要有风洞试验和飞行试验。风洞试验可以考虑气动力影响，但此方法要求将试验对象进行缩比设计，缩比模型与真实结构存在一定的差别，且由于风洞洞壁与支架的干扰气动力难免失真；此外对于高速、热环境等情况，风洞试验模拟费用昂贵且实施困难。飞行试验可以完全模拟试验对象的真实工作环境，但试验的条件受限、费用高且风险大，飞机一旦在空中发生颤振，会在几秒甚至更短的时间内解体，飞行员几乎没有处置时间，逃脱概率基本为零。
地面颤振模拟试验就是一种可以有效弥补传统试验不足的、极具生命力的颤振研究方法。地面试验以飞行器地面颤振试验系统为研究对象，以多学科设计优化理论研究为核心，密切结合飞行器地面颤振试验系统的工程特点，突破等效试验建模方法、多点分布式气动力建模与控制、颤振试验一体化检测方法等关键技术，着力解决飞行器颤振气动力模型难实现、多点激振力无法精确控制、颤振试验结果无法反复回放等问题，提高总体设计水平。
航空界、力学界虽然较早对避免颤振的问题进行了研究，但目前的研究还是初级阶段，没有形成一个系统的理论方法体系；现有的方法缺乏飞行器等价地面颤振试验方法和评价；特别是现有技术方法难以描述飞行器在不同飞行速度、大气密度、气流环境、不同温度等气动力和强度变化影响下的复杂颤振模型，使得颤振地面试验研究难以有工程化进展。
技术实现要素：
为了克服现有技术不能有效表达气动力和强度变化影响下复杂颤振模型的问题，本发明提供了一种飞行器颤振分析网格模型拉盖尔建模方法，该方法在飞行器机体轴系选择多个网格点，在不同飞行速度、大气密度、气流环境、不同温度等气动力和强度变化影响下按照机体轴系分解方法表示复杂颤振网格模型，根据建立该模型的要求提出安装传感器和数据、图像记录要求，通过有效颤振飞行试验获取数据，通过气流传感器测量值获得激励函数，采用拉盖尔函数对振动变量进行逼近和等效描述，按照辨识方法同时确定了机体轴系坐标网格点处三个轴向振动方程求解，解决了现有技术不能有效表达气动力和强度变化影响下复杂颤振模型的技术问题。
本发明解决其技术问题采用的技术方案是，一种飞行器颤振分析网格模型拉盖尔建模方法，其特征包括以下步骤：
步骤1：以飞行器机体轴系oxyz分析复杂颤振模型，在机体轴系选取n个网格点:xi,yi,zi,i＝1,2,…,n,振动时网格点动态三轴位置分量x(y,z,t),y(x,z,t),z(x,y,t)为时间t和其它两轴位置的函数，为了便于表达，以xix(y,z,t)为例，下标i＝1,2,…,n为网格点标号，下标第二个字母x,y,z分别表示振动在机体轴系oxyz的三个轴分量，为了简化问题，考虑第i＝1,2,…,n个网格点的x轴方向振动时，x(y,z,t)＝xix(t),y(x,z,t)＝yix(x),z(x,y,t)＝zix(x)，考虑第i＝1,2,…,n个网格点的y轴方向振动时，x(y,z,t)＝xiy(y),y(x,z,t)＝yiy(t),z(x,y,t)＝ziy(y)，考虑第i＝1,2,…,n个网格点的z轴方向振动时，x(y,z,t)＝xiz(z),y(x,z,t)＝yiz(z),z(x,y,t)＝ziz(t)；
在网格点邻域内建立的近似模型为：
式中，x[xix(t),yix(x),zix(x),θt,t]为在机体轴系坐标xi,yi,zi网格点的邻域内x轴向振动函数，axi[xix(t),yix(x),zix(x),θt]、bxi[xix(t),yix(x),zix(x),θt]为x轴向振动方程的结构系数函数，xix(t),yix(x),zix(x)分别为在机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n处x轴向振动时对应于xi,yi,zi的变化值；y[xiy(y),yiy(t),ziy(y),θt,t]为在机体轴系坐标xi,yi,zi网格点的邻域内y轴向振动函数，ayi[xiy(y),yiy(t),ziy(y),θt]、byi[xiy(y),yiy(t),ziy(y),θt]为y轴向振动方程的结构系数函数，xiy(y),yiy(t),ziy(y)分别为在机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n处y轴向振动时对应于xi,yi,zi的变化值；z[xiz(z),yiz(z),ziz(t),θt,t]为在机体轴系坐标xi,yi,zi网格点的邻域内z轴向振动函数，azi[xiz(z),yiz(z),ziz(t),θt]、bzi[xiz(z),yiz(z),ziz(t),θt]为z轴向振动方程的结构系数函数，xiz(z),yiz(z),ziz(t)分别为在机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n处z轴向振动时对应于xi,yi,zi的变化值；ui(xi,yi,zi,θt,t)为在xi,yi,zi网格点的等效激励函数，t为时间；θt＝[tihmafziρ]t为参数向量，ti表示xi,yi,zi网格点的温度，h为飞行高度，ma为马赫数，fzi为xi,yi,zi网格点的气流环境影响，ρ为大气密度；
步骤2：对应步骤1的机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n，安装微型温度传感器，x、y、z轴向气流、位置和振动传感器，特别是在机翼上下方和所有舵面两边安装微型x、y、z轴向气流、位置和振动传感器，同时在机身加装大于1000帧/秒的图像传感器记录观测机翼翼尖、所有舵面的振动幅值和频率；飞机机载传感器记录时间、飞行高度、马赫数，大气密度；
步骤3：将飞行器到达给定高度和马赫数后颤振试验的过程表达成有效颤振飞行试验，有效颤振飞行试验数据采样时间为tk＝0,ts,2ts,…,nts，ts为记录数据的采样周期，n+1为有效颤振飞行试验的总采样次数；通过颤振飞行试验获得离散时间tk＝0,ts,2ts,…,nts时刻机体轴系xi,yi,zi,i＝1,2,…,n网格点的xix(tk)、yiy(tk)、ziz(tk)和θt的测试值；
步骤4：根据机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n，安装微型x、y、z轴向气流传感器，特别是在机翼上下方和所有舵面两边安装微型x、y、z轴向气流传感器，确定tk＝0,ts,2ts,…,nts时刻机体轴系xi,yi,zi,i＝1,2,…,n的激励函数
对x[xix(t),yix(x),zix(x),θt,t]、y[xiy(y),yiy(t),ziy(y),θt,t]、z[xiz(z),yiz(z),ziz(t),θt,t]分别采用给定函数逼近，得到：
且
以及
步骤5：令：
以及
可将(1)式描述成：
令
式中：
px＝[px(0)px(1)…px(mx-1)px(mx)],py＝[py(0)py(1)…py(my-1)py(my)],
mx、my、mz为对应于
为m＝mx,my,mz阶拉盖尔正交多项式的递推形式，可得
式中，
令
axi(xix,θt)＝paxξx(xix),bxi(xix,θt)＝pbxξx(xix),ayi(yiy,θt)＝payξy(yiy),byi(yiy,θt)＝pbyξy(yiy),
azi(ziz,θt)＝pazξz(ziz),bzi(ziz,θt)＝pbzξz(ziz)，
式中：
pax＝[ax(0)ax(1)…ax(mx-1)ax(mx)],pbx＝[bx(0)bx(1)…bx(mx-1)bx(mx)],
pay＝[ay(0)ay(1)…ay(my-1)ay(my)],pby＝[by(0)by(1)…by(my-1)by(my)],
paz＝[az(0)az(1)…az(mz-1)az(mz)],pbz＝[bz(0)bz(1)…bz(mz-1)bz(mz)],
可得
或写成
以(3)式第一项为例，对
根据步骤3和步骤4得到的xix(tk)、yiy(tk)、ziz(tk)和
式中，
进而可得：
本发明的有益结果是：在飞行器机体轴系选择多个网格点，考虑不同飞行速度、大气密度、气流环境、不同温度等气动力和强度变化影响下按照机体轴系分解方法表示复杂颤振网格模型，根据建立该模型的要求提出安装传感器和数据、图像记录要求，通过有效颤振飞行试验获取数据，通过气流传感器测量值获得激励函数，通过气流传感器测量值获得激励函数，采用拉盖尔函数对振动变量进行逼近和等效描述，按照辨识方法同时确定了机体轴系坐标网格点处三个轴向振动方程求解，从而给出了完整的复杂颤振模型网格模型建模技术方案，解决了现有技术不能有效表达气动力和强度变化影响下复杂颤振模型的技术问题。
下面结合具体实例对本发明作详细说明。
具体实施方式
步骤1：以飞行器机体轴系oxyz分析复杂颤振模型，在机体轴系选取n个网格点:xi,yi,zi,i＝1,2,…,n,振动时网格点动态三轴位置分量x(y,z,t),y(x,z,t),z(x,y,t)为时间t和其它两轴位置的函数，为了便于表达，以xix(y,z,t)为例，下标i＝1,2,…,n为网格点标号，下标第二个字母x,y,z分别表示振动在机体轴系oxyz的三个轴分量，为了简化问题，考虑第i＝1,2,…,n个网格点的x轴方向振动时，x(y,z,t)＝xix(t),y(x,z,t)＝yix(x),z(x,y,t)＝zix(x)，考虑第i＝1,2,…,n个网格点的y轴方向振动时，x(y,z,t)＝xiy(y),y(x,z,t)＝yiy(t),z(x,y,t)＝ziy(y)，考虑第i＝1,2,…,n个网格点的z轴方向振动时，x(y,z,t)＝xiz(z),y(x,z,t)＝yiz(z),z(x,y,t)＝ziz(t)；
在网格点邻域内建立的近似模型为：
式中，x[xix(t),yix(x),zix(x),θt,t]为在机体轴系坐标xi,yi,zi网格点的邻域内x轴向振动函数，axi[xix(t),yix(x),zix(x),θt]、bxi[xix(t),yix(x),zix(x),θt]为x轴向振动方程的结构系数函数，xix(t),yix(x),zix(x)分别为在机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n处x轴向振动时对应于xi,yi,zi的变化值；y[xiy(y),yiy(t),ziy(y),θt,t]为在机体轴系坐标xi,yi,zi网格点的邻域内y轴向振动函数，ayi[xiy(y),yiy(t),ziy(y),θt]、byi[xiy(y),yiy(t),ziy(y),θt]为y轴向振动方程的结构系数函数，xiy(y),yiy(t),ziy(y)分别为在机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n处y轴向振动时对应于xi,yi,zi的变化值；z[xiz(z),yiz(z),ziz(t),θt,t]为在机体轴系坐标xi,yi,zi网格点的邻域内z轴向振动函数，azi[xiz(z),yiz(z),ziz(t),θt]、bzi[xiz(z),yiz(z),ziz(t),θt]为z轴向振动方程的结构系数函数，xiz(z),yiz(z),ziz(t)分别为在机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n处z轴向振动时对应于xi,yi,zi的变化值；ui(xi,yi,zi,θt,t)为在xi,yi,zi网格点的等效激励函数，t为时间；θt＝[tihmafziρ]t为参数向量，ti表示xi,yi,zi网格点的温度，h为飞行高度，ma为马赫数，fzi为xi,yi,zi网格点的气流环境影响，ρ为大气密度；
步骤2：对应步骤1的机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n，安装微型温度传感器，x、y、z轴向气流、位置和振动传感器，特别是在机翼上下方和所有舵面两边安装微型x、y、z轴向气流、位置和振动传感器，同时在机身加装大于1000帧/秒的图像传感器记录观测机翼翼尖、所有舵面的振动幅值和频率；飞机机载传感器记录时间、飞行高度、马赫数，大气密度；
步骤3：将飞行器到达给定高度和马赫数后颤振试验的过程表达成有效颤振飞行试验，有效颤振飞行试验数据采样时间为tk＝0,ts,2ts,…,nts，ts为记录数据的采样周期，n+1为有效颤振飞行试验的总采样次数；通过颤振飞行试验获得离散时间tk＝0,ts,2ts,…,nts时刻机体轴系xi,yi,zi,i＝1,2,…,n网格点的xix(tk)、yiy(tk)、ziz(tk)和θt的测试值；
步骤4：根据机体轴系坐标网格点xi,yi,zi,i＝1,2,…,n，安装微型x、y、z轴向气流传感器，特别是在机翼上下方和所有舵面两边安装微型x、y、z轴向气流传感器，确定tk＝0,ts,2ts,…,nts时刻机体轴系xi,yi,zi,i＝1,2,…,n的激励函数
对x[xix(t),yix(x),zix(x),θt,t]、y[xiy(y),yiy(t),ziy(y),θt,t]、z[xiz(z),yiz(z),ziz(t),θt,t]分别采用给定函数逼近，得到：
且
以及
步骤5：令：
以及
可将(1)式描述成：
令
式中：
px＝[px(0)px(1)…px(mx-1)px(mx)],py＝[py(0)py(1)…py(my-1)py(my)],
mx、my、mz为对应于
为m＝mx,my,mz阶拉盖尔正交多项式的递推形式，可得
式中，
令
axi(xix,θt)＝paxξx(xix),bxi(xix,θt)＝pbxξx(xix),ayi(yiy,θt)＝payξy(yiy),byi(yiy,θt)＝pbyξy(yiy),
azi(ziz,θt)＝pazξz(ziz),bzi(ziz,θt)＝pbzξz(ziz)，
式中：
pax＝[ax(0)ax(1)…ax(mx-1)ax(mx)],pbx＝[bx(0)bx(1)…bx(mx-1)bx(mx)],
pay＝[ay(0)ay(1)…ay(my-1)ay(my)],pby＝[by(0)by(1)…by(my-1)by(my)],
paz＝[az(0)az(1)…az(mz-1)az(mz)],pbz＝[bz(0)bz(1)…bz(mz-1)bz(mz)],
可得
或写成
以(3)式第一项为例，对
根据步骤3和步骤4得到的xix(tk)、yiy(tk)、ziz(tk)和
式中，
进而可得：