LARGE EDDY SIMULATION OF FLOW OVER ELLIPTIC CYLINDER ARRAY IN SQUARE CONFIGURATION AT SUBCRITICAL REYNOLDS NUMBERS

Flow over arrays of cylinders of different cross-sections has been an essential area of research interest for a long time. For the estimation of flow characteristics around arrays of cylinders, numerical simulation of the flow past four elliptic cylinders in a square arrangement for a range of spacing ratios is analysed in this paper. Three spacing ratios considered in the present study are 3.45, 4.14 and 5.17. The results from the numerical simulation are compared with the experimental findings obtained previously using square cylinder arrays. The fluid flow modelling is performed by applying three-dimensional LES (Large-eddy simulation) with commercial software ANSYS Fluent 19R1. The results from the simulation for elliptic cylinder arrays involve both quantitative and qualitative analyses in the form of various patterns of flow, drag and lift coefficients, St (Strouhal number) and PSD (Power Spectral Density) plots. The Strouhal number values found for cylinder 1 and cylinder 2 in case of elliptic cylinder arrays are 0.642 and 0.703. It is observed that the force coefficients encountered by the cylinders are moderately varying for different spacing ratios.


INTRODUCTION
Flows past arrays of elliptical and circular cross-section cylinders is a significant research area in ocean engineering. The elliptic cross-section is categorized as a fundamental cross-section in consideration of its possibility to acquire a flat plate and a circular cross-section by altering its axis ratio. Also, elliptic cylinders possess the advantage over circular cylinders by having small values of drag coefficient than the circular cylinders and a little wake region as well. Various experimental and numerical studies have been performed considering the different arrangement of cylinder arrays such as inline, staggered and tandem.
Modi and Dikshit [1] experimentally studied the aerodynamics of a single elliptic cylinder for a Reynolds number 68000 by varying the attack angle and axis ratio. Lugt and Haussling [2] found that flow around the elliptic cylinder is quite complicated. Still, some common flow characteristics were similar to circular cylinders such as the formation of different types of bubbles single as well as multiple. There was the formation of some kind of alleyway which depends upon the vortices and the bubble formation. In their study, they found steady eddies which vary with the angle considered in the study. Ota and Nishiyama [3] carried out experimental research on flow over two elliptic cylinders in a tandem arrangement at an axis ratio of 1:3. The attack angles taken for the study range from the values of 0 to 90°. The value of Reynolds number chosen for the study was in the subcritical regime. The angle of attack was varied, and it was found that for different attack angles, the value of Cp (pressure coefficient) was different. Ota [4] studied heat transfer between elliptic cylinders in tandem. The major and minor axis ratio, in this case, was 1:2, and the air was taken as the working fluid for the study. The Reynolds number used was 15000 and 80000 according to the major axis length.
Jackson [5] observed that vortex shedding starts from the third cylinder. When the flow moves towards downstream cylinders, an increment in the vortex shedding was found in the study. Mittal and Balachander [6] found the location of the separation layer and the start of the first vortex formation of fluid flow around a cylinder. Wong et al. [7] obtained the simulation for the stream flowing over two square cylinders at high Reynolds number. Nair and Sengupta [8] studied the flow over elliptic cylinder Journal of Thermal Engineering, Research Article, Vol. 7, No. 1, pp. 204-219, January, 2021 205 numerically by applying DNS (Direct Numerical Simulation) and obtained a comparison of the results for Reynolds number 3000.
Castiglia and Balabani [9] conducted both numerical and experimental studies on arrays of elliptic cylinders in the subcritical regime with an axis ratio of 1:2. They applied the Smagorinsky turbulence model for performing the three-dimensional simulation. In their work, it was observed that the shedding of vortices starts after the first cylinder. Franke and Franke [10] used the Reynolds number of 15000 for LES (Large-eddy simulation), and the spacing ratio used was 1.5. The bistable flow around the cylinder was observed in the investigation at this spacing ratio. Kim and Sengupta [11] noted that at spacing ratio 3.5, the patterns of flow are different and consist of two phases. The extent of occurrence of a vortex varies with the Reynolds number and the spacing ratio.
Agrawal et al. [12] performed simulation at a low value of Reynolds number for flow past two square cross-section cylinders arranged beside each other by applying the Lattice Boltzmann method. They observed in-phase and antiphase regimes. Liu and Cui [13] studied 3-D (three-dimensional) flow numerically, past two side-by-side arranged circular cross-section cylinders at Reynolds number 200. Ibrahim and Gomaa [14] investigated numerically and experimentally the flow characteristics with heat transfer past arrays of the elliptic cylinder. The range 5600 to 40000 for Reynolds number was used. They performed experiments by considering air as a working fluid and used the model for the turbulent flow as k-Ɛ RNG (Re-Normalisation Group) model. They found that the heat transfer coefficient per unit of pumping power is higher when the angle of attack is 0 º and lower when it is 90º.
Lam and Zou [15] numerically and experimentally investigated the flow over circular crosssection cylinder arrays for a spacing ratio, which varies from 1.5 to 5.0, and Reynolds number ranges from 11000 to 20000. Kumar and Vengadesan [16] used Large-eddy simulation to investigate the flow interface between two different size square cylinders which were arranged in a side-by-side arrangement. The value of 50,000 was used for Reynolds number. The vortex shedding was higher for the smaller cylinder than for the larger cylinder. Alawadhi [17] studied flow past elliptic cylinder arrays in an inline arrangement with an inclination angle. The inclination angle increases from 0° to 90°. No appearance of the vortex shedding was observed till 3rd cylinder, but there was instability in the flow. Vortex shedding was observed from the 4th cylinder onwards, the shear layer is separated from the 3rd cylinder roll up and impinges on the 4th cylinder. The value of the Nusselt number is more significant for the cylinder having a larger inclination angle. It was observed that with an increase in the inclination angle, there is a substantial increase in Cf (force coefficient). Peng et al. [18] investigated flow over two elliptic cylinders arranged beside each other. They analysed the effect of variations in the various parameters of the stream, such as the attack angle and spacing ratio. Liu et al. [19] studied experimentally flow past arrays of four square cylinders in a square arrangement at subcritical Reynolds number. They found the values of drag and lift coefficient for the various spacing ratios and angles of attack.
Mahmood et al. [20] presented a hybrid approach to solve the three-dimensional fluid flow over five staggered tubes in the direction of the flow. The values of Reynolds numbers chosen for the investigation was 100 and 300. The values of the Prandtl number and pitch-to-diameter ratio taken for the inquiry was 0.7 and 0.5, respectively. The five tubes were arranged in the equilateral triangle pattern. Sudarma et al. [21] studied the combustion process numerically inside a bluff body and compared the results for different turbulence models. Sheikholeslami et al. [22] simulated nanofluid flow through a pipe consisting of twisted tape turbulators. The exergy variations in the flow were investigated. Kariman et al. [23] carried out energy and exergy analysis to optimize the water desalination system, which can be very beneficial on an industrial scale. Kariman et al. [24] introduced accurate modelling of a new desalination system with brine tank and carried out energy and economic analysis for the system.
Hoseinzadeh et al. [25] investigated the heat transfer of laminar and turbulent pulsating nanofluid flow in a two-dimensional channel. Sheikholeslami et al. [26] simulated steady turbulent flow of copper oxide nanofluid by utilizing the k-ε model. They found that exergy loss is inversely proportional to the increase of pumping power. Sheikholeslami et al. [27] carried out modelling of turbulent heat transfer of homogeneous nanofluid due to inserting double twisted tapes. Sarafraz et al. [28] reported the results of the experiments conducted on the convective heat transfer of graphene nanoplatelets dispersed in water-ethylene glycol. Sarafraz et al. [29] did experiments for estimating thermal performance and heat transfer of n-pentane-acetone and n-pentane-methanol mixtures inside a gravityassisted thermosyphon heat pipe. The optimum heat pipe tilt angle was 55° for the binary mixtures, while it was 65° for the pure liquids. Laouira et al. [30] investigated heat transfer phenomena numerically inside a horizontal channel with an open trapezoidal enclosure subjected to a heat source of different lengths. It was observed that both the local and the average Nusselt numbers increase as the local heat source length increases. Gourari et al. [31] studied numerically natural convection flow of water between two coaxial cylinders. Results were determined for the laminar flow at different inclination angle. The best heat transfer was obtained for the inclination angle 90°. Salih et al. [32] developed an adaptive mesh refinement algorithm to simulate the fluid flow around a thin object and determine the computational cost. Baghban et al. [33] did an experimental and modelling study of heat transfer performance of water-SiO2 nanofluid in quadrangular cross-section channels. Hoseinzadeh et al. [34] performed exergeoeconomic, and optimization analyzes for a model for seawater desalination. Abadi et al. [35] conducted a numerical and experimental study on the energy efficiency of a regenerative heat and mass exchanger by utilizing the counter-flow Maisotsenko cycle. Turan [36] numerically analysed for laminar mixed convection in a square cross-section cylinder. Various flow contours were observed for the study.
It can be seen from the above discussion that minimal numbers of studies have been performed on arrays of elliptic cylinders. For investigating the variation in the force coefficients, a square crosssection cylinder studied previously has been changed to elliptic cross-section keeping the same hydraulic diameter. The present study aims to fill this gap to some extent by numerically simulating and characterizing the flow over elliptic cylinder arrays having a square flow configuration using large-eddy simulation. For validation purpose, the experimental work by Liu et al. [19] is simulated numerically by applying LES to square cylinder arrays in a square configuration with commercial software ANSYS Fluent 19R1.

METHODOLOGY
The numerical simulation can be performed by LES, RANS (Reynolds Average Navier Stokes) and DNS (Direct Numerical Simulation). In the present work, LES is utilized to perform the numerical simulation. LES, a space filtering technique in CFD (Computational Fluid Dynamics), directly computes the large-scale eddies while modelling the smaller-scale swirls. The computational cost taken by LES is less as compared to DNS.
Moreover, LES is more accurate and give better results than RANS because large eddy contains most of the energies and LES directly solve them. LES is a numerical technique to solve turbulent flow problems. A filter function is used in LES to differentiate between the large scales and small scales. The characteristic filter cut-off width usually denoted by Δ in LES. All eddies larger than Δ are directly solved, while smaller eddies are modelled. In 3-D computations with grid cells having a length, width and height as Δx, Δy and Δz respectively, the cut-off width is given by Equation (1). (1) The filtered momentum equations are given in Equation (2).
Here, denotes the rate of change and convective fluxes of filtered momentum.
, this term indicates the gradient and diffusive fluxes of filtered energy, respectively.
, means the divergence term caused by the filtering operation. 207

Computational Domain
The computational domain for the numerical simulation of the flow around four elliptic cylinder arrays in a square configuration with the computational domain dimensions is shown in Fig. 1. The hydraulic diameter of the elliptic cylinder is taken as 0.029 m, which is same as the hydraulic diameter of the square cylinder considered in experimental work by Liu et al. [19]. The major axis of the elliptic cylinder is taken as 0.0222 m and the length of the minor axis, b is taken as 0.0111 m (i.e. the axis ratio is equal to 1:2). The spacing ratio in square configuration elliptic cylinder arrays is L/De = 3.45 for case 1, 4.14 for case 2 and 5.17, for example, 3. The spacing ratios were selected equal to the values considered in the previous findings of square-cylinder arrays in the square configuration by Liu et al. [19]. Air at temperature 298.15 K is used as the working fluid. The computational domain for case 1, i.e. for spacing ratio 3.45 is shown in Fig. 1. The freestream velocity for all the three cases is selected as 23.7 m/s corresponding to the value of Re = 45800.

Meshing
Meshing is one of the essential steps in obtaining an accurate boundary layer solution. In the present work, meshing is done in ICEM CFD (Integrated Computer Engineering and Manufacturing), which provides better control over grid generation and creates structured mesh. The structured mesh is designed because it has fewer convergence issues and gives better results for less complicated geometries. Figure 2 shows the quality of the mesh used for the elliptic cylinder arrays.
To resolve the separation of the boundary layer and the wake region, mesh cells used in this research were refined near the cylinder wall. The value of y + = 1 is used for the first cell height, which is placed away from the cylinder walls. Ten per cent stretching is given after the first cell height. In the numerical simulation, 2531160 nodes and 2384880 cells are used for the computational domain.

Boundary Conditions
Suitable boundary conditions are applied to the boundaries of the computational flow domain for achieving the actual conditions of the flow. Outflow boundary conditions and velocity inlet are used at the outlet and inlet in the flow domain. For the sidewalls, the symmetry boundary condition is applied. Noslip boundary condition is applied on the cylinder walls, and the periodic boundary condition is used in the span-wise direction.

Solution Procedure
The processing operations available in simulation settings are serial processing and parallel processing. In this study, parallel processing is applied to increase the computations speed of the simulation. The simulation type used is three-dimensional, double-precision, unsteady, implicit and pressure based. Large-eddy simulation with dynamic Smagorinsky-Lilly model is used to perform threedimensional simulations. In this study, the SIMPLE segregated algorithm is used. Least squares cellbased technique is utilised for gradient evaluation. This technique offers more flexibility about the order of accuracy achieved. Pressure interpolation is done with the help of a second-order interpolation scheme. For obtaining the precision of the transient formulation bounded second-order implicit scheme is utilized. ANSYS Fluent provides two methods for solution initialization, i.e. Standard initialization and Hybrid initialization. Among these, Standard initialization method is used for solution initialization in this study which is simple and takes less initialization time. The number of iterations used per time step is fixed to 20 iterations so that errors can be reduced. The solution is initialised with the value of the velocity at the inlet for the velocity inlet boundary condition. The results obtained from the three-dimensional Largeeddy simulations are presented in the next section.

Model Validation
This section includes the results of LES applied to square cylinder arrays and elliptic cylinder arrays. To validate our numerical simulation, model validation is done for spacing ratio (L/D) = 3.45 and Reynolds number = 45800. The results obtained by applying LES for square cylinder arrays are compared with the experimental results of Liu et al. [19]. The vortex shedding frequency, drag and lift coefficient are compared with the experimental findings. Figure 3 shows the comparison of the PSD (Power spectral density) plot for the present values obtained by LES and the experimental data of Liu et al. [19]. Table 1 shows the values of the drag and the lift coefficients for the four cylinders which are arranged in the square configuration in the clockwise direction. From Table 1, it can be observed that the values of the numerical simulation investigated are in good agreement with the mean drag and lift coefficient values of the published experimental data. From the PSD plot, it is seen that the vortex shedding frequency is approximate. 94 Hz obtained by Liu et al. [19] and 104 Hz, is captured in the present work. Table 1 shows the comparison of average drag and lifts coefficient values and the Strouhal numbers obtained by LES in the current work with the published experimental data collected by Liu et al. [19]. From the PSD plot in Fig. 3 and drag, lift and Strouhal number values in Table 1, it can be observed that the results of three dimensional LES match well with the experimental results.

RESULTS AND DISCUSSION
Results comprise of quantitative data obtained from the large-eddy simulation of flow past arrays of the elliptic cylinder for different spacing ratios. The plot between PSD and Vortex shedding frequency for spacing ratio 3.45 for cylinder 1 and 2 are shown in Fig. 4. (a) and (b) respectively. Frequency plot is obtained from FFT ( Fast Fourier Transform ) analysis of time sequence data. The frequency of vortex shedding is lower for the cylinders placed upstream as compared to downstream cylinders. The increase in vortex shedding frequency is because of the early separation of the boundary layer due to the impingement of the vortices coming from the upstream cylinders.  Values of the mean drag coefficient and lift coefficient for different spacing ratios (L/De) between the elliptic cylinders that vary from 3.45 to 5.17 are shown in Table 3. From the data, we can see that the drag coefficient values on the upstream side of the cylinder are more significant than the downstream of the cylinder because the flow initially impinges on the upstream cylinder that reduces the drag effect on the cylinder that is placed in the wake region of the upstream cylinders. Still, the mean lift coefficient of the downstream cylinder is higher than the upstream cylinder. It is due to the uninterrupted formation of periodic vortices, and cylinders gain more fluctuating forces. From Table 3, it can be observed that the drag experienced by the downstream elliptic cylinder's arrays is less than the drag on the upstream cylinders. Table 3 is showing average drag and lift coefficients values at different spacing ratios for elliptic cylinder arrays. The Variation of the mean drag coefficient for four elliptic cylinders at different spacing ratios is shown in Fig. 5 by plotting the values obtained from the LES (Large-eddy simulation) for elliptic cylinder arrays. At a spacing ratio of 3.45, the value of the mean drag coefficient for cylinder 1 and cylinder 4 is more as compared to the cylinders 2 and 3 respectively. The initial impact of the flowing fluid on the upstream cylinders increases the drag forces on the cylinders. It can be observed from the plot in Fig. 5 that as the spacing ratio increases from 3.45 to 4.14, the drag force increases, but when the spacing ratio is further increased to 5.17, then the drag on the cylinder starts decreasing. The Variation of the mean lift coefficient for four elliptic cylinders at different spacing ratios is shown in Fig. 6 by plotting the values obtained from the LES (Large-eddy simulation) for elliptic cylinder arrays. It can be observed from the plot in Fig. 6 that as the spacing ratio increases from 3.45 to 4.14 the lift force on the cylinder decreases but when the spacing ratio is increased further from 4.14 to 5.17 then the lift force on the cylinder starts rising due to the enhancement of the shedding of the vortices in the flow field. The Variation of the drag coefficient with flow time for four elliptic cylinders at a spacing ratio of 3.45 is shown in Fig. 7 by plotting the values obtained from the LES (Large-eddy simulation) for elliptic cylinder arrays. It can be observed from Fig. 7 that initially when the numerical simulation starts, the drag coefficient value is fluctuating more but as the numerical simulation proceeds further the values of the drag coefficients start converging with the flow time.  The drag coefficient variation with flow time for four elliptic cylinders at 4.14 spacing ratio is shown in Fig. 9 by plotting the values obtained from the LES (Large-eddy simulation) for elliptic cylinder arrays. Lift coefficient variation with flow time for four elliptic cylinders at a spacing ratio of 4.14 is shown in Fig. 10 by plotting the values obtained from the LES (Large-eddy simulation) for elliptic cylinder arrays. The drag coefficient variation with flow time for four elliptic cylinders at 5.17 spacing ratio is shown in Fig. 11 by plotting the values obtained from the LES (Large-eddy simulation) for elliptic cylinder arrays.  Eddy viscosity affects the velocity profile in a turbulent flow. It has been observed in our simulation that as the spacing ratio between the cylinder increases, the intensity of eddy viscosity increases. Figure 14 shows the eddy viscosity for a spacing ratio of 3.45. From the eddy viscosity contours, it can be observed that the wake regions are different for the upstream and the downstream cylinders in the flow direction. The reason for this is the change in the velocity caused after striking of the incoming fluid on the elliptic cylinder walls. For different values of the Reynolds number, the contours of the turbulent viscosity vary. This variation is according to the change in the velocity in the flow direction.

Figure 14. Eddy viscosity contours
As the flow with a turbulence intensity of 0.7% comes in contact with the elliptic cylinder, the velocity becomes zero at the stagnation point. Position of stagnation point is symmetric on both the cylinders that are on the upstream side, but asymmetrically located stagnation points are obtained on the downstream cylinders. Figure 15 shows the stream-wise velocity component contours for a spacing ratio of 3.45.  Figure 16 shows the contours of the velocity component in the transverse direction. In the wake region, we are getting positive and negative velocity, as shown in Fig. 16, due to the formation of clockwise and anti-clockwise vortices. As the flow proceeds, in the downstream of the flow, the strength of the vortices decreases.

CONCLUSION
Based on the investigation on an array of four elliptic cylinders in the square configuration using commercial CFD software for different spacing ratios and subcritical Reynolds number, the following conclusions have been drawn:  It is observed that the Large-eddy simulation model using dynamic Smagorinsky-Lilly turbulence model contributes to more appropriate predictions of the force coefficients (drag and lift) for each chosen spacing ratio and flow velocity.  The values obtained indicate that this configuration of elliptic cylinder arrays experiences lower drag forces than the square cylinder arrays.  At higher spacing ratio, the drag coefficient decreases as compared to the drag coefficient at the minimum spacing ratio. This reduction is due to the drop in the strength of the vortices.  The lift coefficient increases at large spacing ratio due to more fluctuations in the flow.  Velocity contours are observed in the streamwise as well as transverse direction. It is seen from the contours that free shear layers separated from the upstream cylinders roll up into mature vortices in between the cylinders and impinges on the downstream cylinders.  The PSD plots are plotted using FFT (Fast Fourier Transform) for the model validation and elliptic cylinder arrays. The vortex shedding frequency obtained was 104 Hz for model validation and 525 Hz for the elliptic cylinder arrays investigated.  It is observed that as the spacing ratio between the cylinders increases, the intensity of eddy viscosity increases.  The value of the mean drag and lift coefficients at different spacing ratios were observed to be moderately varying at the subcritical Reynolds number used for the numerical investigation. The obtained results for the model validation are observed to be in synchronism with the experimental published results. Moreover, the comparison of the simulated results with elliptic cylinder arrays using Large-eddy simulation proves that the bluff body, which is streamlined to the flow gives smoother flow with smaller drag forces. In the present work, Large-eddy simulation of flow past arrays of four cylinders in a square configuration is performed. However, there can be various other configurations possible, such as a hexagonal or an octagonal arrangement. Cylinders of the different cross-section can be arranged in arrays in such type of arrangements. Numerical simulation can be performed for different spacing ratios and at different Reynolds numbers. According to the desired accuracy of the results and availability of the computing resources, simulation models can be utilized. It can be regarded as a significant future scope of research to understand the flow characteristics better in case of arrays of cylinders. Greek symbols ρ air density, kg/m 3 ν kinematic viscosity of air, m 2 /s