“数值传热学”中贴体网格与非结构化网格上的有限容积法对比教学

宇 波， 禹国军， 韩东旭， 李敬法， 王 艺， 孙东亮

（1.北京石油化工学院机械工程学院，北京 102617；2.上海海事大学商船学院，上海 201306；3.中国石油大学（北京）机械与储运工程学院，北京 102249）

0 引言

流动与传热数值计算中多涉及不规则的计算区域，因此不规则区域的流动与传热数值计算理论是流动传热仿真软件开发的重要基础理论，其中贴体网格和非结构化网格有限容积法是核心内容。然而，教学经验表明，初学者不易准确区分贴体网格和非结构化网格，认为不规则的网格都是非结构化网格，对两者的有限容积法计算理论（如控制方程的离散与求解等）更是困惑。此外，由于两种网格的定义、内涵、有限容积理论等知识点在教材及教学过程中又较为分散，学生难以对非规则区域有限容积法体系进行整体把握，因此在分别完成非结构化有限容积法各教学内容之后，进行一次完整理论体系涉及的各知识要点的综合辨析教学将会是一个有效的手段。综合辨析教学可以让学生更容易厘清两者的异同、掌握两者对应的控制方程离散方法、理解两者对物理问题的适应性，实现对该理论体系理解和应用能力的升华，从而为将来进行流动传热数值仿真软件的自主研发奠定基础。

1 贴体网格和非结构化网格的差异

推荐采用图1 所示的两套网格系统来区分贴体网格和非结构化网格。为了讲清楚贴体网格的定义“贴体网格是一种通过建立物理域边界点与计算域边界点的映射关系，从而得到物理域与计算域内部节点间的对应关系，将不规则的物理域转化成规则的计算域，然后由计算域上规则的结构化网格根据对应关系转换到原物理域上的网格”［1］。以图1（a）来进行解释：为了生成图1（a1）所示的不规则六边形区域的网格，将该不规则区域的边界映射到图1（a2）所示的正方形边界上，然后生成计算域上的规则的结构化网格，最后根据边界的映射关系将计算域上的网格转换到原物理域上，得到图1（a1）所示的贴体网格。通过该例，可以让学生更容易理解为什么贴体网格本质上是结构化网格：贴体网格是结构化网格的一个数学变换，其具有统一的拓扑结构。对于非结构化网格，建议从网格生成的过程来理解其定义：即非结构化网格是一种直接生成适应物理域几何形状的不规则的网格，网格点生成的顺序不存在几何上的连续性，即相邻的网格点的生成顺序不一定是连续的。以图1（b）所示的非结构化三角形网格为例进行说明：对于图1（b）所示的某个三角形单元，无法直接给出其邻点的关系，因为网格点的生成顺序是无规律的，因此邻点的关系须在计算时进行判断。通过这样的讲述，学生就很容易理解非结构化网格的本质特征“网格不具有统一的拓扑结构”。通过以上两套网格的对比给出结论：贴体网格和非结构化网格的本质差异就在于结构化网格节点具有统一的拓扑结构，而非结构化网格节点拓扑结构杂乱。

图1 贴体网格和非结构化网格

2 贴体网格和非结构化网格的有限容积法实施过程差异

课堂教学中可采用表1 来对比介绍基于贴体网格和非结构化网格的有限容积法实施过程的差异。通过该表，让学生从宏观上理解两种方法的差异及造成差异的原因，即：由于网格坐标系的差异，两种方法所采用的控制方程的形式、控制方程及边界条件的离散均存在较大的差异。基于贴体网格的有限容积法在与物理域不同的计算域上进行控制方程的离散求解，因此控制方程及边界条件需要转化到计算域对应的坐标系，而非结构化网格则直接在物理域上进行离散求解，无须进行坐标变换，但由于非结构化网格不具有固定的邻点关系，不易直接计算偏导数，因而采用矢量形式的控制方程［1-3］。

表1 基于两种网格系统的有限容积法整体差异对比

3 贴体网格和非结构化网格有限容积法控制方程的异同及变换要点

以下式所示的二维对流扩散方程为例，对比介绍基于这两种网格的有限容积法所采用的控制方程的差异及方程转换的要点。

式中：ρ为流体密度；x、y为空间坐标分量；t为时间坐标分量；u、v分别为物理平面上x和y方向的速度分量；φ为通用变量；Γφ为扩散系数；Sφ为源项。

3.1 贴体网格有限容积法控制方程及转换要点

（1）控制方程变换。贴体网格对应的计算域所基于的坐标系不同于直角坐标系，因此须对式（1）所示的适用于直角坐标系的控制方程进行坐标变换后方能使用。初学者对于贴体坐标系下的控制方程的形式较为陌生，建议以图2 所示的转换过程进行教学。通过该图，让学生掌握坐标变换的本质，即通过一定的坐标变换法则将物理域上的控制方程转换为计算域上的控制方程，从而在规则的计算域上进行离散求解。而该转换法则［1］就是方程转换的关键。接着举例介绍转换法则的使用方法：图2 中ψ为待求导的物理量或物理量的组合，例如式（1）中的ρuφ和Γφ（∂φ/∂x）等均可作为ψ代入图中所示的相应坐标变换法则中进行转换。二阶偏导数的转换等价于作2 次一阶偏导数的转换，例如将Γφ（∂φ/∂x）作为ψ值进行第1 次坐标变换后得到了下式所示的表达式，其中含有φx（即φ对x的一阶偏导），再将其进行一次坐标变换即可［4-6］。此外，通过此图也让学生了解为何只有跟空间坐标有关的项需要进行坐标变换，非稳态项无须进行变换，源项中只需变换跟坐标有关的部分。

图2 贴体网格控制方程的坐标变换方法

式中：J为雅克比系数；（）x，（）ξ，（）η分别表示括号里的变量（或变量组合）对坐标x、ξ和η的偏导数。

（2）边界条件的变换。通过边界条件的统一形式aφ+b（∂φ/∂n）=c来引出边界条件也应进行坐标变换，原因是其中的一阶导数项跟坐标有关。然后以图3 所示的左边界为例介绍一阶导数的变换方法。

图3 边界条件的处理

由图3 可知，对于物理域左边界上任意一点P，其切线方向即为η坐标在当地的方向，按照基矢量变换原理［1］可得：gη=xηex+yηey。向学生解释清楚这一点后，按图4 所示的流程介绍如何由该切向基矢量推导得到左边界的外法向导数。接着让学生举一反三，推导得到右边界、上边界和下边界的外法向导数的表达式。由此，学生就比较容易掌握不同边界的边界条件在贴体坐标系下的表达式的来龙去脉。

图4 左边界外法向导数的推导思路

3.2 非结构化网格有限容积法采用的控制方程及边界条件

与贴体网格不同，非结构化网格基于直角坐标系，因此，进行有限容积法离散时控制方程无须进行坐标变换。但是，应向学生讲清楚非结构化网格有限容积法采用矢量形式的控制方程，其原因是：非结构化网格不存在统一的网格节点拓扑结构，难以直接给出统一的偏导数的离散格式，因此采用如下式所示的矢量形式的控制方程［1-3］，即

式中，u为速度矢量。

接着介绍边界条件处理的特点：边界条件通用表达式aφ+b（∂φ/∂n）=c中的一阶偏导数通过与边界相邻的单元p0的梯度计算得到，即（∂φ/∂n）=∇φp0·n（其中n为由p0指向相邻节点的方向向量）。

4 两种网格控制方程离散过程的差异

通过图5 介绍两种方法的计算网格系统及网格参数，让学生了解贴体网格系统跟直角坐标系结构化网格系统无本质区别，而非结构化网格系统与前面所学的结构化网格系统有本质的区别。以此引出贴体网格离散方法与结构化网格完全一致，而非结构化网格的离散是一种全新的方法，本次课将重点学习。

图5 贴体网格和非结构化网格以及计算节点布置

通过下式所示的对流项的ξ方向分量在贴体网格的离散让学生理解贴体网格有限容积法方程离散方法跟结构化网格有限容积法完全一样，即控制容积上的分部积分；不同之处在于贴体网格有限容积法的方程离散是在与物理域不同的计算域上进行［4-6］。其公式为

通过以下两个步骤介绍非结构化网格有限容积法控制方程的离散方法：①结合高斯散度定理［7-8］，将式（3）在控制容积V上进行积分可得式（5）所示的积分形式的控制方程；②以扩散项为例，按式（6）进行离散。同理可写出其他项的离散形式［9-10］。

式中：A和A分别为控制容积界面的面积和面积矢量；j为界面的编号。

让学生观察式（6）所示的离散格式以及前面所讲的边界条件的表达式，让他们发现梯度计算在非结构化网格有限容积法离散中是核心。接着介绍求解梯度常采用的Green-Gauss 梯度法［11-13］及最小二乘法［14-15］，具体方法用图6 进行对照式讲解。

图6 两种梯度求解法的思路对比

5 两种网格有限容积法的适应性

采用表2 对比介绍两种网格有限容积法的差异，让同学们从整体上把握二者的特征及适应性。

表2 贴体网格与非结构化网格有限容积法综合对比

结合表2 的分析，可以给出以下建议：对一些几何形状相对规则、流动结构较为复杂又需要高精度的问题，建议用贴体网格法，而对于不规则性较大的一般性流动与传热问题，建议使用非结构化网格。

6 结语

不规则区域问题是流动与传热数值计算普遍遇到的问题，基于贴体网格和非结构化网格有限容积法是解决此类问题普遍采用的方法。作为商业软件的使用者，应当了解该数值计算理论，作为肩负历史重任的未来潜在软件开发者，更是要踏实掌握该根基性的理论。虽然该内容一直是数值传热学教学的难点，但只要从本文所述的5 个方面进行对比式教学，可以让学生在对比学习中深入了解两种方法的本质，准确掌握两种方法的技术细节，为将来成为流动传热仿真软件开发人才、实现计算模拟软件国产化、突破技术封锁奠定基础。

·名人名言·

大学的荣誉，不在它的校舍与人数；而在于它一代一代人的质量。

——柯南特