# Directional Diffusivity of UDS Source

 Register Blogs Members List Search Today's Posts Mark Forums Read September 4, 2003, 02:34 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  September 4, 2003, 02:40 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  September 4, 2003, 03:02 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 Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post [swak4Foam] swak4foam building problem GGerber OpenFOAM Community Contributions 54 April 24, 2015 16:02 hjasak OpenFOAM Post-Processing 69 April 24, 2008 01:24 jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51 Greg Perkins FLUENT 8 October 20, 2000 12:40 Greg Perkins FLUENT 0 October 11, 2000 03:43

All times are GMT -4. The time now is 13:42.