Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2011, Article ID 712372, 20 pages
http://dx.doi.org/10.1155/2011/712372
Research Article

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

1Department of Engineering Mechanics, Kunming University of Science and Technology, Yunnan, Kunming 650051, China
2School of Engineering, University of Aberdeen, Aberdeen AB24 3UE, UK

Received 24 January 2011; Revised 24 April 2011; Accepted 25 May 2011

Academic Editor: Sergio Preidikman

Copyright © 2011 Wang Wenquan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [68]. 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 [1215]. 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.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𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖;𝜐𝑡=𝐶𝑠Δ22𝑆𝑖𝑗𝑆𝑖𝑗=𝐶𝑠Δ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𝑢𝑖𝜕𝑥𝑗𝜕𝑥𝑗𝜕𝜕𝑥𝑗𝑇𝑖𝑗,(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|||𝑆|||𝑆𝑖𝑗,(2.9) and the SGS eddy viscosity 𝜐sgs in the dynamic subgrid-scale model are defined by𝜐sgs=2𝐶Δ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|||𝑆|||𝑆𝑖𝑗Δ2|||𝑆|||𝑆𝑖𝑗.(2.13) Lilly [22] solved the system of equations given by (2.12) using the least squares minimization, obtaining the following expression for 𝐶:1𝐶=2𝐿𝑖𝑗𝑀𝑖𝑗𝑀𝑖𝑗𝑀𝑖𝑗.(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/2sgs𝑆𝑖𝑗𝑆𝑖𝑗𝑥𝑦𝑧,(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:𝑥𝑖=𝑥𝑖𝐿,𝑢𝑖=𝑢𝑖𝑢0,𝑡=𝑡𝑢0𝐿,𝑝𝑖=𝑝𝑖𝜌𝑢02,(2.16) where 𝐿 is the characteristic length and 𝑢0 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=𝜌𝜕𝑃𝜕𝑥𝑖+1𝜕Re2𝑢𝑖𝜕𝑥𝑗𝜕𝑥𝑗+𝜕𝜕𝑥𝑗1Resgs𝜕𝑢𝑖𝜕𝑥𝑗,(2.18) where Re=𝐿𝑢0/𝜈 and Resgs=𝐿𝑢0/𝜈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=𝑢𝑖𝑛𝑢𝑗𝑛𝜕𝑥𝑗𝜕𝑝𝑛𝜕𝑥𝑖+1𝜕Re2𝑢𝑖𝑛𝜕𝑥𝑗𝜕𝑥𝑗+𝜕𝜕𝑥𝑗1Resgs𝜕𝑢𝑖𝑛𝜕𝑥𝑗.(3.1)

Step 2. 𝑢𝑖𝑛+1/2𝑢𝑖𝑛𝜕Δ𝑡/2=𝑢𝑖𝑛+1/3𝑢𝑗𝑛+1/3𝜕𝑥𝑗𝜕𝑝𝑛𝜕𝑥𝑖+1𝜕Re2𝑢𝑖𝑛+1/3𝜕𝑥𝑗𝜕𝑥𝑗+𝜕𝜕𝑥𝑗1Resgs𝜕𝑢𝑖𝑛+1/3𝜕𝑥𝑗.(3.2)

Step 3. 𝑢𝑖𝑢𝑖𝑛𝜕Δ𝑡=𝑢𝑖𝑛+1/2𝑢𝑗𝑛+1/2𝜕𝑥𝑗+1𝜕Re2𝑢𝑖𝑛+1/2𝜕𝑥𝑗𝜕𝑥𝑗+𝜕𝜕𝑥𝑗1Resgs𝜕𝑢𝑖𝑛+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=𝐴𝑛𝑖𝑗𝑢𝑗𝑛𝐵𝑖𝑗𝑝𝑗𝑛1𝑆Re𝑖𝑗𝑢𝑗𝑛+1Re𝑅𝑖𝑛1Resgs𝑆𝑖𝑗𝑢𝑗𝑛+1Resgs𝑅𝑖𝑛.(3.5)

For Step 2, 𝑀𝑖𝑗𝑢𝑗𝑛+1/2𝑢𝑗𝑛Δ𝑡/3=𝐴𝑛+1/3𝑖𝑗𝑢𝑗𝑛+1/3𝐵𝑖𝑗𝑝𝑗𝑛1𝑆Re𝑖𝑗𝑢𝑗𝑛+1/3+1Re𝑅𝑖𝑛+1/31Resgs𝑆𝑖𝑗𝑢𝑗𝑛+1/3+1Resgs𝑅𝑖𝑛+1/3.(3.6)

For Step 3, 𝑀𝑖𝑗𝑢𝑗𝑛+1/2𝑢𝑗𝑛Δ𝑡/3=𝐴𝑛+1/2𝑖𝑗𝑢𝑗𝑛+1/21𝑆Re𝑖𝑗𝑢𝑗𝑛+1/2+1Re𝑅𝑖𝑛+1/21Resgs𝑆𝑖𝑗𝑢𝑗𝑛+1/2+1Resgs𝑅𝑖𝑛+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 as2𝑝𝑛+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.

tab1
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 𝑢+=𝑢1/𝑢𝜏, 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].

712372.fig.001
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.

712372.fig.002
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.

712372.fig.003
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.

tab2
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.

712372.fig.004
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.

712372.fig.005
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.

712372.fig.006
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.

712372.fig.007
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×105, and the attack angle of guide vane relative the blade is 0°.

fig8
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×103𝑇 (𝑇=𝐿/𝑈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)/(𝜌𝑈2ref), 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𝜏𝑤/𝜌𝑈2ref, 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.

fig9
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. 𝑆1 is the pressure surface and 𝑆4 is the suction surface, and 𝑆2 and 𝑆3 are formed by 5° rotation of 𝑆1 and 𝑆4 around the hydroturbine shaft, respectively. Near the pressure side (𝑆1 and 𝑆2), 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 (𝑆3 and 𝑆4) 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.

712372.fig.0010
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.

712372.fig.0011
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.

712372.fig.0012
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.

fig13
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.

fig14
Figure 14: Contours of isosurface of spanwise vorticity in blade passage Blue: ≤−120 s1, Red: ≥60 s1.

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.

References

  1. C. M. Klaij, Space-time DG-FEM for complex flows, Ph.D. thesis, University of Twente Enschede, 2006. View at Zentralblatt MATH
  2. P. Moin, “Advances in large eddy simulation methodology for complex flows,” International Journal of Heat and Fluid Flow, vol. 24, no. 5, pp. 710–720, 2002. View at Publisher · View at Google Scholar · View at Scopus
  3. S. Conway, D. Caraeni, and L. Fuchs, “Large eddy simulation of the flow through the blades of a swirl generator,” International Journal of Heat and Fluid Flow, vol. 21, no. 5, pp. 664–673, 2000. View at Publisher · View at Google Scholar · View at Scopus
  4. L. X. Zhang, Y. Guo, and W. Q. Wang, “Large eddy simulation of turbulent flow in a true 3D Francis hydro turbine passage with dynamical fluid-structure interaction,” International Journal for Numerical Methods in Fluids, vol. 54, no. 5, pp. 517–541, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  5. P. Sagaut, Large Eddy Simulation for Incompressible Flows: An Introduction, Scientific Computation, Springer, New York, NY, USA, 2001.
  6. M. Tyagi and S. Acharya, “Large eddy simulation of turbulent flows in complex and moving rigid geometries using the immersed boundary method,” International Journal for Numerical Methods in Fluids, vol. 48, no. 7, pp. 691–722, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  7. W. Q. Wang, L. X. Zhang, Y. Yan, and Y. Guo, “Turbulent flow simulation using LES with dynamical mixed one-equation subgrid models in complex geometries,” International Journal for Numerical Methods in Fluids, vol. 63, no. 5, pp. 600–621, 2010. View at Publisher · View at Google Scholar · View at Scopus
  8. K. Mahesh, G. Constantinescu, and P. Moin, “A numerical method for large-eddy simulation in complex geometries,” Journal of Computational Physics, vol. 197, no. 1, pp. 215–240, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  9. J. N. Reddy and D. K. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press, London, UK, 2nd edition, 2001.
  10. T. L. Popiolek, A. M. Awruch, and P. R. F. Teixeira, “Finite element analysis of laminar and turbulent flows using LES and subgrid-scale models,” Applied Mathematical Modelling, vol. 30, no. 2, pp. 177–199, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  11. D. L. Young, T. I. Eldho, and J. T. Chang, “Large eddy simulation of turbulent flows in external flow field using three-step FEM-BEM model,” Engineering Analysis with Boundary Elements, vol. 30, no. 7, pp. 564–576, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  12. R. Codina, “Stabilized finite element approximation of transient incompressible flows using orthogonal subscales,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 39-40, pp. 4295–4321, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  13. P. A. B. Sampaio, P. H. Hallak, A. L. G. A. Coutinho, and M. S. Pfeil, “A stabilized finite element procedure for turbulent fluid-structure interaction using adaptive time-space refinement,” International Journal for Numerical Methods in Fluids, vol. 44, no. 6, pp. 673–693, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  14. A. E. Tejada-Martínez and K. E. Jansen, “On the interaction between dynamic model dissipation and numerical dissipation due to streamline upwind/Petrov-Galerkin stabilization,” Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 9–11, pp. 1225–1248, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  15. J. Hoffman and C. Johnson, “A new approach to computational turbulence modeling,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 23-24, pp. 2865–2880, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  16. A. N. Brooks and T. J. R. Hughes, “Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations,” Computer Methods in Applied Mechanics and Engineering, vol. 32, no. 1–3, pp. 199–259, 1982. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  17. R. Lohner, K. Morgan, and O. C. Zienkiewicz, “The solution of nonlinear hyperbolic equations systems by the finite element method,” International Journal for Numerical Methods in Fluids, vol. 4, no. 11, pp. 1043–1063, 1984. View at Publisher · View at Google Scholar · View at Scopus
  18. C. B. Jiang and M. Kawahara, “The analysis of unsteady incompressible flows by a three-step finite element method,” International Journal for Numerical Methods in Fluids, vol. 16, no. 9, pp. 793–811, 1993. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  19. A. Keating, U. Piomelli, and K. Bremhorst, “Large-eddy simulation of heat transfer downstream of a backward-facing step,” Journal of Turbulence, vol. 5, pp. 20–46, 2004. View at Publisher · View at Google Scholar · View at Scopus
  20. D. You and P. Moin, “A dynamic global-coefficient subgrid-scale model for large-eddy simulation of turbulent scalar transport in complex geometries,” Physics of Fluids, vol. 21, Article ID 045109, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  21. M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, “A dynamic subgrid-scale eddy viscosity model,” Physics of Fluids A, vol. 3, no. 7, pp. 1760–1765, 1991. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  22. D. K. Lilly, “A proposed modification of the Germano subgrid-scale closure method,” Physics of Fluids A, vol. 4, no. 3, pp. 633–635, 1992. View at Publisher · View at Google Scholar · View at Scopus
  23. H. A. Van Der Vorst, “BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of non-symmetric linear system,” SIAM Journal on Scientific and Statistical Computing, vol. 13, pp. 631–644, 1992. View at Publisher · View at Google Scholar
  24. R. D. Moser, J. Kim, and N. N. Mansour, “Direct numerical simulation of turbulent channel flow up to Res = 590,” Physics of Fluids, vol. 11, no. 4, pp. 943–945, 1999. View at Publisher · View at Google Scholar · View at Scopus
  25. Y. Bazilevs, V. M. Calo, J. A. Cotrell, T. J. R. Hughes, A. Reali, and G. Scovazzi, “Multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows,” Tech. Rep. 07-16, ICES, 2007. View at Google Scholar · View at Zentralblatt MATH
  26. Y. Morinishi and O. V. Vasilyev, “A recommended modification to the dynamic two-parameter mixed subgrid scale model for large eddy simulation of wall bounded turbulent flow,” Physics of Fluids, vol. 13, no. 11, pp. 3400–3410, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  27. J. C. Vogel and J. K. Eaton, “Combined heat transfer and fluid dynamic measurements downstream of a backward facing step,” Journal of Heat Transfer, vol. 107, no. 4, pp. 922–929, 1985. View at Publisher · View at Google Scholar · View at Scopus
  28. J. C. Vogel, Heat transfer and fluid mechanics measurements in the turbulent reattaching flow behind a backward-facing step, Ph.D. thesis, Stanford University, 1984.
  29. E. W. Adams, J. P. Johnston, and J. K. Eaton, “Experiments on the structure of turbulent reattaching flow,” Tech. Rep. MD-43, Thermosciences Division, Department of Mechanical Engineering, Stanford University, 1984. View at Google Scholar
  30. C. D. Pierce and P. Moin, “Method for generating equilibrium swirling inflow conditions,” The American Institute of Aeronautics and Astronautics Journal, vol. 36, no. 7, pp. 1325–1327, 1998. View at Publisher · View at Google Scholar · View at Scopus
  31. I. Orlanski, “A simple boundary condition for unbounded hyperbolic flows,” Journal of Computational Physics, vol. 21, no. 3, pp. 252–269, 1976. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  32. K. Akselvoll and P. Moin, “Large eddy simulation of turbulent confined coannular jets and turbulent flow over a backward facing step,” Tech. Rep. TF-63, Department of Mechanical Engineering, Stanford University, 1995. View at Google Scholar