🔬
Solver Purpose & Physical Scope
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.
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 464 times
📄 Source File:
porous_media_cond.f90
📁 calcul/Conduction /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 porous_media_cond.f90 -o porous_media_cond
2. Execution with input.txt redirection:
porous_media_cond < input.txt
📄 Sample input.txt File Structure
Sample Data:
1.4 0.026 0.4 100.0 25.0 0.05 0.01
Parameter Description:
Solid conductivity k_solid [W/mK] Fluid conductivity k_fluid [W/mK] Porosity (void fraction) epsilon Hot temperature T_hot [C] Cold temperature T_cold [C] Sample thickness L [m] Cross-section area A [m2]