Numerical simulation of 3 D backward facing step flows at various Reynolds numbers

The work deals with the numerical simulation of 3D turbulent flow over backward facing step in a narrow channel. The mathematical model is based on the RANS equations with an explicit algebraic Reynolds stress model (EARSM). The numerical method uses implicit finite volume upwind discretization. While the eddy viscosity models fail in predicting complex 3D flows, the EARSM model is shown to provide results which agree well with experimental PIV data. The reference experimental data provide the 3D flow field. The simulations are compared with experiment for 3 values of Reynolds number.


Introduction
The numerical simulation of 3D turbulent flows in applications requires a high degree of modelling the turbulence in the framework of averaged Navier-Stokes equations.The model should provide acceptable accuracy so that used computational resources are not wasted due to its errors.The commonly used eddy-viscosity turbulence models are not entirely acceptable for 3D simulations due to the fact that they can not correctly capture secondary flows.The class of explicit algebraic Reynolds stress (EARSM) models gives promising results at a cost not much higher than that of eddy-viscosity model.Along with this more accurate mathematical model, similar degree of accuracy is required from numerical approximations.The higher order upwind finite volume methods are standard and has been used by some of present authors in implicit form for steady as well as unsteady simulations [1,2].The higher accuracy can be expected from more rigorous formulation of boundary conditions inherent in finite element methods (FEM).The authors published results [3] of the streamline upwind/pressure stabilizing Petrov-Galerkin (SUPG/PSPG) FEM applied to 2D backward facing step flows and compared results with present finite volume method.
The finite volume method is also applied to 3D simulation of turbulent flow over perpendicular backward facing step where the PIV experimental data were kindly provided by their authors Uruba and Jonáš [4,5].The published simulation extends the older results presented earlier [3].

Mathematical model
We consider the incompressible viscous laminar or turbulent flow in 3D Cartesian coordinates.The governing equations are the Reynolds averaged Navier-Stokes equations a e-mail: petr.louda@fs.cvut.cz(RANS) which formally differ from laminar model by the Reynolds stress terms.The RANS equations can be written where u = u(x, t) denotes mean velocity vector, p = p(x, t) mean kinematic pressure, , ν the constant kinematic viscosity, t time, S = 1 2 (∇u + ∇ T u) is the mean strain rate tensor and T is the Reynolds stress tensor.
The Reynolds stress tensor T can be modelled by an eddy-viscosity model where with I being unity tensor, ν t eddy-viscosity and k turbulent energy.In this work, the k-ω model in the SST variant by Menter [6] is considered.The more general constitutive relation has been proposed by Wallin and Johansson [7] in the form of the explicit algebraic Reynolds stress model (EARSM), where with Ω being mean rotation rate tensor and τ turbulent time scale.The turbulent scale are predicted e.g. by a k-ω system.In this work the Hellsten proposal is adopted [8].
The turbulent scales as k and τ are modelled by twoequation k-ω system in the form where the production terms P k = ν T S(u) : S(u), and P ω = f (P k , k, ω), and C D is the cross-diffusion term.The k-ω system can be integrated up to the wall where k = 0 and ω is sufficiently large in order to obtain viscous sublayer.

Numerical method
The artificial compressibility method [9] is used to obtain equation for pressure.It consists of adding a pressure time derivative term a −2 ∂p/∂t to the continuity equation, where a is selectable positive parameter, making the inviscid part of the system of equations hyperbolic.The parameter a in this work is chosen equal to the maximum inlet velocity.This value ensures good convergence to steady state but is not large enough to make the transient accurate in time.Therefore it is suitable for steady flows only.The discretization is done by a cell-centered finite-volume method (FVM) with hexahedral finite volumes (7) where the unknown cell average W i, j,k ≡ (p, u 1 , u 2 , u 3 , k, ω) i. j,k and F I , F V are inviscid and viscous flux through boundaries of finite volume (i, j, k).The ΔV denotes the size of the finite volume.The convective terms are discretized using the third order van Leer upwind interpolation in the direction of grid lines.The viscous terms are approximated using second order central approximation on dual finite volume grid of octahedrons constructed over each face of the primary finite volume.A pressure diffusion is added to the right hand side of the modified continuity equation in order to improve the pressure-velocity coupling [10].The Eq. ( 6) is integrated in time by the backward Euler (implicit) scheme where the steady residual Rez(W) n+1 i, j,k is linearized by the Newton method with respect to (i, j, k)-finite volume and its six neighbours with highest weights.The linearized system then has a block-diagonal structure, where the size of blocks equals the number of coupled RANS equations, i.e. 4 in 3D.The system is solved by a block relaxation method.For more details see [1].The turbulence model equations are solved in the same manner, however, decoupled from the RANS equations.The decoupling saves memory and computational time since there is no need for matrix-vector multiplication if the source terms Jacobi matrix is chosen diagonal.

Computational results
The numerical simulations correspond to PIV experiments by Uruba and Jonáš [4] and Uruba [5] where backward facing step in relatively narrow channel was investigated.The experiment can provide time averages in different planes which can be combined into 3D time averaged flow-field.
Three values of free-stream velocity are considered: U e = {4.5, 10.9, 21.7} m/s.The corresponding Reynolds numbers based on step height H = 25 mm and U e are 7 200, 17 600 and 34 400 respectively.The width of the channel equals 4H.
The simulation is carried out by finite volume method with the EARSM turbulence model and with SST model for comparison.The height of the channel is truncated to 5H (upstream of the step) with symmetry boundary condition on the top plane.The length of the domain is 4H and 25H upstream and downstream of the step respectively.The grid has 96 steps in transverse direction, with total number of approx.1.2 × 10 6 hexahedral finite volumes.The grid is highly refined near walls with minimum finite volume thickness 2 × 10 −4 H.The time step for the implicit scheme was Δt ≈ 2.5 × 10 −3 H/U.
The experimental data extend around 6H downstream of the step and this portion of flow field is shown in following figures.The isolines of velocity are printed in red for positive values and in blue for the rest.The increment of isolines is 0.02, where 1 corresponds to free-stream velocity.
The figures 1 and 2 show horizontal and vertical component of velocity in the center-plane.The computation differs by the presence of secondary vortex in the corner of recirculation zone.It might not be sufficiently steady to appear in experimental data, and the RANS approach in principle is least reliable in cases at the edge of statistical steadiness, where the unsteadiness is not sufficiently random.
The next figures show isolines of velocity in horizontal planes behind the step at different heights above the bottom wall.Only half of the span is shown, so that the centerplane is lower boundary of the figures and on the left and top there is wall.Experimental and computational results are shown.The vertical component of velocity is shown in Figs. 3 and 4 respectively in planes 0.2H and 0.5H above the bottom wall.The longitudinal component of velocity is shown in Figs. 5 and 6.The computational results agree with experiment in positions of local maxima of velocities and also the trends with changing Reynolds number seems captured well.
The figure 7 shows longitudinal velocity very close (1.25 mm) the bottom wall.The secondary recirculation in the corner is present in simulations whereas is smaller in experiment where it disappears completely at highest Reynolds number.The down-stream boundary of the main recirculation zone moves slightly closer to the step with increasing Re which is opposite trend than appearing from measurement.
It has been remarked that an eddy-viscosity model predicts unreliably secondary flows.Here we compare results of the SST model with the EARSM model at Re = 17 600.The results of SST in the center-plane seem acceptable, see Figs. 8, 9, only the SST model typically predicts stronger flow in the recirculation zone.However the next figures 10,11,12,13 show that the SST model gives qualitatively different results.Finally the boundary of main recirculation zone, as shown in figure 14, is longer with SST model and has wrong shape near side walls.

Conclusions
The paper presented results of numerical simulation of turbulent flow over backward-facing step in a narrow channel.The turbulence was modelled by the EARSM turbulence model.The results are compared with PIV measurement at 3 Reynolds numbers.The comparison in terms of mean velocity components shows good agreement qualitatively as well as quantitatively, except for region close to step corner.It should be noted however that also measurement in this region is less reliable because of lower degree of statistical steadiness.Further it has been shown that an eddy-viscosity model provides wrong results considering secondary flow in recirculation zone, although the length of recirculation zone is predicted with acceptable accuracy.