calc_pg_grad.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CALC_PG_GRAD                                            !
!  Purpose: Calculate cell centered pressure force exerted on the      !
!           particles in the cell by the gas/fluid phase               !
!           Note that P_force is evaluated as -dp/dx                   !
!                                                                      !
!  Notes: This pressure force only needs to be calculated once during  !
!         the DEM loop (at the beginning) since the gas/fluid phase    !
!         is essentially static at that point (i.e., gas field is not  !
!         updated during DEM loop                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_pg_grad

! Model B momentum equation
      use run, only: model_b
! Particle volume.
      use discretelement, only: pvol
! Gas pressure force by fluid cell
      use discretelement, only: p_force
! Particle drag force
      use discretelement, only: drag_fc
! Flag for 3D simulatoins.
      use geometry, only: do_k
! Loop bounds for fluid grid
      USE compar, only: ijkstart3, ijkend3
! Flags for cyclic BC with pressure drop 
      use geometry, only: cyclic_x_pd, cyclic_y_pd, cyclic_z_pd
! Specified pressure drop
      use bc, only: delp_x, delp_y, delp_z
! Domain length
      use geometry, only: xlength, ylength, zlength
! Gas phase pressure
      use fldvar, only: p_g

      use discretelement, only: max_pip, pijk, des_explicitly_coupled
      use particle_filter, only: filter_cell
      use particle_filter, only: filter_weight
      use particle_filter, only: des_interp_on

      use functions, only: fluid_at

      use functions, only: is_normal

! Global Parameters:
!---------------------------------------------------------------------//
! Double precision values.
      use param1, only: zero

      implicit none

! Loop counters: Particle, fluid cell, neighbor cells
      INTEGER :: NP, IJK, LC
! Interpolation weight
      DOUBLE PRECISION :: WEIGHT
! Interpolated gas phase quanties.
      DOUBLE PRECISION :: lPF(3)
! Loop bound for
      INTEGER :: LP_BND
! mean pressure gradient for the case of periodic boundaries
      DOUBLE PRECISION :: cPG(3)
!......................................................................!

! Calculate the gas phase pressure gradient. (dP/dx)
      CALL calc_grad_des(p_g, p_force)


! Add in cyclic BC pressure drop.
      cpg(1) = merge(delp_x/xlength, zero, cyclic_x_pd)
      cpg(2) = merge(delp_y/ylength, zero, cyclic_y_pd)
      cpg(3) = merge(delp_z/zlength, zero, cyclic_z_pd)

      DO ijk=ijkstart3, ijkend3
         p_force(:,ijk) = cpg - p_force(:,ijk)
      ENDDO


      IF(des_explicitly_coupled .AND. .NOT.model_b) THEN

! Loop bounds for interpolation.
         lp_bnd = merge(27,9,do_k)

! Calculate the gas phase forces acting on each particle.

!$omp  parallel do default(none) &
!$omp  private(NP,lPF,lc,ijk,weight) &
!$omp  shared(MAX_PIP,DES_INTERP_ON,LP_BND,p_force,drag_fc, &
!$omp     filter_cell,filter_weight,pijk,pvol)
         DO np=1,max_pip

            IF(is_normal(np)) THEN
               IF(.NOT.fluid_at(pijk(np,4))) cycle

               IF(des_interp_on) THEN
                  lpf = zero
                  DO lc=1,lp_bnd
                     ijk = filter_cell(lc,np)
                     weight = filter_weight(lc,np)
                     lpf = lpf + p_force(:,ijk)*weight
                  ENDDO
               ELSE
                  lpf = p_force(:,pijk(np,4))
               ENDIF

! Include gas pressure and gas-solids drag
               drag_fc(np,:) = drag_fc(np,:) + lpf*pvol(np)
            ENDIF

         ENDDO
      ENDIF

      RETURN
      END SUBROUTINE calc_pg_grad