🔍 Search & Filter All 178+ Engineering Calculators →

💻 Fortran Source Code Library

We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.

Binary Diffusivity Estimator

Core Numerical Engine in Fortran 90 • 69 total downloads

binary_diffusivity.f90 
! =========================================================================
! Source File: binary_diffusivity.f90
! =========================================================================

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

Solver Description

Estimate binary diffusion coefficient D_AB for gas pairs using Chapman-Enskog and dilute solutes in liquids using Wilke-Chang.

Key Numerical Methods & Architecture

  • Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
  • Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
  • Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.

🛠️ Local Compilation

To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).

Compilation Command:

gfortran -O3 binary_diffusivity.f90 -o binary_diffusivity

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

binary_diffusivity < input.txt

📥 Downloads & Local Files

Preview of the required input file (input.txt):

! Diffusion Model (1=Gas, 2=Liquid)
1
! Temperature T [C]
25.0
! Pressure P [atm] (Gas only)
1.0
! Solute Molecular Weight MA [g/mol]
18.015
! Solvent Molecular Weight MB [g/mol]
28.97
! Solute Lennard-Jones Parameter sigmaA [A]
2.641
! Solvent Lennard-Jones Parameter sigmaB [A]
3.711
! Solute Energy Parameter epsA/k [K]
809.1
! Solvent Energy Parameter epsB/k [K]
78.6
! Solute Molar Volume VA [cm3/mol] (Liquid only)
18.0
! Solvent Viscosity muB [cP] (Liquid only)
1.0
! Association Factor phi (Liquid only)
1.0
! Liquid Density rho [kg/m3] (Liquid only)
1.184
📦 Download Offline Bundle (ZIP) 22 dl 📄 Download input.txt 25 dl 📥 Download binary_diffusivity.f90 22 dl