我做的是一個 用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;    /* 賦值移動速度給宏的參數，乘以時間步長即為此時間步發生的位移 */ 
}