RasmahRasmah

Radiation Transport: Discrete Ordinates, Scattering, and View Factors

using Rasmah
using LinearAlgebra

Overview

Radiation carries heat by photons, and in a participating medium the radiative intensity obeys an integro-differential transport equation. Rasmah solves it three complementary ways: the discrete-ordinates (S_N) method for directional resolution, anisotropic scattering phase functions, and the view factors of surface-to-surface exchange. Each is verified against its classical closed form below.

Discrete ordinates: the S_N slab

A gray, purely-absorbing slab of optical thickness τ has the exact transmitted flux 2 E_3(τ), where E_3 is the third exponential integral. The S_N method resolves directions with a Gauss–Legendre quadrature and converges to that value as the number of ordinates grows:

τ = 1.0
exact = 2 * e3(τ)
println("exact 2·E_3(τ) = ", round(exact, digits=6))
for n in (2, 4, 8, 16)
    println("S_$n flux = ", round(discrete_ordinates_slab_flux(τ, n), digits=6))
end
exact 2·E_3(τ) = 0.219384
S_2 flux = 0.20429
S_4 flux = 0.210985
S_8 flux = 0.219359
S_16 flux = 0.219406

Anisotropic scattering

Scattering redistributes radiation without absorbing it. The Henyey–Greenstein phase function Φ(μ) = (1−g²)/(4π(1+g²−2gμ)^{3/2}) is normalized and forward-peaked for g > 0:

# normalization 2π ∫ Φ dμ = 1
println("normalization g=0.5 = ", round(phase_normalization(henyey_greenstein, 0.5), digits=6))
# forward scattering (g > 0): Φ(1) > Φ(-1)
println("Φ(1) = ", round(henyey_greenstein(1.0, 0.5), digits=4), "   Φ(-1) = ", round(henyey_greenstein(-1.0, 0.5), digits=4))
normalization g=0.5 = 1.0
Φ(1) = 0.4775   Φ(-1) = 0.0177

Surface-to-surface view factors

Between surfaces, the view factor F_{1→2} is the fraction of radiation leaving 1 that strikes 2. For two identical parallel rectangles it has the Howell-catalogue closed form, and it satisfies reciprocity:

println("F(1,1) = ", round(parallel_rectangle_view_factor(1.0, 1.0), digits=6), "  (parallel squares)")
println("diffuse element → disk: F = ", differential_to_disk_view_factor(1.0, 1.0), "  (= R²/(R²+d²))")
println("reciprocity: F21 = ", view_factor_reciprocity(1.0, 2.0, 0.4), "  (= A1/A2 · F12)")
F(1,1) = 0.199825  (parallel squares)
diffuse element → disk: F = 0.5  (= R²/(R²+d²))
reciprocity: F21 = 0.2  (= A1/A2 · F12)

Next steps