des_granular_temperature.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DES_GRANULAR_TEMPERATURE                                C
!  Purpose: Calculate the DES granular temperature                     C
!                                                                      C
!                                                                      C
!  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
!  Reviewer:                                          Date:            C
!  Revision: For parallel processing modifications are made to         C
!            accomodate ghost cells  (Pradeep G.)                      C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE des_granular_temperature

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE discretelement
      USE fldvar, only: u_s, v_s, w_s, theta_m
      USE fun_avg
      USE functions
      USE geometry
      USE indices
      USE mpi_utility
      USE param
      USE param1
      USE physprop, only: mmax
      USE run
      USE sendrecv
      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! indices
      INTEGER :: I, J, K, IJK
!
      INTEGER :: M, LL
! counter for no. of particles in phase m in cell ijk
      INTEGER :: NP_PHASE(dimension_3, dimension_m)
! temporary variable for mth phase granular temperature in cell ijk
      DOUBLE PRECISION :: TEMP(dimension_3, dimension_m)
! accounted for particles
      INTEGER :: PC
! squared particle velocity v.v
      DOUBLE PRECISION :: SQR_VEL, SQR_ROT_VEL
!-----------------------------------------------

! Calculate a local species granular temperature for current instant of
! time.  Note that the following calculation of species granular
! temperature employs a fluctuating particle velocity that is defined
! as the difference between a particles velocity and the corresponding
! local mean velocity of that particles species evaluated at the same
! instant of time.
!----------------------------------------------------------------->>>
! The following calculations are performed on the 'fluid' grid
      temp(:,:) = zero
      np_phase(:,:) = 0
      pc = 0
      DO ll = 1, max_pip
! skipping particles that do not exist
         IF(is_nonexistent(ll)) cycle
         pc = pc + 1
! skipping ghost particles
         IF(is_ghost(ll) .or. is_entering_ghost(ll) .or. is_exiting_ghost(ll)) cycle

         i = pijk(ll,1)
         j = pijk(ll,2)
         k = pijk(ll,3)
         ijk = pijk(ll,4)

         m = pijk(ll,5)
         np_phase(ijk,m) = np_phase(ijk,m) + 1

         temp(ijk,m) = temp(ijk,m) + &
            (des_vel_new(ll,1)-u_s(ijk,m))**2
         temp(ijk,m) = temp(ijk,m) + &
            (des_vel_new(ll,2)-v_s(ijk,m))**2
         IF(do_k) THEN
            temp(ijk,m) = temp(ijk,m) + &
               (des_vel_new(ll,3)-w_s(ijk,m))**2
         ENDIF

         IF(pc .EQ. pip) EXIT
      ENDDO

! loop over all fluid cells
      DO ijk = ijkstart3, ijkend3
         IF(fluid_at(ijk)) THEN

            DO m = mmax+1,des_mmax+mmax
               IF (np_phase(ijk,m) > 0 ) THEN
                  theta_m(ijk,m) = temp(ijk,m)/&
                     dble(dimn*np_phase(ijk,m))
               ELSE
                  theta_m(ijk,m) = zero
               ENDIF
            ENDDO
         ENDIF
      ENDDO


! Calculate global quantities: granular temperature, kinetic energy,
! potential energy and average velocity at the current instant of time
!----------------------------------------------------------------->>>

! initialization for calculations
      des_ke = zero
      des_rote = zero
      des_pe = zero
      des_vel_avg(:) = zero

! Calculate global average velocity, kinetic energy &
! potential energy
      pc = 0
      DO ll = 1, max_pip
! skipping ghost particles and particles that don't exist
         IF(is_nonexistent(ll)) cycle
         pc = pc + 1
         IF(is_ghost(ll) .OR. is_entering_ghost(ll) .OR. is_exiting_ghost(ll)) cycle

         sqr_vel = zero
         sqr_rot_vel = zero
         DO i = 1, dimn
            sqr_vel = sqr_vel + des_vel_new(ll,i)**2
         ENDDO

         DO i = 1, merge(1,3,no_k)
            sqr_rot_vel = sqr_rot_vel + omega_new(ll,i)**2
         ENDDO


         des_ke = des_ke + pmass(ll)/2.d0 * sqr_vel
! Calculation of rotational kinetic energy
         des_rote = des_rote +                                         &
            (0.4d0*pmass(ll)*des_radius(ll)**2)/2.d0 * sqr_rot_vel
         des_pe = des_pe + pmass(ll)*dble(abs(grav(2)))*&
            des_pos_new(ll,2)
         des_vel_avg(:) =  des_vel_avg(:) + des_vel_new(ll,:)

         IF(pc .EQ. pip) EXIT
      ENDDO

!Calculating total number of particles in the entire domain
!Needed for correct average velocities and granular temp.
      CALL global_all_sum(pip-ighost_cnt,tot_par)
      CALL global_all_sum(des_vel_avg(1:dimn))

!J.Musser changed PARTICLES TO PIP
      IF(tot_par > 0) des_vel_avg(:) = des_vel_avg(:)/dble(tot_par)

! The following quantities are primarily used for debugging/developing
! and allow a quick check of the energy conservation in the system.
! In their current form they are best applied to monodisperse cases.
! Calculate x,y,z components of global energy & granular temperature
      global_gran_energy(:) = zero
      global_gran_temp(:)  = zero
      pc = 0
      DO ll = 1, max_pip

! skipping ghost particles and particles that don't exist
         IF(is_nonexistent(ll)) cycle
         pc = pc + 1
         IF(is_ghost(ll) .OR. is_entering_ghost(ll) .OR. is_exiting_ghost(ll)) cycle

         global_gran_energy(:) = global_gran_energy(:) + &
            0.5d0*pmass(ll)*(des_vel_new(ll,:)-des_vel_avg(:))**2
         global_gran_temp(:) = global_gran_temp(:) + &
            (des_vel_new(ll,:)-des_vel_avg(:))**2

         IF(pc .EQ. pip) EXIT
      ENDDO

      CALL global_all_sum(global_gran_temp)
      CALL global_all_sum(global_gran_energy)
      CALL global_all_sum(des_ke)
      CALL global_all_sum(des_pe)
      CALL global_all_sum(des_rote)

      IF(tot_par > 0) global_gran_energy(:) =                          &
         global_gran_energy(:)/dble(tot_par)
      IF(tot_par > 0) global_gran_temp(:) =                            &
         global_gran_temp(:)/dble(tot_par)

      RETURN

      END SUBROUTINE des_granular_temperature

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CALC_DES_BEDHEIGHT                                      !
!  Purpose: Calculate an average bed height for each solids phase      !
!           present                                                    !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!

      SUBROUTINE calc_des_bedheight

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE des_bc
      USE discretelement
      USE fldvar
      USE functions
      USE geometry
      USE indices
      USE param, only: dimension_m
      USE param1, only: zero
      USE physprop
      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! particle index
      INTEGER :: L
! accounted for particles
      INTEGER :: PC
! solids phase no.
      INTEGER :: M
! ijk indices
      INTEGER :: J, IJK
! average height of fluid cell
      DOUBLE PRECISION :: hcell
! height of particle (y-position)
      DOUBLE PRECISION :: hpart
! volume fraction of phase M in fluid cell
      DOUBLE PRECISION :: EP_SM
! tmp variables for calculations
      DOUBLE PRECISION :: tmp_num(dimension_m), tmp_den(dimension_m)
!-----------------------------------------------

! calculation of bed height following the formulation of Goldschmidt et
! al., Powder Technology 138, 2003, 135-159
      tmp_num(:) = zero
      tmp_den(:) = zero
      DO ijk = ijkstart3, ijkend3
         j = j_of(ijk)
         DO m = mmax+1, des_mmax+mmax
            IF(rop_s(ijk,m) > zero) THEN
               hcell = 0.5d0*(yn(j)+yn(j-1))
               ep_sm = ep_s(ijk,m)
               tmp_num(m) = tmp_num(m) + ep_sm*hcell*vol(ijk)
               tmp_den(m) = tmp_den(m) + ep_sm*vol(ijk)
            ENDIF
         ENDDO
      ENDDO
! calculate avg height for each phase
      IF (pip >0) bed_height(:) = tmp_num(:)/tmp_den(:)

! alternative method to calculating bed height (turned off atm)
      IF(.false.) THEN
      tmp_num(:) = zero
      tmp_den(:) = zero

      pc = 0
      DO l = 1, max_pip
! skipping ghost particles and particles that don't exist
         IF(is_nonexistent(l)) cycle
         pc = pc + 1
         IF(is_ghost(l) .OR. is_entering_ghost(l) .OR. is_exiting_ghost(l)) cycle

         m = pijk(l,5)
         hpart = des_pos_new(l,2)
         tmp_num(m) = tmp_num(m) + hpart
         tmp_den(m) = tmp_den(m) + 1
         IF(pc .EQ. pip) EXIT
      ENDDO
      ! calculate avg height for each phase
      IF (pip >0) bed_height(:) = tmp_num(:)/tmp_den(:)
      ENDIF

      RETURN

      END SUBROUTINE calc_des_bedheight