🔬
Solver Purpose & Physical Scope
Calculate convective mass transfer coefficients, Sherwood numbers, Schmidt numbers, and mass transfer fluxes for various geometries and flow regimes.
📂 Discipline: Masstransfer
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 493 times
📄 Source File:
convective_mass_transfer.f90
📁 calcul/MassTransfer /
convective_mass_transfer.f90
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 convective_mass_transfer.f90 -o convective_mass_transfer
2. Execution with input.txt redirection:
convective_mass_transfer < input.txt
📄 Sample input.txt File Structure
Sample Data:
2 1 1.0 0.05 0.5 2.0 1.2 1.8e-5 2.5e-5 0.0 1.0 28.97 5.0
Parameter Description:
Geometry (1=Flat Plate, 2=Circular Pipe, 3=Sphere, 4=Single Cylinder) Flow Regime (1=Laminar, 2=Turbulent, 3=Mixed) Characteristic Length L [m] Pipe Inner Diameter D [m] Plate Width W [m] Free-stream/Mean Velocity U [m/s] Fluid Density rho [kg/m3] Fluid Dynamic Viscosity mu [Pa-s] Binary Diffusivity Dab [m2/s] Bulk Solute Concentration C_bulk [mol/m3] Surface Solute Concentration C_surface [mol/m3] Solute Molecular Weight [g/mol] Pipe / Plate Length [m]