Fluent案例：利用VOF+凝固熔化+UDF实现激光焊熔池模拟
案例描述
利用VOF+凝固熔化+UDF对激光焊熔池过程的模拟，该过程中包括金属材料的相变、激光热源加热、热源移动等物理变化过程。通过UDF完整地模拟了激光焊接过程中的热源加热、熔池流动、反冲压力效应、表面张力变化和浮力效应等关键物理现象。构建模型如下图所示，模型的正面为对称面，其他均为壁面，焊接的起始点为（-3.5,0,0），激光热源的移动速度为0.01m/s（文末附案例文件获取方式）。
步骤
1.导入网格模型
导入模型网格文件，如下图所示：
2.通用设置
采用压力基求解器，瞬态求解，重力方向为z的负方向。
3. 模型选择
开启能量方程，粘性模型选择层流，开启凝固和熔化模型并保持默认设置。
4. UDF加载
添加udf.c源文件，点击编译并加载。
UDF完整地模拟了激光焊接过程中的热源加热、熔池流动、反冲压力效应、表面张力变化和浮力效应等关键物理现象。注意：UDF定义焊接速度的过程中，速度的定义并非热源点的移动速度，而是计算域的移动速度。
// 定义常数宏//激光热源函数DEFINE_SOURCE(laser_source,c,t,dS,eqn)// 定义热源项，Fluent UDF 标准格式{real xc[ND_ND],x,y,z,time,Q,r,H,source; // 声明变量：坐标、时间、功率、半径、深度、源项Thread *sec_th; // 声明第二相线程指针sec_th=THREAD_SUB_THREAD(t,1); // 获取第二相线程（多相流中的气相）time=RP_Get_Real("flow-time"); // 获取当前流动时间C_CENTROID(xc, c, t); // 获取当前网格单元的中心坐标x=xc[0]+xx0-v*time; // 计算移动热源的x坐标（考虑焊接速度，速度的定义并非热源点的移动速度，而是计算域的移动速度）y=xc[1]; // y坐标z=xc[2]; // z坐标H=0.003; /* 热源深度*/ // 定义热源作用深度r=sqrt(x*x+y*y); // 计算点到热源中心的水平距离Q=n*P; // 计算有效功率 Q = 效率 × 功率if(r<=R0&&z<=0&&z>=-H) // 判断条件：在热源半径内且在深度范围内{// 高斯旋转体热源模型source=9*Q/(PI*R0*R0*H*0.95)*exp(-6/(log(H/(-z))*R0*R0)*r*r);dS[eqn]=0; // 源项对变量的导数为0（常数源项）}else{source=0; // 在热源范围外，源项为0dS[eqn]=0; // 导数为0}return source; // 返回热源项值}//VOF梯度存储函数DEFINE_ADJUST(store_VOF_gradient,domain)// 定义调整函数，在每个时间步执行{Thread *t; // 线程指针Thread *ppt; // 主相线程指针Thread **pt; // 线程指针数组cell_t c; // 网格单元标识符int phase_domain_index=1; // 相域索引Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,phase_domain_index); // 获取子域// 分配存储空间用于VOF梯度计算Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);// 重构VOF标量场Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);// 计算VOF梯度Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG, Vof_Deriv_Accumulate);// 循环遍历所有线程mp_thread_loop_c (t,domain,pt){if (FLUID_THREAD_P(t)) // 如果是流体线程{ppt = pt[phase_domain_index]; // 获取主相线程begin_c_loop (c,t) // 开始单元循环{// 将VOF梯度分量存储到UDM中C_UDMI(c,t,0) = C_VOF_G(c,ppt)[0]; // x方向梯度C_UDMI(c,t,1) = C_VOF_G(c,ppt)[1]; // y方向梯度C_UDMI(c,t,2) = C_VOF_G(c,ppt)[2]; // z方向梯度}end_c_loop (c,t) // 结束单元循环}}// 释放存储空间Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);}//X方向反冲压力源项DEFINE_SOURCE(x_recoil_source,c,t,dS,eqn){real xc[ND_ND],x,y,z,r,time,T,xfen,recoil,x_source; // 声明变量Thread *sec_th;sec_th=THREAD_SUB_THREAD(t,1); // 获取第二相线程T=C_T(c,t); // 获取当前单元温度time=RP_Get_Real("flow-time"); // 获取当前时间C_CENTROID(xc, c, t); // 获取单元中心坐标x=xc[0]+xx0-v*time; // 计算移动坐标y=xc[1];z=xc[2];r=sqrt(x*x+y*y); // 计算到热源中心距离// 计算x方向界面法向分量（归一化）xfen=C_UDMI(c,t,0)/sqrt(C_UDMI(c,t,0)*C_UDMI(c,t,0)+C_UDMI(c,t,1)*C_UDMI(c,t,1)+C_UDMI(c,t,2)*C_UDMI(c,t,2)+0.00000001);if(r<=0.001) // 在反冲压力作用半径内{// 条件：温度达到汽化点且在气液界面区域if(T>=T_SAT&&C_VOF(c,sec_th)>=0.1&&C_VOF(c,sec_th)<=0.9) //反冲压力必须施加在第二相上，0表示第一相-空气相{// 反冲压力计算（基于温度的经验公式）recoil=3.35*Re*exp(12.37*(T-T_SAT)/T);x_source=recoil*xfen; // x方向反冲压力分量C_UDMI(c,t,3)=x_source; // 存储到UDM3}elsex_source=0; // 不满足条件时源项为0}else{x_source=0; // 在作用半径外源项为0}return x_source*Ree; // 返回x方向反冲压力源项}//Y方向反冲压力源项DEFINE_SOURCE(y_recoil_source,c,t,dS,eqn){real xc[ND_ND],x,y,z,r,time,T,yfen,recoil,y_source;Thread *sec_th;sec_th=THREAD_SUB_THREAD(t,1);T=C_T(c,t);time=RP_Get_Real("flow-time");C_CENTROID(xc, c, t);x=xc[0]+xx0-v*time;y=xc[1];z=xc[2];r=sqrt(x*x+y*y);// 计算y方向界面法向分量yfen=C_UDMI(c,t,1)/sqrt(C_UDMI(c,t,0)*C_UDMI(c,t,0)+C_UDMI(c,t,1)*C_UDMI(c,t,1)+C_UDMI(c,t,2)*C_UDMI(c,t,2)+0.00000001);if(r<=0.001){if(T>=T_SAT&&C_VOF(c,sec_th)>=0.1&&C_VOF(c,sec_th)<=0.9){recoil=3.35*Re*exp(12.37*(T-T_SAT)/T);y_source=recoil*yfen; // y方向反冲压力分量C_UDMI(c,t,4)=y_source; // 存储到UDM4}elsey_source=0;}else{y_source=0;}return y_source*Ree; // 返回y方向反冲压力源项}//Z方向反冲压力源项DEFINE_SOURCE(z_recoil_source,c,t,dS,eqn){real xc[ND_ND],x,y,z,r,time,T,zfen,recoil,z_source;Thread *sec_th;sec_th=THREAD_SUB_THREAD(t,1);T=C_T(c,t);time=RP_Get_Real("flow-time");C_CENTROID(xc, c, t);x=xc[0]+xx0-v*time;y=xc[1];z=xc[2];r=sqrt(x*x+y*y);// 计算z方向界面法向分量zfen=C_UDMI(c,t,2)/sqrt(C_UDMI(c,t,0)*C_UDMI(c,t,0)+C_UDMI(c,t,1)*C_UDMI(c,t,1)+C_UDMI(c,t,2)*C_UDMI(c,t,2)+0.00000001);if(r<=0.001){if(T>=T_SAT&&C_VOF(c,sec_th)>=0.1&&C_VOF(c,sec_th)<=0.9){recoil=3.35*Re*exp(12.37*(T-T_SAT)/T);z_source=recoil*zfen; // z方向反冲压力分量C_UDMI(c,t,5)=z_source; // 存储到UDM5}elsez_source=0;}else{z_source=0;}return z_source*Reez; // 返回z方向反冲压力源项（不同系数）}//表面张力属性定义DEFINE_PROPERTY(surf_tension,c,t){real surf;int curr,i,j;real temp1=C_T(c,t); // 获取当前单元温度if(temp1>=930) // 温度 ≥ 930Ksurf=1.0-0.0003*(temp1-930); // 线性减小else if(temp1>=2730) // 温度 ≥ 2730Ksurf=0.41; // 保持最小值elsesurf=1.0; // 低温时保持最大值return surf; // 返回表面张力系数}//浮力源项DEFINE_SOURCE(fuli, cell, thread, dS, eqn){real source;real p, g, B, T, Tl;T = C_T(cell, thread); // 获取温度p = 2700; // 密度g = -9.81; // 重力加速度B = 0.0001; // 热膨胀系数Tl = 930; // 液相线温度if (T > Tl) // 温度高于液相线source = -p * g*B*(T - Tl); // 计算浮力源项elsesource = 0; // 固态时无浮力dS[eqn] = 0; // 导数为0return source; // 返回浮力源项}
5. 材料参数设置
添加并定义金属材料，在流体材料中修改空气材料的参数并另存为金属材料，再对相关参数进行修改。其中密度、比热、导热率、粘度参数为分段线性函数。
其中密度、比热、导热率、粘度参数定义如下：
①密度根据温度定义4个点：300K，2700kg/m3；820K，2700kg/m3；930K，2400kg/m3；5000K，2400kg/m3；
②比热根据温度定义4个点：300K，871J/(kg K)；820K，871J/(kg K)；930K，1060J/(kg K)；5000K，1060J/(kg K)；
③导热率根据温度定义4个点：300K，238W/(m K)；820K，238W/(m K)；930K，100W/(m K)；5000K，100W/(m K)；
④粘度根据温度定义2个点：930K，0.0016kg/(m s)；5000K，0.0016kg/(m s)。
添加金属材料后，定义VOF模型，表面张力系数选择UDF中定义的surf_tension。
6.单元区域条件设置
在单元区域设置中勾选源项，并分别添加源项的数量，包括XYZ三个方向上的反冲压力源项、Z方向上的浮力源项和一个能量源项。
7.边界条件设置
定义底部壁面的传热边界条件：
定义四周壁面的传热边界条件：
8.初始化
新建一个区域标记，用于后续的局部初始化。
初始化选择标准初始化，温度为300K，第二相的体积为0。局部初始化将标记出的区域第二相的体积分数设置为1。
9.运行计算
时间步数和时间步长定义如下：
温度分布云图
空气相体积分布云图
