CN106055860A - 一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法 - Google Patents

Abstract

一种基于Newmark方法的内埋武器舱舱门气动弹性动力响应数值计算方法，将内埋武器舱打开和关闭过程中折叠舱门的动态响应视为两部分的叠加：一部分为电机驱动舱门机构运动所带来的刚体位移；另一部分为舱门结构在气动力作用下产生的弹性变形。在本方法中，将活塞理论作为气动弹性分析中的气动力模型，利用气动刚度矩阵和气动阻尼矩阵对气动力进行描述，借鉴四连杆机构的运动原理来分析舱门机构的运动规律，然后通过修正气动刚度矩阵和气动阻尼矩阵的方法来计及舱门的刚体运动对气动力产生的影响，根据Newmark方法进行舱门结构气动弹性动力学方程的数值求解。本方法考虑了武器舱打开和关闭动态过程中舱门的刚体运动对气动力产生的影响，提高了折叠舱门气动弹性动力响应分析的精确性，为折叠舱门气动弹性动力响应的数值计算提供了新思路。

Description

一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力 响应数值计算方法

技术领域

本发明涉及内埋武器舱舱门气动弹性动力响应分析领域，特别涉及一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法。

背景技术

气动弹性如今在航空航天领域有着重要的应用，特别是在一些刚度小、速度高的飞行器设计过程中更需要对气动弹性问题进行分析和计算。气动弹性动力响应问题是分析结构在惯性力、弹性力和气动力三者作用下的动态响应历程。气动弹性动力响应问题的分析方法可分为频域方法和时域方法两类。频域方法应用频域气动力模型，求解动力学方程的频域形式，得到结构响应的功率谱；时域方法则采用时域非定常或准定常气动力理论，在时域内求解动力学方程，得到随时间变化的结构动力学响应历程。

活塞理论是一种常用的时域准定常气动力理论，适用于翼剖面厚度较薄且飞行马赫数较高的情况，因此在壁板类结构的气动弹性问题中得到了广泛的应用。活塞理论认为翼面上某一点的扰动对其他点产生的影响十分微弱，故可以忽略这种影响，认为翼型上某一点的压力只与该点处的下洗速度有关，如同活塞在一圆管道运动，活塞上作用的压力只与活塞的运动速度有关。在来流马赫数在的范围内时，适用一阶活塞理论：

$\mathrm{\Δ}p=-\frac{2q}{\sqrt{{\mathrm{Ma}}^2-1}}(\frac{\∂w}{\∂x}+\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\frac{1}{V}\frac{\∂w}{\∂t})---\left(16\right)$

当翼型厚度对气动力产生的影响不可忽略时，应当采用二阶活塞理论进行气动力计算；当来流马赫数达到Ma>5的高超声速范围时，应该采用非定常的高阶活塞理论来计算气动力。

Newmark法属于积分类型的动力数值分析方法，是一种时域分析方法。其核心思想是对时间步距内加速度的分布做出适当的假设，然后通过积分获得速度反应和位移反应的表达式，进而求得步距末点的反应值。最常用的两种Newmark法为平均常加速度法和线性加速度法：平均常加速度法相当于假定在各步距内加速度为常数，且取值为时间步始、末两点加速度的平均值；线性加速度法则假设步距内加速度服从线性分布。针对如下单自由度体系的运动方程：

$m\stackrel{\·\·}{u}\left(t\right)+c\stackrel{\·}{u}\left(t\right)+ku=p\left(t\right)---\left(17\right)$

其中m、c和k分别为体系的质量、阻尼和刚度特性，t为时间，和u(t)分别为体系的加速度、速度和位移，p(t)为外载荷。将时间计算域[0,t_end]离散化得到各离散点t₀,t₁,t₂,…,t_n，并记时间步距△t_i＝t_i+1-t_i。在已知t＝t_i及其以前时刻的全部响应的情况下，为了得到t＝t_i+1时刻的动力响应，Newmark法假设：

$\begin{array}{c}\stackrel{\·}{u}\left(t_{i+1}\right)=\stackrel{\·}{u}\left(t_i\right)+\left(1-\mathrm{\γ}\right){\mathrm{\Δt}}_i\stackrel{\·\·}{u}\left(t_i\right)+{\mathrm{\γ\Δt}}_i\stackrel{\·\·}{u}\left(t_{i+1}\right)\\ u\left(t_{i+1}\right)=u\left(t_i\right)+{\mathrm{\Δt}}_i\stackrel{\·}{u}\left(t_i\right)+\left(\frac{1}{2}-\mathrm{\β}\right){\left({\mathrm{\Δt}}_i\right)}^2\stackrel{\·\·}{u}\left(t_i\right)+\mathrm{\β}{\left({\mathrm{\Δt}}_i\right)}^2\stackrel{\·\·}{u}\left(t_{i+1}\right)\end{array}---\left(18\right)$

其中γ和β是控制数值分析精度和稳定性的参数，当满足β≥0.25,γ≥0.25(0.5+β)²时，Newmark法是无条件稳定的。当取时为平均常加速度模式，当取时为线性加速度模式。

就飞行器的投放物而言，其携带形式主要有外挂式、保形式和内埋式。武器舱内埋可以叫嚣超音速飞行时的阻力，从而提高飞行器的飞行性能，并且能够减少雷达反射面进而提高飞行器的隐身能力。根据现代飞行器的发展趋势，机动性、隐身性等性能必不可少。因此，内埋式武器舱在现代飞行器的设计中得到了十分广泛的应用。如今，内埋式武器舱舱门的动弹性问题的研究也有相当的基础，然而这些研究大多关注于舱门打开或关闭状态下的气动弹性分析。针对武器舱开关门动态过程中的气动弹性问题的研究目前相对少见。

发明内容

本发明要解决的技术问题为：填补内埋武器舱折叠舱门气动弹性研究的空白，提供一种针对武器舱门开关门动态过程的折叠舱门气动弹性动力响应数值计算方法，该方法考虑了武器舱开关门过程中，电机驱动舱门机构运动造成的气动力变化，并利用Newmark方法求解气动弹性运动方程，实现了内埋武器舱开关门过程中的折叠舱门气动弹性动力响应分析，利用工程算法进行气动力计算，提高了计算效率，节省了计算时间。

本发明解决上述技术问题采用技术方案为：一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法，包括以下步骤：

(1)根据内埋武器舱折叠舱门的结构特点，在商业软件中建立折叠舱门的几何模型，折叠舱门的两部分分别称为大舱门和小舱门，大舱门一边与机身铰接，另一边与小舱门铰接，小舱门再通过两根支撑杆连接到机身上，为方便进行步骤(7)和步骤(8)中的插值，舱门主体用面进行建模；

(2)将折叠舱门等效为一个四连杆机构，大、小舱门分别视为四连杆机构中的连杆，根据四连杆机构的运动原理分析舱门机构的运动规律，由给出的大舱门转动角速度的变化历程得到小舱门转动角速度和平动速度的变化历程；

(3)在商业软件中对步骤(1)中建立的舱门几何模型进行网格划分，舱门主体由壳单元构成，并赋予相应的材料属性，然后利用商业软件的功能提取出舱门有限元模型的总体刚度矩阵和总体质量矩阵以及结点坐标；

(4)为适应气动力计算，建立舱门的升力面模型，升力面与舱门所在平面重合，并对升力面进行网格划分，并提取升力面结点坐标；

(5)根据步骤(4)得到的升力面结点坐标，并结合大气密度、马赫数飞行参数生成初始总体气动刚度矩阵K_0q和初始总体气动阻尼矩阵C_0q，为步骤(7)中的结构动力学分析做准备，初始状态下总体气动刚度矩阵和总体气动阻尼矩阵的具体形式可根据气动力模型，结合Hamilton原理进行推导得到；

(6)根据当前时刻舱门做刚体运动的平动速度和转动角速度，分析舱门与来流之间的相对速度，进而得到舱门的攻角α和气流偏角并对步骤(5)中的初始总体气动刚度矩阵K_0q和初始总体气动阻尼矩阵C_0q进行修正，得到当前时刻下的总体气动刚度矩阵K_q和总体气动阻尼矩阵C_q；

(7)根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，将步骤(9)得到的舱门结构当前时刻的结点位移向量w和结点速度向量通过样条插值方法插值到升力面上，得到升力面的结点位移向量w_q和结点速度向量并根据下式计算气动力：

式中F_q为作用于升力面结点的气动载荷列向量；

(8)根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，通过样条插值方法将步骤(7)中的气动载荷列向量插值到舱门有限元模型上，结合步骤(3)得到的结构总体刚度矩阵、总体质量矩阵和步骤(6)得到的总体气动刚度矩阵和总体气动阻尼矩阵，并根据Newmark方法进行舱门气动弹性动力学方程的求解，得到舱门结构下一时刻的响应，包括结点的位移向量w、速度向量和加速度向量

(9)判断当前时刻是否达到设置的结束时间t_end，即是否满足：

t≥t_end (20)

若不满足，则转到步骤(6)，时间步增加1，继续进行下一时间步的舱门气动弹性响应分析；若满足，则认为已完成0～t_end时间段内舱门的气动弹性响应计算，根据步骤(8)得到的舱门结构各时刻的气动弹性响应，输出舱门在0～t_end时间段内的气动弹性动力响应历程，实现内埋武器舱舱门气动弹性动力响应数值计算。

所述步骤(2)中，将舱门机构等效为一个四连杆机构，根据四连杆机构的运动原理计算舱门机构的运动规律，将支撑杆与机身的铰接点视为坐标原点，x方向水平向右，电机驱动大舱门转动，若任意时刻大舱门与水平方向的夹角θ₂及其转动角速度已知，支撑杆的运动规律可由下面两式计算：

$\begin{array}{c}{x}_{D}cos\mathrm{\θ}-2l_1{y}_{D}sin\mathrm{\θ}-2l_1{l}_3{\mathrm{cos\θcos\θ}}_2-2l_1{l}_3{\mathrm{sin\θsin\θ}}_2\\ =l_2^2-l_4^2-l_3^2-l_1^2\end{array}---\left(21\right)$

$\begin{array}{c}-2l_3{x}_D{\mathrm{sin\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_{2}2l_3{y}_D{\mathrm{cos\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_2+2l_1{x}_D\sin \mathrm{\θ}\stackrel{\·}{\mathrm{\θ}}-2l_1{y}_D\cos \mathrm{\θ}\stackrel{\·}{\mathrm{\θ}}+2l_1{l}_3{\mathrm{sin\θcos\θ}}_2\stackrel{\·}{\mathrm{\θ}}\\ +2l_1{l}_3{\mathrm{cos\θsin\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_2-2l_1{l}_3{\mathrm{cos\θsin\θ}}_2\stackrel{\·}{\mathrm{\θ}}-2l_1{l}_3{\mathrm{sin\θcos\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_2=0\end{array}---\left(22\right)$

其中l₁为支撑杆长度，l₂为小舱门宽度，l₃为大舱门宽度，(x_D,y_D)为大舱门与机身铰接点的坐标，θ为支撑杆与水平方向的夹角，为支撑杆的转动角速度。

根据下式求解小舱门与水平方向的夹角θ₁：

${\mathrm{\θ}}_1=arctan\frac{l_{1}sin\mathrm{\θ}-y_D-l_3{\mathrm{sin\θ}}_2}{l_{1}cos\mathrm{\θ}-x_D-l_3{\mathrm{cos\θ}}_2}---\left(23\right)$

θ₁对时间求导即可得到小舱门的转动角速度小舱门的平动速度可以由大舱门的转动角速度来确定：

$\begin{array}{c}{v}_x=-{\stackrel{\·}{\mathrm{\θ}}}_2{l}_2{\mathrm{sin\θ}}_2\\ v_y={\stackrel{\·}{\mathrm{\θ}}}_2{l}_2{\mathrm{cos\θ}}_2\end{array}---\left(24\right)$

所述步骤(5)具体实现如下：采用准定常气动力理论中的一阶活塞理论作为气动力模型，作用在舱门上的压强为p，无穷远处来流的静压为p_∞，它们之间的差值△p可表示为：

$\mathrm{\Δ}p=p-p_{\mathrm{\∞}}=-\frac{2q}{\sqrt{{\mathrm{Ma}}^2-1}}(\frac{\∂w}{\∂x}+\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\frac{1}{V}\frac{\∂w}{\∂t})---\left(25\right)$

其中Ma代表来流马赫数，V代表来流速度，ρ_∞为无穷远处来流的大气密度，q代表来流动压，并满足w为舱门横向振动位移，坐标系x轴沿来流方向，t代表时间；

将活塞理论与Hamilton原理结合，并引入有限元相关知识，可以得到单元上外力虚功为：

$\begin{array}{c}\mathrm{\δ}U=\∫\∫\mathrm{\λ}(\frac{\∂N^T}{\∂x}{w}_e+\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\frac{1}{V}{N}^T{\stackrel{\·}{w}}_e)N^T{\mathrm{\δw}}_{e}dxdy\\ ={{\mathrm{\δw}}_e}^T\∫\∫\mathrm{\λ}\left(N\frac{\∂N^T}{\∂x}\right)w_{e}dxdy+{{\mathrm{\δw}}_e}^T\∫\∫\frac{\mathrm{\λ}}{V}\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\left({\mathrm{NN}}^T\right){\stackrel{\·}{w}}_{e}dxdy\end{array}---\left(26\right)$

N为单元形函数组成的行向量，w_e为单元结点位移列向量，为其关于时间的导数，x和y为单元局部坐标系的坐标，λ为常数，并且满足因此可以定义初始单元气动刚度矩阵K_0qe和初始单元气动阻尼矩阵C_0qe为如下形式：

$\begin{array}{c}{K}_{0qe}=\mathrm{\λ}\∫\∫N^T\frac{\∂N}{\∂x}dxdy\\ C_{0qe}=\frac{\mathrm{\λ}}{V}\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\∫\∫N^{T}Ndxdy\end{array}---\left(27\right)$

然后结合单元的结点坐标，对初始单元气动刚度矩阵K_0qe和初始单元气动阻尼矩阵C_0qe进行等参变换。假设四个结点的坐标为(x_i,y_i)i＝1,2,3,4。那么单元局部坐标系的坐标x和y可以表示为：

$\left\{\begin{array}{c}x\\ y\end{array}\right\}=\left[\begin{array}{cccccccc}{N_1}^{\′}& 0& {N_2}^{\′}& 0& {N_3}^{\′}& 0& {N_4}^{\′}& 0\\ 0& {N_1}^{\′}& 0& N_2& 0& {N_3}^{\′}& 0& {N_4}^{\′}\end{array}\right]\left\{\begin{array}{c}{x}_1\\ y_1\\ x_2\\ y_2\\ x_3\\ y_3\\ x_4\\ y_4\end{array}\right\}---\left(28\right)$

形状函数N′_i的具体形式如下：

$\begin{array}{c}{N_1}^{\′}=\frac{1}{4}(1-\mathrm{\ξ})(1-\mathrm{\η}),{N_2}^{\′}=\frac{1}{4}(1+\mathrm{\ξ})(1-\mathrm{\η})\\ {N_3}^{\′}=\frac{1}{4}(1+\mathrm{\ξ})(1+\mathrm{\η}),{N_4}^{\′}=\frac{1}{4}(1-\mathrm{\ξ})(1+\mathrm{\η})\end{array}---\left(29\right)$

其中ξ和η为等参单元的自然坐标，取值范围为-1到1。等参变换后的初始单元气动刚度矩阵K′_0qe和初始单元气动阻尼矩阵C′_0qe为如下形式：

$\begin{array}{c}{K}_{0qe}^{\′}=\mathrm{\λ}\∫\∫N^T(\frac{\∂N}{\∂\mathrm{\ξ}}\frac{\∂\mathrm{\ξ}}{\∂x}+\frac{\∂N}{\∂\mathrm{\η}}\frac{\∂\mathrm{\η}}{\∂x})\left|J\right|d\mathrm{\ξ}d\mathrm{\η}\\ C_{0qe}^{\′}=\frac{\mathrm{\λ}}{V}\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\∫\∫N^{T}N\left|J\right|d\mathrm{\ξ}d\mathrm{\η}\end{array}---\left(30\right)$

其中雅各比矩阵J的形式如下：

$J=\left[\begin{array}{cc}\frac{\∂x}{\∂\mathrm{\ξ}}& \frac{\∂y}{\∂\mathrm{\ξ}}\\ \frac{\∂x}{\∂\mathrm{\η}}& \frac{\∂y}{\∂\mathrm{\η}}\end{array}\right]---\left(31\right)$

将等参变换后的初始单元气动刚度矩阵K′_0qe和初始单元气动阻尼矩阵C′_0qe进行组装得到初始总体气动刚度矩阵K_0q和初始总体气动阻尼矩阵C_0q。

所述步骤(6)中，对初始总体气动刚度矩阵和初始总体气动阻尼矩阵进行修正为：

根据当前时刻舱门相对于来流的攻角α和气流偏角对活塞理论进行修正，气动载荷计算公式如下：

并且在舱门结构上附加大小为ρ_∞c_∞V sinα的横向气动载荷，c_∞为无穷远处来流的声速；对等参变换后的初始单元气动刚度矩阵K′_0qe和初始单元气动阻尼矩阵C′_0qe进行修正后，得到当前时刻下的单元气动刚度矩阵K′_qe和单元气动阻尼矩阵C′_qe为如下形式：

当前时刻下的单元气动刚度矩阵K′_qe和单元气动阻尼矩阵C′_qe组装后即可得到当前时刻下的总体刚度矩阵K_q和总体气动阻尼矩阵C_q。

本发明的有益效果是：本发明提供了内埋武器舱折叠舱门气动弹性动力响应数值计算的新思路，在折叠舱门的气动弹性动力响应分析中，考虑电机驱动舱门机构运动造成的气动力变化，通过修正气动刚度矩阵和气动阻尼矩阵的方法来对这种变化进行定量化表征，提高了内埋武器舱开关门过程中折叠舱门的气动弹性动力响应数值计算的精确性；同时，利用工程算法进行气动力计算，提高了计算效率，节省了计算时间。

附图说明

图1为折叠舱门简化模型示意图；

图2为折叠舱门简化为四连杆机构示意图；

图3为大舱门转动角速度和角度变化曲线图；

图4为小舱门转动角速度和角度变化曲线图；

图5为小舱门上某点处横向位移响应曲线图；

图6为小舱门各时刻横向变形对比图；

图7为本发明方法实现流程图。

具体实施方式

如图7所示，本发明提出了一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法，包括以下步骤：

(1)建立舱门简化模型如图1所示，大舱门长2米，宽0.8米，小舱门长2米，宽0.5米，支撑杆长度为1.04米，舱门主体用面进行建模，舱门内侧的筋板和肋用线进行建模。

(2)将折叠舱门简化为一个四连杆机构，如图2所示，大舱门对应的连杆长度与大舱门宽度相同，小舱门对应的连杆长度与小舱门宽度相同，支撑杆对应的连杆长度与支撑杆长度相同。坐标原点为图2中A点，x轴方向沿水平方向，支撑杆与水平方向夹角设为θ，小舱门和大舱门与水平方向的夹角分别为θ₁和θ₂，以逆时针方向为正。本实例中计算内埋武器舱开门过程中小舱门的气动弹性动力响应，给出大舱门转动角速度服从梯形分布，其转动角速度和角度变化曲线如图3所示，小舱门的平动速度根据大、小舱门铰接点的运动速度得到，借鉴四连杆机构的运动原理可以计算出小舱门转动角速度的变化历程，如图4所示，左图给出了小舱门转动角速度随时间的变化曲线，右图给出了小舱门与水平方向夹角随时间的变化曲线。

(3)将小舱门几何模型导入商业软件MSC.PATRAN中并完成有限元模型的建立，舱门主体用壳单元进行建模，材料设置为复合材料层合板，筋板和肋用一维梁单元进行建模，材料设置为钢。将小舱门有限元模型的总体刚度矩阵和总体质量矩阵输出到结果文件中，并编写程序进行读取。

(4)在商业软件MSC.PATRAN中建立升力面，升力面与小舱门所在平面重合，然后对升力面进行网格划分，并提取升力面网格的结点坐标。

(5)给出飞行高度为11千米，飞行马赫数为2，此高度的大气密度为0.364千克/立方米，声速为295米/秒，利用权利要求3中的方法，结合升力面的结点坐标，根据单元等参变换的思想，建立单元气动刚度矩阵和单元气动阻尼矩阵，然后组装总体气动刚度矩阵和气动阻尼矩阵。

(6)根据小舱门当前时刻的平动速度和转动角速度，计算该状态下小舱门的攻角α和气流偏角然后对气动刚度矩阵和气动阻尼矩阵进行修正，并计算小舱门上的横向附加载荷。

(7)将小舱门当前时刻的结点位移列向量和结点速度列向量插值到升力面上，并进行气动力计算，得到作用于升力面结点的气动载荷列向量。

(8)通过插值方法将气动力加载到小舱门有限元模型的结点上，本发明实例中大舱门的刚度远大于小舱门的刚度，因此在小舱门的气动弹性分析中可以将大舱门视为刚体，小舱门与支撑杆和大舱门的铰接点可视为简支约束，并利用Newmark方法计算小舱门下一时刻的结点位移、速度和加速度响应。

(9)判断当前时刻是否达到设置的结束时间t_end，即是否满足：

t≥t_end (34)

若不满足，则转到步骤(6)，时间步增加1，继续进行舱门的气动弹性响应分析。

重复步骤(6)～(9)，直至完成0～t_end时间段内小舱门的气动弹性响应计算，根据小舱门结构各时刻的气动弹性响应，可以得到小舱门在0～t_end时间段内的气动弹性动力响应历程，包括小舱门上所有结点在0～t_end时间段内各自由度的位移变化情况。根据小舱门上某一点处各时刻的横向位移可以绘制图5所示曲线，从图中可以判断出该点在整个响应历程中出现的最大位移为0.0878667毫米以及出现最大位移的时间为0.972秒；根据小舱门各结点在t＝0,1,2,3s时的横向位移，可以绘制图6所示图像，从图中可以看出小舱门的整体变形情况。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制；其可扩展应用于内埋武器舱折叠舱门气动弹性设计领域，凡采用等同变换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。

Claims (4)

1.一种基于Newmark方法的内埋武器舱舱门气动弹性动力响应数值计算方法，其特征在于实现步骤如下：

(1)根据内埋武器舱折叠舱门的结构特点，在商业软件中建立折叠舱门的几何模型，折叠舱门的两部分分别称为大舱门和小舱门，大舱门一边与机身铰接，另一边与小舱门铰接，小舱门再通过两根支撑杆连接到机身上，为方便进行步骤(7)和步骤(8)中的插值，舱门主体用面进行建模；

(2)将折叠舱门等效为一个四连杆机构，大、小舱门分别视为四连杆机构中的连杆，根据四连杆机构的运动原理分析舱门机构的运动规律，由给出的大舱门转动角速度的变化历程得到小舱门转动角速度和平动速度的变化历程；

(3)在商业软件中对步骤(1)中建立折叠舱门的几何模型进行网格划分，舱门主体由壳单元构成，并赋予相应的材料属性，然后利用商业软件的功能提取出舱门有限元模型的总体刚度矩阵和总体质量矩阵以及结点坐标；

(4)为适应气动力计算，建立舱门的升力面模型，升力面与舱门所在平面重合，并对升力面进行网格划分，并提取升力面结点坐标；

(5)根据步骤(4)得到的升力面结点坐标，并结合大气密度、马赫数飞行参数生成初始总体气动刚度矩阵K_0q和初始总体气动阻尼矩阵C_0q，为步骤(7)中的结构动力学分析做准备，初始状态下总体气动刚度矩阵和总体气动阻尼矩阵的具体形式可根据气动力模型，结合Hamilton原理进行推导得到；

(6)根据当前时刻舱门做刚体运动的平动速度和转动角速度，分析舱门与来流之间的相对速度，进而得到舱门的攻角α和气流偏角并对步骤(5)中的初始总体气动刚度矩阵K_0q和初始总体气动阻尼矩阵C_0q进行修正，得到当前时刻下的总体气动刚度矩阵K_q和总体气动阻尼矩阵C_q；

(7)根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，将步骤(9)得到的舱门结构当前时刻的结点位移向量w和结点速度向量通过样条插值方法插值到升力面上，得到升力面的结点位移向量w_q和结点速度向量并根据下式计算气动力：

$F_q=K_q{w}_q+C_q{\stackrel{\·}{w}}_q---\left(1\right)$

式中F_q为作用于升力面结点的气动载荷列向量；

(8)根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，通过样条插值方法将步骤(7)中的气动载荷列向量插值到舱门有限元模型上，结合步骤(3)得到的结构总体刚度矩阵、总体质量矩阵和步骤(6)得到的总体气动刚度矩阵和总体气动阻尼矩阵，并根据Newmark方法进行舱门气动弹性动力学方程的求解，得到舱门结构下一时刻的响应，包括结点的位移向量w、速度向量和加速度向量

(9)判断当前时刻是否达到设置的结束时间t_end，即是否满足：

t≥t_end (2)

若不满足，则转到步骤(6)，时间步增加1，继续进行下一时间步的舱门气动弹性响应分析；若满足，则认为已完成0～t_end时间段内舱门的气动弹性响应计算，根据步骤(8)得到的舱门结构各时刻的气动弹性响应，输出舱门在0～t_end时间段内的气动弹性动力响应历程，实现内埋武器舱舱门气动弹性动力响应数值计算。

2.根据权利要求1所述的一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法，其特征在于：所述步骤(2)中，将舱门机构等效为一个四连杆机构的过程为，根据四连杆机构的运动原理计算舱门机构的运动规律，将支撑杆与机身的铰接点视为坐标原点，x方向水平向右，电机驱动大舱门转动，若任意时刻大舱门与水平方向的夹角θ₂及其转动角速度已知，支撑杆的运动规律由下面两式计算：

$\begin{array}{c}{x}_{D}cos\mathrm{\θ}-2l_1{y}_{D}sin\mathrm{\θ}-2l_1{l}_3{\mathrm{cos\θcos\θ}}_2-2l_1{l}_3{\mathrm{sin\θsin\θ}}_2\\ =l_2^2-l_4^2-l_3^2-l_1^2\end{array}---\left(3\right)$

$\begin{array}{c}-2l_3{x}_D{\mathrm{sin\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_2+2l_3{y}_D{\mathrm{cos\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_2+2l_1{x}_D\sin \mathrm{\θ}\stackrel{\·}{\mathrm{\θ}}-2l_1{y}_D\cos \mathrm{\θ}\stackrel{\·}{\mathrm{\θ}}+2l_1{l}_3{\mathrm{sin\θcos\θ}}_2\stackrel{\·}{\mathrm{\θ}}\\ +2l_1{l}_3{\mathrm{cos\θsin\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_2-2l_1{l}_3{\mathrm{cos\θsin\θ}}_2\stackrel{\·}{\mathrm{\θ}}-2l_1{l}_3{\mathrm{sin\θcos\θ}}_2{\stackrel{\·}{\mathrm{\θ}}}_2=0\end{array}---\left(4\right)$

其中l₁为支撑杆长度，l₂为小舱门宽度，l₃为大舱门宽度，(x_D,y_D)为大舱门与机身铰接点的坐标，θ为支撑杆与水平方向的夹角，为支撑杆的转动角速度；

根据下式求解小舱门与水平方向的夹角θ₁：

${\mathrm{\θ}}_1=\arctan \frac{l_{1}sin\mathrm{\θ}-y_D-l_3{\mathrm{sin\θ}}_2}{l_{1}cos\mathrm{\θ}-x_D-l_3{\mathrm{cos\θ}}_2}---\left(5\right)$

θ₁对时间求导即可得到小舱门的转动角速度小舱门的平动速度可以由大舱门的转动角速度来确定：

$\begin{array}{c}{v}_x=-{\stackrel{\·}{\mathrm{\θ}}}_2{l}_2{\mathrm{sin\θ}}_2\\ v_y={\stackrel{\·}{\mathrm{\θ}}}_2{l}_2{\mathrm{cos\θ}}_2\end{array}---\left(6\right).$

3.根据权利要求1所述的一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法，其特征在于：所述步骤(5)具体实现如下：采用准定常气动力理论中的一阶活塞理论作为气动力模型，作用在舱门上的压强为p，无穷远处来流的静压为p_∞，它们之间的差值Δp表示为：

$\mathrm{\Δ}p=p-p_{\mathrm{\∞}}=-\frac{2q}{\sqrt{{\mathrm{Ma}}^2-1}}(\frac{\∂w}{\∂x}+\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\frac{1}{V}\frac{\∂w}{\∂t})---\left(7\right)$

其中Ma代表来流马赫数，V代表来流速度，ρ_∞为无穷远处来流的大气密度，q代表来流动压，并满足w为舱门横向振动位移，坐标系x轴沿来流方向，t代表时间；

将活塞理论与Hamilton原理结合，并引入有限元相关知识，得到单元上外力虚功为：

$\begin{array}{c}\mathrm{\δ}U=\∫\∫\mathrm{\λ}\left(\frac{\∂N^T}{\∂x}{w}_e+\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\frac{1}{V}{N}^T{\stackrel{\·}{w}}_e\right)N^T{\mathrm{\δw}}_{e}dxdy\\ ={{\mathrm{\δw}}_e}^T\∫\∫\mathrm{\λ}\left(N\frac{\∂N^T}{\∂x}\right)w_{e}dxdy+{{\mathrm{\δw}}_e}^T\∫\∫\frac{\mathrm{\λ}}{V}\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\left({\mathrm{NN}}^T\right){\stackrel{\·}{w}}_{e}dxdy\end{array}---\left(8\right)$

N为单元形函数组成的行向量，w_e为单元结点位移列向量，为其关于时间的导数，x和y为单元局部坐标系的坐标，λ为常数，并且满足定义初始单元气动刚度矩阵K_0qe和初始单元气动阻尼矩阵C_0qe为如下形式：

$\begin{array}{c}{K}_{0qe}=\mathrm{\λ}\∫\∫N^T\frac{\∂N}{\∂x}dxdy\\ C_{0qe}=\frac{\mathrm{\λ}}{V}\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\∫\∫N^{T}Ndxdy\end{array}---\left(9\right)$

然后结合单元的结点坐标，对初始单元气动刚度矩阵K_0qe和初始单元气动阻尼矩阵C_0qe进行等参变换，假设四个结点的坐标为(x_i,y_i)i＝1,2,3,4，单元局部坐标系的坐标x和y表示为：

$\left\{\begin{array}{c}x\\ y\end{array}\right\}=\left[\begin{array}{cccccccc}{N_1}^{\′}& 0& {N_2}^{\′}& 0& {N_3}^{\′}& 0& {N_4}^{\′}& 0\\ 0& {N_1}^{\′}& 0& {N_2}^{\′}& 0& {N_3}^{\′}& 0& {N_4}^{\′}\end{array}\right]\left\{\begin{array}{c}{x}_1\\ y_1\\ x_2\\ y_2\\ x_3\\ y_3\\ x_4\\ y_4\end{array}\right\}---\left(10\right)$

形状函数N′_i的具体形式如下：

$\begin{array}{c}{N_1}^{\′}=\frac{1}{4}(1-\mathrm{\ξ})(1-\mathrm{\η}),{N_2}^{\′}=\frac{1}{4}(1+\mathrm{\ξ})(1-\mathrm{\η})\\ {N_3}^{\′}=\frac{1}{4}(1+\mathrm{\ξ})(1+\mathrm{\η}),{N_4}^{\′}=\frac{1}{4}(1-\mathrm{\ξ})(1+\mathrm{\η})\end{array}---\left(11\right)$

其中ξ和η为等参单元的自然坐标，取值范围为-1到1，等参变换后的初始单元气动刚度矩阵K′_0qe和初始单元气动阻尼矩阵C′_0qe为如下形式：

$\begin{array}{c}{K}_{0qe}^{\′}=\mathrm{\λ}\∫\∫N^T(\frac{\∂N}{\∂\mathrm{\ξ}}\frac{\∂\mathrm{\ξ}}{\∂x}+\frac{\∂N}{\∂\mathrm{\η}}\frac{\∂\mathrm{\η}}{\∂x})\left|J\right|d\mathrm{\ξ}d\mathrm{\η}\\ C_{0qe}^{\′}=\frac{\mathrm{\λ}}{V}\frac{{\mathrm{Ma}}^2-2}{{\mathrm{Ma}}^2-1}\∫\∫N^{T}N\left|J\right|d\mathrm{\ξ}d\mathrm{\η}\end{array}---\left(12\right)$

其中雅各比矩阵J的形式如下：

$J=\left[\begin{array}{cc}\frac{\∂x}{\∂\mathrm{\ξ}}& \frac{\∂y}{\∂\mathrm{\ξ}}\\ \frac{\∂x}{\∂\mathrm{\η}}& \frac{\∂y}{\∂\mathrm{\η}}\end{array}\right]---\left(13\right)$

将等参变换后的初始单元气动刚度矩阵K′_0qe和初始单元气动阻尼矩阵C′_0qe进行组装得到初始总体气动刚度矩阵K_0q和初始总体气动阻尼矩阵C_0q。

4.根据权利要求1所述的一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法，其特征在于：所述步骤(6)中，对初始总体气动刚度矩阵和初始总体气动阻尼矩阵进行修正为：

根据当前时刻舱门相对于来流的攻角α和气流偏角对活塞理论进行修正，气动载荷计算公式如下：

并且在舱门结构上附加大小为ρ_∞c_∞V sinα的横向气动载荷，c_∞为无穷远处来流的声速；对等参变换后的初始单元气动刚度矩阵K′_0qe和初始单元气动阻尼矩阵C′_0qe进行修正后，得到当前时刻下的单元气动刚度矩阵K′_qe和单元气动阻尼矩阵C′_qe为如下形式：

当前时刻下的单元气动刚度矩阵K′_qe和单元气动阻尼矩阵C′_qe组装后即可得到当前时刻下的总体刚度矩阵K_q和总体气动阻尼矩阵C_q。