Numerical Analysis of Flow Around a Cylinder in Critical and Subcritical Regime

: Modeling the wind flow around cylindrical buildings is one of the problems within urban physics. Despite the simple geometry of the cylinder, it is an interesting physical phenomenon. Par ‐ tial knowledge of flow field properties can be found in the literature, but in terms of their use for practical tasks, the data are still incomplete. The authors performed a numerical analysis of the flow around the smooth cylinder in the subcritical and critical regime for Reynolds numbers in the range of Re = 2.3 × 10 3 to 4 × 10 5 . Turbulent flow was solved using LES model and the numerical solution was compared with available data from experiments or standard. Analysis of the mean stream ve ‐ locity showed the elongation of the core of the wake with decreasing Re. The pressure coefficient evaluation showed a big difference between its distribution in the subcritical and critical regime. In the subcritical regime, a significant increase in the minimum value and a shift of the extreme close to the axis of the cylinder is proven. The results of the drag coefficient confirm a significant decrease in the transition from subcritical to critical regime, which is indicated in the cited experiments.


Introduction
Civil engineering and architectural engineering are strongly associated with the rapidly developing applied scientific discipline of urban physics, which offers a wide range of areas of interest [1]. Researchers investigate the complex relationship between spatial composition and building typology on the one hand and thermal and climatic conditions within and between buildings on the other hand in [2].
A holistic view of the system in the building industry, in terms of sustainability, takes into consideration environmental, economic, cultural, and social issues. The use of computational fluid dynamics (CFD) can significantly help to design a structure in specific situations, evaluate the requirements for load-bearing capacity and reliability, and verify its properties. The common goal of both sectors, urban and civil engineering, is to design ecological and energy efficient buildings, as well as to adapt to the natural and cultural environment.
Flow around a cylindrical object is one of the frequently solved problems; the simple geometry of the cylinder is an interesting physical phenomenon. It concerns various types of structures such as cooling towers, chimneys [3], buildings of circular shape (Figure 1), pylons for cable cars, bridge structures [4][5][6], offshore structures, air-cooled heat exchanges [7], storage tanks, and other industrial buildings [8] and their structural components. More detailed knowledge of the flow field and the effects of wind flow on objects finds its application not only in civil engineering [9], but also in wind engineering [10]. Problems of wind flow can relate to the layout of buildings and the shapes of the buildings themselves [11,12], as well as the type of cladding [3], the shape of balconies [13], and the like. , (b) design of a complex of residential buildings [13].
The nature of the wind flow is defined by the Reynolds number (1), which is a dimensionless parameter representing the ratio of inertia force to viscous force in a flow [15].
where u [m•s −1 ] represents the flow velocity,  m 2 s −1  is the kinematic viscosity of the running air, and D [m] is the diameter of the cylinder.
The Reynolds number affects the value of the drag coefficient cd. It defines the degree of drag force that acts in a direction that is opposite the relative flow velocity (horizontal direction). This is an important quantity in the dimensioning of structures of circular cross-section, and it defines the degree of loading of the structure due to wind. The course of cd depending on Re obtained from experimental measurements [16] is shown in Figure  2. It is clear from the graph that there is a significant decrease in cd in the region of the transition between subcritical and critical regimes and that the research is not fully clarified in this zone. There is a change in the boundary layer of the cylinder which significantly increases the complexity of numerical simulations. The lack of experimental data is associated with the frequent problem of achieving high Re numbers in the wind tunnel. Although partial results of experiments for some of the Re values for subcritical and critical regions can be found in the scientific literature, they do not give a sufficiently detailed picture of the cd in this region.
Another fundamental quantity used in civil engineering is the pressure coefficient cp, which defines the distribution of the pressure load on the cylinder [18][19][20]. The boundary layer and the structure of the near wake behind a cylinder have always attracted attention for theoretical reasons and consequently for practical applications [7,[21][22][23].
The objective of this work is to examine the suitability of using LES turbulent model for the effective calculation of turbulent characteristics of the flow around a cylindrical object and to contribute to the addition of information about the flow field properties. There are investigated drag coefficient cd, pressure coefficient cp [18][19][20], and velocity profile in the wake behind the cylinder [7,[21][22][23], these quantities are essential for wind engineering. Specifically, the flow around the cylinder in this research is done for Re = 2.3 × 10 3 to 4 × 10 5 , which is an interval that includes the subcritical and critical region [16]. Numerical simulations are carried out in ANSYS Fluent software by using high performance computers of the National Supercomputer Center IT4Innovations and are verified on the basis of available experimental data [24][25][26][27][28][29]. The ability to obtain relevant and reasonably reliable information regarding the course of cd in the critical area by the method of numerical modeling would be of great importance for further specification of wind load on structures.

Task Description
Isothermal flow around a smooth cylinder with the nature of the flow for eight different Re numbers, Re 2.3 × 10 3  4 × 10 5 ), is modeled according to Table 1. In all cases the cylinder diameter D = 0.1 m is identical, the Reynolds number change is ensured by the change in velocity of the flow. The basic parameters for the calculation are given in Table  2.

Numerical Model and Boundary Conditions
The task is solved in the academic version of ANSYS Fluent software (version 2020 R2) using the numerical Large eddy simulation turbulence model (LES model). The large eddy method is a simulation technique based on filtering a flow field into macro-and micro-scale eddy structures. Large-scale eddy structures are simulated directly. Turbulent structures of microscales, which are generally isotropic, are expressed using subgrid-scale models. In addition, these small vortices contribute little to the momentum transfer (and to heat transfer in anisotropic tasks), therefore they are expressed by parameterization schemes embedded in the equations for large vortices.
Filtration of the continuity equation (2) and the momentum equation (3) yields initial relations for the mathematical description of the present isotropic process by LES method.
Continuity equation: momentum equation (Navier-Stokes): where  [kg•m −3 ] is the density of the flowing medium, see Table 2, represents the time average of the velocity components, is stress tensor, ̂ is the time average of the static pressure, Pa is the tensor of the residual subgrid stress, which arises due to the filtration of large vortices. The left sides in (2) to (3) describe macrovortex structures (variables denoted by the canopy), the subgrid term for microstructures is expressed on the right side of the Navier-Stokes equation (3) and is expressed using subgrid-scale models in which subgrid turbulent viscosity is defined. Individual subgrid-scale models differ from each other in the description of turbulent viscosity. The presented problem is solved by wall-modeled large eddy simulation (WMLES) subgrid model [30], in which subgrid turbulent eddy viscosity is calculated with the use of a hybrid length scale where is the wall distance, is the strain rate, 0.4187, and 0.2 are constants, and is the normal to the wall inner scaling. The LES model is based on a modified grid scale to account for the grid anisotropies in wall-modeled flows: where ℎ is the maximum edge length for a cell, ℎ is the wall-normal grid spacing, and 0.15 is a constant.
The boundary conditions for all simulations are given in Table 3. symmetry: there is no friction on the side walls, the sides of the computational area do not affect the longitudinal velocity, the normal velocity and the flow of quantities across the border are zero upper and lower horizontal walls: wall: corresponds to wind tunnel conditions

Meshing
Unstructured tetrahedral mesh ( Figure 3) with hexahedral prismatic cells ( Figure 4) covering the boundary layer was used for numerical simulations. External dimensions of the mesh are 7 × 1.8 ×0.2 m (x × y × z), which corresponds to multiples of the cylinder diameter 70 D × 18 D × 2 D ( Figure 5). The dimensions of the computational domain were chosen to approximately match to the boundary conditions in the compared experiments. Because the experimental data come from different authors, different countries, and different periods, the computational area is not a model of a specific wind tunnel but tries to respect the general principles of the modeling of fluids. These are mainly: giving the fluid enough space to be affected as little as possible by the shear forces at the wind tunnel walls, and to comply with the recommended maximum of blockage ratio (ratio of the windward area of the object to the cross-section of the tunnel) of approximately 5%.
In total, two variants of the mesh with different number and average size of cells were used for the reason of the rather big range between minimal and maximal Re and because of the higher requirements of the LES model on the resolution of the mesh for high velocities.
The mesh for the lower velocities (u0 = 0.35-15 m•s −1 ) consisted of 1,138,960 cells. It has height of the first cell at the wall 1 × 10 −5 m, the hexahedral prismatic layer has 20 layers with the total thickness 0.0035 m. Prismatic layer is followed by tetrahedral cells that increasingly grow outward from the wall of the cylinder from the size 0.0015 to 0.03 m. The mesh for the higher velocities (u0 = 21-60 m•s −1 ) consisted of 9,111,680 cells. These parameters are different compared to the previously described mesh: height of the first cell at the wall 3.5 × 10 −6 m, 40 layers of the prismatic layer, the size of tetrahedral cells grows from the size 0.00075 to 0.015 m.
Both grids meet the condition for near wall modeling y + ≤ 1. This is a dimensionless quantity, dependent on the type of flow and other parameters, given by the formula: where  is density of the flowing medium [kg•m −3 ], yp is distance of the first cell point from It follows from (6) that, in keeping y + ≤ 1, there is an inverse relation between the first cell size at the wall and the velocity of the flow, which is why the mesh for all Re is not identical. To maintain the condition y + ≤ 1, it is necessary to create a very small first cell near the wall for higher velocities.

Flow Field Characteristics
Based on the available experimental data, it can generally be assumed that for the flow around the cylinder up to Re = 1.4 × 10 5 , the distance of the minimum mean flow velocity ux in the wake increases from the axis of the cylinder with decreasing Re number. The length of the so-called core of the wake also increases with decreasing Re (core of the wake means an area of swirling wake with a negative value of ux). This is obvious from the illustrative pictures of numerical simulations in Figure 6. Different velocity field distributions for two different flow regimes: subcritical, Re = 4 × 10 3 and critical, Re = 1.4 × 10 5 are shown in Figure 6a,b. The root mean square (RMS) of mean stream velocity fluctuation is shown in Figure 6c,d for the same flow regimes. Based on the streamwise velocity in the wake behind a cylinder axis, the velocity field characteristics are evaluated in this paper. The streamwise velocity is defined as a dimensionless quantity, so-called normalized mean stream velocity, and is expressed by the ratio ux/u0. The distance from the cylinder axis is defined by the ratio x/D, where x/D = 0 applies to the cylinder axis. The normalized mean stream velocity profiles in the wake at the level of the cylinder axis obtained for similar Re from physical measurements, as well as from numerical simulations in subcritical and critical regimes, are recorded in Figure 7. Figure 7. Normalized mean stream velocity for selected Re; CFD and experimental data (Frohlich [22], Beudan [21], Khashehchi [7], Chao Fu [23], Cantwell [18]), subcritical and critical regimes.

Normalized Mean Stream Velocity in the Subcritical Region
Experimental data for smooth cylinders in the subcritical region are taken from [21] and [22] for Re = 3.9 × 10 3 and from the newer [7] for Re = 4 × 10 3 or [23] for Re = 5 × 10 3 . Figure 7 shows CFD data for Re numbers close to the experiments. The minimum value of the normalized mean stream velocity in all cases is within the range ux/u0(−0,2; −0.3) and its distance of the cylinder axis is x/D = 2. The length of the core of the wake (the transition from negative to positive values ux/u0) is in these cases about x/D = 2.5. Only the source [22] presents a different course of ux/u0. However, the area in the immediate vicinity of the cylinder remains unclear (approximately x/D to 1.2). According to [23], the stream velocity is also negative in this part, while all numerical simulations show positive ux values. In [7], the results for Re = 4 × 10 3 are presented up to the distance x/D = 1.5 and the values in the immediate vicinity of the cylinder are not presented here. However, the profile data ux/u0 for the lower Re numbers, which are given in [7], prove the difficulty of the objective description of the velocity field and indicate the possibility of a short range with positive values ux in the immediate vicinity behind the cylinder.
Another situation is from the point of view of determining ux at a greater distance from the axis of the cylinder and determining its maximum value. In this case, both numerical simulations coincide (Re = 2. 3 × 10 3 and Re = 4 × 10 3 ). At the distance x/D = 3.5, where the data from [23] end, the normalized mean stream velocity values based on CFD calculations are in the range 0.57-0.7, while [23] presents ux/u0 = 0.46 and [7] even ux/u0 = 0.40. It can be said that all numerical simulations for Re = 4 × 10 3 show for the distance x/D ≥ 3.5 a better agreement with [22], when the normalized mean stream velocity is about the value ux/u0 = 0.7. Figure 7 have shown a shortening of the core of the wake with increasing Re number. The minimum value of normalized mean stream velocity decreases to the value ux/u0 = −0.31 and their distance from the axis of the cylinder by assumption is shortened by increasing Re up to the distance x/D ≈ 1. However, for flows with Re = 1.4 × 10 5 the experimental data [18,22] differ from CFD calculations (Figure 7).

CFD simulations in
From the point of view of defining the maximum value of ux at a greater distance from the cylinder axis, the numerical simulations for the flow with Re = 1.4 × 10 5 are compared with experimental data at a distance approximately x/D = 4.5 at the value ux/u0 = 0.75. For flows with the highest presented Re in CFD, the mean stream velocity stabilizes at a distance x/D ≈ 3 with the value ux/u0 = 0.81.

Pressure Load on the Cylinder Circumference-Pressure Coefficient cp
The pressure load on the cylinder circumference for all presented simulations was evaluated using a dimensionless cp coefficient. It is given by the ratio of the static pressure pi and the dynamic pressure pdyn relative to the reference point. The reference point was set 0.2 m behind the entrance to the area, which is 0.8D in front of the axis of the cylinder, when the flow field is not yet affected by the flowing obstacle.
Pressure coefficient for a constant air density at isothermal flow  is defined where pref is the static pressure at the reference point [Pa], pci is the resulting static pressure on the cylinder surface at the i-point [Pa], uref is the mean streamwise velocity at the reference point [m•s −1 ]. All data for the verification of numerical simulations obtained from experimental research are shown in Figure 8. They are taken mainly from [19], where the author focuses on the position of the minimum cp (maximum compressive load) for Re(1.3 × 10 2 ; 2.1 × 10 5 ). Results of these experiments fall into both the subcritical and critical regimes. Further results from the experiments used for comparison in this article fall into the critical regime and they are taken from [29] for Re = 1.5 × 10 5 and from [27] for Re = 1.5 × 10 5 .  (Norberg [19], Tani [29], James [27]).
One of the other goals of this work is to compare the pressure coefficient distribution cp for different Re with available experimental results. Due to the complexity of describing the problem for higher Re and due to the limitation of the number of available relevant experimental results, the resulting analyzes of the pressure coefficient cp in this article are divided separately for the subcritical and critical regime.

Pressure Coefficient cp-Subcritical Region
It can be seen in Figure 8 and Figure 9 that in the subcritical region, the value of the minimum of the pressure coefficient is in the range of cp,min ≈ (−1.15 to −1.25) both in the results of the experiments and in the results of the performed simulations. Very slight differences depending on Re are seen for cp distribution for Re > 8 × 10 3 . However, the difference between the minimum value of cp for Re = 3 × 10 3 and Re = 8 × 10 3 is evident. The CFD results also correspond to these results, there is a different minimum of cp around Re = 3 × 10 3 and Re = 2 × 10 4 .

Pressure Coefficient cp-Critical
Region, Re ≥ 1 × 10 5 As mentioned above, from the available experimental data, the description of the flow field in the critical region is not yet fully understood. The complexity of the situation is proved by the results of the mean pressure distribution results obtained from the experimental research in Figures 8 and 10, which differ significantly from each other, and the CFD results also differ from them. There is also an interesting comparison of cp with standard results [31]. The minimum position for Re = 5 × 10 5 according to [31] and CFD is   (Tani [29], Norberg [19], James [27]) and standard [31].

Drag Coefficient cd
In the present work the drag coefficient cd was examined. The results of the numerical simulations are compared mostly with the values from the literature [16,17] in the graph in Figure 2. The graph shows both a significant decrease in cd during the transition from subcritical to critical region and a partially unexplained value in part of the interval. When examining the value of the drag coefficient, the complexity of the issue in the transition from the subcritical to the critical region has been again demonstrated.
It is clear from Figure 2 and Figure 11 that the value of cd in this flow regime is not fully clarified. In Table 4 and Figure 11 there are synoptic values of the cited experimental results, as well as the numerical simulations solved in this study.  Figure 11. Relationship of drag coefficient cd and Reynolds number. Data from [16,[31][32][33], and Table 4.

Discussion
The paper is focused on the numerical simulation of the flow around the cylinder in the range Re2.3 × 10 3  4 × 10 5 ) using the LES model. It focuses on the selected flow characteristics that are necessary for sustainable architecture.
Velocity profile in the wake behind the cylinder confirms the elongation of the core of the wake with decreasing Re, which is a finding that finds application mainly in research on the issue of heat transfer conditions. The pressure coefficient evaluation confirmed a big difference between its distribution in the subcritical and critical region. According to Figure 8, this is a clear complexity of the experimental analysis for high Re. CFD results differ significantly in the critical area compared to (mutually different-inconsistent) experimental data, but they approach standard calculations. The distribution of the compressive load on the circumference of the cylinder showed in the critical regime a significant increase in the minimum values and a shift of the extreme close to the axis of the cylinder (min ≈ 85°). This approached the values given in the standard [31].
The drag coefficient is known to be sensitive to the transition between regimes and from the obvious differences between the curves in Figure 11, it is clear that the turbulent flow in the region of the transition from the subcritical to the critical region is a phenomenon that is not easy to analyze experimentally, which also applies to simplified mathematical models or incorrectly set LES model parameters. In the presented work, it was managed to capture the decrease in cd, which corresponds to the trend of available experiments [16,32,33] and the curves given in the standard [31].
The results of the presented numerical simulations with the selected calculation domain, computational mesh, and setting of calculation parameters showed that in the subcritical and critical modes it is possible to capture the trend of the observed flow properties.
Although various recommendations for performing CFD calculations of bluff bodies can be found in the literature, a comprehensive methodology of the procedure that would lead to reliable results in civil engineering is still not available. The presented work complements the existing research and contributes at least a little to the progress of the methodology of mathematical modeling in the critical regime. A valuable finding from this study is that the parameters used in the presented calculation lead to results close to the experimental and standard values. Unfortunately, they cannot yet be considered as a generally optimal calculation setting for this type of tasks, which could then be reliably used in practical calculations of wind load of buildings of more complex shapes. In further research, the authors will focus on performing calculations with other subgrid LES models Re [-] and try to perform a more demanding numerical simulation on a domain with significantly greater thickness to better capture the spatial character of vortex structures, which may be one of the factors influencing the results.