[Cooper24]

I am trying to model laser melting of metal particles(Selective Laser Melting). The laser will move over the particles and melt them. There is an inert gas domain above the particles and hence the 2 phase flow VOF is used.
I have created the setup. I have written the UDF for the energy source term for laser. Now, the laser will be applied at the interface of metal particles and the gas. When the laser melts the particles, volume fraction will change and hence the cells in which the laser will be applied. Can someone please help with the UDF for tracking the interfacial cells for the application of laser? Any help is much appreciated.
Till now, I have been using this code but I am not sure if it is correct or not as I am not getting correct results.

#include "udf.h"
#include "sg_mphase.h"
#include "mem.h"
#include "sg_mem.h"
#include "math.h"
#include "flow.h"
#include "unsteady.h"
#include "metric.h"

#define A 0.4 // Absorption coefficient
#define P 200 // Laser power
#define R 40e-6 // spot radius
#define v 0.5 // scan speed of laser
#define h 25 // Heat transfer coefficient
#define Ta 298 // Ambient air temperature
#define s 5.67e-8 // Stefan Boltzmann constant
#define e 0.5 // Emmisivity
#define Pi 3.1415926535
#define Ts 1658 // Solidus temperature
#define Tl 1723 // Liquidus temperature
#define x0 -400e-6 // Initial x position of the laser
#define y0 0.0 // Intiial y position of the laser
#define Lv 7.45e6 // Latent heat of Vaporisation
#define Tv 3090 // Evaporation Temperature
#define Rg 8.314 // Universal Gas constant
#define M 0.05593 // Molar mass
#define Pa 101325 // Atmospheric pressure
#define sigma 1.6 // Surface tension coefficient
#define sigmaT 0.8e-3 // Temperature gradient of surface tension
#define domain_ID 3 // Domain ID of metal substrate

#define rhog 1.662 // density of argon
#define Cpg 520.64 // specific heat of argon



DEFINE_ADJUST(adjust_gradient_heat, domain)
{
Thread *t;
Thread **pt;
cell_t c;
int phase_domain_index = 1.0;
Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,phase_domain_index);
{
Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_N ULL);
Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);
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))
{
Thread *ppt = pt[phase_domain_index];

begin_c_loop (c,t)
{
C_UDMI(c,t,0) = C_VOF_G(c,ppt)[0];
C_UDMI(c,t,1) = C_VOF_G(c,ppt)[1];
C_UDMI(c,t,2) = C_VOF_G(c,ppt)[2];
C_UDMI(c,t,3) = 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)); // magnitude of gradient of volume fraction
C_UDMI(c,t,4) = C_UDMI(c,t,0)/C_UDMI(c,t,3); // nx -> x gradient of volume fraction divided by magnitude of gradient of volume fraction
C_UDMI(c,t,5) = C_UDMI(c,t,1)/C_UDMI(c,t,3); // ny -> y gradient of volume fraction divided by magnitude of gradient of volume fraction
C_UDMI(c,t,6) = C_UDMI(c,t,2)/C_UDMI(c,t,3); // nz -> z gradient of volume fraction divided by magnitude of gradient of volume fraction
}
end_c_loop (c,t)
}
Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NU LL);
}

DEFINE_SOURCE(heat_source, c, t, dS, eqn) // The name of the UDF is heat_source
{
Thread *pri_th;
Thread *sec_th;

real source;
real x[ND_ND], time; // Define face centroid vector, time
time = RP_Get_Real("flow-time"); // Acquire time from Fluent solver
C_CENTROID(x, c, t); // Acquire the cell centroid location
real T = C_T(c,t);

pri_th = THREAD_SUB_THREAD(t, 0);
sec_th = THREAD_SUB_THREAD(t, 1);

real alpha = C_VOF(c,sec_th); // cell volume fraction
real gamma;

if (T>293.0 && T<1658.0)
{
gamma = 0;
}
else if (T>=1658.0 && T<=1723.0)
{
gamma = (T-Ts)/(Tl-Ts);
}
else if (T>1723.0)
{
gamma = 1;
}


real Pv = 0.54*Pa*exp((Lv*M*(T-Tv))/(Rg*T*Tv));
real mv = (0.82*M*Pv)/(sqrt(2*Pi*M*Rg*T));

real rhos = 7900; // density of solid SS316
real rhol = 7433 + 0.0393*T - 0.00018*pow(T,2); // density of liquid SS316
real rhom = rhol*gamma + rhos*(1-gamma); // density of SS316
real rho = alpha*rhom + rhog*(1-alpha); // density of cell containing metal and gas
real Cps = 462 + 0.134*T; // specific heat of solid SS316
real Cpl = 775; // specific heat of liquid SS316
real Cpm = Cpl*gamma + Cps*(1-gamma); // specific heat of SS316
real Cp = alpha*Cpm + Cpg*(1-alpha); // specific heat of cell containing metal and gas

real factor = (2*rho*Cp)/(rhom*Cpm + rhog*Cpg);

real r = sqrt(pow(x[0]-x0-v*time,2.0) + pow(x[1]-y0,2.0));

if(C_VOF(c,sec_th)>0.05 && C_VOF(c,sec_th)<1)
{
if(C_T(c,t) < 3090)
{
source = (((2*A*P)/(Pi*R*R))*exp((-2*(r*r))/(R*R)) - h*(T-Ta) - s*e*(pow(T,4) - pow(Ta,4)))*C_UDMI(c,t,3)*factor;
dS[eqn] = 0.0;
}
else if(C_T(c,t) >= 3090)
{
source = (((2*A*P)/(Pi*R*R))*exp((-2*(r*r))/(R*R)) - h*(T-Ta) - s*e*(pow(T,4) - pow(Ta,4)) - Lv*mv)*C_UDMI(c,t,3)*factor;
dS[eqn] = 0.0;
}
}
else
{
source = 0.0;
dS[eqn] = 0.0;
}
return source;
}