🔬
Solver Purpose & Physical Scope
Evaluate Thiele modulus, effectiveness factor, and observable Weisz-Prater criterion for diffusion with chemical reaction in catalyst pellets.
📂 Discipline: Masstransfer
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 236 times
📄 Source File:
diffusion_reaction.f90
📁 calcul/MassTransfer /
diffusion_reaction.f90
program diffusion_reaction
implicit none
integer :: geom_type, reaction_order, i, n_points, iostat_val
double precision :: Rpellet, Lslab, De, k, Cs, rhoCat, targetPhi
double precision :: Lc, phi, eta, wp, rateIntrinsic, rateObs, phi_i, eta_slab_i, eta_cyl_i, eta_sph_i, wp_i
character(len=80) :: geomName, regime
read(*,*,iostat=iostat_val) geom_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid geometry input.'
stop
end if
read(*,*,iostat=iostat_val) reaction_order
read(*,*,iostat=iostat_val) Rpellet
read(*,*,iostat=iostat_val) Lslab
read(*,*,iostat=iostat_val) De
read(*,*,iostat=iostat_val) k
read(*,*,iostat=iostat_val) Cs
read(*,*,iostat=iostat_val) rhoCat
read(*,*,iostat=iostat_val) targetPhi
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all diffusion-reaction inputs.'
stop
end if
if (De <= 0.0d0 .or. k < 0.0d0 .or. Cs <= 0.0d0) then
write(*,*) 'ERROR: De and Cs must be positive; k must be non-negative.'
stop
end if
select case(geom_type)
case(1)
if (Lslab <= 0.0d0) then
write(*,*) 'ERROR: Slab half-thickness must be positive.'
stop
end if
geomName = 'Slab / Flat Pellet'
Lc = Lslab
case(2)
if (Rpellet <= 0.0d0) then
write(*,*) 'ERROR: Cylinder radius must be positive.'
stop
end if
geomName = 'Long Cylinder'
Lc = Rpellet
case(3)
if (Rpellet <= 0.0d0) then
write(*,*) 'ERROR: Sphere radius must be positive.'
stop
end if
geomName = 'Sphere'
Lc = Rpellet
case default
write(*,*) 'ERROR: Geometry must be 1 slab, 2 cylinder, or 3 sphere.'
stop
end select
if (reaction_order == 1) then
phi = Lc*sqrt(k/De)
rateIntrinsic = k*Cs
else
phi = Lc*sqrt(max(k,0.0d0)/(De*Cs))
rateIntrinsic = k
end if
eta = effectiveness(geom_type, phi)
rateObs = eta*rateIntrinsic
wp = rateObs*Lc**2/(Cs*De)
if (phi < 0.3d0 .and. wp < 0.3d0) then
regime = 'Kinetic-control likely; internal diffusion small'
else if (phi < 3.0d0 .or. wp < 1.0d0) then
regime = 'Mixed kinetic/diffusion regime'
else
regime = 'Strong internal diffusion limitation likely'
end if
write(*,'(A)') '============================================================'
write(*,'(A)') ' DIFFUSION WITH CHEMICAL REACTION ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Geometry Name = ', trim(geomName)
write(*,'(A,I8)') ' Reaction Order Model = ', reaction_order
write(*,'(A,ES12.4,A)') ' Characteristic Length = ', Lc, ' m'
write(*,'(A,ES12.4,A)') ' Effective Diffusivity = ', De, ' m2/s'
write(*,'(A,ES12.4)') ' Rate Constant k = ', k
write(*,'(A,ES12.4,A)') ' Surface Concentration Cs = ', Cs, ' mol/m3'
write(*,*)
write(*,'(A)') '--- REACTION-DIFFUSION RESULTS -----------------------------'
write(*,'(A,ES12.4)') ' Thiele Modulus (phi) = ', phi
write(*,'(A,ES12.4)') ' Effectiveness Factor (eta) = ', eta
write(*,'(A,ES12.4)') ' Intrinsic Surface Rate = ', rateIntrinsic
write(*,'(A,ES12.4)') ' Observable Rate = ', rateObs
write(*,'(A,ES12.4)') ' Weisz-Prater Criterion = ', wp
write(*,'(A,A)') ' Regime Assessment = ', trim(regime)
write(*,*)
write(*,'(A)') '--- EFFECTIVENESS PROFILE -----------------------------------'
write(*,'(A)') ' phi eta_slab eta_cylinder eta_sphere WP_if_sphere'
write(*,'(A)') ' -------------------------------------------------------------------'
n_points = 80
do i=0,n_points-1
phi_i = 0.01d0 * (1000.0d0)**(dble(i)/dble(n_points-1))
eta_slab_i = effectiveness(1, phi_i)
eta_cyl_i = effectiveness(2, phi_i)
eta_sph_i = effectiveness(3, phi_i)
wp_i = eta_sph_i*phi_i**2
write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
phi_i, eta_slab_i, eta_cyl_i, eta_sph_i, wp_i
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' First-order Thiele modulus: phi = Lc*sqrt(k/De).'
write(*,'(A)') ' Slab eta = tanh(phi)/phi.'
write(*,'(A)') ' Sphere eta = 3/phi^2*(phi*coth(phi)-1).'
write(*,'(A)') ' Cylinder eta approximated by interpolation between slab and sphere behavior.'
write(*,'(A)') ' Observable Weisz-Prater: C_WP = r_obs*Lc^2/(Cs*De).'
contains
double precision function effectiveness(g, ph)
implicit none
integer, intent(in) :: g
double precision, intent(in) :: ph
double precision :: cothv, slab, sph
if (ph < 1.0d-8) then
effectiveness = 1.0d0
return
end if
slab = tanh(ph)/ph
if (ph > 50.0d0) then
cothv = 1.0d0
else
cothv = cosh(ph)/sinh(ph)
end if
sph = 3.0d0/(ph*ph)*(ph*cothv - 1.0d0)
if (g == 1) then
effectiveness = slab
else if (g == 3) then
effectiveness = sph
else
effectiveness = 0.5d0*(slab+sph)
end if
if (effectiveness > 1.0d0) effectiveness = 1.0d0
if (effectiveness < 0.0d0) effectiveness = 0.0d0
end function effectiveness
end program diffusion_reaction
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 diffusion_reaction.f90 -o diffusion_reaction
2. Execution with input.txt redirection:
diffusion_reaction < input.txt
📄 Sample input.txt File Structure
Sample Data:
3 1 0.002 0.002 1.0e-9 0.01 2.0 1200.0 1.0
Parameter Description:
Geometry (1=Slab, 2=Cylinder, 3=Sphere) Reaction Order (0=Zero, 1=First) Pellet Radius R [m] Slab Half-Thickness L [m] Effective Diffusivity De [m2/s] Reaction Rate Constant k Surface Concentration Cs [mol/m3] Catalyst Bulk Density rhoCat [kg/m3] Reference Target Thiele Modulus