program convective_mass_transfer implicit none integer :: geom_type, flow_type, i, n_points, iostat_val double precision :: L, Dpipe, width, velocity, rho, mu, Dab, CA_bulk, CA_surface, MWA, lengthPipe double precision :: Lc, area, nu, Re, Sc, Sh, km, NA, mass_flux, total_mdot, delta_m, dC double precision :: x, Rex, Shx, kmx, NAx, filmx, localLc character(len=80) :: geom_name, regime_name, corr_name double precision, parameter :: PI = 3.1415926535897932384626433832795d0 read(*,*,iostat=iostat_val) geom_type if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid geometry type input.' stop end if read(*,*,iostat=iostat_val) flow_type read(*,*,iostat=iostat_val) L read(*,*,iostat=iostat_val) Dpipe read(*,*,iostat=iostat_val) width read(*,*,iostat=iostat_val) velocity read(*,*,iostat=iostat_val) rho read(*,*,iostat=iostat_val) mu read(*,*,iostat=iostat_val) Dab read(*,*,iostat=iostat_val) CA_bulk read(*,*,iostat=iostat_val) CA_surface read(*,*,iostat=iostat_val) MWA read(*,*,iostat=iostat_val) lengthPipe if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all convective mass transfer parameters.' stop end if if (rho <= 0.0d0 .or. mu <= 0.0d0 .or. Dab <= 0.0d0 .or. MWA <= 0.0d0) then write(*,*) 'ERROR: Fluid properties, diffusivity, and molecular weight must be positive.' stop end if if (velocity < 0.0d0) then write(*,*) 'ERROR: Velocity cannot be negative.' stop end if select case (geom_type) case (1) if (L <= 0.0d0) then write(*,*) 'ERROR: Plate length L must be positive.' stop end if if (width <= 0.0d0) width = 1.0d0 geom_name = 'External Flat Plate' Lc = L area = L * width case (2) if (Dpipe <= 0.0d0) then write(*,*) 'ERROR: Pipe diameter must be positive.' stop end if if (lengthPipe <= 0.0d0) lengthPipe = 1.0d0 geom_name = 'Circular Pipe Internal Flow' Lc = Dpipe area = PI * Dpipe * lengthPipe case default write(*,*) 'ERROR: Invalid geometry. Use 1 flat plate or 2 pipe.' stop end select nu = mu / rho Re = rho * velocity * Lc / mu Sc = mu / (rho * Dab) dC = CA_surface - CA_bulk call sherwood_average(geom_type, flow_type, Re, Sc, Sh, regime_name, corr_name) km = Sh * Dab / Lc NA = km * dC mass_flux = NA * MWA / 1000.0d0 total_mdot = mass_flux * area if (km > 0.0d0) then delta_m = Dab / km else delta_m = 0.0d0 end if write(*,'(A)') '============================================================' write(*,'(A)') ' CONVECTIVE MASS TRANSFER ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A,A)') ' Geometry Name = ', trim(geom_name) write(*,'(A,I2)') ' Geometry Type Code = ', geom_type write(*,'(A,A)') ' Flow Regime = ', trim(regime_name) write(*,'(A,A)') ' Correlation = ', trim(corr_name) write(*,'(A,ES12.4,A)') ' Characteristic Length = ', Lc, ' m' write(*,'(A,ES12.4,A)') ' Transfer Area = ', area, ' m2' write(*,*) write(*,'(A)') '--- FLUID / DIFFUSION PROPERTIES ---------------------------' write(*,'(A,ES12.4,A)') ' Density (rho) = ', rho, ' kg/m3' write(*,'(A,ES12.4,A)') ' Dynamic Viscosity (mu) = ', mu, ' Pa.s' write(*,'(A,ES12.4,A)') ' Kinematic Viscosity = ', nu, ' m2/s' write(*,'(A,ES12.4,A)') ' Binary Diffusivity = ', Dab, ' m2/s' write(*,'(A,ES12.4,A)') ' Velocity = ', velocity, ' m/s' write(*,'(A,ES12.4,A)') ' Species MW_A = ', MWA, ' g/mol' write(*,*) write(*,'(A)') '--- DIMENSIONLESS GROUPS ------------------------------------' write(*,'(A,ES12.4)') ' Reynolds Number (Re) = ', Re write(*,'(A,ES12.4)') ' Schmidt Number (Sc) = ', Sc write(*,'(A,ES12.4)') ' Sherwood Number (Sh) = ', Sh write(*,*) write(*,'(A)') '--- MASS TRANSFER RESULTS -----------------------------------' write(*,'(A,ES12.4,A)') ' Mass Transfer Coeff (k_m) = ', km, ' m/s' write(*,'(A,ES12.4,A)') ' Molar Flux (N_A) = ', NA, ' mol/m2.s' write(*,'(A,ES12.4,A)') ' Mass Flux = ', mass_flux, ' kg/m2.s' write(*,'(A,ES12.4,A)') ' Total Mass Rate = ', total_mdot, ' kg/s' write(*,'(A,ES12.4,A)') ' Film Thickness = ', delta_m, ' m' write(*,'(A,ES12.4,A)') ' Concentration Difference = ', dC, ' mol/m3' write(*,*) write(*,'(A)') '--- LOCAL MASS TRANSFER PROFILE -----------------------------' write(*,'(A)') ' x [m] Re_x Sh_x k_x [m/s] N_A [mol/m2.s] film [m]' write(*,'(A)') ' ---------------------------------------------------------------------------' n_points = 50 do i = 1, n_points if (geom_type == 1) then x = L * dble(i) / dble(n_points) localLc = max(x, 1.0d-12) Rex = rho * velocity * localLc / mu call sherwood_local_plate(flow_type, Rex, Sc, Shx) kmx = Shx * Dab / localLc else x = lengthPipe * dble(i) / dble(n_points) localLc = Dpipe Rex = Re call sherwood_average(geom_type, flow_type, Re, Sc, Shx, regime_name, corr_name) ! Simple entrance enhancement fading with x/D for visualization. if (Dpipe > 0.0d0) Shx = Shx * (1.0d0 + 0.25d0 * exp(-x/(10.0d0*Dpipe))) kmx = Shx * Dab / Dpipe end if NAx = kmx * dC if (kmx > 0.0d0) then filmx = Dab / kmx else filmx = 0.0d0 end if write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') x, Rex, Shx, kmx, NAx, filmx end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' write(*,'(A)') ' Re = rho u Lc / mu; Sc = mu/(rho D_AB); k_m = Sh D_AB/Lc.' write(*,'(A)') ' Flat plate laminar average: Sh_L = 0.664 Re_L^0.5 Sc^(1/3).' write(*,'(A)') ' Flat plate turbulent average: Sh_L = (0.037 Re_L^0.8 - 871) Sc^(1/3).' write(*,'(A)') ' Pipe turbulent Dittus-Boelter mass analogy: Sh = 0.023 Re^0.8 Sc^n.' write(*,'(A)') ' Film thickness: delta_m = D_AB/k_m.' contains subroutine sherwood_average(g, mode, Re_in, Sc_in, Sh_out, regime, corr) implicit none integer, intent(in) :: g, mode double precision, intent(in) :: Re_in, Sc_in double precision, intent(out) :: Sh_out character(len=80), intent(out) :: regime, corr double precision :: n_exp n_exp = 1.0d0/3.0d0 if (g == 1) then if (mode == 3 .or. (mode == 1 .and. Re_in < 5.0d5)) then regime = 'Laminar external boundary layer' corr = 'Flat plate average Sh = 0.664 Re^0.5 Sc^(1/3)' Sh_out = 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0) else regime = 'Turbulent / transition external boundary layer' corr = 'Flat plate average Sh = (0.037 Re^0.8 - 871) Sc^(1/3)' Sh_out = (0.037d0 * max(Re_in,0.0d0)**0.8d0 - 871.0d0) * Sc_in**(1.0d0/3.0d0) if (Sh_out < 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0)) & Sh_out = 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0) end if else if (mode == 3 .or. (mode == 1 .and. Re_in < 2300.0d0)) then regime = 'Laminar internal flow' corr = 'Pipe fully developed laminar Sh = 3.66' Sh_out = 3.66d0 else regime = 'Turbulent internal flow' corr = 'Dittus-Boelter mass analogy Sh = 0.023 Re^0.8 Sc^0.33' Sh_out = 0.023d0 * max(Re_in,0.0d0)**0.8d0 * Sc_in**n_exp end if end if if (Sh_out < 0.0d0) Sh_out = 0.0d0 end subroutine sherwood_average subroutine sherwood_local_plate(mode, Rex, Sc_in, Shx) implicit none integer, intent(in) :: mode double precision, intent(in) :: Rex, Sc_in double precision, intent(out) :: Shx if (mode == 3 .or. (mode == 1 .and. Rex < 5.0d5)) then Shx = 0.332d0 * sqrt(max(Rex,0.0d0)) * Sc_in**(1.0d0/3.0d0) else Shx = 0.0296d0 * max(Rex,0.0d0)**0.8d0 * Sc_in**(1.0d0/3.0d0) end if if (Shx < 0.0d0) Shx = 0.0d0 end subroutine sherwood_local_plate end program convective_mass_transfer