calc_dif_g.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_DIF_g                                              C
!  Purpose: Calculate the effective diffusivity of fluid phase         C
!                                                                      C
!                                                                      C
!  Comments:                                                           C
!  This routine will not be called if dif_g is defined                 C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_dif_g()

! Modules
!---------------------------------------------------------------------//
      USE param1, only: undefined
      USE physprop, only: dif_g0, dif_g
      USE sendrecv, only: send_recv
! invoke user defined quantity
      USE usr_prop, only: usr_difg, calc_usr_prop
      USE usr_prop, only: gas_diffusivity
      IMPLICIT NONE
!---------------------------------------------------------------------//

      IF (usr_difg) THEN
         CALL calc_usr_prop(gas_diffusivity,lm=0)
      ELSEIF (dif_g0 == undefined) THEN
! unncessary check but included for clarity
         CALL calc_default_dif_gas
      ENDIF

      CALL send_recv(dif_g, 2)

      RETURN
      END SUBROUTINE calc_dif_g


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_DEFAULT_DIF_GAS                                    C
!  Purpose: Compute the default value for diffusivity of each gas      C
!  species; each species is assigned the same value with the base      C
!  value corresponding to CO2 in N2 at T=298K and P~=1atm.             C
!                                                                      C
!  Author:M. Syamlal                                  Date: 13-FEB-98  C
!                                                                      C
!  Revision: Include dilute mixture approximation for calculation of   C
!  multicomponent diffusion coefficients                               C
!  Author:N. Reuge                                    Date: 11-APR-07  C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!  Fuller relation - to account for influence of gas temperature       C
!     and pressure                                                     C
!  Curtiss-Hirschfelder, Wilke & Blanc - dilute mixture approximation  C
!     for multicomponent diffusion. Valid if the mass fraction of the  C
!     carrier species > 0.9                                            C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_default_dif_gas

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use fldvar, only: rop_g, t_g, p_g, x_g
      use functions, only: fluid_at
      use param1, only: zero, one
      use physprop, only: nmax, dif_g
      use scales, only: unscale_pressure
      use toleranc, only: zero_x_gs
      IMPLICIT NONE

! Local Variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK
! Species index
      INTEGER :: N, L
! Binary diffusion coefficient
      DOUBLE PRECISION :: Dab(nmax(0),nmax(0))
! Reference temperature and pressure for each species diffusion
! coefficient
      DOUBLE PRECISION :: Tg_ref(nmax(0))
      DOUBLE PRECISION :: Pg_ref(nmax(0))
! Intermediate calculation to determine weighted average diffusion
! coefficient
      DOUBLE PRECISION :: lDab, Sum_XgLoDab, Sum_XgL, lXgN
!---------------------------------------------------------------------//

      CALL set_binary_dab_gas(dab, tg_ref, pg_ref)

!!$omp  parallel do private(ijk) &
!!$omp& schedule(dynamic,chunk_size)
      DO n = 1, nmax(0)
         DO ijk = ijkstart3, ijkend3
            IF (fluid_at(ijk)) THEN

               IF (nmax(0) == 1) THEN
! Only 1 species present in gas phase
                  ldab = dab(n,n)

               ELSE
! Dilute mixture approximation for multi-component diffusion

                  sum_xglodab = zero
                  sum_xgl = zero
                  DO l = 1, nmax(0)
                     IF (l /= n) THEN
                        sum_xglodab = sum_xglodab+&
                                     x_g(ijk,l)/dab(n,l)
! Sum_XgL = 1-XgN
                        sum_xgl = sum_xgl + x_g(ijk,l)
                     ENDIF
                  ENDDO
! it is possible for lXgN to evaluate <0
                  lxgn = one-sum_xgl

                  IF (lxgn > zero_x_gs .AND. &
                     sum_xgl > zero_x_gs) THEN
! i.e. when cell is not only species N
! If this criteria is too strict (i.e. when XgN->1), then this section
! may be evaluated when the other XgN do not carry significant value
! (i.e. are effectively zero). As a result, lDab may be evaluated and
! become unrealistically large. This may happen when happen when
! 1-XgN != sum_XgL. Generally, sum_XgL should equal 1-XgN but they may
! differ due to the numerical evaluation of each Xg. If XgN->1, then
! use both sum_XgL and 1-XgN to determine whether calculations should
! proceed. Given the noted limitation of this approximation then an
! additional criteria should probably be added to only evaluate when
! sum_xgL <0.1.
                     IF (sum_xglodab > zero) THEN
! for numerical reasons use Sum_XgL rather than ONE-X_g(IJK,N)
                        ldab = sum_xgl/sum_xglodab
                     ELSE
! this should not occur...
                        ldab = dab(n,n)
                     ENDIF
                  ELSE
! Address case when the mass fraction of the Nth species is nearly 1.
                     ldab = dab(n,n)
                  ENDIF
               ENDIF

! Influence of gas temperature and gas pressure from Fuller relation
               dif_g(ijk,n) = rop_g(ijk)*ldab* &
                           (t_g(ijk)/tg_ref(n))**1.75* &
                           pg_ref(n)/unscale_pressure(p_g(ijk))
            ELSE
               dif_g(ijk,n) = zero
            ENDIF
         ENDDO
      ENDDO

      RETURN
      END SUBROUTINE calc_default_dif_gas


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_BINARY_DAB_GAS                                      C
!  Purpose: Set binary diffusion coefficients for all cross species    C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!  Bird, Stewart and lightfoot (1960)                                  C
!  Reid, Prausnitz and Poling (1987)                                   C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_binary_dab_gas(lDab, lTg_ref, lPg_ref)

! Modules
!---------------------------------------------------------------------//
      use physprop, only: nmax
      use run, only: units
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Binary diffusion coefficient
      DOUBLE PRECISION, INTENT(OUT) :: lDab(nmax(0),nmax(0))
      DOUBLE PRECISION, INTENT(OUT) :: lTg_ref(nmax(0))
      DOUBLE PRECISION, INTENT(OUT) :: lPg_ref(nmax(0))

! Local Variables
!---------------------------------------------------------------------//
! Species index
      INTEGER :: N, L
!---------------------------------------------------------------------//


! Default gas diffusion coefficient
! Bird, Stewart, and Lightfoot (1960) -- CO2--N2 at 298.2 K
      DO n = 1,nmax(0)
         DO l = n,nmax(0)
!         DO L = N+1, NMAX(0)
!            IF (N /= L) THEN
! assign the diaganol case for when nmax = 1 or when only species N present
! in given cell.
               ldab(n,l) = 0.165d0       !cm^2/s
               ldab(l,n) = ldab(n,l)
!            ENDIF
         ENDDO
         ltg_ref(n) = 298.2d0
         lpg_ref(n) = 1.01d6             !dyne
      ENDDO

! Gas diffusion coefficients for 3 species system: SiH4, H2 & N2
! Calculated using relation derived from Chapman and Enskog's theory
! of gases - Reid, Prausnitz and Poling (1987)
! Binary diffusion coefficient SiH4/H2 at 873 K
!      lDab(1,2) = 3.78      ! cm^2/s
!      lDab(2,1) = lDab(1,2)
! Binary diffusion coefficient SiH4/N2 at 873 K
!      lDab(1,3) = 1.02      ! cm^2/s
!      lDab(3,1) = lDab(1,3)
! Binary diffusion coefficient H2/N2 at 873 K
!      lDab(2,3) = 4.52      ! cm^2/s
!      lDab(3,2) = lDab(2,3)
!      DO N = 1,NMAX(0)
!         lTg_ref(N) = 873.0
!         lPg_ref(N) = 1.01e6         ! dyne
!      ENDDO

      IF(units == 'SI') THEN
         DO n = 1,nmax(0)
            DO l = n,nmax(0)
!               IF (N /= L) THEN
! assign the diaganol case nmax = 1 or when only species N present
! in given cell.
                  ldab(n,l) = ldab(n,l)*0.0001d0   !m^2/s
                  ldab(l,n) = ldab(n,l)
!               ENDIF
            ENDDO
            lpg_ref(n) = lpg_ref(n)/10.d0          !Pa
         ENDDO
      ENDIF

      RETURN
      END SUBROUTINE set_binary_dab_gas