🔬
Solver Purpose & Physical Scope
For grey, diffuse, and opaque surfaces, radiation exchange is solved using a network of thermal resistances. The surface resistance represents the emissivity barrier:
📂 Discipline: Radiation
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 364 times
📄 Source File:
gray_radiation.f90
📁 calcul/Radiation /
gray_radiation.f90
program gray_radiation
implicit none
! Inputs
integer :: mode
integer :: geom_type
double precision :: T1_C, T2_C
double precision :: A1, A2
double precision :: eps1, eps2
double precision :: F12
! Shields parameters
integer :: n_shields
double precision, dimension(5) :: eps_s
double precision, dimension(5) :: Ts_K, Ts_C
! Constants
double precision, parameter :: SIGMA = 5.670374d-8
! Calculated properties
double precision :: T1_K, T2_K
double precision :: Eb1, Eb2
double precision :: R_s1, R_sp, R_s2, R_total
double precision :: J1, J2, Q12
double precision :: Q_no_shields, Q_with_shields, R_shields, reduction
double precision :: temp_term
integer :: i, iostat_val
! Read inputs
read(*,*,iostat=iostat_val) mode
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid calculation mode.'
stop
end if
if (mode == 1) then
! Mode 1: Two-surface enclosure
read(*,*,iostat=iostat_val) geom_type
read(*,*,iostat=iostat_val) T1_C
read(*,*,iostat=iostat_val) T2_C
read(*,*,iostat=iostat_val) A1
read(*,*,iostat=iostat_val) A2
read(*,*,iostat=iostat_val) eps1
read(*,*,iostat=iostat_val) eps2
read(*,*,iostat=iostat_val) F12
else if (mode == 2) then
! Mode 2: Radiation shields
read(*,*,iostat=iostat_val) T1_C
read(*,*,iostat=iostat_val) T2_C
read(*,*,iostat=iostat_val) A1
read(*,*,iostat=iostat_val) eps1
read(*,*,iostat=iostat_val) eps2
read(*,*,iostat=iostat_val) n_shields
do i = 1, 5
read(*,*,iostat=iostat_val) eps_s(i)
end do
else
write(*,*) 'ERROR: Invalid calculation mode.'
stop
end if
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read inputs.'
stop
end if
! Basic validation
T1_K = T1_C + 273.15d0
T2_K = T2_C + 273.15d0
if (T1_K <= 0.0d0 .or. T2_K <= 0.0d0) then
write(*,*) 'ERROR: Temperatures must be above absolute zero.'
stop
end if
if (A1 <= 0.0d0) then
write(*,*) 'ERROR: Area 1 must be positive.'
stop
end if
if (eps1 <= 0.0d0 .or. eps1 > 1.0d0 .or. eps2 <= 0.0d0 .or. eps2 > 1.0d0) then
write(*,*) 'ERROR: Surface emissivities must be in range ]0.0, 1.0].'
stop
end if
Eb1 = SIGMA * T1_K**4
Eb2 = SIGMA * T2_K**4
! Select Mode for logic execution
if (mode == 1) then
! Mode 1: Two-surface gray enclosures
select case (geom_type)
case (1)
! Parallel plates
if (F12 <= 0.0d0 .or. F12 > 1.0d0) then
write(*,*) 'ERROR: View factor must be between 0.0 and 1.0.'
stop
end if
case (2, 3)
! Concentric cylinders/spheres => F12 = 1.0
F12 = 1.0d0
if (A1 >= A2) then
write(*,*) 'ERROR: Inner area A1 must be less than outer area A2.'
stop
end if
case default
write(*,*) 'ERROR: Invalid geometry type selected.'
stop
end select
! Resistor network resistances
R_s1 = (1.0d0 - eps1) / (eps1 * A1)
R_sp = 1.0d0 / (A1 * F12)
R_s2 = (1.0d0 - eps2) / (eps2 * A2)
R_total = R_s1 + R_sp + R_s2
! Heat exchange
Q12 = (Eb1 - Eb2) / R_total
! Radiosities
J1 = Eb1 - Q12 * R_s1
J2 = Eb2 + Q12 * R_s2
! Output Results
write(*,'(A)') '============================================================'
write(*,'(A)') ' TWO-SURFACE GRAY ENCLOSURE THERMAL SOLVER'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,F12.2,A)') ' T1 (Blackbody Node Eb1) = ', T1_C, ' deg-C'
write(*,'(A,F12.2,A)') ' T2 (Blackbody Node Eb2) = ', T2_C, ' deg-C'
write(*,'(A,ES12.4,A)') ' Eb1 = ', Eb1, ' W/m2'
write(*,'(A,ES12.4,A)') ' Eb2 = ', Eb2, ' W/m2'
write(*,*)
write(*,'(A)') '--- RESISTOR NETWORK VALUES --------------------------------'
write(*,'(A,ES12.4,A)') ' Surface Resistance R1 = ', R_s1, ' 1/m2'
write(*,'(A,ES12.4,A)') ' Space Resistance R12 = ', R_sp, ' 1/m2'
write(*,'(A,ES12.4,A)') ' Surface Resistance R2 = ', R_s2, ' 1/m2'
write(*,'(A,ES12.4,A)') ' Total Resistance R_tot = ', R_total, ' 1/m2'
write(*,*)
write(*,'(A)') '--- RADIOSITY NODES AND FLUX -------------------------------'
write(*,'(A,ES12.4,A)') ' Radiosity Node J1 = ', J1, ' W/m2'
write(*,'(A,ES12.4,A)') ' Radiosity Node J2 = ', J2, ' W/m2'
write(*,'(A,ES12.4,A)') ' Net Heat Transfer Rate = ', Q12, ' W'
write(*,*)
else if (mode == 2) then
! Mode 2: Multi-shield Radiation Screens
if (n_shields < 1 .or. n_shields > 5) then
write(*,*) 'ERROR: Number of shields must be between 1 and 5.'
stop
end if
! Calculate resistance without shields
R_total = (1.0d0/eps1 + 1.0d0/eps2 - 1.0d0) / A1
Q_no_shields = (Eb1 - Eb2) / R_total
! Calculate resistance of shields
R_shields = 0.0d0
do i = 1, n_shields
if (eps_s(i) <= 0.0d0 .or. eps_s(i) > 1.0d0) then
write(*,*) 'ERROR: Shield emissivity must be in range ]0, 1].'
stop
end if
R_shields = R_shields + (2.0d0 / eps_s(i) - 1.0d0) / A1
end do
Q_with_shields = (Eb1 - Eb2) / (R_total + R_shields)
reduction = 100.0d0 * (1.0d0 - Q_with_shields / Q_no_shields)
! Calculate shield temperatures sequentially
! Q = A1 * sigma * (Ts_{i-1}^4 - Ts_i^4) / (1/eps_s(i-1) + 1/eps_s(i) - 1)
do i = 1, n_shields
if (i == 1) then
temp_term = T1_K**4 - (Q_with_shields / (A1 * SIGMA)) * (1.0d0/eps1 + 1.0d0/eps_s(1) - 1.0d0)
else
temp_term = Ts_K(i-1)**4 - (Q_with_shields / (A1 * SIGMA)) * (1.0d0/eps_s(i-1) + 1.0d0/eps_s(i) - 1.0d0)
end if
if (temp_term < 0.0d0) then
Ts_K(i) = 0.0d0
else
Ts_K(i) = temp_term**0.25d0
end if
Ts_C(i) = Ts_K(i) - 273.15d0
end do
! Output Results
write(*,'(A)') '============================================================'
write(*,'(A)') ' RADIATION SHIELDS SYSTEM SOLVER'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,F12.2,A)') ' T1 (Hot Surface) = ', T1_C, ' deg-C'
write(*,'(A,F12.2,A)') ' T2 (Cold Surface) = ', T2_C, ' deg-C'
write(*,'(A,I2)') ' Number of Shields = ', n_shields
write(*,'(A,ES12.4,A)') ' Q without Shields = ', Q_no_shields, ' W'
write(*,'(A,ES12.4,A)') ' Q with Shields = ', Q_with_shields, ' W'
write(*,'(A,F12.2,A)') ' Heat Transfer Reduction = ', reduction, ' %'
write(*,*)
write(*,'(A)') '--- SHIELD TEMPERATURE DISTRIBUTION ------------------------'
write(*,'(A)') ' Shield # Emissivity Temp [K] Temp [deg-C]'
write(*,'(A)') ' ----------------------------------------------------------'
do i = 1, n_shields
write(*,'(I5,10X,F6.3,9X,F8.2,7X,F8.2)') i, eps_s(i), Ts_K(i), Ts_C(i)
end do
write(*,*)
end if
end program gray_radiation
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 gray_radiation.f90 -o gray_radiation
2. Execution with input.txt redirection:
gray_radiation < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 1 600.0 100.0 1.0 1.0 0.8 0.2 1.0
Parameter Description:
Mode (1=Two-Surface Enclosure, 2=Radiation Shields) Enclosure Geometry (1=Parallel Plates, 2=Concentric Cylinders, 3=Concentric Spheres) Temperature T1 [°C] Temperature T2 [°C] Area A1 [m2] Area A2 [m2] Emissivity eps1 Emissivity eps2 View Factor F12