💻 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.
Porous Media Effective Conductivity
Core Numerical Engine in Fortran 90 • 28 total downloads
porous_media_cond.f90
! =========================================================================
! Source File: porous_media_cond.f90
! =========================================================================
program porous_media_cond
implicit none
integer :: i, j
double precision :: ks,kf,eps,Th,Tc,Lth,A
double precision :: dT,kpar,kser,kHSu,kHSl,kEMT,kgeo,keff
double precision :: Q,Rth,es,f1,f2,df,dk
read(*,*) ks; read(*,*) kf; read(*,*) eps
read(*,*) Th; read(*,*) Tc; read(*,*) Lth; read(*,*) A
dT=abs(Th-Tc)
kpar=eps*kf+(1d0-eps)*ks
kser=1d0/(eps/kf+(1d0-eps)/ks)
if(abs(kf-ks)>1d-12) then
kHSu=ks+eps/(1d0/(kf-ks)+(1d0-eps)/(3d0*ks))
kHSl=kf+(1d0-eps)/(1d0/(ks-kf)+eps/(3d0*kf))
else; kHSu=ks; kHSl=ks; endif
kgeo=ks**(1d0-eps)*kf**eps
kEMT=(kpar+kser)/2d0
do i=1,200
f1=(1d0-eps)*(ks-kEMT)/(ks+2d0*kEMT)+eps*(kf-kEMT)/(kf+2d0*kEMT)
df=-(1d0-eps)*3d0*ks/(ks+2d0*kEMT)**2-eps*3d0*kf/(kf+2d0*kEMT)**2
if(abs(df)<1d-30) exit
dk=-f1/df; kEMT=kEMT+dk
if(kEMT<1d-10) kEMT=1d-10
if(abs(dk)<1d-12) exit
enddo
write(*,'(A)') '============================================'
write(*,'(A)') ' POROUS MEDIA EFFECTIVE CONDUCTIVITY'
write(*,'(A)') '============================================'
write(*,'(A)') ''
write(*,'(A)') '--- INPUTS ---'
write(*,'(A,F12.4,A)') ' k_solid = ',ks,' W/mK'
write(*,'(A,F12.6,A)') ' k_fluid = ',kf,' W/mK'
write(*,'(A,F10.4)') ' Porosity epsilon = ',eps
write(*,'(A,F10.2,A)') ' T_hot = ',Th,' C'
write(*,'(A,F10.2,A)') ' T_cold = ',Tc,' C'
write(*,'(A,F10.4,A)') ' Thickness L = ',Lth,' m'
write(*,'(A,F10.4,A)') ' Area A = ',A,' m2'
write(*,'(A,F10.2)') ' k_solid/k_fluid ratio = ',ks/kf
write(*,'(A)') ''
write(*,'(A)') '--- MODEL RESULTS ---'
write(*,'(A)') ' Model k_eff[W/mK] Q[W] R_th[K/W]'
write(*,'(A)') ' -------------------------------------------------------'
Q=kpar*A*dT/Lth; Rth=Lth/(kpar*A)
write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Parallel (Voigt) ',kpar,Q,Rth
Q=kser*A*dT/Lth; Rth=Lth/(kser*A)
write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Series (Reuss) ',kser,Q,Rth
Q=kHSu*A*dT/Lth; Rth=Lth/(kHSu*A)
write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Hashin-Shtrik Up ',kHSu,Q,Rth
Q=kHSl*A*dT/Lth; Rth=Lth/(kHSl*A)
write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Hashin-Shtrik Lo ',kHSl,Q,Rth
Q=kEMT*A*dT/Lth; Rth=Lth/(kEMT*A)
write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'EMT (Bruggeman) ',kEMT,Q,Rth
Q=kgeo*A*dT/Lth; Rth=Lth/(kgeo*A)
write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Geometric Mean ',kgeo,Q,Rth
write(*,'(A)') ''
write(*,'(A)') '--- POROSITY SWEEP ---'
write(*,'(A)') ' eps k_par k_ser k_HSu k_HSl k_EMT k_geo'
write(*,'(A)') ' ---------------------------------------------------------------------'
do i=1,25
es=0.01d0+0.98d0*dble(i-1)/24d0
kpar=es*kf+(1d0-es)*ks
kser=1d0/(es/kf+(1d0-es)/ks)
if(abs(kf-ks)>1d-12) then
kHSu=ks+es/(1d0/(kf-ks)+(1d0-es)/(3d0*ks))
kHSl=kf+(1d0-es)/(1d0/(ks-kf)+es/(3d0*kf))
else; kHSu=ks; kHSl=ks; endif
kgeo=ks**(1d0-es)*kf**es
keff=(kpar+kser)/2d0
do j=1,200
f1=(1d0-es)*(ks-keff)/(ks+2d0*keff)+es*(kf-keff)/(kf+2d0*keff)
f2=-(1d0-es)*3d0*ks/(ks+2d0*keff)**2-es*3d0*kf/(kf+2d0*keff)**2
if(abs(f2)<1d-30) exit; keff=keff-f1/f2
if(keff<1d-10) keff=1d-10; if(abs(f1)<1d-12) exit
enddo
write(*,'(2X,F6.3,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5)') es,es*kf+(1d0-es)*ks,kser,kHSu,kHSl,keff,kgeo
enddo
write(*,'(A)') ''
write(*,'(A)') '--- CORRELATIONS ---'
write(*,'(A)') ' Parallel (Voigt): k=eps*kf+(1-eps)*ks (upper bound)'
write(*,'(A)') ' Series (Reuss): 1/k=eps/kf+(1-eps)/ks (lower bound)'
write(*,'(A)') ' Hashin-Shtrikman: tightest bounds for isotropic composites'
write(*,'(A)') ' EMT/Bruggeman: self-consistent effective medium'
write(*,'(A)') ' Geometric: k=ks^(1-eps)*kf^eps'
write(*,'(A)') ' Ref: Kaviany, Heat Transfer in Porous Media, Ch.3'
write(*,'(A)') ' Incropera Ch.3, Nield & Bejan Ch.2'
end program porous_media_cond
Solver Description
Calculates effective thermal conductivity ($k_{eff}$) of multi-phase porous media. Computes Parallel (Voigt upper limit), Series (Reuss lower limit), Hashin-Shtrikman isotropic bounds, self-consistent EMT/Bruggeman formulation, and Geometric Mean model.
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 porous_media_cond.f90 -o porous_media_cond
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
porous_media_cond < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Solid conductivity k_solid [W/mK]
1.4
! Fluid conductivity k_fluid [W/mK]
0.026
! Porosity (void fraction) epsilon
0.4
! Hot temperature T_hot [C]
100.0
! Cold temperature T_cold [C]
25.0
! Sample thickness L [m]
0.05
! Cross-section area A [m2]
0.01
1.4
! Fluid conductivity k_fluid [W/mK]
0.026
! Porosity (void fraction) epsilon
0.4
! Hot temperature T_hot [C]
100.0
! Cold temperature T_cold [C]
25.0
! Sample thickness L [m]
0.05
! Cross-section area A [m2]
0.01