set_fluidbed_p.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_FLUIDBED_P                                          C
!  Purpose: Set the pressure field inside the bed assuming a fluidized C
!           bed with gravity acting the -ve y-direction                C
!                                                                      C
!  Author: M. Syamlal                                 Date: 21-JAN-92  C
!  Reviewer:M. Syamlal, S. Venkatesan, P. Nicoletti,  Date: 29-JAN-92  C
!           W. Rogers                                                  C
!                                                                      C
!  Revision Number: 1                                                  C
!  Purpose: Modifications for including cylindrical geometry           C
!  Author: M. Syamlal                                 Date:  6-MAR-92  C
!  Reviewer: M. Syamlal                               Date: 11-DEC-92  C
!  Revision Number: 2                                                  C
!  Purpose: Set pressure drop for cyclic boundary condition w/         C
!           pressure drop                                              C
!  Author: M. Syamlal                                 Date: 29-APR-94  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced: BC_DEFINED, BC_TYPE, IC_P_g, BC_P_g           C
!                        EP_g, MW_MIX_G, RO_g0, T_g,                   C
!                        SMAX, ROP_s,                                  C
!                        DX, DY, DZ, BFY_G, DELP_X, DELP_Y, DELP_Z,    C
!                        DO_I, DO_J, DO_K, IMIN1, KMIN1, JMIN1, IMAX1, C
!                        IMAX2, JMAX1, JMAX2, KMAX1, KMAX2             C
!  Variables modified: P_g                                             C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE set_fluidbed_p

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE bc
      USE bodyforce
      USE compar
      USE constant
      USE discretelement
      USE eos, ONLY: eosg
      USE exit, only: mfix_exit
      USE fldvar
      USE functions
      USE funits
      USE geometry
      USE ic
      USE indices
      USE machine, only: start_log, end_log
      USE mpi_utility
      USE param
      USE param1
      USE physprop
      USE scales
      USE sendrecv
      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! indices
      INTEGER :: I, J, K, IJK, M
! Local loop counter
      INTEGER :: L
! Gas pressure at the axial location j
      DOUBLE PRECISION :: PJ
! Bed weight per unit area
      DOUBLE PRECISION :: BED_WEIGHT
! Total area of a x-z plane
      DOUBLE PRECISION :: AREA
! x-z plane area of one cell
      DOUBLE PRECISION :: dAREA
! Average pressure drop per unit length
      DOUBLE PRECISION :: DPoDX, DPoDY, DPoDZ
!-----------------------------------------------

! If any initial pressures are unspecified skip next section
! calculations.
      DO l = 1, dimension_ic
         IF (ic_defined(l)) THEN
            IF (ic_p_g(l) == undefined) GOTO 60
            pj = ic_p_g(l)
         ENDIF
      ENDDO

! Here the pressure in each cell is determined from a specified pressure
! drop across the domain length. This section requires that the pressure
! is already defined in all initial condition regions (otherwise this
! section would be skipped)
! ---------------------------------------------------------------->>>
      IF (do_i .AND. delp_x/=undefined) THEN
         dpodx = delp_x/xlength
         pj = pj - dpodx*half*(dx(imax1)+dx(imax2))
         DO i = imax1, imin1, -1
            pj = pj + dpodx*half*(dx(i)+dx(i+1))
            DO k = kmin1, kmax1
               DO j = jmin1, jmax1
! Bound Checking
                  IF(.NOT.is_on_mype_owns(i,j,k)) cycle
                  IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
                  ijk = funijk(i,j,k)
                  IF (fluid_at(ijk)) p_g(ijk) = scale_pressure(pj)
               ENDDO
            ENDDO
         ENDDO
      ENDIF

      IF (do_j .AND. delp_y/=undefined) THEN
         dpody = delp_y/ylength
         pj = pj - dpody*half*(dy(jmax1)+dy(jmax2))
         DO j = jmax1, jmin1, -1
            pj = pj + dpody*half*(dy(j)+dy(j+1))
            DO k = kmin1, kmax1
               DO i = imin1, imax1
! Bound Checking
                  IF(.NOT.is_on_mype_owns(i,j,k)) cycle
                  IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
                  ijk = funijk(i,j,k)
                  IF (fluid_at(ijk)) p_g(ijk) = scale_pressure(pj)
               ENDDO
            ENDDO
         ENDDO
      ENDIF

      IF (do_k .AND. delp_z/=undefined) THEN
         dpodz = delp_z/zlength
         pj = pj - dpodz*half*(dz(kmax1)+dz(kmax2))
         DO k = kmax1, kmin1, -1
            pj = pj + dpodz*half*(dz(k)+dz(k+1))
            DO j = jmin1, jmax1
               DO i = imin1, imax1
! Bound Checking
                  IF(.NOT.is_on_mype_owns(i,j,k)) cycle
                  IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
                  ijk = funijk(i,j,k)
                  IF (fluid_at(ijk)) p_g(ijk) = scale_pressure(pj)
               ENDDO
            ENDDO
         ENDDO
      ENDIF
! ----------------------------------------------------------------<<<
      GOTO 100   ! pressure in all intial condition region cells was defined

   60 CONTINUE   ! pressure in an initial condition region cell was undefined


! ---------------------------------------------------------------->>>
! Search for an outflow boundary condition where pressure is specified
      pj = undefined
      DO l = 1, dimension_bc
         IF (bc_defined(l) .AND. bc_type_enum(l)==p_outflow) pj = bc_p_g(l)
      ENDDO

      IF (pj == undefined) THEN
! either a PO was not specified and/or a PO was specified but not the
! pressure at the outlet
         IF (ro_g0 /= undefined) THEN
! If incompressible flow set P_g to zero
            DO ijk = ijkstart3, ijkend3
               IF (fluid_at(ijk)) p_g(ijk) = zero
            ENDDO
            GOTO 100

         ELSE   ! compressible case

! Error condition -- no pressure outflow boundary condition is specified
! if a case is compressible and pressure in any of the initial
! conditions regions is unspecified, then a PO is effectively required
! (i.e., is specifies a bc_p_g).
            CALL start_log
            IF(dmp_log)WRITE (unit_log, 1000)
            CALL mfix_exit(mype)
         ENDIF
      ENDIF


! Set an approximate pressure field assuming that the pressure drop
! balances the weight of the bed, if the initial pressure-field is not
! specified
      DO j = jmax2, jmin1, -1

! Find the average weight per unit area over an x-z slice
         bed_weight = 0.0
         area = 0.0
         DO k = kmin1, kmax1
            DO i = imin1, imax1
               IF(.NOT.is_on_mype_owns(i,j,k)) cycle
               IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
               ijk = funijk(i,j,k)
               IF (fluid_at(ijk)) THEN
                  IF (coordinates == 'CARTESIAN') THEN
                     darea = dx(i)*dz(k)
                  ELSE IF (cylindrical) THEN
                     darea = dx(i)*x(i)*dz(k)
                  ENDIF
                  area = area + darea
                  IF (ro_g0 == undefined) THEN
                     bed_weight = bed_weight - dy(j)*bfy_g(ijk)*ep_g(ijk)*eosg(&
                        mw_mix_g(ijk),pj,t_g(ijk))*darea
                  ELSE
                     bed_weight = bed_weight - dy(j)*bfy_g(ijk)*ep_g(ijk)*ro_g0&
                        *darea
                  ENDIF
! This code is turned off for DEM runs until the value of rop_s can be
! ensured valid values for a DEM run at this point in the code.
                  IF (.NOT.discrete_element) THEN
                     DO m = 1, smax
                        bed_weight = bed_weight - dy(j)*bfy_s(ijk,m)*rop_s(ijk,m)*&
                           darea
                     ENDDO
                  ENDIF  ! end if (.not.discrete_element)
               ENDIF  ! end if (fluid_at(ijk))
            ENDDO    ! end do loop (i=imin1,imax1)
         ENDDO    ! end do loop (k=kmin1,kmax1)

! Global Sum
         call global_all_sum(bed_weight)
         call global_all_sum(area)
         IF (area /= 0.0) bed_weight = bed_weight/area

         pj = pj + bed_weight
         DO k = kmin1, kmax1
            DO i = imin1, imax1
               IF(.NOT.is_on_mype_owns(i,j,k)) cycle
               IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
               ijk = funijk(i,j,k)
               IF(fluid_at(ijk).AND.p_g(ijk)==undefined)p_g(ijk)=scale_pressure(pj)
            ENDDO    ! end do (i=imin1,imax1)
         ENDDO   ! end do (k = kmin1,kmax1)
      ENDDO   ! end do (j=jmax2,jimn1, -1)
! end setting an undefined pressure in an initial condition region
! ----------------------------------------------------------------<<<

  100 CONTINUE

      call send_recv(p_g,2)

      RETURN

 1000 FORMAT(/1x,70('*')//' From: SET_FLUIDBED_P'/' Message: Outflow ',&
         'pressure boundary condition (P_OUTFLOW) not found.',/&
         'All the initial pressures (IC_P_g) or at least one P_OUTFLOW',/&
         'condition need to be specified',/1x,70('*')/)

      END SUBROUTINE set_fluidbed_p