
[Sponsors] 
May 17, 2021, 07:31 
Obtaining fluid properties using UDF

#1 
New Member
Join Date: Jan 2018
Posts: 29
Rep Power: 5 
Hello ,
I am trying to run one of those cases, where the properties from outlet from few surfaces are assigned as inlet to another. In this setup, I intend to calculate total pressure and velocity on the surface, and assigned to another surface with the change in friction factor value. However I am not sure if I have set the the total pressure and velocity with the correct formula. The massflow rates have been giving values in the correct neighborhood. however the other two are currently problemsome. TLR Could anyone point out, how can I obtain the total pressure and velocity at a surface using UDF? As I am using a constant density assumption, I am not obtaining this with F_R(f,t). Dia, mu, L are also constants. Code:
DEFINE_PROFILE(pressure_outlet, thread, position) { /* Averages static pressure on downstream interface inflow BC and applies it at the upstream outflow interface */ int ID; face_t f; real NV_VEC(area); real Amagnitude, pressure, press_avg; real A_sum, A_glob, velocity_glob, density_sum, density_glob ; real mdot_sum, mdot_glob, Flux, mflux_glob , reynolds, velocity, velocity_sum, friction_factor; real PA_glob; face_t f_upstream; Domain* domain_upstream; domain_upstream = Get_Domain(1); real density = 0.9462; int zone_ID = THREAD_ID(thread); /* set the ID of the downstream interface */ if (zone_ID == out1) { ID = in1; } else if (zone_ID == out2) { ID = in2; } else if (zone_ID == out3) { ID = in3; } else if (zone_ID == out4) { ID = in4; } else { ID = 10000; } if (ID!=10000) { Thread* thread_upstream = Lookup_Thread(domain_upstream, ID); mdot_sum = 0.0; A_sum = 0.0; velocity_sum = 0.0; reynolds = 0.0; begin_f_loop(f_upstream, thread_upstream) if PRINCIPAL_FACE_P(f_upstream, thread_upstream) { F_AREA(area, f_upstream, thread_upstream); Amagnitude += NV_MAG(area); Flux += F_FLUX(f_upstream, thread_upstream); pressure += F_P(f_upstream, thread_upstream); velocity += F_FLUX(f_upstream, thread_upstream)/(density*NV_MAG(area)); /*pressure = F_P(f_upstream, thread_upstream) + 0.5 * density * pow(velocity,2);*/ /*PA_sum += pressure * NV_MAG(area);*/ /*velocity_sum += velocity;*/ } end_f_loop(f_upstream, thread_upstream) mdot_glob = PRF_GRSUM1(Flux * Amagnitude); A_glob = PRF_GRSUM1(Amagnitude); mflux_glob = mdot_glob / A_glob; PA_glob = PRF_GRSUM1(pressure*Amagnitude); velocity_sum = PRF_GRSUM1(velocity*Amagnitude); velocity_glob = velocity_sum / A_glob; /* AWA velocity */ reynolds = density * fabs(velocity_glob) * Dia / mu; friction_factor = 1/pow((1.8*log10(reynolds)1.5),2); press_avg = (PA_glob/A_glob) + (0.5*density*pow(velocity_glob,2))  (friction_factor*(L/Dia)*(0.5*density*pow(velocity_glob,2)))*frelax; } else { press_avg = 0.0; } begin_f_loop(f, thread) { F_PROFILE(f, thread, position) = press_avg; } end_f_loop(f, thread) Message0("reynolds at %d is %g\n", zone_ID, reynolds); Message0("ff at at %d is %g\n", zone_ID, friction_factor); Message0("velocity %d is %g\n", zone_ID, velocity_glob); Message0("calculated_p at %d is %g\n", zone_ID, press_avg); Message0("\n"); } Outputs looks as follows Code:
reynolds at 51 is 7.5611e+07 ff at at 51 is nan velocity 51 is 290.34 calculated_p at 51 is nan reynolds at 50 is 2.11107e+08 ff at at 50 is nan velocity 50 is 810.635 calculated_p at 50 is nan reynolds at 49 is 1.34733e+08 ff at at 49 is nan velocity 49 is 517.363 calculated_p at 49 is nan reynolds at 48 is 8.88475e+07 ff at at 48 is nan velocity 48 is 341.167 calculated_p at 48 is nan 

May 17, 2021, 20:43 

#2 
Senior Member
Alexander
Join Date: Apr 2013
Posts: 1,723
Rep Power: 28 
to get ff you are using nothing else but reynolds number, which is defined.
so the problem here could be in log10() try add #include <math.h> in header also you are using Amagnitude Flux pressure velocity which are not defined before summation, fix it. they are not zeros by default. compile your code
__________________
best regards ****************************** press LIKE if this message was helpful 

May 18, 2021, 05:44 

#3 
New Member
Join Date: Jan 2018
Posts: 29
Rep Power: 5 
Thank you for your response AlexenderZ,
Yes, I do have all that you have mentioned, was just showing the DEFINE_Profile bit Code:
#include "udf.h" #include <math.h> #define in1 36 #define in2 35 #define in3 34 #define in4 33 #define out1 51 #define out2 50 #define out3 49 #define out4 48 #define frelax 0.0001 #define L 0.672 #define Dia 0.006 #define mu 2.18e5 #define PI 3.14159 #define density 0.9462 /* Step 1: integrate the mass flux at every outlet and set value on corresponding BC Step 2: integrate the static pressure at every inlet and set value on corresponding BC */ DEFINE_PROFILE(massflux_inlet, thread, position) { /* Integrates mass flux through upstream outflow interface Sets the mass flux profile on the downstream interface */ int ID; /* the ID of zone with upstream boundary */ face_t f; real area[ND_ND]; real Amagnitude, Flux, mflux_glob; real mdot_sum, mdot_glob; real A_sum, A_glob; face_t f_upstream; Domain* domain_upstream; /* domain is declared as a variable */ domain_upstream = Get_Domain(1); /* returns fluid domain pointer */ /* Find the mass flux through the upstream interface outlet */ /* set the ID of the upstream outlet interface surface */ int zone_ID = THREAD_ID(thread); if (zone_ID == in1) { ID = out1; } else if (zone_ID == in2) { ID = out2; } else if (zone_ID == in3) { ID = out3; } else if (zone_ID == in4) { ID = out4; } else { ID = 10000; } if (ID != 10000) { Thread* thread_upstream = Lookup_Thread(domain_upstream, ID); Flux = 0.0; Amagnitude = 0.0; begin_f_loop(f_upstream, thread_upstream) if PRINCIPAL_FACE_P(f_upstream, thread_upstream) { F_AREA(area, f_upstream, thread_upstream); Amagnitude += NV_MAG(area); Flux += F_FLUX(f_upstream, thread_upstream); } end_f_loop(f_upstream, thread_upstream) mdot_glob = PRF_GRSUM1(Flux*Amagnitude); A_glob = PRF_GRSUM1(Amagnitude); mflux_glob = mdot_glob / A_glob; } else { mflux_glob = 0.0; } begin_f_loop(f, thread) { F_PROFILE(f, thread, position) = mflux_glob; } end_f_loop(f, thread) Message0("Amagnitude is %f\n", Amagnitude); Message0("massflowrate at %d is %f\n", zone_ID, mflux_glob); } DEFINE_PROFILE(pressure_outlet, thread, position) { /* Averages static pressure on downstream interface (in:*) inflow BC and applies it at the upstream outflow interface (out:*) */ int ID; /* the ID of zone with downstream boundary */ face_t f; real area[ND_ND]; real Amagnitude, pressure, press_avg; real PA_sum, PA_glob; real A_sum, A_glob, velocity_glob, density_sum, density_glob ; real mdot_sum, mdot_glob, Flux, mflux_glob , reynolds, velocity, velocity_sum, friction_factor; face_t f_downstream; Domain* domain_downstream; domain_downstream = Get_Domain(1); /* set the ID of the downstream interface */ int zone_ID = THREAD_ID(thread); /* applied to */ if (zone_ID == out1) { ID = in1; } else if (zone_ID == out2) { ID = in2; } else if (zone_ID == out3) { ID = in3; } else if (zone_ID == out4) { ID = in4; } else { ID = 10000; } if (ID!=10000) { Thread* thread_downstream = Lookup_Thread(domain_downstream, ID); Amagnitude = 0.0; Flux = 0.0; pressure = 0.0; begin_f_loop(f_downstream, thread_downstream) if PRINCIPAL_FACE_P(f_downstream, thread_downstream) { F_AREA(area, f_downstream, thread_downstream); Amagnitude += NV_MAG(area); Flux += F_FLUX(f_downstream, thread_downstream); pressure += F_P(f_downstream, thread_downstream); } end_f_loop(f_downstream, thread_downstream) mdot_glob = PRF_GRSUM1(Flux*Amagnitude); A_glob = PI * 0.25 * pow(Dia,2); mflux_glob = mdot_glob / A_glob; PA_glob = PRF_GRSUM1(pressure*Amagnitude); velocity_glob = mflux_glob / (density*A_glob); reynolds = density * fabs(velocity_glob) * Dia / mu; friction_factor = 1.0/pow((1.8*log10(reynolds)1.5),2); press_avg = (PA_glob/A_glob) + (0.5*density*pow(velocity_glob,2))  (friction_factor*(L/Dia)*(0.5*density*pow(velocity_glob,2))); } else { press_avg = 0.0; } begin_f_loop(f, thread) { /*velocity = mflux_glob/ (density*(PI*0.25*pow(Dia,2)));*/ F_PROFILE(f, thread, position) = press_avg; } end_f_loop(f, thread) Message0("\n"); Message0("Amagnitude is: %0.12f\n", Amagnitude); Message0("A_glob is %f\n", A_glob); Message0("PA_glob is: %f\n", PA_glob); Message0("massflowrate at in%d: is %f\n", zone_ID, mflux_glob); Message0("reynolds at %d is %f\n", zone_ID, reynolds); Message0("ff at at %d is %f\n", zone_ID, friction_factor); Message0("velocity %d is %f\n", zone_ID, velocity_glob); Message0("calculated_p at %d is %f\n", zone_ID, press_avg); Message0("\n"); } the output for which is something like. Here the main concern is that despite following the manual example for obtaining total area, I still get a zero, and perhaps due to it, the other values suffers. Would you have a suggesstion to correct Amagnitude? Idon't see why this should yeild 0. Code:
Amagnitude is 0.000000000000 massflowrate at 36 is 0.000223 Amagnitude is 0.000000000000 massflowrate at 35 is 0.000283 Amagnitude is 0.000000000000 massflowrate at 34 is 0.000173 Amagnitude is 0.000000000000 massflowrate at 33 is 0.000176 Amagnitude is: 0 A_glob is 0.000028 PA_glob is: 0.538942 massflowrate at in51: is 0.000707 reynolds at 51 is 6883.657060 ff at at 51 is 0.034191 velocity 51 is 26.432700 calculated_p at 51 is 18125.940007 Amagnitude is: 0 A_glob is 0.000028 PA_glob is: 0.738375 massflowrate at in50: is 0.001392 reynolds at 50 is 13549.274851 ff at at 50 is 0.028366 velocity 50 is 52.028146 calculated_p at 50 is 23326.701371 Amagnitude is: 0 A_glob is 0.000028 PA_glob is: 0.009078 massflowrate at in49: is 0.001392 reynolds at 49 is 13549.275805 ff at at 49 is 0.028366 velocity 49 is 52.028150 calculated_p at 49 is 2466.904770 Amagnitude is: 0 A_glob is 0.000028 PA_glob is: 0.401181 massflowrate at in48: is 0.000925 reynolds at 48 is 9008.356102 ff at at 48 is 0.031680 velocity 48 is 34.591377 calculated_p at 48 is 15631.361372 

May 19, 2021, 20:36 

#4 
Senior Member
Alexander
Join Date: Apr 2013
Posts: 1,723
Rep Power: 28 
you are not making global summation for area
you've tried, but didn't use it: was Code:
A_glob = PRF_GRSUM1(Amagnitude); Code:
Amagnitude = PRF_GRSUM1(Amagnitude);
__________________
best regards ****************************** press LIKE if this message was helpful 

May 24, 2021, 11:13 

#5 
New Member
Join Date: Jan 2018
Posts: 29
Rep Power: 5 
Tried this and thank you for your suggestion.
But no joy! 

May 24, 2021, 21:47 

#6 
Senior Member
Alexander
Join Date: Apr 2013
Posts: 1,723
Rep Power: 28 
show here fixed code and log
__________________
best regards ****************************** press LIKE if this message was helpful 

Tags 
udf 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Error in Two phase (condensation) modeling  adilsyyed  CFX  15  June 24, 2015 19:42 
Sample UDF for Characterizing the TubeSide 2Phase Fluid  thodij  Fluent UDF and Scheme Programming  0  June 22, 2015 15:09 
Changing the properties of a fluid  TheForumLord  FLUENT  2  June 12, 2014 01:25 
energy eqn + constant fluid properties  cfx_user  Main CFD Forum  1  March 6, 2013 07:13 
change particle properties using udf?  ljp  FLUENT  0  April 1, 2010 15:12 