*2.2. Solving Temperature Distribution with FDM*

The temperature distribution inside the computational domain is obtained by discretizing the energy equation with FDM by using the standard central difference scheme in space and the forward difference scheme in time [37]. After discretization, the equation for finding the temperature value at a grid point in new time level is (*q* in the below expression represents the heat source term due to the constant temperature BC on the hexagonal edges)

$$T\_{i,j}^{\text{new}} = T\_{i,j}^{\text{old}} + (\text{RHS}^{\text{old}} + q\_{i,j}^{\text{old}}) \Delta t\_{\prime} \tag{5}$$

with

$$\begin{array}{lcl} \text{RHS}^{\text{old}} = & -4\alpha T\_{i,j}^{\text{old}} + \left(\alpha - \frac{1}{2}\mu\_{i+1,j}^{\text{old}}\right)T\_{i+1,j}^{\text{old}} + \left(\alpha + \frac{1}{2}\mu\_{i-1,j}^{\text{old}}\right)T\_{i-1,j}^{\text{old}} \\ & + \left(\alpha - \frac{1}{2}\nu\_{i,j+1}^{\text{old}}\right)T\_{i,j+1}^{\text{old}} + \left(\alpha + \frac{1}{2}\nu\_{i,j-1}^{\text{old}}\right)T\_{i,j-1}^{\text{old}} \end{array} \tag{6}$$

In the above equation, α is the thermal diffusivity, and the subscripts *i* & *j* represent lattice grid indices in *x* − & *y*− directions, respectively.

### *2.3. Evaluation of ffl*(*<sup>x</sup>*, *t*) *and q*(*<sup>x</sup>*, *t*) *with SPM*

In SPM, a smoothed profile function (also termed as concentration function or indicator function), φ*k*(**<sup>x</sup>**, *t*), is used to identify the solid regions [13] (here, *k* is the index value for the hexagonal cylinders; *k* = 1 for the inner hexagon and *k* = 2 outer one). The equation for φ*k*(**<sup>x</sup>**,*<sup>t</sup>*) is defined in such a way that φ*k* = 0 in the fluid region, φ*k* = 1 in the solid region, and φ*k* smoothly varies from 0 to 1 at the fluid-solid interface. Here, the following equation is used to evaluate φ*k*(**<sup>x</sup>**, *t*) of each hexagon

$$
\phi\_k(\mathbf{x}, t) = f(\mathbf{d}\_k(\mathbf{x}, t)),
\tag{7a}
$$

$$f(\mathbf{d}\_k(\mathbf{x}, t)) = \begin{cases} 0 & r < -\xi\_k/2 \\ \frac{1}{2} \left( \sin \left( \pi \frac{\mathbf{d}\_k(\mathbf{x}, t)}{\xi\_k} \right) + 1 \right) & |r| < \xi\_k/2 \\ 1 & r > \xi\_k/2 \end{cases},\tag{7b}$$

where **d***k*(**<sup>x</sup>**,*<sup>t</sup>*) and ξ*k* are the signed normal distance function to the solid surface (here, edges of two hexagons) and the interface thickness of each hexagon, respectively. Unless otherwise mentioned, throughout this work, the values for the interface thickness for the two hexagons are chosen as, ξ1 = ξ2 ≡ 0.5. The following equation is used for finding **d***k*(**<sup>x</sup>**, *t*) of each hexagon

$$\mathbf{d}\_{k}(\mathbf{x},t) = \max\{\frac{L\_{\mathbf{x},k}}{2}, \max\left[\left(\left|\mathbf{x} - X\_{\varepsilon,k}\right|\sin(30^{\circ}) + \left|y - Y\_{\varepsilon,k}\right|\cos(30^{\circ})\right), \left(\left|y - Y\_{\varepsilon,k}\right| - \frac{L\_{yk}}{2}\right)\right] \tag{8}$$

*Lx*,*<sup>k</sup>* and *Ly*,*<sup>k</sup>* in the above equation are the distance between two corners and two opposite sides of each hexagon, respectively (*Ly*,<sup>1</sup> = *Lin* and *Ly*,<sup>2</sup> = *Lout* are the sizes of the two hexagons), and *Xc*,*<sup>k</sup>* and *Yc*,*<sup>k</sup>* are the center positions of hexagons, in *x*− and *y*− directions, respectively. The smoothly distributed concentration field of two hexagons, φ(**<sup>x</sup>**, *t*), is obtained by adding the φ*k*(**<sup>x</sup>**, *t*) values of two hexagons

$$\phi(\mathbf{x},t) = \sum\_{k=1}^{2} \phi\_k(\mathbf{x},t). \tag{9}$$

The body force term for enforcing no-slip BC on hexagonal edges, **<sup>f</sup>**fl(**<sup>x</sup>**,*<sup>t</sup>*) in Equation (1), is evaluated by using

$$\mathbf{f}^{\mathrm{fl}}(\mathbf{x},t) = \left[\mathbf{u}\_{\mathrm{P}}(\mathbf{x},t) - \mathbf{u}(\mathbf{x},t)\right] \phi(\mathbf{x},t) / \Delta t,\tag{10}$$

**<sup>u</sup>**p(**<sup>x</sup>**,*<sup>t</sup>*) in the above equation is the velocity field for the solid regions, which is zero as the two hexagonal cylinders are stationary. Similarly, the heat source term for treating the constant temperature BC, *q*(**<sup>x</sup>**, *t*) of Equation (5), can be obtained by

$$q(\mathbf{x},t) = \left[T\_{\mathbb{P}}(\mathbf{x},t) - T(\mathbf{x},t)\right] \phi(\mathbf{x},t) / \Delta t,\tag{11}$$

where *<sup>T</sup>*p(**<sup>x</sup>**, *t*) is the hexagonal cylinders temperature field, which is evaluated by using

$$\phi(\mathbf{x},t)T\_{\mathbb{P}}(\mathbf{x},t) = \sum\_{k=1}^{2} \phi\_{\mathbb{k}}(\mathbf{x},t)T\_{\mathbb{k}}(t),\tag{12}$$

where *Tk*(*t*) is the temperature of each hexagon; *<sup>T</sup>*1(*t*) = *Tin* ≡ 1 for hot inner cylinder and *<sup>T</sup>*2(*t*) = *Tout* ≡ 0 the cold outer cylinder.
