CN104133933B - 一种高超声速飞行器热环境下气动弹性力学特性分析方法 - Google Patents

Abstract

本发明涉及一种气动弹性分析方法，提供了一种高超声速飞行器热环境下气动弹性力学特性分析方法。使用活塞理论对高超声速飞行器全机有限元模型进行频域非定常气动力计算，在此基础上考虑了在高超声速热环境中，模型受到的气动热效应，忽略气动热输入对气动力输入和弹性力输入的弱耦合效应，仅考虑气动热载荷作用下的结构温度的稳态特性，在进行定常气动力求解之后采用参考焓法求得飞行器表面的热流密度，进而计算出飞行器表面的稳态温度分布，并求得此时结构的实际等效刚度矩阵，并采用工程方法求解出临界颤振速度。本发明解决了气动热环境下超声速飞行器气动弹性分析问题，通过对颤振速度的分析从而提高高超声速飞行器气动弹性性能。

Description

一种高超声速飞行器热环境下气动弹性力学特性分析方法

技术领域

本发明涉及一种航空航天领域的气动弹性分析方法，具体地说，是高超声速飞行器考虑热效应情况下的气动弹性分析方法。

背景技术

在气动弹性学科的研究初期，通常采用的方法是以结构和流体动力学模型为基础进行线性化假设，颤振边界问题归根到底是非自伴系统模型的特征值求解。非自伴系统模型基于颤振运动微分方程而构建，对于位于颤振临界点附近的简谐振动特性分析也能获得可靠的结果。该方法的不足之处是，需要对模型进行线性化假设，系统中若带有非线性的气动弹性问题，就显无能为力了。伴随着飞行器理论的发展，结构上采用高升阻比轮廓、机身采用智能材料、运用智能变形技术，气流三维效应、结构间的相互干扰等非线性气动力因素，使得材质、负荷、装配缝隙和外形柔性形变等结构非线性效应变得更加突出。

高超声速气动热弹性问题是典型的多场耦合问题，通过对此一次求解而得出结果在工程实践中可望而不可及。在气动热环境下，原本由气动力、惯性力和弹性力三者相互作用的耦合关系，又会加入气动热耦合，使这四者相互影响形成不同程度的耦合作用关系。高超声速气动热弹性力学问题求解的关键之一就是获得适合于气动弹性分析的非定常气动力。非线性的引入，使得非定常气动力求解方法的使用受到极大限制，比如：基于势流理论和模态叠加法只能局限于时域求解微分方程，获得系统响应曲线，从而对系统的稳定性进行判别。时域分析的好处是可以让模型脱离结构/气动非线性的限制，进行多场耦合状态下的时域仿真，并且能够结合工作模态分析中的最小二乘法曲线拟合、小波变换法、分段滑动、自回归滑动平均法等具体手段，从而确定系统频率、阻尼系数等气动弹性的重要临界参数。但是计算过程复杂、工作量大，尤其在时域方法中做CFD/CSD的直接耦合，对计算机的配置提出了很高的要求，限制了其在工程实践中的使用范围。

从国内外研究资料现状可以发现，对大高超声速飞行器进行热气动弹性分析方面的研究工作还有待发掘，利用现有商用有限元软件，结合工程计算方法，对具有真实翼型的高超声速飞行器全机模型进行热气动弹性分析对于高超声速飞行器的发展有着很强的现实意义，怎样使用热颤振分析理念，对高超声速飞行器进行热气动弹性计算分析已成为了一个新的命题。

发明内容

本发明针对上述现有技术状况而设计提供了一种高超声速飞行器热环境下气动弹性力学特性分析方法，其目的是解决在气动热环境下，对高超声速飞行器进行气动弹性分析，判断高超声速飞行器气动弹性性能满足技术要求与否，进而提高高超声速飞行器气动弹性性能。

本发明所采用的技术方案包括以下步骤：

(1)建立高超声速飞行器全机几何模型和有限元模型，包括以下内容：

(a)在三维造型软件中建立高超声速飞行器全机几何模型，几何模型采用有利于机身和超燃冲压发动机实现一体化布局设计的二维升力体外形，同时机身和全动平尾采用了双楔形薄翼型，垂尾采用的是梯形翼型，其中最易发生颤振的平尾部分与机身由转轴相连；

(b)在将几何模型读入MSC.Patran之后，建立有限元模型并对模型进行有限元网格划分，根据全机的组成部件分别处理，机身蒙皮、垂尾表面蒙皮和平尾表面蒙皮等结构部分选择壳单元建模，并采用MSC.Patran中的四边形单元进行自动划分，机身内部和平尾内部采用梁结构，并赋予梁单元属性，飞机头部、垂尾前缘和平尾前缘部分采用实体单元建模，使用MSC.Patran中Match Parasoild Faces/Neighbor Solid List命令的智能控制功能可以很好的保障各部分之间单元的连续性并使用MSC.Patran中的四面体单元进行自动划分；

(c)在材料的选用上，机身主要选用了高强度、高密度、耐高温的金属材料结构，外层用氧化铝隔热层覆盖以保护其内部结构，其中在机头和平尾前部都采用了C-C复合材料以耐受高温，机头C-C材料后方是Densalloy 180钨合金，同时对刚度和密度进行了适当折减，在逐一对各个部件划分有限单元并赋予材料属性的基础上，对各部件加入集中质量，使各个部件的质量特性更加合理，系统及燃油质量等，用集中质量卡(CONM2)施加于相应质心上，并用MPC元约束在主盒段上，全动平尾部分采用的是Haynes 230镍合金蒙皮和梁结构骨架，用BEAM单元建立刚轴大梁来模拟水平安定面和升降舵的刚度特性，CONM2单元模拟结构质量，水平安定面与升降舵之间采用约束X、Y、Z三个方向自由度的MPC来模拟铰链连接，采用BAR、ROD单元来模拟作动器和摇臂；

(2)针对高超声速飞行器全机有限元模型，求解其在初始状态下振动模态：根据模型的对称性，取半机身，根部固支，通过MSC.Nastran求解序列SOL103对模型在初始状态下进行振动模态分析，再用MSC.Patran进行后处理分析得到的前n阶模态；

(3)计算模型表面受到的热流密度场分布：

机翼前缘部分受到的气动加热现象最为明显，高超声速飞行器薄型机翼机身的前缘部分和翼面部分可以忽略厚度因素看做平板来进行计算：

参考粘性系数μ^*通过萨特兰表达式求解得出：

参考密度ρ^*通过状态方程表达式求解得出：

同时可以得到参考雷诺数，如下：

斯坦顿数通过雷诺比拟式求解得出：

最终高超声速飞行器表面热流通过以下公式求得：

Q_aero＝St^*ρ^*V_∞c_ρ(T_r-T_ω)

沿着机身x方向将机身分成若干等份，假设每份所处的热流密度一定，则热流场为若干个不相互连续的离散温度值构成，在初始环境温度一定时，计算得到的热流场分布；

(4)将热流密度场加载到有限元模型上，计算其有限元模型的稳态温度场分布：考虑流经飞行器表面的气流，飞行器表面任一点的τ_ω等价与相应的可压流的τ_ω，不可压缩气流的温度升高至一个给定的参考温度值T^*：

T^*＝0.5T_ω+0.22T_r+0.28T_∞

其中：

T^*表示飞行器表面温度，

表示边界层外缘温度，由气动力计算部分参数可以得到，

T_r为恢复温度，

根据算出的热流密度场，在MSC.Patran中将其加载在有限元模型表面，同时加上适当的辐射值加以平衡，通过MSC.Nastran求解序列SOL153对其进行温度分布计算，得到其有限元模型的稳态温度场分布；

(5)根据稳态温度场分布情况得到稳态热应力变形后的结构参数：

令结构初始单元线性刚度矩阵为K₀，为热应力作用生成的多余刚度为K_σ，则结构受热效应作用后的总的应力刚度矩阵等效为

K＝K₀+K_σ

则在热效应作用下的结构振动表达式为

式中，M表示质量矩阵，

表示振型，

ω为振频，

此时的动力学方程即描述了为结构在热载荷作用下结构振动特性，即简化为求上式的广义特征值问题，根据得到的温度分布情况，进一步的通过MSC.Nastran求解序列SOL153对对结构进行静态分析，得到稳态热应力变形后的结构参数，包括刚度矩阵，此时的刚度矩阵即为考虑热效应下的等效刚度矩阵；

(6)建立气动网格模型：应用MSC Flightloads中的气动弹性模块，将模型划分成数目适合的气动分区，每个气动分区又会分成数目适合的气动片，获得气动网格参数；

(7)根据结构的刚度矩阵等结构参数和划分的气动网格参数，使用活塞理论进行线性的频域非定常气动力计算：

初始状态下的气体压力大小为p_∞，密度为ρ_∞，声速大小为a_∞，假设在等熵条件下

其中，γ为气体比热比。当地流音速大小为

两边同时积分可以得到

无穷远处的p＝p_∞，且此时v_∞＝0，得到

得扰动压力为

三阶扰动压力表达式：

假设x为顺气流方向，则

当结构模型为上下对称时，上式则被简化为

表示上表面的下洗

表示下表面的下洗

令

得到气体压差的表达式为

从MSC.Patran中导出结构的刚度矩阵等结构参数和划分的气动网格参数，并编程实现活塞理论计算其非定常气动力；

(8)根据(5)得到的稳态热应力变形后的结构参数和(7)得到的频域非定常气动力，使用工程计算方法进行颤振计算，得到颤振速度：

根据有限元法建立结构动力学方程为

式中：x表示节点自由度，

F_N表示气动力对应的等效节点力。

把结构动力学方程变换到模态坐标系下，其表示如下：

式中，F表示等效节点气动力对应的等效模态气动力，节点自由度x与模态自由度q相互转化关系方程表达式为

x＝Φ^Tq＝∑φ_iq_i

式中

Φ＝[φ₁，φ₂，φ₃...]

通过MSC.Nastran计算得到模型相关的振型参数，根据活塞理论气动力计算方法求解得到模态坐标系下的气动力，其中第i模态时，模态力的大小可以表示为

式中：φ_i(x，y)表示的是第i阶连续模态，是通过φ_i插值计算得到的，

F_A(x，y)表示的是连续气动力，

φs_i(x，y)表示的是第j个气动有限单元上的第i阶连续模态，其中略去去了下标j以简化表达式，

Fs_A(x，y)表示的是第j个气动有限单元上的连续气动力，其中略去去了下标j以简化表达式，

ws(x，y，t)表示的是第j个气动有限单元上的下洗速度。由定义我们可以得到；

序号为j气动有限单元有四个节点构成，他们按照(Nd1，Nd2，Nd3，Nd4)的顺序逆时针分布，这四个节点的节点坐标可以表示为(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，(x₄，y₄)，第i阶模态对应的节点在z方向的位移可以用φ_i来表示，令其取值分别为(此处采用线性插值，因而忽略斜率改变)，选择的插值函数为：

将x，y投影到r，s坐标下的变换矩阵表达式为：

为简便计算对公式进行进一步化简，令

在任一气动有限单元中为常数，气动力可以进一步简化，假设

其中气动力系数在任一气动有限单元片中为定值，无需进行积分运算，对于模态载荷

对式中的部分进行运算，

过程如下：

因为划分的有限单元数量庞大，我们选用简单高效的高斯积分法，同时符合3阶精度要求，选用四点进行积分，其中：

r₁＝0.57735，r₂＝-0.57735，s₁＝0.57735，s₂＝-0.57735，

通过MSC.Nastran计算得到的函数f′(x，y)是一个与厚度有关的量，因为选用的气动面单元是平面，则气动面单元的函数表达式为

f(x，y)＝a+bx+cy

式中参变量系数可以通过拟合的方法求出，其中

当气动有限单元片为三角形单元时，通过假设第三个节点域第四个节点重合来处理，进行颤振计算后得到颤振速度以及对应的颤振频率。

本发明的有益效果是：

1)研究了高超声速飞行器热气动弹性分析的理论基础，探讨了包括高超声速非定常气动力的工程计算方法等常用的研究方法，为进一步探索提供了理论上的指导。

2)建立了高超声速飞行器几何模型和有限元模型，作为后续进行热颤振计算和气动弹性分析的研究对象。

3)实现了一种求解高超声速飞行器热气动弹性分析过程中气动力、气动热和颤振速度的计算方法，并以有限元模型为研究对象，使用该方法可以得到其考虑热效应情况下的颤振计算结果。

通过以下的描述并结合附图，本发明将变得更加清晰，这些附图用于解释本发明的实施例。

附图说明

图1是本发明一种高超声速飞行器热环境下气动弹性力学特性分析方法的基本流程的示意图；

图2和图3分别是本发明步骤(1)建立的某型高超声速飞行器全机几何模型和有限元模型图；

图4是本发明步骤(2)求解某型高超声速飞行器在初始状态下的前六阶振动模态的振型图；

图5是本发明步骤(3)对某型高超声速飞行器全机进行气动热计算，得到热流场分布图；

图6是本发明步骤(4)将热流密度场加载到某型高超声速飞行器全机有限元模型上，计算得到的稳态温度场分布图；

图7和图8分别是本发明步骤(6)某型高超声速飞行器全机模型的气动分区及气动网格划分；

图9和图10分别是本发明步骤(8)为计算得到的某型高超声速飞行器全机颤振V-g和V-f曲线图。

具体实施方式

方法实施例：以某高超声速飞行器全机气动弹性一体化的分析为例，并针对一种高超声速飞行器热环境下气动弹性力学特性分析方法的基本流程的示意图如图1所示，说明利用一种高超声速飞行器热环境下气动弹性力学特性分析方法分析高超声速飞行器气动弹性性能的具体实施方法。

(1)建立高超声速飞行器全机几何模型和有限元模型，包括以下内容：

(a)基于某型高超声速飞行器实际尺寸轮廓与布局形式在三维造型软件CATIA中建立全机几何模型，几何模型采用有利于机身和超燃冲压发动机实现一体化布局设计的二维升力体外形，同时机身和全动平尾采用了双楔形薄翼型，垂尾采用的是梯形翼型，其中最易发生颤振的平尾部分与机身由转轴相连；

(b)在将某型高超声速飞行器几何模型读入MSC.Patran之后，建立有限元模型所示并对模型进行有限元网格划分，根据全机的组成部件分别处理，机身蒙皮、垂尾表面蒙皮和平尾表面蒙皮等结构部分选择壳单元建模，并采用MSC.Patran中的四边形单元进行自动划分，机身内部和平尾内部采用梁结构，并赋予梁单元属性，飞机头部、垂尾前缘和平尾前缘部分采用实体单元建模，使用MSC.Patran中Match Parasoild Faces/Neighbor SolidList命令的智能控制功能可以很好的保障各部分之间单元的连续性并使用MSC.Patran中的四面体单元进行自动划分；

(c)在材料的选用上，机身主要选用了高强度、高密度、耐高温的金属材料结构，外层用氧化铝隔热层覆盖以保护其内部结构，其中在机头和平尾前部都采用了C-C复合材料以耐受高温，机头C-C材料后方是Densalloy180钨合金，同时对刚度和密度进行了适当折减，在逐一对各个部件划分有限单元并赋予材料属性的基础上，对各部件加入集中质量，使各个部件的质量特性更加合理，系统及燃油质量等，用集中质量卡(CONM2)施加于相应质心上，并用MPC元约束在主盒段上，全动平尾部分采用的是Haynes230镍合金蒙皮和梁结构骨架，用BEAM单元建立刚轴大梁来模拟水平安定面和升降舵的刚度特性，CONM2单元模拟结构质量，水平安定面与升降舵之间采用约束X、Y、Z三个方向自由度的MPC来模拟铰链连接，采用BAR、ROD单元来模拟作动器和摇臂，建立的某型高超声速飞行器全机几何模型和有限元模型图分别如图2和图3所示；

(2)针对高超声速飞行器全机有限元模型，求解其在初始状态下振动模态：根据模型的对称性，取半机身，根部固支，通过MSC.Nastran求解序列SOL103对模型在初始状态下进行振动模态分析，再用MSC.Patran进行后处理分析得到的前n阶模态，对应的前六阶模态的频率为：

模态	第一阶	第二阶	第三阶	第四阶	第五阶	第六阶
							频率ω(Hz)	37.596	77.909	84.968	93.229	118.18	178.68

对应的前六阶模态的振型图如图4所示；

(3)计算模型表面受到的热流密度场分布：

机翼前缘部分受到的气动加热现象最为明显，高超声速飞行器薄型机翼机身的前缘部分和翼面部分可以忽略厚度因素看做平板来进行计算：

参考粘性系数μ^*通过萨特兰表达式求解得出：

参考密度ρ^*通过状态方程表达式求解得出：

同时可以得到参考雷诺数，如下：

斯坦顿数通过雷诺比拟式求解得出：

最终高超声速飞行器表面热流通过以下公式求得：

Q_aero＝St^*ρ^*V_∞c_ρ(T_r-T_ω)

沿着机身x方向将机身分成若干等份，假设每份所处的热流密度一定，则热流场为若干个不相互连续的离散温度值构成，在初始环境温度一定时，计算得到的热流场。本例沿着机身x方向将机身分成20等份，假设每份所处的热流密度一定，则热流场为20个不相互连续的离散温度值构成，在初始环境温度为300K时，计算得到的热流场分布如图5所示，从该图可以看出热流密度沿机身x向从183935J/(m²·s)到140000J/(m²·s)逐渐降低，飞机头部为所受热流量最大的地方；

(4)将热流密度场加载到有限元模型上，计算其有限元模型的稳态温度场分布：考虑流经飞行器表面的气流，飞行器表面任一点的τ_ω等价与相应的可压流的τ_ω，不可压缩气流的温度升高至一个给定的参考温度值T^*：

T^*＝0.5T_ω+0.22T_r+0.28T_∞

其中：

T^*表示飞行器表面温度，

表示边界层外缘温度，由气动力计算部分参数可以得到，

T_r为恢复温度，

根据算出的热流密度场，在MSC.Patran中将其加载在有限元模型表面，同时加上适当的辐射值加以平衡，通过MSC.Nastran求解序列SOL153对其进行温度分布计算，得到其有限元模型的稳态温度场分布如图6所示，机身前舱段温度、垂尾与平尾温度约为1600℃，机身温度相对低一些，约为819℃；

(5)根据稳态温度场分布情况得到稳态热应力变形后的结构参数：

令结构初始单元线性刚度矩阵为K₀，为热应力作用生成的多余刚度为K_σ，则结构受热效应作用后的总的应力刚度矩阵等效为

K＝K₀+K_σ

则在热效应作用下的结构振动表达式为

式中，M表示质量矩阵，

表示振型，

ω为振频，

此时的动力学方程即描述了为结构在热载荷作用下结构振动特性，即简化为求上式的广义特征值问题，根据得到的温度分布情况，进一步的通过MSC.Nastran求解序列SOL153对对结构进行静态分析，得到稳态热应力变形后的结构参数，包括刚度矩阵，此时的刚度矩阵即为考虑热效应下的等效刚度矩阵，此时的前六阶固有振动频率与之前为考虑热因素情况下的对比表为：

随着温度场平均温度的升高，各阶刚度呈下降趋势；

(6)建立气动网格模型：应用MSC Flightloads中的气动弹性模块，将模型划分成数目适合的气动分区，每个气动分区又会分成数目适合的气动片，获得气动网格参数，本例中沿着展向划分为8个气动分区，每个气动分区又会分成6个的气动片，而平尾和垂尾作为单独的气动分区，沿展向分别划分为6条气动片。模型的气动分区及气动网格划分分别如图7和图8所示；

(7)根据结构的刚度矩阵等结构参数和划分的气动网格参数，使用活塞理论进行线性的频域非定常气动力计算：

初始状态下的气体压力大小为p_∞，密度为ρ_∞，声速大小为a_∞，假设在等熵条件下

其中，γ为气体比热比。当地流音速大小为

两边同时积分可以得到

无穷远处的p＝p_∞，且此时v_∞＝0，得到

得扰动压力为

三阶扰动压力表达式：

假设x为顺气流方向，则

当结构模型为上下对称时，上式则被简化为

表示上表面的下洗

表示下表面的下洗

令

得到气体压差的表达式为

从MSC.Patran中导出结构的刚度矩阵等结构参数和划分的气动网格参数，并编程实现活塞理论计算其非定常气动力；

(8)根据(5)得到的稳态热应力变形后的结构参数和(7)得到的频域非定常气动力，使用工程计算方法进行颤振计算，得到颤振速度：

根据有限元法建立结构动力学方程为

式中：x表示节点自由度，

F_N表示气动力对应的等效节点力。

把结构动力学方程变换到模态坐标系下，其表示如下：

式中，F表示等效节点气动力对应的等效模态气动力，节点自由度x与模态自由度q相互转化关系方程表达式为

x＝Φ^Tq＝∑φ_iq_i

式中

Φ＝[φ₁，φ₂，φ₃...]

通过MSC.Nastran计算得到模型相关的振型参数，根据活塞理论气动力计算方法求解得到模态坐标系下的气动力，其中第i模态时，模态力的大小可以表示为

式中：φ_i(x，y)表示的是第i阶连续模态，是通过φ_i插值计算得到的，

F_A(x，y)表示的是连续气动力，

φs_i(x，y)表示的是第j个气动有限单元上的第i阶连续模态，其中略去去了下标j以简化表达式，

Fs_A(x，y)表示的是第j个气动有限单元上的连续气动力，其中略去去了下标j以简化表达式，

ws(x，y，t)表示的是第j个气动有限单元上的下洗速度。由定义我们可以得到：

序号为j气动有限单元有四个节点构成，他们按照(Nd1，Nd2，Nd3，Nd4)的顺序逆时针分布，这四个节点的节点坐标可以表示为(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，(x₄，y₄)，第i阶模态对应的节点在z方向的位移可以用φ_i来表示，令其取值分别为(此处采用线性插值，因而忽略斜率改变)，选择的插值函数为：

将x，y投影到r，s坐标下的变换矩阵表达式为：

为简便计算对公式进行进一步化简，令

在任一气动有限单元中为常数，气动力可以进一步简化，假设

其中气动力系数在任一气动有限单元片中为定值，无需进行积分运算，对于模态载荷

对式中的部分进行运算，

过程如下：

因为划分的有限单元数量庞大，我们选用简单高效的高斯积分法，同时符合3阶精度要求，选用四点进行积分，其中：

r₁＝0.57735，r₂＝-0.57735，s₁＝0.57735，s₂＝-0.57735，

通过MSC.Nastran计算得到的函数f′(x，y)是一个与厚度有关的量，因为选用的气动面单元是平面，则气动面单元的函数表达式为

f(x，y)＝a+bx+cy

式中参变量系数可以通过拟合的方法求出，其中

当气动有限单元片为三角形单元时，通过假设第三个节点域第四个节点重合来处理，进行颤振计算后得到颤振马赫数为4.16，其对应的颤振频率287.5Hz，其耦合形式未发生改变，仍是平尾部分翼面扭转振型和挥舞振型的耦合。其颤振V-g和V-f曲线图分别如图9和图10所示。

Claims (1)

1.一种高超声速飞行器热环境下气动弹性力学特性分析方法，其特征在于，包括以下步骤：

(1)建立高超声速飞行器全机几何模型和有限元模型，包括以下内容：

(a)在三维造型软件中建立高超声速飞行器全机几何模型，几何模型采用有利于机身和超燃冲压发动机实现一体化布局设计的二维升力体外形，同时机身和全动平尾采用了双楔形薄翼型，垂尾采用的是梯形翼型，其中最易发生颤振的平尾部分与机身由转轴相连；

(b)在将几何模型读入MSC.Patran之后，建立有限元模型并对模型进行有限元网格划分，根据全机的组成部件分别处理，机身蒙皮、垂尾表面蒙皮和平尾表面蒙皮结构部分选择壳单元建模，并采用MSC.Patran中的四边形单元进行自动划分，机身内部和平尾内部采用梁结构，并赋予梁单元属性，飞机头部、垂尾前缘和平尾前缘部分采用实体单元建模，使用MSC.Patran中Match Parasoild Faces/Neighbor Solid List命令的智能控制功能可以很好的保障各部分之间单元的连续性并使用MSC.Patran中的四面体单元进行自动划分；

(c)在材料的选用上，机身主要选用了高强度、高密度、耐高温的金属材料结构，外层用氧化铝隔热层覆盖以保护其内部结构，其中在机头和平尾前部都采用了C-C复合材料以耐受高温，机头C-C材料后方是Densalloy180钨合金，同时对刚度和密度进行了适当折减，在逐一对各个部件划分有限单元并赋予材料属性的基础上，对各部件加入集中质量，使各个部件的质量特性更加合理，系统及燃油质量，用集中质量卡CONM2施加于相应质心上，并用MPC元约束在主盒段上，全动平尾部分采用的是Haynes230镍合金蒙皮和梁结构骨架，用BEAM单元建立刚轴大梁来模拟水平安定面和升降舵的刚度特性，CONM2单元模拟结构质量，水平安定面与升降舵之间采用约束X、Y、Z三个方向自由度的MPC来模拟铰链连接，采用BAR、ROD单元来模拟作动器和摇臂；

(2)针对高超声速飞行器全机有限元模型，求解其在初始状态下振动模态：根据模型的对称性，取半机身，根部固支，通过MSC.Nastran求解序列SOL103对模型在初始状态下进行振动模态分析，再用MSC.Patran进行后处理分析得到的前n阶模态；

(3)计算模型表面受到的热流密度场分布：

机翼前缘部分受到的气动加热现象最为明显，高超声速飞行器薄型机翼机身的前缘部分和翼面部分可以忽略厚度因素看做平板来进行计算：

参考粘性系数μ^*通过萨特兰表达式求解得出：

${\mathrm{\μ}}^{*}={\left(\frac{T^{*}}{288.15}\right)}^{1.5}\frac{398.55}{T^{*}+110.4}\×1.7894\×{10}^{-5}$

参考密度ρ^*通过状态方程表达式求解得出：

${\mathrm{\ρ}}^{*}=\frac{p_{\mathrm{\∞}}}{{\mathrm{RT}}^{*}}$

同时可以得到参考雷诺数，如下：

${\mathrm{Re}}_x^{*}=\frac{{\mathrm{\ρ}}^{*}{V}_{\mathrm{\∞}}x}{{\mathrm{\μ}}^{*}}$

斯坦顿数通过雷诺比拟式求解得出：

${\mathrm{St}}^{*}=\frac{c_f^{*}}{2}\frac{1}{{\left(\mathrm{Pr}\right)}^{2/3}}$

最终高超声速飞行器表面热流通过以下公式求得：

Q_aero＝St^*ρ^*V_∞c_ρ(T_r-T_ω)

沿着机身x方向将机身分成若干等份，假设每份所处的热流密度一定，则热流场为若干个不相互连续的离散温度值构成，在初始环境温度一定时，计算得到热流场分布；

(4)将热流密度场加载到有限元模型上，计算其有限元模型的稳态温度场分布：考虑流经飞行器表面的气流，飞行器表面任一点的τ_ω等价与相应的可压流的τ_ω，不可压缩气流的温度升高至一个给定的参考温度值T^*：

T^*＝0.5T_ω+0.22T_r+0.28T_∞

其中：

T^*表示飞行器表面温度，

表示边界层外缘温度，由气动力计算部分参数可以得到，

T_r为恢复温度，

根据算出的热流密度场，在MSC.Patran中将其加载在有限元模型表面，同时加上适当的辐射值加以平衡，通过MSC.Nastran求解序列SOL153对其进行温度分布计算，得到其有限元模型的稳态温度场分布；

(5)根据稳态温度场分布情况得到稳态热应力变形后的结构参数：令结构初始单元线性刚度矩阵为K₀，热应力作用生成的多余刚度为K_σ，则结构受热效应作用后的总的应力刚度矩阵等效为

K＝K₀+K_σ

则在热效应作用下的结构振动表达式为

式中，M表示质量矩阵，

表示振型，

ω为振频，

此时的动力学方程即描述了结构在热载荷作用下结构振动特性，即简化为求上式的广义特征值问题，根据得到的温度分布情况，进一步的通过MSC.Nastran求解序列SOL153对结构进行静态分析，得到稳态热应力变形后的结构参数，包括刚度矩阵，此时的刚度矩阵即为考虑热效应下的等效刚度矩阵；

(6)建立气动网格模型：应用MSC Flightloads中的气动弹性模块，将模型划分成数目适合的气动分区，每个气动分区又会分成数目适合的气动片，获得气动网格参数；

(7)根据结构的刚度矩阵结构参数和划分的气动网格参数，使用活塞理论进行线性的频域非定常气动力计算：

初始状态下的气体压力大小为p_∞，密度为ρ_∞，声速大小为a_∞，假设在等熵条件下

$\frac{p}{p_{\mathrm{\∞}}}={\left(\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_{\mathrm{\∞}}}\right)}^{\mathrm{\γ}}$

其中，γ为气体比热比，当地流音速大小为

$a^2=\mathrm{\γ}\frac{p}{\mathrm{\ρ}}$

$p^{(-\frac{\mathrm{\γ}+1}{2\mathrm{\γ}})}dp=\frac{\mathrm{\γ}}{a_{\mathrm{\∞}}}{p_{\mathrm{\∞}}}^{\left(\frac{\mathrm{\γ}-1}{2\mathrm{\γ}}\right)}du$

两边同时积分可以得到

$\frac{2\mathrm{\γ}}{\mathrm{\γ}-1}{p}^{\left(\frac{\mathrm{\γ}+1}{2\mathrm{\γ}}\right)}=\frac{\mathrm{\γ}}{a_{\mathrm{\∞}}}{p_{\mathrm{\∞}}}^{\left(\frac{\mathrm{\γ}-1}{2\mathrm{\γ}}\right)}{v}_{\mathrm{\∞}}+C$

无穷远处的p＝p_∞，且此时v_∞＝0，得到

$C=\frac{2\mathrm{\γ}}{\mathrm{\γ}-1}{p_{\mathrm{\∞}}}^{\left(\frac{\mathrm{\γ}-1}{2\mathrm{\γ}}\right)}$

得扰动压力为

$\frac{p}{p_{\mathrm{\∞}}}={(1+\frac{\mathrm{\γ}-1}{2}\frac{v_n}{a_{\mathrm{\∞}}})}^{\frac{2\mathrm{\γ}}{\mathrm{\γ}-1}}$

三阶扰动压力表达式：

$p-p_{\mathrm{\∞}}=p_{\mathrm{\∞}}\[\mathrm{\γ}\frac{v_n}{a_{\mathrm{\∞}}}+\frac{\mathrm{\γ}(\mathrm{\γ}+1)}{4}{\left(\frac{v_n}{a_{\mathrm{\∞}}}\right)}^2+\frac{\mathrm{\γ}(\mathrm{\γ}+1)}{12}{\left(\frac{v_n}{a_{\mathrm{\∞}}}\right)}^3\]$

假设x为顺气流方向，则

$v_n=\frac{\∂Z(x,y,t)}{\∂t}+V\frac{\∂Z(x,y,t)}{\∂x}$

当结构模型为上下对称时，上式则被简化为

表示上表面的下洗

表示下表面的下洗

令

得到气体压差的表达式为

从MSC.Patran中导出结构的刚度矩阵结构参数和划分的气动网格参数，并编程实现活塞理论计算其非定常气动力；

(8)根据(5)得到的稳态热应力变形后的结构参数和(7)得到的频域非定常气动力，使用工程计算方法进行颤振计算，得到颤振速度：根据有限元法建立结构动力学方程为

$M_N\stackrel{\·\·}{x}+K_{N}x=F_N$

式中：x表示节点自由度，

F_N表示气动力对应的等效节点力，

把结构动力学方程变换到模态坐标系下，其表示如下：

$\stackrel{\·\·}{q}+Kq=F$

式中，F表示等效节点气动力对应的等效模态气动力，节点自由度x与模态自由度q相互转化关系方程表达式为

x＝Φ^Tq＝∑φ_iq_i

式中

Φ＝[φ₁，φ₂，φ₃...]

通过MSC.Nastran计算得到模型相关的振型参数，根据活塞理论气动力计算方法求解得到模态坐标系下的气动力，其中第i模态时，模态力的大小可以表示为

$F_i={\∫}_s{\mathrm{\φ}}_i(x,y)F_A(x,y,t)ds=\underset{surface}{\mathrm{\Σ}}\underset{\sec tionj}{\∫}{\mathrm{\φs}}_i(x,y){\mathrm{Fs}}_A(x,y,t)ds$

式中：φ_i(x，y)表示的是第i阶连续模态，是通过φ_i插值计算得到的，

F_A(x，y)表示的是连续气动力，

φs_i(x，y)表示的是第j个气动有限单元上的第i阶连续模态，其中略去去了下标j以简化表达式，

Fs_A(x，y)表示的是第j个气动有限单元上的连续气动力，其中略去去了下标j以简化表达式，

ws(x，y，t)表示的是第j个气动有限单元上的下洗速度，由定义我们可以得到：

$ws(x,y,t)=\underset{i}{\mathrm{\Σ}}{\mathrm{\φs}}_i(x,y)q_i$

序号为j气动有限单元有四个节点构成，他们按照Nd1，Nd2，Nd3，Nd4的顺序逆时针分布，这四个节点的节点坐标依次表示为(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，(x₄，y₄)，第i阶模态对应的节点在z方向的位移可以用φ_i来表示，令其取值分别为选择的插值函数为：

$\begin{array}{c}x=\frac{1}{4}(1-r)(1-s)x_1+\frac{1}{4}(1+r)(1-s)x_2\\ +\frac{1}{4}(1+r)(1+s)x_3+\frac{1}{4}(1-r)(1+s)x_4\end{array}$

$\begin{array}{c}y=\frac{1}{4}(1-r)(1-s)y_1+\frac{1}{4}(1+r)(1-s)y_2\\ +\frac{1}{4}(1+r)(1+s)y_3+\frac{1}{4}(1-r)(1+s)y_4\end{array}$

$\begin{array}{c}{\mathrm{\φs}}_i=\frac{1}{4}(1-r)(1-s)z_1^i+\frac{1}{4}(1+r)(1-s)z_2^i\\ +\frac{1}{4}(1+r)(1+s)z_3^i+\frac{1}{4}(1-r)(1+s)z_4^i\end{array}$

将x，y投影到r，s坐标下的变换矩阵表达式为：

$\begin{array}{cc}\left(\begin{array}{c}\frac{\∂}{\∂r}\\ \frac{\∂}{\∂s}\end{array}\right)=\left(\begin{array}{cc}\frac{\∂x}{\∂r}& \frac{\∂y}{\∂r}\\ \frac{\∂x}{\∂s}& \frac{\∂y}{\∂s}\end{array}\right)\left(\begin{array}{c}\frac{\∂}{\∂x}\\ \frac{\∂}{\∂y}\end{array}\right)& J=\left(\begin{array}{cc}\frac{\∂x}{\∂r}& \frac{\∂y}{\∂r}\\ \frac{\∂x}{\∂s}& \frac{\∂y}{\∂s}\end{array}\right)\end{array}$

为简便计算对公式进行进一步化简，令

$\frac{\∂ws}{\∂t}=\underset{i}{\mathrm{\Σ}}{\mathrm{\φs}}_i(x,y){\stackrel{\·}{q}}_i$

$\begin{array}{c}\frac{\∂ws}{\∂x}=\underset{i}{\mathrm{\Σ}}{q}_i\[\begin{array}{cc}\frac{\∂r}{\∂x}& \frac{\∂s}{\∂x}\end{array}\]\left[\begin{array}{c}\frac{\∂{\mathrm{\φs}}_i(x,y)}{\∂r}\\ \frac{\∂{\mathrm{\φs}}_i(x,y)}{\∂s}\end{array}\right]\\ =\underset{i}{\mathrm{\Σ}}{q}_i(\frac{\∂r}{\∂x}\frac{\∂{\mathrm{\φs}}_i(x,y)}{\∂r}+\frac{\∂s}{\∂x}\frac{\∂{\mathrm{\φs}}_i(x,y)}{\∂s})\end{array}$

在任一气动有限单元中为常数，气动力可以进一步简化，假设

$Cp=(\frac{2p_{\mathrm{\∞}}\mathrm{\γ}}{a_{\mathrm{\∞}}}+\frac{p_{\mathrm{\∞}}M\mathrm{\γ}(1+\mathrm{\γ})}{a_{\mathrm{\∞}}}\frac{\∂f(x,y)}{\∂x})$

其中气动力系数在任一气动有限单元片中为定值，无需进行积分运算，对于模态载荷

$F_i=\underset{surface}{\mathrm{\Σ}}\underset{\sec tionj}{\∫}{\mathrm{\φs}}_i(x,y){\mathrm{Fs}}_A(x,y,t)ds,$

对式中的部分进行运算，

过程如下：

因为划分的有限单元数量庞大，我们选用简单高效的高斯积分法，同时符合3阶精度要求，选用四点进行积分，其中：

r₁＝0.57735，r₂＝-0.57735，s₁＝0.57735，s₂＝-0.57735，

通过MSC.Nastran计算得到的函数f′(x，y)是一个与厚度有关的量，因为选用的气动面单元是平面，则气动面单元的函数表达式为

f(x，y)＝a+bx+cy

式中参变量系数可以通过拟合的方法求出，其中

$\frac{\∂f(x,y)}{\∂x}=b$

当气动有限单元片为三角形单元时，通过假设第三个节点域第四个节点重合来处理，进行颤振计算后得到颤振速度以及对应的颤振频率。