💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Gray Enclosures & Shields
Core Numerical Engine in Fortran 90 • 49 total downloads
gray_radiation.f90
! =========================================================================
! Source File: 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
Solver Description
For grey, diffuse, and opaque surfaces, radiation exchange is solved using a network of thermal resistances. The surface resistance represents the emissivity barrier:
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 gray_radiation.f90 -o gray_radiation
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
gray_radiation < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Mode (1=Two-Surface Enclosure, 2=Radiation Shields)
1
! Enclosure Geometry (1=Parallel Plates, 2=Concentric Cylinders, 3=Concentric Spheres)
1
! Temperature T1 [°C]
600.0
! Temperature T2 [°C]
100.0
! Area A1 [m2]
1.0
! Area A2 [m2]
1.0
! Emissivity eps1
0.8
! Emissivity eps2
0.2
! View Factor F12
1.0
1
! Enclosure Geometry (1=Parallel Plates, 2=Concentric Cylinders, 3=Concentric Spheres)
1
! Temperature T1 [°C]
600.0
! Temperature T2 [°C]
100.0
! Area A1 [m2]
1.0
! Area A2 [m2]
1.0
! Emissivity eps1
0.8
! Emissivity eps2
0.2
! View Factor F12
1.0