🔬
Solver Purpose & Physical Scope
Participating gases (such as $CO_2$ and $H_2O$) emit and absorb radiation in discrete spectral bands. Unlike solids, they do not act as black or gray bodies. The total gas mixture emissivity ($\varepsilon_g$) is modeled using Leckner's bivariate correlations (1972) as a function of temperature ($T_g$) and partial pressure-length parameter ($p_g L_e$):
📂 Discipline: Radiation
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 228 times
📄 Source File:
gas_radiation.f90
📁 calcul/Radiation /
gas_radiation.f90
program gas_radiation
implicit none
! Inputs
double precision :: Tg_C, Ts_C
double precision :: pc, pw, L, Pt
! Internal temperatures
double precision :: Tg_K, Ts_K
double precision :: tg, ts
! Leckner matrices
double precision :: c(0:3, 0:2), d(0:3, 0:2)
! Emissivity variables
double precision :: L_cm, pL_c, pL_w, pL_sum
double precision :: x_g, y_c, y_w
double precision :: eps_0c, eps_0w
double precision :: Pec, Pew, Cc, Cw
double precision :: eps_c, eps_w
double precision :: zeta, overlap_eps, eps_total
! Absorptivity variables
double precision :: pL_c_s, pL_w_s, pL_sum_s
double precision :: x_s, y_c_s, y_w_s
double precision :: eps_0c_s, eps_0w_s
double precision :: alpha_c, alpha_w, overlap_alpha, alpha_total, temp_term
integer :: i, j, iostat_val
! Initialize Leckner matrices
! CO2 Emissivity Matrix (c_ij)
c(0,0) = -3.5430d0; c(0,1) = 1.6370d0; c(0,2) = -0.5050d0
c(1,0) = 0.3800d0; c(1,1) = -0.2760d0; c(1,2) = 0.1060d0
c(2,0) = -0.1180d0; c(2,1) = 0.0540d0; c(2,2) = -0.0160d0
c(3,0) = 0.0170d0; c(3,1) = -0.0070d0; c(3,2) = 0.0020d0
! H2O Emissivity Matrix (d_ij)
d(0,0) = -2.2118d0; d(0,1) = -1.1987d0; d(0,2) = 0.035596d0
d(1,0) = 0.85667d0; d(1,1) = 0.93048d0; d(1,2) = -0.14391d0
d(2,0) = -0.10833d0; d(2,1) = -0.17156d0; d(2,2) = 0.045915d0
d(3,0) = -0.016333d0; d(3,1) = 0.01988d0; d(3,2) = -0.005830d0
! Read inputs
read(*,*,iostat=iostat_val) Tg_C
read(*,*,iostat=iostat_val) Ts_C
read(*,*,iostat=iostat_val) pc
read(*,*,iostat=iostat_val) pw
read(*,*,iostat=iostat_val) L
read(*,*,iostat=iostat_val) Pt
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read gas calculation parameters.'
stop
end if
! Basic validations
Tg_K = Tg_C + 273.15d0
Ts_K = Ts_C + 273.15d0
if (Tg_K <= 0.0d0 .or. Ts_K <= 0.0d0) then
write(*,*) 'ERROR: Temperatures must be above absolute zero.'
stop
end if
if (pc < 0.0d0 .or. pw < 0.0d0 .or. Pt <= 0.0d0) then
write(*,*) 'ERROR: Gas pressures must be positive.'
stop
end if
if ((pc + pw) > Pt) then
write(*,*) 'ERROR: Sum of partial pressures cannot exceed total pressure.'
stop
end if
if (L <= 0.0d0) then
write(*,*) 'ERROR: Mean beam length must be positive.'
stop
end if
! Units conversion: path length to cm
L_cm = L * 100.0d0
pL_c = pc * L_cm
pL_w = pw * L_cm
pL_sum = (pc + pw) * L_cm
tg = Tg_K / 1000.0d0
x_g = log(tg)
! =======================================================
! 1. EMISSIVITY CALCULATIONS (at gas temperature Tg)
! =======================================================
! Base Emissivity CO2
y_c = log(max(1.0001d0, pL_c))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + c(i, j) * (x_g**i) * (y_c**j)
end do
end do
eps_0c = exp(temp_term)
! Base Emissivity H2O
y_w = log(max(1.0001d0, pL_w))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + d(i, j) * (x_g**i) * (y_w**j)
end do
end do
eps_0w = exp(temp_term)
! Pressure corrections
Pec = Pt + 0.28d0 * pc
Cc = 1.0d0 + ((Pec - 1.0d0) / (Pec - 1.0d0 + 0.28d0)) * 0.15d0 * &
(1.0d0 - tanh(log10(max(1.0001d0, pL_c))))
Pew = Pt + 0.85d0 * pw
Cw = 1.0d0 + ((Pew - 1.0d0) / (Pew - 1.0d0 + 0.5d0)) * 0.35d0 * &
(1.0d0 - tanh(log10(max(1.0001d0, pL_w))))
eps_c = eps_0c * Cc
eps_w = eps_0w * Cw
! Overlap correction (zeta)
overlap_eps = 0.0d0
if ((pc + pw) > 1.0d-8) then
zeta = pw / (pc + pw)
if (pL_sum >= 1.0d0) then
overlap_eps = zeta * (10.7d0 / (10.7d0 + 101.0d0 * zeta) - 0.0089d0 * zeta**10.4d0) * &
(log10(pL_sum))**2.76d0
end if
end if
eps_total = eps_c + eps_w - overlap_eps
if (eps_total < 0.0d0) eps_total = 0.0d0
! =======================================================
! 2. ABSORPTIVITY CALCULATIONS (for surface at Ts)
! =======================================================
ts = Ts_K / 1000.0d0
x_s = log(ts)
! CO2 Absorptivity (evaluated at Ts and scaled pL)
pL_c_s = pc * L_cm * (Ts_K / Tg_K)
y_c_s = log(max(1.0001d0, pL_c_s))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + c(i, j) * (x_s**i) * (y_c_s**j)
end do
end do
eps_0c_s = exp(temp_term)
alpha_c = Cc * eps_0c_s * (Tg_K / Ts_K)**0.65d0
! H2O Absorptivity
pL_w_s = pw * L_cm * (Ts_K / Tg_K)
y_w_s = log(max(1.0001d0, pL_w_s))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + d(i, j) * (x_s**i) * (y_w_s**j)
end do
end do
eps_0w_s = exp(temp_term)
alpha_w = Cw * eps_0w_s * (Tg_K / Ts_K)**0.45d0
! Overlap correction for absorptivity
overlap_alpha = 0.0d0
pL_sum_s = (pc + pw) * L_cm * (Ts_K / Tg_K)
if ((pc + pw) > 1.0d-8) then
zeta = pw / (pc + pw)
if (pL_sum_s >= 1.0d0) then
overlap_alpha = zeta * (10.7d0 / (10.7d0 + 101.0d0 * zeta) - 0.0089d0 * zeta**10.4d0) * &
(log10(pL_sum_s))**2.76d0
end if
end if
alpha_total = alpha_c + alpha_w - overlap_alpha
if (alpha_total < 0.0d0) alpha_total = 0.0d0
! =======================================================
! OUTPUT DATA
! =======================================================
write(*,*) '============================================================'
write(*,*) ' LECKNER MIXTURE GAS EMISSION ENGINE'
write(*,*) '============================================================'
write(*,*)
write(*,'(A,F12.2,A)') ' Gas Temperature (Tg) = ', Tg_C, ' deg-C'
write(*,'(A,F12.2,A)') ' Wall Temperature (Ts) = ', Ts_C, ' deg-C'
write(*,'(A,F12.4,A)') ' Path Length (L) = ', L, ' m'
write(*,'(A,F12.4,A)') ' Total Pressure (Pt) = ', Pt, ' bar'
write(*,'(A,F12.4,A)') ' CO2 Partial Pressure = ', pc, ' bar'
write(*,'(A,F12.4,A)') ' H2O Partial Pressure = ', pw, ' bar'
write(*,*)
write(*,*) '--- GAS EMISSIVITY PROPERTIES ------------------------------'
write(*,'(A,F12.6)') ' CO2 Base Emissivity = ', eps_0c
write(*,'(A,F12.6)') ' CO2 Pressure Coeff Cc = ', Cc
write(*,'(A,F12.6)') ' CO2 Corrected Emiss = ', eps_c
write(*,'(A,F12.6)') ' H2O Base Emissivity = ', eps_0w
write(*,'(A,F12.6)') ' H2O Pressure Coeff Cw = ', Cw
write(*,'(A,F12.6)') ' H2O Corrected Emiss = ', eps_w
write(*,'(A,F12.6)') ' Mixture Overlap Corr = ', overlap_eps
write(*,'(A,F12.6)') ' Total Gas Emissivity = ', eps_total
write(*,*)
write(*,*) '--- GAS ABSORPTIVITY PROPERTIES ----------------------------'
write(*,'(A,F12.6)') ' CO2 Absorptivity = ', alpha_c
write(*,'(A,F12.6)') ' H2O Absorptivity = ', alpha_w
write(*,'(A,F12.6)') ' Mixture Overlap Corr = ', overlap_alpha
write(*,'(A,F12.6)') ' Total Gas Absorptivity = ', alpha_total
write(*,*)
end program gas_radiation
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 gas_radiation.f90 -o gas_radiation
2. Execution with input.txt redirection:
gas_radiation < input.txt
📄 Sample input.txt File Structure
Sample Data:
1000.0 300.0 0.10 0.15 2.0 1.0
Parameter Description:
Gas temperature Tg [°C] Wall temperature Ts [°C] CO2 partial pressure pc [bar] H2O partial pressure pw [bar] Mean beam length L [m] Total pressure Pt [bar]