Decomposition of the drag force in steady and oscillatory flow through a hexagonal sphere pack

Abstract

We investigate steady and oscillatory flow through a hexagonal close-packed arrangement of spheres in the framework of the volume-averaged momentum equation. We quantify the friction and pressure drag based on a direct numerical simulation dataset. Using the pressure decomposition of Graham (J. Fluid Mech., vol. 881, 2019), the pressure drag can be further split up into an accelerative, a viscous and a convective contribution. For the accelerative pressure, a closed-form expression can be given in terms of the potential flow solution. We investigate the contributions of the different drag components to the volume-averaged momentum budget and their Reynolds number scaling. For steady flow, we find that the friction and viscous pressure drag are proportional to $Re$ at low Reynolds numbers and scale with $Re^{1.4}$ for high Reynolds numbers. This is close to the steady laminar boundary layer scaling. For the convective pressure drag, we find a cubic scaling at low and a quadratic scaling at high Reynolds numbers. The Reynolds stresses have a minor contribution to the momentum budget. For oscillatory flow at low and medium Womersley numbers, the amplitudes of the drag components are similar to the steady cases at the same Reynolds number. At high Womersley numbers, the drag components behave quite differently and the friction and viscous pressure drag are relatively insignificant. The drag components are not in phase with the forcing and the superficial velocity; the phase lag increases with the Womersley number. This suggests that new models beyond the current quasisteady approaches need to be developed.

1. Introduction

In this contribution we investigate the behaviour of the drag in steady and oscillatory flow through a hexagonal sphere pack. The sphere pack can be decomposed into triply periodic unit cells of size $\ell$ and is assumed to have a spatial extent $L \gg \ell$ (figure 1a). When the flow is driven by pressure or velocity variations on the macroscale, for example by a pressure wave of wavelength ${O}(L)$, the flow is locally almost periodic (Ene & Sanchez-Palencia 1975). The pore-scale flow is described in terms of the velocity $\boldsymbol {u}$ and the pressure $P$. The large-scale flow is described in terms of the macroscopic pressure gradient $\boldsymbol {f}$, which is defined such that the pore-scale pressure deviation $p=P-\boldsymbol {f}\boldsymbol {\cdot }\boldsymbol {x}$ is a periodic function, and in terms of the superficial velocity $\left \langle \boldsymbol {u}\right \rangle _{s}$, which is defined as the volume average of the pore-scale velocity $\boldsymbol {u}$ over the unit cell

(1.1)\begin{equation} \left\langle\boldsymbol{u}\right\rangle_{s} = \frac{1}{V} \int_{V_{f}} \boldsymbol{u}\,\mathrm{d}V . \end{equation}

Here, $V$ is the volume of the unit cell and $V_{f}$ is the fluid volume within the unit cell; the porosity $\epsilon$ is defined as the ratio $V_{f}/V$. Note that depending on the flow regime, the unit cell has to be chosen larger than the primitive unit cell of the geometry (Agnaou, Lasseux & Ahmadi 2016). Here, the unit cell contains four primitive unit cells (figure 1b). In the limit $\ell /L\to 0$, the superficial velocity is governed by the continuity equation

(1.2)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left\langle\boldsymbol{u}\right\rangle_{s} = 0 \end{equation}

and a local relation between the superficial velocity and the macroscopic pressure gradient that follows from the solution of the Navier–Stokes equations on the unit cell (Ene & Sanchez-Palencia 1975). The pore-scale velocity $\boldsymbol {u}$ and the pore-scale pressure deviation $p$ are regarded as triply periodic fields on the unit cell. The relation between the superficial velocity and the macroscopic pressure gradient can also be expressed by the volume-averaged Navier–Stokes equations, which are obtained by averaging equation (2.1b) over the unit cell. By Gauss’ theorem and the periodic boundary conditions, the integrals over the open pore areas cancel; thus we obtain

(1.3)\begin{equation} \rho\frac{\partial\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\partial t}=\underbrace{-\frac{1}{V}\int_{A_{fs}} p\,\boldsymbol{n}\,\mathrm{d}A}_{\textit{pressure drag}} \underbrace{-\frac{1}{V} \int_{A_{fs}} \boldsymbol{\tau}_{w}\,\mathrm{d}A}_{\textit{friction drag}} - \epsilon\,\boldsymbol{f} . \end{equation}

In this equation, the symbol $\boldsymbol {\tau }_{w}=\mu \,(\boldsymbol {\nabla }\otimes \boldsymbol {u})^{\mathrm {T}}\vert _{w}\boldsymbol {\cdot }\boldsymbol {n}$ represents the wall shear stress vector, and $A_{fs}$ denotes the fluid–solid interface, i.e. the surface of the spheres. Note that the force exerted by the fluid onto the spheres also contains a contribution from the macroscopic pressure gradient. While (1.2) and (1.3) have been derived assuming a periodic porous medium, they can also be obtained for non-periodic porous media by the volume-averaging theory of Whitaker (1986, 1996) if the pore scale, the averaging scale and the macroscale are sufficiently separated. A comparison between the homogenisation approach outlined above and the volume-averaging approach can be found in Davit et al. (2013). The pressure drag and the friction drag terms appearing in (1.3) are unclosed with respect to $\left \langle \boldsymbol {u}\right \rangle _{s}$ and $\boldsymbol {f}$. In general, they can be obtained only by direct numerical simulation (DNS) of the pore-scale flow. The aim of modelling is to replace the solution to the pore-scale flow problem by an explicit relationship between $\boldsymbol {f}$ and $\left \langle \boldsymbol {u}\right \rangle _{s}$.

Figure 1. Conceptual sketch of the volume approach for the hexagonal sphere pack. (a) The sphere pack consists of triply periodic unit cells. The periodic flow $\boldsymbol {u}$, $p$ inside the unit cell is driven by the macroscopic pressure gradient $\boldsymbol {f}$. (b) The simulation domain of the hexagonal sphere pack consists of four primitive unit cells (one of which is highlighted in yellow).

In this work, we investigate this relationship and consider the macroscopic pressure gradient as a known quantity. The pore-scale flow is computed numerically in a triply periodic domain of the hexagonal sphere pack (shown in figure 1b) for a constant and a sinusoidally oscillating forcing. The pore-scale flow then depends on two dimensionless numbers that are formed with the sphere diameter $d$, the density $\rho$ and the kinematic viscosity $\nu$: The Hagen number $Hg=\vert \boldsymbol {f}\vert d^3/(\rho \nu ^2)$ describes the magnitude of the macroscopic pressure gradient relative to the viscous forces. In other work, the Hagen number is referred to as the pressure-gradient-based Reynolds number (Ene & Sanchez-Palencia 1975; Firdaouss, Guermond & Le Quéré 1997; Iervolino, Manna & Vacca 2010; Lasseux, Valdés-Parada & Bellet 2019). The Womersley number is defined as $Wo=\sqrt {\varOmega d^2/\nu }$ where $\varOmega$ is the angular frequency of the forcing. The Womersley number is proportional to the ratio of the sphere diameter to the Stokes layer thickness $\sqrt {2\nu /\varOmega }$ and thus determines the region that is affected by the wall friction via diffusive transport. For each parameter combination, a Reynolds number $Re=\vert \left \langle \boldsymbol {u}\right \rangle _{s}\vert d/\nu$ for steady flow or

(1.4)\begin{equation} Re=\limsup_{t\to\infty}\frac{\left\vert\left\langle\boldsymbol{u}\right\rangle_{s}\right\vert d}{\nu} \end{equation}

for oscillatory flow results from solving the pore-scale flow problem.

Next, we briefly summarise some important findings on the flow resistance behaviour of steady flow. In this case there are closed-form expressions which allow the pore-scale problem to be bypassed. Figure 2 shows the drag coefficient defined according to (Macdonald et al. 1979) as

(1.5)\begin{equation} F_k' = \frac{\left\vert\,\boldsymbol{f}\right\vert d}{\rho\!\left\langle u\right\rangle_{s}^2} = \frac{Hg}{Re^2} \end{equation}

as a function of the Reynolds number for the DNSs of Sakai & Manhart (2020) of steady flow through the hexagonal sphere pack. For very small Reynolds numbers, the superficial velocity depends linearly on the macroscopic pressure gradient. This relationship is described by Darcy's law (dotted line) which can be written in dimensionless form as

(1.6)\begin{equation} Hg = \frac{d^2}{K} Re \quad\text{or}\quad F_k' = \frac{d^2}{K} Re^{{-}1}, \end{equation}

where $K$ denotes the permeability. For Reynolds numbers $\lesssim$10, Mei & Auriault (1991) derived a cubic correction to Darcy's law of the form

(1.7)\begin{equation} Hg = \frac{d^2}{K} Re + \hat{b}\,Re^3 \quad\text{or}\quad F_k' = \frac{d^2}{K} Re^{{-}1} + \hat{b}\, Re \end{equation}

for isotropic porous media (solid line). Firdaouss et al. (1997) derived the same law under the condition that ‘if the pressure gradient is reversed, the seepage velocity should also be reversed with no change in modulus’. They supported their derivation with numerical simulations of two-dimensional porous media flow and further demonstrated that (1.7) is satisfied for several classical experimental datasets up to Reynolds numbers of approximately $16$. Hill, Koch & Ladd (2001) confirmed the theory of Mei & Auriault (1991) for numerical simulations of flow through regular and random sphere packs. For higher Reynolds numbers, the drag is commonly described in terms of the Forchheimer equation (Forchheimer 1901) which is composed of a linear and a quadratic term

(1.8)\begin{equation} Hg = a\,Re + b\,Re^2 \quad\text{or}\quad F_k' = a\,Re^{{-}1} + b . \end{equation}

It should be noted that the Forchheimer equation is not consistent with (1.7). Ergun (1952) proposed empirical correlations for the coefficients $a$ and $b$ based on packed bed experiments. The correlations were further refined by, for example, Macdonald et al. (1979) who aggregated multiple experimental datasets. When the flow becomes turbulent, a change of slope of the resistance curve occurs and a different set of coefficients $a'$, $b'$ must be determined (Fand et al. 1987; Burcharth & Andersen 1995).

Figure 2. Drag coefficient in steady flow through a hexagonal sphere pack together with Darcy's law, the correction of Mei & Auriault (1991) and the modified Ergun equation (Macdonald et al. 1979).

A major difficulty in describing unsteady and oscillatory flow is that the drag force does not depend solely on the instantaneous Reynolds number, but is generally a function of the history of the flow. For example, figure 3 shows the instantaneous drag force (i.e. the first two terms on the right-hand side of (1.3)) as a function of the instantaneous Reynolds number for two of our oscillatory flow simulations. It can be seen in figure 3(a) that for the low frequency case LF5 the instantaneous drag is mostly close to the drag observed in a steady flow at the same instantaneous Reynolds number. Conversely, for the medium frequency case MF5 (figure 3b) the instantaneous drag differs significantly from the drag observed in a steady flow at the same instantaneous Reynolds number. Moreover, a hysteresis loop can be observed which suggests a different behaviour of the flow in the acceleration and deceleration phases of the cycle. For linear oscillatory flow the models of Johnson, Koplik & Dashen (1987) and Pride, Morgan & Gangi (1993) provide an accurate description of the history-dependent drag. The models are formulated in the frequency domain as a non-rational transfer function (‘dynamic permeability’). For nonlinear oscillatory flow, the drag has been described by a Forchheimer-type expression (see (1.8)) and variations thereof with frequency- or time-dependent coefficients (van Gent 1993; Hall, Smith & Turcke 1995). However, the inherent assumption of the model is that the nonlinear drag is a function of the instantaneous Reynolds number. As discussed in the above, this assumption is not valid for some flow configurations. Furthermore, we have shown in our previous work (Unglehrt & Manhart 2022a) that at medium and large Womersley numbers the nonlinear parts of the flow can be out of phase with the superficial velocity.

Figure 3. Comparison of the relation between the instantaneous drag force and the superficial velocity for steady and oscillatory flow: (a) LF5 ($Re=158, Wo=10$); (b) MF5 ($Re=157, Wo=31.62$).

Therefore, the objective of the present contribution is to identify and quantify the drag generation processes in steady and oscillatory flow through a sphere pack. In particular, the analysis is guided by the following questions: How large is the contribution of the pressure drag and the friction drag? What effects contribute to the pressure drag? What is the effect of turbulence? How do these contributions scale with the dimensionless numbers governing the flow? What is the phase of these contributions in oscillatory flow?

We adapt the pressure decomposition of Graham (2019) to unsteady incompressible flow through a periodic porous medium. Based on the Poisson equation for the pressure, the pressure is decomposed into three different components: the first component is a reaction force to the imposed macroscopic pressure gradient; the second component represents the displacement of the flow from the wall due to viscosity; and the third component represents the pressure drag induced by vorticity and dissipation. The resulting decomposition of the volume-averaged Navier–Stokes equations is similar to the approach of Aghaei-Jouybari et al. (2022) and is also closely related to various decompositions of the force on a moving body (Quartapelle & Napolitano 1983; Howe 1989; Yu 2014; Li & Wu 2018; Menon & Mittal 2021). We comment on the relationship between this decomposition and the theory of Johnson et al. (1987) for linear oscillatory porous media flow in Appendix B.2. We then investigate DNS datasets for laminar oscillatory flow (Unglehrt & Manhart 2022a) and steady flow (Sakai & Manhart 2020) through a hexagonal sphere pack. In order to establish a baseline, we apply this decomposition to nonlinear steady flow and linear oscillatory flow. We then proceed to analyse nonlinear oscillatory flow. We investigate the evolution of the drag components over the cycle depending on the Reynolds and Womersley number, with particular focus on the Reynolds number scaling of the drag components. Finally, we discuss the implications of our results for the physical understanding and for the modelling of unsteady flow in porous media.

2. Theory

2.1. Mathematical notation

The equations are written in vector notation according to the ISO 80000-2:2019 standard. In particular, $\boldsymbol {\nabla }=\boldsymbol {e}_i \,\partial (.)/\partial x_i$ denotes the Nabla operator, where $\boldsymbol {e}_i$ are the Cartesian unit vectors; $\boldsymbol {a}\,{\cdot }\,\boldsymbol {b}=a_i b_i$ denotes the inner product, $\boldsymbol{\mathsf {A}}:\boldsymbol{\mathsf {B}}=A_{ij} B_{ij}$ denotes the double inner product, $(\boldsymbol {a}\otimes \boldsymbol {b})_{ij}=a_i b_j$ denotes the outer product and $(\boldsymbol {a}\times \boldsymbol {b})_{i}=\epsilon _{ijk} a_j b_k$ denotes the cross product.

2.2. Differential equations

The flow in the pore space is governed by the incompressible Navier–Stokes equations

(2.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

(2.1b)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{u}\otimes\boldsymbol{u}\right) ={-}\frac{1}{\rho}\boldsymbol{\nabla} p +\nu \Delta\boldsymbol{u} + \frac{1}{\rho}\boldsymbol{f} . \end{gather}$$

The flow is driven by a constant force $\boldsymbol {f}=f_x\,\boldsymbol {e}_x$ or a sinusoidal force $\boldsymbol {f}=f_x\sin \varOmega t \,\boldsymbol {e}_x$ which is constant in space and represents to the macroscopic pressure gradient. On the spheres, the velocity satisfies no-slip and impermeability boundary conditions and for both the velocity $\boldsymbol {u}$ and the deviation pressure $p$ triply periodic boundary conditions are imposed.

2.3. Decomposition of the pressure

In this section, we recall the decomposition of the pressure of Graham (2019) that forms the basis of the discussion in the rest of the article. We start from the Poisson equation for the pressure, which can be derived by taking the divergence of the momentum equation (2.1b):

(2.2a)\begin{equation} \Delta p ={-}\rho\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{u}\otimes\boldsymbol{u}\right) = 2 \rho Q, \end{equation}

where $Q=-\frac {1}{2}(\boldsymbol {\nabla }\otimes \boldsymbol {u}):(\boldsymbol {\nabla }\otimes \boldsymbol {u})^{\mathrm {T}}$ is the second invariant of the velocity gradient tensor (Chong, Perry & Cantwell 1990) that is frequently used for vortex identification (Hunt, Wray & Moin 1988; Dubief & Delcayre 2000). The pressure satisfies periodic boundary conditions at the open domain boundaries and the Neumann boundary condition

(2.2b)\begin{equation} \boldsymbol{\nabla} p \boldsymbol{\cdot}\boldsymbol{n} = \mu\Delta\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{n}+ \boldsymbol{f}\boldsymbol{\cdot}\boldsymbol{n} \end{equation}

at solid walls where $\mu =\rho \nu$ is the dynamic viscosity. The boundary condition can be obtained by projecting the Navier–Stokes equations onto the normal $\boldsymbol {n}$ and using the no-slip and no-penetration conditions for $\boldsymbol {u}$. Note that this boundary condition is not required to solve for the pressure, but it is a property of any sufficiently smooth solution (Sani et al. 2006). Thus, the pressure has three sources with a generally different scaling: the macroscopic pressure gradient; the viscous force; and the convective force. The additive decomposition of Graham (2019) separates these different scalings and results in the following three boundary value problems:

(2.3a) \begin{align} &\text{differential equation} \quad\qquad \text{wall boundary condition}\nonumber\end{align}

(2.3a) \begin{align} &\Delta p^{(a)} = 0, \boldsymbol{\nabla} p^{(a)}\boldsymbol{\cdot}\boldsymbol{n} = \boldsymbol{f}\boldsymbol{\cdot}\boldsymbol{n}, \end{align}

(2.3b) \begin{align} &\Delta p^{(v)} = 0, \boldsymbol{\nabla} p^{(v)}\boldsymbol{\cdot}\boldsymbol{n} = \mu\Delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}, \end{align}

(2.3c) \begin{align} &\Delta p^{(c)} = 2\rho Q, \boldsymbol{\nabla} p^{(c)}\boldsymbol{\cdot}\boldsymbol{n} = 0. \end{align}

By summing up the equations, the pressure Poisson equation and the boundary condition are recovered. The accelerative pressure $p^{(a)}$ counterbalances the wall-normal component of the macroscopic pressure gradient $\boldsymbol {f}$ and therefore ensures that the force field acts tangentially to the wall. This effect is illustrated in figure 4 for flow around a cylinder. Note that Graham (2019) defined the accelerative pressure in terms of the acceleration of a moving body in a stationary frame of reference. Upon changing to a comoving frame of reference, the body becomes stationary and a fictitious force appears in the momentum equation. By identifying the acceleration of the body with $-\boldsymbol {f}/\rho$, we have adapted the decomposition to the present setting. The viscous pressure $p^{(v)}$ arises from unbalanced viscous stresses at the wall (Graham 2019). In particular, we show in Appendix A.1 that the boundary condition of the viscous pressure is given by the divergence of the wall shear stress. Finally, as the $Q$-invariant is equal to the difference between the rotation rate magnitude and the strain rate magnitude (Dubief & Delcayre 2000), the convective pressure $p^{(c)}$ is caused by vortical ($Q>0$) and dissipative ($Q<0$) flow features.

Figure 4. Illustration of the effect of the pressure component $p^{(a)}$. (a) External force field $\boldsymbol {f}$. (b) Projected force field $\boldsymbol {f}-\boldsymbol {\nabla } p^{(a)}$ (blue) and $-\boldsymbol {\nabla } p^{(a)}$ (red).

For turbulent flow, we follow Aghaei-Jouybari et al. (2022) and take the Reynolds average of the pressure decomposition. Then, the mean convective pressure $\bar {p}\vphantom {p}^{(c)}$ contains contributions from the mean velocity $\bar {\boldsymbol {u}}$ and the Reynolds stress tensor:

(2.4) \begin{align} \Delta \bar{p}\vphantom{p}^{(c)} = 2 \rho \bar{Q} ={-}\rho\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\cdot}(\overline{\boldsymbol{u}\otimes\boldsymbol{u}}) ={-}\rho\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\cdot}(\overline{\boldsymbol{u}}\otimes \bar{\boldsymbol{u}})-\rho\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\cdot}( \overline{\boldsymbol{u}^\prime\otimes\boldsymbol{u}^\prime}). \end{align}

In analogy to the terminology for the dissipation rate, we refer to the former contribution as ‘direct’ convective pressure $\bar {p}\vphantom {p}^{(d)}$ and to the latter as ‘turbulent’ convective pressure $\bar {p}\vphantom {p}^{(t)}$.

2.4. Decomposition of the pressure drag

In this section, we decompose the pressure drag in the volume-averaged momentum equation (1.3) into the components due to the accelerative pressure $p^{(a)}$, the viscous pressure $p^{(v)}$ and the convective pressure $p^{(c)}$. An auxiliary potential field allows the pressure drag components to be directly expressed in terms of the sources in the boundary value problems (2.3). First, we define an auxiliary potential $\boldsymbol {\varPhi }$ which satisfies the Laplace equation $\Delta \boldsymbol {\varPhi }=0$ with periodic boundary conditions and $(\boldsymbol {\nabla }\otimes \boldsymbol {\varPhi })^{\mathrm {T}}\boldsymbol {\cdot }\boldsymbol {n}=\boldsymbol {n}$ at solid walls (Batchelor 2000, (6.4.11)). This auxiliary potential also forms the basis of other force decompositions (Quartapelle & Napolitano 1983; Howe 1989; Yu 2014; Li & Wu 2018; Menon & Mittal 2021). Then, we apply Green's second identity to the pressure $p$ and the components of the auxiliary potential $\boldsymbol {\varPhi }$:

(2.5)\begin{equation} \int_{V_{f}}\boldsymbol{\varPhi}\,\Delta p \,\mathrm{d}V = \underbrace{\int_{V_{f}}p \,\Delta \boldsymbol{\varPhi} \,\mathrm{d}V}_{{=}0} + \int_{\partial V_{f}} \boldsymbol{\varPhi} (\boldsymbol{\nabla} p \boldsymbol{\cdot}\boldsymbol{n})\, \mathrm{d}A - \int_{\partial V_{f}} p \,(\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi})^{\mathrm{T}} \boldsymbol{\cdot}\boldsymbol{n}\,\mathrm{d}A . \end{equation}

By the definition of the auxiliary potential, its Laplacian is zero and its wall-normal gradient can be replaced by the normal vector. Furthermore, the integrals over the pore areas cancel due to the periodic boundary conditions on $\boldsymbol {\varPhi }$ and $p$. We get

(2.6)\begin{equation} \int_{V_{f}}\boldsymbol{\varPhi}\,\Delta p \,\mathrm{d}V = \int_{A_{fs}} \boldsymbol{\varPhi}(\boldsymbol{\nabla} p \boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d}A - \int_{A_{fs}} p\,\boldsymbol{n}\,\mathrm{d}A . \end{equation}

Finally, we insert the boundary value problems (2.3) and we obtain the components of the pressure drag force per unit volume

(2.7a)$$\begin{gather} -\boldsymbol{f}_{p}^{(a)}:={-}\frac{1}{V}\int_{A_{fs}} p^{(a)}\,\boldsymbol{n}\,\mathrm{d}A ={-}\frac{1}{V}\int_{A_{fs}} \boldsymbol{\varPhi}\left(\boldsymbol{f}\boldsymbol{\cdot}\boldsymbol{n}\right)\mathrm{d}A, \end{gather}$$

(2.7b)$$\begin{gather}-\boldsymbol{f}_{p}^{(v)}:={-}\frac{1}{V}\int_{A_{fs}} p^{(v)}\,\boldsymbol{n}\,\mathrm{d}A ={-} \frac{1}{V} \int_{A_{fs}} \boldsymbol{\varPhi}\,(\mu\Delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d}A, \end{gather}$$

(2.7c)$$\begin{gather}-\boldsymbol{f}_{p}^{(c)}:={-}\frac{1}{V}\int_{A_{fs}} p^{(c)}\,\boldsymbol{n}\,\mathrm{d}A = \frac{1}{V}\int_{V_{f}}\boldsymbol{\varPhi}\,2 \rho Q \,\mathrm{d}V . \end{gather}$$

The auxiliary potential $\boldsymbol {\varPhi }$ can be considered as an analogue of the influence line in structural mechanics and represents the effect of a pressure source on the integral pressure drag. Vice versa, the components $\varPhi _x$, $\varPhi _y$ and $\varPhi _z$ of the auxiliary potential can be seen as the pressure fields in response to a unit source $\boldsymbol {e}_x$, $\boldsymbol {e}_y$ or $\boldsymbol {e}_z$ distributed uniformly over the surface. Note that the auxiliary potential $\boldsymbol {\varPhi }$ is defined up to a constant, but both $Q$ and $\Delta \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n}$ have zero mean for a domain with periodic and no-slip boundary conditions (see Soria, Ooi & Chong (1997) and Appendix A.2, respectively) and the constant does not affect the result. For simplicity, we constrain $\boldsymbol {\varPhi }$ to have zero mean. Figure 5 shows the distribution of the $x$-component of this potential in the hexagonal sphere pack. We can see that $\varPhi _x$ is an antisymmetric function with respect to the planes $x=0, x=d/2, x=d, \ldots$, and has periodicity $d$ in the $x$-direction. The auxiliary potential is largest at the wall and takes its extreme values near the contact points of the spheres.

Figure 5. Auxiliary potential $\varPhi _x$ in the hexagonal sphere pack (a) in a three-dimensional view and (b) in the plane $\sqrt {3}/3\, y-\sqrt {6}/3\,z = 0$ with the mirror planes $x=0$, $x=d/2$ and $x=d$ of the hexagonal sphere pack. The auxiliary potential $\varPhi _x$ is antisymmetric with respect to these mirror planes. The colours range from $-0.15\,d$ (blue) to $0.15\,d$ (red).

In the following, we briefly discuss the pressure drag components in (2.7). With the (dimensionless) tensor of virtual inertia

(2.8)\begin{equation} \boldsymbol{\mathsf{A}} = \frac{1}{V_{s}} \int_{A_{fs}} \boldsymbol{\varPhi}\otimes\boldsymbol{n} \,\mathrm{d}A \end{equation}

defined in Batchelor (2000, (6.4.15)), we can rewrite the accelerative pressure drag as

(2.9)\begin{equation} -\frac{1}{V}\int_{A_{fs}} p^{(a)}\,\boldsymbol{n}\,\mathrm{d}A ={-}\left(\frac{1}{V}\int_{A_{fs}} \boldsymbol{\varPhi}\otimes\boldsymbol{n}\,\mathrm{d}A \right) \boldsymbol{\cdot} \boldsymbol{f} ={-} \left(1-\epsilon\right) \boldsymbol{\mathsf{A}} \boldsymbol{\cdot}\boldsymbol{f}. \end{equation}

Consequently, the accelerative pressure drag directly counteracts the macroscopic pressure gradient. The tensor $-(1-\epsilon )\boldsymbol{\mathsf {A}}$ is equivalent to the hydrodynamic drag tensor ${\lambda }_{\infty }$ introduced by Lafarge (2009, p. 159) based on the work of Johnson & Sen (1981). Moreover, we demonstrate in Appendix B.1 that the tensor of virtual inertia $\boldsymbol{\mathsf {A}}$ can be related to the well-known ‘high-frequency limit of the dynamic tortuosity’ $\alpha _{\infty }$ by Johnson et al. (1987). As the tensor of virtual inertia can be precomputed for a given geometry, no further model is necessary to describe the accelerative pressure drag. The viscous pressure drag is a weighted surface integral of the source term of the viscous pressure $p^{(v)}$. As demonstrated in Appendix A.3, the viscous pressure drag can be reformulated in terms of the wall shear stress as

(2.10)\begin{equation} -\frac{1}{V}\int_{A_{fs}} p^{(v)}\,\boldsymbol{n}\,\mathrm{d}A = \frac{1}{V}\int_{A_{fs}} (\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi})^{\mathrm{T}} \boldsymbol{\cdot} \boldsymbol{\tau}_{w}\,\mathrm{d}A \end{equation}

for the present boundary conditions. As both the friction drag and the viscous pressure drag are integrals of the wall shear stress with only geometry-dependent weights, these terms should have the same scaling. The convective pressure drag term (2.7c) has also been referred to as ‘$Q$-induced force’ (Aghaei-Jouybari et al. 2022). Due to the fore–aft antisymmetry of the auxiliary potential $\varPhi _x$ in the hexagonal sphere pack (figure 5), drag can only be produced from the part of the distribution of the $Q$-invariant that is antisymmetric with respect to the fore–aft symmetry.

Finally, we can insert the decomposition (2.7) into the volume-averaged momentum equation (1.3):

(2.11)\begin{align} \rho\frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t}= \underbrace{-\frac{1}{V}\int_{A_{\mathrm{fs}}} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\right)^{\mathrm{T}} \boldsymbol{\cdot}\boldsymbol{\tau}_{w}\,\mathrm{d}A}_{\textit{friction and viscous pressure drag}} + \underbrace{\frac{1}{V}\int_{V_{f}}\boldsymbol{\varPhi}\,2 \rho Q \,\mathrm{d}V}_{\textit{convective pressure drag}} + \underbrace{\vphantom{\int_{V_{f}}}\left[\epsilon\boldsymbol{\mathsf{I}}- \left(1-\epsilon\right) \boldsymbol{\mathsf{A}}\right] \boldsymbol{\cdot}\boldsymbol{f}}_{\textit{effective forcing}} . \end{align}

In this form, the drag in the porous media flow is separated into a surface contribution due to the viscous term and a volume contribution due to the convective term of the Navier–Stokes equations. It can also be seen that only a fraction of the macroscopic pressure gradient acts onto the flow. In the remainder of the paper, we apply the pressure drag decomposition to a DNS dataset of steady and oscillatory flow in a hexagonal sphere pack. We also show in the Appendices B.2 and B.3 how the drag terms in (2.11) can be used to rederive the results of Johnson et al. (1987) for linear oscillatory flow at high Womersley numbers and to generalise the theory of Mei & Auriault (1991) to oscillatory flow at low Reynolds numbers, respectively.

3. Methodology

3.1. Description of the flow solver

The simulation dataset used in this paper was obtained using our in-house code MGLET (Manhart, Tremblay & Friedrich 2001). It employs a block-structured Cartesian grid with a staggered arrangement of variables (Harlow & Welch 1965) on which the incompressible Navier–Stokes equations are discretised by means of a finite volume method with second-order central approximations. A third-order low-storage Runge–Kutta method is used for the time integration. In every stage a projection step is performed to make the stage velocity divergence-free. This requires the solution of a discrete Poisson problem for a correction pressure. The no-slip boundary conditions on the spheres are enforced by a second-order accurate ghost-cell immersed boundary method (Peller et al. 2006; Peller 2010). In this approach, the velocity field in the interface cells is approximated by a linear least-squares interpolant that satisfies the no-slip boundary condition. From this the specific volume fluxes are computed for the convective velocities, whereas the point values are computed for the convected velocities. The convective velocities are made divergence-free by a cell-by-cell iterative correction that is coupled to the pressure correction in the field. As a result the immersed boundary method is mass conserving.

3.2. Description of the porous medium geometry

A hexagonal close-packed arrangement of spheres (simply referred to as hexagonal sphere pack) is considered as a porous medium geometry. It is triply periodic with the lattice vectors $d\,\boldsymbol {e}_x$, $1/2\,d\,\boldsymbol {e}_x+\sqrt {3}/2\,d\,\boldsymbol {e}_y$ and $2\sqrt {6}/3\,d\,\boldsymbol {e}_z$, and the primitive unit cell (figure 1b) contains two spheres of diameter $d$ that are placed at the locations $(0,0,0)\,d$ and $(1/2,2\sqrt {3}/3,\sqrt {6}/3)\,d$ (Conway & Sloane 1999, p. 114). The sphere pack has a porosity $\epsilon =1-{\rm \pi} /(3\sqrt {2})=0.26$ which is identical to the porosity of the cubic close-packing studied for example in Hill et al. (2001), Hill & Koch (2002) and He et al. (2019). The hexagonal close-packing arrangement possesses a total number of 24 symmetries (Cockroft 1999, space group 194), for example mirror symmetries about the planes $x=1/2\,d$ and $z=\sqrt {6}/3\,d$.

3.3. Description of the simulation database

The drag decomposition will be evaluated for a collection of DNSs of steady and oscillatory flow through a hexagonal close-packed arrangement of spheres (Conway & Sloane 1999, p. 114). The simulation parameters of the steady and the oscillatory cases are summarised in tables 1 and 2, respectively. Figure 6 shows the distribution of the simulated cases in the $Hg$–$Wo$ parameter space. Note that since oscillatory flow in the quasisteady limit ($Wo\to 0$) is in equilibrium at every instant, its behaviour is identical to the corresponding steady flow.

Table 1. Simulation parameters of the steady cases and root mean square of the pressure drag decomposition residual. The simulations marked with ${\dagger}$ were recomputed at a finer grid resolution compared with Sakai & Manhart (2020). The value $N_{samples}$ denotes the number of snapshots that were collected during the averaging time $T_{avg}$. The residual $r=f_{px}-f_{px}^{(a)}-f_{px}^{(v)}-f_{px}^{(c)}$ of the pressure drag decomposition was computed for each snapshot.

Table 2. Simulation parameters of the oscillatory cases and root mean square of the pressure drag decomposition residual. The simulations marked with ${\dagger}$ were taken from Unglehrt & Manhart (2022a). The residual $r=f_{px}-f_{px}^{(a)}-f_{px}^{(v)}-f_{px}^{(c)}$ of the pressure drag decomposition was computed for each snapshot.

Figure 6. Parameter space for the hexagonal sphere pack. The blue crosses represent the steady simulations in Sakai & Manhart (2020) that correspond to the limit $Wo\to 0$. The open circles denote the simulations in Unglehrt & Manhart (2022a) of laminar oscillatory flow and the red filled circles represent simulations of transitional and turbulent oscillatory flow. The dashed line separates the linear regime on the left-hand side from the nonlinear regime on the right-hand side (Unglehrt & Manhart 2022a).

For all cases, we used a triply periodic simulation domain with the lengths $L_x=2\,d$, $L_y = \sqrt {3}\,d$ and $L_z = 2\sqrt {6}/3\,d$ in the $x$-, $y$- and $z$-directions, respectively. The simulation domain thus contains four primitive unit cells. Choosing an appropriate domain size is essential since the flow may otherwise be constrained to a periodic state far from what would be observed in larger domains. Since laminar flow has the same periodicity as the geometry, it would be sufficient to consider one unit cell. Using multiple unit cells allows us to observe the breaking of this periodicity, which is an indicator of a transitional or turbulent flow state. In these regimes, it is plausible that structures spanning multiple pores could form. However, in their study of turbulent flow through a cubic-close packed array of spheres at $Re=222$, $370$ and $740$, He et al. (2019) found that ‘[$\ldots$] the integral scales for all Reynolds numbers studied in this work are much smaller than the particle diameter and thus the unit cell domain showed little variation in statistics compared to a larger domain’. Consequently, we would expect only minor changes if the domain size were increased. The findings of Agnaou et al. (2016) further support this view; they observed that the critical Reynolds number for the onset of unsteady flow in arrays of cylinders is essentially independent of the domain size if the porosity is small ($\epsilon \lesssim 0.45$).

The steady cases are based on the transient flow simulations by Sakai & Manhart (2020). They classified their flow cases into linear (L), steady nonlinear (SNL), unsteady nonlinear (UNL) and turbulent (T) regimes. For large times, the linear and steady nonlinear cases resulted in a constant flow field, whereas the unsteady nonlinear and turbulent cases resulted in a temporally fluctuating velocity field. The low-Reynolds-number cases were recomputed on a finer grid in order to reduce the errors in the evaluation of the pressure decomposition. Thus, all simulations of steady flow used in the present paper were performed using a resolution of 320 cells per sphere diameter (cpd). The high-Reynolds-number simulations UNL2–T4 were continued in order to collect instantaneous flow fields for a statistical evaluation of the mean and turbulent drag components. When the case UNL1 was continued up to a time $t\left \langle u \right \rangle _{i}/d=10.5$, the chaotic oscillations changed into a decaying harmonic oscillation which indicates that the flow converges to a steady state.

The oscillatory cases are based on the simulations in Unglehrt & Manhart (2022a) of linear and nonlinear laminar oscillatory flow. The cases are grouped according to their Womersley number into the low frequency regime (LF) at $Wo=10$, the medium frequency regime (MF) at $Wo=31.62$ and the high frequency regime (HF) at $Wo=100$ and numbered consecutively from $1$ to $4$ with increasing Hagen number. We performed additional simulations of oscillatory flow (cases LF5, LF6, MF5, MF6 and HF5) that were classified as transitional or turbulent based upon their symmetry behaviour (Unglehrt & Manhart 2022b). The cases HF4 and HF5 were computed at a resolution of 768 cpd and the other oscillatory flow cases were computed at a resolution of 384 cpd. We found in Unglehrt & Manhart (2022a) that the cases LF1, LF2, MF1, MF2, HF1, HF2 show effectively linear behaviour and the cases LF3, MF3, HF3 and LF4, MF4, HF4 exhibit nonlinear effects of comparable strength, respectively. For all simulations, the time series of the volume-averaged velocity as well as instantaneous velocity and pressure fields were saved. For the simulations LF5, LF6, MF5, MF6 and HF5, which are transitional or turbulent, only the instantaneous values are discussed as it is computationally expensive to obtain converged statistics for an oscillatory flow.

In Sakai & Manhart (2020) and Unglehrt & Manhart (2022a), the relationship between the imposed pressure gradient and the superficial velocity in the steady and linear oscillatory cases was validated with results from the literature and the grid resolution was determined based on a grid study. The grid convergence of the new cases LF5, LF6, MF5, MF6 and HF5 was assessed based on simulations with coarser grids at resolutions of 48, 96 and 192 cpd (Appendix C). We found that the differences in the cycle-averaged kinetic energy and in the maximum amplitude of the superficial velocity between the finest and the second finest resolution were less than $1.8\,\%$ for all cases.

3.4. Calculation of the terms in the decomposition

In this section, we describe the details of the evaluation of the terms in the drag decomposition from our simulation data. Moreover, we quantify the errors introduced by the decomposition and the statistical errors. First, the pressure drag components were determined from snapshots of the flow fields. The accelerative pressure drag $\boldsymbol {f}_{p}^{(a)}$ could be calculated in closed form using the tensor of virtual inertia (2.8) obtained from the auxiliary potential. The viscous pressure drag $\boldsymbol {f}_{p}^{(v)}$ was calculated in the form of (2.7b). To obtain the second derivative $\Delta \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n}$ at the surface, the wall normal velocity $v$ was interpolated to a point at wall distance $h=1.5\Delta x$. The value of $\Delta \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n}\vert _{w}$ was then calculated using a Taylor expansion of the wall normal velocity profile

(3.1)\begin{equation} v(y) = \underbrace{\vphantom{\left.\frac{\partial v}{\partial y}\right\vert_{w}} \left.v\right\vert_{w}}_{{=}0} + \underbrace{\left.\frac{\partial v}{\partial y}\right\vert_{w}}_{{=}0} y + \left.\frac{\partial^2 v}{\partial y^2}\right\vert_{w} \frac{y^2}{2} = \left.\Delta\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{n}\right\vert_{w} \frac{y^2}{2} \end{equation}

that satisfies the no-slip, impermeability and incompressibility conditions. The convective pressure drag $\boldsymbol {f}_{p}^{(c)}$ was determined by the volume integral (2.7c). As the integrand $\boldsymbol {\varPhi } Q$ can take large positive and negative values, the numerical evaluation of the integral is a delicate task. Due to the symmetry of the hexagonal sphere pack, flow in the positive and negative $x$-direction should behave the same. To enforce this behaviour, we made the values of the auxiliary potential $\varPhi _x$ antisymmetric with respect to the mirror planes $x=0$, $x=d/2$, $x=d$, etc. (figure 5) by setting $\varPhi _x(x):= [\varPhi _x(x)-\varPhi _x(2d-x)]/2$. Note that this is unnecessary in the continuous setting due to the identities $\left \langle Q \right \rangle _{s}=0$ and (A7), which are, however, not perfectly satisfied in the discrete sense. In addition, the $Q$-invariant was formulated as the divergence of the convective term in order to be consistent with the projection method used in our flow solver. The interface cells were not included in the integration as $Q=0$ at no-slip walls.

We determined the residual of the pressure drag decomposition with respect to the pressure drag force that was directly computed from the instantaneous pressure fields. In tables 1 and 2, we report the root mean square residuals over all snapshots; for the steady cases L4, L6, SNL1, SNL2 and SNL4 we report only the residual at the final time. It can be seen that the balance is closed with satisfactory accuracy considering that the total and viscous pressure drag terms have been computed at a ghost-cell immersed boundary. The residual of the decomposition increases with the Womersley number; this can be explained by the formation of boundary layers that increase the error in the evaluation of the source term for the viscous pressure especially near the contact points of the spheres. A higher residual is also observed for the transitional and turbulent cases.

Second, we determined the friction drag using the volume-averaged momentum balance (1.3). The pressure drag term was computed directly from the instantaneous pressure fields and the superficial acceleration was obtained from the derivative of the time series of the superficial velocity $\left \langle u \right \rangle _{s}$.

Third, for the line plots of the oscillatory cases in § 4 the snapshot values in the last period of each simulation were shifted such that the abscissae $\varphi :=\varOmega t$ lie in $[0,2{\rm \pi} ]$. Since the sinusoidal behaviour of the linear cases was misrepresented by a piecewise linear curve due to the relatively low number of samples (table 2), we used a Fourier series interpolation of the snapshot values. For the cases LF5, LF6, MF5 and MF6 a piecewise cubic interpolation (Akima 1974) was used due to high frequency fluctuations during parts of the cycle.

Finally, we averaged the snapshot values for the steady chaotic and turbulent cases UNL2, T1, T2, T3 and T4. For the cases L4–UNL1 we used only the final flow field of the simulation. To decompose the time-averaged convective pressure drag into its direct and turbulent contributions (see § 2.3), we determined the direct convective pressure drag from the time-averaged velocity field and then computed the turbulent convective pressure drag from the difference between the total and the direct contribution. The time-averaged velocity field was estimated from the snapshots. The number of samples is given in table 1. Since our simulation domain contains eight repetitions of the same pore geometry (Unglehrt & Manhart 2022a), we included shifted copies of every instantaneous field into the average. This led to a nominal increase of the sample size by a factor of eight.

We estimated the statistical error for each drag component with the Student's t-distribution. In all cases the $95\,\%$ confidence interval of the sample average had a half-width smaller than $0.75\,\%$ of the average value. While the underlying assumption of a Gaussian distribution of the sample values was not satisfied for some of the cases, we nevertheless expect that the statistical error has in a similar order of magnitude.

3.5. Calculation of the auxiliary potential field

As the ghost cell immersed boundary method in MGLET (see 3.1) is tailored towards flow with no-slip boundary conditions, we computed the auxiliary potential field with the finite element method (FEM) using the FEniCS solver framework (Logg, Mardal & Wells 2012). We employed uniform meshes of linear tetrahedral elements with resolutions up to 384 cpd. From the numerical solution for the auxiliary potential $\boldsymbol {\varPhi }$, we obtained the tensor of virtual inertia

(3.2)\begin{equation} \boldsymbol{\mathsf{A}} = \begin{bmatrix} 0.1345 & 0 & 0 \\ 0 & 0.1345 & 0 \\ 0 & 0 & 0.1329 \end{bmatrix} \end{equation}

where the off-diagonal terms are numerically zero. Furthermore, we computed the length scale tensor $\boldsymbol{\mathsf {L}}$ defined in § B.2,

(3.3)\begin{align} \boldsymbol{\mathsf{L}} &= 2\left[\epsilon\,\boldsymbol{\mathsf{I}}-\left(1-\epsilon\right)\boldsymbol{\mathsf{A}}\right] \boldsymbol{\cdot}\left[\frac{1}{V}\int_{A_{fs}} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes \boldsymbol{\varPhi}\right)^{\mathrm{T}} \boldsymbol{\cdot} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\right) \mathrm{d}A \right]^{{-}1} \nonumber\\ &= \begin{bmatrix} 0.05886 & 0 & 0 \\ 0 & 0.05922 & 0 \\ 0 & 0 & 0.06011 \end{bmatrix}\,d \end{align}

where the off-diagonal elements are numerically zero, too.

The hexagonal sphere pack is isotropic in the $x$–$y$ plane and possesses the same arrangement of spheres as the face-centred cubic sphere pack. Therefore, we can compare our results with the values of Chapman & Higdon (1992) who give a value $1/F=1.612\times 10^{-1}$ for the ‘electrical formation factor’ $F$, corresponding to a value $\alpha _{\infty } = \epsilon /F = 1.61$ for the high-frequency limit of the dynamic tortuosity, and a value $\varLambda = 0.062\,d$ for the length scale defined by Johnson et al. (1987). From (3.2) and (3.3), we obtain the values

(3.4)$$\begin{gather} \alpha_{\infty} = \left(1-\frac{1-\epsilon}{\epsilon}\frac{A_{11}+A_{22}}{2}\right)^{{-}1} = 1.622 , \end{gather}$$

(3.5)$$\begin{gather}\varLambda = \frac{L_{11}+L_{22}}{2} = 0.05904\,d, \end{gather}$$

which show a satisfactory agreement with the results of Chapman & Higdon (1992).

Figure 7 shows the convergence of $\alpha _{\infty }$ and $\varLambda$ over the resolution, which was successively doubled starting from 12 cpd. At intermediate resolutions, we observe a second-order convergence for $\alpha _{\infty }$ and a first-order convergence for $\varLambda$. The value of $\boldsymbol{\mathsf {L}}$ is uncertain as we expect the velocity potential to behave as ${O}(r^{\sqrt {2}-1})$ close to the contact point, leading to a singular velocity (Cox & Cooker 2000). Consequently, we observe a decrease in the rate of convergence. Nevertheless, we consider the numerical solution for the auxiliary potential $\boldsymbol {\varPhi }$ at a resolution of 384 cpd as well converged.

Figure 7. Mesh convergence of the auxiliary potential solution. We give the difference in the high-frequency limit of the dynamic tortuosity $\alpha _{\infty }$ and the length scale $\varLambda$ relative to their values at a resolution of 384 cpd.

4. Results

In this section, we apply the decomposition of the pressure drag (2.7) to our DNS dataset of flow through a hexagonal sphere pack. First, we analyse the steady flow (§ 4.1) and linear oscillatory flow cases (§ 4.2). These represent the quasisteady limit $Wo\to 0$ and the small amplitude limit $Re\to 0$ and serve as a baseline for discussing of the effects of the Reynolds number and the Womersley number in nonlinear oscillatory flow. We then analyse the nonlinear oscillatory flow data (§ 4.3).

4.1. Stationary flow

In this section, we discuss the decomposed drag of our DNS dataset for steady nonlinear flow. In particular, we analyse the dependence of the different drag components on the Reynolds number. Figure 8 shows the contributions of the drag components to the Reynolds-averaged momentum budget in the $x$-direction:

(4.1) \begin{equation} \frac{1}{\epsilon f_x} \begin{bmatrix}\underbrace{\rho \frac{\mathrm{d}\!\left\langle\bar{u}\right\rangle_{s}}{\mathrm{d}t}}_{{=}0} +\ \bar{f}_{px}^{(a)} +\ \bar{f}_{px}^{(v)} \underbrace{+\ \bar{f}_{px}^{(d)}+\bar{f}_{px}^{(t)}}_{=\bar{f}_{px}^{(c)}}+ \ \bar{f}_{\tau_{w}x}\end{bmatrix} = 1. \end{equation}

Since we have divided the momentum equation by the magnitude of the macroscopic pressure gradient $\epsilon f_x$, the terms represent the fraction of the total drag for each drag component. The accelerative pressure drag $\boldsymbol {f}_{p}^{(a)}$ is a pure function of the macroscopic pressure gradient and the geometry due to its definition in (2.3a); its relative contribution to the total stress balance has a value of $38.4\,\%$ independent of the Reynolds number. The viscous pressure drag $\boldsymbol {f}_{p}^{(v)}$ and the friction drag both decrease with the Reynolds number. At low Reynolds numbers the friction drag is approximately twice as large as the viscous pressure drag. For Reynolds numbers above $36$, the ratio between the terms remains almost constant around $1.7$. The direct convective pressure drag $\bar {f}_{px}^{(d)}$ caused by the time-averaged velocity field starts at zero and increases with the Reynolds number. It overtakes the friction and pressure drag at a Reynolds number of approximately $250$. The drag $\bar {f}_{px}^{(t)}$ caused by the Reynolds stresses is non-zero only for the unsteady nonlinear and turbulent cases. Its share increases with the Reynolds number and reaches $6\,\%$ of the total drag at the highest Reynolds number (which is $22\,\%$ of the direct convective pressure drag).

Figure 8. Drag components in steady flow normalised with the imposed macroscopic pressure gradient $\epsilon f_x$.

In order to investigate the scaling of the drag components with $Re$, we form a friction factor-like quantity by normalising the drag with $\rho$, $\left \langle u \right \rangle _{s}$ and $d$. The result is shown in figure 9. For small Reynolds numbers, especially between the cases L4 and L6, the viscous pressure drag coefficient and the friction drag coefficient decrease with $1/Re$, indicating a linear dependence of these drag components on the Reynolds number. The convective pressure drag coefficient increases proportionally to $Re$, corresponding to a cubic dependence of the drag on $Re$. These observations are consistent with the theory of Mei & Auriault (1991). For large Reynolds numbers (T1–T4), the friction drag coefficient and the viscous pressure drag coefficient approach a scaling with exponents $-0.63$ and $-0.61$, respectively. This is very close to the classical laminar boundary layer scaling $1/\sqrt {Re}$ of the friction coefficient (dashed line). The direct convective drag due to the mean velocity field shows a nearly perfect scaling with $Re^2$ for Reynolds numbers between $91$ and $263$, as indicated by a constant drag coefficient. For higher Reynolds numbers, the direct convective pressure drag coefficient shows a slight decrease. There is no clear scaling for the turbulent convective pressure drag. Although we see neither a quadratic scaling of the convective pressure drag nor a linear scaling of the friction and viscous pressure drag in the steady nonlinear regime ($Re=10$–59), the total drag can be described by the Forchheimer equation (1.8), i.e. the sum of a linear and a quadratic term (Sakai & Manhart 2020).

Figure 9. Drag components in steady flow normalised with $\frac {1}{2}\rho \!\left \langle u \right \rangle _{s}^2/d$. The black lines represent different scalings with the Reynolds number: $1/Re$ (dotted), $1/\sqrt {Re}$ (dashed), $1$ (dash-dotted) and $Re$ (solid). The scaling line for $Re$ was anchored at the case L4, the scaling lines for $1/\sqrt {Re}$ were anchored at the case T4, and the scaling line for $1$ was set to the mean value of the cases UNL1–T2.

4.2. Linear oscillatory flow

In this section, we present the results of the drag decomposition for linear oscillatory flow and compare them with theoretical results from the literature. In particular, we discuss the cases LF1 and LF2 at $Wo=10$, MF1 and MF2 at $Wo=31.62$, and HF1 and HF2 at $Wo=100$; all of which have been shown to exhibit linear behaviour in Unglehrt & Manhart (2022a).

The theoretical behaviour of linear oscillatory flow is well understood (Landau & Lifshits (1987, pp. 83f); Batchelor (2000, pp. 353f); Lafarge (2009)) and is summarised below. The velocities and forces are directly proportional to the magnitude of the macroscopic pressure gradient, $\epsilon \,f_x$; the velocities and forces normalised by $\epsilon \,f_x$ depend only on the Womersley number. At low frequencies ($Wo\to 0$), the velocity is in phase with the forcing and is governed by the steady Stokes equations. At high frequencies, the flow has a boundary layer structure: the bulk flow is irrotational and has a phase lag of $90^\circ$ with respect to the forcing, and the amplitude of the bulk flow decreases as $Wo^{-2}$. Near the wall, the flow behaves like the Stokes boundary layer for which the wall shear stress is history dependent and advances the outer flow velocity by $45^\circ$ (Schlichting & Gersten 2017, p. 142). The superficial velocity can be predicted by Darcy's law or the unsteady Darcy equation (Zhu & Manhart 2016) for low frequencies and by the asymptotics of Johnson et al. (1987) for high frequencies. The well-known model of Johnson et al. (1987) blends these asymptotes and predicts the response of the superficial velocity with good accuracy (Chapman & Higdon 1992; Unglehrt & Manhart 2022a). Please note that the asymptotics of Johnson et al. (1987) can be directly obtained from the drag decomposition (2.11) and the Stokes boundary layer solution (see Appendix B.2). This calculation suggests that the viscous pressure drag and the friction drag have the same time dependence as $Wo\to \infty$.

In the following, we address the question of which processes take up the momentum that is supplied to the flow by the macroscopic pressure gradient. To this end, we rearrange the volume-averaged momentum equation like in (4.1):

(4.2)\begin{equation} \frac{1}{\epsilon f_x} \left[\rho \frac{\mathrm{d}\!\left\langle u\right\rangle_{s}}{\mathrm{d}t} +f_{px}^{(a)} +f_{px}^{(v)}+ f_{px}^{(c)} +f_{\tau_{w}x}\right] = \sin(\varOmega t). \end{equation}

Figure 10 displays the terms of this equation over the course of one period of oscillation ($\varphi :=\varOmega t \mod 2{\rm \pi}$) for the simulations LF1, MF1 and HF1. We observe that the acceleration term increases with the Womersley number whereas the viscous pressure and friction drag decrease with the Womersley number. By definition, the accelerative pressure drag remains constant at $38.4\,\%$ of the macroscopic pressure gradient. At $Wo=10$ more than half of the drag is caused by friction and the viscous pressure. On the other hand, at $Wo=100$ most of the pressure drag is caused by the accelerative pressure and the contributions of the friction and viscous pressure drag decrease. Table 3 summarises the relative amplitudes and the phase lag of the different terms with respect to the macroscopic pressure gradient. It can be seen that both quantities are in line with the theoretical expectations and reflect the change of the velocity field from a Stokes flow to a potential flow with thin boundary layers.

Figure 10. Drag components in linear flow normalised with the amplitude of the imposed macroscopic pressure gradient $\epsilon f_x$ for $Wo=10$ (LF1), $Wo=31.62$ (MF1) and $Wo=100$ (HF1).

Table 3. Relative amplitude and phase lag of the acceleration, the accelerative pressure drag, the friction drag and the viscous pressure drag with respect to to the macroscopic pressure gradient $\epsilon \,f_x\sin (\varOmega t)$ in linear flow. The limits $Wo\to 0$ and $Wo\to \infty$ correspond to Stokes flow (case L4) and potential flow, respectively. Note that the accelerative pressure drag $f_p^{(a)}$ always has a relative amplitude of $38.4\,\%$; the convective pressure drag $f_p^{(c)}$ is negligible in linear flow.

The convective pressure has almost no contribution to the force balance. As in the steady state, the convective pressure drag exhibits a cubic scaling with the Reynolds number. This is demonstrated by the collapse of the suitably normalised $f_{px}^{(c)}$ curves for LF1 and LF2, MF1 and MF2, and HF1 and HF2 in figure 11. The relative intensity of the convective pressure drag decreases strongly with the Womersley number. The cubic scaling follows from the drag decomposition when the symmetries of the flow in the hexagonal sphere pack are taken into account (see Appendix B.3). For $Wo=10$ we can observe a saddle point at the zero crossing of the convective pressure drag, which is consistent with a $\left \langle u \right \rangle _{s}^3$ behaviour of the convective pressure drag. For $Wo=31.62$ and $Wo=100$, this saddle point is absent.

Figure 11. Convective pressure drag normalised with $\rho (\max \left \langle u \right \rangle _{s})^{3}/\nu$ corresponding to the scaling of Mei & Auriault (1991).

4.3. Nonlinear oscillatory flow

In this section, we analyse the simulations of nonlinear oscillatory flow. The momentum budgets for the weakly nonlinear cases LF3, MF3 and HF3 are not shown, as they differ only slightly from the linear regime. However, it can be seen in figure 11 that for these cases the convective pressure drag deviates from the cubic Reynolds number scaling.

For the strongly nonlinear cases, figures 12, 13 and 14 show the terms of the momentum equation for $Wo=10$, $Wo=31.62$ and $Wo=100$, respectively. Like in the previous section, the forces are normalised with the amplitude $\epsilon f_x$ of the macroscopic pressure gradient (cf. (4.2)) such that all terms sum up to $\sin (\varOmega t)$ and the accelerative pressure drag appears as $0.384 \sin (\varOmega t)$.

Figure 12. Drag components in nonlinear flow at $Wo=10$ normalised with the amplitude of the imposed macroscopic pressure gradient $\epsilon f_x$ for $Re=77$ (LF4), $Re=158$ (LF5) and $Re=306$ (LF6).

Figure 13. Drag components in nonlinear flow at $Wo=31.62$ normalised with the amplitude of the imposed macroscopic pressure gradient $\epsilon f_x$ for $Re=73$ (MF4), $Re=157$ (MF5) and $Re=297$ (MF6).

Figure 14. Drag components in nonlinear flow at $Wo=100$ normalised with the amplitude of the imposed macroscopic pressure gradient $\epsilon f_x$ for $Re=252$ (HF4) and $Re=468$ (HF5).

At the lowest Womersley number (figure 12), the acceleration is very small compared with the drag forces and the drag components are mostly in phase with the macroscopic pressure gradient. Hence, the flow can be considered quasisteady. The acceleration shows a distinct non-sinusoidal behaviour due to the nonlinear relationship between the macroscopic pressure gradient and the superficial velocity. The convective pressure drag shows a short plateau at the zero crossings; the duration of the plateau decreases with the Reynolds number. As the Reynolds number increases, the friction drag and the viscous pressure drag decrease whereas the convective pressure drag increases. The amplitudes of these components agree well with the results of the steady cases (figure 8). For the cases LF5 and LF6 we can observe fluctuations in the acceleration and in the convective pressure drag, while the friction drag and the viscous pressure drag do not show any fluctuations. These fluctuations could be attributed to vortex shedding and the transition to turbulence.

At the intermediate Womersley number (figure 13), the acceleration is significantly larger than at the lower Womersley number. The amplitudes of the friction drag, viscous pressure drag and convective pressure drag are comparable to $Wo=10$, but the phases lag behind the macroscopic pressure gradient. The convective pressure drag is close to zero during the acceleration phase of each half-cycle, the duration of which decreases with increasing Reynolds number. This behaviour is similar to the plateaus observed at $Wo=10$. When the acceleration reaches its maximum, the convective pressure drag starts to increase; the acceleration goes to zero and changes its sign. Consequently, the maximum convective pressure drag occurs later than the maximum of the superficial velocity. This is consistent with the observations in Unglehrt & Manhart (2022a) that the maximum kinetic energy of the nonlinear part of the velocity field is delayed with respect to the maximum of the superficial velocity.

At the highest Womersley number (figure 14), the acceleration is the dominant term in the momentum balance. The friction drag and the viscous pressure drag are much smaller than for the other Womersley numbers and have approximately the same magnitude as for linear flow at the same Womersley number. Furthermore, they are shifted in phase with respect to the macroscopic pressure gradient. For the case HF4, the convective pressure drag has a relative magnitude of $8\,\%$ and a nearly sinusoidal waveform; for the case HF5, the magnitude increases to $24\,\%$ and the waveform becomes triangular. The phase lag between the convective pressure drag and the macroscopic pressure gradient decreases with increasing Reynolds number. Remarkably, the triangular waveform of the convective pressure drag can also be observed at low Reynolds numbers (figure 11).

In the following, we investigate the high-Reynolds-number scaling of the friction drag and the viscous and convective pressure drag components. In particular, do the scalings observed in steady flow extend to oscillatory flow? For this analysis we construct different normalisations for the drag components based on the sphere diameter $d$, the density $\rho$, the kinematic viscosity $\nu$ and the cycle maximum of the superficial velocity $\max \left \langle u \right \rangle _{s}$. For the inertial scaling, the convective pressure drag $f_{px}^{(c)}$ is normalised with $\rho (\max \left \langle u \right \rangle _{s})^2/d$, and for the steady laminar boundary layer scaling, the friction drag $f_{\tau _{{w}}}$ and the viscous pressure drag $f_{px}^{(v)}$ are normalised with $\rho \sqrt {\nu }\,(\max \left \langle u \right \rangle _{s})^{3/2}/d^{3/2}$.

Figures 15 and 16 show the friction drag and the viscous pressure drag in the steady laminar boundary layer scaling. At $Wo=10$, the curves of the viscous pressure drag collapse for the cases LF5 and LF6. We do not observe a collapse of the friction drag, but the curves are close. At $Wo=31.62$, we find an excellent agreement of the friction drag amplitude with the steady boundary layer scaling for the cases MF5 and MF6. The normalised amplitudes of the viscous pressure drag also agree with the scaling, but the shape of the curves is different between the cases. At $Wo=100$, we do not observe a collapse of the friction drag and the viscous pressure drag in the steady boundary layer scaling.

Figure 15. Friction drag normalised with $\rho \sqrt {\nu } (\max \left \langle u \right \rangle _{s})^{3/2}/d^{3/2}$ corresponding to a steady laminar boundary layer scaling.

Figure 16. Viscous pressure drag normalised with $\rho \sqrt {\nu } (\max \left \langle u \right \rangle _{s})^{3/2}/d^{3/2}$ corresponding to a steady laminar boundary layer scaling.

Figure 17 shows the convective pressure drag in the inertial normalisation. We observe similar amplitudes of the convective pressure drag at $Wo=10$ and $Wo=31.62$. Moreover, the normalised amplitude of the cases LF6 ($Re=307$) and MF6 ($Re=298$) is consistent with the normalised amplitude of the sum of the direct and turbulent convective pressure drag for the cases T2–T4 in the same Reynolds number range ($Re=263$–$354$). However, we do not observe a collapse of the curves at neither Womersley number and thus we cannot confirm the inertial scaling of the convective pressure drag for the oscillatory cases. At $Wo=100$, we do not observe an inertial scaling in the present range of Reynolds numbers ($Re\leqslant 465$). A striking feature in figure 17 is the phase behaviour at $Wo=31.62$. While at low Reynolds numbers the convective pressure drag is approximately $70^\circ$ out of phase with the forcing, the phase shift decreases with increasing Reynolds number. At $Wo=100$, we can also observe a variation of the phase shift, but no clear trend can be identified.

Figure 17. Convective pressure drag normalised with $\rho (\max \left \langle u \right \rangle _{s})^{2}/d$ corresponding to an inertial scaling.

5. Discussion

In this section, we interpret our results with regard to the dynamics of the pore-scale flow. We then discuss the implications of our findings for model descriptions of unsteady porous media flow.

5.1. Steady flow

For steady flow we observed that the direct convective pressure drag due to the time-averaged velocity field scales approximately with $Re^2$ for high Reynolds numbers ($Re=140$–$350$); the friction drag and the viscous pressure drag scale with $Re^{2\unicode{x2013}0.6} = Re^{1.4}$ for $Re=200$–$350$. Dybbs & Edwards (1984) conducted experiments of steady flow through a hexagonal sphere pack. They reported the emergence of boundary layers and an ‘inertial core flow’ between $Re=1$ and $10$. A consistent flow pattern has been observed in the DNSs (Sakai & Manhart 2020). Similarly, for a simple cubic sphere pack Horton & Pokrajac (2009) put forward a conceptual division of the velocity field into a high speed ‘core flow’ and low speed regions near the spheres. Using dye visualisations, Wegner, Karabelas & Hanratty (1971) obtained the skin friction line pattern in a face-centred cubic sphere pack. In a follow-up study, Karabelas, Wegner & Hanratty (1973) hypothesised the presence of boundary layers between the attachment points and the separation lines along the spheres. A simple boundary layer calculation based on a pressure profile resulted in an approximate agreement with the experimental data. Furthermore, Jolls & Hanratty (1969) electrochemically measured the mass transfer rate and the wall shear stress over a sphere inside a packed bed of porosity $\epsilon =0.41$ at Reynolds numbers between $5$ and $1120$. ‘With the exception of the very rearward portion of the spheres the effect of Reynolds number on the local mass transfer rate and on the local shear stress is what is predicted by boundary layer theory for isolated spheres. This would seem to suggest that flow over most of the surface of the sphere could be described by a three-dimensional boundary layer flow.’

Our results seem to support this conceptual picture in that the observed scaling of the friction drag and viscous pressure drag are consistent with the $Re^{3/2}$ scaling predicted by laminar boundary layer theory under the assumption of a Reynolds number independent core flow. The nearly quadratic scaling of the direct convective pressure drag indeed suggests that the time-averaged core flow varies only weakly with the Reynolds number. Furthermore, He et al. (2019, figures 2 and 3) and Sakai & Manhart (2020, figure 15) found that the turbulent kinetic energy is concentrated in the large pores and is low near the walls and where the time-averaged velocity is high. This substantiates the hypothesis of a laminar boundary layer even in the ‘turbulent’ flow regime.

Future research should attempt to confirm the applicability of the boundary layer concept to the present flow configuration based on velocity profiles or the local momentum budget. The presence of laminar boundary layers would allow us to extrapolate the viscous drag to higher Reynolds numbers and would also imply a scaling for the heat and mass transfer in the vicinity of the wall (Karabelas, Wegner & Hanratty 1971; Schlichting & Gersten 2017, ch. 9). This could be important, for example, in the design of chemical reactors. It would also be interesting to extend the present analysis to higher Reynolds numbers to investigate the scaling of the turbulent convective pressure drag.

Given that the observed low-Reynolds-number behaviour agrees with the theory of Mei & Auriault (1991) for isotropic porous media and that the experiments in disordered packed beds point to a quadratic scaling of the drag (Macdonald et al. 1979) and a boundary layer scaling of the friction drag (Jolls & Hanratty 1969) at high Reynolds numbers, we expect the present scalings to carry over qualitatively also to other kinds of sphere packs.

5.2. Oscillatory flow

For oscillatory flow at $Wo=10$ we found that the amplitudes and scalings of the different drag components are very similar to the steady case. In the cases LF5 and LF6 some fluctuations can be observed in the convective pressure drag and in the acceleration (and thus the superficial velocity); the friction drag and the viscous pressure drag show only small traces of these fluctuations (figure 12). This further supports the above hypothesis that the laminar boundary layers are only weakly influenced by inertial and turbulent effects.

At $Wo=31.62$, the amplitudes of the drag components are still close to the steady values, but the phases differ considerably from the lower Womersley number. As the Reynolds number increases, the friction drag and viscous pressure drag become increasingly in phase with the macroscopic pressure gradient (figures 10 and 13). Since the Womersley number is relatively high and since the friction and viscous pressure drag approach a steady laminar boundary layer scaling for higher Reynolds numbers, we explain this behaviour using the boundary layer concept. Generally, the boundary layer thickness can be estimated as $\delta \propto \sqrt {\nu t_{B}}$ where $t_{B}$ is the time that a fluid particle spends inside the boundary layer (Schlichting & Gersten 2017, p. 141). In an accelerating flow, $t_{B}$ is just the elapsed time $t$ since the start of the boundary layer formation. When the time reaches the convection time $d/\left \langle u \right \rangle _{s}$, the boundary layer starts to become steady and its thickness is $\delta \propto \sqrt {\nu d/\left \langle u \right \rangle _{s}}$ or $\delta /d \propto Re^{-1/2}$. In this case, the drag is in phase with the superficial velocity. If the period of oscillation is shorter than the convection time, the flow never becomes steady and the boundary layer thickness is $\delta \propto \sqrt {\nu /\varOmega }$ or $\delta /d\propto Wo^{-1}$. In this case, the boundary layer flow is essentially linear and the drag is out of phase with the superficial velocity (cf. § 4.2). When the Womersley number is fixed, the Reynolds number determines if the boundary layer flow reaches a quasisteady state. The process outlined above can be seen in the case MF5 (figure 18a). In the acceleration phase, the boundary layer is thinner than in the steady case; consequently, the drag is higher than in the steady case. Then, the boundary layer growth reaches the steady state value and during the deceleration, the boundary layer remains quasisteady. Thus, the drag coincides with the steady state curve. For the convective pressure drag (figure 18b) we observe a non-sinusoidal time evolution with a plateau around the zero crossings and a high magnitude in between. The shape and phase of the waveform vary considerably with the Reynolds number (figure 17).

Figure 18. Comparison of the relation between the instantaneous drag components and the superficial velocity for steady and oscillatory flow in the case MF5 ($Re=157$, $Wo=31.62$). (a) Sum of friction and viscous pressure drag; (b) convective pressure drag.

In order to extend our understanding of the convective pressure drag, we look at the instantaneous velocity fields of the case MF5 at the beginning and at the end of the steep increase of the convective pressure drag (the times are highlighted by the markers in figure 13). At the first time ($\varphi =0.28{\rm \pi}$), the flow has an instantaneous Reynolds number of $85$ and the convective pressure drag in the $x$-direction is $-3\,\%$ of the instantaneous macroscopic pressure gradient ( $f_{px}^{(c)}/(\epsilon f_x )=-0.03 \sin (0.28{\rm \pi} )$). At the second time ($\varphi =0.52{\rm \pi}$), the instantaneous Reynolds number is at its peak value $157$ and the convective pressure drag in the $x$-direction is $-27\,\%$ of the instantaneous macroscopic pressure gradient ( $f_{px}^{(c)}/(\epsilon f_x )=-0.26 \sin (0.52{\rm \pi} )$). It can be seen in figure 19 that at the beginning of the increase the distribution of the velocity magnitude is roughly fore–aft symmetric with respect to the planes $x=d/2$ and $x=3d/2$. Since a symmetric velocity field has a symmetric distribution of the $Q$-invariant, which is then multiplied with the antisymmetric auxiliary potential $\varPhi _x$, a relatively low convective pressure drag is produced. On the other hand, a non-symmetric velocity magnitude distribution can be observed at the end of the increase of the convective pressure drag. The zero contour of the streamwise velocity component ($u=0$) indicates that the latter field exhibits a large separation region behind the contact points in the oblique cut plane. The comparison of the two velocity fields shown in figure 19 suggests that the steep increase in convective pressure drag is caused by the emergence of the flow separation regions. The plateaus near the zero crossings of the convective pressure drag could thus be seen as attached flow whereas the parts of the cycle with a large convective pressure drag would correspond to separated flow.

Figure 19. Instantaneous velocity magnitude $\vert \boldsymbol {u}\vert$ and $u=0$ contour of the case MF5 at the times marked in figure 13. The colours are normalised with respect to the instantaneous superficial velocity. (a) Beginning of the steep increase of the convective pressure drag ($\varphi =0.28{\rm \pi}, Re(t)=85$). (b) End of the steep increase of the convective pressure drag ($\varphi =0.52{\rm \pi}, Re(t)=157$).

At $Wo=100$, the drag components do not follow the same scalings as at the lower Womersley numbers and clear phase differences between the drag components can be observed. A possible explanation for these discrepancies is that at low Womersley numbers the boundary layers are quasisteady if the Reynolds number is high enough, whereas the boundary layers do not become steady at the highest Womersley number in the considered Reynolds number range. The convective pressure drag has an almost triangular waveform at low Reynolds numbers (figure 11) and at high Reynolds numbers (figure 17). This qualitatively different behaviour of the convective pressure drag in comparison with the lower Womersley numbers could be understood if one assumes a finite formation time for the drag producing structures. Then, at $Wo=10$ the formation time would be small compared with the period of oscillation, resulting in a small phase lag of the convective pressure drag. At $Wo=31.62$, the formation time would be relatively large compared with the period of oscillation (similar to the duration of the plateaus at the zero crossings), resulting in a larger phase lag of the convective pressure drag. Figure 17(b) suggests that the formation time would decrease with increasing Reynolds numbers. Finally, at $Wo=100$ the frequency of oscillation is so high that the formation and destruction in subsequent half-cycles overlaps in time. Thus, the plateau would disappear.

5.3. Implications for modelling

We have presented a new form (2.11) of the volume-averaged momentum equation for a spatially constant macroscopic pressure gradient $\boldsymbol {f}$ where we can express the drag in terms of the wall shear stress and the second invariant of the velocity gradient tensor. The auxiliary potential $\boldsymbol {\varPhi }$ (and derived from it the tensor of virtual inertia $\boldsymbol{\mathsf {A}}$) only depends on the geometry of the porous medium. In this formulation, the components of the pressure drag with a viscous scaling, an inertial scaling and a direct proportionality to the macroscopic pressure gradient are separated. We have shown in Appendix B.2 how this form of the volume-averaged momentum equation can be used to directly derive the asymptotic drag behaviour at high Womersley numbers of Johnson et al. (1987).

For steady flow, we found that the friction drag and the viscous pressure drag depend linearly on $Re$ at low Reynolds numbers and scale with $Re^{1.4}$ at high Reynolds numbers. The convective pressure drag scales with $Re^3$ at low Reynolds numbers and with $Re^2$ at high Reynolds numbers. At low Reynolds numbers, these results are in line with Darcy's law (1.6) and its correction (1.7) by Mei & Auriault (1991). However, the Forchheimer equation (1.8) is incompatible with the low-Reynolds-number behaviour of the convective pressure drag and with the high-Reynolds-number behaviour of the friction drag and of the viscous pressure drag.

In nonlinear oscillatory flow at $Wo=10$ the drag components show the same scaling as in steady flow. Moreover, the momentum balance indicates that the flow is quasisteady. This flow can thus be modelled by extending the steady state drag law with an acceleration term (Zhu et al. 2014; Zhu & Manhart 2016). At $Wo=31.62$ the Reynolds number scalings of the drag components are similar to the steady case, but the drag components are out of phase with the superficial velocity (figure 18). To model the friction and viscous pressure drag, a promising approach could be to blend the parametrisation of Johnson et al. (1987) with the $Re^{3/2}$ behaviour of the laminar boundary layer. As the convective pressure drag cannot be expressed as a function of the instantaneous superficial velocity alone and, furthermore, scales with $Re^3$ at low Reynolds numbers, it seems necessary to think beyond the traditional parametrisation in terms of $\left \langle u \right \rangle _{s}^2$. In particular, we could observe a smaller hysteresis between the convective pressure drag and a time-lagged superficial velocity or the instantaneous kinetic energy. For the simulation cases at $Wo=100$ no clear high-Reynolds-number scalings could be identified; thus, further research is required in this direction. As a starting point for the development of improved models, we provide the time series of the superficial velocity and the drag components from our simulations as supplementary material available at https://doi.org/10.1017/jfm.2023.798.

Finally, our decomposition provides a new point of view on the time constant in the volume-averaged momentum equation. In most model equations for unsteady porous media flow, the resistance of the bulk flow to acceleration has been incorporated with the ad hoc addition of a ‘virtual mass coefficient’ (Sollitt & Cross 1972; Burcharth & Andersen 1995) or ‘acceleration coefficient tensor’ (Nield 1991) to the volume-averaged momentum equation. This was done by analogy to the added-mass effect in inviscid flow. For example, Nield (1991) suggested an unsteady extension to Darcy's law (1.6),

(5.1)\begin{equation} \rho \boldsymbol{\mathsf{C}}_{a} \boldsymbol{\cdot}\frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t}={-} \frac{\mu}{K}\left\langle\boldsymbol{u}\right\rangle_{s}+\boldsymbol{f} \end{equation}

where the acceleration coefficient tensor is assumed to be of the form $\boldsymbol{\mathsf {C}}_{a}=\epsilon ^{-1}\boldsymbol{\mathsf {I}}+\boldsymbol{\mathsf {N}}$; the tensor $\boldsymbol{\mathsf {N}}$ representing ‘the contribution from “fractures” ’. The volume-averaged momentum equation (2.11) can also be brought to such a form by multiplying the equation with $\boldsymbol{\mathsf {C}}_{a}:=[\epsilon \boldsymbol{\mathsf {I}}- (1-\epsilon ) \boldsymbol{\mathsf {A}}]^{-1}$. Then, the accelerative pressure drag is absorbed into the prefactor of the acceleration and all other drag terms are rescaled:

(5.2)\begin{equation} \rho \boldsymbol{\mathsf{C}}_{a} \boldsymbol{\cdot}\frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t}= \boldsymbol{\mathsf{C}}_{a} \boldsymbol{\cdot}\left[-\frac{1}{V}\int_{A_{fs}} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\right)^{\mathrm{T}} \boldsymbol{\cdot}\boldsymbol{\tau}_{w}\,\mathrm{d}A + \frac{1}{V}\int_{V_{f}}\boldsymbol{\varPhi}\,2 \rho Q \,\mathrm{d}V\right] +\boldsymbol{f} . \end{equation}

The term $-\mu /K\,\left \langle \boldsymbol {u}\right \rangle _{s}$ in (5.1) can be identified as a parametrisation of the first term on the right-hand side of (5.2) with the Darcy expression for the drag. Our decomposition thus gives a new interpretation to the ‘virtual mass’ in a porous medium in terms of the accelerative pressure drag, which possesses a clear physical meaning also for viscous flow. As discussed in Appendix B.1, this definition of the acceleration coefficient reduces to the ‘high-frequency limit of the dynamic tortuosity’ by Johnson et al. (1987) in the isotropic case, i.e. $\boldsymbol{\mathsf {C}}_{\mathrm {a}}=\alpha _{\infty }/\epsilon \,\boldsymbol{\mathsf {I}}$.

6. Conclusion

In this paper, we studied the behaviour of the drag force in steady and oscillatory flow through a hexagonal sphere pack. Based on the pressure decomposition of Graham (2019) we derived a new form of the volume-averaged momentum equation in which the pressure drag force is split into three contributions. The accelerative pressure drag is a reaction force directly proportional the macroscopic pressure gradient. It prevents the macroscopic pressure gradient from accelerating the fluid normal to the wall. The viscous pressure drag results from unbalanced viscous stresses and can be expressed to a weighted integral of the wall shear stress. The convective pressure drag can be expressed as a weighted volume integral of the $Q$-invariant of the velocity gradient tensor representing effects like vortices, shear layers and flow separation.

Using this decomposition, the drag law for high Womersley numbers (Johnson et al. 1987) and the $Re$ dependence of the drag for low Reynolds numbers could be derived using relatively simple arguments (see §§ B.2 and B.3). Moreover, we could provide a new theoretical basis for the virtual mass coefficient commonly employed in models for unsteady porous media flow (see § 5.3).

We then applied the drag decomposition to a DNS dataset of steady and oscillatory flow through a hexagonal sphere pack. We investigated the contributions of the different drag terms to the volume-averaged momentum budget. The accelerative pressure drag is proportional to the macroscopic pressure gradient and thus has a fixed contribution of $38.4\,\%$ to the momentum budget. For steady flow, the remaining drag is dominated by the friction and viscous pressure drag at low Reynolds numbers and by the convective pressure drag at high Reynolds numbers. For the considered Reynolds numbers, the Reynolds stresses only have a minor effect on the drag. For oscillatory flow at low and medium Womersley numbers, the friction drag, viscous pressure drag and convective pressure drag have a similar magnitude as in the steady case. At high Womersley numbers, the friction and viscous pressure drag are significantly smaller than in the steady case. Thus, the drag at high Womersley numbers is made up mostly by the accelerative and the convective pressure drag. An important feature of the drag in oscillatory flow is that the drag components are not in phase with the body force and the superficial velocity. The phase differences increase with the Womersley number.

We investigated the Reynolds number scalings of the friction drag, the viscous pressure drag and the convective pressure drag. In the steady case, the friction and viscous pressure drag are proportional to $Re$ at small Reynolds numbers and scale with $Re^{1.4}$ for Reynolds numbers between $200$ and $350$. The convective pressure drag of the time-averaged velocity field scales with $Re^3$ up to a Reynolds number of $10$ and with $Re^2$ for $Re=140\unicode{x2013}350$. For oscillatory flow, the same amplitude scalings can be observed at $Wo=10$ and $Wo=31.62$, whereas no clear high $Re$ scaling could be found for the cases at $Wo=100$.

These scalings support the picture of Dybbs & Edwards (1984) who divided the flow at higher Reynolds numbers into an inertial core flow and viscous boundary layers, where we linked the former with the convective pressure drag and the latter with the friction and viscous pressure drag. The visualisation of instantaneous velocity fields suggests that the convective pressure drag in the hexagonal sphere pack is caused by large flow separations. Moreover, the clear scalings of the friction and viscous pressure drag and of the convective pressure drag indicate that the inertial core and the boundary layers are only weakly affected by the turbulence for $Re=200\unicode{x2013}350$.

In future work, the present theory for periodic porous media could be extended to non-periodic porous media. This might be realised by rewriting the identity

(6.1)\begin{equation} \left\langle\boldsymbol{\varPhi}\,\Delta P\right\rangle_{s} = \left\langle\boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{\nabla}P\otimes\boldsymbol{\varPhi})\right\rangle_{s}-\left\langle\boldsymbol{\nabla}\boldsymbol{\cdot}[(\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi})P]\right\rangle_{s} \end{equation}

with the spatial averaging theorem (Whitaker 1985); together with the volume-averaged Navier–Stokes equations (Whitaker 1986, 1996) a generalisation of (2.11) would be obtained.