UDF for tracking interface in 2 phase VOF method.
 

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

It's easier to help if you say which results you expected, and which results you got.

Hello Sir, thanks for your reply. Actually, I am new to UDF writing and wrote this code taking help from different sources.
Now, there are 2 things I need:
1. Gradient of Volume Fraction to convert surface sources to volumetric sources. To get volume fraction gradient, I have used the above DEFINE_ADJUST which I found in the forum. Hope it is right.
2. Next is I want to apply my source terms at the interfacial cells. For this I am taking a range of volume fraction, 0.05 < C_VOF < 1, for the secondary phase. The cells having this volume fraction of secondary phase are to be considered the interfacial cells.
Now, I want to write a UDF to identify these interfacial cells. From, what I have learnt it can be done by comparing the z coordinates of the cell centroids and the source will be applied in the cells with z coordinates in which C_VOF 1. I am unable to figure out how to implement this.
Please correct me if I am wrong anywhere.

Maybe my question was not very clear.

You wrote the following:
So you get incorrect results, in other words something in your results was different than what you expected.

My question: what is that something?

What did you see in the results that was different than what you expected? Temperature? Pressure? Velocities? How did you conclude that your results were incorrect?

Sorry for the ambiguous reply. Actually the laser should be applied only at the interfacial cells and the other source terms also. So, the laser should only move at the top surface of the solid block. Now, during the post processing, the laser is moving, the contours are right and the temperature values are also in the suitable range.
The issue is that the laser is applied continuously the the thickness of the block. When I create a plane below the top surface, I can see a laser moving on that plane. I have attached 2 images also.
Now, as I have applied the sources using the condition for volume fraction of secondary phase, C_VOF(c,sec_th), I think these cells should be at the interface only and that's what I need help with.

I have attached the image of the particle bed which I want to model actually. But before that I have been trying just with a solid block. I want to move the laser over the particles as you can see in the image.

The "real" solution is probably very complex, but there is a simple hack that might work...

if(C_VOF(c,sec_th)>0.05 && C_VOF(c,sec_th)<1)

Change that to

if(C_VOF(c,sec_th)>0.05 && C_VOF(c,sec_th)<0.95)

And you might have a "good enough" result.

You might need to play with the values (0.05 and 0.95) until it looks like what you want.

Thanks for your suggestion sir. I have read in a paper how to track the interfacial cells. It was written that z coordinates of the cell centroids should be stored along with the cell volume fraction. Next the cells which have a volume fraction less than 1, their z coordinates should be checked and if the cells below these have a volume fraction = 1, then these are the interfacial cells.
Can you suggest how to implement this? I know that I can store cell centroids and cell volume fractions in UDMs using DEFINE_ADJUST. I am unable to think of the comparison of z coordinates part.

That is the complex "real" solution that I mentioned.
I could probably implement it, but it would take me a lot of time. I would reserve 2 work days for this including testing (excluding documenting). That's more time than I'm willing to invest in a stranger's problem, I'm afraid.

So the above "hack" is the best you get from me. I hope it's useful.

Thanks for your reply sir. I completely understand and I would not want you to invest your time in my problem. Just wanted to check if that is the correct way. I understand that it is complex and will see if I am able to do it.

Hello Sir, taking help from different sources, I tried to write a UDF for tracking the interfacial cells. The UDF compiles without any error but the simulation throws segmentation fault and FLUENT crashes at the first iteration itself. Can you please tell if my approach is right or not? I need this very badly and would be really grateful. Here is my UDF for tracking interfacial cells:

DEFINE_ADJUST(adjust_interface, domain)
{
Thread *t;
Thread **pt;
cell_t c;
face_t f;
face_t ftop;
cell_t c0;
cell_t c1;
cell_t ca;
cell_t cb;

int phase_domain_index = 1;

real zunit[ND_ND];
real ZDOT;

real es[ND_ND];
real ds;
real A[ND_ND];
real A_by_es;
real dr0[ND_ND];
real dr1[ND_ND];
mp_thread_loop_c(t,domain,pt)
if (FLUID_THREAD_P(t))
{
Thread *ppt = pt[phase_domain_index];

begin_c_loop(c,t) // This will loop over all the cells in the domain
{
real xc[ND_ND];
C_CENTROID(xc,c,t);
NV_D(zunit, = ,0,0,1);

int n;
int top;
C_UDMI(c,t,7) = xc[2];
C_UDMI(c,t,8) = C_VOF(c,ppt);

c_face_loop(c, t, n) // This will loop over all the faces of a cell and n is the local face index number
{
f = C_FACE(c,t,n);
INTERIOR_FACE_GEOMETRY(f,t,A,ds,es,A_by_es,dr0,dr1 );
ZDOT = NV_DOT(es,zunit);
ca = F_C1(f,t);
cb = F_C0(f,t);

if(ZDOT == 1 && 0.05 < C_VOF(cb,ppt) < 1 && C_VOF(ca,ppt) == 0)
{
top = n;
}

ftop = C_FACE(c,t,top); // ftop is the global face index number. C_FACE converts local face index number, n to global face index number, ftop
}

c1 = F_C1(ftop,t); // c1 is the index of face's neighbouring C1 cell
c0 = F_C0(ftop,t); // c0 is the index of face's neighbouring C0 cell

/* This is the volume fraction condition */
if (C_VOF(c1,ppt) = 0.0 && 0.05 < C_VOF(c0,ppt) < 1)
{
C_UDMI(c,t,9) = c0;
C_UDMI(c,t,10) = c1;
}
}
end_c_loop(c,t)
}
}


I tried to print lines also for checking and the error comes at the line If(FLUID_THREAD_P(t))

That's probably the mistake. "if(ZDOT == 1 && 0.05 < C_VOF(cb,ppt) < 1 && C_VOF(ca,ppt) == 0)"

It maybe so "if(ZDOT == 1 && C_VOF(cb,ppt)>0.05&&C_VOF(cb,ppt) < 1 && C_VOF(ca,ppt) == 0)"

you can use this expression as a condition in C, check syntax

Thanks for your help. I will change that. Although the UDF compiles without any error. The issue is when I start the simulation. It throws a mpt_accept error. I can change the syntax but I am more concerned about the logic I am using to track the interfacial cells. I am new to writing UDFs and notsure whether my logic is right.

Thanks for your help. I will change that although it compiles without any error. The issue is I am not confident about the logic I have used for tracking the interfacial cells. I am new to writing UDFs and even if the logic is right, I am not confident whether I have written the code correctly. Pleae give some advice in this regard if you can.

Thanks

Were you aware that you need to allocate at least 11 UDMs for your code to work?

Thanks for your reply. Yes I was aware of that and infact I have allocated 12 UDMs in fluent.