望结构外部绕流学习者解决难题啊

1987guochunping

我在fluent的udf动网格中嵌入振动方程（用纽马克法表示出来的），在flunet中提取出升力带入振动方程，解出每时刻的位移和速度，但我在fluent中采集出来的升力曲线很不理想，论文中体现的是正弦曲线那样。望各位指点指点。下面是我编的程序：
     #include"stdio.h"
#include"udf.h"
#include"unsteady.h"
static real v_prev=0.0;
#define k  760.832
#define M  1
#define C  0.005
#define ep0 1e-10
DEFINE_CG_MOTION(shuwan,dt,vel,omega,time,dtime)
{               Thread *t=DT_THREAD(dt);
                 face_t f;
                             real NV_VEC(A);
                                 real force,dv;
                                 real K;
                                 real P;
                                 real dy;
                                 real a;
                                 real disp,loc;
                                 real dforce;
                                 real li;
                                 real ep;
                                 NV_S(vel,=,0.0);
                                 NV_S(omega,=,0.0);
                                 if(!Data_Valid_P())
                                         return;
                                  force=0.0;
                                 begin_f_loop(f,t)
                                 {  
                                     F_AREA(A,f,t);
                                         force+=F_P(f,t)*A[1];
                                 }
                                 end_f_loop(f,t)
              read_loc_v_prev_file(&loc,&v_prev);
                                      ep=1;
                                          disp=0.0;
                      dforce=0.0;
                                          K=0.0;
                     li=force-dforce;
                                         dforce=force;
                              
                                 while(ep>ep0)
                                         { a=(force-C*v_prev-K*disp)/M;
                     K=k+6*M/(0.5*dtime*dtime)+3*C/dtime;
                     P=li+M*(v_prev/(0.5*dtime)+a)+C*(v_prev-0.5*dtime*a);
                                         dy=P/K;
                                         dv=dy/dtime-v_prev+0.5*dtime*a;
                                         v_prev=dy/dtime+0.5*dtime*a;
                                         disp+=dy;
                                         loc+=v_prev*dtime;
                     ep=abs(disp-loc);
                                         }
                               
                                     Message("\nUDF value:time=%f,y_velo=%f,force=%f\n",time,v_prev,force);
                                       
                        write_loc_v_prev_file(loc,v_prev);
                               
#if PARALLEL
                                         host_to_node_real_1(v_prev);
#endif           
                                         vel[1]=v_prev;

                                        }