[NEEPU]

#include "udf.h"
#define T_w 268.15
#define T0 273.15
#define L_f 3.349*pow(10,5)
#define V_w 5
#define L_e 2.5*pow(10,6)
#define c_a 1009
#define T_s 273.15
#define C_i 2261
#define C_w 4187
#define X1 6.75*pow(10,-2)
#define z 60
#define u_a 14.8*pow(10,-6)
#define p_a 1293
#define u_c 17.9*pow(10,-6)
#define R_v 461.1
#define k 0.0237
#define r 0.0015
#define m_r 15

FILE *fp1;
FILE *fp2;
#define LWC 0.6
enum
{
CEN_X,
CEN_Y,
f_pres,
NUM_OF_HITS,
num,
beta,
m_in,
m_out,
m_va,
m_so,
m_im
};
int i;
face_t f1;
real C2[ND_ND];
int num_in_date;
int num_col=52;
DEFINE_ON_DEMAND(face_num_pres_max)
{
Domain *domain = Get_Domain(1);
int ID = 3;
int j = 0;
real x[ND_ND];
face_t f;
int c ;
real max;
real a[300];
FILE *fp0;
Thread *t = Lookup_Thread(domain, ID);
fp0 = fopen("position.txt", "a");
begin_f_loop(f, t)
{
F_CENTROID(x, f, t);
F_UDMI(f, t, CEN_X) = x[0];
F_UDMI(f, t, CEN_Y) = x[1];
F_UDMI(f, t, f_pres) = F_P(f, t);
a[i] = F_UDMI(f, t, f_pres);
i = i + 1;
F_UDMI(f, t, num) = i;
printf("x=%f y=%f i=%d n=%f\n", F_UDMI(f, t, CEN_X), F_UDMI(f, t, CEN_Y), i, F_UDMI(f, t, num));
fprintf(fp0, "x = %f y = %f i = %d n=%d\n", F_UDMI(f, t, CEN_X), F_UDMI(f, t, CEN_Y), i,F_UDMI(f,t,num));
}
end_f_loop(f, t)
for (max = a[0], c = 0; c < i; c++)
{
if (max < a[c])
max = a[c];
}
begin_f_loop(f, t)
{
if (F_UDMI(f, t, f_pres) == max)
{
F_CENTROID(C2, f, t);
f1 = f;
F_CENTROID(x, f1, t);
printf("x[0]=%f x[1]=%f", x[0], x[1]);
fprintf(fp0, "max x[0]=%f x[1]=%f\n", x[0], x[1]);


}

}
end_f_loop(f, t)
fclose(fp0);

}
DEFINE_DPM_EROSION(SJ, p, t, f, normal, alpha, Vmag, mdot)
{
real x[ND_ND];
real A[ND_ND], area;
real A1, B1, B2;
F_CENTROID(x, f, t);
num_in_date = num_in_date + 1;
F_UDMI(f, t, NUM_OF_HITS) = num_in_date;
F_AREA(A, f, t);
area = NV_MAG(A);
A1 = area;
B1 = F_UDMI(f, t, NUM_OF_HITS);
B2 = (B1*M_PI * 1000 * pow(P_DIAM(p) / 2, 3)) * 4 / (0.0006*area * 5 * 3 * CURRENT_TIME);
F_UDMI(f, t, beta) = B2;
printf("B1=%f B2=%f A1=%f x[0]=%f x[1]=%f d=%f i=%d\n", B1, B2, A1, x[0], x[1],P_DIAM(p),F_UDMI(f,t,num));
fp1 = fopen("position_collision.txt", "a");
fprintf(fp1, "B1=%f B2=%f A1=%f x[0]=%f x[1]=%f d=%f i=%d\n", B1, B2, A1, x[0], x[1], P_DIAM(p), F_UDMI(f, t, num));
fclose(fp1);
}
DEFINE_EXECUTE_AT_END(icing)
{
int i2;
R_eq = r * (1 + m_r / 2);
int i1;
int j;
real m_out1;
real m_va1;
real m_so1;
real m_in1;
real A[ND_ND];
real area;
real Pr, Re, Gr, Nu_f, Nu_p, h_c, r_c, P_ss, P_sw,A_p,n,f_so,m_im1;
Domain *domain = Get_Domain(1);
Thread *t=Lookup_Thread(domain,3);
face_t f;
i2 = F_UDMI(f1, t, num);
i1 = F_UDMI(f1, t, num);
F_UDMI(f1, t, m_in) = 0;
Pr = (u_a*c_a) / k;
Re = (2 * R_eq*p_a*V_w) / u_a;
if (Re >= 40 && Re <= 4000)
{
A_p = 0.6830;
n = 0.4660;
}
else if (Re >= 4000 && Re <= 40000)
{
A_p = 0.1930;
n = 0.6180;
}
else if (Re >= 40000 && Re <= 400000)
{
A_p = 0.0266;
n = 0.8050;
}
Gr = (9.81*(T_s - T_w)*pow(2 * R_eq, 3)) / (pow(V_w, 2)*(273.15 + (T_s + T_w) / 2));
Nu_f = 0.48*pow(Gr*Pr, 3);
Nu_p = A_p * pow(Re, n)*pow(Pr, 1 / 3);
h_c = ((Nu_f + Nu_p)*k) / (2 * R_eq);
r_c = 1 - (pow(u_c / V_w, 2))*(1 - pow(Pr, 1 / 3));
P_ss = 611.01109 - 44.48156*5 + 1.4188*pow(5, 2) -0.02389*pow(5, 3);
P_sw = 611.01109;
F_CENTROID(A, f, t);
area = NV_MAG(A);
m_va1 = (h_c*(P_ss / T_s - P_sw / T_w))*CURRENT_TIME / (461.1*p_a*c_a);
m_im1= F_UDMI(f1, t, beta)*LWC*V_w*area*CURRENT_TIME;
F_UDMI(f1, t, m_im) = m_im1;
F_UDMI(f1, t, m_va) = m_va1;
m_in1 = F_UDMI(f1, t, m_in);
f_so = 1 - (m_im1*((T_w - T0) + L_f + pow(V_w, 2) / 2) + m_in1 * (L_f + pow(V_w, 2)) - m_va1 * (2 * L_f + L_e + pow(V_w, 2) / 2) + h_c * area*r_c + pow(V_w, 2) / 2 * c_a - (T_s - T_w)*(h_c*area + k / sqrt(M_PI*X1*z))) / ((m_im1 + m_in1)*(L_f + pow(V_w, 2) / 2));
F_UDMI(f1, t, m_out) = ((1 - f_so)*F_UDMI(f1, t, m_im) - m_va1)*CURRENT_TIME;
m_out1 = F_UDMI(f1, t, m_out) / 2;
m_so1 = f_so * F_UDMI(f1, t, m_im);
F_UDMI(f1, t, m_so) = m_so1;


for (j = 0; j < num_col / 2; j++)
{
begin_f_loop(f, t)
{
if (F_UDMI(f, t, num) == i2 + 1)
{
F_CENTROID(A, f, t);
area = NV_MAG(A);
F_UDMI(f, t, m_in) = m_out1;
m_in1 = F_UDMI(f, t, m_in);
m_im1 = F_UDMI(f, t, beta)*LWC*V_w*area*CURRENT_TIME;
F_UDMI(f, t, m_im) = m_im1;
f_so = 1 - (m_im1*((T_w - T0) + L_f + pow(V_w, 2) / 2) + m_in1 * (L_f + pow(V_w, 2)) - m_va1 * (2 * L_f + L_e + pow(V_w, 2) / 2) + h_c * area*r_c + pow(V_w, 2) / 2 * c_a - (T_s - T_w)*(h_c*area + k / sqrt(M_PI*X1*z))) / ((m_im1 + m_in1)*(L_f + pow(V_w, 2) / 2));
F_UDMI(f, t, m_out) = ((1 - f_so)*(F_UDMI(f, t, m_im)+F_UDMI(f,t,m_in)) - m_va1)*CURRENT_TIME;
m_out1 = F_UDMI(f, t, m_out);
F_UDMI(f, t, m_va) = m_va1;
m_so1 = (f_so * (m_im + m_in));
F_UDMI(f, t, m_so) = m_so1;
}
}
end_f_loop(f, t)
i2 = i2 + 1;
}
}
DEFINE_EXECUTE_AT_END(icing1)
{
R_eq = r * (1 + m_r / 2);
int i1;
int j;
real m_out1;
real m_va1;
real m_in1;
real m_so1;
face_t f;
real A[ND_ND];
real area;
real Pr, Re, Gr, Nu_f, Nu_p, h_c, r_c, P_ss, P_sw, A_p, n, f_so,m_im1;
Domain *domain = Get_Domain(1);
Thread *t = Lookup_Thread(domain, 3);

i1 = F_UDMI(f1, t, num);
F_UDMI(f1, t, m_in) = 0;

Pr = (u_a*c_a) / k;
Re = (2 * R_eq*p_a*V_w) / u_a;
if (Re >= 40 && Re <= 4000)
{
A_p = 0.6830;
n = 0.4660;
}
if (Re >= 4000 && Re <= 40000)
{
A_p = 0.1930;
n = 0.6180;
}
if (Re >= 40000 && Re <= 400000)
{
A_p = 0.0266;
n = 0.8050;
}
Gr = (9.81*(T_s - T_w)*pow(2 * R_eq, 3)) / (pow(V_w, 2)*(273.15 + (T_s + T_w) / 2));
Nu_f = 0.48*pow(Gr*Pr, 3);
Nu_p = A_p * pow(Re, n)*pow(Pr, 1 / 3);
h_c = ((Nu_f + Nu_p)*k) / (2 * R_eq);
r_c = 1 - (pow(u_c / V_w, 2))*(1 - pow(Pr, 1 / 3));
P_ss = 611.01109 + 44.48156*(-5) + 1.4188*pow(5, 2) + 0.02389*pow(-5, 3);
P_sw = 611.01109;
F_CENTROID(A, f, t);
area = NV_MAG(A);
m_va1 = (h_c*(P_ss / T_s - P_sw / T_w))*CURRENT_TIME / (461.1*p_a*c_a);
m_im1 = F_UDMI(f1, t, beta)*LWC*V_w*area*CURRENT_TIME;
F_UDMI(f1, t, m_im) = m_im1;
F_UDMI(f1, t, m_va) = m_va1;
m_in1 = F_UDMI(f1, t, m_in);
f_so = 1 - (m_im1*((T_w - T0) + L_f + pow(V_w, 2) / 2) + m_in1 * (L_f + pow(V_w, 2)) - m_va1 * (2 * L_f + L_e + pow(V_w, 2) / 2) + h_c * area*r_c + pow(V_w, 2) / 2 * c_a - (T_s - T_w)*(h_c*area + k / sqrt(M_PI*X1*z))) / ((m_im1 + m_in1)*(L_f + pow(V_w, 2) / 2));
F_UDMI(f1, t, m_out) = ((1 - f_so)*F_UDMI(f,t,m_im) - m_va1)*CURRENT_TIME;
m_out1 = F_UDMI(f1, t, m_out) / 2;
m_so1 = f_so * F_UDMI(f1, t, m_im);
F_UDMI(f1, t, m_so) = m_so1;
for (j = 0; j < num_col/ 2; j++)
{
begin_f_loop(f, t)
{
if (F_UDMI(f, t, num) == i1 - 1)
{
F_CENTROID(A, f, t);
area = NV_MAG(A);
F_UDMI(f, t, m_in) = m_out;
m_in1 = F_UDMI(f, t, m_in);
m_im1 = F_UDMI(f, t, beta)*LWC*V_w*area*CURRENT_TIME;
F_UDMI(f, t, m_im) = m_im1;
f_so = 1 - (m_im1*((T_w - T0) + L_f + pow(V_w, 2) / 2) + m_in1 * (L_f + pow(V_w, 2)) - m_va1 * (2 * L_f + L_e + pow(V_w, 2) / 2) + h_c * area*r_c + pow(V_w, 2) / 2 * c_a - (T_s - T_w)*(h_c*area + k / sqrt(M_PI*X1*z))) / ((m_im1 + m_in1)*(L_f + pow(V_w, 2) / 2));

F_UDMI(f, t, m_out) = ((1 - f)*F_UDMI(f,t,m_im) - m_va1)*CURRENT_TIME;
m_out1 = F_UDMI(f, t, m_out);
F_UDMI(f, t, m_va) = m_va1;
m_so1 = f_so * (m_im + m_in);
F_UDMI(f, t, m_so) = m_so1;
}
}
end_f_loop(f, t)
i1 = i1 - 1;
}
}
DEFINE_ON_DEMAND(icing_sum)
{
Domain *d;
Thread *t = Lookup_Thread(d, 3);
face_t f;
fp2 = fopen("position_icing_sum.txt", "a");
begin_f_loop(f, t)
{
fprintf(fp2, "m_im=%f m_in=%f m_out=%f m_so=%f m_va=%f x=%f y=%f i=%d\n", F_UDMI(f, t, m_im), F_UDMI(f, t, m_in), F_UDMI(f, t, m_out), F_UDMI(f, t, m_so), F_UDMI(f, t, m_va), F_UDMI(f, t, CEN_X),F_UDMI(f, t, CEN_Y), F_UDMI(f, t, num));
}
end_f_loop(f, t)
fclose(fp2);
}