🔬
Solver Purpose & Physical Scope
Calculate thermal contact resistance and conductance using Mikic-Rohsenow correlations. Evaluate micro-asperity interface gaps with and without Thermal Interface Materials.
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 309 times
📄 Source File:
thermal_contact_resistance.f90
📁 calcul/Conduction /
thermal_contact_resistance.f90
program thermal_contact_resistance
implicit none
! Input variables
integer :: material_pair_idx
real(8) :: k1, k2, Hc
real(8) :: Ra1, Ra2
real(8) :: P
integer :: fluid_tim_idx
real(8) :: k_tim_input, t_tim_input
real(8) :: temp_val
real(8) :: area_val, Q_val
! Computed intermediate variables
real(8) :: ks, sigma, m1, m2, m, Y, hc_s, kg_air, kg_dry, hc_g_dry
real(8) :: h_joint_dry, R_double_prime_dry, R_c_dry, delta_T_dry
real(8) :: kg_tim, hc_g_tim, h_joint_tim, R_double_prime_tim, R_c_tim, delta_T_tim
real(8) :: k_tim_act, t_tim_act
! Read inputs
read *, material_pair_idx
if (material_pair_idx == 7) then
read *, k1
read *, k2
read *, Hc
else if (material_pair_idx == 1) then
! Al-Al
k1 = 180.0d0
k2 = 180.0d0
Hc = 950.0d0
else if (material_pair_idx == 2) then
! Cu-Cu
k1 = 390.0d0
k2 = 390.0d0
Hc = 700.0d0
else if (material_pair_idx == 3) then
! Steel-Al
k1 = 50.0d0
k2 = 180.0d0
Hc = 950.0d0
else if (material_pair_idx == 4) then
! Steel-Steel
k1 = 50.0d0
k2 = 50.0d0
Hc = 1500.0d0
else if (material_pair_idx == 5) then
! SS-SS
k1 = 15.0d0
k2 = 15.0d0
Hc = 2900.0d0
else if (material_pair_idx == 6) then
! Cu-Al
k1 = 390.0d0
k2 = 180.0d0
Hc = 700.0d0
end if
read *, Ra1
read *, Ra2
read *, P
read *, fluid_tim_idx
read *, k_tim_input
read *, t_tim_input
read *, temp_val
read *, area_val
read *, Q_val
! 1. Harmonic thermal conductivity
ks = 2.0d0 * k1 * k2 / (k1 + k2)
! 2. Surface roughness (Ra1, Ra2 in microns)
sigma = sqrt(Ra1**2 + Ra2**2)
! 3. Surface slopes (m1, m2 from Lambert & Fletcher relation)
m1 = 0.076d0 * (Ra1**0.52d0)
m2 = 0.076d0 * (Ra2**0.52d0)
m = sqrt(m1**2 + m2**2)
! 4. Solid contact conductance (hc_s in W/m2.K)
! hc_s = 1.25 * ks * (m / (sigma * 10^-6)) * (P / Hc)^0.95
hc_s = 1.25d0 * ks * (m / (sigma * 1.0d-6)) * ((P / Hc)**0.95d0)
! 5. Mean gap separation (Y in microns)
! Y = 1.53 * sigma * (P / Hc)^-0.097
Y = 1.53d0 * sigma * ((P / Hc)**(-0.097d0))
! 6. Dry contact computations (Reference for comparison, or if selected)
! gas conductivity for air as a function of temperature: k_air = 0.0243 + 7.1e-5 * T
kg_air = 0.0243d0 + 7.1d-5 * temp_val
! Determine gas conductivity for dry reference: Air or Vacuum
if (fluid_tim_idx == 2) then
kg_dry = 0.0d0 ! Vacuum
else
kg_dry = kg_air ! Air (default dry reference)
end if
hc_g_dry = kg_dry / (Y * 1.0d-6)
h_joint_dry = hc_s + hc_g_dry
R_double_prime_dry = 1.0d0 / h_joint_dry
! Area-based resistance (area is in cm2, convert to m2)
R_c_dry = R_double_prime_dry / (area_val * 1.0d-4)
delta_T_dry = Q_val * R_c_dry
! 7. Actual configuration computations (TIM or Dry)
if (fluid_tim_idx == 1) then
! Air (Dry Air)
h_joint_tim = h_joint_dry
R_double_prime_tim = R_double_prime_dry
k_tim_act = kg_air
t_tim_act = Y
else if (fluid_tim_idx == 2) then
! Vacuum (Dry Vacuum)
kg_dry = 0.0d0
h_joint_tim = hc_s
R_double_prime_tim = 1.0d0 / hc_s
k_tim_act = 0.0d0
t_tim_act = Y
else
! Thermal Grease or Pad
k_tim_act = k_tim_input
t_tim_act = t_tim_input
! unified TIM model
if (t_tim_act > Y) then
R_double_prime_tim = ((t_tim_act - Y) * 1.0d-6) / k_tim_act + &
1.0d0 / (hc_s + k_tim_act / (Y * 1.0d-6))
else
R_double_prime_tim = 1.0d0 / (hc_s + k_tim_act / (Y * 1.0d-6))
end if
h_joint_tim = 1.0d0 / R_double_prime_tim
end if
R_c_tim = R_double_prime_tim / (area_val * 1.0d-4)
delta_T_tim = Q_val * R_c_tim
! Print detailed results report
print *, '=================================================='
print *, ' THERMAL CONTACT RESISTANCE CALCULATION REPORT'
print *, '=================================================='
print *, ''
print *, 'INPUT CHARACTERISTICS:'
print *, '--------------------------------------------------'
print '(A,F10.2,A,F10.2,A)', ' Material Conductivities (k1, k2): ', k1, ' / ', k2, ' W/(m.K)'
print '(A,F10.2,A)', ' Harmonic Mean Conductivity (k_s): ', ks, ' W/(m.K)'
print '(A,F10.2,A)', ' Microhardness Softer Solid (H_c): ', Hc, ' MPa'
print '(A,F10.3,A,F10.3,A)', ' Surface Roughnesses (Ra1, Ra2): ', Ra1, ' / ', Ra2, ' um'
print '(A,F10.3,A)', ' Effective Roughness (sigma): ', sigma, ' um'
print '(A,F10.4,A,F10.4,A)', ' Profile Slopes (m1, m2): ', m1, ' / ', m2, ' rad'
print '(A,F10.4,A)', ' Effective Interface Slope (m): ', m, ' rad'
print '(A,F10.3,A)', ' Apparent Contact Pressure (P): ', P, ' MPa'
print '(A,F10.3,A)', ' Force Ratio (P/H_c): ', P/Hc, ' '
print '(A,F10.2,A)', ' Contact Temperature: ', temp_val, ' C'
print '(A,F10.2,A)', ' Contact Area: ', area_val, ' cm2'
print '(A,F10.2,A)', ' Dissipated Heat Power (Q): ', Q_val, ' W'
print *, ''
print *, 'INTERFACE MICRO-GEOMETRY:'
print *, '--------------------------------------------------'
print '(A,F10.4,A)', ' Mean Separation Gap (Y): ', Y, ' um'
print '(A,ES12.4,A)', ' Solid-Solid Spot Conductance (hc_s):', hc_s, ' W/(m2.K)'
print *, ''
print *, 'THERMAL PERFORMANCE COMPARISON (DRY VS TIM):'
print *, '--------------------------------------------------'
print *, '1. DRY CONTACT REFERENCE (AIR INTERFACE)'
print '(A,ES12.4,A)', ' Gas Conductivity (k_air): ', kg_air, ' W/(m.K)'
print '(A,ES12.4,A)', ' Gas Gap Conductance (hc_g_dry): ', hc_g_dry, ' W/(m2.K)'
print '(A,ES12.4,A)', ' Total Dry Conductance (h_c): ', h_joint_dry, ' W/(m2.K)'
print '(A,ES12.4,A)', ' Contact Resistance (R"_c, dry): ', R_double_prime_dry, ' m2.K/W'
print '(A,F12.4,A)', ' Joint Resistance (R_c, dry): ', R_c_dry, ' C/W'
print '(A,F10.2,A)', ' Interface Temp Drop (delta_T): ', delta_T_dry, ' C'
print *, ''
print *, '2. ACTUAL CONTACT CONFIGURATION'
if (fluid_tim_idx == 1) then
print *, ' Mode: Dry Air'
else if (fluid_tim_idx == 2) then
print *, ' Mode: Dry Vacuum'
else if (fluid_tim_idx == 3) then
print *, ' Mode: Thermal Grease'
print '(A,F10.3,A)', ' TIM Conductivity (k_TIM): ', k_tim_act, ' W/(m.K)'
print '(A,F10.2,A)', ' TIM Bond Thickness (t_TIM): ', t_tim_act, ' um'
else if (fluid_tim_idx == 4) then
print *, ' Mode: Thermal Pad'
print '(A,F10.3,A)', ' TIM Conductivity (k_TIM): ', k_tim_act, ' W/(m.K)'
print '(A,F10.2,A)', ' TIM Thickness (t_TIM): ', t_tim_act, ' um'
end if
print '(A,ES12.4,A)', ' Total Active Conductance (h_c): ', h_joint_tim, ' W/(m2.K)'
print '(A,ES12.4,A)', ' Contact Resistance (R"_c, tim): ', R_double_prime_tim, ' m2.K/W'
print '(A,F12.4,A)', ' Joint Resistance (R_c, tim): ', R_c_tim, ' C/W'
print '(A,F10.2,A)', ' Interface Temp Drop (delta_T): ', delta_T_tim, ' C'
print *, ''
print *, 'SUMMARY COMPARISON INDICATORS:'
print *, '--------------------------------------------------'
print '(A,F10.2,A)', ' TIM Resistance Reduction Factor: ', R_double_prime_dry / R_double_prime_tim, ' x'
print '(A,F10.2,A)', ' TIM Temperature Drop Reduction: ', (delta_T_dry - delta_T_tim), ' C'
print *, '=================================================='
end program thermal_contact_resistance
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 thermal_contact_resistance.f90 -o thermal_contact_calc
2. Execution with input.txt redirection:
thermal_contact_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 1.5 2.0 2.0 3 2.5 30.0 60.0 10.0 100.0
Parameter Description:
Material pair (1=Al-Al, 2=Cu-Cu, 3=Steel-Steel, etc., 7=Custom) Roughness Ra1 [µm] Roughness Ra2 [µm] Contact pressure P [MPa] TIM fluid (1=Air, 2=Helium, 3=Grease, etc.) TIM thermal conductivity k [W/m-K] TIM thickness t [µm] Contact temperature [°C] Joint area [cm2] Heat transfer rate Q [W]