EI、Scopus 收录
中文核心期刊

基于自适应法向射线加密的笛卡尔网格流体仿真方法

邓露, 宋伦继, 罗灿炎, 高鹤, 毕林

邓露, 宋伦继, 罗灿炎, 高鹤, 毕林. 基于自适应法向射线加密的笛卡尔网格流体仿真方法. 力学学报, 2025, 57(6): 1-10. DOI: 10.6052/0459-1879-25-073
引用本文: 邓露, 宋伦继, 罗灿炎, 高鹤, 毕林. 基于自适应法向射线加密的笛卡尔网格流体仿真方法. 力学学报, 2025, 57(6): 1-10. DOI: 10.6052/0459-1879-25-073
Deng Lu, Song Lunji, Luo Canyan, Gao He, Bi Lin. Cartesian grid fluid simulation method based on adaptive normal ray refinement. Chinese Journal of Theoretical and Applied Mechanics, 2025, 57(6): 1-10. DOI: 10.6052/0459-1879-25-073
Citation: Deng Lu, Song Lunji, Luo Canyan, Gao He, Bi Lin. Cartesian grid fluid simulation method based on adaptive normal ray refinement. Chinese Journal of Theoretical and Applied Mechanics, 2025, 57(6): 1-10. DOI: 10.6052/0459-1879-25-073
邓露, 宋伦继, 罗灿炎, 高鹤, 毕林. 基于自适应法向射线加密的笛卡尔网格流体仿真方法. 力学学报, 2025, 57(6): 1-10. CSTR: 32045.14.0459-1879-25-073
引用本文: 邓露, 宋伦继, 罗灿炎, 高鹤, 毕林. 基于自适应法向射线加密的笛卡尔网格流体仿真方法. 力学学报, 2025, 57(6): 1-10. CSTR: 32045.14.0459-1879-25-073
Deng Lu, Song Lunji, Luo Canyan, Gao He, Bi Lin. Cartesian grid fluid simulation method based on adaptive normal ray refinement. Chinese Journal of Theoretical and Applied Mechanics, 2025, 57(6): 1-10. CSTR: 32045.14.0459-1879-25-073
Citation: Deng Lu, Song Lunji, Luo Canyan, Gao He, Bi Lin. Cartesian grid fluid simulation method based on adaptive normal ray refinement. Chinese Journal of Theoretical and Applied Mechanics, 2025, 57(6): 1-10. CSTR: 32045.14.0459-1879-25-073

基于自适应法向射线加密的笛卡尔网格流体仿真方法

详细信息
    通讯作者:

    毕林, 研究员, 主要研究方向为高空复杂流动机理研究与应用. E-mail: bzbaby1010@163.com

  • 中图分类号: O354.1

CARTESIAN GRID FLUID SIMULATION METHOD BASED ON ADAPTIVE NORMAL RAY REFINEMENT

  • 摘要: 在黏性流动模拟中, 近壁流动在物面法向方向上的速度梯度远大于切线方向, 呈现出显著的各向异性特征. 传统的各向同性笛卡尔网格方法在捕捉边界层流动细节时面临网格数量剧增和计算效率下降的挑战. 针对这一问题, 提出了一种基于自适应法向射线加密(adaptive normal ray refinement, ANRR)的笛卡尔网格黏性流体仿真方法. 该方法的核心要义在于, 根据物面切线方向上的角度变化程度来自适应生成法向射线种子点, 在法向射线附近进行网格加密, 以精确捕捉边界层流动特征, 同时在射线之间采用较粗糙的网格过渡, 从而在保证计算精度的前提下有效减少整体网格数量. 然后利用插值构建了射线间的高效信息传递技术, 确保流场求解过程中的准确性. 最后, 对层流平板、低雷诺数圆柱绕流以及NACA0012翼型绕流等典型算例进行模型验证. 结果表明, 与传统几何自适应加密网格方法相比, ANRR网格方法在边界层流动区域显著减少了网格规模, 在保持高精度的同时提升了计算效率, 为自适应笛卡尔网格的黏性流动问题高效求解提供了新的解决方案.
    Abstract: In viscous flow simulations, the velocity gradient in the normal direction of the wall is much larger than that in the tangential direction, presenting a significant anisotropic feature. Traditional isotropic Cartesian grid methods face the challenges of a sharp increase in the number of grids and a decline in computational efficiency when capturing the details of boundary layer flows. To address this issue, this paper proposes a Cartesian grid viscous fluid simulation method based on Adaptive Normal Ray Refinement (ANRR). The core of this method lies in adaptively generating normal ray seed points according to the degree of angle change in the tangential direction of the surface, and performing grid refinement near the normal rays to accurately capture the characteristics of boundary layer flows. Meanwhile, coarser grids are used for transition between rays, effectively reducing the overall number of grids while maintaining computational accuracy. Then, an efficient information transfer technology between rays is constructed through interpolation to ensure the accuracy of the flow field solution process. Finally, typical cases such as laminar flow over a flat plate, low Reynolds number flow around a circular cylinder, and flow around an NACA0012 airfoil are used for model verification. The results show that compared with the traditional geometric adaptive grid refinement method, the ANRR grid method significantly reduces the grid scale in the boundary layer flow region, improves computational efficiency while maintaining high accuracy, and provides a new solution for the efficient solution of viscous flow problems using adaptive Cartesian grids.
  • 近年来, 基于笛卡尔网格的流体仿真方法发展迅猛, 已成为计算流体力学(computational fluid dynamics, CFD)在工程领域中应用的重要工具[1-5]. 笛卡尔网格是一种基于笛卡尔坐标系的网格系统, 其所有网格单元的边或面与笛卡尔坐标系的各轴相互正交. 这种网格系统能够以非贴体的方式生成空间网格, 对物体表面的依赖程度较低, 易于实现网格的自动生成. 另外, 笛卡尔网格通常采用叉树结构, 这种独特的树状存储方式天然适用于动态自适应过程[6-10]. 有学者认为[11], 笛卡尔网格方法是实现CFD第三次技术突破的关键因素之一.

    然而, 对于黏性流动, 尤其是高雷诺数湍流问题, 物面法向方向上流动梯度远大于切线方向, 呈现各向异性特征, 但笛卡尔网格因其各向同性特点, 需要较大规模的网格来求解边界层流动. 再者对于湍流问题, 物面附近的第一层网格尺度要求Δy + < 1, 带来“网格灾难”问题, 是制约笛卡尔网格工程实用的主要问题之一[12-13].

    针对上述问题, 国际上尝试采用不同的方法进行研究[14]. 比如各向异性自适应笛卡尔网格, Wang等[15]提出将流场变量梯度或者二阶导数作为各向异性自适应的判据, 实现了在各向异性特征明显的区域进行特定方向的加密或粗化. 还有许多学者[16-18]发展了各向异性自适应笛卡尔网格方法用于减少网格量. 但是各向异性笛卡尔网格下流场变量值的计算将变得复杂, 且面临计算鲁棒性问题. 也有学者采用外部为笛卡尔网格, 物面附近为结构/非结构网格的混合/重叠网格技术[19-26]来规避这个问题, 但这种方法需要多套网格之间的信息交互, 带来额外的插值误差, 且破坏了纯粹笛卡尔网格自动化及自适应优势. 针对湍流问题, 有学者采用壁面函数方法[27-31], 将壁面附近第一层网格尺度放宽到对数律区, 显著减少了计算资源的需求, 提高了模拟的效率. 然而, 目前的壁面函数方法多用于平衡湍流模拟, 对于非平衡湍流模拟能力欠缺. 此外, Ruffin等[32]提出了一种法向射线加密(normal ray refinement, NRR)方法. NRR方法间歇地在物面放置法向射线, 每一条法向射线附近聚集了很多精细的网格, 而在两条法向射线之间则由粗网格隔开, 这样可以大幅减少附面层流动模拟的网格需求. 在此基础上, Zaki等[33]在NRR求解迭代过程中对法向射线长度自适应进行了测试, Arslanbekov等[34]提出一种新的法向射线放置方法, Bopp等[35]实现了NRR方法在三维流动模拟中的应用.

    本文在Ruffin等[32]的基础上, 拟依据物面切线角度变化自动生成法向射线种子点, 实现法向射线几何自适应放置, 发展射线间通信技术, 然后建立起一种基于自适应法向射线加密(adaptive normal ray refinement, ANRR)的笛卡尔网格流体仿真方法. 并用层流平板、低雷诺数圆柱绕流以及NACA0012翼型绕流等典型算例进行考核验证. 以期为自适应笛卡尔网格的黏性流动问题高效求解提供新思路.

    本文使用有限差分法来求解直角坐标系下(横轴为x轴, 纵轴为y轴)的二维可压缩Navier-Stokes方程. 通过无量纲化处理, 最终得到的守恒型方程可表示为以下形式

    $$ \frac{{\partial {\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{E}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{F}}}}{{\partial y}} = \frac{{\partial {{\boldsymbol{E}}_{\text{v}}}}}{{\partial x}} + \frac{{\partial {{\boldsymbol{F}}_{\text{v}}}}}{{\partial y}} $$ (1)

    式中, U为守恒通量, EF为无黏流项; EvFv为黏性流项. 在后续数值算例计算中, 采用浸入边界-虚拟单元法处理物面边界[9]; 空间离散对流项采用Roe-MUSCL 三阶(Roe近似Riemann解为代表的通量差分分裂方法[36], MUSCL: monotonic upstream-centered scheme for conservation[37]) 5点模板格式离散, 黏性项采用二阶中心差分格式; 时间推进采用显式Runge-Kutta方法配合当地时间步长离散.

    笛卡尔网格得益于叉树结构, 叶子节点的可扩展性使得网格自适应易于实现. 叉树数据结构本质上是一种分层数据结构, 各层之间的划分依据是网格的层级大小. 在二维情况下对应四叉树结构, 如图1所示. 自适应笛卡尔网格的初始均匀网格由层级为0的网格单元构成, 是叉树结构的根节点, 并且这些单元没有父单元. 没有子单元的网格单元被称为叶子单元, 其具体层级取决于加密次数. 在各向同性的笛卡尔网格中, 每次加密网格时, 子单元的边长会减半, 同时网格的层级增加1. 子单元的位置排列顺序遵循特定的规则(逆时针方向)进行编号(son1 ~ son4).

    图  1  二维四叉树数据结构及网格图
    Figure  1.  Two dimensional quadtree data structure and grid diagram

    在ANRR框架中, 具体的网格生成步骤如下.

    步骤1: 读入外形数据文件;

    步骤2: 根据计算域大小和初始网格间距h生成均匀的初始背景网格;

    步骤3: 将网格根据是否与物面相交分成流场中的外部网格、相交网格、以及物面内部网格;

    步骤4: 根据需要进行块加密, 在某一计算域加密所有网格到提前设定好的层级, 并进行光顺处理, 使得相邻网格的层次之差不能超过1;

    步骤5: 根据物面切线角度变化, 绕物体表面自适应放置法向射线种子点, 形成自适应法向射线加密网格, 具体方法见第3节.

    图2所示, 物面内部单元为白色无填充区域, 法向射线沿着物面在外部网格和相交网格中放置. 具体来说, 物面会生成法向射线种子点, 从种子点出发, 法向射线延伸到预定的长度(与几何自适应加密网格方法的最大加密层级法向方向上所有网格总高度一致); 与法向射线相交的初始网格会加密到提前设定好的网格层级; 同时, 加密后的网格会拓宽到一定程度, 以便于后续插值的使用. 即, 物面法向射线沿伸的地方放置非常精细的网格, 在两条法向射线之间则由较为粗的网格组成. 这也符合边界层流动在法向上较切向上变化显著的规律.

    图  2  法向射线加密
    Figure  2.  Normal ray refinement

    我们根据物面外形设计了自动生成法向射线种子点. 生成方法如下.

    判据1: 如图3所示, 将黑色线段视为一部分物面(示意图, 物面可以是曲折的, 即黑色线段可以是黑色曲线); 固定一个长度记为Δs, 从A点出发沿着绿色线段到达B点的长度为Δs; 从A点到B点物面切线变化的角度记为Δθ. THR是一个给定的阈值, 如果ΔθTHR, 且Δθ≤90°, 则在B点生成法向射线种子点.

    图  3  B点生成种子点
    Figure  3.  Generate seed points at point B

    判据2: 沿用上述的定义, 如果ΔθTHR, 将Δs乘以一个给定的系数Coef, 绿色线段延伸到C点, 然后在C点生成法向射线种子点, 如图4所示.

    图  4  C点生成种子点
    Figure  4.  Generate seed points at point C

    判据3: 物面外形是由点点连线构成的. 如图5所示, 当3个外形点(绿色小圆)形成了尖角. 即当Δθ≥90°, 分别在尖角两边(蓝色线段)上的绿色小圆处生成法向射线种子点(红色虚线为延伸出的法线); 且在角尖处生成法向射线种子点(紫色虚线为延伸出的法线), 它的方向是两条蓝色线段上的绿色小圆连接而成的黑色线段的法向. 图6中红色小圆处就是判据三法向射线加密的情况示意.

    图  5  尖角处生成种子点
    Figure  5.  Generate seed points at sharp corners
    图  6  尖角处法向射线加密
    Figure  6.  Refinement of normal ray at sharp corners

    通过判据1和判据2可知, 在物面曲率大的位置, 所生成的法向射线间隔较短; 而在曲率小的位置, 所生成的法向射线间隔较长. 图7展示了这种自适应法向射线生成特征.

    图  7  自适应法向射线生成特征
    Figure  7.  Adaptive normal ray generation features

    定义了3种类型的计算区域[32], 如图8所示. 在图8中, 蓝色区域表示法线高度外的流体网格, 红色区域表示法向射线上的精细网格, 绿色区域表示法线高度内的粗网格, 位于两条法向射线之间. 为了将精细网格流场信息从一条法向射线传输到下一条, 实现法向射线间通信, 对绿色单元定义了指针(PTR). 图8中, 对于黑色三角所代表的绿色网格单元来说, 它的指针#1指向距离最近的法向射线上的精细网格a(后续会利用到指向网格的流场信息); 指针#2指向第二近的法向射线上的精细网格b; 指针#3指向网格ab构成的射线上的网格c. 同时, 网格bc在同一条法向射线上存在间隔. 这些指向的精细单元(图8中橙色圆点)与黑色三角所代表的绿色单元距离物面的距离是基本相同的. 同时, 为了保证插值的准确性, 指向的精细单元不在法向射线边缘. 特殊情况, 指针指向的位置可能位于不同精细单元共同的边上, 如图9所示, 黑色三角为被插值单元, 黑色圆点是指向位置; 这时取周围最近的4个单元的流场信息的均值替代指针指向的精细单元的流场信息, 即蓝色正方形框内4个单元.

    图  8  插值指针
    Figure  8.  Interpolation pointer
    图  9  指针指向的特殊情况
    Figure  9.  The special case pointed to by the pointer

    图10所示是计算程序所使用的网格计算模板, 即使出现悬挂网格(计算中心网格周围的网格层级水平和本身不同), 仍将周围的数据利用最小二乘法插值到计算模板中. 对于上述3种计算区域来说, 计算方法如下.

    图  10  计算模板
    Figure  10.  Calculation template

    (1) 对于所有位于法线高度外的网格(蓝色单元)使用有限差分法;

    (2) 粗网格(绿色单元)的流场信息用指针指向的两个相邻法向射线上的3个精细单元的流场信息进行插值替换. 精细网格单元的计算模板延伸到粗网格区域中和绿色单元重叠的位置被称为NRR虚拟单元[13], 即图11中绿色单元, 它的流场信息同样被插值替换.

    图  11  NRR虚拟单元
    Figure  11.  NRR ghost cell

    (3) 因为黏性区域中的红色精细网格单元足够小, 直接利用计算模板中的数据使用有限差分法求解.

    无论是绿色粗网格区域的插值还是虚拟单元的插值, 都使用同样的插值公式. 插值点的状态矢量可以通过二项式插值方法[32]确定

    $$ {{{\boldsymbol{V}}}}''_{P} = {{\boldsymbol{C}}}_{1} + {{\boldsymbol{C}}}_{2}{d}_{cP} + {{\boldsymbol{C}}}_{3}{d}_{cP}^{2} $$ (2)

    其中

    $$ {{\boldsymbol{C}}_1} = {{\boldsymbol{V}} _c} $$ (3)
    $$ {{\boldsymbol{C}}_2} = \frac{{{{{\boldsymbol{V}} }_b} - {{{\boldsymbol{V}} }_c}}}{{{d_{cb}}}} $$ (4)
    $$ {{\boldsymbol{C}}_3} = \frac{{{{{\boldsymbol{V}} }_a} - {{\boldsymbol{C}}_1} - {{\boldsymbol{C}}_2}{d_{ca}}}}{{d_{ca}^2}} $$ (5)

    式中, 见图8, P单元为被插值单元(黑色三角), abc单元分别是指针#1、#2和#3指向的单元; ${\boldsymbol{V}} $是某个流场信息(速度、压力和密度等); d为距离, dcP则表示c单元和P单元之间的距离, 其余类似.

    为增强稳定性, 不允许在指针指向的精细单元ab之间引入新的最大值或最小值. 因此, 对单元P处的插值状态矢量使用以下表达式

    $$ {{{\boldsymbol{V}}}}_{p} = \left\{\begin{aligned} &{{{\boldsymbol{V}}}}''_{P}\;\; {\mathrm{if}}\;\; {{{\boldsymbol{V}}}}_{P\mathrm{min}} < {{{\boldsymbol{V}}}}''_{P} < {{{\boldsymbol{V}}}}_{P\mathrm{max}} \\ &{{{\boldsymbol{V}}}}'_{P}\;\; {\mathrm{if}}\;\; {{{\boldsymbol{V}}}}''_{P}\leqslant {{{\boldsymbol{V}}}}_{P\mathrm{min}}\;\; {\mathrm{or}}\;\; {{{\boldsymbol{V}}}}''_{P}\geqslant {{{\boldsymbol{V}}}}_{P\mathrm{max}}\end{aligned}\right. $$ (6)

    其中

    $$ {{\boldsymbol{V}}_{P\max }} = \max \left( {{{\boldsymbol{V}}_a},{{\boldsymbol{V}}_b}} \right) $$ (7)
    $$ {{\boldsymbol{V}}_{P\min }} = \min \left( {{{\boldsymbol{V}}_a},{{\boldsymbol{V}}_b}} \right) $$ (8)
    $$ {\boldsymbol{V}} '_P = \dfrac{{\dfrac{{{{{\boldsymbol{V}} }_a}}}{{{d_{Pa}}}} + \dfrac{{{{{\boldsymbol{V}} }_b}}}{{{d_{Pb}}}}}}{{\dfrac{1}{{{d_{Pa}}}} + \dfrac{1}{{{d_{Pb}}}}}} $$ (9)

    首先采用自适应法向射线加密网格方法求解层流平板边界层. 基于来流参数的单位雷诺数Re = 105, 马赫数Ma = 0.2. 入口为自由来流, 出口为压力出口, 顶部为自由来流边界条件. 计算域大小设置为[−0.48,1.28] × [−0.08,1.20], 其中平板底面x∈[−0.48,0.0]为对称滑移壁面, x∈[0.0,1.28]为无滑移壁面. 初始网格尺寸为h = 0.04, 设置无滑移壁面最大加密层数为7层, 滑移壁面最大加密层数为6层. 几何自适应加密网格方法[9]的网格数为134320, 网格图如图12(a)所示; 自适应法向射线加密网格方法的网格数为82324, 网格图如图12(b)所示.

    图  12  层流平板算例的网格图对比
    Figure  12.  Comparison of grid diagrams for laminar flow flat plate examples

    图13所示, 对比了自适应法向射线加密(ANRR)网格方法和传统几何自适应(geometric adaptive, GeoA)加密网格方法以及Blasius平板解在x = 0.565处的切向速度(U)分布的实验结果. 在图中可以看到, ANRR网格方法得到的实验结果和Blasius平板解一致. 为了进一步验证ANRR网格框架计算的准确性, 还对表面摩擦阻力系数进行了对比, 如图14所示, 可以看到三者吻合良好.

    图  13  层流平板算例在x = 0.565处的U速度对比
    Figure  13.  Comparison of U-velocity at x = 0.565 for laminar flow flat plate example
    图  14  层流平板算例表面摩擦系数对比
    Figure  14.  Comparison of surface friction coefficients for laminar flow flat plate examples

    计算条件设置如下: 圆径D = 1, 来流马赫数Ma = 0.2. 计算域大小设置为[−24,64] × [−34,34]. 初始网格尺寸为h = 2, 设置最大加密层数为7层. 为了捕捉尾部产生的周期漩涡, 在x∈[−3,15], y∈[−3,3]设置了一个5层加密的加密块(注: 圆柱中心位于(0,0)). 几何自适应的流场网格依据到壁面的距离逐渐加密, 网格数量为58244, 网格图如图15(a)所示; 自适应法向射线的流场网格数为48592, 网格图如图15(b)所示. 除来流面设置为速度入口边界条件外, 计算域的其他面设置为远场无反射边界条件, 圆柱壁面周围采用无滑移壁面边界条件.

    图  15  圆柱绕流算例的网格图对比
    Figure  15.  Comparison of grid diagrams for flow around a cylinder calculation examples

    利用ANRR方法分别计算了基于直径的雷诺数Re = 10, 20, 40, 60, 80和100的流动. 在图16中可以看到Re≤60的Ma流场和流线图, 观察其流动特征, 在低雷诺数下(Re≤40时), 流体平行且有序地流动, 流动几乎保持在层流状态. 流线沿着圆柱表面流动, 在圆柱周围, 流线呈现出对称的特征. 在圆柱的后侧, 流动沿表面贴附, 不会立即分离, 形成一个尾迹区, 随着雷诺数的增加, 在尾迹区的漩涡逐渐被拉伸. 当Re = 60时, 流动开始不稳定, 流线不再是平行且规则的, 而是出现了涡流分离和再附现象, 尾迹摆动明显.

    图  16  Re≤60的圆柱绕流流场和流线图
    Figure  16.  Flow field and streamline diagram around a cylinder with Re≤60

    对于Re = 100的流动, 将几何自适应加密网格和法向射线加密网格计算的流场在推进过程中进行了对比, 如图17所示. 流场速度云图分布并无太大区别.

    图  17  Re = 100时圆柱绕流算例的速度云图对比
    Figure  17.  Comparison of velocity cloud diagrams for flow around a cylinder with Re = 100

    将ANRR方法计算的不同雷诺数下的平均阻力以及S数(非定常时间变化尺度的无量纲数)与Tritton[38]和Zhong等[39]的计算结果进行了对比. 表1数据是Re = 10, 20和40的平均阻力系数的比较. 图18分别展示了Re = 80和100时ANRR方法和GeoA方法计算所得升力系数的对比, 图中所示, 升力系数呈正弦波动, 幅值和频率彼此接近, 但由于二者收敛速度不同, 图像在初期的升力曲线不重合. 随着雷诺数的变大, 升力系数也逐渐增加, 具有周期性, 也可据此计算St数. 表2数据是Re = 60, 80和100的平均阻力系数(Cda)和St数的比较. 总体而言, 本文的数值结果与参考文献中的实验结果呈现出较好的一致性, 表明了ANRR方法流场求解的准确性.

    表  1  Re = 10, 20和40时的圆柱绕流平均阻力系数的比较
    Table  1.  Comparison of average drag coefficients of flow around a cylinder at Re = 10, 20 and 40
    Re Tritton Zhong et al GeoA ANRR
    10 2.926 3.015 2.964 2.996
    20 2.103 2.115 2.131 2.144
    40 1.605 1.557 1.589 1.595
    下载: 导出CSV 
    | 显示表格
    图  18  Re = 80和100时圆柱绕流升力系数对比
    Figure  18.  Comparison of lift coefficients around flow in a cylinder with Re = 80 and 100
    表  2  Re = 60, 80和100时的圆柱绕流平均阻力系数以及St数的比较
    Table  2.  Comparison of the average drag coefficient and St number of flow around a cylinder around Re = 60, 80 and 100
    Re Tritton Zhong et al. GeoA ANRR
    Cda Cda St Cda St Cda St
    60 1.398 1.423 0.137 7 1.392 0.126 7 1.450 0.135 0
    80 1.316 1.375 0.152 3 1.403 0.154 5 1.391 0.151 2
    100 1.271 1.349 0.165 3 1.376 0.167 2 1.369 0.162 8
    下载: 导出CSV 
    | 显示表格

    本小节计算NACA0012翼型上的低雷诺数黏性流动, 机翼弦长C为1. 设置流动参数: Ma = 0.8, α = 10°(攻角), Re = 500; 计算域大小设置为[−8,12] × [−10,10], 初始网格尺寸为h = 0.1, 最大加密层数为6层. 在x∈[−0.2,1.2], y∈[−0.3,0.3]区域设置了一个3层加密的加密块. 几何自适应加密网格方法网格如图19(a)所示, 网格数为116552; 自适应法向射线加密网格如图19(b)所示, 网格数为71 488.

    图  19  NACA0012算例的网格图对比
    Figure  19.  Comparison of grid diagrams for NACA0012 calculation example

    将GeoA和ANRR网格计算的速度流场进行比较, 如图20所示. 图中可以看到, 流场和流场线分布基本相同.

    图  20  NACA0012算例的速度云图对比
    Figure  20.  Comparison of velocity cloud diagrams for NACA0012 calculation example

    同时, 将计算数据和Jawahar等[40]的非结构化网格及其论文中提到的结构三角网格(STRI-V), 以及GAMM研讨会中记录的可压缩Navier-Stokes计算结果[41]进行了对比. 如图21所示, 比较了ANRR方法和Jawahar等得到的解, 展示了表面压力系数的对比, 可以看到表面压力系数吻合良好, ANRR与参考文献相比, 在x = 0.1的上表面附近出现最大偏差, 但不超过7%.

    图  21  NACA0012算例表面压力系数对比
    Figure  21.  Comparison of surface pressure coefficient for NACA0012 calculation example

    表3展示了本文计算数据和参考文献压差阻力系数(Cdp)、摩擦阻力系数(Cdf)、总阻力系数(Cd)和总升力系数(Cl)的对比, 它们的数值彼此接近, 并且完全在GAMM研讨会报告的范围内.

    表  3  NACA0012翼型算例的升力和阻力系数的比较
    Table  3.  Comparison of lift and drag coefficients for NACA0012 airfoil example
    Parameter Cdp Cdf Cd Cl
    Jawahar et al. 0.15287 0.12439 0.27726 0.50231
    STRI-V 0.14930 0.12286 0.27216 0.49394
    GeoA 0.14981 0.12318 0.27299 0.46503
    ANRR 0.15007 0.12238 0.27245 0.45484
    GAMM 0.243 ~ 0.286 8 0.414 5 ~ 0.517
    下载: 导出CSV 
    | 显示表格

    在本文中, 开发和测试了一种基于自适应法向射线加密的笛卡尔网格流体仿真方法. 利用物面切线角度变化的3个判据, 自适应生成了法向射线种子点, 以此生成的法向射线加密网格具有几何适应性. 在此基础上搭建了法向射线加密网格框架, 构造了射线间通信方法以求解Navier-Stokes方程. 为了验证ANRR框架下的流场求解器的可行性, 通过层流平板、低雷诺数圆柱绕流以及NACA0012翼型绕流等典型算例测试了该方法. 与几何自适应加密网格方法相比, 笛卡尔网格的数量可以大幅减少, 并且得到了准确鲁棒性的数值结果, 与参考文献吻合. 下一步, 拟将该方法推广到三维, 实现其工程实用化.

  • 图  1   二维四叉树数据结构及网格图

    Figure  1.   Two dimensional quadtree data structure and grid diagram

    图  2   法向射线加密

    Figure  2.   Normal ray refinement

    图  3   B点生成种子点

    Figure  3.   Generate seed points at point B

    图  4   C点生成种子点

    Figure  4.   Generate seed points at point C

    图  5   尖角处生成种子点

    Figure  5.   Generate seed points at sharp corners

    图  6   尖角处法向射线加密

    Figure  6.   Refinement of normal ray at sharp corners

    图  7   自适应法向射线生成特征

    Figure  7.   Adaptive normal ray generation features

    图  8   插值指针

    Figure  8.   Interpolation pointer

    图  9   指针指向的特殊情况

    Figure  9.   The special case pointed to by the pointer

    图  10   计算模板

    Figure  10.   Calculation template

    图  11   NRR虚拟单元

    Figure  11.   NRR ghost cell

    图  12   层流平板算例的网格图对比

    Figure  12.   Comparison of grid diagrams for laminar flow flat plate examples

    图  13   层流平板算例在x = 0.565处的U速度对比

    Figure  13.   Comparison of U-velocity at x = 0.565 for laminar flow flat plate example

    图  14   层流平板算例表面摩擦系数对比

    Figure  14.   Comparison of surface friction coefficients for laminar flow flat plate examples

    图  15   圆柱绕流算例的网格图对比

    Figure  15.   Comparison of grid diagrams for flow around a cylinder calculation examples

    图  16   Re≤60的圆柱绕流流场和流线图

    Figure  16.   Flow field and streamline diagram around a cylinder with Re≤60

    图  17   Re = 100时圆柱绕流算例的速度云图对比

    Figure  17.   Comparison of velocity cloud diagrams for flow around a cylinder with Re = 100

    图  18   Re = 80和100时圆柱绕流升力系数对比

    Figure  18.   Comparison of lift coefficients around flow in a cylinder with Re = 80 and 100

    图  19   NACA0012算例的网格图对比

    Figure  19.   Comparison of grid diagrams for NACA0012 calculation example

    图  20   NACA0012算例的速度云图对比

    Figure  20.   Comparison of velocity cloud diagrams for NACA0012 calculation example

    图  21   NACA0012算例表面压力系数对比

    Figure  21.   Comparison of surface pressure coefficient for NACA0012 calculation example

    表  1   Re = 10, 20和40时的圆柱绕流平均阻力系数的比较

    Table  1   Comparison of average drag coefficients of flow around a cylinder at Re = 10, 20 and 40

    Re Tritton Zhong et al GeoA ANRR
    10 2.926 3.015 2.964 2.996
    20 2.103 2.115 2.131 2.144
    40 1.605 1.557 1.589 1.595
    下载: 导出CSV

    表  2   Re = 60, 80和100时的圆柱绕流平均阻力系数以及St数的比较

    Table  2   Comparison of the average drag coefficient and St number of flow around a cylinder around Re = 60, 80 and 100

    Re Tritton Zhong et al. GeoA ANRR
    Cda Cda St Cda St Cda St
    60 1.398 1.423 0.137 7 1.392 0.126 7 1.450 0.135 0
    80 1.316 1.375 0.152 3 1.403 0.154 5 1.391 0.151 2
    100 1.271 1.349 0.165 3 1.376 0.167 2 1.369 0.162 8
    下载: 导出CSV

    表  3   NACA0012翼型算例的升力和阻力系数的比较

    Table  3   Comparison of lift and drag coefficients for NACA0012 airfoil example

    Parameter Cdp Cdf Cd Cl
    Jawahar et al. 0.15287 0.12439 0.27726 0.50231
    STRI-V 0.14930 0.12286 0.27216 0.49394
    GeoA 0.14981 0.12318 0.27299 0.46503
    ANRR 0.15007 0.12238 0.27245 0.45484
    GAMM 0.243 ~ 0.286 8 0.414 5 ~ 0.517
    下载: 导出CSV
  • [1]

    Buning PG, Nichols RH. OVERFLOW2 training class-afternoon session//10th Symposium on Overset Composite Grids & Solution Technology, Moffett Field, 2010

    [2]

    Spurlock WM, Aftosmis MJ, Nemec M. Cartesian Mesh Simulations for the third AIAA sonic boom prediction workshop. Journal of Aircraft, 2021, 59(3): 708-724

    [3] 陈坚强, 袁先旭, 涂国华等. 计算流体力学2035愿景. 北京: 科学出版社, 2023 (Chen Jianqiang, Yuan Xianxu, Tu Guohua, et al. Computational Fluid Dynamics 2035 Vision in China. Beijing: Science Press, 2023 (in Chinese)

    Chen Jianqiang, Yuan Xianxu, Tu Guohua, et al. Computational Fluid Dynamics 2035 Vision in China. Beijing: Science Press, 2023 (in Chinese)

    [4]

    Anderson GR, Aftosmis MJ, Nemec M. Cart3D simulations for the second AIAA sonic boom prediction workshop. Journal of Aircraft, 2019, 56(3): 896-911 doi: 10.2514/1.C034842

    [5]

    Mani M, Dorgan AJ. A perspective on the state of aerospace computational fluid dynamics technology. Annual Review of Fluid Mechanics, 2023, 55(1): 431-457 doi: 10.1146/annurev-fluid-120720-124800

    [6]

    Zeeuw De. A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations. [PhD Thesis]. Ann Arbor: University of Michigan, 1993

    [7]

    Lourenço M A, Padilla EL. An octree structured finite volume based solver. Applied Mathematics And Computation, 2020, 365: 1-28

    [8] 孟爽. 各向异性自适应笛卡尔网格方法研究及在高速列车中的应用. [博士论文]. 长沙: 中南大学, 2023 (Meng Shuang. Research on anisotropic adaptive Cartesian grid method and its application to high-speed trains. [PhD Thesis]. Changsha: Central South University, 2023 (in Chinese)

    Meng Shuang. Research on anisotropic adaptive Cartesian grid method and its application to high-speed trains. [PhD Thesis]. Changsha: Central South University, 2023 (in Chinese)

    [9] 罗灿炎. 浸入式笛卡尔网格数值方法及高速列车气动性能应用研究. [博士论文]. 长沙: 中南大学, 2024 (Luo Canyan. Immersed boundary Cartesian grid numerical method and its application in high-speed train aerodynamics. [PhD Thesis]. Changsha: Central South University, 2024 (in Chinese)

    Luo Canyan. Immersed boundary Cartesian grid numerical method and its application in high-speed train aerodynamics. [PhD Thesis]. Changsha: Central South University, 2024 (in Chinese)

    [10]

    Meng S, Zhou D, Yuan X, et al. Enhanced strategy for adaptive Cartesian grid generation with arbitrarily complex 3D geometry. Advances in Engineering Software, 2022, 174: 103304 doi: 10.1016/j.advengsoft.2022.103304

    [11]

    Nakahashi K. Aeronautical CFD in the age of petaflops-scale computing: From unstructured to Cartesian meshes. European Journal of Mechanics B/Fluids, 2013, 40: 75-86 doi: 10.1016/j.euromechflu.2013.02.005

    [12]

    Chawner JR, Dannenhoffer J, Taylor NJ. Geometry, mesh generation, and the CFD 2030 vision//46th AIAA Fluid Dynamics Conference, Washington, 2016

    [13]

    Verzicco R. Immersed boundary methods: Historical perspective and future outlook. Annual Review of Fluid Mechanics, 2023, 55(1): 129-155 doi: 10.1146/annurev-fluid-120720-022129

    [14] 赵宁, 刘剑明, 田琳琳等. 可压缩流动问题笛卡尔网格模拟方法研究进展与展望. 力学学报, 2025, 57(2): 285-314 (Zhao Ning, Liu Jianming, Tian Linlin, et al. Progress and prospects of cartesian mesh simulation methods for compressible flow problems. Chinese Journal of Theoretical and Applied Mechanics, 2025, 57(2): 285-314 (in Chinese) doi: 10.6052/0459-1879-24-374

    Zhao Ning, Liu Jianming, Tian Linlin, et al. Progress and prospects of cartesian mesh simulation methods for compressible flow problems. Chinese Journal of Theoretical and Applied Mechanics, 2025, 57(2): 285-314 (in Chinese) doi: 10.6052/0459-1879-24-374

    [15]

    Wang ZJ, Chen RF. Anisotropic solution-adaptive viscous Cartesian grid method for turbulent flow simulation. AIAA Journal, 2002, 40(10): 1969-1978

    [16]

    Li K, Wu ZN. Nonet-cartesian grid method for shock flow computations. Journal of Scientific Computing, 2003, 20(3): 303-329

    [17]

    Capizzano F. A compressible flow simulation system based on Cartesian grids with anisotropic refinements//45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, 2007

    [18]

    Keats WA, Lien FS. Two-dimensional anisotropic cartesian mesh adaptation for the compressible Euler equations. International Journal for Numerical Methods in Fluids, 2010, 46(11): 1099-1125

    [19]

    Steger JL, Dougherty FC, Benek JA. A chimera grid scheme//ASME Mini-Symposium on Advances in Grid Generation, Houston, 1982

    [20]

    Karman LJS. Splitflow-A 3D unstructured Cartesian/prismatic grid CFD code for complex geometries//33rd Aerospace Sciences Meeting and Exhibit, Reno, 1995

    [21] 常兴华, 王年华, 马戎等. 并行重叠/变形混合网格生成技术及其应用. 气体物理, 2019, 4(6): 12-21 (Chang Xinghua, Wang Nianhua, Ma Rong, et al. Dynamic hybrid mesh generator coupled with overset and deformation in parallel environment. Physics of Gases, 2019, 4(6): 12-21 (in Chinese)

    Chang Xinghua, Wang Nianhua, Ma Rong, et al. Dynamic hybrid mesh generator coupled with overset and deformation in parallel environment. Physics of Gases, 2019, 4(6): 12-21 (in Chinese)

    [22]

    Ueno Y, Ochi A. Airframe noise prediction using Navier-Stokes code with Cartesian and boundary-fitted layer meshes//25th AIAA/CEAS Aeroacoustics Conference, Delft, 2019

    [23]

    Chesshire G, Henshaw WD. Composite overlapping meshes for the solution of partial differential equations. Journal of Computational Physics, 1990, 90: 1-64 doi: 10.1016/0021-9991(90)90196-8

    [24]

    Sitaraman J, Floros M, Wissink A, et al. Parallel domain connectivity algorithm for unsteady flow computations using overlapping and adaptive grids. Journal of Computational Physics, 2010, 229(12): 4703-4723 doi: 10.1016/j.jcp.2010.03.008

    [25] 韩少强, 宋文萍, 韩忠华等. 高速共轴刚性旋翼非定常流动高精度数值模拟. 航空学报, 2024, 45(9): 177-196 (Han Shaoqiang, Song Wenping, Han Zhonghua, et al. High-accuracy numerical-simulation of unsteady flow over high-speed coaxial rigid rotors. Acta Aeronautica et Astronautica Sinica, 2024, 45(9): 177-196 (in Chinese)

    Han Shaoqiang, Song Wenping, Han Zhonghua, et al. High-accuracy numerical-simulation of unsteady flow over high-speed coaxial rigid rotors. Acta Aeronautica et Astronautica Sinica, 2024, 45(9): 177-196 (in Chinese)

    [26] 刘周, 周伟江. 适于黏性计算的自适应笛卡儿网格生成及其应用. 航空学报, 2009, 30(12): 2280-2287 (Liu Zhou, Zhou Weijiang. Adaptive viscous Cartesian grid generation and application. Acta Aeronautica et Astronautica Sinica, 2009, 30(12): 2280-2287 (in Chinese) doi: 10.3321/j.issn:1000-6893.2009.12.007

    Liu Zhou, Zhou Weijiang. Adaptive viscous Cartesian grid generation and application. Acta Aeronautica et Astronautica Sinica, 2009, 30(12): 2280-2287 (in Chinese) doi: 10.3321/j.issn:1000-6893.2009.12.007

    [27]

    Lee JD. Development of an efficient viscous approach in a Cartesian grid framework and application to rotor-fuselage interaction. [PhD Thesis]. Atlanta: Georgia Institute of Technology, 2006

    [28] 沈志伟, 赵宁, 胡偶. 可压缩黏性流动笛卡尔网格虚拟单元方法研究. 空气动力学学报, 2014, 32(6): 748-754 (Shen Zhiwei, Zhao Ning, Hu Ou. Numerical research of Cartesian based ghost cell method for compressible viscous flows. Acta Aerodynamica Sinica, 2014, 32(6): 748-754 (in Chinese)

    Shen Zhiwei, Zhao Ning, Hu Ou. Numerical research of Cartesian based ghost cell method for compressible viscous flows. Acta Aerodynamica Sinica, 2014, 32(6): 748-754 (in Chinese)

    [29]

    Liu X, Yang B, Ji C, et al. Research on the turbine blade vibration based on the immersed boundary method. Journal of Fluids Engineering, 2018, 140(6): 061402 177

    [30]

    Yang B, Song M, Zhu G. Research on the ghost cell immersed boundary method for compressible flow. Processes, 2024, 12: 1182 doi: 10.3390/pr12061182

    [31] 罗灿炎, 毕林, 徐晶磊等. 笛卡尔网格下不同湍流模型的壁面函数方法研究. 工程力学, 2024, 41(8): 11-22 (Luo Canyan, Bi Lin, Xu Jinglei, et al. Study on wall function method of different turbulence models based on cartesian grid. Engineering Mechanics, 2024, 41(8): 11-22 (in Chinese)

    Luo Canyan, Bi Lin, Xu Jinglei, et al. Study on wall function method of different turbulence models based on cartesian grid. Engineering Mechanics, 2024, 41(8): 11-22 (in Chinese)

    [32]

    Ruffin SM, Sekhar S. A normal ray refinement technique for Cartesian-grid based Navier–Stokes solvers. International Journal of Computational Fluid Dynamics, 2012, 26(4): 231-246 doi: 10.1080/10618562.2012.691970

    [33]

    Zaki M, Ruffin SM. Conservation and grid adaptation enhancements to a normal ray refinement technique for cartesian-grid based navier-stokes solvers//50th Aiaa Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Nashville, 2012

    [34]

    Arslanbekov R, Kolobov V, Ruffin SM, et al. Implementation and evaluation of normal ray refinement technique in adaptive cartesian framework//42nd Aiaa Fluid Dynamics Conference and Exhibit, New Orleans, 2013

    [35]

    Bopp MS, Dement DC, Ruffin SM, et al. An improved object-oriented cartesian grid framework implementing three dimensional normal ray refinement//32nd AIAA Applied Aerodynamics Conference, Atlanta, 2014

    [36]

    Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. Journal of Computational Physics, 1997, 135: 250-258 doi: 10.1006/jcph.1997.5705

    [37]

    Leer BV. Towards the ultimate conservation difference scheme V. A second-order sequel to Godunov’s Method. Journal of Computational Physics, 1979, 32: 101-136

    [38]

    Tritton DJ. Experiments on the flow past a circular cylinder at low reynolds numbers. Journal of Fluid Mechanics, 1959, 6(4): 547-567 doi: 10.1017/S0022112059000829

    [39]

    Zhong M, Zou S, Pan D et al. A simplified discrete unified gas kinetic scheme for incompressible flow. Physics of Fluids, 2020, 32: 093601 doi: 10.1063/5.0021332

    [40]

    Jawahar P, Kamath H. A high-resolution procedure for euler and Navier-Stokes computations on unstructured grids. Journal of Computational Physics, 2000, 164(1): 165-203 doi: 10.1006/jcph.2000.6596

    [41]

    Dervieux A, Rizzi A, Van Leer B, et al. Numerical simulation of compressible Euler flows. Biotechnology Progress, 1989, 30(3): 523–534

图(21)  /  表(3)
计量
  • 文章访问数:  0
  • HTML全文浏览量:  0
  • PDF下载量:  0
  • 被引次数: 0
出版历程
  • 网络出版日期:  2025-04-16
  • 刊出日期:  2025-06-17

目录

/

返回文章
返回