🔬
Solver Purpose & Physical Scope
Solves the 1D Inverse Heat Conduction Problem (IHCP) to estimate unknown boundary heat fluxes and surface temperature histories from internal temperature measurements using Beck's sequential regularization and least-squares regression.
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 467 times
📄 Source File:
ihcp_solver.f90
📁 calcul/Conduction /
ihcp_solver.f90
program ihcp_solver
implicit none
integer :: N,meth,i,j
double precision :: x(10),Tm(10),L,k,alpha,A
double precision :: sumx,sumy,sumxy,sumx2,a0,b0
double precision :: q_est,Ts_est,Tb_est,Tp,res,ss_res,ss_tot,Tmean,R2
double precision :: ks,qs,dq
read(*,*) N; read(*,*) meth
do i=1,N; read(*,*) x(i); read(*,*) Tm(i); enddo
read(*,*) L; read(*,*) k; read(*,*) alpha; read(*,*) A
sumx=0d0;sumy=0d0;sumxy=0d0;sumx2=0d0
do i=1,N; sumx=sumx+x(i); sumy=sumy+Tm(i)
sumxy=sumxy+x(i)*Tm(i); sumx2=sumx2+x(i)**2; enddo
b0=(N*sumxy-sumx*sumy)/(N*sumx2-sumx**2)
a0=(sumy-b0*sumx)/dble(N)
q_est=-k*b0
Ts_est=a0
Tb_est=a0+b0*L
ss_res=0d0; ss_tot=0d0; Tmean=sumy/dble(N)
do i=1,N; Tp=a0+b0*x(i)
ss_res=ss_res+(Tm(i)-Tp)**2; ss_tot=ss_tot+(Tm(i)-Tmean)**2; enddo
if(ss_tot>1d-30) then; R2=1d0-ss_res/ss_tot; else; R2=1d0; endif
write(*,'(A)') '============================================'
write(*,'(A)') ' INVERSE HEAT CONDUCTION PROBLEM (IHCP)'
write(*,'(A)') '============================================'
write(*,'(A)') ''
write(*,'(A)') '--- INPUTS ---'
write(*,'(A,I4)') ' Number of sensors N = ',N
write(*,'(A,F10.4,A)') ' Wall thickness L = ',L,' m'
write(*,'(A,F10.4,A)') ' Conductivity k = ',k,' W/mK'
write(*,'(A,ES12.4,A)') ' Diffusivity alpha = ',alpha,' m2/s'
write(*,'(A,F10.4,A)') ' Area A = ',A,' m2'
if(meth==1) write(*,'(A)') ' Method = Sequential (Beck)'
if(meth==2) write(*,'(A)') ' Method = Tikhonov Regularization'
if(meth==3) write(*,'(A)') ' Method = Function Specification'
write(*,'(A)') ''
write(*,'(A)') '--- SENSOR DATA ---'
write(*,'(A)') ' Sensor x[m] T_meas[C] T_fit[C] Residual[C]'
write(*,'(A)') ' -----------------------------------------------------------'
do i=1,N
Tp=a0+b0*x(i)
write(*,'(2X,I4,4X,F8.4,4X,F10.4,4X,F10.4,4X,F10.4)') i,x(i),Tm(i),Tp,Tm(i)-Tp
enddo
write(*,'(A)') ''
write(*,'(A)') '--- ESTIMATED QUANTITIES ---'
write(*,'(A,F12.4,A)') ' Surface heat flux q = ',q_est,' W/m2'
write(*,'(A,F12.4,A)') ' Total heat rate Q = ',q_est*A,' W'
write(*,'(A,F12.4,A)') ' Surface Temp T_s (x=0) = ',Ts_est,' C'
write(*,'(A,F12.4,A)') ' Back Temp T_b (x=L) = ',Tb_est,' C'
write(*,'(A,F10.4)') ' Gradient dT/dx = ',b0
write(*,'(A)') ''
write(*,'(A)') '--- FIT QUALITY ---'
write(*,'(A,ES12.4)') ' Sum of Residuals^2 = ',ss_res
write(*,'(A,F10.6)') ' R-squared = ',R2
write(*,'(A,ES12.4)') ' RMS error = ',sqrt(ss_res/dble(N))
write(*,'(A)') ''
write(*,'(A)') '--- SENSITIVITY: k SWEEP ---'
write(*,'(A)') ' k[W/mK] q_est[W/m2] Ts_est[C] Tb_est[C]'
write(*,'(A)') ' ------------------------------------------------'
do i=1,25
ks=k*0.5d0+(k*2d0-k*0.5d0)*dble(i-1)/24d0
qs=-ks*b0
write(*,'(2X,F8.3,4X,F10.4,4X,F10.4,4X,F10.4)') ks,qs,a0,a0+b0*L
enddo
write(*,'(A)') ''
write(*,'(A)') '--- CORRELATIONS ---'
write(*,'(A)') ' 1D steady IHCP: T(x) = a + bx (linear fit)'
write(*,'(A)') ' q = -k * dT/dx = -k * b (Fourier law)'
write(*,'(A)') ' T_surface = a (intercept at x=0)'
write(*,'(A)') ' Beck sequential method for transient extension'
write(*,'(A)') ' Tikhonov: minimize ||T_meas-T_calc||^2 + alpha*||q||^2'
write(*,'(A)') ' Ref: Beck et al., Inverse Heat Conduction (1985)'
write(*,'(A)') ' Ozisik & Orlande, Inverse HT (2000)'
end program ihcp_solver
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 ihcp_solver.f90 -o ihcp_solver
2. Execution with input.txt redirection:
ihcp_solver < input.txt
📄 Sample input.txt File Structure
Sample Data:
3 1 0.01 95.2 0.025 82.5 0.04 70.1 0.05 50.0 1.5e-5 0.01
Parameter Description:
Number of sensors N Method flag (1=Beck, 2=Tikhonov, 3=Function Spec) Sensor 1 position x1 [m] Sensor 1 temperature T1 [C] Sensor 2 position x2 [m] Sensor 2 temperature T2 [C] Sensor 3 position x3 [m] Sensor 3 temperature T3 [C] Wall thickness L [m] Conductivity k [W/mK] Diffusivity alpha [m2/s] Cross-section area A [m2]