尖扭转前的侧视图，化是翼尖 扭转后的侧视图；
[0034] 图7是不同方法的升力曲线对比图；
[0035] 图8是优化后的非均匀面元网格分布图；
[0036] 图9是机翼最大位移和总升力变化曲线，其中，9a是机翼最大位移变化曲线，9b是 机翼总升力变化曲线；
[0037] 图10是本发明的流程图。图中：
[0038] 1.迎角；2.升力/N;3.C抑曲线；4.未修正面元法曲线点修正面元法曲线；6.迭 代次数；7.结构最大位移/mm;8.升力/N;9.修正面元法曲线；10.C抑曲线；11.前梁腹板； 12.后梁腹板；13.翼肋。
【具体实施方式】
[0039] 本实施例是NACA0012对称翼型的M6机翼的面元修正与网格预先自适应计算方 法。
[0040] 所述NACA0012对称翼型的M6机翼根部完全固支。机翼参数如下：翼根弦长为 0.8139m，翼尖弦长为0. 4573m，展长为1. 1963m，前缘后掠角为30。，后缘后掠角为15. 8。， 无扭转角，蒙皮厚度为0. 002m，梁腹板厚度为0. 0015m，翼肋厚度为0. 0015m，材料为侣合 金，弹性模量E= 70GPa，泊松比y= 0. 30,密度为P= 2. 7X103kg/m3。图3为机翼的整 体有限元模型，其中，（a)代表从机翼翼根到翼尖的展向方向，化）代表从机翼前缘到后缘 的弦向方向。图4为机翼前梁腹板11、后梁腹板12和翼肋13的有限元模型。采用S角形 和四边形壳单元进行有限元建模，机翼由翼肋、梁腹板和蒙皮组成。11个翼肋沿展向均匀 分布。梁腹板包括前梁腹板和后梁腹板，分别位于距离机翼前缘27.0%和63.5%处。弦 向各网格节点位置分别为距离前缘0. 0 %、5. 0 %、10. 0 %、27. 0 %、45. 5 %、63. 5 %、82. 0和 100%处，展向各网格节点位置与翼肋位置一致。
[0041] 本实施例的具体实施过程是：
[0042] 步骤1，计算不同迎角下机翼部件原刚体气动载荷数据。
[0043] 通过CFD软件计算不同迎角下机翼部件原刚体的气动载荷数据。将机翼部件作为 刚性几何体划分CFD计算的空间网格，所述的空间网格采用非结构网格。翼面采用=角形 单元划分，最小尺度为1mm。附面层用=棱柱网格划分，共20层，附面层之外的空间网格采 用四面体单元。从附面层最内层开始网格尺度W1. 2倍的体积比率增长。CFD数据由FLUENT 计算，采用禪合压力基求解器。空气假设为理想气体，溫度为300K，参考压力为海平面大气 压，端流模型采用S-A模型，其余参数选择默认设置，机翼的CFD表面网格见图5。计算亚音 速时机翼迎角在0°到20°之间的C抑数据，迎角间隔为r，提取机翼物面网格节点上的 压力分布数据，得到机翼物面网格节点上压力分布的数据库。本实施例中，马赫数取0. 6,机 翼迎角范围为0°到10。。
[0044] 步骤2,计算变形刚体气动载荷数据。
[0045] 在机翼升力曲线的非线性段选择一个迎角aW，将步骤1中计算得到的该迎角时 的气动载荷加载到机翼上，利用NASTRAN软件计算得到翼尖的扭转角。用CAD软件重新几 何造形W使翼尖扭转相同的角度且使扭转角从翼根到翼尖线性变化。再按照步骤1中的方 式重新划分CFD计算网格并实施常规CFD数值计算，提取机翼变形后各翼面网格节点上的 压力分布，得到机翼的变形刚体气动载荷数据。
[0046] 本实施例中，取aW为8°，计算得到的翼尖扭转角为向下扭转4°，翼尖扭转前 后的CFD网格如图6所示，其中，6a是机翼翼尖扭转前的侧视图，化是翼尖扭转后的侧视 图。
[0047] 步骤3,对升力面的弦平面进行面元网格参数初始化。
[0048] 上述的面元网格参数即指网格加密方法的参数。通用的代数辅助网格加密方法有 =次多项式插值法、等差数列法、等比数列法、拉伸函数法。本实施例采用能够实现多种加 密方式的=次多项式插值法，需要两条=次曲线分别对机翼展向和弦向进行网格点布置， 分别记为展向S次曲线和弦向S次曲线。
[0049] 给定展向和弦向的网格点数量范围和每条S次曲线两端的斜率值作为约束条件， 需初始化的变量为弦向和展向的网格点数量，即待优化变量。
[0050]S次曲线能够实现在需要的位置处进行网格加密。图2给出了（a)、化）和（C)S 条典型的=次曲线，图中对横坐标和纵坐标均进行了归一化处理；所述归一化处理指将横 坐标最大值和纵坐标最大值均缩放为1。因确定一条=次曲线需要四个约束，其中的两个约 束为对横坐标和纵坐标进行归一化后曲线两端的位置化0)和（1，1)，另外两个约束为所 述给定的=次曲线两端的斜率；通过所述四个约束能够完全确定=次曲线的形状。本实施 例中，如图2所示，S条给定的S次曲线（a)、化）、（c)在（0,0，）点处的斜率分别为2. 6、1、 0. 1，在（1，1)点的斜率分别为0. 1、1、2. 6。
[0051] 步骤4,用几何算法对机翼弦平面进行面元网格划分。
[0052] 根据给定的=次曲线参数及展向和弦向网格点数量对机翼弦平面进行面元网格 节点布置，即划分机翼弦平面的面元网格。
[0053] 在对展向布置网格节点时，先对展向S次曲线的横轴按照展向节点数量进行等距 划分，此时横轴上各点在S次曲线上对应的纵坐标即表示此节点的展向比例位置，且横轴 正方向代表从机翼翼根到翼尖的方向；在对弦向布置网格节点时，先对弦向=次曲线的横 轴按照弦向节点数量进行等距划分，此时横轴上各点在=次曲线上对应的纵坐标即表示此 节点的弦向比例位置，且横轴正方向代表从机翼前缘到后缘的方向。
[0054] 对于弦向S次曲线，若需要在机翼前缘加密，则应使曲线在化0)点的斜率值减 小；对于展向S次曲线，若需要在翼尖加密，则应使曲线在（1，1)点的斜率值减小；若要采 用均匀划分形式，则只需要将（〇,〇)、（1，1)处的斜率值均设为1即可。
[00巧]步骤5,面元法计算气动力矩阵。
[0056] 通过ZAER0软件生成气动力矩阵。具体是，将步骤4得到的网格分布输入到ZAER0 软件，计算给定飞行状态下的气动力影响系数矩阵[AIC]。上述飞行状态与步骤1中计算 CFD数据的飞行状态一致。本实施例中，飞行状态为0. 6马赫数，海平面大气压，溫度300K。
[0057] 步骤6,采用分段斜率修正面元法对气动载荷分布进行修正，获得修正后的载荷分 布。
[0058] 利用步骤1得到的定常刚体CFD计算数据对面元法的气动力影响系数矩阵[AIC] 实施精细修正。
[0059] 修正前的面元法求解气动力公式为
[0060] {CJ= [AIC] {W} (1) 其中，[AI幻为空气动力影响系数矩阵，{W}为无量纲下洗列阵，{Cp}为各面元上的压 差系数列阵。
[0061] 面元修正法基本思路是将面元法计算的气动压差系数与C抑计算数据进行对比 修偏，获得气动力影响系数矩阵的修正因子矩阵[WT]，并将修正因子矩阵前乘于气动力影 响系数矩阵[AIC]，此方法称为斜率修正法，如下式。
[0062] [AIC]*=[WTJ[AIC]似
[0063] 式中，[WTi]为修正因子矩阵。[AlCr为修正气动影响系数矩阵。
[0064] 斜率修正法将定常压差系数随迎角的变化率作为计算依据来计算修正因子矩阵 [WTi]，见公式（3)、（4)。公式（4)中[WTi] 11代表CFD定常气动力数据随迎角的变化率与面 元法定常气动力数据随迎角的变化率的比值。对于进入非线性段的气动力系数修正，显然 相邻迎角越相近，修正效果越好，即高精度面元修正。
[006引{Cp} = [WTJ [AIC] {W} (3)
[0066]
(4)
[0067] 其中，[CpJigu。。表示第i个面元网格上C抑压差系数关于迎角的差分数值；[Cp。] 。。。。1为第i个网格面元法计算得出的压差系数关于迎角的差分数值；[WT1] 11表示割线修正 因子矩阵的对角元素第i项，每个面元对应的修正因子构成了对角修正矩阵[WTi] ;a1、a2 表示两相邻的不同飞行迎角。
[0068]由式（4)，似给定两个不同迎角的C抑气动力系数即能够得到精确的修正结果，实 际运没有意义，本实施例的目的是将此修正因子应用于任意飞行迎角下的弹性变形条件， 运可理解为各剖面的迎角不同。于是，用尽可能多个不同迎角下的C抑数据进行修正，所得 修正因子无疑是高精度的，但修正因子将会成为随迎角变化的非线性关系。为此，本实施例 采用公式（5)进行分段线化高精度面元修正：
[006引
（5)
[0070] 式中，n为迎角范围划分的子区间数目，[WTji为第i个迎角区间采用斜率修正法 得到的修正因子矩阵，[Cor]i为各面元在第i个迎角区间内的下洗速度值占总下洗速度的 比例系数对角矩阵。
[0071] 步骤7,根据式（6)计算CFD数据与修正面元法结果的误差范数，判断是否达到收 敛条件。若是，则得到最优网格分布，迭代过程结束；反之采用自适应模拟退火算法更新网 格分布参数，更新后转至步骤4,重复步骤4~步骤6;在重复步骤4~步骤6的过程中不断 获得新的修正后的载荷分布。对得到的新的修正后的载荷分布通过本步骤计算C抑数据与 修正面元法结果的误差范数，判断是否达到收敛条件