Fluent UDF造波源程序集
法一：边界造波法
程序一：inlet.c
#include "udf.h" /* #include
#include
#include
#define HH 0.5 /*waver hight*/ /*不同波形需要修改的波形参数*/
#define LL 0.95 /*波长*/
#define g 9.81
#define pi 3.1415925
#define d 0.5 /*water deepth */
DEFINE_PROFILE(x_velocity,thread,index)
{
real kk = 2.0*pi/LL;
real c = sqrt(g*tanh(kk*d)/kk); /*参见一般性公式的推导*/
real TT = LL/c;
real x[ND_ND]; /* this will hold the position vector */
real y = 0;
real s = 0;
real ct = 0; /*相位角*/
face_t f;
real t = CURRENT_TIME;
real u = 0;
t=RP_Get_Real("flow-time");
begin_f_loop(f,thread) /* loops over all faces in the thread passed
in the DEFINE macro argument */ {
F_CENTROID(x,f,thread);
y = x[1];
s = y+d;
ct = kk*(x[0]-c*t); /*参见一般性公式的推导*/
if(y < 0.5*HH*sin(ct)) /*水面以下，其中，0.5*HH*sin(ct)为波面方程*/ u = pi*HH*cosh(kk*s)*sin(ct)/(TT*sinh(kk*d)); /*x方向速度分量公
式，参见一般性公式的推导*/ else u = 0.0; /*水面以上流体单位速度矢量的x方向分量*/
F_PROFILE(f,thread,index) = u;
}
end_f_loop(f,thread)
}
DEFINE_PROFILE(y_velocity,thread,index)
{
real kk = 2.0*pi/LL;
real c = sqrt(g*tanh(kk*d)/kk);
real TT = LL/c;
real x[ND_ND]; /* this will hold the position vector */
real y = 0;
real s = 0;
real ct = 0;
face_t f;
real t = CURRENT_TIME;
real v = 0;
t=RP_Get_Real("flow-time");
begin_f_loop(f,thread)
{
F_CENTROID(x,f,thread);
y = x[1];
s = y+d;
ct = kk*(x[0]-c*t);
if(y < 0.5*HH*sin(ct))
v = pi*HH*sinh(kk*s)*cos(ct)/(TT*sinh(kk*d));
else v = 0.0;
F_PROFILE(f,thread,index) = v;
}
end_f_loop(f,thread)
}
DEFINE_PROFILE(voffactor,thread,index)
{
real kk = 2.0*pi/LL;
real c = sqrt(g*tanh(kk*d)/kk);
real TT = LL/c;
real x[ND_ND]; /* this will hold the position vector */ real y = 0;
real s = 0;
real ct = 0;
face_t f;
real t = CURRENT_TIME;
t=RP_Get_Real("flow-time");
begin_f_loop(f,thread)
{
F_CENTROID(x,f,thread);
y = x[1];
s = y+d;
ct = kk*(x[0]-c*t);
if(y < 0.5*HH*sin(ct))
F_PROFILE(f,thread,index) = 1.0;
else
F_PROFILE(f,thread,index) = 0.0;
}
end_f_loop(f,thread)
程序二：wave.c （此程序同程序一大致相同）
#include "udf.h" /* must be at the beginning of every UDF you write new9 case */ real AA=0.5; /*waver amplitude*/
real LL=0.95; /*不同波形需要修改的波形参数*/
real TT=0.78;
real pi=3.1415926
real kk=2.0*pi/TT;
real ww=2.0*pi/LL;
real h=0.5; /*water deepth */
real ux=1.0; /*此为何变量？*/
DEFINE_PROFILE(x_velocity,thread,index)
{
real x[ND_ND]; /* this will hold the position vector */
real y;
face_t f;
real t = CURRENT_TIME;
real u;
begin_f_loop(f,thread) /* loops over all faces in the thread passed
in the DEFINE macro argument */ {
F_CENTROID(x,f,thread);
y = x[1];
if(y<(AA*cos(kk*x[0]-ww*t))) /*(-a*sin(w*t)+2.0为入口的波面随时间的变化*/
/*u=a*w*(exp(k*(y+h))+exp(-k*(y+h)))*-sin(w*t)/(exp(k*h)-exp(-k*h))*/
u=9.8*AA*kk/ww*cosh(kk*(y+h))*cos(kk*x[0]-ww*t)/cosh(kk*h);
else u =0.0; /*水面以上流体单位速度矢量的x方向分量*/
F_PROFILE(f,thread,index)=u;
}
end_f_loop(f,thread)
}
DEFINE_PROFILE(y_velocity,thread,index)
{
real x[ND_ND]; /* this will hold the position vector */
real y;
face_t f;
real t = CURRENT_TIME;
real v;
t=RP_Get_Real("flow-time");
begin_f_loop(f,thread) /* loops over all faces in the thread passed
in the DEFINE macro argument */ {
F_CENTROID(x,f,thread);
y = x[1];
if(y<(AA*cos(kk*x[0]-ww*t)))
v=9.8*AA*kk/ww*sinh(kk*(y+h))*sin(kk*x[0]-ww*t)/cosh(kk*h);
else v=0.0;
F_PROFILE(f,thread,index)=v;
}
end_f_loop(f,thread)
}
DEFINE_PROFILE(voffactor,thread,index)
{
real x[ND_ND]; /* this will hold the position vector */ real y;
face_t f;
real t = CURRENT_TIME;
begin_f_loop(f,thread) /* loops over all faces in the thread passed
in the DEFINE macro argument */ {
F_CENTROID(x,f,thread);
y = x[1];
if(y<(AA*cos(kk*x[0]-ww*t)))
F_PROFILE(f,thread,index) = 1.0;
else F_PROFILE(f,thread,index) = 0.0;
}
end_f_loop(f,thread)
}
法二：推波板造波法（动网格）
程序三：pushboard.c
#include
#include"udf.h"
#define T 5.8 /*周期*/
#define S 1.04 /*冲程*/
#define PI 3.1415926
DEFINE_CG_MOTION(pushboard, dt, cg_vel, cg_omega, time, dtime) {
real u = 0;
real ww = 0;
ww = 2*PI/T;
if(time <= 2*T)
u = ww*S*time*cos(ww*time)/(4*T);
else
u = ww*S*cos(ww*time)/2;
cg_vel[0] = u;
}
法三：摇板造波法
程序四：11.c （三维造波）
#include"udf.h"
#define k 0.77
#define w 2.75
#define A 0.059 /*振幅*/
#define l 1.6 /*波长*/
DEFINE_CG_MOTION(yaoban,dt,vel,omega,time,dtime)
{
vel[0]=0.0; /*x方向的速度*/
vel[1]=0.0;
vel[2]=0.0;
omega[0]=0.0;
omega[1]=0.0;
omega[2]=atan(A/l)*w*cos(w*time); /*z方向的角频率*/ }
程序五：b.c
#include "udf.h"
static real G=9.81;
static real VLA_M=100;
static real VLA_F=0.0;
static real VLA_V=0.0;
DEFINE_CG_MOTION(vla,dt,vel,omega,time,dtime)
{
Thread *t;
Domain *d;
real dv,CG[ND_ND],force[2],moment[2];
if(!Data_Valid_P())
return;
NV_S(vel,=,0);
NV_S(omega,=,0);
d=THREAD_DOMAIN (DT_THREAD ((Dynamic_Thread *)dt));
t=DT_THREAD(dt);
NV_S (CG, =, 0.0);
Compute_Force_And_Moment(d,t,CG,force,moment,FALSE);
dv = dtime * (force[1]-G*VLA_M)/ VLA_M;
VLA_V+= dv;
vel[1]=VLA_V;
VLA_F=force[1];
Message("time = %f F(y) = %f V(y) = %f\n", time, force[1], VLA_V); }