#### Abstract

An innovative computational model is presented for the large eddy simulation (LES) of multidimensional unsteady turbulent flow problems in complex geometries. The main objectives of this research are to know more about the structure of turbulent flows, to identify their three-dimensional characteristic, and to study physical effects due to complex fluid flow. The filtered Navier-Stokes equations are used to simulate large scales; however, they are supplemented by dynamic subgrid-scale (DSGS) models to simulate the energy transfer from large scales toward subgrid-scales, where this energy will be dissipated by molecular viscosity. Based on the Taylor-Galerkin schemes for the convection-diffusion problems, this model is implemented in a three-dimensional finite element code using a three-step finite element method (FEM). Turbulent channel flow and flow over a backward-facing step are considered as a benchmark for validating the methodology by comparing with the direct numerical simulation (DNS) results or experimental data. Also, qualitative and quantitative aspects of three-dimensional complex turbulent flow in a strong 3D blade passage of a Francis turbine are analyzed.

#### 1. Introduction

An accurate prediction of turbulent flow inside/over a complex geometry has been one of the most important issues in fluid mechanics. Although formulation of mathematical models to simulate numerically such complex flows is a challenging task, many researches have been developed and some reliable results have been obtained; for example, Klaij [1] developed a space-time discontinuous Galerkin Finite Element Method (DG-FEM) for complex flows, Moin [2] and Conway et al. [3] simulated incompressible flows in complex geometries using large eddy simulation (LES), Zhang et al. [4] investigated the flow through the blades of a swirl generator using LES, and Sagaut [5] simulated the turbulent flow in a true 3D Francis hydroturbine passage considered fluid-structure interaction.

Flows with high Reynolds numbers, where the influence of the different turbulence scales must be taken into account, cannot be solved by direct numerical simulation (DNS) due to the large amount of data and unknowns involved in the computational solution. It is only suitable to be used in low Reynolds number flow and for a relative simple flow passage until now. In present, turbulent flows were often governed by the Reynolds Averaged Navier-Stokes (RANS) equations in many industrial fields and solved for the mean flow field together with a chosen turbulence closure model. This approach is based on the separation of the instantaneous value of a specific flow variable in its mean value and fluctuations with respect to this mean value. The well-known Reynolds stress components are originated substituting mean values and fluctuations of the variables in the conservation equations. But the RANS equations have more unknowns than equations, and, for this reason, it is necessary to use closure models to define the Reynolds stress components, in which there are many artificial parameters so that the applied fields are limited.

Alternatively, large eddy simulation (LES) may be used to analyze complex turbulent flows [6β8]. Using a grid filter, eddies of different scales are separated of the large eddies and subgrid-scales. Large eddies are associated to the low flow frequencies, and they are originated by the domain geometry and the boundaries. Subgrid-scales (SGSs) are associated to high frequencies and they have an isotropic and homogeneous behavior, maintaining their independence with respect to the main stream. Then, in LES the large eddies are simulated directly, whereas SGSs are simulated using closure models. It cannot only improve the computational accurate comparing with the traditional RANS model, but also reduce the demanding of the computational resource relative to DNS.

At present, numerous researchers have explored the different applications of LES with finite element method (FEM) [9] for example, Popiolek et al. [10] adopts the finite element analysis of laminar and turbulent flows using LES and subgrid-scale models simulating the two- and three-dimensional flows in a lid-driven cavity and over a backward-facing step and Young et al. [11] numerically simulated the turbulent flows in external field based on the LES with the Smagorinskyβs SGS eddy viscosity model using three-step FEM-BEM model. Also, various forms of FEM formulations are widely available in some of the literature [12β15]. Commonly used Galerkin schemes have limitations to effectively deal with the convective terms, and hence other forms of FEM like Petrov-Galerkin formulations [16] and Taylor-Galerkin schemes [17] were developed. Jiang and Kawahara [18] developed a three-step finite element formulation for the solution of an unsteady incompressible viscous flow based on the Taylor-Galerkin scheme. Of these works, it seems that the three-step finite element formulation for LES is more practical and provide possible routes to the solution of turbulence transient modelling problem in complex geometries.

In many engineering fields, accurate predictions of complex transient turbulent flows, which may be defined as a three-dimensional flow with highly disordered, intermittent, and rotational fluid motion, become more and more important. Thus, the applicability of DSGS model to complex transient turbulent flow combined with three-step finite element method should be examined. A very well-known benchmark for validating the methodology is turbulent channel flow and flow over a backward-facing step [19, 20]. Therefore, in this work, using three-step finite element method, we apply the DSGS model to LES of turbulent channel flow and turbulent flow over a backward-facing step to examine their performances in transient flow. Furthermore, a three-dimensional complex transient flow in a blade passage of a Francis turbine is simulated.

#### 2. Description of DSGS Model

Applying a grid filter to the continuity and momentum equations, respectively, the following filtered expressions for a Newtonian incompressible fluid are obtained: where and are the filtered velocity components and pressure, respectively, is the molecular kinematic viscosity, and is the components of the SGS stress tensor, which is given by In (2.1)β(2.3), the overbar indicates filtered quantities and represent components of the large turbulence scales, while in (2.3) represents components of the small turbulence scales or subgrid-scales.

The eddy viscosity model is frequently used to represent the effects of SGS in LES. is assumed as a nonlinear function of the strain rate, and it may be written as follows: where the filtered strain rate component and the eddy viscosity are given by In (2.5), is the Smagorinskyβs constant, which has values varying between 0.1 and 0.2 and is the grid filter width (representing a length scale). Usually, in three-dimensional flows , but in the finite element context may be taken as the cubic root of the element volume.

In the dynamic subgrid-scale model, formulated first by Germano et al. [21] and modified after by Lilly [22], values of have variations in space and time. These values may be calculated in a systematic way, without any interference of the user. The dynamic subgrid-scale model is characterized by two filtering processes: in the first one, using the grid filter, the filtered expressions are given by (2.1)-(2.2), where the SGS Reynolds stress was included. In the second filtering process, a test filter is applied, and the corresponding equation are given by where indicates the application of the second filtering process using a test filter, and is component of a stress tensor, with its component given by Taking into account (2.4) and (2.5), may be expressed as and the SGS eddy viscosity in the dynamic subgrid-scale model are defined by where .

Using (2.3) and (2.8), it is obtained that: From (2.4), (2.5), (2.9), and (2.11), the following system of equations is obtained: where Lilly [22] solved the system of equations given by (2.12) using the least squares minimization, obtaining the following expression for : This model has some important characteristics. (a) The eddy viscosity is equal to zero in laminar flows. (b) The eddy viscosity may take negative values, simulating the energy transfer from small to large scales (the backscatter phenomenon). (c) The model has an appropriated asymptotic behaviour near the solid boundaries.

The computation reveals that the local coefficient often yields a highly oscillating eddy viscosity field including a significant partition with negative values, which destabilizes the numerical calculation. To this end, a homogeneous coefficient is often used in the practical simulations. is computed with the requirement that the SGS dissipation of the resolved kinetic energy in the whole computational domain remains the same as with the local coefficient , that is, where denotes space averaging over the entire domain. is used to define the SGS turbulent viscosity, , in (2.10) and in the momentum equations.

The following dimensionless forms of the variables are used, namely: where is the characteristic length and is the maximum flow velocity. Substituting (2.9) into (2.7) and using (2.10) and (2.16), now (2.6) and (2.7) can be written as (after dropping the asterisk from the dimensionless variables for brevity from now on) where and .

#### 3. The Finite Element Algorithm

In the present model as mentioned earlier, a coupled three-step FEM approach is used in the solution of the governing differential equations in the LES formulation. The numerical formulation is briefly described in this section.

In the present model, the mass momentum (2.2) is approximated using an explicit three-step FEM based on a Taylor series expansion in time and standard Galerkin FEM in space, or the so-called Taylor-Galerkin method as proposed by Jiang and Kawahara [18]. Applying the three-step FEM scheme to (2.18), the three-step scheme can be obtained as far as time discretization is concerned.

*Step 1. *

*Step 2. *

*Step 3. * are the apparent velocities from which the velocities in the present time step can be derived as,
The superscripts and on the variables and represent these values of and at time and , respectively.

Applying the classical Galerkin method for space discretization, the following matrix expressions are obtained for (3.1)β(3.4), respectively

*For Step 1, *

* For Step 2, *

* For Step 3, *

Equation (3.4) can be discretized as where in which , , and are the shape functions.

Then the flow is analyzed by the following algorithm: (1) Determine with (3.5). (2) Determine with (3.6). (3) Calculate with (3.8). (4) Determine with (3.7).

We can gain the solution of only if the is got. We implement the divergence to (3.4) and using the incompressible conditions at time , the pressure Poisson equation is derived to correct the velocity equation as In the present model, the pressure poisson equation is solved using Bi-stable Conjugate Gradient method [23].

#### 4. Numerical Applications

##### 4.1. Turbulent Channel Flow

Turbulent flow through a plane channel has been widely considered as a benchmark for validating turbulence models and numerical methods. Reynolds numbers of 590 based on the channel half height and friction velocity are considered. The computational parameters for large-eddy simulations are summarized in Table 1. Coarse-mesh (), medium (, and fine-mesh () statistics were collected over 50, 20, and 15 flow-through times, respectively. One flow-through time is the domain length in the streamwise direction divided by the bulk velocity.

Figure 1 shows plots of the mean streamwise velocity for different meshes in plus units, as a function of , where , and, , represents the average of space-time. The mean velocity is in excellent agreement with these DNS results of Moser et al. [24], and all results are in agreement with the well-known law of the wall. In the mean velocity profile, the upward shift in the log-layer compared with the DNS data a little overpredicts the mean streamwise velocity in the core region of the channel. Our experience has shown that the addition of dissipation, by including a traditional LES subgrid-scale model such as the dynamic Smagorinsky model, leads to higher values of the mean velocity in the core region and in some cases leads to an overprediction of the mean velocity, which is also in agreement with the results of [25, 26].

The root mean squares (RMS) of velocities in our numerical model are plotted in Figure 2. It is defined as , , and . The RMS velocity profiles are slower to converge to the DNS results of Van Der Vorst [23] than mean velocity profiles; furthermore, our numerical model leads to an overprediction of the peak value and an little under-prediction of the peak and values. In whole, the RMS velocity is agreement with these DNS results of [24] well.

##### 4.2. Flow over a Backward-Facing Step

###### 4.2.1. Computational Setup

The simulation was carried out in the same configuration as the experiments of Vogel and Eaton [27]. Further data on these experiments are presented in [28, 29]. A schematic layout of the simulation domain is shown in Figure 3. The channel expansion ratio was 1.25, with a Reynolds number of 28β000 (based on the freestream velocity and step height, ). The experiment was carried out with an inflow condition consisting of two developing boundary layers separated by a relatively undisturbed core. These boundary layers had a measured thickness of . The total domain size used for the computations was , which included an entry length of . A grid containing and a finer grid containing nodes were used, which was stretched in the wall-normal and streamwise directions using hyperbolic tangent functions to cluster grid points at the step edge and in the wall boundary layers. The grid stretching can be observed in Figure 3.

Due to the need to supply a time-varying turbulent inflow condition, a time-series obtained from a separate periodic channel flow simulation was used at the inflow plane. A forcing method [30] was used to force the periodic channel flow simulation to match the experimental results for the mean and fluctuations of streamwise velocity. A convective boundary condition [31] was used at the exit plane. Statistics were averaged in the homogeneous spanwise direction and over 80 (for the grid) or 60 (for the finer grid) flow-through times.

###### 4.2.2. Results and Discussion

Table 2 summarizes the time-averaged mean reattachment lengths obtained for the two grid simulations, the simulation of Akselvoll and Moin [32] and the experimental data. The present results using the dynamic model containing nodes and the results obtained by Akselvoll and Moin [32] are within the estimated experimental error bounds, while the present results using nodes are just outside these bounds.

The computed coefficient of friction along the lower wall is compared to the experimental results in Figure 4. The coefficient of friction is defined as , where is the shear stress at the wall and is the freestream velocity. There are no known differences between the two grid simulations. The results show a similar agreement with the experimental results as the simulations by Akselvoll and Moin [32]: good agreement upstream of and from reattachment to , but poor agreement in both the recirculation zone and downstream of . The reason for the poor agreement downstream of is probably the effect of the outflow boundary condition. In the recirculation zone, it is unclear why all the LES simulations predict a larger negative value of . Akselvoll and Moin [32] noted that this could be caused by either the inflow generation method used or by inadequate grid resolution in this region. Results obtained using the finer grid are also shown in Figure 4 and indicate that the agreement with the experimental results is not improved with grid refinement. Another cause could be the limited spanwise extent of the domain, however, some preliminary simulations in wider domains did not result in improved results.

Figure 5 shows the mean streamwise velocity at a number of locations downstream of the step. Both grid simulations show generally good agreement with the experimental results with the major differences occurring downstream of reattachment and in the recirculation region. The coarse grid results show a slightly stronger reversed flow in the recirculation zone comparing with the finer grid, but both grid simulations show generally a little lower of the experimental results. Near reattachment region both grids results show smaller velocities than the experiments, which it is attributed to the flow gaining momentum as it passes downstream (due to side-wall boundary-layer growth) in the experiments, however, downstream of reattachment region, the LES results of finer grid is showed that velocities are good agreement with the experimental results in this region. The primary discrepancy with the experimental results is in the recirculation region (at = 3.2 and 5.9), where subgrid-scale model predicts a stronger backflow in the boundary layer, linked to the prediction of a larger coefficient of friction in this region, which is caused by the no slip wall condition and the boundary layer thickness decreasing on the backward step wall attacked by the upstream flow, also, the physical dissipation of SGS model and numerical dissipation may be affected by the results for the complex turbulent flow in the recirculation region.

Profiles of the RMS of the resolved streamwise velocity fluctuations are shown in Figure 6. Overall, the agreement of the simulation results with the measurements is good except for underprediction of the velocity fluctuation in recirculation region. Downstream of reattachment, there is some over-prediction of the streamwise velocity fluctuations with finer grid, however, the results obtained on the coarse grid show, in general, a lower value of streamwise velocity fluctuations, especially near the walls, which is probably caused by insufficient grid resolution and the physical dissipation of SGS model.

The resolved Reynolds stress at different streamwise sections are shown in Figure 7. Though the coarse-mesh result tends to over-predict the peak Reynolds stress, overall not much difference is seen in this quantity for different meshes. The position of the peak Reynolds stress in the normalwise is between and , which is similar to fully developed incompressible channel flow.

##### 4.3. The Turbulent Flow in a Francis Turbine Passage

###### 4.3.1. Computational Setup

The numerical example is a runner blade of a test Type-HLA551-LJ-43 Francis hydroturbine model. The computational domain including the distributor (stay vanes and guide vanes) domain and the runner domain is shown in Figure 8, and consists of one stay vane, one guide vane, and a runner blade. The distributor computational domain corresponds to an interblade channel that is bounded upstream by a cylindrical patch A-A and downstream by a conical patch B-B. The distributor inlet section corresponds to the spiral casing outlet section, while the outlet section is conventionally considered to be the distributor-runner interface. The runner computational domain also corresponds to an inter-blade channel that is bounded upstream by a conical patch (wrapped on the same conical surface as the distributor outlet), then across the runner middle axis C-C, and is extended downstream up to the draft tube inlet of radius D-D. For separation of the made-up turbulence that is caused by the strong 3D skewing passage, the geometry is normalized with half of the external diameter of the runner blade, βm. The mean chord dimension of the runner blade passage (streamwise dimension) is prescribed as and the mean blade-pitch dimension as , and the mean spanwise dimension of the blade passage is prescribed as . The mean coming flow velocity at the inlet of the runner domain (B-B) passage is defined as reference velocity βm/s. The flow Reynolds number is defined as , and the attack angle of guide vane relative the blade is 0Β°.

**(a) The configuration of a blade**

**(b) The whole configuration of the flow passage**

The meshes around the runner blade surfaces are first generated, and then the other domains are calculated in turn. The tetrahedral fluid element number of the whole model is 1,957,828 and 11,160 nodes are distributed over one side of the blade surfaces (averaging 90 points in the streamwise direction and 124 points in the spanwise direction).

The computational boundary conditions are shown in Figure 8, and are presented as follows. A uniform velocity field that is normal at the inflow section is imposed on the distributor inlet section, the turbulence quantities (the turbulence intensity is 6% and turbulence length scale is βm, where is the hydroradius of the distributor inlet section) are prescribed on the inflow section of the distributor, the free outflow condition is specified on the runner outlet (draft tube inlet), the periodic conditions are imposed on the pitchwise periodic boundary, and no-slip wall conditions are imposed on the stay vane, guide vane, runner blade, and distributor upper and lower rings as well as on the crown and band surfaces of the runner blade, respectively.

A strategy of changeable time interval is adopted in simulating progress for the sake of convergence and accuracy. The maximum time step is limited less than ( as the passing period of the blade passage). The grid turbulence simulation has progressed to 20β and collection of the sampling data for time statistics starts from 4β.

###### 4.3.2. Numerical Results

The time-averaged static pressure coefficients, , along the suctionand pressure sides of the blade are shown in Figure 9(a), where is the mean static pressure in flowing fileds, is the mean static pressure in the draft tube inlet section (D-D section in Figure 8). The time-averaged skin friction coefficients, , along the suction and pressure sides of the blade are shown in Figure 9(b), where is the wall shear stress. On the suction side of the blade, a short region of acceleration is seen at the leading zone of the passage, and after a short decelerated from to , a high pressure gradient is kept until the trailing edge of the blade. The flow direction may be changed due to reverse pressure gradient, which leads an intense negative pressure zones and a high wall-shear stress on the suction surface. At the same time, it can be seen that the pressure distributions are more regular on the pressure side, although largely different along the , , and lines on the suction side, which is demonstrated that the flows become more complex near the suction surface affected by the high 3D curved blade.

**(a)**

**(b)**

Figure 10 shows the static pressure distribution based on the time-averaged Root-Mean-Square (RMS) on the sections along the pitch-wise direction. is the pressure surface and is the suction surface, and and are formed by 5Β° rotation of and around the hydroturbine shaft, respectively. Near the pressure side ( and ), the pressure fluctuation is up to maximum near the leading of the blade passage, which is gradual decreasing along the streamline and up to minimum at the trailing edge of the blade passage. The contours of the RMS static pressure on the suction surfaces ( and ) show that local pressure gradient is larger than that on the pressure surface; moreover, the local low pressure region mainly lies in the suction side, which suggested that the cavitations more easily come into being on the suction surface.

Figure 11 shows the velocity magnitude distributions based on time-averaged RMS at the sections. The velocity fluctuation gradually increase along the streamline in the blade passage on the pressure side; however, it is not all true on the suction side. Some contour lines self-closed are exhibited on these blade sections, which is the reason why there is a strongly swirling flow structure in this zone. The remarkable difference of the velocity and pressure distributions based on the RSM between the suction and pressure sides also implies that the flow patterns are greatly affected by the distorted wakes from the upstream due to the attack angle of guide vane and the curvature of the blade wall itself.

Figure 12 shows the evolution of computational SGS kinetic energy of some classic points along the streamline on blade passage. It can be seen that the subgrid kinetic energy is large difference at different space points; moreover, the subgrid kinetic energy is also large fluctuation as time even though at the same space point, which is suggested that the difference of energy transport for different eddy scales at different time. Simultaneously, Figure 12 also shows the difference of SGS kinetic energy at different space points, which is verifies that the difference of eddy scales induced by skew blade.

Figure 13 shows the evolution of computational SGS stress of some classic points along the streamline in blade passage Numbered A, B, and C in Figure 8(a), where . It is shown that the peak of both the normal ones and the shear stress in point is higher than point and . It is also shown that the SGS stress distributions in space and time are large difference, which shows strong turbulent characters in complex blade passage.

(a) |

(b) |

(c) |

(d) |

(e) |

(f) |

Figure 14 shows the instantaneous isosurface of spanwise vorticity in the blade passage at different time. These figures clearly show the swirling flow structures. A strong clockwise spanwise vortex are viewed in the leading zone, then they are gradually elongated as time and the vortexes are then broken, a few vortex pairs rotating at clockwise and counterclockwise directions come into being along the band zone. Firstly, the big swirling vortex controls the main flow structure in the leading zone. With the collective spanwise vortex evolving continuously towards the downstream, the distances between vortex pairs became longer its intensity is gradually decreasing, and the flow becomes finally instability.

(a) |

(b) |

(c) |

(d) |

#### 5. Conclusion

A computational method is presented to apply DSGS model for the LES of unsteady turbulent flow problems in complex geometries. Based on the Taylor-Galerkin schemes for the convection-diffusion problems, this model is implemented in a three-dimensional finite element code using a three-step FEM. Qualitative and quantitative aspects of three-dimensional turbulent flow in a channel, turbulent flow over a backward-facing step and turbulent flow in a blade passage of a Francis turbine are simulated. The statistic characteristics of pressure, velocity, Reynolds stress and skin friction velocity, and so forth, are present. The spatiotemporally vortex structures are also shown in detail. The present numerical results are well in agreement with the DNS results or experimental data, which is verified that the three-step FEM for LES with DSGS model is validity for transient turbulent flow in complex geometries.

#### Acknowledgments

The authors thank the National Natural Science Foundation of China (NSFC) (Grants no. 11002063 and 50839003), the Natural Science Foundation of Yunnan Province of China (Grants no. 2008GA027 and 2009ZC035M), and the Kunming University of Science and Technology (Grant no. KKZ3200906003) for financial support of this research.