基于Fluent的激光熔覆熔池流场仿真——VOF耦合融化凝固模型+源项UDF、二维模型、自由界面

激光熔覆熔池流场仿真是笔者本科的毕设课题，趁着这些天有时间整理了一下，因为毕业后就没有接触Fluent这个软件了，文中如有错误还请指正。

一、前言

废话不多说，咱就直接步入正题，首先解释一下几个问题。

为什么仿真是二维模型，而不用三维模型？

不是笔者不想用三维模型，实在是电脑太拉了（i5 10代  16G  1050ti），这样的配置跑一次0.5秒的仿真要计算10个小时。上三维模型的话，计算量指数增长，内存和CPU扛不住。除非有条件能租到服务器，不然还是老老实实二维仿真。

如何学习Fluent？

笔者学习Fluent时，看的是B站这个Up主的教程，讲解得很是详细。

怂管木觉兽的个人空间-怂管木觉兽个人主页-哔哩哔哩视频 (bilibili.com)https://space.bilibili.com/503748872

 本篇博客为以下论文的简化介绍，更多详细的内容可查看论文：
激光内送粉变姿态熔覆非水平熔池流场的数值模拟 - 中国知网 (cnki.net)https://kns.cnki.net/kcms2/article/abstract?v=2F6201taHdf5I3Mhnq_1SFWmmO_lD5xDytAOzzH2EvYu21EJrWWq98ajiw30lx_ZuPyg2d8Rh-220iAne0g0PpurfplSQNaj-9IFje238uXlneH6Lqu4Wc2igpDrPgYlgHU8lJeJUFs=&uniplatform=NZKPT&language=CHS

 二、模型简介

如下图所示，本文选取一个特征面（图中红框）作为研究对象，分析该截面处的熔池流场。

首先基于FLUENT软件的离散相模型(DPM)对喷嘴流场进行了计算，获取金属粉末的质量分布关系，并基于此结合激光热源模型，基于多相流 (VOF)耦合凝固/熔化模型，通过施加质量源项模拟同步送粉，施加能量源项模拟激光扫描加热，对熔池流场进行了计算，进而获取熔池流场数据。下图是笔者论文中的研究路线。

看了笔者的模型可能会有以下疑问：

为什么不直接用DPM-to-VOF模型，而是拆分为两个模型？

原因是DPM-to-VOF必须是三维模型，考虑到笔者个人电脑的算力，只能妥协使用二维模型，并将激光熔覆过程分步骤计算。

激光热源的施加为什么要通过热源源项，而不直接耦合辐射模型？

在多相流计算中，只有使用欧拉模型才能耦合辐射模型，由于欧拉模型对算力要求大，操作难度高，索性笔者使用VOF模型，并使用源项来施加激光热源。

总之，如果有服务器跑仿真，读者可以不用这些妥协让步，还原更加真实的熔池流场。

三、喷嘴流场模型

3.1 模型假设

基于FLUENT的DPM模型和k-ω SST湍流模型，使用UDF设定壁面边界温度和判定DPM粒子状态，建立气/粉射流仿真计算模型。为提高计算效率做出以下假设：

1. 金属粉末粒子在气/粉射流中所占的体积分数很低，忽略粉末颗粒间的粘性及碰撞，只考虑射流对粒子的单向耦合，同时忽略激光脉冲对金属粉末的冲击作用。
2. 假设金属粉末具有相同的直径0.045mm（325目）且在同轴喷嘴入口处均匀分布。金属材料比热容、导热率、密度不随温度发生变化（比热容：502J/(kg*K)、导热率：20.9W/(m*K)、密度：7800kg/m3）。
3. 颗粒进入熔池中不考虑液体飞溅，生成等离子体等。颗粒与基体或堆积层等区域作用时，均按照弹性碰撞处理，当粒子速度小于逃逸速度，粒子被熔池捕获，反之视为逃逸。

3.2 几何模型和网格划分

模型及网格划分如下图所示，使用Mesh软件的混合网格对流体域进行网格划分，全局网格尺寸0.025mm，壁面网格设置边界层，采用对称网格计算一半流体域以提高仿真效率。喷嘴使用的是课题组研制的光内同轴送粉喷头，读者根具自己的喷嘴建立模型即可。

3.3 模型建立

下面笔者将讲解如何设置Fluent仿真模型，因为太久没有使用Fluent软件了，如有讲错的地方，还请指正。 

（1）通用设置

导入几何模型后，在通用界面求解器类型选择压力基，速度格式选择绝对，时间选择稳态，2D空间选择平面，添加重力。

（2）模型设置

考虑到激光熔覆是一个金属加热融化的过程，激光产生的热量会影响周围空气流动，进而影响喷嘴金属粉末的分布，因此如下图所示，开启能量方程。

如图所示，使用k-ω SST湍流模型，该模型有更广的适用范围和更高的准确精度，尤其是在壁面附近。其参数选择默认设置。

DPM模型设置如下，在喷射源中创建新的group类型，使用惰性粒子，其中流的数量表示仿真粒子的个数，数值越大，仿真的粒子越多，结果也更加准确。仿真材料的参数读者按实际情况设置即可，笔者在本模型中假设金属粉末具有相同的直径0.045mm（325目），且在轴喷嘴入口处均匀分布。

 

（3）材料设置

材料的话，也没有什么特别介绍，就按自己的实际情况设置。

（4）单元区域条件设置

单元域的话默认使用空气即可，也无需做别的更改。

（5）边界条件设置

• 入口设置为速度型入口，气流速度和热量按实际情况设置，DPM中的离散相BC类型选择escape；
• 出口为压力型出口，表压和热量按实际情况设置，DPM中的离散相BC类型选择escape；
• 底部壁面需要添加激光热源，以及判断金属粉末的捕获和逃逸情况，笔者在此展开介绍。

壁面热源的添加可以参考笔者的该文章的介绍。

Fluent UDF教程——壁面温度设定，实现动态高斯热源的施加，DEFINE_PROFILE宏讲解_udf 宏-CSDN博客https://blog.csdn.net/weixin_48501028/article/details/127224390?spm=1001.2014.3001.5502

粉末是否被熔池捕获可以通过比较以下两个力的大小：粉末与基体间碰撞产生的弹力Fr ，由表面张力所引起的Fa ，弹力使粉末离开熔池，黏附力使粉末被熔池吸收。当弹力Fr大于黏附力Fa，粉末将会反弹（reflect）；当弹力Fr小于黏附力Fa，粉末将会被熔池吸收（trap）。

若对于已选定的金属粉末则颗粒粒径是确定的，其他参数也是确定的，所以粉末是否被熔池捕获取决于粉末颗粒的运动速度。

通过宏DEFINE_DPM_BC可以对边界进行设置，采用return返回参数，有以下三种类型：PATH_ACTIVE继续追踪，PATH_ABORT因为错误而终止，PATH_END程序逻辑上的停止。粉末捕获/逃逸的判定程序的逻辑如下图所示：

具体的UDF代码如下：

#include "udf.h"
DEFINE_DPM_BC(bc_reflect,tp,t,f,f_normal,dim)
/*tp被追踪颗粒的数据集f颗粒撞击的面的indexf_normal由颗粒指向壁面单位向量dim维度
*/
{real alpha; 				//颗粒与壁面之间的夹角real vn=0; 				//法向速度real nor_coeff = 1; 			//法向恢复系数real tan_coeff = 0.3; 		//切向恢复系数real normal[3]; 			//存放颗粒指向壁面的法向向量int i, idim = dim; 			//获取维度real NV_VEC(x);			//声明定义一个维度相关的数组可以用for (i = 0; i<idim; i++)normal[i] = f_normal[i]; 	//获取由颗粒指向壁面单位向量if (TP_TYPE(tp) == DPM_TYPE_INERT) {alpha = M_PI / 2. - acos(MAX(-1., MIN(1., NV_DOT(normal, TP_VEL(tp)) / MAX(NV_MAG(TP_VEL(tp)), DPM_SMALL))));/* 判断粒子到达壁面边界 */if ((NNULLP(t)) && (THREAD_TYPE(t) == THREAD_F_WALL)) F_CENTROID(x, f, t);for (i = 0; i<idim; i++) vn += TP_VEL(tp)[i] * normal[i];if(vn<8)		//小于逃逸速度，粒子被捕获return PATH_END; 	//停止计算else{for (i = 0; i<idim; i++)TP_VEL(tp)[i] -= vn*normal[i]; 	//求解切向速度for (i = 0; i<idim; i++)TP_VEL(tp)[i] *= tan_coeff; 	//撞击后切向速度for (i = 0; i<idim; i++)TP_VEL(tp)[i] -= nor_coeff*vn*normal[i]; 	//撞击后法向速度for (i = 0; i<idim; i++)TP_VEL0(tp)[i] = TP_VEL(tp)[i]; 	//更新速度return PATH_ACTIVE; 		//继续追踪}}	
}

（6）求解器设置

求解器的方法和控制使用默认的设置即可：

  

 初始化笔者使用混合初始化，并不需要过多的讲究，运行计算相关的设置如下。

3.4 模型的求解及分析

对于计算域的速度场可以通过以下方式获取，这里笔者就不展开赘述。

 而粒子流场，可以通过Fluent后处理器Sample Trajectory获取基体上表面边界处粒子的数据信息。

之后对数据进行分析处理，最终获得金属粉末落入熔池前的分布曲线。

 四、激光热源模型

在激光熔覆数值模拟过程中，热源模型对模拟结果起着至关重要作用，激光束的形状、能量分布、能量密度都会对成形工艺有重要影响。常用的激光热源的模型有高斯热源模型、半椭球体热源模型、双椭球热源模型、三维椎体热源模型等。

不过本文模型的热源是直接通过热源公式施加。例如高斯热源的表达式如下所示：

五、激光熔覆流场仿真

5.1 相关知识点介绍

5.1.1 多相流模型

VOF(Volume of Fluid)方法是一种在固定的欧拉网格下的界面追踪方法，是通过研究网格单元中流体和网格体积比来确定自由面，追踪流体的变化[[i]]。基本原理为：计算每个控制单元内互不相融的两相的体积占比，计算第一相(Primary Phase)的体积分数α，第二相的体积分数则为1-α。当α为1表示该单元为第一相，为0表示第二相，介于0-1时表示该处为固液界面。

5.1.2 熔化凝固模型

5.2 模型假设

 经过漫长的前期准备，我们总算可以仿真计算流场了！！！

由于激光熔覆是一个复杂的瞬态过程，涉及热传导、对流、熔化凝固等诸多物理化学过程。为减小激光熔覆数值仿真过程中的操作难度和计算量，对激光熔覆作合理的简化假设是很有必要的：

1. 考虑表面张力、重力和浮力对熔池流体流动的驱动作用，忽略准直气流、送粉气流及激光脉冲对熔池和熔覆层的冲击作用，忽略熔池的蒸发作用。

2. 开启熔化凝固模型将计算域分为固相、液相和固液两相区。熔池中的液相视为不可压缩性且为层流流动的液体。固液两相区视为连续多孔介质材料，液相和固相分别视为孔隙度为1和0的多孔介质。

3. 熔覆材料为各向同性，金属的导热系数、定压热容、粘度、表面张力系数随温度变化而改变，其他参数（密度、凝固潜热、固相线温度、液相线温度等）均为定值。

4. 忽略粉末颗粒在落入熔池前，吸收激光束的能量升高的温度，忽略粉末落入熔池时的动能。

5. 设置环境温度和工件的初始温度为293K，大气压力为101325Pa。

5.3 几何模型与网格设置

本文使用Design Modele软件完成二维模型建立。设置计算域尺寸为4mm×3.5mm，初始时计算域划分为空气计算域和316L基体，定义顶部边界为速度型入口，空气域两侧边界为压力型出口，基体计算域边界为壁面。使用Mesh软件的四边型网格进行划分，设置全局网格尺寸0.02mm。模型如下图所示：

5.4 模型建立

（1）通用设置

导入几何模型后，在通用界面求解器类型选择压力基，速度格式选择绝对，时间选择瞬态，2D空间选择平面，添加重力，可以通过设置重力分量，进而模拟基体的倾斜情况。

（2）模型设置

在VOF计算过程中界面模型有Sharp、Sharp/Dispersed、Dispersed三种模式，对于激光熔覆，熔覆材料和空气相界面明确但相界面流动过程非明确，选用Sharp/Dispersed有利于计算收敛以及获得更清晰的相界面。此外VOF模型中有Explicit和Implicit两种计算模式，Explicit更适合瞬态模型的计算，因此本文选择Explicit模式进行仿真。设置第一相(Primary Phase)为316L钢，第二项(Secondary Phase)为空气，开启体积力和表面张力计算。同时，在求解方法中开启VOF限速处理提高计算收敛。

仿真设计激光热源，当然需要开启能量方程。

将熔池的流动假设为层流流动。

开启融化凝固模型，融化凝固模型的本质是控制材料的孔隙率，当孔隙率1视为流体，孔隙率0视为固体。

（3）材料设置

本文基体和金属粉末材料选用316L钢，316L材料性质如下表所示：

 

如下图所示，本文的材料分为固相的和液相的金属材料，两者设置的材料属性并无区别。

 （4）单元区域条件设置

单元域需要开启源项和多孔介区域。

通过宏DEFINE_SOURCE在计算域定义能量源项实现环形激光热源扫描加热：

具体的实现代码如下，笔者使用的为中空激光热源，读者按实际情况修改即可。

# include"udf.h"
# include "sg_mphase.h"
DEFINE_SOURCE(heat_flux,c,t,ds,eqn)
{real xc[ND_ND];real source_heat_air;real x,y,v,time,r,ze,a,p,pi,espl;Thread *pri_th,*sec_th;pri_th=THREAD_SUB_THREAD(t,0); 		//指向首相的指针sec_th=THREAD_SUB_THREAD(t,1); 		//第二相指针C_CENTROID(xc,c,t); 						//获取坐标值v=0.006; 			//扫描速度pi=3.14159; a=0.3; 				//激光吸收率p=800; 				//激光功率r=0.001; 				//光斑半径ze=0.0005; 			//离焦量espl=0; 				//能量峰位置系数time=RP_Get_Real("flow-time");			//获取时间/* 判断熔覆层边界，在边界处施加环形激光热源 */if(C_VOF(c,pri_th)>0.5 && C_VOF(c,pri_th)<0.95 && time<=0.3333 ){source_heat =(2*a*p)/(pi*(r*r+2*r*ze*0.5774))*exp(-2*pow(pow(pow(0.001-v*time,2)+pow(xc[0],2),0.5)-(ze*0.5774+espl*r),2)/(r*r))*6000;ds[eqn]=0;}else{source_heat =0;ds[eqn]=0;}return source_heat;
}

 此外呢，还需要添加质量源项，模拟金属粉末施加的过程，

# include"udf.h"
# include "sg_mphase.h"
DEFINE_SOURCE(den_flux,c,t,ds,eqn)
{real xc[ND_ND];real source_den;real time,fit_line,x,pi,a,p;a = 700000; 			//修正系数p = 7800; 			//金属密度pi = 3.14159265;C_CENTROID(xc,c,t); 		//获取坐标/* 金属粉末落入熔池前的分布曲线 */x = xc[0]*1000;fit_line=(-448.05*pow(x,8)+105.375*pow(x,6)- 9.6225*pow(x,4)+0.41*pow(x,2)+0.5587)*0.001;if(fit_line<0)fit_line = 0;time=RP_Get_Real("flow-time");		//仿真时间/* 在熔覆层边界处施加质量源项 */if(time>=0.0833 && C_VOF(c,t)>0.5 && C_VOF(c,t)<1 && time<0.25){source_den=p*a*fit_line*(1-C_VOF(c,t));ds[eqn]=0;}else{source_den=0;ds[eqn]=0;}return source_den;
}

（5）边界条件设置

以热对流形式定义模型边界处的金属内部换热和空气对流换热，并定义模型顶部与空气接触面存在环境辐射换热。

（6）求解器设置 

由于激光熔覆热影响区间小，加热及冷却时间快，并且本文关注点为熔池流场流动和形貌演变，当温度低于固相线温度时熔池停止流动，因此设置较小的计算域和仿真时间即可满足激光熔覆熔池仿真要求。设置瞬态仿真步长0.0001s，时间步数5000步，计算0.5s的激光熔覆仿真过程。在0s-0.3333s由UDF施加能量源项模拟激光热源扫过截面，0.08333s-0.25s期间施加质量源项模拟同步送粉的过程，0.3333s-0.5s设置为冷却过程。

设置好了之后就可以进行一个漫长又痛苦的等待了，模型参数要是不正确，很有可能会计算不收敛，因此，要有很大的耐心。 

5.5 模型的求解及分析

 ​​​​熔覆层界面变化仿真结果如下：

熔池流场仿真结果如下：

但实际中激光熔覆是一个复杂的三维瞬态问题；并且因为VOF模型的局限性，未对变姿态下熔池流场进行追踪分析。若为进一步完善提高激光熔覆熔池流场仿真的真实性及准确性，可采用DPM-to-VOF建立三维模型。

以上是笔者激光熔覆熔池流场的仿真过程，如有错误还请指正。