Eddy Viscosity Assumption Research Articles

Are random forests better suited than neural networks to augment RANS turbulence models?

Machine-learning (ML) techniques have bloomed in recent years, especially in fluid mechanics applications. In this paper, we trained, validated and compared two types of ML-based models to augment Reynolds-averaged Navier–Stokes (RANS) simulations. The methodology was tested in a series of flows around bumps, characterized by different levels of flow separation and curvatures. Initially, the ML-based models were trained in three configurations presenting attached flow, small and moderate separation and tested in two configurations presenting incipient and large separation. The output quantity of the machine-learning model is the turbulent viscosity as done in Volpiani et al. (2022). The new models based on artificial neural networks (NN) and random forest (RF) improved the results if compared to the baseline Spalart–Allmaras model, in terms of velocity field and skin-friction profiles. We noted that NN has better extrapolation properties than RF, but the skin-friction distribution can present small oscillations when using certain input features. These oscillations can be reduced if the RF model is employed. One major advantages of RF is that raw quantities can be given as input features, avoiding normalization issues (such as division by zero) and allowing a larger number of universal inputs. At the end, we propose a mixed NN-RF model that combines the strengths of each method and, as a result, improves considerably the RANS prediction capability, even for a case with strong separation where the Boussinesq hypothesis (and therefore the eddy-viscosity assumption) lacks accuracy.

Numerical and experimental study of two-vehicle overtaking process

The overtaking process between two simplified Ahmed body models is studied using a set of static experiments to provide a comparison basis for the numerical study. Both steady-state (static) and transient (dynamic) simulations are performed. There are notable differences of the drag and side force coefficients between the predicted and experimental results in the wake effects on the overtaking model as well as under a situation when two models are parallel for the static tests. These differences are attributed to the SST [Formula: see text] model, which is a two-equation Reynolds-averaged Navier-Stokes equation (RANS) turbulence model with an assumption of isotropic eddy viscosity, employed in the simulation. A more advanced turbulence model that simulates large-eddy organized motion in a turbulent flow and which uses an anisotropic eddy viscosity formulation in the near-wall flow subregion would increase the prediction accuracy for the overtaking simulation. A comparison of the results for the steady-state and the transient simulations at various relative velocity ratios ([Formula: see text]) shows that the quasi-steady approach is suited to simulate the dynamic two-vehicle overtaking process when the condition [Formula: see text] is met.

Effect of constant eddy viscosity assumption on optimization using a discrete adjoint method

This paper presents a thorough study of the effect of the Constant Eddy Viscosity (CEV) assumption on the optimization of a discrete adjoint-based design optimization system. First, the algorithms of the adjoint methods with and without the CEV assumption are presented, followed by a discussion of the two methods’ solution stability. Second, the sensitivity accuracy, adjoint solution stability, and Root Mean Square (RMS) residual convergence rates at both design and off-design operating points are compared between the CEV and full viscosity adjoint methods in detail. Finally, a multi-point steady aerodynamic and a multi-objective unsteady aerodynamic and aeroelastic coupled design optimizations are performed to study the impact of the CEV assumption on optimization. Two gradient-based optimizers, the Sequential Least-Square Quadratic Programming (SLSQP) method and Steepest Descent Method (SDM) are respectively used to draw a firm conclusion. The results from the transonic NASA Rotor 67 show that the CEV assumption can deteriorate RMS residual convergence rates and even lead to solution instability, especially at a near stall point. Compared with the steady cases, the effect of the CEV assumption on unsteady sensitivity accuracy is much stronger. Nevertheless, the CEV adjoint solver is still capable of achieving optimization goals to some extent, particularly if the flow under consideration is benign.

Structural aspects of the attached turbulent boundary layer flow over a hill

Three-dimensional turbulent boundary layers under strong pressure gradients and curvature are the rule in real-world flow applications but are typically not well predicted by turbulence models due to isotropic eddy viscosity or equilibrium assumptions. Validation-quality data in complex 3D flows is necessary for continued efforts to improve simulation accuracy. The Benchmark Validation Experiments for RANS/LES Investigations (BeVERLI) Hill bump model, designed specifically for validation experiments, was tested in the Virginia Tech Stability Wind Tunnel to collect validation experiment data on the three-dimensional (3D) boundary layer flow over a 3D hill. Laser Doppler velocimetry measurements on the bump model were used to study the mean flow and turbulence structure and evaluate the impact of pressure gradient and curvature upon the total shear stress in the boundary layer and evaluate the impact of pressure gradient and curvature upon the total shear stress behavior in the near-wall region. From analysis of the BeVERLI Hill flow, including the boundary layer just upstream of the hill, and comparison with the 3D flow around a wing-body junction of Ölçmen et al. (Exp Fluids 31:219–228, 2001), it is shown that none of the stations studied exhibit a constant shear stress region over any significant region of the boundary layer.Graphical abstract

Multi-Row Turbomachinery Aerodynamic Design Optimization by an Efficient and Accurate Discrete Adjoint Solver

This paper proposes an approach that combines manual differentiation (MD) and automatic differentiation (AD) to develop an efficient and accurate multi-row discrete adjoint solver. In this approach, the structures of adjoint codes generated using an AD tool are first analyzed. Then, the AD-generated codes are manually adjusted to reduce memory and CPU time consumption. This manual adjustment is performed by replacing the automatically generated low-efficient differentiated codes with manually developed ones. To demonstrate the effectiveness of the proposed approach, the single-stage transonic compressor–NASA Stage 35 and the 1.5-stage Aachen turbine–are used. The solution information exchange at a rotor-stator/stator-rotor interface is achieved by a conservative, non-reflective, and robust discrete adjoint mixing plane method. The results show that the discrete adjoint solver developed by hybrid automatic and manual differentiation is more economical in computational cost than that developed purely by an AD tool and has higher sensitivity accuracy than the adjoint solver with the constant eddy viscosity (CEV) assumption. Moreover, the multi-row turbomachinery design optimizations can be efficiently performed by the discrete adjoint solver developed by the hybrid automatic and manual differentiation.

Bayesian parameter estimation of SST model for shock wave-boundary layer interaction flows with different strengths

The Shock Wave-Boundary Layer Interaction (SWBLI) flow generated by compression corner widely occurs in engineering. As one of the primary methods in engineering, the Reynolds Averaged Navier-Stokes (RANS) methods usually cannot correctly predict strong SWBLI flows. In addition to the defects of the eddy viscosity assumption, the uncertainty of the closure coefficients in RANS models often significantly impacts the simulation results. This study performs parametric sensitivity analysis and Bayesian calibration on the closure coefficients of the Menter k-ω Shear-Stress Transport (SST) model based on the SWBLI with different strengths. Firstly, the parametric sensitivity on prediction results is analyzed using the Sobol index. The results indicate that the Sobol indices of wall pressure and skin friction exhibited opposite fluctuation trends with the increase of SWBLI strength. Then, the Bayesian uncertainty quantification method is adopted to obtain the posterior probability distributions and Maximum A Posteriori (MAP) estimates of the closure coefficients and the posterior uncertainty of the Quantities of Interests (QoIs). The results indicate that the prediction ability for strong SWBLI of the SST model is significantly improved by using the MAP estimates, and the relative errors of QoIs are reduced dramatically.

Improved Dissipation Rate Equation for Swirling and Rotating Pipe Flows

The nature of turbulent swirling and rotating flow in a straight pipe is investigated using a family of near-wall two-equation models. Specifically, the viability of three different near-wall two-equation models is assessed. These models are asymptotically consistent near the wall. The first two models, one with isotropic and another with anisotropic eddy viscosity invoked, solved a dissipation rate equation with no explicit correction made to account for swirl and flow rotation. The third model assumes an isotropic eddy viscosity but solves an improved dissipation rate equation that has explicit corrections made to account for swirl and flow rotation. Calculations of turbulent flows in the swirl number range 0.25 - 1.3 with and without a central recirculation region reveal that, with the exception of the third model, neither one of the other two models can replicate the mean field at the swirl numbers tested. Furthermore, taking stress anisotropy into account also fails to model swirl effect correctly. Significant improvements can be realized from the third model, which is based on an improved dissipation rate equation and the assumption of isotropic eddy viscosity. The predicted mean flow and turbulence statistics correlate well with measurements at low swirl. At high swirl, the two-equation model with an improved dissipation rate equation can still be used to model swirling and rotating pipe flows with a central recirculation region. However, its simulation of flows without a central recirculation region is not as satisfactory.

Optimal Design of Computational Fluid Dynamics: Numerical Calculation and Simulation Analysis of Windage Power Losses in the Aviation

Based on the theory of computational fluid dynamics (CFD), with the help of the Fluent software and the powerful parallel computing capability of the super cloud computer, the single-phase flow transient simulation calculation of the windage power loss of the engagement spiral bevel gear pair (SBGP) was performed. The two-equation SST k-ω turbulence model based on the assumption of eddy viscosity was adopted, which was improved from the standard k-ε model combined with the Wilcox k-ω model. The SST k-ω turbulence model inherited the respective advantages of the Wilcox k-ω model in the near-wall region and the k-ε model in the free shear layer and could more accurately describe the resistance and separation effect of the gear tooth surface on the airflow. The simulation analyzed the airflow characteristics around SBGP and the mechanism of the windshield to reduce the windage loss of the gear. It also studied the influence of the windshield clearance and opening size on the windage power loss. Then the orthogonal experimental analysis method was adopted to perform numerical simulation analysis. The windage torque was studied under different clearance values between the windshield and the gear tooth surface, as well as the large end and the small end. The variance analysis was performed on the numerical simulation data. The results showed that when the windshield clearance value was 1 mm and the engagement opening was 30°, the windage torque was the smallest, and the effect of reducing the windage power loss was the best. According to the changes in the pressure, velocity, and turbulent kinetic energy cloud diagram of the flow field in the reducer during multi-group simulation tests, the local optimal windshield configuration was obtained, which provided a method for further research on the multi-objective optimization of the windshield and the windage loss of the gear pair under the oil–gas two-phase flow and also provided a reference for the practical engineering application of the windshield.

Mean Current Profile over Rippled-Beds in the Presence of Non-Breaking Waves and Analysis of Its Influencing Factors

Classical eddy viscosity model deviates from the actual mean current profiles, when calculating the mean current profiles over rippled-beds in the presence of non-breaking waves, owing to the neglect of the enhancement of the wave boundary layer thickness by ripples and the wave-induced shear stress (the radiation stress and the wave Reynolds stress). Considering these shortcomings, a semi-empirical one-dimensional vertical (1DV) model is presented in this study. The present model was obtained using the two-dimensional Navier–Stokes equations and eddy viscosity assumptions, which differ from those of previous researchers, while a top-to-bottom sequence was adopted to calculate the mean current profiles. Empirical formulae were derived from the laboratory measurements and used in the present model to accurately predict the wave boundary layer thickness and bed roughness. The present model is in satisfactory agreement with the data from laboratory experiments. The factors influencing the mean current profiles were analyzed also. The wave-induced second-order shear stresses were found to be the principal reason for the deviations of the mean current profiles in the near-surface layer; as the influencing factors of wave-induced shear stress, the intensity of the wave relative to the current, the angle between the wave and current, and the size of ripples can also have a non-negligible effect on the mean current profiles.

ILES of melt turbulent convection with conjugated heat transfer in the crucible and gas flow for Czochralski silicon crystal growth system

Results of Implicit Large Eddy Simulation (ILES) of buoyancy, rotation, and surface tension driven turbulent melt convection during Czochralski (Cz) silicon crystal growth are presented. Convective and conductive heat transfer in the melt is coupled with heat transfer in surrounding crystal, crucible, and gas flow. Distributions of resolved turbulence kinetic energy (TKE), Reynolds stress tensor, production terms and dissipation rate in the silicon melt are discussed and compared to published results. We have found significant anisotropy of the melt velocity fluctuations along the melt free surface and the crucible. Conjugated heat transfer in the melt flow and in the crucible has provided important results of the temperature fluctuation intensity along the liquid/solid interface. Our ILES results have also shown that there are significant regions in the melt with negative values of the shear production term. To estimate the quality of calculations using a conventional turbulence model for a Cz configuration, we have compared our ILES data to respective results obtained within the Boussinesq eddy viscosity assumption. As well we have evaluated how precise are available approximations of the production terms and dissipation rate of TKE for Si melt turbulent convection.

Wave-Modified Ekman Current Solutions for the Time-Dependent Vertical Eddy Viscosity

This study aimed to highlight a general lack of clarity regarding the scale of the temporal averaging implicit in Ekman-type models. Under the assumption of time and depth-dependent eddy viscosity, we present an analytical Fourier series solution for a wave-modified Ekman model. The depth dependence of eddy viscosity is based on the K-Profile Parameterization (KPP) scheme. The solution reproduces major characteristics of diurnal variation in ocean velocity and shear. Results show that the time variability in eddy viscosity leads to an enhanced mean current near-surface and a decrease in the effective eddy viscosity, which finally results in an intensified near-surface shear and wakes a low-level jet flow. Rectification values are dominated by the strength of diurnal mixing, and partly due to the nonlinear depth dependence of the eddy viscosity.

Coupling Mechanism of Rotating Casing Effect and Impeller Structure of Roto-Jet Pump

An increase in internal pressure of a Roto-Jet pump is due to combined action of its impeller and rotating casing. Internal pressure of the pump chamber was determined according to outlet pressure of the impeller, and the influence of the rotary casing effect was ignored. To study the combined action mechanism of the rotating casing effect and impeller structure on the Roto-Jet pump, we used the open test bed of the Roto-Jet pump and four model pumps with impellers of different structures as research objects. We also conducted a comprehensive experimental study on the coupling mechanism between the rotating casing effect and impeller structure. Numerical calculation was performed to avoid the assumption of isotropic eddy viscosity, Reynolds stress linear pressure-strain model is selected, and the numerical calculation results are compared with the experimental results to verify its credibility. The results show that the rotating casing effect has multiple functions to reduce the friction loss of the disc, improve the pressure distribution inside the rotating casing, and increase the pump head. All scheme, pressure, and velocity fluctuations occur in the upstream and wake regions of the collecting pipe and the energy loss is concentrated in the upstream region. The difference in velocity distribution inside the collecting pipe is small and negligible. As long as the impeller and rotating casing continue to rotate synchronously, the liquid shear velocity at the same coordinate position of each scheme remains unchanged, and the liquid rotation angular velocity in the rotating casing is approximately 75% of the rotating casing which conforms to the rigid motion law. In the same scheme, the coefficient of uneven velocity inside the rotating casing gradually increases along the radial direction. The closer to the axis, the faster is the decrease of the peak tangential velocity and the velocity tends to be uniform. The size, shape, and position of the vortex core inside the rotating casing change constantly with various schemes. The distribution of vortex cores varies under each scheme. The front cover and rotating casing have a serious effect on the vortex core. The extremely poor analysis of the test results shows that the performance of the Roto-Jet pump is better when the closed impeller is rotated in synchronisation with the rotating casing. The advantages and disadvantages of each blade type can be determined according to the situation. The research results can exhibit the influence mechanism of the Roto-Jet pump shell effect. Selection of an impeller structure provides a reference.

Eddy viscosity from bottom Ekman veering profiles

ADCP profiles near the sea bed, at a 130 m depth continental shelf, were used to estimate the mean veering of the velocity vectors using a complex correlation coefficient. An analytic expression for the veering vertical gradient (∂θ/∂z) as function of turbulent Ekman layer thickness (δE) was developed and, using constant eddy viscosity assumptions, sea bottom stress and additional turbulent boundary layer features were evaluated. Further statistics on the time series yields median values for δE of about 13 m, eddy viscosity Az = 4.2×10−3 m2s−1 and median bottom stress of 9.2×10−2 Pa, thus leading to high Reynolds numbers. Very good agreement is obtained between observations and bottom Ekman balance solutions according to determination coefficients assessed (R≥0.99).

On a Proper Tensor-Diffusivity Model for Large-Eddy Simulation of Buoyancy-Driven Turbulence

In this work, we aim to shed light to the following research question: can we find a nonlinear tensorial subgrid-scale (SGS) heat flux model with good physical and numerical properties, such that we can obtain satisfactory predictions for buoyancy-driven turbulent flows? This is motivated by our findings showing that the classical (linear) eddy-diffusivity assumption, \(\varvec{q}^{eddy} \propto \nabla \overline{T}\), fails to provide a reasonable approximation for the actual SGS heat flux, \(\varvec{q}= \overline{\varvec{u}T} - \overline{\varvec{u}} \overline{T}\): namely, a priori analysis for air-filled Rayleigh-Bénard convection (RBC) clearly shows a strong misalignment. In the quest for more accurate models, we firstly study and confirm the suitability of the eddy-viscosity assumption for RBC carrying out a posteriori tests for different models at very low Prandtl numbers (liquid sodium, \(Pr=0.005\)) where no heat flux SGS activity is expected. Then, different (nonlinear) tensor-diffusivity SGS heat flux models are studied a priori using DNS data of air-filled (\(Pr=0.7\)) RBC at Rayleigh numbers up to \(10^{11}\). Apart from having good alignment trends with the actual SGS heat flux, we also restrict ourselves to models that are numerically stable per se and have the proper cubic near-wall behavior. This analysis leads to a new family of SGS heat flux models based on the symmetric positive semi-definite tensor \(\mathsf {G}\mathsf {G}^{T}\) where \(\mathsf {G}\equiv \nabla \overline{\varvec{u}}\), i.e. \(\varvec{q}\propto \mathsf {G}\mathsf {G}^{T}\nabla \overline{T}\), and the invariants of the \(\mathsf {G}\mathsf {G}^{T}\) tensor.

Stochastic flow approach to model the mean velocity profile of wall-bounded flows.

There is no satisfactory model to explain the mean velocity profile of the whole turbulent layer in canonical wall-bounded flows. In this paper, a mean velocity profile expression is proposed for wall-bounded turbulent flows based on a recently proposed stochastic representation of fluid flows dynamics. This original approach, called modeling under location uncertainty, introduces in a rigorous way a subgrid term generalizing the eddy-viscosity assumption and an eddy-induced advection term resulting from turbulence inhomogeneity. This latter term gives rise to a theoretically well-grounded model for the transitional zone between the viscous sublayer and the turbulent sublayer. An expression of the small-scale velocity component is also provided in the viscous zone. Numerical assessments of the results are provided for turbulent boundary layer flows, pipe flows and channel flows at various Reynolds numbers.

An inconsistence-free integral-based dynamic one- and two-parameter mixed model

ABSTRACTThe development of the dynamic procedure as well as the deeper understanding of the link between filtering, modelling and numerics, allowed Large Eddy Simulation (LES) to make great progresses during the last years. Among several modelling approaches, the scale-similar-based modelling is based on the observation that the smallest resolved scales are the most active in the interaction with the unresolved ones. Owing to the low dissipation introduced by the scale-similar models (SSMs), the coupling with the eddy-viscosity model is often used in the so-called mixed models. Dynamic version of mixed models is historically based on the application of the test-filtering on the differential form of the filtered momentum equation. Such an approach is used for both the one and the two-coefficients mixed models. The use of the differential form of the filtered equations produces the well-known mathematical inconsistence caused by the need to extract arbitrarily the model functions out of filtering. It is known that, along with the eddy viscosity assumption, the magnitude of the Germano identity error (GIE) is strongly influenced. The mathematical inconsistence in the extraction of the dynamic eddy viscosity coefficient was recently superseded by using the new integral-based formulation. However, owing to the intrinsic limits of the Smagorinsky model, also in those results, the GIE is still remarkable therefore, the present paper presents a new formulation to the integral-based dynamic procedure for both one and two-coefficients mixed models (IDMM). The original contributions of the present paper can be summarised: (1) A theoretical analysis comparing the spectral errors for the differential and integral-based SSM, assessing that the errors are less relevant for the integral form; (2) The implementation of one and two parameters IDMM for the simulation of turbulence in a plane channel flow, assessing the reduction of the GIE and the good behaviour of the statistics that are compared with those of the other LES codes used in the LESinItaly project.

Extending the applicability of RANS turbulence closures to the simulation of transitional flow around hydrofoils at low Reynolds number

Growing interest in marine autonomous systems such as underwater gliders and autonomous underwater vehicles (AUVs) has led to the increase of research on hydrofoil performance at low Reynolds numbers, bringing forward the issue of simulating laminar to turbulent transition using computational fluid dynamics (CFD). Most conventional turbulence closures in practical Reynolds-averaged Navier-Stokes (RANS) simulations are incapable of predicting transition with excessive turbulence generation near the leading edge, deriving from their eddy viscosity assumption. Having the potential of inherently predicting transition by solving the near wall flow accurately, three Low-Re RANS closures with anisotropy consideration are selected and evaluated for their ability to capture both separation-induced transition and natural transition. The empirical γ−Reθ transition model is also assessed in this work. Comparison with published experiments shows the capability of the four turbulence closures to capture separation-induced transition accurately. Predicting natural transition remains more challenging: two Low-Re closures exhibit the potential to capture natural transition only at medium grid finesse levels, whereas the γ−Reθ model evidences challenges in reproducing the correct flow physics. The lag elliptic blending model provides the most accurate solution owing to its consistent grid convergence, high accuracy for predicting the hydrodynamic coefficients, and ability to reproduce the correct transition physics.

Stochastic representation of the Reynolds transport theorem: Revisiting large-scale modeling

We explore the potential of a formulation of the Navier-Stokes equations incorporating a random description of the small-scale velocity component. This model, established from a version of the Reynolds transport theorem adapted to a stochastic representation of the flow, gives rise to a large-scale description of the flow dynamics in which emerges an anisotropic subgrid tensor, reminiscent to the Reynolds stress tensor, together with a drift correction due to an inhomogeneous turbulence. The corresponding subgrid model, which depends on the small scales velocity variance, generalizes the Boussinesq eddy viscosity assumption. However, it is not anymore obtained from an analogy with molecular dissipation but ensues rigorously from the random modeling of the flow. This principle allows us to propose several subgrid models defined directly on the resolved flow component. We assess and compare numerically those models on a standard Green-Taylor vortex flow at Reynolds numbers Re = 1600, Re= 3000 and Re = 5000. The numerical simulations, carried out with an accurate divergence-free scheme, outperform classical large-eddies formulations and provides a simple demonstration of the pertinence of the proposed large-scale modeling.

Compound channel flow with a longitudinal transition in hydraulic roughness over the floodplains

Flows in a compound open-channel (two-stage geometry with a main channel and adjacent floodplains) with a longitudinal transition in roughness over the floodplains are experimentally investigated in an 18 m long and 3 m wide flume. Transitions from submerged dense vegetation (meadow) to emergent rigid vegetation (wood) and vice versa are modelled using plastic grass and vertical wooden cylinders. For a given roughness transition, the upstream discharge distribution between main channel and floodplain (called subsections) is also varied, keeping the total flow rate constant. The flows with a roughness transition are compared to flows with a uniformly distributed roughness over the whole length of the flume. Besides the influence of the downstream boundary condition, the longitudinal profiles of water depth are controlled by the upstream discharge distribution. The latter also strongly influences the magnitude of the lateral net mass exchanges between subsections, especially upstream from the roughness transition. Irrespective of flow conditions, the inflection point in the mean velocity profile across the mixing layer is always observed at the interface between subsections. The longitudinal velocity at the main channel/floodplain interface, denoted $$U_{int}$$ , appeared to be a key parameter for characterising the flows. First, the mean velocity profiles across the mixing layer, normalised using $$U_{int}$$ , are superimposed irrespective of downstream position, flow depth, floodplain roughness type and lateral mass transfers. However, the profiles of turbulence quantities do not coincide, indicating that the flows are not fully self-similar and that the eddy viscosity assumption is not valid in this case. Second, the depth-averaged turbulent intensities and Reynolds stresses, when scaled by the depth-averaged velocity $$U_{d,int}$$ exhibit two plateau values, each related to a roughness type, meadow or wood. Lastly, the same results hold when scaling by $$U_{d,int}$$ the depth-averaged lateral flux of momentum due to secondary currents. Turbulence production and magnitude of secondary currents are increased by the presence of emergent rigid elements over the floodplains. The autocorrelation functions show that the length of the coherent structures scales with the mixing layer width for all flow cases. It is suggested that coherent structures tend to a state where the magnitude of velocity fluctuations (of both horizontal vortices and secondary currents) and the spatial extension of the structures are in equilibrium.

Turbulence modeling for flows in wall bounded porous media: An analysis based on direct numerical simulations

Various models are available for simulating turbulent flows in porous media. Models based on the eddy viscosity assumption are often adopted to close the Reynolds stress term. In order to validate the assumptions behind such turbulence models, we studied the dynamics of macroscopic momentum and turbulence kinetic energy in porous media flows by utilizing Direct Numerical Simulation (DNS). The generic porous matrix is composed of regularly arranged spheres. The resulting periodic porous medium is bounded by two walls. The DNS analyses with a Lattice Boltzmann method were performed for various values of the applied pressure gradient, pore size to channel width ratio, and porosity. The DNS results were averaged over time and volume to obtain macroscopic results. The results show that the macroscopic shear Reynolds stress in all Representative Elementary Volumes (REVs), independent of their location, is negligibly small, although the mean velocity gradient takes nonzero values near the wall. The turbulence kinetic energy production rate is generally balanced by the dissipation rate in each REV. The DNS results support a zero-equation turbulence model that accounts for the fact that turbulent structures are restricted in size by the pore scale. The DNS results also suggest that the Brinkman term, which expresses the diffusion of momentum, has an important effect near the wall where the gradient of the shear stress is large. Therefore, the Brinkman term should be taken into account in the macroscopic momentum equation as a component of the total drag. A preliminary macroscopic model for calculating turbulent porous media flows has been proposed and compared with our DNS results.