set_constprop.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_constprop                                           C
!  Purpose: This routine serves two purposes:                          C
!    1) initializes various variables everywhere in the domain with    C
!       a zero value. the real need for this is unclear. undefined     C
!       may be a better approach...                                    C
!    2) if defined, sets physical properties to their specified        C
!       constant value in the fluid domain. cannot set in flow         C
!       boundaries or later checks will also complain (may be          C
!       overly strict check)                                           C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_constprop

! Modules
!-----------------------------------------------
      USE param1, only: zero, half, one, undefined

      use discretelement, only: des_mmax

      USE fldvar, only: ro_g
      USE fldvar, only: ro_s, d_p
      use fldvar, only: p_s

      USE visc_g, only: mu_gt, epmu_gt, lambda_gt, eplambda_gt
      USE visc_s, only: mu_s, epmu_s, lambda_s, eplambda_s, lambda_s_c
      USE visc_s, only: ep_star_array
      USE visc_s, only: ep_g_blend_start, ep_g_blend_end

      USE physprop, only: nmax, mmax, smax
      USE physprop, only: ro_s0, ro_g0, d_p0
      USE physprop, only: mu_g0, mu_g, mu_s0
      USE physprop, only: mw_avg, mw_mix_g
      USE physprop, only: c_pg0, k_g0, c_pg, k_g, dif_g0, dif_g
      USE physprop, only: c_ps0, k_s0, c_ps, k_s, dif_s0, dif_s
      USE physprop, only: cv

      USE constant, only: ep_s_max_ratio, d_p_ratio, ep_s_max, m_max
      use constant, only: ep_star

      USE run, only: call_dqmom
      USE run, only: yu_standish, fedors_landel
      USE run, only: kt_type_enum, ia_2005, gd_1999, gtsh_2012
      USE run, only: blending_stress, sigm_blend, tanh_blend

      USE drag, only: f_gs, f_ss

      use kintheory, only: mu_sm_ip, mu_sl_ip, xi_sm_ip, xi_sl_ip
      use kintheory, only: fnu_s_ip, ft_sm_ip, ft_sl_ip
      use kintheory, only: kth_sl_ip, knu_sm_ip, knu_sl_ip, kvel_s_ip
      use kintheory, only: ed_ss_ip, edvel_sl_ip, edt_s_ip, edvel_sm_ip
      use kintheory, only: a2_gtsh, xsi_gtsh

      use mms, only: use_mms

      USE compar, only: ijkstart3, ijkend3
      use functions, only: wall_at, fluid_at

      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! indices
      INTEGER :: IJK, M, I, J
      DOUBLE PRECISION :: old_value, DP_TMP(smax)
!-----------------------------------------------

! First, initialize certain transport coefficients, physical
! properties, and other variable types everywhere in the
! domain to zero. Then, set these to their specified value
! (if defined) in fluid cells. Some are also set in flow cells.
! NOTE: DO NOT simply zero existing field variables.
      ro_g = zero
      mu_g = zero
      mu_gt = zero
      epmu_gt = zero
      lambda_gt = zero
      eplambda_gt = zero
      k_g = zero
      c_pg = zero
      mw_mix_g = zero
      dif_g = zero
      f_gs = zero

      ro_s = zero
      d_p = zero     ! this is done in init_fvars
      mu_s = zero
      epmu_s = zero
      lambda_s_c = zero
      lambda_s = zero
      eplambda_s = zero
! not defining p_s in flow cells can become a problem for certain
! cases as gradients in P_s are used (e.g., when a PO is west of a
! MI...?). Simply assigning it a zero value may not be the best
! approach but it is currently used....
      p_s = zero
      k_s = zero
      c_ps = zero
      dif_s = zero
      f_ss = zero

! Set the flag for recalculating gas viscosity.
!      RECALC_VISC_G = (MU_g0==UNDEFINED .OR. L_SCALE0/=ZERO .OR.&
!                       K_EPSILON .OR. ISHII)

! Set default value for virtual mass coefficient
      cv = half

! Variables specific to various kinetic theory models
      IF (kt_type_enum == ia_2005) THEN
         mu_sm_ip = zero
         mu_sl_ip = zero
         xi_sm_ip = zero
         xi_sl_ip = zero
         fnu_s_ip = zero
         ft_sm_ip = zero
         ft_sl_ip = zero
         kth_sl_ip = zero
         knu_sm_ip = zero
         knu_sl_ip = zero
         kvel_s_ip = zero
         ed_ss_ip = zero
         edvel_sl_ip = zero
      ENDIF

      IF (kt_type_enum == ia_2005 .OR. &
          kt_type_enum == gd_1999 .OR.  &
          kt_type_enum == gtsh_2012) THEN
         edt_s_ip = zero
         edvel_sm_ip = zero
      ENDIF

      IF(kt_type_enum == gtsh_2012) THEN
         a2_gtsh = zero
         xsi_gtsh = zero
      ENDIF

! Set specified constant physical properties values
      DO ijk = ijkstart3, ijkend3

         IF (.NOT.wall_at(ijk)) THEN
! Fluid and inflow/outflow cells: FLAG < 100
            IF (ro_g0 /= undefined) ro_g(ijk) = ro_g0
            IF (c_pg0 /= undefined) c_pg(ijk) = c_pg0
            IF (mw_avg /= undefined) mw_mix_g(ijk) = mw_avg
         ENDIF

         IF (fluid_at(ijk)) THEN
! Strictly Fluid cells: FLAG = 1
            IF (mu_g0 /= undefined) THEN
               mu_g(ijk) = mu_g0
               mu_gt(ijk) = mu_g0
               epmu_gt(ijk) = mu_g0
               lambda_gt(ijk) = -(2.0d0/3.0d0)*mu_g0
               eplambda_gt(ijk) = -(2.0d0/3.0d0)*mu_g0
            ENDIF
            IF (k_g0 /= undefined) k_g(ijk) = k_g0
            IF (dif_g0 /= undefined) dif_g(ijk,:nmax(0)) = dif_g0
         ENDIF

         IF (use_mms) THEN
            IF (mu_g0 /= undefined) THEN
               mu_g(ijk) = mu_g0
               mu_gt(ijk) = mu_g0
               epmu_gt(ijk) = mu_g0
               lambda_gt(ijk) = -(2.0d0/3.0d0)*mu_g0
               eplambda_gt(ijk) = -(2.0d0/3.0d0)*mu_g0
            ENDIF
            IF (k_g0 /= undefined) k_g(ijk) = k_g0
            IF (dif_g0 /= undefined) dif_g(ijk,:nmax(0)) = dif_g0
         ENDIF

      ENDDO

! ensure ro_s(ijk,m) is assigned to ro_s0(m) or ro_s(np) so that the
! function ep_s works for discrete phases. might be able to move this
! to set_ic_dem and set_ic_mppic. also required d_p(ijk,m) for hybrid
! use.  note check_solids_common_all ensures d_p0 is set for all
! solids defined also either ro_s0 must be set or ro_s0
      DO m = 1, mmax+des_mmax
         DO ijk = ijkstart3, ijkend3
            IF(.NOT.wall_at(ijk)) THEN
! Fluid and inflow/outflow cells: FLAG < 100
               IF (ro_s0(m) /= undefined) ro_s(ijk,m) = ro_s0(m)
               IF (c_ps0(m) /= undefined) c_ps(ijk,m) = c_ps0(m)
               IF (d_p0(m) /= undefined) d_p(ijk,m) = d_p0(m)
            ENDIF

            IF (fluid_at(ijk)) THEN
! Strictly fluid cells: FLAG = 1
               IF (mu_s0(m) /= undefined) THEN
                  mu_s(ijk,m) = mu_s0(m)
                  epmu_s(ijk,m) = mu_s0(m)
                  lambda_s(ijk,m) = (-2./3.)*mu_s(ijk,m)
                  eplambda_s(ijk,m) = (-2./3.)*mu_s(ijk,m)
               ENDIF
               IF (k_s0(m) /= undefined) k_s(ijk,m) = k_s0(m)
               IF (dif_s0(m) /= undefined) dif_s(ijk,m,:nmax(m)) = dif_s0(m)
            ENDIF

            IF (use_mms) THEN
               IF (mu_s0(m) /= undefined) THEN
                  mu_s(ijk,m) = mu_s0(m)
                  epmu_s(ijk,m) = mu_s0(m)
                  lambda_s(ijk,m) = (-2./3.)*mu_s(ijk,m)
                  eplambda_s(ijk,m) = (-2./3.)*mu_s(ijk,m)
               ENDIF
               IF (k_s0(m) /= undefined) k_s(ijk,m) = k_s0(m)
               IF (dif_s0(m) /= undefined) dif_s(ijk,m,:nmax(m)) = dif_s0(m)
            ENDIF

! set ep_star_array to user input ep_star in all cells.
            ep_star_array(ijk) = ep_star
! initializing blending stress parameters
            IF(blending_stress.AND.tanh_blend) THEN
               ep_g_blend_start(ijk) = ep_star_array(ijk) * 0.99d0
               ep_g_blend_end(ijk)   = ep_star_array(ijk) * 1.01d0
            ELSE IF(blending_stress.AND.sigm_blend) THEN
               ep_g_blend_start(ijk) = ep_star * 0.97d0
               ep_g_blend_end(ijk) = ep_star * 1.01d0
            ELSE
               ep_g_blend_start(ijk) = ep_star_array(ijk)
               ep_g_blend_end(ijk)   = ep_star_array(ijk)
            ENDIF

         ENDDO   ! end loop over ijk
      ENDDO   ! end loop over MMAX


! Initializing parameters needed if a correlation is used to compute
! ep_star: initializing the indexing system.
      IF(yu_standish .OR. fedors_landel) THEN
         DO m = 1, smax
            IF(ep_s_max(m) == undefined) ep_s_max(m) = one-ep_star
         ENDDO

         IF (.NOT.call_dqmom) THEN

! refer to Syam's dissertation
            IF (smax == 2) THEN
               ep_s_max_ratio(1,2) = ep_s_max(1)/ &
                  (ep_s_max(1)+(1.-ep_s_max(1))*ep_s_max(2))
            ENDIF

! initialize local variables
            DO i = 1, smax
               dp_tmp(i) = d_p0(i)
               m_max(i) = i
            ENDDO

! Rearrange the indices from coarsest particles to finest to be
! used in CALC_ep_star. Done here because it may need to be done
! for auto_restart
            DO i = 1, smax
               DO j = i, smax
                  IF(dp_tmp(i) < dp_tmp(j)) THEN
                     old_value = dp_tmp(i)
                     dp_tmp(i) = dp_tmp(j)
                     dp_tmp(j) = old_value
                  ENDIF
               ENDDO
            ENDDO

            DO i = 1, smax
               DO j = 1, smax
                  IF(dp_tmp(i) == d_p0(j) .AND. d_p0(i) .NE. d_p0(j)) THEN
                     m_max(i) = j
                  ENDIF
               ENDDO
            ENDDO
         ENDIF    ! if .not. call_dqmom
      ELSE   ! if .not. Yu-standish or Fedors-Landel
         ep_s_max(:) = zero
         ep_s_max_ratio(:,:) = zero
         d_p_ratio(:,:) = zero
         m_max(:) = 0
      ENDIF

      RETURN
      END SUBROUTINE set_constprop