|
[Sponsors] |
September 4, 2003, 02:34 |
Directional Diffusivity of UDS Source
|
#1 |
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 |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
[swak4Foam] swak4foam building problem | GGerber | OpenFOAM Community Contributions | 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 |