🔬
Solver Purpose & Physical Scope
Size packed-bed gas absorption columns using Henry's law equilibrium, Kremser theoretical-stage estimate, and the NTU-HTU height method.
📂 Discipline: Masstransfer
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 344 times
📄 Source File:
absorption_packed_bed.f90
📁 calcul/MassTransfer /
absorption_packed_bed.f90
program absorption_packed_bed
implicit none
integer :: i, n_points, iostat_val
double precision :: y_in, y_out_t, x_in, m, G, L, KGa, area, packingFactor
double precision :: LG, Aabs, x_out, y_star_bottom, y_star_top, removal
double precision :: NTU, HTU, Z, Nstage, Gflux, y, x, ystar, drive, ntu_cum, dy, y_hi, y_lo
double precision :: minL, pinch_top, pinch_bottom
read(*,*,iostat=iostat_val) y_in
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid inlet gas composition input.'
stop
end if
read(*,*,iostat=iostat_val) y_out_t
read(*,*,iostat=iostat_val) x_in
read(*,*,iostat=iostat_val) m
read(*,*,iostat=iostat_val) G
read(*,*,iostat=iostat_val) L
read(*,*,iostat=iostat_val) KGa
read(*,*,iostat=iostat_val) area
read(*,*,iostat=iostat_val) packingFactor
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all absorber inputs.'
stop
end if
if (y_in <= 0.0d0 .or. y_in >= 1.0d0 .or. y_out_t < 0.0d0 .or. y_out_t >= y_in) then
write(*,*) 'ERROR: Require 0 <= y_out < y_in < 1.'
stop
end if
if (x_in < 0.0d0 .or. x_in >= 1.0d0 .or. m <= 0.0d0) then
write(*,*) 'ERROR: Liquid inlet composition and Henry slope are invalid.'
stop
end if
if (G <= 0.0d0 .or. L <= 0.0d0 .or. KGa <= 0.0d0 .or. area <= 0.0d0) then
write(*,*) 'ERROR: G, L, KGa, and column area must be positive.'
stop
end if
if (packingFactor <= 0.0d0) packingFactor = 1.0d0
LG = L/G
Aabs = L/(m*G)
x_out = x_in + G/L * (y_in - y_out_t)
y_star_top = m*x_in
y_star_bottom = m*x_out
pinch_top = y_out_t - y_star_top
pinch_bottom = y_in - y_star_bottom
if (pinch_top <= 0.0d0 .or. pinch_bottom <= 0.0d0) then
write(*,*) 'ERROR: Pinch detected. Increase L/G or reduce target removal.'
stop
end if
removal = (y_in - y_out_t)/y_in * 100.0d0
NTU = integrate_ntu(y_in, y_out_t, x_in, m, G, L)
Gflux = G/area
HTU = Gflux / KGa * packingFactor
Z = NTU * HTU
if (abs(Aabs - 1.0d0) < 1.0d-8) then
Nstage = (y_in - y_out_t) / max(1.0d-12, (y_out_t - y_star_top))
else
! Practical Kremser-style estimate using end-point driving force ratio.
Nstage = log((y_in - y_star_bottom)/(y_out_t - y_star_top)) / log(Aabs)
if (Nstage < 0.0d0) Nstage = abs(Nstage)
end if
if (Nstage < 0.0d0 .or. Nstage /= Nstage) Nstage = NTU
minL = m*G*(y_in - y_out_t)/max(1.0d-12, y_in - m*x_in)
write(*,'(A)') '============================================================'
write(*,'(A)') ' ABSORPTION COLUMN DESIGN - PACKED BED ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,ES12.4)') ' Gas Inlet y_in = ', y_in
write(*,'(A,ES12.4)') ' Outlet Gas y_out = ', y_out_t
write(*,'(A,ES12.4)') ' Lean Liquid x_in = ', x_in
write(*,'(A,ES12.4)') ' Henry Slope (m) = ', m
write(*,'(A,ES12.4,A)') ' Gas Molar Flow (G) = ', G, ' mol/s'
write(*,'(A,ES12.4,A)') ' Liquid Molar Flow (L) = ', L, ' mol/s'
write(*,'(A,ES12.4,A)') ' Volumetric K_Ga = ', KGa, ' mol/m3.s'
write(*,'(A,ES12.4,A)') ' Column Area = ', area, ' m2'
write(*,*)
write(*,'(A)') '--- EQUILIBRIUM / OPERATING LINE ---------------------------'
write(*,'(A,ES12.4)') ' L/G Ratio = ', LG
write(*,'(A,ES12.4)') ' Absorption Factor = ', Aabs
write(*,'(A,ES12.4)') ' Outlet Liquid x_out = ', x_out
write(*,'(A,ES12.4)') ' y*_top = ', y_star_top
write(*,'(A,ES12.4)') ' y*_bottom = ', y_star_bottom
write(*,'(A,ES12.4,A)') ' Minimum Liquid Flow = ', minL, ' mol/s'
write(*,*)
write(*,'(A)') '--- DESIGN RESULTS ------------------------------------------'
write(*,'(A,ES12.4)') ' Theoretical Stages = ', Nstage
write(*,'(A,ES12.4)') ' Overall Gas NTU = ', NTU
write(*,'(A,ES12.4,A)') ' Overall Gas HTU = ', HTU, ' m'
write(*,'(A,ES12.4,A)') ' Column Height = ', Z, ' m'
write(*,'(A,ES12.4,A)') ' Removal Efficiency = ', removal, ' percent'
write(*,*)
write(*,'(A)') '--- COLUMN PROFILE ------------------------------------------'
write(*,'(A)') ' z/H y_gas x_liq y_star driving NTU_cum'
write(*,'(A)') ' ---------------------------------------------------------------------------'
n_points = 60
ntu_cum = 0.0d0
do i = 0, n_points-1
y = y_out_t + (y_in - y_out_t)*dble(i)/dble(n_points-1)
x = x_in + G/L*(y - y_out_t)
ystar = m*x
drive = y - ystar
if (i > 0) then
y_hi = y_out_t + (y_in - y_out_t)*dble(i)/dble(n_points-1)
y_lo = y_out_t + (y_in - y_out_t)*dble(i-1)/dble(n_points-1)
dy = y_hi - y_lo
ntu_cum = ntu_cum + dy / max(1.0d-12, drive)
end if
write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
dble(i)/dble(n_points-1), y, x, ystar, drive, ntu_cum
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Henry law: y* = m x.'
write(*,'(A)') ' Operating line: x = x_in + (G/L)(y - y_out).'
write(*,'(A)') ' NTU_OG = integral dy/(y - y*) from y_out to y_in.'
write(*,'(A)') ' HTU_OG = G_molar_flux / K_Ga; Z = NTU_OG * HTU_OG.'
write(*,'(A)') ' Kremser-style stage estimate based on absorption factor A = L/(mG).'
contains
double precision function integrate_ntu(yin, yout, xin, mh, Gf, Lf)
implicit none
double precision, intent(in) :: yin, yout, xin, mh, Gf, Lf
integer :: j, N
double precision :: y, x, yst, f, sumv, dy, w
N = 2000
dy = (yin - yout)/dble(N)
sumv = 0.0d0
do j = 0, N
y = yout + dy*dble(j)
x = xin + Gf/Lf*(y - yout)
yst = mh*x
f = 1.0d0/max(1.0d-12, y - yst)
if (j == 0 .or. j == N) then
w = 0.5d0
else
w = 1.0d0
end if
sumv = sumv + w*f
end do
integrate_ntu = sumv*dy
end function integrate_ntu
end program absorption_packed_bed
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 absorption_packed_bed.f90 -o absorption_packed_bed
2. Execution with input.txt redirection:
absorption_packed_bed < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.10 0.02 0.0 1.25 1.0 2.2 0.08 0.25 1.0
Parameter Description:
Inlet gas solute mole fraction y_in Target outlet gas mole fraction y_out Lean solvent inlet mole fraction x_in Henry slope m = y*/x Gas molar flow G [mol/s] Liquid molar flow L [mol/s] Volumetric gas-side coefficient K_Ga [mol/m3-s] Column cross-sectional area A [m2] Packing safety factor