MODELLING OF EFFECT OF INFLOW TURBULENCE ON LARGE EDDY SIMULATION OF BLUFF BODY FLOWS

A random flow generation (RFG) algorithm for a previously established large eddy simulation (LES) code is successfully incorporated into a finite element fluid flow solver to generate the required inflow/initial turbulence boundary conditions for the LES computations of viscous incompressible turbulent flow over a nominally twodimensional circular cylinder at Reynolds number of 140,000. The effect of generated turbulent inflow boundary conditions on the near wake flow and the shear layer and on the prediction of integral flow parameters is studied based on long time average results. No-slip velocity boundary function is used but wall effects are taken into consideration with a near wall modelling methodology based on van Driest Damping approach. The numerical results obtained from simulations are compared with each other and with the experimental data for different turbulent inflow boundary conditions to assess the functionality of the RFG algorithm for the present LES code and hence its influence on the vortex shedding mechanism and the resulting flow field predictions.


INTRODUCTION
The turbulent flow field around a circular cylinder is a great practical importance for many engineering applications such as hydrodynamic loading on ocean marine piles and risers, bridges.Due to the occurrence of turbulence in the boundary layer the computational studies of the vortex shedding from the cylinder may require turbulence modelling [1].In the literature, there are a wide variety of numerical studies based on different types of turbulence modelling ranging from Reynolds averaged Navier-Stokes equations (RANS) based turbulence models [1] to LES [2][3] to predict the turbulent vortex shedding from a bluff body.Tutar and Holdo [3] conduct various transient simulations of the turbulent flow around a nominally two-dimensional circular cylinder and conclude that LES calculations yield much more realistic flow field prediction than the RANS based turbulence model calculations even for two-dimensional mesh systems including moderate spatial resolutions.In order to generate inflow turbulence, two groups of approaches in the numerical studies may be used.One is to conduct the auxiliary simulations of turbulent flow fields using LES approach [3] and to store the time series of fluctuating velocity components for inflow boundary conditions of the main simulations.In this approach, the velocity field extracted from a plane near the domain exit is basically rescaled and then is reintroduced as a boundary condition at the inlet.This approach will not be applicable to cases where scale laws are non-trivial.Another group of approach is to artificially generate time series of random velocity fluctuations by performing an inverse Fourier transform for prescribed spectral densities.There are several numerical studies with varying degree of success in this approach [4,5].This method may require a fairly lengthy development section as a transient region due to the fact that the synthetic velocity field generated by the random fluctuation approach may lack both turbulent structure and non-linear energy transfer.To eliminate the disadvantage of the method, Celik et al. [6] suggest a relatively simple random flow generation (RFG) algorithm, which can be used to prescribe inlet condition as well as initial conditions for spatially developing inhomogeneous anisotropic turbulent flows.The algorithm is based on the study of Kraichnan [7] and is successfully applied to particle dispersion modelling by Smirnov et al. [8].This is the approach used in the present study.

THE GOVERNING EQUATIONS
For the turbulent flow computations the filtered Navier-Stokes equations and the continuity equation for an incompressible fluid can be written as follows: ( ) ( ) where ρ is the fluid density, µ is the fluid dynamic viscosity, i u is the filtered value of the velocity, p is the filtered value of the pressure, ij τ is the sub-grid scale stress defined by To close the filtered equations ij τ is modelled using a viscous analogy as below where sgs k is the sub-grid scale kinetic energy, ij S is the strain rate tensor and t µ is the sub-grid scale viscosity term calculated by the Smagorinsky model [9] given as with l is the characteristic sub-grid length scale related to the with of the filter, ∆ used as The first term on the right-hand side of Eq. ( 6), the Smagorinsky constant s C , is chosen in a range from 0.10 to 0.30.The second term, the filter width, ∆ , as an indication of the characteristic length scale separates large and small scale eddies from each other and can considered to be an average cell size.It is calculated for two-dimensional elements in the present finite element method (FEM) based code [14] due to Tutar and Holdo [3] as Therefore the sub-grid length is calculated directly from the local grid size and the grid size distribution is thus very important for the present sub-grid scale (SGS) model.The calculation of the sub-grid viscosity is performed using a subroutine that returns an array of values at each integration point in the computational domain at each time step.

RANDOM FLOW GENERATION (RFG) ALGORITHM
The present RFG algorithm is a modified form of the previously developed technique [10] and is incorporated into a finite element code [10] with which the LES code is applied to generate a realistic inflow field.The RFG technique can be used to prescribe inlet conditions as well as initial conditions for spatially developing inhomogeneous, anisotropic turbulent flows [6].The fluctuating component is derived from the turbulence intensity and length-scales, provided by the turbulence model or via empirical relations.The procedure adopted here consists of the following steps.

Given an anisotropic velocity correlation tensor
of a turbulent flow field find an orthogonal transformation tensor ij a that would diagonalize ij r ( ) As a result of this step both ij a and n c become known functions of space.Coefficients play the role of turbulent fluctuating velocities ( ) in the new coordinate system produced by transformation tensor ij a .
2. Generate a transient flow-field in a three -dimensional domain the modified method of Kraichnan [7] ( ) (16) The procedure above takes as input the correlation tensor of the original flow field ij r and the information on length and time scales of turbulence.These quantities can be obtained from a steady state RANS simulations or experimental data.The outcome of the procedure is a time dependent flow field ( ) , and it is, to a high degree, a divergence -free inhomogeneous anisotropic flow field as shown by Smirnov et al. [8].The stress tensor of the generated flow field is equal to ij r .Due to the application of Eq. ( 11), spatial and temporal variations of i u follow Gaussian distributions with characteristic length and time scales of l ,τ ; however other distributions can be used to simulate different turbulence spectra.

MESH CONFIGURATION AND NUMERICAL DETAILS
A stationary two-dimensional circular cylinder with a non-dimensional unit diameter, D is situated in the centre of the vertical plane and is exposed to a timeaveraged free stream velocity, U.The Reynolds number based on the free stream velocity and the cylinder diameter is 140,000.Inlet (upstream), upper and lower boundaries are extended laterally to minimize the effects of the boundaries on the cylinder.The flow domain is also extended long enough downstream to eliminate the far field effects on the near wake and to produce full development of the vortex street.Then simulations are carried out for the mesh resolution of 200 x 175 and 250 x 215 and 335 x 300 grid points in the (x, y) direction to study the effect of mesh resolution on the flow results.As seen in Fig. 1 the computational domain contains grid points, which are non-uniformly spaced in the x-and y-directions.Very fine mesh resolution is used surrounding the cylinder vicinity.The mesh size is increased with the distance from the cylinder surface and the location of the first grid point from the cylinder wall is chosen to be 0.0015D.All simulations are carried out using a general purpose finite element solver [10].The whole computational domain is subdivided into a finite number of small two-dimensional four-noded finite elements.The non-linear equation systems resulting from the finite element discretization of the flow equations are solved using a segregated solution algorithm with a second-order trapezoid time integration scheme in order that the advection terms may have a higher order of accuracy and the solution may be stable.Non-dimensional time step, * t (= D t U ) where t is the time, is chosen at a value of 0.001.The iterative convergence criteria is set to Cylinder surface: No-slip velocity boundary condition is prescribed at the cylinder surface.A near wall methodology, which comprises the no-slip velocity boundary condition with a modified form of damped mixing length model, is used.

RESULTS AND DISCUSSION
The numerical simulations are performed in conjunction with the RFG algorithm for different turbulent inflow boundary conditions i.e. varying degree of turbulence at a Reynolds number of 140,000.At this Reynolds number the flow is in the sub-critical flow regime and the reported wind tunnel inflow turbulence of which isotropic turbulence level, u I is 0.6 percent and length scale of turbulence D l t of 0.02 [11].The simulations are conducted in two-dimensional manner as the turbulent vortex shedding from the cylinder is considered to be two-dimensional in the mean.The

Primary Flow Features
Basic flow features such as primary and secondary separations, vortices developing from the cylinder surface and the well known von Karman vortex street with periodic vortex shedding are observed for each simulation.For brevity, only those obtained from a simulation conducted for a mesh resolution of 335 x 300 grid points at turbulence intensity of 0.6 percent (Case 3) are presented in this section.Figure 2 shows the instantaneous streamline patterns for Case 3 at a different nondimensional time instant of 2.5 in the near flow field of the cylinder before the instability develops for the present flow.The secondary vortices which are successfully reproduced by the present LES code give rise to the turbulent momentum exchange towards free shear layers, Tutar and Holdo [3].

Influence of Mesh Resolution
The effect of mesh resolution on the present calculations is tested for three different mesh systems (Cases 1-3) containing 200 x 175, 250 x 215 and 335 x 300 grid points in the (x, y) direction with isotropic inflow turbulence level of 0.6 percent.Figure 5 illustrates the time-averaged streamwise velocity distributions in the near wake along a constant x/D = 1 position for all three mesh systems and that of experimental data of Cantwell and Coles [11].The time averaging is applied to the LES calculations for a period of 15 shedding cycles at the end of the each simulation.In general, it is seen in Fig. 3 that the agreement with the experimental data is satisfactory in consideration of two-dimensional simulations except for small differences in the regions of peak velocity and in the very near wake region near the centreline.LES calculations with each mesh system slightly underestimate the velocity profile in the free shear layer region while overestimating it in the recirculation zone.The higher streamwise velocity on the centreline of the cylinder indicates that the recirculation length in the LES calculations is much smaller than the experimental one and this can be attributed to two-dimensional simulations.The deviations from the experimental data decrease with the increase in mesh resolution.As the mesh resolution is increased, the vortex structures near the separation point becomes larger and the time averaged location of the separation point moves further downstream as noted in Table 1 which illustrates the improvement of pressure profiles and hence the time averaged base pressure and drag coefficient values with the higher mesh resolution.Fig. 2.The instantaneous streamline contours in the near wake of the cylinder for Case 3 at non-dimensional of 2.5

Influence of Inflow Turbulence
The influence of turbulent inflow boundary conditions on the present flow is extensively studied with varying degree of isotropic turbulence.With this regard, numerical simulations are performed in conjunction with the RFG technique to generate the isotropic inflow turbulence of which the turbulence level is varied from 0.6 to 6 percent.Figures 4(a) and (b) show the time history of non-dimensional x and y velocity components at the inflow boundary at two different inflow locations for varying degree of inflow turbulence.Figures 4(a) and (b) clearly demonstrate the effect of turbulence on flow at the inflow boundary.The frequency and amplitude of the velocity fluctuations increase as the inflow turbulence level is increased from 0.6 to 6 percent.The spatial variations (i.e.phase lag) observed in the history of each velocity component for each simulation case can be attributed to spatial inhomogenity that is satisfied at the inflow boundary for realistic flow simulations.The present RFG technique therefore ensures the spatially evolving isotropic turbulence along the inflow boundary.The corresponding one-dimensional energy spectrums of y velocity component for different inflow turbulence conditions at the inflow boundary and at two more selected locations along the centreline of the cylinder are computed.The energy spectrums are non-dimensionalized with respect to the total energy in each case, and are plotted with respect to the non-dimensional vortex shedding frequency U D f , where f is the frequency.Figures 5(a) and (b) show the energy spectrums at two different locations along the centreline of the cylinder for three cases (Cases 3-5).The energy spectrums are consistent with the instantaneous flow visualizations, where the presence of small scales is clearly observed even at this very low isotropic turbulence intensity value of 0.6 percent.Figure 5(a) shows the energy spectrum at a point P(0, 7D) at the inflow boundary for isotropic inflow turbulence of 0.6 and 6 percents.The energy spectrum is confined to a narrow band for Case 3 and a broader banded energy spectrum is obtained for Case 4. Both cases can capture the inertial sub-range as observed in Fig. 5(a).The energy spectra for the y velocity component appear to exhibit an inertial sub-range about a half decade of inertial sub-range from about non-dimensional frequency of 2 to 10 where energy spectra exhibit a slope close to -5/3 predicted by Kolmogorov theory.The energy spectrum determined at point P(8D,7D) behind the cylinder can be also seen in Fig. 5(b).The broader banded distribution of the energy at this in-wake location compared to that at the selected front stagnation point may be attributed to the occurrence of the near wake turbulence and the associated large irregular wake velocity fluctuations.Figure 5(b) clearly shows dominant frequencies at non-dimensional vortex shedding frequencies with different magnitudes for each case.It is clearly seen in Fig. 5(b) that when the isotropic turbulence intensity is increased from 0.6 to 6 percent (Case 4) the energy spectrum shows more energy at the larger frequencies in the inertial subrange and the distribution of the energy is much broader banded.The effect of inflow turbulence on the time averaged pressure distributions over the cylinder surface is illustrated in Fig. 6.As clearly seen in Fig. 6, the position of the minimum time averaged pressure coefficient, min P C moves slightly rearwards on the cylinder as the inflow turbulence increases for all turbulence inflow data generated by the RFG technique.The minimum pressure decreases and the base pressure increases with the increasing turbulence intensity.This is mainly due to the direct effect of the inflow turbulence on the boundary layer leading to changes in the position of the separation lines.

CONCLUSIONS
The following conclusions can be withdrawn from the present study: 1) The turbulent inflow data generated by the proposed RFG algorithm improves the flow resolution in the near wake and hence the predictions for the integral parameters.The RFG technique promotes the boundary layer characteristics, pressure distribution and other integral parameters with increasing level of inflow turbulence.
2) The position of the separation points shift rearwards, the base pressure increases and the minimum pressure, time averaged and fluctuating drag and fluctuating lift coefficients decrease with increasing inflow turbulence level at the investigated Reynolds number.The present model predicts this trend in accordance with the experimental observations.
3) The turbulent inflow data generated by the RFG technique generally broadens the vortex stretching frequencies but does not disrupt the vortex shedding mode at the investigated sub-critical Reynolds number of 140,000 with increasing turbulence intensity.
4) The differences in the energy spectra between simulations with and without RFG applications of different level of isotropic inflow turbulence are found to become smaller in the near wake of the cylinder compared to that of inflow boundary.

3 .
where l ,τ are the length and time-scale of turbulence respectively, ijk ε is the permutation tensor used in vector product operation, and distribution with mean M and the standard deviation σ .Numbers n j k , n ω represent a sample of n wavenumber vectors and frequencies of the modelled turbulence spectrum Apply scaling and orthogonal transformation to the flow field i v generated in the previous step to obtain a new flow field i u ( ) ( )

Table 1
Overview of calculated integral parameters at Re = 140,000