program binary_diffusivity implicit none integer :: mode, i, n_points, iostat_val double precision :: T_C, T_K, P_atm, MA, MB, sigmaA, sigmaB, epsA, epsB double precision :: VA, mu_cP, phi, rho, D_cm2s, D_m2s, Sc double precision :: sigmaAB, epsAB, Tstar, OmegaD, Tprof_C, Tprof_K, Dprof_cm2s, Dprof_m2s character(len=80) :: method_name read(*,*,iostat=iostat_val) mode if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid model selection input.' stop end if read(*,*,iostat=iostat_val) T_C read(*,*,iostat=iostat_val) P_atm read(*,*,iostat=iostat_val) MA read(*,*,iostat=iostat_val) MB read(*,*,iostat=iostat_val) sigmaA read(*,*,iostat=iostat_val) sigmaB read(*,*,iostat=iostat_val) epsA read(*,*,iostat=iostat_val) epsB read(*,*,iostat=iostat_val) VA read(*,*,iostat=iostat_val) mu_cP read(*,*,iostat=iostat_val) phi read(*,*,iostat=iostat_val) rho if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all binary diffusivity parameters.' stop end if T_K = T_C + 273.15d0 if (T_K <= 0.0d0) then write(*,*) 'ERROR: Temperature must be above absolute zero.' stop end if if (MA <= 0.0d0 .or. MB <= 0.0d0) then write(*,*) 'ERROR: Molecular weights must be positive.' stop end if if (mode == 1) then if (P_atm <= 0.0d0 .or. sigmaA <= 0.0d0 .or. sigmaB <= 0.0d0 .or. epsA <= 0.0d0 .or. epsB <= 0.0d0) then write(*,*) 'ERROR: Gas pressure, sigma, and epsilon/k inputs must be positive.' stop end if method_name = 'Chapman-Enskog gas diffusion' call gas_diffusivity(T_K, P_atm, MA, MB, sigmaA, sigmaB, epsA, epsB, D_cm2s, sigmaAB, epsAB, Tstar, OmegaD) D_m2s = D_cm2s * 1.0d-4 Sc = 0.0d0 else if (mode == 2) then if (VA <= 0.0d0 .or. mu_cP <= 0.0d0 .or. phi <= 0.0d0) then write(*,*) 'ERROR: Liquid molar volume, viscosity, and association factor must be positive.' stop end if method_name = 'Wilke-Chang liquid diffusion' D_cm2s = 7.4d-8 * sqrt(phi * MB) * T_K / (mu_cP * VA**0.6d0) D_m2s = D_cm2s * 1.0d-4 if (rho > 0.0d0) then Sc = (mu_cP * 1.0d-3) / (rho * D_m2s) else Sc = 0.0d0 end if sigmaAB = 0.0d0; epsAB = 0.0d0; Tstar = 0.0d0; OmegaD = 0.0d0 else write(*,*) 'ERROR: Invalid model. Use 1 for gas or 2 for liquid.' stop end if write(*,'(A)') '============================================================' write(*,'(A)') ' BINARY DIFFUSIVITY ESTIMATOR ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A,A)') ' Method = ', trim(method_name) write(*,'(A,F12.4,A)') ' Temperature = ', T_C, ' deg-C' write(*,'(A,F12.4,A)') ' Temperature (K) = ', T_K, ' K' write(*,'(A,F12.6,A)') ' Pressure (atm) = ', P_atm, ' atm' write(*,*) write(*,'(A)') '--- SPECIES / MOLECULAR INPUTS -----------------------------' write(*,'(A,F12.4,A)') ' Molecular Weight A = ', MA, ' g/mol' write(*,'(A,F12.4,A)') ' Molecular Weight B = ', MB, ' g/mol' if (mode == 1) then write(*,'(A,F12.4,A)') ' sigma_A = ', sigmaA, ' Angstrom' write(*,'(A,F12.4,A)') ' sigma_B = ', sigmaB, ' Angstrom' write(*,'(A,F12.4,A)') ' epsilon/k_A = ', epsA, ' K' write(*,'(A,F12.4,A)') ' epsilon/k_B = ', epsB, ' K' write(*,'(A,F12.4,A)') ' sigma_AB = ', sigmaAB, ' Angstrom' write(*,'(A,F12.4,A)') ' epsilon/k_AB = ', epsAB, ' K' write(*,'(A,F12.4)') ' Reduced Temperature T* = ', Tstar write(*,'(A,F12.4)') ' Collision Integral = ', OmegaD else write(*,'(A,F12.4,A)') ' Solute Molar Volume VA = ', VA, ' cm3/mol' write(*,'(A,F12.4,A)') ' Solvent Viscosity = ', mu_cP, ' cP' write(*,'(A,F12.4)') ' Association Factor phi = ', phi write(*,'(A,F12.4,A)') ' Liquid Density = ', rho, ' kg/m3' end if write(*,*) write(*,'(A)') '--- DIFFUSIVITY RESULTS -------------------------------------' write(*,'(A,ES12.4,A)') ' Binary Diffusivity (D_AB) = ', D_m2s, ' m2/s' write(*,'(A,ES12.4,A)') ' D_AB [cm2/s] = ', D_cm2s, ' cm2/s' if (mode == 2 .and. Sc > 0.0d0) write(*,'(A,ES12.4)') ' Schmidt Number (Sc) = ', Sc write(*,*) write(*,'(A)') '--- PARAMETRIC PROFILE --------------------------------------' write(*,'(A)') ' T_C [C] P [atm] D_AB [m2/s] D_AB [cm2/s]' write(*,'(A)') ' ----------------------------------------------------------' n_points = 40 do i = 1, n_points Tprof_C = max(-250.0d0, T_C - 50.0d0 + 100.0d0 * dble(i-1) / dble(n_points-1)) Tprof_K = Tprof_C + 273.15d0 if (mode == 1) then call gas_diffusivity(Tprof_K, P_atm, MA, MB, sigmaA, sigmaB, epsA, epsB, Dprof_cm2s, sigmaAB, epsAB, Tstar, OmegaD) else Dprof_cm2s = 7.4d-8 * sqrt(phi * MB) * Tprof_K / (mu_cP * VA**0.6d0) end if Dprof_m2s = Dprof_cm2s * 1.0d-4 write(*,'(F9.3,2X,F10.5,2X,ES12.4,2X,ES12.4)') Tprof_C, P_atm, Dprof_m2s, Dprof_cm2s end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' if (mode == 1) then write(*,'(A)') ' Chapman-Enskog: D_AB = 0.001858 T^(3/2) sqrt(1/MA+1/MB)/(P sigma_AB^2 Omega_D) [cm2/s]' write(*,'(A)') ' sigma_AB = (sigma_A + sigma_B)/2; epsilon_AB/k = sqrt(epsilon_A/k * epsilon_B/k)' write(*,'(A)') ' Omega_D uses the Neufeld et al. collision-integral fit.' else write(*,'(A)') ' Wilke-Chang: D_AB = 7.4E-8 sqrt(phi MB) T/(mu_B VA^0.6) [cm2/s]' write(*,'(A)') ' Intended for dilute solutes diffusing in liquid solvents.' end if contains subroutine gas_diffusivity(TK, Patm, M1, M2, sig1, sig2, eps1, eps2, Dcm, sigAB, epsAB, Tst, Om) implicit none double precision, intent(in) :: TK, Patm, M1, M2, sig1, sig2, eps1, eps2 double precision, intent(out) :: Dcm, sigAB, epsAB, Tst, Om sigAB = 0.5d0 * (sig1 + sig2) epsAB = sqrt(eps1 * eps2) Tst = TK / epsAB Om = collision_integral(Tst) Dcm = 0.001858d0 * TK**1.5d0 * sqrt(1.0d0/M1 + 1.0d0/M2) / (Patm * sigAB**2 * Om) end subroutine gas_diffusivity double precision function collision_integral(Ts) implicit none double precision, intent(in) :: Ts collision_integral = 1.06036d0 / Ts**0.15610d0 + 0.19300d0 / exp(0.47635d0*Ts) + & 1.03587d0 / exp(1.52996d0*Ts) + 1.76474d0 / exp(3.89411d0*Ts) end function collision_integral end program binary_diffusivity