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

Abstract

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. Partial 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³ to 4 × 10⁵. 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 velocity 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.

1. 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.

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].

$$Re= \frac{u·D}{\upsilon },$$

(1)

where u [m·s⁻¹] represents the flow velocity, ν [m²⋅s⁻¹] 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 c_d. 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 c_d 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 c_d 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 c_d in this region.

Another fundamental quantity used in civil engineering is the pressure coefficient c_p, 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 c_d, pressure coefficient c_p [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³ to 4 × 10⁵, 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 c_d in the critical area by the method of numerical modeling would be of great importance for further specification of wind load on structures.

2. Methods

2.1. Task Description

Isothermal flow around a smooth cylinder with the nature of the flow for eight different Re numbers, Re ∈ (2.3 × 10³; 4 × 10⁵), is modeled according to

Table 1. Selected Re and corresponding velocities.

Subcritical	Re = 2.3 × 103	Re = 4 × 103	Re = 2 × 104
	u0 = 0.35 m·s−1	u0 = 0.6 m·s−1	u0 = 3 m·s−1
Critical	Re = 1 × 105	Re = 1.4 × 105	Re = 2 × 105	Re = 3 × 105	Re = 4 × 105
	u0 = 15 m·s−1	u0 = 21 m·s−1	u0 = 30 m·s−1	u0 = 45 m·s−1	u0 = 60 m·s−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. Basic calculation parameters.

Cylinder	D = 0.1 m
Flowing medium—air	density (constant) kinematic viscosity dynamic viscosity	ρ = 1.225 kg·m−3 ν = 1.5 × 10−5 m2·s−1 μ = ν · ρ = 1.8 × 10−5 kg·(m·s)−1

.

2.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:

$$\frac{\partial \rho }{\partial \mathrm{t}}+\frac{\partial \left(\rho \widehat{u_i}\right)}{\partial {\mathrm{x}}_i}=0,$$

(2)

momentum equation (Navier–Stokes):

$$\frac{\partial \left(\rho {\widehat{u}}_i\right)}{\partial \mathrm{t}}+\frac{\partial \left(\rho {\widehat{u}}_i{\widehat{u}}_j\right)}{\partial {\mathrm{x}}_j}-\frac{\partial \left({\sigma }_{ij}\right)}{\partial {\mathrm{x}}_j}+\frac{\partial \widehat{p}}{\partial {\mathrm{x}}_i}=\frac{\partial {\tau }_{ij}}{\partial {\mathrm{x}}_j},$$

(3)

where ρ [kg·m⁻³] is the density of the flowing medium, see Table 2, t [s] is time, ${\widehat{u}}_{i,j}$ [m·s⁻¹] represents the time average of the velocity components, ${\sigma }_{ij}$ is stress tensor, $\widehat{p}$ is the time average of the static pressure, ${\tau }_{ij}$ [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

$${\mathsf{\nu }}_t=\min \left[{\left(\kappa d_w\right)}^2,{\left(C_{Smag}\mathsf{\Delta }\right)}^2\right]\cdot S\cdot \{1-\exp [-{\left(y^{+}/25\right)}^3]\},$$

(4)

where $d_w$ is the wall distance, $S$ is the strain rate, $\kappa =0.4187$, and $C_{Smag}=0.2$ are constants, and $y^{+}$ 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:

$$\Delta =\min \left(\max \left(C_w\cdot d_w;C_w\cdot h_{max},h_{wn}\right);h_{max}\right),$$

(5)

where $h_{max}$ is the maximum edge length for a cell, $h_{wn}$ is the wall-normal grid spacing, and $C_w=0.15$ is a constant.

The boundary conditions for all simulations are given in

Table 3. Boundary conditions.

Boundary	Type of Boundary Condition
entry into computing area:	velocity inlet
output from computing area:	pressure outlet
vertical side walls:	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

.

2.3. 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 (u₀ = 0.35–15 m·s⁻¹) consisted of 1,138,960 cells. It has height of the first cell at the wall 1 × 10⁻⁵ 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 (u₀ = 21–60 m·s⁻¹) 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⁻⁶ 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:

$$y^{+}=\frac{\rho ·{\mathrm{y}}_p·u^{*}}{\mu },$$

(6)

where ρ is density of the flowing medium [kg·m⁻³], y_p is distance of the first cell point from the wall [m], u* is friction velocity [m·s⁻¹], μ is dynamic viscosity [kg·(m·s)⁻¹].

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.

3. Results

3.1. 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⁵, the distance of the minimum mean flow velocity u_x 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 u_x). 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³ and critical, Re = 1.4 × 10⁵ 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 u_x/u₀. 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.

3.1.1. 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³ and from the newer [7] for Re = 4 × 10³ or [23] for Re = 5 × 10³. 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 u_x/u₀∈(−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 u_x/u₀) is in these cases about x/D = 2.5. Only the source [22] presents a different course of u_x/u₀.

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 u_x values. In [7], the results for Re = 4 × 10³ 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 u_x/u₀ 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 u_x in the immediate vicinity behind the cylinder.

Another situation is from the point of view of determining u_x 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³ and Re = 4 × 10³). 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 u_x/u₀ = 0.46 and [7] even u_x/u₀ = 0.40. It can be said that all numerical simulations for Re = 4 × 10³ show for the distance x/D ≥ 3.5 a better agreement with [22], when the normalized mean stream velocity is about the value u_x/u₀ = 0.7.

3.1.2. Normalized Mean Stream Velocity in the Critical Region

CFD simulations in 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 u_x/u₀ = −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⁵ the experimental data [18,22] differ from CFD calculations (Figure 7).

From the point of view of defining the maximum value of u_x at a greater distance from the cylinder axis, the numerical simulations for the flow with Re = 1.4 × 10⁵ are compared with experimental data at a distance approximately x/D = 4.5 at the value u_x/u₀ = 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 u_x/u₀ = 0.81.

3.2. Pressure Load on the Cylinder Circumference—Pressure Coefficient c_p

The pressure load on the cylinder circumference for all presented simulations was evaluated using a dimensionless c_p coefficient. It is given by the ratio of the static pressure p_i and the dynamic pressure p_dyn 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

$$c_p=\frac{p_i}{p_{dyn}}=\frac{p_{ci}-p_{ref}}{1/2·\rho ·u_{ref}^2},$$

(7)

where p_ref is the static pressure at the reference point [Pa], p_ci is the resulting static pressure on the cylinder surface at the i-point [Pa], u_ref is the mean streamwise velocity at the reference point [m·s⁻¹].

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 c_p (maximum compressive load) for Re ∈ (1.3 × 10²; 2.1 × 10⁵). 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⁵ and from [27] for Re = 1.5 × 10⁵.

One of the other goals of this work is to compare the pressure coefficient distribution c_p 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 c_p in this article are divided separately for the subcritical and critical regime.

3.2.1. Pressure Coefficient c_p—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 c_p,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 c_p distribution for Re > 8 × 10³. However, the difference between the minimum value of c_p for Re = 3 × 10³ and Re = 8 × 10³ is evident. The CFD results also correspond to these results, there is a different minimum of c_p around Re = 3 × 10³ and Re = 2 × 10⁴.

The position of c_p,min of the presented Re coincides (φ_min ≈ 70°) and with increasing Re is shifting slightly nearer to the vertical axis of the cylinder (about 5°, φ_min ≈ 75°). For the lowest Re, presented in the subcritical region (Re = 3 × 10³), the results of the pressure coefficient of the experiment and the CFD (Re = 2.3 × 10³) are close and have the value of c_p,min ≈ −1.15. All other experimental results for higher Re in this region are c_p,min ≈ −1.25, which is again in agreement with the results of numerical simulations.

On the leeward side of the cylinder, the results of experiments and numerical simulations slightly differ. The best agreement with the experimental results is shown by the simulation for the value Re = 2.3 × 10³ (experiment Re = 3 × 10³), distributions are almost identical and the value c_p ≈ −0.93 is on the whole back side of the cylinder in the range φ ∈(90°,180°). In general, no values of c_p on the leeward side exceed the value c_p,min ≈ −1 in numerical simulations. Conversely, the pressure coefficient recorded in the experiments on the back of the cylinder (φ > 110°) showed with the increase in Re also increase in compressive load up to the value of c_p ≈ −1.2.

3.2.2. Pressure Coefficient c_p—Critical Region, Re ≥ 1 × 10⁵

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 Figure 8 and Figure 10, which differ significantly from each other, and the CFD results also differ from them. There is also an interesting comparison of c_p with standard results [31]. The minimum position for Re = 5 × 10⁵ according to [31] and CFD is φ_min ≈ 85°, however, they differ in the value. CFD calculations tend to overestimate the extreme pressure on the cylinder, for Re = 2.4 × 10⁵ is c_p,min,CFD = −2.6, while c_{p,min,standard} = −2.2. The position of the c_p minimum in experiments for different Re is approximately φ_min ≈ 65°, but the values of c_p,min differ significantly from each other. The CFD results are close to standard on the leeward side for φ > 135°.

3.3. Drag Coefficient c_d

In the present work the drag coefficient c_d 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 c_d 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 c_d in this flow regime is not fully clarified. In

Table 4. Values of drag coefficient c_d for selected Re; CFD.

cd—Subcritical Regime			cd—Critical Regime
Re = 2.3 × 103	Re = 4 × 103	Re = 2 × 104	Re = 1.4 × 105	Re = 2 × 105	Re = 3 × 105	Re = 4 × 105
1.05	1.02	1.02	0.90	0.82	0.45	0.25

and Figure 11 there are synoptic values of the cited experimental results, as well as the numerical simulations solved in this study.

4. Discussion

The paper is focused on the numerical simulation of the flow around the cylinder in the range Re ∈ (2.3 × 10³; 4 × 10⁵) 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 c_d, 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 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.