Proceedings of the 10th World Congress on Mechanical, Chemical, and Material Engineering (MCM'24)
Barcelona, Spain –August 22-24, 2024
Paper No. HTFF 269
DOI: 10.11159/htff24.269

Solving the Radiation Transfer Equation in Participating Media Using
Physics Informed Neural Networks
Pratibha Biswal1, Jetnis Avdijaj1, Alessandro Parente1,2, Axel Coussement1

pratibha.biswal@ulb.be
Aero-Thermo-Mechanics Department, Université Libre de Bruxelles, Belgium
Brussels Institute for Thermal-Fluid Systems and Clean Energy (BRITE), Université Libre de Bruxelles and Vrije
Universiteit Brussel, 1050 Brussels, Belgium
2
WEL Research Institute, Avenue Pasteur 6, 1300 Wavre, Belgium
1

Abstract - The radiative transfer equation (RTE) serves as a fundamental framework for modeling the propagation of electromagnetic
waves through a medium. Traditionally, solving the RTE has been challenging and computationally intensive. In this work, a physicsinformed neural network (PINN) model is used to solve the 1D radiative transfer equation. The PINN approach integrates physical laws
into the neural network training process, offering a novel way to address the computational complexities of RTE solution. The resultsfrom PINN model are validated against results from previous studies. Findings for various extinction coefficient are presented
demonstrating the efficacy and accuracy of the PINN approach. This work contributes to the theoretical understanding of the RTE andhighlights the potential of PINNs to enhance and streamline numerical methods in this domain.
Keywords: Thermal Radiation, Radiation Transfer Equation (RTE), Scattering, Participating media, Physics-informed neural network(PINN), artificial neural network (ANN)

1. Introduction

Radiative heat transfer in participating media, such as furnace gases, the atmosphere, and clouds, involves absorption,
emission, and scattering phenomena [1]. Thermal radiation interacting with a medium leads to energy absorption, reducing
transmitted energy, while scattering redirects radiation in multiple directions, causing out-scattering and in-scattering.
Scattering can be isotropic or anisotropic, influenced by factors like temperature, composition, and spectral properties. In
high-temperature environments like furnaces, combustion byproducts such as carbon dioxide and water vapor significantly
affect radiation absorption and scattering, with their properties being spectrally dependent and temperature-sensitive.
Additionally, soot further complicates these radiation interactions.
Traditionally, radiative transfer equation (RTE) solvers adopt either physics-based (stochastic) or deterministic
methods. Methods like Monte-Carlo [2] excel in parallel computing but face challenges with numerically handling
optically thick media and integrating with other physics, such as fluid mechanics. Consequently, we often turn to numericalmethods that account for the spectral properties of gas species and particulate matter, as well as their interactions within
combustion chambers. Methods like this such as finite volume [3] and finite element [4] methods incur significant
computational costs, rendering the overall process slow, cumbersome, and expensive. In addition, mesh-based methods(e.g discrete ordinate) are very sensitive to the computational dimension and suffer from the curse of dimensionality, given
the high dimensionality of RTE. In order to address these challenges, there is a growing need for RTE models that require
fewer computational resources and less time.
In recent computational trends, artificial neural networks (ANNs) are increasingly favoured over elaborate physicalmodels due to their minimal computational requirements [5]. One promising approach involves substituting exhaustive
RTE solutions with ANN prediction models in strongly scattering media. ANNs, however, often require large amounts ofdata and can struggle with generalization, higher computational costs and less accurate predictions. Recent breakthroughs
like physics-informed neural networks (PINNs) show potential for problems governed by partial differential equations [6].
Unlike pure data-driven ANN, PINNs embed differential equations into the training process [7]. These mesh-free models
selects random discrete points in the computational region (or take data from simulations/experiments), making them lesssensitive to dimensionality issues especially in RTE problems [8].

HTFF 269-1

In this work, the one-dimensional radiative transfer equation (RTE) is solved using a Physics-Informed Neural
Network (PINN) model. The physical principles of the RTE are integrated into the neural network framework, allowing for
efficient handling of the complexities of radiative transfer. The solutions obtained from the PINN model are validated
against previous studies, and results are presented in terms of radiation intensity and flux for various extinction
coefficients. This work not only enhances the theoretical understanding of the RTE but also suggests a fast and efficientmethod for solving thermal radiative transfer problems in scattering media, demonstrating the significant potential of
PINNs in this domain.
(a)

(b)

Figure 1 Radiation intensity I vs μ at (a) x=0 and (b) x=1 involving different Gaussian quadrature schemes.

2. Mathematical modelling[9].

2.1. Radiation transport equation (RTE), scattering phase function and boundary conditions
Generalized steady radiation transfer equation (RTE) for a system with participating medium can be written as follows

κs,λ ❑Ω ∙ ∇ I λ ( r , Ω )+ β λ I λ ( r , Ω ) =κ a , λ I b , λ [ T ( r ) ] +∫ I ¿ ¿, Ω ' ) Φ λ ( Ω∙ Ω' ) d Ω' (1)4 π 4π λHere, I λ ( r , Ω ) is specific intensity of radiation at position 𝑟 in direction Ω and wavelength 𝜆. β λ is the extinction
coefficient at wavelength 𝜆 which accounts for both absorption and scattering. κ a , λ and κ s , λ are absorption and scattering
coefficients, respectively at wavelength 𝜆. I b , λ [ T ( r ) ] is source term, typically representing the blackbody radiation
intensity at temperature 𝑇(𝑟). Φ λ ( Ω∙ Ω' ) is scattering phase function, which describes the angular distribution of scatteredradiation. In the left side of Eq. (1), first term represents rate of change of specific radiation intensity due to spatial
variation. Second term of LHS represents the extinction term due to absorption and scattering. In the right-hand side of Eq.
(1), the first term represents the radiation emission by the medium which can be described by a blackbody at temperatureT(r). The most important term for radiation in scattering medium is the last term at right hand side. This term represents
the scattering function that denotes radiation scattering from all directions. This term makes the equation [Eq. (1)] an
integro-differential equation. For monochromatic radiation in one-dimensional geometry involving non-absorbing medium,
the RTE can be expressed as follows:

HTFF 269-2

μκ 1
d I (x ,μ)+ βI ( x , μ )= s ∫ I ¿ ¿, μ ' ) Φ ( μ ∙ μ ' ) dμ'
dx
2 −1

(2)

In this equation, I ( x , μ ) is monochromatic specific intensity at position x in direction 𝜇. In this case, μ represents the
direction cosine (cos(θ)), representing the direction of radiation with respect to the x-axis. It ranges from −1 to 1. Theκs 1
spatial axis, x ranges from 0 to 1. Integral term in Eq. (2),∫ I ¿ ¿ , μ ' ) Φ ( μ ∙ μ ' ) dμ' integrates the contributions of
2 −1

radiative intensity from all directions 𝜇′ to the direction 𝜇, weighted by the phase function Φ(𝜇𝜇′). In one-dimension, thescattering phase function can be written as follows.

Φ(μ∙μ' ) =2π

1∫ Φ ( Ω∙ Ω ' ) d φ
2π 0

(3)

Here, Φ(Ω⋅Ω′) is the full phase function depending on the solid angles Ω and Ω′. φ is the azimuthal angle, integratingover 2𝜋 to account for all possible scattering directions in the plane perpendicular to the initial direction. The scattering
phase function can be expanded in terms of Legendre Polynomials, 𝑃𝑚(𝜇) to provide an orthogonal basis for the function
defined on the interval [−1,1]. The scattering phase function Φ (cosΘ ) describes the angular distribution of scatteredradiation and can be expressed in terms of expansion coefficient (Am) and Legendre polynomials (𝑃𝑚) asMΦ ( Ω∙ Ω' ) = Φ (cosΘ ) = 1+ ∑ Am Pm (cosΘ )

(4)

m=1

Figure 2 Schematic representation of the architecture of the physics informed neural network (PINN) used to solve
the RTE.
For this problem with one dimension,

HTFF 269-3

M

Φ ( μ ∙ μ ' ) = 1+ ∑ Am Pm ( μ ) Pm ( μ ' ¿ ¿

(5)

m=1

The final form of RTE in terms of Legendre polynomial approximation can be written as follows,

(

)

M
κs 1
d I (x ,μ)μ+ βI ( x , μ )= ∫ I ¿ ¿, μ ' ) 1+ ∑ A m Pm ( cos Θ) dμ '
dx
2 −1
m=1

(6)

In order to solve Eq. (6), the following boundary conditions are applied.
Incoming radiation intensity: I ( 1 , μ )=1 for μ ∈ (0,1]Outgoing radiation intensity: I ( 0 , μ ) =0 for μ ∈ [−1,0)2.2. Quadrature rules
Angular discretization is a key factor that dictates the numerical solution accuracy of radiative heat transfer problems.
The selection of a proper angular quadrature is often essential for an efficient solution. In methods like PINN, thecomputation is significantly affected by angular space quadrature. The solid angular space is two dimensional anddescribed by the zenith angle θ and azimuthal angle φ . For the one-dimensional case, the radiative intensity is only afunction of zenith angle θ due to axisymmetry. Five different methods are tested in this work. Gauss–Legendre, Gauss–
Labotto, Gauss-Chebyshev, Gauss-Hermite and Gauss-Kronrod are explored based on corresponding quadrature points andweights for integrating over the interval [−1, 1]. These works are compared with an analytical solution scheme as presented
in Figure 1. For a quadrature size of 10, Gauss-Legendre method provided better results compared to other methods.
Gauss-Legendre quadrature rule takes the following form
1

n

−1

i=1

∫ f ( μ ) d μ ≈ ∑ wi f ( μi ) (7)Here, n is the number of sample points used in the approximation; 𝑤𝑖 are the quadrature weights; μi are the roots of thei Legendre polynomial.
th

HTFF 269-4

Figure 3 Comparison of radiation intensity on μ-x plane as obtained in an earlier work of finite element method results (Pontaza andReddy [10]; Left panel) and present work (Right panel)). Figure from left panel is reprinted with permission from Elsevier.

2.3. Training points
As is customary in supervised learning, we need to generate or obtain data to train the network. Generally,
experimental or simulation data can be used for detailed PINN studies. In this work, low-discrepancy sequences are
considered for training set pertaining to the simplicity of the domain. Interior and boundary collocation points are
established with the 1D domain (0<x<1). These data will be the interior training points. These points are needed to be
strategically placed to capture the behavior of the system accurately. Each point represents a specific instance in space andintensity.

HTFF 269-5

Figure 4 Comparison of results of present work with an earlier work (Hu et al. [11]). (a) Distribution of radiation intensity withμ at x = 0 and (b) distribution of radiation intensity with μ at x = 1. (c) Radiation heat flux vs x. Figures in the left panel are reprinted
with permission from Elsevier.

Different methods are explored to understand the effects on the radiation intensity in the 1D domain. Each sampling
method offers a different approach to selecting points within the 1D domain. We focus on sequence sampling method, thatinvolves generating points in a deterministic, ordered sequence. These points are systematically distributed within the

HTFF 269-6

domain, ensuring a more uniform coverage. In this work, Sobol, Latin Hypercube Sampling (LHS), random and Haltonsequences are considered for tests and the results are not shown in this article for brevity. For the considered cases, allmethods other than random sampling showed promising results. For subsequent calculations, Sobol sequence is used.
2.4. Physics informed neural networks (PINNs)
A deep feedforward neural network is used transforming inputs into outputs through multiple layers of neurons. Eachlayer consists of affine-linear maps and scalar non-linear activation functions. The network has an input layer, an output
layer, and multiple hidden layers. In order to solve Eq. (2), the equation is approximated the DNN, which takes spatial ( x)and angular ( μ) variables as inputs. The outputs the approximated solution Ip , where p denotes the neural network
parameters including the weights and biases of neural networks. Figure 2 illustrated the architecture of the PINN. The NN
contains two input layers for x and μ. After detailed testing of network hyperparameters, 6 hidden layers with 24 neurons
each are considered. The tanh activation function is used for all neurons in the hidden layers after comparing with otheravailable activation functions.
The partial differential operator ∂/∂x is implemented by using automatic differentiation using PyTorch. The neuralnetwork parameters p are obtained by optimizing the loss function. Loss function includes residual of the RTE andboundary conditions, that makes the neural network dependent on physical governing equations.
For interior collocation points, the residual is defined as follows:Rκ 1(x , μ )+βI ( x , μ) − s ∫ I ¿ ¿¿,∫ ¿→ μ d I dx
2 −1

(

M

)

μ’) 1+ ∑ A m Pm ( cos Θ) d μ' −f (8)m=1

For the boundary conditions, the residuals are

Rbc 1 → I ( 1 , μ )−1 and Rbc 2 → I ( 0 , μ )−0(9)The objective is to minimize these residuals simultaneously to determine the weights and biases. The network, RTE
residual, boundary residuals, and loss functions are evaluated followed by an optimization algorithm (Adam) to obtain
optimal network parameters. The network is trained iteratively, adjusting the weights and biases to minimize the loss
function.

3. Benchmark and validation

Results from this work are validated with earlier physics based numerical simulation works. Figure 3 shows thecomparison of radiation intensity on spatial-angular ( μ− x) plane with an earlier work by Pontaza and Reddy [10]. They
have carried out 1D simulations for unit thickness using least square (LS) finite element method (FEM). The spatialangular distributions of the radiation intensity for the current work involving a PINN and the earlier work involving FEM
are in excellent agreement.
Further validation has been carried out with the detailed work done by Hu et al., [11]. Exit distributions of radiative
intensity (I−) at x = 0 and (I +) at x = 1 are plotted in Figure 4 (a) and (b). The results based on different methods angulardiscretization are in good agreement with slight discrepancy. Note that, LS denotes the least square method used by
Pontaza and Reddy [10]. CSM denotes collocation spectral method. The results are also compared with analytical solution
presented by Cengel and Ozisik [12].

4. Results and discussions

In the case of isotropic scattering, a linear space-dependent scattering coefficient, κ s=x, and unit extinction coefficient,
β=1, are considered. The corresponding results are presented earlier with the benchmark results for a unit extinction
coefficient. In thermal radiation, the extinction coefficient represents the ability of a medium to absorb and scatter radiation

HTFF 269-7

as it passes through. It characterizes how much radiation is absorbed or scattered per unit distance traveled through thematerial. In combustion, the extinction coefficient represents how radiation interacts with combustion products like gases,soot, and particles. In transparent media like air or water, β tends to be relatively low, reflecting minimal interaction ofradiation with the medium. However, in semi-transparent materials such as smoke or steam, β can increase substantially
due to higher levels of absorption and scattering by suspended particles or molecules.

Figure 5 Effect of extinction coefficient, β on (a-c) radiation intensity contour on x−μ plane [(a) β = 0.6, (b) β = 0.8 and (c) β=1] and (d) radiation flux.

The extinction coefficient contributes to this attenuation by determining the rate at which radiation is absorbed along
its path. The effect of the extinction coefficient on radiation intensity with μ at x=0 and x=1 is plotted but not shown tomaintain brevity. At lower values of β, radiation intensity is highest as the absorption into the medium is minimal. On theother hand, at larger β values, radiation intensity is minimum as the medium absorbs significant amount of radiation.
Common to all cases, radiation intensity decreases exponentially with distance as it passes through a medium due toabsorption and scattering processes. A clearer representation with the evolution of radiation intensity on the x-μ plane ispresented in Figure 5(a-c). A higher extinction coefficient indicates stronger absorption, resulting in a greater reduction in

HTFF 269-8

the intensity of radiation as it penetrates deeper into the medium. This is evident from the lower values of radiation
intensity for larger β.
A high extinction coefficient means that the material strongly absorbs thermal radiation, resulting in lower radiation
intensity. The material may heat up as it absorbs thermal radiation, leading to a higher temperature gradient across its
thickness or volume. If the material cannot efficiently dissipate the absorbed heat, it may accumulate thermal energy,
leading to further heating. On the other hand, a low extinction coefficient implies weaker absorption of thermal radiationby the material. Thermal radiation emitted by a hot object can more readily pass through the material. With less absorptionof thermal radiation, the temperature gradient across the material may be more uniform. This nature is clearly illustrated interms of radiation flux for various β in Figure 5(d). In environments where thermal radiation plays a significant role,materials with low extinction coefficients may contribute to improved thermal comfort. However, for applications likefurnace combustion, control and optimization of the extinction coefficient is necessary.

5. Conclusion

This study focuses on radiative heat transfer in participating media, applicable for various industrial andenvironmental applications, Solution of radiation transport equation (RTE) based on utilization of a Physics-Informed
Neural Network (PINN) model is explored. Through benchmarking and validation against established numerical methods,
PINN demonstrated the effectiveness and accuracy for the considered problems. Our future work is focused on 2D and 3D
geometries with higher dimensional systems. Detailed investigation into angular quadrature rules and sampling methods
are presented. For an isotropic medium, the significance of extinction coefficients on radiation intensity and flux is
explained. Findings of this work is pivotal. This study provides a kick-start approach to tackle high dimensionality curse of
RTE problems. Fast and reliable PINN based solution provides better theoretical understanding of radiative heat transfer inscattering media.

Acknowledgements

Authors acknowledge support and contributions of the corresponding funding agency, REPowerEU & Plan deRelance et de Resilience (PRR), which have been instrumental in the process of this research.

References[1] P. J. Coelho, “Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer
problems in participating media”, Journal of Quantitative Spectroscopy and Radiative Transfer, 2014, vol. 145, pp.
121-146.
[2] J. R. Howell, K. J. Daun, “The past and future of the Monte Carlo method in thermal radiation transfer”, Journal of
Heat Transfer, 2021, pp. 143, Art. No. 100801[3] S.A.N. Heugang, H.T.K. Tagne, F. B. Pelap. Discrete transfer AND finite volume methods FOR highly anisotropicallyscattering IN radiative heat analysis, Journal of Computational and Theoretical Transport, 2020, vol. 49, pp. 195-214[4] L.H. Liu, L. Zhang, H.P. Tan, “Finite element method for radiation heat transfer in multi-dimensional graded index
medium”, Journal of Quantitative Spectroscopy & Radiative Transfer, 2006, vol. pp. 436–445[5] Amit Kumar Yadav, S.S. Chandel, “Solar radiation prediction using Artificial Neural Network techniques: A review”,
Renewable and Sustainable Energy Reviews, 2014, vol. 33, pp. 772-781[6] S. Cuomo, V.S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, “Scientific machine learning throughphysics–informed neural networks: Where we are and what’s next.” Journal of Scientific Computing, 2022, vo. 92,
Art. No. 88.
[7] Q. A. Huhn, M. E. Tano, J. C. Ragusa, Physics-informed neural network with Fourier features for radiation transport inheterogeneous media, Nuclear Science and Engineering, 2023 vol. 197, pp. 2484-2497[8] S. Mishra, R. Molinaro, Physics informed neural networks for simulating radiative transfer Journal of Quantitative
Spectroscopy & Radiative Transfer, 2021, vol. 270, Art. No. 107705[9] M.F. Modest , Radiative Heat Transfer, Academic press, 2013.

HTFF 269-9

[10] J.P. Pontaza, J.N. Reddy, Least-squares finite element formulations for one-dimensional radiative transfer, Journal of
Quantitative Spectroscopy & Radiative Transfer, 2005, vol. 95, pp. 387-406[11] Z. Hu, H. Tian, B. Li, W. Zhang, Y. Yin, M. Ruan, D. Chen. Incident energy transfer equation and its solution bycollocation spectral method for one-dimensional radiative heat transfer, Journal of Quantitative Spectroscopy &Radiative Transfer, 2017, vol. 200, pp. 163–172[12] Y.A. Cengel, M.N. Ozisik, Radiation transfer in an anisotropically scattering plane-parallel medium with spacedependent albedo w(x): JQSRT, 1985, vol. 34, pp. 263–270.

HTFF 269-10