drag_gs_des1.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CALC_GS_DES1                                            !
!  Author: J.Musser                                   Date: 21-NOV-14  !
!                                                                      !
!  Purpose: This routine is called from the DISCRETE side to calculate !
!  the gas-based forces acting on each particle using interpolated     !
!  values for gas velocity, gas pressure, and gas volume fraction.     !
!                                                                      !
!  Notes:                                                              !
!                                                                      !
!   F_gp is obtained from des_drag_gp subroutine is given as:          !
!    F_GP = beta*VOL_P/EP_s where VOL_P is the particle volume.        !
!                                                                      !
!  The drag force on each particle is equal to:                        !
!    D_FORCE = beta*VOL_P/EP_s*(Ug - Us) = F_GP *(Ug - Us)             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE drag_gs_des1

! Modules
!---------------------------------------------------------------------//
! Gas phase volume fraction
      use fldvar, only: ep_g
! Gas phase velocities
      use fldvar, only: u_g, v_g, w_g
! Function to deterine if a cell contains fluid.
      use functions, only: fluid_at
      use functions, only: is_normal
! Flag for 3D simulatoins.
      use geometry, only: do_k
! Size of particle arrays on this processor.
      use discretelement, only: max_pip
! Flag: The fluid and discrete solids are explicitly coupled.
      use discretelement, only: des_explicitly_coupled
! IJK of fluid cell containing particles center
      use discretelement, only: pijk
! Drag force on each particle
      use discretelement, only: f_gp
! Particle velocity
      use discretelement, only: des_vel_new
! Total forces acting on particle
      use discretelement, only: fc
! Gas pressure force by fluid cell
      use discretelement, only: p_force
! Particle volume.
      use discretelement, only: pvol
! Particle drag force
      use discretelement, only: drag_fc
! Flag for MPPIC runs.
      use mfix_pic, only: mppic
! Flag to use implicit drag for MPPIC
      use mfix_pic, only: mppic_pdrag_implicit
! Double precision values.
      use param1, only: zero
      use param, only: dimension_3
! Flag to use interpolation
      use particle_filter, only: des_interp_on
! Interpolation cells and weights
      use particle_filter, only: filter_cell, filter_weight
! Model B momentum equation
      use run, only: model_b
      use param, only: dimension_3
      use compar
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
!     Cell-center gas velocities.
      DOUBLE PRECISION :: UGC(dimension_3)
      DOUBLE PRECISION :: VGC(dimension_3)
      DOUBLE PRECISION :: WGC(dimension_3)
! Loop counters: Particle, fluid cell, neighbor cells
      INTEGER :: NP, IJK, LC
! Interpolation weight
      DOUBLE PRECISION :: WEIGHT
! Interpolated gas phase quantities and particle velocity.
      DOUBLE PRECISION :: lEPg, VELFP(3), lPF(3), vel_p(3)
! Drag force acting on each particle.
      DOUBLE PRECISION :: D_FORCE(3)
! Flag for Model A momentum equation
      LOGICAL :: MODEL_A
! Loop bound for filter
      INTEGER :: LP_BND
!......................................................................!

! Set flag for Model A momentum equation.
      model_a = .NOT.model_b
! Loop bounds for interpolation.
      lp_bnd = merge(27,9,do_k)

! Calculate the cell center gas velocities.
      CALL calc_cell_center_gas_vel(u_g, v_g, w_g, ugc, vgc, wgc)

! Calculate the gas phase forces acting on each particle.

!$omp parallel default(none) private(np,lepg,velfp,ijk,weight,lpf,d_force,vel_p)    &
!$omp          shared(max_pip,des_interp_on,lp_bnd,filter_cell,filter_weight, &
!$omp          ep_g,pijk,des_vel_new,f_gp,mppic,ugc,vgc,wgc,p_force,    &
!$omp          des_explicitly_coupled,drag_fc,mppic_pdrag_implicit,fc,model_a,pvol)
!$omp do
      DO np=1,max_pip
         IF(.NOT.is_normal(np)) cycle

         lepg = zero
         velfp = zero
         lpf = zero

! Calculate the gas volume fraction, velocity, and pressure force at
! the particle's position.
         IF(des_interp_on) THEN
            DO lc=1,lp_bnd
               ijk = filter_cell(lc,np)
               weight = filter_weight(lc,np)
! Gas phase volume fraction.
               lepg = lepg + ep_g(ijk)*weight
! Gas phase velocity.
               velfp(1) = velfp(1) + ugc(ijk)*weight
               velfp(2) = velfp(2) + vgc(ijk)*weight
               velfp(3) = velfp(3) + wgc(ijk)*weight
! Gas pressure force.
               lpf = lpf + p_force(:,ijk)*weight
            ENDDO
         ELSE
            ijk = pijk(np,4)
            lepg = ep_g(ijk)
            velfp(1) = ugc(ijk)
            velfp(2) = vgc(ijk)
            velfp(3) = wgc(ijk)
            lpf = p_force(:,ijk)
         ENDIF

! Avoid drag calculations in cells without fluid (cut-cell)
         IF(.NOT.fluid_at(pijk(np,4))) THEN
            drag_fc(np,:) = 0.0d0
            f_gp(np) = 0.0d0

! For explicit coupling, use the drag coefficient calculated for the
! gas phase drag calculations.
         ELSEIF(des_explicitly_coupled) THEN
            drag_fc(np,:) = f_gp(np)*(velfp - des_vel_new(np,:))
         ELSE

! Calculate the drag coefficient.
            vel_p = des_vel_new(np,:)
            CALL des_drag_gp(np, vel_p, velfp, lepg)

! Calculate the gas-solids drag force on the particle
            IF(mppic .AND. mppic_pdrag_implicit) THEN
! implicit treatment of the drag term for mppic
               d_force = f_gp(np)*velfp
            ELSE
! default case
               d_force = f_gp(np)*(velfp - des_vel_new(np,:))
            ENDIF

! Update the contact forces (FC) on the particle to include gas
! pressure and gas-solids drag
            fc(np,:) = fc(np,:) + d_force(:)
            IF(model_a) fc(np,:) = fc(np,:) + lpf*pvol(np)

         ENDIF

      ENDDO
!$omp end parallel

      RETURN
      END SUBROUTINE drag_gs_des1

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DRAG_GS_GAS1                                            !
!  Author: J.Musser                                   Date: 21-NOV-14  !
!                                                                      !
!                                                                      !
!  Purpose: This routine is called from the CONTINUUM. It calculates   !
!  the scalar cell center drag force acting on the fluid using         !
!  interpolated values for the gas velocity and volume fraction. The   !
!  The resulting sources are interpolated back to the fluid grid.      !
!                                                                      !
!  NOTE: The loop over particles includes ghost particles so that MPI  !
!  communications are needed to distribute overlapping force between   !
!  neighboring grid cells. This is possible because only cells "owned" !
!  by the current process will have non-zero weights.                  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE drag_gs_gas1

! Flag: The fluid and discrete solids are explicitly coupled.
      use discretelement, only: des_explicitly_coupled
! Size of particle array on this process.
      use discretelement, only: max_pip
! IJK of fluid cell containing particles center
      use discretelement, only: pijk
! Drag force on each particle
      use discretelement, only: f_gp
! Particle velocity
      use discretelement, only: des_vel_new
! Contribution to gas momentum equation due to drag
      use discretelement, only: drag_bm
! Scalar cell center total drag force
      use discretelement, only: f_gds
! Gas phase volume fraction
      use fldvar, only: ep_g
! Gas phase velocities
      use fldvar, only: u_g, v_g, w_g
      use fldvar, only: u_go, v_go, w_go
      use functions, only: fluid_at
      use functions, only: is_normal
! Volume of scalar cell.
      use geometry, only: vol
! Flag for 3D simulatoins.
      use geometry, only: do_k
! Flag for MPPIC runs
      use mfix_pic, only: mppic
! Statical weight of each MPPIC parcel
      use mfix_pic, only: des_stat_wt
! Double precision values.
      use param1, only: zero, one
! Flag to use interpolation
      use particle_filter, only: des_interp_on
! Interpolation cells and weights
      use particle_filter, only: filter_cell, filter_weight
! Array size for fluid variables
      use param, only: dimension_3
! MPI function for collecting interpolated data from ghost cells.
      use sendrecvnode, only: des_collect_gdata
! MPI wrapper for halo exchange.
      use sendrecv, only: send_recv

      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
!     Cell-center gas velocities.
      DOUBLE PRECISION :: UGC(dimension_3)
      DOUBLE PRECISION :: VGC(dimension_3)
      DOUBLE PRECISION :: WGC(dimension_3)
! Loop counters: Particle, fluid cell, neighbor cells
      INTEGER :: NP, IJK, LC
! Interpolation weight
      DOUBLE PRECISION :: WEIGHT
! Interpolated gas phase quantities and particle velocity
      DOUBLE PRECISION :: lEPg, VELFP(3), vel_p(3)
! Loop bound for filter
      INTEGER :: LP_BND
! Drag force (intermediate calculation)
      DOUBLE PRECISION :: lFORCE
! Drag sources for fluid (intermediate calculation)
      DOUBLE PRECISION :: lDRAG_BM(3)
!......................................................................!

! Initialize fluid cell values.
      f_gds = zero
      drag_bm = zero

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

! Calculate the cell center gas velocities.
      CALL calc_cell_center_gas_vel(u_g, v_g, w_g, ugc, vgc, wgc)

! Calculate the gas phase forces acting on each particle.

!$omp parallel default(none) private(np,lepg,velfp,ijk,weight,ldrag_bm,lforce,vel_p) &
!$omp          shared(max_pip,des_interp_on,lp_bnd,filter_cell,filter_weight, &
!$omp          ep_g,pijk,des_vel_new,f_gp,vol,des_stat_wt,mppic,drag_bm,f_gds,ugc,vgc,wgc)
!$omp do
      DO np=1,max_pip

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

         lepg = zero
         velfp = zero

! Calculate the gas volume fraction, velocity, and at the
! particle's position.
         IF(des_interp_on) THEN
            DO lc=1,lp_bnd
               ijk = filter_cell(lc,np)
               weight = filter_weight(lc,np)
! Gas phase volume fraction.
               lepg = lepg + ep_g(ijk)*weight
! Gas phase velocity.
               velfp(1) = velfp(1) + ugc(ijk)*weight
               velfp(2) = velfp(2) + vgc(ijk)*weight
               velfp(3) = velfp(3) + wgc(ijk)*weight
            ENDDO
         ELSE
            ijk = pijk(np,4)
            lepg = ep_g(ijk)
            velfp(1) = ugc(ijk)
            velfp(2) = vgc(ijk)
            velfp(3) = wgc(ijk)
         ENDIF

! This avoids FP exceptions for some ghost particles.
         IF(lepg == zero) lepg = ep_g(pijk(np,4))

! Calculate drag coefficient
         vel_p = des_vel_new(np,:)
         CALL des_drag_gp(np, vel_p, velfp, lepg)

         lforce = f_gp(np)
         IF(mppic) lforce = lforce*des_stat_wt(np)

         ldrag_bm = lforce*des_vel_new(np,:)

         IF(des_interp_on) THEN
            DO lc=1,lp_bnd
               ijk = filter_cell(lc,np)
               weight = filter_weight(lc,np)/vol(ijk)

!$omp atomic
               drag_bm(ijk,1) = drag_bm(ijk,1) + ldrag_bm(1)*weight
!$omp atomic
               drag_bm(ijk,2) = drag_bm(ijk,2) + ldrag_bm(2)*weight
!$omp atomic
               drag_bm(ijk,3) = drag_bm(ijk,3) + ldrag_bm(3)*weight
!$omp atomic
               f_gds(ijk) = f_gds(ijk) + lforce*weight
            ENDDO
         ELSE
            ijk = pijk(np,4)
            weight = one/vol(ijk)

!$omp atomic
            drag_bm(ijk,1) = drag_bm(ijk,1) + ldrag_bm(1)*weight
!$omp atomic
            drag_bm(ijk,2) = drag_bm(ijk,2) + ldrag_bm(2)*weight
!$omp atomic
            drag_bm(ijk,3) = drag_bm(ijk,3) + ldrag_bm(3)*weight

!$omp atomic
            f_gds(ijk) = f_gds(ijk) + lforce*weight
         ENDIF

      ENDDO
!$omp end parallel


! Add in data stored in ghost cells from interpolation. This call must
! preceed the SEND_RECV to avoid overwriting ghost cell data.
      IF(des_interp_on) THEN
         CALL des_collect_gdata(f_gds)
         CALL des_collect_gdata(drag_bm)
      ENDIF

! Update the drag force and sources in ghost layers.
      CALL send_recv(f_gds, 2)
      CALL send_recv(drag_bm, 2)

      RETURN
      END SUBROUTINE drag_gs_gas1

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CALC_CELL_CENTER_GAS_VEL                                !
!  Author: J.Musser                                   Date: 07-NOV-14  !
!                                                                      !
!  Purpose: Calculate the scalar cell center gas velocity. This code   !
!  is common to the DEM and GAS calls for non-interpolated drag        !
!  routines.                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_cell_center_gas_vel(lUg, lVg, lWg, UGC, VGC, WGC)

! Modules
!---------------------------------------------------------------------//
! Fluid grid loop bounds.
      use compar, only: ijkstart3, ijkend3
! Flags and correction factors for cut momentum cells.
      use cutcell, only: cut_u_treatment_at, theta_ue, theta_ue_bar
      use cutcell, only: cut_v_treatment_at, theta_vn, theta_vn_bar
      use cutcell, only: cut_w_treatment_at, theta_wt, theta_wt_bar
! Functions to average momentum to scalar cell center.
      use fun_avg, only: avg_x_e, avg_y_n, avg_z_t
! Functions to lookup adjacent cells by index.
      use functions, only: im_of, jm_of, km_of
! Function to deterine if a cell contains fluid.
      use functions, only: fluid_at
! Flag for 3D simulatoins.
      use geometry, only: do_k
      use indices, only: i_of
! Double precision parameters.
      use param, only: dimension_3
      use param1, only: zero

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
      DOUBLE PRECISION, INTENT(IN) :: lUg(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: lVg(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: lWg(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: UGC(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: VGC(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: WGC(dimension_3)
! Local variables:
!---------------------------------------------------------------------//
! Indices of adjacent cells
      INTEGER :: IJK, IMJK, IJMK, IJKM
!......................................................................!


! Calculate the cell center gas velocity components.
      DO ijk=ijkstart3, ijkend3
         IF(fluid_at(ijk)) THEN
            imjk = im_of(ijk)
            IF(cut_u_treatment_at(imjk)) THEN
               ugc(ijk) = (theta_ue_bar(imjk)*lug(imjk) +              &
                  theta_ue(imjk)*lug(ijk))
            ELSE
               ugc(ijk) = avg_x_e(lug(imjk),lug(ijk),i_of(ijk))
            ENDIF

            ijmk = jm_of(ijk)
            IF(cut_v_treatment_at(ijmk)) THEN
               vgc(ijk) = (theta_vn_bar(ijmk)*lvg(ijmk) +              &
                  theta_vn(ijmk)*lvg(ijk))
            ELSE
               vgc(ijk) = avg_y_n(lvg(ijmk),lvg(ijk))
            ENDIF

            IF(do_k) THEN
               ijkm = km_of(ijk)
               IF(cut_w_treatment_at(ijkm)) THEN
                  wgc(ijk) = (theta_wt_bar(ijkm)*lwg(ijkm) +           &
                     theta_wt(ijkm)* lwg(ijk))
               ELSE
                  wgc(ijk) = avg_z_t(lwg(ijkm),lwg(ijk))
               ENDIF
            ELSE
               wgc(ijk) = zero
            ENDIF
         ELSE
            ugc(ijk) = zero
            vgc(ijk) = zero
            wgc(ijk) = zero
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_cell_center_gas_vel