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

本文主要是介绍基于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粒子状态，建立气/粉射流仿真计算模型。为提高计算效率做出以下假设：

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

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 熔化凝固模型

FLUENT的熔化凝固模型应用了焓-孔隙度技术，将计算域视为多孔介质材料。多孔介质是指内部含有众多空隙的固体材料，当孔隙率1为流体，孔隙率0为固体，孔隙率介于0-1可以认为是糊状区。由此可通过设置单元域条件开启多孔区域选项模拟基体的凝固与熔化。

5.2 模型假设

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

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

考虑表面张力、重力和浮力对熔池流体流动的驱动作用，忽略准直气流、送粉气流及激光脉冲对熔池和熔覆层的冲击作用，忽略熔池的蒸发作用。
开启熔化凝固模型将计算域分为固相、液相和固液两相区。熔池中的液相视为不可压缩性且为层流流动的液体。固液两相区视为连续多孔介质材料，液相和固相分别视为孔隙度为1和0的多孔介质。
熔覆材料为各向同性，金属的导热系数、定压热容、粘度、表面张力系数随温度变化而改变，其他参数（密度、凝固潜热、固相线温度、液相线温度等）均为定值。
忽略粉末颗粒在落入熔池前，吸收激光束的能量升高的温度，忽略粉末落入熔池时的动能。
设置环境温度和工件的初始温度为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设置为冷却过程。