🔬
Solver Purpose & Physical Scope
The view factor $F_{ij}$ is the fraction of radiation leaving surface $i$ that is directly intercepted by surface $j$. The general definition is given by the double area integral:
📂 Discipline: Radiation
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 449 times
📄 Source File:
view_factors.f90
📁 calcul/Radiation /
view_factors.f90
program view_factors
implicit none
! Inputs
integer :: geom_type
double precision :: dim1, dim2, dim3
! Constants
double precision, parameter :: PI = 3.141592653589793d0
! Calculated properties
double precision :: F12, F21
double precision :: X, Y, S, term1, term2, term3, term4, term5, ln_term
double precision :: r1, r2, A1, A2
integer :: iostat_val
! Read inputs
read(*,*,iostat=iostat_val) geom_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid geometry type.'
stop
end if
read(*,*,iostat=iostat_val) dim1
read(*,*,iostat=iostat_val) dim2
read(*,*,iostat=iostat_val) dim3
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read dimension parameters.'
stop
end if
! Basic validations
if (dim1 <= 0.0d0 .or. dim2 <= 0.0d0 .or. dim3 <= 0.0d0) then
write(*,*) 'ERROR: All geometric dimensions must be strictly positive.'
stop
end if
F12 = 0.0d0
F21 = 0.0d0
select case (geom_type)
case (1)
! Aligned parallel rectangles: dim1 = width (W), dim2 = height (H), dim3 = separation (L)
X = dim1 / dim3
Y = dim2 / dim3
term1 = log(((1.0d0 + X**2) * (1.0d0 + Y**2)) / (1.0d0 + X**2 + Y**2)) / 2.0d0
term2 = X * sqrt(1.0d0 + Y**2) * atan(X / sqrt(1.0d0 + Y**2))
term3 = Y * sqrt(1.0d0 + X**2) * atan(Y / sqrt(1.0d0 + X**2))
term4 = X * atan(X)
term5 = Y * atan(Y)
F12 = (2.0d0 / (PI * X * Y)) * (term1 + term2 + term3 - term4 - term5)
F21 = F12 ! Symmetric area
case (2)
! Perpendicular rectangles with common edge: dim1 = height (H), dim2 = length (L), dim3 = common edge (W)
X = dim1 / dim3
Y = dim2 / dim3
term1 = X * atan(1.0d0/X)
term2 = Y * atan(1.0d0/Y)
term3 = sqrt(X**2 + Y**2) * atan(1.0d0 / sqrt(X**2 + Y**2))
ln_term = ((1.0d0 + X**2)*(1.0d0 + Y**2)) / (1.0d0 + X**2 + Y**2) * &
( (X**2 * (1.0d0 + X**2 + Y**2)) / ((1.0d0 + X**2)*(X**2 + Y**2)) )**(X**2) * &
( (Y**2 * (1.0d0 + X**2 + Y**2)) / ((1.0d0 + Y**2)*(X**2 + Y**2)) )**(Y**2)
term4 = log(ln_term) / 4.0d0
F12 = (1.0d0 / (PI * X)) * (term1 + term2 - term3 + term4)
! Reciprocity: A1 = H*W, A2 = L*W => F21 = F12 * A1 / A2 = F12 * H / L
F21 = F12 * dim1 / dim2
case (3)
! Coaxial parallel disks: dim1 = radius 1 (R1), dim2 = radius 2 (R2), dim3 = separation (L)
r1 = dim1 / dim3
r2 = dim2 / dim3
S = 1.0d0 + (1.0d0 + r2**2) / (r1**2)
F12 = 0.5d0 * (S - sqrt(S**2 - 4.0d0 * (dim2/dim1)**2))
! Reciprocity: A1 = PI*R1^2, A2 = PI*R2^2 => F21 = F12 * R1^2 / R2^2
F21 = F12 * (dim1 / dim2)**2
case (4)
! Concentric cylinders (infinite): dim1 = inner radius (r1), dim2 = outer radius (r2), dim3 = length (ignored)
if (dim1 >= dim2) then
write(*,*) 'ERROR: Inner radius r1 must be less than outer radius r2.'
stop
end if
F12 = 1.0d0
F21 = dim1 / dim2
case (5)
! Concentric spheres: dim1 = inner radius (r1), dim2 = outer radius (r2), dim3 = ignored
if (dim1 >= dim2) then
write(*,*) 'ERROR: Inner radius r1 must be less than outer radius r2.'
stop
end if
F12 = 1.0d0
F21 = (dim1 / dim2)**2
case default
write(*,*) 'ERROR: Invalid geometry code.'
stop
end select
! Output results
write(*,'(A)') '============================================================'
write(*,'(A)') ' RADIATION VIEW FACTOR SOLVER'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,I2)') ' Geometry Code = ', geom_type
write(*,'(A,F12.4)') ' Dimension 1 (d1) = ', dim1
write(*,'(A,F12.4)') ' Dimension 2 (d2) = ', dim2
write(*,'(A,F12.4)') ' Dimension 3 (d3) = ', dim3
write(*,*)
write(*,'(A)') '--- CALCULATED VIEW FACTORS --------------------------------'
write(*,'(A,F12.6)') ' View Factor F1->2 = ', F12
write(*,'(A,F12.6)') ' View Factor F2->1 = ', F21
write(*,*)
write(*,'(A)') '--- RECIPROCITY VERIFICATION -------------------------------'
if (geom_type == 1) then
write(*,'(A)') ' Symmetric geometry: A1 = A2, F12 = F21.'
else if (geom_type == 2) then
write(*,'(A,F12.6)') ' Area ratio A1/A2 (H/L) = ', dim1 / dim2
write(*,'(A,F12.6)') ' F21 calculated = ', F21
else if (geom_type == 3) then
write(*,'(A,F12.6)') ' Area ratio A1/A2 = ', (dim1 / dim2)**2
write(*,'(A,F12.6)') ' F21 calculated = ', F21
else if (geom_type == 4) then
write(*,'(A,F12.6)') ' Area ratio A1/A2 (r1/r2)= ', dim1 / dim2
write(*,'(A,F12.6)') ' F21 calculated = ', F21
else if (geom_type == 5) then
write(*,'(A,F12.6)') ' Area ratio A1/A2 = ', (dim1 / dim2)**2
write(*,'(A,F12.6)') ' F21 calculated = ', F21
end if
end program view_factors
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 view_factors.f90 -o view_factors
2. Execution with input.txt redirection:
view_factors < 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]