PAPER • OPEN ACCESS

# A numerical method based on domain decomposition to solve coupled conduction-radiation physics using parallel computing within large porous media

To cite this article: Atin Kumar et al 2021 *J. Phys.: Conf. Ser.* 2116 012057

View the [article online](#) for updates and enhancements.

You may also like

- A domain decomposition preconditioning for an inverse volume scattering problem
Carlos Borges and George Biros

- Particle image velocimetry analysis with simultaneous uncertainty quantification using Bayesian neural networks
Mia C Morrell, Kyle Hickmann and Brandon M Wilson

- Comparing Macroscale and Microscale Simulations of Porous Battery Electrodes
Kenneth Higa, Shao-Ling Wu, Dilworth Y. Parkinson et al.

UNITED THROUGH SCIENCE & TECHNOLOGY

<figure><img src="image_1.png" /></figure>

The Electrochemical Society
Advancing solid state & electrochemical science & technology

248th
ECS Meeting
Chicago, IL
October 12-16, 2025
Hilton Chicago

<figure><img src="image_2.png" /></figure>

Science +
Technology +
YOU!

SUBMIT
ABSTRACTS by
March 28, 2025

SUBMIT NOW
# A numerical method based on domain decomposition to solve coupled conduction-radiation physics using parallel computing within large porous media

Atin Kumar¹,³, Jérôme Vicente², Jean-Vincent Daurelle², Yann Favennec¹ and Benoit Rousseau¹

¹LTeN, CNRS, UMR 6607, Université de Nantes, Nantes, France

²IUSTI, CNRS, UMR 7343, Aix-Marseille Université, Marseille, France

³Mersen, Amiens, France

E-mail: atin.kumar@mersen.com

**Abstract.** A domain decomposition approach is developed to solve coupled conductive-radiative heat transfer within highly porous materials. In this work, a Kelvin-cell foam with five cells in each direction which has $\approx 15.6 \times 10^6$ of voxels is considered. The coupled heat transfer is solved using the finite volume method where deterministic ray tracing is used to calculate radiative exchange. The temperature distribution is computed and cross-validated with the distribution obtained using a commercial software STAR-CCM+.

## 1. Introduction

Highly porous media possess many interesting properties, making them suitable for applications such as energy conservation, thermal protection systems, etc. Heat transfer through such materials may involve all three modes and its accurate estimation through numerical [1, 2] or experimental methods [3] is a challenging task. Designing and optimizing heat transfer through such complicated porous media can be leveraged by appropriate determination of effective heat transfer within the volume and optimizing relevant contributing textural parameters such as shape, size, porosity, etc. Till date, several methods are developed to approximate coupled heat transfer at local scale but they are limited to small and simple geometries as the computation over large and realistic geometries is computationally expensive.

Here, we present an approach based on domain decomposition to estimate effective heat transfer by coupling conduction and radiation in participating media (an opaque solid surrounded by vacuum). The work is based on the voxel-based finite volume method (FVM) which uses ray tracing for calculation of the radiative exchange factor. Using this approach, a Kelvin-cell foam which has $\approx 15.6 \times 10^6$ voxels is studied and results are shown and discussed.

## 2. Mathematical formulation

### 2.1. Governing equations

Coupled conductive-radiative heat transfer at steady state can be solved by energy balance as:

$$ \nabla \cdot (Q_{\text{cond}} + Q_{\text{rad}}) = 0 \quad \forall x \in \Omega \subset \mathbb{R}^3, \qquad (1) $$

<figure><img src="image_3.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
where $Q$ denotes the thermal energy, subscripts 'cond' and 'rad' denote conduction and radiation, respectively, $\Omega = \Omega_s \cup \Omega_v$ denotes the whole volume, subscripts 's' and 'v' denote solid and void, respectively, and $\boldsymbol{x} = (x, y, z)$ denotes the spatial location. The conductive heat energy is denoted as:

$$Q_{\text{cond}} = -\boldsymbol{\lambda}_s A \nabla T(\boldsymbol{x}) = 0 \quad \forall \boldsymbol{x} \in \Omega_s \subset \mathbb{R}^3, \qquad (2)$$

where $\boldsymbol{\lambda} = (\lambda_x, \lambda_y, \lambda_z)^T$ denotes the thermal conductivity vector, $A$ denotes the cross-sectional area perpendicular to the direction of heat flow, and $T$ is the temperature at spatial location $\boldsymbol{x}$. Further, the net radiative heat energy, adapted from eq. (8.14) in [4], can be represented as:

$$Q_{\text{rad}} = \sigma \epsilon n^2 A(\boldsymbol{x}) T^4(\boldsymbol{x}) - \int_{A'} \sigma \epsilon A(\boldsymbol{x}') T^4(\boldsymbol{x}') \mathcal{F}_{\boldsymbol{x}' \to \boldsymbol{x}} d\boldsymbol{A}' \quad \forall \boldsymbol{x} \in \Omega_v \subset \mathbb{R}^3, \qquad (3)$$

where $\boldsymbol{x}$ and $\boldsymbol{x}'$ denote emitting and reflecting surfaces, respectively, $T$ denotes the temperature, $A$ denotes the surface area, $\sigma$ denotes the Stefan-Boltzmann constant, $\epsilon$ denotes the total hemispherical emissivity of the surface, and $\mathcal{F}_{\boldsymbol{x}' \to \boldsymbol{x}}$ denotes the radiative exchange factor. To solve eq. (2) and eq. (3), following boundary conditions (BC) are imposed (for schematic representation, cf. fig 1):

* Dirichlet BC for eq. (2) is prescribed at hot and cold boundaries as:

$$T(x=0) = T_{\text{hot}}, \quad T(x=L_x) = T_{\text{cold}}; \qquad (4)$$

* Neumann BC for eq. (2) is prescribed on side faces to avoid any heat loss as:

$$\nabla T \cdot \boldsymbol{n}(y=0 \ \& \ y=L_y) = 0; \quad \nabla T \cdot \boldsymbol{n}(z=0 \ \& \ z=L_z) = 0; \qquad (5)$$

* The emitting BC is prescribed at fixed temperature boundaries. Either reflecting or absorbing (surface at $T=0$ K) BC is prescribed at all four side walls;

* Also, the reflecting BC is applied at the interface with reflection where $\rho = 1 - \epsilon$.

## 2.2. Numerical methods

The domain decomposition approach [5] is based on division of the entire volume into subvolumes of almost equal size (cf. fig 1). It requires calculation of the equivalent thermal conductivity vector ($\boldsymbol{\lambda}_{i,j}$) of each subvolume and radiative exchange factor between subvolumes ($\mathcal{F}_{i,j \to m,n}$) $\forall i, j, m, n = 1, 2, ..., N$ where $N$ denotes the number of subvolumes.

A voxel-based FVM is used to solve the steady-state heat conduction equation by performing the energy balance on each subvolume (resp. control volume (CV) at discrete level to calculate $\boldsymbol{\lambda}_{i,j}$) [1, 6]. The temperature distribution throughout the domain is calculated using the second order accurate finite difference stencil formula as:

$$T_{i,j}^{n+1} = (w_{i+1,j} T_{i+1,j}^n + w_{i-1,j} T_{i-1,j}^n + w_{i,j+1} T_{i,j+1}^n + w_{i,j-1} T_{i,j-1}^n)/w_{i,j}^t, \quad i, j = 1, 2, ..., N, \qquad (6)$$

where $T$ is the temperature of the subvolume (resp. CV) corresponding to attached index, $w$ is the weight of exchange with the neighbour subvolume (resp. CV) and $w_{i,j}^t = w_{i+1,j} + w_{i-1,j} + w_{i,j+1} + w_{i,j-1}$ is the total weight. The $w_{i+1,j}$ (i.e., in x-direction) is calculated as $2/\left[\frac{dx/2}{\lambda_x(i,j)} + \frac{dx/2}{\lambda_x(i+1,j)}\right]$, where $\lambda_x$ is the conductivity (resp. CV) and $dx$ is the length of the subvolume (resp. CV) in x-direction. The development of the method is relatively simple, the code is cross-validated with several conventional test cases.

Furthermore, the ray tracing approach based on a fast voxel traversal algorithm [7] is used to compute the radiative exchange factor. Rays are emitted from the center of each triangle,
<figure><img src="image_4.png" /><figcaption>Figure 1: Schematic representation of a porous medium, placed between solid plates and decomposed into subvolumes, with prescribed boundary conditions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2: Kelvin-cell foam generated in genMat and divided into 5 × 5 × 5 subvolumes.</figcaption></figure>

towards the void phase, in prescribed number of directions within the hemisphere, and they are traced until they are fully absorbed or leave the system. Once $\mathcal{F}_{(i,j)\to(m,n)}$ is computed, the temperature distribution throughout the domain can be computed using the discretized form of eq. (3) as:

$$ \sum_{\substack{i,j=1 \\ i,j \neq m,n}}^{N} h_{r,(i,j)\to(m,n)} \mathcal{F}_{(i,j)\to(m,n)}(T_{i,j} - T_{m,n}) = 0 \quad \forall m, n = 1, 2, \dots, N, \quad (7) $$

where $h_{r,(i,j)\to(m,n)}$ is a constant and is equal to $\epsilon\sigma A_{i,j} (T_{i,j}^2 + T_{m,n}^2) (T_{i,j} + T_{m,n})$ and $A$ denotes the total surface area of all triangles. Generally, computation of $\mathcal{F}_{(i,j)\to(m,n)}$ is long and can be memory intensive if one records the exchange between each triangle [1]. It would then be nearly impossible to consider a large mesh with millions of surfaces. Using domain decomposition approach, we calculate the radiative exchange among all the triangles within the volume but store the radiative exchange factor between subvolumes. Finally, the coupling is performed over a simplified geometry having subvolumes equal to number of subdivisions.

### 3. Numerical experiments

The 3D geometry (or 3D image) used for the computation is generated using an in-house software genMat developed in C++ and Qt [8]. The geometry was placed between two walls whose faces were maintained at fixed temperatures ($T_{\text{hot}} = 1800$ K and $T_{\text{cold}} = 1200$ K) to generate strong temperature gradient. The 3D geometry is composed of 250 voxels in each direction giving total number of voxels $\approx 15.6 \times 10^6$. The foam has a porosity of 90% (cf. fig 2.). We assumed that the solid phase is opaque and optically smooth with a local thermal conductivity ($\lambda_s$) equal to $0.3 \text{ W m}^{-1} \text{ K}^{-1}$ and a fixed emmisivity ($\epsilon$) equal to 0.9.

The plane-averaged temperature profile of a Kelvin-cell foam for both absorbing and reflecting BC is shown in fig 3. The plane-averaged temperature $\bar{T} = \frac{1}{A} \int_{A_s} T(x = x_c)dA$ is the average of temperature over a slice positioned at $x = x_c$, where $A$ is the cross-sectional area perpendicular to the direction of heat flow. The number of points at which the temperature is computed is equal to the number of divisions in the direction of heat flow. The results obtained using our approach were compared with the results produced using the STAR-CCM+. STAR-CCM+ is also based on FVM and uses ray tracing to calculate the view factor.

From fig 3, it is evident that the results are in good agreement with those of STAR-CCM+, hence cross-validating our approach. For absorbing boundaries, the temperature in the middle
<figure><img src="image_6.png" /><figcaption>Figure 3: Plane-averaged temperature ($\bar{T}$) profile from conduction and conduction-radiation computation for Absorbing (left) and Reflecting (right) boundary conditions</figcaption></figure>

goes lower than those fixed at hot and cold boundaries. This occurs due to the loss of heat due to radiation (as the temperature at the side walls is 0 K) and gives a significant deviation from pure conductive profile. For reflecting boundaries, the temperature does not go beyond those fixed at the inlet and outlet. Further, we noticed asymmetric behavior of temperature profile as in [1]. The average temperature increases near the cold boundary as high intensity radiation coming from the hot region will be absorbed here and vice versa occurs near the hot boundary.

Furthermore, the number of rays (or directions) used for exchange factor computation may have a significant effect on correct estimation of radiative heat transfer. For absorbing boundary, the effect is negligible because most of the rays striking the side faces will be absorbed permanently. However, this is not the case for reflecting boundaries as the reflected rays can travel further before being absorbed. Hence, it is necessary for reflecting boundaries to consider maximum possible directions to account for proper exchange of radiation, as shown in fig 3.

## 4. Conclusion

In this work, a method based on domain decomposition was established to solve the coupled conductive-radiative heat transfer within Kelvin-cell foam. The coupling requires $\lambda_{i,j}$ and $F_{(i,j)→(m,n)}$ whose computation is independent of the temperature and simplifies the coupling efforts. The results obtained for reflecting and absorbing boundaries are in good agreement with those computed using STAR-CCM+. It was found that the effect of number of rays on results may be significant depending on the BC and geometry, and should be used carefully.

This approach is a new milestone in dealing with large and complex geometries composed of millions of voxels due to low computational effort. In future, it will be extended to solve bigger and more complex geometries such as fibrous materials composed of randomly distributed fibers.

## References

[1] Perraudin D Y S and Haussener S 2017 *International Journal of Heat and Mass Transfer*. **112** 387-400

[2] Badri M A, Favennec Y, Jolivet P and Rousseau B 2020 *Finite Elements in Analysis and Design*. **178**

[3] Coquard R, Rochais D and Baillis D 2009 *International Journal of Heat and Mass Transfer*. **52** 4907-4918

[4] Modest M F 2013 *Radiative heat transfer* Academic press

[5] Kumar A 2021 *Textural Control of Coupled Radiative and Conductive Heat Transfer in Fibrous Media at Very High Temperatures* University of Nantes

[6] Mendes M A A Ray S and Trimis D 2013 *International Journal of Heat and Mass Transfer*. **66** 412-422

[7] Amanatides J and Woo A 1987 *Eurographics*. **87** 3-10

[8] Guévelou S, Rousseau B, Domingues G and Vicente J 2017 *J. Quant. Spectrosc. Radiat. Transfer*. **189** 329-338