12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

SYMMETRY BREAKING AND TURBULENCE IN OSCILLATORY FLOW
THROUGH A HEXAGONAL SPHERE PACK

Lukas Unglehrt
Professorship of Hydromechanics
Technical University of Munich
Arcisstr. 21, 80335 Munich, Germany
lukas.unglehrt@tum.deMichael Manhart
Professorship of Hydromechanics
Technical University of Munich
Arcisstr. 21, 80335 Munich, Germany
michael.manhart@tum.de

ABSTRACT2016) or on the processes at the interface between the porous
medium and a free flow (Shigematsu et al., 2018). On the other
hand, numerical simulations of oscillatory flow through twodimensional porous media configurations were performed byGraham & Higdon (2002), Iervolino et al. (2010) and Kardgar& Jafarian (2021). In our upcoming work (Unglehrt & Manhart, 2022), the onset of nonlinear effects was investigated inlaminar oscillatory flow through a hexagonal sphere pack bynumerical simulation.

The appearance of turbulence in oscillatory flow througha hexagonal sphere pack was investigated by means of direct
numerical simulation. The Reynolds numbers lie between 26.9
and 297 and the Womersley number is 31.62. We characterised the flow state based on instantaneous velocity fields andthe time series of the superficial velocity. The velocity fields
were decomposed into symmetric and antisymmetric components with respect to the symmetries of the sphere pack. Basedon the temporal evolution of the kinetic energy of the antisymmetric components, it is possible to distinguish between laminar, transitional and turbulent flow. Finally, the velocity fields
were decomposed with respect to the average over the realisations that can be obtained by applying the symmetries of
the hexagonal sphere pack. For the turbulent flow case, a significant scale separation between the average and fluctuating
velocity fields is observed.

The flow state in oscillatory flow depends on two dimensionless parameters: the Hagen number Hg representing theamplitude and the Womersley number Wo representing the frequency of the forcing. The Reynolds number Re is a uniquefunction of the Hagen and Womersley number. Therefore, either of Hagen or Reynolds number can be used to characterise
the flow state. Based on a model equation, Gu & Wang (1991)
estimated which flow regime might be expected for a given
combination of Re and Wo.
In the present contribution, we investigate the onset of
turbulence in oscillatory flow through an idealised porous
medium geometry by means of direct numerical simulation. A
visual inspection of the instantaneous velocity fields suggeststhat the simulation database includes both laminar and turbulent simulations. However, it is unclear how to objectively
quantify the flow state. Commonly, turbulence is characterisedby temporal and spatial properties of the velocity fluctuations
about the mean, e.g. the frequency or wavenumber spectraor the spatial two-point correlations. In the present case, thedefinition of a fluctuation is not trivial since the flow is unsteady and inhomogeneous in all spatial directions. One could
perform a phase-dependent time average, but this comes at alarge computational expense since the flow needs to be integrated over many cycles. As our objective is to merely differentiate between laminar and turbulent flow, we consider thiscost to be disproportionate. On the other hand, especially athigher Womersley numbers analyses of the full velocity havethe problem that the mean flow does not vary slowly comparedto the turbulence and that the boundary and shear layers are notlarge compared to the turbulent scales. Therefore, the objective of the present work is to find quantitative evidence thatallows to discern the laminar and the turbulent flow state.

INTRODUCTION
Oscillatory porous media flow occurs in wave-induced
transport processes in coral reefs (Lowe et al., 2008) or in
the sediment bed; moreover, it could be of interest to enhance solute transport in chemical reactors (Crittenden et al.,2005). Other applications include for example wave-induced
flow through rubble-mound breakwaters (van Gent, 1993; Hallet al., 1995). In order to describe these transport processes, it
is essential to know whether the pore scale flow is turbulent,
as turbulence is associated with strong mixing.
Starting with Dybbs & Edwards (1984), the microscale
behaviour of stationary turbulent flow through regular andrandom sphere packs has been studied numerically (e.g. Hill
& Koch, 2002; He et al., 2019; Sakai & Manhart, 2020)and experimentally (e.g. Patil & Liburdy, 2012, 2013). Fourflow regimes are commonly distinguished (Dybbs & Edwards,
1984) depending on the Reynolds number: i) linear flow, ii)
steady nonlinear flow, iii) unsteady nonlinear flow and iv)
chaotic / turbulent flow. The transition to turbulent flow generally occurs for Reynolds numbers around 120 (Fand et al.,1987). In regular sphere packs, the emergence of chaotic andturbulent flow can be related to symmetry breaking bifurcations that occur around a Reynolds number of 100 (Hill &Koch, 2002) or for 138 < Re < 209 (Sakai & Manhart, 2020).
Experimental studies of turbulent oscillatory porous media flow have been performed with a focus on the behaviourof bulk flow quantities (van Gent, 1993; Hall et al., 1995;
Losada et al., 1995; Pamuk & Özdemir, 2014; Bağcı et al.,First, we investigate the behaviour of the streamwise
and cross-streamwise components of the superficial volumeaveraged velocity. Second, we analyse the breaking of thesymmetries that are imposed onto the flow by the sphere packgeometry. Third, we propose to average the flow over an ensemble of realisations that is generated using the symmetries

112th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
Table 1: Simulation parameters.

of the sphere pack geometry. This allows us to separate thevelocity fields into an average and fluctuating part based on
which we can distinguish laminar, transitional and turbulent
cases.

METHODOLOGY
Numerical method
We performed direct numerical simulation of the incompressible Navier-Stokes equations with our in-house code
MGLET that is based on Cartesian block-structured grids. The
spatial discretisation in MGLET uses a second-order central finite volume scheme with a staggered arrangement of variables
and the temporal discretisation employs a third-order explicit
Runge-Kutta method. The no-slip boundary conditions on thespheres is enforced by a discrete-forcing immersed boundary method that is described in Peller et al. (2006); Peller(2010). The momentum fluxes near the immersed boundary
are evaluated using a linear least-squares ghost-cell interpolation/extrapolation approach. The conservation of mass is enforced in every substep of the time integration scheme using aflux correction procedure in the interface cells and by solving
a Poisson equation for a correction pressure in the field.

WoRe

number of cycles

MF3

105.5

31.62

26.9

3†

MF4

106

31.62

74.0

4† / 1.32‡

MF5

106.5

31.62

157

6.3975

MF6

107

31.62

297

1.63

‡ This simulation was recomputed with a more finely resolved

triangle representation of the sphere pack geometry to match theother simulations.

The flow was computed with grid resolutions of 48, 96,
192 and 384 cells per diameter (cpd), resulting in a total number of 88 million fluid cells at 384 cpd. The error of the numerical solution was estimated based on the space-time L2 -normof the velocity and on the oscillation amplitude of the superficial velocity huis which is used to form the Reynolds number.
For all simulations the relative difference between the finest
and the second finest resolution is less than 1.3%.
The resolution requirements come from the wall boundary layers and possibly the turbulence in the bulk. The characteristic thickness
p of the oscillatory boundary layer can be
estimated as δ = 2ν/Ω which we resolve with 17 cells. On
the other hand, He et al. (2019) employed a grid resolution of
250 (cpd) to simulate stationary turbulent flow up to Re = 741.

As a porous medium, we choose a hexagonal closepacked arrangement of spheres of diameter d. It has a porosity of 0.259 which is the lower limit for packings of equal
spheres. The flow is described by the incompressible NavierStokes equations with an an oscillatory volume force fx sin Ωt
in the x–direction. Initially, the flow is at rest.
The problem is governed by two independent parameters:
the Hagen number Hg = fx d 3 /(ρν 2 ) represents the ratio of
the amplitude of the applied volumepforce to the viscous forces
and the Womersley number Wo = Ωd 2 /ν represents the ratio of the sphere diameter to the thickness of the oscillatory
Stokes boundary layer. The Reynolds number Re is basedon the sphere diameter d and the amplitude of the superficial
volume-averaged velocity after the decay of the transient. The
superficial velocity is defined as
1V

Hg† These simulations are presented in (Unglehrt & Manhart, 2022).

Study design

huis :=case

RESULTS
Instantaneous velocity fields
Figure 1a–d displays a section through the velocity field
perpendicular to the main flow direction for our √simulations.
As the Hagen number is increased by powers of 10, the velocity field develops pronounced flow structures in the large
pores, resembling those in transient flow described by Sakai &Manhart (2020). Due to the regular layout of the structures, wewould consider the cases MF3–MF5 as laminar. At the highest Reynolds number, the velocity field is asymmetrical andirregular vortical structures can be made out. Intuitively, we
would consider the case MF6 (figure 1d) as turbulent. In thefollowing, we discuss quantitative evidence that allows us toobjectively decide on the flow state.

Z

u dV
Vf

with the fluid volume Vf and the total volume V .
We consider flow at a Womersley number of 31.62 andfour values of the Hagen number. In linear flow, the selectedvalue of the Womersley number lies at the transition between
the low and the high frequency regime that are governed by
Stokes flow and potential flow, respectively. Consequently, oursimulations belong to the mid frequency (MF) regime. The
simulation parameters are reported in table 1.

Superficial velocity
In this section, we discuss the behaviour of the superficial
velocity as the flow state changes from laminar to (presumably) turbulent. The time series of the streamwise superficial
velocity (figure 1e–h) remains relatively smooth throughout
the changes in the velocity field. This is in contrast to the stationary case for which Hill & Koch (2002) and Sakai & Manhart (2020) reported (quasi-) periodic or chaotic oscillationsof the streamwise superficial velocity. At the lowest Reynolds
number, huis has a time lag with respect to the forcing and theamplitude is 25% below the steady state amplitude. As the
Reynolds number is increased, the phase lag reduces and thepeak amplitude approaches the steady state value that would be
attained for the same Hagen number. At the highest Reynolds
number, we observe a plateau rather than a pronounced peak.

Domain size and grid spacing
The simulation domain consists of two unit cells in thex–direction and one√unit cell√in the y- and z–direction. It thus
has an extent 2 d × 3 d × 2 6/3 d. For this size, the domain
contains two spheres in every lattice direction. He et al. (2019)used a domain of the same volume for their direct numericalsimulations of turbulent flow in a face-centred cubic spherepack. They state that ”the unit cell domain showed little variation in statistics compared to a larger domain”.

212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022(a) MF3 (Re = 26.9) †(b) MF4 (Re = 74.0) †(c) MF5 (Re = 157)(d) MF6 (Re = 297)(e) MF3 (Re = 26.9) †(f) MF4 (Re = 74.0) †(g) MF5 (Re = 157)(h) MF6 (Re = 297)

† simulations from (Unglehrt & Manhart, 2022)
Figure 1: (a)–(d) Velocity magnitude in the y-z-plane perpendicular to the flow direction. The velocity fields correspond to

the time of the maximum superficial velocity. (e)–(f) Forcing (dotted red line) and superficial velocity (blue line) over one
period. The superficial velocity is normalised with the steady state value (i.e. for a constant pressure gradient) obtained
from the amplitude of the forcing via the relations given by Sakai & Manhart (2020). The black symbol marks the timeof the snapshot presented in (a)–(d).
Quantification of symmetry breaking

The cross-streamwise components of the superficial velocity, hvis and hwis , are orthogonal to the imposed volume
force. For the cases MF3 and MF4, these components are verysmall as a consequence of the symmetries of the sphere packimprinted on the flow. On the other hand, for the cases MF5
and MF6 significant nonzero values of the cross-streamwise
components can be observed. Figure 2 shows a plot of hvis
and hwis for the latter simulations. The curves start out at theorigin and spread out into the plane. For the case MF5, the amplitude of hvis is approximately 30 times larger than the amplitude of hwis and approximately 400 times smaller than theamplitude of huis . A clustering of points at the origin indicatesthat the cross-streamwise components grow in time. For thecase MF6, the amplitude of hvis and hwis is approximately 60
times smaller than the amplitude of huis . Both cases shown infigure 2 exhibit irregular orbits, indicating chaotic behaviour.
Moreover, the cross-streamwise components of the superficial
velocity are only weakly correlated with the streamwise superficial velocity (represented by the colour of the curves).
While for MF6 the observed chaotic behaviour of thecross-streamwise components is in line with irregular structure of the velocity field shown in figure 1d, the visual impression of the velocity field of the case MF5 (figure 1c) does
not suggest any chaotic behaviour in this case. This could
be explained with the relatively low magnitude of the crossstreamwise components.
In conclusion, the streamwise superficial velocity changes
gradually with Reynolds number and we could not identify
features in the time series that could be used to clearly discern
laminar and turbulent flow. Suprisingly, for higher Reynolds
numbers the time and value of the maximum huis coincidewith the results of a quasi-steady approximation. On the other
hand, the cross-streamwise components are very small for thecases MF3 and MF4 whereas they are nonzero and exhibit
chaotic behaviour for MF5 and MF6.

In this section, we investigate the symmetry breaking of
the flow which was already indicated by the cross-streamwise
components of the superficial velocity.1 Based on the works of
Hill & Koch (2002) and Sakai & Manhart (2020), the symmetry breaking can be considered a prerequisite of turbulent flowin the sphere pack.
When a forcing is applied along the x–direction, laminarflow through the hexagonal sphere pack exhibits four symmetries: the flow is invariant with respect to
1. a translation by d in the x–direction (Tx )◦2. a translation by d at a 60
√ angle to the x–direction (Txy )3. a reflection around z = 6/3 d (S√z)√
4. a rotation by π about the axis y = 3/3, z = 6/6 (Rx )Figure 3 displays the result of these transformations applied tothe sphere pack inside the simulation domain. Considering theperiodic boundary conditions of the domain, the original andthe transformed configurations of the spheres are congruent.
In order to investigate the possible breaking of these symmetries, we decompose the instantaneous velocity fields into acomponent usym that is symmetric with respect to one of thesymmetries, and a corresponding antisymmetric component
uanti . The kinetic energy of the antisymmetric component,
hkanti is = 12 ρu2anti s , measures the violation of the symmetry
under consideration. Figure 4 displays the temporal evolutionof hkanti is for the various symmetries. For the cases MF3 andMF4 we observe that all symmetries are satisfied to an accuracy of more than 10−7 . The amplitude of the antisymmetriccomponents remains stationary. This indicates that the flow

1 Please note the following connection between the symmetries andthe cross-streamwise components of the superficial velocity: If theflow is symmetric under the rotation Rx , then hvis = hwis = 0 and if
the flow is symmetric under the reflection Sz , then hwis = 0.

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

(a) Translation Tx(b) Translation Txy(c) Reflection Sz(d) Two-fold rotation Rx(a) MF5 (Re = 157)

Figure 3: Symmetries of laminar flow in a hexagonal

sphere pack due to a volume force along the x–direction.
The black box represents the position of the sphere packbefore the application of the symmetry operation.
(b) MF6 (Re = 297)
ble or chaotic dynamics at Re = 157 and Re = 297. A direct inspection of the spatial distribution of the antisymmetriccomponents is however not informative as in general, turbulent motion is neither symmetric nor antisymmetric. On theother hand, the mean velocity field of a turbulent flow is usually symmetric. Therefore, we decompose the velocity field
into a part that satisfies all of the symmetries

Figure 2: Cross-streamwise components of the superfi-

cial velocity; the curves are coloured by the streamwise
component of the superficial velocity. The irregular orbits indicate chaotic flow.

is stable with respect to the symmetry breaking and supports
the hypothesis that the flow is laminar in these cases. For thecase MF5 a simultaneous exponential growth of the energy of
all antisymmetric components can be observed. After approximately 5 cycles, the kinetic energy of the Rx -antisymmetricpart saturates at 0.3% of the maximum total kinetic energy (cf.
figure 4d). This suggests that (i) there exists a linear instability mechanism that facilitates the growth of the antisymmetricpart and (ii) a nonlinear self-interaction of the instability occurs which limits its growth. However, this does not necessarily imply that the flow becomes turbulent. Finally, for thecase MF6 a fast growth of the antisymmetric components canbe observed. The kinetic energy of the antisymmetric components with respect to Rx , Sz and Txy peaks between 9 and
10% of the cycle maximum of the total kinetic energy.
The initial values of hkanti is turned out to depend on theaccuracy of the triangle representation of the spheres, fromwhich the interpolation stencils and coefficients in the immersed boundary method are generated. There exist twosources of the initial symmetry breaking perturbations in ourcode: asymmetric interpolants at the immersed boundary andthe residual of the pressure correction. A visual inspection of
the antisymmetric parts of the velocity field shows that only
the physically meaningful parts of these perturbations are amplified.

u=
1(I + Rx ) (I + Tx ) I + Txy (I + Sz ) u
16(1)

where I is the identity, and a residual component u0 = u − uwhich is in general neither symmetric nor antisymmetric withrespect to any symmetry.
Alternatively, this decomposition can be interpreted as anaverage over the ensemble of 16 velocity fields that is generated by the action of the symmetry group. Such a symmetry
average was used by Sirovich (1987); Sirovich & Park (1990)
to enlarge the sample size for a proper orthogonal decomposition. Please note that while equation (1) can be understood as
an ensemble average, it may not converge to the same result asa time or phase average. For example, Srikanth et al. (2021)
investigated a persistent symmetry breaking phenomenon in a
porous medium for which the time-averaged flow is asymmetrical. As there is an equal probability for the flow to deviate toeither side, the ensemble mean flow would be symmetrical.
The symmetry group average (1) underestimates the energy of the turbulent fluctuations compared to the ensemble
mean. This can be demonstrated as follows. As the flow
problem is symmetric, the ensemble mean E[u] is symmetric and invariant with respect to Rx , Tx , Txy and Sz . For
example, Rx E[u] = E[u] and E[Rx u] = E[u]. Hence, the ensemble mean does not change under the symmetry group average, i.e. E[u] = E[u] and E[u] = E[u], and the expectation of
the fluctuation about the symmetry group average is zero, i.e.
E[u0 ] = E[u − u] = 0. Therefore u0 is a part of the fluctuation
about the ensemble mean.

Symmetry group averaging
In the preceding sections, it was demonstrated that theantisymmetric components of the velocity field show unsta-

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
Figure 5 displays the spatial distribution of the magnitudeof u and u0 . For the case MF5 the symmetry group averaged
velocity field u looks identical to the instantaneous field in figure 1c. This is due to the low energy of the antisymmetriccomponents (cf. figure 4). The residual component u0 features
symmetric patterns that are arranged near the high-velocity
features in u. For the case MF6, the symmetry group averaged velocity field contains symmetric versions of the highintensity regions near the contact points and inside the large
pores that can be identified in figure 1d. The residual velocity
field u0 does not contain these features; instead many irregularand small scale vortices can be identified.
(a) Translation Tx

CONCLUSION
We performed direct numerical simulation of oscillatory
flow through a hexagonal sphere pack at Wo = 31.62 and four
Reynolds numbers. We investigated instantaneous velocity
fields, the time series of the streamwise and cross-streamwise
superficial velocity and the breaking of the symmetries of theflow. Moreover, the velocity field was averaged over the ensemble of realisations generated from the instantaneous velocity fields by applying the symmetries of the sphere pack,and the spatial distributions of the average and fluctuation field
were discussed. The resulting average field approximates theensemble mean of the flow and the resulting fluctuations are asubset of the fluctuations about the ensemble mean.
At Re = 26.9 and Re = 74.0 the cross-streamwise components of the superficial velocity are nearly zero and the flow
shares all symmetries of the sphere pack. Consequently, thesecases can be considered as laminar. On the other hand, thesymmetries of the flow are broken for Re = 157 and Re = 297
and the cross-streamwise superficial velocity show nonzero
values and chaotic behaviour. At Re = 157 the antisymmetric parts of the velocity field with respect to the symmetriesof the sphere pack grow exponentially in time, indicating thata linear instability is present in the flow. Eventually, the antisymmetric part with respect to the rotation symmetry (cf.
figure 3d) saturates at 0.3% of the total kinetic energy and asimilar behaviour could be expected for the other symmetries.
The fluctuations with respect to the symmetry group average
have an ordered appearance, corroborating that this flow caseis not yet turbulent. Consequently, we categorise it as a transitional flow. Finally, at Re = 297 the symmetries are broken
rapidly and the antisymmetric parts of the velocity field peakat approximately 10% of the total kinetic energy. The fluctuations about the symmetry group average consist of disorderd
vortical structures and a scale separation can be observed between the fluctuating and the average field. Together with the
chaotic behaviour of the cross-streamwise superficial velocity,
this strongly supports the view that this case exhibits turbulent
flow.
Further investigations may be necessary to assess the effect of the domain size and the dependency of the results on
the Womersley number. Moreover, our findings for the case at
Re = 157 suggest that a Floquet-type linear instability analysisof the flow could be interesting.

(b) Translation Txy(c) Reflection Sz(d) Two-fold rotation Rx
Figure 4: Volume-averaged kinetic energy of the an-tisymmetric components.
MF3 (Re = 26.9),
MF4 (Re = 74.0),
MF5 (Re = 157), andMF6(Re = 297). The kinetic energy hkanti is is normalisedwith the maximum total kinetic energy of the last cycle.

ACKNOWLEDGEMENTS
The authors gratefully acknowledge the financial supportof the DFG under grant no. MA2062/13-1. Computing time
was granted by the Leibniz Supercomputing Centre on its
Linux-Cluster.

5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022(a) |u| for MF5 (Re = 157)(b) 30 |u0 | for MF5 (Re = 157)(c) |u| for MF6 (Re = 297)(d) |u0 | for MF6 (Re = 297)

Figure 5: Magnitude of the symmetry group average velocity u and the corresponding fluctuation u0 in the y-z-plane

perpendicular to the flow direction. The velocity fields correspond to figure 1c–d.
REFERENCES

ties with and without waves: A synthesis of porous media
and canopy flow approaches. Limnology and Oceanography
53 (6), 2668–2680.
Pamuk, M. T. & Özdemir, M. 2014 Hydrodynamic Aspects ofOscillating Flow Through Porous Media Consisting of Steel
Spheres. In Progress in Exergy, Energy, and the Environment (ed. I. Dincer, A. Midilli & H. Kucuk), pp. 633–640.
Cham: Springer International Publishing.
Patil, V. A. & Liburdy, J. A. 2012 Flow Structure Identificationin Unsteady Flow in Porous Media. In Proceedings of theAsme Fluids Engineering Division Summer Meeting, 2012,
Vol 1, Pts a and B, Symposia, pp. 775–781. New York:Amer Soc Mechanical Engineers.
Patil, V. A. & Liburdy, J. A. 2013 Flow structures and theircontribution to turbulent dispersion in a randomly packed
porous bed based on particle image velocimetry measurements. Physics of Fluids 25 (11), 113303.
Peller, N. 2010 Numerische Simulation turbulenter Strömungen mit Immersed Boundaries. PhD thesis, Technische Universität München, München.
Peller, N., Le Duc, A., Tremblay, F. & Manhart, M. 2006 Highorder stable interpolations for immersed boundary methods. International Journal for Numerical Methods in Fluids
52 (11), 1175–1193.
Sakai, Y. & Manhart, M. 2020 Consistent Flow Structure Evolution in Accelerating Flow Through Hexagonal Sphere
Pack. Flow, Turbulence and Combustion 105 (2), 581–606.
Shigematsu, T., Nakajo, S. & Okada, Y. 2018 An experimental
study on the interaction between oscillatory flow and idealized porous bed. Journal of Coastal Research pp. 981–985.
Sirovich, L. 1987 Turbulence and the Dynamics of Coherent
Structures Part II: Symmetries and Transformations. Quarterly of Applied Mathematics XLV, 573–82.
Sirovich, L. & Park, H. 1990 Turbulent thermal convection ina finite domain: Part I. Theory. Physics of Fluids A: Fluid
Dynamics 2 (9), 1649–1658.
Srikanth, V., Huang, C.-W., Su, Timothy S. & Kuznetsov, A. V.
2021 Symmetry breaking of turbulent flow in porous media
composed of periodically arranged solid obstacles. Journalof Fluid Mechanics 929.
Unglehrt, L. & Manhart, M. 2022 Onset of nonlinearity in oscillatory flow through a hexagonal sphere pack. Submitted
to Journal of Fluid Mechanics.
van Gent, M. R. A. 1993 Stationary and oscillatory flow
through coarse porous media. Communications on hydraulic and geotechnical engineering, No. 1993-09 .

Bağcı, Ö., Arbak, A., De Paepe, M. & Dukhan, N. 2016 Experimental results for oscillatory water flow in 10-ppi metal
foam at low-frequencies. Journal of Physics: Conference
Series 745, 032040.
Crittenden, B. D., Lau, A., Brinkmann, T. & Field, R. W.
2005 Oscillatory flow and axial dispersion in packed bedsof spheres. Chemical Engineering Science 60 (1), 111–122.
Dybbs, A. & Edwards, R. V. 1984 A New Look at Porous Media Fluid Mechanics — Darcy to Turbulent. In Fundamentals of Transport Phenomena in Porous Media (ed. J. Bear
& M. Y. Corapcioglu), NATO ASI Series, vol. 82, pp. 199–
256. Dordrecht: Springer Netherlands.
Fand, R. M., Kim, B. Y. K., Lam, A. C. C. & Phan, R. T.
1987 Resistance to the Flow of Fluids Through Simple andComplex Porous Media Whose Matrices Are Composed ofRandomly Packed Spheres. Journal of Fluids Engineering
109 (3), 268–273.
Graham, D. R. & Higdon, J. J. L. 2002 Oscillatory forcing offlow through porous media. Part 2. Unsteady flow. Journalof Fluid Mechanics 465.
Gu, Z. & Wang, H. 1991 Gravity waves over porous bottoms.
Coastal Engineering 15 (5-6), 497–524.
Hall, K. R., Smith, G. M. & Turcke, D. J. 1995 Comparison of oscillatory and stationary flow through porous media.
Coastal Engineering 24 (3-4), 217–232.
He, X., Apte, S. V., Finn, J. R. & Wood, B. D. 2019 Characteristics of turbulence in a face-centred cubic porous unit cell.
Journal of Fluid Mechanics 873, 608–645.
Hill, R. J. & Koch, D. L. 2002 The transition from steady to
weakly turbulent flow in a close-packed ordered array ofspheres. Journal of Fluid Mechanics 465, 59–97.
Iervolino, M., Manna, M. & Vacca, A. 2010 Pulsating Flowthrough Porous Media. In Turbulence and Interactions(ed. E. H. Hirschel, W. Schröder, K. Fujii, W. Haase,
B. Leer, M. A. Leschziner, M. Pandolfi, J. Periaux, A. Rizzi,
B. Roux, Y. I. Shokin, M. Deville, T.-H. Lê & P. Sagaut), ,vol. 110, pp. 167–173. Berlin, Heidelberg: Springer Berlin
Heidelberg.
Kardgar, A. & Jafarian, A. 2021 Numerical simulation of turbulent oscillating flow in porous media. Scientia Iranica
28 (2), 743–756.
Losada, I. J., Losada, M. A. & Martı́n, F. L. 1995 Experimental
study of wave-induced flow in a porous structure. Coastal
Engineering 26 (1-2), 77–98.
Lowe, R. J., Shavit, U., Falter, James L., Koseff, J. R. &Monismith, S. G. 2008 Modeling flow in coral communi-

6