基于修正活塞理论的壁板非线性热气动弹性分析

郭静 王建民 张忠 李海波 高博

（北京强度环境研究所 可靠性与环境工程技术重点实验室，北京100076）

0 引言

高马赫数飞行的航天飞行器，气动力、热、结构多物理场相互耦合，易引发热气动弹性稳定性问题[1]。新型复合材料的广泛应用，使得壁板结构易出现极限环振荡、混沌等复杂非线性热气动弹性响应形式，影响飞行安全，亟需开展壁板结构非线性热气动弹性分析方法研究[2-3]。

大变形引起的几何非线性和非定常气动力的高效预示是非线性热气动弹性分析应考虑的主要问题。壁板颤振、超声速颤振和热颤振广泛采用简单、高效的活塞理论预示非定常气动力[4-8]，活塞理论也在应用中不断改进和发展[9]，先后提出了Van Dyke活塞理论[10]、当地流活塞理论[11]等。杨智春等[12]以超音速气流中二维受热壁板为对象，开展了非线性热颤振分析，给出了其基本型和二次失稳型两种颤振模式。夏巍等[13]用一阶活塞理论模拟壁板受到的气动力，分析了超音速气流中受热壁板后屈曲状态下的稳定性。杨翊仁等[14-16]研究了二维、三维壁板在超音速气流作用下各种复杂的热弹性颤振动力学现象，给出了二维壁板发生动态Hopf分岔的条件[17]。夏巍[18]、Mei[19]推导了考虑热效应的颤振有限元方程，研究了几何非线性引起的壁板分岔与混沌响应。

本文充分利用Nastran软件SOL400中的非线性有限元求解程序，采用DMAP语言对其进行二次开发，将温度场分析和Van Dyke修正活塞理论非定常气动力计算引入SOL400非线性瞬态响应直接求解模块，编制非线性热气动弹性分析程序，并采用PCL语言在Patran软件中开发嵌入式前后处理界面，将Van Dyke修正活塞理论非定常气动力计算的高效性与SOL400非线性瞬态响应计算的高效性相结合，建立了工程实用的非线性热气动弹性分析工具。以复合材料壁板为研究对象，进行典型温度工况下的非线性热气动弹性分析。结果表明，本文方法可对壁板结构LCO，周期，非周期和混沌等复杂非线性热气动弹性响应进行有效预示。

1 Van Dyke修正活塞理论

根据经典活塞理论，翼面上某点压力与该点下洗速度之间满足如下关系[20]

式中，p为翼面上某点的压力，p∞、ρ∞、a∞分别为无穷远来流压力，密度和声速，W为翼面下洗速度，γ为气体绝热比，空气通常取 =1.4γ。当下洗速度小于音速时，活塞理论可根据麦克劳林级数保留项阶次称为一阶、二阶或三阶活塞理论。

针对马赫数较小的超声速工况，采用Van Dyke二阶理论修正后的压力系数Cp可表示为

式中，M为马赫数，a为当地音速，则活塞理论可转化为

2 基于Patran/Nastran的非线性热气动弹性分析程序

矩阵形式的热气动弹性运动方程表达式为

式中，[M ]，[K]分别为质量矩阵和刚度矩阵，A(u)为结点力列向量，由Van Dyke修正活塞理论求得。刚度矩阵[K]由小位移线性刚度矩阵[k0]，几何刚度矩阵[kσ]和初位移矩阵[kL]叠加表示，[K ]= [k0]+ [kσ]+[kL]。刚度矩阵[K]和结点力列向量A(u)均为位移向量{u}的非线性函数，须采用直接时域积分方法进行非线性热气动弹性分析。本文基于Nastran软件中的非线性瞬态响应直接求解模块SOL400中的非线性有限元求解程序，采用DMAP语言对其进行二次开发，将温度场分析和Van Dyke修正活塞理论非定常气动力计算引入SOL400非线性瞬态响应直接求解模块，编制了非线性热气动弹性分析程序，将Van Dyke修正活塞理论求解非定常气动力的高效性和SOL400求解结构非线性瞬态响应的高效性充分结合，避免了非线性瞬态响应复杂有限元程序的编制，技术路线如图 1所示。

图1 非线性热气动弹性分析技术方案Fig.1 Technical scheme of nonlinear thermoaeroelastic analysis

SOL 400引入非定常气动力的DMAP二次开发技术路线如图 2所示，并采用PCL语言在Patran软件中开发非线性热气动弹性分析嵌入式前后处理界面，如图 3所示，建立了工程实用的非线性热气动弹性分析工具。

图2 SOL 400引入非定常气动力DMAP二次开发技术方案Fig.2 Technical scheme for SOL 400 secondary development by introducing unsteady aerodynamics with DMAP language

图3 非线性热气动弹性分析界面Fig.3 Nonlinear thermoaeroelastic analysis interface

3 数值算例

针对石墨-环氧复合材料壁板标准算例，壁板长度、宽度、厚度分别为0.381m、0.3048m、1.3e-3m，位移边界为四边固支，壁板力学性能参数和有限元模型分别如表 1和图 4所示。

表1 壁板力学性能参数Table1 Mechanical property parameters of panel

图4 壁板有限元模型Fig.4 Finite element model of panel

本文采用与文献标准算例相同的Van Dyke修正活塞理论计算非定常气动力，进行了复合材料壁板典型温度工况的非线性热气动弹性分析。取壁板展向50%站位3/4弦线位置处结点的位移和速度响应与标准算例进行对比验证，maxW表示最大位移振幅。

1）常温工况

常温工况取空气密度为1.225kg/m3，马赫数为2.0，采用本程序分析了不同飞行速度下极限环幅值随无量纲动压的变化规律，如图 5所示。由图 5可知，本文计算得到的无量纲临界颤振动压为275，文献中无量纲临界颤振动压为250，误差不超过10%，且随着无量纲动压增大，极限环幅值增大，本文与文献结果变化趋势一致。研究结果表明，本程序适用于复合材料壁板结构常温工况气动弹性分析。

图5 壁板极限环幅值随无量纲动压变化规律Fig.5 Variation of limit cycle amplitude of panel with dimensionless dynamic pressure

2）均匀温度场工况

均匀温度场工况取空气密度为1.225kg/m3，马赫数为2.0，考虑四种均匀温度场，温升分别为ΔT=Tcr,2Tcr,3Tcr,4Tcr,屈 曲 临 界 温 升 为Tcr= 21.8°C。采用本程序计算得到不同飞行速度下的非线性热气动弹性响应如图 6-图 9所示，分别对应了LCO，周期，非周期和混沌响应。非线性热气动弹性响应无量纲最大幅值与文献对比如表2所示。

表2 非线性热气动弹性响应无量纲最大幅值（ max/W h）对比 Table2 Comparison of non-dimensional maximum amplitudes of nonlinear thermoaeroelastic responses

图6 ΔT = Tcr，LCOFig.6 ΔT = Tcr , Limit cycle oscillation LCO

图7 ΔT = 2Tcr ，周期Fig.7 ΔT = 2Tcr, periodic vibration

图8 ΔT = 3Tcr ，非周期Fig.8 ΔT = 3Tcr, aperiodic vibration

图9 ΔT = 4Tcr ，混沌Fig.9 ΔT = 4Tcr, chaos

由分析可知，本程序可以对壁板结构LCO，周期，非周期和混沌等复杂非线性热气动弹性响应进行有效预示，且本文与文献结果一致性很好，误差不超过6%，表明针对Nastran SOL 400多学科非线性求解模块中的非线性瞬态响应求解器的DMAP二次开发流程是正确的。

4 结论

本文采用Van Dyke修正活塞理论分析非定常气动力，针对商用软件Nastran SOL 400多学科非线性求解模块中的非线性瞬态响应求解器，增加温度场和非定常气动力求解模块，采用PCL和DMAP语言，基于Patran/Nastran二次开发平台，编制了非线性热气动弹性分析程序，以复合材料壁板为研究对象，分析了其在典型温度工况下的非线性气动弹性响应，标准算例对比验证了本文方法的正确性。

研究结果表明，常温工况下，本文计算得到的临界颤振动压与文献结果的误差不超过10%，且变化趋势一致，且可以对复合材料壁板结构LCO，周期，非周期和混沌等复杂非线性热气动弹性响应进行有效预示。