Hi Guys

I've noticed a lot of people are asking how to add source terms of diffusion of UDSs to either UDS eqns or maybe other eqns.

The source is d/dx( T d(phi)/dx)) where phi is a UDS scalar and T = T(x,y,z) is a vector for the directional diffusivity.

Find below my code for calculating this. Its not complete (some helper functions are missing) and I'm not going to explain it in detail. But it should give you guys some idea.

The basic principle is to go back to the finite volume conservation of fluxes for the entire volume. Fluxes thus depend on the gradient across faces times face properties - which are calculated as the average of the neighbour cells. If you don't understand please read the book by Patankar - Nuemrical heat Transfer 1980.

So here it is:-

/* ---------------------------------------------------

Directional Diffusivity Source Term

--------------------------------------------------- */ real Directional_Diffusivity_Srce(cell_t cell,Thread *thread,real dS[],int eqn,

int uds, real Uds_Diff_Coeff[], int p) {

real source=0.0;

real A[ND_ND];

real x0[ND_ND], x1[ND_ND];

real dx[ND_ND], dxp[ND_ND], Diff_Coeff[ND_ND];

face_t face,fparent;

cell_t c0,c1;

Thread *t0,*t1;

Thread *f_thread, *tparent;

real thi, dens, mag_dx;

int i,numbf,nf,maxf;

/* ---- source is for directional diffusivity */

dS[eqn] = 0.0;

/* --- zero source at first iteration ! */

if (CURRENT_ITER > 0)

{

begin_loop_dimensions(i)

Diff_Coeff[i] = Uds_Diff_Coeff[i] - C_UDSI_DIFF(cell,thread,uds);

/* ---- loop over all faces - adding diffusion source */

source = 0.0;

dS[eqn] = 0.0;

/* ---- loop over faces !*/

c_face_loop(cell,thread,numbf)

{

face = C_FACE(cell,thread,numbf);

f_thread = C_FACE_THREAD(cell,thread,numbf);

fparent = face;

tparent = f_thread;

/* ---- any child faces ? */

nf=0;

maxf = (PARENT_FACE_THREAD_P(tparent) ? F_NKIDS(fparent,tparent) : 1);

/* ---- loop over all faces */

while (nf < maxf)

{

/* ---- parent thread ? */

if (PARENT_FACE_THREAD_P(tparent))

{

face = F_CHILD(fparent,tparent,nf);

f_thread = F_CHILD_THREAD(fparent,tparent,nf);

}

nf++;

/* ---- face fluxes are what??? */

F_AREA(A,face,f_thread);

/* ---- interior face thread! */

if (THREAD_TYPE(f_thread) == THREAD_F_INTERIOR)

{

c0 = F_C0(face,f_thread);

if (c0 != cell)

{

c0 = F_C1(face,f_thread);

t0 = F_C1_THREAD(face,f_thread);

c1 = F_C0(face,f_thread);

t1 = F_C0_THREAD(face,f_thread);

Neg_Vector(A,A);

}

else

{

c1 = F_C1(face,f_thread);

t1 = F_C1_THREAD(face,f_thread);

t0 = F_C0_THREAD(face,f_thread);

}

thi = Dot_Product(A,Diff_Coeff);

dens = FACE_PROPERTY(CELL_SEC_PHASE_APP_DENSITY,c0,t0,c1, t1,p);

C_CENTROID(x0,c0,t0);

C_CENTROID(x1,c1,t1);

Subtract_Vector(dx,x1,x0);

Positive_Unit_Vector(dxp,dx);

mag_dx = Dot_Product(dx,dxp);

if (fabs(mag_dx) > MIN_REAL_NUMBER)

{

source += (thi*dens)*(C_UDSI(c1,t1,uds) - C_UDSI(c0,t0,uds))/mag_dx;

if (((thi*dens)/mag_dx) > 0.0) dS[eqn] -= (thi*dens)/mag_dx;

#ifdef SRM_DEBUG_FACEFLUX

if (cell == 1208)

{

Message("DD: c0 %d c1 %d, flux %e dens %e u0 %e u1 %e dx %e uf %e \n",c0,c1,

(thi*dens)*(C_UDSI(c1,t1,uds) - C_UDSI(c0,t0,uds))/mag_dx,

dens,C_UDSI(c0,t0,uds),C_UDSI(c1,t1,uds),

mag_dx,

thi*(C_UDSI(c1,t1,uds) - C_UDSI(c0,t0,uds))/mag_dx);

} #endif

}

}

/* ---- wall face thread */

/*

if (THREAD_TYPE(f_thread) == THREAD_F_WALL)

{

thi = fabs(Dot_Product(A,Diff_Coeff))*

(1.0 - Phase[p].Porosity)*Phase_Face_Density_(face,f_thread,p);

C_CENTROID(x0,cell,thread);

F_CENTROID(x1,face,f_thread);

Subtract_Vector(dx,x1,x0);

Positive_Unit_Vector(dxp,dx);

mag_dx = Dot_Product(dx,dxp);

if (fabs(mag_dx) > MIN_REAL_NUMBER)

{

source += thi*(F_UDSI(face,f_thread,uds) - C_UDSI(cell,thread,uds))/mag_dx;

if ((thi/mag_dx) > 0.0) dS[eqn] -= thi/mag_dx;

}

}

*/

}

}

#ifdef SRM_DEBUG_FACEFLUX

if (cell == 1208)

{

Message("DD: cell %d, sum %e, sum/dv %e \n",cell,source,source/CELL_VOLUME(cell,thread));

} #endif

source = source/CELL_VOLUME(cell,thread);

dS[eqn] = dS[eqn]/CELL_VOLUME(cell,thread);

}

else

source = 0.0;

return source; }

There are a couple of vector routines I wrote - which do pretty much what they say. Dot_Product, Subtract_Vector and Positive_Unit_Vector compute what they say.

Good Luck

Greg

Actually lookign more closely the source calculates

d/dx(rho_app * T * d(phi)/dx)

where rho_app = an apparent density. For your applications this can be replaced with 1 - ie. ignore.

Greg

One other thing about this code - I know it works, but since I wrote is for FLUENT 5.x the udf documentation has imprved a lot and you might be able to make it more efficient by using internal FLUENT Macros etc.

For example, there are FLUENT macros for the vector routines, which I have written (not shown) and called.

Also it may be useful to rewrite it using the _RG gradients suggested by Kim (below). This will speed up tinsg alot, since the astute reader will notice for neighbour cells my routine re-calculates the gradient on each side of the face for each cell. Thus the computations are much greater than they need be if these are stored. Anyway if you make improvements, please post them so I can use them too!!!

Greg