1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: CALC_DIF_g C 4 ! Purpose: Calculate the effective diffusivity of fluid phase C 5 ! C 6 ! C 7 ! Comments: C 8 ! This routine will not be called if dif_g is defined C 9 ! C 10 ! C 11 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 12 SUBROUTINE CALC_DIF_G() 13 14 ! Modules 15 !---------------------------------------------------------------------// 16 USE param1, only: undefined 17 USE physprop, only: dif_g0, dif_g 18 USE sendrecv, only: send_recv 19 ! invoke user defined quantity 20 USE usr_prop, only: usr_difg, calc_usr_prop 21 USE usr_prop, only: gas_diffusivity 22 IMPLICIT NONE 23 !---------------------------------------------------------------------// 24 25 IF (USR_Difg) THEN 26 CALL CALC_USR_PROP(Gas_Diffusivity,lm=0) 27 ELSEIF (Dif_g0 == UNDEFINED) THEN 28 ! unncessary check but included for clarity 29 CALL CALC_DEFAULT_DIF_GAS 30 ENDIF 31 32 CALL send_recv(DIF_G, 2) 33 34 RETURN 35 END SUBROUTINE CALC_DIF_G 36 37 38 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 39 ! C 40 ! Subroutine: CALC_DEFAULT_DIF_GAS C 41 ! Purpose: Compute the default value for diffusivity of each gas C 42 ! species; each species is assigned the same value with the base C 43 ! value corresponding to CO2 in N2 at T=298K and P~=1atm. C 44 ! C 45 ! Author:M. Syamlal Date: 13-FEB-98 C 46 ! C 47 ! Revision: Include dilute mixture approximation for calculation of C 48 ! multicomponent diffusion coefficients C 49 ! Author:N. Reuge Date: 11-APR-07 C 50 ! C 51 ! C 52 ! Literature/Document References: C 53 ! Fuller relation - to account for influence of gas temperature C 54 ! and pressure C 55 ! Curtiss-Hirschfelder, Wilke & Blanc - dilute mixture approximation C 56 ! for multicomponent diffusion. Valid if the mass fraction of the C 57 ! carrier species > 0.9 C 58 ! C 59 ! C 60 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 61 SUBROUTINE CALC_DEFAULT_DIF_GAS 62 63 ! Modules 64 !---------------------------------------------------------------------// 65 use compar, only: ijkstart3, ijkend3 66 use fldvar, only: ROP_g, T_g, P_g, X_g 67 use functions, only: fluid_at 68 use param1, only: zero, one 69 use physprop, only: NMAX, Dif_g 70 use scales, only: unscale_pressure 71 use toleranc, only: zero_x_gs 72 IMPLICIT NONE 73 74 ! Local Variables 75 !---------------------------------------------------------------------// 76 ! Indices 77 INTEGER :: IJK 78 ! Species index 79 INTEGER :: N, L 80 ! Binary diffusion coefficient 81 DOUBLE PRECISION :: Dab(NMAX(0),NMAX(0)) 82 ! Reference temperature and pressure for each species diffusion 83 ! coefficient 84 DOUBLE PRECISION :: Tg_ref(NMAX(0)) 85 DOUBLE PRECISION :: Pg_ref(NMAX(0)) 86 ! Intermediate calculation to determine weighted average diffusion 87 ! coefficient 88 DOUBLE PRECISION :: lDab, Sum_XgLoDab, Sum_XgL, lXgN 89 !---------------------------------------------------------------------// 90 91 CALL SET_BINARY_DAB_GAS(Dab, Tg_ref, Pg_ref) 92 93 !!$omp parallel do private(ijk) & 94 !!$omp& schedule(dynamic,chunk_size) 95 DO N = 1, NMAX(0) 96 DO IJK = IJKSTART3, IJKEND3 97 IF (FLUID_AT(IJK)) THEN 98 99 IF (NMAX(0) == 1) THEN 100 ! Only 1 species present in gas phase 101 lDab = Dab(N,N) 102 103 ELSE 104 ! Dilute mixture approximation for multi-component diffusion 105 106 SUM_XgLoDab = ZERO 107 SUM_XgL = ZERO 108 DO L = 1, NMAX(0) 109 IF (L /= N) THEN 110 SUM_XgLoDab = SUM_XgLoDab+& 111 X_g(IJK,L)/Dab(N,L) 112 ! Sum_XgL = 1-XgN 113 SUM_XgL = SUM_XgL + X_g(IJK,L) 114 ENDIF 115 ENDDO 116 ! it is possible for lXgN to evaluate <0 117 lXgN = ONE-SUM_XgL 118 119 IF (lXgN > ZERO_X_GS .AND. & 120 SUM_XgL > ZERO_X_GS) THEN 121 ! i.e. when cell is not only species N 122 ! If this criteria is too strict (i.e. when XgN->1), then this section 123 ! may be evaluated when the other XgN do not carry significant value 124 ! (i.e. are effectively zero). As a result, lDab may be evaluated and 125 ! become unrealistically large. This may happen when happen when 126 ! 1-XgN != sum_XgL. Generally, sum_XgL should equal 1-XgN but they may 127 ! differ due to the numerical evaluation of each Xg. If XgN->1, then 128 ! use both sum_XgL and 1-XgN to determine whether calculations should 129 ! proceed. Given the noted limitation of this approximation then an 130 ! additional criteria should probably be added to only evaluate when 131 ! sum_xgL <0.1. 132 IF (SUM_XgLoDab > ZERO) THEN 133 ! for numerical reasons use Sum_XgL rather than ONE-X_g(IJK,N) 134 lDab = SUM_XgL/SUM_XgLoDab 135 ELSE 136 ! this should not occur... 137 lDab = Dab(N,N) 138 ENDIF 139 ELSE 140 ! Address case when the mass fraction of the Nth species is nearly 1. 141 lDab = Dab(N,N) 142 ENDIF 143 ENDIF 144 145 ! Influence of gas temperature and gas pressure from Fuller relation 146 DIF_G(IJK,N) = ROP_G(IJK)*lDab* & 147 (T_g(IJK)/Tg_ref(N))**1.75* & 148 Pg_ref(N)/UNSCALE_PRESSURE(P_g(IJK)) 149 ELSE 150 DIF_G(IJK,N) = ZERO 151 ENDIF 152 ENDDO 153 ENDDO 154 155 RETURN 156 END SUBROUTINE CALC_DEFAULT_DIF_GAS 157 158 159 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 160 ! C 161 ! Subroutine: SET_BINARY_DAB_GAS C 162 ! Purpose: Set binary diffusion coefficients for all cross species C 163 ! C 164 ! C 165 ! Literature/Document References: C 166 ! Bird, Stewart and lightfoot (1960) C 167 ! Reid, Prausnitz and Poling (1987) C 168 ! C 169 ! C 170 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 171 SUBROUTINE SET_BINARY_DAB_GAS(lDab, lTg_ref, lPg_ref) 172 173 ! Modules 174 !---------------------------------------------------------------------// 175 use physprop, only: NMAX 176 use run, only: units 177 IMPLICIT NONE 178 179 ! Dummy arguments 180 !---------------------------------------------------------------------// 181 ! Binary diffusion coefficient 182 DOUBLE PRECISION, INTENT(OUT) :: lDab(NMAX(0),NMAX(0)) 183 DOUBLE PRECISION, INTENT(OUT) :: lTg_ref(NMAX(0)) 184 DOUBLE PRECISION, INTENT(OUT) :: lPg_ref(NMAX(0)) 185 186 ! Local Variables 187 !---------------------------------------------------------------------// 188 ! Species index 189 INTEGER :: N, L 190 !---------------------------------------------------------------------// 191 192 193 ! Default gas diffusion coefficient 194 ! Bird, Stewart, and Lightfoot (1960) -- CO2--N2 at 298.2 K 195 DO N = 1,NMAX(0) 196 DO L = N,NMAX(0) 197 ! DO L = N+1, NMAX(0) 198 ! IF (N /= L) THEN 199 ! assign the diaganol case for when nmax = 1 or when only species N present 200 ! in given cell. 201 lDab(N,L) = 0.165D0 !cm^2/s 202 lDab(L,N) = lDab(N,L) 203 ! ENDIF 204 ENDDO 205 lTg_ref(N) = 298.2d0 206 lPg_ref(N) = 1.01D6 !dyne 207 ENDDO 208 209 ! Gas diffusion coefficients for 3 species system: SiH4, H2 & N2 210 ! Calculated using relation derived from Chapman and Enskog's theory 211 ! of gases - Reid, Prausnitz and Poling (1987) 212 ! Binary diffusion coefficient SiH4/H2 at 873 K 213 ! lDab(1,2) = 3.78 ! cm^2/s 214 ! lDab(2,1) = lDab(1,2) 215 ! Binary diffusion coefficient SiH4/N2 at 873 K 216 ! lDab(1,3) = 1.02 ! cm^2/s 217 ! lDab(3,1) = lDab(1,3) 218 ! Binary diffusion coefficient H2/N2 at 873 K 219 ! lDab(2,3) = 4.52 ! cm^2/s 220 ! lDab(3,2) = lDab(2,3) 221 ! DO N = 1,NMAX(0) 222 ! lTg_ref(N) = 873.0 223 ! lPg_ref(N) = 1.01e6 ! dyne 224 ! ENDDO 225 226 IF(UNITS == 'SI') THEN 227 DO N = 1,NMAX(0) 228 DO L = N,NMAX(0) 229 ! IF (N /= L) THEN 230 ! assign the diaganol case nmax = 1 or when only species N present 231 ! in given cell. 232 lDab(N,L) = lDab(N,L)*0.0001D0 !m^2/s 233 lDab(L,N) = lDab(N,L) 234 ! ENDIF 235 ENDDO 236 lPg_ref(N) = lPg_ref(N)/10.D0 !Pa 237 ENDDO 238 ENDIF 239 240 RETURN 241 END SUBROUTINE SET_BINARY_DAB_GAS 242