Heat transfer to proximal cylinders in hypersonic flow

Uncontrolled atmospheric entry often results in fragmentation and partial or complete ablation of the entry body. The fragments formed during the entry can interact, altering each other's aerothermal environment. Such interactions occur in meteor entries¹ and are known to affect their trajectories and ground dispersion.² Moreover, discarded spacecraft components and whole satellites³ also fragment during reentry, and their demise is also likely to be affected by the fragment interactions. Motivated by these applications, several studies of proximal body interactions in hypersonic flow have emerged in the last 15 years. These studies have focused on aerodynamic forces, and aerothermal heating has not been investigated, despite its importance in reentry body demise. However, studies^4–6 have investigated extreme localized heating due to shock-interactions, which can occur over complex surfaces or in certain proximal body arrangements. Nevertheless, these studies do not contribute significantly toward the understanding of the aerothermal environment during atmospheric entry demise.

One of the first studies on proximal body interactions by Laurence et al.⁷ experimentally and computationally investigated the forces on binary proximal cylinders and spheres in fixed positions, specifically when the second body was entirely within the primary body's shock layer. First, they developed an analytical theory based on the blast wave analogy and validated it experimentally and numerically. They found that the secondary body experiences not only a drag force but also a lift force. This lift force is exclusively attractive toward the second body if its diameter is larger than one-sixth of the primary body. This follows the intuition of the secondary body being pulled behind the primary body and becoming trapped in its wake. More recently, Marwege et al.⁸ developed a data-driven method for calculating forces in binary proximal body interactions with an accuracy within a few percent of numerical simulations.

Another early study by Laurence et al.⁹ investigated the dynamic separation of two, initially attached, spheres. They classified three main behaviors depending on the size and position of the detaching body. First, expulsion from the bow shock layer, which occurs with increasing size ratios and detachment position toward the front stagnation point. It results in the secondary body being permanently excluded from the primary body's bow shock. Second is shock surfing,¹⁰ where at a critical angle for a given radii ratio, the detached body follows the primary body's bow shock downstream. Third is entrainment in bow shock, where moving the detachment position toward the rear centerline increases the probability of the secondary body being entrained within the primary body's bow shock.

Moreover, several recent studies expand on the initial studies. For example, Park et al.¹¹ experimentally investigated the separation behavior of cubes and rough spheres. They noted that the body shape, surface roughness, and rotational motion affected the separation velocity. Furthermore, Park and Park,¹² Whalen and Laurence¹³ experimentally studied sphere clusters as large as 36 spheres and found that larger clusters separate faster than smaller ones. Similarly, Sousa et al.,¹⁴ Butler et al.¹⁵ experimentally and computationally investigated the dynamics of a spherical body shedding from a hypersonic ramp.

The literature suggests the possibility of persisting proximal body interactions during reentry, as the fragmented bodies could be trapped in the parent body's bow shock during entry. Despite this, aerothermal heating in proximal bodies during atmospheric entry has not been investigated. As a first study focusing on heat transfer, this work aims to investigate the heat transfer to proximal circular cylinders in frozen and two-dimensional hypersonic flow.

This paper is organized as follows: Sec. II reviews the relevant flow physics; Sec. III describes the governing equations and computational method; Sec. IV investigates the heat transfer with binary proximal bodies; Sec. V investigates the heat transfer to multiple proximal bodies; and finally, Sec. VI concludes the study.

Hypersonic blunt body wakes consist of a rotational inviscid shock layer, a viscous inner wake and various shock and expansion waves. Many studies have investigated blunt body hypersonic near wakes, see the introduction in Hinman and Johansen,¹⁶ whereas less is known about far hypersonic wakes.^17–19 The key features and length scales in hypersonic wakes are illustrated in Fig. 1. The bow shock around the blunt body drastically slows the flow to subsonic speeds around the front of the cylinder. A boundary layer starts from the front stagnation point and remains attached for larger turn angles compared to subsonic flows due to the favorable pressure gradient created by the expansion fan. However, adverse pressure gradients eventually form and separate the boundary layer, resulting in a separation shock—also known as a lip shock.²⁰ The separated flow converges toward the wake centerline at the reattachment point, where some flow is directed back to create a counter-rotating vortex pair near the base. With an increasing Reynolds number, the re-circulation may have more than one vortex pair. The flow that continues away from the body along the centerline creates a reattachment shock and forms a viscous inner wake which may transition to turbulence and mix with the inviscid shock layer downstream. The flow field is symmetrical along the centerline for symmetric blunt bodies, but the symmetry breaks as the wake transitions to turbulence.

Wake Reynolds number (Re_w) governs the wake topology,

where h is the wake base height, the shortest distance between the separation point and the centerline. Semi-analytical expressions for calculating the dividing streamline density (ρ_d), temperature and viscosity (μ_d) can be found in Eqs. 3.1–3.7 from Hinman and Johansen.¹⁶ Although a more convenient scaling parameter is the free-stream Reynolds number, $Re∞$⁠, which is based on free-stream flow properties and body length scale. It is not useful in describing the wake behavior as it does not account for the flow features that occur upstream of the wake. Figure 2 shows the relationship between the two Reynolds numbers and free-stream Mach number, $M∞$⁠, without thermochemical effects. It shows that generally $Rew≪Re∞$⁠, and that increasing $M∞$ for a given $Re∞$⁠, reduces Re_w due to increasing temperature and viscosity. Moreover, it outlines four flow regimes: A is diffusion-dominated, essentially Stokes flow; B—where the re-circulation region changes from diffusion to convection-dominated; C—where the wake is convection dominated and maybe unsteady; and D—where the wake is convection dominated and will transition to a turbulent wake.

Stagnation-point heat transfer is governed by the exact solution of compressible boundary layer equations.²¹ For frozen flow around a cylinder, the stagnation heat flux is

where q₀ is the stagnation-point heat flux, Pr is the Prandtl number, ρ_e is the boundary layer edge density, μ_e is the boundary layer edge viscosity, h_aw is adiabatic wall enthalpy, h_w is the wall enthalpy, and $due/dξ$ is the velocity gradient in the wall tangent direction (ξ) at the boundary layer edge. A similar expression with chemically reacting flows is not shown here as current work assumes frozen flow. The velocity gradient can be approximated using Newtonian theory,²² as in Eq. (3), where R is the nose radius, $P∞$ is the free-stream pressure, and P_e and ρ_e are boundary layer edge pressure and density.

There are many simplified forms of Eq. (2) and are reviewed by Tauber,²³ where they state a simple expression for cylinders:

They also state that the distribution of surface heat transfer is approximately

where θ is the position angle on the cylinder surface. Assuming $haw≫hw$⁠, Eq. (4) can also be written in terms of non-dimensional flow parameters, free-stream temperature and the length scale,

where C is a constant which includes free-stream fluid properties, namely, specific heats' ratio, gas constant and viscosity. Expression (6) scales with $Re∞1/2$ and $M∞5/2$ assuming the wall temperature is much colder than the stagnation temperature.

Shock interaction heating occurs when an impinging shock interacts with the body's bow shock near its surface. These interactions can be decomposed into six canonical flow patterns, as initially described by Edney.⁵ Types I, II, and V are associated with a shock-boundary-layer interaction; type III is characterized by an attaching free shear layer; type IV is characterized by an impinging or grazing supersonic jet; and type VI by an expansion-fan-boundary-layer interaction. Only type VI results in a reduction of surface gradients and pressure, all other types result in an increase in local pressure and surface gradients. Most severe heating generally occurs when a shock or a shear layer impinges on the body surface and can lead to $10×$ nominal heating or even higher.⁶ Typically, this increase in heating also correlates with an increase in surface pressure of a similar magnitude.

Semi-empirical heat transfer correlations have been developed for each type of interaction by Keyes and Hains.⁶ The exact equations are not shown here for brevity, but the correlations suggest that all types of interactions strongly depend on the impinging and bow shock angles at the intersection; as well as the local surface inclination and interaction-specific flow length scale.

This work uses a Navier–Stokes computational fluid dynamics solver, originally designed for modeling reacting and dense compressible fluids.^24,25 More recently, the solver has been further developed with a focus on modeling hypersonic flows.²⁶ It has been coupled with BoxLib,²⁷ a structured Cartesian adaptive mesh refinement massively parallel framework and a ghost-point forcing immersed boundary method. The solver has been validated over a range of test cases with a focus on surface properties predictions (pressure, velocity gradients, and temperature gradients). The tests include stand-off and surface pressure predictions up to sub-orbital Mach numbers (⁠$M ≲ 30$⁠) in Euler flow. In addition to pressure, heat transfer and skin friction predictions over compression ramp with shock–shock and shock–boundary layer interactions in viscous flow. The solver is well suited for resolving a large range of flow scales, often found in hypersonic flows, and modeling flows around complex geometries.

The solver computes the Navier–Stokes equations conservation laws as follows:

where ρ is the density, u is the velocity vector, p is the pressure, e_t is the total energy, I is the identity matrix, $∇$ is the gradient operator, h_t is the total enthalpy, $τ$ is the shear stress tensor, μ is the dynamic viscosity, q is the thermal conduction flux vector, and λ is the thermal conductivity. The Euler fluxes are calculated via a fourth-order accurate central-skew-like conservative difference method with artificial dissipation shock-capturing.²⁸ 

The viscous fluxes are computed with second-order accurate standard central difference. The time integration is with a local time-stepping strategy in the AMR framework—where the time step is halved, and the number of steps doubled, for each refinement level to maintain a constant Courant–Friedrichs–Lewy (CFL) number on all refinement levels in uniform flow. Equation (7) is time-marched with an explicit third-order total variation diminishing Runge–Kutta scheme.²⁹ In closing Eq. (7), the following assumptions are made: N₂ gas free-stream, frozen flow, Sutherland's law viscosity, constant specific heat capacity (c_p = 1040 kJ/kg K), and constant Prandtl number (Pr = 0.7).

Previous numerical studies focusing on proximal body separation^7,9 have used computational methods similar to the current work: combining a high-order shock-capturing method with embedded ghost-cell boundary conditions and parallel dynamically adaptive structured Cartesian meshes.

An expression for the transverse velocity (V_t) of proximal spherical bodies was first presented by Passey and Melosh,² which was later verified further by dynamical separation studies of Laurence et al.⁹ and Park et al.¹¹ The expression reads

where C is a constant, r₁ and r₂ are radii of the two spherical bodies, $ρ∞$ is free-stream density, $V∞$ is free-stream velocity, and ρ_b is the density of the detached body. The order of magnitudes for each variable considering Earth entries are as follows. The constant C, as discussed by Laurence et al.,⁹ is $10−2<O(C)<1$⁠. The number of fragments are assumed to be small, so O(n) = 1. The fragmented bodies are assumed to be smaller than the primary body but not too small, which means $10−1<(r1/r2)<1$⁠. Finally, the ratio of the free-stream density to body density is not trivial. Spacecraft components and satellites are usually metallic, and common meteorites³⁰ are stony, stony-iron, or iron. A reasonable assumption is to consider a metallic body, so $O(ρb)=103$ kg/m³. Passey and Melosh² suggest that breakup altitudes can be as low as 6 km and as high as the mesosphere (approximately 80 km), so $10−5<O(ρ∞)<1$ kg/m³. Therefore, $10−6<O(Vt/V∞)<10−1$⁠; this suggests that the separation velocity varies over orders of magnitude and is much slower than free-stream velocity. Therefore, the flow adjusts very quickly around the bodies and neglecting body dynamics can be justified as a first approximation.

In this work, two identical bodies with 0.1 m diameter are considered. The body positions are non-dimensionalized using flow length scales, instead of the body length scale, as it allows a meaningful comparison of positions across different Mach numbers. The positions in real space (x, y) and non-dimensional space (⁠$x̂,ŷ$⁠) are illustrated in Fig. 3. The figure also shows a hatched region representing a collision area where the second body cannot be placed. The second body's position is systematically varied in the primary body's near wake, $0≤x̂≤2$ and $0≤ŷ≤2$⁠. The reference length scale in the y direction is the shock radius (y_s) and in the x direction is the rear stagnation-point distance (l_s) from the body center. The different positions are investigated at Mach numbers 2, 4, and 8. These Mach numbers are lower than typical atmospheric entry flows, which could be up to Mach 30, where thermochemical relaxation and radiative heat transfer have a significant effect on the overall heat transfer. Including these physical phenomena is beyond the scope of the current work, and therefore, the conditions investigated are representative of low altitude conditions. The free-stream Reynolds number is fixed at 10 000, and the free-stream temperature is 300 K. Therefore, the flow is completely characterized by the Mach number. The wall temperature is 300 K for Mach 2 and 500 K for the other two conditions. The wake Reynolds numbers are less than 100 for the free-stream Reynolds and Mach numbers, suggested by Fig. 2, and laminar flow behavior is expected.

A mesh refinement study is conducted at the three Mach numbers, from which the meshes are summarized in Table I. The Mach 2 mesh has more number of computational points than the Mach 4 because a larger domain is required due to the wider shock layer. The mesh requirement for accurate stagnation-point heat transfer prediction can be estimated using Eq. (4). Figures 4(a) and 4(b) show the mesh in a fully developed flow field, with refinement around shock and cylinder boundary at M = 4, and with refinement around the cylinder boundary only at M = 8. The mesh for M = 2 is similar and not shown for brevity. At M = 8, the higher levels of refinement (3, 4, and 5) are added around the front of the cylinder using criteria based on the wall pressure coefficient. The simulations are initialized at the free-stream temperature and pressure with free-stream Mach number everywhere except for a small annulus around the cylinder, where the Mach number is zero. The domain boundary conditions, considering Fig. 4, are free-stream inflow on the left boundary and zero-gradient for all other boundaries. For the immersed boundary, no-slip, constant wall temperature and zero pressure gradient conditions are imposed. The simulations are computed for 10 flow times, where a flow time is defined as the domain length divided by the free-stream velocity. The maximum CFL number is around 0.5. Re-gridding occurs every 0.5 flow times. The flow length scales required for the non-dimensionalization are computed from the mesh refinement results and are listed in  Appendix A.

The non-dimensional surface properties, pressure (C_p), skin friction (C_f), and heat transfer (C_h) coefficients are defined as follows:

Figure 5 shows surface properties down-sampled to 100 surface elements and averaged with a three-point moving filter. The pressure field does not vary significantly with increasing resolution at all Mach numbers. The reference (Ref) heat transfer profile is from Eq. (5), which are known to include around 25% error when compared to ground shock-tube and ballistic range experiments.²³ The stagnation-point heat transfer predictions from all meshes are within this error. Table II suggests that at M = 2, the stagnation heat transfer is overpredicted by 10%. At M = 4 and 8, the stagnation heat transfer is underpredicted by around 10%. Another observation is that the heat transfer profile qualitatively agrees but is offset around a factor of the stagnation-point heat transfer error, and it does not agree well beyond $3π/8$ as expected. For all Mach numbers, the medium mesh resolution gives results within 3% of the fine mesh results; therefore, it is deemed acceptable for binary bodies computations.

Initial studies suggested secondary body arrangements with $ŷ=2$ were too far to observe any effect on each other. So, despite the illustration in Fig. 3, the study is limited to $ŷ≤1$⁠. Eight arrangements of the bodies are considered, in (⁠$x̂,ŷ$⁠) coordinates and in order from 1 to 8 are {(1,0), (2,0), (0,0.5), (1,0.5), (2,0.5), (0,1), (0.5,1), and (1,1)}. Simulations are conducted using these arrangements and the numerical setup is identical to the medium mesh simulations from the refinement study in Sec. IV A, but with two immersed bodies. The meshing criteria around both bodies are the same. The simulations are advanced for 10 flow-through times; this duration is sufficient for the surface gradients to be developed close to a steady state, as shown in Fig. 6. The figure shows the heat load variation of primary and secondary bodies with non-dimensional time. The heat load is more unsteady at lower Mach numbers but only for a couple of arrangements. Specifically, at M = 2 in arrangements 2 and 4, and at M = 4 in arrangements 2 and 5. Therefore, the surface properties for these cases are time-averaged over 5–10 flow times, but for all other cases, the surface properties are taken at their steady state.

The typical mesh contains around $[4.6,1.5,5.6]×106$ points for Mach numbers [2, 4, 8]. In all the meshes, more than 80% of computational points are on the base level. Each binary body arrangement costs around [980, 170, 1590] CPU hours to compute. The M = 2 mesh is almost the same size as the M = 8 mesh as the lower Mach number domain is twice as large and the base mesh resolution is the same for both Mach numbers.

Figures 7–9 show similar flow features to Fig. 1, namely bow shock, lip shock, and reattachment shocks. In arrangements 1 and 2 for all Mach numbers, the secondary body does not have a bow shock when directly behind the primary body. In arrangements 3–8, both bodies have interacting bow shocks. In arrangements 3, 6, and 7 at M = 2, 3, and 6 for M = 4 and only in arrangement 6 at M = 8; the bow shocks combine and form a single smooth bow shock, or they result in interactions which can lead to unsteady wakes with oscillating reattachment shocks and shocklet shedding. Counter-intuitively, the wake becomes more steady at increasing Mach numbers. This is because the wake behavior is governed by the wake Reynolds number which decreases with increasing free-stream Mach number for a fixed free-stream Reynolds number, as shown in Fig. 2.

The results show Edney-type shock interactions but mostly far from the body surface. They are known to occur with near-body shock interactions, as discussed in Sec. II B. However, most of the binary body arrangements are such that the secondary body is positioned fully inside the shock layer of the primary body. This is the case at M = 2 and M = 4 for all arrangements; and at M = 8 with arrangements 1, 2, 6, and 7. Therefore, in these arrangements, there is no impinging shock and bow shock interaction. Arrangements 3, 4, 5, and 8 at M = 8 are close to producing Edney-type interactions. However, the secondary body positions are such that shock interactions occur far from the body. Small changes in relative body position may produce Edney type IV or type I.

Figures 10–12 show the surface pressure, skin friction and heat transfer coefficients. Both C_p and C_f profiles imply the forces acting on the bodies, and C_h implies the total heat load. The results show that:

• Often C_p is symmetrical, which means there is no net lift and only a drag force on the body; whereas an asymmetrical C_p profile suggests the presence of a lift force and a drag force. Furthermore, the C_p distribution is smooth, peaks at the stagnation point, and is close to zero behind the cylinder.

• Often C_f is antisymmetric. This suggests zero net torque acting on the body and only drag force due to viscosity. C_f magnitude is $100×$ smaller than C_p and its contribution to drag is negligible. Whereas an asymmetrical profile suggests a net torque acting on the body. Moreover, C_f tends to zero at the stagnation point and peaks when surface tangents are aligned with the free-stream, and its local maxima and minima are due to boundary layer separations and attachments.

• A symmetrical C_h profile around the body centerline suggests similar heating on both sides of the cylinder, whereas, in an asymmetrical profile, one side is heated more than the other. C_h peaks at the stagnation point and is negligible on a surface with zero projected area (shadow surfaces) in the flow direction. The profile's asymmetry increases with increasing Mach number while the peak heat transfer and heat load are also visibly different with the arrangements.

Figures 13(a) and 13(c) show the normalized heat load flux and the peak heat flux on both bodies in all arrangements. The results are conservatively presented to two significant figures considering mesh convergence in Sec. IV A. Overall, the normalized heat load and peak heat flux can increase up to 1.2 and decrease to 0.1. Generally, the magnitude of variation in heat load and peak heat augmentation is observed to be insensitive to the Mach number. A similar magnitude of variation is also likely to be the case for even higher Mach numbers, and thermochemically frozen flow, as the wake topology remains similar.

When the secondary body is directly behind the primary body, it is significantly thermally shielded. The shielding is due to a reduction in heat transfer on the front half of the cylinder (⁠$−π/2<θ<π/2$⁠). This occurs in arrangements 1 and 2 for all Mach numbers, but can also be observed in arrangements 4 and 5 at M = 8. The heat transfer is reduced on the front half of the cylinder due to a weaker flow impingement around the stagnation region because of shielding from the leading body. This pressure shielding can also be observed in arrangement 5 at M = 8.

On the other hand, the increase in heat load is due to an increase in heat transfer on the sides of the cylinder (⁠$π/4<|θ|<3π/4$⁠), because of the hot shock layer flow being forced between the two bodies. As the secondary body moves toward the primary body's bow shock, as in arrangements 3 and 6, the heat transfer to both bodies increases.

Another observation is that arrangements 1 and 2 at M = 2 generate larger heat loads and peak heating than at M = 4 and M = 8. This is because of mixing in the shock layer due to a higher wake Reynolds number at M = 2, the mixing was also observed with the numerical Schlieren in Sec. IV B 1.

Following the investigation of heat transfer to binary proximal bodies, a natural first step is to consider multiple bodies, also referred to as a cluster. Many of the parameters which affect the bodies' surface properties are the same as in binary bodies, namely the flow Mach number and Reynolds number. However, with multiple bodies, the number of bodies in a cluster and the size and position of the bodies widen the parameter space. In this study, the focus is only on the effect of body positions on heat transfer. Only one flow condition with free-stream at M = 4 and Re = 10 000; and clusters of five equal-sized bodies at the specified flow conditions are considered.

The clusters are generated by randomly placing cylinder centers while ensuring no overlap with each other, and that the bodies are located within the computational domain. Eight positions are selected, see arrangements in  Appendix B. The variety in the spacial arrangement of the bodies in the cluster is ensured: cases 7 and 8 are closely arranged; 1 and 2 are medium distance apart; 5 and 6 are loosely arranged; and 3 and 4 are arranged roughly linearly—horizontally and vertically. The numerical setup is identical to the M = 4 binary body study in Sec. IV B. Figure 14 shows the instantaneous mesh for cases 2 and 5.

Let $hij(θ)$ be the local time varying heat transfer for body i in a given cluster at time j,

where N_b is the number of bodies in a cluster, N_t is the number of samples at different times, Q_ij is the heat load for body i and time j, $Q¯i$ is the time-averaged heat load for body i in a cluster, $⟨Q¯⟩$ is the time- and cluster-averaged heat load and $σ{Q¯}$ is the standard deviation of the time-averaged heat load over bodies in a cluster. The heat transfer coefficient C_h is averaged similarly.

Figure 15 shows the instantaneous numerical Schlieren image from the eight cases, where the five bodies are labeled according to Table IV. The primary bow shock shape varies dramatically depending on the body arrangements, and highly unsteady wakes are generated with shock interactions and shocklet shedding. The unsteadiness is significant compared to binary body arrangements in Fig. 8 at the same free-stream conditions.

Figure 16 shows the cluster-averaged unsteady evolution heat load and heat load for each body. The flow field observed in Fig. 15 results in unsteady heat transfer to many bodies. A time-average over five flow times is deemed to be a reasonable approximation of the steady-state heat load. However, in some cases, such as Body D in case 5, the current sampling frequency of 0.5 flow time may not be sufficient for an accurate time-average. This is unlikely to change the qualitative trends.

Figure 16 indicates that the maximum time-averaged heat load on a given body is around 1.4 (body 8A), compared to a factor of 1.2 in the binary bodies configuration. While the minimum heat load is around 0.1 (body 3B), same as in the binary bodies. The figure also shows that some bodies experience significant unsteadiness. Case 5 has the largest heat load fluctuations (0.2 peak-to-peak). In contrast, all bodies in cases 4 and 7 experience a steady heat load. In the remaining cases—1, 2, 3, 6, and 8—only some bodies experience an unsteady heat load. Bodies C and E in case 6 experience unsteady heat load whereas bodies A, B, and E undergo steady heating. Finally, an important observation is that the normalized cluster and time-averaged heat load approximately varies between 1.2 and 0.4.

Figure 17 shows the time-averaged heat transfer coefficient for all cases. The normalized maximum peak heat flux is around 1.4 larger (on body 2C), and the minimum is 0.2 (case 2E). This suggests a large variation in peak heating occurs on bodies in a given cluster. Comparatively, binary body configurations produced peak heating in the range of 1.2–0.1. All peaks occur between $−π/2<θ<π/2$⁠, on the front half of the cylinders, as expected. Large spatial variations occur in case 8; whereas cases 1, 2, 4, and 7 have medium variations; and 3 and 6 have small variations. Overall, the additional heating on the sides and back of the cylinders is the cause of the increased heat load.

Figure 18 shows the time- and cluster-averaged coefficient of heat transfer, $⟨C¯h⟩$⁠, profiles around the cylinder for all the cases considered. All cases have their average stagnation-point heat transfer reduced by around 10%–50% compared to an isolated body. This is because some bodies are shielded from the direct shocked flow, for example, some bodies in cases 1, 2, 3, and 7. These shielded bodies have a reduced heat load around the stagnation region. However, the heat transfer around the sides and the back, for $π/4 ≲ |θ|<π$⁠, is increased in many cases. This increase not only offsets the reduction in heat load around the stagnation region but also results in a net increased heat load when compared to an isolated body.

Figure 19 shows the magnitude of normalized cluster-average heat load, which can also be inferred from Fig. 16. It shows that bodies on average, in cases 4, 5 and 8, receive a larger heat load than an isolated body (⁠$⟨Q¯⟩/Q1>1$⁠). In case 5, the unsteadiness and mixing in the wake are responsible for the increased heat load, as suggested by Fig. 16. However, in case 4, where the flow is steady, additional heating is due to the flow being forced to pass between the cylinders—as suggested by Figs. 16 and 17. In contrast, case 8 has increased heating on the side and back, as the hot gases from the stagnation region of the combined bow shock linger in the small spaces in between the bodies—as suggested by Fig. 17. Overall, the increase in time- and cluster-averaged heat load seems to be caused by different mechanisms. On the other hand, cases 1, 2, 3, and 7 are subject to a reduced heat load (⁠$⟨Q¯⟩/Q1<1$⁠). This occurs because some bodies are shielded from the direct shock layer flow and are in the cooler viscous wake of another body, as seen in Fig. 15. Such shielded bodies have reduced heat transfer around their stagnation region, as shown in Fig. 18. Moreover, the heat load in case 6 is on average the same as an isolated body (⁠$⟨Q¯⟩/Q1≈1$⁠), as the increase in heat transfer to the leading bodies is negated by the reduced heat transfer to the trailing bodies. Finally, bodies in cases 7 and 8 are arranged in a visually similar manner, but they differ in their normalized average heat load: 1.2 and 0.85, respectively. Figures 16 and 17 suggest that case 8 is more unsteady and receives larger heating on the sides and back than case 7.

Figure 19 also shows the scaling of normalized time- and cluster-averaged heat load with cluster geometrical properties. Simple properties like the standard deviation of the body positions in both coordinate directions give good correlations with current cases. More complex properties like projected area in the flow direction, and for a collection of points (body centers) Hamiltonian path, Minkowski distance, did not give improved correlations. The standard deviation of the coordinate parallel to the flow vector (⁠$σ{x}$⁠) and perpendicular to the flow vector (⁠$σ{y}$⁠) correlate well with the time- and cluster-averaged heat load. However, the ratio $σ{x}/σ{y}$ correlates the best.

Overall, the correlations suggest that increasing $σ{x}$ decreases the time- and cluster-averaged heat load. In contrast, $σ{y}$ is directly proportional to time- and cluster-averaged heat load. Importantly, it suggests that when $σ{x}/σ{y}≳1$⁠, the average heat load is less than compared to an isolated body; in other words, thermal shielding occurs. Alternatively, when $σ{x}/σ{y}≲1$⁠, the average heat load is more than an isolated body; in other words, thermal amplification occurs.

Note, these correlations are only valid for closely packed clusters with $σ{x} ≲ 3D$ and $σ{y} ≲ 3D$⁠. For $σ{y}≫3D$⁠, meaning the bodies in the cluster are far apart and are essentially isolated. Therefore, the current linearly increasing heat load behavior cannot be sustained for all values of $σ{y}/D$⁠.

Figure 20 shows the time-averaged heat load standard deviation in a cluster, correlated with four parameters. The strongest (negative) correlation is observed with $σ{y}$ and the weakest positive correlation is with $σ{x}$⁠. Note, increasing $σ{y}$ cannot decrease $σ{Q¯}$ beyond 0. Intuitively, the heat load variance must go to zero as the cluster size increases and the bodies become isolated. This means that for $σ{y}≳3D$⁠, there must be a change in trend. Moreover, $σ{x}/σ{y}$ shows a positive correlation; and $⟨Q¯⟩/Q1$ shows a negative correlation. This means that heated clusters have a smaller variation of heat load among their bodies than cooled clusters.

This paper studied the heat transfer in proximal cylinders in the near wake at Mach numbers 2, 4, and 8, with constant wall temperature, under the calorically perfect gas assumption. First, binary bodies are considered, where the secondary body's positions are systematically varied in the primary body's near wake. The key findings are as follows:

• The heat load and peak heat transfer can be augmented for either one or both proximal bodies by $+20%$ to −90% of an isolated body. The magnitude of variation in the heat load and the peak heat transfer seems insensitive to the Mach numbers investigated. Moreover, this trend is expected to continue in flows with larger Mach numbers than investigated, in frozen flow and given a similar free-stream Reynolds number.

• Minimum heat load and minimum peak heat transfer (maximum thermal shielding) occur when the secondary body is in the viscous wake of the primary body. The thermal shielding is due to a reduction in heat transfer to the front half of the cylinder due to a weaker flow impingement, because of flow shielding from the leading body.

• Maximum heat load and peak heat flux occur when the secondary body approaches the primary body's bow shock and stagnation region. The increase in heat load is due to increased heat transfer around the sides of the cylinder (⁠$π/4<|θ|<3π/4$⁠) due to the hot shock layer gases being forced to pass between the two bodies.

Second, multiple bodies near each other are considered. Eight clusters with five randomly arranged cylinders are selected, and the cluster-average heat transfer and standard deviation are described and found to correlate well with cluster geometrical properties. The key findings are:

• The maximum time-averaged heat load on a given body in a cluster can vary between around $+40%$ to −90% compared to an isolated body. The increase in heat load is mainly from additional heating to the sides and the back of the cylinders, whereas the decrease is mainly due to reduced heat transfer on the front half of the cylinder. Moreover, flow fields can be highly unsteady, resulting in unsteady heat transfers.

• The time- and cluster-averaged heat load in a cluster varies between +20% and −60% of an isolated body. The change in heat load among the clusters is mainly due to the change in heat transfer around the sides and the back. The average heat load around the stagnation region is also reduced in most clusters, as the following bodies are shielded by the leading bodies from the direct flow. The increased heat transfer to the sides and back is caused by: unsteadiness and mixing in the wake; the flow being forced in between the cylinders; or the hot flow from the stagnation region lingering in the small spaces in between the bodies.

• The time- and cluster-averaged heat load in a cluster shows a strong negative correlation with the ratio of the standard deviation of body coordinates in the direction parallel to the free-stream velocity to the deviation of body coordinates in the perpendicular direction. In other words, clusters thin in the direction perpendicular to the free-stream velocity and long in the direction parallel to the free-stream velocity have their heat load reduced (thermally shielded) while thick and short clusters have an increased heat load.

• The time-averaged heat load standard deviation in a cluster shows a negative correlation with the magnitude of the time- and cluster-averaged heat load. This means that heated clusters have a smaller variation of heat load among their bodies than cooled clusters.

Despite the current limitation of this study to two-dimensional flows, the trends and intuition gained in the current work are likely to apply to three-dimensional flows, at least qualitatively. Future works should investigate heat transfer to proximal bodies in three-dimensions and the effect of thermochemical relaxation, which could have large effects on heat transfer. Also, the effect of the wake Reynolds number on the heat transfer, particularly in the far wake, should also be investigated.

The authors would like to thank and acknowledge P. Bruce for his comments on the work. The authors also acknowledge Engineering and Physical Sciences Research Council (EPSRC) (No. EPSRC-DTP 2018/19) and Mechanical Engineering Department at Imperial College London for the funding.

The authors have no conflicts to disclose.

Monal Patel: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (lead). Salvador Navarro-Martinez: Funding acquisition (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Flow length scales, estimated from the mesh refinement study, are summarized in Table III. These values are used to calculate the positions of the second cylinder in real space from non-dimensional space, as illustrated in Fig. 3.

The body positions for each case in Sec. V are listed in Table IV.