找回密码
 注册
查看: 2155|回复: 0

不规则波造波udf,帮忙看看有什么问题

[复制链接]
发表于 2013-3-31 20:39:55 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?注册

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;
}
您需要登录后才可以回帖 登录 | 注册

本版积分规则

快速回复 返回顶部 返回列表