[求助]大家帮忙看一下这个流固耦合的UDF

xinxinxiangsun

仿照zhoushuaihd做的风振涡致振动的例子，弹簧振子体系的运动UDF，但是在计算的一开始圆柱就跑到很远的地方去了，导致负网格，大家帮忙看一下是什么问题，非常感谢
/*coupling of fluid and solid,single degree equation, revised by Cao*/
#include "udf.h"
#include "math.h"
#include "stdio.h"
static real vel_prev;
static real disp_prev;
static real acel_prev;
static real x_cg;
/*static real body_centroid[3]={0.275,0.275,0.0}; */
DEFINE_CG_MOTION(wall_udf,dt,vel,omega,time,dtime)
{
Thread *t;
face_t f;
real NV_VEC (A);
real force, dv;
real beta,gama,m,k,damp_ratio,wn,wd,c,kkk;
real b0,b1,b2,b3,b4,b5,b6,b7,P,delt;
real disp_old,vel_old,acel_old,force_old;
real aa[4];
int i,j;
FILE *out,*fp1;

aa[0]=0.0;
aa[1]=0.0;
aa[2]=0.0;
aa[3]=0.0;
/*make sure the suitable of the output file*/
#if !RP_NODE
if((fp1=fopen("iternation values","r"))==NULL)
{
  Message(" Cannot open iternation values file...creating one!\n");
  if((fp1=fopen("iternation values","w"))==NULL)
  {
   Message(" Cannot create iternation values file...dying!\n");
  }
  else
  {
   fprintf(fp1,"%e %e %e %e",aa[0],aa[1],aa[2],aa[3]);
   fclose(fp1);
   fp1=fopen("iternation values","r");
  }
}
#endif
/*get irrelavite value for iternation*/
fscanf(fp1,"%e %e %e %e",&disp_old,&vel_old,&acel_old,&force_old);
fclose(fp1);
/*define the form of the output*/
if(acel_old==0.000000e+000)
{
  out=fopen("E:\\CFD\\my procedure\\v=1.50 ξ=0.0002\\result.txt","w");
  fprintf(out,"%10.8s  %10.5s  %10.5s  %10.5s  %10.5s\n","time","force","disp","vel","acel");
  fclose(out);
}

/*reset velocites*/
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);
if(!Data_Valid_P())
  return;
/*get the thread pointer for which this motion is defined*/
t = DT_THREAD (dt);
x_cg = DT_CG(dt)[1];
/*computer presure force on body by looping through all faces*/
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)
    beta=0.25;
gama=0.50;
m=0.000562*10;
k=1.0*10;
damp_ratio=0.0002;
wn=sqrt(k/m);
wd=2*damp_ratio*wn;
delt=dtime;
c=2*damp_ratio*wn*m;
/*calculate the integration constants*/
b0=1.0/(beta*delt*delt);
b1=gama/(beta*delt);
b2=1.0/(beta*delt);
b3=1.0/(2.0*beta)-1.0;
b4=gama/beta-1.0;
b5=delt*0.5*(gama/beta-2.0);
b6=delt*(1.0-gama);
b7=delt*gama;
/* acel_old=(force-c*vel_old-k*disp_old)/m; */
kkk=k+b0*m+b1*c;
P=force+(b0*disp_old+b2*vel_old+b3*acel_old)*m+(b1*disp_old+b4*vel_old+b5*acel_old)*c;
disp_prev=P/kkk;
/*Message ("b0 = %e, b2 = %e, b3 = %e, b6 = %e, b7 = %e\n",b0,b2,b3,b6,b7);*/
/*Message ("disp_prev = %e, disp_old = %e, vel_old = %e, acel_old = %e\n",disp_prev,disp_old,vel_old,acel_old);*/
acel_prev=b0*(disp_prev-disp_old)-b2*vel_old-b3*acel_old;
vel_prev=vel_old+b6*acel_old+b7*acel_prev;
/*Calculate the new CG position*/
x_cg=x_cg+vel_prev*dtime;
vel[1]=vel_prev;
Message ("time = %e, force = %e, disp = %e, vel = %e, acel = %e\n",time, force, disp_prev, vel_prev, acel_prev);
if ((fp1=fopen("iternation values","w"))==NULL)
{
  Message("cannot open iternation values file\n");
}
fprintf(fp1,"%e %e %e %e",disp_prev,vel_prev,acel_prev,force);
fclose(fp1);
out=fopen("E:\\CFD\\my procedure\\v=1.50 ξ=0.0002\\result.txt","a+");
fprintf(out,"%f, %e, %e, %e, %e\n",time, force, disp_prev, vel_prev, acel_prev);
fclose(out);
}
DEFINE_CG_MOTION(fluid_udf,dt,vel,omega,time,dtime)
{
vel[1]=vel_prev;
}

xinxinxiangsun

自己先顶一下，我看论坛上很多人在探讨怎么做涡致振动的问题，如果能够把这个UDF的问题解决好，相信会有很大的帮助的，高手帮助看一下

xinxinxiangsun

UDF采用的是纽马克方法进行的求解

wangsheng37

你在编写的时候考虑过无因次量没有？比如说位移无因次量和真实位移有个比例关系，如果你用无因次量代替真实值相当于将真实值扩大了好多，便会出现负体积。由于你没有给出这些物理量的定义，我也不好判断是不是这样。我目前也在做这个方向，希望能和您一起交流。我的qq是;1047858311