Simulation of two-dimensional high reynolds number wake flows with the use of filtered lagrangian navier-stokes equations

TECHNICAL PAPERS

FLUID MECHANICS

^[e-mail], Universidade Federal do Paraná, Department of Environmental Engineering 80240-030 Curitiba, PR, Brazil

^[e-mail], Universidade Federal do Rio de Janeiro, Inst. Alberto. Luiz Coimbra 21945-970 Rio de Janeiro, RJ, Brazil

^[e-mail], Universidade Federal da Bahia, Instituto de Física 40170-280 Salvador, BA, Brazil

ABSTRACT

This work presents a method for numerical simulation of viscous fluid flows based on the use of filtered Lagrangian Navier-Stokes equations. It is based on the hypothesis that the small scale motions are responsible for the homogenization of the physical properties in infinitesimal fluid parcels through a mixing process, mathematically expressed as a spatial filter, which implicitly works as an averaging operator. When the filter was applied to the Navier-Stokes equations for incompressible flows, in the Lagrangian frame of reference, the functional form of the resulting equations was preserved, and the number of unknowns remained the same as the number of available equations, therefore circumventing the need of an explicit subgrid scale model. The method was validated by performing two-dimensional simulations of high Reynolds number wake flows. It was verified that the main large scale characteristics of the flows were properly replicated.

Keywords: filtered Lagrangian Navier-Stokes equations, wake flow, two-dimensional turbulence, semi-Lagrangian method, turbulent flows

Introduction

The Navier-Stokes (N-S) equations give a complete description of fluid flows. Although those equations are seemingly simple, there are no analytical solutions for them when applied to real flows. Therefore, it is necessary the use of numerical methods to obtain solutions for such cases.

Several different approaches have been developed to obtain numerical solutions for the N-S equations when applied to high Reynolds number (Re) fluid flows. Those based on the Reynolds-Averaged Navier-Stokes (RANS) equations have to deal with the closure problem, which occurs due to the non-linear character of the original N-S equations (Stull, 1995). Usually the averaged system uses the concept of eddy viscosity to model the terms that contain correlations of disturbance quantities (Wilcox, 2000). Another approach is the Direct Numerical Simulation (DNS), where all scales of motion of the fluid flow are represented (Moin and Mahesh, 1998). An intermediate approach is the Large-Eddy Simulation (LES), where, after application of a suitable filtering operator, the larger (filtered) scales of motion are simulated explicitly, whereas the smaller (residual) scales are modeled. The modeling of the smaller scales is also related to the closure problem since the filtering of the N-S equations produces subgrid terms that are parameterized, basically by using the eddy viscosity concept (Lesieur and Métais, 1996). Other approaches of the LES methodology use numerical methods to implicitly perform the filtering of the governing equations and the modeling of the residual stresses (Pope, 2000). One of these is the monotonically integrated LES (MILES) method (Fureby and Grinstein, 2002), where the implicit subgrid scale model is provided by a non-linear high-frequency filter built into a finite volume discretization.

Some turbulence models use the Lagrangian frame of reference for the development of closure schemes. The LANS-α model of Chen et al. (1999) employs mean values of two velocity vectors, calculated in the Lagrangian frame of reference, to model the effects of the small scale motions on the larger scales. In the transilient matrix method of Stull (1995), interactions of vortices of different spatial scales moving in the fluid are the basis for the development of the turbulent closure scheme. Probability Density Function methods (Pope, 2000), on the other hand, are stochastic, and use statistical properties of the smaller scale motions calculated in the Lagrangian frame of reference to develop the turbulence model. A characteristic shared by those three methods is that they seek solutions to the N-S equations written in the Eulerian frame of reference using closure schemes developed in the Lagrangian frame of reference.

This paper presents a method based on the use of filtered Lagrangian N-S equations for numerical simulation of incompressible turbulent fluid flows. The method is based on the physical interpretation of the small scale motions in fluid flows as being responsible for the homogenization of physical properties in the smaller scales, represented by an implicit spatial filter. Another feature of the method is its computational efficiency related to the use of the Lagrangian frame of reference, which provides high numerical stability.

In this paper, the problem will be treated in two dimensions only and it is not expected that all the features of real flows will be exactly reproduced. However, it is anticipated that the main large scale characteristics of the flows will be suitably reproduced by the numerical experiments.

In Section 2 the problem of two-dimensional high Reynolds number wake flow is presented. The main physical features of the problem, its analytical details, as well as experimental results are reviewed. That problem will be used to test the numerical method proposed in this study. Section 3 presents the mathematical model and introduces the filtered Lagrangian Navier-Stokes equations method. The procedure to obtain the numerical model with the use of the method is also described in that Section. The details of the implementation of the numerical model for the simulations of turbulent wake flows and the results of the numerical experiments are presented in Section 4. A discussion of the numerical method and its interpretation from a physical perspective is presented in Section 5, followed by the conclusions of this study, in Section 6.

Nomenclature

a = fluid particle displacement in horizontal direction, m

B _λ = filter function

b = fluid particle displacement in vertical direction, m

F = auxiliary Lagrangian notation

F_S = GCI safety factor

f(ξ) = normalized velocity deficit function

h = generic grid spacing, m

J = number of grid point in domain interior

j = grid point index

L = plate length, m

m = asymptotic order of discretization error

p = pressure, Pa

p_k = pressure boundary value, Pa

q₁ = value of a generic dependent variable on a fine grid

q₂ = value of a generic dependent variable on a coarse grid

q = generic dependent variable

= time averaged generic dependent variable

Re = Reynolds number

= wake Reynolds number

r = grid refinement ratio

t = time, s

= free stream horizontal velocity component, m/s

u = horizontal velocity component, m/s

u_i = generic velocity component, m/s

u_s = horizontal velocity deficit, m/s

v = vertical velocity component, m/s

W₀ = wake similarity parameter

x = horizontal coordinate, m

x_i = generic coordinate direction, m

y = vertical coordinate, m

y_c = ordinate of wake flow central axis, m

y_1/2 = wake boundary layer half-thickness, m

Greek Symbols

Δt = time step, s

Δx = horizontal grid space, m

Δy = vertical grid space, m

Δ₀ = wake similarity parameter

α = auxiliary integration variable

= GCI error estimate

= relative error estimate

ζ = auxiliary integration variable

η = enstrophy, 1/s²

= wake momentum thickness, m

λ = half filter width, m

v = kinematic viscosity, m²/s

ξ = dimensionless wake vertical coordinate

ρ = density, kg/m³

= vorticity, 1/s

Description of the Problem

A complete presentation of the two-dimensional turbulent wake problem can be found in Pope (2000). A summary of the main aspects of the problem is given.

Consider a stationary two-dimensional flow that is symmetric about the x-axis. A plane wake is formed when a uniform stream, with velocity in the x direction, flows along a solid body (wake generator).

The characteristic convective velocity is the free stream velocity , and the characteristic velocity deficit u_s is defined as

where is the ordinate of the flow central axis and is the time average of the horizontal velocity component u. The boundary layer half-thickness is the value of the ordinate where the local velocity deficit equals one half of the characteristic velocity deficit. It is defined by the following expression:

The expected behavior is for the wake thickness to increase and the velocity deficit to decrease along the x direction.

Since there are two characteristic velocities and varies in space the flow cannot be exactly self-similar. Nevertheless, it becomes asymptotically self-similar in the far-wake region, where .

We set a dimensionless coordinate in the direction normal to the main stream . The normalized velocity deficit function is defined as

The average longitudinal velocity component will be therefore

Using the two-dimensional stationary boundary layer equations and the definitions aforementioned, Pope (2000) presents the following solution for the universal dimensionless velocity deficit profile

Wygnanski, Champagne and Marasli (1986) carried out a systematic experimental study on two-dimensional, small deficit, turbulent wakes to verify the universality of the velocity deficit profiles in this type of flow. It was observed that the solution in Eq. (5) overestimates the average velocity at the outer limit of the wake. Based on experimental results the authors proposed the following expression for :

Wygnanski, Champagne and Marasli (1986) also defined the following similarity parameters:

where x₀ is the virtual origin and is the wake momentum thickness given by

In this work the virtual origin was arbitrarily set equal to zero, because it is two orders of magnitude smaller than the characteristic length scale of the wake region.

The wake Reynolds number is defined as

The product of and can be written as

where and are the normalizing velocity and length scales, respectively. Pope (2000) shows that the product is independent of x in the far-wake region, which is the region sufficiently away from the wake generator. Therefore, it is expected that W₀Δ₀ will be constant in that region. In the experiments performed by Wygnanski, Champagne and Marasli (1986), ranged between 3% and 15% in the far-wake, and it was observed that W₀Δ₀ was approximately constant in that region. In their experiments Re_è ranged from 1360 to 5900 along the wake.

The Filtered Lagrangian Navier-Stokes Equations Method

Mathematical model and filtering procedure

The mathematical model used in this problem is composed of the two-dimensional N-S equations for an incompressible Newtonian fluid flow, without external forces, in a Lagrangian frame of reference

and the two-dimensional continuity equation

where is the material derivative operator, u_i are the velocity components along the x_i directions, p is pressure, is the fluid density and v is the kinematic viscosity of the fluid.

The mathematical model uses the Navier-Stokes equations written in the Lagrangian frame of reference. To obtain the filtered equations we define the following auxiliary notation for two-dimensional flows

Coordinates (x,y) in Eqs. (13) and (14) refer to the position of an infinitesimal fluid parcel moving in the Lagrangian frame of reference, at time , with 2a and 2b being the fluid parcel displacements from the departure point at time to the final position at time , along the x and y directions respectively, and Δt being a discrete time interval.

For the implementation of the method, the total derivative of u_i in Eq. (11) will be discretized in the Lagrangian frame of reference with the use of the auxiliary notation of Eqs. (13) and (14):

In LES low-pass filtering of the N-S equations is explicitly performed in order to allow for the resulting fluid velocity field to be suitably resolved on a relatively coarse grid. The filtering operation is introduced after the velocity field is decomposed into resolved and subgrid components. The application of the filter on the non-linear advection terms of the decomposed Eulerian N-S equations yields additional unknowns, the stress tensors. The number of unknowns of the filtered system is greater than the number of equations. Therefore, these new unknowns have to be modeled to obtain the solution of the system of equations.

The filtered Lagrangian Navier-Stokes approach proposed in this study employs a spatial box filter of the form (Sorbjan, 1989)

where q is a dependent variable, x is length, B_λ = 1 is the filter function, 2λ is the filter width and ζ is an auxiliary variable.

We apply the box filter operator (16) to the differential-difference equation (15). Considering that ρ and v are constants, the conservation of a constant and linearity properties of the filter, and the fact that the operator is uniform, so that filtering and differentiation commute (Pope, 2000), then the filtered momentum equation is

Application of the filtering operator to the continuity Eq. (12) yields

Equations (17) and (18) are the filtered Lagrangian equations. Notice that the variables were not decomposed into resolved and subgrid components before the application of the filtering operator. The filtered equations refer only to the resolved scales of motion and they have the same functional form of the original N-S equations, except for the filtering operator. It is important to notice that the resulting system of equations has the same number of unknowns as the number of equations. Therefore, they make a closed system. Differently from LES, there are no additional unknowns in the filtered system. Consequently, there is no need to include additional closure equations. That is possible because there are no explicit non-linear terms in the original equations written in the Lagrangian frame of reference. In that frame of reference non-linear effects are implicit in the material derivative of the velocity components. Now it remains to define a procedure to implement the filtering operator.

Fureby and Grinstein (2002) claim that finite differences, volume and element methods filter the N-S equations over cells of the discrete computational domain. We use this idea to analyze the properties of the box filter (16). According to Sorbjan (1989), the box filter works as a spatial averaging operator. Wilcox (2000) also states that the values of the properties of the flow at grid points in a numerical simulation represent average values. He analyzed the approximation by centered finite differences for the first derivative of a continuous variable q(x) in a grid with spacing h between the grid points, and concluded that

It can be seen that the term between brackets has the same functional form of the box filter in Eq. (16).

The expression above states that the centered finite difference approximation can be interpreted as an operator that produces the spatial derivative of the mean value of q(x). The averaging implicit in the finite difference approximation also filters out almost all scales smaller than 2h (Sorbjan, 1989). Based on this property, it will not be necessary the actual calculation of the filtering operator by using Eq. (16). Instead, the spatial derivatives will be approximated by finite differences, and it will be assumed that the filtering operation implicit in the discretization will yield the mean value of the variables in discrete fluid parcels represented by the grid cells. This procedure is analogous to the implicit filtering of the finite volume discretization of the MILES method (Fureby and Grinstein, 2002).

One of the fundamental assumptions of this method is that the flow must be represented in the Lagrangian frame of reference. Consequently, the numerical solution of the filtered N-S equations requires the use of that frame of reference for consistency. To abide by this requirement we will employ the semi-Lagrangian method.

The semi-Lagrangian method was first proposed by Sawyer (1963) to numerically solve weather forecasting problems. Robert (1981) introduced a variation into Sawyer's semi-Lagrangian method where in each time step the future positions of the fluid particles are coincident with the grid points. In this variation, the departure point of the fluid particles is determined by moving the particle along a backward-in-time trajectory. With this change it became possible to have the advantage of a regular spatial representation of the domain, characteristic of Eulerian methods, as well as high computational efficiency due to the numerical stability of Lagrangian methods. A comprehensive description of the semi-Lagrangian method can be found in Durran (1999).

As stated before, an important advantage of the semi-Lagrangian method is the possibility of representation of the fluid properties at the same points of a fixed grid at all times. Therefore, it is possible to obtain time series of a variable at a grid point . Consequently, time averages of that variable can also be obtained as if the Eulerian frame of reference had been employed.

Numerical model

The semi-Lagrangian three-time-level scheme associated to the semi-implicit algorithm proposed by Robert et al. (1985) will be used in this work. As mentioned before, the problems will be restricted to two dimensions. We start from the filtered Lagrangian N-S equations (17) and continuity equation (18) in component form in two dimensions (x,y), where the symbol is omitted for simplicity

In the auxiliary notation of Eqs. (13) and (14), the position of a fluid parcel at time (t+Δt) will be coincident with a grid point, and displacements 2a and 2b from the departure point along the x and y directions, respectively, will be calculated iteratively (Robert, 1981). In general, the departure point of a fluid particle will not be coincident with a grid point position. Therefore, it will be necessary to use spatial interpolations to obtain the values of the variables at time .

Equations (20) and (21) are discretized implicitly for the pressure terms and explicitly for the viscous terms. Using the notation of Eqs. (13) and (14) we obtain

The discretization of pressure and diffusion terms of Eq. (17), which yields Eqs. (23) and (24), is valid due to the properties of the filter and the linearity of those terms. The same applies for the continuity equation (18).

The continuity equations for times and are

After some manipulation, Eqs. (23) and (24) can be rewritten as

Next, Eqs. (27) and (28) are differentiated with respect to x and y, respectively. The resulting equations are added and, with the use of Eqs. (25) and (26), the following elliptic equation for pressure is obtained

Equations (27)-(29) constitute the numerical model that will be employed for the prediction of u, v and p. Again, notice that these equations form a closed system, and no additional equations besides the filtered ones (20, 21 and 22) were needed.

Implementation of the Numerical Model

Computational domain geometry and boundary conditions

The computational domain simulates a smooth flat plate (the wake generator body) immersed in an initially uniform flow at zero incidence angle. The plate is placed at the flow centerline with respect to the y axis. The grid is uniform in the x and y directions with grid spaces and respectively, with . There is a short region upstream of the plate leading edge to allow for the adjustment of the flow to the presence of the plate. The inflow boundary therefore is not coincident with the plate leading edge. The wake region extends downstream of the plate trailing edge.

Along the inflow boundary Those values were kept constant during the simulation. In the upper and lower boundaries The vertical velocity component was adjusted along the simulation using the gradient boundary condition of Mason and Sykes (1979) to control the reflection of waves at the boundary back to the domain interior. Along the outflow boundary u was also computed with the use of the gradient boundary condition, and v was arbitrarily set equal to zero. The pressure along the external boundary was set to an arbitrary reference value p_k = 1.0 Pa. At the external boundaries the normal components of the Laplacians were arbitrarily set equal to zero and the along-boundary components were approximated by centered finite differences.

Discretization in space will be made in a collocated grid where u, v and p are calculated at the same grid point. Time discretization of the total derivative and pressure terms will use the three-time-level semi-Lagrangian scheme aforementioned, and the Euler scheme will be employed for the viscous terms. Spatial derivatives will be approximated by centered finite differences in the domain interior and by one-sided finite differences at the boundaries. Spatial interpolations in the domain interior will be local, with the use of cubic polynomials in the x and y directions. At the external boundaries interpolations will be cubic along the boundary and quadratic in the direction normal to the boundary.

The simulation configuration parameters were set in such a way that a long region was left downstream of the plate in order to obtain sufficiently large, to characterize a high Re wake flow comparable to the experiments of Wygnanski, Champagne and Marasli (1986).

Numerical simulations

The objectives of the numerical experiments were the following:

- to evaluate the ability of the method to replicate the main features of turbulent wake flows;

- to study the effects of resolution on the simulation results;

- to verify the ability of the method to simulate some universality properties of turbulent wake flows under different flow conditions;

- to assess the sensitivity of the method to different physical characteristics of the fluid represented in the simulations; and

- to verify possible effects of numerical diffusion, caused by spatial interpolations on the representation of physical viscous effects of the flow.

In order to verify the numerical stability of the simulations, two global parameters were used: the domain averaged kinetic energy (per unit mass) KE and the domain averaged entropy η, defined by

where u_j and v_j are the instantaneous velocity components at each grid point (j) of the interior of the computational domain, Ω_j is the vorticity at each grid point of the interior of the computational domain and J is the number of grid points in the interior of the computational domain. These variables were calculated at each time step, and they were chosen because, in the case of occurrence of numerical instability, they present significant variations during the simulation. As stated before, one of the advantages of the semi-Lagrangian method is that the values of the variables are referred to the same positions in space at each time step, since the arrival positions of the fluid particles are always coincident with grid points. For that reason, it is possible to obtain time or space averages of the variables referring to the same grid points (i.e. the same computational domain) along the time of simulation. It is important to notice that the velocity components used to obtain KE and η are the filtered (or resolved) ones, since those are the variables present in the numerical model, expressed by Eqs. (27)-(29).

Seven numerical simulations of high Re wake flows downstream of a flat plate were configured with the parameters presented in Table 1. For experiments 1, 2, 3 and 4 the values of ρ and v characterize air at 293 K, and for experiments 5 and 6 the fluid simulated is water also at 293 K. The domain dimensions were 25.0 m in the x direction and 0.2 m in the y direction. Experiments 1 and 2 had the same configuration parameters, except for the grid spacing, where in experiment 2 (high resolution) they were respectively half the size of those used in experiment 1 (low resolution). The same applies for experiments 3/4 (low/high resolution), and 5/6 (low/high resolution), respectively.

In all cases, the wake Reynolds number Re_è had values comparable to those of the experiments of Wygnanski, Champagne and Marasli (1986) (Table 1), characteristic of turbulent wakes. The simulated flow time in all experiments was 160 s.

It is important to point out that the filtering operator (16), which is implicitly applied through the finite differencing discretization in space, produces some smoothing of the variables in time too, since small scale fluctuations are averaged out in the limit of the grid cell. Nevertheless, the resolved part of the motion still varies in time. Figure 1 shows the time evolution of the u velocity component between 45 s and 60 s, at position x = 7.75 m, at the wake centerline ordinate, for experiment 1. It can be seen the presence of high frequency disturbances in the flow.

Figure 2 shows the temporal evolution of the domain averages of enstrophy and kinetic energy per unit mass for experiment 1. Both variables had sharp variations at the beginning of the simulation, followed by a quasi-stationary regime. It can be seen that the simulation remained numerically stable. In all experiments, the domain averages of kinetic energy and enstrophy had similar behavior. For that reason, only the results for experiment 1 are presented.

In order to analyze the behavior of the large scale features of the wake, grid point time averages of the flow variables were computed with data recorded in the time interval of the quasi-stationary regime for each experiment. Those grid point time averages of the variables are identified by an overbar (e.g. means the grid point time average of u).

Figure 3 presents the fields for experiments 1 through 6, from the plate trailing edge to the outflow boundary. The vertical representation was limited to the region where significant variations in were observed (about 30% of the domain height). The most external isotachs define the positions where = 99% , and they can be used to indicate the outer limits of the wakes. As expected, the thicknesses of the wakes increase and the velocity deficit u_s decreases in the x direction, and the wakes are symmetric about the flow central axes.

The corresponding low and high resolution simulations produced wakes with similar widths. However, the regions near the central axes presented some differences in the velocity fields. In all cases, the high resolution simulations reduced the velocity gradient along the x direction. That is, there was a more effective diffusion effect of the longitudinal velocity component along the x direction. Experiments 3/4, which represent higher turbulence values than experiments 1/2, presented more significant differences between the high and low resolution simulations. The differences are also significant in experiments 5/6, where the simulated fluid was water.

In section 2 it was seem that in the far wake region the product W₀Δ₀ is expected to be constant. Longitudinal profiles of W₀Δ₀ obtained in experiments 1 through 6 are presented in Fig. 4. In all cases, at the beginning of the wake regions there is a sharp decrease in W₀Δ₀. Qualitatively, from the figure, we can see that for experiments 1/2, downstream of x = 7.5 m approximately, W₀Δ₀ varies slowly in the x direction and asymptotically tends to a constant value. That is the expected behavior in the far wake region.

The same occurs for experiments 3/4, except that the asymptotic trend in W₀Δ₀ appears after x = 12.5 m, approximately. On the other hand, for experiments 5/6, W₀Δ₀ starts tending to a constant value around x = 5.0 m. In all experiments the higher resolution simulations produced smoother transitions in the value of W₀Δ₀ along the x direction than the low resolution simulations, but, in general, the differences are not too significant.

Four positions in the final part of the wake were chosen to verify if the simulations would be able to suitably replicate universal characteristics of the turbulent far wake, Those positions were located at positions corresponding to 78% (x = 19.40 m), 82% (x = 20.61 m), 87% (x = 21.83 m) and 92% (x = 23.04 m) of the domain length. They were chosen because they are in the region where parameter W₀Δ₀ is almost constant, which characterizes the far wake region. That is the region where we expect both self-similar and universal behavior of the wake velocity deficit profile f(ξ).

Figure 5 shows the simulated profiles of f(ξ) at the four aforementioned positions, together with the empirical profile proposed by Wygnanski, Champagne and Marasli (1986). By and large there was a good superposition of the profiles, indicating that the simulations were able to reproduce the self-similar characteristics of f(ξ) in the far wake, for different flow regimes and different fluids. The profiles also closely matched the empirical function of Eq. (6). The high resolution simulations produced better profiles near the wake edges than the low resolution ones. In the inner part of the wake all simulations replicated the velocity deficit profiles equally well.

The method based on the use of filtered Lagrangian N-S equations is implemented with the use of an implicit box filter, through a finite difference spatial discretization, which causes the filtering of the fluid properties in the smaller scales. Consequently, it is expected that short wavelength (high frequency) features will be filtered out. It is interesting to verify this filtering effect on the simulated flow.

Figure 6 shows energy spectra of the u and v velocity components, respectively, at positions corresponding to 75% of the wake length (x = 19.3 m for experiments 1/2, x = 19.7 m for experiments 3/4 and x = 19.2 m for experiments 5/6), at the ordinate of the flow central axis. According to Kraichnan (1967), the two-dimensional turbulence spectrum can present two inertial subranges: one in the large scales, where there is inverse energy cascade, and other in the small scales, where there is direct vorticity cascade. The graph of those regions would appear with slopes -5/3 and -3, respectively, in log-log representation. In Fig. 6, in the high frequency regions of the spectra, the graphs for experiments 1/2 and 3/4 are fairly similar, and they show slopes close to -3, spanning ranges of approximately one-half to three-quarters of a decade. For each case, the high frequency parts of the spectra are similar. In the intermediate range, with frequencies just below the inertial subrange, however, the graphs present some differences for the high and low resolution simulations, with the high resolution ones producing greater values of energy for those frequencies. The same occurs in experiments 5/6, but in those simulations the resolution effects were more significant in the spectrum of the u velocity component. Nevertheless, all experiments produced some inertial subrange, which is better defined for the v velocity component when water was the simulated fluid. Hence, despite the filtering effect of the finite difference scheme, and the relatively coarse resolution of the grids, the simulations were able to reproduce processes related to the production, decay and transference of two-dimensional turbulence vorticity in the smaller scales.

It is known that spatial interpolations in the semi-Lagrangian method introduce dissipation into a numerical simulation. It can be therefore questioned if the turbulent mixing and viscous dissipation processes reproduced by the simulations using the method presented in this paper are mostly due to numerical dissipation effects introduced by the interpolations. An experiment was made to clarify this point.

Simulation 7 was run for the exact same conditions of experiment 1, but in this case the viscosity term of the N-S equations was "switched-off" at time t = 83 s, by setting the kinematic viscosity (actually 10^-30), all other conditions remaining unchanged. This new flow would be therefore virtually inviscid. If the simulated small scale effects (turbulent mixing and viscous dissipation) are mostly due to numerical dissipation of the spatial interpolations, it should not be expected a significant change in the characteristics of the simulated flow. On the other hand, in case a remarkable change occurs in the conditions of the numerical simulation, it will demonstrate that the turbulent and viscous processes simulated by the method are more likely to be related to actual physical responses to the viscosity term in the filtered Lagrangian N-S equations, not to numerical diffusion introduced by the interpolation procedure.

Figure 7 presents the temporal evolutions of the domain averages of enstrophy and kinetic energy, and of the u velocity component at the first grid point downstream of the flat plate trailing edge at the ordinate of the flow central axis (that grid point would be in the "shadow" of the plate with respect to the flow) for simulation 7. It is apparent that as v is set H" 0 there was an immediate response in the simulation, with increase both in the kinetic energy and the enstrophy, and reduction of u to zero immediately downstream of the plate. Those effects are expected since there was neither viscous nor turbulent drag acting on the flow. Hence, the results of simulation 7 strongly suggest that the small scale turbulent effects reproduced by the simulation are actual physical processes related to the viscous term in the filtered Lagrangian N-S equations and not produced by numerical dissipation due to spatial interpolation of the semi-Lagrangian method.

Errors Estimation

Grid convergence analysis was performed by using the GCI discretization error estimator of Roache (1994). In GCI, estimate of the discretization error of variable q produced by the numerical solution on a fine grid ε is obtained by the following expression:

where q₁ and q₂ are the solutions obtained on the fine and coarse grids, respectively, F_S is a safety factor of 3 for general applications, r is the refinement ratio between the fine and coarse grids, and m is the asymptotic order of the discretization error.

Relative errors between the fine grid and coarse grid numerical solutions were also estimated, with the use of the following expression:

In this study, the grid refinement ratio between corresponding fine grid experiments (2, 4 and 6) and coarse grid experiments (1, 3 and 5) was r = 2. Centered finite differences were employed in the domain interior therefore m = 2. The variables selected to analyze the discretization errors were grid point values of at 11 positions, five along the wake centerline, and three on each side of the wake center, symmetrically positioned about the wake centerline. Values of were chosen because it is a critical variable for the normalization of all parameters that characterize the turbulent wake.

Tables 2 through 4 present the values of at the selected points and the corresponding GCI error estimates and relative errors for experiments 1/2, 3/4 and 5/6, respectively.

From the results presented in the tables, we can see that larger errors occurred in the centerline of the wake (y = 0.1 m), what is expected, since that is the region where the largest disturbances occur in the flow. Also, the errors decrease in the centerline as we move to the far wake region (x > 18.0 m). Therefore, we can expect that the self-similar and universal characteristics of the wake flow will be appropriately represented, since the errors are relatively small in that region. We can also see that at the points on the sides of the wake center the errors are approximately symmetric about the center, that they are remarkably smaller than those at the wake center, and also that they slightly increase as we move towards the far wake region.

Also, comparing the data of Tables 2 and 3 it can be seen that the differences between the values obtained at corresponding grid points for experiments 1/2 and 3/4 (both simulating air flows) are somewhat larger than the errors of each pair of simulations. Therefore, the numerical model was able to produce conditions that were numerically different in each pair of experiments, reflecting the different turbulent levels (i.e. Re numbers) of each flow. In other words, if the errors were larger than the differences between each pair of flow conditions we would not be sure that the numerical model was able to reproduce different physical flow conditions. We can say the same when we compare the results of Tables 2 and ⁴ as well as 3 and ⁴. The analysis shows that characteristics of different fluids (air for experiments 1 through 4 and water for experiments 5 and 6) were able to be captured by the simulations, although there is a significant overlap in the Reynolds numbers ranges of the different experiments.

Errors of the normalized velocity profiles were also analyzed. For each experiment, root mean square errors (RMS) of the normalized velocity deficit profiles simulated at each of the four selected x positions in the far wake region (Fig. 5) were calculated, using the semi-empirical values obtained by Eq. (6) as reference. Table 5 summarizes these results.

Discussion

The results of the experiments show that it is possible to simulate several aspects of turbulent wake flows with the use of the filtered Lagrangian N-S equations. Now it is necessary to propose a physical interpretation for the method.

According to the transilient turbulence theory for the atmospheric boundary layer developed by Stull (1995), the larger eddies are responsible for the transport of fluid along finite distances whereas the smaller eddies mix fluid properties in small scales. Starting from this idea the method presented in this paper assumes that the mixing in small scales causes homogenization of the fluid properties in those scales. Still focusing on small scales, one can suppose that the fluid particles are continuously interacting and exchanging properties among them. Then, it can be considered that infinitesimal fluid parcels will have their physical properties homogenized after each instantaneous interaction. The outcome of this homogenization process is that the physical properties in an infinitesimal fluid parcel will be averaged out, with the resulting averages being determined from the values the properties had in each fluid particle before the interactions. Notice that when the mean value is used to represent the result of the interaction of several fluid particles, all those particles that participated in the interaction process can be replaced by only one. That is possible because a certain physical property will have the same (mean) value in all fluid particles due to the mixing and homogenization. This substitution reduces the number of degrees of freedom of the system that represents the fluid flow (Stevens and Lenschow, 2001).

In incompressible Newtonian fluid flows the advection terms in the momentum equations are non-linear. Those terms appear in the equations when they are written in the Eulerian frame of reference. On the other hand, when the equations are written in the Lagrangian frame of reference the non-linearities are not explicit, because they are (implicitly) included in the total derivative.

The physical interpretation of the averaging in the Lagrangian frame of reference is that it represents the homogenization of properties of a fluid parcel due to interaction (mixing) of its several constituent particles. After the interactions, a certain physical property in the fluid parcel will assume a mean value computed from the value the property had in each of the individual particles before the interactions.

Notice that in the conceptual interpretation of the method it was neither necessary to invoke the ergodic hypothesis, nor to assume the existence of a spectral gap, as is usually required when applying the Reynolds averaging (Stull, 1995; Wilcox, 2000). Another important point is that during the averaging process we refer to a fluid parcel defined by a material surface. That assumption guarantees that we will be always referring to the same set of fluid particles during the infinitesimal time interval independently of any deformation process undergone by the infinitesimal fluid parcel.

Due to the assumption of homogenization of the fluid properties, the existence of homogeneity and isotropy after the particles interaction is implicit. However, these conditions will be limited to the characteristic length scale of the filter.

As stated before, Eqs. (27)-(29) constitute a closed numerical model for the prediction of u, v and p. Therefore, since in the resulting filtered system the number of unknowns is equal to the number of equations available, there is no need to employ an explicit turbulence model. Any non-linear processes in the flow, such as those related to turbulent effects, will not be explicitly represented in the filtered system, since in the Lagrangian frame of reference non-linear effects are implicitly represented in the material derivative.

Conclusion

This work presented a method based on filtered Lagrangian N-S equations for numerical simulation of boundary layer flows. The method is based on the hypothesis that the small scale motions are responsible for the homogenization of physical properties in infinitesimal fluid parcels through a mixing process. The mixing is mathematically represented by a filter that works as an averaging operator that, when applied to a fluid physical property following a streamline, yields the Lagrangian mean of that property. When the operator was applied to the Navier-Stokes equations for an incompressible flow written in the Lagrangian frame of reference, the functional form of the equations did not change, and the number of unknowns remained the same as the number of equations available.

The study focused on two-dimensional flows only, and its objective was to verify if the main characteristics of a high Re wake formed downstream of a solid body could be suitably replicated by the numerical simulations using the filtered Lagrangian N-S equations. The study on turbulent wakes of Wygnanski, Champagne and Marasli (1986) was used as reference to assess the accuracy of the numerical simulations. The results obtained were satisfactory since the simulated profiles closely matched the universal reference profile, and they clearly showed to be self-similar.

The energy spectra of the u and v components at selected grid points were also examined. The presence of a direct vorticity cascade along the flow central axis, expected in two-dimensional turbulent flows, was observed.

Also, it has been demonstrated in the study that the method allows for the replication of actual physical responses to the viscous term in the N-S equations, and that the simulated turbulent effects are not just a fortuitous product of numerical dissipation.

It is important to point out the numerical stability of the method. The use of the semi-Lagrangian scheme associated to the semi-implicit algorithm of Robert et al. (1985) provided unconditional linear stability and prevented the occurrence of aliasing.

Based on the results stated above the conclusion is that, for the flow problem focused on in this study, the underlying hypothesis of the method based on the use of the filtered Lagrangian N-S equations has been demonstrated to be correct. Although more work is needed to validate the method under other flow conditions, the good performance shown by the approach indicates that it may be competitive with other numerical methods presently applied to simulating boundary layer flows.

Paper accepted August, 2010.

Technical Editor: Francisco Ricardo Cunha (Footnotes)

• Aris, R., 1962 "Vectors, Tensors, and the Basic Equations of Fluid Mechanics", Prentice-Hall Inc., Englewood Cliffs.
• Chen, S., Holm, D.D., Margolin, L.G. and Zhang R., 1999, "Direct numerical simulation of the Navier-Stokes alpha model", Physica D, Vol. 133, pp. 66-83.
• Durran, D.R., 1999, "Numerical Methods for Wave Equations in Geophysical Fluid Dynamics", Springer-Verlag, New York.
• Fureby, C. and Grinstein, F.F., 2002, "Large eddy simulation of high-Reynolds-number free and wall-bounded flows", J Comp Phys, Vol. 181, pp. 68-97.
• Kraichnan, R.H., 1967, "Inertial ranges in two-dimensional turbulence", Phys Fluids, Vol. 10, pp. 1417-1423.
• Lesieur, M. and Métais, E., 1996, "New trends in large-eddy simulations of turbulence", Annu Rev Fluid Mech, Vol. 28, pp. 45-82.
• Mason, P.J. and Sykes, R.I., 1979, "Three-dimensional numerical integrations of the Navier-Stokes equations for flow over surface-mounted obstacles", J Fluid Mech, Vol. 91, pp. 433-450.
• Moin, P. and Mahesh, K., 1998, "Direct numerical simulation: a tool in turbulence research", Annu Rev Fluid Mech, Vol. 30, pp. 539-578.
• Pope, S.B., 2000, "Turbulent Flows", Cambridge University Press, Cambridge.
• Roache, P.J., 1994, "Perspective: A method for uniform reporting of grid refinement studies", J of Fluids Engineering, Vol. 116, pp. 405-413.
• Robert, A., 1981, "A stable numerical integration scheme for the primitive meteorological equations", Atmos-Ocean, Vol. 19, pp. 35-46.
• Robert, A., Yee, T.L. and Ritchie, H., 1985, "A semi-Lagrangian and semi-implicit numerical integration scheme for multilevel atmospheric models", Mon Wea Rev, Vol. 113, pp. 388-394.
• Sawyer, J.S., 1963, "A semi-Lagrangian method of solving the vorticity advection equation", Tellus, Vol. 15, pp. 336-342.
• Sorbjan, Z., 1989, "Structure of the Atmospheric Boundary Layer", Prentice-Hall, Inc., Englewood Cliffs.
• Stevens, B. and Lenschow, D.H., 2001, "Observations, experiments, and large eddy simulation", Bull Amer Meteor Soc, Vol. 82, pp. 283-294.
• Stull, R.B., 1995, "An Introduction to Boundary Layer Meteorology", Kluwer Academic Publishers, The Netherlands.
• Wilcox, D.C., 2000, "Turbulence Modeling for CFD", 2^nd ed., DCW Industries, Inc., Anaheim.
• Wygnanski, I., Champagne, F. and Marasli, B., 1986, "On the large-scale structures in two-dimensional, small-deficit, turbulent wakes", J Fluid Mech, Vol. 168, pp. 31-71.

Simulation of two-dimensional high reynolds number wake flows with the use of filtered lagrangian navier-stokes equations

Ricardo Carvalho de Almeida^I; José Luis Drummond Alves^II; Clemente Augusto S. Tanajura^III

Publication Dates

• Publication in this collection
  02 May 2011
• Date of issue
  Mar 2011