
[Sponsors] 
UDF for tracking interface in 2 phase VOF method. 

LinkBack  Thread Tools  Search this Thread  Display Modes 
August 25, 2021, 13:43 
UDF for tracking interface in 2 phase VOF method.

#1 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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 40e6 // 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.67e8 // Stefan Boltzmann constant #define e 0.5 // Emmisivity #define Pi 3.1415926535 #define Ts 1658 // Solidus temperature #define Tl 1723 // Liquidus temperature #define x0 400e6 // 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.8e3 // 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("flowtime"); // 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 = (TTs)/(TlTs); } else if (T>1723.0) { gamma = 1; } real Pv = 0.54*Pa*exp((Lv*M*(TTv))/(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*(1gamma); // density of SS316 real rho = alpha*rhom + rhog*(1alpha); // 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*(1gamma); // specific heat of SS316 real Cp = alpha*Cpm + Cpg*(1alpha); // specific heat of cell containing metal and gas real factor = (2*rho*Cp)/(rhom*Cpm + rhog*Cpg); real r = sqrt(pow(x[0]x0v*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*(TTa)  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*(TTa)  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; } 

August 26, 2021, 05:54 

#2 
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 
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". 

August 26, 2021, 06:59 

#3 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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. 

August 26, 2021, 12:03 

#4  
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 
Maybe my question was not very clear.
You wrote the following: Quote:
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". 

August 26, 2021, 12:47 

#5 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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. 

August 26, 2021, 12:52 

#6 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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.


August 26, 2021, 13:18 

#7 
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 
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". 

August 26, 2021, 13:38 

#8 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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. 

August 26, 2021, 14:03 

#9 
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 
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". 

August 26, 2021, 14:29 

#10 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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.


October 22, 2021, 12:47 

#11 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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)) 

October 24, 2021, 22:01 

#12  
Member
xingangzheng
Join Date: Jul 2009
Posts: 37
Rep Power: 16 
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:


October 25, 2021, 00:38 

#13 
Senior Member
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34 
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 

October 25, 2021, 11:47 

#14 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
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.


October 25, 2021, 11:49 

#15  
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
Quote:
Thanks 

October 25, 2021, 14:15 

#16 
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 
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". 

October 26, 2021, 05:42 

#17 
Member
Join Date: Jul 2020
Location: India
Posts: 65
Rep Power: 5 
Thanks for your reply. Yes I was aware of that and infact I have allocated 12 UDMs in fluent.


Tags 
interface tracking, laser melting, volume of fluid 
Thread Tools  Search this Thread 
Display Modes  


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  STARCCM+  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 