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