🔬
Solver Purpose & Physical Scope
Calculate conduction shape factors (S), thermal resistance, and multi-dimensional heat flow rates for buried pipes, eccentric cylinders, floor heating arrays, and wall edges.
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 553 times
📄 Source File:
conduction_shape_factors.f90
📁 calcul/Conduction /
conduction_shape_factors.f90
program conduction_shape_factors
implicit none
! Pi constant
real(8), parameter :: pi = 3.14159265358979d0
! Inputs
integer :: config_idx
real(8) :: L_val, z_val, r_val, r1_val, r2_val, w_val, D_val, D1_val, D2_val, e_val, s_val, Di_val, Do_val, ri_val, ro_val
real(8) :: k_val
real(8) :: T1_val, T2_val
! Outputs
real(8) :: S_factor, R_th, Q_rate
integer :: err_flag
character(len=100) :: err_msg
character(len=50) :: config_name
! Inline acosh function helper
! acosh(x) = log(x + sqrt(x^2 - 1.0)) for x >= 1.0
real(8) :: x_temp
err_flag = 0
err_msg = ""
S_factor = 0.0d0
! Read configuration index
read *, config_idx
! Read global parameters: k, T1, T2
read *, k_val
read *, T1_val
read *, T2_val
select case (config_idx)
case (1)
! 1. Buried horizontal cylinder in semi-infinite medium
config_name = "Buried Horizontal Cylinder"
read *, L_val ! Length [m]
read *, z_val ! Depth [m]
read *, r_val ! Radius [m]
if (L_val <= 0.0d0 .or. z_val <= 0.0d0 .or. r_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (z_val < r_val) then
err_flag = 1
err_msg = "Burial depth (z) must be greater than or equal to cylinder radius (r)."
else
x_temp = z_val / r_val
S_factor = 2.0d0 * pi * L_val / log(x_temp + sqrt(x_temp**2 - 1.0d0))
end if
case (2)
! 2. Buried vertical cylinder in semi-infinite medium
config_name = "Buried Vertical Cylinder"
read *, L_val ! Length [m]
read *, D_val ! Diameter [m]
if (L_val <= 0.0d0 .or. D_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (L_val < 2.0d0 * D_val) then
err_flag = 1
err_msg = "Formula valid for L >= 2D (long cylinder/well)."
else
S_factor = 2.0d0 * pi * L_val / log(4.0d0 * L_val / D_val)
end if
case (3)
! 3. Two parallel buried cylinders in infinite medium
config_name = "Two Parallel Cylinders"
read *, L_val ! Length [m]
read *, w_val ! Spacing between centers [m]
read *, r1_val ! Radius 1 [m]
read *, r2_val ! Radius 2 [m]
if (L_val <= 0.0d0 .or. w_val <= 0.0d0 .or. r1_val <= 0.0d0 .or. r2_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (w_val <= r1_val + r2_val) then
err_flag = 1
err_msg = "Centerline spacing (w) must be greater than r1 + r2 to prevent overlap."
else
x_temp = (w_val**2 - r1_val**2 - r2_val**2) / (2.0d0 * r1_val * r2_val)
if (x_temp < 1.0d0) then
err_flag = 1
err_msg = "Invalid parameters for acosh calculation."
else
S_factor = 2.0d0 * pi * L_val / log(x_temp + sqrt(x_temp**2 - 1.0d0))
end if
end if
case (4)
! 4. Isothermal sphere in semi-infinite medium
config_name = "Buried Isothermal Sphere"
read *, r_val ! Radius [m]
read *, z_val ! Depth [m]
if (r_val <= 0.0d0 .or. z_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (z_val < r_val) then
err_flag = 1
err_msg = "Burial depth (z) must be greater than or equal to sphere radius (r)."
else
S_factor = 4.0d0 * pi * r_val / (1.0d0 - r_val / (2.0d0 * z_val))
end if
case (5)
! 5. Cylinder centered in square solid
config_name = "Cylinder in Square Solid"
read *, L_val ! Length [m]
read *, w_val ! Width of square side [m]
read *, D_val ! Cylinder diameter [m]
if (L_val <= 0.0d0 .or. w_val <= 0.0d0 .or. D_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (w_val < D_val) then
err_flag = 1
err_msg = "Square side width (w) must be greater than cylinder diameter (D)."
else
S_factor = 2.0d0 * pi * L_val / log(1.08d0 * w_val / D_val)
end if
case (6)
! 6. Eccentric cylinders
config_name = "Eccentric Concentric Cylinders"
read *, L_val ! Length [m]
read *, D1_val ! Inner cylinder outer diameter [m]
read *, D2_val ! Outer cylinder inner diameter [m]
read *, e_val ! Eccentricity (offset) [m]
if (L_val <= 0.0d0 .or. D1_val <= 0.0d0 .or. D2_val <= 0.0d0 .or. e_val < 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive, eccentricity non-negative."
else if (D2_val <= D1_val + 2.0d0 * e_val) then
err_flag = 1
err_msg = "Outer diameter D2 must be larger than D1 + 2*e to prevent wall collision."
else
x_temp = (D1_val**2 + D2_val**2 - 4.0d0 * e_val**2) / (2.0d0 * D1_val * D2_val)
if (x_temp < 1.0d0) then
err_flag = 1
err_msg = "Invalid parameters for acosh calculation."
else
S_factor = 2.0d0 * pi * L_val / log(x_temp + sqrt(x_temp**2 - 1.0d0))
end if
end if
case (7)
! 7. Edge of two adjoining walls
config_name = "Edge of Two Adjoining Walls"
read *, D_val ! Length of edge [m]
if (D_val <= 0.0d0) then
err_flag = 1
err_msg = "Edge length (D) must be positive."
else
S_factor = 0.54d0 * D_val
end if
case (8)
! 8. Corner of three walls
config_name = "Corner of Three Walls"
read *, L_val ! Wall thickness [m]
if (L_val <= 0.0d0) then
err_flag = 1
err_msg = "Wall thickness (L) must be positive."
else
S_factor = 0.15d0 * L_val
end if
case (9)
! 9. Disk on semi-infinite surface
config_name = "Isothermal Disk on Surface"
read *, r_val ! Disk radius [m]
if (r_val <= 0.0d0) then
err_flag = 1
err_msg = "Radius (r) must be positive."
else
S_factor = 4.0d0 * r_val
end if
case (10)
! 10. Hemisphere on surface
config_name = "Hemisphere on Surface"
read *, r_val ! Radius [m]
if (r_val <= 0.0d0) then
err_flag = 1
err_msg = "Radius (r) must be positive."
else
S_factor = 2.0d0 * pi * r_val
end if
case (11)
! 11. Row of cylinders in semi-infinite solid
config_name = "Row of Cylinders (Floor Heating)"
read *, L_val ! Length [m]
read *, s_val ! Center-to-center spacing [m]
read *, D_val ! Cylinder diameter [m]
read *, z_val ! Burial depth [m]
if (L_val <= 0.0d0 .or. s_val <= 0.0d0 .or. D_val <= 0.0d0 .or. z_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (s_val <= D_val) then
err_flag = 1
err_msg = "Spacing (s) must be greater than cylinder diameter (D)."
else if (z_val < D_val / 2.0d0) then
err_flag = 1
err_msg = "Burial depth (z) must be at least half the diameter."
else
x_temp = (s_val / (pi * D_val)) * sinh(2.0d0 * pi * z_val / s_val)
if (x_temp <= 0.0d0) then
err_flag = 1
err_msg = "Logarithm argument is non-positive."
else
S_factor = 2.0d0 * pi * L_val / log(x_temp)
end if
end if
case (12)
! 12. Concentric cylinders
config_name = "Concentric Cylinders (1D Radial)"
read *, L_val ! Length [m]
read *, Di_val ! Inner diameter [m]
read *, Do_val ! Outer diameter [m]
if (L_val <= 0.0d0 .or. Di_val <= 0.0d0 .or. Do_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (Do_val <= Di_val) then
err_flag = 1
err_msg = "Outer diameter (Do) must be greater than inner diameter (Di)."
else
S_factor = 2.0d0 * pi * L_val / log(Do_val / Di_val)
end if
case (13)
! 13. Concentric spheres
config_name = "Concentric Spheres (1D Radial)"
read *, ri_val ! Inner radius [m]
read *, ro_val ! Outer radius [m]
if (ri_val <= 0.0d0 .or. ro_val <= 0.0d0) then
err_flag = 1
err_msg = "Dimensions must be positive values."
else if (ro_val <= ri_val) then
err_flag = 1
err_msg = "Outer radius (ro) must be greater than inner radius (ri)."
else
S_factor = 4.0d0 * pi * ri_val * ro_val / (ro_val - ri_val)
end if
case default
err_flag = 1
err_msg = "Invalid configuration index chosen."
end select
! Calculate thermal resistance and heat flow rate if valid
if (err_flag == 0) then
if (k_val <= 0.0d0) then
err_flag = 1
err_msg = "Thermal conductivity (k) must be positive."
else
R_th = 1.0d0 / (S_factor * k_val)
Q_rate = S_factor * k_val * (T1_val - T2_val)
end if
end if
! Write output
print *, '=================================================='
print *, ' 2D CONDUCTION SHAPE FACTOR CALCULATION REPORT'
print *, '=================================================='
print *, ''
print *, 'CONFIGURATION DETAILS:'
print *, '--------------------------------------------------'
print *, ' Selected Profile: ', trim(config_name)
print *, ' Index Code: ', config_idx
print *, ''
if (err_flag /= 0) then
print *, '❌ ERROR DETECTED:'
print *, ' ', trim(err_msg)
print *, '=================================================='
stop
end if
print *, 'THERMO-PHYSICAL INPUTS:'
print *, '--------------------------------------------------'
print '(A,F10.4,A)', ' Conductivity of Solid (k): ', k_val, ' W/(m.K)'
print '(A,F10.2,A)', ' Boundary Temperature 1 (T1): ', T1_val, ' C'
print '(A,F10.2,A)', ' Boundary Temperature 2 (T2): ', T2_val, ' C'
print '(A,F10.2,A)', ' Temperature Difference (dT): ', T1_val - T2_val, ' C'
print *, ''
print *, 'CALCULATED THERMAL METRICS:'
print *, '--------------------------------------------------'
print '(A,F12.6,A)', ' Conduction Shape Factor (S): ', S_factor, ' m'
print '(A,F12.6,A)', ' Conduction Thermal Resistance (R):', R_th, ' C/W'
print '(A,F12.2,A)', ' Conduction Heat Flow Rate (Q): ', Q_rate, ' W'
print *, '=================================================='
end program conduction_shape_factors
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 conduction_shape_factors.f90 -o shape_factor_calc
2. Execution with input.txt redirection:
shape_factor_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 1.5 60.0 20.0 5.0 1.2 0.3
Parameter Description:
Configuration (1-12) Thermal conductivity k [W/m-K] Temperature T1 [°C] Temperature T2 [°C] Length L [m] Depth/Distance z [m] Radius r [m]