Finite Element Analysis of Turbulent Flows Using LES and Dynamic Subgrid-Scale Models in Complex Geometries

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:𝜕𝑢_𝑗𝜕𝑥_𝑗𝜕=0,(2.1)𝑢_𝑖𝜕𝑡+𝜕𝜕𝑥_𝑗𝑢_𝑖𝑢_𝑗=−1𝜌𝜕𝑝𝜕𝑥_𝑖+𝜐𝜕²𝑢_𝑖𝜕𝑥_𝑗𝜕𝑥_𝑗−𝜕𝜕𝑥_𝑗𝜏_𝑖𝑗,(2.2) 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 𝜏_𝑖𝑗=𝑢_𝑖𝑢_𝑗−𝑢_𝑖𝑢_𝑗+𝑢_𝑖𝑢^𝑗+𝑢_𝑗𝑢^𝑖+𝑢^𝑖𝑢^𝑗.(2.3) 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:𝜏_𝑖𝑗=−2𝜐_𝑡𝑆_𝑖𝑗,(2.4) where the filtered strain rate component 𝑆_𝑖𝑗 and the eddy viscosity 𝜈_𝑡 are given by𝑆_𝑖𝑗=12𝜕𝑢_𝑖𝜕𝑥_𝑗+𝜕𝑢_𝑗𝜕𝑥_𝑖;𝜐^𝑡=𝐶_𝑠Δ²2𝑆_𝑖𝑗𝑆_𝑖𝑗=𝐶_𝑠Δ²|||𝑆|||.(2.5) 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 Δ=(Δ𝑥Δ𝑦Δ𝑧)^1/3, 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𝜕𝑢_𝑗𝜕𝑥_𝑗𝜕=0,(2.6)𝑢_𝑖𝜕𝑡+𝜕𝜕𝑥_𝑗𝑢_𝑖𝑢_𝑗=−1𝜌𝜕𝑃𝜕𝑥_𝑖+𝜐𝜕²𝑢_𝑖𝜕𝑥_𝑗𝜕𝑥_𝑗−𝜕𝜕𝑥_𝑗𝑇_𝑖𝑗,(2.7) 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𝑇_𝑖𝑗=𝑢_𝑖𝑢_𝑗−𝑢_𝑖𝑢_𝑗−𝑢^𝑖𝑢^𝑗.(2.8) Taking into account (2.4) and (2.5), 𝑇_𝑖𝑗 may be expressed as 𝑇_𝑖𝑗=−2𝐶Δ²|||𝑆|||𝑆_𝑖𝑗,(2.9) and the SGS eddy viscosity 𝜐_sgs in the dynamic subgrid-scale model are defined by𝜐_sgs=2𝐶Δ²|||𝑆|||,(2.10) where ⟨Δ⟩=2Δ.

Using (2.3) and (2.8), it is obtained that:𝐿_𝑖𝑗=𝑇_𝑖𝑗−𝜏_𝑖𝑗=𝑢_𝑖𝑢_𝑗−𝑢_𝑖𝑢_𝑗.(2.11) From (2.4), (2.5), (2.9), and (2.11), the following system of equations is obtained:𝐿_𝑖𝑗=𝑇_𝑖𝑗−𝜏_𝑖𝑗=−2𝐶𝑀_𝑖𝑗,(2.12) where𝑀_𝑖𝑗=Δ²|||𝑆|||𝑆_𝑖𝑗−Δ²|||𝑆|||𝑆_𝑖𝑗.(2.13) Lilly [22] solved the system of equations given by (2.12) using the least squares minimization, obtaining the following expression for 𝐶:𝐶=−12𝐿^𝑖𝑗𝑀_𝑖𝑗𝑀_𝑖𝑗𝑀_𝑖𝑗.(2.14) 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 𝐶_hom(𝑡) is often used in the practical simulations. 𝐶_hom(𝑡) 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,2𝐶(𝑥,𝑦,𝑧,𝑡)Δ𝑘^1/2sgs𝑆_𝑖𝑗𝑆_𝑖𝑗_𝑥𝑦𝑧=2𝐶_hom(𝑡)Δ𝑘^1/2_sgs𝑆_𝑖𝑗𝑆_𝑖𝑗_𝑥𝑦𝑧,(2.15) where ⟨⋅⟩_𝑥𝑦𝑧 denotes space averaging over the entire domain. 𝐶_hom(𝑡) is used to define the SGS turbulent viscosity, 𝜐_sgs, in (2.10) and in the momentum equations.

The following dimensionless forms of the variables are used, namely:𝑥^∗𝑖=𝑥^𝑖𝐿,𝑢_𝑖^∗=𝑢_𝑖𝑢₀,𝑡^∗=𝑡𝑢⁰𝐿,𝑝_𝑖^∗=𝑝_𝑖𝜌𝑢₀₂,(2.16) 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)𝜕𝑢_𝑗𝜕𝑥_𝑗𝜕=0,(2.17)𝑢_𝑖𝜕𝑡+𝜕𝜕𝑥_𝑗𝑢_𝑖𝑢_𝑗=−1𝜌𝜕𝑃𝜕𝑥_𝑖+1Re𝜕²𝑢_𝑖𝜕𝑥_𝑗𝜕𝑥_𝑗+𝜕𝜕𝑥_𝑗Re1_sgs𝜕𝑢_𝑖𝜕𝑥_𝑗,(2.18) where Re=𝐿𝑢₀/𝜈 and Re_sgs=𝐿𝑢₀/𝜈_sgs.

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. 𝑢_𝑖^𝑛+1/3−𝑢_𝑖^𝑛Δ𝑡/3=−𝜕𝑢_𝑖^𝑛𝑢_𝑗^𝑛𝜕𝑥_𝑗−𝜕𝑝^𝑛𝜕𝑥_𝑖+1Re𝜕²𝑢_𝑖^𝑛𝜕𝑥_𝑗𝜕𝑥_𝑗+𝜕𝜕𝑥_𝑗Re1_sgs𝜕𝑢_𝑖^𝑛𝜕𝑥_𝑗.(3.1)

Step 2. 𝑢_𝑖^𝑛+1/2−𝑢_𝑖^𝑛Δ𝑡/2=−𝜕𝑢_𝑖^𝑛+1/3𝑢_𝑗^𝑛+1/3𝜕𝑥_𝑗−𝜕𝑝^𝑛𝜕𝑥_𝑖+1Re𝜕²𝑢_𝑖^𝑛+1/3𝜕𝑥_𝑗𝜕𝑥_𝑗+𝜕𝜕𝑥_𝑗Re1_sgs𝜕𝑢_𝑖^𝑛+1/3𝜕𝑥_𝑗.(3.2)

Step 3. 𝑢_𝑖^∗−𝑢_𝑖^𝑛Δ𝑡=−𝜕𝑢_𝑖^𝑛+1/2𝑢_𝑗^𝑛+1/2𝜕𝑥_𝑗+1Re𝜕²𝑢_𝑖^𝑛+1/2𝜕𝑥_𝑗𝜕𝑥_𝑗+𝜕𝜕𝑥_𝑗Re1_sgs𝜕𝑢_𝑖^𝑛+1/2𝜕𝑥_𝑗.(3.3)⟨𝑢_𝑖⟩^∗ are the apparent velocities from which the velocities in the present time step can be derived as, 𝑢_𝑖^𝑛+1=𝑢_𝑖^∗−Δ𝑡𝜕𝑝^𝑛+1𝜕𝑥_𝑖.(3.4) The superscripts 𝑛,𝑛+1/3,𝑛+1/2 and 𝑛+1 on the variables 𝑢_𝑖 and 𝑝 represent these values of 𝑢_𝑖 and 𝑝 at time 𝑛,𝑛+1/3,𝑛+1/2 and 𝑛+1, respectively.

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

For Step 1, 𝑀_𝑖𝑗𝑢_𝑗^𝑛+1/3−𝑢_𝑗^𝑛Δ𝑡/3=−𝐴^𝑛𝑖𝑗𝑢_𝑗^𝑛−𝐵_𝑖𝑗𝑝_𝑗^𝑛−1Re𝑆^𝑖𝑗𝑢_𝑗^𝑛+1Re𝑅_𝑖^𝑛−1Re_sgs𝑆_𝑖𝑗𝑢_𝑗^𝑛+1Re_sgs𝑅_𝑖^𝑛.(3.5)

For Step 2, 𝑀_𝑖𝑗𝑢_𝑗^𝑛+1/2−𝑢_𝑗^𝑛Δ𝑡/3=−𝐴^{𝑛+1/3𝑖𝑗}𝑢_𝑗^𝑛+1/3−𝐵_𝑖𝑗𝑝_𝑗^𝑛−1Re𝑆^𝑖𝑗𝑢_𝑗^𝑛+1/3+1Re𝑅_𝑖^𝑛+1/3−1Re_sgs𝑆_𝑖𝑗𝑢_𝑗^𝑛+1/3+1Re_sgs𝑅_𝑖^𝑛+1/3.(3.6)

For Step 3, 𝑀_𝑖𝑗𝑢^∗𝑗^𝑛+1/2−𝑢_𝑗^𝑛Δ𝑡/3=−𝐴^{𝑛+1/2𝑖𝑗}𝑢_𝑗^𝑛+1/2−1Re𝑆^𝑖𝑗𝑢_𝑗^𝑛+1/2+1Re𝑅_𝑖^𝑛+1/2−1Re_sgs𝑆_𝑖𝑗𝑢_𝑗^𝑛+1/2+1Re_sgs𝑅_𝑖^𝑛+1/2.(3.7)

Equation (3.4) can be discretized as𝑀_𝑖𝑗𝑢_𝑗^𝑛+1=𝑀_𝑖𝑗𝑢^∗𝑗−Δ𝑡𝐵_𝑖𝑗𝑝_𝑗^𝑛+1,(3.8) where 𝑀_𝑖𝑗=_Ω𝑁_𝑖𝑁_𝑗𝑑Ω,𝐴^𝑛𝑖𝑗=_Ω𝑁_𝑖𝑁^𝑘𝑢_𝑘^𝑛𝜕𝑁_𝑗𝜕𝑥_𝑘𝑑Ω,𝐵^𝑖𝑗=_Ω𝑁_𝑖𝜕𝑁_𝑗𝜕𝑥_𝑘𝑆𝑑Ω,_𝑖𝑗=_Ω𝜕𝑁_𝑖𝜕𝑥_𝑘𝜕𝑁_𝑗𝜕𝑥_𝑘𝑑Ω,𝑅_𝑖^𝑛=_𝜕Ω𝑁_𝑖𝜕𝑢^𝑛𝜕𝑛𝑑𝑆,𝑅_𝑖^𝑛+1/3=_𝜕Ω𝑁_𝑖𝜕𝑢^𝑛+1/3𝜕𝑛𝑑𝑆,𝑅_𝑖^𝑛+1/2=_𝜕Ω𝑁_𝑖𝜕𝑢^𝑛+1/2𝜕𝑛𝑑𝑆,(3.9) in which 𝑁_𝑖, 𝑁_𝑗, and 𝑁_𝑘 are the shape functions.

Then the flow is analyzed by the following algorithm: (1) Determine 𝑢^𝑖𝑛+1/3 with (3.5). (2) Determine 𝑢^𝑖𝑛+1/2 with (3.6). (3) Calculate 𝑝^𝑛+1 with (3.8). (4) Determine 𝑢^𝑖𝑛+1 with (3.7).

We can gain the solution of 𝑢^𝑖𝑛+1 only if the 𝑝^𝑛+1 is got. We implement the divergence to (3.4) and using the incompressible conditions 𝜕⟨𝑢_𝑖⟩^𝑛+1/𝜕𝑥_𝑖=0 at time 𝑛+1, the pressure Poisson equation is derived to correct the velocity equation as∇²𝑝^𝑛+1=1Δ𝑡𝜕𝑢_𝑖^∗𝜕𝑥_𝑖.(3.10) 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 (64×96×64), medium (128×193×128), and fine-mesh (192×257×192) 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.

Table 1 

Simulation parameters for the three channel numerical simulations.

Figure 1 shows plots of the mean streamwise velocity for different meshes in plus units, ⟨𝑢⟩⁺ as a function of 𝑦⁺, where ⟨𝑢⟩⁺=⟨𝑢₁⟩/𝑢_𝜏, and, 𝑦⁺=(1−|𝑦/ℎ|)𝑢_𝜏/𝜈, ⟨⋅⟩ 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].

Figure 1 

Profiles of the mean streamwise velocity in turbulent channel flow at Re=590.

The root mean squares (RMS) of velocities in our numerical model are plotted in Figure 2. It is defined as 𝑢_rms=⟨𝑢^1𝑢^1⟩, 𝑣_rms=⟨𝑢^2𝑢^2⟩, and 𝑤_rms=⟨𝑢^3𝑢^3⟩. 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 𝑢_rms value and an little under-prediction of the peak 𝑣_rms and 𝑤_rms values. In whole, the RMS velocity is agreement with these DNS results of [24] well.

Figure 2 

Profiles of the RMS velocity fluctuations in turbulent channel flow at Re=590.

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 𝛿/ℎ=1.1. The total domain size used for the computations was 22ℎ×5ℎ×3ℎ, which included an entry length of 2ℎ. A grid containing 144×96×96 and a finer grid containing 192×128×128 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.

Figure 3 

Computational geometry with computational grid.

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 144×96×96 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 192×128×128 nodes and the results obtained by Akselvoll and Moin [32] are within the estimated experimental error bounds, while the present results using 144×96×96 nodes are just outside these bounds.

Table 2 

Mean reattachment lengths.

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 𝐶_𝑓=2𝜏_𝑤/𝜌𝑈^2𝑐, 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 𝑥/ℎ=2 and from reattachment to 𝑥/ℎ=16, but poor agreement in both the recirculation zone and downstream of 𝑥/ℎ=16. The reason for the poor agreement downstream of 𝑥/ℎ=16 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 4 

Coefficient of friction along the lower wall.

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.

Figure 5 

Mean streamwise velocity profiles.

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.

Figure 6 

RMS of resolved streamwise velocity fluctuations.

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 𝑦/ℎ=0.8 and 𝑦/ℎ=1.0, which is similar to fully developed incompressible channel flow.

Figure 7 

Resolved Reynolds stress at different streamwise sections.

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, 𝑅_ref=0.225 m. The mean chord dimension 𝐿 of the runner blade passage (streamwise dimension) is prescribed as 𝐿/𝑅_ref=0.724 and the mean blade-pitch dimension as 𝑦/𝑅_ref=0.46, and the mean spanwise dimension of the blade passage is prescribed as 𝑧/𝑅_ref=2.06. The mean coming flow velocity at the inlet of the runner domain (B-B) passage is defined as reference velocity 𝑈_ref=0.68 m/s. The flow Reynolds number is defined as Re=𝑈_ref𝐿/𝜈^𝑓=1.1×10⁵, 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

Figure 8 

The geometric configuration, computational domain, and boundary conditions.

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 𝑙=0.07𝑅=0.003 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 2.0×10⁻³𝑇 (𝑇=𝐿/𝑈_ref 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, 𝐶_𝑝=2(𝑝−𝑝_out)/(𝜌𝑈²_ref), along the suctionand pressure sides of the blade are shown in Figure 9(a), where 𝑝 is the mean static pressure in flowing fileds, 𝑝_out is the mean static pressure in the draft tube inlet section (D-D section in Figure 8). The time-averaged skin friction coefficients, 𝐶_𝑓=2𝜏_𝑤/𝜌𝑈²_ref, 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 𝑥′/𝐿=0.1 to 𝑥′/𝐿=0.25, 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 9 

Averaged pressure (a) and skin friction coefficient distribution (b) on blade surface.

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 10 

Static pressure distributions based on RMS (unit: Pa).

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 11 

Velocity magnitude distribution based on RMS (unit: m/s).

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 12 

Evolution of computational subgrid kinetic energy on two sides of 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 𝑥^_A/𝐿=0.2,𝑥^_B/𝐿=0.5,𝑥^_C/𝐿=0.8). 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 13 

Evolution of computational SGS stress.

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) 𝑡/𝑇=18.48

(b) 𝑡/𝑇=18.89

(c) 𝑡/𝑇=19.3

(d) 𝑡/𝑇=19.71

Figure 14 

Contours of isosurface of spanwise vorticity in blade passage Blue: ≤−120 s−1, Red: ≥60 s−1.

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.