，直至满足收敛条件。
[0072] 本实施例中，收敛条件为连续20步最优解不变。
[0073]面元网格自适应优化的目标函数为式化），即C抑数据与修正面元数据的误差。
[0074]
(6 )
[00巧]其中，N是计算C抑数据时的迎角数量，本实施例中为11，L、MX、MY分别表示机翼 的气动力、俯仰力矩W及滚转力矩。所述1-MX、MY的下标CFD均表示CFD计算结果，由CFD 计算得到的压力系数在翼面积分得到；所述LMX、MY的下标panel表示面元修正法的计算 结果，由步骤6中的压差系数Cp在弦平面上积分得到。
[0076] 本实施例采用自适应模拟退化算法来优化网格分布，使误差范数最小化。自适应 模拟退化算法的基本过程为：把每种组合状态S看成某一物质系统的微观状态，而将其对 应的目标函数C(s)看成该物质系统在状态S下的内能；用控制参数T类比溫度，让T从 一个足够高的值慢慢下降，对每个T用Metropolis抽样法在计算机上模拟该体系的热平 衡态，即对当前状态S作随机扰动W产生一个新状态s'，如果C(s'）<C(s)则接受s'为 下一状态，否则W概率eWs'WsWT接受。经过一定次数（Markov链长）的捜索，认为系 统在此溫度T下达到平衡，则降低溫度T再进行捜索，直到满足结束条件（"RA0Jian.A StatisticalPhysicsStudyonSpatialInhomogeneousSystemswithSmallNumberof Particles巧巧unan化iversity, 2013.")。在本实施例中，组合状态s就是展向网格数量 和弦向网格数量组成的向量，C(s)为式（6)定义的目标函数值。
[0077] 图7中给出了分别采用面元法、斜率修正面元法及C抑=种方法计算得到的刚性 M6机翼在0. 6马赫数下不同迎角的机翼总升力。其中C抑数据为步骤1中得到的原刚体气 动载荷数据。
[007引由图7可W看出，迎角从0°变化到10°时，M6机翼气动力基本呈线性增长。未 修正面元法曲线4和C抑曲线3的误差随着迎角增大逐渐变大，运是由于该方法不能考虑 机翼附面层、粘性等影响导致的。由修正面元法曲线5可W看出，经过斜率修正面元修正法 计算的气动升力与CFD结果吻合良好，能够准确反映气动力随迎角变化的非线性特征。运 验证了分段线化面元修正方法的有效性和精度。
[0079] 网格优化时，分别取机翼展向与弦向网格点数范围为巧，1引及巧，2引两种情况， 并采用均匀网格划分及非均匀网格划分两种方法。对于非均匀的网格，曲线起点和终点斜 率分别设为0.4和0.5。优化结果如表1及表2所示。表中优化结果第一列为展向网格点 数量，第二列为弦向网格点数量。表中的A1表示总升力相对误差，A2表示总滚转力矩相 对误差，A3表示总俯仰力矩的相对误差。
[0080] 由表中结果可W看出，通过面元网格的自适应优化，目标函数最小值对应的并非 最细网格情形。通过均匀网格划分与加密网格的比较可知，适当加密机翼前缘及翼尖网格 可W较大地提高俯仰力矩及滚转力矩的计算精度，并且采用非均匀网格划分所需网格节点 更少，运可W提高气动力计算效率。对于气动弹性优化等需要反复气动力计算的领域来讲， 其效果非常可观。
[0081] 表1均匀网格划分结果
[0082]
[0083] 表2非均匀网格划分结果
[0084]
[0085]为验证本发明应用于静气动弹性分析时的效率和精度，在此对机翼在4°迎角下的 几何变形进行分析比较，面元网格分布采用表2的第二个结果，对应图8。分别用CFD和高 精度面元修正法与结构有限元程序计算禪合，反复迭代计算结构变形与气动力。机翼的气 动力总和W及结构最大变形随迭代步变化曲线如图9所示。由修正面元法曲线9和CFD曲 线10可W看出，C抑与高精度面元修正法均迭代4步收敛，且收敛结果一致。弹性气动力 较刚性机翼损失5. 6%。计算时间方面C抑迭代共耗时4小时，修正面元法仅耗时2分钟。 高精度面元修正法准确计算了弹性变形后的机翼气动力，与CFD结果几乎一致，运说明优 化后的网格结合高精度面元修正法在计算弹性变形后机翼气动力方面是准确可靠的。
【主权项】
1. 一种面元修正与网格预先自适应计算方法，其特征在于， 步骤1，计算不同迎角下机翼部件原刚体气动载荷数据； 将机翼部件作为刚性几何体划分CFD计算的空间网格；所述的空间网格采用非结构网 格；计算亚音速时机翼迎角在0°到20°之间的CFD数据，迎角间隔为Γ ; 提取机翼物面网格节点上的压力分布数据，得到各机翼物面网格节点上压力分布，从 而确定不同迎角下机翼部件原刚体气动载荷数据； 步骤2,计算变形刚体气动载荷数据； 在机翼升力曲线的非线性段选择一个迎角a raf，将步骤1中计算得到的该迎角时的气 动载荷加载到机翼上，利用NASTRAN软件计算得到翼尖的扭转角；用CAD软件重新几何造形 以使翼尖扭转相同的角度且使扭转角从翼根到翼尖线性变化；再按照步骤1中的方式重新 划分CFD计算网格并实施常规CFD数值计算，提取机翼变形后各翼面网格节点上的压力分 布，得到机翼的变形刚体气动载荷数据； 步骤3,对升力面的弦平面进行面元网格参数初始化； 采用通用的三次多项式插值法，通过展向三次曲线和弦向三次曲线分别对机翼展向和 弦向进行网格点布置； 对所述展向三次曲线和弦向三次曲线分别进行归一化；以给定的展向网格点数量范 围、弦向网格点数量范围和每条三次曲线两端的斜率值作为约束条件，需初始化的变量为 弦向和展向的网格点数量； 步骤4,用几何算法对机翼弦平面进行面元网格分划； 根据给定的三次曲线参数及展向网格点数量和弦向网格点数量对机翼弦平面进行面 元网格节点布置； 步骤5,面元法计算气动力矩阵； 通过ZAERO软件生成气动力矩阵；具体是，将步骤4得到的网格分布输入到ZAERO软 件，计算给定飞行状态下的气动力影响系数矩阵[AIC];上述飞行状态与步骤1中计算CFD 数据的飞行状态一致； 步骤6,采用分段斜率修正面元法对气动载荷分布进行修正，获得修正后的载荷分布； 利用步骤1得到的定常刚体CFD计算数据对面元法的气动力影响系数矩阵[AIC]实施 修正，具体是采用公式（5)进行分段线化高精度面元修正：15) 式中，η为迎角范围划分的子区间数目，[WTJ1为第i个迎角区间采用斜率修正法得到 的修正因子矩阵，[Cori1为各面元在第i个迎角区间内的下洗速度值占总下洗速度的比例 系数对角矩阵； 步骤7,计算CFD数据与修正面元法结果的误差范数，判断是否达到收敛条件； 若是，则得到最优网格分布，迭代过程结束；反之采用自适应模拟退火算法更新网格分 布参数，更新后转至步骤4,重复步骤4~步骤6 ;在重复步骤4~步骤6的过程中不断获得 新的修正后的载荷分布；对得到的新的修正后的载荷分布通过本步骤计算CFD数据与修正 面元法结果的误差范数，判断是否达到收敛条件，直至满足收敛条件。2. 如权利要求1所述面元修正与网格预先自适应计算方法，其特征在于，在确定不同 迎角下机翼部件原刚体气动载荷数据时，翼面采用三角形单元划分，最小尺度为1mm ;附面 层用三棱柱网格划分，共20层，附面层之外的空间网格采用四面体单元；从附面层最内层 开始网格尺度以1. 2倍的体积比率增长；CFD数据由FLUENT计算，采用耦合压力基求解器； 空气假设为理想气体，温度为300K，参考压力为海平面大气压，湍流模型采用S-A模型，其 余参数选择默认设置。3. 如权利要求1所述面元修正与网格预先自适应计算方法，其特征在于，在对展向布 置网格节点时，先对展向三次曲线的横轴按照展向节点数量进行等距划分，此时横轴上各 点在三次曲线上对应的纵坐标即表示此节点的展向比例位置，且横轴正方向代表从机翼翼 根到翼尖的方向；在对弦向布置网格节点时，先对弦向三次曲线的横轴按照弦向节点数量 进行等距划分，此时横轴上各点在三次曲线上对应的纵坐标即表示此节点的弦向比例位 置，且横轴正方向代表从机翼前缘到后缘的方向。4. 如权利要求1所述面元修正与网格预先自适应计算方法，其特征在于，对于弦向三 次曲线，若需要在机翼前缘加密，则应使曲线在（0,0)点的斜率值减小；对于展向三次曲 线，若需要在翼尖加密，则应使曲线在（1，1)点的斜率值减小；若要采用均匀划分形式，则 只需要将（〇, 〇)、（1，1)处的斜率值均设为1即可。
【专利摘要】一种面元修正与网格预先自适应计算方法，用多组不同迎角下的CFD气动力载荷数据对低阶面元法进行分段线化修正，同时采用自适应模拟退火算法对面元计算网格的分布进行优化。通过对面元网格的自适应优化提高了修正面元法的精度，弥补了传统修正面元法依赖于网格分布的缺点。优化后的网格在应用于修正面元法时不仅保持了面元法计算效率高的优点，还能保证机翼整体受力接近CFD数据的精度，有效提高气动弹性优化迭代设计过程中的精度和效率。本发明得到的气动载荷数据与CFD计算结果精度误差在2％以内，且有效延伸到气动载荷随迎角变化的非线性段，提高了结构刚度参数变化过程中的气动载荷计算效率。
【IPC分类】G06F17/50, G06T17/30
【公开号】CN105183996
【申请号】CN201510582632
【发明人】孙秦, 刘祥, 刘琳颖, 贾欢
【申请人】西北工业大学
【公开日】2015年12月23日
【申请日】2015年9月14日