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