High Accurate Finite Volume Method for Large Eddy Simulation of Complex Turbulent Flows

High Accurate Finite Volume Method for Large Eddy Simulation of Complex Turbulent Flows

Lan Xu a,d , Guixiang Cui b[1], Chunxiao Xu b, Zhishi Wang c, Zhaoshun Zhang b and Naixiang Chen a

a Department of thermal energy engineering, Tsinghua university, Beijing, 100084, China

b Department of Engineering Mechanics, Tsinghua university, Beijing, 100084, China

c Department of Civil and Environmental Engineering, Faculty of Science and Technology,

University of Macau , Macao , China

d College of Science, Donghua University, Shanghai, 200051, China

______

Abstract

This paper proposes a finite volume method with compact fourth order accuracy scheme for large eddy simulation (LES). Two-dimensional lid-driven cavity flow and a flow over an oscillating plate are used as examples to verify both the accuracy and convenience of the proposed scheme. A turbulent channel flow and a turbulent flow over a backward facing step are numerically tested for its effectiveness by LES with dynamic Smagorinsky subgrid model. Immersed boundary method (IBM) is applied in this paper to deal with flows with complex configuration so that the boundary condition on the rigid wall can be satisfied well. A curved channel flow and a flow around an airfoil of NACA0012 are computed with immersed boundary method, and comparison with experimental data is also made, showing that the higher accurate finite volume method for LES is proved to be a promising numerical method.

Keywords: large eddy simulation, finite volume method, compact scheme, immersed boundary method

______

1. Introduction

Turbulence is the most difficult problem in the fluid mechanics, and it has not yet fully solved, it appears everywhere even in nano flows[1-3]. Large eddy simulation is believed to be a potential prediction method for complex turbulent flows in foreseeable future[4], however, a number of problems are still needed to be further resolved before it can be used in practice. On the physical aspect the subgrid model is the most significant issue in LES and great efforts have been made in the development of LES since it was proposed in 1970’s. A number of subgrid models are available, among which the dynamic Smagroinsky model is considered to be the best one, which will be used in this paper. Wall model for subgrid stress is another important issue in practical use of LES for complex wall-bounded turbulent flows [5]. A number of wall models have been proposed, e.g., equilibrium layer model, two-layer model and others. Although the wall model is not a fully resolved issue in LES, the equilibrium layer model based on the logarithm law in near-wall turbulence is easy to use in numerical calculation and will be used in the paper.

On contrast to the physical issues the numerical issues are more serious in the practical development of LES. Since LES is an unsteady solver it needs powerful computer for long time calculation to obtain reliable statistics. Both accuracy and efficiency of the numerical method are the major keys to practical use of LES. Comparing with other various numerical solutions, the finite volume method (FVM) [6,7] is more robust and has better conservative property of mass, momentum and energy transfer and also it is easier to accommodate the wall subgrid model. However the conventional FVM, e.g., SIMPLE with up-wind scheme, is in low accuracy and is not suitable for LES. As we know that the subgrid stress is proportional to the square of filter length D, which is in the same order of magnitude as the grid mesh length[8]. When a lower accuracy numerical scheme is used, the numerical errors may exceed the subgrid stress and one can not obtain a real LES. To keep the advantage of FVM we try to use higher order accuracy scheme in the interpolation between volume and surface averages of flow variables. The Padé finite volume compact method [9], which is successfully applied in numerical computation of two dimensional Navier-Stokes equation [10], is extended to LES of three dimensional unsteady flow. To effectively resolve wall-bounded turbulence, the non-uniform grids is needed in the wall region and we deduce an higher order accuracy FVM compact scheme on non-uniform grids. The proposed method is tested in a three-dimensional lid-driven cavity flow, a turbulent channel flow and a turbulent flow over a backward facing step with satisfaction. To deal with flows with complex configuration the immersed boundary method (IBM) is used to satisfy the boundary condition on the rigid wall. A curved channel flow and a flow around an airfoil of NACA0012 are computed successfully with immersed boundary method.

The paper is arranged as follows. In section 2 a high accurate scheme is formulated, and numerical verification is made; In section 3 the accuracy and effectives of the proposed scheme are numerically illustrated; Applications of the numerical method to the simulation of complex turbulent flows are presented in section 4. Conclusions are summarized in section 5.

2. Formulation of the numerical method

2.1. Governing equations of LES

The governing equations of large eddy simulation can be obtained by filtering Navier-Stokes equation which can be written as:

/ (1)
/ (2)

in which variables with upper bar stand for the filtered quantities or so-called resolved scale turbulent quantities and is the subgrid stress. We apply eddy viscosity type model in the paper:

/ (3)

in which and the eddy viscosity will be determined by dynamical procedure[11].

2.2. The spatial discretization of the equations

The conservative laws of mass and momentum for an element of incompressible fluid can be written as follows

/ (4)
/ (5)

where the Cartesian coordinates are used with x-axis in west-east direction, y-axis in south-north direction and z-axis in bottom-top direction. are the mass fluxes on the element surfaces , f stands for the velocity components and m is the fluid viscosity. The subscripts e, w, n, s, t, b in equation (5) denote the surface of the element and P is the center of the element as illustrated in Figure 1 and S is the external force terms in momentum equation, such as the pressure gradient.

(a) the control element (b) the grids

Figure 1 Illustration for the non-staggered grids used in computation

We use non-staggered grids in computation since it is easy in programming. To avoid the decoupling between velocity and pressure in non-staggered grids we use momentum interpolation with fourth order compact scheme. The accuracy of FVM depends on the discretization formula used in connection between primary variables at grid points W, E, P and fluxes on element surfaces e, w as shown in Figure 1 (b). We use the fourth order Padé type compact scheme proposed by Kobayashi (1999) as the discretization formula. The one dimensional formulae for interior points in x direction can be written as follows

/ (6)

with

in which is the average of f on the volume, i.e. the variable at center of the element, is the average of f on the element surfaces, i.e. the variables at center of surfaces e and w; are mesh lengths of the grids; tDx is a shift operator defined as follows

/ (7)

Similar formula can be derived for the boundary points as follows

/ (8)
/ (9)

with

When the formulas on uniform grids are as follows

/ (10)
/ (11)
/ (12)

The derivatives of variable f can also be discretized with fourth order accuracy, for instance in interior points

/ (13)

with

and on boundary

/ (14)
/ (15)

with

The formulas on uniform grids are as follows

/ (16)
/ (17)
/ (18)

For brevity all interpolation formulas can be written in matrix form as：

/ (19)
/ (20)

Inserting the compact interpolation into equation (4) and (5) we have the semi-discrete form of equation (2).

/ (21)

in which and b is a source term as in usual FVM.

2.3. Time advancement

Fourth order Runge-Kutta integration is used in time marching, the governing equation is written as

/ (22)

in which C and D are the convective and diffusive operators respectively and discretized by the compact scheme as stated above. Gp stands for the pressure gradient which is solved by coupling the momentum equation with continuity equation as follows.

/ (23)

in which M stands for the operator of continuity equation and , with

/ (24)
/ (25)
/ (26)
/ (27)

2.4. The solution of pressure

The pressure is calculated by continuity equation. To avoid the decoupling between velocity and pressure on non-staggered grids we use momentum interpolation with fourth order accuracy to determine the velocity at the surface of the FVM cell. For instance the momentum equation in x direction in the second step of Runge-Kutta integration is

/ (28)

in which ,, du is geometry parameter for FVM cell.

From equation (16) the momentum interpolation can be obtained as:

/ (29)

and velocity on the cell surface can be calculated by the following equation

/ 30

After the solution of cell surface velocity the pressure will be computed by strong implicit procedure (SIP) [12].

2.5. Immersed boundary method

The immersed boundary method (IBM) was originally introduced in the pioneering work of Peskin (1972) [13]. The basic idea of IBM lies in the definition of solid (either moving or stationary) boundaries. Traditionally most computational methods use complicated boundary fitted grids to define the geometrical configuration and the flow regions. The IBM actually simulates the presence of solid bodies by means of suitably defined body forces applied on the momentum equations. The N-S equations allow the specification of such a body force fi, inserted as a source term i.e.,

/ (31)

The body force is computed at every time step, so that the velocity field on an arbitrary surface is driven to a specified value. In general, that surface can move and does not necessarily have to coincide with a grid line. As shown by Mohd-Yusof (1997) [14], considering the time-discretized version of Eq. (31), one can write

/ (32)

where the term RHS contains all the pressure gradient, advection and diffusion terms. In order to drive the velocity field at the next time step , to the desired level , it is sufficient to formulate the source term f of the Eq. (32) as

/ (33)

which is imposed appropriately in the discretized form of the conservation equations through the boundary conditions of the the solid walls. In general the boundary of the region where we want is not coincide with a coordinate line. So the value of fi on the cell node closest to the surface but outside the body is linearly interpolated.

3. The check of numerical accuracy

3.1. Two-dimensional lid-driven cavity flow with gravity

The real accuracy of proposed method is checked by a cavity flow with gravity at Re=1. The flow is driven by upper cover with particular motion and an analytical solution can be obtained for the flow field inside cavity. The details of the flow case can be found in Pereira (2001)[10]. The accuracy of proposed method is presented in Table 1 and 2 with non-uniform and uniform grids respectively. The errors are the maximum deviations between the numerical solution and analytical results and the accuracy is estimated by the exponent of power law. It is clear that the proposed method is of fourth order accuracy approximately and error is smaller on non-uniform grids than that on uniform grids with same size.

Table 1

Errors and accuracy of proposed method with non-uniform grids

Grid points / Longitudinal velocity / Vertical velocity / Pressure
errors / accuracy / errors / accuracy / errors / accuracy
8×8
16×16
32×32
64×64 / 3.24E-3
2.59E-4
2.18E-5
1.60E-6 / —
3.54
3.45
3.70 / 7.32E-3
6.13E-4
4.37E-5
3.52E-6 / —
3.45
3.75
3.54 / 2.95E-2
2.51E-3
2.12E-4
1.63E-5 / —
3.44
3.44
3.60

Table 2

Errors and accuracy of proposed method with uniform grids

Grid points / Longitudinal velocity / Vertical velocity / Pressure
errors / accuracy / errors / accuracy / errors / accuracy
8×8
16×16
32×32
64×64 / 3.43E-3
2.48E-4
1.86E-5
1.43E-6 / —
3.79
3.74
3.70 / 6.71E-3
4.72E-4
3.46E-5
2.61E-6 / —
3.83
3.77
3.73 / 4.55E-2
3.43E-3
2.74E-4
2.21E-5 / —
3.73
3.65
3.63

3.2. Oscillating plate

A flow above an oscillating plate is calculated in order to check the accuracy of the method for unsteady flows. This coordinates are set up with the x-axis along the plate, and the y-axis normal to it. The velocity of the plate is prescribed as and the exact solution can be obtained for this problem. The errors between our computational results and exact solutions are listed in Table 3. And the accuracy is estimated by the exponent of power law. It is clear that the method has approximately fourth order accuracy.

Table 3

Errors and accuracy of proposed method with different numerical resolutions

Grid points / Longitudinal velocity / Vertical velocity / Pressure
errors / accuracy / errors / accuracy / errors / accuracy
25×25
50×50
100×100 / 8.90E-3
6.21E-4
4.43E-5 / —
3.84
3.81 / 7.41E-3
5.43 E-4
3.82E-5 / —
3.77
3.83 / 1.11E-2
1.03E-3
8.09E-5 / —
3.43
3.67

4. Applications

4.1. Case 1 Turbulent channel flow

The first testing case is a turbulent flow through channel; the geometry is illustrated in Figure 2. The Reynolds number is and the half width of channel H=1.

Figure 2 The geometry of turbulent flow in channel

The streamwise, normal and spanwise directions are , , respectively，and ,, are the velocity components in correspondent directions. Computational domain is composed of a rectangular box with the normal height equaling 2, streamwise length 4p and spanwise width 2p. The periodic conditions are posed in both streamwise and spanwise directions; and a wall model of power law (Werner, 1991) [15] is given at the channel wall such that