🔬
Solver Purpose & Physical Scope
Calculate transient mass diffusion profiles, fractional saturation, and total mass transfer over time for standard geometries.
📂 Discipline: Masstransfer
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 314 times
📄 Source File:
transient_diffusion_solid.f90
📁 calcul/MassTransfer /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 transient_diffusion_solid.f90 -o transient_diffusion_solid
2. Execution with input.txt redirection:
transient_diffusion_solid < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.01 0.0 0.0 1e-9 1e-4 0.0 1.0 3600.0 1000.0
Parameter Description:
Geometry (1=Slab, 2=Cylinder, 3=Sphere) Half-Thickness/Radius L [m] Position x/r [m] Area A [m2] Diffusivity D [m2/s] Convective Mass Transfer Coeff hm [m/s] Initial Concentration C0 [mol/m3] Ambient/Surface Concentration Cs [mol/m3] Time t [s] Solid Density rho_s [kg/m3]