[Cooper24]

I am trying to model laser melting and solidification of a block. I am using VOF with solidification/melting in Fluent. I have written a UDF for energy source. The image attached shows the equation I am trying to incorporate. The other image shows the domain I am using and BCs. My UDF compiles without any error but when I try to run the simulation, Fluent crashes and window closes. I am new to UDF writing and may committing some mistake.
I would be really helpful if anyone could provide some insight.
Here is my UDF

#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 80e-6 // spot radius
#define v 0.1 // 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 -100e-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.005593 // Molar mass
#define Pa 101325 // Atmospheric pressure
#define domain_ID 3 // Domain ID of metal substrate


DEFINE_INIT(volume_fraction, mixture_domain)
{
int phase_domain_index = 1;
cell_t c;
Thread *t;
Domain *subdomain;
real xc[ND_ND];
sub_domain_loop(subdomain, mixture_domain, phase_domain_index)
{
if(DOMAIN_ID(subdomain) == 3)
thread_loop_c(t,subdomain)
{
begin_c_loop_all(c,t)
{
C_CENTROID(xc,c,t);
if(xc[0] > -0.5e-3 && xc[0] < 0.5e-3 && xc[1] > -0.25e-3 && xc[1] < 0.25e-3 && xc[2] < 0 && xc[2] > -0.3e-3)
C_VOF(c,t) = 1.0;
else
C_VOF(c,t) = 0.0;
}
end_c_loop_all(c,t)
}
}
}



DEFINE_ADJUST(adjust_gradient, 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));
}
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);
real alpha = C_VOF(c,t); // cell volume fraction
real gamma = C_LIQF(c,t); // cell liquid fraction

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 rhog = 1.6228; // density of argon
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 Cpg = 520.64; // specific heat of argon
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);


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

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*(pow(x[0]-x0-v*time,2.0) + pow(x[1]-y0,2.0)))/(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*(pow(x[0]-x0-v*time,2.0) + pow(x[1]-y0,2.0)))/(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;
}

Maybe there is some issue with the gradient of volume fraction term and the smearing out factor multiplied at the end. Any help is much appreciated.
Thanks in advanceHeat_Source_Eqn.jpg

Domain.png

[AlexanderZ]

my compiler show the following error

line(98) : error C2275: 'real' : illegal use of this type as an expression

[Cooper24]

Thank you very much for your reply. I am not very good at coding and new to writing UDF. I directly compiled it in Fluent and it compiled without any error. But when I tried to run it, Fluent crashed with SIGSEGV fault.
I am really relieved that you replied as I am really facing difficulty with this being new to the field.
Can you suggest if the UDF is right for modeling moving laser source and storing gradient of volume fraction? And how to check validity of the code even if it compiles in Fluent without error?

Thanks in advance.

[pakk]

If you are new to writing UDFs, maybe you got the same error but did not see it because you don't know where to look... See my signature below.

[AlexanderZ]

user-defined-> press build -> copy here console output

[Cooper24]

Thank you for your reply. The UDF compiles without error and after I will share the screenshot of the console. The issue I am facing is that the "factor" which I multiplied in the source term at the last is creating some problem as you can seen in my original post. It's a factor to be multiplied as shown in the equation I am tried to add in the image. That factor depends on the property of two phases and volume fraction and liquid fraction in the cell.
If I remove the factor, the following UDF runs fine and simulation also starts although it starts diverging after some time.
This is the working UDF without the factor:

#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 80e-6 // spot radius
#define v 0.1 // 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 -100e-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.005593 // Molar mass
#define Pa 101325 // Atmospheric pressure
#define domain_ID 3 // Domain ID of metal substrate


DEFINE_INIT(volume_fraction, mixture_domain)
{
int phase_domain_index = 1;
cell_t c;
Thread *t;
Domain *subdomain;
real xc[ND_ND];
sub_domain_loop(subdomain, mixture_domain, phase_domain_index)
{
if(DOMAIN_ID(subdomain) == 3)
thread_loop_c(t,subdomain)
{
begin_c_loop_all(c,t)
{
C_CENTROID(xc,c,t);
if(xc[0] > -0.5e-3 && xc[0] < 0.5e-3 && xc[1] > -0.25e-3 && xc[1] < 0.25e-3 && xc[2] < 0 && xc[2] > -0.3e-3)
C_VOF(c,t) = 1.0;
else
C_VOF(c,t) = 0.0;
}
end_c_loop_all(c,t)
}
}
}



DEFINE_ADJUST(adjust_gradient, 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));
}
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);

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));

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

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*(pow(x[0]-x0-v*time,2.0) + pow(x[1]-y0,2.0)))/(R*R)) - h*(T-Ta) - s*e*(pow(T,4) - pow(Ta,4)))*C_UDMI(c,t,3);
dS[eqn] = 0.0;
}
else if(C_T(c,t) > 3090)
{
source = (((2*A*P)/(Pi*R*R))*exp((-2*(pow(x[0]-x0-v*time,2.0) + pow(x[1]-y0,2.0)))/(R*R)) - h*(T-Ta) - s*e*(pow(T,4) - pow(Ta,4)) - Lv*mv)*C_UDMI(c,t,3);
dS[eqn] = 0.0;
}
}
else
{
source = 0.0;
dS[eqn] = 0.0;
}
return source;
}


And this is what I want to add in the source term:

real alpha = C_VOF(c,t); // cell volume fraction
real gamma = C_LIQF(c,t); // cell liquid fraction

real rhog = 1.6228; // density of argon
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 Cpg = 520.64; // specific heat of argon
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);

This "factor", I want to multiply in the source term at last.

[Cooper24]

I will share the screenshot

[Cooper24]

This is the message I get when I build the UDF