CN112580241A - 一种基于结构降阶模型的非线性气动弹性动响应分析方法 - Google Patents

Abstract

本发明属于结构动力学和气动弹性力学分析领域，特别涉及一种基于结构降阶模型的非线性气动弹性动响应分析方法。本方法采用对输入测试载荷及其对应结构变形进行回归分析的方法求解给定形式的结构非线性刚度系数，利用非线性刚度系数构建大变形结构动力学方程，在此基础上利用模态组合乘比例系数的方式构建测试载荷，通过构建位移残量基函数准确恢复机翼展向位移；最后结合结构降阶模型给出的动力学方程及非定常曲面涡格法搭建几何非线性气动弹性阵风响应求解流程。本发明兼顾求解精度、计算效率与复杂模型适用性，可应用于航空航天类飞行器几何非线性气动弹性阵风响应分析计算。

Description

一种基于结构降阶模型的非线性气动弹性动响应分析方法

技术领域

本发明属于结构动力学和气动弹性力学分析领域，特别涉及一种基于结构降阶模型的非线性气动弹性动响应分析方法。

背景技术

气动弹性力学作为应用力学的分支，主要研究气动力、弹性力与惯性力之间耦合的问题。在气动力作用下弹性结构发生振动与变形，结构弹性运动反过来又会导致气动力的大小及分布发生变化，这种相互作用会带来多种多样的气动弹性现象，包括气动弹性变形、颤振、阵风响应等。在有动力飞行的历史发展进程中，气动弹性现象扮演了重要角色，是保证飞行器安全、提高飞行器飞行性能所必须考虑的关键因素。

现代飞行器设计对于飞行性能的要求不断提高，以高空长航时无人机、飞翼布局无人机为代表的大柔性飞行器，自20世纪90年代以来逐渐成为各航空大国的研究重点。该类飞行器能够以更灵活的方式不间断的执行诸如军事侦察、民用科学探索及通信中继等任务，具有远超传统飞行器的飞行高度及续航能力，具有广阔的应用前景。由于其飞行性能的要求，大柔性飞行器通常利用较大的机翼展弦比来提高升阻比，同时在结构设计汇总增加复合材料比例。然而这些因素会导致这类大柔性飞行器相对结构重量低、柔性大，气动力加载引起的结构变形较大，带来几何非线性气动弹性问题。

当飞行器结构变形较大或结构上加载了较大的初始载荷时，结构动力学特性按非线性关系依赖于结构瞬时变形及应力状态。在力学机理上表现为结构应变与位移呈现非线性关系、受力平衡呈现非线性关系，而应力与应变的本构关系依然是线性的。在气动弹性分析中，考虑几何非线性因素后，气动弹性问题即成为几何非线性气动弹性问题。此时，线性小变形假设不再适用。一方面，大变形下其气动力分布情况与线性有明显不同；另一方面，结构动力学特性也会随着结构变形及载荷的改变而发生变化。这两方面特性导致几何非线性因素对于气动弹性分析具有重要影响，成为制约大柔性飞行器设计的关键问题。阵风响应问题是气动弹性分析中的重要问题。阵风不仅会激起飞机的刚体运动，还会激发结构的弹性振动，由此产生的动态响应会降低飞机的飞行品质，影响驾驶员操纵，甚至引发安全事故，直接影响飞行安全和飞行包线范围。而大柔性飞机由于其自身的特点使得其对于阵风扰动更加敏感。

结构建模作为气动弹性分析中的核心问题，其重要性不言而喻，准确的大变形结构建模方法是几何非线性气动弹性分析的基础。目前，气动弹性研究中通常采用大变形梁理论，难以应用于实际工程模型，而工程分析通常使用现有的商业有限元软件，模型自由度数大、非线性阶数高、收敛性差、计算效率低且难以对结构动力学特性进行机理性研究。非线性结构降阶模型能够较好的平衡求解效率、计算精度及复杂模型适用性，是一种适用于几何非线性气动弹性分析的结构模型。本申请人在中国专利CN108052772A“一种基于结构降阶模型的几何非线性静气动弹性分析方法”，CN108363843A“一种基于结构降阶模型的几何非线性静气动弹性全机配平方法”，以及文章“Static aeroelastic analysisincluding geometric nonlinearites based on reduced order model”Xie Changchuan等，《Chinese Journal of Aeronautics》,第638-650页，20170216”中，将非线性结构降阶模型引入几何非线性静气动弹性分析，构建静气动弹性分析框架，但无法处理动力学响应问题。另外，本申请人在文章“Flutter and gust response analysis of a wing modelincluding geometric nonlinearities based on a modified structural ROM”An Chao等，《Chinese Journal of Aeronautics》,第48-63页，20190731”中将非线性结构降阶模型引入颤振与阵风分析，但对于较为重要的机翼展向位移恢复中采用了计算展向位移正交基函数来表达展向位移的方式。对于复杂机翼模型而言，展向位移正交基函数并不便于求解。同时，构建降阶模型所用的测试载荷为在一定变形下的气动力载荷，计算量较大，且造成生成的结构降阶模型适用性降低。

发明内容

为了克服现有几何非线性动气动弹性阵风响应计算中结构非线性有限元计算的不足，本发明基于非线性结构降阶模型及曲面气动力进行机翼几何非线性气动弹性阵风响应计算，对现有建立的结构降阶模型中的机翼展向位移恢复方法及测试载荷选取方法进行改进，并将其应用于动气动弹性响应分析中，提供一种高效准确的针对复杂机翼模型的几何非线性气动弹性阵风响应分析方法，建立机翼几何非线性阵风响应计算流程。

特别地，本发明核心设计思想可描述如下：采用对输入测试载荷及其对应结构变形进行回归分析的方法求解给定形式的结构非线性刚度系数，利用非线性刚度系数构建大变形结构动力学方程，在此基础上利用模态组合乘比例系数的方式构建测试载荷，通过构建位移残量基函数准确恢复机翼展向位移；最后结合结构降阶模型给出的动力学方程及非定常曲面涡格法搭建几何非线性气动弹性阵风响应求解流程。综上，本发明在结构建模时综合考虑了气动力随动效应及展向位移，并且在气动建模时综合考虑了大变形及三维效应。本发明所提供的几何非线性气动弹性阵风响应分析方法兼顾求解精度，计算效率与复杂模型适用性，可应用于航空航天类飞行器几何非线性气动弹性阵风响应分析计算。

为实现上述目的，本发明提供了一种基于结构降阶模型的非线性气动弹性动响应分析方法，包括如下步骤：

S1：建立飞行器机翼模型或全机模型的有限元模型，所述有限元模型包括多个节点、多个梁单元以及多个集中质量单元，其中至少一部分节点上加载随动力，飞行器机翼的根部六自由度固支；

S2:确定所述有限元模型的结构模态，包括结构垂直弯曲、扭转及水平弯曲方向线性模态，其中每个模态由其固有频率和模态刚度表征；

S3:采用若干阶结构模态作为非线性结构降阶模型的模态基底，建立考虑包含几何非线性大变形作用的结构动力学方程：

方程(1)对应第m阶结构模态并应用了Einstein求和约定，m＝1,…,N，N为参与降阶的结构模态数目；其中，M_m对应第m阶结构模态的广义质量项；K_m对应第m阶结构模态的广义刚度项，数值已知；

为非线性刚度系数，数值未知，上标(2)和(3)表示对应的广义坐标多项式阶数，下标对应广义坐标的模态阶次；q_m,q_n,q_l,q_p分别对应第m,n,l,p阶结构模态的物理广义坐标，n＝1,…,N，l＝1,…,N，p＝1,…,N；F_m对应第m阶结构模态的广义力；

取结构动力学方程(1)的静力学部分，并将已知的线性刚度部分K_mq_m移至方程右侧得到：

S4:确定非线性结构降阶模型建立所需的测试载荷及测试变形，具体过程为：

选择所述有限元模型的弯曲及扭转线性模态，弯曲线性模态共nb阶，扭转线性模态共nt阶，组合得到如下测试载荷：

其中，

和

分别为对应第i阶弯曲线性模态和第j阶扭转线性模态的加权系数，φ_ib,i＝1,…,nb和φ_jt,j＝1,…,nt分别为选定的第i阶弯曲线性模态振型及第j阶扭转线性模态振型，下标b及t分别表示模态为弯曲线性模态和扭转线性模态；对于第i阶线性弯曲模态振型取Na个不同的加权系数进行组合，对于第j阶线性扭转模态振型取Nc个不同的加权系数进行组合。

将NT组测试载荷以随动力的形式加载于所述有限元模型，并进行非线性有限元计算，得到NT组测试变形；

S5:利用逐步回归分析方法建立非线性结构降阶模型，具体过程为：

将NT组测试载荷和测试变形输入方程(2)，给定回归分析问题：

其中，上标1,2,…,NT表示测试载荷编号；

利用逐步回归法对回归分析问题方程(4)进行回归分析，求解非线性刚度系数，得到非线性结构降阶模型；

S6:利用位移残量基函数恢复非线性结构降阶模型的准确位移，所述位移残量基函数为

Ψ＝RQ⁺ (6)

其中，Ψ为位移残量基函数；R＝[R¹,R²,R³,…,R^NT]为位移残量矩阵，位移残量R^j＝X^j-Φq^j，j＝1,2,...,NT，X^j为第j组测试载荷对应的测试位移向量，Φ为线性模态矩阵，q^j为第j组测试载荷对应的物理广义坐标；Q⁺为矩阵Q的广义逆，Q为表征位移残量的基函数广义坐标矩阵：

其中，Q^Nm，Nm＝1,2,…,NT为第Nm组测试载荷的位移残量基函数的广义坐标行向量；Q^Nm中共有M个位移残量基函数广义坐标：

其中，

表示第Nm组测试载荷的第j个位移残差基函数广义坐标，

对应第Nm组载荷的非线性结构降阶模型中的共N个结构模态广义坐标，位移残差基函数的广义坐标为降阶模型中的结构模态广义坐标逐个的二阶多项式相乘形式，即：

位移残差基函数广义坐标个数M与非线性结构降阶模型中的结构模态广义坐标个数N关系有：

S7:利用非线性结构降阶模型结合非定常曲面气动力及曲面样条插值方法确定飞行器机翼的几何非线性阵风响应。

进一步，步骤S7具体过程为：

1)输入结构气动模型数据，并给定工况条件，所述工况条件包括风速和攻角；

2)建立飞行器机翼初始结构的气动网格模型，计算结构气动插值矩阵；

3)给定迭代计算的时间步增量及终止时间步数；

4)给定飞行器机翼的阵风模型及频率信息后，利用非定常曲面涡格法求解对象的曲面非定常气动力；

5)进行力插值，利用非线性结构降阶模型计算飞行器机翼的结构变形；

6)根据步骤5)中的结构变形更新结构位移及速度，作为下一计算时间步的瞬态响应初始条件；

7)进行位移插值，更新气动网格模型作为下一计算时间步的气动力计算初始条件，进入下一计算时间步进行计算；

8)重复步骤5)～7)直至达到终止时间步数，得到飞行器机翼在步骤1)给定的工况条件下的一段时间内的阵风响应计算结果。

本发明的有益效果：

1)本发明给出基于非线性结构降阶模型的非线性气动弹性阵风响应分析方法，与传统的基于非线性有限元方法的气动弹性分析相比有效减小了计算资源消耗，提高了计算效率；

2)在建立非线性结构降阶模型中，本发明采用线性模态作为降阶基底，提高了对模型的适用性；

3)本发明利用模态组合乘比例系数的方式生成测试载荷，相对采用气动力作为测试载荷的方法计算量更低，操作简单；

4)本发明采用求位移残量基函数的方式表达展向位移，相比计算展向位移正交基的方式而言，模型适用性更好，而对于复杂模型而言，计算展向位移正交基往往很难实现。

附图说明

图1是本发明的基于结构降阶模型的非线性气动弹性动响应分析方法流程图；

图2是本发明的求解几何非线性阵风响应问题流程图；

图3是本发明实施例的机翼模型示意图；

图4是本发明实施例的机翼的有限元模型示意图；

图5是本发明实施例的机翼的线性弯曲及扭转模态振型图；

图6是本发明实施例的机翼气动力模型示意图。

具体实施方式

如图1所示，本发明的基于结构降阶模型的非线性气动弹性动响应分析方法包括如下步骤:

S1：计算初始化

建立机翼模型的有限元模型，划分有限元模型的结构有限元网格，并定义结构有限元网格的单元属性(包括梁单元的截面尺寸及法线方向)和材料特性，设置计算边界条件，即在需要固支的机翼根部边界点上进行固支约束；

S2：计算结构模态

非线性结构降阶模型基于结构的线性模态建立，在建立非线性结构降阶模型前，计算结构垂直弯曲、扭转及水平弯曲方向线性模态；

S3：确定结构动力学方程

考虑包含几何非线性大变形作用的结构动力学方程：

方程(1)对应第m阶结构模态并应用了Einstein求和约定，m＝1,…,N，N为参与降阶的模态数目；其中，M_m对应第m阶结构模态的广义质量项；K_m对应第m阶结构模态的广义刚度项，数值已知；

为非线性刚度系数，数值未知，上标(2)和(3)表示对应的广义坐标多项式阶数，下标对应广义坐标的模态阶次；q_m,q_n,q_l,q_p分别对应第m,n,l,p阶结构模态的物理广义坐标，n＝1,…,N，l＝1,…,N，p＝1,…,N；F_m对应第m阶结构模态的广义力。

将结构动力学方程(1)取静力学部分如方程(2)：

将方程(2)中的已知的线性刚度部分移至方程右侧，可得：

方程(2)即可化为回归问题进行分析，进而求解未知的非线性刚度系数

S4：确定非线性结构降阶模型建立需要的测试载荷及测试变形

测试载荷的给定是结构降阶模型建立的关键部分，测试载荷合适与否直接影响回归分析的准确性及结构降阶模型的合理性。本步骤具体过程如下：

选择有限元模型的弯曲及扭转线性模态，弯曲线性模态共nb阶，扭转线性模态共nt阶，给出组合得到的测试载荷如下：

其中，

分别为对应第i阶弯曲模态和第j阶扭转线性模态的加权系数，φ_ib(i＝1,…,nb),φ_jt(j＝1,…,nt)为选定的第i阶线性弯曲及第j阶扭转线性模态振型，下标b及t分别表示模态为弯曲及扭转模态；对于第i阶线性弯曲模态振型可取Na个不同的加权系数进行组合，对于第j阶线性扭转模态振型可取Nc个不同的加权系数进行组合。

将NT组测试载荷以随动力的形式加载于有限元模型，利用MSC.Nastran有限元软件计算机翼变形，调整比例系数大小使得最终得到的NT组机翼变形在感兴趣的非线性变形范围内，作为建立结构降阶模型需要的测试变形。

S5：利用逐步回归分析方法建立非线性结构降阶模型

假设提供NT组测试载荷，对应有NT组测试位移，将测试载荷与对应的测试位移投影至模态空间，下面给出模态空间广义坐标与物理空间位移关系：

其中，q_m,M_m分别为第m阶模态对应的物理广义坐标及广义质量项，

为物理空间下的结构质量矩阵，X为物理空间的测试位移向量。

由回归分析得到方程(3)左侧的非线性刚度系数，回归问题可表达为：

模态对应的物理广义坐标q上标表示静态测试载荷编号。对方程(6)给出的回归问题进行逐步回归分析即可求解未知的非线性刚度系数。有利地，相比最小二乘方法，本发明所采用的逐步回归方法可以避免回归共线性问题，保证模型非线性刚度系数数值准确可靠。

S6：准确位移恢复

对于大柔性机翼而言，低阶线性模态一般为弯曲模态，扭转模态及水平弯曲模态，其模态分量不能表征因为弯曲导致的展向位移，步骤S5中利用低阶线性模态建立的结构降阶模型自然也无法恢复这一部分位移，导致结构变形计算不准确。

为此，在本步骤将因弯曲导致的展向位移以位移残量表征，对于第j组测试载荷及其对应位移，给出位移残量表达式：

R^j＝X^j-Φq^j,j＝1,2,...,NT (7)

其中，X^j为第j组测试载荷对应的物理空间的测试位移向量，Φ为线性模态矩阵，Φ＝{φ₁,φ₂,…φ_j,…,φ_N}，φ_j为参与降阶的第j阶结构模态。考虑所有NT组测试载荷情况，得到位移残量矩阵为：

R＝[R¹,R²,R³,…,R^NT] (8)

给定表征位移残量的基函数广义坐标矩阵Q：

其中，Q^Nm，Nm＝1,2,…,NT为第Nm组测试载荷的位移残量基函数的广义坐标行向量；Q^Nm中共有M个位移残量基函数广义坐标：

其中，上标Nm对应第Nm组载荷求得的数据，

表示第Nm组载荷对应的第j个位移残差基函数广义坐标，

对应第Nm组载荷的非线性结构降阶模型中的共N个结构模态广义坐标。则表征位移残量的基函数表达为：

位移残差基函数广义坐标个数与非线性结构降阶模型中的结构模态广义坐标个数关系有：

因此，位移残量基函数表示为

Ψ＝RQ⁺ (12)

其中，Ψ为位移残量基函数；Q⁺为矩阵Q的广义逆。

上述位移残量基函数Ψ可以准确表征展向位移及其他由于使用线性模态进行位移恢复所导致的位移恢复损失，特别是展向位移。因此结构的完整变形u可表征为：

u＝qΦ+rΨ (13)

其中r_j＝p_j,j＝1,2,…,M。p_j为式(10)中对应结构模态广义坐标二阶多项式相乘形式的位移残量基函数广义坐标。

S7：求解几何非线性气动弹性阵风响应问题

得到非线性结构降阶模型，利用降阶模型结合非定常曲面气动力及曲面样条插值方法可以进行大柔性机翼几何非线性阵风响应的计算。大柔性机翼的变形与气流的作用力有关，而气动力又与机翼在空中的形状和位置有关，两者是互相影响相互制约的非线性关系，本发明采用分析迭代计算的方法得到大柔性机翼的阵风响应，具体过程如下：

1)输入结构气动模型数据并给定工况条件，包括风速及攻角；

2)建立机翼初始结构的气动网格模型，计算结构气动插值矩阵；

3)给定迭代计算时间步增量及终止时间步数；

4)给定阵风模型及频率信息后，利用非定常涡格法求解曲面非定常气动力；

5)进行力插值，利用非线性结构降阶模型计算机翼的结构变形；

6)根据5)中的结构变形更新结构位移及速度，作为下一计算时间步的瞬态响应初始条件；

7)进行位移插值，更新气动网格模型作为下一时间步气动力计算初始条件，进入下一时间步计算；

8)重复5)～7)过程直至达到规定的终止时间步数。所有时间步计算完毕后，即得到大柔性机翼在该工况下一段时间内的阵风响应计算结果。几何非线性阵风响应分析求解流程如图2所示。

下面结合附图和实施例进一步描述本发明，应该理解，以下所述实施例旨在便于对本发明的理解，而对其不起任何限定作用。

本实施例采用大柔性机翼，翼型选用NACA0015对称翼型，翼根完全固支；机翼参数为：翼根弦长100mm，翼尖弦长100mm，展长1000mm，后掠角0Deg，无扭转角；机翼主梁采用矩形截面，截面长35mm，截面宽1.5mm，位于机翼50％弦长处；机翼外形由8个轻木盒段来维持，每个盒段由轻木骨架与热缩膜蒙皮组成，为减小由于弯曲变形使得盒段相互挤压而产生对机翼刚度的影响，盒段与主梁单点连接并且相邻盒段之间留有一定间隙。特别地，为了调整机翼的动力学特性和结构变形，在翼尖设置了配重杆，配重杆总长200mm。机翼模型总重0.613kg，该机翼模型示意图如图3所示。

本实施例的基于结构降阶模型的几何非线性静气动弹性分析方法具体过程如下：

第一步：计算初始化

利用有限元软件MSC.Nasrtan建立机翼的有限元模型，其中，采用CBEAM单元模拟主梁及翼肋，采用CONM2单元模拟集中质量，FORCE1添加随动力结构载荷，根部六自由度固支。本实施例的有限元模型信息如表1所示，有限元模型示意图如图4所示。

表1有限元模型信息

模型节点数	213
		模型梁单元数	212
模型集中质量数	65
		加载随动力的节点数	51

第二步：计算结构模态

有限元模型建立完成后，利用Nastran软件SOL103线性模态求解器计算结构线性模态，本实施例计算前5阶模态及第8阶模态二阶扭转模态，结构线性模态信息如表2所示，振型信息如图5所示。前5阶模态中，第1阶模态为垂直一弯模态，如图5中(a)所示，固有频率1.179Hz，模态刚度5.488E+01，固有频率较低，模型柔度较大；第2阶模态为垂直二弯模态，如(b)所示，固有频率7.724Hz，模态刚度2.355E+03；第3阶模态为垂直三弯模态，如(c)所示，固有频率22.19Hz，模态刚度1.944E+04；第4阶模态为一阶扭转模态，如(d)所示，固有频率22.95Hz，模态刚度2.080E+04；第5阶模态为水平一弯模态，如(e)所示，固有频率27.47Hz，模态刚度2.978E+04。第8阶模态为二阶扭转模态，如(h)所示，固有频率95.72Hz，模态刚度3.617E+05。

表2结构线性模态信息

模态阶数	模态名称	固有频率/Hz	模态刚度
				1	垂直一弯	1.179	5.488E+01
2	垂直二弯	7.724	2.355E+03
				3	垂直三弯	22.19	1.944E+04
4	一阶扭转	22.95	2.080E+04
				5	水平一弯	27.47	2.978E+04
8	二阶扭转	95.72	3.617E+05

第三步：确定结构动力学方程

结构动力学方程如方程(1)所示，本实施例中采用前5阶线性模态及第8阶模态共6阶模态作为模态基底建立非线性结构降阶模型，方程(1)中已知项有如下数值：

M₁＝1M₂＝1M₃＝1

M₄＝1M₅＝1M₆＝1

K₁＝5.488E+01K₂＝2.355E+03K₃＝1.944E+04

K₄＝2.080E+04K₅＝2.978E+04K₆＝3.617E+05

第四步：确定非线性结构降阶模型建立需要的测试载荷及测试变形

本实施例采用前5阶线性模态及第8阶模态共6阶模态作为结构分析模态，包括：垂直一弯{φ₁}、垂直二弯{φ₂}、垂直三弯{φ₃}、一阶扭转{φ₄}、水平一弯{φ₅}、二阶扭转{φ₆}，φ₁，φ₂，φ₃，φ₄，φ₅，φ₆分别为前5阶线性模态及第8阶模态的振型。测试载荷选定为在由一阶弯曲模态和一阶扭转模态组合出的作用力。对应一阶弯曲模态的模态加权系数

依次取值0、100、200、300、400、500，对应一阶扭转模态的模态加权系数

依次取值0、3、6、9、12、15。

将每组测试载荷加载在有限元模型上并利用Nastran软件SOL106几何非线性静力学分析求解器计算，得到270组测试变形。翼端最大垂直变形243mm，接近机翼翼展的25％。

第五步：利用逐步回归分析方法建立非线性结构降阶模型

将270组测试载荷与测试变形输入结构动力学方程(1)，给定回归分析问题如方程(6)式所示，在MATLAB中利用逐步回归方法对问题(6)进行回归分析，通过回归分析求解非线性刚度系数后，得到非线性结构降阶模型。

第六步：准确位移恢复

得到非线性结构降阶模型后，利用位移残量基函数恢复模型准确位移，本实施例考虑所有270组测试载荷情况求解位移残量基函数，选择6阶模态情况下，位移残量基函数共21阶，位移残量基函数求解及位移恢复如方程(12)(13)所示。

第七步：利用非线性结构降阶模型求解几何非线性阵风响应问题

在10m/s，12m/s，14m/s三种风速状态，攻角3°工况下，对正弦离散阵风条件下的机翼阵风响应进行分析计算，阵风频率取3Hz及4Hz。

以风速10m/s为例：

1)向机翼的有限元模型输入机翼的结构气动模型数据，给定风速10m/s，攻角3°；

2)机翼的初始结构没有变形，根据机翼的有限元模型建立长1000mm，宽100mm的矩形气动面，并将该气动面共划分160个气动网格建立气动网格模型，如图6所示，其中沿机翼展向40个网格，沿机翼弦向4个网格，网格大小一致，分布均匀；之后利用MATLAB语言计算结构气动插值矩阵；

3)给定迭代时间步增量0.0025s，计算总时间5s；

4)利用MATLAB语言编写程序进行曲面气动力的求解，气动力方法选用定常曲面涡格法，在160个网格内布置涡环基本解。求解坐标系规定为：x轴沿来流方向，y轴水平向右，z轴右手定则确定，在涡格内布置涡环单元。涡环单元由四段等强度直线涡首尾相接而成，翼面自由涡由后缘涡格拖出，平行于来流方向。第一次计算气动力时，机翼气动面处于初始状态即未变形状态，之后给定阵风模型及频率信息，计算阵风响应；

5)进行力插值，利用非线性结构降阶模型计算机翼气动面的结构变形；

6)根据5)中的结构变形更新结构位移及速度，作为下一时间步的瞬态响应初始条件；

7)进行位移插值，更新气动网格模型作为下一时间步气动力计算初始条件；

8)进入下一时间步计算，重复5)～7)过程直至5s时域响应计算结束。

将本发明方法计算结果与现有的基于非线性有限元方法的几何非线性气动弹性阵风响应计算结果进行对比，验证计算准确性。表3给出了两种方法下翼尖垂向位移响应幅值的对比结果，从结果看两者误差很小。同时，在计算效率上，模型阶数低，相比基于非线性有限元方法的几何非线性气动弹性分析而言具有较大优势。

表3两种方法下翼尖垂向位移响应幅值计算结果对比

风速/m/s	阵风频率	非线性降阶模型结果/mm	有限元方法结果/mm
				10	3Hz	19.2	18.9
10	4Hz	11.3	11.2
				12	3Hz	25.1	24.7
12	4Hz	13.0	13.2
				14	3Hz	33.1	32.0
14	4Hz	17.2	16.8

综上，本发明在计算大柔性飞行器几何非线性气动弹性阵风响应中能够兼顾计算效率，求解精度及复杂模型适用性，在实际飞行器设计应用中具有重要意义。

本发明未详细描述内容为本领域技术人员公知内容。

对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以对本发明的实施例做出若干变型和改进，这些都属于本发明的保护范围。

Claims (2)

1.一种基于结构降阶模型的非线性气动弹性动响应分析方法，其特征在于，包括如下步骤：

S1：建立飞行器机翼模型或全机模型的有限元模型，所述有限元模型包括多个节点、多个梁单元以及多个集中质量单元，其中至少一部分节点上加载随动力，飞行器机翼的根部六自由度固支；

S2:确定所述有限元模型的结构模态，包括结构垂直弯曲、扭转及水平弯曲方向线性模态，其中每个模态由其固有频率和模态刚度表征；

S3:采用若干阶结构模态作为非线性结构降阶模型的模态基底，建立考虑包含几何非线性大变形作用的结构动力学方程：

方程(1)对应第m阶结构模态并应用了Einstein求和约定，m＝1,…,N，N为参与降阶的结构模态数目；其中，M_m对应第m阶结构模态的广义质量项；K_m对应第m阶结构模态的广义刚度项，数值已知；

为非线性刚度系数，数值未知，上标(2)和(3)表示对应的广义坐标多项式阶数，下标对应广义坐标的模态阶次；q_m,q_n,q_l,q_p分别对应第m,n,l,p阶结构模态的物理广义坐标，n＝1,…,N，l＝1,…,N，p＝1,…,N；F_m对应第m阶结构模态的广义力；

取结构动力学方程(1)的静力学部分，并将已知的线性刚度部分K_mq_m移至方程右侧得到：

S4:确定非线性结构降阶模型建立所需的测试载荷及测试变形，具体过程为：

选择所述有限元模型的弯曲及扭转线性模态，弯曲线性模态共nb阶，扭转线性模态共nt阶，组合得到如下测试载荷：

其中，

和

分别为对应第i阶弯曲线性模态和第j阶扭转线性模态的加权系数，φ_ib,i＝1,…,nb和φ_jt,j＝1,…,nt分别为选定的第i阶弯曲线性模态振型及第j阶扭转线性模态振型，下标b及t分别表示模态为弯曲线性模态和扭转线性模态；对于第i阶线性弯曲模态振型取Na个不同的加权系数进行组合，对于第j阶线性扭转模态振型取Nc个不同的加权系数进行组合。

将NT组测试载荷以随动力的形式加载于所述有限元模型，并进行非线性有限元计算，得到NT组测试变形；

S5:利用逐步回归分析方法建立非线性结构降阶模型，具体过程为：

将NT组测试载荷和测试变形输入方程(2)，给定回归分析问题：

其中，上标1,2,…,NT表示测试载荷编号；

利用逐步回归法对回归分析问题方程(4)进行回归分析，求解非线性刚度系数，得到非线性结构降阶模型；

S6:利用位移残量基函数恢复非线性结构降阶模型的准确位移，所述位移残量基函数为

Ψ＝RQ⁺ (6)

其中，Ψ为位移残量基函数；R＝[R¹,R²,R³,…,R^NT]为位移残量矩阵，位移残量R^j＝X^j-Φq^j，j＝1,2,...,NT，X^j为第j组测试载荷对应的测试位移向量，Φ为线性模态矩阵，q^j为第j组测试载荷对应的物理广义坐标；Q⁺为矩阵Q的广义逆，Q为表征位移残量的基函数广义坐标矩阵：

其中，Q^Nm，Nm＝1,2,…,NT为第Nm组测试载荷的位移残量基函数的广义坐标行向量；Q^Nm中共有M个位移残量基函数广义坐标：

其中，

表示第Nm组测试载荷的第j个位移残差基函数广义坐标，

对应第Nm组载荷的非线性结构降阶模型中的共N个结构模态广义坐标，位移残差基函数的广义坐标为降阶模型中的结构模态广义坐标逐个的二阶多项式相乘形式，即：

位移残差基函数广义坐标个数M与非线性结构降阶模型中的结构模态广义坐标个数N关系有：

S7:利用非线性结构降阶模型结合非定常曲面气动力及曲面样条插值方法确定飞行器机翼的几何非线性阵风响应。

2.根据权利要求1所述的方法，其特征在于，步骤S7具体过程为：

1)输入结构气动模型数据，并给定工况条件，所述工况条件包括风速和攻角；

2)建立飞行器机翼初始结构的气动网格模型，计算结构气动插值矩阵；

3)给定迭代计算的时间步增量及终止时间步数；

4)给定飞行器机翼的阵风模型及频率信息后，利用非定常曲面涡格法求解对象的曲面非定常气动力；

5)进行力插值，利用非线性结构降阶模型计算飞行器机翼的结构变形；

6)根据步骤5)中的结构变形更新结构位移及速度，作为下一计算时间步的瞬态响应初始条件；

7)进行位移插值，更新气动网格模型作为下一计算时间步的气动力计算初始条件，进入下一计算时间步进行计算；

8)重复步骤5)～7)直至达到终止时间步数，得到飞行器机翼在步骤1)给定的工况条件下的一段时间内的阵风响应计算结果。

Images