# 使用Gmsh为OpenFoam轴对称模拟准备结构化网格

OpenFoam在底层只接受三维网格。虽然也具有进行二维仿真、轴对称仿真的能力，但实际上是通过高层处理在三维网格上实现的。这样一来，在Gmsh画网格过程中就需要一些技巧。

其中，实现二维模拟相对简单，只需要画出来一个二维网格，然后通过`Extrude`拉伸一层，形成一块由三棱柱（来自三角网格）或四棱柱（来自四边形网格）拼接而成的薄板即可。仿真时可以将薄板的两个大面的边界条件类型设置为`empty`。

{{ :science:openfoam:2d_unstructure_for_openfoam.png?400 | 用于OpenFoam二维仿真的非结构网格}}

[[https://www.youtube.com/watch?v=C_sCHpLaz3Y|这是一个关于制作二维OpenFoam网格的视频教程]]

轴对称网格的要求则要复杂一些，包括以下：

+ 先画出二维截面，而后旋转扫略为约$5°$的楔形体
+ 楔形体关于$z=0$平面对称
+ 楔形体两边网格点对称
+ 计算时，把楔形体两边两个对称面的边界条件类型设置为`wedge`

这些要求也不算十分复杂。[[https://www.youtube.com/watch?v=MUGnlYncqck|这里]]就有一个利用轴对称做自由射流中的轮胎面的流场的例子。

但是，简单是对于非结构网格而言的。如果想要结构化网格，Gmsh的一个bug或者feature会带来一些麻烦：**Gmsh中的`Rotate`命令应用于`Surface`对象时，会改变组成`Surface`的那些`Line`的编号。**

Gmsh中创建结构网格的方法是，利用`Transfinite`特征，规定某一条线上有多少网格点，继而指定每一块结构化网格的四个角点——角点之间可以有一条或多条线，但所有线被细分的**段数**应该与网格对面的两个角点间的段数相同。

不管怎么说，`Transfinite`信息必须被赋予到线上，而且这一步必须在`Extrude`之前完成。`Extrude`之后，几何信息就成了楔形体，没法再区分出六面体来了。

同时，为了保证楔形体关于$z=0$平面对称，必须首先在$z=0$平面上画图，然后把平面向一个方向旋转$2.5°$，然后回马枪扫略$5°$，[[https://www.youtube.com/watch?v=MUGnlYncqck|上面提到的视频教程]]也是这么干的。不幸的是，如果你对做好了`Transfinite`信息的`Surface`这么干，最后出来的都是非结构网格。

分析原因，大概是因为，`Transfinite`信息只是一个“记录信息”的动作，命令里面写`Transfinite`的时候，不是说结构化网格马上就画好了，而是把这些信息记录到线上去，等待后面真的画网格的时候用。但是，`Rotate`命令应用于`Surface`对象时，好像是把原来的线全删掉了，然后在新的位置重新建立了一些线（观察到ID都不一样），而那些与线的ID绑定的`Transfinite`信息就这样丢失了。

Gmsh甚至专门加了一个`Geometry.CopyMeshingMethod`变量，在`Rotate`命令之前将其置为`1`可以复制`Transfinite`信息到旋转后的线，仅对`Builtin`内核有效。但实际用的时候发现，加了这个会报错说“仅能复制自动生成的`Transfinite`信息”，画网格过程倒是也能继续，但画的过程也会继续零星报错，最后画出来东西有没有错不得而知，我没花时间尝试。

兜兜转转一圈下来，貌似最好就是不要用`Rotate`命令。我尝试了一下不旋转直接扫略，确实画出来的楔形体网格很漂亮，但它不符合$z=0$平面对称的条件啊。理论上也可以在画2D截面时就把他画在负向旋转过$2.5°$的平面上，但是这样搞的话算坐标的工作显然累死人。

然后开始想，能不能先画好网格再整体旋转$2.5°$？理论上很容易实现，但没发现相关命令。那能不能先旋转`Surface`再写`Transfinite`信息？理论上可以，但新的线条ID全是动态生成的，写起来也很累。

那可否将旋转操作再向前推动，画完点就把他们旋转过去？这样初次画出的线就是旋转过的。测试一下，居然可行！点的旋转不影响点的ID，并且还挺方便的，可以直接用`:`符号批量旋转。

下面提供一个例子：

<code text structured-mesh.geo>
// 不依赖于Builtin内核
SetFactory("OpenCASCADE");

lc = 0.00025;

// 在x-y平面画点
Point(1) = {-0.01, 0.02, 0, lc};
Point(2) = {0, 0.02, 0, lc};
Point(3) = {0, 0.025, 0, lc};
Point(4) = {0.08, 0.025, 0, lc};
Point(5) = {0.09, 0.035, 0, lc};
Point(6) = {0.1, 0.035, 0, lc};
Point(7) = {0.125, 0.01, 0, lc};
Point(8) = {0.127, 0.01, 0, lc};
Point(9) = {0.129, 0.0105, 0, lc};
Point(10) = {0.129, 0, 0, lc};
Point(11) = {-0.01, 0, 0, lc};
Point(12) = {0, 0, 0, lc};
Point(13) = {0.08, 0, 0, lc};
Point(14) = {0.1, 0, 0, lc};

// 马上旋转以上所有点
Rotate {{1,0,0},{0,0,0},2.5*Pi/180.0}
{
	Point{1:14};
}

// 在旋转过的点上定义线
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 9};
Line(9) = {9, 10};
Line(10) = {10, 14};
Line(11) = {14, 13};
Line(12) = {13, 12};
Line(13) = {12, 11};
Line(14) = {11, 1};

// tune for structured mesh begain
Line(15) = {2, 12};
Line(16) = {4, 13};
Line(17) = {6, 14};
Transfinite Curve {14, 15} = 80 Using Progression 1;
Transfinite Curve {16, 17, 9} = 120 Using Progression 1;
Transfinite Curve {2} = 41 Using Progression 1;
Transfinite Curve {1, 13} = 20 Using Progression 1;
Transfinite Curve {3, 12} = 200 Using Progression 1;
Transfinite Curve {11} = 99 Using Progression 1;
Transfinite Curve {4, 5} = 50 Using Progression 1;
Transfinite Curve {10} = 250 Using Progression 1;
Transfinite Curve {6} = 192 Using Progression 1;
Transfinite Curve {7, 8} = 30 Using Progression 1;
Curve Loop(101) = {14, 1, 15, 13};
Plane Surface(201) = {101};
Curve Loop(102) = {15, -12, -16, -3, -2};
Plane Surface(202) = {102};
Curve Loop(103) = {16, -11, -17, -5, -4};
Plane Surface(203) = {103};
Curve Loop(104) = {17, -10, -9, -8, -7, -6};
Plane Surface(204) = {104};
Transfinite Surface {201} = {11, 1, 2, 12};
Transfinite Surface {202} = {12, 3, 4, 13};
Transfinite Surface {203} = {13, 4, 6, 14};
Transfinite Surface {204} = {14, 6, 9, 10};
Recombine Surface {201, 202, 203, 204};
// tune for structured mesh end

// 旋转扫略
fluid[] = Extrude {{1,0,0},{0,0,0},-5*Pi/180.0}
{
	Surface{201}; Surface{202}; Surface{203}; Surface{204};
	Layers{1};
	Recombine;
};

Physical Surface("wedge0") = {201, 202, 203, 204};
Physical Surface("wedge1") = {208, 212, 216, 221};
Physical Surface("inlet") = {205};
Physical Surface("inlet-wall") = {206, 211};
Physical Surface("propellant") = {210};
Physical Surface("after-burner-wall") = {215, 214, 220};
Physical Surface("nozzle-wall") = {219, 218};
Physical Surface("outlet") = {217};

Physical Volume("fluidVolume") = {1, 2, 3, 4};

</code>

{{ :science:openfoam:2d_structure_for_openfoam.png?400 | 用于OpenFoam轴对称仿真的结构网格}}