/*UDF that nucleation of deopet*/
#include "udf.h"
#define PI 3.14159 /* 定义π*/
#define K 1.38e-23 /*定义玻尔兹曼常数 */
#define M 2.99e-26 /*定义水分子质量*/
#define R 461.5 /*定义水蒸气气体常数*/
#define f_v 0.4 /*相对湿度*/
#define q_c 1.0 /*凝结系数*/
#define h_f 2417.8 /*凝结潜热*/
/* Virial State Function Coefficient */
#define v_B_a 10000.0
#define v_B_a1 0.0015
#define v_B_a2 -0.000942
#define v_B_a3 -0.0004882
#define n_e 2.781828
/* Virial State Function Coefficient */
#define cp_a1 46.0
#define cp_a2 1.47276
#define cp_a3 8.3893e-4
#define cp_a4 -2.19989e-7
#define cp_a5 2.46619e-10
#define cp_a6 -9.70466e-14
real S1;G;J;pv1;ps1;T;P;sur_tension;density_c;density_v;thet;r_star;G_star;gama;/*不饱和度，吉布斯焓，成核率，水蒸气分压力，水蒸气饱和压力，温度，压力*/
/* Coefficient B in the Virial State Function*/
double Virial_B(double T)
{
double v_B = 0.;
double v_B_tao = 0.;
v_B_tao = 1500 / T;
v_B = v_B_a1/(1+T/v_B_a)+v_B_a2*pow(n_e, v_B_tao)*pow((1-pow(n_e, -v_B_tao)), 2.5)*pow(v_B_tao, 0.5)+v_B_a3*v_B_tao;
return v_B;
}
/* B1=dBdt in the definition of specific heat, enthalpy and entropy*/
double Virial_dBdt(double T)
{
double V_dBdt = 0.;
double v_B_tao = 0.;
v_B_tao = 1500 / T;
V_dBdt = 0.7323/pow(T,2)+0.0182418*pow(n_e,v_B_tao)*(T+3000)*pow((1-pow(n_e,-v_B_tao)),2.5)
*pow((1/T),2.5)+136.813*pow((1-pow(n_e,-v_B_tao)),1.5)*pow((1/T),2.5)-15./pow((T+10000),2);
return V_dBdt;
}
/* B2=d2B/dt2 in the definition of specific heat, enthalpy and entropy*/
double Virial_ddBddt(double T)
{
double V_ddBddt = 0.;
double v_B_tao = 0.;
v_B_tao = 1500 / T;
V_ddBddt = -1.4646/pow(T,3)-202220.*pow((1-pow(n_e,-v_B_tao)),1.5)*pow((1/T),4.5)-307830.*pow(n_e,-v_B_tao)
*pow((1-pow(n_e,-v_B_tao)),0.5)*pow((1/T),4.5)-410.441*pow((1-pow(n_e,-v_B_tao)),1.5)*pow((1/T),3.5)
+pow(n_e,v_B_tao)*(-82088.1*pow((1/T),4.5)-164.176*pow((1/T),3.5)-0.0273627*pow((1/T),2.5))
*pow((1-pow(n_e,-v_B_tao)),2.5)+30/pow((T+10000),3);
return V_ddBddt;
}
double cp=0.;
double cp0=0.;
double ag=11.16;
v_B = Virial_B(T);
dv_Bdt = Virial_dBdt(T); /* dB/dt */
dv_Bdt2 = Virial_ddBddt(T); /* ddB/ddt */
c_p0 = cp_a1/T+cp_a2+cp_a3*T+cp_a4*pow(T,2)+cp_a5*pow(T,3)+cp_a6*pow(T,4);
c_p = cp0+R*((1-ag*T)*(v_B-dv_Bdt)-dv_Bdt2)*density_v;
c_v= cp0-R*((1+(2*dv_Bdt+dv_Bdt2)))*density_v;
/*计算成核率，液滴生长速度，临界半径*/
DEFINE_ADJUST(chenghelv, d)
{
Thread *t;
cell_t c;
thread_loop_c (t,d)
{
begin_c_loop_all(c,t)
{
T = C_T(c,t);
P = C_P(c,t);
ps1 = -2.88934 + 0.37966 * T - 0.00925 * pow(T, 2) + 1.54e-4 * pow(T, 3);/*饱和压力拟合*/
density_v = -0.00134 + 0.00136 * T - 3.27e-5 * pow(T, 2) + 7.76e-7 * pow(T, 3);/*蒸气密度拟合*/
density_c = 1000.07904 + 0.01231 * T - 0.00586 * pow(T,2) + 1.6e-5 * pow(T,3);/*水密度拟合*/
sur_tension = 1.07*pow((1-T/647.3),1.258)*0.2538*(1-(1-T/647.3));/*表面张力拟合*/
thet=2*h_f /R*T*(h_f/R*T-0.5);/*修正因子*/
S1 = P / ps1;
pv1 = P * f_v;
r_star=2* sur_tension /density_c*R*T* log(S1) ;/*临界半径*/
G_star=4/3* PI*(pow(r_star,2))* sur_tension;/*临界自由能*/
G = 16 / 3 * PI * pow((M / (density_c * log(S1) * K * T)),2) * pow(sur_tension, 3); /*计算吉布斯自由焓*/
C_UDMI(c,t,0) =q_c /(1+thet) *sqrt(2 / PI * sur_tension) * pow(M ,-3/2) * pow(density_v, 2) / density_c * exp( -G / K * T);/*成核率*/
end_c_loop_all(c,t)
}
}
}
评论4