|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?注册
x
#include "udf.h"
/*#include<stdio.h> (输入输出)*/
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define g 9.81
#define pi 3.1415926
#define Ts 1.39 /*周期*/
#define Hs 0.035 /*波高*/
#define d 0.52 /*水深*/
#define dx 0.01 /*网格长度*/
DEFINE_SOURCE(irregular,c,t,dS,eqn) /*质量源不规则造波*/
{
real x[ND_ND];
real y=x[1];
real s=y+d;
real SW[100];
real eps[100];
real K[100];
real w[100];
real Tp=Ts/0.93447;
real wp=2*pi/Tp;
real a;
real b;
real m;
real irr_source;
real tt=RP_Get_Real("flow-time"); /*返回当前计算时间*/
int i;
C_CENTROID(x,c,t); /*返回单元格重心的x值*/
irr_source=0.0; /*初始化源为0*/
for(i=0;i<=99;i++)
{
w=2.6+0.092*i; /*划分wi*/
if(w<=wp) /*J谱*/
SW=0.2189*pow(Hs,2)*pow(Tp,-4)*pow((w+0.046)/2/pi,-5)*exp(-1.25*pow(Tp*(w+0.046)/2/pi,-4))*pow(3.3,exp(-((w+0.046)/wp-1)*((w+0.046)/wp-1)/0.0098))/2/pi;
else
SW=0.2189*pow(Hs,2)*pow(Tp,-4)*pow((w+0.046)/2/pi,-5)*exp(-1.25*pow(Tp*(w+0.046)/2/pi,-4))*pow(3.3,exp(-((w+0.046)/wp-1)*((w+0.046)/wp-1)/0.0162))/2/pi;
a=0.0; b=15.0; m=(a+b)/2;
while(fabs(w*w-g*m*tanh(m*d))>0.00001 && fabs(a-b)>0.00001) /*二分法解色散方程求Ki*/
{
if((w*w-g*m*tanh(m*d))*(w*w-g*b*tanh(b*d))<0) a=m;
if((w*w-g*a*tanh(a*d))*(w*w-g*m*tanh(m*d))<0) b=m;
m=(a+b)/2;
}
K=m;
srand((unsigned)time(NULL));
eps=((double)rand())/RAND_MAX*2*pi; /*相位角θ在0~2π之间取随机数*/
irr_source=irr_source+2*pow(2*SW*0.092,0.5)*w*cosh(K*s)/sinh(K*d)*cos((w+0.046)*tt+eps)/dx;
}
dS[eqn] = 2/dx;
return irr_source;
} |
|