Numerical analysis of flow around a cylinder in critical and subcritical regime

sustainability Article Numerical Analysis of Flow Around a Cylinder in Critical and Subcritical Regime Ivan Kološ * , Vladimíra Michalcováand Lenka Lausová   Citation: Kološ, I.; Michalcová, V.; Lausová, L. Numerical Analysis of Flow Around a Cylinder in Critical and Subcritical Regime. Sustainability 2021,13, 2048. https://doi.org/ 10.3390/su13042048 Academic Editors: Marc A. Rosen and JoséAlvarez Received: 30 November 2020 Accepted: 10 February 2021 Published: 14 February 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Department of Structural Mechanics, Faculty of Civil Engineering, VŠB-Technical University of Ostrava, 708 00 Ostrava, Czech Republic; [email protected] (V.M.); [email protected] (L.L.) *Correspondence: [email protected] 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 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 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. Keywords: CFD; LES; circular cylinder; subcritical and critical regime; drag coefficient; lift coefficient; pressure coefficient; normalized mean stream velocity 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 – 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. Sustainability 2021,13, 2048. https://doi.org/10.3390/su13042048 https://www.mdpi.com/journal/sustainability Sustainability 2021,13, 2048 2 of 13 Figure 1. Buildings of circular cross-section: (a) [14], (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 ]. Re =u·D υ, (1) 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 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 crosssection, 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. Figure 2. Variation of cdand Flow Transitions for Single Cylinder Flow [16,17]. 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 cdin 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–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–23]. Sustainability 2021,13, 2048 3 of 13 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 – 20 ], and velocity profile in the wake behind the cylinder [ 7 , 21 – 23 ], these quantities are essential for wind engineering. Specifically, the flow around the cylinder in this research is done for Re = 2.3 ×103to 4 ×105 , 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 – 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 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. Table 1. Selected Re and corresponding velocities. Subcritical Re = 2.3 ×103Re = 4 ×103Re = 2 ×104 u0= 0.35 m·s−1u0= 0.6 m·s−1u0=3m·s−1 Critical Re = 1 ×105Re = 1.4 ×105Re = 2 ×105Re = 3 ×105Re = 4 ×105 u0= 15 m·s−1u0= 21 m·s−1u0= 30 m·s−1u0= 45 m·s−1u0= 60 m·s−1 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−5m2·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 macroand 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: ∂ρ ∂t+∂(ρˆ ui) ∂xi =0, (2) Sustainability 2021,13, 2048 4 of 13 momentum equation (Navier–Stokes): ∂(ρˆ ui) ∂t+∂ρˆ uiˆ uj ∂xj −∂σij ∂xj +∂ˆ p ∂xi =∂τij ∂xj , (3) where ρ [kg · m −3 ] is the density of the flowing medium, see Table 2,t[s] is time, ˆ ui,j [m · s −1 ] represents the time average of the velocity components, σij is stress tensor, ˆ p is the time average of the static pressure, τ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 νt=minh(κdw)2,CSmag∆2i·S· {1−exp[−y+/253]}, (4) where dw is the wall distance, S is the strain rate, κ= 0.4187, and CSmag = 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: ∆=min(max(Cw·dw;Cw·hmax,hwn);hmax), (5) where hmax is the maximum edge length for a cell, hwn is the wall-normal grid spacing, and Cw=0.15 is a constant. The boundary conditions for all simulations are given in Table 3. 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 × 2D(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%. Sustainability 2021,13, 2048 5 of 13 Figure 3. Unstructured tetrahedral mesh—the size of the cells progressively decreases towards the wall. Figure 4. Dimensions of cells in the prismatic layer at the wall of the cylinder for Re ≥1×105. Figure 5. Dimensions of the computational domain. 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 = 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 (u 0 = 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: y+=ρ·yp·u∗ µ, (6) Sustainability 2021,13, 2048 6 of 13 where ρ is density of the flowing medium [kg · m −3 ], y p is distance of the first cell point from the wall [m], u* is friction velocity [m·s−1], µis dynamic viscosity [kg·(m·s)−1]. 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 5 , 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 3 and critical, Re = 1.4 ×105 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. Figure 6. Velocity field for subcritical and critical regimes. Mean stream velocity: ( a )Re=4 × 10 3 ;u 0 = 0.6 m · s −1 ; (b) Re = 1.4 ×105 ;u 0 =3m · s −1 . RMS of mean stream velocity fluctuations: ( c )Re=4 × 10 3 ;u 0 = 0.6 m · s −1 ; (d) Re = 1.4 ×105;u0=3m·s−1. 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 0 .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. Sustainability 2021,13, 2048 7 of 13 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. 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 3 and from the newer [ 7 ] for Re = 4 × 10 3 or [ 23 ] for Re = 5 ×103 . Figure 7shows 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∈ ( − 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 0 )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 u x values. In [ 7 ], the results for Re = 4 ×103 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 0 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 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 u x /u 0 =0.46 and [ 7 ] even u x /u 0 =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. 3.1.2. Normalized Mean Stream Velocity in the Critical Region CFD simulations in Figure 7have 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 = − 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 ×105the 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 5 are compared with experimental data at a distance approximately x/D = 4.5 at the value u x /u 0 = 0.75. Sustainability 2021,13, 2048 8 of 13 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. 3.2. 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 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.8Din 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 cp=pi pdyn =pci −pre f 1/2·ρ·u2 re f , (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−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 c p (maximum compressive load) for Re ∈ (1.3 × 10 2 ; 2.1 ×105 ). 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 ×105and from [27] for Re = 1.5 ×105. Figure 8. Mean pressure distribution—experimental data; subcritical and critical regime (Norberg [19], Tani [29], James [27]). 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 cp—Subcritical Region It can be seen in Figures 8and 9that 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 3 . However, the difference Sustainability 2021,13, 2048 9 of 13 between the minimum value of c p 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 c p around Re = 3 ×103and Re = 2 ×104. Figure 9. Mean pressure distribution, subcritical regime, CFD and experimental data (Norberg [ 19 ]). 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 3 ), the results of the pressure coefficient of the experiment and the CFD (Re = 2.3 × 10 3 ) are close and have the value of c p,min ≈ − 1.15. All other experimental results for higher Re in this region are cp,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 3 (experiment Re = 3 × 10 3 ), 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 cp≈ −1.2. 3.2.2. Pressure Coefficient cp—Critical Region, Re ≥1×105 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 8and 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 5 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 5 is c p,min,CFD = − 2.6, while cp,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◦.