CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > FLUENT

Directional Diffusivity of UDS Source

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   September 4, 2003, 02:34
Default Directional Diffusivity of UDS Source
  #1
Greg Perkins
Guest
 
Posts: n/a
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
  Reply With Quote

Old   September 4, 2003, 02:40
Default Re: Directional Diffusivity of UDS Source
  #2
Greg Perkins
Guest
 
Posts: n/a
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
  Reply With Quote

Old   September 4, 2003, 03:02
Default Re: Directional Diffusivity of UDS Source
  #3
Greg Perkins
Guest
 
Posts: n/a
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
  Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
swak4foam building problem GGerber OpenFOAM Installation 54 April 24, 2015 16:02
DxFoam reader update hjasak OpenFOAM Post-Processing 69 April 24, 2008 01:24
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51
UDF Scalar Code: HT 1 Greg Perkins FLUENT 8 October 20, 2000 12:40
UDFs for Scalar Eqn - Fluid/Solid HT Greg Perkins FLUENT 0 October 11, 2000 03:43


All times are GMT -4. The time now is 11:22.