Yoshiyuki Sakai TUM School of Engineering and Design Technical University of Munich Arcisstr. 21, 80333 Munich, Germany yoshiyuki.sakai@tum.de Michael Manhart TUM School of Engineering and Design Technical University of Munich Arcisstr. 21, 80333 Munich, Germany michael.manhart@tum.de
ABSTRACT
We analysed a high-fidelity DNS dataset of accelerating flow through a porous medium consisting of uniform spheres in hexagonal close-packed arrangement. The fluid was accel- erated from rest by a constant pressure gradient until the super- ficial volume-averaged velocity attained its steady-state ampli- tude.  Our focus is on the generation mechanism of the over- shoot of the superficial streamwise velocity that appears just before the steady-state amplitude is reached. The relative mag- nitude of the velocity overshoot appears to be Reynolds num- ber dependent, and can be larger than 5% of the steady-state value. The current study suggests that the velocity overshoot is a consequence of enhanced intercomponent energy transfer carried out by pressure, which redistributes streamwise energy into the orthogonal directions, and it mainly takes place in the wake regions behind the sphere contact points.
BACKGROUND
Transient porous media flow establishes an important class of flow phenomena. There are numerous technical and environmental flow problems which fall into this class, e.g.: flow through plant canopies in atmospheric boundary layers; coral colonies under the influence of ocean waves (Lowe et al., 2005, 2008); or the interaction of a turbulent boundary layer with the flow and transport in the upper layer of the soil or a snow layer. In our recent contributions, we studied accelerating flow through hexagonal sphere-pack porous media by means of direct numerical simulations (Sakai & Manhart, 2020; Sakai et al., 2021). Figure 1 depicts the sphere-pack geometry. We observed that during the development, superficial (i.e. volume-averaged including solid phase) streamwise velocity (u) ς overshoots before reaching to the steady-state value when the Reynolds number based on the steady-state ( u ) s , sphere di- ameter d and the fluid kinematic viscosity ν —which will be referred to as Re steady — is larger than 10 (see figure 2). The relative intensity of the overshoot with respect to the steady- state level depends on Resteady , and it could be larger than 5%, therefore this phenomenon is technically relevant (see fig- ure 3). Conversely if the Reynolds number is lower, { u } s re- laxes to its steady-state without any overshooting (see the L6 case in figure 2).  A similar observation was made by Hill & Koch (2002) with a different sphere pack geometry, however the underlining mechanics was remained unexplored. There- fore, it is of our interest to study such mechanism particu- larly because the commonly-used unsteady flow models (e.g. unsteady Darcy, Darcy-Forshheimer equations) cannot predict the overshoot. Moreover, it was also demonstrated that the ex- act thin boundary layer asymptotics for the porous media flow by Johnson et al. (1987) does not incorporate such overshoot (Sakai et al., 2021). In Sakai et al. (2021), it was found that the the appearance time of the velocity overshoot asymptotically converges to a constant multiple of the inviscid dimensionless time τ inv : 1 Figure 2:  Streamwise superficial velocity development, normalised by the inviscid scales. Colored symbols im- which is based on the intrinsic pressure gradient (∂ x ( p ) ; ) driv- ing the flow, the sphere diameter d and the fluid density ρ . From the flow topological viewpoint, this inviscid time-scale was shown to control the wake formation process inside the accelerating porous media flow, and more specifically, the de- velopment of the flow separations emerging behind the sphere contact points. In the current contribution, we investigate the generation mechanism of the velocity overshoot in the accelerating porous media flow by means of the volume-averaged energy equation. The a-posteriori application of such energy equation to our DNS dataset (Sakai & Manhart, 2020) allows us to quantify the exchange of energy between the three spatial directions. Finally, our special attention will be paid upon the possible link between the velocity overshoot and the pore-scale wake formation, which appear to be controlled by the same univer- sal time-scale.
NUMERICAL METHODS
We employ our in-house flow solver MGLET, which is based on a 2nd-order central finite-volume spatial discretisa- tion, and a low-storage variant of a three-step Runge-Kutta method (Williamson, 1980), to integrate in time the incom- pressible Navier-Stokes equations for velocity and pressure, ViZ.
where u = [ u , v , w ] is velocity vector, p , and ρ are pressure, fluid density respectively. Those variables are stored in a Cartesian grid with stag- gered arrangement.  Complex geometries are represented by a second-order mass-conserving immersed boundary method (Peller et al. , 2006; Peller, 2010; Sakai & Manhart, 2020). The pressure computation is decoupled from the velocity computa- tion by Chorin's projection method (Chorin, 1968).  Conse- quently, a Poisson equation needs to be solved for the pressure for each Runge-Kutta sub-step by means of Strongly Implicit Procedure (SIP) (Stone, 1968). The code is MPI-parallel based on a conventional domain decomposition approach, and capable of a satisfactory scaling up to O (10 ² ) parallel processes, thanks to the continuous per- formance optimisation being applied to the code (Sakai et al. , 2019; Sakai & Manhart, 2021).
DNS DATASET
In this contribution, we analyse a high-fidelity dataset consisting of 16 DNS cases, which was first published in Sakai & Manhart (2020). We consider porous media flow that is rep- resented by a triply-periodic numerical box of [ L x , L y , L z ] /d = [2, √3, 2 √ 6], which is filled by uniform spheres being packed in a hexagonal arrangement (cf.  figure 1).  The porosity of the sphere-pack is at the close-pack limit of ε ≈ 0.26. A rest- ing fluid was accelerated at a constant streamwise pressure gradient in x-direction (∂x(p);) until steady-state being es- tablished. The applied pressure gradient, dynamic viscosity (μ = ρν), and the sphere diameter ( d ) are fixed, whereas ρ is adjusted to achieve desired Reynolds numbers.  We employed equidistant grid points distribution in x , y , z -directions, with 162 finite-volume cells per diameter up to Resteady = 48, and 324 points per diameter between 59 ≤ Resteady ≤ 347.  This highly adequate grid resolution was adopted to resolve the very thin potential-flow type boundary layers forming on the sphere surface when the flow is suddenly accelerated from the rest. The range of Resteady covered by the current dataset is 5.62 × 10 −7 ≤ Resteady ≤ 347, spanning over linear, steady non- linear, unsteady nonlinear and turbulent flow regimes. In this proceeding, we study closely the UNL2 (Unsteady Non-Linear) and the T4 (Turbulent) cases from the dataset, which correspond to Resseady = 138 and 347 respectively. Those cases are well above the critical point where the velocity overshoot appears (see figure 3), yet their Reynolds numbers being sufficiently separated so that the generality of our find- ings could be tested to the limit of the database.
SUPERFICIAL VELOCITY EVOLUTION
Figure 4 shows the temporal evolution of superficial ve- locity components (solid lines), as well as the square-root of the energy components (dashed lines). Both in the case UNL2 and T4, { v } ç and { w } ç are few order of magnitude smaller than the predominant component, {u} s . This is due to the reflec- tional and the rotational symmetries of the pore geometry be- ing imprinted to the ν and w fields which cancel the most of their contributions under the volume-averaging operation, al- though those symmetries are weakly broken already at those Reynolds numbers (Sakai & Manhart, 2020).  Note when the flow Reynolds number is smaller than unity (i.e.  linear flow regime), the above symmetries are fully preserved so that { v }。 and {w}。are zero to the limit of the numerical accuracy. Con- versely, ⟨ u ⟩ ς is finite for all Re steady although there also ex- ists a d -periodic translational symmetry in x -direction. This is the consequence of the pressure gradient being applied in the streamwise direction destroying this symmetry. In Sakai & Manhart (2020) it was shown that reverse flow emerges for Resteady > 10. Consequently, the flow recircula- tion regions contribute negatively to {u} ς , whilst it does not for the corresponding root-mean-square value  ( uu ) s . This dif- ference is the reason why the level of  ( uu ) s is considerably higher than {u}s (see figure 4). For a similar reason, √ {vv}s and √{ww} s are also significantly larger than {v} s and {w} s . Through the application of the intrinsic pressure gradi- ent in x-direction the streamwise kinetic energy grows rapidly and peaks out at around t/ τ inv = 1.5 in both cases. Subse- quently, the energy components in the perpendicular directions (i.e. {vv} s and {ww} s } peak out around t/ t inv = 1.7, with rela- tively steep gains between 1.0 < t / Tinv < 1.7. Whilst the { vv } s and 〈 ww 〉 s are climbing up to their peaks, 〈 uu 〉 s drops its level noticeably. Based on the above observation, we hypothesise that an intercomponent energy transfer from the primary flow direc- tion to the directions perpendicular to the flow is responsible for the appearance of the superficial velocity overshoot. In the following section, the above hypothesis is evaluated by means of the superficial energy budget equations.
SUPERFICIAL ENERGY BUDGET
Let us suppose the streamwise superficial velocity over- shoot as a consequence of the energy transfer to other spatial directions, or dissipation, or in combination. Consequently, we consider the following component-wise superficially-averaged energy transport equations: where p ′ is the fluctuations around ⟨ p ⟩ i, i.e. p = ⟨ p ⟩ i + p ′ , whereas the subscript ( i ) implies i being a free index, i.e. no summation over the index. Note that the above transport equa- tions are the component-wise variant of the one discussed in Zhu et al. (2014), which consists of advective , pressure , energy input and viscous transport terms (see annotations in equation 4). Note that the viscous term consists of energy diffusion and dissipation, whereas the energy input term is only active in x -direction in the current configuration.  As in Zhu et al. (2014), we expect the following by definition: the advective terms should vanish under periodic boundary and divergence- free conditions; similarly the viscous diffusion vanishes due to the periodicity, whereas viscous dissipation is active as a sink term. Finally the pressure term can be further decomposed into pressure diffusion and inter-component transfer terms, as:
inter-comp (5) Note that the pressure diffusion term vanishes also due to the periodic boundary condition, whereas the inter-component term acts as both source and sink terms that allow energy in different spatial directions to be exchanged and re-distributed. Those superficial energy transport terms were evaluated a-posteriori based on a series of flow snapshots: for the UNL2 case, a series of 61 snapshots were saved which were sepa- rated by 쓰 ≈ 0.04; whereas for the T4 case we collected 181 snapshots with energy transport terms. The kinetic energy is injected through the pressure input term (the dotted line). The advective terms (the solid lines) are virtually zero due to the periodicity, which is in accordance with our expectation. Conversely, finite lev- els of the viscous terms (the solid-dashed lines) can be ob- served, which should be exclusively from the energy dissipa- tion process in the respective components. Similarly, the pres- sure terms (the dashed lines) consist exclusively of the inter- component energy terms, transferring energy from the primary direction x to the orthogonal directions in y and z, according to their respective signs. Since any disturbance in the pressure field propagates at the infinite speed in incompressible flow, the pressure trans- port terms in the respective directions peak out simultane- ously, which is slightly later than t/τ inv = 1.5 for both cases. The magnitude of the pressure terms increases rapidly between t / Tinv = 1.2 and 1.6, and this time range coincides with when the pore-scale vortical structures break down (see Fig. 9 and 12 in Sakai & Manhart (2020)), directly after the Re-independent consistent pore-scale flow structure evolution emerges. More- over, it coincides with the aforementioned time range when we observed rapid increases in {vv} s and {ww} s (cf. figure 4). Finally, notice a delay between the growth of the pressure re-distribution and the viscous dissipation terms.  This could be attributed to the time required for the aforementioned fine- scale structures in the bulk flow region to build up, which con- tribute to the enhanced energy dissipation. In the subsequent section, we further investigate the mechanism of the pressure inter-component energy transfer by evaluating the spatial distributions of the terms inside the pores.
RELEVANT PORE-SCALE FLOW
To examine the mechanism of the rapid growth of the pressure transport terms, their pore-scale spatial distributions in the UNL2 case are depicted at two instances: t / Tiny = 0.6 in figure 6, and t/τ inv = 1.5 in figure 7. Note that the earlier instance (t/Tinv = 0.6) is in the middle of the (u)s accelera- tion, whereas the later instance ( t / τ inv = 1.5) is approximately at tpeak- A local symmetry plane defined byֲlected for this visual study. The outline of the pore geometry is defined by four spheres at the corners, in addition to a contact point of a pair of spheres being aligned in the direction perpen- dicular to the depicted plane.  In figure 6 and figure 7, such con- tact point is visible at the centre of the plane. The primary flow direction is from left to right. Consequently, the flow first en- ters the larger pore, before it is split into two high-momentum jets by the contact point.  Those split jets flow through very narrow gaps besides the contact point, then are merged inside the larger pore space behind the contact point.  Note that this split-merge sequence repeats itself due to the periodic bound- ary condition. Finally, only one half of the simulation domain in x-direction is shown in those figures to improve the visibil- ity. At t/Tiny = 0.6 (cf. figure 6), the flow is entirely attached to the sphere surfaces (not shown here) and only minor en- ergy transfer exists, which is necessary to guide the attached flow through the pore geometry.  It is evident that the energy is drawn from the streamwise to the perpendicular directions directly in front of, as well as behind, the sphere contact point. Conversely, the energy is fed back to x-direction on either side of the contact point where a pair of very narrow pores is lo- cated. Under this condition, the spatial distribution of the pres- sure transport terms appear to be quite antisymmetric in x- direction (i.e.  left-right), in addition to the reflectional sym- metry around the horizontal line through the touching point (top-bottom). This apparent antisymmetry is only slightly bro- ken in a way that: there exists a stagnation point only in front of the contact point; the local extrema are moderately shifted downstream; and the amplitude of the peaks are slightly larger in the upstream of the contact point than the downstream in all spatial components. This quasi-antisymmetric distribution of the energy transfer terms means that a large portion of the contribution is cancelled out under the volume-averaging op- eration. Figure 7 shows the corresponding fields from t / Tinv = 1.5. At this instance, a massive flow separation and recirculation have emerged in the downstream of the contact point (again not shown here). Since the inertial influence is so high that the high-momentum jets coming out of the narrow gaps besides the contact point can no longer trace the pore geometry and be merged in the downstream.  Instead, the jets remain being separated and impinge on to the sphere surfaces on the side, where the strongest level of negative energy transfer can be observed in the x-component (cf. figure 7a). As the consequence of this flow regime change, the afore- mentioned antisymmetry in the spatial distribution of the pres- sure transport terms has been destroyed entirely.  This com- plete disappearance of the antisymmetry results a rapid in- crease in the net energy transfer rate from the streamwise to the perpendicular directions. Evidently, the energy level in the streamwise direction drops due to the enhanced re-distribution, and the difference between the peak to the subsequent equilib- rium levels forms the aforementioned velocity overshoot. 12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12) Osaka, Japan (Online), July 19–22, 2022
SUMMARY
In this contribution, we consider accelerating porous me- dia flow by means of fully-resolved DNS simulations, and we investigate the mechanism of the velocity overshoot observed in the superficial velocity development for sufficiently high Reynolds number.  For all Re , the intrinsic pressure gradient feeds energy into the streamwise component, whereas only un- der a strong influence of inertia —i.e. high enough Re when the velocity overshoot emerges— the net inter-component en- ergy transfer from the streamwise into the perpendicular com- ponents becomes significant.  Through the detailed exam- ination inside the pores, it was demonstrated that the en- hanced inter-component energy transfer is the consequence of flow separation behind the contact points.  Finally, this rapid pressure-induced energy transfer takes place between t / Tiny = 1.0 and 1.7, where Tiny is inviscid time-scale.