CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > ANSYS > FLUENT > Fluent UDF and Scheme Programming

UDF for tracking interface in 2 phase VOF method.

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 25, 2021, 13:43
Default UDF for tracking interface in 2 phase VOF method.
  #1
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
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;
}
Cooper24 is offline   Reply With Quote

Old   August 26, 2021, 05:54
Default
  #2
Senior Member
 
Join Date: Nov 2013
Posts: 1,926
Rep Power: 23
pakk will become famous soon enough
It's easier to help if you say which results you expected, and which results you got.
__________________
"The UDF library you are trying to load (libudf) is not compiled for parallel use on the current platform" is NOT the error after compiling. It is the error after loading. To see compiler errors, look at your screen after you click "build".
pakk is offline   Reply With Quote

Old   August 26, 2021, 06:59
Default
  #3
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
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.
Cooper24 is offline   Reply With Quote

Old   August 26, 2021, 12:03
Default
  #4
Senior Member
 
Join Date: Nov 2013
Posts: 1,926
Rep Power: 23
pakk will become famous soon enough
Maybe my question was not very clear.

You wrote the following:
Quote:
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.
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?
__________________
"The UDF library you are trying to load (libudf) is not compiled for parallel use on the current platform" is NOT the error after compiling. It is the error after loading. To see compiler errors, look at your screen after you click "build".
pakk is offline   Reply With Quote

Old   August 26, 2021, 12:47
Default
  #5
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
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.
Attached Images
File Type: png keyhole3.png (112.9 KB, 8 views)
File Type: png FFF2.png (145.5 KB, 6 views)
Cooper24 is offline   Reply With Quote

Old   August 26, 2021, 12:52
Default
  #6
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
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.
Attached Images
File Type: jpg particlebed1.jpg (83.9 KB, 4 views)
Cooper24 is offline   Reply With Quote

Old   August 26, 2021, 13:18
Default
  #7
Senior Member
 
Join Date: Nov 2013
Posts: 1,926
Rep Power: 23
pakk will become famous soon enough
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.
__________________
"The UDF library you are trying to load (libudf) is not compiled for parallel use on the current platform" is NOT the error after compiling. It is the error after loading. To see compiler errors, look at your screen after you click "build".
pakk is offline   Reply With Quote

Old   August 26, 2021, 13:38
Default
  #8
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
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.
Cooper24 is offline   Reply With Quote

Old   August 26, 2021, 14:03
Default
  #9
Senior Member
 
Join Date: Nov 2013
Posts: 1,926
Rep Power: 23
pakk will become famous soon enough
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.
__________________
"The UDF library you are trying to load (libudf) is not compiled for parallel use on the current platform" is NOT the error after compiling. It is the error after loading. To see compiler errors, look at your screen after you click "build".
pakk is offline   Reply With Quote

Old   August 26, 2021, 14:29
Default
  #10
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
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.
Cooper24 is offline   Reply With Quote

Old   October 22, 2021, 12:47
Default
  #11
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
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))
Cooper24 is offline   Reply With Quote

Old   October 24, 2021, 22:01
Default
  #12
Member
 
xingangzheng
Join Date: Jul 2009
Posts: 30
Rep Power: 14
ustbdynamic is on a distinguished road
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)"



Quote:
Originally Posted by Cooper24 View Post
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))
ustbdynamic is offline   Reply With Quote

Old   October 25, 2021, 00:38
Default
  #13
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 1,788
Rep Power: 28
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
you can use this expression as a condition in C, check syntax
Code:
 C_VOF(cb,ppt)>0.05&&C_VOF(cb,ppt) < 1
__________________
best regards


******************************
press LIKE if this message was helpful
AlexanderZ is offline   Reply With Quote

Old   October 25, 2021, 11:47
Default
  #14
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
Quote:
Originally Posted by ustbdynamic View Post
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)"
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.
Cooper24 is offline   Reply With Quote

Old   October 25, 2021, 11:49
Default
  #15
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
Quote:
Originally Posted by AlexanderZ View Post
you can use this expression as a condition in C, check syntax
Code:
 C_VOF(cb,ppt)>0.05&&C_VOF(cb,ppt) < 1
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
Cooper24 is offline   Reply With Quote

Old   October 25, 2021, 14:15
Default
  #16
Senior Member
 
Join Date: Nov 2013
Posts: 1,926
Rep Power: 23
pakk will become famous soon enough
Were you aware that you need to allocate at least 11 UDMs for your code to work?
__________________
"The UDF library you are trying to load (libudf) is not compiled for parallel use on the current platform" is NOT the error after compiling. It is the error after loading. To see compiler errors, look at your screen after you click "build".
pakk is offline   Reply With Quote

Old   Yesterday, 05:42
Default
  #17
Member
 
Join Date: Jul 2020
Location: India
Posts: 36
Rep Power: 3
Cooper24 is on a distinguished road
Thanks for your reply. Yes I was aware of that and infact I have allocated 12 UDMs in fluent.
Cooper24 is offline   Reply With Quote

Reply

Tags
interface tracking, laser melting, volume of fluid

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Phase Field Method & Diffuse Interface Modeling in FOAM holger_marschall OpenFOAM Programming & Development 13 December 19, 2019 00:54
How to write my own phase change with VOF method yamifm0f STAR-CCM+ 0 September 27, 2018 07:17
Question about adaptive timestepping Guille1811 CFX 25 November 12, 2017 17:38
two phase model with CICSAM vof method hilllike Main CFD Forum 6 February 21, 2014 14:31
Radiation interface hinca CFX 15 January 26, 2014 17:11


All times are GMT -4. The time now is 12:47.