usr1_des.f

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: URS1_DES                                               !
!                                                                      !
!  Purpose: This routine is called within the discrete phase time loop !
!  after the source terms have been calculated but before they are     !
!  applied. The user may insert code in this routine or call user      !
!  defined subroutines.                                                !
!                                                                      !
!  This routine is called from the time loop, but no indicies (fluid   !
!  cell or particle) are defined.                                      !
!                                                                      !
!  Author: J.Musser                                   Date: 06-Nov-12  !
!                         


!  RAJATH KANTHARAJ, 
!  Purdue University, Spring 2020

!  Comments:
!  Stokes drag is implemented in this subroutine with the solution to 2D squeeze of Newtonian liquid
!  Top wall -- constant velocity 
!  Side wall -- NOT needed                                                             
!                                                                      
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE USR1_DES

      Use discretelement
      Use des_rxns
      Use des_thermo
      Use run
      Use indices  
      Use usr
      Use comp_mean_fields_mod, only: comp_mean_fields ! This will call calc_epg_des to update gas/solids volume fraction after reassigning radius
      Use fldvar
      Use compar  
      Use drag ! Call function to calculate Cd  
      Use functions, only: is_nonexistent, fluid_at, is_entering_ghost, is_exiting_ghost, is_ghost 
      Use exit  , only: mfix_exit 
      Use constant, only: pi
      Use vtk, only: vtk_time, vtk_dt 
      Use geometry
      Use physprop ! Fluid viscosity, MU_G is defined here  
      use vtp, only: write_vtp_file  
      Use derived_types, only: pic    
      Use param1  
  
      IMPLICIT NONE
      INTEGER :: LL, calls_write ! particle index LL, frequency flag to write .txt file      
      INTEGER, save :: flag_call = 1, flag_txt ! first time function call flag, TXT file write frequency flag  
      DOUBLE PRECISION :: OVERLAP_top ! Particle-wall overlaps with side/top walls 
      DOUBLE PRECISION :: F_part_top ! Force on particles by side/top walls
      DOUBLE PRECISION :: P_top, F_top, F_exer !TOP wall velocity, pressure and force on it exerted by particles 
      DOUBLE PRECISION :: v, top_particle, B
      INTEGER :: max_rad_loc, max_part_rad 

      ! Save top, side wall velocities, positions for future function calls to avoid corruption
      DOUBLE PRECISION, parameter :: vel_top = -10.0E-3, npart=2.0
      DOUBLE PRECISION, save :: F_supp ! Force limit (N) for
                                                ! squeezing         

      DOUBLE PRECISION, save :: x_star, pos_top, max_diap, Xf_new  
      DOUBLE PRECISION, DIMENSION(MAX_PIP,3) :: dist_ij_1 ! Center-center distance vector        
      DOUBLE PRECISION :: dist_ij_2, dist_ij ! Other variables to store distance 


      DOUBLE PRECISION :: pos_change, tol_change, tol_force ! simulation end criterion when top
                                     ! plate has reached a steady-state        
      DOUBLE PRECISION :: mass_top, A_top, F_net_top ! TOP wall mass, cross-section area, net force on it
      DOUBLE PRECISION :: vel_top_old 
      DOUBLE PRECISION :: pos_top_old 
      DOUBLE PRECISION :: KN_DES_wall, KT_DES_wall, Z_domain, X_domain, mean_dia
      DOUBLE PRECISION :: VRELx, VRELy, vflx, vfly, alph, time_drag
      DOUBLE PRECISION, parameter :: vflz = 0.0, P_supp=500.0E3 !500 kPa pressure   
      DOUBLE PRECISION :: F_fluid_top, F_ratio, P_sqz, vtop_x 
      
      DOUBLE PRECISION, DIMENSION(MAX_PIP) :: top_overlap, part_dist
      DOUBLE PRECISION, DIMENSION(3) :: F_drag
      DOUBLE PRECISION :: P_part, P_fluid, F_total, pmass_rk 
      double precision :: V_FP_mag, Re_p, beta_drag ! Relative particle-fluid velocity magnitude, Particle Reynold's number, drag coefficient
      DOUBLE PRECISION, PARAMETER :: rho_fl=1000.0, eta_fl = 20.0 ! Fluid viscosity
      DOUBLE PRECISION :: min_diam, mass_one, mass_effs, damp_cft, tcoll
      DOUBLE PRECISION :: vol_particle, vol_frac_rk, vol_sum_mfix
      DOUBLE PRECISION :: eps_new, epg_new, vol_ijk  
  
        ! Local variables
        !---------------------------------------------------------------------//
        ! general i, j, k indices
              INTEGER :: I, J, K, IJK, cur_ijk
              INTEGER :: II, JJ, KK
              INTEGER :: IPJK, IJPK, IJKP,&
                         IPJPK, IPJKP, IJPKP, IPJPKP
        ! indices used for interpolation stencil (unclear why IE, JN, KTP are
        ! needed)
              INTEGER :: IW, IE, JS, JN, KB, KTP
        ! i,j,k indices of the fluid cell the particle resides in minus 1
        ! (e.g., shifted 1 in west, south, bottom direction)
              INTEGER, DIMENSION(3) :: PCELL
        ! order of interpolation set in the call to set_interpolation_scheme
        ! unless it is re/set later through the call to
        ! set_interpolation_stencil
              INTEGER :: ONEW
        ! particle number index, used for looping
              INTEGER :: NP, nindx
        ! constant whose value depends on dimension of system
        ! avg_factor=0.125 (in 3D) or =0.25 (in 2D)
        ! avg_factor=0.250 (in 3D) or =0.50 (in 2D)
              DOUBLE PRECISION :: AVG_FACTOR
        DOUBLE PRECISION :: desposnew(3), velfp(3), vel_new(3)

   
      Z_domain = NZ2 - SZ1 ! Domain length along Z 
      X_domain = EX2 - WX1 ! Domain length along X  
      x_star = 0.5*X_domain ! Used to translate origin for fluid
                            ! velocity equation

      F_top = 0.0  
      mass_top = 1.0E-5  
      tol_change = 1.0E-10  
      tol_force = 1.0E-12  
      calls_write = 50 ! .txt file write frequency 
             
      open(UNIT=77, FILE='wall_dynamics.txt', ACCESS='append')
      open(UNIT=95, FILE='blt_plate_pos.txt', ACCESS='append')  
      open(UNIT=101, FILE='Culprit_particle.txt', ACCESS='append')
      open(UNIT=107, FILE='interparticle_overlap.dat', ACCESS='append')  
      open(UNIT=197, FILE='local_particle_vfrac.txt', ACCESS='append')  
              
      top_overlap = 0.0   
      
      IF(flag_call == 1) THEN
        F_supp = 345.0E3 * Z_domain * X_domain ! Ultimate squeeze force based on 100kPa pressure

        flag_call = 0
        flag_txt = 49 ! Write wall positions to text file every calls_write number of times 
      
        min_diam = 2.0*MINVAL(DES_RADIUS, MASK=DES_RADIUS .GT. 0) ! Minimum diameter
        write(*,*) "Minimum particle diameter: ", min_diam
        write(*,*) "Particle density: ", RO_S0(1)
        mass_one = PI/6.0 * min_diam**3.0 * RO_S0(1) ! Only one particle phase M = 1
        write(*,*) "Mass of one particle: ", mass_one
        mass_effs = mass_one * mass_one/(mass_one + mass_one)
        write(*,*) "Coefficient of restitution: ", DES_EN_INPUT(1)
        damp_cft = 2*SQRT(KN*mass_effs) * ABS(log(DES_EN_INPUT(1))) &
                   / SQRT(PI*PI + log(DES_EN_INPUT(1))**2)
        write(*,*) "Effective mass: ", mass_effs

        tcoll = PI/(SQRT(KN/mass_effs - &
                   (damp_cft/mass_effs)**2/4.0))
        write(*,*) "Collision time: ", tcoll
        DTSOLID = 1.0/50.0 * tcoll ! 3.0E-8 !1.0/50.0 * t_coll ! Modify the DEM time step based on minimum
                       ! particle diameter
        write(*,*) "Corrected time step: ", DTSOLID
  
        max_diap = MAXVAL(2*DES_RADIUS)
        write(*,*) 'Max diameter: ', max_diap
        write(*,*) "Min particle radius (microns): ", min_diam*1e6
        
        write(77,*) "Top_overlapR      Top_wall_position(mm)" //&  
                    "  positionChange  P_top_particles" //& 
                    "  F_top_part  " //&
                    " max_particle_position_top"
 
        max_rad_loc = MAXLOC(DES_POS_NEW(:,1), DIM=1) 
        alph = MAXVAL(DES_POS_NEW(:,1)) + DES_RADIUS(max_rad_loc) &
                          - x_star
        Xf_new = alph
 
        max_part_rad = MAXLOC(DES_POS_NEW(:,2), DIM=1)
        pos_top_old = MAXVAL(DES_POS_NEW(:,2)) + MAXVAL(DES_RADIUS) !DES_RADIUS(max_part_rad)

        write(*,*) "Top wall initially at ", pos_top_old/1.0E-3, " mm"
        write(*,*) "Top wall initial velocity: ", vel_top*1000, & 
                                                           " mm/s"
      ELSE
        pos_top_old = pos_top

        ! Velocity of mid-point i.e., y = b/2 and x = x_o,f
        vtop_x = -6.0*abs(vel_top)/(pos_top_old) * (-1d0/4d0 * Xf_new)
        !write(*,*) "Velocity at mid-point (mm/s): ", vtop_x*1000

        ! New farthest fluid X position
        Xf_new = Xf_new + vtop_x * DTSOLID
        alph = Xf_new

       !! ---------------- DRAG FORCE -------------------------- !!
       ! Use MFIX subroutines to compute local volume fraction !
       ! Then, use Ding and Gidaspow formulation to calculate drag force 
       ! Add drag force to interparticle contact forces  

       DO ijk = ijkstart3,ijkend3
         if(.not.fluid_at(ijk) .or. pinc(ijk).eq.0) cycle
           i = i_of(ijk)
           j = j_of(ijk)
           k = k_of(ijk)

          ! generally a particle may not exist in a ghost cell. however, if the
          ! particle is adjacent to the west, south or bottom boundary, then pcell
          ! may be assigned indices of a ghost cell which will be passed to
          ! set_interpolation_stencil
           pcell(1) = i-1
           pcell(2) = j-1
           pcell(3) = merge(1, k-1, NO_K)
          
           
           CALL COMP_MEAN_FIELDS ! Compute local solid/fluid volume fraction

           MU_G(IJK) = eta_fl     
           RO_G(IJK) = rho_fl
           
           IF(pos_top < YN(j) .AND. & 
              pos_top > YN(j)-DY(j)) THEN
              vol_ijk = DX(i)*DZ(k) * & 
                       ( pos_top - (YN(j)-DY(j)) )

             eps_new = EP_S(IJK,1)*VOL(IJK)/vol_ijk
             epg_new = ONE - eps_new 
           ELSE
             vol_ijk = VOL(IJK)   
             EPG_NEW = EP_G(IJK)        
           ENDIF

          
           vol_particle = 0.0
           vol_sum_mfix = 0.0
          
           DO nindx = 1,PINC(IJK)
              NP = PIC(ijk)%p(nindx)
              !write(*,*) "Particle number NP: ", NP
              ! skipping indices that do not represent particles and ghost particles
              if(is_nonexistent(np)) cycle
              if(is_ghost(np).or.is_entering_ghost(np).or.is_exiting_ghost(np)) cycle

              desposnew(:) = des_pos_new(np,:)
              vol_particle = vol_particle + &
                             pi/6.0*(2*DES_RADIUS(NP))**3
              vol_sum_mfix = vol_sum_mfix + PVOL(NP)      
         
             ! ENDIF  
              !call DRAG_INTERPOLATION(gst_tmp,vst_tmp,desposnew,velfp,weight_ft)

              ! Calculate the particle centered drag coefficient (F_GP) using the
              ! particle velocity and the interpolated gas velocity.  Note F_GP
              ! obtained from des_drag_gp subroutine is given as:
              !    f_gp=beta*vol_p/eps where vol_p is the particle volume.
              ! The drag force on each particle is equal to:
              !    beta(u_g-u_s)*vol_p/eps.
              ! Therefore, the drag force = f_gp*(u_g - u_s)
              VEL_NEW(:) = DES_VEL_NEW(NP,:)

              vflx = -6.0*abs(vel_top)/(pos_top_old)**3d0 * & 
                      (desposnew(1)-x_star) * (desposnew(2)**2d0 & 
                                       - pos_top_old*desposnew(2))

              vfly = 6.0*abs(vel_top)/(pos_top_old)**3d0 * & 
                         (1.0/3.0*desposnew(2)**3d0 - & 
                          1.0/2.0*pos_top_old*desposnew(2)**2d0)

              
              VELFP = (/vflx, vfly, vflz /) ! vflz=0; Z-velocity not used in drag force   
               
              CALL DES_DRAG_GP(NP, VEL_NEW, VELFP, EPG_NEW)

              F_drag(1:3) = F_GP(NP)*(VELFP-VEL_NEW)
              !! Ignore Z-component of drag force !!  
                      
              FC(NP,:3) = FC(NP,:3) + F_drag  
              !write(*,*) "Total force: ", FC(NP, :3)
              !read(*,*)    
           ENDDO
            vol_frac_rk = vol_particle/vol_ijk
            IF(MOD(IJK,2)==0 .AND. MOD(flag_txt, 200)==0) THEN    
               write(197,*) vol_particle, vol_sum_mfix, &
                            vol_ijk, & 
                            vol_frac_rk, EP_S(IJK,1)
            ENDIF
        ENDDO

      ENDIF 
        
  
        KN_DES_wall = KN
        F_top = 0.0
! SIDE, TOP WALL dynamics loop through all particles for current time step
      DO LL = 1, MAX_PIP
       IF(IS_NONEXISTENT(LL)) CYCLE 
        
         OVERLAP_top = DES_RADIUS(LL)-(pos_top_old - DES_POS_NEW(LL,2)) ! y-position
         top_overlap(LL) = OVERLAP_top/DES_RADIUS(LL)
          
        IF(OVERLAP_top > 0) THEN
          F_part_top = -KN_DES_wall*OVERLAP_top ! Negative sign indicates that force on particle acts along -Y

          ! Update particle force array
          FC(LL,2) = FC(LL,2) + F_part_top
          
          ! Update total force acting on wall due to particles
          F_top = F_top - F_part_top ! This force acts along +Y
        ENDIF 

      ENDDO

  
      flag_txt = flag_txt + 1  
      
      ! A_top = pos_side_old * Z_domain ! X*Z  
        A_top = X_domain * Z_domain

        pos_top = pos_top_old + DTSOLID*vel_top ! Wall is coming down, so pos_top is decreasing in magnitude
        pos_change = ABS(pos_top_old - pos_top)/pos_top_old

         !! ---- Calculate squeeze force, pressure ---- !!
              ! Fluid force
        !F_fluid_top = 8.0*eta_fl*abs(vel_top)/pos_top**3 &
        !                 * Z_domain * alph**3
 
        ! Particle force is F_top (see above)
        ! Total squeeze force & pressure
       ! F_total = F_top + F_fluid_top   
       ! P_sqz = F_total/(2*alph * Z_domain)  

        !P_part = F_top/(2*alph * Z_domain)
        !P_fluid = F_fluid_top/(2*alph * Z_domain)   
        
         ! FLAG to check how far down the top wall has been displaced 
           ! If top wall has come down such that there is only about 5
           ! particle diameters across, then stop squeeze simulations
        mean_dia = 2*sum(DES_RADIUS)/size(DES_RADIUS)

           !! ---- Calculate squeeze force, pressure ---- !!
              ! Fluid force
        F_fluid_top = 8.0*eta_fl*abs(vel_top)/pos_top**3 &
                         * Z_domain * alph**3
 
        ! Particle force is F_top (see above)
        ! Total squeeze force & pressure
        F_total = F_top + F_fluid_top   
        P_sqz = F_total/(2*alph * Z_domain)  

        P_part = F_top/(2*alph * Z_domain)
        P_fluid = F_fluid_top/(2*alph * Z_domain)
 
        IF(pos_top < npart*max_diap) THEN ! Pressure-controlled squeezing
           ! Exit function and complete squeeze simulation
           write(95,*) "User-specified squeeze pressure (kPa): ", & 
                                                P_supp/1.0E3
           write(95,*) "Total squeeze force (N): ", F_total
           write(95,*) "Total squeeze pressure (kPa): ", P_sqz/1.0E3
           write(95,*) "Particle pressure (kPa): ", P_part/1.0E3
           write(95,*) "Fluid pressure (kPa): ", P_fluid/1.0E3
           write(95,*) "Max particle diameter (m): ", max_diap
           write(95,*) "BLT is ", pos_top*1.0E6, " um ", & 
                "at squeeze pressure of ", P_supp/1.0E3, " kPa"
           write(95,*) "Top plate position-to-particle ratio: ", &
                                pos_top/max_diap
           write(95,*) "Change in top plate position: ", pos_change     
           !read(*,*)
           max_rad_loc = MAXLOC(DES_POS_NEW(:,2), DIM=1)
           top_particle = MAXVAL(DES_POS_NEW(:,2)) + &
                           DES_RADIUS(max_rad_loc)     

           CALL WRITE_VTP_FILE(1,0) ! Write final VTP file for ParaView analysis
     
           write(107,*) top_overlap
           write(107,*) " "
           write(77, *) MAXVAL(top_overlap), pos_top, pos_change, &
                     P_sqz, F_total, top_particle

           call mfix_exit(0)
         ENDIF
  
      IF(MOD(flag_txt, calls_write) == 0) THEN 
        max_rad_loc = MAXLOC(DES_POS_NEW(:,2), DIM=1)
        top_particle = MAXVAL(DES_POS_NEW(:,2)) + & 
                       DES_RADIUS(max_rad_loc)     
       
        write(77, *) MAXVAL(top_overlap), pos_top, pos_change, &
                     P_sqz, F_total, top_particle
                           
      ENDIF
      
      IF(MOD(flag_txt, 200)==0) THEN
        write(107,*) MAXVAL(top_overlap)
        write(107,*) " "
      ENDIF        

      RETURN
      END SUBROUTINE USR1_DES