UDFs for modelling heat transfer in a porous medium

sutherlandshire
sutherlandshire Member Posts: 1 **
edited November 2023 in Fluids

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

Answers

Categories