CFD网格雕刻从入门到精通(1):HelloWorld

摘要

目前所有的数值模拟研究和应用中都离不开网格(mesh), 而且网格的质量直接影响这模拟结果的质量。从今天开始我将带来基于OpenFOAM（后面简称为OF）的动力学模拟网格剖分教程，OpenFOAM支持各种各样的网格格式及转换。我最开始喜欢用 Gmsh ,但是gmsh还是有一定的局限性，在网格质量控制方面不太好实现，而OF对网格的正交性要求还是挺好的。所以逐渐转向用OF自带的网格剖分工具 snappyHexMesh，其实它更应该理解为网格雕刻工具。详情参见 ESI官方文档 .

snappyHexMesh的基本思想

snappyHexMesh的基本思想就是将一个几何模型(用stl文件表达)放在一个基本的正交网格里面，然后经过一定的算法和参数约束将模型的表面雕刻出来，其内部基本上还是正交网格。 这个背景正交网格一般应该比几何模型的范围大，能将几何模型包住，用blockMesh工具生成即可。所以第一步应该使用blockMesh生成背景网格。snappyHexMesh的网格质量控制参数非常多，总共一百多个，所以其功能还是挺强大的而且使用自由度也是挺大。

雕刻流程

用snappyHexMesh进行网格雕刻的工作流程主要分为以下几步：

1. 生成背景网格（或基准网格）

2. 几何体的定义

3. 生成笛卡尔网格

4. 雕刻网格

5. 在表面附近添加层

6. 检查网格质量

图 21 snappyHexMesh的工作流程

这几个流程的示意图如图 21所示。

生成背景网格

背景网格必须符合三个条件：纯六面体(一般用立方体规则网格)网格；网格比例应该接近于1；必须至少有一个cell与STL表面相交。

几何模型(STL文件)

STL文件可以通过很多建模软件获得，比如blender，我这里介绍一种用gmsh生成stl文件的方法(gmsh在建模方面的方便之处不用多讲)。 Gmsh 可以将几何体导出为多种格式，包括STL和基于CAD的step文件。导出为stl的命令为 gmsh -0 xxx.geo -o xxx.stl -format stl 后面的format参数也可以不用指定，会根据输出文件名的后缀名自动判断是什么格式。 也可以导出为step文件(这是 OpenCASCADE )的一种模型描述文件，然后也可以基于 OpenCASCADE 库写一个小程序(比如 step2stl )将step文件转换为stl文件，或者用其他的CAD软件进行进一步编辑几何体。这一步就这么搞定了！ 注意: STL文件通常都是放在 constant/triSurface 目录下的。

snappyHexMesh

snappyHexMesh工具读取 system 目录下的 snappyHexMeshDict 字典文件，这个文件里面给定了具体的网格质量和边界等控制的参数。 雕刻完成之后的网格会写入 constant/polyMesh 目录形成标准的OF网格文件。

snappyHexMesh的基本流程可以分为七步：

1. 生成背景六面体网格

2. 在特征棱处进行单元分裂

单元分裂的更多细节控制在 snappyHexMeshDict 字典文件中的 castellatedMeshControls 子字典中。而STL几何体的特征棱和特征表面信息可以使用 surfaceFeatureExtract 工具进行提取，提取完成后会在 STL文件同目录下生成对应的 .eMesh 文件。

3. 表面上的单元分裂

棱上的单元分裂和细化之后，将对表面单元进行分裂和细化。相应的参数是在 castellatedMeshControls 子字典中通过 refinementMeshControls 关键词指定。

4. 单元移除

完成特征棱和表面上的单元分裂操作之后，就开始进行单元移除。STL的几何体现在处于背景网格中，那么最终需要那一部分作为计算的网格呢？是STL几何体实体还是背景网格减去STL几何体所剩余的部分（STL几何体表面将作为边界，而内部是空的）？这个判断是通过 locationInMesh 关键词指定的三维坐标进行的，最终保留的网格是此点所在的体积。这个关键词也是在 castellatedMeshControls 子字典中。

5. 特定区域的单元分裂 refinementRegions 在 castellatedMeshControls 子字典中指定细化区域。

6. 贴合表面

移除指定区域的单元以后，表面贴合的点会组成一个贴合表面网格。而关于控制贴合过程的参数在 snapControls 子字典中指定。

7. 边界层网格

贴合操作完成之后的网格基本上可以用于模拟计算了，但是沿着表面依然有一些不规则的单元。可以通过 addLayersControls 子字典在表面附近创建与表面平行的层状网格过渡，既可以很好的描述STL几何体本来的形状，又可以得到较好的网格质量。

snappyHexMesh命令和文件解析

OpenFOAM中的很多工具都是通过相应的字典文件读取参数的，snappyHexMesh也不例外，必须的字典文件有: snappyHexMeshDict,当然了一般也会有blockMesh用于生成背景网格，surfaceFeatureExtractDict用于提取STL表面特征，meshQualityDict指定网格质量等。基本的文档结构列举如下：

.
├── constant
│   └── triSurface
│       └── surfacemesh.stl
├── run.sh
└── system
   ├── blockMeshDict
   ├── controlDict
   ├── fvSchemes
   ├── fvSolution
   ├── meshQualityDict
   ├── snappyHexMeshDict
   └── surfaceFeatureExtractDict

snappyHexMeshDict的顶层参数

• castellatedMesh : 关键词 (必须) [true`|:code:`false]

• snap : 关键词 (必须) [true`|:code:`false]

• addLayers : 关键词 (必须) [true`|:code:`false]

• mergeTolerance : 关键词 (必须) [1e-6 (数字)]

• debug : 关键词 (可选) [3 (数字：0-4)]

• geometry : 子字典 (必须)

• castellatedMeshControls : 子字典 (必须)

• snapControls : 子字典 (必须)

• addLayersControls : 子字典 (可选)

• meshQualityControls : 子字典 (可选)

子字典示例：

1. geometry子字典

geometry
{
   sphere.stl // STL 文件名
   {
      type triSurfaceMesh;
      regions
      {
            secondSolid // STL文件中的几何体区域名称
            {
               name mySecondPatch; // patch名称，可用于指定边界条件
            }
      }
   }

   box1x1x1  // 自定义立方体区域1
   {
      type   searchableBox; // 通过对角线两个点的坐标定义
      min    (1.5 1 -0.5);
      max    (3.5 2 0.5);
   }

   sphere2  // 自定义区域2（球形区域）
   {
      type   searchableSphere; // 通过中心和半径定义
      centre (1.5 1.5 1.5);
      radius 1.03;
   }
};

2. castellatedMeshControls子字典

常用参数有：

• locationInMesh : 三维坐标向量，指定哪个区域被网格化。比如 (0 2 0)

• maxLocalCells : 在细化的时候，每个处理器处理的最大单元个数。比如 1e+6

• maxGlobalCells : 网格细化过程中全局网格单元最大个数。比如 2e+6

• minRefinementCells : 如果大于等于细化网格单元个数，表面细化则停止。比如 0

• maxLoadUnbalance : 在网格细化的时候，最大处理器不平衡指标，当此值为0时达到最理想平衡。比如 0.1

• nCellsBetweenLevels : 不同细化级别之间的单元缓冲层数。比如 1

• resolveFeatureAngle : 将最大细化级别应用于可以看到其角度超过此值的相交单元。 比如 30

• allowFreeStandingZoneFaces : 是否允许生成独立的区域面。 [true`|:code:`false]

• features : 需要细化的特征列表

• refinementSurfaces : 需要细化**表面**的字典

• refinementRegions : 需要细化**区域**的字典

castellatedMeshControls
{
   maxLocalCells 100000;
   maxGlobalCells 2000000;
   minRefinementCells 0;
   maxLoadUnbalance 0.10;
   nCellsBetweenLevels 2;
   features
   (
      {
            //可以用surfaceFeatureExtract工具处理得到
            file "surfacemesh.eMesh";
            level 0;
      }
   );
   refinementSurfaces
   {
      banana_stlSurface
      {
            level (2 2);
      }
   }
   resolveFeatureAngle 60; //为了避免曲面处过度细化
   planarAngle 30;
   refinementRegions
   {
      refinementBox
      {
            mode inside;   // inside;
            levels ((1E15 1));  //1E15 4
      }
   }
   locationInMesh (4.0 0.1 0.1);
   allowFreeStandingZoneFaces true;
}

3. snapControls子字典

常用参数

• nSmoothPatch : 在找到相应表面之前的patch光滑迭代的次数。比如 3

• tolerance : 曲面特征点或边所吸引点的距离与局部最大边长的比率。比如 4.0

• nSolveIter : 网格位移松弛迭代次数。比如 30

• nRelaxIter : 捕捉过程的松弛迭代的最大次数。比如 5

• nFeatureSnapIter : 特征边缘捕捉迭代次数。比如 10

• implicitFeatureSnap : 通过采样表面来检测（仅几何）特征。比如 [false`|:code:`true]

• explicitFeatureSnap : 使用castellatedMeshControls指定的特征 比如 [true`|:code:`false]

• multiRegionFeatureSnap : 使用explicitFeatureSnap时检测多个表面之间的特征。 [false`|:code:`true]

4. addLayersControls子字典

关键参数

• layers : 层的字典

• relativeSizes : 层厚度是相对于外层的未变形单元大小还是绝对值？[true`|:code:`false]

• expansionRatio : 层网格的膨胀系数。比如 1.0

• finalLayerThickness : 离壁面最远的层的厚度，相对或绝对同样是通过relativeSizes关键词指定的。比如 1

• firstLayerThickness : 离壁面最近的的层的厚度，相对或绝对同样是通过relativeSizes关键词指定的。比如 0.3

• thickness : 所有层的总厚度。比如 0.3

• minThickness : 所有层的最小总厚度，低于该厚度的表面将不被挤出。比如 0.1

• nGrow :

• featureAngle : 大于此值，则表面不被基础。 比如 60

• maxFaceThicknessRatio : 表面不被挤压的表面厚度比，适用于弯曲的单元。比如 0.5

• nSmoothSurfaceNormals : 曲面法线的平滑迭代次数。比如 1

• nSmoothThickness : 表面patch上的平滑层厚度。比如 10

• minMedialAxisAngle : 用于拾取中轴点的角度。比如 90

• maxThicknessToMedialRatio : 在厚度与中间距离之比较大的地方减少层的生长。比如 0.3

• nSmoothNormals : 内部网格运动方向的平滑迭代次数。比如 3

• nRelaxIter : 捕捉松弛迭代的最大次数。比如 5

• nBufferCellsNoExtrude : 0

• nLayerIter : 层的附加迭代总次数。比如 50

• nRelaxedIter : meshQuality 子字典中的 relaxed 被使用后的最大迭代次数。比如 20

5. meshQualityControls子字典

• maxNonOrtho : 被允许的最大的非正交度数，180表示禁用。比如 65

• maxBoundarySkewness : 许的最大边界面偏度，小于0表示被禁用。比如 20

• maxInternalSkewness : 许的最大内部面偏度，小于0表示被禁用。比如 4

• maxConcave : 允许的最大凹度，180表示禁用。比如 80

• minFlatness : 最小投影面积与实际面积之比，-1被禁用。比如 0.5

• minVol : 金字塔最低体积，较大的负数(e.g. -1e30)则被禁用。比如 1e-13

• minTetQuality : 由面心和可变基点最小分解三角形以及单元中心形成的四面体的最低质量,较大的负数(e.g. -1e30)则被禁用。比如 1e-13

• minArea : 最小面的面积，小于零表示被禁用。比如 -1

• minTwist : 面的最小扭曲，小于-1表示被禁用。比如 0.05

• minDeterminant : 最小归一化单元的行列式（这是个什么鬼指标）,1表示六面体，小于等于0表示错误单元。比如 0.001

• minFaceWeight : 0->0.05。比如 0.05

• minVolRatio : 0->1.0。比如 0.01

• minTriangleTwist : >0则对Fluent兼容。一般设置为 -1.

• nSmoothScale : 错误分布迭代次数。比如 4

• errorReduction : 缩小误差点处的后退位移量。比如 0.75

• nSmoothScale :

应用效果

OpenFOAM技术内幕:热物理模型中的能量场——thermo.he()   硬石膏沉淀对海底热液系统中金属运移的作用