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

摘要

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

权利要求

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(jkz)为工程软件计算提供的相应离散减缩频率的原始广义气动力系数矩阵；

是归一化有理近似平方误差和的无约束化目标函数优化计算在相应离散减缩频率kz下的结果。

4.根据权利要求3所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：所述归一化有理近似平方误差和为:其中，下标ij与广义气动力系数矩阵的第i行及第j列元素相对应；上标T表示矩阵转置； 为归一化系数， 为原始广义气动力系数矩阵第i行及第j列元素的实部与虚部；

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

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

5.根据权利要求1所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：步骤3)中所述总误差加权拟合系数作为总目标函数，如下式：其中，n为广义气动力系数矩阵的列数，m为广义气动力系数矩阵的列数，l为选定的离散减缩频率数，wij为对矩阵第i行第j列元素的加权系数，可取1；Fij(kp)，Gij(kp)分别为减缩频率为离散值kp时原始广义气动力系数矩阵各复数元素的实部与虚部，项为对二次误差函数有理近似形式实施线性优化后所得的广义气动力系数矩阵各复数元素的实部与虚部；

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

0＜Lj＜βj＜Uj，(j＝1，…，N-2)

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

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

各νk(k＝i或j)为将不等式约束转换为等式约束的松弛变量，在优化迭代运算中可运用目标函数关于νk的导数为零条件而被消除，故不作为新变量。

6.根据权利要求5所述的气动伺服弹性系统非定常气动载荷的有理近似优化方法，其特征在于：所述总目标函数可写作两项纯实数的计算方式：其中，wij为加权系数，可仅选1；上标DLS为有理近似拟合所得各矩阵元素与二次误差函数有理近似拟合所得各矩阵元素的相加算子，由该函数可得关于滞后根的偏导数算式为：

其中，

将上式简记为向量形式：

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

说明书

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

技术领域

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

背景技术

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

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

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

发明内容

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

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

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

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

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

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

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

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

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

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

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

[0016]

[0017] 其中，jk是纯虚数，k为实数域内的减缩频率变量， 为取二次误差函数为有理近似形式的复数拟合系数矩阵，βj(j＝1,…,N-2)称为滞后根，m为含滞后根分式项的下标，工程中N的取值一般不大于6。

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

[0019]

[0020] 其中，Q(jk)为工程软件计算提供的原始广义气动力系数矩阵； 是归一化有理近似平方误差和的无约束化目标函数优化计算所得结果。

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

[0022]

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

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

[0025]

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

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

[0028] 进一步地：步骤3)中所述总误差加权拟合系数作为总目标函数，如下式：其中，n为广义 气动 力系数 矩阵的 列 数 ，m 为广 义气 动力 系数 矩阵的 列数 ，l 为

[0029] 选定的离散减缩频率数，wij为对矩阵第i行第j列元素的加权系数，可取1；Fij(kp)，Gij(kp)分别为相应各减缩频率离散值kp时所得原始广义气动力系数矩阵各复数元素的实部与虚部， 项为对二次误差函数有理近似形式实施线性优化后所得的广义气动力系数矩阵各复数元素的实部与虚部 ；

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

[0031] 0＜Lj＜βj＜Uj，(j＝1，…，N-2)

[0032] 其中，Lj和Uj分别为第j个滞后根βj取值的上下限,

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

[0034]

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

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

[0037]

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

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

[0040]

[0041] 其中，

[0042]

[0043]

[0044] 将上式简记为向量形式：

[0045]

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

[0047] 本发明的有益效果为：

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

附图说明

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

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

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

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

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

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

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

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

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

[0058] 其中：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六阶对称弯曲模态。

具体实施方式

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

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

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

[0062]

[0063]

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

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

[0066] 式(1)中的Mqq、Cqq及Kqq阵当采用模态质量归一化的结构振型模态时，前m阶结构振型模态阵为

[0067]

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

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

[0070]

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

[0072] 假定结构模型的一致质量矩阵为Ms，广义质量阵Mqq和Mqδ的表达式如下Mqq＝FqTMsFq＝In×n Mqδ＝FqTMsFδ   (5)

[0073] 同理，假定结构的总刚度矩阵为Ks，结构的广义刚度矩阵Kqq为

[0074]

[0075] 其中，diag表示对角矩阵，即非对角元素均为零，ωi为第i阶正交模态的圆频率，ωi2为系统特征方程的特征值，刚体模态的对应值为零。广义阻尼阵Cqq为下述形式[0076] Cqq＝diag(2ξ1ω1,2ξ2ω2,…,2ξmωm)   (7)其中，ξi为第i阶固有振型的模态阻尼系数。

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

[0078]

[0079]

[0080] 上两式中的GH和Gδ分别为机体结构及操纵面表面有限元网格与其弦平面面元网格间的几何插值转换矩阵； 为沿来流方向的导数算子，i为单位虚数，k＝ωb/U∞称为减缩频率，ω为气流受扰的圆频率，b为升力面半弦长，U∞为来流速度。

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

[0082]

[0083]

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

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

[0086] Q＝F(k)+jG(k)   (12)

[0087] 其中，j为单位虚数，即j2＝-1。

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

[0089]

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

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

[0092]

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

[0094]

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

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

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

[0098]

[0099] 式中，Dz＝β2z+k2，z＝1,2,…,N-2；{x}＝{a0ij,a1ij,…,a(N-2)ij}，Fij(k)和Gij(k)分别为式(12)广义气动力系数矩阵第i行第j列复数元素的实部与虚部。

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

[0101]

[0102] 式中，k2通常取接近于零但大于零的某个小数值。

[0103] 3.约束某个拟合系数为零(如拟合系数a2ij为零),即

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

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

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

[0107]

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

[0109]

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

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

[0112]

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

[0114]

[0115] 式中， 其他符号见式(19)～(20)。

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

[0117]

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

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

[0120]

[0121] 其中，Q(jkz)为工程软件计算提供的相应离散减缩频率的原始广义气动力系数矩阵； 为由式(23)拟合优化所得结果。

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

[0123]

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

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

[0126]

[0127] 式中，n为广义气动力系数矩阵的列数，m为广义气动力系数矩阵的列数，l为选定的离散减缩频率数，wij为对矩阵第i行第j列元素的加权系数，可取1；Fij(kp)，Gij(kp)分别为减缩频率为离散值kp时原始广义气动力系数矩阵各复数元素的实部与虚部，项为对二次误差函数有理近似形式实施线性优化后所得的广义
气动力系数矩阵各复数元素的实部与虚部。为后续描述方便，称式(26)为广义气动力拟合系数的总误差函数。

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

[0129]

[0130] 其中， m＝3，…，N。

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

[0132] 0

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

[0134]

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

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

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

[0138] f＝f(β1,β2,...,βN-2)   (30)

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

[0140]

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

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

[0143]

[0144] 其中，

[0145]

[0146]

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

[0148]

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

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

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

[0152] 实施例

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

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

[0155] 表1机翼平面尺寸参数

[0156] 位置 长度(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

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

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

[0159]

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

[0161] 表2模型材料性能参数

[0162]

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

[0164] 编号 f(Hz) k 编号 f(Hz) k 编号 f(Hz) k1 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

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

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

[0167]

[0168]

[0169] 机翼模型所取最高阶模态对应的减缩频率值为1.618，故将最大减缩频率kmax值选为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。

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

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

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

[0173]

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

[0175]

[0176]

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

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

[0179]N β1 β2 β3 β4 β5 β6
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

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

[0181] 非线性系数个数 线性优化误差 非线性优化误差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

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

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

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

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

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