💻 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.
Mass Diffusion — Fick's Law & Stagnant Film
Core Numerical Engine in Fortran 90 • 26 total downloads
mass_diffusion.f90
! =========================================================================
! Source File: mass_diffusion.f90
! =========================================================================
program mass_diffusion
implicit none
! Inputs
integer :: mode ! 1 = Equimolar Counter-Diffusion (ECD), 2 = Diffusion in Stagnant Gas (Stefan)
double precision :: T_C, T_K ! Temperature in C and Kelvin
double precision :: P_kPa, P_Pa ! Pressure in kPa and Pascals
double precision :: D_AB ! Diffusivity in m2/s
double precision :: L ! Film thickness in m
double precision :: P_A1_kPa, P_A1_Pa ! Boundary 1 partial pressure of solute A in kPa and Pa
double precision :: P_A2_kPa, P_A2_Pa ! Boundary 2 partial pressure of solute A in kPa and Pa
double precision :: M_A ! Molecular weight of solute A in g/mol
! Constants
double precision, parameter :: R = 8.314462618d0 ! Universal gas constant, J/(mol.K)
! Physical properties
double precision :: C_total ! Total concentration, mol/m3
double precision :: y_A1, y_A2 ! Mole fractions at boundaries
double precision :: P_B1, P_B2, P_Blm ! Partial pressures of stagnant gas B
double precision :: N_A ! Molar flux of A, mol/m2.s
double precision :: n_A_mass ! Mass flux of A, g/m2.s
double precision :: h_m ! Mass transfer coefficient, m/s
! Local profile variables
integer :: i, n_points
double precision :: x, dx, local_yA, local_CA, local_PA
character(len=50) :: mode_name
integer :: iostat_val
! Read inputs
read(*,*,iostat=iostat_val) mode
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid mode input (must be 1 or 2).'
stop
end if
read(*,*,iostat=iostat_val) T_C
read(*,*,iostat=iostat_val) P_kPa
read(*,*,iostat=iostat_val) D_AB
read(*,*,iostat=iostat_val) L
read(*,*,iostat=iostat_val) P_A1_kPa
read(*,*,iostat=iostat_val) P_A2_kPa
read(*,*,iostat=iostat_val) M_A
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all mass diffusion parameters.'
stop
end if
! Validations
if (mode /= 1 .and. mode /= 2) then
write(*,*) 'ERROR: Mode must be 1 (Equimolar Counter-Diffusion) or 2 (Stagnant Film).'
stop
end if
if (T_C < -273.15d0) then
write(*,*) 'ERROR: Temperature must be above absolute zero.'
stop
end if
if (P_kPa <= 0.0d0) then
write(*,*) 'ERROR: Pressure must be positive.'
stop
end if
if (D_AB <= 0.0d0) then
write(*,*) 'ERROR: Diffusivity must be positive.'
stop
end if
if (L <= 0.0d0) then
write(*,*) 'ERROR: Film thickness must be positive.'
stop
end if
if (P_A1_kPa < 0.0d0 .or. P_A1_kPa > P_kPa) then
write(*,*) 'ERROR: Boundary 1 partial pressure must be between 0 and total pressure.'
stop
end if
if (P_A2_kPa < 0.0d0 .or. P_A2_kPa > P_kPa) then
write(*,*) 'ERROR: Boundary 2 partial pressure must be between 0 and total pressure.'
stop
end if
if (M_A <= 0.0d0) then
write(*,*) 'ERROR: Molecular weight must be positive.'
stop
end if
! Conversions
T_K = T_C + 273.15d0
P_Pa = P_kPa * 1000.0d0
P_A1_Pa = P_A1_kPa * 1000.0d0
P_A2_Pa = P_A2_kPa * 1000.0d0
! Total molar concentration (ideal gas law C = P / RT)
C_total = P_Pa / (R * T_K)
y_A1 = P_A1_Pa / P_Pa
y_A2 = P_A2_Pa / P_Pa
! Core Calculations
if (mode == 1) then
mode_name = "Equimolar Counter-Diffusion"
! Molar flux N_A = D_AB * (C_A1 - C_A2) / L
N_A = D_AB * C_total * (y_A1 - y_A2) / L
h_m = D_AB / L
P_Blm = 0.0d0 ! Not applicable
else
mode_name = "Diffusion through a Stagnant Gas B"
! Boundary partial pressures of stagnant B
P_B1 = P_Pa - P_A1_Pa
P_B2 = P_Pa - P_A2_Pa
! Log-mean partial pressure of B
if (abs(P_B1 - P_B2) < 1.0d-5 * P_Pa) then
P_Blm = (P_B1 + P_B2) / 2.0d0
else
P_Blm = (P_B2 - P_B1) / log(P_B2 / P_B1)
end if
! Molar flux N_A = D_AB * P * (P_A1 - P_A2) / (R * T * L * P_Blm)
N_A = D_AB * P_Pa * (P_A1_Pa - P_A2_Pa) / (R * T_K * L * P_Blm)
h_m = D_AB * P_Pa / (L * P_Blm)
end if
! Convert to mass flux (g/m2.s)
n_A_mass = N_A * M_A
! Output results
write(*,'(A)') '============================================================'
write(*,'(A)') ' MASS DIFFUSION CALCULATION ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Calculation Mode = ', trim(mode_name)
write(*,'(A,I2)') ' Mode Code = ', mode
write(*,'(A,F12.2,A)') ' Temperature = ', T_C, ' deg-C'
write(*,'(A,F12.2,A)') ' Total Pressure = ', P_kPa, ' kPa'
write(*,'(A,ES12.4,A)') ' Binary Diffusivity D_AB = ', D_AB, ' m2/s'
write(*,'(A,F12.4,A)') ' Film Thickness (L) = ', L, ' m'
write(*,'(A,F12.2,A)') ' Molecular Weight M_A = ', M_A, ' g/mol'
write(*,*)
write(*,'(A)') '--- BOUNDARY CONDITIONS -------------------------------------'
write(*,'(A,F12.4,A)') ' Boundary 1 Partial P_A1 = ', P_A1_kPa, ' kPa'
write(*,'(A,F12.4,A)') ' Boundary 2 Partial P_A2 = ', P_A2_kPa, ' kPa'
write(*,'(A,F12.6)') ' Boundary 1 Mole Frac yA1 = ', y_A1
write(*,'(A,F12.6)') ' Boundary 2 Mole Frac yA2 = ', y_A2
write(*,*)
write(*,'(A)') '--- ENGINE RESULTS ------------------------------------------'
write(*,'(A,F12.4,A)') ' Total Concentration C = ', C_total, ' mol/m3'
if (mode == 2) then
write(*,'(A,F12.4,A)') ' Log-Mean Pressure P_Blm = ', P_Blm/1000.0d0, ' kPa'
end if
write(*,'(A,F12.6,A)') ' Mass Transfer Coeff (hm) = ', h_m, ' m/s'
write(*,'(A,ES12.4,A)') ' Molar Flux of A (N_A) = ', N_A, ' mol/m2.s'
write(*,'(A,ES12.4,A)') ' Mass Flux of A (n_A) = ', n_A_mass, ' g/m2.s'
write(*,*)
! ============================================
! LOCAL PROFILE DATA (along path coordinate x)
! ============================================
write(*,'(A)') '--- LOCAL PROFILE ALONG FILM --------------------------------'
write(*,'(A)') ' x [m] y_A [-] C_A [mol/m3] P_A [kPa]'
write(*,'(A)') ' -----------------------------------------------------------'
n_points = 40
dx = L / dble(n_points)
! Include x = 0 initially
write(*,'(F8.4,4X,F10.6,4X,F12.4,4X,F10.4)') 0.0d0, y_A1, y_A1 * C_total, P_A1_kPa
do i = 1, n_points
x = dble(i) * dx
if (mode == 1) then
! Linear profile for ECD
local_yA = y_A1 - (y_A1 - y_A2) * (x / L)
else
! Non-linear exponential-like profile for stagnant gas
if (abs(y_A1 - y_A2) < 1.0d-7) then
local_yA = y_A1
else
local_yA = 1.0d0 - (1.0d0 - y_A1) * ((1.0d0 - y_A2)/(1.0d0 - y_A1))**(x/L)
end if
end if
local_CA = local_yA * C_total
local_PA = local_yA * P_kPa
write(*,'(F8.4,4X,F10.6,4X,F12.4,4X,F10.4)') x, local_yA, local_CA, local_PA
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
if (mode == 1) then
write(*,'(A)') " Fick's First Law for Equimolar Counter-Diffusion (linear gradient)."
else
write(*,'(A)') " Stefan's Diffusion Law (diffusion through stagnant gas B, logarithmic profile)."
end if
write(*,'(A)') ' Total concentration computed from ideal gas equation of state.'
end program mass_diffusion
Solver Description
Models steady 1D binary gas diffusion across a flat film. It supports Equimolar Counter-Diffusion (ECD) where species A and B diffuse in opposite directions at equal rates (linear concentration gradient), and Stefan's diffusion where A evaporates through stagnant component B (logarithmic mole fraction profile). Computes total concentration, log-mean stagnant pressure, mass transfer coefficient, and molar/mass fluxes.
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 mass_diffusion.f90 -o mass_diff
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
mass_diff < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Diffusion Mode (1=Equimolar Counter-Diffusion, 2=Stagnant Gas B)
2
! System Temperature [°C]
25.0
! Total Pressure [kPa]
101.325
! Binary Diffusivity D_AB [m2/s]
2.6e-5
! Film Thickness L [m]
0.1
! Boundary 1 Solute Partial Pressure [kPa]
3.17
! Boundary 2 Solute Partial Pressure [kPa]
0.0
! Solute Molecular Weight [g/mol]
18.015
2
! System Temperature [°C]
25.0
! Total Pressure [kPa]
101.325
! Binary Diffusivity D_AB [m2/s]
2.6e-5
! Film Thickness L [m]
0.1
! Boundary 1 Solute Partial Pressure [kPa]
3.17
! Boundary 2 Solute Partial Pressure [kPa]
0.0
! Solute Molecular Weight [g/mol]
18.015