|
|
for(m=0;m<save;m++)
{
tnow+=dts;//当前时间值
//-----------------接触力和矩初始化--------------------------------------------
for(i=0;i<nu;i++){fxp=fyp=fzp=0;mx=my=mz=0;}
//-----------------计算空隙率--------------------------------------------------
for(i=0;i<l1;i++)for(j=0;j<m1;j++)for(k=0;k<n1;k++){voi[j][k]=1.0;cell[j][k]=0;}
for(i=0;i<nu;i++)
{
printf("%d.",m);
cx=(int)(px/dx)+1;cz=(int)(pz/dz)+1;cy=1;//cy=(int)(py/dy)+1;//当前颗粒所在网格
printf("%d,%d,%d,%f,%f,%f\n",cx,cy,cz,px,py,pz);
voi[cx][cy][cz]-=(4*pi*r[pk]*r[pk]*r[pk]/3/vol);//pk:i颗粒的种类
ini=cell[cx][cy][cz];cpn[cx][cy][cz][ini]=i;cell[cx][cy][cz]++;//cell[][][]该网格内的颗粒数 cpn[][][][]第i个颗粒在该网格的序号为ini
printf("ini=%d",ini);
if(cell[cx][cy][cz]==(1+CM))
{
printf("cell_n is too full!\n该网格内颗粒数=%d\tcx=%d\tcy=%d\tcz=%d\n",cell[cx][cy][cz],cx,cy,cz);
exit(0);
}
}
// for(i=8;i<13;i++)for(j=1;j<3;j++)for(k=8;k<15;k++)voi[j][k]=(float)0.3;
//-----------------碰撞计算----------------------------------------------------
ac();
//-----------------计算流体力--------------------------------------------------
if(tnow>tt)//tt为投颗粒过程,开始进风时间
{
fluidf();
for(i=0;i<nu;i++)
{
fxp+=fxf;
fyp+=fyf;
fzp+=fzf;
}
}
} |
|