一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法
【专利摘要】一种基于Newmark方法的内埋武器舱舱门气动弹性动力响应数值计算方法，将内埋武器舱打开和关闭过程中折叠舱门的动态响应视为两部分的叠加：一部分为电机驱动舱门机构运动所带来的刚体位移；另一部分为舱门结构在气动力作用下产生的弹性变形。在本方法中，将活塞理论作为气动弹性分析中的气动力模型，利用气动刚度矩阵和气动阻尼矩阵对气动力进行描述，借鉴四连杆机构的运动原理来分析舱门机构的运动规律，然后通过修正气动刚度矩阵和气动阻尼矩阵的方法来计及舱门的刚体运动对气动力产生的影响，根据Newmark方法进行舱门结构气动弹性动力学方程的数值求解。本方法考虑了武器舱打开和关闭动态过程中舱门的刚体运动对气动力产生的影响，提高了折叠舱门气动弹性动力响应分析的精确性，为折叠舱门气动弹性动力响应的数值计算提供了新思路。
【专利说明】
一种基于Newmark方法的内埋武器舱折叠舱门气;动弹性动力 响应数值计算方法
技术领域
[0001 ]本发明涉及内埋武器舱舱门气动弹性动力响应分析领域，特别涉及一种基于 Newmark方法的内埋武器舱折叠舱门气动弹性动力响应数值计算方法。
【背景技术】
[0002] 气动弹性如今在航空航天领域有着重要的应用，特别是在一些刚度小、速度高的 飞行器设计过程中更需要对气动弹性问题进行分析和计算。气动弹性动力响应问题是分析 结构在惯性力、弹性力和气动力三者作用下的动态响应历程。气动弹性动力响应问题的分 析方法可分为频域方法和时域方法两类。频域方法应用频域气动力模型，求解动力学方程 的频域形式，得到结构响应的功率谱;时域方法则采用时域非定常或准定常气动力理论，在 时域内求解动力学方程，得到随时间变化的结构动力学响应历程。
[0003] 活塞理论是一种常用的时域准定常气动力理论，适用于翼剖面厚度较薄且飞行马 赫数较高的情况，因此在壁板类结构的气动弹性问题中得到了广泛的应用。活塞理论认为 翼面上某一点的扰动对其他点产生的影响十分微弱，故可以忽略这种影响，认为翼型上某 一点的压力只与该点处的下洗速度有关，如同活塞在一圆管道运动，活塞上作用的压力只 与活塞的运动速度有关。在来流马赫数在七~ 5的范围内时，适用一阶活塞理论：
[0004]
(16)
[0005] 当翼型厚度对气动力产生的影响不可忽略时，应当采用二阶活塞理论进行气动力 计算；当来流马赫数达到Ma>5的高超声速范围时，应该采用非定常的高阶活塞理论来计算 气动力。
[0006] Newmark法属于积分类型的动力数值分析方法，是一种时域分析方法。其核心思想 是对时间步距内加速度的分布做出适当的假设，然后通过积分获得速度反应和位移反应的 表达式，进而求得步距末点的反应值。最常用的两种Newmark法为平均常加速度法和线性加 速度法:平均常加速度法相当于假定在各步距内加速度为常数，且取值为时间步始、末两点 加速度的平均值;线性加速度法则假设步距内加速度服从线性分布。针对如下单自由度体 系的运动方程：
[0007] mu{t)-\-cii{i) + ku = p(t) (17)
[0008] 其中m、c和k分别为体系的质量、阻尼和刚度特性，t为时间，/《/；)、A⑴和u(t)分别 为体系的加速度、速度和位移，P(t)为外载荷。将时间计算域[(^te3nd]离散化得到各离散点 t〇，tl，t2,…，tn，并记时间步距Ati = ti + l_ti。在已知t = ti及其以前时刻的全部响应的情况 下，为了得至1J t = ti+i时刻的动力响应，Newmark法假设：
[0009] :(18:)
[0010]其中γ和β是控制数值分析精度和稳定性的参数，当满足於0.25, γ》0.25(0.5+ β)2时，Newmark法是无条件稳定的。当Φ
时为平均常加速度模式，当取
时为线性加速度模式。
[0011]就飞行器的投放物而言，其携带形式主要有外挂式、保形式和内埋式。武器舱内埋 可以叫嚣超音速飞行时的阻力，从而提高飞行器的飞行性能，并且能够减少雷达反射面进 而提高飞行器的隐身能力。根据现代飞行器的发展趋势，机动性、隐身性等性能必不可少。 因此，内埋式武器舱在现代飞行器的设计中得到了十分广泛的应用。如今，内埋式武器舱舱 门的动弹性问题的研究也有相当的基础，然而这些研究大多关注于舱门打开或关闭状态下 的气动弹性分析。针对武器舱开关门动态过程中的气动弹性问题的研究目前相对少见。
【发明内容】
[0012] 本发明要解决的技术问题为:填补内埋武器舱折叠舱门气动弹性研究的空白，提 供一种针对武器舱门开关门动态过程的折叠舱门气动弹性动力响应数值计算方法，该方法 考虑了武器舱开关门过程中，电机驱动舱门机构运动造成的气动力变化，并利用Newmark方 法求解气动弹性运动方程，实现了内埋武器舱开关门过程中的折叠舱门气动弹性动力响应 分析，利用工程算法进行气动力计算，提高了计算效率，节省了计算时间。
[0013] 本发明解决上述技术问题采用技术方案为：一种基于Newmark方法的内埋武器舱 折叠舱门气动弹性动力响应数值计算方法，包括以下步骤：
[0014] (1)根据内埋武器舱折叠舱门的结构特点，在商业软件中建立折叠舱门的几何模 型，折叠舱门的两部分分别称为大舱门和小舱门，大舱门一边与机身铰接，另一边与小舱门 铰接，小舱门再通过两根支撑杆连接到机身上，为方便进行步骤(7)和步骤(8)中的插值，舱 门主体用面进行建模；
[0015] (2)将折叠舱门等效为一个四连杆机构，大、小舱门分别视为四连杆机构中的连 杆，根据四连杆机构的运动原理分析舱门机构的运动规律，由给出的大舱门转动角速度的 变化历程得到小舱门转动角速度和平动速度的变化历程；
[0016] (3)在商业软件中对步骤（1)中建立的舱门私.几何模型进行网格划分，舱门主体 由壳单元构成，并赋予相应的材料属性，然后利用商业软件的功能提取出舱门有限元模型 的总体刚度矩阵和总体质量矩阵以及结点坐标；
[0017] (4)为适应气动力计算，建立舱门的升力面模型，升力面与舱门所在平面重合，并 对升力面进行网格划分，并提取升力面结点坐标；
[0018] (5)根据步骤(4)得到的升力面结点坐标，并结合大气密度、马赫数飞行参数生成 初始总体气动刚度矩阵Ko q和初始总体气动阻尼矩阵CQq，为步骤(7)中的结构动力学分析做 准备，初始状态下总体气动刚度矩阵和总体气动阻尼矩阵的具体形式可根据气动力模型， 结合Hami Iton原理进行推导得到；
[0019] (6)根据当前时刻舱门做刚体运动的平动速度和转动角速度，分析舱门与来流之 间的相对速度，进而得到舱门的攻角α和气流偏角炉，并对步骤(5)中的初始总体气动刚度矩 阵Ko q和初始总体气动阻尼矩阵Coq进行修正，得到当前时刻下的总体气动刚度矩阵Kq和总体 气动阻尼矩阵C q;
[0020] (7)根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，将步骤(9)得 到的舱门结构当前时刻的结点位移向量w和结点速度向量#通过样条插值方法插值到升力 面上，得到升力面的结点位移向量Wq和结点速度向量 1^^，并根据下式计算气动力：
[0021]
(19)
[0022]式中Fq为作用于升力面结点的气动载荷列向量；
[0023] (8)根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，通过样条插 值方法将步骤(7)中的气动载荷列向量插值到舱门有限元模型上，结合步骤(3)得到的结构 总体刚度矩阵、总体质量矩阵和步骤(6)得到的总体气动刚度矩阵和总体气动阻尼矩阵，并 根据Newmark方法进行舱门气动弹性动力学方程的求解，得到舱门结构下一时刻的响应，包 括结点的位移向量w、速度向量W和加速度向量#;
[0024] (9)判断当前时刻是否达到设置的结束时间t^，即是否满足：
[0025] t ^ tend (20)
[0026] 若不满足，则转到步骤(6)，时间步增加1，继续进行下一时间步的舱门气动弹性响 应分析;若满足，则认为已完成〇~W1时间段内舱门的气动弹性响应计算，根据步骤⑻得 到的舱门结构各时刻的气动弹性响应，输出舱门在〇~^ nd时间段内的气动弹性动力响应历 程，实现内埋武器舱舱门气动弹性动力响应数值计算。
[0027]所述步骤(2)中，将舱门机构等效为一个四连杆机构，根据四连杆机构的运动原理 计算舱门机构的运动规律，将支撑杆与机身的铰接点视为坐标原点，X方向水平向右，电机 驱动大舱门转动，若任意时刻大舱门与水平方向的夹角Θ2及其转动角速度為已知，支撑杆 的运动规律可由下面两式计算：
[0030] 其中1:为支撑杆长度，I2为小舱门宽度，I3为大舱门宽度，（xD，y D)为大舱门与机身 铰接点的坐标，θ为支撑杆与水平方向的夹角j为支撑杆的转动角速度。
[0031] 根据下式求解小舱门与水平方向的夹角θ1:
[0032]
(23)
[0033] Q1对时间求导即可得到小舱门的转动角速度小舱门的平动速度可以由大舱门 的转动角速度4来确定：
[0034]
(24)
[0035] 所述步骤(5)具体实现如下：采用准定常气动力理论中的一阶活塞理论作为气动 力模型，作用在舱门上的压强为p，无穷远处来流的静压为!)"，它们之间的差值Λρ可表示 为：
[0036]
Ρ5)
[0037] 其中Ma代表来流马赫数，V代表来流速度，Ρ~为无穷远处来流的大气密度，q代表来 流动压，并满5
w为舱门横向振动位移，坐标系X轴沿来流方向，t代表时间；
[0038] 将活塞理论与Hami Iton原理结合，并引入有限元相关知识，可以得到单元上外力虚功为：
[0040] N为单元形函数组成的行向量，I为单元结点位移列向量，乾!为其关于时间的导
[0039] (26) 数，X和y为单元局部坐标系的坐标，λ为常数，并且满足
因此可以定义初始单 元气动刚度矩阵Koqe和初始单元气动阻尼矩阵Coqe为如下形式：
[0041 ]
(27)
[0042] 然后结合单元的结点坐标，对初始单元气动刚度矩阵Koqe3和初始单元气动阻尼矩 阵Coqe进行等参变换。假设四个结点的坐标为( Xl，yi)i = l，2,3,4。那么单元局部坐标系的坐 标X和y可以表示为：
[0043]
[0044] 形状函数M i的具体形式如下：
[0045]
Pi))
[0046] 其中ξ和η为等参单元的自然坐标，取值范围为-1到1。等参变换后的初始单元气动 刚度矩阵f Oqe和初始单元气动阻尼矩阵C Oqe为如下形式：
[0047] (:3〇)
[0048]
[0049]
[0050]将等参变换后的初始单元气动刚度矩阵C _和初始单元气动阻尼矩阵(/(^进行 组装得到初始总体气动刚度矩阵Koq和初始总体气动阻尼矩阵C0q。
[0051] 所述步骤(6)中，对初始总体气动刚度矩阵和初始总体气动阻尼矩阵进行修正为：
[0052] 根据当前时刻舱门相对于来流的攻角α和气流偏角免对活塞理论进行修正，气动载 荷计算公式如下：
[0053]
[0054]并且在舱门结构上附加大小为P~c~V sina的横向气动载荷，c"为无穷远处来流的 声速;对等参变换后的初始单元气动刚度矩阵f _和初始单元气动阻尼矩阵C _进行修正 后，得到当前时刻下的单元气动刚度矩阵f @和单元气动阻尼矩阵(/qe3为如下形式：
[00551
[0056] 当前时刻下的单元气动刚度矩阵和单元气动阻尼矩阵Cqe组装后即可得到当 前时刻下的总体刚度矩阵K q和总体气动阻尼矩阵Cq。
[0057] 本发明的有益效果是:本发明提供了内埋武器舱折叠舱门气动弹性动力响应数值 计算的新思路，在折叠舱门的气动弹性动力响应分析中，考虑电机驱动舱门机构运动造成 的气动力变化，通过修正气动刚度矩阵和气动阻尼矩阵的方法来对这种变化进行定量化表 征，提高了内埋武器舱开关门过程中折叠舱门的气动弹性动力响应数值计算的精确性；同 时，利用工程算法进行气动力计算，提高了计算效率，节省了计算时间。
【附图说明】
[0058] 图1为折叠舱门简化模型示意图；
[0059] 图2为折叠舱门简化为四连杆机构示意图；
[0060] 图3为大舱门转动角速度和角度变化曲线图；
[0061] 图4为小舱门转动角速度和角度变化曲线图；
[0062] 图5为小舱门上某点处横向位移响应曲线图；
[0063] 图6为小舱门各时刻横向变形对比图；
[0064]图7为本发明方法实现流程图。
【具体实施方式】
[0065]如图7所示，本发明提出了一种基于Newmark方法的内埋武器舱折叠舱门气动弹性 动力响应数值计算方法，包括以下步骤：
[0066] (1)建立舱门简化模型如图1所示，大舱门长2米，宽0.8米，小舱门长2米，宽0.5米， 支撑杆长度为1.04米，舱门主体用面进行建模，舱门内侧的筋板和肋用线进行建模。
[0067] (2)将折叠舱门简化为一个四连杆机构，如图2所示，大舱门对应的连杆长度与大 舱门宽度相同，小舱门对应的连杆长度与小舱门宽度相同，支撑杆对应的连杆长度与支撑 杆长度相同。坐标原点为图2中A点，X轴方向沿水平方向，支撑杆与水平方向夹角设为Θ，小 舱门和大舱门与水平方向的夹角分别为Θ#ΡΘ 2，以逆时针方向为正。本实例中计算内埋武器 舱开门过程中小舱门的气动弹性动力响应，给出大舱门转动角速度服从梯形分布，其转动 角速度和角度变化曲线如图3所示，小舱门的平动速度根据大、小舱门铰接点的运动速度得 到，借鉴四连杆机构的运动原理可以计算出小舱门转动角速度的变化历程，如图4所示，左 图给出了小舱门转动角速度随时间的变化曲线，右图给出了小舱门与水平方向夹角随时间 的变化曲线。
[0068] (3)将小舱门几何模型导入商业软件MSC. PATRAN中并完成有限元模型的建立，舱 门主体用壳单元进行建模，材料设置为复合材料层合板，筋板和肋用一维梁单元进行建模， 材料设置为钢。将小舱门有限元模型的总体刚度矩阵和总体质量矩阵输出到结果文件中， 并编写程序进行读取。
[0069] (4)在商业软件MSC.PATRAN中建立升力面，升力面与小舱门所在平面重合，然后对 升力面进行网格划分，并提取升力面网格的结点坐标。
[0070] (5)给出飞行高度为11千米，飞行马赫数为2,此高度的大气密度为0.364千克/立 方米，声速为295米/秒，利用权利要求3中的方法，结合升力面的结点坐标，根据单元等参变 换的思想，建立单元气动刚度矩阵和单元气动阻尼矩阵，然后组装总体气动刚度矩阵和气 动阻尼矩阵。
[0071] (6)根据小舱门当前时刻的平动速度和转动角速度，计算该状态下小舱门的攻角α 和气流偏角免，然后对气动刚度矩阵和气动阻尼矩阵进行修正，并计算小舱门上的横向附加 载荷。
[0072] (7)将小舱门当前时刻的结点位移列向量和结点速度列向量插值到升力面上，并 进行气动力计算，得到作用于升力面结点的气动载荷列向量。
[0073] (8)通过插值方法将气动力加载到小舱门有限元模型的结点上，本发明实例中大 舱门的刚度远大于小舱门的刚度，因此在小舱门的气动弹性分析中可以将大舱门视为刚 体，小舱门与支撑杆和大舱门的铰接点可视为简支约束，并利用Newmark方法计算小舱门下 一时刻的结点位移、速度和加速度响应。
[0074] (9)判断当前时刻是否达到设置的结束时间t^，即是否满足：
[0075] t ^ tend (34)
[0076] 若不满足，则转到步骤(6)，时间步增加1，继续进行舱门的气动弹性响应分析。
[0077] 重复步骤(6)~(9)，直至完成0~tend时间段内小舱门的气动弹性响应计算，根据 小舱门结构各时刻的气动弹性响应，可以得到小舱门在〇~Wi时间段内的气动弹性动力响 应历程，包括小舱门上所有结点在0~时间段内各自由度的位移变化情况。根据小舱门 上某一点处各时刻的横向位移可以绘制图5所示曲线，从图中可以判断出该点在整个响应 历程中出现的最大位移为0.0878667毫米以及出现最大位移的时间为0.972秒;根据小舱门 各结点在t = 0，l，2,3s时的横向位移，可以绘制图6所示图像，从图中可以看出小舱门的整 体变形情况。
[0078]以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制；其可扩展应 用于内埋武器舱折叠舱门气动弹性设计领域，凡采用等同变换或者等效替换而形成的技术 方案，均落在本发明权利保护范围之内。
【主权项】
1. 一种基于Newmark方法的内埋武器舱舱门气动弹性动力响应数值计算方法，其特征 在于实现步骤如下： (1) 根据内埋武器舱折叠舱门的结构特点，在商业软件中建立折叠舱门的几何模型，折 叠舱门的两部分分别称为大舱门和小舱门，大舱门一边与机身铰接，另一边与小舱门铰接， 小舱门再通过两根支撑杆连接到机身上，为方便进行步骤(7)和步骤(8)中的插值，舱门主 体用面进行建模； (2) 将折叠舱门等效为一个四连杆机构，大、小舱门分别视为四连杆机构中的连杆，根 据四连杆机构的运动原理分析舱门机构的运动规律，由给出的大舱门转动角速度的变化历 程得到小舱门转动角速度和平动速度的变化历程； (3) 在商业软件中对步骤（1)中建立折叠舱门的几何模型进行网格划分，舱门主体由壳 单元构成，并赋予相应的材料属性，然后利用商业软件的功能提取出舱门有限元模型的总 体刚度矩阵和总体质量矩阵以及结点坐标； (4) 为适应气动力计算，建立舱门的升力面模型，升力面与舱门所在平面重合，并对升 力面进行网格划分，并提取升力面结点坐标； (5) 根据步骤(4)得到的升力面结点坐标，并结合大气密度、马赫数飞行参数生成初始 总体气动刚度矩阵K〇q和初始总体气动阻尼矩阵C Qq，为步骤(7)中的结构动力学分析做准 备，初始状态下总体气动刚度矩阵和总体气动阻尼矩阵的具体形式可根据气动力模型，结 合Hamilton原理进行推导得到； (6) 根据当前时刻舱门做刚体运动的平动速度和转动角速度，分析舱门与来流之间的 相对速度，进而得到舱门的攻角a和气流偏角於，并对步骤(5)中的初始总体气动刚度矩阵 K〇q和初始总体气动阻尼矩阵CQq进行修正，得到当前时刻下的总体气动刚度矩阵Kq和总体气 动阻尼矩阵C q; (7) 根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，将步骤(9)得到的 舱门结构当前时刻的结点位移向量w和结点速度向量脅通过样条插值方法插值到升力面 上，得到升力面的结点位移向量Wq和结点速度向量并根据下式计算气动力：式中Fq为作用于升力面结点的气动载荷列向量； (8) 根据步骤(3)中的舱门结点坐标和步骤(4)中的升力面结点坐标，通过样条插值方 法将步骤(7)中的气动载荷列向量插值到舱门有限元模型上，结合步骤(3)得到的结构总体 刚度矩阵、总体质量矩阵和步骤(6)得到的总体气动刚度矩阵和总体气动阻尼矩阵，并根据 Newmark方法进行舱门气动弹性动力学方程的求解，得到舱门结构下一时刻的响应，包括结 点的位移向量w、速度向量#和加速度向量# (9) 判断当前时刻是否达到设置的结束时间t^，即是否满足： t^tend (2) 若不满足，则转到步骤(6 )，时间步增加1，继续进行下一时间步的舱门气动弹性响应分 析;若满足，则认为已完成〇~时间段内舱门的气动弹性响应计算，根据步骤(8)得到的 舱门结构各时刻的气动弹性响应，输出舱门在〇~tmd时间段内的气动弹性动力响应历程， 实现内埋武器舱舱门气动弹性动力响应数值计算。2. 根据权利要求1所述的一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力 响应数值计算方法，其特征在于:所述步骤(2)中，将舱门机构等效为一个四连杆机构的过 程为，根据四连杆机构的运动原理计算舱门机构的运动规律，将支撑杆与机身的铰接点视 为坐标原点，x方向水平向右，电机驱动大舱门转动，若任意时刻大舱门与水平方向的夹角 0 2及其转动角速度成已知，支撑杆的运动规律由下面两式计算：其中11为支撑杆长度，12为小舱门宽度，13为大舱门宽度，（XD，yD)为大舱门与机身铰接 点的坐标，9为支撑杆与水平方向的夹角，4为支撑杆的转动角速度； 根据下式求解小舱门与水平方向的夹角01:0:对时间求导即可得到小舱门的转动角速度為:，小舱门的平动速度可以由大舱门的转 动角速度疼来确定：3. 根据权利要求1所述的一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力 响应数值计算方法，其特征在于:所述步骤(5)具体实现如下:采用准定常气动力理论中的 一阶活塞理论作为气动力模型，作用在舱门上的压强为P，无穷远处来流的静压为1)~，它们 之间的差值A p表示为：其中Ma代表来流马赫数，V代表来流速度，口~为无穷远处来流的大气密度，q代表来流动 压，为舱门横向振动位移，坐标系x轴沿来流方向，t代表时间； 将活塞理论与Hamilton原理结合，并引入有限元相关知识，得到单元上外力虚功为：N为单元形函数组成的行向量，^为单元结点位移列向量，#。为其关于时间的导数，x和 y为单元局部坐标系的坐标，A为常数，并且满足 ，定义初始单元气动刚度矩阵 K〇qe和初始单元气动阻尼矩阵C〇qe为如下形式： 然后结合单元的结点坐标，对初始单元气动刚度矩阵Ko#和初始单元气动阻尼矩阵Co# 进行等参变换，假设四个结点的坐标为(^，71)1 = 1，2,3,4，单元局部坐标系的坐标1和7表 示为：形状函数M :的具体形式如下：其中I和n为等参单元的自然坐标，取值范围为-1到1，等参变换后的初始单元气动刚度 矩阵f oqe和初始单元气动阻尼矩阵C oqe为如下形式：其中雅各比矩阵J的形式如下：将等参变换后的初始单元气动刚度矩阵f Qqe和初始单元气动阻尼矩阵咖进行组装 得到初始总体气动刚度矩阵Koq和初始总体气动阻尼矩阵C0q。4.根据权利要求1所述的一种基于Newmark方法的内埋武器舱折叠舱门气动弹性动力 响应数值计算方法，其特征在于:所述步骤(6)中，对初始总体气动刚度矩阵和初始总体气 动阻尼矩阵进行修正为： 根据当前时刻舱门相对于来流的攻角a和气流偏角伊对活塞理论进行修正，气动载荷计 算公式如下：并且在舱门结构上附加大小Sina的横向气动载荷，c"为无穷远处来流的声 速;对等参变换后的初始单元气动刚度矩阵foqe和初始单元气动阻尼矩阵(/咖进行修正 后，得到当前时刻下的单元气动刚度矩阵f @和单元气动阻尼矩阵qe5为如下形式：当前时刻下的单元气动刚度矩阵K'e和单元气动阻尼矩阵C'ji装后即可得到当前时 刻下的总体刚度矩阵Kq和总体气动阻尼矩阵Cq。
【文档编号】G06F19/00GK106055860SQ201610286167
【公开日】2016年10月26日
【申请日】2016年5月3日
【发明人】邱志平, 张泽晟, 王晓军, 耿新宇, 蔡逸如, 郑宇宁, 姜南
【申请人】北京航空航天大学