模拟阀底弹簧的一个UDF，还没调试，有几个有问题想请教

求是创新

以下就是模拟阀弹簧的UDF，感觉有几个地方肯定有问题，其中一个就是我想定义整个下底面为动网格，这样的话用 F_CENTROID 宏找Y坐标轴的时候我不确定找的是阀底面的哪个面（我本意是想找最中心的那个面的Y轴值以便于求解弹簧伸缩量），阀底模型图请见附件，另外请看一下还有其他问题么，如real数据格式的话小数点后取6位是否可行？

#include "udf.h"
static real v_prev = 0.0;
static real vector[3] = {0,1,0};
static real f = 120.211498;
static real k = 53.0022;
DEFINE_CG_MOTION(spring,dt,vel,omega,time,dtime)
{
    Thread *t;
    face_t f;
    real NV_VEC(A);
    real force, dv;
    real x[ND_ND];
    real y = 0.0;
    /* 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);
    /* compute pressure 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) * NV_DOT(A,vector);
      }
    end_f_loop(f,t)
    /* compute change in velocity, that is, dv = F * dt / mass  velocity update using explicit Euler formula */

    F_CENTROID(x,f,t);
    y = 0.0161 - x[1];
    dv = dtime * ( force + y * k ) / 0.1;
    v_prev += dv;
    Message ("time = %f, x_vel = %f, force = %f\n", time, v_prev,
      force);
    /* set x-component of velocity */
    vel[1] = v_prev;
}

[ 本帖最后由 求是创新 于 2014-3-7 10:08 编辑 ]