求指点：dpm模型中加入与颗粒位置相关的体积力

madongdong

我的udf程序加载进去，每次都提示这个错误FLUENT received fatal signal (ACCESS_VIOLATION)；不知道问题出在哪儿了？能否帮我看看程序哪里有问题，万分感谢！
#include "udf.h"
#include "dpm.h"
#include "surf.h"
#define a 1
#define b 7
#define U 3000
#define D 0.2
#define R1 0.00015
#define T0 2000   
#define TSTART 0.0

float q_function(float M)
{
  float y,k;
  k=(b-a)/(b+2*a);
  y=((12*3.14*a*b*R1*R1*U)/((a+2*b)*D*(1-k*M)));
  return(y);
}

float E_function(float M)
{
  float y,k;
  k=(b-a)/(b+2*a);
  y=(3*b*U)/((b+2*a)*(1-k*M)*D);
  return(y);
}


DEFINE_DPM_BODY_FORCE(p_b_force,p,i)
  {
    cell_t c=RP_CELL(&p->cCell);
    Thread *t=RP_THREAD(&p->cCell);
    int j=0,num0=0,num1=0;
    float pos[4]={0};
    Particle *pi;
    float bforce;
    float q,E,s,M,bforce_x=0,bforce_y0=0,bforce_y=0,bforce_z=0;


    begin_particle_cell_loop(pi,c,t)
    {
      num0++;
    }
    end_particle_cell_loop(pi,c,t)

    if(P_TIME(p)>TSTART)
    {
      if(num0>0)
     {

        begin_particle_cell_loop(pi,c,t)
        {
          for(j=0;j<3;j++)
          {
            pos[j]=P_POS(pi)[j]-P_POS(p)[j];
            s=sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
        
                                                                                                                                                                                                                                                                                                                                                                            if(s<=3*R1)                                                               
            {
              num1=num1+1;
              M=2*R1/(s*s*s);
              q=q_function(M)*(P_TIME(p)/(P_TIME(p)+T0));
              E=E_function(M)-q/(4*3.14*a*R1*R1);
              bforce=E*q;
              bforce_y0+=bforce*pos[1]/s;
              bforce_y=bforce_y0/num1;
            }
          }
        }
        end_particle_cell_loop(pi,c,t)
      }

      if(i==0)
      {
        bforce=bforce_x;
      }
      else if(i==1)
      {
        bforce=bforce_y;
      }
    }
    else bforce=bforce_z;
    return(bforce/P_MASS(p));
  }

classic1573

请问您这个问题解决了吗？我也是同样的问题，谢谢指点