💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Diffusion with Chemical Reaction
Core Numerical Engine in Fortran 90 • 52 total downloads
diffusion_reaction.f90
! =========================================================================
! Source File: 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
Solver Description
Evaluate Thiele modulus, effectiveness factor, and observable Weisz-Prater criterion for diffusion with chemical reaction in catalyst pellets.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 diffusion_reaction.f90 -o diffusion_reaction
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
diffusion_reaction < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry (1=Slab, 2=Cylinder, 3=Sphere)
3
! Reaction Order (0=Zero, 1=First)
1
! Pellet Radius R [m]
0.002
! Slab Half-Thickness L [m]
0.002
! Effective Diffusivity De [m2/s]
1.0e-9
! Reaction Rate Constant k
0.01
! Surface Concentration Cs [mol/m3]
2.0
! Catalyst Bulk Density rhoCat [kg/m3]
1200.0
! Reference Target Thiele Modulus
1.0
3
! Reaction Order (0=Zero, 1=First)
1
! Pellet Radius R [m]
0.002
! Slab Half-Thickness L [m]
0.002
! Effective Diffusivity De [m2/s]
1.0e-9
! Reaction Rate Constant k
0.01
! Surface Concentration Cs [mol/m3]
2.0
! Catalyst Bulk Density rhoCat [kg/m3]
1200.0
! Reference Target Thiele Modulus
1.0