1 Introduction

Incompressible turbulent flows around a single square cylinder (SC) are basically taken as a crucial academic problem in the society of applied mathematics and mechanics,since in fact there exist long-time interests of deeply perceiving the characteristics and properties of the turbulent SC flows.A brief literature review of the works done before the year 1990 was given previously^[1],which was prepared for the project of larger mechanics problem in ocean engineering financially supported by national science foundation of china (NSFC) in the latter of 1980s.Here,we trace the research history from two lines before presenting the objective of this paper. 1.1 Experiments

The earlier encompassment was reported in 1966 by Vickery^[2],where the influence of incom- ing stream turbulence intensity on the lift and drag on a long cylinder of square cross-section was explored.Similarly,some studies were done to seek the turbulence effect on drag of rectan- gular cylinders^{[3, 4, 5, 6]}.Some measurements of SC wake were carried out by Okajima^[3] and Davis et al.^[7].The flapping shear layer formed by flow separation from the forward corner of an SC was studied experimentally by Lyn and Rodi^[8].By virtue of one-component laser-Doppler velocimetry,they also investigated the associated recirculation region on the sidewall formed in flow separation from the forward corner of an SC.It was suggested that the recirculation effect can be a possible explanation for some flow properties.

On the other hand,Gu and Sun^[9] carried out wind tunnel tests to measure the interference between two circular cylinders (CC) in staggered arrangement at high subcritical Reynolds numbers (i.e.,2.2 × 10⁵ and 3.3 × 10⁵),and found three main flow patterns.They also mea- sured the flows around three parallel CC in equ_ilateral-triangular arrangements by wind tunnel tests^[10].By measuring pressure distributions on three CC at Reynolds number 5.5 × 10⁴ and visualizing flow regimes at Reynolds number 1.4×10⁴,when the spacing ratio is over the range 1.7≤6 N/d 6≤2.2,where N and d are respectively the distance between the centers of adjacent cylinders and the cylinder diameter,there are three basic interference types: interference of proximity,shear layer reattachment,and wake.As the incident angle varies from zero to 60°, these flow patterns can be subdivided further.

Recently,using dye- and laser-induced fluorescence visualizations,the spanwise and stream- wise vortex structures were captured by Luo et al.^[11].It was found that the critical Reynolds numbers were in the range from 160 to 200,at which mode transition from the mode A to B occurs.For the flow pattern mode A or B,the spanwise wavelengths are slightly different from that observed in CC flow.Highlighting the effect of incoming flow turbulence,after particle imaging velocimetry (PIV) measurements,So et al.^[12] found that the incoming flow turbulence increases the vortex shedding frequency in the post-lock-in region but does not change the system natural frequency.The so-called lock-in point is the point at which the frequency of vortex-shedding frequency coincides with the system natural frequency.

For wake flows,many concluding remarks can be found in recent literatures.For instance, to study the vortical structures and the interactions in the turbulent wakes of three side-by-side circular CC,Zhou^[13] measured the turbulent wake using various techniques,including PIV, hot-wire anemometry,laser Doppler anemometry,and laser-illuminated flow visualization.The cylinder centre-to-centre spacing T/d is arranged to be 1.5,at which both wide and narrow wakes exist.It can be observed that the base pressures of the upper and lower cylinders are identical but smaller than that of the central cylinder.Consequently,the gap flows between cylinders are deflected towards the upper and lower cylinders,respectively,forming one wide wake behind the central cylinder and one narrow wake on each side of the wide wake.It was indicated that while the vortical structures in the narrow wakes may be generated from the shear layer separation from the upper or lower cylinder,those in the wide wake may originate from the shear layer instability.

For the near wake around a finite-length SC,with one end mounted on a flat plate and the other free,an experimental study was carried out by Wang and Zhou^[14].The cylinder aspect ratio,or height-to-width ratio H/d ranges from 3 to 7.The Reynolds number based on d and the free-stream velocity is 9 300.Using the hot-wire anemometry,laser Doppler anemometry, and PIV,measurements were carried out mainly in a closed-loop low-speed wind tunnel.It was reported that the instantaneous flow structure around the cylinder is arch-type regardless of H/d,consisting of two spanwise vortical “legs”,one on each side of the cylinder,and their connection or “bridge” is near the free end.

For Reynolds number in the range of (2 500,6.0 × 10⁴),Alam et al.^[15] measured the CC wake flows using two tripwires attached at azimuthal angle alpha = ±10° to 70° measured from the forward stagnation point.Based on the measurements and those in literature,five flow regimes (A - E) are classified.In particular,it was proposed that the bi-stable phenomenon (regime D) has its possible connection to the mechanism of well reported rain-wind-induced cable vibration.

A systematic experimental study of the flow past two side-by-side square cylinders was carried out by Alam et al.^[16].Using different techniques,including hot wires,load cell,PIV and laser-induced fluorescence flow visualization,the flow experiments were performed at Re = 4.7 × 10⁴ and a changeable cylinder centre-to-centre spacing ratio (T/d) ranged from 1.02 to 6.00.Based on the flow structure and the Strouhal number (St),four distinct flow regimes and their corresponding T/d ranges are identified.The results also unveil a crucial role played by the gap flow and its passage geometry in contributing to the differences observed.

In addition to the significant insights from experiments,the onset of instability in the flow around an SC with linear stability analysis was studied by Kelkar and Patankar^[17].It was reported that the critical Reynolds number is inferred from the variation of the rate of growth of the perturbations with the Reynolds number.For unsteady flow at a Reynolds number above the critical value,the calculated Strouhal number is very close to its experimentally observed value.Furthermore,a stability analysis revealed that there is a critical Reynolds number beyond which the SC wake flow turns to three dimensional regime,as reported by Robichaux et al.^[18]. On the other hand,for CC wake flow,the corresponding vortex dynamics has been reviewed by Williamson^[19].Both extrinsic and intrinsic features are described in detail. 1.2 Simulations

Using different turbulence models,numerical simulations of vortex shedding past a free- standing SC at Re = 2.2 × 10⁴ were done by Bosch and Rodi^[20].Using wall functions,the standard k - ε model was compared with a modification suggested by Kato and Launder^[21]. It was shown that the main quantities of engineering interest can be predicted reasonably well with k -ε type models.The application of the KatoõLaunder modification can bring a major improvement.

A direct numerical simulation (DNS) of three-dimensional flow around an SC at Reynolds numbers ranged from 150 to 500 was carried out by Sohankar et al.^[22].In the DNS,the grid is non-staggered,and the pressure Poisson equation is solved by a multi-grid approach,implying that some numerical disturbances should have been unexpectedly introduced into the DNS software,for which,some explanations can be found in the previous work^{[23, 24]}.

Similarly,to predict the appearance of three-dimensional phenomena in the flow past an SC numerically,the DNS results for the scenarios of Reynolds numbers from 125 to 500 were reported by Saha et al.^[25].It was found that the values of the Strouhal number and the time-averaged drag-coefficient are closely associated with each other and can reflect the spatial structure of the wake.It can be noted that their results are based on the MAC algorithm of Harlow and Welch^[26].Noticeably,in the numerical approximation,both the convective and diffusive terms are dealt with the central difference scheme.

To understand air pollutant dispersion among high-density,high-rise bu_ildings,Niu and Zhu^[27] reported the numerical results of three-dimensional flows around two identical square cylinders in staggered arrangements at Reynolds number 250 and zero incident angle,when the angle between the incoming velocity vector and the line connecting the centers of the cylinders is 45°.

Niu et al.^[28] reported the numerical results of convective heat transfer from two identical square cylinders (TISC) arranged at close proximity in a uniform cross air-flow at Re=250.It can be found that depending on the TISC arrangement,a bifurcation sometimes occurs,with the bifurcated lower Strouhal number value corresponding to the shedding from the outer shear layer,and the higher one corresponding to the shedding from the inner shear layer.

The numerical scheme in the previous studies^{[22, 25, 27]} can be further improved by using a 4th order upwind scheme to deal with the nonlinear convective terms in the governing equations. Actually,in a staggered grid system,numerical approximation by interpolation should adhere to the law of Taylor expansion,otherwise,the discrepancy in numerical solution can occur,because in numerical scheme design the implicitly-introduced artificial-viscosity due to truncation error is comparatively high.

In general,for turbulent flow simulations,three strategies: the phenomenological models, the DNS,and the large eddy simulations (LES) are widely used.The first strategy was reviewed by Hanjalic^[29].It was based on Reynolds averaging,but the lacking of universal property for the model coefficients implies that these coefficients examined in some benchmark situations might be unsu_itable to other relatively complicated engineering problems.

The DNS,as reported in Refs.^{[30, 31, 32]},in the traditional meaning,should use extremely fine grids so that the total spatial grid numbers should be as large as Re^2.25,where Re is the Reynolds number.This indicates that super computers must be used in DNS of fully developed turbulent flows.

The LES,as reported in Refs.^{[33, 34]},uses the grid-filtering to explicitly simulate the large motion,while the small eddy effect on the large eddy motion should be modeled.The LES works reported previously^{[35, 36, 37, 38, 39, 40, 41, 42]} indicate that this is an appropriate strategy for predicting turbulent flows.

For the LES,a series of sub-grid scale (SGS) models were developed^{[35, 43, 44, 45, 46, 47, 48]},among which are the earlier SGS model of Smagorinsky^[35],the dynamic SGS model given by Germano et al.^[43],a modification of dynamic SGS model^[44],the SGS model based on the Kolmogorov equation for resolved turbulence^{[45, 46]},the Vreman model based on algebraic theory,and the recent SGS model based on dynamic estimation of Lagrangian time scale developed by Verma and Mahesh^[48].However,it should be noted that the SGS models reported so far have not ascertained under what flow conditions the LES can lapse into the DNS.

It should be mentioned that even in LES the unclosed stresses are usually modeled by SGS-viscosity models,alternatively,to smooth the dynamics of the Navier-Stokes equations, another way is to directly modify the non-linear convective terms^[49].Thus,models based on these regulations occur^{[50, 51, 52, 53, 54, 55]}.

Turbulent flow is intermittent in time and space,indicating that at a local point the flow may be swirling or non-swirling,depending on if the swirling strength value is non-zero,or zero. The SGS model described in this paper can extent reflect this flow property.The results of the LES are compared with existing experimental data. 1.3 Objective

This paper presents the swirling-strength based on the LES results of the incompressible zero-incident flows around a single SC at low Reynolds numbers in the range from 10³ to 10⁴. The objective is to explore the SC flow characteristics,primarily related to three aspects: (i) the Reynolds number dependence,(ii) the distribution of (t-z) averaged vorticity component and sub-grid viscosity ratio,and (iii) the distribution of the (t-z) averaged Kolmogorovmicro-scale. 2 Governing equations

At higher Reynolds numbers,turbulent effect is significant,so the SGS effect on large scale motions should be carefully considered.In previous LES works^{[36, 37, 38, 39, 40, 41, 42]},sub-grid stress is usually assumed to be the product of sub-grid viscosity and filtered strain rate,which is just an analogy of the constitutive relationship for laminar Newtonian flu_id flows.This indicates that these SGS models have not satisfactorily ascertained at what Reynolds number they are applicable, because the sub-grid viscosity still exists,even when the flow is switched to a laminar regime and the swirl strength is generally zero.

In the earlier work of Smagorinsky^[35],the sub-grid viscosity is assumed to be proportional to the modular of local strain rate of filtered velocity field.But unfortunately,using the three- dimensional velocity fields based on flow instability analysis^[56],numerical tests revealed that the peak values of the modular of strain rate occurs at regions without swirling.To overcome this ambigu_ity,considering that turbulent flows have coherent structures and strong intermittency, here we propose a relatively fresh SGS model,where the sub-grid viscosity is supposed to be proportional to the FSI and local swirling strength,

Here,C_μ (= 0.09) is an artificially-defined constant,λ_ci is the swirling strength of filtered velocity gradient ∇u,the FSI,f_I,is the ratio of swirling strength λ_ci to the magnitude of the complex eigenvalue λ (= λ_cr + iλ_ci) of ∇u,i.e., where δ is the length scale,defined by harmonic-average of grid intervals, where δ_i denotes the grid interval in the x_i direction.If the filtered velocity gradient ∇u is denoted by its eigenvalue λ must satisfy the characteristic equation where with |A| representing the determinant of A.If the turbulent flow is not affected by the Coriolis force due to coordinate rotation,at a local point if the roots of (5) are real,and the swirling strength equals zero,it indicates that the vortices certainly have not arrived at the point, implying that the local flow regime is temporally laminar.On the other hand,when the local flow is in a turbulent regime,it indicates that there are vortices at the local point,suggesting that the roots must have two conjugate complex eigenvalues denoted by λ_cr ± iλ_ci,here,i (= ) is the unit of imaginary number.The iso-surfaces of swirling strength were used to illustrate vortices in turbulent flows previously^{[57, 58, 59]}.

Because the sub-grid viscosity ν_s has a factor f_I,the decaying factor of the length scale (1-exp(-y⁺/25)) near the wall as used in previous study^[36] is abolished in the present study, where y⁺ = yu_τ/ν,y denotes the normal distance to the wall,with u_τ denoting the mean friction velocity and ν representing the kinematic viscosity of fluid.

Since the swirling strength can occur in laminar flow regime when there is the effect of Cori- olis force due to coordinate rotation,in order to see the peculiarity of the present SGS model,a comparison of the present model with that of Vreman^[47] is therefore given in Figs.1(a)-1(b), where two top-views of the contours of sub-grid viscosity are illustrated.It can be seen that the sub-grid viscosity based on the present model is zero if the flow is non-swirling,while sub-grid viscosity based on the Vreman model does not hold this peculiarity.

Consider the turbulent SC flows in a Cartesian coordinate system,as schematically shown in Fig. 2,in which x is the horizontal coordinate,while y and z denote the vertical and spanwise directions.The origin (O) is allocated at the left-bottom corner of the SC.The spanwise coordinate is given merely by the symbol z.Let the flu_id kinematic viscosity be ν,for the incoming flow velocity u_in,the Reynolds number of the SC flow can be defined by Re = u_ind/ν.

We select u₀ = u_in and the SC cross-sectional side length d as the velocity and length scales,hence the time and pressure scales should be t₀ = d/u₀ and ρu₀ ² .Here,ρ is the flu_id density.The dimensionless governing equations for the incompressible SC flows can be written as follows:

For simplicity of presentation,here we have omitted the over bar for filtering of variables.The normalized total viscosity γ₁ can be represented by R denotes a force depending on viscosity gradient,

The solutions of (7) and (8) are sought under appropriate boundary and initial conditions. For the boundary conditions (BC) on the SC walls,non-slip BC is used,which implies

while for the BC at the outlet,similar to the treatment published elsewhere^{[22, 25]},we use the Orlanski^[60] type In the normalized form,we choose u_c = 1.At the inlet section,we choose However,on the lower and upper side boundaries of the computational domain,we use The initial condition is given by The periodical condition is used on the spanwise boundaries,which can be simply expressed as where B is the spanwise (z) length of the SC (see Fig. 2). 3 Numerical method 3.1 Solution procedure

The governing equations (7)-(8) of the SC flow problem were discretized by a finite difference method in a staggered grid system,where the convective terms were treated using the 4th order upwind scheme^[61].It can be noted that by making some proper changes,the existing numerical methods^{[24, 26, 62, 63, 64, 65, 66]} should also be applicable.

The solution procedure is based on the accurate projection algorithm PmIII developed by Brown et al.^[67].The detailed steps of the procedure were described elsewhere^[68].Now just a brief description is described below.

Let the intermediate velocity vector,the pressure potential,and the time level be u,φ,and n,respectively.Define H = (u · ∇)u - R,then let

We can calculate u by and calculate pressure p by where the pressure potential φ must satisfy the Poisson equation and the terms at the level of (n + 1/2) are calculated explicitly using a second order Euler prediction-correction scheme in temporal.As mentioned above,the nonlinear convective terms in the governing equations are spatially discretized by the 4th upwind finite difference scheme given in Ref.^[61],with the viscous diffusion terms by the second-order central difference scheme.

The pressure potential Poisson’s equation is solved by the approximate factorization one (AF1) method^[69] at first,and then by the stabilized bi-conjugate gradientmethod (Bi-CGSTAB) given by van der Vorst^[70] to improve the solution accuracy.In the AF1 application,for the SC wake flow simulation,the numerical experience reported elsewhere^[1] is useful for the nu- merical stability of pressure potential iteration.The criterion for pressure potential iteration was chosen so that the relative error defined previously^[71] should be less than 3 × 10^-8.The implicit second-order Crank-Nicolson for the right-hand side diffusion terms are used in spatial discretization. 3.2 Scheme improvement

Considering its intrinsic potential in numerical simulation^[24],the staggered grid system is used.we have therefore arranged u_ijk,v_ijk,and w_ijk at the points P₁(x_-1/2+i,y_j,z_k), P₂(x_i,y_-1/2+j,z_k),and P₃(x_i,y_j ,z_-1/2+k),respectively,with pressure potential φijk set at point P(x_i,y_j,z_k).Based on this stagger arrangement,the non-linear convective terms in the momentum equations are spatially discretized with a 4th order upwind scheme,in which those terms at the level o(ε4) are neglected.To briefly describe this scheme,we just give an expression for uux in the u-momentum equation.Adhering to the upwind idea,omitting the subscripts of j and k,for uux at P₁,when u(x_-1/2+i,y_j,z_k) = u_i > 0,we have

otherwise, The scheme coefficients can be obtained by the Taylor expansion,for u_i > 0, where δ₁ = +_x-1/2+i - x_-1/2+i-1,δ₂ = +_x-1/2+i - x_-1/2+i-2.On the other hand,when u_i≤6 0,we have where δ₁ = -x_-1/2+i+x_-1/2+i+1,δ₂ = -x_-1/2+i+x_-1/2+i+2.For vuy and wuz,i.e.,the cross-convective terms in the u-momentum equation,some descriptions were reported elsewhere^[61]. Similar numerical treatments can be encompassed for the v- and w-momentum equations.It is necessary to mention that to comparative-rigorously adhere to the upwind concept,the value of u_x⁽³⁾ must be yielded at x_-1/2+i-1/2 for u_i > 0,or at x_-1/2+i+1/2 for u_i < 0.Correspondingly, some Taylor-expansion-based modification for the scheme coefficients (A,B,C,and D) is requ_ired. In particular,for numerical treatment of the divergence of intermediate velocity vector ∇·u and the ∇²φ at the time level n,the influences of heterogeneous stagger grid (i.e.,φ_i is arranged between u_i and u_i+1) are taken into account to improve the accuracy of spatial discretization. 4 Results and discussion

In terms of the numerical method described above,the LES of the single SC flows as shown schematically in Fig. 2 was encompassed using an Inter-Cored personal computer with a memory 3.2 Gb and a CPU frequency 3.30 GHz.Correspondingly,the geometrical parameters for computational domain were given in Table 1.Each case requ_ires a CPU time of about 68 hours.Classified by Re₃ = Re/10³,6 cases were studied,as shown in Table 2.For each case, the spanwise length of the SC (B) is set as 11.

The grid in z -direction is uniform,with grid number N_k=41.While the grids in x- and y- directions are non-uniform,the grid-numbers are 191 and 145,respectively.The grid-partition uses a geometrical-series algorithm.The grid,in principle,should be denser near the SC walls. That means the grids near SC corners must be denser.The assignment of grids on the SC section walls has used two key parameters named as the amplifying-factor(1/0.88),the series termination number (14).The total number of vertical grid-line passing through the SC cross- sectional side is 34,as indicated in the range of x ∈ ^{[0, 1]} in Fig. 3(a).

The grid heterogeneity near the SC can be seen in Figs.3(a) and 3(b).In the bottom and top sides of the SC,the y-grid is assigned with an amplifying-factor 1/0.85 and a series- termination number 17.While for x-grid assignment,the relevant values are 1/0.828 and 17 for the region upwind the SC,but 1/0.935 and 47 for the region downstream the SC.Using the geometrical-series algorithm for grid partition,it seen that the finest grid distance to each SC wall is 4.165 × 10^-3.Such a grid arrangement leads to a time step Δt = 2.125 × 10^-3,which is smaller than those reported previously^{[22, 25, 27]}.

Just in the xy-plane,to show the distribution of Kolmogorov micro-scales,which are ob- tained by time-spanwise (t-z) averaging,we use a subscript “m” to represent (t-z) averaged variables,and a superscript (′) to represent the root mean square values (RMS) of the variables. Especially,the subscript “m” was omitted for (t - z) average lift (CL) and drag (CD),just for the convenience of representation.For Re ∈ [10³,10⁴],the vortex motion can result in relatively complicated flow patterns,which should be explored carefully and appropriately.It is noted for the 6 cases,the grid arrangement and the time step are unchanged.In post-processing of the LES output data,to remove the effect of initial condition,merely the data in the time period of t ∈ ^{[68, 170]} are used.The aspects of flow induced forces,and flow patterns are reported and discussed in the following subsections. 4.1 Flow characteristics

For turbulent SC flows,the characteristics can be recognized by observing the CD and CL coefficients,their mean and RMS values.The RMS values should be corresponding to the time averaging for the temporal data train obtained from the CD and CL coefficients based on z- average.Because the SC flows have fixed vortex separation points,the vortex shedding should be different from that in a CC flows.

To illustrate the Reynolds number dependence,the (t - z) averaged CL,CD,(ν_sr)pm,their relevant RMS values denoted by CL′,CD′ and (ν_sr)′ _p,and the Strouhal number are given in Table 2,Figs.4(a)-4(b),and Figs.5(a)-5(b).As shown in the third column of Table 2,for Re₃ = Re/10³ ∈ [2.5,10],the mean CD is about 1.825 (±0.052),with the maximum of relative deviation being 2.83%.This means that over the range of Re₃ ∈ [2.5,10],the drag coefficient is almost certainly Reynolds-number independent.For Re₃ ∈ [1,2.5],also given in the third column of Table 2,it can be seen that the present λ_ci-based LES indicates that with the increase of Reynolds number,the mean CD increases gradually from 1.605 to 1.775.Furthermore,as seen in the second column of Table 2,the mean CL is approximately zero.Also,this Reynolds- number dependence of mean CD and CL,can be observed in Fig. 4(a).

It should be mentioned that for Re = 500 (Re₃ = 0.5),the mean CD is 1.84,or 2.17 as respectively reported in the DNS^{[22, 25]}.This suggests even if the numerical results are based on DNS,the difference still exists due to the different solution procedure.On the other hand,the earlier measurements^{[4, 5, 6]} indicated the mean drag coefficient CD is around 2.0.As reported previously^[1],for turbulent flow around a single SC,the mean CD value depends on the turbulence intensity of incoming flow,an increase of the intensity can cause a decease of the mean CD.

The Re dependence of relevant RMS values denoted by CL′,CD′,and (ν_sr)′ _p can be seen from the 5th,6th,and 7th columns of Table 2.Different from the work of Bosch and Rodi^[20], these values are obtained on the basis of the current three-dimensional (3D) LES,where the sub-grid viscosity is calculated in terms of (1).Corresponding to the values given by the 5th and 6th columns of Table 2,the CL′ and CD′ are plotted as functions of Re₃,as seen in Fig. 4(b). Coarsely seen,in the low Re number range of Re₃ ∈ ^{[1, 10]},the CL′ is approximately 0.2,with CD′ having a value around 0.1.

It is noted that the peak value of viscosity ratio (ν_sr) is sought in the whole 3D computational domain.Because the turbulent SC flow is intrinsically dynamic,the peak value point must be varied temporally.The mean peak value denoted by (ν_sr)pm is merely based on t-average.From Fig. 4(c),it can be seen that the mean peak value of viscosity ratio (ν_sr)pm and its RMS value exhibit growing trend with Reynolds number,but the growing rate for (ν_sr)pm is larger.For instance,at Re₃ = 10,the (ν_sr)pm value is 17.16 with its RMS value equal to 3.673,as can be seen from the 7th line of Table 2.

Crucially,the Strouhal number has an evident independence of Reynolds number,as seen in 8th column of Table 2.For Re₃ ∈ ^{[1, 10]},St is about 0.112 8 with a relative deviation less than 4.5%,which can also be seen clearly from Figs.5(a)-5(d),in which the power spectra at Re₃ = 1,2,2.5,and 10 for parts 5(a)-5(d) are obtained by discrete Hilbert transform as done previously^[72].This Strouhal number is very close to the St number 0.124 for the case of Re = 2.2 × 10⁴ (or Re₃ = 22) as predicted by Bosch and Rodi^[20],when the two-dimensional flow calculation was based on the two-layer (TL) approach with Standard k -ε equation when the inlet viscosity ratio (r_μ = ν_t/ν) was set as 100.Note ν_t is the eddy viscosity which depends on kinetic energy k and its dissipation rate ε.As reported^[20],the Strouhal number measured by Bearman and Trueman^[4] is 0.123 for the SC flow at Re ≈ 5 × 10⁴ (or Re₃ ≈ 50) when the inlet flow turbulence level T u (= was less than 1.2%.For Re = 500,the previous DNS^{[22, 25]} indicated that the Strouhal number value is 0.122,or 0.116 respectively,which shows that the present LES has predicted the Strouhal number agreed qu_ite well with that on the basis of DNS of the single SC flow at Re = 500.

The relevant evolutions of lift and drag at Re₃ = 2.5 can be seen in Figs.6(a)-6(b),from which a visual comparison with c-DNS can be made.Obviously,for Re₃ = 2.5,the SC flow is turbulent,the use of λ_ci-based sub-grid viscosity model has not brought large discrepancies of the evolution curves of z-averaged lift and drag coefficients.

This comparison again provides numerical-solution evidence for the Chinese Society of Me- chanics to believe that the nonstandard analysis theory of turbulence^{[73, 74, 75]},is to some extent correct,and can be applied to get a higher computational efficiency as reported elsewhere^{[61, 68]}. This nonstandard analysis theory allows the use of c-DNS to obtain relatively accurate results, unnecessary to follow the DNS grid number requ_irement: be at the level of Re^2.25.

The properties of the peak value of viscosity ratio can be ascertained by reading in the evolution given by Fig. 7,from which we can see that both the mean peak value of viscosity ratio (ν_sr)pm and its RMS value have growing tendency,as seen directly in Fig. 4(c).The values of (ν_sr)pm given by Fig. 7 implies that there exists a consistency when those given by Table 2 were compared just for the cases of Re₃ = 2.5 and 10.

This consistency can also be evidently seen from Figs.8(a)-8(b).As given by Table 2,it is seen that at Re₃ = 2.5,CD′ ≈ 0.1; while at Re₃ = 10,CD′ = 0.128 8.In Fig. 8(b),it easily seen that the oscillation of CD in the period of t ∈ (140,170) does have a larger fluctuation magnitude.

4.2 Flow patterns

The flow patterns vary temporally and spatially^{[76, 77, 78]}.The illustration of the turbulent flow patterns is useful to deepen the understanding of the zero-incident incompressible flows around a single SC at low Reynolds numbers located in the range Re ∈ [10³,10⁴].

When t = 170,the contours of vorticity component ω₁ at x = 3,6,9,and 12 for Re₃ = 2.5 are shown by Figs.9(a)-9(d),i.e.,the top-layer,with the relevant ω1 contours for Re₃ = 10 given by Figs.9(e)-9(h),i.e.,the bottom-layer.A comparison between the layers of Fig. 9 reveals that the secondary flow pattern are Re dependent.While a comparison between parts 9(a)-9(d) in the top-layer or parts 9(e)-9(h) in the bottom-layer indicates that the secondary flow pattern is x-dependent.

There is an intrinsic correlation between the secondary flows (see Fig. 9) and the main stream flow patterns shown by contours of ω₃ as seen in Figs.10(a)-10(b).Note that the calculation of flow vorticity components has utilized the advantage of staggered grid arrangement.For instance,the (ω₃)_ijk stored at the point (x_-1/2+i,y_-1/2+j ,z_k) can be predicted more conve-niently.

4.3 (t - z) averaged variable distributions

Even though at the peak value of viscosity ratio can change its spatial position as a result of large eddy motion in the 3D computation domain,exploring where the peak value could appear should be useful for understanding the characteristics of turbulent SC flows at low Reynolds numbers.Therefore,the (t-z) averaged viscosity ratio ν_sr were shown by the contours labeled by 0.05,0.5,1,1.25,1.5,as seen in Figs.11(a)-11(b).In red colored region,the (t - z) averaged viscosity ratio ν_sr is beyond 1.5.While in the blue colored region,its value is less than 0.05.As seen in Figs.11(a)-11(b),for Re₃ = 2.5 and 10,the peak occur possibly in the two downstream regions,i.e.,E_Ⅰ {x ∈ (2,4); y ∈ (-1.5,-0.5)},and E_Ⅱ {x ∈ (2,4); y ∈ (1.5,2.5)}, where the coherent interaction of large eddies are more intensive.

The swirling-strength based LES of turbulent SC flows has predicted the Kolmogorov micro- scale on the basis of estimating the kinetic energy dissipation rate,similar to the work re- ported previously^[68].Let λu be the Taylor micro-scale,with respect to the work of turbulence introduction^[79],we can express the Kolmogorov miscro-scale η_u as below

where k = is the time and z-averaged turbulent kinetic energy,with u′ _i being the instan- taneous fluctuation of velocity component in the x_i direction,while ε is the dissipation rate of k given by with s_ij denoting the fluctuation of strain rate of velocity field.

The (t - z) averaged Kolmogorov micro-scale field is shown in Figs.12(a)-12(b) for two Reynolds numbers labeled by Re₃ = 2.5 and 10.It seen that the micro-scale is certainly small, but even smaller in the near SC wall region.The logarithm of Kolmogorov micro-scale to the base of 10 (log(η_u)) is less than -2.6 in the near SC wall regions.With respect to the definition of η_u given by 25,this property of log(η_u) distribution indicates that near the SC solid walls, the dissipation rate of turbulent kinetic energy are larger because of the flu_id-solid interaction. From Figs.12(a)-12(b),there exists an asymmetry of distribution of log(η_u),suggesting that in the SC flow the vortex shedding from the two fixed points are not equ_ivalent,which may lead to asymmetrical distribution of calculated k and its dissipation rate ε^[80].

4.4 Instantaneous iso-surfaces

The iso-surfaces of FSI in the computational sub-domain {x ∈ (0,7),y ∈ (-2,3),z ∈ (5,7)} are shown in Figs.13 (a) and 13(c),for t = 170.As seen in (2),the iso-surfaces of FSI are intrinsically dependent on the vortex motion,which can be quantified by swirling strength λ_ci^[57] .The FSI gives rise to some laminated structure involved with vortex shedding,in addition to relatively small vortex structure in the SC wake.The distribution of FSI iso-surfaces is dependent on the Reynolds number,since the FSI iso-surfaces shown by Fig. 13 (a) has some differences from that shown by Fig. 13 (c).

The iso-surfaces of viscosity ratio ν_sr are shown in Fig. 13 (b) and 13(d).They to some extent reflect the reliability of (t - z) averaged distribution given by Figs.11(a) and 11(b). 5 Conclusions

A swirling-strength based LES was developed for zero-incident turbulent flows around a single SC at low Reynolds numbers ranged from 10³ to 10⁴.The LES is based on a finite difference method which uses the Taylor-expansion to improve the accuracy of discretization. For the case of Re = 2 500,a comparison with c-DNS was made,so that some numerical evidence for the nonstandard analysis theory of turbulence^{[73, 74, 75]} can be found.The present LES results have the following findings:

(Ⅰ) The Strouhal number has an evident independence of the Reynolds number.It has a value around 0.112 8,with a relative deviation less than 4.5%.The mean drag coefficient is 1.825,also almost independent of Reynolds number for Re ∈ [2 500,10⁴],with a maximum relative deviation 2.83%.

(Ⅱ) There exist coherent vortical structures in the zero-incident turbulent SC flows at the low Reynolds numbers.The peak value of sub-grid viscosity ratio and its RMS values grow with the Reynolds number.The two downstream regions where the coherent interaction of large eddies are more intensive are found possible for the peak to occur.

(Ⅲ) The logarithm of Kolmogorov micro-scale to the base of 10 (log(η_u)) is less than -2.6 in the near SC wall regions.This indicates that near the SC solid walls,the dissipation rate of turbulent kinetic energy are larger.The instantaneous FSI exhibits some laminated structure involved with vortex shedding,in addition to relatively small vortex structure in the SC wake.

(Ⅳ) The nonstandard analysis theory of turbulence is to some extent correct. 6 Numerical test of sub-grid viscosity

To compare the present SGS model with that proposed by Vreman^[47],we choose a three- dimensional flow field based on the analysis of thermal instability of a layer of flu_id heated from below with the effect of rotation^[56],the steady flow field has the following velocity components:

In the numerical test,we simply set a_x = a_y = a/,a = π,W₀ = 1,and the Taylor number ,where Ω is the angular speed of rotation around the z-axis, d₀ is the depth of the flu_id layer,ν is the fluid kinematic viscosity.The domain is given by { x ∈ (0,2),y ∈ (0,2),z ∈ (-1,1)},which is partitioned by 51 × 51 × 51 grids.The grids are distributed uniformly in each direction.The calculated contours of the sub-grid viscosity are shown in Figs.1(a)-1(b).