初学者在进行相关流体仿真时，小球入水往往是一个非常经典的案例。重要的是，许多研究者曾对小球入水进行了大量的实验研究，保留了丰富的实验成果。这使得人们在进行入水领域的仿真工作时，往往利用小球入水来进行数值算法的验证。本文将利用ANSYS Fluent中的动网格局部重构法（Local Remeshing），动态层铺法（Dynamic Layering）以及嵌套网格（Overset mesh）三种手段对小球入水进行简单的仿真模拟。
几何模型
真实模型是三维的，这里为了方便起见，采用二维模型考虑。模型尺寸如图所示。单位为mm。重心坐标(0,0,0)，惯性矩Ixx = 20431.71mm^4，Iyy= 20431.71mm4，Izz=40863.42mm4，面积A=506.71mm^2。
对于平面模型计算，fluent默认拉伸厚度为1m，假设小球为钢质小球，密度7800kg/m^3，因此可以计算小球质量mass=7800*506.71e-6=3.952kg。因此UDF采用DEFINE_SDOF_PROPERTIES宏，给质量及三方向惯性矩赋值即可。将以下代码保存为sphere.cpp文件。
1 |
|
局部重构法（Local Remeshing）
在使用非结构网格的区域上一般采用弹簧光顺模型进行动网格划分，但是如果运动边界的位移远远大于网格尺寸，则采用弹簧光顺模型可能导致网格质量下降，甚至出现体积为负值的网格，或因网格畸变过大导致计算不收敛。为了解决这个问题，Fluent 在计算过程中将畸变率过大，或尺寸变化过于剧烈的网格集中在一起进行局部网格的重新划分，如果重新划分后的网格可以满足畸变率要求和尺寸要求，则用新的网格代替原来的网格，如果新的网格仍然无法满足要求，则放弃重新划分的结果。
网格划分
在ICEM CFD中划分网格，根据几何模型创建点线，如图所示。
根据四条矩形边线创建面，然后以小球边线为划分线划分整个面，之后删除小球内部的面，如图所示。
点击创建体，根据矩形对角的两个点进行创建体BODY。接下来创建各个part，分别将上边界（OUTLET），左右边界（SIDEWALL），下边界（DOWNWALL），小球边界（SPHEREWALL），流场面（FLUID）创建为part，如图所示。
设置全局网格参数如下，生成全三角形计算网格。
各个part的网格划分参数如下。
点击compute mesh创建面网格，得到粗略的网格划分如图，可见，小球外层应有明显的边界层网格。
查看网格质量如图。
点击smooth elements globally，进行网格光顺，参数如下。
光顺后的网格质量如下，网格质量在0.5以上。由于这个例子只做简单的教程，因而不进一步完善网格。
仿真计算
点击保存,注意路径不能有中文。点击OUTPUT，选择求解器Fluent，点击保存，选择刚才的prj，弹出选项框，可修改命名，保存为msh。点击保存并输出msh文件。Fluent读入msh文件，scale网格设定合适的单位。将网格尺寸设为mm，显示尺寸设为m。如图所示。
设置VOF多相流模型，（models-Multiphase:多相设置，流体中有两相或以上的介质），设置body force formulation为Implicit Body Force。
Viscous表示粘性方程设置，（可选择无粘流动、层流流动、S-A方程湍流模型、k-e两方程湍流模型、k-Ω两方程湍流模型、雷诺应力模型、离散涡模型、大涡模型等）。选择湍流模型为标准k-epsilon模型，采用增强壁面函数Near-Wall。材料数据库中添加water-liquid（Materials-Fluid-Air-Fluent Database-Copy）。设置水为主相，空气为第二相(Multiphase-phases-phase1-water-phase2-air)。
本例中需要设置的边界条件为出口边界。该边界设置为压力出口，静压为0。设置第二相空气相的体积分数为1。压力出口边界点击operating conditions，将参考压力边界设置为y=0.1m(即上边界所处的y轴坐标)，如图所示。
点击initialization对流场进行初始化设置，在弹出的框中下方点击initialize。之后右键cell registers→new→field variable（19前的版本选择iso value）
在弹出的窗口中，选择field value of 分别为cell info 和 cell element type，type选择为cells in range，iso-min和iso-max均设置为3，点击save/display。这一步的主要目的是选择出小球的边界层网格，如图所示。
点击separate→cells，进行网格剖分，如图所示。
现在点击cell zone conditions，可以发现，多了一个fluid:008。之后在菜单栏中的surface点击create→zone，在左侧选择fluid:008，点击create。
点击菜单栏中的user-defined→functions→complied，编译并load此UDF。
激活动网格。在Dynamic mesh面板中，激活smoothing与Remeshing方法，同时勾选Six DOF选项。点击mesh methods 下的settings，选择spring layer。然后点击remeshing，点击default，然后点击ok。
关闭mesh method settings面板，进入six dof setting面板。勾选Write Motion History，这样运动数据即可写入文本文件。设置完毕后关闭6DOF设置面板。Create创建运动区域，如图所示，指定运动区域为spherewall，勾选six dof options中的on选项。点击create。
之后指定运动区域为fluid:008，勾选six dof options中的on选项和passive选项。点击create。
之后右键cell registers→new→region，创建参数如下。
设置水的静压梯度。点击菜单栏中的user-defined→field functions→custom，将静压表达式P=rou×g×h写入，由于在这个例子中水面所处坐标为y=-0.24m，因此表达式写为998.2×9.8×(-0.24-y)，点击define即可定义该场函数，如图所示。
设置完毕后，保存case与data文件。点击initialization，之后点击patch，设置参数如图。点击patch。
设置水的静压梯度。设置参数如下，勾选field function，选择之前定义的场函数，对划分好的region进行patch(也就是对水域进行静压梯度设置)。
点击results→graphics→contours，新建一个云图，选择各相的体积分数云图。
如果显示如图，则patch无误。
另外，再新建一个压力云图，如果显示如图，则水的静压梯度设置无误。
点击calculation activities→autosave，设置每100步自动保存。
点击calculation activities→solution animations，设置每100步自动保存动画文件。
设置时间步长time size为0.0001s，时间步数time step为4000，即计算0.4s情况。设置内迭代数量为20，然后进行计算。
仿真结果
可见入水空泡壁面非常粗糙，但是大致给出了入水空泡的扩张和颈缩状态，这主要是因为网格本身密度不够以及网格变形的问题。如果要得到更加精细的结果，可以对水面处的网格和小球运动区域的网格进行加密处理，另外，调整网格变形的参数也能得到很好的效果。
动态层铺法（Dynamic Layering）
对于棱柱型网格区域（六面体和或者三棱柱），可以应用动态层铺法。动态层铺法的中心思想是根据紧邻运动边界网格层高度的变化，添加或者减少动态层，即在边界发生运动时，如果紧邻边界的网格层高度增大到一定程度，就将其划分为两个网格层；如果网格层高度降低到一定程度，就将紧邻边界的两个网格层合并为一个层。
几何模型
几何模型示意图如图：
网格划分
在ICEM CFD中划分网格，需要进行两套网格划分然后进行网格拼接，这里均采用结构化网格划分方式。最需要注意的一点是要对各个边界进行清晰命名，这样在设置interface时就不会弄混。
侧边网格的最终划分及part设置如图，并save mesh as保存为sidemesh.uns。
变形网格的划分需要多点步骤，其part设置为：
点击blocking→create block，选择Type为2D planar，将part命名为slipmesh，点击apply。
点击blocking→split block，创建（±50，±50）四个几何点，并采用点划分的方式对block进行划分。
点击blocking→split block，点击ogrid block，点击select blcoks，选择小球周围的block，并中键确定，点击apply，完成ogrid block创建，如图所示。
点击blocking→associate，完成edge与curve之间的关联，在对小球边界进行关联时，需要勾选project vertices，然后点击apply即可，如图所示。
点击blocking→Delete block，点击小球内部的block并删除。
点击blocking→Pre-Mesh Params，点击Edge Params，外围的边线节点分布数量均与其长度相等，切均匀分布，并勾选copy parameters。内部四条斜边的节点分布数量为58，且由内向外由密集到稀疏，参数设置如下：
点击Pre-Mesh，进行网格划分如图，并save mesh as保存为slipmesh.uns。
打开slipmesh.uns，点击file→mesh→open mesh，导入sidemesh.uns，点击merge，点击covert to unstruct mesh，随即导出网格layermesh.msh。
仿真计算
Fluent读入msh文件，scale网格设定合适的单位。将网格尺寸设为mm，显示尺寸设为m。选择基于压力的瞬态求解器，开启重力。点击菜单栏中的user-defined→functions→complied，编译并load之前保存好的UDF。
设置VOF多相流模型，设置body force formulation为Implicit Body Force。
选择湍流模型为标准k-epsilon模型，采用增强壁面函数。材料数据库中添加water-liquid。设置水为主相，空气为第二相。
在boundary conditions将四个交界面设置为interface，顶端边界均设置为压力出口，静压为0。设置第二相空气相的体积分数为1。压力出口边界点击operating conditions，将参考压力边界设置为y=0.1m(即上边界所处的y轴坐标)，如图所示。
在mesh interface中用Manual create对交界面进行配对，如图所示。
激活动网格。在Dynamic mesh面板中，激活smoothing与Layering方法，同时勾选Six DOF选项。点击mesh methods 下的settings，选择spring layer。然后点击Layering，采用默认设置，然后点击ok。进入six dof setting面板，勾选Write Motion History，这样运动数据即可写入文本文件。设置完毕后关闭6DOF设置面板。
创建运动与变形区域，如图所示，指定运动区域为spherewall，勾选six dof options中的on选项。在mesh option 中设置层高为0.5mm，点击create。之后指定运动区域为slipmesh，勾选six dof options中的on选项和passive选项。点击create。接下来设置变形区域的顶端边界和底端边界为stationary类型，在mesh option 中设置层高为0.5mm，点击create。如图所示。
之后右键cell registers→new→region，创建参数如下。
设置水的静压梯度。点击菜单栏中的user-defined→field functions→custom，将静压表达式P=rou×g×h写入，由于在这个例子中水面所处坐标为y=-0.24m，因此表达式写为998.2×9.8×(-0.24-y)，点击define即可定义该场函数，如图所示。
设置完毕后，保存case与data文件。点击initialization，之后点击patch，设置参数如图所示，从而选出水域所在的网格。
设置水的静压梯度。设置参数如下，勾选field function，选择之前定义的场函数，对划分好的region进行patch(也就是对水域进行静压梯度设置)。
点击results→graphics→contours，新建一个云图，选择各相的体积分数云图。如果显示如图，则水域patch无误。(由于网格划分的问题，初始气液界面处存在一些起伏，这主要是因为O-BLOCK划分造成的)
另外，再新建一个压力云图，如果显示如图，则水的静压梯度设置无误。
点击calculation activities→autosave，设置每100步自动保存。点击calculation activities→solution animations，设置每100步自动保存动画文件。设置时间步长time size为0.0001s，时间步数time step为4000，即计算0.4s情况。设置内迭代数量为20，然后进行计算。
仿真结果
下面分别给出了0.26s，0.32s，0.34s时刻下的体积分数云图，可见小球入水后形成了一个较长的入水空泡，并在液面处形成了明显的飞溅，空泡在不断拉长的过程中先是扩张然后颈缩，并在入水深度约1/2处完成了空泡深闭合。0.26s时，已经能观察到空泡趋于表面闭合，但是很不明显。同时，空泡的最大直径仅略大于小球直径，可见空泡的扩张也非常弱。在0.34s时甚至观察到了空泡多次闭合的现象。可见相对于网格重构法的入水空泡壁面，动态层铺方法计算的壁面更加光滑，对于自由边界的行为也有着较好的捕捉。
嵌套网格（Overset mesh）
嵌套网格（overset）又称重叠网格，该网格由两部分组成：包含主要流体区域的背景网格（background），围绕部件的构建的单独的组件网格（component）。此时，部件网格的外边界指定为overset（特有的边界类型），其他的边界和通常的网格则没有区别。相比于传统的MDM动网格方式，嵌套网格有以下优势：
更简单的网格划分方式（边界更为简单）
当算例需要更改网格布局的时候，仅需要添加/删除部件网格即可
当涉及到边界运动的时候，通常仅需要部件网格独立进行运动即可
不会出现负体积
网格质量不随边界移动而变化
不需要思考复杂的边界运动逻辑
总的来讲，嵌套网格的最大特点，就是让流体仿真变得更为简单和快捷，提高流体工程师的工作效率。
几何模型
背景网格（background）如图所示：
组件网格（component）如图所示：
网格划分
在ICEM CFD中划分网格，这里对背景网格和组件网格均采用结构化网格划分方式。由于背景网格划分十分简单，因此这里只对前景网格的O-BLOCK做一点记录。
在ICEM中建立好几何模型后（记得删除小球内部的面）进行part创建，将外围4条边命名为OVERSETWALL，将小球边界命名为SPHEREWALL，如图所示。
点击blocking→create block，选择Type为2D planar，将part命名为frontmesh，点击apply。
点击blocking→split block，点击ogrid block，点击select blcoks，选择屏幕中的block，并中键确定，如图所示。
点击apply，完成ogrid block创建，如图所示。
点击blocking→associate，完成edge与curve之间的关联，在对小球边界进行关联时，需要勾选project vertices，然后点击apply即可，如图所示。
点击blocking→Delete block，点击小球内部的block并删除，如图所示。
点击blocking→Pre-Mesh Params，点击Edge Params，外围4条边节点分布数量均为100，切均匀分布，并勾选copy parameters。内部四条斜边的节点分布数量为58，且由内向外由密集到稀疏，参数设置如下：
点击Pre-Mesh，进行网格划分如图。
右键Pre-Mesh，点击covert to unstruct mesh，随即导出前景网格frontmesh.msh。此外，对背景网格进行划分，并保存为background.msh。
仿真计算
Fluent读入background.msh文件，点击Append→append case file，读入frontmesh.msh文件，scale网格设定合适的单位。将网格尺寸设为mm，显示尺寸设为m。选择基于压力的瞬态求解器，开启重力。点击菜单栏中的user-defined→functions→complied，编译并load之前保存好的UDF。
设置VOF多相流模型，设置body force formulation为Implicit Body Force。
选择湍流模型为标准k-epsilon模型，采用增强壁面函数。材料数据库中添加water-liquid。设置水为主相，空气为第二相。
点击cell zone conditions→frontmesh，在弹出的窗口中勾选Mesh Motion，如图所示。
本例中需要设置的边界条件为出口边界。该边界设置为压力出口，静压为0。设置第二相空气相的体积分数为1。
将oversetwall边界条件type设为overset，其余边界条件均为壁面。
双击overset interfaces，弹出的窗口中命名为oversetinterface，选择background和frontmesh，点击create即可。
激活动网格。在Dynamic mesh面板中，激活smoothing方法，同时勾选Six DOF选项。进入six dof setting面板，勾选Write Motion History，这样运动数据即可写入文本文件。设置完毕后关闭6DOF设置面板。创建运动区域，如图所示，指定运动区域为spherewall，勾选six dof options中的on选项。点击create。
之后指定运动区域为frontmesh，勾选six dof options中的on选项和passive选项。点击create。
之后右键cell registers→new→region，创建参数如下。
设置水的静压梯度。点击菜单栏中的user-defined→field functions→custom，将静压表达式P=rou×g×h写入，由于在这个例子中水面所处坐标为y=-0.24m，因此表达式写为998.2×9.8×(-0.24-y)，点击define即可定义该场函数，如图所示。
设置完毕后，保存case与data文件。点击initialization，之后点击patch，设置参数如动态层铺一节所示，从而选出水域所在的网格并patch水域的静压梯度。
点击results→graphics→contours，新建一个云图，选择各相的体积分数云图。如果显示如图，则patch无误。
另外，再新建一个压力云图，如果显示如图，则水的静压梯度设置无误。
点击calculation activities→autosave，设置每100步自动保存。点击calculation activities→solution animations，设置每100步自动保存动画文件。设置时间步长time size为0.0001s，时间步数time step为4000，即计算0.4s情况。设置内迭代数量为20，然后进行计算。
仿真结果
下面分别给出了0.26s，0.32s，0.34s时刻下的体积分数云图，可见小球入水后形成了一个较长的入水空泡，并在液面处形成了明显的飞溅，空泡在不断拉长的过程中先是扩张然后颈缩，并在入水深度约1/2处完成了空泡深闭合。在这个例子中，空泡的最大直径接近于两倍的小球直径，可见空泡的扩张速度较高，但是空泡表面闭合现象没有观察到，0.34s时可以观察到初生的沃辛顿射流。可见相对于网格重构法的入水空泡壁面，重叠网格方法计算的壁面更加光滑，对于自由边界的行为也有着较好的捕捉。
