How can we model forced diffusion in binary mixtures (single phase)?
|
ANALYSIS
J = J_ordinary + J_forced J1_ordinary = -rho D12 Grad yi J1_forced = -rho^2 * w1 * w2 * y1 * y2 * D12 * (F1 - F2) / (P * W * W) where w1 = mol wt 1 w2 = mol wt 2 y1 = mass fract 1 y2 = 1-y1 W = mol wt mixture F1 = Force Vector on species 1 (N/kg) F2 = Force Vector on species 2 (N/kg) P = pressure This is solved by UDF species equation source is equal to S = - DIV (J1_forced) Problem is solved using UDS: UDS(0) = y1 UDS(1) = xcomp of J1_forced UDS(2) = y-comp of J1_forced UDS(3) = z-comp of J1_forced In Define Adjust function, we * compute J_forced. * compute grad of components of J_forced (nine terms stored in memory)( * add appropriate components to get divergence: div J = d(J1x)/dx + D(J1y)/dy + d(J1z)/dz Source term = Div J is stored in UDM(0) For boundary condition, at a wall: J.n = 0 J = J_ord + J_forced, thus Flux(1) = J1_ord . n = - J_forced . n Where J1_forced is evaluated at the wall So far only wall boundary conditions are in the UDF. Others can be added as necessary. /******************************************************************** UDF to compute forced diffusion in FLUENT 6.1 by William Wangard, Ph.D. Fluent, Inc. All rights reserved, 2003. Fluent, Inc. does not guarantee the rigor of this User defined function BINARY mixtures only. UDSs are used to store species: UDS(0) = species 0 mass fraction species 1 mass fraction = 1 - UDS(0) UDS(1) = x component of forced diffusion flux of species 1 UDS(2) = y component .... UDS(3) = z component .... UDM(0) = DIV (forced diffusion flux) NOTE the variable FAC is a control to underrelax the body force, since it may be large. Adjust it as necessary to converge the solver Recommended settings, * Body-force weighted pressure scheme * Body force underrelax factor = 0.1 * Mom underrleax factor = 0.1 * Species 0 underrelax factor = 0.9 Convergence criteria for UDS 0 should be < 1E-6. This UDF should work in PARALLEL, but it has not been tested. Address other questions to ww@fluent.com Thank you. **********************************************************************************/ #include "udf.h" #include "sg.h" #define NSPECIES 2 /* D12 of Oxygen in air is 2E-5 m2/s */ /* In gui, define Diffusivity as D12 * rho_air */ #define D12 2E-5 /* Underrelaxation factor for body forces */ #define FAC 0.01 real wt[NSPECIES] = {32.0, 28.0134}; real b0[3] = {10000., 10000., 0.}; real b1[3] = {0., 0., 0.}; enum { M = NSPECIES - 1, N_UDS_REQ = M + ND_ND }; enum { N_UDM_REQ = ND_ND }; DEFINE_PROFILE(flux_1,tf,i) { Thread* tc; cell_t c; face_t f; real NV_VEC(normal), NV_VEC(area), Amag; real p,w,x,flux; real NV_VEC(d1), NV_VEC(d2); begin_f_loop(f,tf) { F_AREA(area,f,tf); Amag = NV_MAG(area); NV_VS(normal,=,area,/,Amag); c = F_C0(f,tf); tc = F_C0_THREAD(f,tf); p = ABS_P(C_P(c,tc),op_pres); w = 1.0 / ( F_UDSI(f,tf,0) / wt[0] + (1.- F_UDSI(f,tf,0)) / wt[1]); x = C_R(c,tc) * C_R(c,tc) * wt[0] * wt[1] * F_UDSI(f,tf,0) * (1.0 - F_UDSI(f,tf,0)) * D12 / (p * w * w); NV_VV(d1,=,b0,-,b1); NV_VS(d2,=,d1,*,x); flux = NV_DOT(d2,normal); /* Message("normal = %10.3e, %10.3e flux = %10.3e n", normal[0], normal[1], flux);*/ F_PROFILE(f,tf,i) = flux; } end_f_loop(f,tf); } DEFINE_ADJUST(compute_forced_diffusion,d) { real w; int n, id, ns = 0; Thread* tc; face_t f; cell_t c; Thread* tf; real p; real x; if (n_uds < N_UDS_REQ) Internal_Error("You have defined not enough UDSsn"); if (sg_udm < N_UDM_REQ) Internal_Error("Not enough UDMs defined.n"); /* Put */ thread_loop_c(tc,d) { begin_c_loop(c,tc) { p = ABS_P(C_P(c,tc),op_pres); w = 1.0 / ( C_UDSI(c,tc,0) / wt[0] + (1.-C_UDSI(c,tc,0)) / wt[1]); x = (C_R(c,tc) * C_R(c,tc) * wt[0] * wt[1] * C_UDSI(c,tc,0) * (1.-C_UDSI(c,tc,0)) * D12 ) / ( p * w * w ); for(n=0; n<ND_ND; n++) C_UDSI(c,tc,M+n) = x * (b0[n] - b1[n]); } end_c_loop(c,tc); } thread_loop_f(tf,d) { if (THREAD_TYPE(tf) == THREAD_F_WALL) { begin_f_loop(f,tf) { id = THREAD_ID(tf); c = F_C0(f,tf); tc = F_C0_THREAD(f,tf); p = ABS_P(C_P(c,tc),op_pres); w = 1.0 / ( F_UDSI(f,tf,0) / wt[0] + (1.- F_UDSI(f,tf,0)) / wt[1]); x = C_R(c,tc) * C_R(c,tc) * wt[0] * wt[1] * F_UDSI(f,tf,0) * (1.0 - F_UDSI(f,tf,0)) * D12 / (p * w * w); for(n=0; n<ND_ND; n++) F_UDSI(f,tf,M+n) = x * (b0[n] - b1[n]); /* Message("Thread %i, y0 = %10.3e y1 = %10.3e p = %10.3e w = %10.3e x = %10.3en", id, F_UDSI(f,tf,0), (1.- F_UDSI(f,tf,0)), p, w, x); */ } end_f_loop(f,tf); } } /* Compute Gradient of UDSs */ if (sg_uds) { for(n=M+0; n<M+ND_ND; n++) { MD_Alloc_Storage_Vars(d, SV_UDSI_RG(n),SV_UDSI_G(n),SV_NULL); Scalar_Reconstruction(d, SV_UDS_I(n), -1, SV_UDSI_RG(n), NULL); Scalar_Derivatives(d, SV_UDS_I(n), -1, SV_UDSI_G(n), SV_UDSI_RG(n), NULL); } } /* Now, UDS(0) = y_0 (as always, y_1 = 1.0 - y_0) UDS(1) = jf_0_x UDS(2) = jf_0_y UDS(3) = jf_0_z C_UDSI_G(1)[0] = d(jf_0_x)/dx C_UDSI_G(1)[1] = d(jf_0_x)/dy C_UDSI_G(1)[2] = d(jf_0_x)/dz C_UDSI_G(2)[0] = d(jf_0_y)/dx C_UDSI_G(2)[1] = d(jf_0_y)/dy C_UDSI_G(2)[2] = d(jf_0_y)/dz C_UDSI_G(2)[0] = d(jf_0_z)/dx C_UDSI_G(2)[1] = d(jf_0_z)/dy C_UDSI_G(2)[2] = d(jf_0_z)/dz In species, eqn, source = -div (j0) for eqn 0 source = UDM(0) = -( d(jf_0_x)/dx + d(jf_0_y)/dy + d(jf_0_z)/dz) */ thread_loop_c(tc,d) { begin_c_loop(c,tc) { C_UDMI(c,tc,0) = 0.; for(n=0; n<ND_ND; ++n) C_UDMI(c,tc,0) -= C_UDSI_G(c,tc,M+n)[n]; } end_c_loop(c,tc); } } DEFINE_SOURCE(uds_source,c,t,dS,eqn) { real source; dS[eqn] = 0.0; source = C_UDMI(c,t,0); return source; } DEFINE_SOURCE(x_momentum_source,c,t,dS,eqn) { real source, fac; dS[eqn] = 0.0; fac = FAC; source = C_R(c,t) * (C_UDSI(c,t,0) * b0[0] + (1.-C_UDSI(c,t,0)) * b1[0]); return source; } DEFINE_SOURCE(y_momentum_source,c,t,dS,eqn) { real source, fac; dS[eqn] = 0.0; fac = FAC; source = C_R(c,t) * (C_UDSI(c,t,0) * b0[1] + (1.-C_UDSI(c,t,0)) * b1[1]); return source; } DEFINE_SOURCE(z_momentum_source,c,t,dS,eqn) { real source, fac; dS[eqn] = 0.0; fac = FAC; source = C_R(c,t) * (C_UDSI(c,t,0) * b0[2] + (1.-C_UDSI(c,t,0)) * b1[2]); return source; } |
||
|