Penalty Methods for Turbulent Flows Interacting with Obstacles

Chapter 10 - Companion Computer Programs

The fast and accurate simulation of laminar and turbulent incompressible flows is crucial in science and engineering. Low-order numerical strategies, particularly finite-volume (FV) methods based on cell-centred and vertex-centred approaches, remain essential in CFD, especially in industrial settings, due to their favourable trade-off between computational cost and accuracy. However, these methods face challenges such as stabilising convective-dominated flows, handling velocity-pressure coupling, and managing numerical flux reconstruction on distorted and stretched grids. The face-centred finite volume (FCFV) method emerges as an alternative to standard FV methods. Derived from a mixed formulation of the discontinuous Galerkin method, FCFV avoids flux reconstruction at cell faces, making its accuracy and convergence nearly insensitive to grid quality. It also satisfies the Ladyzhenskaya-Babuˇska- Brezzi (LBB) condition without special treatment for velocity-pressure coupling. This work showcases the FCFV method for simulating laminar and turbulent incompressible flows for the first time. The formulation is based on the Reynolds-averaged Navier-Stokes (RANS) equations coupled with the negative Spalart-Allmaras (SA) model. Three new convective stabilisations inspired by Riemann solvers are proposed, along with monolithic and staggered solution strategies for the RANS-SA system and two relaxation strategies for pseudo-time marching. A new hybrid pressure FCFV formulation is also introduced to improve the FCFV accuracy in laminar flow simulations by enriching the pressure space through relaxed compressibility conditions. Both FCFV and hybrid pressure FCFV achieve first-order con- vergence of velocity, velocity gradient tensor, and pressure, accurately predicting engineering quantities such as drag and lift on structured and unstructured meshes. By avoiding gradient reconstruction, these methods are less sensitive to mesh quality, even on highly distorted grids. Numerical benchmarks for laminar and turbulent, steady and transient cases assess the performance, accuracy, and robustness of the proposed methodologies. Implemented in Fortran 90, these methods lay the foundation for the integration of FCFV techniques within the CFD community. The work concludes with a detailed discussion of the Fortran 90 FCFV code implementation.

In this paper, the effects of adiabatic and isothermal conditions on the statistics in compressible turbulent channel flow are investigated using direct numerical simulation (DNS). DNS of two compressible turbulent channel flows (Cases 1 and 2) are performed using a mixed Fourier Galerkin and B-spline collocation method. Case 1 is compressible turbulent channel flow between isothermal walls, which corresponds to DNS performed by Coleman et al. (1995). Case 2 is the flow between adiabatic and isothermal walls. The flow of Case 2 can be a very useful framework for the present objective, since it is the simplest turbulent channel flow with an adiabatic wall and provides ideal information for modelling the compressible turbulent flow near the adiabatic wall. Note that compressible turbulent channel flow between adiabatic walls is not stationary if there is no sink of heat. In Cases 1 and 2, the Mach number based on the bulk velocity and sound speed at the isothermal wall is 1.5, and the Reynolds number based on the bulk density, bulk velocity, channel half-width, and viscosity at the isothermal wall is 3000. To compare compressible and incompressible turbulent flows, DNS of two incompressible turbulent channel flows with passive scalar transport (Cases A and B) are performed using a mixed Fourier Galerkin and Chebyshev tau method. The wall boundary conditions of Cases A and B correspond to those of Cases 1 and 2, respectively. Case A corresponds to the DNS of Kim & Moin (1989). In Cases A and B, the Reynolds number based on the friction velocity, the channel half-width, and the kinematic viscosity is 150. The mean velocity and temperature near adiabatic and isothermal walls for compressible turbulent channel flow can be explained using the non-dimensional heat flux and the friction Mach number. It is found that Morkovin's hypothesis is not applicable to the near-wall asymptotic behaviour of the wall-normal turbulence intensity even if the variable property effect is taken into account. The mechanism of the energy transfers among the internal energy, mean and turbulent kinetic energiesis investigated, and the difference between the energy transfers near isothermal and adiabatic walls is revealed. Morkovin's hypothesis is not applicable to the correlation coeffcient between velocity and temperature fluctuations near the adiabatic wall.

Two solution methods for compressible flows and their application in Large-Eddy Simulation (LES) of turbulent flows are discussed. After a summary of the basic equations and subgrid scale models, the applied discretization methods, solution schemes and boundary conditions with their implementation on high-performance computers are presented. The two discretization methods for the advective terms used in the simulations are a second-order mixed central-upwind discretization based on a AUSM formulation and a higher-order compact finite difference scheme. First, results of the second-order scheme are shown for the case of a turbulent channel and pipe flow, which are only weakly influenced by subgrid scale models, but are in good agreement with DNS and experimental data. Higher-order compact schemes are then used for both laminar and turbulent flows and the results are compared to those from simulations with the second-order scheme. It is found that for the laminar flow and the LES of plane turbulent jets the higher-order compact scheme is not of principal advantage. This is primarily related to the sensitivity of the compact scheme to stretched and curvilinear grids, which requires explicit filtering steps to avoid unphysical oscillations in the solution. Further results are then presented, which are obtained with the computationally cheaper second-order scheme. The results include turbulent flows in pipe bends with different Reynolds and Dean numbers. The solutions show swirl switching with a distinct low frequency peak. Finally, turbulent flows in non-trivial geometries are considered. This includes a first step towards the LES of the flow in an internal combustion engine with four valves and a realistic geometry.

A novel finite volume discretization method for advection–diffusion systems on stretched meshes

Incompressible laminar and turbulent flows over flat plates and airfoils have been investigated numerically and experimentally in unsteady flow conditions. Important differences were found between laminar and turbulent flat plate flows over a wide range of oscillation frequencies. Also, the importance of unsteady effects on laminar boundary layers was found to diminish rapidly with increasing longitudinal pressure gradients, whereas turbulent separation on airfoils was significantly affected by oscillatory motion when the incidence approached the stall angle. The calculated hysteresis in turbulent separation followed in a qualitative sense the well-known trends of dynamic stall delay and reattachment. However, the numerical analysis failed to indicate some of the important features of dynamic stall observed in the present experiment and in previous studies.

This paper presents a new family of semi-implicit hybrid finite volume/finite element schemes on edge-based staggered meshes for the numerical solution of the incompressible Reynolds-Averaged Navier–Stokes (RANS) equations in combination with the k−ε turbulence model. The rheology for calculating the laminar viscosity coefficient under consideration in this work is the one of a non-Newtonian Herschel–Bulkley (power-law) fluid with yield stress, which includes the Bingham fluid and classical Newtonian fluids as special cases. For the spatial discretization, we use edge-based staggered unstructured simplex meshes, as well as staggered non-uniform Cartesian grids. In order to get a simple and computationally efficient algorithm, we apply an operator splitting technique, where the hyperbolic convective terms of the RANS equations are discretized explicitly at the aid of a Godunov-type finite volume scheme, while the viscous parabolic terms, the elliptic pressure terms and the stiff algebraic source terms of the k−ε model are discretized implicitly. For the discretization of the elliptic pressure Poisson equation, we use classical conforming P1 and Q1 finite elements on triangles and rectangles, respectively. The implicit discretization of the viscous terms is mandatory for non-Newtonian fluids, since the apparent viscosity can tend to infinity for fluids with yield stress and certain power-law fluids. It is carried out with P1 finite elements on triangular simplex meshes and with finite volumes on rectangles. For Cartesian grids and more general orthogonal unstructured meshes, we can prove that our new scheme can preserve the positivity of k and ε. This is achieved via a special implicit discretization of the stiff algebraic relaxation source terms, using a suitable combination of the discrete evolution equations for the logarithms of k and ε. The method is applied to some classical academic benchmark problems for non-Newtonian and turbulent flows in two space dimensions, comparing the obtained numerical results with available exact or numerical reference solutions. In all cases, an excellent agreement is observed.

An accurate and reliable solution procedure is presented for solving the two-dimensional, compressible, unsteady boundary layer equations. The procedure solves the governing equations in a coupled manner using a fully implicit finite-difference numerical algorithm. Several unsteady compressible and incompressible laminar flows are considered. Example results for two unsteady incompressible turbulent flows are also included. An algebraic mixing length closure model is used for the turbulent flow calculations. The computed results compare favorably with experimental data and available analytical/numerical solutions.

We give a brief introduction to research on adaptive computational methods for laminar compressible and incompressible flows and then focus on computability and adaptivity for turbulent incompressible flow, where we present a framework for adaptive finite element methods with duality‐baseda posteriorierror control for chosen output quantities of interest. We show in concrete examples that outputs such as mean values in time of drag and lift of a bluff body in a turbulent flow are computable to a tolerance of a few percent, for a simple geometry using some hundred thousand mesh points and for complex geometries using some million mesh points.

A method of efficiently computing turbulent compressible flow over complex two‐dimensional configurations is presented. The method makes use of fully unstructured meshes throughout the entire flow field, thus enabling the treatment of arbitrarily complex geometries and the use of adaptive meshing techniques throughout both viscous and inviscid regions of the flow field. Mesh generation is based on a locally mapped Delaunay technique in order to generate unstructured meshes with highly stretched elements in the viscous regions. The flow equations are discretized using a finite element Navier‐Stokes solver, and rapid convergence to steady state is achieved using an unstructured multigrid algorithm. Turbulence modelling is performed using an inexpensive algebraic model, implemented for use on unstructured and adaptive meshes. Compressible turbulent flow solutions about multiple‐element aerofoil geometries are computed and compared with experimental data.

The main objective of this study is to clarify the effect of thermal wall boundary conditions on turbulence statistics and structures in a compressible turbulent flow. This work is an extension of Morinishi et al. (J. Fluid Mech. vol. 502, 2004, p. 273), who performed DNS of compressible turbulent channel flow between adiabatic and isothermal walls at Mach number $M\,{=}\,1.5$ (Case 2). We address the question of whether the modification of turbulence statistics is attributable to the effect of the adiabatic wall boundary condition or the effect of the increase of wall temperature caused by the adiabatic wall boundary condition. New DNS of the compressible turbulent channel flow between isothermal walls with the wall temperature difference at the Mach number $M=1.5$ (Case 1) and DNS of the corresponding incompressible turbulent flow with passive scalar transport (Case I) are performed. The present study shows that the mean temperature profile near the high-temperature wall for Case 1 has an additional maximum due to the friction work, while such an additional maximum does not appear for Cases 2 and I. The additional maximum leads to a corresponding near-wall maximum of temperature fluctuations. We find the direction of energy transfer due to pressure work near the adiabatic wall for Case 2 being opposite to that near the isothermal wall to be due to the effect of the high-temperature wall, not to the effect of the adiabatic wall. These findings are explained by using the budgets of internal energy and temperature variance transport equations.

This paper presents results from a series of numerical studies of various geometry riblets in both laminar and turbulent external flows, which complement results already reported by colleagues at IMST, ONERA-CERT, and UMIST, for V and L-groove riblets in laminar internal and external flows. Although these earlier studies suggested that a 2–3% reduction could be achieved in the latter case, questions remained as to the sensitivity of this finding to further grid refinement, particularly as no drag reduction was indicated for fully developed internal flows. New, higher resolution, CRAY 2 V-groove riblet computations, employing both orthogonal cartesian and conformal grids, have tended to confirm an optimum 2.5% drag reduction for riblets scaling on the boundary layer thickness (rather than wall units as for turbulent flow), and computed mean velocity profiles agreed well with LDA experimental data. Conformal mesh computations for semi-circular U, L, and intermediate form, riblets indicated a similar drag reduction in laminar flow, but also revealed important differences in the spanwise distribution of Cf across these riblets compared to the V-groove type. Subsequent turbulent computations for V and L-groove riblets, employing a conformally transformed mixing length, resulted in mean velocity profiles which again compared well with the limited available turbulent flow data. Parametric computations revealed that the maximum drag reduction was considerably overpredicted for the V-groove riblets, but the optimum performance of the L-groove riblets and the variation in drag reduction with h+ were rather better predicted. Any discrepancies could be associated with secondary motions set-up inside the riblets, since these may be responsible for additional momentum transport, and hence higher Cf, but were not accounted for in the modelling.

Accelerate the convergence of turbulent flows simulation: A novel progressive locally power-law preconditioning method

Numerical modeling of unsteady turbulent separated flows of incompressible fluid around vertical axis wind turbines (VAWT) with Darrieus and Savonius rotors was performed on the basis of the Reynolds-averaged Navier-Stokes equations, closed by the SpalartAllmaras turbulence model. The numerical solution was obtained using original homemade program code that implements an implicit finite volume algorithm based on the artificial compressibility method. Verification of developed CFD code has been carried out by simulation of several typical testing problems, namely: laminar and turbulent flow around steady and unsteady cylinder, subcritical and supercritical flow around steady and oscillating NACA 4412, NACA 0015 airfoils, aerodynamics of the turbulent flow around a 30P30N multi-element airfoil in cruise, takeoff and landing configurations, turbulent flow around tractor-trailer. The results of Darrieus and Savonius rotors simulation, physical and computational experiments with different number and geometry of blades are presented in the form of detailed flow visualization and reconstruction for two- and three-blade Darrieus rotor. The analysis of flow field near rotors was carried out and the main stages of vortex development were identified and analyzed. The results of three series of computational experiments on aerodynamic of two and three-bucket Savonius rotor were obtained, presented and discussed.

A gas-kinetic scheme coupled with SST model for turbulent flows