CN110674599A - 气动伺服弹性系统非定常气动载荷的有理近似优化方法 - Google Patents

Abstract

本发明公开了气动伺服弹性系统非定常气动载荷的有理近似优化方法通过已有商用工程数值计算软件获取分析结构的两类原始数据，在此基础上，运用现代工程分析公认的非定常气动力有理近似表达形式，建立了一种线性项系数及非线性项系数两类变量的两级交替双重迭代数值优化算法。两级双重优化迭代直至达到获得高精度的有理近似拟合式中的最佳系数，提高了飞机选定结构的气动伺服弹性系统建模分析中非定常气动载荷历程的预计精度，极大降低有理近似拟合误差，从而为飞机气动伺服弹性系统的动响应、动稳定性以及主动弹性控制提供高精度载荷历程数值。

Description

气动伺服弹性系统非定常气动载荷的有理近似优化方法

技术领域

本发明涉及气动弹性力学领域，尤其涉及气动伺服弹性系统非定常气动载荷 的有理近似优化方法。

背景技术

飞机按操纵指令偏转操纵面，使飞机改变原飞行轨迹，即形成飞机的各种机 动飞行动作。飞机机动飞行时在其升力面上形成空气动力载荷历程，又称非定 常气动力，这种变化的载荷历程必然导致结构的弹性振动；反之，弹性振动会 对非定常气动力产生影响，两者耦合最不利的情况是导致升力面振幅自激放大， 轻者飞机难以控制，重者导致结构破坏。因此，工程结构设计中，非定常气动 载荷历程预计是弹性动响应、动稳定性以及主动弹性控制性能评估的前提技术 工作。

尽管现代基于计算流体力学CFD与计算结构力学CSD耦合的数值分析技术 日益成熟，但这类数值技术有两类缺陷使其难以在工程结构迭代设计过程中得 以广泛应用：其一，虽CFD与CSD耦合的数值分析精度高，但数值建模过程复 杂，完整飞机的模型计算网格规模巨大，目前的常规高性能数值计算技术仍难 适应于大型结构的非定常气动载荷历程预计及其动响应分析，效率极低；其二， 当结构柔性较大，飞行机动会导致结构振动的变形量较大，此种情形下需要CFD 与CSD耦合的动网格更新技术，目前该类算法技术仍难于工程应用。飞机结构 设计过程中，更多需要的是数值建模效率高、迭代计算耗时相对少的数值分析 技术与工具，以至于适当牺牲分析精度。

针对上述飞机设计过程中大型结构的非定常气动载荷预计及其高效动响应 或稳定性分析等问题，现行工程技术常用的方法是结构模态缩阶法和小扰动线 性非定常气动载荷法的数值计算综合。这类工程数值模型分析方法中，非定常 气动载荷是利用离散频域内偶极子面元数值算法，计算获取飞机升力面上各面 元压力点相对于某一振动频率的非定常气动力影响系数矩阵，并通过与结构固 有模态的插值关系计算获得作用于结构上的广义非定常气动力系数矩阵。在此 基础上，对各离散频率下的广义非定常气动力系数矩阵给予固定形式的有理数 值拟合逼近，以获得连续变化的广义非定常气动力系数据矩阵。文献 1“Application of aeroservoelastic modeling using minimum-state unsteadyaerodynamic approximations[J].Journal of guidance,control,and dynamics,1991,14(6):1267-1276”提出了一种最小状态有理近似拟合算法，但该算法的约束经验 性过强，拟合效果并不理想；文献2“Method for aerodynamic unsteady forces timecalculations on an F/A-18aircraft[J].Aeronautical Journal,2008,112(1127):27-32” 在具体飞机的非定常气动载荷历程计算中采用了Roger型最小二乘算法；文献 3“ZAERO User’s Manual.Version8.2,ZONA Technology Inc.March 2008”实现了 Roger型最小二乘算法的工程算法程序的工程应用。但这些算法均固定了气动力 滞后根，使广义非定常气动力系数矩阵各元素在多离散频率点上的拟合效果并 不理想。

发明内容

本发明针对飞机设计过程中非定常气动载荷预计应用中的上述不足，提出了 气动伺服弹性系统非定常气动载荷的有理近似优化方法。

为了解决上述技术问题，本发明采用以下技术方案：

气动伺服弹性系统非定常气动载荷的有理近似优化方法，包括如下步骤：

1)通过已有商用工程数值计算软件获取分析气动及结构的原始数据：构建 航空器结构动力学分析的有限元网格模型，计算并选取一定数量的结构 固有振型模态矩阵；构建气动布局的气动载荷面元网格，确定选取的多 个离散减缩频率值，计算获取相应各减缩频率下的非定常气动力影响系 数矩阵，结合两类原始数据计算可得相应各离散减缩频率的非定常广义 气动力系数矩阵；

2)预设初始参数：初始滞后根、初始Hessian矩阵和收敛精度控制值，将初 始参数以及步骤1)中的各离散减缩频率下的广义气动力系数矩阵数据送 入近似有理子程序进行计算，返回二次误差函数的有理近似拟合数据；

3)计算广义气动力系数矩阵的有理近似拟合数据，结合步骤1)中的各离散 减缩频率下广义气动力系数矩阵数据和步骤2)中二次误差函数的有理近 似拟合数据，计算总误差加权拟合函数，并结合滞后根的上下限约束， 形成无约束化目标函数，计算无约束化目标函数的函数值和导数值；

4)按数学标准方法计算拟牛顿优化搜索方向，并沿搜素方向进行一维线性 搜索计算步长因子；

5)计算新的滞后根，并将新的滞后根以及各离散减缩频率下广义气动力系 数的矩阵数据一起送入近似有理拟合子程序中，调用所述子程序计算得 到二次误差函数的近似有理拟合数据；

6)根据步骤5)中二次误差函数的有理近似拟合数据计算新的广义气动力系 数矩阵系数，并计算总误差函数；

7)将步骤6)中迭代的总误差函数与前一次迭代的总误差函数相减，取绝对 值；判断绝对值是否小于预设的收敛精度控制值：是，则停止优化迭代， 否则，用BFGS标准数学公式构造新的Hessian矩阵,令迭代计数加1,转 至步骤3)继续迭代优化。

进一步地：步骤2)中二次误差函数的近似有理拟合形式为：

其中，jk是纯虚数，k为实数域内的减缩频率变量，

(i＝0,…,N)为取 二次误差函数为有理近似形式的复数拟合系数矩阵，β_j(j＝1,…,N-2)称为滞后 根，m为含滞后根分式项的下标，工程中N的取值一般不大于6。

进一步地：所述二次误差函数为：

其中，Q(jk)为工程软件计算提供的原始广义气动力系数矩阵；

是归一化有理近似平方误差和的无约束化目标函数优化计算所得结果。

进一步地：所述归一化有理近似平方误差和为:

其中，下标ij与广义气动力系数矩阵的第i行及第j列元素相对应；上标T表示 矩阵的转置；

为归一化系数，

为原始广义气动力系数矩阵第i行及第j列元素的实部与虚部；

上式为广义气动力系数矩阵各元素的有理近似拟合的矩阵乘积形式； D_m＝β_m ²+k_i ²(m＝1,…，N-2)，当行号小于l时，i取矩阵行号，当行号大于l时， i取行号减l，l为用户选定的离散减缩频率数；右端列向量为有理近似各矩阵中 的第i行第j列元素，为待定系数；

{C}_ij为用户提高有理近似的针对性所选择的线性约束的组合系数，下标ij 的含义同前，C_0ij为用户约束的常数项；{λ}_ij为拉氏乘子列向量。

进一步地：步骤3)中所述外层迭代是将所有减缩频率下各矩阵元素的拟合 误差加权总和作为总目标函数，如下式：

其中，n为广义气动力系数矩阵的列数，m为广义气动力系数矩阵的列数，l为 选定的离散减缩频率数，w_ij为对矩阵第i行第j列元素的加权系数，可取1；F_ij(k_p)， G_ij(k_p)分别为相应各减缩频率离散值k_p时所得原始广义气动力系数矩阵各复数 元素的实部与虚部，

项为对二次误差函数有理近似形式实施 线性优化后所得的广义气动力系数矩阵各复数元素的实部与虚部。

外层非线性滞后根优化的约束条件设置为各滞后根的上下限：

0＜L_j＜β_i＜U，(j＝1，…，N-2)

其中，L_j和U_j分别为第j个滞后根β_j取值的上下限,

结合拉格朗日乘子法，含所述滞后根上下限约束的归一化有理近似平方误 差和的无约束化目标函数作为外层优化迭代的目标函数：

其中，f如前述，

为外层无约束化目标函数中第i个滞后根上下限约束的拉氏乘子；各ν_k(k＝i或j)为将不等式约束转换为等式约束的松弛变量，在优化迭代 运算中可运用目标函数关于ν_k的导数为零条件而被消除，故不作为新变量。

进一步地：所述总目标函数可写作两项纯实数的计算方式：

其中，w_ij为加权系数，可仅为1；上标DLS为有理近似拟合所得各矩阵元素与 二次误差函数有理近似拟合所得各矩阵元素的相加算子，

由该函数可得关于滞后根的导数算式为：

其中，

将上式简记为向量形式：

其中,[·]^T表示行向量的转置，即为列向量。

本发明的有益效果为：

本发明采用了一种两类变量的非线性交替优化迭代算法，建立了精确求导 公式，实现了多变量的高效导数计算；借助于近似二阶导的数值计算技术，形 成拟牛顿优化搜索方向，继而实现两类变量的交替优化，完成关注频域范围内 更高精度的非定常气动载荷预计。该方法基于离散频率下的广义非定常气动力 系数矩阵数据，在关心的频段范围内，对多点频率下的广义非定常气动力系数 矩阵元素实施固定有理近似格式的两级交替双重优化迭代拟合计算，可极大降 低数值拟合误差，从而为飞机气动伺服弹性系统的动响应、动稳定性以及主动 弹性控制提供高精度载荷历程数值。

附图说明

图1是气动伺服弹性系统非定常气动载荷的有理近似优化方法流程图；

图2是近似有理拟合子程序流程图；

图3是一机翼平面几何及骨架结构图；

图4是机翼结构有限元网格及离散质量模型图；

图5是机翼结构固有弹性对称振型模态；

图6是机翼弦平面气动面元网格模型图；

图7是机翼弦平面气动面元网格的典型振型插值图；

图8广义非定常升力系数有理拟合曲线；

图9广义非定常俯仰系数有理拟合曲线；

其中：3a为机翼平面几何形状；3b为机翼结构内部骨架；4a为机翼结构有 限元网格；4b为机翼结构离散质量分布；5a为模态1一阶对称弯曲模态；5b为 模态2二阶对称弯曲模态；5c为模态3三阶对称弯曲模态；5d为模态4四阶对 称弯曲模态；5e为模态5一扭三弯耦合对称模态；5f为模态6五阶对称弯曲模 态；5g为模态7二阶对扭转模态；5h为模态8六阶对称弯曲模态；5i为模态9 三阶对称弯曲模态；7a为模态7二阶对称扭转模态；7b模态8六阶对称弯曲模 态。

具体实施方式

本发明气动伺服弹性系统非定常气动载荷的有理近似优化方法，计算流程 见图1所示，通过已有商用工程数值计算软件获取分析结构的两类原始数据： 其一，计算获取选定离散频域下的非定常气动力影响系数矩阵数据库；其二， 完成飞机结构集体及操纵面的低阶固有频率及位移振型计算，得到结构固有振 型模态矩阵。在此基础上，运用现代工程分析公认的非定常气动力有理近似表 达形式，建立了一种线性项系数及非线性项系数两类变量的两级交替双重迭代 数值优化算法：在归一化误差平方和以及拉氏乘子转化的无约束化目标函数基 础上，建立了二次误差函数；内迭代在非线性项系数固定条件下，对拟合误差 实施双重线性项系数的最佳拟合计算；外迭代则对非定常气动力有理近似式中 的非线性项系数进行基于拟牛顿方向的优化搜索。两级双重优化迭代直至达到 获得高精度的有理近似拟合式中的最佳系数。本发明技术对于获取工程关注频 域范围内的非定常气动载荷高精度预计具有重要的工程应用价值及意义。

本发明中主要推导方法包括以下步骤：

步骤1:原始数据及理论计算关系。针对一固定翼飞机的气动布局及结构布 置，运用工程商用数值分析软件构建其气动布局的气动载荷计算面元网格模型 以及结构动力学分析的有限元网格模型，调用结构固有振动品质的特征值求解 器，即可获得式(1)中用广义模态坐标描述的飞机结构动力学运动方程左端项； 调用非定常气动载荷求解器，可以获得式(2)中的非定常气动力影响系数矩阵 [AIC]。

式(1)中，q为列向量，对应选定的飞机刚体模态及弹性模态的广义坐标， 也即各振型模态在振动过程中的时变参与系数，其上的“·”表示对其的时间导数； δ为操纵面的刚体偏转模态坐标；ρ为大气密度；V为飞机飞行速度；M_qq和M_qδ 为机体广义质量及控制面耦合质量矩阵，C_qq为广义阻尼矩阵，K_qq为机体结构 广义刚度矩阵；Q_qq和Q_qδ分别为飞机升力面振动和操纵面偏转引起的广义非定 常气动力系数矩阵。

式(2)中，Δp表示面元各网格压力点处的压强分布列向量，维数为J×1；J 为面元网格数；ρ为大气密度，V为飞行速度；[AIC]矩阵为飞机的固定升力面 和可动操纵面在某一振动频率下的联合非定常气动力影响系数矩阵，维数为 J×J，是本发明的原始数据矩阵；dw_q、dw_δ分别为飞机升力面和操纵面振动时在 相应面元网格控制点处形成的无量纲下洗速度向量，维数均为J×1，在小变形振 动条件下，与各振型模态位移呈线性关系。

式(1)中的M_qq、C_qq及K_qq阵当采用模态质量归一化的结构振型模态时，前 m阶结构振型模态阵为

其中，n为结构有限元模型的总自由度个数，m为所取的模态阶次，f_ij表示结构 模型节点第i个自由度的第j阶模态振型值。对于飞行状态，需注意飞机刚体运 动模态的选取。本发明以飞机纵向运动为例，刚体模态仅选取沉浮及俯仰模态。 飞机的沉浮模态可令各节点的向上运动位移为1即得；俯仰模态选绕飞机质心 抬头转动1个弧度单位所形成的模型各节点运动位移。以飞机质心为原点，顺 航向为飞机x轴的正向，则俯仰刚体振动模态的运动位移为f(x,y,z)＝-x，且对 应的广义模态坐标即为其时变俯仰角。

操纵面刚体偏转模态是指操纵面绕其转轴旋转单位角度形成的操纵面模型 节点位移，操纵面以外的模型节点位移为零。假定飞机有l个操纵面，相应的刚 体偏转模态矩阵为

其中，n为全部机体结构的节点位移自由度，当节点不在操纵面上，则其节点各 自由度的模态位移为零。

假定结构模型的一致质量矩阵为M_s，广义质量阵M_qq和M_qδ的表达式如下

M_qq＝F_q ^TM_sF_q＝I_n×n M_qδ＝F_q ^TM_sF_δ (5)

同理，假定结构的总刚度矩阵为K_s，结构的广义刚度矩阵K_qq为

其中，diag表示对角矩阵，即非对角元素均为零，ω_i为第i阶正交模态的圆频 率，ω_i ²为系统特征方程的特征值，刚体模态的对应值为零。广义阻尼阵C_qq为 下述形式

C_qq＝diag(2ξ₁ω₁,2ξ₂ω₂,…,2ξ_mω_m) (7)

其中，ξ_i为第i阶固有振型的模态阻尼系数。

为计算相应于各阶振型模态的面元各网格控制点处的相对下洗速度，需将结 构振型转移到气动力面元网格上，再通过式(8)和(9)计算可得用结构广义 坐标表达的时变相对下洗速度dw_q和dw_δ向量：

上两式中的G_H和G_δ分别为机体结构及操纵面表面有限元网格与其弦平面面元 网格间的几何插值转换矩阵；

为沿来流方向的导数算子，i为单位虚数， k＝ωb/U_∞称为减缩频率，ω为气流受扰的圆频率，b为升力面半弦长，U_∞为来 流速度。

联立式(2)、(8)和(9)计算可得式(1)右端相应于某气流扰动减缩频 率k的两广义非定常气动力系数矩阵：

其中，S是对角阵，是面元各网格的面积值。

观察式(10)及(11)，矩阵元素均为复数。为描述方便，可将上两式统一 简记为：

Q＝F(k)+jG(k) (12)

其中，j为单位虚数，即j²＝-1。

步骤2：广义非定常气动载荷矩阵的有理近似形式。式(12)是小扰动速度 势线化微分方程数值计算结果的形式描述，需完成包括奇异积分在内的大量数 值计算，通常由工程专业计算软件对离散减缩频率实施计算，如Zaero软件。实 际上，非定常气动载荷的历程特性与升力面的结构振动特性相耦合，即减缩频 率k并不能事先确定，故首先需采用数值拟合方法将其连续化。本发明的非线 性两级交替双重优化方法采用如下公认的非定常气动载荷有理近似描述形式：

式中，

即为式(12)的近似值；jk是纯虚数，k为实数域内的减缩频率变量； A_i(i＝0,…,N)为待定的复数形式的拟合系数矩阵，与式(12)的维数相同； β_j(j＝1,…,N-2)称为滞后根，工程中N的取值一般不大于6。

另外，式(13)中的各A_i项关于有理分式是线性的，故称为线性系数项； 而滞后根β_j是有理分式的非线性项，称为非线性系数。实际拟合优化计算时需将 式(13)变换为与式(12)相同的形式，即：

将式(14)复数的实部与虚部写成列向量形式为

其中，

为有理近似的广义气动力系数矩阵中第i行及第j列复元素的实部与虚部；右端列向量是各线性系数矩阵A_i(i＝0,…,N)的第i 行第j列元素，为待定变量。右端矩阵系数为减缩频率k及非线性系数的矩阵函 数，D_m＝β² _m+k²(m＝1,…，N-2)，当给定k值及滞后根时，此矩阵是确定的，即 当事先给定滞后根，则式(15)的拟合优化只是线性代数方程组的求解问题。

步骤3：含约束的最小二乘拟合基本算法。为提高有理近似形式的针对性， 工程上多采用以下三类约束：

1.指定某个减缩频率k的拟合值与理论值相等，包括实部及虚部相等，即

式中，D_z＝β² _z+k²，z＝1,2,…,N-2；{x}＝{a_0ij,a _1ij,…,a_(N-2)ij}，F_ij(k)和G_ij(k)分别 为式(12)广义气动力系数矩阵第i行第j列复数元素的实部与虚部。

2.约束k＝0的拟合斜率与理论值相等，即

式中，k₂通常取接近于零但大于零的某个小数值。

3.约束某个拟合系数为零(如拟合系数a_2ij为零),即

[0 0 1 0 … 0]{x}＝0 (18)

以上三种约束或保证函数值及拟合斜率的精度，或限制给定的拟合系数。

工程中常用的一种含约束最小二乘基本算法是在选定滞后根条件下的平方 和误差最小优化算法。设在工程关注的频段范围内选定了l个离散减缩频率值， 由工程数值计算软件计算得到的各离散频率下广义非定常气动力系数矩阵按式 (12)形式为

同理，式(14)各离散减缩频率下有理近似拟合的矩阵形式可记为

式中，右端的矩阵项是式(15)中取所选定的离散减缩频率所得的具体结果， D_m＝β_m ²+k_i ²(m＝1,…，N-2)，当行号小于l时，i取矩阵行号，当行号大于l时， i取行号减l，l为用户选定的离散减缩频率数；右端列向量为有理近似各矩阵中 的第i行第j列元素，为待定系数。

为后叙方便起见，将式(16)～(18)的等式约束简记为

由式(19)和(20)，工程中通常应用归一化的有理近似平方误差和，即：

式中，

其他符号见式(19)～(20)。

结合拉格朗日乘子法，归一化有理近似平方误差和的无约束化目标函数为

式中，ij对应于广义其动力学系数矩阵的第i行第j列元素，上标T表示矩阵的 转置，{λ}_ij为拉氏乘子向量，其他符号见式(19)～(21)。

步骤4：二次误差双重线性拟合优化的内迭代算法。应用式(23)的常规线 性拟合优化算法的技术效果并不理想。为此，本发明在两级交替迭代优化算法 中建立了一个二次误差函数，并作为内迭代的线性拟合运算部分。二次误差函 数是建立在式(23)的拟合优化结果基础上，可写作：

其中，Q(jk)为工程软件计算提供的原始广义气动力系数矩阵；为由式(23)拟合优化所得结果。

取二次误差函数仍为有理近似形式，即：

于是，运用各离散的减缩频率值，且固定的气动力滞后根不变，则由式(24) 可得二次误差函数的离散值，将其作为原始输入数据，同样可实施形同式(23) 的优化目标函数，应用最小二乘算法即得线性方程组求解问题。

步骤5：非线性滞后根的外迭代优化算法。实际计算发现，非线性滞后根对 拟合结果的影响甚为显著。为此，本发明在前述双重线性优化内迭代算法的基 础上，提出以下针对滞后根的非线性外层迭代算法。外层迭代是将所有减缩频 率下各矩阵元素的拟合误差加权总和作为总目标函数，如下式

式中，n为广义气动力系数矩阵的列数，m为广义气动力系数矩阵的列数，l为 选定的离散减缩频率数，w_ij为对矩阵第i行第j列元素的加权系数，可取1；F_ij(k_p)， G_ij(k_p)分别为减缩频率为离散值k_p时原始广义气动力系数矩阵各复数元素的实 部与虚部，

项为对二次误差函数有理近似形式实施线性优 化后所得的广义气动力系数矩阵各复数元素的实部与虚部。为后续描述方便， 称式(26)为广义气动力拟合系数的总误差函数。

另需注意，式(26)的

项是对式(25)实施线性双重 优化后的结果。参考式(13)和(14)的形式，结合式(25)这两项可表示为

其中，

经验可知，当非线性滞后根的取值处于所研究问题的缩减频率范围内时， 拟合结果更好。为此，外层非线性滞后根优化问题的约束条件设置为各滞后根 的上下限：

0<L_j＜β_j＜U_j,(j＝1,…,N-2) (28)

其中，L_j和U_j分别为第j个滞后根β_j取值的上下限，

结合式(26)和(28)，应用式(23)的优化算法构造形式，可得外层优化 迭代的无约束化目标函数

式中，f如式(26)，

为外层无约束化目标函数中第i个滞后根上下限约束的 拉氏乘子；各ν_k(k＝i或j)为将不等式约束转换为等式约束的松弛变量，优化迭 代运算中可运用目标函数关于ν_k的导数为零条件而被消除，故不作为新变量。

步骤6：两级交替双重迭代优化算法。本发明提出的非线性两级交替双重迭 代优化算法流程框图见图1所示。其外层优化迭代计算采用标准数学方法的拟 牛顿方向与一维线性搜索的成熟优化迭代方法，其中的关键步骤在于总误差函 数关于各气动力滞后根β_i的一阶导数值计算问题。本发明提出了针对复数形式 的有理近似函数一阶导解析计算方法，算法如下：

当给定非线性参数β₁,β₂,…,β_N-2时，由步骤3和4可知，各线性系数矩阵由 两次线性方程组的求解唯一确定，这也意味着线性系数是依赖于非线性系数的 非独立变量。因此，描述总误差函数的式(26)可形式地表达为

f＝f(β₁,β₂,...,β_N-2) (30)

实际上，式(26)因隐含复杂复数运算而难于实施导数计算。为此，本发 明变化其求和方式，将式(30)式等价地改写为两项纯实数计算方式，即：

其中，w_ij为加权系数，可仅为1；上标DLS为归一化有理近似拟合所得的各矩 阵元素与二次误差函数有理近似拟合所得的各矩阵元素相加算子。

式(31)的最后一个等式说明总误差函数是滞后根的函数。为实施式(29) 的优化计算，需计算式(31)关于滞后根的导数。由导数算子的线性性质及链 式法则，本发明给出关于式(31)式的滞后根导数算式为：

其中，

为方便起见，将上式简记为向量形式：

其中，符号[·]^T表示行向量的转置，即为列向量。

按本发明上述步骤建立的无约束化目标函数及其导数计算公式，应用标准 数学方法即可完成总误差函数的极小化运算。关于本发明的完整计算流程见图1 所示，下面将针对具体机翼结构，结合附图和实施例对本发明的实施步骤及技 术效果做进一步说明。

基于上述技术方案，本发明的有益效果是：本发明提高了飞机结构气动伺 服弹性系统建模分析中非定常气动载荷历程的预计精度，本发明采用了一种高 精度的两类变量两级交替双重优化迭代算法，该方法基于离散频率下的广义非 定常气动力系数矩阵数据，在关心的频段范围内，对多点频率下的广义非定常 气动力系数矩阵元素实施固定有理近似格式的两级交替双重优化迭代拟合计 算，可极大降低数值拟合误差，从而为飞机气动伺服弹性系统的动响应、动稳 定性以及主动弹性控制提供高精度载荷历程数值。

实施例

如图1至8所示，本实施例是对一大展弦比机翼薄壁结构的广义非定常气 动力系数矩阵进行处理的过程，结合附图和实施例对本发明进一步说明。

步骤Ⅰ，从CATIA软件中取出机翼结构的CAD数字模型数据。无人机几 何外形及内部骨架结构如图2所示，可划分为中央翼、左右外翼，分割点为A3L 及A3R，图中标注了关键点的位置。该模型关于中面A0-B0对称，表1给出了 右半机翼尺寸参数，其半展长为8900mm。

表1机翼平面尺寸参数

位置	长度(mm)	位置	长度(mm)
				A0-B0	937.00	B31R-B32R	272.00
A1R-B12R	2141.00	A4R-B41R	427.00
				A2R-B2R	1055.00	B41R-B42R	227.00
A3R-B31R	566.00	A5R-B5R	440.00

步骤Ⅱ，采用Hypermesh软件对步骤Ⅰ的机翼CAD模型进行有限元网格剖 分，如图3所示。注意，本算例采用了较小尺度的板壳和梁两类单元的有限元 网格剖分方法。这类结构剖分方法对单元剖分的尺寸没有严苛要求，密度可大 些；结构质量采用集中质量方法建模，即将每个结构单元的分布质量均分到有 限元模型的节点上，但注意需将一个结构框格内的质量要均摊到有梁缘条、翼 肋支持的结构节点上，以获得高质量的结构整体位移振型解。

步骤Ⅲ，完成结构有限元模型参数的输入，包括结构材料及质量密度，如 表2所示。将完整的机翼结构模型运用Nastran软件进行结构固有振型计算。图 4为该机翼模型前25阶弹性模态中的9个对称变形模态云图。除此弹性模态外， 利用刚体模态方法引入了沉浮、俯仰两个刚体模态，故后续的开环系统降阶方 法共选用结构的前11阶固有位移振型模态。表3为9个弹性模态对应固有频率 f及其对应的减缩频率k值。减缩频率k的计算式为：

其中，机翼参考弦长REFC为924mm，马赫数Ma为0.6，假定飞行高度为1000m， 当地音速为V_a＝336.4m/s。

表2模型材料性能参数

表3机翼固有频率及相应的减缩频率

编号	f(Hz)	k	编号	f(Hz)	k	编号	f(Hz)	k
									1	3.11E+00	0.044203	4	4.19E+01	0.596063	7	9.02E+01	1.282887
2	3.11E+00	0.044203	5	4.86E+01	0.692093	8	1.03E+02	1.47131
									3	2.54E+01	0.361366	6	7.07E+01	1.005685	9	1.14E+02	1.618

步骤Ⅳ，在机翼弦平面上剖分气动力系数计算的面元网格，如图5所示。 各区域网格划分形式如表4所示，网格均为顺气流均匀划分。表4中区域1～6 见图2所注，左右机翼各划分为六个区域，其区域及面元划分均关于中面对称。

表4机翼弦平面各分区面元划分个数

机翼模型所取最高阶模态对应的减缩频率值为1.618，故将最大减缩频率 k_max值选为1.65。为获取广义非定常气动力系数矩阵的拟合公式，利用ZAERO 计算出0.6Ma下减缩频率在0.0～1.65范围内的23个非定常气动力系数矩阵 [AIC]，所取离散减缩频率k值分别为0.0、0.0001、0.05、0.1、0.15、0.20、0.25、 0.3、0.4、0.45、0.5、0.6、0.7、0.8、0.9、0.95、1.1、1.2、1.3、1.4、1.5、1.6、 1.65。

步骤Ⅴ，应用三维样条函数插值方法，将机翼结构的11个位移振型插值到 机翼弦平面的面元网格上，第7和8阶的结构弹性变形插值结果如图6所示， 其他未示出。

步骤Ⅵ，由ZAERO软件的ZONA6模块计算0.6Ma下23个减缩频率k对 应的非定常气动力系数矩阵[AIC]。结合结构有限元模型参数，利用式(10)和 (11)求解机体模态和控制面模态对应的广义非定常气动力系数矩阵Q_qq、Q_qδ。 本算例模型中含有11个结构模态，即2个刚体和9个弹性模态；另有2个控制 面模态。面元模型中左、右机翼控制面分别建立，此处两个控制面模态分别对 应左、右翼面的操纵面面模态。将对应一个减缩频率下的广义非定常气动力系 数矩阵组合成[Q_qq(jk)|Q_qδ(jk)]，为11×13的复数矩阵。

步骤Ⅶ，选择本算例的拟合约束条件为：减缩频率k＝0处，约束拟合值等 于原始计算结果；取减缩频率k＝0.0001处，约束拟合值的虚部等于实际计算值 的虚部。按ZAERO软件手册给出的如下经验公式计算滞后根，表5给出了N取 1～6的值。按最小二乘法计算式(23)得线性方程组解。

表5 ZAERO软件计算不同N值的滞后根

N	β<sub>1</sub>	β<sub>2</sub>	β<sub>3</sub>	β<sub>4</sub>	β<sub>5</sub>	β<sub>6</sub>
							1	0.70125
2	0.311667	1.246667
							3	0.175313	0.70125	1.577813
4	0.1122	0.4488	1.0098	1.7952
							5	0.077917	0.311667	0.70125	1.246667	1.947917
6	0.057245	0.22898	0.515204	0.915918	1.431122	2.060816

步骤Ⅷ，按附图1所示的算法流程设计两级交替的双重迭代优化流程，并 将步骤Ⅶ计算所得的有理近似拟合矩阵作为初始数据，取收敛精度ε＝10^-4，权重 系数w_ij＝1；同样取N为1～6，最终计算得到的滞后根见表6。步骤Ⅶ与步骤Ⅷ 计算所得的总误差结果对比如表7所示。

表6本发明非线性方法计算不同N值的滞后根

N	β<sub>1</sub>	β<sub>2</sub>	β<sub>3</sub>	β<sub>4</sub>	β<sub>5</sub>	β<sub>6</sub>
							1	0.450792
2	0.290515	1.467621
							3	0.136502	0.967423	1.018347
4	0.122447	1.226503	1.529251	1.693386
							5	0.119598	1.23527	1.44566	1.503459	1.5668133
6	0.11648	1.012977	1.187568	1.238286	1.2931153	1.3499809

表7线性优化与本发明非线性优化的误差值比较

非线性系数个数	线性优化误差	非线性优化误差
			1	8.690468	6.134711
2	5.222127	3.275647
			3	4.59832	1.651599
4	4.00438	0.9940788
			5	3.67175	0.6367089
6	2.696798	0.4750968

从步骤Ⅵ到步骤Ⅷ，可完成相应各离散减缩频率值下广义非定常气动力系 数的气动物理计算，以及线性最小二乘有理近似拟合计算和本发明的非线性两 级交替双重有理近似优化计算。对比表5和表6可以看出，本发明算法所得的 优化滞后根分布特性明显较线性拟合计算的更为集中；且表7说明优化的滞后 根明显改进了广义非定常气动力系数的有理近似拟合效果，且采用3个滞后根 优化的计算总误差已小于固定6个滞后根的线性计算总误差。图7给出了三种 计算所得的[Q_qq，Q_qδ]矩阵中第一行各元素的复数值，对应各模态的广义非定常 升力系数；图8给出了三种计算所得的[Q_qq，Q_qδ]矩阵中第二行各元素的复数值， 对应各模态的广义非定常俯仰力矩系数。各图中曲线REAL代表直接利用 ZAERO软件计算出的离散点对(Ma,k)的理论值，APP01、APP04及APP06表 示滞后根个数分别为1、4和6时的线性最小二乘算法近似结果；而APP04-1表 示本发明方法取N等于4的非线性拟合优化结果。从图中可以得到三点结论：

1.随着滞后根的增多，线性拟合效果会变好。但对于原始数据变化比较复 杂的元素，如弹性模态4、5、7、8、9，即使滞后根个数为6时，拟合效果仍不 理想。

2.对比APP04、APP04-1可知，取相同滞后根的个数，对滞后根实施非线 性优化的结果明显优于固定滞后根的线性最小二乘拟合效果，体现了滞后根对 拟合结果的敏感性；即使原始数据变化比较复杂的元素，本发明算法总能使拟 合结果逼近于原始理论值，肯定了算法对高度非线性的拟合能力。

3.对比APP04-1、APP06可知，本发明算法下4个优化滞后根的拟合结果 已优于6个固定滞后根的线性拟合结果，也即本发明算法将高拟合精度、低阶 气动力模型两个目标集于一体。

由此可见，本发明算法虽其复杂度高于线性最小二乘拟合算法，但对数值 拟合精度的提高非常明显，且算法在现代常规个人电脑上的执行并无难度。可 为现代固定翼飞机结构设计迭代中实施弹性动响应、动稳定分析以及弹性主动 控制模拟评估提供精度可靠的有效工程应用方法。

Claims (6)

1.气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：包括如下步骤：

1)通过已有商用工程数值计算软件获取分析气动及结构的原始数据：构建航空器结构动力学分析的有限元网格模型，计算并选取一定数量的结构固有振型模态矩阵；构建气动布局的气动载荷面元网格，确定选取的多个离散减缩频率值，计算获取相应各减缩频率下的非定常气动力影响系数矩阵，结合两类原始数据计算可得相应各离散减缩频率的非定常广义气动力系数矩阵；

2)预设初始参数：初始滞后根、初始Hessian矩阵和收敛精度控制值，将初始参数以及步骤1)中的各离散减缩频率下的广义气动力系数矩阵数据送入近似有理子程序进行计算，返回二次误差函数的有理近似拟合数据；

3)计算广义气动力系数矩阵的有理近似拟合数据，结合步骤1)中的各离散减缩频率下广义气动力系数矩阵数据和步骤2)中二次误差函数的有理近似拟合数据，计算总误差加权拟合系数，并结合滞后根的上下限约束，形成无约束化目标函数，计算无约束化目标函数的函数值和导数值；

4)按数学标准方法计算拟牛顿优化搜索方向，并沿搜素方向进行一维线性搜索计算步长因子；

5)计算新的滞后根，并将新的滞后根以及各离散减缩频率下广义气动力系数的矩阵数据一起送入有理近似拟合子程序中，调用所述子程序计算得到二次误差函数的有理近似拟合数据；

6)根据步骤5)中二次误差函数的有理近似拟合数据计算新的广义气动力系数矩阵系数，并计算总误差函数；

7)将步骤6)中迭代的总误差函数与前一次迭代的总误差函数相减，取绝对值；判断绝对值是否小于预设的收敛精度控制值，是，则停止优化迭代；否则，用BFGS标准数学公式构造新的Hessian矩阵,令迭代计数加1,转至步骤3)继续迭代优化。

2.根据权利要求1所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：步骤2)中二次误差函数的有理近似拟合形式为：

其中，jk是纯虚数，k为实数域内的减缩频率变量，

为取二次误差函数为有理近似形式的复数拟合系数矩阵，β_j(j＝1,…,N-2)称为滞后根，m为含滞后根分式项的下标，工程中N的取值一般不大于6。

3.根据权利要求2所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：所述二次误差函数为：

其中，Q(jk)为工程软件计算提供的原始广义气动力系数矩阵；

是归一化有理近似平方误差和的无约束化目标函数优化计算所得结果。

4.根据权利要求3所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：所述归一化有理近似平方误差和为:

其中，下标ij与广义气动力系数矩阵的第i行及第j列元素相对应；上标T表示矩阵转置；

为归一化系数，

为原始广义气动力系数矩阵第i行及第j列元素的实部与虚部；

上式为广义气动力系数矩阵各元素的有理近似拟合的矩阵及向量乘积形式；D_m＝β_m ²+k_i ²(m＝1,…，N-2)，当行号小于l时，i取矩阵行号，当行号大于l时，i取行号减l，l为用户选定的离散减缩频率数；右端列向量为有理近似各矩阵中的第i行第j列元素，为待定系数；

{C}_ij为用户提高有理近似的针对性所选择的线性约束的组合系数，下标ij的含义同前，C_0ij为用户约束的常数项；{λ}_ij为拉氏乘子列向量。

5.根据权利要求1所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：步骤3)中所述外层迭代是将所有减缩频率下各矩阵元素的拟合误差加权总和作为总目标函数，如下式：

其中，n为广义气动力系数矩阵的列数，m为广义气动力系数矩阵的列数，l为选定的离散减缩频率数，w_ij为对矩阵第i行第j列元素的加权系数，可取1；F_ij(k_p)，G_ij(k_p)分别为减缩频率为离散值k_p时原始广义气动力系数矩阵各复数元素的实部与虚部，

项为对二次误差函数有理近似形式实施线性优化后所得的广义气动力系数矩阵各复数元素的实部与虚部,

外层非线性滞后根优化的约束条件设置为各滞后根的上下限：

0＜L_j＜β_j＜U_j，(j＝1，…，N-2)

其中，L_j和U_j分别为第j个滞后根β_j取值的上下限,

结合拉格朗日乘子法，含所述滞后根上下限约束的归一化有理近似平方误差和的无约束化目标函数作为外层优化迭代的目标函数：

其中，f如前述，λ_i ^OL为外层无约束化目标函数中第i个滞后根上下限约束的拉氏乘子；各ν_k(k＝i或j)为将不等式约束转换为等式约束的松弛变量，在优化迭代运算中可运用目标函数关于ν_k的导数为零条件而被消除，故不作为新变量。

6.根据权利要求5所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：所述总目标函数可写作两项纯实数的计算方式：

其中，w_ij为加权系数，可仅选1；上标DLS为有理近似拟合所得各矩阵元素与二次误差函数有理近似拟合所得各矩阵元素的相加算子，

由该函数可得关于滞后根的偏导数算式为：

其中，

将上式简记为向量形式：

其中，符号[·]^T表示行向量的转置，即为列向量。

Images