LES - IB analysis of a flow in channel with an adverse pressure gradient.

This paper presents results of Large Eddy Simulation (LES) of a flow in a channel with an adverse pressure gradient. The applied computational code (SAILOR) is based on a high-order compact finite difference scheme on half-staggered meshes and is combined with an Immersed Boundary (IB) method. The IB method is implemented using the direct forcing approach in a simplified stepwise variant. The results are verified using the experimental data and the general agreement is good. The observed discrepancies are small close to the inlet section and increase along the channel length where strong flow separation occurs.


Introduction
Finite volume methods are the most popular in industry applications. Their flexibility and possibility of the use of various type unstructured grids allow carrying out simulations in very complex geometries. Unfortunately, a generation process of high quality unstructured meshes is often very difficult task. The level of complexity increases when one want to generate an unstructured mesh for Large Eddy Simulation, as in this case a required large number of the computational cells makes this process very time consuming. Additionally, considering cases with moving elements with dynamically modified meshes the need to generate a new mesh every time step complicates and extends the overall solution times significantly. A relatively simple solution of these problems consist in the use of an Immersed Boundary (IB) method [1], which allows performing the computations in geometries that do not exactly match the grid nodes, as shown in Fig. 2. In this work we combine the IB method with a high-order LES code and apply the LES-IB approach to study a problem of flow separation in a channel with an adverse pressure gradient.
Understanding of the separation phenomenon is an important and ambitious task. The regions of separation include wide ranges of turbulent high-energy structures influencing on turbulence production near the walls, which makes the flow in the separation region very difficult to analyse experimentally. In this view the LES method not only complements the measurements but also gives deep insight into every detail of the flow. The LES requires the use of non-dissipative high order methods [2] and a compact finite difference method used in the applied code (SAILOR) belongs to this family. The weak point of these methods is that they can be applied in cases with rather simple flow domains. To avoid this limitation in the present work we use LES-IB approach, which, as mentioned above enables performing simulations in almost arbitrary domain shapes.

Governing equations and solution algorithm
For incompressible Newtonian fluid flow in the LES framework the momentum and continuity equations are given as: where the bar symbol denotes spatially filtered variables [3]; u -velocity components, Ppressure, ρ,ν,ν t -density, viscosity and turbulent viscosity. The termf IBM i is a body-force term, which is used in IB method based on a direct forcing approach [4]. The spatial discretization is performed using 6th order compact difference scheme on half-staggered meshes. The time integration is performed with 2nd order Adams-Bashfort/Adams-Moulton scheme with the projection method for pressure-velocity coupling. The subgrid viscosity is modeled by the model proposed be Vreman [5]. The correctness of the applied discretization and temporal integration methods and their combination with the IB method was confirmed in laminar and turbulent flow problems [6,7].

Solution algorithm
High-order compact finite differences are used for spatial discretization which are based on half staggered meshes. In that approach all components of the velocity are in the same points in space and pressure nodes are shifted half a cell-size in respect to them. More details can be found in [8].
The solution algorithm is schematically presented in Figure 1. It can be seen that the IB method acts four times during one algorithm loop. Its influence is expressed in the equation 2 byf IBM i terms. In the present calculations the simplest stepwise IB approach is used. The schematic view of that method is presented in figure 2. In this approach the solid body is determined by the locations of computational nodes. In the nodes inside the body (Ω b ) velocity is imposed by equalizing it to the object velocityū Ω , in the nodes outside the body (Ω f ) all velocity components are computed using the governing equations. It is easy to see that the shape of the body has a stepwise appearance (dark gray color in figure 2).

Sub-grid modelling
An appropriate sub-grid modelling is a very important aspect of LES. In wall bounded flows it is required that the turbulent viscosity vanishes in near wall regions. Correctness of the results obtained using a model proposed in [5] within the LES-IB approach was confirmed in the simulations of a flow turbulent in a converging-diverging channel [6]. According to [5] the turbulent viscosity is calculated using the formula: where the model constant C is equal to 0.025.

Immersed Boundary method
The forcing terms used in the equation 2 are formulated using an assumed velocity of the immersed object, u D i , and they take the forms: in the corrector step. The forcing terms (6) and (8) are used in the Navier-Stokes equations, whereas (7) and (9) in the velocity correction steps. Using the above formulas, it turns out that in the immersed body region the computed velocity, e.g. u * * i , is equal to the assumed velocity u D i . Hence, it is easy to see that there is no need to compute the forcing terms using the equations (6-9), but one can directly impose the velocity u D i on u * i , u * * i , u n+1 i .

Domain and mesh configuration
A computational domain analyzed in this paper is shown in the figure 3. The flow is from the left to right side and the adverse pressure gradient is caused by a curvature of the upper wall and the grid mounted in the outlet section. The complexity of the flow in this seemingly simple configuration is illustrated in the figure 3 by the contours of instantaneous values of the streamwise velocity component and the iso-surfaces of Q-parameter. The Q-parameter is very good indicator of the vortex structures [9] and is defined as Q = 1 2 (S ij S ij − Ω ij Ω ij ), where S ij and Ω ij are symmetric and antisymmetric parts of the velocity gradient tensor. The walls and the solid region modelled by the IB method are marked by the violet color. The width of domain and the height of inlet section are equal to 0.32 m. In the z direction the periodic boundary conditions are applied. At the inlet, the mean velocity profile of the streamwise velocity component corresponding to experimental data is imposed with 5% noise. In the experiment, in order to increase the total pressure inside the channel a grid with 40% clearance was used at the outlet section. In the performed computations this grid was mimicked using the IB method. The mean Reynolds number based on the friction velocity on the lower wall is equal Re τ = 3000, whereas the Reynolds number based on the mean inlet velocity and the height of the channel is equal to Re = 2 × 10 5 , approximately. The sizes of the domain and the mesh parameters are given in Table 1. In the directions y and z the mesh is uniform, in the direction x the mesh is defined in the following way. In the region near the bottom wall (0 -0.32 m) the mesh is compacted. The y + of the first node is equal 1.5, in the physical coordinates it is 0.0005 m. In the region of the curved wall the mesh is uniform. The smooth transition between these regions was obtained by distributing the mesh nodes using a high-order polynomial functions.      regions in the boundary layer close to the bottom and upper walls are apparent. The sub-filter viscosity is calculated in a classical way, i.e. without taking into account the use of the IB method. However, it can be seen that the sub-grid viscosity inside the solid region is virtually equal to zero. This is due the fact that the velocity gradients inside the solid body are very small and as these values are used to calculate the eddy viscosity (see formula 4) it becomes almost zero as well. On the other hand, as one could expect large values of the eddy viscosity occurs in the separation region, where also the vorticity has the biggest value and the flow is regarded as the most turbulent.

Comparison with experimental data
The comparison with experimental data are presented for three locations from the inlet, i.e. at 0.4, 0.8 and 1.2 m. The time averaging procedure started when the flow in the channel was fully developed and lasted for 6 seconds. Figures 7, 8 and 9 show the profiles of streamwise velocity and the symbols correspond to the experimental data. In general the results obtained using the LES-IB method are close to the experimental points. As can be seen the best results are for the first line (0.4 m from inlet) and they get worse along the channel length. This is probably due to the fact that in the location of first line the channel is straight and the influence of the wall curvature on the solution is insignificant. On the other hand, in the regions where the separation occur on the upper and lower wall the discrepancies are clearly visible. Most likely, the number of points used in the present studies is insufficient to correctly capture the physics of the flow in the presence of separation. The grid sensitivity studies are planned for future, however, the present results have already proven that applying LES-IB allows analysing the flow separation phenomena caused by the adverse pressure gradients.

Conclusions
This paper presented the results of LES computations of the flow in the channel with adverse pressure gradient. The applied computation code was based on high-order compact methods in combination with the Immersed Boundary method. The simulation conditions and domain shape corresponded to experimental configuration and the obtained results were in acceptable agreement with the experimental data. The discrepancies close to the inlet section were small and increased along the channel length. It is likely, that this could result from to coarse computational mesh or not exact representation of the experimental inlet velocity fluctuations.
In the computations the inlet fluctuations were modelled by a white noise. This usually causes damping of the turbulent intensity downstream, which in the present case could result in inaccurate prediction of the separation point. Future research will focus on this subject and on mesh sensitivity studies.