多孔介质-UDF讨论

[wangzhixiong08]

/************************************************** **********************/ /* Porous Media Source Code */ /* */ /* By Greg Perkins */ /* CANCES, Australian Technology Park */ /* Ph: 02 9318 0004, Fax: 02 9319 2328 */ /* Email: perkinsg@cances.atp.com.au */ /* */ /* Started: 10-11-99 */ /* */ /* Revised: 22-11-99 */ /************************************************** **********************/
 
/* ---- These routines implement a porous media model (same as standard
 
Fluent model) that can vary spatially and temporily.
 
To use in Fluent compile and link source library as described
 
in the Fluent UDF manual. Then use the udf in the source terms
 
for momentum equations */
 
#include "udf.h" #include "sg.h"
 
/* ---- use a user defined scalar to track material properties of the domain */ enum { PERM, N_REQUIRED_UDS };
 
/* --- misc defines */ #define PERM_EXP_VOID 35.0 #define PERM_EXP_MATERIAL -20.0 #define PERM_ADJ_FALSE 0 #define PERM_ADJ_TRUE 1 #define UDS_PERM(alpha) log10(alpha) /* --- permeability function for scalar-0 */ #define WEIGHT 1.0e20 /* --- weighting co-eff in linearized eqn */ #define WEIGHT2 1.0e20 #define INLET_VELOCITY 10.0 /* --- inlet velocity in m/s */
 
/* ---- define criteria for various possible configurations */ #define MAT_CHANNEL(x,y,z) ((y<=0.25) || (y>=0.75)) #define MAT_STEP(x,y,z) ( ((x<=1.0)&&(y<=0.5)) || ((x>1.0)&&(y<=0.25)) ) #define MAT_ONECHANNEL(x,y,z) (y<=0.5) #define MAT_CHANNEL2(x,y,z) ( ((x<=1.0)&&(y<=0.5)) || ((x>1.0)&&(y<=0.25)) || (y>=0.95) )
 
/* ------------------------------------------------------------------------
 
Material Properties
 
This routine returns the local properties of the material (alpha,C2) for
 
a given location in the domain (x,y,z).
 
------------------------------------------------------------------------ */ void Material_Properties(real x, real y, real z, real *alpha, real *C2) {
 
/* --- use this routine to return the permeability, alpha and the
 
co-efficient C2, for each location x,y,z in the domain. No
 
Modifications needed for 2D. For time dependent porous media
 
use the RP_Get_Real("flow time") function to find out t (secs) */
 
if MAT_STEP(x,y,z)
 
{
 
*alpha = pow(10.0,PERM_EXP_MATERIAL);
 
*C2 = 0.0;
 
}
 
else
 
{
 
*alpha = pow(10.0,PERM_EXP_VOID);
 
*C2 = 0.0;
 
} }
 
/* ------------------------------------------------------------------------
 
X_Momentum_Source
 
This routine returns the source term for the X-momentum term for each
 
control volume in the domain. The local properties are obtained by
 
calling Material_Properties
 
------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Xmom,cell,thread,dS,eqn) {
 
real x[ND_ND];
 
real alpha, C2, constant1, constant2, Ux, source;
 
/* --- determine x,y,z co-ordinates of cell */
 
C_CENTROID(x,cell,thread);
 
/* --- determine local properties */
 
Material_Properties(x[0], x[1], x[2], &alpha, &C2);
 
/* --- determine constants 1,2 */
 
constant1 = C_MU_L(cell,thread)/alpha;
 
constant2 = 0.5 * C_R(cell,thread) * C2;
 
/* --- determine x-velocity */
 
Ux = C_U(cell,thread);
 
source = -(constant1*Ux + constant2 * fabs(Ux) * Ux);
 
dS[eqn] = -(constant1 + 2 * constant2 * fabs(Ux)); /* XXX CHECK */
 
return source; }
 
/* ------------------------------------------------------------------------
 
Y_Momentum_Source
 
This routine returns the source term for the Y-momentum term for each
 
control volume in the domain. The local properties are obtained by
 
calling Material_Properties
 
------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Ymom, cell, thread, dS, eqn) {
 
real x[ND_ND];
 
real alpha, C2, constant1, constant2, Uy, source;
 
/* --- determine x,y,z co-ordinates of cell */
 
C_CENTROID(x,cell,thread);
 
/* --- determine local properties */
 
Material_Properties(x[0], x[1], x[2], &alpha, &C2);
 
/* --- determine constants 1,2 */
 
constant1 = C_MU_L(cell,thread)/alpha;
 
constant2 = 0.5 * C_R(cell,thread) * C2;
 
/* --- determine y-velocity */
 
Uy = C_V(cell,thread);
 
source = -(constant1*Uy + constant2 * fabs(Uy) * Uy);
 
dS[eqn] = -(constant1 + 2 * constant2 * fabs(Uy)); /* XXX CHECK */
 
return source; }
 
/* ------------------------------------------------------------------------
 
Z_Momentum_Source
 
This routine returns the source term for the Z-momentum term for each
 
control volume in the domain. The local properties are obtained by
 
calling Material_Properties
 
------------------------------------------------------------------------ */ DEFINE_SOURCE(SRCE_Zmom, cell, thread, dS, eqn) {
 
real x[ND_ND];
 
real alpha, C2, constant1, constant2, Uz, source;
 
/* --- determine x,y,z co-ordinates of cell */
 
C_CENTROID(x,cell,thread);
 
/* --- determine local properties */
 
Material_Properties(x[0], x[1], x[2], &alpha, &C2);
 
/* --- determine constants 1,2 */
 
constant1 = C_MU_L(cell,thread)/alpha;
 
constant2 = 0.5 * C_R(cell,thread) * C2;
 
/* --- determine z-velocity */
 
Uz = C_W(cell,thread);
 
source = -(constant1*Uz + constant2 * fabs(Uz) * Uz);
 
dS[eqn] = -(constant1 + 2 * constant2 * fabs(Uz)); /* XXX CHECK */
 
return source; }

[l.j刘侃]

注意格式还有先调试一下吧 这个问题其实最难的是最后的各个方向上的分量处理
而如果是平面各向同性的话 难点在于处理最后源项中针对的方向长度问题
注意源项的单位 你好像并没有关心