国家数值风洞(NNW)工程中的黏性自适应笛卡尔网格方法研究进展

陈浩，袁先旭，王田天，周丹，赵宁，唐志共，毕林,*

1. 中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000

2. 中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

3. 中南大学 交通运输工程学院 轨道交通安全教育部重点实验室，长沙 410075

4. 南京航空航天大学 航空学院，南京 210016

随着计算机技术和计算流体力学学科的飞速发展，数值模拟在湍流研究中扮演着越来越重要的角色。然而，由于湍流流动具有非定常、多尺度、强扩散与耗散等复杂特性，且目前面临的几何外形越来越复杂[1-2]，数值模拟要在认识、解决、应用湍流问题上发挥更重要的作用，还需攻克高质量计算网格生成的难题。

目前，以多块分区结构网格[3-5]、三角形(四面体)贴体非结构网格[6-7]等为基础发展的湍流流动数值模拟方法，在提高湍流认知水平、解决湍流工程应用问题等方面发挥了重要作用。然而，这些网格类型对于湍流这样的复杂流动，要么需要花费大量的人力资源，自动化程度不高，要么网格分布难以满足复杂湍流流动结构动态、精细捕捉要求，自适应能力不强，逐渐成为制约湍流精细化研究的瓶颈问题之一。

相对于传统的分区结构网格和贴体非结构网格，自适应笛卡尔网格[8-11]不依赖于物面网格直接生成空间网格，具有网格生成自动化程度高、复杂外形适应性好、非定常/多尺度等流动结构捕捉能力强等优势，适用于湍流等复杂流动问题的仿真模拟。因此自适应笛卡尔网格方法具有重要的研究意义和广阔的应用前景。

笛卡尔网格也有其劣势，主要是非贴体性带来的物面边界处理困难。对于非贴体的笛卡尔网格来说，网格与物面边界相交，在处理黏性流动问题时面临较大的局限性，尤其是对于边界层内大纵横比流动，雷诺数越高，网格规模越大。因此，通过笛卡尔网格方法实现三维空间效应强、旋涡尺度跨度大的非定常复杂流动问题的高精度求解具有一定难度。针对这种非贴体性，不同学者提出了不同的边界处理方式，主要有：切割单元法和浸入边界法。

切割单元法使笛卡尔网格具有贴体网格的特性，在处理时需要考虑网格被物面切割后的形状和分布。该方法能使物面处满足守恒性的要求，一般结合有限体积方法进行计算，应用较为广泛[12-13]。商业软件CART3D就是基于切割单元法进行处理非贴体物面边界的。中国空气动力研究与发展中心的肖涵山等[14]对切割单元法进行了研究和发展，并成功应用于机弹分离问题的模拟仿真。虽然切割单元法在一些计算中取得了良好的结果，但是这种方法会产生形状多样、大小不一的网格单元，这样不仅会限制时间步长，还会增加数值解的非物理震荡，进而增加求解难度，限制了其在黏性流动问题中的应用。

浸入边界法[15-16]则不需要对在物面附近的单元进行特殊处理，保留了笛卡尔网格原有形状，故此方法有一定的便捷性和灵活性。该方法最初由Peskin[16]在1972年提出，用于模拟人类心脏中的血液流动。这种方法的思想是：不直接考虑物面的存在，而是把物体作为边界条件嵌入到流场中，对边界的处理则是将边界模化成流动控制方程中的的力源项。Mohd[17]将浸入边界法与B-Spline法结合来获取高阶精度，并且通过用时间离散导数去掉时间步长的限制，大大扩展了浸入边界法的使用范围。Dadone和Grossman[18-20]在浸入边界法的基础上，提出了虚拟单元方法(Ghost Cell Method)，该方法是在Dadone提出的物面曲率修正技术(Curvature Corrected Symmetry Technique)的基础上演化而来，最初应用于贴体网格，之后拓展到笛卡尔网格，他成功地将这种方法应用于二维和三维无黏流动问题的数值模拟，并在之后结合数值壁面函数，实现了对于高雷诺数可压缩黏性流动问题的模拟。在上述工作的基础上，Lee和Ruffin[21]将浸入边界方法推广到黏性问题的流场计算中，提出了一种新的壁面函数，并针对湍流问题进行了模拟研究，取得了较好的效果。国内在浸入边界法和自适应笛卡尔网格技术相结合的研究领域成果则相对较少，刘剑明[22]、胡偶[23]对自适应笛卡尔网格技术和浸入边界法进行了发展，在黏性可压缩流动问题中作了应用研究，并扩展到了高雷诺数黏性流动。陈浩等[24]将上述技术方法进一步扩展至超声速和高超声速黏性流动问题，在二维典型算例中得到了可靠的结果。上述研究所面临的外形相对简单，对于复杂构型的扩展应用仍然有待开展。

为了保证求解精度，笛卡尔网格在自适应加密过程会产生大量网格。尤其当求解三维黏性流动问题时，为降低非物理振荡、保证模拟精度，非贴体的笛卡尔网格会在边界层区域进行不计成本的加密，使得相同模拟精度下，笛卡尔网格明显大于贴体类网格，甚至会有量级上的差距。巨大的网格量会导致计算效率低下以及存储量过大的难题。在应对上述问题时，不同学者形成了两种研究思路：

1) 在网格技术上做改进，主要包括：混合笛卡尔网格技术、各向异性自适应技术以及法向射线加密技术等。

其中，混合笛卡尔网技术相对使用较多，将贴体结构网格和笛卡尔网格的优势相结合，通过在近壁使用贴体网格来避免笛卡尔网格边界处理的困难。网格生成与空气动力研究国际权威普林斯顿大学Baker[25]认为混合网格技术是易用度和黏性模拟精度兼顾较好的方法。Charlton等[26]较早提出了内部贴体网格和外部直角坐标网格相结合的混合网格思想。之后，Delanaye[27]在Beger[28]与Aftosim等[29]的工作基础上发展了贴体结构与笛卡尔网格相混合的Navier-Stokes(N-S)求解器。目前，国外已发展了多个基于笛卡尔/结构网格混合的知名网格技术求解器，例如，NASA著名的OVERFLOW[30]求解器就是使用笛卡尔背景网格耦合贴体结构网格的形式。此外，由美国国防部支持开发的旋翼计算平台HELIOS，是基于笛卡尔/贴体非结构网格混合的形式。Wang和Chen[31]提出的 “黏性自适应笛卡尔网格方法” 也是混合网格方法，其是在切割单元方法基础上，通过投影方式构造近壁贴体结构网格。之后，Fujimoto等[32]继续对Wang的“黏性自适应笛卡尔方法”进行改进，使之能够处理包括缝隙、凹槽、凸起等结构的更为复杂的几何外形流动。

相比于国外形成的较为成熟的网格求解技术与求解程序，国内在基于笛卡尔网格的混合网格技术方面的研究相对有限，相关成果也较少。肖涵山等[33]在其发展自适应笛卡尔网格技术基础上，借鉴“投影”方法的思想，实现了混合笛卡尔网格的自动生成，同时构建了Euler数值求解器，形成了飞行器气动力快速预测和评估工具。张来平等[34]在混合网格方法上进行了较为深入的研究，涉及到了混合笛卡尔网格方法的相关工作。他们针对不可压非定常流动问题，在近壁面采用贴体结构网格，远场采用笛卡尔网格，中间用非结构网格连接过渡，建立了动态混合网格方法。田书玲等[35-37]发展了混合笛卡尔网格技术，并对网格生成、混合重叠网格装配和动态重叠网格算法开展了深入研究。借鉴传统重叠网格思想，胡偶[23]、沈志伟[38]在其所建立的自适应混合笛卡尔网格方法基础上，开展了非定常旋涡主导流动问题的数值模拟与流动机理分析。以上工作显著促进了国内混合笛卡尔网格技术的发展，部分成果已应用于相对复杂的流动问题研究。

各向异性笛卡尔网格技术是利用流场中很多区域是各向异性的特点，发展与之相匹配的网格加密方法，实现在保证网格质量的同时明显降低网格规模。传统的笛卡尔网格方法多是采用各项同性加密的方法，事实上，流场中的很多区域是各向异性的，如附面层、激波间断、剪切层等，在这些区域如果还采用各向同性的加密方法，势必会造成网格点的大量浪费。因此，很多学者将各向异性加密技术引入到笛卡尔网格方法当中。Ham等[39]针对不可压缩层流非定常流动问题，基于非结构数据结构，建立了各向异性笛卡尔网格加密和粗化方法，显著减少了网格数目和内存占用。在此基础上，Keats和Lien[40]将该方推广到可压缩黏性流动问题的模拟当中，经验证，他们所发展的方法在降低网格数目的同时能够保证计算精度。Wang和Chen[31]提出的 “黏性自适应笛卡尔网格方法”正是采用了基于2N叉树数据结构的各向异性流场解自适应方法和“投影”方法相结合，实现了对于湍流黏性流动问题较高精度的模拟。Wu和Li[41]也对各向异性的网格加密技术进行了研究和发展，他们最初是用各向异性网格方法捕捉黏性边界层中的固有各向异性特性，在经修正之后，在捕捉斜激波和正激波方面取得了较好的效果。Wu和Li[42]分别将各向同性和各向异性的网格加密技术应用到笛卡尔网格方法中，进行对比研究，验证了所发展的各向异性解自适应方法加密笛卡尔网格可以有效地捕捉流场中激波、接触间断等关键区域。Sang和Shi[43]基于全叉树(Omni-Tree)数据结构发展了自适应笛卡尔网格方法，相对于传统的仅支持各向同性加密的八叉树数据结构，该数据结构同时支持各向同性和各向异性两种加密方式，通过对民用飞机这种复杂外形开展的模拟仿真和对比研究，证明其所发展的各向异性笛卡尔网格技术能够在一定程度上减少网格数目和提高计算效率。

法向加密技术则是借鉴贴体结构网格大纵横比的特点，通过沿壁面法向配置的射线，调整附面层内的笛卡尔网格的粗细分布，减少附面层内网格数目。这种方法最初由Ruffin等[44]在2012年提出，他们针对边界层高精度模拟需要大量网格的问题，提出边界层法向加密(Normal-Ray Refinement)的概念，并建立了物面法向射线的分布方式以及不同射线之间加密网格的数据交互方法，极大减少了网格点数目和内存占用，节省了计算时间。但是这种方法需要人为地指定法向加密的位置，在非加密区其壁面处理精度较低，一般用于求解层流黏性流动问题，对于边界层特征很复杂的湍流问题的高精度模拟仍然需要发展。

2) 引入通过简化方法刻画近壁流动的思想，降低壁面网格尺度要求，以达到减少网格量的目的，主要包括近壁模型(Wall Model)和壁面函数(Wall Function)方法。

其中，近壁模型方法通过对湍流模型进行修正，使其能够对于近壁黏性影响区域进行求解[45-46]。壁面函数方法则是考虑到边界层存在解析解这一特征，通过构造一维代数模型来代替黏性底层的直接求解，使得对壁面法向第1层网格的y+需求放宽至50～100，甚至200，大大降低计算资源的使用[21]。相比近壁模型方法，计算量相对较小，且更容易实现，在贴体类网格方法中得到了广泛应用。

国外研究人员将近壁模型和壁面函数方法在笛卡尔网格框架下进行应用，定义湍流壁面边界条件，从而降低边界层内计算网格数量。Kalitzin和Iaccarino[47]提出了一种基于壁面律的浸入边界方法处理高雷诺数湍流问题。Lee和Ruffin[21]通过壁面函数模型耦合k-ω湍流模型定义湍流壁面边界条件，并用于高雷诺数问题的数值模拟，包括二维平板、翼型以及三维旋翼机的机身旋翼相互作用等。Dadone和Grossman[18-20]在其无黏虚拟单元方法的基础上，提出在壁面处构造数值壁面函数结合虚拟单元的方法，求解高雷诺数流动问题。Capizzano[48]则采用双层壁面模型修正浸入边界网格交接面上湍流变量值的方法，进行高雷诺数可压缩流动问题的模拟研究。Bernardini等[49]基于浸入边界方法，通过壁面模型耦合SA模型，用于延迟脱体涡模拟(Delayed Detached Eddy Simulation, DDES)中，对于高雷诺数下壁面分离的湍流流动开展了模拟研究，获得了高精度、可靠的数据结果。这些学者的研究工作对基于笛卡尔网格技术进行湍流问题的模拟起到了极大的推进作用。

国内有关笛卡尔网格方法在高雷诺数黏性流动问题中的研究和应用起步较晚，成果也较少。刘剑明[22]基于自适应笛卡尔技术，发展了针对可压缩Navier-Stokes方程的径向基函数浸入边界方法，结合Menter 剪切应力输运(SST)两方程湍流模型，对于高雷诺数下的圆柱和NACA0012翼型等简单外形的绕流问题进行了研究，取得了一定的成果。胡偶[23]针对SSTk-ω模型湍流模型首次构造了一种壁面函数——虚拟单元方法(Wall Function-Ghost Cell Method, WF-GCM)，实现了湍流壁面边界重构在“非贴体”笛卡尔网格上的应用。他的研究对象比较简单，雷诺数量级也不是很高，对于湍流边界层的流动特征研究较少，有进一步发展的需求和空间。

虽然国内外学者对于黏性笛卡尔网格技术已经开展了部分研究工作，促进了笛卡尔网格方法在复杂黏性流动问题中的应用。但是由于自适应笛卡尔网格技术模拟黏性流动难度较大，且对于计算资源要求较高，国内外对于黏性自适应笛卡尔网格方法的研究相对其他网格技术而言仍然偏少，尤其是国内，尚处于起步阶段，公开的学术成果极少。因此，目前对于黏性自适应笛卡尔网格方法国内外尚缺少系统性的综述介绍，亟待相关工作进行补充。

国家数值风洞(NNW)工程关注黏性自适应笛卡尔网格方法的发展潜力和应用前景，资助了若干研究课题，力争突破关键技术瓶颈，促进实现黏性笛卡尔网格技术的工程实用化。本文聚焦国家数值风洞工程，总结介绍了黏性笛卡尔网格方法相关课题的研究进展，重点从笛卡尔网格高效生成技术、笛卡尔网格自适应技术以及黏性物面边界处理等3个方面展开。

1 网格高效生成技术

1.1 笛卡尔网格数据结构

网格的数据结构是指用来存储如位置、物理参数、邻居关系等网格节点信息的一种结构形式，它是进行数值计算的前提和基础。合适的网格数据结构可以有效减小内存的需求并提高访问速度，为网格的高效生成提供支撑。

目前笛卡尔网格的数据结构主要分为以下3种：结构式、非结构式和叉树式[48-52]。结构网格由于其相邻网格之间存在自然的联系，只需要采用I、J、K的顺序对数据进行存储即可，使用起来简单且易于编程实现，但灵活性较差，自适应困难。非结构网格在存储自身网格流场信息的基础上，还要存储其邻居网格单元的编号及其与自身网格相邻边的相对关系信息，导致其数据结构复杂，存储量大，访问效率相对较低。叉树数据结构是一种分层式的数据结构，依据网格大小划分出不同的层次，由父单元即上一层单元决定其网格的位置，同样的，该单元也可以划分成数个子单元，并通过指针来建立不同层级之间父子单元的访问和联系。对于自适应笛卡尔网格而言，网格的加密过程实质上是对空间沿着垂直坐标方向进行等分(一维为二等分，二维为四等分，三维为八等分)，因此传统的叉树数据结构在自适应加密和粗化、网格调用方面具有很大的优势，但是其所存在的问题也比较突出：查询网格单元的邻居单元所需要的搜索计算量大，例如最差情况需要从叶子节点回溯到根节点，然后从根节点的邻居节点搜索到叶子节点；在网格进行自适应后需重新确定每个网格单元的邻居关系；存储网格信息所需要的存储量较大。此外，叉树数据结构还存在一些无用的内存分配，例如，叉树数据结构的每个单元都有子单元指针，但是这些指针对于叶子单元而言并不会被用到，这就会造成内存的浪费[28]。

为了解决上述传统叉树数据结构的问题，本文作者团队在笛卡尔网格存储时采用了叉树数据结构的改进形式FTT数据结构[51-53]，如图1[52]所示。它将叉树数据结构中叶子节点闲置的指针利用起来，指向邻居单元的父单元，从而构建快速访问邻居单元的“线程”，实现在不同检索操作中数据信息的快速访问，提供了高效、易并行的叉树信息访问路径，同时节省数据节点，降低叉树信息存储的内存占用,大大提高计算效率。

图1 叉树数据结构 vs FTT数据结构[52]Fig.1 Fork tree data structure vs FTT data structure[52]

总体而言，采用FTT数据结构进行邻居单元的查找，平均遍历的数据层数是传统的邻居查找方式的1/4，这样在搜索邻居上可以节约大约75%的时间。在内存占用方面，采用传统的叉树数据结构，要存储完上述信息，每个单元需要17 words的内存，而如果采用FTT方式，每个单元需要2 words的内存，相对于传统方式明显降低[52]。

1.2 物面网格数据结构

在生成网格的过程中，网格类型的判断占据了大部分时间。为了高效判断网格类型，Bonet和Peraire[54]采用ADT(Alternating Digital Tree)数据结构来存储物面网格，以缩小可能相交的物面单元检测范围，从而提高相交判断的效率。在该方法中，每个ADT节点都可作为父节点或子节点对应一个笛卡尔空间，父节点包含其所对应所有子节点的笛卡尔空间[55]。ADT数据结构利用以上的层次包含关系，只要父节点不与目标空间相交，则该父节点所包含的子节点均将在查询中被忽略，从而大大减少查询次数，提高了搜索效率。然而传统的ADT所构建的二叉树的平衡性往往受物面单元的分布密度影响，导致查询效率不稳定。为此，本文作者团队采用其改进版本KDT(K-Dimensional Tree)结构，其更容易构造平衡优质的二叉树[56-57]。

KDT是一种多维的二叉树结构，KD树的每一层都和某个维度相关，每个节点都是一个k维节点。KDT的叉树在构建时，只需要将当前层的相关点进行排序后取其中值点，随后确定经过该点与相关维度垂直的平面,用来剖分空间的搜索区域，与ADT数据结构在区域等分点进行划分不同，KDT可以保证多层叉树的平衡性并减少树的深度，从而便于查找单元信息。如图2所示[51]，对于图2(a) 中所示的物面网格，KDT叉树数据结构相比于ADT叉树数据结构更加平衡，其相应的操作复杂度也明显降低。对于物面网格较少的情况，宜采用ADT来存储物面单元；当物面网格量较多时，则宜采用KDT。这是因为，KDT每进行一次空间划分，均需要进行一次排序操作，当网格量较大时，树的深度将会比ADT小很多，此时这种排序的优势才会凸显。

图2 随机分布的网格点和ADT、KDT数据结构 [51]Fig.2 Randomly distributed grids and ADT、KDT data structure [51]

1.3 相交判断方法

笛卡尔网格的相交判断是开展自适应和浸入边界处理的前提。基于网格相对于物面的位置关系，可将网格划分为外部网格、内部网格以及相交网格单元。实践表明，相交判断占据了笛卡尔网格生成的大量时间，并且是网格生成成败的关键因素。

Aftosmis等[58]以及Dadone和Grossman[20]采用简单且有效的射线相交方法进行网格点的类型判断。该方法通过计算由目标点生成的射线与物面的交点个数，来判断网格点与物面的相对位置：相交数目为奇数，网格点在物面内部；相交数目为偶数，网格点在物面外部。然而在使用射线相交方法时存在局限性，即可能存在射线恰好经过外形尖点的情况，导致判断的结果出现错误。

针对上述传统射线相交方法存在的问题，本文作者团队提出了一种新的鲁棒的判断算法，与射线相交方法相配合，可以准确识别射线经过外形拐点的不同情形，称为奇异性判断算法[52]。这种方法通过判断相应的线段或面与射线的相对位置从而更加准确地判断目标网格点的类型。

图3[52]显示的是该方法在二维情况中的应用，在图3(a)中，射线R以待判断的目标点T为起始点，经过离散线段之间的尖点P，两条以P为端点的线段分别为PA和PB。PL和PR均在射线的同一侧，因此目标点为物面外的点。在图3(b) 中，尖点处的线段PA和PB在射线R的不同侧，则目标点为物面内的点。其中，r1为以P为起点，方向与R相同的矢量，r2则为方向相反的矢量。PA对应的向量l1(对应外法向n1)和射线的反向矢量以及PB对应的向量l2(对应外法向n2)与射线R的夹角分别为α和β。

图3 射线经过拐点[52]Fig.3 Ray passing inflection point [52]

针对三维外形，可通过切割平面方法建立目标射线与相交线的相对位置关系，从而将三维问题转化为二维问题。以图4[29]为例，在射线R与切割线PA和PB组成的切割平面内，即可按照前面的二维问题处理方式来识别射线与离散三角形的相对位置关系。其中，nPA和nPB分别为切割线PA和PB对应的平面内的外法向，R为过目标点T的射线。

图4 通过切割面将三维问题转化为二维问题[29]Fig.4 Conversion of 3D problem to 2D problem by cutting plane passing ray[29]

对于网格在物面内外的判断，Aftosmis等提出了染色算法[58]，以快速确定与物面不相交的网格单元类型。这种算法通过查询目标网格单元相邻非物面相交单元的类型，直接给予其相同的网格类型，目标单元就像被邻居单元“染色”，故称为染色算法。肖涵山等[33]提出了一种与之相类似的“推进查找”方法，以判断网格单元是否在物面内部。该方法首先将与物面相交的单元删除，这样便切断了外部单元和物面内部单元的联系，然后将流场外边界上的单元标记为外部单元，将与之相连的也标记为外部单元并进行推进，从而快速准确地判定出内外单元。

本文作者团队基于上述方法思想，结合FTT数据结构，发展了适用于自适应笛卡尔网格的染色算法[52]。如图5所示，首先在外部和内部分别找一种子单元，并标记为OutCell和InCell，接下来对其邻居单元进行判断，只要是不与物面相交的邻居单元便标记为与种子单元相同的网格类型，从而快速确定网格单元是在物面内部还是外部。从表1中可以看出，使用改进后的染色算法相比于传统的判别方式，所用的时间大大减少。

表1 染色算法 vs 传统判别方法

由于空间笛卡尔网格单元都是标准的AABB(Axis-Aligned Bounding Box)[59]，因此另一种笛卡尔网格单元与三角形的相交判断方法可归结为AABB与三角形的相交判定，如图6所示。该种方法基于独立轴理论，此理论可表述为当两个凸多面体A和B满足以下两个条件之一则可判定二者不相交：①A和B可以被平行于A或B的任意一个面的法向的轴分离；②A和B可以沿着由A与B的边叉乘形成的轴分离。在具体实现过程中，需要对13根轴进行测试：前3根轴为AABB的面法向轴(e1,e2,e3)，用来判定三角形的AABB与笛卡尔网格的AABB是否有重叠，如果没有，则不相交；第4根轴为三角形的面法向轴(n)，用来判定AABB是否与三角形所在平面相交；最后9根轴为AABB的3个面法向轴分别与三角形的三条边叉乘形成的轴。以上13根轴的测试均通过，则表明没有分离轴存在，即AABB与三角形不相交；否则，只要有分离轴存在，即可结束测试，AABB与三角形相交。通过这种方法可以对待判断网格单元进行快速筛选，减小一定的计算量。

图6 AABB与三角形相交判断示意图Fig.6 Diagram of intersection judgment between AABB and triangle

1.4 并行加速技术

笛卡尔网格在复杂构型/流动问题中应用时存在网格规模过大的难题，为了保证网格生成的效率，有必要开展大规模并行计算技术的应用。事实上，高可扩展性、高负载平衡的并行流场算法和实现仍是目前国际上的一个难点，是CFD模拟中提高计算效率的重要环节[60]。若在计算分区过程中无法做到负载平衡，则并行效率会急剧下降。在非结构网格和混合网格并行计算方法中广泛使用的网格分区软件如METIS和Zoltan等，为了保持分区的负载平衡，每次笛卡尔网格进行解自适应后都需要重新对网格进行分区，这在处理非定常流动问题时，由于笛卡尔网格需要频繁地进行自适应，会在分区上耗费大量的时间。针对以上问题和现状，当下需发展高效的适用于自适应笛卡尔网格的并行策略。

本文作者团队针对自适应笛卡尔网格下高性能并行计算技术的应用开展了深入研究。在流场计算时基于空间填充曲线(Z型排序是Z曲线，如图7所示[61]，逆时针排序是Hilbert曲线)按照排序规则将所有子节点连接起来形成一维链表型数据结构；然后，采用网格分区软件METIS的思想，直接对一维链表式数据进行分区。此外，在流场计算过程中，随着解自适应的进行，该方法还便于实现动态的网格分区，很好地保障计算负载平衡。同时，前述网格自适应加密时遵循2:1加密准则，即加密的网格单元和其周围单元的加密层级最多不超过1，结合网格单元快速ADT搜索算法，可开发最小化的分区交接面信息交换算法，实现分区网格之间信息的快速传递。针对三维等熵涡问题对动态分区效果、负载平衡以及并行效率等方面展开了技术验证，如图8所示。可以看出随着涡的运动，网格进行了自适应并实现了自动分区，与此同时也较好地保持了物理区域的连续性，能够有效地减少进程之间的通讯开销，进而提升计算效率。随着进程数的增加，并行效率保持情况良好，拥有可行的并行规模扩展前景。

图7 排序方案和Z型空间填充曲线示意图[61]Fig.7 Schematic diagram of sorting scheme and Z-type space filling curves[61]

图8 三维涡运动问题网格分区及计算时间占比Fig.8 Grid partition and calculation time of 3D vortex motion problem

2 网格自适应技术

2.1 各向同性自适应

上述工作为构建面向任意外形的自适应笛卡尔网格生成方法奠定了基础。接下来根据几何和流场特征进行网格的自适应加密或粗化。

由于初始网格一般尺寸较大，无法准确识别模型的几何信息和流场的特征结构信息，需要建立相应的加密判据。对于基于几何特征信息进行的加密，一般称为几何自适应。其基本思路是：通过射线相交或染色算法判断出相交网格后，对于与物面相交的单元进行标记和初步加密一定的次数。为了保证物面与流场之间的网格过渡均匀，一般对过渡层网格也进行加密。之后，通过识别几何的曲率、厚度等特征信息，对精细结构进行进一步加密，以保证模型信息的高保真度。最后，对于网格质量进行检测和优化，即对于加密不合理、影响叉树结构平衡性、质量不高等区域进行排查和修正，直至达到需要的高质量网格。例如，对需要加密的网格进行合适的细化、对相邻网格层次差大于1的网格进行降层等。基于构建的几何自适应方法，针对多种复杂构型开展了算例测试，以导弹和翼身组合体为例，如图9所示[52]。

图9 几何自适应网格[52]Fig.9 Geometrical adaptive mesh[52]

笛卡尔网格在流场结构自适应方面具有天然优势。通过流场自适应，可以提高对流动特征结构的分辨率，降低数值耗散，提高流动模拟的精度，同时也能够使网格分布更加合理。为了精细捕捉流动中出现的涡旋、激波、剪切层等特征结构，需要构建适用性好、精度可靠的判据。为此，本文作者团队发展了一种综合判据，即速度散度和旋度相结合的判断方法[52]。该判据既可以通过速度的散度捕捉激波，又可以基于速度的旋度捕捉涡旋和剪切层等特征，相对于单一判据具有明显的优势。针对典型算例测试了构建的自适应方法，得到了较好的效果。流场自适应结果如图10 所示。

图10 超声速流动自适应网格Fig.10 Adaptive mesh based on supersonic flow

2.2 各向异性自适应

各向异性自适应即笛卡尔网格在加密过程中会存在多种加密方式，从而生成形状、大小不一的子单元。以二维各向异性为例，一个父单元有4种加密方式，分别为在x、y两个方向同时加密,在x、y其中一个方向加密以及在x、y两个方向均不加密。加密示意图如图11所示。

图11 二维各向异性笛卡尔网格加密形式Fig.11 2D anisotropic Cartesian mesh refinement form

在各向异性自适应加密过程中，会产生方向、大小不同的相邻网格，因此邻居单元查询算法比各向同性更为复杂。为便于显示，以二维为例说明可能出现的邻居情况，如图12所示。为了防止相邻单元层级差过大导致解不稳定的情况，本文作者团队采用相邻单元层级差不大于1的加密/粗化原则，从根本上解决各向异性笛卡尔网格邻居查询和网格生成质量问题。

图12 各向异性笛卡尔网格邻居分布示意图Fig.12 Schematic diagram of anisotropic Cartesian mesh neighborhood distribution

在网格自适应过程中，先依据给定的网格尺寸构造初始粗网格，以粗网格为基础根据需要进行加密。为了减少加密边界的复杂程度，首先将网格单元进行一定程度的各向同性加密；然后在物面曲率变化较大的区域进行各向异性加密，若单元对应的物面斜率在i方向上分量和临近单元在相同方向上的斜率分量相差大于设定阈值，则此单元在i方向上进行加密。至于流场解自适应判据，若继续采用各向同性笛卡尔网格解自适应判据会出现失效的情况。因此，采用基于速度旋度和散度以及压力梯度的综合判据，以实现不同流场结构(激波、漩涡等)的细致捕捉。

在进行流场自适应之前，首先需要解决插值模板点的问题。由于各向异性网格邻居排列组合形式多样，为了提高结果精度，插值模板点的选取需要对特定情况进行特殊处理。以图13为例，当网格i为Type=1的Case 1情况时，i+1模板点的流场值记为其所在单元的值，而Case 2、Case 3以及Case4的i+1模板点的流场值需要周围网格单元插值得到。图14为球柱流场解自适应结果，图例中的ρ代表无量纲的密度分布。从图中可以看出，各向异性笛卡尔网格可较好捕捉激波和尾流涡结构，同时，在激波附近，网格分布合理，过渡均匀。网格数目较各向同性笛卡尔网格可降低15%～30%，内存占用下降10%～25%，网格生成效率提高20%～30%。

图13 部分插值模板点选取示意图Fig.13 Schematic diagram of partial interpolation template points selection

图14 各向异性流场解自适应结果Fig.14 Adaptive results of anisotropic flow field solutions

2.3 法向射线加密

针对传统自适应笛卡尔网格方法在面对高雷诺数黏性流动问题时，存在网格数目过大的问题，Ruffin等[44]借鉴贴体结构网格大纵横比的特点，通过在壁面配置射线来确定加密区域，射线间则用粗网格过渡，构建了一种拟结构化的网格生成方法，即法向射线加密技术，提供了一种在可接受的网格数量下模拟高雷诺数黏性流动的笛卡尔方法。实际上，在边界层法向上布置高分辨率网格相对于沿壁面方向布置同等密度的网格会拥有更高的计算收益。法向射线加密技术则是基于该思想，在法向射线上布置足够细密的网格单元，在沿壁面方向就可以使用比较粗的网格。本文作者团队对于该方法进行了发展和实现，如图15所示，采用法向射线技术相应的网格数目比传统在边界层均匀加密的笛卡尔网格数目少很多，其网格总体分布与与贴体结构网格类似，利于在可接受的网格规模下捕捉边界层特征。

图15 自适应下笛卡尔网格和法向射线比较Fig.15 Adaptive Cartesian grid and normal ray comparison

法向射线加密技术分为3部分，分别是射线区域划分、射线间信息交互、计算插值与修正。射线区域的划分如图16所示，将整个计算域分为4个区域。4个区域中，蓝色和灰色区域并不参与流场求解。蓝色网格的流场值通过射线插值得到，通过应用射线间信息交互技术，构建了基于指针方法的相邻射线数据交互方法，并根据这些流场信息可靠度高的射线单元进行流场计算，射线间的粗网格的流场信息基于插值方法获得，实现了在较少法向射线的情况下高精度模拟边界层流动，大幅度降低了网格规模。

图16 法向射线加密技术的计算区域划分Fig.16 Region division of normal ray refinement technology

3 黏性物面边界处理

由于笛卡尔网格的非贴体特性，与物面相交会产生大小不同、形状各异的切割单元。与传统的贴体结构网格相比，这种单元形式不利于边界层流动特征的捕捉，容易引起边界层流动数据的非物理振荡。为了减小这种数值误差，若通过不计成本的加密方式，还会带来网格规模过大的问题。因此，发展了两种处理方法。

一种方法是采用重叠网格技术[62]，将贴体结构网格和笛卡尔网格的优势相结合，在物面边界附近提前生成贴体的结构或非结构网格，而在流场其他区域生成笛卡尔网格。通过两套网格的衔接，将笛卡尔网格的非贴体物面边界问题转移到两套网格的交界面处。

另一种方法是在笛卡尔网格框架下引入通过边界层近似理论得到的壁面函数方法，以降低数值方法对网格尺度的依赖、减小边界层内计算网格数量。

3.1 重叠网格技术

基于重叠网格开展数值计算需要网格间的数据传递和交互，在由粗网格贡献单元向细网格插值(或者双方长细比相差很大)时，都会引入较大数值误差。因此贡献单元与插值点所在网格单元应尽量保持一致的形状和尺度。基于笛卡尔网格的自适应优势，本文作者团队发展了重叠区尺度匹配的笛卡尔网格自适应技术，如图17所示。在初始笛卡尔网格的基础上，依据背景网格与贴体网格的尺寸特征进行自适应加密处理。通过自适应匹配，保证边界附近网格尺度的一致性，使流场在混合网格边界处光滑过渡，避免由于背景直角坐标网格与贴体结构网格在边界处的网格尺度相差太大而引起数值耗散。

图17 重叠网格交界面自适应笛卡尔网格Fig.17 Adaptive Cartesian mesh at hybrid mesh interfaces

贡献单元搜索[63]即“寻点”，是确定网格点在其他网格中的位置、查找网格贡献单元的统称。在网格的生成过程中需要进行大量的“寻点”操作，准确、高效的贡献单元搜索方法是重叠网格的关键技术之一。对笛卡尔背景网格和近壁贴体网格的重叠边界贡献单元搜索过程中，可结合笛卡尔网格正交规则的几何特征这一优势，实现搜索的高效、准确以及存储的低需求。笛卡尔网格贡献单元搜索如图18所示，以目标笛卡尔网格格心为中心，以网格的空间特征尺寸为参考，构建空间大小合适的包围盒，其中P1和P2表示包围盒两对角顶点；以空间坐标为衡量尺度，搜索该包围盒内所含有的结构网格节点；统计包围盒内结构网格节点数目，计算与目标笛卡尔网格格心的距离。若网格节点数满足贡献单元所需数目，依据距离参数选取所需结构网格节点；若网格节点数少于贡献单元所需数目，则继续扩大包围盒范围，直到可选贡献单元数目满足条件为止。以包围盒中所选结构网格节点为贡献单元，结合相应插值算法，获得笛卡尔网格流动参数。

图18 笛卡尔网格贡献单元搜索Fig.18 Cartesian grids donor cell search

通过上述发展的重叠区网格匹配技术和数据交互技术，与附面层贴体结构网格生成方法相结合[64]，构建了基于近壁贴体结构网格/笛卡尔网格的重叠网格技术。同时开展了真实管道列车的算例验证工作，初步考核了所发展的重叠笛卡尔网格方法的有效性和可靠性(图19)。

图19 光顺外形管道列车Fig.19 Smooth shape pipe train

3.2 壁面函数方法

传统的湍流壁面函数基于贴体网格提出，求解过程中直接应用壁面边界条件，这在非贴体网格下难以满足。非贴体网格的单元与物面的不一致性，将传统壁面函数直接推广至非贴体网格存在一定困难。

2002年2月，中国矿业大学北京东校区（其图书馆以下简称“北馆”）并入中国传媒大学。经过多年建设与深耕，中国传媒大学以新闻传播学、戏剧与影视学、信息与通信工程为龙头，逐渐形成文、工、管、经、法、理多学科协调发展的新闻传媒特色高校，为传媒行业的台前幕后培养专门人才，所有专业人才培养与课程设计都与传媒行业紧密结合。由于两馆原来服务的读者对象不同，藏书重点和收藏范围也不一样。笔者所在图书馆在这方面所做的规划和实践过程如下：

本文作者团队基于Spalding[65]提出的壁面函数模型，对该模型在笛卡尔网格下进行适用性修正。该模型格式较简单，且经过很多学者的考核验证，具有较高的可信度。针对可压缩湍流黏性流动问题的计算，为解决边界层非线性结构特征，引入虚拟点相对物面的参考点来取代无黏情况下的对称点。同时基于壁面函数修正湍流壁面边界条件，建立了面向笛卡尔网格的壁面函数方法，并开展了圆柱(图20，Cp为圆柱表面压力系数分布，θ为沿圆柱周向角度值)、翼型(图21)等典型二维算例的考核验证工作。通过计算结果对比，证明了采用壁面函数方法可以有效降低壁面网格尺度要求，并在网格较粗时达到良好的模拟精度。

图20 二维圆柱算例Fig.20 Two dimensional cylindrical calculation case

图21 NACA0012翼型算例Fig.21 NACA0012 airfoil calculation case

胡偶在其前期的工作基础上，继续对基于SSTk-ω湍流模型的Spalding壁面函数模型进行发展[23]。为了使参考点在壁面附近均匀分布以及提高壁面附近的湍流模拟精度和数值插值精度，提出了一种探针型浸入边界法(Ghost Cell Probe Immersed Boundary Method), 如图22所示，再结合针对SA和SST湍流模型的新型常微分化(ODE)壁面模型定义黏性壁面条件，形成耦合虚拟单元-壁面模型的湍流壁面边界条件重构方法。此外，除常用的Spalding等提出的壁面函数外，在现有浸入边界-虚拟单元方法体系下，耦合了Knopp[66]提出的自适应壁面函数模型(Adaptive Wall Function)，能够根据不同湍流模型自适应地调整经验公式，扩展了壁面函数模型的适用范围。针对三维问题，基于自适应混合笛卡尔网格体系进行了壁面函数理论研究，并通过典型的ONERA M6翼型绕流问题初步验证了理论方法，如图23所示，为后续移植到三维自适应笛卡尔网格体系做了技术积累。

图22 探针型浸入边界法示意图Fig.22 Schematic diagram of probe immersed boundary method

图23 ONERA M6 网格及绕流模拟结果Fig.23 ONERA M6 grid and flow simulation results

目前，针对非贴体网格下的边界条件的主要处理办法是，通过线性重构或插值方式获得壁面边界的流动变量。然而在湍流边界层内，不合理的线性重构得到的边界条件难以满足湍流情况下的边界层分布特性，容易诱发非物理的数值求解振荡。另一方面，常见的基于贴体网格发展的湍流模型(如SA、SST等)要求积分到壁面，对于湍流边界层内的网格具有较高的精度要求，并且在壁面边界及近壁区存在相匹配的近壁处理手段(Near-wall Treatment)，例如一方程SA模型使用壁面距离，k-ω模型使用Damping Function和Yap修正项等。在非贴体笛卡尔网格条件下，这也带来了匹配性问题，即在笛卡尔网格下特殊构造的壁面函数在参数形态上难以无缝对接原有成熟的湍流模型。这种匹配问题会带来计算误差，污染边界层求解精度。因此非贴体笛卡尔网格下鲁棒、高精度的湍流壁面函数是关键问题，有待进一步深入研究。

本文作者团队发展了一种适用于边界层全域的壁面函数。针对T3B平板边界层算例，如图24所示，其中New_wf为基于新型壁面函数方法得>到的计算结果，Full_grd为全网格下的摩阻计算结果，Exp代表实验参照结果。可以看出，基于SST模型，第1层网格高度在y+=50处；预测的摩阻系数Cf与y+(1)=0.01的全网格摩阻相当，具有较高的精度。后续将开展其在笛卡尔网格框架下的适用性应用研究。

图24 沿平板边界层流向摩擦阻力分布Fig.24 Friction drag distribution along flow direction of plate boundary layer

4 结 论

本文总结了国家数值风洞工程黏性自适应笛卡尔网格方法的研究进展，涵盖了网格生成、网格自适应方法和物面边界处理技术等内容，涉及南京航空航天大学、中南大学、中国空气动力研究与发展中心等单位在笛卡尔网格技术优化、非贴体笛卡尔网格物面边界处理等方面深入开展的工作，主要得出了以下结论：

1) 网格生成技术方面，分别在网格数据结构、物面单元检索技术、网格单元类型判断方法、高性能并行计算技术等开展了应用或改进工作，提高了三维自适应笛卡尔网格生成的效率和鲁棒性，为开展大规模网格量下复杂构型黏性流动问题的模拟研究奠定了基础。

2) 网格自适应技术方面，分别对各向同性自适应方法、各向异性自适应方法、法向射线加密方法进行了实现，构建了高效可靠的自适应判据，可优化网格分布，并有效降低网格规模，为复杂外形/流动问题的高效、高精度模拟丰富了技术手段。

3) 物面边界处理方面，分别对贴体式的重叠笛卡尔网格方法、非贴体式的浸入边界方法进行了研究与发展，提供了黏性物面边界模拟的不同技术手段，通过壁面函数方法的引入和适用性改进，初步形成了可有效降低湍流模拟的网格规模、保证非贴体黏性物面模拟精度的技术方法，为基于自适应笛卡尔网格方法开展高雷诺数黏性流动高效、高精度的模拟提供了有力的技术支撑。

本文对黏性自适应笛卡尔网格方法的国内外研究进展进行了简要梳理，重点对国家数值风洞工程黏性笛卡尔网格技术课题进展进行了阶段性总结介绍。在后续工作中，将在国家数值风洞工程的支撑下，直面黏性笛卡尔网格技术的难点和痛点，继续开展笛卡尔网格方法基础理论的深化研究和工程推广应用，为国家数值风洞工程丰富复杂构型和复杂流动模拟能力提供技术储备。