Hello, I am unsure if this is the place to ask this, however I am really struggling to create UDFs to fit a numerical model which must be applied as UDFs in ANSYS Fluent to model a phase change material in a porous medium. I have attached the mathematical model I am trying to follow as well as my current code. I am specifically stuck on adding the energy source term, however I am also facing difficulty in adding an interfacial area. I believe my source terms in FxFy SxSy are correct however I am not definite on this. Any help would be greatly appreciated.
#include "udf.h"
#define Tsolidus 300.0
#define Tliquidus 304.0
#define L_f 250 //Latent heat of fusion
#define UDM_RHO_F_PREV 0 // UDM index for storing previous rho_f
#define UDM_LAMBDA_PREV 1 // UDM index for storing previous lambda
DEFINE_EXECUTE_ON_LOADING(init_values, udf_name)
{
Thread *t;
cell_t c;
Domain *d = Get_Domain(1); // Get domain using Fluent API, assuming single-phase flow
/* Loop over all cells to initialize UDM values */
thread_loop_c(t, d)
{
begin_c_loop_all(c, t)
{
C_UDMI(c,t,UDM_RHO_F_PREV) = 1000.0; //initializing previous rho_f
C_UDMI(c,t,UDM_LAMBDA_PREV) = 0.0; //initializing previous lambda
}
end_c_loop_all(c, t)
}
}
DEFINE_SOURCE(energy_sourceN, c, t, dS, eqn)
{
real lambda, lambda_prev, V[ND_ND], rho_f, Sl, d_lambda_dt, grad_lambda_dot_V, T, rho_f_prev;
T = C_T(c, t);
rho_f = C_R(c, t);
C_CENTROID(V, c, t);
rho_f_prev = C_UDMI(c, t, UDM_RHO_F_PREV); // Fetch the previous value of rho_f from UDM
lambda_prev = C_UDMI(c, t, UDM_LAMBDA_PREV); // Fetch the previous value of lambda from UDM
// Calculate lambda
if (T < Tsolidus)
lambda = 0.0;
else if (T > Tliquidus)
lambda = 1.0;
else
lambda = (T - Tsolidus) / (Tliquidus - Tsolidus);
// Calculate d(lambda*rho_f*L_f)/dt
d_lambda_dt = ((rho_f * lambda * L_f * 0.85) - (rho_f_prev * lambda_prev * L_f * 0.85)) / CURRENT_TIMESTEP;
// The gradient calculation requires more attention. Below is just a placeholder, and you need to implement the correct gradient calculation.
grad_lambda_dot_V = (C_R_G(c,t)[0] * V[0] + C_R_G(c,t)[1] * V[1]) * lambda * L_f;
Sl = d_lambda_dt + grad_lambda_dot_V;
dS[eqn] = 0; // Adjust the source term linearization if necessary
// Store the current time step values of rho_f and lambda to UDM
C_UDMI(c, t, UDM_RHO_F_PREV) = rho_f;
C_UDMI(c, t, UDM_LAMBDA_PREV) = lambda;
return Sl;
}
DEFINE_SOURCE(F_xmomentum, c, t, dS, eqn)
{
real mu_f, rho_f, K, C, magnitude_V, force;
real velocity[ND_ND];
velocity[0] = C_U(c,t); // u-component
velocity[1] = C_V(c,t); // v-component
/* Fetching effective viscosity and density from Fluent */
mu_f = C_MU_EFF(c, t); // Dynamic Viscosity of fluid
rho_f = C_R(c,t); // Density of fluid
K = 1.452827637854882e-04; // Define or fetch the value of K
C = 0.058424315707417; // Define or fetch the value of C
magnitude_V = sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1]);
force = (mu_f/K + rho_f*C*magnitude_V/sqrt(K)) * (velocity[0]);
dS[eqn] = 0; // source term linearization, adjust if necessary
return force;
}
DEFINE_SOURCE(F_ymomentum, c, t, dS, eqn)
{
real mu_f, rho_f, K, C, magnitude_V, force;
real velocity[ND_ND];
velocity[0] = C_U(c,t); // u-component
velocity[1] = C_V(c,t); // v-component
/* Fetching effective viscosity and density from Fluent */
mu_f = C_MU_EFF(c, t); // Dynamic Viscosity of fluid
rho_f = C_R(c,t); // Density of fluid
K = 1.452827637854882e-04; // Define or fetch the value of K
C = 0.058424315707417; // Define or fetch the value of C
magnitude_V = sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1]);
force = (mu_f/K + rho_f*C*magnitude_V/sqrt(K)) * (velocity[1]);
dS[eqn] = 0; // source term linearization, adjust if necessary
return force;
}
DEFINE_SOURCE(S_xmomentum,c,t,dS,eqn)
{
real velocity[ND_ND], Am, lambda, T, source;
velocity[0] = C_U(c,t); // u-component
velocity[1] = C_V(c,t); // v-component
Am = 10000; // Mushy Zone Constant
T = C_T(c, t);
if (T < Tsolidus)
lambda = 0.0;
else if (T > Tliquidus)
lambda = 1.0;
else
lambda = (T - Tsolidus) / (Tliquidus - Tsolidus);
source = (Am * ((1.0 - lambda)*(1.0 - lambda)) /(lambda*lambda*lambda+0.001))*velocity[0] ;
return source;
}
DEFINE_SOURCE(S_ymomentum,c,t,dS,eqn)
{
real velocity[ND_ND], Am, lambda, T, source;
velocity[0] = C_U(c,t); // u-component
velocity[1] = C_V(c,t); // v-component
Am = 10000; // Mushy Zone Constant
T = C_T(c, t);
if (T < Tsolidus)
lambda = 0.0;
else if (T > Tliquidus)
lambda = 1.0;
else
lambda = (T - Tsolidus) / (Tliquidus - Tsolidus);
source = (Am * ((1.0 - lambda)*(1.0 - lambda)) /(lambda*lambda*lambda+0.001))*velocity[1] ;
return source;
}
DEFINE_SOURCE(h_sf,c,t,dS,eqn)
{
real dl, epsilon, mu, rho, nu, pr, rho_pcm, u, v, mu_pcm, k_pcm, Re_d, h_sf;
dl = 0.0194; //Diameter Matlab
epsilon = 0.85; //Porosity Matlab
mu = C_MU_L(c,t); // Laminar viscosity
rho = C_R(c,t); // Density
nu = mu/rho; // Kinematic Viscosity
pr = nu/rng_alpha(1.,mu,mu); // Prandtl Number
rho_pcm = C_R(c,t);
u = C_U(c,t);
v = C_V(c,t);
mu_pcm = C_MU_L(c,t);
k_pcm = C_K_L(c,t);
Re_d = (rho_pcm * sqrt(u*u + v*v) *dl) / (epsilon * mu_pcm);
if (Re_d > 0 && Re_d <= 40)
h_sf = (0.76 * pow(Re_d, 0.4) * pow(pr, 0.37) * k_pcm)/dl;
else if (Re_d > 40 && Re_d <= 1000)
h_sf = (0.52 * pow(Re_d, 0.5) * pow(pr, 0.37) * k_pcm)/dl;
else if (Re_d > 1000 && Re_d <= 20000)
h_sf = (0.26 * pow(Re_d, 0.6) * pow(pr, 0.37) * k_pcm)/dl;
else
Error("Re_d value out of range!\n");
exit(1);
return h_sf;
}
