On the turbulent flow past a realistic open-cell metal foam

Abstract Turbulence is investigated in the lee of an open-cell metal foam layer. In contrast to canonical grids, metal foams are locally irregular but statistically isotropic. The solid matrix is characterised by two lengths, the ligament thickness $d_f$ and the pore diameter $d_p$. A direct numerical simulation is conducted on a realistic metal foam geometry for which $d_f/d_p = 0.14$ and the porous layer thickness is five times the pore diameter. The Reynolds number based on the pore size is ${\textit {Re}}_{d_p} = 4000$, corresponding to a Taylor-scale Reynolds number ${\textit {Re}}_\lambda \approx 80$. Closer to the foam than two pore diameters, the pressure and turbulent transports of turbulent kinetic energy are non-negligible. In the same region, ${\textit {Re}}_\lambda$ undergoes a steep decrease whereas the dissipation coefficient $C_{\epsilon }$ increases like ${\textit {Re}}_\lambda ^{-1}$. At larger distances from the porous layer, the classical grid turbulence situation is recovered, where the mean advection of turbulent kinetic energy equals dissipation. This entails a power-law decay of turbulent quantities and characteristic lengths. The decaying exponents of integral, Taylor and Kolmogorov scales are close to one-half, indicating that the turbulence simulated here differs from Saffman turbulence. Analysis of the scaling exponents of structure functions and the decorrelation length of dissipation reveals that small-scale fluctuations are weakly intermittent.

grids have attracted a lot of interest because of the specific type of turbulence generated. A remarkable increase in the Reynolds number based on the Taylor scale with respect to usual passive grids was noticed by Seoud & Vassilicos (2007). Further studies investigating the decaying turbulence downstream of a set of multiscale grids (Krogstad & Davidson 2012) and a square-element fractal grid (Hearst & Lavoie 2014) have revealed that, while the region close to the grids can be characterised by residual inhomogeneity and is grid-dependent, in the far field -where development is accomplished -flow characteristics are in accordance with classical grid turbulence measurements.
The in-depth study of turbulence generated by fractal grids has revealed a peculiar behaviour of the dissipation coefficient C = L/u 3 (here is the rate of dissipation.) in a region close to the turbulence-generating grid but in conjunction with energy spectra, which follow the −5/3 slope for a wide range of wavenumbers. This behaviour has been described as a breakdown of the classical dissipation scaling and is observed in wind tunnel experiments of grid-generated turbulence of different geometries (Mazellier & Vassilicos 2010;Isaza, Salazar & Warhaft 2014;Valente & Vassilicos 2014;Mora et al. 2019) and also in direct numerical simulation (DNS) data of decaying turbulence in a periodic box (Goto & Vassilicos 2016). The breakdown consists in a behaviour of C at the initial stages of the decay that depends upon the inlet Reynolds number Re M and the local Reynolds number Re λ as follows: C ∼ Re 1/2 M /Re λ . This implies that L/λ ∼ Re 1/2 M along the direction of Re λ decay (Valente & Vassilicos 2012).
The turbulent flow behind an open-cell metal foam has never been investigated before. Similar to regular grids, the open-cell metal foam geometry is characterised by two main length scales, i.e. the mean pore diameter d p and the mean ligament thickness d f . In contrast to regular grids, open cells are arranged randomly in space and their morphology, based on a polyhedral frame, is never exactly repeated. This generates a structure that is highly irregular and anisotropic at the pore scale but statistically isotropic at the macroscale. In addition, the metal foam layer investigated here has a thickness larger than a single ligament or a pore. Moreover, the solidity of metal foams, measured as 1 − ε, where ε represents grid porosity, is very different from the solidities typical of grids employed for grid turbulence, which generally range between 0.30 and 0.45. In high-porosity metal foams, this value is typically lower than 0.10 (Calmidi & Mahajan 2000).
In this paper, a description is reported of turbulence downstream of a high-porosity open-cell metal foam. After a qualitative introduction to the flow investigated in § 3.1, the degree of homogeneity and isotropy of turbulence is investigated in § 3.2. The power-law behaviour of turbulent kinetic energy k and its dissipation rate are assessed in § 3.3 and § 3.4, respectively. The streamwise variations of turbulent length scales are considered in § 3.5. In § 3.6 the behaviour of the dissipation rate coefficient is investigated in detail, also in the context of recent research on the topic. In § 3.7 the discussion is about whether or not Saffman turbulence is achieved for the present flow configuration. The multiscale nature of the flow is examined by means of one-dimensional spectra in § 3.8. Then § 3.9 investigates high-order statistics, the role of intermittency and the decorrelation length of dissipation. Budget terms of turbulent kinetic energy and of velocity variances are reported in § 3.10, where the main mechanisms for the transport of energy and fluctuations are described in the vicinity of the solid structure and in the fully developed region. Final remarks are drawn in § 4.

Mathematical formulation and flow configuration
The system of equations solved numerically comprises the mass conservation and the Navier-Stokes equations for incompressible flows: complemented with appropriate boundary conditions. The subscripts i and j take the values 1, 2 and 3 to denote the streamwise direction, x, and the two cross-flow directions, y and z, respectively. All the variables are made non-dimensional by the velocity at the inlet U ∞ and the mean pore diameter of the metal foam d p . The Reynolds number based on the unperturbed velocity and the mean pore diameter is Re d p = 4000. The governing equations (2.1) are solved using the high-order finite-difference method implemented in Incompact3d (Laizet & Lamballais 2009). Sixth-order compact schemes are used for spatial discretisation (Lele 1992), and time integration is done by the third-order Adams-Bashforth scheme. The velocity field is evaluated on a Cartesian grid with uniform spacing along the three directions, and pressure is defined on a staggered grid. Pressure-velocity decoupling is accomplished by a fractional-step method, which determines the divergence-free velocity field by solving a Poisson equation. The Poisson problem is tackled in the spectral space by using a modified wavenumber formalism, which allows for any kind of boundary conditions for the velocity field in the physical space. Inflow/outflow boundary conditions are enforced along the streamwise direction and periodic boundary conditions are set along the cross-flow directions to represent statistical homogeneity. A uniform velocity field is prescribed at the inlet (U ∞ , 0, 0), while at the outlet the velocity is determined by a convection equation: The convection velocity c is calculated at each time step as the mean between the maximum and minimum values of the streamwise velocity component at the outlet. The representation of the intricate metal foam geometry is achieved via an immersed boundary method (IBM) based on direct forcing (Gautier, Laizet & Lamballais 2014). The IBM enforces the no-slip condition at the solid walls while preserving the simplicity of the finite-difference schemes applied to the Cartesian grid (Laizet & Lamballais 2009). Further details on the numerical methodology employed here are provided in Corsini & Stalio (2020).
A sketch of the computational domain is displayed in figure 1. The extents of the domain along the streamwise and the two cross-flow directions are L x = 45d p and L y = L z = 11.25d p , respectively. The origin of the coordinate system is located at the centre of the downstream face of the porous matrix; thus x = 0 describes the most upstream cross-flow plane where the fluid is not in contact with the solid phase. The thickness of the metal foam layer in the streamwise direction is 5d p and it spans the whole domain in the cross-flow directions. It is placed at a 5d p distance from the inlet section; this avoids interference between the upstream boundary and the solid matrix. The computational domain is discretised by n x = 3073 grid nodes in the streamwise direction and n y = n z = 768 in the cross-flow directions. The spatial resolution is sufficiently fine to ensure that x = y = z 2η for x 5. Close to the porous layer (0 < x < 5), where dissipation is larger, x = y = z 5η. In the above comparisons, the Kolmogorov microscale η is calculated a posteriori from its definition. The time step is kept at t = 0.001d p /U ∞ during the simulation, which in terms of the Kolmogorov time scale τ η yields t 0.033τ η . This corresponds to a Courant-Friedrichs-Lewy number CFL < 0.3. Statistical quantities are computed by averaging in time and along the homogeneous y and z directions. Gathering of statistics begins after one flow-through time from the start of the simulation. In order to obtain well-converged statistics, the time interval of collection is T = 225d p /U ∞ , and the three-dimensional snapshots of the velocity and pressure field are sampled at equal time intervals of T = 4.5d p /U ∞ .

Metal foam geometry
The problem of the computer modelling of an intricate metal foam porous structure has been tackled in different ways. Pore-scale morphology can be reconstructed through X-ray tomography (Piller et al. 2013) or generated mathematically assuming an ideal cell geometry based on a virtual sample of regular polyhedra (Boomsma, Poulikakos & Ventikos 2003). In this study, where geometric periodicity is a key feature of the numerical representation and irregularity of the foam is a requirement, a third approach is adopted: the open-pore cellular structure is generated synthetically through a numerical algorithm, developed by August et al. (2015). Besides the excellent realism and periodicity, one further favourable feature of synthetic metal foams is the possibility to tune their porosity and permeability. Thanks to a diffuse interface representation of the phase-field approach (August et al. 2015), the thickness of ligaments can be easily adjusted. Figure 2 shows the details of a couple of cells of the synthetic structure.
The synthetic metal foam structure used in the simulation is characterised by the geometrical features listed in table 1. Both d p and d f are calculated by spatial averages and thus represent the mean pore diameter and the mean ligament thickness of the metal foam sample. The value of grid porosity set, ε = 0.92, is representative of high-porosity metal foams (Calmidi & Mahajan 2000). Based on typical sizes of open-cell aluminium foams, an inflow velocity of U ∞ = 15 m s −1 is obtained at Re d p = 4000, assuming air at standard conditions as the working fluid. Table 1 reports experimental conditions from previous wind tunnel studies on regular planar grids. Figure 2. Cells of the aluminium foam generated algorithmically (August et al. 2015).
A sample of the metal foam geometry with superimposed computational points is displayed in figure 2 of Corsini & Stalio (2020). Staircase patterns of the immersed boundaries approximate the rounded borders of the solid region. This referenced picture and the computed ratio between average ligament diameter and grid spacing d f / x ≈ 10 suggest that the ligaments are discretised by an adequate number of grid points. Figure 3(a) shows the streamwise component of the instantaneous velocity field in one of the snapshots collected. The uniform free stream of the inflow is disrupted by the irregularly arranged ligaments of the solid structure and velocity fluctuations arise within the porous matrix. Vortices of different orientations are shed from the ligaments and a wake is formed. The largest perturbations are observed close to the downstream edge of the porous matrix. The wakes originated by the ligaments develop in a non-uniform fashion and interact at variable lengths. The smaller wakes are seen to disappear after a couple of pore diameters, whereas larger wakes stemming from ligament clumps meet at a further distance from the foam. The larger wakes also last in time, as revealed by the streaks in the time-averaged velocity field U t shown in figure 3(b).

Homogeneity and isotropy
The approximation to statistical homogeneity in the cross-flow directions in grid turbulence is known to depend on the grid geometry and the Reynolds number. While, for regular grids, experiments suggest that the flow becomes nearly homogeneous for x/M > 40 (Comte-Bellot & Corrsin 1966;Mohamed & Larue 1990), for fractal grids, homogeneity is usually retrieved further downstream (Hearst & Lavoie 2014) and, for example, Valente & Vassilicos (2011) The distribution of U t in cross-flow planes is shown in figure 4 for six streamwise positions at increasing distance from the metal foam. Solid lines mark regions where the percentage of variation of U t relative to U ∞ has magnitude greater than 10 %, while dashed lines encompass regions of magnitude greater than 5 %. These are seen to gradually shrink along x. While inhomogeneity is in part to be ascribed to the limited size of the sample composed only by the collection of snapshots in time, for x > 20 their extent is   Krogstad & Davidson (2010) and Kitamura et al. (2014) on the turbulence generated by regular grids are also included. Here M, d and σ denote the mesh width, the rod thickness and the solidity of the grid, respectively. Asterisk * denotes quantities expressed in dimensional form. Quantities in parentheses indicate values that have been deduced from other quantities provided in the same work. still appreciable. In the x = 30 station, the time-averaged velocity U t varies between 0.86 and 1.12. The isotropy level of the large scales of motion can be investigated through the ratio of root-mean-square (r.m.s.) velocity fluctuations along orthogonal directions. In the present case, the fluctuations of the x-component of velocity are the largest: in the developed region, indicators u rms /w rms and u rms /v rms oscillate within the interval (1.5, 1.6) about a mean of 1.55 for u rms /v rms and 1.56 for u rms /w rms , where the difference is finally due to the size of the sample employed in the simulations. In previous grid turbulence measurements (Kurian & Fransson 2009;Krogstad & Davidson 2010;Kitamura et al. 2014), the observed isotropy indicators are in general smaller than those measured here. Kitamura et al. (2014), who also collected experimental results by other authors, report u rms /w rms < 1.2 in all cases; similar values are also reported in the lee of fractal grids (Hurst & Vassilicos 2007;Gomes-Fernandes, Ganapathisubramani & Vassilicos 2012). More details about large-scale isotropy measures downstream of the present metal foam are reported in Corsini & Stalio (2020). Figure 5 displays the streamwise evolution of the variance of the three components of velocity u i u i as well as the turbulent kinetic energy k . Very close to the foam and for x < 1, velocity fluctuations are observed to remain constant. The negative slope increases gradually in x until fluctuations exhibit a power-law decay that persists until the end of the computational domain. As demonstrated, for example, in Tennekes & Lumley (1972), this is expected in the region where advection and dissipation of the turbulent kinetic energy become the only non-negligible terms in the transport equation of k ; see § 3.10.

Decay of velocity fluctuations
Power-law parameters are sought in the form through a numerical procedure. In (3.1), A is the multiplicative coefficient, x 0 is the virtual origin and n is the decay exponent. As n is positive, the power law has a vertical asymptote (and a singularity) at the virtual origin x = x 0 . As the parameters in (3.1) depend greatly upon the interval of sampling data considered, also the interval limits are determined inside the fitting procedure. A similar approach has been applied in Hearst & Lavoie (2014).
In the present work, a developed region I d = [x min , x max ] is employed, where the right border of the interval is kept fixed to x max = 30.0, clear of possible -yet not evidentoutflow condition effects, while x min is discretely varied in the interval [0.015, 24.7] to seek an x min coordinate that ensures the best fit. The coordinate x min will be taken as the start of the developed region. For each selection of a value for x min , the virtual origin x 0 is discretely varied within I 0 = [0.015, x min ], as the singularity is not supposed to belong to I d . The intervals of variation of both x min and x 0 are discretised by the same subdivision as the computational mesh. Parameters A and n are determined through a least-squares fit. Deviations between computed data and fitting laws are then calculated as the Euclidean norm of the error divided by the number of data points: In (3.2), δ j represents the difference between computed statistics and the least-squares fitted power-law approximation of u 2 rms at the jth point of I d ; and N d is the number of uniformly spaced points in the data fit region I d . The (x 0 , x min ) couple which ensures the smallest deviation from computed data provides the final A and n coefficients. This procedure leads to the results given in table 2 with error distribution as in figure 6. In the region of parameters investigated, only one minimum is found, which is located far from the boundaries of the region investigated. The dependence on porosity of the power-law exponents obtained is shown to be only weak in the Appendix.
In order to check the dependence of the results from the fitting method employed, a procedure from the literature is also applied to the same set of data; this alternative technique is that utilised by Hearst & Lavoie (2014). In the application of the method to the present case, the virtual origin x 0 and lower bound x min are varied within I 0 = [0, 4] and [4, 14.7], respectively. For both u 2 rms and k , the process converges after three iterations. Applied to u 2 rms , this fitting procedure provides x 0 = 0.610 and x min = 6.78, which yield the estimate forñ u = 1.12. In the case of k , the values obtained are x 0 = 0.360, x min = 5.02 andñ k = 1.14. The results obtained by the method proposed by Hearst & Lavoie (2014) are only marginally different from those obtained as described above and reported in table 2.
Comte-Bellot & Corrsin (1966) found 1.15 n 1.29 for regular grids, while according to Mohamed & Larue (1990) n ≈ 1.3. More recently, Krogstad & Davidson (2010) found n = 1.13 ± 0.02. Besides the parameters x 0 , A and n in (3.1), the procedure provides the coordinate x min of the start of the developed region as the left boundary of the interval I d . All the subsequent fittings in this work will be carried out on I d = [7.98, 30.0].

Turbulence decay rate
As opposed to experimental studies, where rate of dissipation needs to be evaluated using the frozen turbulence assumption or isotropy (3.5), in DNS studies can be computed directly from its definition, is the fluctuating rate of strain. Figure 7 displays the spatial evolution of , together with the least-squares fitted power law in the form ∼ a x h . The coefficients that fit the data over I d = [7.98, 30.0] are a = 0.217 and h = −2.20.
The scaling of dissipation can also be set in relation to the scaling of k in § 3.3. As also shown quantitatively in § 3.10, in the developed region of a statistically steady, high-Reynolds-number flow, one has Equation (3.4) suggests that the decay exponent in this case should equal h = −(n k + 1).
In the present study, h = −2.20 and n k = 1.14 are calculated. The percentage difference between −(n k + 1) and h is below 3 %. Dissipation is compared in figure 7 to the same quantity computed under the hypothesis of isotropic turbulence: (3.5) The close similarity between and iso is only in apparent contradiction with isotropy ratios larger than 1.5 reported in § 3.2. From the demonstration by Taylor (1935), it appears that hypotheses less strict than isotropy are required for equation (3.5) to hold. The weaker set of hypotheses hold true in the present case; see Corsini & Stalio (2020).

Length scales
The length scales examined for the characterisation of the turbulence generated by a metal foam are the Kolmogorov scale η, the Taylor microscale λ and the integral scales. All the length scales are computed directly from their definitions. Their distribution within the developed region is approximated by a power law of the distance x. The range of variation of each length scale along I d is reported in table 1 together with data from the literature on classical grid turbulence.

Kolmogorov scale
The Kolmogorov scale is defined through dissipation . It is predicted from the power-law expression for (derived from the expression for k ) that η can be represented by a function of the form η ∼ a η x s , where s equals −h/4 and finally s = (n k + 1)/4. The decay exponent for k computed here, n k = 1.14, gives s = 0.54. The power-law approximation is displayed in figure 8 Figure 9 displays the streamwise distribution of the Taylor microscale, defined by Taylor-scale values in a few measurement points  (see their table 4); their data fit a power law of exponent c = 0.53. It is predicted in the study by George (1992) that the Taylor scale of homogeneous isotropic turbulence increases in time with t 1/2 and, for grid turbulence, λ grows as x 1/2 in the laboratory frame. More recently, Kurian & Fransson (2009) reported that λ increases approximately like the square-root of the streamwise coordinate. In the present study, fitting a power-law approximation of the form λ ∼ a λ x c over I d gives a λ = 0.0577 and c = 0.52.  Results from the present investigation, and not reported for brevity, show that the difference between (λ/η) 2 / √ 15 and Re λ is less than 3 % over the whole computational domain.

Integral scales
The autocorrelation coefficient is defined as the ratio between the autocorrelation function of separation r = re j and the autocorrelation function for r = 0: The autocorrelation coefficients along y of the streamwise ρ 11 (x, re 2 ) and the cross-flow ρ 22 (x, re 2 ) velocity components are reported in figure 8 of Corsini & Stalio (2020). A distinction is made between the transverse and longitudinal correlations depending on the relative direction of velocity components and separation vector. Transverse correlations built with streamwise velocity fluctuations have longer tails with respect to longitudinal correlations built with spanwise fluctuations. Integral scales are defined as integrals over r of autocorrelation coefficients: (3.9) These depend upon the coordinates along non-homogeneous directions as well as on the direction of separation e j . In practice, since the computational domain has finite boundaries, the integral scales are calculated here as the distance over which the autocorrelation function decreases from 1 to 1/e, where e is Euler's number. For the calculation of integral scales, correlations are assumed to decay exponentially. Variable Power law Multiplicative coefficient Exponent Table 3. Parameters of the power-law functions of the form f (x) ∼ ax b fitting the turbulent quantities analysed.
Integral scales for the present case are based on the streamwise or the cross-flow velocity components. Given the inhomogeneity of the streamwise direction, separation is set in the cross-flow directions e 2 and e 3 . Thus, the integral scale based on the streamwise velocity is always transverse and will be denoted by L. The integral scales based on cross-flow velocity components can be either longitudinal or transverse and are denoted by L g and L t , respectively. As statistically L 11 (x, e 2 ) = L 11 (x, e 3 ) = L(x), L 22 (x, e 2 ) = L 33 (x, e 3 ) = L g (x) and L 33 (x, e 2 ) = L 22 (x, e 3 ) = L t (x), these equalities have been exploited in the calculation of integral scales. Figure 11 displays the integral scales. The power-law approximations L ∼ a L x q have the coefficients reported in table 3; the exponents computed are very close to the L ∼ x 1/2 behaviour predicted by Wang & George (2002). Also, in the case of higher porosity presented in the Appendix, the power-law exponent is close to 1/2.
Notice that L is one order of magnitude smaller than the domain size, the transverse length of which is 11.25d p ; this suggests that the imposed lateral periodic boundary conditions can satisfactorily represent cross-flow homogeneity in the present case. The evolution of the Reynolds number based on the integral scale, defined as Re L ≡ Lu rms /ν, along the x-axis is displayed in figure 9 of Corsini & Stalio (2020).

Dissipation rate coefficient
In high-Reynolds-number turbulent flows away from solid walls, the dissipation rate can be scaled on the integral length scale and velocity fluctuations through an order-one constant: (3.10) Figure 12 shows that, in the present case, after an initial steep increase for x < 2, the dissipation rate coefficient C based on L fluctuates over I d between 0.45 and 0.50, where its spatial average isC = 0.483. This value is very close to values reported in Pearson, Krogstad & van de Water (2002) for shear turbulence at different Re λ numbers. With regard to the initial steep increase, in recent years Vassilicos and coworkers have observed that, close to the turbulence-generating grid, there is a region characterised by spectra that closely match the −5/3 power law, in conjunction with an increase in C . In the hypothesis of = iso , combining the definition (3.10) with equation (3.5) leads to As only small variations of the ratio between length scales L/λ are observed for given Reynolds number of the mesh size, C is seen to increase like Re −1 λ . This behaviour is reported for both fractal and regular grids (Valente & Vassilicos 2012). Equation (3.11) is studied for the present case in figure 13, where the logarithmic plot emphasises the C (Re λ ) behaviour close to the metal foam. Corresponding to the initial steep decrease in Re λ shown in figure 10, C is seen to increase like Re −1 λ . In the same region, a well-defined −5/3 energy spectra behaviour is observed over a broader wavenumber range than the fully developed region; see the inset of figure 15 in § 3.8. The situation described in Valente & Vassilicos (2012) is thus observed in the present case.
The streamwise evolution of C experiences a transition between the Re −1 λ behaviour and a region where the variations in C are much smaller; see figures 12 and 13. This transition is about x = 2. The turbulent kinetic energy budgets reported in § 3.10 indicate that, downstream of this location, the turbulent transport terms become negligible and the mean advection of k equals dissipation. On the contrary, for x < 2, this equality does not hold, and the variations of C suggest a non-equilibrium condition between energy at the large scales and dissipation. It should be noted that the present results confirm the predictions reported by Tennekes & Lumley (1972) in their (3.2.29) and (3.2.30), where transition is expected at streamwise distances from the grid that are much larger than the integral scale. In the present configuration, x = 2 corresponds to x ≈ 6L. In the theory by Richardson and Kolmogorov, the constancy of C requires that turbulence is at a high Reynolds number and far from solid walls. While the small variations of C in the fully developed region can be attributed to the Reynolds number, which is not very high, the steep increase in C in the vicinity of the porous matrix (x < 2) is to be ascribed to the vicinity of the solid filaments and ultimately to non-negligible turbulent transport terms in the budget equation of k .

Is grid turbulence Saffman turbulence?
In recent articles (Krogstad & Davidson 2010;Kitamura et al. 2014) it was discussed whether grid turbulence can be considered to be of the Saffman type. Both Krogstad & Davidson (2010) and Kitamura et al. (2014) conclude that grid turbulence is Saffman turbulence.
The theory by Saffman (1967) describes the decay of homogeneous turbulence as u 2 rms = KC 2/5 t −6/5 , L = K C 1/5 t 2/5 , (3.12a,b) where K and K are constants and C is expressed by an invariant integral. Both equations (3.12a,b) have to hold for turbulence to be of the Saffman type. This is sometimes expressed by the requirement that u 2 rms L 3 = const. during decay, but the latter is a necessary, not sufficient, condition for Saffman turbulence.
The exponents reported in tables 2 and 3 suggest that turbulence investigated in the present study is not of the Saffman type. Figure 14 provides graphical confirmation for this conclusion. In the fully developed region of grid turbulence, as well as in the case investigated here, the advection of turbulent kinetic energy is almost perfectly balanced by dissipation: (3.13) (see § 3.4). Setting n k = n u = n, (3.13) combined with the definition of C in (3.10) and the hypotheses L ∼ x q and C ∼ x f leads to the following relation: As is apparent from (3.12a,b) and (3.14), grid turbulence generated at n = 6/5 is not of the Saffman type unless C stays constant during kinetic energy decay (f = 0).

Spectral scaling and energy transfer
One-dimensional spectra E ii (κ j ) are obtained from the discrete Fourier transform of the two-point velocity correlation functions along e j , where κ j is the jth component of the wavenumber vector. Only j = 2 and j = 3 are used here because the streamwise direction is not homogeneous. Depending upon the pairing between velocity and wavenumber components, three distinct spectra can be calculated: streamwise spectrum E s = E 11 (κ j ), While, in general, spectra at different stages of the decay could be expected to scale with different reference quantities, it is observed here that, as factors ηu 2 η , λu 2 rms and Lu 2 rms evolve at very similar rates in the streamwise direction (see table 3), any of those can be used equivalently. Spectra scaled by λu 2 rms and evaluated at increasing distance along the x-axis are displayed in figure 15. Figure 16 compares the three types of spectra (E s , E g and E t ) at a coordinate x = 20. For κ 2 η > 0.2, the streamwise and transverse spectra almost coincide, which is confirmation that anisotropy is at large scales only; see Mohamed & Larue (1990) and Corsini & Stalio (2020). Only a narrow range in the wavenumber space (0.025 < κ 2 η < 0.1) is noticed where E s exhibits a −5/3 behaviour. Figure 16 includes the longitudinal spectra E 11 (κ 1 ) of Comte-Bellot & Corrsin (1971) for two regular grids of different mesh size and Re λ values in the same range as the present case (Re λ = 65 and 41, compared to the present Re λ = 81). Close agreement is observed between measurements in turbulence generated by classical grids and turbulence simulated in the high-porosity metal foam for κ 2 η > 0.1. The present results match canonical grid turbulence spectra for κ 2 η > 0.1, the −5/3 law is not observed for wavenumbers κ 2 η > 0.1, and local isotropy is observed for κ 2 η > 0.2. That is the range of scales where dissipation becomes non-negligible. The distinction between scales containing the bulk of the energy from those responsible for dissipation is done by considering the energy spectrum E(κ) and the dissipative spectrum D(κ) = 2νκ 2 E(κ).

The turbulence spectrum E(κ) is obtained by
where the summation convention is applied. The spectrum in (3.15) removes the directional information from both the velocities and the Fourier modes, as E(κ) is given as a function of the wavenumber magnitude κ = |κ|. The kinetic energy cumulated at wavenumbers lower than κ is the complement to k (0,κ) is indicated by k (κ,∞) . The corresponding quantities for the dissipation are obtained in the same fashion using D(κ) and are indicated by (0,κ) and (κ,∞) . In figure 17 the fraction of cumulative turbulent kinetic energy at wavenumbers higher than κ and the fraction of cumulative dissipation at wavenumbers lower than κ at a distance from the foam x = 20 are depicted as functions of κη and of the corresponding wavelength /η = 2π/κη. The peak of the energy spectrum occurs at κη ≈ 0.02 (see figure 16) but it may be observed that the main fraction of the kinetic energy (k (κ,∞) = 0.1k) is contained in the range of wavenumbers up to κη ≈ 0.25, one decade further. The peak of dissipation is roughly at κη ≈ 0.25 (which corresponds to the maximum derivative in (0,κ) / ). The contribution to dissipation reaches (0,κ) = 0.9 at κη ≈ 0.7. The bulk of turbulent kinetic energy is contained in motions of length scales > 25η ≈ 1 2 L; this range could be viewed as the energy-containing range. On the other hand, the dissipation is effective at the length scales 10η < < 50η. The overlap between k (κ,∞) and (0,κ) reveals that energy starts to be dissipated at a length scale where the energy content is non-negligible. 3.9. Structure functions In this section, the local structure of turbulence is investigated at x = 25 by analysing the scaling properties of the structure functions with separation r. The general definition is given by (3.17) Because of turbulence decay along x, two different structure functions are identified in this work: longitudinal structure functions (δu g ) p = (δu j,j ) p and transverse structure functions (δu t ) p = (δu i,j ) p , where i = 2, j = 3 or i = 3, j = 2. According to the first, original theory by Kolmogorov (1941), for high Reynolds numbers when the separation lies in the inertial subrange η r L, the moments of velocity difference (δu g ) p take a universal form that depends only on through the following scaling property: In the refined similarity theory (Kolmogorov 1962), the intermittency effects are taken into account by rewriting expression (3.18) in terms of r , the dissipation averaged over a volume of linear dimension r, and assuming that r has a log-normal distribution, where μ is the exponent of the dissipation autocorrelation function, again for inertial range separations (Monin & Yaglom 1975). Given the Reynolds number Re λ ≈ 80 of the present case, the study is conducted using the extended self-similarity (ESS) observation by Benzi et al. (1993). The pth moments that the effect generated by intermittent structures on small-scale turbulence in the moderate-Reynolds-number range is less intense than at high Reynolds numbers. It is noticed that the method based on sixth-order structure functions, (3.21), displays a decay in the x-direction that is not expected nor recovered in the calculation of μ directly from the definition (3.20).
The decorrelation scaler is defined as the length of decorrelation of the instantaneous dissipation , i.e. the separation scale r where (x + r) (x) / 2 becomes unitary with a 1 % error. In figure 21 the development ofr along x is compared to that of the integral scale L. The decorrelation scale can be considered as a very large length scale, which in homogeneous isotropic turbulence depends on the Reynolds number and also the intermittency characteristics of the flow.

Turbulent energy budgets
The equation governing the transport of k in a statistically steady case can be written in symbols as where A is the contribution by mean advection and T p , T t and D v represent pressure, turbulent and diffusive transports; P and ˜ stand for the production and pseudo-dissipation rate of turbulent kinetic energy, defined as (3.23) As the difference − ˜ is seldom important (Pope 2000), in this section is used. In grid turbulence cases, where mean velocity gradients are zero and cross-flow directions y and z are homogeneous, the production term vanishes: In addition, for the Reynolds number computed here, the molecular contribution to transport is negligible. The resulting expression for the budget equation (3.22) thus becomes where the single terms are distributed along the x-axis as depicted in figure 22. Two regions can be identified in figure 22: a near-field region where the three conservative terms redistribute the kinetic energy along x, and a far-field region where x pressure and turbulent transports become negligible. Thus dissipation is solely balanced by turbulent kinetic energy advection; see figure 22(a) and its inset. A possible downstream boundary for the near-field region can be set at x = 2, where the turbulent and pressure transports become 100 times smaller than advection and dissipation. The near-field region can be subdivided into two subregions. By considering the ligament diameter d f as the length scale, a region very close to the metal foam is identified for x * /d f < 0.5, where turbulent transport is larger than the advective one. This region is characterised by the sustainment of turbulence by means of streamwise fluctuations, which is counter-balanced by dissipation, as pressure and advective transports account for negligible shares. The second subregion extends for 0.5 < x * /d f < 14, where turbulent kinetic energy is provided through pressure and advective mechanisms while dissipation and turbulent transport drain k ; see figure 22(a).
Data gathered by Norberg (1998) indicate that the vortex formation length f for the present Reynolds number based on the mean diameter of the filaments, Re d f = 560, lies in the range 1.5d f < f < 2d f ; see Bloor, Gerrard & Lighthill (1966) for other possible definitions of the vortex formation length. Figure 22(a) displays that this range matches the streamwise location where all the transport terms in (3.25) peak. Therefore, the present data suggest that kinetic energy is originated in the second part of the near-field region, where transport mechanisms show local peaks. Turbulent kinetic energy k is drained from here by the turbulent transport, which transfers turbulent fluctuations upstream, as indicated by the negative sign of the turbulent flux (see figure 22b), feeding the region very close to the metal foam. On the other hand, pressure and advective mechanisms are found to provide turbulent kinetic energy in the entire region 0.5 < x * /d f < 14, and the signs of the relative fluxes indicate that velocity fluctuations are transferred downstream. At the coordinate x * /d f = 14 (x = 2), fluxes up and uk become constant in the streamwise direction, leading to the budget reported in (3.4) for x > 2, which is typical in developed decaying flows; see figure 22(a) and its inset, which reports the budget terms in the I d region.
In order to separately assess the behaviour of streamwise and cross-flow velocity fluctuations, the transport equations for velocity variances are analysed. The budget of streamwise velocity variance reads (3.26) where, from left to right, the first three terms, respectively, represent the advective, turbulent and pressure transports of u 2 , p∂u/∂x is the pressure strain term and u stands for the u-variance pseudo-dissipation rate, u = 2((∂u/∂x j ) (∂u/∂x j ))/Re d p . The diffusive transport is again neglected with respect to the other terms. Equations similar to (3.26) can be derived for the transverse velocity components, v and w. As the flow is isotropic on y-z planes, statistics for the transverse velocity component are obtained by averaging results along y and z. For conciseness, the transverse velocity component is indicated by v and its direction by y, (3.27) where the terms have the same interpretation as in (3.26), and the pseudo-dissipation rate of the v-variance is computed as v = 2((∂v/∂x j ) (∂v/∂x j ))/Re d p . Note that, with respect to equation (3.25), budgets of velocity component variances include the pressure strain terms. As will be shown in the following, the role of these terms is to redistribute kinetic energy among the velocity components and thus it is responsible for the 'return to isotropy' in homogeneous turbulence. Such a phenomenon has been studied for several decades as it is involved in second-order turbulence models; see, for example, the works by Rotta (1951) and Lumley & Newman (1977). Pressure strain terms are not included in (3.25) as their sum S u + 2S v vanishes because of incompressibility. Figures 23 and 24 show the velocity variance budgets and fluxes along the streamwise coordinate. It appears that transport mechanisms are more intense in the u 2 budget with respect to v 2 , while dissipations u and v are almost equal. The behaviour of the pressure strain terms reveals that turbulent energy is drained from the streamwise velocity fluctuations, as S u represents a sink term in the u 2 budget, and provided to cross-flow velocity fluctuations, where S v acts as a source; see figures 23(a) and 24(a). The role of S u and S v does not change along the whole streamwise extension of the domain, as shown by the insets in figures 23(a) and 24(a). This occurs because, as reported in § 3.2, flow anisotropy is conserved in the flow domain considered. In the decaying region I d , the pressure strain terms maintain an intensity comparable to the advective and dissipation terms, suggesting that the isotropic condition would have been reached in a longer computational domain.
The profiles in figures 23 and 24 show that the u-variance terms behave very similarly to the budgets and fluxes of turbulent kinetic energy reported in figure 22. In particular, the peaks of turbulent, advective and pressure transports are located at the same distance from the metal foam, comparable to the vortex formation length. This suggests that velocity fluctuations are more intensely triggered along the streamwise direction with respect to the transverse ones. On the other hand, the v-variance terms peak slightly downstream with respect to terms in the u 2 and k budgets, and the negative peak of turbulent transport is not very intense (see figure 24a), but it drains enough energy to sustain cross-flow fluctuations in the region just downstream of the metal foam. Another difference with respect to the u 2 terms is the positive -yet not that intense -turbulent flux very close A

Conclusions
An analysis is presented of the turbulent flow behind a synthetic metal foam layer of thickness equal to five times the mean pore diameter. Unlike classical grid turbulence geometries, the metal foam ligaments are variably oriented, unevenly spaced and in general less ordered. Similar to classical grids, metal foams are mainly characterised by two length scales i.e. the pore diameter and the ligament thickness. The analysis encompasses one single Reynolds number Re d p = 4000 and ε = 0.92 but the effect of a different porosity (ε = 0.97) on selected quantities is addressed in the Appendix.
The analysis of the turbulent kinetic energy budget suggests that turbulence is triggered in the region close to the metal foam, approximately at a distance equal to the vortex formation length indicated by Norberg (1998), x * /d f ≈ 2. In the near-field region, x < 2, turbulent kinetic energy is distributed by transport mechanisms; the turbulent transport moves k upstream, sustaining fluctuations in close proximity to the metal foam, while pressure and advective transports provide velocity fluctuations to larger x-coordinates. At x = 2, pressure and turbulent transports are found to be negligible with respect to the other terms, entailing for x > 2 the expected behaviour of grid turbulence, where viscous dissipation is balanced by the mean advection of k . Budgets of streamwise and cross-flow velocity variances indicate that fluctuations along x are the most intensely triggered and pressure strain terms act to redistribute turbulent energy towards the isotropic condition, which however is not observed due to the limited domain extension.
In that same near-field region, Re λ decreases steeply and the dissipation rate coefficient C approximates C ∝ Re −1 λ while L/λ remains almost constant. These observations, in conjunction with the typical −5/3 slope in the turbulent power spectra, represent a typical behaviour already observed in the literature (Valente & Vassilicos 2012). Given that this occurs very close to the porous matrix, where the hypotheses for the Richardson cascade are not verified, a considerable variation in C is not interpreted here as an anomalous behaviour.
The developed region is defined here as the region where u rms decays following a power law and starts at x min = 7.98. Besides k , a number of relevant quantities, like the integral length scale, the Taylor microscale and the Kolmogorov scale of the flow, are observed to follow a power-law behaviour. The exponents obtained from a least-squares fit are in many cases very close to values predicted and measured in classical grid turbulence experiments. Given n k = 1.14 and q = 0.52 (where n k and q are the exponents of the power laws for k and L), the constancy of u 2 rms L 3 does not hold in the developed region simulated here, and thus the turbulence does not follow the theory by Saffman.
Structure functions of different orders are calculated at a fixed position x = 25 in the fully developed region. The extended self-similarity is used to calculate the scaling exponents. The results are interpreted in the light of the refined similarity hypothesis. The intermittency exponent μ computed directly from its definition is seen to behave uniformly within the computational domain. It reveals that intermittency is not very intense at the moderate Reynolds number of this study. The decorrelation length of dissipationr is larger than the integral scale. It depends on the Reynolds number as well as the intermittency characteristics of the flow. body of the paper. The two foams have equal pore size but different ligament thickness, i.e. the foam with larger porosity is characterised by a thinner ligament, d f < d f . The comparison is conducted by means of DNS performed with the same numerical parameters as outlined in § 2 except for the computational grid, which consists of n x = 2049 and n y = n z = 512 grid nodes. The porous matrix of higher porosity generates a turbulent field where fluctuations are less intense in all the spatial directions. Application of the fitting procedure outlined in § 3.3 on the high-porosity foam provides decay exponents for u 2 rms and k of n u = 1.05 and n k = 1.08, respectively. The results obtained on the coarse grid for the ε = 0.92 case yield n u = 1.12 and n k = 1.13; see table 2 for comparisons against results on the fine mesh. The power-law decay of the turbulent kinetic energy begins more upstream than in the ε = 0.92 case; see figure 25.
Inside the developed region, the Kolmogorov and the longitudinal integral length scales evaluated for ε = 0.97 behave according to power laws with growing rates estimated by the exponents s = 0.54 and q = 0.47, respectively (see table 3 for notation). The lower-porosity case exhibits decay exponents along the x-direction of s = 0.56 and q = 0.51 (see table 3 for comparisons against results on the fine mesh) but the Kolmogorov and the integral scales are different in magnitude, as displayed by figure 26(a) and (b). Associated with the smaller kinetic energy content of the flow and a smaller dissipation rate observed at ε = 0.97, the turbulence induced by the higher-porosity foam is characterised by a larger Kolmogorov scale and a smaller integral scale. Figure 27 displays the effects of porosity on the terms of the turbulent kinetic energy budget. The transport mechanisms of k for ε = 0.97 are the same as described in § 3.10 for ε = 0.92 but are characterised by an upstream shift of the transport peaks with respect to that case. This shift is consistent with the reduction of the vortex formation length associated with a smaller ligament thickness, as indicated by Norberg (1998) and reported in § 3.10.
In summary, while the results obtained for ε = 0.92 show small quantitative differences with respect to ε = 0.97, no significant discrepancies are observed between the two cases. An upstream shift in the peaks of k transports characterises the higher-porosity case, which in turn implies a reduced vortex formation length and the achievement of power-law decays at shorter distances from the porous matrix. The size of such displacement is apparently very small.