基于伴随方程方法的非结构网格自适应技术及应用

李 立，白 文，梁益华

(1.中国航空工业集团公司第六三一研究所，陕西 西安 710068;2.中国航空研究院航空数值模拟技术研究应用中心，北京 100012)

0 引言

随着技术进步，计算流体力学(CFD)已成为飞行器气动力设计的例行工具和日常应用手段。但计算网格生成一直是制约CFD技术发展的一个瓶颈。通常对于复杂的计算问题，整个数据准备的70%以上工作量都被用于生成计算网格。近年来，非结构网格(尤其是非结构/结构混合网格)因其对复杂几何外形的适应能力和良好的自动生成能力，在复杂外形流场模拟中的应用受到越来越多的重视［1］。本文围绕非结构自适应网格生成技术开展研究，其目标是改进非结构网格技术在复杂流动物理流场的模拟能力。

网格自适应技术的核心是构造各种自适应探测器。不同的探测器对应不同的网格自适应方法［2］。目前，主要有两类构造自适应探测器的典型方法。一类是基于流动特征的传统方法，如采用流场变量的差量、梯度等作为探测器。该类方法普遍认为［3］，流场中梯度变化大的位置对流场计算精度的影响越大。主流商业软件CFX、Fluent等所采用的策略都属此类。一类是基于敏感性分析的新兴方法，如本文研究的基于流场伴随方程的网格自适应探测器。与传统方法相比，基于伴随方程的网格自适应方法是一种基于全局后验误差估计的方法，直接建立了流场误差与目标输出(如升力、阻力和力矩等)误差的关联关系，可有效避免因错误判断加密区域导致的网格过度加密［3］，提高计算效率和目标输出的计算精度［4-5］，从而有效削减气动数值模拟对计算网格的敏感性。

本文对基于伴随方程的网格自适应技术的基本原理进行了研究，构造了基于流场伴随方程的网格自适应探测器(伴随自适应探测器)，利用局部加密网格的办法，建立并实现了基于Euler方程的无粘伴随自适应能力。对于粘性流动问题，提出首先应用无粘计算和伴随自适应从基准网格自动生成具有较好网格分辨率的非结构自适应网格，然后采用层推进法自动生成粘性网格的求解策略。通过对具体算例开展的数值试验，表明本文建立的方法是有效的。

1 数值方法

1.1 流场主控制方程的计算方法

出发方程为Euler方程或雷诺平均的N-S方程，应用格点格式的有限体积方法离散，控制体采用对偶方法生成，空间离散采用标准的Jameson中心差分格式，时间推进采用三阶显式Runge-Kutta格式。另外，在实际计算中，均采用三重多重网格W循环的聚合多重网格算法进行计算加速。对于粘性计算，湍流模型采用Menter SST k-ω 模型。

1.2 流场伴随方程的计算方法

流场伴随方程的定义式如下［4-5］:

其中，λ即伴随变量，R为流场控制方程残差，Q为流场解变量，f为目标输出函数(如升力、阻力、力矩等)。显然，当流场解已知，伴随变量就可以通过求解式(1)得到。为了解的稳定性和计算过程的鲁棒性，本文采用求解流场主控制方程的类似办法，通过在伴随方程中添加虚拟时间项进行时间推进求解［7］:

这里，Rλ(λ)为伴随方程的残差，

实际求解时，对式(2)，与流场主控制方程一样，采用格点格式的有限体积方法离散，空间离散采用标准Jameson中心差分格式，时间推进采用三阶Runge-Kutta显式格式，并采用多重网格加速。

1.3 伴随自适应探测器

通过伴随变量可直接建立目标输出残差与流场残差之间的关系［4］:

式中，(λ0，Q0)表示初始网格上的伴随解和流场解;(λ*，Q*)表示理论解或非常密网格上的解，一般无法直接给出;R(Q0)表示粗网格上流场控制方程的残差;表示在非常密网格上伴随方程的残差。

由式(4)不难分析，影响目标输出的误差项可分为两部分。一部分是可利用现有网格直接计算的绝对误差项(λ0)TR(Q0)，而剩下的一部分是需在更密网格上计算才能得到的修正误差项。这表明，如对修正误差项加以估计，就可判断计算网格中哪些位置对目标输出有更大的影响，进而可用于建立新型的网格自适应探测器，即伴随自适应探测器。本文中，伴随自适应探测器定义如下:

通过估计或判断该值大小就可标记出流场中哪些位置的网格最需要加密。

1.4 自适应网格生成策略

自适应网格生成的一般性策略和方法在文献［2，8］有较详细的介绍。本文采用局部加密网格的办法(又称h-方法)来进行自适应网格生成。h-方法通过细划网格单元来实现自适应网格生成，通常分为点插入、粗网格细分和网格重构等三类方法。本文采用的方法属于粗网格细分法。该方法依赖于自适应探测器所标记的边。根据边的标记情况，对于四面体，共有三种细分网格的模式，分别为二分模式、四分模式和八分模式。本文自适应网格生成策略如下:首先通过伴随自适应探测器对网格单元进行标记，然后根据网格单元的标记情况，采用二分模式进行各向同性加密。对四面体单元，其方法是连接四面体最长边的中点和该边所对的2个点，从而将四面体分成2个子四面体。为了利用式(5)定义的自适应探测器对网格单元或边进行标记，公式中的理论解或非常密网格上的解首先必须用从原始粗网格的解(λ0，Q0)往一个全局加密网格上插值得到的解(λh，Qh)代替。这里，全局加密网格采用八分模式进行网格划分来自动生成。为了简单，文中采用分段线性插值来计算全局加密网格的解。通过循环扫描粗网格的每个单元所包含的全局加密网格的网格点并由式(5)来估计当地实际的修正误差:

自适应准则定义如下:用户给定全局误差阀值ε，按照误差等分布的原则，判断当地误差与平均全局误差的比值η=N·εk/ε是否大于1，如该条件成立，则对所在的单元进行标记。这里，N表示现有粗网格的网格点数。

值得指出，采用二分模式对四面体网格进行自适应剖分可能会导致自适应后的网格包含悬挂点或非四面体单元，并且还有可能使自适应网格的网格质量比较差，因此采取措施消除悬挂点和对网格进行光顺非常必要［9］。推荐将自适应后的网格导入通用网格生成软件进行必要的处理。

1.5 粘性网格生成

为了捕捉边界层内的粘性流动，必须在沿物面法向的一定厚度内生成各向异性网格，即需将原始的纯四面体非结构网格转换成混合网格。混合网格中可能包含的网格单元类型有四面体、棱柱、六面体或金字塔单元。有多种方法可以实现从四面体网格自动生成混合网格，本文推荐采用层推进方法［10］。其思路是在物面附近按期望的边界层厚度通过网格局部重构产生贴体的棱柱单元，而在这一厚度之外的大部分网格单元保持不变，过渡单元采用四面体或金字塔单元填充。

2 程序校验和方法验证

2.1 RAE2822跨声速翼型计算

RAE2822跨声速翼型是一个超临界翼型，在全世界得到广泛研究，其试验结果被用于大量流场解算器的验证［11］。这里，选取的计算状态为 Case 9:M∞=0.734，攻角α=2.79°。初始计算网格如图1(a)所示，共有15852个网格点，31388个三角形网格单元。图1(b)给出由初始网格得到的马赫数云图。可以看到，该初始网格较好地捕捉到流场的激波结构。选用升力系数作为目标输出函数进行自适应，设定全局误差阀值为0.1。自适应后的网格如图1(c)所示，共有20763个网格点，41141个三角形单元。可以看到，伴随自适应加密的区域主要在翼型的前缘、上表面和激波区，这和主要流动的特征非常匹配，也正是我们希望加密的区域。按伴随自适应探测器的基本原理，这同时还表明，RAE2822翼型的前缘、上表面和激波区对于准确计算升力有主要贡献。

2.2 ONERA-M6跨声速机翼计算

图1 RAE2822跨声速翼型的伴随自适应结果Fig.1 Effect of mesh adaptation for RAE2822 airfoil

图2 ONERA-M6跨声速机翼的伴随自适应结果Fig.2 Effect of mesh adaptation for ONERA-M6 wing

图3 自适应前后ONERA-M6机翼物面沿展向压力系数分布的对比Fig.3 Comparison of spanwise pressure coefficient distributions between before and after adaptation with experimental data for ONERA-M6 wing

与RAE2822翼型一样，ONERA-M6机翼也是一个经典的验证算例［12］。这里，选取的计算状态为:M∞=0.84，攻角α=3.06°。初始计算网格为纯四面体网格，包含52144个网格点，296517个四面体单元，其中，机翼物面分布的网格点数为5459。图2(a)给出了初始的机翼物面网格。可以看到，初始的物面网格点在机翼表面分布比较均匀，除在翼尖、机翼前缘位置从几何保形的角度作了处理外，没有从流动特征方面作任何特殊考虑。同样，选用升力系数作为目标输出函数进行自适应，设定全局误差阀值为0.1。自适应后的四面体网格共有128386个点，740366个四面体单元，其中机翼物面分布的网格点数为8041。图2(b)给出自适应后物面网格的示意图。从自适应的物面网格点分布可以很清晰看到一个呈现λ形的加密区域。这也正好对应本计算状态下ONERA-M6机翼表面应呈现的激波结构。从图中可看出，通过自适应，网格点对机翼前后缘也进行了加密。图2(c)给出由自适应网格得到的ONERA-M6机翼物面马赫数云图，计算得到的机翼表面 激波非常清晰。图3比较了自适应前后计算得到的机翼物面沿展向压力系数分布。从压力系数(尤其在靠近翼根的位置)的比较，可以看到，对于ONERA-M6机翼，采用伴随自适应还可在一定程度上提高计算的精度。

3 复杂流场计算应用

这一节讨论伴随自适应技术在复杂流场计算中的应用。选择的考核算例是第二届国际涡流试验项目(VFE-2)中的尖前缘65°三角翼模型［13］，主要考察伴随自适应技术在提高基于非结构网格的CFD技术预测大攻角复杂涡流场的能力。VFE-2项目是美国和欧洲在2004年联合发起的国际合作项目，采用NASA Langley研究中心提供的不同前缘钝度的三角翼标准模型作为研究对象，旨在确认和评估当前CFD方法在复杂涡流场的预测能力和技术状态。关于该模型的介绍和具体描述参见文献［14］。本文选取的计算状态是:M∞=0.4，攻角α=20.3°。对于粘性计算，基于平均气动弦长的雷诺数Re=2.0×106。为了减小计算量，实际计算中，采用半翼展模型开展计算研究。

3.1 无粘计算结果

初始网格采用纯四面体网格，包含64835个网格点，365681个四面体单元，其中，物面上分布的的网格点为9464个。对于初始计算网格，除了从几何保形方面考虑对三角翼前后缘适当进行加密控制外，没有进行任何特殊处理。图4(a)上半部分给出了三角翼初始的物面网格，可以看到，物面上点的分布几乎是均匀的。图4(a)下半部分是由该初始计算网格采用无粘计算得到的物面压力系数分布云图，给出了一个很明显的低压区。这一低压区正是三角翼上表面前缘涡扫过的痕迹，说明由初始计算网格可捕捉到涡流场的存在。仍选用升力系数作为目标输出函数进行伴随自适应，全局误差阀值设为0.1。自适应后的四面体网格共有139682个网格点，766948个四面体单元，其中，物面上分布的网格点为13333个。图4(b)上半部分给出了自适应后的三角翼物面网格示意图。很明显，自适应加密的物面网格点主要分布在机翼翼根、机翼前后缘以及计算捕捉到的前缘涡所扫过的区域。图4(b)下半部分给出了由自适应网格计算得到的物面压力系数分布云图。与上图比较，不难看到，自适应后得到的低压区更加清晰，这表明捕捉到的前缘涡的强度有所增加。

图4 VFE-2尖前缘65°三角翼模型自适应前后结果的对比Fig.4 Comparison of mesh and pressure coefficient distribution on the wall surface between before and after adaptation for VFE-2 delta wing with sharp leading-edge

图5通过沿流向方向截取两个典型站位x-c=0.4及x-c=0.6的空间网格，对比了自适应前后空间网格的变化。可以看到，对于此算例，空间网格主要在靠近机翼上表面的位置(尤其是翼根和主涡的附近)进行加密，而在机翼下表面附近网格几乎没有任何加密。并且，在流向方向的不同站位，网格的加密位置和密度随主涡涡核的位置和强度有所变化，但并不完全集中在涡核周围。这与基于流动特征的涡自适应或熵自适应方法观察到的有明显差别［2］。

图5 尖前缘65°三角翼模型自适应前后两个典型站位的空间网格对比Fig.5 Comparison of mesh resolution in space from two typical sections for VFE-2 delta wing with sharp leading-edge

图6比较了自适应前后采用无粘计算得到的尖前缘65°三角翼在典型站位的物面压力系数分布。可以看到，自适应前后，采用无粘计算在预测涡核的位置和主涡吸力峰强度方面都有很大改善，但与试验相比，整体上偏差都比较大。这也从侧面印证了人们长期形成的一个计算经验，即采用Euler方程进行无粘计算不足以很好地预测三角翼大攻角涡流场这类存在大的分离和以涡流为主导的复杂流场。

图6 自适应前后尖前缘65°三角翼模型在典型站位的物面压力系数分布与试验值的对比(无粘计算)Fig.6 Comparison of pressure coefficient distributions between before and after adaptation with experimental data for VFE-2 delta wing with sharp leading-edge(in inviscid computation)

表1 两套粘性计算网格的网格信息统计Table1 Information for two meshes in viscous computations

3.2 粘性计算结果

在粘性计算中，粘性计算网格分别在初始四面体网格和自适应后的四面体网格基础上生成。为了保证第一层网格的Y+≈1，边界层内第一层网格距离物面的距离按0.01进行估计。在实际计算中，取第一层网格的实际间距为5.0×10-6，网格拉伸比为1.2，总的附面层网格层数为30层。两套网格的详细信息参见表1。

为了比较计算效果，图7首先给出利用两套粘性计算网格得到的物面极限流线和空间流场解的对比。计算中所有参数设置在前后保持一致。从物面极限流线的结果看，由自适应局部加密的粘性网格得到的解不仅观察到主再附线(说明涡流场存在)，而且还清晰观察到二次再附线，而由初始粘性网格得到的解仅能观察到主再附线，二次再附线不明显。并且，两者主再附线的位置明显有差别，由初始粘性网格得到的主再附线更靠近翼尖的位置。这表明，由自适应局部加密的网格得到的主涡位置会明显沿展向翼根靠近。从空间流场解的结果进行分析，图中采用沿流向方向不同站位的总压比等值线和绕涡核的流线的方式清晰表现出涡流场的存在，但由自适应局部加密的粘性网格得到的解明显观察到二次涡，并且主涡核的位置明显更靠近翼根。因此，两者的分析结论完全一致，即采用伴随自适应局部加密后，粘性计算进行涡捕捉的能力明显增强。另外，还容易注意到，由初始粘性网格计算得到的解还提前观察到涡破裂。这与试验观察到的情况不太吻合［15］。

图7 尖前缘65°三角翼模型的物面极限流线和空间流场解对比(粘性计算)Fig.7 Comparison of flow solutions between before and after adaptation for VFE-2 delta wing with sharp leading-edge(in viscous computation)

图8比较了粘性计算得到的物面压力系数分布。与无粘计算仅给出三个典型站位的结果不同，图8给出了在试验中所测量的所有站位数据的比较。其中，试验结果包含两类，一类由NASA Langley中心的跨音速风洞在雷诺数Re=6.0×106下进行试验得到，一类由德国风洞协会的跨音速风洞在雷诺数Re=2.0×106下进行试验得到，二者有较小的雷诺数效应差别。可以看到，由初始粘性网格得到的各站位上表面压力系数分布与试验值相差甚远，主要表现在压力系数吸力峰位置明显沿展向朝翼尖位置偏移，而吸力峰峰值比试验值偏小;与之相较，由自适应局部加密的粘性网格除了在x/c=0.2这一站位外，在其他四个站位都给出与试验符合较好的结果。自适应局部加密对计算结果的改进非常明显。在x/c=0.2这一站位，推测计算结果没有明显改善的原因是，两套粘性网格在机翼前端附近的位置保持了基本一致的网格分辨率，而这一网格分辨率不足以保证在x/c=0.2站位给出与试验符合良好的计算结果。可能的解决办法是在这一区域附近进行手动加密以保证在机翼前端都有足够的网格分布。

图8 尖前缘65°三角翼模型在典型战位的物面压力系数分布与试验值的对比(粘性计算)Fig.8 Comparison of pressure coefficient distributions with experimental data for VFE-2 delta wing with sharp leading-edge(in viscous computation)

4 结论

对基于伴随方程的网格自适应技术的基本原理进行了研究，构造和实现了基于流场伴随方程的非结构网格自适应探测器，采用局部加密网格的方法，建立了基于Euler方程的无粘伴随自适应能力。对粘性计算，提出首先采用无粘计算和伴随自适应生成网格分辨率好的初始网格，然后采用层推进方法生成粘性计算网格的求解策略。采用二维和三维经典算例对方法进行了验证，并应用本文提出的求解策略对VFE-2尖前缘三角翼大攻角涡流场进行了数值模拟。计算结果表明，采用自适应技术能提高计算精度，尤其对三角翼大攻角涡流场这类存在大的分离和以涡主导的复杂流场，采用本文策略，可明显改善预测结果，有效提高预测准确度。粘性计算时，在网格生成阶段，采用Euler无粘计算进行伴随自适应生成具有较好网格分辨率的初始四面体网格还可在一定程度上减小计算代价，提高计算效率。

［1］张来平，王志坚.三维复杂外形的非结构网格自动生成技术与应用［J］.计算物理，1999，16(5):552-558.

［2］BAI W，QIU Z，LI L.Recent efforts to establish adaptive hybrid grid computing capability at ACTRI［J］.Computational Fluid Dynamics journal，2007，7(4):438-449.

［3］WARREN G P，ANDERSON J T，THOMAS J T，KRIST S L.Grid convergence for adaptive methods［R］.AIAA paper 1991-1592，1991.

［4］PARK M A.Adjoint-based three-dimensional error prediction and grid adaptation［R］.AIAA paper 2002-3286，2002.

［5］NIELSEN E J.Adjoint-based algorithms for adaptation and design optimization on unstructured grids［A］.The 3rd East-West High-Speed Flow-field Conference［C］.Beijing，China，2005.

［6］LEE-RAUSCH E A，PARK M A，JONES W T，HAMMOND D P，NIELSEN E J.Application ofparalleladjoint-based error estimation and anisotropic grid adaptation for three-dimensionalaerospace configurations［A］.The 23rd AIAA Applied Aerodynamics Conference［C］.Toronto，Ontario，2005.

［7］VENDITTI D A.Grid adaptation for functional outputs of compressible flow simulations［D］.Massachusetts Institute of Technology，Cambridge，MA，2002.

［8］PIRZADEH S Z.An adaptive unstructured grid method by grid subdivision，local remeshing，and grid movement［R］.AIAA paper 1999-3255，1999.

［9］SHAROV D，FUJII K.Three-dimensional adaptive bisection of unstructured grids for transientcompressible flow computations［R］.AIAA paper 1995-1708，1995.

［10］TIMOTHY J B.Mesh generation:art or science［J］.Progress in Aerospace Sciences，2005，41:29-63.

［11］http://www.grc.nasa.gov/www/wind/valid/raetaf/raetaf.html

［12］http://www.grc.nasa.gov/www/wind/valid/m6wing/m6wing.html

［13］UMMEL D，REDEKER G.A new vortex flow experimentfor computer code validation［A］.RTO Symposium on"Advanced Flow Management:Part A-Vortex Flows and High Angle of Attack for Military Vehicles"［C］.Loen，Norway，2001.

［14］CHU J，LUCKRING J M.Experimental surface pressure data obtained on 65°delta wing across Reynolds number and Mach number ranges［R］.NASA TM-1996-4645，1996.

［15］KONRATH R，SCHÖDER A，KOMPENHAS J.Analysis of PIV results obtained for the VFE-265°delta wing configuration at sub-and transonic speeds［A］.24th Applied Aerodynamics Conference［C］.San Francisco，California，2006.