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

权利要求书

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$

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