Add files via upload

src/DEM/ACM_ContactForce_inc.f90

@@ -0,0 +1,131 @@
+    Vel_w=((ri*Rvei+rj*Rvej) .cross. Norm_v)
+    Vrij = Veli-Velj + Vel_w ! Relative velocity
+    vrn  = Vrij .dot. Norm_v
+    if(VelRel_Init(ind)<-1.0_RK)VelRel_Init(ind)=abs(vrn)
+
+    Vij_n= vrn*Norm_v        ! Normal  relative velocity
+    Vij_t= Vrij-Vij_n        ! Tangent relative velocity
+
+    ! Tangential overlap Vector
+    Ovlp_t= TanDelta(ind)
+    normTan1=norm(Ovlp_t)
+    Ovlp_t= Ovlp_t-(Ovlp_t .dot. Norm_v)*Norm_v
+    normTan2=norm(Ovlp_t)
+    if(normTan2>1.0E-20_RK) then
+      Ovlp_t=(normTan1/normTan2)*Ovlp_t
+    else
+      Ovlp_t=zero_r3
+    endif
+    Ovlp_t= (Vij_t*DEM_opt%dt)+Ovlp_t
+
+    ! Computing the normal and tangential contact forces
+    TCollision=Ndt_coll*dtMax
+    select case(DEM_opt%CF_Type)
+    case(ACM_LSD)
+      if(VelRel_Init(ind)>Prop_ij%Vel_Crit) then
+        GPrtcl_HighSt(pid)="Y"
+#ifdef ContactForce_PP
+        GPrtcl_HighSt(pjd)="Y"
+#endif
+        k_n= -Prop_ij%StiffnessCoe_n/(TCollision*TCollision)
+      else
+        k_n=  Prop_ij%StiffnessCoe_n/(TCollision*TCollision)
+        k_n= -max(k_n, Prop_ij%Kn_Grav)
+      endif
+      d_n= -Prop_ij%DampingCoe_n/TCollision
+      fn = k_n*ovrlp +d_n* vrn
+      k_t= -Prop_ij%StiffnessCoe_t/(TCollision*TCollision)
+      d_t= -Prop_ij%DampingCoe_t/TCollision
+      Ftij= k_t*Ovlp_t +d_t*Vij_t
+    case(ACM_nLin)
+      if(VelRel_Init(ind)>Prop_ij%Vel_Crit) then
+        GPrtcl_HighSt(pid)="Y"
+#ifdef ContactForce_PP
+        GPrtcl_HighSt(pjd)="Y"
+#endif
+        k_n= -Prop_ij%StiffnessCoe_n/sqrt(VelRel_Init(ind)*TCollision**5)
+      else
+        k_n=  Prop_ij%StiffnessCoe_n/sqrt(Prop_ij%Vel_Crit*TCollision**5)
+        k_n= -max(k_n,Prop_ij%Kn_Grav)
+      endif
+      d_n= -Prop_ij%DampingCoe_n/TCollision
+      fn = k_n*ovrlp**1.5_RK +d_n* vrn
+      k_t= -Prop_ij%StiffnessCoe_t/(TCollision*TCollision)
+      d_t= -Prop_ij%DampingCoe_t/TCollision
+      Ftij= k_t*Ovlp_t +d_t*Vij_t
+    case(DEM_LSD)
+      if(VelRel_Init(ind)>Prop_ij%Vel_Crit) then
+        GPrtcl_HighSt(pid)="Y"
+#ifdef ContactForce_PP
+        GPrtcl_HighSt(pjd)="Y"
+#endif
+      endif
+      Vel_in=max(VelRel_Init(ind),1.0E-5_RK); Vel_in=Vel_in**0.2_RK
+      k_n= -Prop_ij%StiffnessCoe_n*Vel_in*Vel_in
+      d_n= -Prop_ij%DampingCoe_n*Vel_in
+      fn =  k_n*ovrlp +d_n* vrn
+      k_t= -Prop_ij%StiffnessCoe_t*Vel_in*Vel_in
+      d_t= -Prop_ij%DampingCoe_t*Vel_in
+      Ftij= k_t*Ovlp_t +d_t*Vij_t
+    case(DEM_nLin)
+      if(VelRel_Init(ind)>Prop_ij%Vel_Crit) then
+        GPrtcl_HighSt(pid)="Y"
+#ifdef ContactForce_PP
+        GPrtcl_HighSt(pjd)="Y"
+#endif
+      endif
+      k_n= -Prop_ij%StiffnessCoe_n
+      d_n= -Prop_ij%DampingCoe_n
+      fn =  k_n*ovrlp**1.5_RK +d_n*(ovrlp**0.25_RK)*vrn  ! 2.62, P39
+      k_t= -Prop_ij%StiffnessCoe_t*sqrt(ovrlp)
+      d_t=  0.0_RK  ! No equation is considered for tangential damping yet
+      Ftij= k_t*Ovlp_t +d_t*Vij_t                        ! 272, P44
+    end select
+    Fnij= fn*Norm_v
+
+    ! Coulomb's friction law
+    ft = norm(Ftij)
+    ft_fric = Prop_ij%FrictionCoe_s* abs(fn)
+    if(ft>ft_fric) then
+      ft_fric = Prop_ij%FrictionCoe_k* abs(fn)
+      Ftij =(ft_fric/ft)*Ftij
+      Ovlp_t =(1.0_RK/k_t)*Ftij           ! Caution Here !
+    endif
+
+    ! Computing rolling resistance torque acting on spheres
+    !W_hat = Rvei-Rvej
+    !w_hat_mag= norm(W_hat) ! 2.123, p56
+    !if(w_hat_mag > 1.0E-10_RK) then
+    !  W_hat= W_hat/w_hat_mag
+    !else
+    !  W_hat= zero_r3
+    !endif
+    !if(DEM_opt%CT_Model == CTM_ConstantTorque) then 
+    !  Mrij=(-Prop_ij%FrictionCoe_Roll*abs(fn)*Prop_ij%RadEff)*W_hat  ! 2.122, p56
+    !else
+    !  Mrij=(-Prop_ij%FrictionCoe_Roll*abs(fn)*Prop_ij%RadEff*norm(Vel_w))*W_hat  ! 2.124, p57
+    !endif
+
+    ! Setting the updated contact info pair into the contact list 
+    TanDelta(ind) = Ovlp_t
+
+    ! Updating the contact force and torques of particles i and j
+    Moment= Norm_v .cross. Ftij
+#ifdef ContactForce_PP
+    GPrtcl_cntctForce(pid)= GPrtcl_cntctForce(pid) +(Fnij+Ftij)
+    GPrtcl_torque(pid)= GPrtcl_torque(pid) +ri*Moment!+Mrij
+    GPrtcl_cntctForce(pjd)= GPrtcl_cntctForce(pjd) -(Fnij+Ftij)
+    GPrtcl_torque(pjd)= GPrtcl_torque(pjd) +rj*Moment!-Mrij
+#endif
+#ifdef ContactForce_PPG
+    GPrtcl_cntctForce(pid)= GPrtcl_cntctForce(pid) +(Fnij+Ftij)
+    GPrtcl_torque(pid)= GPrtcl_torque(pid) +ri*Moment!+Mrij
+#endif
+#ifdef ContactForce_PGP
+    GPrtcl_cntctForce(pid)= GPrtcl_cntctForce(pid) -(Fnij+Ftij)
+    GPrtcl_torque(pid)= GPrtcl_torque(pid) +rj*Moment!-Mrij
+#endif
+#ifdef ContactForce_PPFix_W
+    GPrtcl_cntctForce(pid)= GPrtcl_cntctForce(pid) +(Fnij+Ftij)
+    GPrtcl_torque(pid)= GPrtcl_torque(pid) +ri*Moment!+Mrij
+#endif

src/DEM/ACM_Integration.f90

@@ -0,0 +1,94 @@
+module Prtcl_Integration
+  use m_TypeDef
+  use Prtcl_Property
+  use Prtcl_Variables
+  use Prtcl_Parameters
+  use m_Parameters,only:gravity,PrGradData,IsUxConst
+  implicit none
+  private
+  real(RK),parameter,dimension(2):: AB2C = [1.5_RK,-0.5_RK]
+  real(RK),parameter,dimension(3):: AB3C = [23.0_RK,-16.0_RK,5.0_RK]/12.0_RK
+
+  public::Prtcl_Integrate
+#ifdef ObliqueWallTest
+  logical::IsRotate=.false.
+#endif
+contains
+
+  !******************************************************************
+  ! Prtcl_Integrate
+  !******************************************************************
+  subroutine Prtcl_Integrate(iCountACM)
+    implicit none
+    integer,intent(in)::iCountACM
+
+    ! locals
+    integer::pid,itype
+    type(real3)::linVel1,rotVel1,FpLinAcc,FpRotAcc,GravityToTal
+    real(RK)::dt,dth,MassEff,InertiaEff,MassInFluid,TimeIntCoe(2)
+#ifdef ObliqueWallTest
+    real(RK)::rMagnitude,rxDir,ryDir
+#endif
+
+    if(iCountACM==1) then
+      TimeIntCoe(1)=1.0_RK
+      TimeIntCoe(2)=0.0_RK
+    else
+      TimeIntCoe(1)=AB2C(1)
+      TimeIntCoe(2)=AB2C(2)
+    endif
+    GravityToTal=DEM_Opt%Gravity
+    if(IsAddFluidPressureGradient .and. IsUxConst) GravityToTal%x= PrGradData(2)
+
+    dt=DEM_opt%dt;   dth=dt*0.5_RK
+    DO pid =1,GPrtcl_list%nlocal
+      itype  = GPrtcl_pType(pid)
+      InertiaEff = 1.0_RK/PrtclIBMProp(itype)%InertiaEff
+      MassInFluid= PrtclIBMProp(itype)%MassinFluid
+
+      ! linear velocity and position
+      linVel1= GPrtcl_linVel(1,pid)
+      if(GPrtcl_HighSt(pid)=="Y") then   ! Turn off fluid forces for large St collisions
+        MassEff = 1.0_RK/DEMProperty%Prtcl_PureProp(itype)%Mass
+        FpLinAcc= (MassInFluid*MassEff)*GravityToTal
+      else
+        MassEff = 1.0_RK/PrtclIBMProp(itype)%MassEff
+        FpLinAcc= MassEff*GPrtcl_FpForce(pid)+ (MassInFluid*MassEff)*GravityToTal
+      endif
+      GPrtcl_linAcc(1,pid)= MassEff*GPrtcl_cntctForce(pid)
+      GPrtcl_linVel(1,pid)= linVel1 + dt*(FpLinAcc+ TimeIntCoe(1)*GPrtcl_linAcc(1,pid)+ TimeIntCoe(2)*GPrtcl_linAcc(2,pid))
+      GPrtcl_linVel(2,pid)= linVel1
+      GPrtcl_linAcc(2,pid)= GPrtcl_linAcc(1,pid)
+
+      ! rotate position
+      rotVel1= GPrtcl_rotVel(1,pid)
+      GPrtcl_rotAcc(1,pid) = InertiaEff*GPrtcl_torque(pid)
+      FpRotAcc = InertiaEff*GPrtcl_FpTorque(pid)
+      GPrtcl_rotVel(1,pid)=rotVel1 + dt*(FpRotAcc+ TimeIntCoe(1)*GPrtcl_rotAcc(1,pid)+ TimeIntCoe(2)*GPrtcl_rotAcc(2,pid))
+      GPrtcl_rotVel(2,pid)=rotVel1
+      GPrtcl_rotAcc(2,pid)=GPrtcl_rotAcc(1,pid)
+
+#ifdef RotateOnly
+      GPrtcl_linVel(1,pid)=zero_r3
+#endif
+#ifdef ObliqueWallTest
+      GPrtcl_linVel(1,pid)%z=0.0_RK
+      if(.not.IsRotate) then
+        GPrtcl_rotVel(1,pid)=zero_r3
+        rMagnitude= 1.0_RK/sqrt(DEM_Opt%Gravity%x*DEM_Opt%Gravity%x+DEM_Opt%Gravity%y*DEM_Opt%Gravity%y)
+        rxDir= rMagnitude*DEM_Opt%Gravity%x
+        ryDir= rMagnitude*DEM_Opt%Gravity%y
+        rMagnitude=sqrt(GPrtcl_linVel(1,pid)%x*GPrtcl_linVel(1,pid)%x+GPrtcl_linVel(1,pid)%y*GPrtcl_linVel(1,pid)%y)
+        GPrtcl_linVel(1,pid)%x= rMagnitude*rxDir
+        GPrtcl_linVel(1,pid)%y= rMagnitude*ryDir
+      endif
+#endif
+      GPrtcl_PosR(pid)=GPrtcl_PosR(pid)+dth*(linVel1+ GPrtcl_linVel(1,pid))
+#ifdef ObliqueWallTest    
+      if(GPrtcl_PosR(pid)%y<2.0_RK*GPrtcl_PosR(pid)%w) IsRotate=.true.
+#endif
+    ENDDO
+
+  end subroutine Prtcl_Integrate
+
+end module Prtcl_Integration