The topic is the comparison of the accuracy of 2D and 3D LES computations. In this article, turbulent vortex shedding flow past a 2D square cylinder, which is 2D in the mean, is discussed. When the time-averaged flow-field is 2D, it is not unusual for engineering problems to be solved by an unsteady 2D computation because CPU time is reduced greatly in comparison with 3D computations. It is clear that 2D LES is fundamentally incorrect from the view point of turbulence dynamics.
Here, 2D and 3D LES computations are compared, and the accuracy of 2D LES computation is examined. Here, the standard Smagorinsky type SGS model, shown in Table 3, was used. The value of the Smagorinsky constant Cs was 0.1 for 3D computation.
This value was increased by 50 % for 2D computation in comparison with that for 3D computation in order to compensate for the limited diffusion in the span-wise direction. Specifications for 2D and 3D computations are summarized in Table 4. The Reynolds number based on the cylinder width, H, and the velocity at the inflow boundary, U0, is 1.0 x 105. The details of LES computations are described in Note 1.
Mean Surface Pressure (Fig. 4):
Fig. 4 shows distributions of mean surface pressure coefficients <C̅p>t. Here < >t means time averaging and –indicate grid-filtering (cf. Nomenclature). The distribution of the 3D computation is similar to that of the experiment by Ohtsuki et al.(1978) except for the leeward face, while in the case of the 2D computation, the negative value of <C̅p>t becomes too large at the center of the side faces.
RMS Value of Fluctuating Pressure (Fig. 5):
Distributions of the rms pressure coefficients Cprms given by the two computations are compared with the results of experiments by Bearman et al.(1982) and Lee (1975) in Fig. 5. It can be concluded that the 3D computation reproduces the pattern of Cprms distributions in the experiments very well. In the 2D computation, the values of Cprms differ greatly at both the side face and the leeward face. This lack of agreement between the results of the 2D computation and the experiment is not acceptable.
Power Spectrum of Fluctuating Lift and Drag Forces (Fig.6):
Fig. 6 shows the computed and measured spectra of the fluctuating lift force FL and the drag force FD acting on the cylinder. As is shown here, the spectrum curves of FL have a sharp peak in the results of both computations as well as in that of the experiment. The frequency of each peak corresponds to the fluctuation caused by the corresponding the vortex shedding.
The spectrum shapes of FL and FD given by the 3D computation correspond very well with those given by the experiments. On the other hand, the peak frequencies of the fluctuations of FL and FD in the 2D computation do not correspond with those in the experiments, but are located at higher frequencies for both FL and FD. As is shown here, the 2D computation cannot reproduce accurately the surface pressure field and the resulting lift and drag forces.
Difference in the Flow-field Structures near the Side Face (Fig. 7):
Fig.7 compares the power spectra of u,’ based on computations near the side face. This figure reveals significant differences in the flow structures predicted by the 2D and the 3D computations. The spectrum of the 2D computation is confined to a narrow band, and has a sharp peak at the same frequency as that of FL.
On the other hand, the spectrum of the 3D computation extends over a wide band. From this figure, it becomes clear that the energy transfer mechanism throughout a wide spectrum range cannot be reproduced by the 2D computation.
This is hardly surprising when we recall that the vortex stretching mechanism, which is essentially three-dimensional and is not reproduced at all in the 2D computation, plays an important role in the energy cascade around this point. Hereafter, analyses by 3D LES computations are presented.
Comparison of Les with K-Ԑ/DSM Computations and Measurements for the Flow Past 2d Square Cylinder:
In the case of unsteady vortex shedding flows past a 2D square cylinder at high Reynolds numbers, stochastic 3D turbulent fluctuations are superimposed onto the 2D periodic vortex-shedding motion. This is illustrated in Fig. 8 where f is the instantaneous value of a quantity, f̅ the periodic fluctuation, f” the stochastic turbulent fluctuation, <f>t the time averaged value and <f>p the phase averaged value.
Proper testing of the accuracy of a numerical simulation for this type of flow is possible only when time-resolved measurements are available. These measurements give information on the temporal variation of the ensemble-averaged quantities and separate the fluctuations into periodic and stochastic turbulent ones. Lyn et al.(1992, 1994) carried out this kind of measurements on the flow past a square cylinder at Re = 2.2 x 104.
Here, we analyzed the same flow-field by LES (3D) computations. In the computations using RANS models, 3D stochastic turbulent fluctuation f” is averaged out by ensemble-averaging procedure and its effects are modeled, while 2D periodic fluctuation f̅ can be calculated explicitly.
Thus, the only 2D equations need to be solved in the computations using RANS models. Franke and Rodi (1991) carried out the 2D computations using RANS models for the same flow-field. They used the standard k-ԑ model and the standard Differential Second-moment closure Model (DSM) of Launder, Reece and Rodi (1975), with wall corrections to the pressure-strain terms due to Gibson and Launder (1987). For both models, two types of near-wall treatments are applied; wall function and the two-layer approach in which the near-wall region is resolved with a one-equation model due to Norris and Reynolds (1975).
The results given by LES are compared with those given by the RANS model computations of Franke and Rodi as well as with the experimental results. Here, the standard Smagorinsky type subgrid model, shown in Table 3, was applied. The outline of LES computations for this case is summarized in Note 2.
Comparison of LES with k-ԑ/DSM Computations and Experiments:
Values of various integral parameters (e.g., Strouhal number, drag coefficient, etc.) predicted with the various computational methods are compared in Table 5 with experimental values.
In the computations by Franke and Rodi(1991), results with the standard k-ԑ model using wall functions yielded a steady solution and no vortex shedding. Thus this case is not included in Table 5. The two-layer k-ԑ model yielded too low shedding frequency and drag coefficient <CD>t.
The two-layer DSM predicted too high Strouhal number and drag, while the DSM with wall functions predicted these parameters in good agreement with the measurements. The LES computation reproduced slightly lower values, which, however, still agree fairly well with the measurements.
Recently some effort has been made to develop the more refined k-ε models which can be applied to the vortex shedding flows such as the flow-field analyzed here. The results of the modified k-ԑ models, which are also shown in Table 5.
Time averaged velocity <u̅1> (Fig. 9 (a)):
Experimental data of Lyn (1992) and Durao et al.(1988) are also included in Fig. 9. In front of the cylinder, the results are not much influenced by the computation method used, but there are fairly large differences in the wake region. The k-ԑ model overestimates the length of the reverse flow zone considerably.
Both DSM variants predict too short separation, and there is little difference between the results obtained with the two versions of near wall treatment. The LES shows better agreement with the experiments than the DSM computations, although LES still underestimates the length of the separation region.
Kinetic energy of fluctuations <kto>t <ks>t (Fig. 9 (b), (c)):
As shown in Fig. 9 (b), the two-layer k-ԑ model underestimates severely the total fluctuation level behind the cylinder, while the DSM computations give approximately the correct level of distribution of total fluctuations. In front of the cylinder, the standard k-ԑ model yields an unrealistically high fluctuation level (here, only turbulent fluctuations exist) due to excessive turbulence production by this model in stagnation regions, as pointed out by Murakami et.al. (1990, 1992, 1993) and Mochidaetal.(1993).
The DSM does not have this problem. But, Franke and Rodi have shown that all turbulence models tested by them produce too low turbulent fluctuations, <ks>t, in the wake of the cylinder, as shown in Fig. 9(c). This means that the periodic fluctuations are greatly overestimated by the DSM.
As for the LES computation, it predicts the correct distribution of the fluctuations, but the level is too low in general for both <km>t and <ks>t. LES produces a higher turbulent fluctuation level, <ks>t, in downstream region, which would be nearly correct, but is too low by about a factor of 2 in the peak region around x1/H=2. Hence, both total and stochastic turbulent fluctuations are underestimated in the LES simulation. Individually, the LES results are both closer to the measured level than those predicted by the DSM.
Phase averaged velocity < U̅2> (Fig.10):
The agreement of wall function DSM and LES computations with the measurements is fairly good in general. On the other hand, the result of two-layer DSM slightly overestimates the value of <U̅2>p, in comparison with the results of LES and wall function DSM.
Phase averaged turbulent kinetic energy <ks>p (Fig.11):
Here, the LES computation appears to reproduce the observed behavior better than the computations with the DSM. The asymmetry in the contour lines predicted by the DSM is basically the same as in the experiments, but the maxima of <ks>p in the experimental data are located much closer to the centerline than in the computations.
This points to a weakness in the DSM, which appears to predict the centers of the shed vortices considerably farther from the centerline than the experiments indicate. The LES computation appears to predict the maxima of <ks>n closer to the centerline.
Performance of the Modified k- ԑ Models:
Recently, two different types of modified k-ε models were proposed.by Launder and Kato (1993) (hereafter LK model), and Przulj and Younis (1993) (PY model). Details of these models are given in Notes 4 and 5 respectively.
The results of these models are also included in Table 5. The consequence of the modification for Pk in the case of LK model is clear in the values of C̅L. There is an order of magnitude increase in C̅L resulting from the use of modified production as is shown in Table 5. The use of eqn.(8) in Note 4 in place of eqn.(7) gives remarkably closer agreement with the LES. On the other hand, C̅L given from PY model becomes too high in comparison with that of LES.
Fig. 12 (a) compares the time-averaged velocity <u1>t along the centerline. LK model (with Cµ(<S>p)) reproduces the size of the reverse flow zone rather well.
The total kinetic energy of fluctuations <kto>t is illustrated in Fig. 12 (b). <kto>t given from the standard k-ԑ model computation conducted by Launder and Kato (1993) shows excessively large production of turbulent fluctuations in the stagnation region as well as in the result with the standard k-ԑ model by Franke and Rodi(1991) (Fig.9(b), (c)).
Unfortunately this problem still exists in the result of PY model as is indicated in Fig. 13, while the accuracy for the prediction of <kto>t is improved remarkably in the stagnation region by LK model (Fig. 12 (b)). On the other hand, is somewhat underestimated in the wake region by LK model when Cµ is treated as constant.
Significant improvement is attained by varying the value of Cµ in space using eqn.(9) in Note 4 (Fig. 12 (b)). It is rather surprising that distribution of <kto>t given by PY model shows very close agreement with the experiment in the wake region although <kto>t is severely overestimated in the stagnation region, by this model (Fig. 13).
Les Computations of Unsteady Pressure Field around an Oscillating Square Cylinder:
The interactive phenomena between fluid and a bluff-shaped body represent one of the most difficult problems in the field of fluid dynamics. This section deals with the interaction between fluid and body, i.e., the flow-fields around oscillating square cylinders both in cases of forced oscillation and wind-induced free oscillation.
Here, computations based on LES (3D) using the standard Smagorinsky model are presented. A very simple procedure to incorporate the influence of body motion on the flow-field is used. The accuracy of the computations using this procedure is made clear by comparing them with the experimental results of Bearman et al.(1982) and Nakamura et al. (1975).
Outline of Numerical Method for Oscillating Cylinder:
Fig. 14 illustrates the flow conditions treated here. The cylinder oscillates in two different ways. One way is by forced oscillation and another is by wind-induced free oscillation. In forced oscillation, the cylinder is forced to oscillate in a sinusoidal motion in the x2 (lateral) direction.
In the case of free oscillation, the cylinder is spring-mounted in the x2 direction and oscillates by the lift force acting on it. Table 6 lists the cases analyzed one case of the fixed cylinder (case 0), fourteen cases of forced oscillation (cases 10-16, cases 20-26) and four cases of free oscillation (cases 30-33).
In the forced oscillation cases, two different values for the amplitude (A/H= 0.10, 0.25), and nine different values for the reduced velocity (U0/(NH)= 7.0- 12.0) are imposed. In the experiment, the lock-in phenomena occur for all the cases considered herein. As for the cases of free oscillation, the reduced velocity is varied from 9.0 to 12.0, and the Scruton number (Sc= 2mδ/ρH2L3) is set at 1.36. In this range, the self-induced oscillation of the cylinder was observed in the experiment by Nakamura et al.(1975).
The governing equations are given in Tables 7 and 8. The standard Smagorinsky model is used for turbulence modelling of sub-grid scale. In the cases of the fixed cylinder and forced oscillating cylinder, the Reynolds number based on H and U0 is 2.2 x 104, which corresponds to that in the experiment by Bearman et al. (1982). As for the case of the free oscillating cylinder, the Reynolds number is 5 x 104, which approximately corresponds to that in the experiment by Nakamura et al.(1975). Numerical methods and boundary conditions are summarized in NOTE 6.
In the present computation, the whole mesh layout (i.e., the coordinates) is treated to move together with body motion. In order to incorporate the influence of body motion, a technique is applied to simplify the computation.
Instead of moving the body, here the fluctuating inertia force term, -xC, whose magnitude is equal to the acceleration of the body but with the opposite sign, is added to the transport equation of u̅2 as an external force at every grid point of the fluid for every time step, including points at the inflow boundary (cf. Tables 7 and 9). Equivalent form is easily obtained, using the ALE (Arbitrary Lagrangian-Eulerian) method proposed by Hirt et al. (1974), by letting the spatial derivative of xc be 0.
Results and Discussion:
(1) Instantaneous Velocity and Pressure Fields for cases of Positive and Negative Aerodynamic Force Damping (Fig. 16).
The instantaneous velocity and pressure fields around forced oscillating cylinders are illustrated in Fig. 16 for cases of 10 (A/H=0.1, NH/U0 = 1/7.0) and 14 (A/H=0.1, NH/U0 = 1/9.0). The figures show the instantaneous flow-fields when moving velocity, xc, has its peak value in each case. Cylinders are moving to the left side (i.e., upper side in Figs. 16(a) and 16(b)).
In case 10, large negative pressure occurs at the right side of the cylinder (lower side in Fig. 16(a)). On the other hand, it occurs at the left side in case 14 (Fig. 16(b)). This means that the pressure field would damp the cylinder oscillation in case 10 but amplify the oscillation in case 14.
The significant difference between the results of both cases is also observed more clearly in the difference of the phase angles between the cylinder displacement and the pressure acting on the side faces of the cylinder for two cases.
(2) Frequency of Lift Force Acting on the Cylinder. Fig. 17 compares the lift force acting on the fixed cylinder (case 0) and the forced oscillating cylinder (case 23), wherein the imposed body frequency is 0.12 (1/8.5=NH /U0).
The Strouhal number (St)f, which is evaluated from the time history of the lift force, is 0.14 (1/7.3) for the case of the fixed cylinder while it decreases to 0.12 (1/8.5), the value of the imposed body frequency for the case of the oscillating cylinder. This change in the vortex shedding frequency by body motion is known as the lock-in phenomenon, which is reproduced successfully in the present computation. The lock-in phenomenon is also confirmed in all other cases conducted here.
(3) Phase Angle between Body Motion and Fluid Force. The time histories of the displacement of the forced oscillating cylinder and the negative pressure acting at the center of the side face of the cylinder for case 23 are shown in Fig. 18. The phase angle between the two is evaluated as about -30° by Fourier analysis.
In Fig. 19, the phase angles for the other thirteen cases at A/H = 0.10 and 0.25 are also plotted. The predicted phase angles for these cases correspond rather well with the experimental values. When the phase angle is within the range 0° <ф< 180°, the aerodynamic force damping is negative. This means that the lift force would amplify the cylinder oscillation.
As shown in Fig. 19, the phase angle for case 10 (A/H = 0.10, U0/(NH)=7.0) is about -80° which would damp the oscillation, while the angle for case 14 (A/H = 0.10, U0/(NH)=9.0) is about 80° and the aerodynamic force damping is negative. Thus, the oscillation would be amplified by the aerodynamic force in case 14. The aerodynamic force damping of positive and negative observed in the experiment is reproduced well in the computation.
Wind-Induced Free Oscillation:
Fig. 20 shows the oscillation amplitudes of the cylinders in the cases of free oscillation (cases BOSS), where the cylinder is spring- mounted and oscillates within x2 direction by the lift force (cf. Table 8). In the experiment, self (vortex)- induced oscillation starts when the reduced velocity, U0/NH, is about 8 and the oscillation amplitude of the cylinder increases rapidly according to the increase of reduced velocity. The computational results reproduce this tendency very well.
The numerical prediction of turbulent flow-field around a square cylinder is one of the most difficult problems in the field of CFD at this stage, since the flow-field is very complicated, characterized by the impinging, the separation at the frontal corner with steep streamline curvature, the recirculation behind the cylinder and the vortex shedding, etc.
The numerical analyses of vortex shedding flow with various turbulence models are presented here. The standard k-ε model did not succeed in reproducing the vortex shedding at all. However, the modified k-ε model succeeded in reproducing it well, with an accurate prediction of some integral parameters.
DSM reproduced the vortex shedding rather well, but the results show many discrepancies between the experiments, e.g., the overestimation of periodic fluctuation. In general, LES shows best agreement with the experimental results in many respects, even though the standard Smagorinsky SGS model is used. However, recent progress in wind engineering require more accurate SGS models for LES. The dynamic SGS model is expected to become a powerful tool for flow analysis in wind engineering problems in near future.
The improvement of accuracy usually requires the increase of CPU resources and computation instability. Thus a strategy for selecting a turbulence model may be proposed according to the purpose of the analysis and the accuracy required. The relative abilities of various turbulence models are compared in Table 10. LES is most suited to analyze the complicated flow situations usually treated in the applications of wind engineering, but LES needs a great deal of CPU time. Further efforts should be made to develop more reliable turbulence models and more simple and less time-consuming CFD methods.
Outline of 2D and 3D LES computations of flow past fixed square cylinder at Re = 1.0×105:
No velocity fluctuation is imposed at the inflow boundary. The turbulence levels are 0.04 % in the experiment of Bearman (1982) and 0.2 % in the experiment of Ohtsuki (1978), respectively. A staggered grid was adopted and a second-order centered difference scheme was applied for the calculation of spatial derivatives. For time advancement, the Adams-Bashforth scheme was used for both convection and diffusion terms. For boundary conditions at solid walls, the generalized logarithmic law expressed by
was employed in order to estimate the time-averaged wall shear stress, < τw>t. Furthermore, the following relation was used to specify the instantaneous wall shear stress τw.
In the 3D computation, the computational domain covers 32.85H x 11.0H x 2.0H in the stream-wise (x1), lateral (x2) and span-wise (x3) directions, respectively. 99 x 63 x 10 grid points are placed in each direction. For the 2D computation, the computational domain covers the same length of stream-wise and lateral directions, and the same number of grid points are placed in both directions, but the computation has no dimension in the span-wise direction.
The interval of time advancement is 2×10-4 in non-dimensional time scale, which is the scale based on H and U0. Starting from the initial conditions, the computations required 5.0 x 105 time steps (100 in non-dimensional time scale) to reach stable periodic conditions. The computations were then continued further for 2.05 x 106 steps (410 in non-dimensional time scale). Thus statistically reliable averaged values were obtained.
Outline of LES computations of flow past fixed square cylinder at Re=2.2×104:
The value of the Smagorinsky constant Cs was 0.13 for this case. Because the Reynolds number, which is 2.2 x 104, is lower than that in the previous section, the grid scales are multiplied by the Van Driest type wall damping function, 1-exp(-xn+/25), in order to account for the near wall effect.
For the computations conducted here, a second-order centered difference scheme proposed by Piacsek and Williams (1970), which insures spatial conservation of quadratic values, was adopted for the spatial derivatives. For time advancement, the Adams- Bashforth scheme was used for the convection terms and the Crank-Nicolson scheme for the diffusion terms. The interval for time advancement is 1.0 x 103 in non-dimensional time scale based on U0 and H.
The computational domain covers 4.5H upstream of the cylinder, 14.5H downstream and 6.5H on either side of the cylinder. The size of the domain in the plane perpendicular to the cylinder is exactly the same as that in Franke and Rodi(1991) ‘s computation. 104 x 69 grid points were placed in the stream-wise and lateral directions, respectively.
The grid interval hw adjacent to the cylinder wall is 0.022. This value is the same as that for Franke and Rodi’s wall function case. Franke and Rodi set = 0.00125 in the two-layer approach. At the inflow boundary, the approach flow was set to be constant and uniform and no velocity fluctuations were prescribed in the LES computations, while the turbulence level of the approach flow was 2% in Lyn’s experiments and Franke and Rodi’s computations.
At the outflow boundary, zero gradient conditions were imposed. Symmetry conditions were employed for the lateral boundaries, and periodicity conditions were imposed for the boundary planes perpendicular to the cylinder axis.
For the boundary condition at the solid walls, Werner and Wengle(1991) ‘s approach assuming a linear or 1/7 power law distribution of the instantaneous velocity is adopted:
This approach was chosen because the Reynolds number, which is 2.2 x 104, is lower than that for the previous case.
The computations required 1.5 x 105 time steps (150 in non-dimensional time scale) to reach stable periodic conditions. The computations were continued further for 1.0 x 105 steps (100 in non-dimensional time scale) covering 13 shedding periods. The time averaged values were determined by time averaging over these 13 periods as well as by averaging over the span-wise direction.
We could sample data necessary for the <ks>t computation only in the area close to the cylinder. This was due to the limitation of the computer memory.
Modified k- ԑ model proposed by Launder and Kato (1993) (LK model).
Launder and Kato defined dimensionless strain and vorticity parameters:
The production term of turbulent kinetic energy can be written as follows using the definition of <S>p in eqn.(5):
In LK model eqn.(7) is replaced by the following equation:
The values of <S>p become very large in the stagnation region and the very large values of <S>p lead to the excessive level of Pk when the original form of eqn.(7) is used. However, the deformation near a stagnation point is nearly irrotational, i.e., <Ω>p=0.
Thus the replacement Pk=Cµԑ<S>p<Ω>p leads to a marked reduction in turbulence production Pk near the stagnation point while having no effect in a simple shear flow. Further modifications have been proposed by Craft et al. (1993) involving a dependence of Cµ on <S>p.
Kato and Launder (1993) carried out a computation using the model by Craft et al.(1993) in which Cµ is treated as a function of <S>:
Modified k-ԑ model proposed by Przulj and Younis (1993)(PY model).
Pizulj and Younis (1993) added the time-dependent source term to the e equation as follows:
Where q is the mean kinetic energy 1/2<ui>p<ui>p. A new coefficient Ct was introduced; its value was set as 0.26, assigned by computer optimization.
Numerical methods and boundary conditions for LES computations of flow past oscillating square cylinder at Re = 2.2×104.
The size of the computational domain and the boundary conditions are summarized in Table 9. The number of grid points is also indicated in the Table. Grid arrangements are illustrated in Fig. 15. The minimum grid interval adjacent to the solid wall hw was set at 0.025. For the boundary conditions at the solid wall, Werner and Wengle(1991) ‘s approach is adopted.
The interval of time advancement is 1×10-3 in the time scale which is made dimensionless by H and U0. In the case of the fixed cylinder, the computation required 1.0×105 time steps (100 in non-dimensional time scale) to reach stable periodic conditions, starting from the initial conditions.
As for the oscillating cylinder, starting from the result of the fixed cylinder, 1.5×105 steps in the case of forced oscillation and 2.5×105 steps in the case of free oscillation are required to reach stable periodic conditions respectively.
The computation is then continued further for 1.0×105 steps in the case of the fixed cylinder and for 0.7×10s– 1.2×10s steps (70-120 in non-dimensional time scale) covering 10 times of period of the cylinder oscillation in the case of the oscillating cylinder. Thus, statistically reliable averaged values are obtained. The computation for one case takes 20-60 CPU hours on a Fujitsu VP2600 machine.