CFD 模拟二维平板以一个自由度的被动俯仰运动

我做的是一个 用fluent 模拟二维平板在流场中做被动俯仰运动，运动的UDF如下，求解方法为Simple，压力，动量等都采用二阶格式，流场网格的内外采用interface 界面，动网格设置为interface以内的区域和平板一起运动，时间步长取得 也比较小，但是开始计算以后，迭代了几步就开始发散。万能的网友，能不能给我一些建议啊？


          以下是我的UDF（解四阶龙格库塔方程来实现被动运动）


#include "udf.h" 
#include "sg_mem.h" 
#include "dynamesh_tools.h" /* 头文件，不可忽略 Header files  this enables both the definition of DEFINE macros and other ANSYS Fluent-provided macros and functions, and their inclusion in the compilation process*/

#define massr   0.0232 /* 定义全局变量，mass质量通过自主设定的m*求得*/
#define dtm    0.0004 /* 定义全局变量，dtm时间步长，与仿真中的时间步长一致 */
#define freq   42.25 /* freq我这里指定不是频率，是欧米伽的平方，在后面计算用 */

real b_vel=0; /* 初始化速度变量，是一个三维向量，初始化为零 */
real b_cen=0;
real disp;/* 初始化位置变量，是一个三维向量，初始化为零 */
real t; /* 定义一个时间变量t */
FILE *fp; /* 定义一个文件指针*fp */

DEFINE_EXECUTE_AT_END(exe_end)  /* UDF宏，在每一个时间步迭代结束后执行ANSYS Fluent calls the function at the end of an iteration */
{
int i;
real x_cg[3],f_glob[3],m_glob[3]; 
real vx,xn,Vx,Xn; 
real L1,L2,L3,L4;  /* 以上均为初始化局部变量 */

Domain *domain=Get_Domain(1); /* 流场指针，非多相流的话默认 */
Thread *tf=Lookup_Thread(domain,16);  /* 线指针，指向要获取信息的边界的ID，即柱的ID */

for(i=0; i<3; i++) /* 开始一个循环 */
{
x_cg=0;
f_glob=0.0;
m_glob=0.0; /* 此循环作用为初始化局部变量为0 */
}

disp=b_cen;/* 将上一时间步的位移信息赋值给变量 */

Compute_Force_And_Moment(domain,tf,x_cg,f_glob,m_glob,TRUE); /* 此命令为udf固有命令，可查阅到其各个参数的含义，此处为获取目标几何的位移，力和力矩，我们只用到位移和力，均为二维向量 */
  vx=b_vel; /* 将上一时间步的x向速度赋值给vx */
  xn=b_cen; /* 将上一时间步的x向位移赋值给xn */

  L1=m_glob[2]/massr-freq*xn-1.3*vx; 
  L2=m_glob[2]/massr-freq*(xn+vx*dtm/2)-1.3*(vx+L1*dtm/2); 
  L3=m_glob[2]/massr-freq*(xn+vx*dtm/2+dtm*dtm*L1/4)-1.3*(vx+dtm*L2/2); 
  L4=m_glob[2]/massr-freq*(xn+vx*dtm+dtm*dtm*L2/2)-1.3*(vx+dtm*L3);  /* 以上为四阶龙格库塔法解二阶振动微分方程的求解公式，写成C语言形式。注意此处振动微分方程的阻尼为0,因此求解中缺少一项，另外freq是欧米伽的平方 */
  
  Vx=vx+dtm*(L1+2*L2+2*L3+L4)/6; /* 求出此次时间步的x向速度 */
  Xn=xn+dtm*vx+dtm*dtm*(L1+L2+L3)/6;   /* 求出此次时间步的x向位移 */

  b_vel=Vx; /* 赋值 */
  b_cen=Xn; /* 赋值 */

  t=t+dtm; /* 时间叠加 */

  fp=fopen("output.txt","a"; /* 文件指针指向output.txt，打开方式为写入 */
  fprintf(fp,"%5f,%8f,%8f,%8f\n",t,disp,m_glob[2],b_vel);  /* 写入时间，位移，受力，速度 */
  fclose(fp);/* 关闭文件 */
}

DEFINE_CG_MOTION(cylinder,dt,vel,omega,time,dtime) /* UDF宏，作用为移动实体，各参数的具体含义见UDF手册 */
{ 
   NV_S(vel,=,0.0);  /* 初始化移动速度为0 */
   NV_S(omega,=,0.0); /* 初始化移动角速度为0 */
   omega[2]=b_vel;    /* 赋值移动速度给宏的参数，乘以时间步长即为此时间步发生的位移 */ 
}