🔬
Solver Purpose & Physical Scope
Calculate chemical enhancement factor and reactive mass transfer flux for gas absorption into reacting liquids using Hatta number analysis.
📂 Discipline: Masstransfer
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 313 times
📄 Source File:
gas_absorption_reactive.f90
📁 calcul/MassTransfer /
gas_absorption_reactive.f90
program gas_absorption_reactive
implicit none
integer :: model, i, n_points, iostat_val
double precision :: kL, DA, k1, k2, CB, CAi, CAb, H, pA, Ei, area
double precision :: kp, Ha, E, Nphys, Nreact, ratio, totalAbs, driving
double precision :: Ha_i, E_i, Nphys_i, Nreact_i, ratio_i
character(len=90) :: modelName, regime
read(*,*,iostat=iostat_val) model
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid model input.'
stop
end if
read(*,*,iostat=iostat_val) kL
read(*,*,iostat=iostat_val) DA
read(*,*,iostat=iostat_val) k1
read(*,*,iostat=iostat_val) k2
read(*,*,iostat=iostat_val) CB
read(*,*,iostat=iostat_val) CAi
read(*,*,iostat=iostat_val) CAb
read(*,*,iostat=iostat_val) H
read(*,*,iostat=iostat_val) pA
read(*,*,iostat=iostat_val) Ei
read(*,*,iostat=iostat_val) area
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all reactive absorption inputs.'
stop
end if
if (kL <= 0.0d0 .or. DA <= 0.0d0 .or. Ei < 1.0d0) then
write(*,*) 'ERROR: kL, DA must be positive and Ei must be at least 1.'
stop
end if
if (CAi <= 0.0d0 .and. H > 0.0d0 .and. pA > 0.0d0) CAi = H*pA
if (CAi < CAb) then
write(*,*) 'ERROR: Interfacial concentration must be greater than or equal to bulk concentration.'
stop
end if
if (area < 0.0d0) area = 0.0d0
if (model == 1) then
modelName = 'Second-order A+B with pseudo-first-order approximation'
if (k2 < 0.0d0 .or. CB < 0.0d0) then
write(*,*) 'ERROR: k2 and CB must be non-negative.'
stop
end if
kp = k2*CB
else if (model == 2) then
modelName = 'Direct pseudo-first-order reaction in liquid film'
if (k1 < 0.0d0) then
write(*,*) 'ERROR: k1 must be non-negative.'
stop
end if
kp = k1
else
write(*,*) 'ERROR: Model must be 1 second-order pseudo-first-order or 2 first-order.'
stop
end if
driving = CAi - CAb
Ha = sqrt(kp*DA)/kL
E = enhancement(Ha, Ei)
Nphys = kL*driving
Nreact = E*Nphys
if (Nphys > 0.0d0) then
ratio = Nreact/Nphys
else
ratio = 0.0d0
end if
totalAbs = Nreact*area
if (Ha < 0.3d0) then
regime = 'Slow reaction; physical absorption nearly controls'
else if (Ha < 3.0d0) then
regime = 'Moderate reaction enhancement in liquid film'
else if (E >= 0.95d0*Ei) then
regime = 'Near instantaneous-reaction enhancement limit'
else
regime = 'Fast reaction; strong chemical enhancement'
end if
write(*,'(A)') '============================================================'
write(*,'(A)') ' GAS ABSORPTION - REACTIVE SYSTEMS ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Model = ', trim(modelName)
write(*,'(A,ES12.4,A)') ' Liquid Mass Transfer Coeff (kL) = ', kL, ' m/s'
write(*,'(A,ES12.4,A)') ' Liquid Diffusivity D_A = ', DA, ' m2/s'
write(*,'(A,ES12.4,A)') ' Interfacial C_Ai = ', CAi, ' mol/m3'
write(*,'(A,ES12.4,A)') ' Bulk C_Ab = ', CAb, ' mol/m3'
write(*,'(A,ES12.4)') ' Pseudo-first-order kprime= ', kp
write(*,'(A,ES12.4)') ' Instantaneous Limit (Ei) = ', Ei
write(*,*)
write(*,'(A)') '--- REACTIVE ABSORPTION RESULTS -----------------------------'
write(*,'(A,ES12.4)') ' Hatta Number (Ha) = ', Ha
write(*,'(A,ES12.4)') ' Enhancement Factor (E) = ', E
write(*,'(A,ES12.4,A)') ' Physical Absorption Flux = ', Nphys, ' mol/m2.s'
write(*,'(A,ES12.4,A)') ' Reactive Absorption Flux = ', Nreact, ' mol/m2.s'
write(*,'(A,ES12.4)') ' Flux Ratio = ', ratio
write(*,'(A,ES12.4,A)') ' Total Absorption Rate = ', totalAbs, ' mol/s'
write(*,'(A,A)') ' Regime Assessment = ', trim(regime)
write(*,*)
write(*,'(A)') '--- ENHANCEMENT PROFILE -------------------------------------'
write(*,'(A)') ' Ha E N_phys N_react ratio'
write(*,'(A)') ' ----------------------------------------------------------------'
n_points = 80
do i=0,n_points-1
Ha_i = 0.01d0*(1000.0d0)**(dble(i)/dble(n_points-1))
E_i = enhancement(Ha_i, Ei)
Nphys_i = Nphys
Nreact_i = E_i*Nphys_i
if (Nphys_i > 0.0d0) then
ratio_i = Nreact_i/Nphys_i
else
ratio_i = 0.0d0
end if
write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') Ha_i, E_i, Nphys_i, Nreact_i, ratio_i
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Hatta number: Ha = sqrt(kprime*D_A)/kL.'
write(*,'(A)') ' Pseudo-first-order approximation for A+B: kprime = k2*C_B.'
write(*,'(A)') ' Enhancement factor approximation: E = Ha/tanh(Ha), capped by Ei.'
write(*,'(A)') ' Reactive flux: N_A,react = E*kL*(C_Ai-C_Ab).'
contains
double precision function enhancement(Ha_in, Ei_in)
implicit none
double precision, intent(in) :: Ha_in, Ei_in
double precision :: val
if (Ha_in < 1.0d-8) then
val = 1.0d0
else if (Ha_in > 50.0d0) then
val = Ha_in
else
val = Ha_in/tanh(Ha_in)
end if
if (val > Ei_in) val = Ei_in
if (val < 1.0d0) val = 1.0d0
enhancement = val
end function enhancement
end program gas_absorption_reactive
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 gas_absorption_reactive.f90 -o gas_absorption_reactive
2. Execution with input.txt redirection:
gas_absorption_reactive < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 2.0e-4 1.9e-9 0.0 8500.0 100.0 0.020 0.000 0.034 10000.0 50.0 1.0
Parameter Description:
Reaction Model (1=Second-order A+B, 2=First-order A) Physical Mass Transfer Coeff kL [m/s] Liquid Diffusivity DA [m2/s] Pseudo-first-order rate k1 [1/s] Second-order rate k2 [m3/mol-s] Bulk Reactant B Concentration CB [mol/m3] Interfacial A Concentration CAi [mol/m3] Bulk A Concentration CAb [mol/m3] Henry Coefficient H [mol/m3-Pa] Gas Partial Pressure pA [Pa] Instantaneous Enhancement Limit Ei Interfacial Area [m2]