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 Reference Ene and Sanchez-Palencia1975). 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
\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 Reference Agnaou, Lasseux and Ahmadi2016). 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
\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 Reference Ene and Sanchez-Palencia1975). 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
\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 (Reference Whitaker1986, Reference Whitaker1996) 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. (Reference Davit2013). 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}$.
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 Reference Ene and Sanchez-Palencia1975; Firdaouss, Guermond & Le Quéré Reference Firdaouss, Guermond and Le Quéré1997; Iervolino, Manna & Vacca Reference Iervolino, Manna and Vacca2010; Lasseux, Valdés-Parada & Bellet Reference Lasseux, Valdés-Parada and Bellet2019). 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
\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. Reference Macdonald, El-Sayed, Mow and Dullien1979) as
\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 (Reference Sakai and Manhart2020) 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
\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 (Reference Mei and Auriault1991) derived a cubic correction to Darcy's law of the form
\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. (Reference Firdaouss, Guermond and Le Quéré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 (Reference Hill, Koch and Ladd2001) confirmed the theory of Mei & Auriault (Reference Mei and Auriault1991) 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 Reference Forchheimer1901) which is composed of a linear and a quadratic term
\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 (Reference Ergun1952) 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. (Reference Macdonald, El-Sayed, Mow and Dullien1979) 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. Reference Fand, Kim, Lam and Phan1987; Burcharth & Andersen Reference Burcharth and Andersen1995).
Drag coefficient in steady flow through a hexagonal sphere pack together with Darcy's law, the correction of Mei & Auriault (Reference Mei and Auriault1991) and the modified Ergun equation (Macdonald et al. Reference Macdonald, El-Sayed, Mow and Dullien1979).
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 (Reference Johnson, Koplik and Dashen1987) and Pride, Morgan & Gangi (Reference Pride, Morgan and Gangi1993) 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 Reference van Gent1993; Hall, Smith & Turcke Reference Hall, Smith and Turcke1995). 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 Reference Unglehrt and Manhart2022a) that at medium and large Womersley numbers the nonlinear parts of the flow can be out of phase with the superficial velocity.
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 (Reference Graham2019) 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. (Reference Aghaei-Jouybari, Seo, Yuan, Mittal and Meneveau2022) and is also closely related to various decompositions of the force on a moving body (Quartapelle & Napolitano Reference Quartapelle and Napolitano1983; Howe Reference Howe1989; Yu Reference Yu2014; Li & Wu Reference Li and Wu2018; Menon & Mittal Reference Menon and Mittal2021). We comment on the relationship between this decomposition and the theory of Johnson et al. (Reference Johnson, Koplik and Dashen1987) for linear oscillatory porous media flow in Appendix B.2. We then investigate DNS datasets for laminar oscillatory flow (Unglehrt & Manhart Reference Unglehrt and Manhart2022a) and steady flow (Sakai & Manhart Reference Sakai and Manhart2020) 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
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$
$$\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 (Reference Graham2019) 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):
\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 Reference Chong, Perry and Cantwell1990) that is frequently used for vortex identification (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988; Dubief & Delcayre Reference Dubief and Delcayre2000). The pressure satisfies periodic boundary conditions at the open domain boundaries and the Neumann boundary condition
\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. Reference Sani, Shen, Pironneau and Gresho2006). 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 (Reference Graham2019) separates these different scalings and results in the following three boundary value problems:
\begin{align}
&\text{differential equation} \quad\qquad \text{wall boundary condition}\nonumber\end{align}
\begin{align}
&\Delta p^{(a)} = 0, \boldsymbol{\nabla} p^{(a)}\boldsymbol{\cdot}\boldsymbol{n} = \boldsymbol{f}\boldsymbol{\cdot}\boldsymbol{n},
\end{align}
\begin{align}
&\Delta p^{(v)} = 0, \boldsymbol{\nabla} p^{(v)}\boldsymbol{\cdot}\boldsymbol{n} = \mu\Delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n},
\end{align}
\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 (Reference Graham2019) 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 Reference Graham2019). 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 Reference Dubief and Delcayre2000), the convective pressure 
$p^{(c)}$ is caused by vortical (
$Q>0$) and dissipative (
$Q<0$) flow features.
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. (Reference Aghaei-Jouybari, Seo, Yuan, Mittal and Meneveau2022) 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:
\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 Reference Batchelor2000, (6.4.11)). This auxiliary potential also forms the basis of other force decompositions (Quartapelle & Napolitano Reference Quartapelle and Napolitano1983; Howe Reference Howe1989; Yu Reference Yu2014; Li & Wu Reference Li and Wu2018; Menon & Mittal Reference Menon and Mittal2021). Then, we apply Green's second identity to the pressure 
$p$ and the components of the auxiliary potential 
$\boldsymbol {\varPhi }$:
\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
\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
$$\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}$$
$$\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}$$
$$\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 (Reference Soria, Ooi and Chong1997) 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.
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
\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 (Reference Batchelor2000, (6.4.15)), we can rewrite the accelerative pressure drag as
\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 (Reference Lafarge2009, p. 159) based on the work of Johnson & Sen (Reference Johnson and Sen1981). 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. (Reference Johnson, Koplik and Dashen1987). 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
\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. Reference Aghaei-Jouybari, Seo, Yuan, Mittal and Meneveau2022). 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):
\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. (Reference Johnson, Koplik and Dashen1987) for linear oscillatory flow at high Womersley numbers and to generalise the theory of Mei & Auriault (Reference Mei and Auriault1991) 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 Reference Manhart, Tremblay and Friedrich2001). It employs a block-structured Cartesian grid with a staggered arrangement of variables (Harlow & Welch Reference Harlow and Welch1965) 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. Reference Peller, Le Duc, Tremblay and Manhart2006; Peller Reference Peller2010). 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 Reference Conway and Sloane1999, 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. (Reference Hill, Koch and Ladd2001), Hill & Koch (Reference Hill and Koch2002) and He et al. (Reference He, Apte, Finn and Wood2019). The hexagonal close-packing arrangement possesses a total number of 24 symmetries (Cockroft Reference Cockroft1999, 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 Reference Conway and Sloane1999, 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.
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 (Reference Sakai and Manhart2020). 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.
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 (Reference Unglehrt and Manhart2022a). The residual 
$r=f_{px}-f_{px}^{(a)}-f_{px}^{(v)}-f_{px}^{(c)}$ of the pressure drag decomposition was computed for each snapshot.
Parameter space for the hexagonal sphere pack. The blue crosses represent the steady simulations in Sakai & Manhart (Reference Sakai and Manhart2020) that correspond to the limit 
$Wo\to 0$. The open circles denote the simulations in Unglehrt & Manhart (Reference Unglehrt and Manhart2022a) 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 Reference Unglehrt and Manhart2022a).
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. (Reference He, Apte, Finn and Wood2019) 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. (Reference Agnaou, Lasseux and Ahmadi2016) 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 (Reference Sakai and Manhart2020). 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 (Reference Unglehrt and Manhart2022a) 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 Reference Unglehrt and Manhart2022b). 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 (Reference Unglehrt and Manhart2022a) 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 (Reference Sakai and Manhart2020) and Unglehrt & Manhart (Reference Unglehrt and Manhart2022a), 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
\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 Reference Akima1974) 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 Reference Unglehrt and Manhart2022a), 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 Reference Logg, Mardal and Wells2012). 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
\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,
\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 (Reference Chapman and Higdon1992) 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. (Reference Johnson, Koplik and Dashen1987). From (3.2) and (3.3), we obtain the values
$$\begin{gather} \alpha_{\infty} = \left(1-\frac{1-\epsilon}{\epsilon}\frac{A_{11}+A_{22}}{2}\right)^{{-}1} = 1.622 , \end{gather}$$
$$\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 (Reference Chapman and Higdon1992).
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 Reference Cox and Cooker2000). 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.
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:
\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).
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 (Reference Mei and Auriault1991). 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 Reference Sakai and Manhart2020).
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 (Reference Unglehrt and Manhart2022a).
The theoretical behaviour of linear oscillatory flow is well understood (Landau & Lifshits (Reference Landau and Lifshits1987, pp. 83f); Batchelor (Reference Batchelor2000, pp. 353f); Lafarge (Reference Lafarge2009)) 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 Reference Schlichting and Gersten2017, p. 142). The superficial velocity can be predicted by Darcy's law or the unsteady Darcy equation (Zhu & Manhart Reference Zhu and Manhart2016) for low frequencies and by the asymptotics of Johnson et al. (Reference Johnson, Koplik and Dashen1987) for high frequencies. The well-known model of Johnson et al. (Reference Johnson, Koplik and Dashen1987) blends these asymptotes and predicts the response of the superficial velocity with good accuracy (Chapman & Higdon Reference Chapman and Higdon1992; Unglehrt & Manhart Reference Unglehrt and Manhart2022a). Please note that the asymptotics of Johnson et al. (Reference Johnson, Koplik and Dashen1987) 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):
\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.
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).
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.
Convective pressure drag normalised with 
$\rho (\max \left \langle u \right \rangle _{s})^{3}/\nu$ corresponding to the scaling of Mei & Auriault (Reference Mei and Auriault1991).
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)$.
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).
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).
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 (Reference Unglehrt and Manhart2022a) 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.
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.
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.
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 (Reference Dybbs and Edwards1984) 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 Reference Sakai and Manhart2020). Similarly, for a simple cubic sphere pack Horton & Pokrajac (Reference Horton and Pokrajac2009) 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 (Reference Wegner, Karabelas and Hanratty1971) obtained the skin friction line pattern in a face-centred cubic sphere pack. In a follow-up study, Karabelas, Wegner & Hanratty (Reference Karabelas, Wegner and Hanratty1973) 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 (Reference Jolls and Hanratty1969) 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. (Reference He, Apte, Finn and Wood2019, figures 2 and 3) and Sakai & Manhart (Reference Sakai and Manhart2020, 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 Reference Karabelas, Wegner and Hanratty1971; Schlichting & Gersten Reference Schlichting and Gersten2017, 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 (Reference Mei and Auriault1991) for isotropic porous media and that the experiments in disordered packed beds point to a quadratic scaling of the drag (Macdonald et al. Reference Macdonald, El-Sayed, Mow and Dullien1979) and a boundary layer scaling of the friction drag (Jolls & Hanratty Reference Jolls and Hanratty1969) 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 Reference Schlichting and Gersten2017, 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).
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.
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. (Reference Johnson, Koplik and Dashen1987).
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 (Reference Mei and Auriault1991). 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. Reference Zhu, Waluga, Wohlmuth and Manhart2014; Zhu & Manhart Reference Zhu and Manhart2016). 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. (Reference Johnson, Koplik and Dashen1987) 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 Reference Sollitt and Cross1972; Burcharth & Andersen Reference Burcharth and Andersen1995) or ‘acceleration coefficient tensor’ (Nield Reference Nield1991) to the volume-averaged momentum equation. This was done by analogy to the added-mass effect in inviscid flow. For example, Nield (Reference Nield1991) suggested an unsteady extension to Darcy's law (1.6),
\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:
\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. (Reference Johnson, Koplik and Dashen1987) 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 (Reference Graham2019) 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. Reference Johnson, Koplik and Dashen1987) 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 (Reference Dybbs and Edwards1984) 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
\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 Reference Whitaker1985); together with the volume-averaged Navier–Stokes equations (Whitaker Reference Whitaker1986, Reference Whitaker1996) a generalisation of (2.11) would be obtained.
Supplementary material
The time series of the volume-averaged momentum budget terms are provided for all simulation cases. Moreover, the time series of the superficial velocity and kinetic energy components are provided for the cases LF5, LF6, MF5, MF6 and HF5. Supplementary material is available at https://doi.org/10.1017/jfm.2023.798.
Acknowledgements
We would like to express our thanks to Y. Sakai for providing the complete simulation dataset for the steady flow cases. We further thank Y. Sakai, U. Jenssen, J. Brosda and S. Wenczowski for helpful and interesting discussions and we thank H.I. Andersson for reading the manuscript and providing valuable feedback. Finally, we thank the anonymous reviewers for their thorough and constructive criticism.
Funding
The authors gratefully acknowledge the financial support of the DFG under grant no. MA2062/13-1. Computing time was granted by the Leibniz Supercomputing Centre on its Linux-Cluster.
Declaration of interests
The authors report no conflict of interest.
Author contributions
L.U. performed the simulations advised by M.M. L.U. derived the theory of the drag decomposition and performed the evaluation of the simulations. Both authors contributed to reaching conclusions and in writing the paper.
Appendix A. Notes on the viscous pressure drag
This appendix discusses some aspects of the relationship between the viscous pressure and the wall shear stress. In § A.1, we show that the boundary condition of the viscous pressure 
$p^{(v)}$ can be expressed in terms of the wall shear stress divergence. In § A.2, we show that the mean value of the viscous pressure source term is zero for a periodic domain. Finally, in § A.3 we derive an alternative expression for the viscous pressure drag as a weighted integral of the wall shear stress.
A.1. Relationship between the viscous pressure and the wall shear stress divergence
The boundary condition of the viscous pressure (2.3b) can be rewritten using the identity 
$\Delta \boldsymbol {u} = -\boldsymbol {\nabla }\times (\boldsymbol {\nabla }\times \boldsymbol {u})$ and Stokes’ theorem as
\begin{equation} \boldsymbol{\nabla} p^{(v)}\boldsymbol{\cdot}\boldsymbol{n} = \mu\Delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n} ={-}\mu\left[\boldsymbol{\nabla}\times(\boldsymbol{\nabla}\times\boldsymbol{u})\right]\boldsymbol{\cdot}\boldsymbol{n} ={-}\mu\lim_{A\to 0} \frac{1}{A}\int_{\partial A} (\boldsymbol{\nabla}\times\boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{t}\,\mathrm{d}s, \end{equation}
where 
$\boldsymbol {n}$ is the normal vector pointing from the fluid towards the wall, 
$A$ is an small surface patch on the wall and 
$\boldsymbol {t}$ represents the tangent vector on its boundary 
$\partial A$. The vorticity on the wall is related to the wall shear stress by the equation
\begin{equation} \boldsymbol{\omega}_{w} = \left.\boldsymbol{\nabla}\times\boldsymbol{u}\right\vert_{w}= \frac{\boldsymbol{\tau}_{w}}{\mu}\times\boldsymbol{n} \end{equation}
where the cross product expresses a clockwise rotation of the wall shear stress by 
$90^\circ$ around the normal. Figure 20 shows the orientation of the vectors with respect to a single sphere.
(a) Orientation of the normal vector 
$\boldsymbol {n}$, the tangent vector 
$\boldsymbol {t}$ and their cross product 
$\boldsymbol {t}\times \boldsymbol {n}$ with respect to the surface patch 
$A$. The normal vector points from the fluid outside into the sphere. (b) Orientation of the normal vector 
$\boldsymbol {n}$, the wall shear stress 
$\boldsymbol {\tau }_{w}$ and the wall vorticity 
$\boldsymbol {\omega }_{w}$ with respect to the surface.
Using Lagrange's identity for the cross product (Bronstein et al. Reference Bronstein, Semendjajew, Grosche, Ziegler and Ziegler1991, p. 556), we can establish the relation
\begin{equation} \left[(\boldsymbol{\nabla}\times\boldsymbol{u})\times\boldsymbol{n}\right]\boldsymbol{\cdot} \left[\boldsymbol{t}\times\boldsymbol{n}\right]=\left[(\boldsymbol{\nabla}\times\boldsymbol{u})\boldsymbol{\cdot} \boldsymbol{t}\right]\underbrace{\left[\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{n}\right]}_{{=}1}-\underbrace{\left[(\boldsymbol{\nabla} \times\boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{n}\right]}_{{=}0}\underbrace{\left[\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{t}\right]}_{{=}0} = (\boldsymbol{\nabla}\times\boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{t} \end{equation}and by combining the two expressions, we obtain
\begin{equation} (\boldsymbol{\nabla}\times\boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{t} = \left[\left(\frac{\boldsymbol{\tau}_{w}}{\mu}\times\boldsymbol{n}\right)\times\boldsymbol{n}\right]\boldsymbol{\cdot} \left[\boldsymbol{t}\times\boldsymbol{n}\right] ={-}\frac{\boldsymbol{\tau}_{w}}{\mu}\boldsymbol{\cdot}(\boldsymbol{t}\times\boldsymbol{n}) . \end{equation}Finally, we arrive at the boundary condition
\begin{equation} \boldsymbol{\nabla} p^{(v)}\boldsymbol{\cdot}\boldsymbol{n} = \lim_{A\to 0} \frac{1}{A}\int_{\partial A} \boldsymbol{\tau}_{w}\boldsymbol{\cdot}(\boldsymbol{t}\times\boldsymbol{n})\,\mathrm{d}s \end{equation}
where the right-hand side can be understood as the divergence of the wall shear stress, since the vector 
$\boldsymbol {t}\times \boldsymbol {n}$ represents the outward normal vector on the boundary 
$\partial A$ of the surface patch along the wall. Consequently, the viscous pressure 
$p^{(v)}$ is caused by imbalances in the wall shear stress. For example, a stagnation point represents a source of the wall shear stress, hence 
$\boldsymbol {\nabla } p^{(v)}\boldsymbol {\cdot }\boldsymbol {n}>0$ and the viscous pressure increases towards the wall. This is indeed observed in the analytical solution (Graham Reference Graham2019).
A.2. Zero-mean property of the wall-normal friction force for a periodic domain
We apply Gauss's integral theorem to the vector field 
$\Delta \boldsymbol {u}$:
\begin{equation} \int_{\partial V_{f}} \Delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n} \,\mathrm{d}A = \int_{V_{f}} \boldsymbol{\nabla}\boldsymbol{\cdot}(\Delta\boldsymbol{u}) \,\mathrm{d}V = \int_{V_{f}} \varDelta(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})\,\mathrm{d}V = 0 . \end{equation}As the velocity field is periodic, we can further simplify this to
\begin{equation} \int_{A_{fs}} \Delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n} \,\mathrm{d}A = 0, \end{equation}
where 
$A_{fs}$ represents the fluid–solid interface.
A.3. Alternative expression for the viscous pressure drag
We can rewrite viscous pressure drag (2.7b) using the periodic boundary conditions and Gauss's theorem as
\begin{align} -\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 ={-} \frac{1}{V} \int_{\partial V_{f}} \boldsymbol{\varPhi}\,(\mu\Delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d}A \nonumber\\ &={-} \frac{1}{V} \int_{V_{f}} \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \mu\Delta\boldsymbol{u} \otimes\boldsymbol{\varPhi}\right)\mathrm{d}V \nonumber\\ &= \underbrace{\frac{1}{V} \int_{V_{f}} \boldsymbol{\varPhi}\, \mu\varDelta (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}) \,\mathrm{d}V}_{{=}0} - \frac{1}{V} \int_{V_{f}} \mu\Delta\boldsymbol{u}\boldsymbol{\cdot}(\boldsymbol{\nabla} \otimes\boldsymbol{\varPhi})\,\mathrm{d}V. \end{align}
The first term vanishes due to incompressibility. For the second term, we can apply Green's second identity componentwise to move the Laplacian onto the auxiliary potential 
$\boldsymbol {\varPhi }$, which satisfies the Laplace equation. We obtain
\begin{align} -\frac{1}{V} \int_{V_{f}} \mu\Delta\boldsymbol{u}\boldsymbol{\cdot}(\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi})\,\mathrm{d}V & ={-}\frac{1}{V} \int_{A_{fs}} (\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi})^{\mathrm{T}} \boldsymbol{\cdot} \left[\mu\,(\boldsymbol{\nabla}\otimes\boldsymbol{u})^{\mathrm{T}}\boldsymbol{\cdot}\boldsymbol{n}\right]\mathrm{d}A \nonumber\\ &\quad + \frac{1}{V} \int_{A_{fs}} \mu\,\boldsymbol{u} \boldsymbol{\cdot}\left[\boldsymbol{n}\boldsymbol{\cdot}(\boldsymbol{\nabla}\otimes(\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}))\right] \mathrm{d}A. \end{align}
With the no-slip condition 
$\boldsymbol {u}=0$ and the definition of the wall shear stress, we arrive at
\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}
This equation expresses the viscous pressure drag as a weighted integral of the wall shear stress. As the function 
$(\boldsymbol {\nabla }\otimes \boldsymbol {\varPhi })^{\mathrm {T}}$ solely depends on the geometry, we expect that the viscous pressure drag has the same scaling as the wall shear stress and the friction drag.
Appendix B. Asymptotic behaviour of oscillatory flow
This appendix contains a discussion of the behaviour of oscillatory flow in potential flow (§ B.1), at high Womersley numbers (§ B.2) and at low Reynolds numbers (§ B.3). In particular, §§ B.1 and B.2 establish a link between our pressure drag decomposition and the established theory of Johnson et al. (Reference Johnson, Koplik and Dashen1987) for oscillatory porous media flow. Section B.3 generalises the theory of Mei & Auriault (Reference Mei and Auriault1991) and Firdaouss et al. (Reference Firdaouss, Guermond and Le Quéré1997) to oscillatory flow.
B.1. Potential flow
In this section, we derive the potential flow solution in response to a spatially constant macroscopic pressure gradient 
$\boldsymbol {f}$ using the pressure decomposition (2.3). By comparing the boundary value problems for 
$\boldsymbol {\varPhi }$ and 
$p^{(a)}$ (2.3a) we find that 
$p^{(a)} = \boldsymbol {\varPhi }\boldsymbol {\cdot }\boldsymbol {f}$. Since the flow is inviscid, the pressure 
$p^{(v)}$ is zero. It can be shown that for a potential flow the 
$Q$-invariant can be computed as 
$4 Q=-\varDelta |\boldsymbol {u}|^2$. Therefore, we have 
$p^{(c)} = -\frac {1}{2}\rho |\boldsymbol {u}|^2$. Note that 
$p^{(c)}$ satisfies different boundary conditions due to the slip walls where only 
$\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n}=0$. We can now use the momentum equation to determine the velocity:
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}\left(\frac{1}{2}|\boldsymbol{u}|^2\right) ={-}\frac{1}{\rho}\boldsymbol{\nabla}p^{(a)} -\frac{1}{\rho}\boldsymbol{\nabla}p^{(c)} + \frac{1}{\rho} \boldsymbol{f}. \end{equation}
The convective term and 
$\boldsymbol {\nabla }p^{(c)}$ cancel and we are left with
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} ={-}\frac{1}{\rho}\boldsymbol{\nabla}p^{(a)} +\frac{1}{\rho}\boldsymbol{f} = \frac{1}{\rho} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\right) \boldsymbol{\cdot}\boldsymbol{f}. \end{equation}From (B2), the volume-averaged momentum equation in potential flow follows as
\begin{equation} \rho \frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t} = \left\langle\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\right\rangle_{s} \boldsymbol{\cdot} \boldsymbol{f} = \left[\epsilon \boldsymbol{\mathsf{I}} - \frac{1}{V} \int_{V_{f}} \boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\,\mathrm{d}V\right] \boldsymbol{\cdot}\boldsymbol{f} \end{equation}
which we can transform using Gauss's theorem and the periodic boundary conditions of 
$\boldsymbol {\varPhi }$ into
\begin{equation} \rho \frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t} = \left[\epsilon \boldsymbol{\mathsf{I}} - \frac{1}{V} \int_{A_{fs}} \boldsymbol{n}\otimes\boldsymbol{\varPhi}\,\mathrm{d}A\right] \boldsymbol{\cdot}\boldsymbol{f}. \end{equation}With the tensor of added mass (which is symmetric), we can simplify the volume-averaged momentum equation (1.3) to
\begin{equation} \rho\frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t}= \left[\epsilon\,\boldsymbol{\mathsf{I}}-\left(1-\epsilon\right)\boldsymbol{\mathsf{A}}\right] \boldsymbol{\cdot}\boldsymbol{f}. \end{equation}On the other hand, if the porous medium is isotropic, the theory of Johnson et al. (Reference Johnson, Koplik and Dashen1987) gives
\begin{equation} \rho\frac{\alpha_{\infty}}{\epsilon}\,\frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t}= \boldsymbol{f} \end{equation}
in inviscid flow. Consequently, we have 
$\boldsymbol{\mathsf {A}} = {\epsilon }(1-\alpha _{\infty }^{-1})/{(1-\epsilon )}\boldsymbol{\mathsf {I}}$ with the high-frequency limit of the dynamic tortuosity 
$\alpha _{\infty }$ (Johnson et al. Reference Johnson, Koplik and Dashen1987).
B.2. Behaviour in the high Womersley number limit
In this section, we show how the pressure decomposition can be used to derive the high-frequency asymptotics of oscillatory flow in a porous medium (Johnson et al. Reference Johnson, Koplik and Dashen1987). These can be written in the time domain as
\begin{equation} \rho \frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t} ={-} \rho\sqrt{\nu}\,\frac{2}{\varLambda} \int_{0}^{t} \frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}\tau} \,\frac{1}{\sqrt{{\rm \pi}(t-\tau)}}\,\mathrm{d}\tau + \frac{\epsilon}{\alpha_{\infty}}\boldsymbol{f} . \end{equation}We begin the derivation from the volume-averaged momentum equation (1.3) in which we insert the decomposition (2.11) to get
\begin{align} \rho\frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t}={-}\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 + \underbrace{\frac{1}{V}\int_{V_{f}}\boldsymbol{\varPhi}\,2 \rho Q \,\mathrm{d}V}_{{\approx} 0} + \left[\epsilon\boldsymbol{\mathsf{I}}- \left(1-\epsilon\right) \boldsymbol{\mathsf{A}}\right] \boldsymbol{\cdot} \boldsymbol{f}. \end{align}For linear flow, the convective pressure drag can be neglected as it contains the square of the velocity.
In the high-frequency limit, the flow has a boundary layer character and the boundary layer is locally identical to a Stokes boundary layer (Schlichting & Gersten Reference Schlichting and Gersten2017, pp. 352f, pp. 126f). The wall shear stress in the Stokes boundary layer can be written as
\begin{equation} \boldsymbol{\tau}_{w} = \rho\sqrt{\nu} \int_{0}^{t} \left.\frac{\partial \boldsymbol{u}_p}{\partial \tau}\right\vert_{w} \,\frac{1}{\sqrt{{\rm \pi}(t-\tau)}}\,\mathrm{d}\tau, \end{equation}
where 
$\boldsymbol {u}_{p}$ is the potential flow velocity in the core flow. Combining (B2) and (B5), we can establish a one-to-one correspondence between the velocity field in potential flow and its superficial average:
\begin{equation} \rho \frac{\partial \boldsymbol{u}_{p}}{\partial t} = \left(\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\right) \boldsymbol{\cdot} \left[\epsilon\,\boldsymbol{\mathsf{I}}-\left(1-\epsilon\right)\boldsymbol{\mathsf{A}}\right]^{{-}1} \boldsymbol{\cdot}\rho\frac{\mathrm{d}\langle\boldsymbol{u}_{p}\rangle_{s}}{\mathrm{d}t} . \end{equation}With this relation, the wall shear stress can be expressed in terms of the superficial velocity of the potential flow:
\begin{equation} \boldsymbol{\tau}_{w} = \rho\sqrt{\nu} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{\nabla}\otimes\boldsymbol{\varPhi}\right) \boldsymbol{\cdot}\left[\epsilon\,\boldsymbol{\mathsf{I}}-\left(1-\epsilon\right)\boldsymbol{\mathsf{A}}\right]^{{-}1} \boldsymbol{\cdot} \int_{0}^{t} \frac{\mathrm{d}\langle\boldsymbol{u}_{p}\rangle_{s}}{\mathrm{d}\tau} \,\frac{1}{\sqrt{{\rm \pi}(t-\tau)}}\,\mathrm{d}\tau . \end{equation}Using (2.10), we can compute the total viscous drag force as
\begin{equation} -\frac{1}{V}\int_{A_{fs}} (p^{(v)}\,\boldsymbol{n}+\boldsymbol{\tau}_{w})\,\mathrm{d}A ={-} \rho\sqrt{\nu}\,2\,\boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\cdot} \int_{0}^{t} \frac{\mathrm{d}\langle\boldsymbol{u}_{p}\rangle_{s}}{\mathrm{d}\tau} \,\frac{1}{\sqrt{{\rm \pi}(t-\tau)}}\,\mathrm{d}\tau \end{equation}with the tensor
\begin{equation} 2\boldsymbol{\mathsf{L}}^{{-}1} =\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 \boldsymbol{\cdot}\left[\epsilon\,\boldsymbol{\mathsf{I}}-\left(1-\epsilon\right)\boldsymbol{\mathsf{A}}\right]^{{-}1} . \end{equation}Finally, when the potential flow velocity is replaced with the actual fluid velocity, we obtain the volume-averaged momentum equation
\begin{equation} \rho \frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}t} ={-} \rho\sqrt{\nu}\,2\,\boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\cdot} \int_{0}^{t} \frac{\mathrm{d}\!\left\langle\boldsymbol{u}\right\rangle_{s}}{\mathrm{d}\tau} \,\frac{1}{\sqrt{{\rm \pi}(t-\tau)}}\,\mathrm{d}\tau + \left[\epsilon\,\boldsymbol{\mathsf{I}}-\left(1-\epsilon\right)\boldsymbol{\mathsf{A}}\right] \boldsymbol{\cdot}\boldsymbol{f} . \end{equation}Comparing this result with the high-frequency asymptotics of Johnson et al. (Reference Johnson, Koplik and Dashen1987) given in (B7), it is readily apparent that the former is just a tensorial generalisation of the latter.
B.3. Behaviour at finite Reynolds numbers
In this section, we derive the 
$Re^3$ dependence of the first nonlinear correction to the linear drag behaviour from the fore–aft symmetry of the hexagonal sphere pack for oscillatory flow in the 
$x$-direction. The derivation is based on our new representation of the drag in the volume-averaged momentum equation (2.11) and assumes a macroscopic pressure gradient along the 
$x$-direction. This extends the results of Mei & Auriault (Reference Mei and Auriault1991) and Firdaouss et al. (Reference Firdaouss, Guermond and Le Quéré1997) for steady flow at finite Reynolds numbers (see (1.7)) to oscillatory flow.
The viscous pressure drag in the 
$x$-direction is given by (2.10)
\begin{equation} f_{px}^{(v)}={-}\frac{1}{V}\int_{A_{fs}} p^{(v)}\,n_x\,\mathrm{d}A = \frac{1}{V}\int_{A_{fs}} \boldsymbol{\nabla}\varPhi_x\boldsymbol{\cdot} \boldsymbol{\tau}_{w}\,\mathrm{d}A , \end{equation}the friction drag is given by the integral of the wall shear stress
\begin{equation} f_{\tau_{w}x}={-}\frac{1}{V}\int_{A_{fs}} \tau_{{w}x}\,\mathrm{d}A \end{equation}and the convective pressure drag is given by (2.7c)
\begin{equation} -\frac{1}{V}\int_{A_{fs}} p^{(c)}\,n_x\,\mathrm{d}A = \frac{1}{V}\int_{V_{f}} \varPhi_x\,2 \rho Q \,\mathrm{d}V . \end{equation}
The auxiliary potential 
$\varPhi _x$ is fore–aft antisymmetric with respect to the planes 
$x=0, x=d/2, x=d, \ldots$, i.e. an odd function with respect to 
$x$. Therefore, the partial derivative 
${\partial \varPhi _x}/{\partial x}$ is an even function whereas 
${\partial \varPhi _x}/{\partial y}$ and 
${\partial \varPhi _x}/{\partial z}$ are odd functions. Thus, the friction and viscous pressure drag are generated by the even part of 
$\tau _{{w}x}$ and by the odd part of 
$\tau _{{w}y}$ and 
$\tau _{{w}z}$; the convective pressure drag is generated by the odd part of the 
$Q$-invariant. Below we discuss how these parts depend on the Reynolds number.
Like in our previous work (Unglehrt & Manhart Reference Unglehrt and Manhart2022a), we decompose the velocity field into a symmetric part 
$\boldsymbol {u}_{sym}=\frac {1}{2}(\boldsymbol {u}+\mathcal {S}\boldsymbol {u})$ and an antisymmetric part 
$\boldsymbol {u}_{anti}=\frac {1}{2}(\boldsymbol {u}-\mathcal {S}\boldsymbol {u})$ with respect to the fore–aft symmetry
\begin{equation} \mathcal{S}\boldsymbol{u}(\boldsymbol{x},t) = \begin{bmatrix} u(2 d-x,y,z,t) \\ -v(2 d-x,y,z,t) \\ -w(2 d-x,y,z,t) \end{bmatrix} \end{equation}
with respect to the plane 
$x=d$. In laminar flow, the velocity field is 
$d$-periodic in the 
$x$-direction. Consequently, this symmetry operation also expresses the fore–aft symmetries with respect to the planes 
$x=0$, 
$x=d/2$, 
$x=3/2 \,d$ and 
$x=2\,d$.
The wall shear stress points in the direction of the velocity vector near the wall. Consequently, the 
$x$-component of the wall shear stress of the symmetric part 
$\boldsymbol {u}_{sym}$ of the velocity field is an even function whereas the 
$y$- and 
$z$-components are odd functions. The wall shear stress components of the antisymmetric part 
$\boldsymbol {u}_{anti}$ of the velocity field behave in the opposite way. So we find that the friction and viscous pressure drag arise solely from the symmetric part 
$\boldsymbol {u}_{sym}$.
Using the decomposition of the velocity field the 
$Q$-invariant can be rewritten as
\begin{align}
Q = \underbrace{-\tfrac{1}{2}
(\boldsymbol{\nabla}\otimes\boldsymbol{u}_{sym}):(\boldsymbol{\nabla}\otimes\boldsymbol{u}_{sym})^{\mathrm{T}}}_{\textit{fore{--}aft symmetric}}
\underbrace{-\vphantom{\frac{1}{2}}
(\boldsymbol{\nabla}\otimes\boldsymbol{u}_{sym}):(\boldsymbol{\nabla}\otimes\boldsymbol{u}_{anti})^{\mathrm{T}}}_{\textit{fore{--}aft antisymmetric}}
\underbrace{-\tfrac{1}{2}
(\boldsymbol{\nabla}\otimes\boldsymbol{u}_{anti}):(\boldsymbol{\nabla}\otimes\boldsymbol{u}_{anti})^{\mathrm{T}}}_{\textit{fore{--}aft symmetric}}.
\end{align}
We find that the quadratic contributions in 
$\boldsymbol {u}_{sym}$ and 
$\boldsymbol {u}_{anti}$ lead to a fore–aft symmetric distribution of 
$Q$ and hence do not cause any convective pressure drag. On the other hand, the interaction between 
$\boldsymbol {u}_{sym}$ and 
$\boldsymbol {u}_{anti}$ is fore–aft antisymmetric and can cause a convective pressure drag.
For small Reynolds numbers, the velocity field can be described by as the sum of the velocity field in linear flow and corrections proportional to powers of the Reynolds number:
$$\begin{gather} \boldsymbol{u}_{sym} = \boldsymbol{u}_{1|{sym}}\,Re + \boldsymbol{u}_{2|{sym}} \,Re^2 + {O}(Re^3), \end{gather}$$
$$\begin{gather}\boldsymbol{u}_{anti} = \boldsymbol{u}_{1|{anti}}\, Re + \boldsymbol{u}_{2|{anti}}\,Re^2 + {O}(Re^3). \end{gather}$$
Since the velocity field in linear flow is fore–aft symmetric (Unglehrt & Manhart Reference Unglehrt and Manhart2022a), the antisymmetric first-order contribution 
$\boldsymbol {u}_{1|{anti}}$ is zero. The self-interaction of the symmetric first-order contribution 
$(\boldsymbol {u}_{1|{sym}}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {u}_{1|{sym}}$ creates the antisymmetric second-order contribution 
$\boldsymbol {u}_{2|{anti}}$ whereas the symmetric second-order contribution 
$\boldsymbol {u}_{2|{sym}}$ is zero. Then, we have that the symmetric part 
$\boldsymbol {u}_{sym}$ is proportional to 
$Re$ and causes a friction drag and viscous pressure drag proportional to 
$Re$ with a higher-order contribution of order 
$Re^3$. The antisymmetric part 
$\boldsymbol {u}_{anti}$ is proportional to 
$Re^2$ and does not cause any friction and viscous pressure drag. The convective pressure drag arises from the part of the 
$Q$-invariant due to the interaction of 
$\boldsymbol {u}_{sym}$ and 
$\boldsymbol {u}_{anti}$ and is therefore proportional to 
$Re^3$. In conclusion, like in steady flow (Mei & Auriault Reference Mei and Auriault1991) the drag in oscillatory flow at small Reynolds numbers consists of a linear and a cubic part in 
$Re$.
Appendix C. Grid resolution of the simulations
In this appendix, we discuss the grid resolution of the simulation cases LF5, LF6, MF5, MF6 and HF5. For the other oscillatory cases LF1–LF4, MF1–MF4, HF1–HF4 and for the steady cases, convergence with respect to grid resolution was demonstrated in the previous publications, Unglehrt & Manhart (Reference Unglehrt and Manhart2022a) and Sakai & Manhart (Reference Sakai and Manhart2020), respectively.
C.1. Estimate of the required grid resolution
For turbulent flow driven by a constant pressure gradient, the required grid resolution can be estimated following Finn (Reference Finn2013) and He et al. (Reference He, Apte, Finn and Wood2019). It is assumed that a grid spacing in wall units
\begin{equation} {\Delta x}^+=\frac{u_\tau\,{\Delta x}}{\nu}\approx 1\unicode{x2013}3 \end{equation}
is necessary to resolve all scales in the flow, where the friction velocity 
$u_\tau =\sqrt {\langle \tau _{{w}x}\rangle _{A_{fs}}/\rho }$ is defined in terms of the wall shear stress 
$\langle \tau _{{w}x}\rangle _{A_{fs}}$ averaged over the fluid–solid interface. He et al. (Reference He, Apte, Finn and Wood2019) approximate the average wall shear stress as a fraction 
$\beta \approx 0.25$ of the total stress 
$\langle \sigma _{{w}x}\rangle _{A_{fs}}$ which they find from equilibrium as
\begin{equation} \langle \sigma_{{w}x}\rangle_{A_{fs}} = f_x\,\frac{d}{6}\,\frac{\epsilon}{1-\epsilon} . \end{equation}Combining (C1) and (C2), the required grid resolution for turbulent flow at a given Hagen number can be estimated as
\begin{equation} \frac{d}{{\Delta x}}=\frac{1}{{\Delta x}^+ }\sqrt{\frac{\beta}{6}\,\frac{\epsilon}{1-\epsilon}\,Hg} . \end{equation}
For the cases LF6 and MF6, the acceleration is close to zero when the convective pressure drag is large. Therefore, it seems plausible that these cases behave similar to a flow with a constant pressure gradient. Taking a dimensionless grid spacing 
${\Delta x}^+=1$ and setting 
$\beta =0.2$ (which was taken from the momentum budgets in the figures 12 and 13), the estimate results in a required grid resolution of 342 cpd for the cases LF6 and MF6 (
$Hg=10^7$). Consequently, the employed grid resolution of 384 cpd for the cases LF6 and MF6 seems to be sufficient. For the case HF5, the estimate is not applicable, as the flow is far from an equilibrium with the imposed pressure gradient and the wall shear stress is out of phase with the convective pressure drag (figure 14).
C.2. Grid study
In this section, we present a grid study for the cases LF5, LF6, MF5, MF6 and HF5. For each case the simulations were conducted at the resolutions 48 cpd, 96 cpd, 192 cpd and 384 cpd; for the case HF5 an additional simulation at 768 cpd was performed.
For consistency, the discretisation error is assessed using the same procedure as in our previous publication (Unglehrt & Manhart Reference Unglehrt and Manhart2022a). We first consider the Reynolds number based on the maximum superficial velocity in the last cycle and the sphere diameter, which is defined in (1.4). As can be seen in table 4, for every case the Reynolds numbers differ less than 
$1\,\%$ between the two finest grid resolutions. We then consider the space–time 
$L^2$-norm of the velocity field over the last simulated period,
\begin{equation} \left\Vert\boldsymbol{u}\right\Vert_{L^2}^2 = \int_{V_{f}}\int_{T} \left\vert\boldsymbol{u}\right\vert^2\,\mathrm{d}t\,\mathrm{d}V , \end{equation}
corresponding to the signal energy of the velocity field. For all cases the relative difference of the space–time 
$L^2$-norm between the second-finest grid to the finest grid is below 
$1.8\,\%$ (cf. table 4). Consequently, we consider the simulations at the finest grid resolution as converged.
Grid convergence of the velocity field 
$\boldsymbol {u}(\boldsymbol {x},t)$ in steady oscillation. The relative error in 
$\Vert \boldsymbol {u}\Vert _{L^2}^2$ is defined as 
$e_{res} = ( \Vert \boldsymbol {u}_{res}\Vert _{L^2}^2-\Vert \boldsymbol {u}_{384}\Vert _{L^2}^2 ) / \Vert \boldsymbol {u}_{384}\Vert _{L^2}^2$ and as 
$e_{res} = ( \Vert \boldsymbol {u}_{res}\Vert _{L^2}^2-\Vert \boldsymbol {u}_{768}\Vert _{L^2}^2 ) / \Vert \boldsymbol {u}_{768}\Vert _{L^2}^2$ for HF5. The Reynolds number 
$Re$ is defined according to (1.4).
Open access
