💻 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.
Transient Diffusion in Solids
Core Numerical Engine in Fortran 90 • 56 total downloads
transient_diffusion_solid.f90
! =========================================================================
! Source File: transient_diffusion_solid.f90
! =========================================================================
program transient_diffusion_solid
implicit none
integer :: geom_type, i, n_points, n_time, iostat_val
double precision :: L, R, Aref, D, hm, C0, Cs, time, rho_s
double precision :: Lc, Volume, Area_s, Fo, Bi, Ccenter, Csurf, Cavg, Mabs, Mfrac
double precision :: x, xstar, theta, Cx, relpos, tprof, Foprof, fracprof, Mprof
character(len=80) :: geom_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) L
read(*,*,iostat=iostat_val) R
read(*,*,iostat=iostat_val) Aref
read(*,*,iostat=iostat_val) D
read(*,*,iostat=iostat_val) hm
read(*,*,iostat=iostat_val) C0
read(*,*,iostat=iostat_val) Cs
read(*,*,iostat=iostat_val) time
read(*,*,iostat=iostat_val) rho_s
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all transient diffusion parameters.'
stop
end if
if (D <= 0.0d0) then
write(*,*) 'ERROR: Diffusivity D must be positive.'
stop
end if
if (time < 0.0d0) then
write(*,*) 'ERROR: Time cannot be negative.'
stop
end if
if (hm < 0.0d0) then
write(*,*) 'ERROR: Surface mass transfer coefficient cannot be negative.'
stop
end if
select case (geom_type)
case (1)
if (L <= 0.0d0) then
write(*,*) 'ERROR: Slab half-thickness L must be positive.'
stop
end if
if (Aref <= 0.0d0) Aref = 1.0d0
geom_name = 'Plane Wall / Slab'
Lc = L
Volume = 2.0d0 * L * Aref
Area_s = 2.0d0 * Aref
case (2)
if (R <= 0.0d0) then
write(*,*) 'ERROR: Cylinder radius R must be positive.'
stop
end if
if (Aref <= 0.0d0) Aref = 1.0d0
geom_name = 'Long Cylinder'
Lc = R
Volume = PI * R**2 * Aref
Area_s = 2.0d0 * PI * R * Aref
case (3)
if (R <= 0.0d0) then
write(*,*) 'ERROR: Sphere radius R must be positive.'
stop
end if
geom_name = 'Sphere'
Lc = R
Volume = 4.0d0/3.0d0 * PI * R**3
Area_s = 4.0d0 * PI * R**2
case default
write(*,*) 'ERROR: Invalid geometry. Use 1 slab, 2 cylinder, or 3 sphere.'
stop
end select
Fo = D * time / (Lc**2)
if (hm <= 0.0d0) then
Bi = 0.0d0
else
Bi = hm * Lc / D
end if
Ccenter = concentration_at(geom_type, 0.0d0, Fo, Bi, C0, Cs)
Csurf = concentration_at(geom_type, 1.0d0, Fo, Bi, C0, Cs)
Cavg = average_concentration(geom_type, Fo, Bi, C0, Cs)
if (abs(Cs - C0) > 1.0d-30) then
Mfrac = (Cavg - C0) / (Cs - C0)
else
Mfrac = 0.0d0
end if
if (Mfrac < 0.0d0 .and. Cs > C0) Mfrac = 0.0d0
if (Mfrac > 1.0d0 .and. Cs > C0) Mfrac = 1.0d0
Mabs = (Cavg - C0) * Volume
write(*,'(A)') '============================================================'
write(*,'(A)') ' TRANSIENT DIFFUSION IN SOLIDS ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Geometry Name = ', trim(geom_name)
write(*,'(A,I2)') ' Geometry Type Code = ', geom_type
write(*,'(A,ES12.4,A)') ' Characteristic Length = ', Lc, ' m'
write(*,'(A,ES12.4,A)') ' Volume = ', Volume, ' m3'
write(*,'(A,ES12.4,A)') ' Exposed Area = ', Area_s, ' m2'
write(*,*)
write(*,'(A)') '--- DIFFUSION INPUTS ----------------------------------------'
write(*,'(A,ES12.4,A)') ' Diffusivity (D) = ', D, ' m2/s'
write(*,'(A,ES12.4,A)') ' Surface h_m = ', hm, ' m/s'
write(*,'(A,ES12.4,A)') ' Initial Concentration = ', C0, ' units'
write(*,'(A,ES12.4,A)') ' Surface Equilibrium Cs = ', Cs, ' units'
write(*,'(A,ES12.4,A)') ' Elapsed Time = ', time, ' s'
write(*,'(A,ES12.4,A)') ' Solid Density = ', rho_s, ' kg/m3'
write(*,*)
write(*,'(A)') '--- DIMENSIONLESS GROUPS ------------------------------------'
write(*,'(A,ES12.4)') ' Fourier Number (Fo) = ', Fo
write(*,'(A,ES12.4)') ' Biot Number (Bi_m) = ', Bi
write(*,*)
write(*,'(A)') '--- TRANSIENT DIFFUSION RESULTS -----------------------------'
write(*,'(A,ES12.4,A)') ' Center Concentration = ', Ccenter, ' units'
write(*,'(A,ES12.4,A)') ' Surface Concentration = ', Csurf, ' units'
write(*,'(A,ES12.4,A)') ' Average Concentration = ', Cavg, ' units'
write(*,'(A,ES12.4,A)') ' Total Absorbed Mass = ', Mabs, ' kg-equivalent'
write(*,'(A,ES12.4)') ' Mass Fraction Progress = ', Mfrac
if (rho_s > 0.0d0) write(*,'(A,ES12.4,A)') ' Avg Concentration ppmw = ', Cavg / rho_s * 1.0d6, ' ppmw'
write(*,*)
write(*,'(A)') '--- CONCENTRATION PROFILE -----------------------------------'
write(*,'(A)') ' x [m] x* theta C(x,t) relative_pos'
write(*,'(A)') ' ----------------------------------------------------------'
n_points = 61
do i = 0, n_points-1
xstar = dble(i) / dble(n_points-1)
x = xstar * Lc
theta = theta_at(geom_type, xstar, Fo, Bi)
Cx = Cs + theta * (C0 - Cs)
relpos = xstar
write(*,'(F10.6,2X,F10.6,2X,ES12.4,2X,ES12.4,2X,F10.6)') x, xstar, theta, Cx, relpos
end do
write(*,*)
write(*,'(A)') '--- MASS UPTAKE PROFILE -------------------------------------'
write(*,'(A)') ' t [s] Fo fraction M_abs [kg]'
write(*,'(A)') ' ----------------------------------------------------------'
n_time = 60
do i = 0, n_time-1
tprof = time * dble(i) / dble(n_time-1)
if (time <= 0.0d0) tprof = 0.0d0
Foprof = D * tprof / (Lc**2)
if (abs(Cs-C0) > 1.0d-30) then
fracprof = (average_concentration(geom_type, Foprof, Bi, C0, Cs) - C0) / (Cs - C0)
else
fracprof = 0.0d0
end if
Mprof = fracprof * (Cs - C0) * Volume
write(*,'(F10.3,2X,ES12.4,2X,ES12.4,2X,ES12.4)') tprof, Foprof, fracprof, Mprof
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' 1D Fickian diffusion in solids with Fourier-series eigenfunction expansions.'
write(*,'(A)') ' Slab finite Bi: lambda tan(lambda) = Bi; Phi = cos(lambda x*).'
write(*,'(A)') ' Cylinder and sphere use standard radial Fourier solutions for prescribed surface concentration when Bi is large.'
write(*,'(A)') ' For finite Bi in cylinder/sphere, a one-term effective approximation is used for surface-resistance correction.'
contains
double precision function theta_at(g, xs, Fo_in, Bi_in)
implicit none
integer, intent(in) :: g
double precision, intent(in) :: xs, Fo_in, Bi_in
integer :: n, maxn
double precision :: lam, coef, sumv, root
double precision, parameter :: PIc = 3.1415926535897932384626433832795d0
if (Fo_in <= 0.0d0) then
theta_at = 1.0d0
return
end if
if (Bi_in <= 1.0d-12) then
theta_at = 1.0d0
return
end if
maxn = 160
sumv = 0.0d0
if (g == 1) then
do n = 1, maxn
call slab_root(n, Bi_in, root)
coef = 4.0d0 * sin(root) / (2.0d0*root + sin(2.0d0*root))
sumv = sumv + coef * cos(root*xs) * exp(-root*root*Fo_in)
end do
else if (g == 2) then
! Large-Bi radial cylinder approximation based on J0 zeros; finite Bi adjusted by Bi/(Bi+1).
do n = 1, maxn
lam = besselj0_zero_approx(n)
coef = 2.0d0 / (lam * besselj1(lam))
sumv = sumv + coef * besselj0(lam*xs) * exp(-lam*lam*Fo_in)
end do
if (Bi_in < 100.0d0) sumv = 1.0d0 - (1.0d0 - sumv) * Bi_in / (Bi_in + 1.0d0)
else
! Sphere, prescribed surface concentration series.
if (xs <= 1.0d-10) then
do n = 1, maxn
lam = dble(n) * PIc
coef = 2.0d0 * (-1.0d0)**(n+1)
sumv = sumv + coef * exp(-lam*lam*Fo_in)
end do
else
do n = 1, maxn
lam = dble(n) * PIc
coef = 2.0d0 * (-1.0d0)**(n+1)
sumv = sumv + coef * sin(lam*xs)/(lam*xs) * exp(-lam*lam*Fo_in)
end do
end if
if (Bi_in < 100.0d0) sumv = 1.0d0 - (1.0d0 - sumv) * Bi_in / (Bi_in + 1.0d0)
end if
theta_at = max(-0.05d0, min(1.05d0, sumv))
end function theta_at
double precision function concentration_at(g, xs, Fo_in, Bi_in, C0_in, Cs_in)
implicit none
integer, intent(in) :: g
double precision, intent(in) :: xs, Fo_in, Bi_in, C0_in, Cs_in
concentration_at = Cs_in + theta_at(g, xs, Fo_in, Bi_in) * (C0_in - Cs_in)
end function concentration_at
double precision function average_concentration(g, Fo_in, Bi_in, C0_in, Cs_in)
implicit none
integer, intent(in) :: g
double precision, intent(in) :: Fo_in, Bi_in
double precision, intent(in) :: C0_in, Cs_in
integer :: j, N
double precision :: xs, w, integ, theta_avg
N = 800
integ = 0.0d0
do j = 0, N
xs = dble(j)/dble(N)
if (j == 0 .or. j == N) then
w = 0.5d0
else
w = 1.0d0
end if
if (g == 1) then
integ = integ + w * theta_at(g, xs, Fo_in, Bi_in)
else if (g == 2) then
integ = integ + w * theta_at(g, xs, Fo_in, Bi_in) * 2.0d0 * xs
else
integ = integ + w * theta_at(g, xs, Fo_in, Bi_in) * 3.0d0 * xs**2
end if
end do
theta_avg = integ / dble(N)
average_concentration = Cs_in + theta_avg * (C0_in - Cs_in)
end function average_concentration
subroutine slab_root(n, Bi_in, root)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: Bi_in
double precision, intent(out) :: root
integer :: it
double precision :: a, b, c, fa, fc, eps
double precision, parameter :: PIc = 3.1415926535897932384626433832795d0
eps = 1.0d-10
a = (dble(n)-1.0d0)*PIc + eps
b = (dble(n)-0.5d0)*PIc - eps
if (Bi_in > 1.0d8) then
root = (dble(n)-0.5d0)*PIc
return
end if
fa = a*tan(a) - Bi_in
do it = 1, 80
c = 0.5d0*(a+b)
fc = c*tan(c) - Bi_in
if (fa*fc <= 0.0d0) then
b = c
else
a = c
fa = fc
end if
end do
root = 0.5d0*(a+b)
end subroutine slab_root
double precision function besselj0(x)
implicit none
double precision, intent(in) :: x
integer :: k
double precision :: term, sumv
if (abs(x) < 12.0d0) then
term = 1.0d0; sumv = 1.0d0
do k = 1, 80
term = term * (-(x*x)/4.0d0) / (dble(k)*dble(k))
sumv = sumv + term
if (abs(term) < 1.0d-14) exit
end do
besselj0 = sumv
else
besselj0 = sqrt(2.0d0/(3.141592653589793d0*x)) * cos(x - 0.25d0*3.141592653589793d0)
end if
end function besselj0
double precision function besselj1(x)
implicit none
double precision, intent(in) :: x
integer :: k
double precision :: term, sumv
if (abs(x) < 12.0d0) then
term = x/2.0d0; sumv = term
do k = 1, 80
term = term * (-(x*x)/4.0d0) / (dble(k)*dble(k+1))
sumv = sumv + term
if (abs(term) < 1.0d-14) exit
end do
besselj1 = sumv
else
besselj1 = sqrt(2.0d0/(3.141592653589793d0*x)) * cos(x - 0.75d0*3.141592653589793d0)
end if
end function besselj1
double precision function besselj0_zero_approx(n)
implicit none
integer, intent(in) :: n
double precision, parameter :: PIc = 3.1415926535897932384626433832795d0
! McMahon approximation; adequate for profile visualization and series convergence.
besselj0_zero_approx = (dble(n)-0.25d0)*PIc + 1.0d0/(8.0d0*(dble(n)-0.25d0)*PIc)
end function besselj0_zero_approx
end program transient_diffusion_solid
Solver Description
Calculate transient mass diffusion profiles, fractional saturation, and total mass transfer over time for standard geometries.
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 transient_diffusion_solid.f90 -o transient_diffusion_solid
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
transient_diffusion_solid < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry (1=Slab, 2=Cylinder, 3=Sphere)
1
! Half-Thickness/Radius L [m]
0.01
! Position x/r [m]
0.0
! Area A [m2]
0.0
! Diffusivity D [m2/s]
1e-9
! Convective Mass Transfer Coeff hm [m/s]
1e-4
! Initial Concentration C0 [mol/m3]
0.0
! Ambient/Surface Concentration Cs [mol/m3]
1.0
! Time t [s]
3600.0
! Solid Density rho_s [kg/m3]
1000.0
1
! Half-Thickness/Radius L [m]
0.01
! Position x/r [m]
0.0
! Area A [m2]
0.0
! Diffusivity D [m2/s]
1e-9
! Convective Mass Transfer Coeff hm [m/s]
1e-4
! Initial Concentration C0 [mol/m3]
0.0
! Ambient/Surface Concentration Cs [mol/m3]
1.0
! Time t [s]
3600.0
! Solid Density rho_s [kg/m3]
1000.0