|
楼主 |
发表于 2015-1-23 17:26:59
|
显示全部楼层
UDF如下
/***************************************************************************************
Slip_Thermal_BC.c
UDF
for slip velocity(Maxwell)+temperature jump
*****************************************************************************************/
#include "udf.h"
#include "sg.h"
#include <math.h>
/************************************************************
设置前处理命令
**********************************************************/
#define notSCHEMEMFP /*是否设置平均自由程?*/
#define SCHEMETMAC /*TMAC..*/
#define SCHEMEUDRLXCOEFF /*under-relaxation coefficient*/
#define notSCHEMEThAC
#define notSCHEMESpHR
#define SCHEMESIGMASQUARE
#define SCHEMEAMBPRESS
#define UNDERRLX
#define WALLMOTION
#define Boltzmann 1.3806505e-23
#define PI 3.14159265
#define SQRT_2 1.414213562
#ifndef SCHEMESIGMASQUARE
#define sigma_square 1.35305239e-19 /*分子直径平方*/
#endif
#ifndef SCHEMEAMBPRESS
#define ambpress 101325 /*大气压力*/
#endif
#ifndef SCHEMETMAC
#define TMAC 1.0 /*切向冲量调节系数*/
#endif
#ifndef SCHEMEThAC
#define ThAC 1.0 /*热调节系数*/
#endif
#ifndef SCHEMESpHR
#define SpHR 1.4 /*比热比*/
#endif
#ifndef SCHEMEUDRLXCOEFF
#define UDRLXCOEFF 0.2 /*松弛因子*/
#endif
/*
================================================
壁面速度滑移,x方向的程序,带热蠕变项
================================================
*/
DEFINE_PROFILE(temperature_jump,f_thread,index)
{
face_t face;
cell_t cell;
Thread *c_thread;
real Coeff2[ND_ND];
real Prandtl, gamma, temp, normal;
real MeanFreePath=6.8e-8;
#ifdef SCHEMEMFP
MeanFreePath=RP_Get_Real("meanfreepath");
#endif
#ifdef SCHEMEThAC
real ThAC;
ThAC=RP_Get_Real("thac");
#endif
#ifdef SCHEMEUDRLXCOEFF
real UDRLXCOEFF=0.2;
UDRLXCOEFF=RP_Get_Real("udrlxcoeff");
#endif
#ifdef SCHEMESpHR
real SpHR;
SpHR=RP_Get_Real("sphr");
#endif
coord_coeff(f_thread);
begin_f_loop(face,f_thread)
{
cell=F_C0(face,f_thread);
c_thread=THREAD_T0(f_thread);
#ifdef SCHEMEMFP
MeanFreePath = Boltzmann * F_T(face,f_thread)/(sigma_square * (F_P(face,f_thread)+ambpress) * PI * SQRT_2);
#endif
Prandtl=(C_MU_L(cell,c_thread)*C_CP(cell,c_thread))/C_K_L(cell,c_thread);
ND_SET(Coeff2[0],Coeff2[1],Coeff2[2],
F_UDMI(face,f_thread,3),
F_UDMI(face,f_thread,4),
F_UDMI(face,f_thread,5));
normal=NV_DOT(Coeff2,C_T_G(cell,c_thread));
gamma=(2*SpHR)/(SpHR+1);
temp=((2-ThAC)/ThAC) * gamma * MeanFreePath / Prandtl * normal;
#ifdef UNDERRLX
F_PROFILE(face,f_thread,index)=F_T(face,f_thread)*(1-UDRLXCOEFF)+(F_UDSI(face,f_thread,0)+temp)*UDRLXCOEFF;
#else
F_PROFILE(face,f_thread,index)=(F_UDSI(face,f_thread,0)+temp);
#endif
}
end_f_loop(f,thread)
}
DEFINE_PROFILE(maxwell_slip_velocity_x_full,f_thread,index)
{
face_t face;
cell_t cell;
Thread *c_thread;
real slip, thcreep, dveloc;
real normal_slip, tangential_slip, tangential_thcreep;
real Coeff1[ND_ND], Coeff2[ND_ND];
real u[ND_ND];
real TMAC=1;
real MeanFreePath=6.8e-8;
real UDRLXCOEFF=0.2;
real sigma_square=1.0e-19;
real ambpress;
#ifdef SCHEMETMAC
TMAC=RP_Get_Real("tmac");
#endif
#ifdef SCHEMEMFP
MeanFreePath=RP_Get_Real("meanfreepath");
#endif
#ifdef SCHEMEUDRLXCOEFF
UDRLXCOEFF=RP_Get_Real("udrlxcoeff");
#endif
#ifdef SCHEMESIGMASQUARE
sigma_square=RP_Get_Real("sigmasquare");
#endif
#ifdef SCHEMEAMBPRESS
ambpress=RP_Get_Real("ambpress");
#endif
if (RP_Get_Integer("coordeval")==1) coord_coeff(f_thread);
if ((RP_Get_Integer("tempgradeval")>1)&&(Data_Valid_P()))
{
begin_f_loop(face,f_thread)
{
cell=F_C0(face,f_thread);
c_thread=THREAD_T0(f_thread);
/*计算平均自由程*/
#ifndef SCHEMEMFP
MeanFreePath = Boltzmann * F_T(face,f_thread) / (sigma_square*(F_P(face,f_thread)+ambpress)*PI*SQRT_2);
#endif
/*速度坐标存入u[ND,ND]*/
ND_SET(u[0],u[1],u[2],
F_U(face,f_thread),
F_V(face,f_thread),
F_W(face,f_thread));
ND_SET(Coeff1[0],Coeff1[1],Coeff1[2],
F_UDMI(face,f_thread,0),
F_UDMI(face,f_thread,1),
F_UDMI(face,f_thread,2));
ND_SET(Coeff2[0],Coeff2[1],Coeff2[2],
F_UDMI(face,f_thread,3),
F_UDMI(face,f_thread,4),
F_UDMI(face,f_thread,5));
tangential_slip=NVD_DOT(Coeff1,
NV_DOT(Coeff1,C_U_G(cell,c_thread)),
NV_DOT(Coeff1,C_V_G(cell,c_thread)),
NV_DOT(Coeff1,C_W_G(cell,c_thread)));
normal_slip=-1 * NVD_DOT(Coeff1,
NV_DOT(Coeff2,C_U_G(cell,c_thread)),
NV_DOT(Coeff2,C_V_G(cell,c_thread)),
NV_DOT(Coeff2,C_W_G(cell,c_thread)));
slip = ((2-TMAC)/TMAC) * MeanFreePath * (tangential_slip+normal_slip);
tangential_thcreep=NV_DOT(Coeff1,C_T_G(cell,c_thread)) ;
thcreep = 0.75 * C_MU_L(cell,c_thread)/(C_R(cell,c_thread) * C_T(cell,c_thread)) * tangential_thcreep ;
#ifdef WALLMOTION
dveloc = Coeff1[0]* (F_UDSI(face,f_thread,1) + (slip + thcreep)) ;
#else
dveloc = Coeff1[0]* (slip + thcreep) ;
#endif
#ifdef UNDERRLX
dveloc = (1-UDRLXCOEFF) *u[0] + UDRLXCOEFF * dveloc;
#endif
F_PROFILE(face,f_thread,index) = dveloc;
}
end_f_loop(face,f_thread)
}
else
{
begin_f_loop(face,f_thread)
{
ND_SET(u[0],u[1],u[2],
F_U(face,f_thread),
F_V(face,f_thread),
F_W(face,f_thread)) ;
ND_SET(Coeff1[0],Coeff1[1],Coeff1[2],
F_UDMI(face,f_thread,0),
F_UDMI(face,f_thread,1),
F_UDMI(face,f_thread,2)) ;
#ifdef WALLMOTION
F_PROFILE(face,f_thread,index) = Coeff1[0] * F_UDSI(face,f_thread,1);
#else
F_PROFILE(face,f_thread,index) = Coeff1[0] * NV_DOT(u,Coeff1);
#endif
}
end_f_loop(face,f_thread)
}
}
void coord_coeff(Thread *f_thread)
{
face_t face;
cell_t cell;
Thread *c_thread;
real y[ND_ND];
real u[ND_ND], a;
real A[ND_ND];
real dr0[ND_ND], es[ND_ND], ds, A_by_es;
begin_f_loop(face,f_thread)
{
F_CENTROID(y,face,f_thread);
cell=F_C0(face,f_thread);
c_thread=THREAD_T0(f_thread);
BOUNDARY_FACE_GEOMETRY(face,f_thread,A,ds,es,A_by_es,dr0);
ND_SET(u[0],u[1],u[2],
C_U(cell,c_thread),
C_V(cell,c_thread),
C_W(cell,c_thread));
a=NV_MAG(u);
ND_SET(F_UDMI(face,f_thread,0),
F_UDMI(face,f_thread,1),
F_UDMI(face,f_thread,2),
NVD_DOT(u,1,0,0)/a,
NVD_DOT(u,0,1,0)/a,
NVD_DOT(u,0,0,1)/a);
ND_SET(F_UDMI(face,f_thread,3),
F_UDMI(face,f_thread,4),
F_UDMI(face,f_thread,5),
NVD_DOT(A,1,0,0)/A_by_es,
NVD_DOT(A,0,1,0)/A_by_es,
NVD_DOT(A,0,0,1)/A_by_es);
}
end_f_loop(face,f_thread)
}
|
|