Geometrical investigation of forced convective 昀氀ows over staggered arrangement of cylinders employing constructal design

The current numerical work conducted a geometrical investigation of forced convective 昀氀ows over a staggered arrangement of four cylinders applying constructal design. The main objectives are to determine designs maximizing the Nusselt number ( Nu D ) and minimizing the drag coef昀椀cient ( C D ) of the cylinders at three distinct Reynolds numbers ( Re D = 10, 40, and 150) considering Pr = 0.71. The problem involves two constraints: cross-sectional area of the cylinders and the occupation area of the arrangement, and three degrees of freedom: S T / D , S L1 / D and S L2 / D . These represent the ratios between the transversal pitch, the longitudinal pitch between the frontal and intermediate cylinders, and the longitudinal pitch between the intermediate and posterior cylinders over the diameter. The introduction of two longitudinal pitches enabled an exploration of the generation of asymmetrical patterns in the 昀氀ow direction. Results indicated a strong dependency of the arrangement design on the Reynolds number. Moreover, asymmetrical con昀椀gurations generally bene昀椀ted both performance indicators for the present conditions. The most favorable arrangements for the 昀氀uid dynamics were achieved with the lowest S T / D , more concentrated longitudinal pitches for Re D = 10, and asymmetrical pitches for Re D = 40 and 150. Concerning the thermal purpose, the highest magnitude of S T / D and asymmetrical longitudinal pitches led to the highest Nu D . For Re D = 10, the intermediate cylinders are near the upstream one; for Re D = 40 the length of the arrangement increases with the intermediate cylinders near the downstream one, and for Re D = 150, the arrangement occupied all the length of occupation area with slight asymmetry.


Introduction
Convection heat transfer 昀氀ows over arrangement of tubes or cylinders has attracted particular interest due to its accurate representation of various industrial and commercial equipment, including heat exchangers, evaporators, condensers, and other engineering applications as electronic packaging, gas pipelines, and air-conditioning units [1][2][3].Moreover, the complex structures generated in the 昀氀ow over a single cylinder, such as boundary layer detachment, vortex shedding and formation of vortex street, along with the interaction of multiple structures in an arrangement of cylinders, pose a challenging issue for the accurate prediction of their design [4][5][6].
Several studies have been reported in the literature aiming to de昀椀ne strategies to enhance the heat transfer between one cylinder and the surrounding 昀氀ow.These strategies include the consideration of protrusions [7], the incorporation of multiple 昀椀ns in the cylinder periphery [8], or the use of 昀氀exible 昀椀ns [9] to passively control the 昀氀ow.These studies also focus on the understanding the 昀氀ow behavior once the use of protrusions or 昀椀ns dramatically alters the 昀氀ow regime, leading to changes in the vortex shedding that can result in synchronization with the natural frequency of the material (lock-on phenomenon) [8].In addition to circular cylinders and those with protrusions or 昀椀ns, other geometric con昀椀gurations, such as square/rectangular shapes, have also been extensively studied [10].Alternative con昀椀gurations as semicircular cylinder [11], cam-shaped cylinder [12], and elliptical and wave elliptical cylinders [13,14] have also been reported.A comprehensive review on the in昀氀uence of different geometrical con昀椀gurations for a single cylinder over the behavior of 昀氀uid 昀氀ow is described in the work of Lekkala et al. [15].Another strategy to augment the thermal exchange involves placing the cylinder in an enclosure, such as an expansion in a channel [16].Recent advances have been done on the use of magnetohydrodynamics to control the 昀氀uid 昀氀ow and augment the heat transfer rate between the cylinder and the 昀氀uid 昀氀ow [17].Studies on two tubes or bluff bodies under different con昀椀gurations (in line, sideby-side or inclined) had also been reported in the literature.For instance, the work of Sumner [3] provides a detailed explanation about the behavior of the 昀氀uid 昀氀ow over a pair of cylinders under several conditions.In general, the literature has demonstrated that the studies of convective 昀氀ows over single or pair of cylinders remains a relevant issue.
Regarding advancements in the study of cylinder arrangements, investigations have been performed ranging from single rows to staggered and aligned con昀椀gurations with three cylinders or more [18][19][20][21][22].For instance, Ambesi and Kleijn [18] numerically studied forced convection laminar 昀氀ows over a perpendicular single row of cylinders with equidistant and non-equidistant spacing.Authors explored the array's impact on the average Nusselt number for different Reynolds (0.001 ≤ Re D ≤ 600) and Prandtl numbers (0.7 ≤ Pr ≤ 10).Authors noticed a decrease in the average Nusselt number of the array with the decrease of the open frontal area for low Reynolds numbers and equidistant cylinders.Furthermore, for the same frontal area, the non-equidistant array led to lower Nusselt numbers for intermediate Reynolds numbers, with no signi昀椀cant impact on heat transfer rates for low and high Reynolds numbers.Subsequently, Daróczy et al. [19] used computational intelligence to search for the best positions of seven-cylinders arrangement, aiming to maximize the heat transfer rate between the cylinders and the surrounding 昀氀ow.Recently, Pravesh et al. [20] investigated mixed convective non-Newtonian 昀氀ows over inline array of cylinders.The study considered different numbers of Reynolds (1 ≤ Re D ≤ 40) and Prandtl (1 ≤ Pr ≤ 50), as well as varied power-law indexes (0.4 ≤ n ≤ 1.8).The con昀椀guration was assessed using a 昀氀uid volume fraction parameter, representing the 昀氀uid area over the total square area de昀椀ned by (L + D/2) 2 where L is the distance among the cylinders and D the diameter.The investigation revealed an increase in the average Nusselt number with Reynolds, Prandtl and Richardson numbers for shearthinning 昀氀uids (n tending to 0.4), while the opposite was observed for shear-thickening 昀氀uids (n tending to 1.8).Additionally, an increase in volume fraction resulted in a higher Nusselt number for mixed convective 昀氀ows, contrary to observations in forced convective 昀氀ows.Xue et al. [21] examined forced convective laminar 昀氀ows over a bundle with 10 tubes using an equilateral triangle arrangement.The authors evaluated the in昀氀uence of the Reynolds number, longitudinal and transversal pitches over the tube-averaged heat transfer and average friction co-ef昀椀cients.The study focused on varying one parameter while keeping the other two constant, i.e., it is not performed a combination of different pitches.Results indicated that the transversal pitch and Reynolds number had the most signi昀椀cant impact on the thermal and 昀氀uid dynamics performance.Finally, Islam et al. [22] made a numerical study investigating the in昀氀uence of vibration of internal line of cylinders in an initially aligned arrangement with 3 × 3 cylinders.The study also explored the impact of transverse spacing ratio and the diameter on parameters such as mean pressure coef昀椀cient and Nusselt number.The best thermal performance was achieved when the varying cylinders have the lowest transverse spacing investigated, compared to arrangements with the same diameter.
As a 昀椀nite size 昀氀ow system, the design in convection heat transfer over arrangements of cylinders can be guided by the constructal law of design and evolution [23][24][25].The constructal law is a physical principle that posits the design and rhythm in any 昀椀nite-size 昀氀ow system, whether animate or inanimate, evolves toward facilitating the internal currents that 昀氀ow through the thermodynamical system.The constructal design, which is the method of application of the constructal law, has been utilized to demonstrate the designs in nature, society, and engineering applications [23][24][25][26].In the external convection heat transfer 昀氀ows, the method has been applied to investigate the geometry of 昀椀ns that best remove heat from a basement [27], and to direct the con昀椀guration of aligned arrangement of cylinders subjected to natural convection [28], arrays of cylinders in rotation with different sizes under forced convection [29], triangular arrangement of cylinders, and bluff bodies under mixed and forced convection laminar 昀氀ows [30,31].Recent investigations have sought to understand the in昀氀uence of the triangular arrangement of bluff bodies on drag coef昀椀cient and Nusselt number, considering turbulent mixed convection 昀氀ows with Re D = 22,000 and Pr = 0.71 [32].Additionally, a construction function based on velocity and thermal 昀椀elds has been employed to create an arrangement from an elementary con昀椀guration (with one tube) for forced convection laminar 昀氀ows for three different Reynolds numbers (Re D = 10, 50 and 100) and Pr = 0.71 [33].Finally, Cunegatto et al. [34] associated the constructal design and response surface methodology to study the in昀氀uence of arrangements of tubes with different sizes on heat transfer density for convection heat transfer 昀氀ows and with pseudoplastic 昀氀uids.
The literature review evidenced several research works where the constructal design was successfully applied to direct the con昀椀guration of cylinders subjected to convective 昀氀ows.Nevertheless, none of the previously mentioned investigations considered the in昀氀uence of two different longitudinal pitches for the elementary lozenge con昀椀guration of the arrangement of cylinders over both the 昀氀uid dynamic and thermal performance.The present work aims to contribute in this regard by employing the constructal design method to investigate the in昀氀uence of a staggered arrangement of four cylinders on the average drag coef昀椀cient (C D ) and average Nusselt number (Nu D ) in forced convective laminar 昀氀ows.To explore the impact of variable longitudinal pitches in the arrangement design, the problem is de昀椀ned with three degrees of freedom: S T /D (the ratio between the transversal pitch of intermediate cylinders and the diameter of the cylinders), S L1 /D (the ratio between the longitudinal pitch of upstream and intermediate cylinders and the diameter of the cylinders), and S L2 /D (ratio between the longitudinal pitch of the intermediate cylinders and the downstream one and the diameter of the cylinders).Additionally, the effect of the Reynolds number on the best 昀氀uid dynamic and thermal con昀椀gurations is investigated.

Problem description and governing equations
The investigated problem involves an incompressible, laminar, forced convection 昀氀ow of air (Pr = 0.71) over a staggered arrangement of cylinders in a two-dimensional domain, as illustrated in Fig. 1.The domain of Fig. 1 has the following constant dimensions: D = 0.1 m, H = 15D = 1.5 m, L = 35D = 3.5 m, L 1 = 8D = 0.8 m.The 昀氀uid 昀氀ow is analyzed in the unsteady state, with C D and Nu D evaluations performed upon reaching a steady-state condition.The 昀氀ow is induced by a prescribed velocity 昀椀eld at the inlet of the domain, and heat transfer occurs due to the temperature difference between the isothermal cylinders at a higher temperature (T w = 320 K) and the cooler surrounding stream (T ∞ = 300 K).The speci昀椀ed boundary conditions include a constant velocity (u = V ∞ ) with three potential values based on the Reynolds number, and a temperature (T = T ∞ ) imposed on the left-lateral surface of the domain.The four cylinders exhibit constant velocity (u = v = 0 m/s) and temperature (T = T w ).Symmetry conditions are applied to the upper and lower surfaces of the domain, while a constant gauge pressure (p g = 0 Pa) is maintained at the exit, represented by the right surface of the domain.It is worth mentioning that the domain's right surface is placed far enough from the four-cylinder array not to affect the results.Regarding initial conditions, the 昀氀uid is initially at rest (u = v = 0 m/s, T = T ∞ ) at the starting time (t = 0 s).
The modeling of the forced convection problem is given by the following equations of conservation of mass, momentum in x and y directions, and conservation of energy [35]: ∂u ∂t ∂v ∂t ∂T ∂t where x and y are the spatial coordinates (m), u and v are the velocities in x and y directions (m/s), t is the time (s), ρ is the density (kg/m 3 ), p is the pressure (Pa), υ is the kinematic viscosity (m 2 /s), T is the temperature (K) and α is the thermal diffusivity (m 2 /s).

Geometrical investigation with constructal design
The application of the constructal design method in this study is illustrated in Figs. 2 and 3.In the initial steps, Fig. 2, it is crucial to de昀椀ne the 昀氀ow system, including the domain, thermophysical properties, and boundary and initial conditions.Understanding the physical problem is essential to properly de昀椀ne how to facilitate the access to the internal currents (represented by maximizing/minimizing the performance indicators), the constraints, and the degrees of freedom (steps 1 to 5).Simultaneously, it is important to establish the method for solving the physical problem, which involves the numerical solution of Eqs. ( 1)-(4) (step 5b).
In step 3, the objectives are maximizing the average Nusselt number (Nu D ) and minimizing the average drag coef昀椀cient (C D ) in the array.The average Nu D is given by [35]: where the average Nusselt number for each cylinder is given by [36]: And the average C D of the array is given by [1][2]: where i indicates the number of the cylinder (1 ≤ i ≤ 4), * represents the dimensionless variables, n is the normal surface to the cylinder, θ is the angle of the cylinder (rad) where the local Nusselt number is integrated, and F' D is the drag force per unit depth (N/m).
In step 4, two constraints are de昀椀ned for the problem: -the cross-sectional areas of the cylinders -the occupation area of the array where It is worth noting that it is used the similarity principle to ease the convergence of the numerical solution without any loss to the accuracy of the calculus of drag coef昀椀cient and Nusselt number.Consequently, to comprehensively address the problem, a total of 288 different cases were numerically simulated.
Fig. 3 provides a detailed depiction of the three-level 昀氀owchart outlining the geometrical investigation for each Reynolds number (step 6).In the 昀椀rst level, the ratio S T /D is varied while keeping the ratios S L1 / D and S L2 /D constant.For the thermal purpose, the highest Nu D , is the once maximized Nusselt number, (Nu D ) m , and the corresponding S T /D is the once optimized ratio of S T /D for the thermal purpose, (S T /D) o,T .Similarly, for the 昀氀uid dynamic purpose, the lowest C D is the once minimized C D , (C D ) m , and the corresponding S T /D is the once optimized ratio of S T /D for the 昀氀uid dynamic purpose, (S T /D) o,F .For the second level, the process is repeated for different ratios of S L2 /D while maintaining the ratio S L1 /D constant.The maximum Nu D obtained is labeled as (Nu D ) 2m , and the corresponding optimal con昀椀gurations are (S L2 /D) o,T and (S T /D) 2o,T .Simultaneously, the minimum C D is denoted as (C D ) 2m , and the corresponding optimum ratios are (S L2 /D) o,F and (S T /D) 2o,F .The same process is applied in the third level, now varying the ratios of S L1 / D. The maximum Nu D is labeled (Nu D ) 3m , the minimum C D is denoted (C D ) 3m .The optimal ratios are presented in Fig. 3. Additionally, Fig. 3 visually represents potential optimal thermal and 昀氀uid dynamic paths, highlighted in red and blue, respectively.
After de昀椀ning the search space for the geometrical investigation, Eqs. ( 1)-( 4) are solved for each geometrical con昀椀guration (step 7, Fig. 2).Subsequently, the impacts of the degrees of freedom on both 昀氀uid dynamic and thermal performance are analyzed, along with other DOFs (step 8).If new opportunities to facilitate the 昀氀ow system are identi昀椀ed and deemed feasible for implementation, modi昀椀cations to the system can be made based on the principles of the constructal law (step 9).In the present study, the constructal law was applied as a powerful tool to comprehend the in昀氀uence of the design on the arrangement of cylinders.However, the constructal law is a much broader principle and important discussion about the principle can be found in the literature [23][24][25].

Numerical procedures and veri昀椀cation/validation of the computational method
The Eqs. ( 1)-( 4) are numerically solved using the 昀椀nite volume method (FVM), implemented in the software FLUENT, version 14 [36][37][38].The numerical procedures include the utilization of the 2nd order upwind interpolation function for treating the advective terms, the semi-implicit method for pressure linked equations (SIMPLE) to handle For the spatial discretization, hybrid meshes (triangular and rectangular) are generated with the software GMSH (version 3.0.6)[39].A region with structured rectangular mesh, extending up to a diameter of 1.2 m, is implemented around the cylinders to facilitate the capture of velocity and temperature gradients in this area.Fig. 4 illustrates the hybrid mesh and the detail of the re昀椀ned mesh around the cylinders.To de昀椀ne the mesh recommendations, a grid independence study is performed for the case Re D = 150, Pr = 0.71, and S L1 /D = S L2 /D = S T /D = 1.5.The mesh is considered converged when the relative difference between a course and re昀椀ned meshes are lower than 0.5%, as expressed by: where Nu D j is the Nusselt number obtained with the course mesh and Nu D j+1 is the Nusselt number for the next re昀椀ned mesh.Table 1 shows the results of grid independence study performed in the present work, being recommended for continuity of the study a mesh with nearly 145,282 昀椀nite volumes.
Concerning the time-dependent analysis of the problem, for Re D ≤ 40, the structure of 昀氀uid 昀氀ow under forced convection is steady and symmetrical in the upper and lower downstream regions of one isolated cylinder.In this sense, the transient and the steady solutions led to similar 昀椀elds and results for the prediction of C D and Nu D .As Reynolds increases to Re D = 150, it is expected the arise of some unsteady and laminar structures such as von Karman vortices shedding and streets.However, this kind of 昀氀ow structure is meagered for situations where surrounding walls (as the symmetry walls in the present work) have some in昀氀uence on the 昀氀ow around one cylinder or bundle of cylinders.For these cases, it is required higher magnitudes of Reynolds number for the generation of more complex time-dependent structures (as those reported above) or an increase of the height of the channel to a magnitude as high as H = 75D, as recommended in the literature [40].For this work, the height of the computational domain is not augmented to obtain recommendations of the present arrangement for a thermal device equipment or bundle of tubes into channels subjected to any enclosure.Despite the velocity, pressure, and temperature 昀椀elds being more stable and, consequently, more similar to that achieved for the steady-state solution, the later solution led to differences up to around 39% in comparison with the results of literature for the cases with one or four cylinders reported in the literature [41,42] due to inadequate prediction of velocity and temperature 昀椀elds around the cylinder.Consequently, it was necessary to use the unsteady solution for adequate

Table 1
Grid independence study for the case Re D = 150, Pr = 0.71, and S L1 /D = S L2 /D = S T /D = 1.5.prediction of C D and Nu D for Re D = 150.To standardize the mathematical and numerical model throughout the investigation, the unsteady solution was also employed for the other cases of Re D = 10 and 40.

Nu
To de昀椀ne the time-step used in the simulations, a time-step sensibility study was performed for the same case investigated in Table 1 (Re D = 150, Pr = 0.71, S L1 /D = S L2 /D = S T /D = 1.5).Table 2 illustrates the Nusselt number for the arrangement for different time-steps, considering the same criteria used for the grid-independent study.Based on the results of Table 2, a constant time-step of Δt = 5.0 × 10 −3 s with 1000 iterations per time-step is employed in the present simulations.
To verify/validate the computational method, the spatial and timeaveraged Nu D for a single cylinder with Re D = 150 and Pr = 0.71 is compared with Hilpert's correlation [41].Furthermore, the arrangement of four cylinders, as detailed in Tables 1 and 2, is compared with Grimison's correlation [42], and the results are summarized in Table 3.Additionally, the local time-averaged Nu D as a function of the angle for a single cylinder with Re D = 20 and 200 using the present method is compared with the 昀椀ndings of Momose and Kimoto [43] and Bharti et al. [44], as illustrated in Fig. 5.In general, both spatial and time-averaged, and local Nusselt number exhibit good agreement with literature results, showing that the present method is suitable for the geometrical investigation conducted in this study.

Results and discussion
The initial investigation focuses on the impact of the ratio S T /D on the spatial and time-averaged C D (Fig. 6) and Nu D (Fig. 7) for various values of the S L2 /D ratio and the three distinct      Despite that, it is worth mentioning that one additional area of heat exchange is being included in comparison with a triangular arrangement of cylinders, for example, being an advantage in comparison with this previous studied array.Moreover, the effect of the ratio S L2 /D on Nu D,m changes with variation in the Reynolds number.For the lowest Reynolds, it is noticed a proximity between the frontal and intermediate cylinders, seeking to increase the momentum between the cylinders and the posterior cylinder more apart from the group of three cylinders.For the highest Reynolds, the 昀氀ow intensity increases, and the increase of longitudinal pitches becomes advantageous for the thermal performance.Fig. 9(d) shows that the in昀氀uence of the ratio S L2 /D over the (S T /D) o,T depends on the magnitude of the S L1 /D ratio.For S L1 /D = 1.5, an increase from S L2 /D = 1.5 to 2.0 decreases the resistance to 昀氀ow in the array gap, leading to an increase in S T /D.As the ratio S L2 /D increases further, the space among the cylinders increases, and the transversal pitch decreases again to augment the momentum in the region among the 昀椀rst three cylinders.For the other ratios of S L1 /D, once the 昀椀rst cylinder is distant from the intermediate cylinders, the departure of the posterior cylinder makes the distribution of the thermal boundary layer the main responsible for thermal performance.Then, the increase of S L2 /D conducts to an increase of (S T /D) o,T .
In the third level of geometrical investigation for 昀氀uid dynamic performance, the in昀氀uence of the S L1 /D ratio on the twice minimized drag coef昀椀cient, C D,2m , and the corresponding optimal con昀椀gurations, (S L2 /D) o,F and (S T /D) 2o,F , are investigated.Figs.10(a Finally, Figs. 12 and 13 illustrate the impact of Re D on the three times minimized drag coef昀椀cient, C D,3m , Nusselt number three times maximized, Nu D,3m , and the corresponding optimal con昀椀gurations for both performance indicators, (S L1 /D) o , (S L2 /D) 2o and (S T /D) 3o , respectively.Additionally, Fig. 12(b) depicts the pressure and velocity 昀椀elds obtained when the 昀氀ow reaches the steady state for the optimal 昀氀uid dynamics con昀椀gurations, while Fig. 13(b) illustrates the thermal 昀椀elds obtained at the steady state for the optimal thermal con昀椀gurations.
For 昀氀uid dynamic purposes, Fig. 12(a) shows a step decrease in C D,3m from Re D = 10 to Re D = 40, followed by a gradual decrease in the range 40 ≤ Re D ≤ 150.This behavior aligns with expectations, resembling that observed for a single cylinder.Regarding the optimal con昀椀gurations, asymmetrical patterns in the longitudinal direction are evident for all investigated Reynolds numbers for the present conditions, with (S L1 /D) o, F ∕ = (S L2 /D) 2o,F .Moreover, in all range of Re D examined, the upstream pitch (S L1 ) consistently exhibited a higher magnitude than the downstream pitch (S L2 ).This suggests that the use of variable pitches in the 昀氀ow direction can enhance the 昀氀uid dynamic performance of the system.It is also observed that the optimal ratio (S T /D) 3o,T remained constant and at the lowest magnitude irrespective of Re D .In the case Re D = 10, the optimal arrangement was reached for (S L1 /D) o,F = 3.0, (S L2 /D) 2o, F = 1.5, and (S T /D) 3o,F = 1.5, indicating the upstream cylinder is farther from the intermediate cylinders, while the three posterior cylinders are as close to each other as possible.As Re D increases to 40 and 150, the arrangement extends in the longitudinal direction, maintaining the asymmetrical pattern.Furthermore, the best con昀椀gurations are those that effectively distribute the pressure difference between the frontal and posterior regions of the arrangements, as seen in Fig. 12(b).These behaviors align with the principles of constructal law, as systems of higher magnitudes tend to be larger in size than 昀氀ow systems of lower magnitudes, and the design in 昀氀ow systems aims to distribute the imperfections in the most homogeneous manner possible.
For the thermal purposes, Fig. 13(a) illustrates a step increase in Nu D,3m in the range 10 ≤ Re D ≤ 40, followed by a smoother growth for Re D ≥ 40, as expected.Similar to the analysis of C D,3m , the behavior observed here also resembles the in昀氀uence of Re D on Nu D for a single cylinder.Regarding the optimal con昀椀gurations, the ratio (S T /D) 3o,T  this 昀氀ow condition (Re D = 150), or indicating a tendency to obtain symmetrical patterns for higher Reynolds numbers.Despite this, under the present conditions, asymmetrical con昀椀gurations conducted to the best 昀氀uid dynamic and thermal performances, and the use of varied longitudinal pitches along the 昀氀ow stream showed promise as a strategic approach.The results for the thermal 昀椀elds reinforce that 昀氀ow systems of low magnitude tend to generate smaller con昀椀gurations, and as the intensity of the 昀氀ow system increases, large con昀椀gurations emerge.Moreover, the generated designs followed the constructal law principles.Results also indicated that the use of a sort of strategy to modify the arrangement as a function of time and operation conditions can also bene昀椀t the system.
To compare the optimal con昀椀gurations for 昀氀uid-dynamic and thermal purposes to each other, the obtained values of both drag coef昀椀cient and Nusselt number are reported for each optimized arrangement of cylinders.In detail, with Re D = 10, 40, and 150, the optimal con昀椀gurations for 昀氀uid dynamic purposes lead respectively to Nu D = 0.8214, 1.7468, and 3.6737, which are 52.48%,45.90%, and 36.81%lower than the Nusselt values obtained by the corresponding optimal con昀椀gurations for thermal purposes (Nu D = 1.7286, 3.2291, and 5.8140, respectively).On the other hand, with Re D = 10, 40, and 150, the optimal con昀椀gurations for thermal purposes lead respectively to C D = 2.9295, 0.0943, and 0.0041, which are 144.19%,123.46%, and 105.00% lower than the drag coef昀椀cient values obtained by the corresponding optimal con昀椀gurations for 昀氀uid dynamic purposes (C D = 1.1997, 0.0422, and 0.0020, respectively).

Conclusions
The present numerical work explored the con昀椀guration of a staggered arrangement of four cylinders using a constructal design approach.The problem is constrained by the areas of the cylinders and occupation area of the arrangement, and it involves three degrees of freedom (S T /D, S L2 /D, and S L1 /D).This setup enabled an examination The results provided signi昀椀cant insights and recommendations for the design of the arrangement of cylinders concerning their 昀氀uid dynamic and thermal performances, measured by the drag coef昀椀cient (C D ) and Nusselt number (Nu D ).Overall, the lowest and highest magnitudes of transversal pitches (S T /D) lead to the best 昀氀uid dynamic and thermal performances.In particular, comparing the optimal and worst con昀椀gurations for 昀氀uid dynamic (thermal) purposes, differences up to 150% (54%) on the C D (Nu D ) value were obtained by varying the value of S T /D.Notably, the transversal pitch (S T /D) exhibited greater sensitivity in the performance compared to longitudinal pitches, irrespective of the Reynolds number (Re D ) under investigation.Indeed, differences up to 13.6% were obtained in the value of C D,m , and up to 1.2% in the value of Nu D,m , by varying the S L2 /D ratio.In addition, the variations in the value of S L1 / D led to differences of up to 35% in the value of C D,2m , and up to 5.2% in the value of Nu D,2m .
Despite their lower in昀氀uence compared to S T /D, longitudinal pitches made signi昀椀cant contributions to minimizing C D and maximizing Nu D .The results also indicated that employing different pitches (S L1 /D, and S L2 /D) could be an effective strategy for enhancing the 昀氀uid dynamic and thermal performance of the array.This observation was supported by the achievement of asymmetrical optimal con昀椀gurations for all investigated Reynolds numbers.Furthermore, the results showed a strong dependence of the arrangement con昀椀guration on Reynolds , where the highest S T /D led to the best thermal performance, the optimal thermal arrangement spanned over the entire occupation area.In general, the optimal con昀椀gurations were those that distributed imperfections in the most homogeneous manner possible, in accordance with the constructal law.
In future works, it is recommended to explore constraints related to other occupation areas, consider higher Reynolds numbers and other working 昀氀uids, and extend the study to convective 昀氀ows under different conditions, such as the simulation of mixed convective cases.These investigations can provide an even greater comprehension of the in昀氀uence of the design of cylinder' arrays under various scenarios.It is also recommended to develop a sort of transient control of the arrangement to improve its 昀氀uid dynamic and thermal performance.

Declaration of competing interest
The authors declare that they have no known competing 昀椀nancial interests or personal relationships that could have appeared to in昀氀uence the work reported in this paper.

Fig. 1 .
Fig. 1.Problem domain of the staggered arrangement of four cylinders subjected to forced convective 昀氀ows.

Fig. 2 .
Fig. 2. Diagram of the application of constructal design for the geometrical investigation of the arrangement of cylinders subjected to forced convective 昀氀ows.
Re D = 10, 40 and 150.Throughout all cases, a constant ratio S L1 /D = 1.5 is maintained.Spe-ci昀椀cally, Figs.6(a) and 6(b) depict a graph and an isosurface with C D as a function of S T /D for different S L2 /D values at Re D = 10.Similarly, Figs. 6 (c) and 6(d) exhibit similar 昀椀gures for Re D = 40, while Figs.6(e) and 6(f) present the corresponding 昀椀gures for Re D = 150.Figs.7(a) and 7(b) illustrate a graph and an isosurface with Nu D as a function of S T /D for different S L2 /D values at Re D = 10, and Figs.7(c) and 7(d) display analogous plots for Re D = 40.Lastly, Figs.7(e) and 7(f) illustrate similar 昀椀gures for Re D = 150.Fig. 6 consistently demonstrates that an increase in the S T /D ratio results in a proportional increase in C D , regardless of the S L2 /D ratio and Re D values.Therefore, the con昀椀gurations with the smallest transversal pitches, (S T /D) o,F = 1.5, conducted to the optimal 昀氀uid dynamic performance for the cylinder arrangement.When comparing the optimal and worst con昀椀gurations, differences of nearly 150%, 106%, and 63% are obtained, highlighting the sensitivity of the S T /D ratio on the C D .Additionally, some differences in C D magnitudes are observed for different S L2 /D, being more evident in Figs.6(a), 6(c), and 6(e) when S T / D ≥ 2.5.This behavior is corroborated in the isosurfaces presented in Figs.6(b) -(f) for Re D = 10, 40 and 150, respectively.Fig. 7 depicts a similar in昀氀uence of S T /D on Nu D , where Nu D increases with the elevation in S T /D.However, the focus here is on maximizing Nu D .In this context, the optimal con昀椀gurations are those with the intermediate cylinders positioned as far apart as possible, (S T /D) o,T = 5.0.The differences between the best and worst con昀椀gurations are approximately 54%, 51%, and 28%, underscoring the sensitivity of the S T /D ratio on Nu D .Results also reveal variations in Nu D magnitudes for different S L2 /D ratios.For Re D = 10, the differences are more pronounced at the lowest ratios of S T /D.As the Reynolds number increases, particularly for Re D = 150, differences become more apparent in the optimal region of S T /D, S T /D ≥ 2.5.It is also worth mentioning that, for Re D = 150, the optimal S T /D ratio was obtained for (S L2 /D) o,T = 3.0 and (S T /D) 2o,T = 3.5, probably by the balance between the interaction of thermal boundary layers and the momentum magnitude in the gaps around the cylinders.When compared with the triangular arrangement of cylinders studied in Barros et al. [30], the inclusion of an additional cylinder leads to a similar trend in the impact of S T /D on C D and Nu D for forced convective 昀氀ows, with few differences in the optimal con昀椀gurations of S T /D.In the second level of geometrical investigation, Fig. 8 illustrates the in昀氀uence of S L2 /D ratio on the once minimized C D , C D,m .The minimal C D values obtained in Figs.6(a), 6(c), and 6(e) are used to construct the curves for S L1 /D = 1.5 at Re D = 10, 40, and 150, respectively, as shown in Figs.8(a) -(c).The same process used to generate Fig. 6 for S L1 /D = 1.5 is repeated for the different ratios S L1 /D = 3.0 and 4.0, and the results for the optimal con昀椀gurations are also used here.Fig. 8(a) indicates that, for Re D = 10, an increase of S L2 /D leads to an augmentation of C D,m .For S L1 /D = 3.0, a step variation is observed in the range 1.5 ≤ S L2 /D ≤ 2.0, presenting a different behavior than noticed for the other investigated S L1 /D ratios.With an increase in Reynolds number to Re D = 40, Fig. 8(b), the behavior of C D,m as a function of S L2 /D changes, showing a decrease of C D,m with an increase of S L2 /D.At Re D = 150, Fig. 8(c), C D,m is insensitive to changes in S L2 /D for S L1 /D = 1.5 and 4.0, while for S L1 / D = 3.0, a step decrease in C D,m in the range 1.5 ≤ S L2 /D ≤ 2.0 followed by a constant magnitude of C D,m is observed.For all Reynolds number, a step variation of C D,m is perceived in this range, indicating that, in this combined region with S L1 /D = 3.0 and the range 1.5 ≤ S L2 /D ≤ 2.0, the posterior cylinder is under the in昀氀uence of the 昀氀uid dynamic boundary layers of the frontal and intermediate cylinders.It is important to Table 2 Time-step independent study for the case Re D = 150, Pr = 0.71, and S L1 /D = S L2 / D = S T /D = 1.5.

Fig. 5 .
Fig. 5. Nusselt number as function of the angle of the cylinder obtained with the present computational model and results of literature for distinct Reynolds numbers (Re D = 20 and Re D = 200).

Fig. 6 .
Fig. 6.In昀氀uence of the ratio S T /D on drag coef昀椀cient (C D ) for various S L2 /D ratios and distinct Reynolds numbers: a) Re D = 10, b) Re D = 10 (isosurface), c) Re D = 40, d) Re D = 40 (isosurface), e) Re D = 150, f) Re D = 150 (isosurface).(For interpretation of the references to colour in this 昀椀gure, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. In昀氀uence of the ratio S T /D on Nusselt number (Nu D ) for various S L2 /D ratios and distinct Reynolds numbers: a) Re D = 10, b) Re D = 10 (isosurface), c) Re D = 40, d) Re D = 40 (isosurface), e) Re D = 150, f) Re D = 150 (isosurface).(For interpretation of the references to colour in this 昀椀gure, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8. In昀氀uence of the ratio S L2 /D on the once minimized drag coef昀椀cient, C D , m , for various S L1 /D and distinct Re D : a) Re D = 10, b) Re D = 40, c) Re D = 150.

Fig. 9 .
Fig. 9. In昀氀uence of the ratio S L2 /D on the once maximized Nusselt number, Nu D,m , for various S L1 /D and distinct Re D : a) Re D = 10, b) Re D = 40, c) Re D = 150, and d) corresponding optimal con昀椀guration for (S T /D) o,F with Re D = 150.(For interpretation of the references to colour in this 昀椀gure, the reader is referred to the web version of this article.)

Fig. 10 .Fig. 11 .
Fig. 10.In昀氀uence of the ratio S L1 /D on the twice minimized drag coef昀椀cient, C D,2m , and the corresponding optimal con昀椀gurations for distinct Reynolds numbers: a) Re D = 10, b) Re D = 40, c) Re D = 150.(For interpretation of the references to colour in this 昀椀gure, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. In昀氀uence of the Reynolds number on: a) the three times minimized drag coef昀椀cient, (C D,3m ), and the corresponding optimal con昀椀gurations, b) optimal pressure and velocity 昀椀elds.

Fig. 13 .
Fig. 13.In昀氀uence of the Reynolds number on: a) the three times maximized Nusselt number, Nu D,3m , and the corresponding optimal con昀椀gurations, b) optimal temperature 昀椀elds.

Table 3
Comparison of results for Nu D obtained with the present method for one and four cylinders.