CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   FLUENT (http://www.cfd-online.com/Forums/fluent/)
-   -   Scalar problem! (http://www.cfd-online.com/Forums/fluent/75808-scalar-problem.html)

 lig May 5, 2010 15:16

Scalar problem!

Hi, All,

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.

#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)
{
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*/
F_AREA(A,f,t);

/*If face line at domain boundary, use face values;*/
/*If face lines IN the domain, use equations*/
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*/
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 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){

}
else {
}

/*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;
}

 All times are GMT -4. The time now is 03:01.