FLUENT 6 - DEFINE_VECTOR_EXCHANGE_PROPERTY example
|
The use of DEFINE_VECTOR_EXCHANGE_PROPERTY is difficult, because the velocities that are needed for the calculation of the drag coefficient are not in the data fields they are supposed to be. At the time when the UDF is called, C_U etc. do NOT contain what they should contain. Therefore, these data need to be saved to UDM locations first.
Please see also Solution 455, which contains a rather similar example for DEFINE_EXCHANGE_PROPERTY. /* Arbitrary Drag coefficients for mixture multiphase model */ /* Later development stages of this may be found in "jodegas.c"..! */ /* This may FAIL for MOVING reference frames and zones -- not tested!! */ /* This text offers three drag coefficients, marked by mx, my, mz. */ /* They are calculated using the following general form: % % c_D = c_{D0} + 1.0 * 24 / Re + 6.0 / (1. + sqrt(Re)) % | | | % my_cd my_fd my_tt % % Above are shown the RP variables bound to the parameters. */ /* Before using this, please enter the following in FLUENT6: (rp-var-define 'mx-cd 0.7 'real #f) (rp-var-define 'mx-fd 1. 'real #f) ;;; AGM - Felt (rp-var-define 'mx-tt 6. 'real #f) (rp-var-define 'my-cd 0.5 'real #f) (rp-var-define 'my-fd 1. 'real #f) ;;; Air - AGM (rp-var-define 'my-tt 6. 'real #f) (rp-var-define 'mz-cd 10. 'real #f) (rp-var-define 'mz-fd 1. 'real #f) ;;; Air - Felt (rp-var-define 'mz-tt 6. 'real #f) (rp-var-define 'joslip-first-udm-in-use 0 'int #f) Set individual variables using... (rpsetvar 'my_fd 1.) etc. Check variables using... (rpgetvar 'my_fd) etc. ********************************** */ #include "udf.h" #ifdef INTERPRETER Source code level macro naming problem with "INTERPRETER"...! #endif #ifdef PRINT Source code level macro naming problem with "PRINT"...! #endif #ifdef MyError Source code level macro naming problem with "MyError"...! #endif #ifdef STRUCT_REF #define INTERPRETER 1 #define PRINT printf #define MyError PRINT #define Boolean int /* Interpreter doesn't accept char here...:-( */ #else #define INTERPRETER 0 #define PRINT CX_Message #define MyError Error #endif /* Original in mem.h: */ /* #define C_TYPE(c,t)C_STORAGE(c,t,SV_TYPE,uchar_fl *) */ /* This is used in c_face_loop and causes the interpreter to complain: */ /* "non-integer subscript expression: unsigned char" */ /* Therefore, I'm trying a re-definition -- I DON't know if this works! */ #if INTERPRETER #undef C_TYPE #define C_TYPE(c,t)C_STORAGE(c,t,SV_TYPE,int *) #endif /* ================================ */ int joslip_first_udm_in_use = 32000; real mx_cd = 0.; real mx_fd = 0.; real mx_tt = 0.; real my_cd = 0.; real my_fd = 0.; real my_tt = 0.; real mz_cd = 0.; real mz_fd = 0.; real mz_tt = 0.; real drift_relax = 0.; #if INTERPRETER int degas_warnNudm_done = 0; #endif void degas_checkNudm(int minimum_number) { /* Check number of allocated User Defined Memory Locations..: */ #if ! INTERPRETER if (sg_udm < minimum_number) { Error("Please allocate at least %d User Defined Memory Locations!n", minimum_number); } #else if (! degas_warnNudm_done) { PRINT( "Make sure you have allocated at least %d User Defined Memory Locations..n", minimum_number); degas_warnNudm_done = 1; } #endif } DEFINE_ADJUST(adjust_myslip, domain) { Thread *t, *subt; Domain *subd; cell_t c; int phind; joslip_first_udm_in_use = RP_Get_Integer("joslip-first-udm-in-use"); for (subd = DOMAIN_SUB_DOMAIN(domain,(phind=0)); NNULLP(subd); subd = DOMAIN_SUB_DOMAIN(domain,(++phind))) /*..*/ ; degas_checkNudm(ND_ND * phind + joslip_first_udm_in_use); mx_cd = RP_Get_Real("mx-cd"); mx_fd = RP_Get_Real("mx-fd"); mx_tt = RP_Get_Real("mx-tt"); my_cd = RP_Get_Real("my-cd"); my_fd = RP_Get_Real("my-fd"); my_tt = RP_Get_Real("my-tt"); mz_cd = RP_Get_Real("mz-cd"); mz_fd = RP_Get_Real("mz-fd"); mz_tt = RP_Get_Real("mz-tt"); drift_relax = RP_Get_Real("drift/relax"); thread_loop_c(t, domain) { for (subt = THREAD_SUB_THREAD(t,(phind=0)); NNULLP(subt); subt = THREAD_SUB_THREAD(t,(++phind))) /*..*/ { begin_c_loop(c,subt) { ND_SET(C_UDMI(c,t,phind * ND_ND + 0 + joslip_first_udm_in_use), C_UDMI(c,t,phind * ND_ND + 1 + joslip_first_udm_in_use), C_UDMI(c,t,phind * ND_ND + 2 + joslip_first_udm_in_use), C_U(c,subt), C_V(c,subt), C_W(c,subt)); } end_c_loop(c,subt) } } } void m_slip(cell_t cell, Thread *mxtr_thrd, int scnd_clmn, int frst_clmn, real (*fd_function)(real), real *rslt_vect) { Thread *scnd_thrd = THREAD_SUB_THREAD(mxtr_thrd, scnd_clmn); Thread *frst_thrd = THREAD_SUB_THREAD(mxtr_thrd, frst_clmn); /* FiRST is the dispersed (granular) phase, SeCoND is the primary phase. */ /* The numbering refers to the phase columns in the "Interaction" panel. */ /*real diam = generic_property(cell, frst_thrd, %**** DOMAIN_PROPERTY(THREAD_DOMAIN(frst_thrd)), PROP_diameter, 293.); *** */ real diam = C_PHASE_DIAMETER(cell, frst_thrd); real diaQ = diam * diam; real NV_VEC(slipVel), velo, Re, f_d, tau; real ND_VEC(a_x, a_y, a_z); real ND_VEC(*frstU, *frstV, *frstW), ND_VEC(*scndU, *scndV, *scndW); real NV_VEC(torelax); real minus_relax = 1. - drift_relax; if (0. == drift_relax) { MyError("m_slip(): drift_relax may not be zero..!n"); } ND_SET(frstU, frstV, frstW, T_STORAGE_R_XV(mxtr_thrd, SV_UDM_I, frst_clmn * ND_ND + 0 + joslip_first_udm_in_use), T_STORAGE_R_XV(mxtr_thrd, SV_UDM_I, frst_clmn * ND_ND + 1 + joslip_first_udm_in_use), T_STORAGE_R_XV(mxtr_thrd, SV_UDM_I, frst_clmn * ND_ND + 2 + joslip_first_udm_in_use)); ND_SET(scndU, scndV, scndW, T_STORAGE_R_XV(mxtr_thrd, SV_UDM_I, scnd_clmn * ND_ND + 0 + joslip_first_udm_in_use), T_STORAGE_R_XV(mxtr_thrd, SV_UDM_I, scnd_clmn * ND_ND + 1 + joslip_first_udm_in_use), T_STORAGE_R_XV(mxtr_thrd, SV_UDM_I, scnd_clmn * ND_ND + 2 + joslip_first_udm_in_use)); /* Calculate "momentum accelaration" (acceleration due to "convective change") due to flow cross section variation: */ ND_VEC( a_x = NVD_DOT(C_U_G(cell,mxtr_thrd), C_U(cell,mxtr_thrd), C_V(cell,mxtr_thrd), C_W(cell,mxtr_thrd)), a_y = NVD_DOT(C_V_G(cell,mxtr_thrd), C_U(cell,mxtr_thrd), C_V(cell,mxtr_thrd), C_W(cell,mxtr_thrd)), a_z = NVD_DOT(C_W_G(cell,mxtr_thrd), C_U(cell,mxtr_thrd), C_V(cell,mxtr_thrd), C_W(cell,mxtr_thrd))); /* Is the negative sign correct??? No, it wasn't -- cancelled...! */ /* This should be tested by looking at a bubbly multi-phase flow */ /* that experiences some horizontal acceleration (rotation...!) */ /* Add gravity term...: */ if (M_gravity_p) { ND_V(a_x, a_y, a_z, -=, M_gravity); } /* Add unsteady term... */ if (rp_unsteady) { real physical_time_step = RP_Get_Real("physical-time-step"); ND_DS(a_x, a_y, a_z, +=, C_U(cell,mxtr_thrd) - C_U_M1(cell,mxtr_thrd), C_V(cell,mxtr_thrd) - C_V_M1(cell,mxtr_thrd), C_W(cell,mxtr_thrd) - C_W_M1(cell,mxtr_thrd), /, physical_time_step); /* Is the negative sign (in "-=") correct? No, it wasn't -- changed. */ } /* Now, the complete accelaration is in (a_x,a_y,a_z). */ /* Next, calculate the particle relaxation time...: */ /* Calculate [old] slip velocity in order to calculate Reynolds number: */ /* Phase velocities are not available at the time when this runs...! */ /*NV_DD(V, =, C_U(cell, frst_thrd), C_V(cell, frst_thrd), C_W(cell, frst_thrd), %% -, C_U(cell, scnd_thrd), C_V(cell, scnd_thrd), C_W(cell, scnd_thrd)); */ NV_DD(slipVel, =, frstU[cell], frstV[cell], frstW[cell], -, scndU[cell], scndV[cell], scndW[cell]); /* Build magnitude of [old] slip velocity: */ velo = NV_MAG(slipVel); /* Calculated Reynolds number with this [old] slip velocity magnitude: */ Re = RE_NUMBER(C_R(cell, scnd_thrd), velo, diam, C_MU_L(cell, scnd_thrd)); /* Some definitions of drag coefficients..: */ /* Laminar is: C_D = 24. / Re */ /* Turbulent : C_D = const. */ /* we define : f_d = C_D * Re / 24. --- Thus..: */ /* Laminar is: f_d = 1. */ /* Turbulent : f_d ~ Re / 24. */ f_d = (*fd_function)(Re); /*f_d = (mycd + mytt / (1 - sqrt(Re))) * Re / 24. + myfd;*/ /* (mycd, mytt, and myfd are taken from user settings..) */ /* Now calculate the particle relaxation time for the mixture model: */ /* (This has a different definition than the general part. rel. time!) */ tau = (C_R(cell,mxtr_thrd) - C_R(cell,frst_thrd)) * diaQ /*..*/ / (18. * C_MU_L(cell,scnd_thrd) * f_d); /* With this,...: d u_p / d t = 1/tau * (u - u_p) */ NV_DS(torelax, =, a_x, a_y, a_z, *, tau); NV_VS_VS(rslt_vect, =, torelax, *, drift_relax, +, slipVel, *, minus_relax); } real schiller_fd_func(real Re) { if (Re > 1000.) { return 0.44 * Re / 24.; } else { return (1. + 0.15 * pow(Re, 0.687)); } } DEFINE_VECTOR_EXCHANGE_PROPERTY(schi_slip, cell, mxtr_thrd, scnd_clmn, frst_clmn, rslt_vect) { m_slip(cell, mxtr_thrd, scnd_clmn, frst_clmn, schiller_fd_func, rslt_vect); } real unity_fd_func(real Re) { return 1.; } DEFINE_VECTOR_EXCHANGE_PROPERTY(unit_slip, cell, mxtr_thrd, scnd_clmn, frst_clmn, rslt_vect) { m_slip(cell, mxtr_thrd, scnd_clmn, frst_clmn, unity_fd_func, rslt_vect); } real mx_fd_func(real Re) { return (mx_cd + mx_tt / (1 - sqrt(Re))) * Re / 24. + mx_fd; } DEFINE_VECTOR_EXCHANGE_PROPERTY(mxslip, cell, mxtr_thrd, scnd_clmn, frst_clmn, rslt_vect) { m_slip(cell, mxtr_thrd, scnd_clmn, frst_clmn, mx_fd_func, rslt_vect); } real my_fd_func(real Re) { return (my_cd + my_tt / (1 - sqrt(Re))) * Re / 24. + my_fd; } DEFINE_VECTOR_EXCHANGE_PROPERTY(myslip, cell, mxtr_thrd, scnd_clmn, frst_clmn, rslt_vect) { m_slip(cell, mxtr_thrd, scnd_clmn, frst_clmn, my_fd_func, rslt_vect); } real mz_fd_func(real Re) { return (mz_cd + mz_tt / (1 - sqrt(Re))) * Re / 24. + mz_fd; } DEFINE_VECTOR_EXCHANGE_PROPERTY(mzslip, cell, mxtr_thrd, scnd_clmn, frst_clmn, rslt_vect) { m_slip(cell, mxtr_thrd, scnd_clmn, frst_clmn, mz_fd_func, rslt_vect); } |
||
|