我把我的程序粘上，大家看看出再错误的原因，帮帮忙

lantian0606

#include "udf.h"
#include "mem.h"
#include "dynamesh_tools.h"
#include "math.h"
static real v_x = 0.0,v_y=0.0,omega_z=0.0;
DEFINE_CG_MOTION(piston,dt,vel,omega,time,dtime)
{
Thread *t,*tc;
face_t f;
cell_t c;
real NV_VEC(A) ,NV_VEC(CG);
real force_x,force_y,pressure_x,pressure_y,viscosity_x,viscosity_y,moment_z,dv_x,dv_y,domega_z;
real xf[ND_ND],xc[ND_ND];
real dir_differential[3]={0,0,0};
real viscosity_coff,viscosity_stress[3],viscosity_part[3],viscosity_x_part,viscosity_y_part;
real moment_p_z,moment_v_z;
int i;

/* reset velocities */
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);
CG[0]=DT_CG(dt)[0],CG[1]=DT_CG(dt)[1],CG[2]=DT_CG(dt)[2];
if (time<0.00505&&time>=0)
  { vel[0]=0,vel[1]=34.5,omega[2]=0;
    return;
  }
/* compute pressure force , viscosity force and moment on body by looping through all faces */
pressure_x = 0.0,pressure_y=0.0;
moment_p_z=0.0;
viscosity_x=0.0,viscosity_y=0.0;
moment_v_z=0.0;
begin_f_loop(f,t)
{
  F_AREA(A,f,t);
  F_CENTROID(xf,f,t);
  
  c=F_C0(f,t);
  tc=THREAD_T0(t);
  C_CENTROID(xc,c,tc);
  
/* pressure force*/
  pressure_x += -F_P(f,t) * A[0];
  pressure_y += -F_P(f,t) * A[1];

/* viscosity force*/
  
  /*U making*/
  for(i=0;i<3;i++)
    dir_differential[0]+=C_U_G(c,tc)*A/NV_MAG(A)*sqrt(1-A[0]*A[0]/(NV_MAG(A)*NV_MAG(A)));
  /*V making*/
  for(i=0;i<3;i++)
    dir_differential[1]+=C_V_G(c,tc)*A/NV_MAG(A)*sqrt(1-A[1]*A[1]/(NV_MAG(A)*NV_MAG(A)));
  /*W making*/
  for(i=0;i<3;i++)
    dir_differential[2]+=C_W_G(c,tc)*A/NV_MAG(A)*sqrt(1-A[2]*A[2]/(NV_MAG(A)*NV_MAG(A)));
  
  /*compute viscous stress and viscous force of cell faces*/
  viscosity_coff=C_NUT(c,tc);
  for(i=0;i<=3;i++)
    {
     viscosity_stress=dir_differential*viscosity_coff;
     viscosity_part=viscosity_stress*NV_MAG(A);
    }
  /*compute viscous force*/
    viscosity_x_part=viscosity_part[0]*sqrt(1-A[0]*A[0]/(NV_MAG(A)*NV_MAG(A)))-viscosity_part[1]*A[1]/NV_MAG(A)*sqrt(A[0]*A[0]/(A[0]*A[0]+A[2]*A[2]))-viscosity_part[2]*A[2]/NV_MAG(A)*sqrt(A[0]*A[0]/(A[1]*A[1]+A[0]*A[0]));
   
    viscosity_y_part=viscosity_part[1]*sqrt(1-A[1]*A[1]/(NV_MAG(A)*NV_MAG(A)))-viscosity_part[0]*A[0]/NV_MAG(A)*sqrt(A[1]*A[1]/(A[1]*A[1]+A[2]*A[2]))-viscosity_part[2]*A[2]/NV_MAG(A)*sqrt(A[1]*A[1]/(A[1]*A[1]+A[0]*A[0]));
   
    viscosity_x+=viscosity_x_part,viscosity_y+=viscosity_y_part;
  
  /*moment*/
  
   /* pressure force making*/
   moment_p_z+=F_P(f,t) * A[0]*(xf[1]-CG[1])-F_P(f,t) * A[1]*(xf[0]-CG[0]);
   
   /* viscosity force making*/
   moment_v_z+=-viscosity_x_part*(xc[1]-CG[1])+viscosity_y_part*(xc[0]-CG[0]);
}
end_f_loop(f,t)
   
force_x=pressure_x+viscosity_x;
force_y=pressure_y+viscosity_y;
moment_z=moment_p_z+moment_v_z;
/* compute change in velocity and angle velocity, i.e., dv = F * dt / mass
velocity update using explicit Euler formula */
dv_x = dtime * force_x / 39400;
v_x += dv_x;
dv_y = dtime * force_y / 39400;
v_y += dv_y;
domega_z=dtime*moment_z/390000;
omega_z=domega_z;

Message ("\ntime=%f ,x=%f ,y=%f , force_x=%.1f ,force_y=%.1f ,moment_z=%.1f ,Vx=%f ,Vy=%f ,omega_z=%f\n", time,CG[0],CG[1],force_x,force_y,moment_z,v_x,v_y,omega_z);
/* set components of velocity */
vel[0] = v_x,vel[1]=v_y,omega[2]=omega_z;
}
时间步取0.001秒，当运行到0.005秒以后就会出再如下错误：
FLUENT received fatal signal (ACCESS_VIOLATION)
1. Note exact events leading to error.
2. Save case/data under new name.
3. Exit program and restart to continue.
4. Report error to your distributor.
Error Object: ()
请问这是什么错误，原因是什么．
注：本程序想解决的问题是：流场中的一个任意形状物体（对称或非对称），它会受到x向和y向的受力（本程序不考虑Ｚ向）以及绕其质心的力矩，这两个方向的力和此力矩，使它产生x向和y向的速度的变化以及绕其质心角速度的变化，从而使物体运动发生变化。本例采用动网格去做。
大家，高人，快帮帮忙吧。
版主，发这个能给加威望吗，呵呵

lantian0606

版主是个大猛人吧，来，帮我想想，谢谢[br][br][以下内容由 lantian0606 在 2006年07月30日 11:42am 时添加] [br]
快来人帮帮忙啊,急
版主,发动大伙啊

lantian0606

急盼中
另问,如何求压力和粘性力,要程序!!

wangyunfeng

程序很有意思，不过要看懂的话是要有这个方面的经验，，很感谢搂住提供此程序

wangyunfeng

第十一个例子valve的动网格定义中好像有求力的吧，，相信搂住看过的

chbo76

楼主，请将里面涉及到的计算公式注释一下，这样才能帮你分析程序找错误啊。

lantian0606

上面的式主要是计算粘性力,压力和绕质心的力矩(压力和粘性力产生的)的公式,粘性力=粘性系数*垂直于速度方向上的速度梯度.具体公式,我推导过,很多,不知道 对不对.待我整理好,发到网上,望大伙帮忙.我在上面的程序中有简单的注释,大家先可认为我求力的方法是正确的,看看程序有无其它问题.待我把粘性力计算公式发到网上,咱再细研究.本程序能够通过编译,只是计算时出错

lantian0606

我多次试验,发现问题可能在
c=F_C0(f,t);
tc=THREAD_T0(t);
这块地方,其中,c为cell_t型变量,而tc为Thread *型变量
我用这两个语句,是为了找到与运动物体表面紧挨着的单元,来计算速度梯度,和粘性系数.
我只是看着手册写的,对这两句理解不深.
望高手们指点

wangyunfeng

DEFINE_UDS_FLUX(myudsflux,f,t,i)
{
cell_t c0, c1 = -1;
Thread *t0, *t1 = NULL;
real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0;
c0 = F_C0(f,t);
t0 = F_C0_THREAD(f,t);
F_AREA(A, f, t);
if (BOUNDARY_FACE_THREAD_P(t))
{
real dens;
if (NNULLP(THREAD_STORAGE(t,SV_DENSITY)))
dens = F_R(f,t);
else
dens = C_R(c0,t0);
NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, dens);
flux = NV_DOT(psi_vec, A);
}
else
{
c1 = F_C1(f,t);
t1 = F_C1_THREAD(f,t);
NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,C_R(c0,t0));
NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,C_R(c1,t1));
flux = NV_DOT(psi_vec, A)/2.0;
}
return flux;
}
[br][br][以下内容由 wangyunfeng 在 2006年08月03日 09:30am 时添加] [br]
上面的例子是帮助文件中的，，其中就有定义域内部面的uds flux的方法，，就是else后面的，用的就是你说的面的左右临近的单元面的量

lantian0606

&#35;include "udf.h"
&#35;include "mem.h"
&#35;include "dynamesh_tools.h"
static real v_x = 0.0,v_y=0.0,omega_z=0.0;
DEFINE_CG_MOTION(piston,dt,vel,omega,time,dtime)
{
Thread *t,*tc;
face_t f;
cell_t c;
real NV_VEC(A), NV_VEC(CG);
real force_x,force_y,pressure_x,pressure_y,viscosity_x,viscosity_y,viscosity_z,moment_z,dv_x,dv_y,domega_z;
real xf[ND_ND],xc[ND_ND];
real dir_differential[3]={0,0,0};
real viscosity_coff,viscosity_stress[3],viscosity_part[3],viscosity_x_part,viscosity_y_part;
real moment_p_z,moment_v_z;
int i;

/* reset velocities */
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
if (!Data_Valid_P())
return;
if(time<=0.001)
{
  vel[1]=34.5;
  return;
}
/* get the thread pointer for which this motion is defined */
t = DT_THREAD(dt);
CG[0]=DT_CG(dt)[0],CG[1]=DT_CG(dt)[1],CG[2]=DT_CG(dt)[2];
/* compute pressure force , viscosity force and moment on body by looping through all faces */
pressure_x = 0.0,pressure_y=0.0;
moment_p_z=0.0;
viscosity_x=0.0,viscosity_y=0.0;
moment_v_z=0.0;

begin_f_loop(f,t)
{
  F_AREA(A,f,t);
  F_CENTROID(xf,f,t);
  
  c=F_C0(f,t);
  tc=F_C0_THREAD(f,t);
  
  C_CENTROID(xc,c,tc);
  
/* pressure force*/
  pressure_x += -F_P(f,t) * A[0];
  pressure_y += -F_P(f,t) * A[1];
/* viscosity force*/
  
  /*U making*/
  for(i=0;i<3;i++)
    dir_differential[0]+=C_U_G(c,tc)*A/NV_MAG(A)*sqrt(1-A[0]*A[0]/(NV_MAG(A)*NV_MAG(A)));
  /*V making*/
  for(i=0;i<3;i++)
    dir_differential[1]+=C_V_G(c,tc)*A/NV_MAG(A)*sqrt(1-A[1]*A[1]/(NV_MAG(A)*NV_MAG(A)));
  /*W making*/
  for(i=0;i<3;i++)
    dir_differential[2]+=C_W_G(c,tc)*A/NV_MAG(A)*sqrt(1-A[2]*A[2]/(NV_MAG(A)*NV_MAG(A)));
  
  /*compute viscous stress and viscous force of cell faces*/
  viscosity_coff=C_NUT(c,tc);
  for(i=0;i<=3;i++)
    {
     viscosity_stress=dir_differential*viscosity_coff;
     viscosity_part=viscosity_stress*NV_MAG(A);
    }
  /*compute viscous force*/
    viscosity_x_part=viscosity_part[0]*sqrt(1-A[0]*A[0]/(NV_MAG(A)*NV_MAG(A)))-viscosity_part[1]*A[1]/NV_MAG(A)*sqrt(A[0]*A[0]/(A[0]*A[0]+A[2]*A[2]))-viscosity_part[2]*A[2]/NV_MAG(A)*sqrt(A[0]*A[0]/(A[1]*A[1]+A[0]*A[0]));
   
    viscosity_y_part=viscosity_part[1]*sqrt(1-A[1]*A[1]/(NV_MAG(A)*NV_MAG(A)))-viscosity_part[0]*A[0]/NV_MAG(A)*sqrt(A[1]*A[1]/(A[1]*A[1]+A[2]*A[2]))-viscosity_part[2]*A[2]/NV_MAG(A)*sqrt(A[1]*A[1]/(A[1]*A[1]+A[0]*A[0]));
   
    viscosity_x+=viscosity_x_part,viscosity_y+=viscosity_y_part;
  
  /*moment*/
  
   /* pressure force making*/
   moment_p_z+=F_P(f,t) * A[0]*(xf[1]-CG[1])-F_P(f,t) * A[1]*(xf[0]-CG[0]);
   
   /* viscosity force making*/
   moment_v_z+=-viscosity_x_part*(xc[1]-CG[1])+viscosity_y_part*(xc[0]-CG[0]);
}
end_f_loop(f,t)
   
force_x=pressure_x+viscosity_x;
force_y=pressure_y+viscosity_y;
moment_z=moment_p_z+moment_v_z;
/* compute change in velocity and angle velocity, i.e., dv = F * dt / mass
velocity update using explicit Euler formula */
dv_x = dtime * force_x / 39400;
v_x += dv_x;
dv_y = dtime * force_y / 39400;
v_y += dv_y;
domega_z=dtime*moment_z/390000;
omega_z=domega_z;

Message ("\ntime=%f ,x=%f ,y=%f , force_x=%.1f ,force_y=%.1f ,moment_z=%.1f ,Vx=%f ,Vy=%f ,omega_z=%f\n", time,CG[0],CG[1],force_x,force_y,moment_z,v_x,v_y,omega_z);
/* set components of velocity */
vel[0] = v_x,vel[1]=v_y,omega[2]=omega_z;
}
我发现了一些问题,已将
  c=F_C0(f,t);
tc=THREAD_T0(t);
两句改成
  c=F_C0(f,t);
  tc=F_C0_THREAD(f,t);(关于这两句,你看我上传的附件)
经多次调试,发现
程序在计算速度梯度时出错,即在有
C_U_G(c,tc)
C_V_G(c,tc)
C_V_G(c,tc)
的语句的地方可能有问题
你现看看

lantian0606