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

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


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.

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 u i and p are the filtered velocity components and pressure, respectively, υ is the molecular kinematic viscosity, and τ ij is the components of the SGS stress tensor, which is given by τ ij u i u j − u i u j u i u j u j u i u i u j .

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.τ ij is assumed as a nonlinear function of the strain rate, and it may be written as follows: where the filtered strain rate component S ij and the eddy viscosity ν t are given by

2.5
In 2.5 , C s 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 Δ ΔxΔyΔz 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 C s 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 ∂ u j ∂x j 0, 2.6 where • indicates the application of the second filtering process using a test filter, and T ij is component of a stress tensor, with its component given by Taking into account 2.4 and 2.5 , T ij may be expressed as and the SGS eddy viscosity υ sgs in the dynamic subgrid-scale model are defined by 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: where Lilly 22 solved the system of equations given by 2.12 using the least squares minimization, obtaining the following expression for C: 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 C x, y, z, t 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 C hom t is often used in the practical simulations.C hom t 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 C x, y, z, t , that is, where • xyz denotes space averaging over the entire domain.C hom t 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: where L is the characteristic length and u 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 ∂ u j ∂x j 0, 2.17 where Re Lu 0 /ν and Re sgs Lu 0 /ν sgs .

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.
u i * are the apparent velocities from which the velocities in the present time step can be derived as, The superscripts n, n 1/3, n 1/2 and n 1 on the variables u i and p represent these values of u i and p at time n, n 1/3, n 1/2 and n 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,

3.7
where We can gain the solution of u n 1 i only if the p n 1 is got.We implement the divergence to 3.4 and using the incompressible conditions ∂ u i n 1 /∂x i 0 at time n 1, 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 .

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 u τ are considered.The computational parameters for large-eddy simulations are summarized in Table  Figure 1 shows plots of the mean streamwise velocity for different meshes in plus units, u as a function of y , where u u 1 /u τ , and, y 1−|y/h| u τ /ν, • 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 subgridscale 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 u rms u 1 u 1 , v rms u 2 u 2 , and w rms u 3 u 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 u rms value and an little under-prediction of the peak v rms and w rms values.In whole, the RMS velocity is agreement with these DNS results of 24 well.

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, h .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 δ/h 1.1.The total domain size used for the computations was 22h × 5h × 3h, which included an entry length of 2h.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.
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.

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.Adams et al [29] x/h 〈C 〉 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.
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 C f 2τ w /ρU 2 c , where τ w is the shear stress at the wall and U c 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 x/h 2 and from reattachment to x/h 16, but poor agreement in both the recirculation zone and downstream of x/h 16.The reason for the poor agreement downstream of x/h 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 C f .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 x/h 3.2 and 5.9 , where subgridscale 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 u v 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 y/h 0.8 and y/h 1.0, which is similar to fully developed incompressible channel flow.

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 Mathematical Problems in Engineering 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 l 0.07R 0.003 m, where R 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 −3 T T L/U ref as the passing period of the blade passage .The grid turbulence simulation has progressed to 20 T and collection of the sampling data for time statistics starts from 4 T .

Numerical Results
The time-averaged static pressure coefficients, C p 2 p − p out / ρU 2 ref , along the suctionand pressure sides of the blade are shown in Figure 9  flowing fileds, p out is the mean static pressure in the draft tube inlet section D-D section in Figure 8 .The time-averaged skin friction coefficients, C f 2τ w /ρU 2 ref , along the suction and pressure sides of the blade are shown in Figure 9 b , where τ w 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 x /L 0.1 to x /L 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 S c , S m , and S b 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.
Figure 10 shows the static pressure distribution based on the time-averaged Root-Mean-Square RMS on the sections along the pitch-wise direction.S 1 is the pressure surface and S 4 is the suction surface, and S 2 and S 3 are formed by 5 • rotation of S 1 and S 4 around the hydroturbine shaft, respectively.Near the pressure side S 1 and S 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 S 3 and S 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.
Figure 11  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 x A /L 0.2, x B /L 0.5, x C /L 0.8 .It is shown that the peak of both the normal ones and the shear stress in point B is higher than point A and C. 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.
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.

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-

3 . 9 in
which N i , N j , and N k are the shape functions.Then the flow is analyzed by the following algorithm: 1 Determine u n 1/3 i with 3.5 . 2 Determine u n 1/2 i with 3.6 .3 Calculate p n 1 with 3.8 .4 Determine u n 1 i with 3.7 .

Figure 1 :
Figure 1: Profiles of the mean streamwise velocity in turbulent channel flow at Re 590.

Figure 2 :Figure 3 :
Figure 2: Profiles of the RMS velocity fluctuations in turbulent channel flow at Re 590.

Figure 4 :
Figure 4: Coefficient of friction along the lower wall.

Figure 7 :
Figure 7: Resolved Reynolds stress at different streamwise sections.

b
The whole configuration of the flow passage

Figure 8 :
Figure 8: The geometric configuration, computational domain, and boundary conditions.

1 X 1 XFigure 9 :
Figure 9: Averaged pressure a and skin friction coefficient distribution b on blade surface.

2 )Figure 12 :
Figure 12: Evolution of computational subgrid kinetic energy on two sides of blade.

Table 1 :
Simulation parameters for the three channel numerical simulations.