[lig]

Hi, All,

I really need your help, nobody to ask around me.

I am using Eularian method to simulate the particle transport through a porous zone. I use UDS function to modify the scalar equation. The scalar is number concentration. The particle has constant density, so the default scalar equation was divided by density. The gravitational effect was added in the convective part using settling velocity and combine with original convective flux. The porous zone will filtrate the particle and treated as a sink. The flux that particle filtrated by the porous zone was added as source part in the scalar equation.

The following are diffusivity, convective flux and source code. They were able to be compiled, however the result is very odd. The concentration among the porous zone are negative and It is no changed before the porous zone when run the program without input the source term in the porous zone.

Please help me to check if there is something wrong with those code!! Really appreciate!!!

#include "udf.h"

#define mu 1.789e-5 /*viscosity (N s/m^2) */
#define lambda 0.066 /*mean free path (um) */
#define pi 3.1415926 /*constant Pi */
#define K 1.38e-23 /*Boltzmann constant (Nm/K)*/
#define T 293 /*Temperature at 20C (K)*/

DEFINE_DIFFUSIVITY(c_diff_scalar, c, t, i)
{
real dm[10];
dm[0]=0.00032249e-7;
dm[1]=0.00017630e-7;
dm[2]=0.000092080e-7;
dm[3]=0.000058463e-7;
dm[4]=0.000039309e-7;
dm[5]=0.000027884e-7;
dm[6]=0.000019417e-7;
dm[7]=0.000016148e-7;

return (dm[i] + C_MU_T(c, t) / (0.7 * C_R(c, t)));
}


#include "udf.h"

#define Dcyl 0.006 /*cylinder diameter of tree stem&twig (m) */
#define Dleaf 0.001 /*leaf diameter of tree stem&twig (m) */
#define K 1.38e-23 /*Boltzmann constant (Nm/K)*/
#define T 293 /*Temperature at 20C (K)*/
#define mu 1.789e-5 /*viscosity (N s/m^2) */
#define rhoP 1050 /*density of particle (Kg/m^3)*/
#define rhoA 1.225 /*density of air (Kg/m^3)*/
#define pi 3.1415926 /*constant Pi */
#define lambda 0.066 /*mean free path (um) */
#define g 9.81 /*gravity acceleration m/s^2*/
#define Cpol 1.32 /*Polhausen coefficient*/
#define dx 0.16 /*dx=W/n, W=1.6m, n=10 (m)*/

DEFINE_UDS_FLUX(uds_flux_domain, f, t, i)
{
Thread *t0, *t1=NULL;
cell_t c0, c1=-1;
real x[ND_ND];
real NV_VEC(psi_vec), NV_VEC(A); /*spidefault=rhoP*velocity */
real dens, flux;
real Usettle[10];
Usettle[0]=0.00002882;
Usettle[1]=0.00007936; /*particle diameter of scalar_0 (um) */
Usettle[2]=0.00025543;
Usettle[3]=0.00059863;
Usettle[4]=0.00128009;
Usettle[5]=0.00249168;
Usettle[6]=0.00505861;
Usettle[7]=0.00726958;


c0 = F_C0(f,t); /*return the index of a face's neighboring C0 cell*/
t0 = THREAD_T0(t); /*return the cell thread of a given face's adjacent cell C0*/
F_AREA(A,f,t);

/*If face line at domain boundary, use face values;*/
/*If face lines IN the domain, use equations*/
if(! BOUNDARY_FACE_THREAD_P(t)) {
t1 = THREAD_T1(t); /*return the cell thread of a given face's adjacent cell C1, if it exists*/
c1 = F_C1(f,t); /*Get cell on other side of face, i.e. return the cell index for C1, if it exists*/
}
else {
t1 = NULL;
c1 = -1;
}

if(NULL == t1) {
F_CENTROID(x,f,t);
if(x[1]<.00001&&x[0]>-66.00000&&x[0]<110.00000) {
NV_D(psi_vec, =, 0, Usettle[i]);
}
else{
NV_D(psi_vec, =, 0, 0,0);
}

/*NULLP and NNULLP are to check whether storage has been allocated for user-defined scalars*/
/*NULLP returns TRUE if storage is not allocated*/
/*NNULLP returns TRUE if storage is allocated*/
if(NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t);
else dens = C_R(c0,t0);
flux = NV_DOT(psi_vec,A); /*flux through Face*/
}
else {
NV_D(psi_vec, =, 0, Usettle[i]);
dens = 0.5 * (C_R(c0,t0)+C_R(c1,t1));
flux = NV_DOT(psi_vec,A);
}
flux = F_FLUX(f,t)/dens-flux;

return flux;
}



#include "udf.h"

#define dp 15e-6 /*particle diameter (um) */
#define Dcyl 0.006 /*cylinder diameter of tree stem&twig (m) */
#define Dleaf 0.001 /*leaf diameter of tree stem&twig (m) */
#define K 1.38e-23 /*Boltzmann constant (Nm/K)*/
#define T 293 /*Temperature at 20C (K)*/
#define mu 1.789e-5 /*viscosity (N s/m^2) */
#define rhoP 1050 /*density of particle (Kg/m^3)*/
#define rhoA 1.225 /*density of air (Kg/m^3)*/
#define pi 3.1415926 /*constant Pi */
#define lambda 0.066 /*mean free path (um) */
#define g 9.81 /*gravity acceleration m/s^2*/
#define Cpol 1.32 /*Polhausen coefficient*/
#define dx 0.16 /*dx=W/n, W=1.6m, n=10 (m)*/
#define dz 0.1375 /*dz=H/n H=2.2m, n=16, m */

/*Calculated constant by giving a dp*/
#define Cc 1.010 /*Slip corrction factor*/
#define dm 1.61483e-12 /*laminar diffusion coefficient (m^2/s) */
#define Trel 0.000733486 /*relaxation time of the particle (s)*/
#define Usettle 0.00726958 /*settling speed of particle (m/s) */
#define B 1.16206e-7
#define D 68458.70124
/************************************************** ************************************************** ***************/

DEFINE_SOURCE(uds_source7,c,t,dS,eqn)
{

real source;
real x[ND_ND], z, Ucell;

/**********************For collection efficiency of leaves and cylinder of the tree (CEcell) ************************/
/*CE is the total collection efficiency =CEcyl+CEleaf (<=1)at each cell */
real CEcell,CEcell_c,CEcyl,CEleaf, SADcell;
real Ec, Ec_c, El, El_c;/*collection efficiency of each cell due to all deposition mechanisms*/
real Edc, Edl; /*collection efficiency of cylinder&leaf due to DIFFUSION*/
real Eiic, Eiil; /*collection efficiency of cylinder&leaf due to INTERCEPTION+IMPACTION*/
real Eic, Eil; /*collection efficiency of cylinder&leaf due to IMPACTION*/
real Rec, Rel,Stkc, Stkl, Rc, Rl;

Ucell=C_U(c,t);
z=x[1];
Rec=D*Dcyl*Ucell; /*Reynolds number for flow around a cylinder*/
Rel=D*Dleaf*Ucell; /*Reynolds number for flow around a leaf*/
Stkc=Trel*Ucell/Dcyl;/*Stk number for flow around a cylinder*/
Stkl=Trel*Ucell/Dleaf; /*Stk number for flow around a leaf*/
Rc=dp/Dcyl; /*interception parameter for flow around a cylinder*/
Rl=dp/Dleaf; /*interception parameter for flow around a leaf*/


begin_c_loop(c,t)
{
/**********************For collection efficiency of leaves and cylinder of the tree (CEcell) ************************/

/*for cylinder collection efficiency at each cell due to all deposition mechanisms*/
Edc= B/pow((Dcyl*Ucell),0.5);
Eiic = (Stkc/(Stkc+0.8)-(2.56-log10(Rec)-3.2*Rc)/(10*sqrt(Stkc)))*(1+Rc);
Eic = Eiic;
if (Eic >(1+Rc)){
Eiic = pow(Stkc/(Stkc+0.8),2);
}
else if (Eic< 0){
Eiic= pow(Stkc/(Stkc+0.8),2);
}
else{
Eiic= Eic;
}
/*cylinder*/Ec = (Edc + Eiic);
Ec_c=Ec;
if (Ec_c>=1){
Ec=1;
}
else if (z>((3/8)*x[0]+2.5)&&z>((-3/8)*x[0]+1.9)&&z<(-2.1*x[0]-1.46)&&z<(2.1*x[0]+1.9)){
Ec=0;
}
else {
Ec=Ec_c;
}

/*for cylinder collection efficiency at each cell due to all deposition mechanisms*/
Edl = B/pow((Dleaf*Ucell),0.5);
Eiil= (Stkl/(Stkl+0.8)-(2.56-log10(Rel)-3.2*Rl)/(10*sqrt(Stkl)))*(1+Rl);
Eil = Eiil;
if (Eiil >(1+Rl)){
Eiil = pow(Stkl/(Stkl+0.8),2);
}
else if (Eil < 0){
Eiil = pow(Stkl/(Stkl+0.8),2);
}
else{
Eiil = Eil;
}
/*leaf*/El = (Edl + Eiil);
El_c=El;
if (El_c>=1){
El=1;
}
else if (z>((3/8)*x[0]+2.5)&&z>((-3/8)*x[0]+1.9)&&z<(-2.1*x[0]-1.46)&&z<(2.1*x[0]+1.9)){
El=0;
}
else {
El=El_c;
}
/*Surface area density (SAD)of foliage */
if (z <= 2.07){

SADcell = -3.99043*pow(z,3)+12.8994*pow(z,2)-7.58257*z+1.89096;
}
else {
SADcell = -224.925*pow(z,2)+912.217*z-918.31;
}

CEcyl = 1-exp(1-SADcell*Ec*dx);
CEleaf= 1-exp(1-SADcell*El*dx);

/*total collection efficiency at each cell */
CEcell=CEcyl+CEleaf;
CEcell_c=CEcell;
if (CEcell_c>=1){
CEcell=1;
}
else {
CEcell=CEcell_c;
}

source = C_UDSI(c,t,7)*CEcell*(Ucell/dx-Usettle/dz);
dS[eqn] = CEcell*(Ucell/dx-Usettle/dz);
}
end_c_loop(c,t)


return source;
}