possible bugs in template<bool psingle> void JSphCpu::InteractionForcesDEM

ngkhaiching2000

HI Dr Canella and team

I am currently looking into your DEM codes in DualSPHysics 4.4,

I believe there could be a possible bug in the tangential force value?

There should be a "negative" signed in front of 'kn', as the friction force should be opposing to the tangential vector <tx, ty, tz>?

float ft_last=2(kn*float(DemDtForce)-gn)*vt/7;

Hopefully I can get your confirmation on the above?

In the following, I append your way on defining the normal and tangential vectors:

-------------------------------------------------------------------------------------------------------------------------

       const float drx=(psingle? psposp1.x-pspos[p2].x: float(posp1.x-pos[p2].x));

       const float dry=(psingle? psposp1.y-pspos[p2].y: float(posp1.y-pos[p2].y));

       const float drz=(psingle? psposp1.z-pspos[p2].z: float(posp1.z-pos[p2].z));

       const float rr2=drx*drx+dry*dry+drz*drz;

       const float rad=sqrt(rr2);

       const float dvx=velrhop[p1].x-velrhop[p2].x, dvy=velrhop[p1].y-velrhop[p2].y, dvz=velrhop[p1].z-velrhop[p2].z; //vji

       const float nx=drx/rad, ny=dry/rad, nz=drz/rad; //normal_ji        

       const float vn=dvxnx+dvyny+dvz*nz; //vji.nji

        //-Tangential.

        float dvxt=dvx-vn*nx, dvyt=dvy-vn*ny, dvzt=dvz-vn*nz; //Vji_t

        float vt=sqrt(dvxt*dvxt + dvyt*dvyt + dvzt*dvzt);

        float tx=0, ty=0, tz=0; //-Tang vel unit vector.

        if(vt!=0){ tx=dvxt/vt; ty=dvyt/vt; tz=dvzt/vt; }

        float ft_last=2(kn*float(DemDtForce)-gn)*vt/7; //-Elastic frictional string --> ft_elast=2*(kn*fdispl-gn*vt)/7; fdispl=dtforce*vt;