请大侠们帮我看看这段程序吧 我看不懂啊

王建伟

#include "udf.h"
/*
奉献一个动网格的udf。
该udf能实现流场力控制物体转动，同时能够实现入口速度随时间变化

*/
static real Iz=1.28e-8;
static real ph0=1.04719;
static real cg1=2.39e-3;
static real cg2=1.6659e-3;
static real cg3=0.0;

static real UP_PH1=0.0;
static real UP_W1= 0.0;
static real UP_PH2=0.0;
static real UP_W2= 0.0;
static real a0=0.0;
static real a1=0.0;


static real DW_PH1=0.0;
static real DW_W1= 0.0;
static real DW_PH2=0.0;
static real DW_W2= 0.0;
static real b0=0.0;
static real b1=0.0;

DEFINE_CG_MOTION(up, dt, vel, omega, time, dtime)
{
  Thread *t;
  Domain *d;
  real dwp,CG[3],force[3],moment[3];

  NV_S (vel, =, 0.0);
  NV_S (omega, =, 0.0);
  t = DT_THREAD (dt);
  d = THREAD_DOMAIN (DT_THREAD ((Dynamic_Thread *)dt));

  CG[0]=2.39e-3;
  CG[1]=1.6659e-3;
  CG[2]=0.0;

  Compute_Force_And_Moment(d,t,CG,force,moment,FALSE);


          a1=moment[2]/Iz;
          dwp = dtime * ((a1+a0)/2.0);
          a0=a1;

          UP_W1+= dwp;
      UP_PH1+=UP_W1*dtime;
      
  if(UP_PH1>=-1*ph0)
  {
   
    omega[2]=UP_W1;
  }

  if(UP_PH1<-1*ph0)

  {
          UP_W2+= dwp;

          if(UP_W2<=0.0)
          {
          omega[2]=0.0;
          UP_W2=0.0;
          }

          if(UP_W2>0.0)

          {
           UP_PH2+=UP_W2*dtime;

       if(UP_PH2<=1*ph0)
          omega[2]=UP_W2;
       if(UP_PH2>1*ph0)
              omega[2]=0.0;
          }
  }
  Message0("time = %.5e up_dv=%.4e up_PH = %.4e up_W=%.4e M_up = %f.6e\n", time, dwp,UP_PH2, UP_W2, moment[2]);
}

DEFINE_CG_MOTION(down, dt, vel, omega, time, dtime)
{
  Thread *t;
  Domain *d;
  real dwd,CG[3],force[3],moment[3];

  NV_S (vel, =, 0.0);
  NV_S (omega, =, 0.0);
  t = DT_THREAD (dt);
  d = THREAD_DOMAIN (DT_THREAD ((Dynamic_Thread *)dt));

  CG[0]=2.39e-3;
  CG[1]=-1.6659e-3;
  CG[2]=0.0;

  Compute_Force_And_Moment(d,t,CG,force,moment,FALSE);

          b1=moment[2]/Iz;
          dwd = dtime * ((b1+b0)/2.0);
          b0=b1;

          DW_W1+= dwd;
      DW_PH1+=DW_W1*dtime;
     
  if(DW_PH1<=ph0)
  {
          omega[2]=DW_W1;
  }
  if(DW_PH1>ph0)
  {
      DW_W2+= dwd;
          if(DW_W2>=0.0)
          {
      omega[2]=0.0;
          DW_W2=0.0;
          }
      if(DW_W2<0.0)
  {
   
    DW_PH2+=DW_W2*dtime;

        if(DW_PH2>=-1*ph0)
      omega[2]=DW_W2;
    if(UP_PH2<-1*ph0)
          omega[2]=0.0;
  }
  }

  
  Message0("dwd=%.4e dw_PH = %.4e dw_W=%.4e M_dw = %f.6e\n", dwd,DW_PH1, DW_W1, moment[2]);
}


DEFINE_PROFILE(ve_udf, t, i)
{
  face_t f ;

  real vt = RP_Get_Real("flow-time");

if(vt<=0.1)

     begin_f_loop(f,t)
       {
         F_PROFILE(f,t,i) = 0.736*vt ;
       }
    end_f_loop(f,t)

if(vt>0.1&vt<=0.2)

begin_f_loop(f,t)
       {
         F_PROFILE(f,t,i) = 8.832*vt-0.8096 ;
       }
    end_f_loop(f,t)

if(vt>0.2&vt<=0.4)

begin_f_loop(f,t)
       {
         F_PROFILE(f,t,i) = -6.24*vt+2.2816 ;
       }
    end_f_loop(f,t)

}