[lmz]

/*---------------------------------- UDF 头文件及宏定义 ----------------------------------*/
#include "udf.h"

#define UDM_FOULING_RESISTANCE 0 // 污垢热阻 (m²·K/W)
#define UDM_DEPOSITION_RATE 1 // 沉积速率 (kg/m²·s)
#define UDM_FOULING_THICKNESS 2 // 污垢层厚度 (m)

/*---------------------------------- 全局参数定义 ----------------------------------------*/
real C_inf = 0.4; // 近壁颗粒浓度 (kg/m3)
real S_p = 0.8; // 粘附概率 ((无量纲)
real zeta = 1e5; // 污垢粘附强度因子 (N·s/m²)
real rho_f = 2500.0; // 污垢层密度 (kg/m³)
real lambda_f = 1.2; // 污垢导热系数 (W/(m·K))
real dp = 1e-6; // 粒子直径 (m)，需根据案例调整
real kB = 1.380649e-23; // 玻尔兹曼常数 (J/K)
real Ct = 0.55; // 热泳系数(无量纲)
real Sc_t = 0.7; // 湍流施密特数(无量纲)

/*---------------------------------- 计算布朗扩散系数 DB --------------------------------*/
real calculate_DB(cell_t c, Thread* thread) {
// Stokes-Einstein公式: DB = (kB * T) / (6 * π * eta * r)
real r = dp / 2.0; // 粒子半径
real eta = C_MU(c, thread); // 动力粘度 (Pa·s)
real T = C_T(c, thread); // 从单元获取温度 (单位: K)
real DB = (kB * T) / (6.0 * 3.1415926 * eta * r);
return DB; // 单位: m²/s
}

/*---------------------------------- 计算热泳扩散系数 DT --------------------------------*/
real calculate_DT(cell_t c, Thread* thread) {
// 热泳模型: DT = Ct * eta * T * grad_T / (rho_l * lambda_l)
real eta = C_MU(c, thread); // 动力粘度 (Pa·s)
real T = C_T(c, thread); // 从单元获取温度 (单位: K)
real rho_l = C_R(c, thread); // 流体密度 (kg/m³)
real lambda_l = C_K(c, thread); // 流体导热系数 (W/(m·K))

// 获取温度梯度（单位：K/m）
real grad_T[ND_ND], grad_T_mag;
C_T_G(c,thread,grad_T); // 获取温度梯度分量
grad_T_mag = NV_MAG(grad_T); // 温度梯度模长

// 计算DT
real DT = Ct * eta * T * grad_T_mag / (rho_l * lambda_l);
return DT; // 单位: m²/s
}

/*---------------------------------- 计算湍流扩散系数 ε_p ------------------------------*/
real calculate_epsilon_p(cell_t c,Thread* thread) {
// 从k-epsilon模型获取k和ε值
real k = C_K(c, thread); // 湍流动能 (m²/s²)
real epsilon = C_EPSILON(c, thread); // 湍流耗散率 (m²/s³)

// 计算湍流粘度 μ_t = Cμ * ρ * k² / ε
real Cmu = 0.09; // k-epsilon模型常数
real rho = C_R(c, thread); // 流体密度 (kg/m³)
real mu_t = Cmu * rho * k * k / (epsilon + 1e-20); // 避免除以0

// 计算湍流运动粘度 ν_t = μ_t / ρ
real nu_t = mu_t / (rho + 1e-6);

// 湍流扩散系数 ε_p = ν_t / Sc_t
real epsilon_p = nu_t / Sc_t;
return epsilon_p; // 单位: m²/s
}

//添加重力速度计算函数
real u_gravity(cell_t c, Thread* thread) {
real g[ND_ND] = { 0.0, 0.0, 9.81 }; // 示例：Y方向重力

//获取关联单元及流体属性
real eta = C_MU(c, thread); // 动力粘度 (Pa·s)
real rho_l = C_R(c, thread); // 流体密度 (kg/m³)

// 计算密度差
real delta_rho = rho_f - rho_l; // 颗粒与流体密度差

// Stokes沉降速度公式
real u = (dp * dp * delta_rho * g[2]) / (18.0 * eta + 1e-20); // 单位: m/s
return u;
}

/*---------------------------------- 改进的无因次浓度模型 ----------------------------*/
real calculate_C_plus(real y_plus,cell_t c, Thread* thread) {
real nu = C_MU(c, thread) / (C_R(c, thread)+1e-6);
real Sc = nu / calculate_DB(c, thread);

// 近壁区(y+ < 5)线性分布
if (y_plus < 5.0) {
return 1.0 - y_plus / (5.0 * Sc);
}
// 缓冲层(5 < y+ < 30)
else if (y_plus < 30.0) {
real C_log = log(1.0 + Sc * y_plus) / (2.5 * Sc);
return 1.0 - 0.2 - C_log;
}
// 对数区(y+ >= 30)
else {
return (1.0 / Sc_t) * log(y_plus) + 5.0;
}
}

real calculate_dC_plus_dy_plus(real y_plus, cell_t c, Thread* thread) {
real nu = C_MU(c, thread) / (C_R(c, thread) + 1e-6);
real Sc = nu / calculate_DB(c, thread);

// 近壁区导数
if (y_plus < 5.0) {
return -1.0 / (5.0 * Sc);
}
// 缓冲层导数
else if (y_plus < 30.0) {
return -Sc / (2.5 * (1.0 + Sc * y_plus));
}
// 对数区导数
else {
return 1.0 / (Sc_t * y_plus);
}
}

/*---------------------------------- 自定义函数：计算沉积通量 ----------------------------*/
real calculate_J(cell_t c, Thread* thread,face_t f) {
// 获取流体属性
//real T = F_T(f, thread); // 面温度 (K)
real eta1 = C_MU(c, thread); // 动力粘度 (Pa·s)
real rho = C_R(c, thread); // 密度 (kg/m³)
real lambda_l = C_K(c, thread); // 流体导热系数 (W/(m·K))

// 获取温度梯度
real grad_T[ND_ND], grad_T_mag;
C_T_G(c, thread, grad_T); // 获取温度梯度矢量
grad_T_mag = NV_MAG(grad_T); // 温度梯度模长

//计算扩散系数
real DB = calculate_DB(c,thread); // 使用参数T_face和eta
real DT = calculate_DT(c,thread);
real epsilon_p = calculate_epsilon_p(c, thread); // 修正：传递thread指针

// 获取壁面剪切应力
real tau_w_vector[ND_ND]; // ND_ND为维度（2D为2，3D为3）
C_WALL_SHEAR(f, thread, tau_w_vector);// 调用宏获取剪切应力矢量（单位：Pa）
real tau_w_mag = NV_MAG(tau_w_vector); // 计算标量剪切应力（模长）sqrt(tau_x^2 + tau_y^2 [+ tau_z^2])

// 计算摩擦速度m/s
real u_star = sqrt(tau_w_mag / (rho+1e-20));

// 获取流体运动粘度ν
real nu = eta1 / rho;
real y_plus = F_YPLUS(f, thread);

// 计算无因次参数C+及其梯度（假设线性分布）
real C_plus = calculate_C_plus(y_plus, c, thread); // C+ = C/C∞
real dC_plus_dy_plus = calculate_dC_plus_dy_plus(y_plus, c, thread); // 动态梯度模型

// 计算无因次沉积速率u_d+
real term1 = -(DB / nu + epsilon_p / nu) * dC_plus_dy_plus;
real term2 = -(u_gravity(c, thread) / u_star) * C_plus; // 重力沉降项（需定义u_gravity函数）
real term3 = -(DT / (nu * u_star)) * C_plus; // 热泳项（简化模型）
real term4 = -(epsilon_p / u_star) * C_plus; // 湍流扩散项（经验模型）

real u_d_plus = term1 + term2 + term3 + term4;
u_d_plus = max(u_d_plus, 0.0); // 确保非负

// 计算沉积速度u_d
real u_d = u_d_plus * u_star;

// 计算沉积通量J
real J = S_p * u_d * C_inf;
return J; // 单位: kg/(m²·s)
}

/*---------------------------------- 主UDF：更新局部污垢参数 -----------------------------*/
DEFINE_ADJUST(update_local_fouling, domain)
{
//real g[ND_ND] = { 0.0, 0.0, 9.81 }; // 重力加速度 (m/s²)
Thread *thread;
face_t f;

// 获取时间步长（兼容瞬态/稳态）
real dt = NNULLP(RP_Get_Real("physical-time-step")) ? 1.0 : RP_Get_Real("physical-time-step");

// 遍历所有壁面边界
thread_loop_f(thread, domain)
{
if (THREAD_TYPE(thread) == BOUNDARY_THREAD &&
Get_Boundary_Type(thread) == BT_WALL)
{
begin_f_loop(f, thread)
{
cell_t c = F_C0(f, thread);
Thread* cthread = F_C0_THREAD(f, thread);

// 验证单元有效性
if (FLUID_THREAD_P(cthread))
{
real eta = C_MU(c, cthread); // 动力粘度
real rho_l = C_R(c, cthread); // 密度
real lambda_l = C_K(c, cthread); // 导热系数

// 获取剪切应力模长
real tau_w_vector[ND_ND];
C_WALL_SHEAR(f, thread, tau_w_vector);
real tau_w = NV_MAG(tau_w_vector);

// 计算沉积速率和去除速率
real J = calculate_J(c, thread, f);
real m_d = J; // 沉积速率 (kg/(m²·s))
real m_r = tau_w * F_UDMI(f, thread, UDM_FOULING_THICKNESS) / zeta; // 公式(5)

// 更新污垢层厚度（时间步长dt=1s）
real x_f_new = F_UDMI(f, thread, UDM_FOULING_THICKNESS) + (m_d - m_r) * dt;
x_f_new = (x_f_new < 0.0) ? 0.0 : x_f_new; // 确保非负
F_UDMI(f, thread, UDM_FOULING_THICKNESS) = x_f_new;

// 计算局部污垢阻力
real R_f = x_f_new / lambda_f; // 公式(6)
F_UDMI(f, thread, UDM_FOULING_RESISTANCE) = R_f;
F_UDMI(f, thread, UDM_DEPOSITION_RATE) = m_d;
}
}
end_f_loop(f, thread);
}
}
}

/*---------------------------------- UDF初始化 ------------------------------------------*/
DEFINE_ON_DEMAND(init_udm)
{
// 初始化UDM存储为0
Domain* domain = Get_Domain(1);
Thread* thread;
face_t f;

thread_loop_f(thread, domain)
{
// 仅处理壁面边界线程
if (THREAD_TYPE(thread) == BOUNDARY_THREAD &&
Get_Boundary_Type(thread) == BT_WALL)
{
begin_f_loop(f, thread)
{
F_UDMI(f, thread, UDM_FOULING_RESISTANCE) = 0.0;
F_UDMI(f, thread, UDM_DEPOSITION_RATE) = 0.0;
F_UDMI(f, thread, UDM_FOULING_THICKNESS) = 0.0;
}
end_f_loop(f, thread);
}
}
}

The above code debugging is not possible, the main problem is define_adjust and define_on_demand macro definition inside the boundary condition setting, please help the expert teachers to modify and improve, thank you.