1. Introduction
Heave plates are integral components in offshore engineering, particularly for enhancing the stability of structures such as platforms, buoys and floating wind turbines that are exposed to the dynamic forces of ocean waves (Ezoji, Shabakhty & Tao Reference Ezoji, Shabakhty and Tao2022). These broad, flat plates are installed at the base of floating structures to act as hydrodynamic dampers, effectively increasing the added mass and damping properties of the system. This increase often leads to reduced vertical motion, or heave, thereby enhancing the operational stability and safety of these structures. The practical implications of this enhanced stability are substantial, offering improved durability and reliability of offshore operations in challenging marine environments. Furthermore, in the field of renewable energy, heave plates are especially critical in the design of two-body point absorber wave energy converters, where they not only stabilise the device against the unpredictable nature of ocean waves but also influence the energy conversion efficiency by modulating the hydrodynamic interaction between wave forces and the converter itself (Brown & Thomson Reference Brown and Thomson2015; Tamakuwala et al. Reference Tamakuwala, Usman, Tom and Masoud2025). While the importance of heave plates is well recognised, the majority of the hydrodynamic research to date has concentrated on impermeable heave plates (see, e.g. He Reference He2003). Therefore, the dynamics of porous heave plates, which could offer additional benefits (Downie et al. Reference Downie, Wang and Graham2000b ), remains underexplored. Porous plates have the potential to modify fluid flow patterns around the structure to enhance its performance under specific conditions. The interaction between water flow through the pores and the plate adds a layer of complexity absent in impermeable designs, potentially leading to more desired dynamic behaviour.
Among the limited studies that have examined the hydrodynamics of oscillating porous plates, the majority are experimental (Downie et al. Reference Downie, Graham, Hall, Incecik and Nygaard2000a ; Vu, Chenu & Thiagarajan Reference Vu, Chenu and Thiagarajan2004; Chua et al. Reference Chua, Clelland, Huang and Sworn2005; Molin, Remy & Rippol Reference Molin, Remy and Rippol2007; Tao & Dray Reference Tao and Dray2008; Wadhwa, Thiagarajan & Pistani Reference Wadhwa, Thiagarajan and Pistani2009; Wadhwa Reference Wadhwa2011; An & Faltinsen Reference An and Faltinsen2013; Li et al. Reference Li, Liu, Zhao and Teng2013; Mentzoni & Kristiansen Reference Mentzoni and Kristiansen2020; Ezoji et al. Reference Ezoji, Shabakhty and Tao2022), focusing on how amplitude variations affect the force coefficients during heave oscillations. Additionally, a few theoretical investigations utilised potential flow theory, integrated with a quadratic discharge law that governs the pressure drop through the pores, to estimate the force coefficients of perforated oscillating plates (Molin Reference Molin1993, Reference Molin2001; Molin et al. Reference Molin, Remy and Rippol2007; Liu et al. Reference Liu, Li, Li and He2011; Molin Reference Molin2011; An & Faltinsen Reference An and Faltinsen2013). Complementing these, several recent computational studies have simulated the flow around and through perforated plates undergoing heave oscillations (Nokob & Yeung Reference Nokob and Yeung2018; Mentzoni & Kristiansen Reference Mentzoni and Kristiansen2019a , Reference Mentzoni and Kristiansenb , Reference Mentzoni and Kristiansen2020). The aforementioned studies generally find that porous plates provide additional damping at small amplitudes, attributing this effect to flow separation and vortex shedding. However, contrasting these findings, we reported increased damping in our previous theoretical investigation for porous objects, even in the absence of such phenomena (Jafari Kang et al. Reference Jafari Kang, Dehdashti, Vandadi and Masoud2019), see also Sherwood (Reference Sherwood2020). In that study, we analysed the vibrations of porous cylinders at very small amplitudes, employing the Brinkman–Debye–Bueche model for permeable objects. This model accounts for the collective effect of pores through permeability (an indirect measure of porosity) rather than resolving individual pore flows. Notably, we discovered that the additional damping results from an increased phase shift between force and displacement.
The lack of a consensus on the mechanisms responsible for the change in the hydrodynamics of heave plates due to perforation underscores the need for further theoretical studies to elucidate the precise effects of perforation across a wide range of frequencies and amplitudes. To address this issue, here, we theoretically examine the elementary problem of an annular disk oscillating with a small amplitude in an unbounded fluid domain. We select this geometry because it represents the simplest form of a perforated heave plate with clearly defined porosity, thereby eliminating the need for permeability-based continuum modelling to describe the flow through the disk. We employ a combination of the superposition principle, the reciprocal theorem and perturbation expansion to develop semi-analytical formulas for the variations in the added mass and damping coefficients of the annular disk as functions of oscillation frequency and pore radius. These formulations are further simplified in the asymptotic limits of low and high frequencies. Our solution strategy entails breaking down the problem into two more manageable sub-problems: one involving a solid disk vibrating in heave motion, and the other involving pressure-driven flow through a hole in an infinite plate. We apply Tranter’s method to solve the dual integral equations that arise in each sub-problem. Our findings indicate that annular disks indeed provide additional damping (in excess of 60 % in extreme limits) as compared with their solid counterparts over specific ranges of oscillation frequency and pore radius. Additionally, we observe a reduction in the disk’s added mass coefficient when it is made porous, independent of the oscillation frequency. Crucially, the small-amplitude assumption ensures that our predicted variations of the added mass and damping coefficients with pore radius and frequency are not caused by flow separation and vortex shedding, but instead result from changes in the phase relation between force and displacement, as discussed in § 3. For validation, we also solve the governing equations using a second-order finite-element method and demonstrate that our theoretical formulations retain their accuracy across a broad spectrum of parameters.
In what follows, we first describe the problem statement and our solution (§ 2). The results of our calculations are presented and discussed next (§ 3), followed by a short summary and concluding remarks (§ 4).
Figure 1. Schematics of an annular disk with outer radius
$R$
and inner-to-outer radius ratio of
$\varepsilon$
undergoing heave oscillations with amplitude
$A$
at frequency
$\omega$
in the
$z$
direction within an unbounded fluid domain.
2. Problem formulation and solution
Consider an annular disk of outer radius
$R$
and inner-to-outer radius ratio of
$ \varepsilon$
, oscillating sinusoidally with amplitude
$A$
at frequency
$\omega$
within an unbounded, incompressible, Newtonian fluid of density
$\varrho$
and viscosity
$\mu$
(see figure 1). Let
$(\rho ,\theta ,z)$
represent a cylindrical coordinate system positioned at the centre of the disk, with the unit vector
$\boldsymbol e_z$
normal to the broad side of the disk pointing upward, as shown in figure 1.
Suppose the disk’s thickness is negligibly small and the amplitude of oscillations is much smaller than the outer radius of the disk (
$A \ll R$
), such that the convective inertia term in the Navier–Stokes equations can be neglected relative to the unsteady inertia term. This assumption may impose a stricter upper limit on the oscillation amplitude at high frequencies, since the associated error is expected to grow with increasing
$\omega$
(see, e.g. Aureli, Basaran & Porfiri Reference Aureli, Basaran and Porfiri2012). Under these conditions, the dynamics of the fluid flow through and around the disk is governed by the linearised Navier–Stokes equations, which are expressed in dimensionless form below, along with the associated boundary conditions
\begin{align} \eta ^2 \frac {\partial { \tilde {\boldsymbol u} } }{\partial { \tilde t } } = {\boldsymbol{\nabla }} \boldsymbol{\cdot }{ \tilde {\boldsymbol \sigma } } = {\nabla} ^2 { \tilde {\boldsymbol u} } - {\boldsymbol{\nabla }} { \tilde p } \quad \text{and} \quad {\boldsymbol{\nabla }} \boldsymbol{\cdot }{ \tilde {\boldsymbol u} } = 0, \end{align}
with the conditions
\begin{align} { \tilde {\boldsymbol u} } &= \cos { \tilde t } \, { \boldsymbol{e }}_z \quad \text{for} \quad z = 0 \quad \text{and} \quad \varepsilon \leqslant \rho \leqslant 1, \notag \\ { \tilde {\boldsymbol u} } &\to {\boldsymbol{0}} \quad \text{as} \quad \sqrt {\rho ^2 + z^2} \to \infty , \notag \\ { \tilde p } &= 0 \quad \text{for} \quad z = 0 \quad \text{when} \quad 0 \leqslant \rho \lt \varepsilon \quad \text{or} \quad \rho \gt 1. \end{align}
The symbols
$ { \tilde t }$
,
$ { \tilde {\boldsymbol u} }$
,
$ { \tilde p }$
and
$ { \tilde {\boldsymbol \sigma } }$
in the above equations denote the dimensionless time, velocity field, pressure field and stress tensor, respectively. Also,
$\eta ^2 = \varrho R^2 \omega / \mu$
is the oscillatory Reynolds number that measures the time scale of vorticity diffusion relative to the oscillation period. Here, and throughout the article, the length, time, velocity and stress are non-dimensionalised by
$R$
,
$1 / \omega$
,
$A \omega$
and
$\mu A \omega / R$
, respectively.
Considering the sinusoidal nature of the disk’s oscillations, it is mathematically advantageous to express the velocity, pressure and stress as real parts of complex quantities, specifically
$ { \tilde {\boldsymbol u} } = \mathrm{Re} [ {\boldsymbol{u}} \exp ( i { \tilde t } ) ]$
,
$ { \tilde p } = \mathrm{Re}[ p \exp ( i { \tilde t } ) ]$
and
$ { \tilde {\boldsymbol \sigma } } = \mathrm{Re}[ { \boldsymbol \sigma } \exp ( i { \tilde t } ) ]$
, with
$i^2 = -1$
. This transformation simplifies (2.1) and (2.2) to
\begin{align} i \eta ^2 {\boldsymbol{u}} = { \boldsymbol{\nabla }} \boldsymbol{\cdot }{ \boldsymbol \sigma } = {\nabla} ^2 {\boldsymbol{u}} - { \boldsymbol{\nabla }} p \qquad \text{and} \qquad { \boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol{u}} = 0, \end{align}
with the boundary conditions
\begin{align} {\boldsymbol{u}} &= { \boldsymbol{e }}_z \quad \text{for} \quad z = 0 \quad \text{and} \quad \varepsilon \leqslant \rho \leqslant 1, \notag \\ {\boldsymbol{u}} &\to {\boldsymbol{0}} \quad \text{as} \quad \sqrt {\rho ^2 + z^2} \to \infty , \notag \\ p &= 0 \quad \text{for} \quad z = 0 \quad \text{when} \quad 0 \leqslant \rho \lt \varepsilon \quad \text{or} \quad \rho \gt 1. \end{align}
Given the boundary conditions of this boundary-value problem, it follows that the velocity field is symmetric and the pressure field is antisymmetric with respect to the
$z = 0$
plane. This ensures that the normal traction force remains continuous across the disk’s two faces due to the change in the surface normal direction above and below the disk. Therefore, it is sufficient to consider the problem only in the upper half-space (
$z \geqslant 0$
), and express the force exerted on the disk by the fluid in the frequency domain as
\begin{align} {\boldsymbol{F}} = 2 \int _{S_d} { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol \sigma } \; {\textrm { d}S} = - 2 \int _{S_d} p \; {\textrm { d}S} \, { \boldsymbol{e }}_z = F \, { \boldsymbol{e }}_z, \end{align}
where
$ { \boldsymbol n }$
is the unit normal vector and
$S_d$
denotes the upper surface of the annular disk. The dimensionless added mass,
$a$
, and damping coefficients,
$b$
, of the disk are then derived from the following relationship:
\begin{align} b + i a = - \frac {F}{\eta ^2}, \end{align}
where the added mass and damping coefficients are non-dimensionalised by
$\varrho R^3$
and
$\varrho R^3 \omega$
, respectively.
A solution strategy for solving the problem at hand and calculating
$a$
and
$b$
is to apply Hankel transforms to convert (2.3) and (2.4) to triple integral equations involving an unknown function which may be determined using the procedure outlined by Sneddon (Reference Sneddon1966). Here, however, we adopt an alternative approach that involves splitting the original problem into two sub-problems, labelled with subscripts 0 and 1
\begin{align} {\boldsymbol{u}} &= {\boldsymbol{u}}_0 + {\boldsymbol{u}}_1, \end{align}
\begin{align} p &= p_0 + p_1, \end{align}
where
$( {\boldsymbol{u}}_0 , p_0 )$
and
$( {\boldsymbol{u}}_1 , p_1 )$
, respectively, satisfy
\begin{align} i \eta ^2 {\boldsymbol{u}}_0 = { \boldsymbol{\nabla }} \boldsymbol{\cdot }{ \boldsymbol \sigma } _0 = {\nabla} ^2 {\boldsymbol{u}}_0 - { \boldsymbol{\nabla }} p_0 \quad \text{and} \quad { \boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol{u}}_0 = 0, \end{align}
with
\begin{align} {\boldsymbol{u}}_0 &= { \boldsymbol{e }}_z \quad \text{for} \quad z = 0 \quad \text{and} \quad 0 \leqslant \rho \leqslant 1, \notag \\ {\boldsymbol{u}}_0 &\to {\boldsymbol{0}} \quad \text{as} \quad \sqrt {\rho ^2 + z^2} \to \infty , \notag \\ p_0 &= 0 \quad \text{for} \quad z = 0 \quad \text{and} \quad \rho \gt 1, \end{align}
and
\begin{align} i \eta ^2 {\boldsymbol{u}}_1 = { \boldsymbol{\nabla }} \boldsymbol{\cdot }{ \boldsymbol \sigma } _1 = {\nabla} ^2 {\boldsymbol{u}}_1 - { \boldsymbol{\nabla }} p_1 \quad \text{and} \quad { \boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol{u}}_1 = 0, \end{align}
with
\begin{align} {\boldsymbol{u}}_1 &= {\boldsymbol{0}} \quad \text{for} \quad z = 0 \quad \text{and} \quad \varepsilon \leqslant \rho \leqslant 1, \notag \\ {\boldsymbol{u}}_1 &\to {\boldsymbol{0}} \quad \text{as} \quad \sqrt {\rho ^2 + z^2} \to \infty , \notag \\ p_1 &= - p_0 \quad \text{for} \quad z = 0 \quad \text{and} \quad 0 \leqslant \rho \lt \varepsilon , \notag \\ p_1 &= 0 \quad \text{for} \quad z = 0 \quad \text{and} \quad \rho \gt 1. \end{align}
Below, we present a force decomposition that stems from the above superposition, followed by the solution of these sub-problems.
2.1. Force decomposition via reciprocal theorem
Applying the reciprocal theorem (Masoud & Stone Reference Masoud and Stone2019) to the velocity–stress pairs
$({\boldsymbol{u}} , { \boldsymbol \sigma } )$
and
$({\boldsymbol{u}}_0 , { \boldsymbol \sigma } _0)$
yields
\begin{align} \int _S { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol \sigma } \boldsymbol{\cdot }{\boldsymbol{u}}_0 \; {\textrm { d}S} = \int _S { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol \sigma } _0 \boldsymbol{\cdot }{\boldsymbol{u}} \; {\textrm { d}S} = \int _S { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol \sigma } _0 \boldsymbol{\cdot }\left ({\boldsymbol{u}}_0 + {\boldsymbol{u}}_1 \right ) \; {\textrm { d}S}, \end{align}
where
$S$
denotes the surfaces enclosing the upper half of the fluid domain. The integrals over the bounding surface at infinity vanish due to the decay rate of the velocity and stress fields. After implementing the boundary conditions, (2.13) simplifies to
\begin{align} \int _{S_d} { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol \sigma } \; {\textrm { d}S} \boldsymbol{\cdot }{ \boldsymbol{e }}_z = \int _{S_{d,0}} { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol \sigma } _0 \; {\textrm { d}S} \boldsymbol{\cdot }{ \boldsymbol{e }}_z - \int _{S_{d,1}} p_0 \, {\boldsymbol{u}}_1 \; {\textrm { d}S} \boldsymbol{\cdot }{ \boldsymbol{e }}_z, \end{align}
where
$S_{d,0}$
and
$S_{d,1}$
represent areas of
$z = 0$
plane confined to
$0 \leqslant \rho \leqslant 1$
and
$0 \leqslant \rho \lt \varepsilon$
, respectively. Let us define
\begin{align}& F_0 = 2 \int _{S_{d,0}} { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol \sigma } _0 \; {\textrm { d}S} \boldsymbol{\cdot }{ \boldsymbol{e }}_z = - 2 \int _{S_{d,0}} p_0 \; {\textrm { d}S}, \end{align}
\begin{align}&\qquad\quad F_1 =- 2 \int _{S_{d,1}} p_0 \left ( {\boldsymbol{u}}_1 \boldsymbol{\cdot }{ \boldsymbol{e }}_z \right ) \; {\textrm { d}S}. \end{align}
Referring back to (2.5), the total fluid force on the disk can now be expressed as the sum of contributions from the two sub-problems as
\begin{align} F = F_0 + F_1. \end{align}
Consequently, the added mass and damping coefficients can be decomposed as
\begin{align} a &= a_0 + a_1, \end{align}
\begin{align} b &= b_0 + b_1. \end{align}
Here,
$a_0$
and
$b_0$
represent the contributions from the primary sub-problem, while
$a_1$
and
$b_1$
account for the additional effects introduced by the secondary sub-problem.
2.2. First sub-problem
For completeness, we reiterate the solution steps here. By applying Hankel transforms to (2.9), solving the resulting equations using the boundary conditions specified in (2.10) and then performing inverse Hankel transforms, the following expressions can be derived for the velocity and pressure fields:
\begin{align} {\boldsymbol{u}}_0 \boldsymbol{\cdot }{ \boldsymbol{e }}_\rho & = \int _0^\infty k A(k) \big( e^{- k z} - e^{- \xi z}\big) J_1(k \rho ) \; \textrm { d}k, \end{align}
\begin{align} {\boldsymbol{u}}_0 \boldsymbol{\cdot }{ \boldsymbol{e }}_z & = \int _0^\infty k A(k) \left (e^{- k z} - \frac {k}{\xi } e^{- \xi z}\right ) J_0(k \rho ) \; \textrm { d}k, \end{align}
\begin{align} p_0 & = i \eta ^2 \int _0^\infty A(k) \, e^{- k z} \, J_0(k \rho ) \; \textrm { d}k, \end{align}
where
$\xi ^2 = k^2 + i \eta ^2$
,
$J$
is the Bessel function of the first kind, and the unknown function
$A(k)$
satisfies
\begin{align} \int _0^\infty k A(k) \left (1 - \frac {k}{\xi }\right ) J_0(k r) \; \textrm { d}k &= 1 \quad \text{for} \quad 0 \leqslant \rho \leqslant 1, \end{align}
\begin{align} \int _0^\infty A(k) \, J_0(k r) \; \textrm { d}k & = 0 \quad \text{for} \quad \rho \gt 1. \end{align}
Following Tranter’s method (Tranter Reference Tranter1966),
$A(k)$
can be assumed to take the form of
\begin{align} A(k) = \sqrt { \frac {2 k}{\pi } } \, \sum _{m = 0}^\infty \alpha _m \, J_{2m+1/2}(k) \quad \text{for} \quad m = 0, 1, 2 {\cdots} . \end{align}
This representation reduces the above dual integral equations to a system of linear algebraic equations for
$\alpha _m$
\begin{align} \sum _{m = 0}^\infty \alpha _m \int _0^\infty k \left (1 - \frac {k}{\xi } \right ) J_{2m+1/2}(k) \, J_{2n+1/2}(k) \; \textrm { d}k = \delta _{0n} \quad \text{for} \quad n = 0, 1, {\cdots} , \end{align}
with
$\delta _{ij}$
being the Kronecker delta function. Equation (2.26) can be solved numerically using standard approaches. Substituting
$A(k)$
from (2.25) into (2.22) and performing surface integration of the subsequent equation leads to (after some simplifications) the following formulas for the pressure and force on the disk:
\begin{align} \hat {p}_0 &= \frac {2 i \eta ^2}{\sqrt {\pi }} \sum _{m = 0}^\infty \frac {\varGamma (m+1) \, \alpha _m}{\varGamma (m+1/2)} \, _2F_1( m + 1 , 1 / 2 - m ; 1 ; \rho ^2 ), \end{align}
\begin{align} F_0 &= - 2 \int _{S_{d,0}} p_0 \; {\textrm { d}S} = - 8 i \eta ^2 \alpha _0, \end{align}
where
$_2F_1$
is the hypergeometric function and the hat symbol indicates evaluation at the surface of the disk.
2.3. Second sub-problem
This problem describes the fluid flow through and around a stationary annular disk driven by the pressure difference across the disk required to cancel out the contribution from the first sub-problem over the disk’s orifice. Solving for
$( {\boldsymbol{u}}_1 , p_1 )$
without constraints on the pore size of the disk is challenging. Hence, we aim to develop an asymptotic solution in the limit of small
$ \varepsilon$
. As demonstrated in the forthcoming section (§ 3), the approximation error of our solution remains minimal even when
$ \varepsilon$
is as large as
$ 1 / 2$
and, in some cases, beyond.
We use a Taylor series expansion to approximate the pressure distribution on the disk from the prior sub-problem as
\begin{align} \hat {p}_0= \hat {p}_0 \big |_{\rho = 0} + \left .\frac {{\rm d} \hat {p}_0}{{\rm d} \rho } \right |_{\rho = 0} \rho + \frac {1}{2}\left . \frac {{\rm d}^2 \hat {p}_0}{{\rm d} \rho ^2} \right |_{\rho = 0} \rho ^2 + {\cdots} , \end{align}
where
\begin{align} \hat {p}_0 \big |_{\rho = 0} = \frac {2 i \eta ^2}{\sqrt {\pi }} \sum _{m = 0 }^\infty \frac {\varGamma (m+1)}{\varGamma (m+1/2)} \, \alpha _m. \end{align}
Substituting the above relation in (2.12) and then retaining only the first term in the expansion, we arrive at the following equations for
$({\boldsymbol{u}}_1^{(1)} , p_1^{(1)})$
, which are the leading-order contributions to
$({\boldsymbol{u}}_1 , p_1)$
when
$ \varepsilon \ll 1$
:
\begin{align} i \eta ^2 {\boldsymbol{u}}_1^{(1)} = { \boldsymbol{\nabla }} \boldsymbol{\cdot }{ \boldsymbol \sigma } _1^{(1)} = {\nabla} ^2 {\boldsymbol{u}}_1^{(1)} - { \boldsymbol{\nabla }} p_1^{(1)} \quad \text{and} \quad { \boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol{u}}_1^{(1)} = 0 ,\end{align}
with
\begin{align} {\boldsymbol{u}}_1^{(1)} &= {\boldsymbol{0}} \quad \text{for} \quad z = 0 \quad \text{and} \quad \varepsilon \leqslant \rho \leqslant 1, \notag \\ {\boldsymbol{u}}_1^{(1)} &\to {\boldsymbol{0}} \quad \text{as} \quad \sqrt {\rho ^2 + z^2} \to \infty , \notag \\ p_1^{(1)} &= - \hat {p}_0\big |_{\rho = 0} \quad \text{for} \quad z = 0 \quad \text{and} \quad 0 \leqslant \rho \lt \varepsilon , \notag \\ p_1^{(1)} &= 0 \quad \text{for} \quad z = 0 \quad \text{and} \quad \rho \gt 1. \end{align}
The problem posed by (2.31) and (2.32) essentially corresponds to the flow through a pore in an infinite plate, due to the pressure differential
$2 \, \hat {p}_0 \big |_{r=0}$
. This interpretation is substantiated by rescaling the lengths with
$ \varepsilon$
.
Applying the same steps employed to solve the first sub-problem, we obtain
\begin{align} {\boldsymbol{u}}_1^{(1)} \boldsymbol{\cdot }{ \boldsymbol{e }}_\rho & = \int _0^\infty k B(k) \big(e^{-k z} - e^{-\xi z}\big) J_1(k \rho ) \; \textrm { d}k, \end{align}
\begin{align} {\boldsymbol{u}}_1^{(1)} \boldsymbol{\cdot }{ \boldsymbol{e }}_z & = \int _0^\infty k B(k) \left (e^{-k z} - \frac {k}{\xi } \, e^{-\xi z}\right ) J_0(k \rho ) \; \textrm { d}k, \end{align}
\begin{align} p_1^{(1)} & = i \eta ^2 \int _0^\infty B(k) \, e^{-k z} \, J_0(k \rho ) \; \textrm { d}k, \end{align}
where the unknown function
$B(k)$
satisfies the dual integral equations
\begin{align} \int _0^\infty B(k) \, J_0(k \rho ) \; \textrm { d}k & = i p^\star \quad \text{for} \quad 0 \leqslant \rho \leqslant \varepsilon , \end{align}
\begin{align} \int _0^\infty k B(k) \left (1 - \frac {k}{\xi }\right ) J_0(k \rho ) \; \textrm { d}k & = 0 \quad \text{for} \quad \rho \gt \varepsilon , \end{align}
where
$p^\star = \hat {p}_0 \big |_{\rho = 0} / \eta ^2$
. Following Tranter’s method (Tranter Reference Tranter1966) once again, we represent
$B(k)$
as a series of Bessel functions
\begin{align} B(k) = i p^\star \sqrt {\frac {2 k \varepsilon }{\pi }} \left [ k \left (1 - \frac {k}{\xi }\right ) \right ]^{-1} \sum _{m = 0}^\infty \beta _m \, J_{2m+1/2} (k \varepsilon ). \end{align}
This formulation inherently satisfies the velocity boundary condition on the disk. Upon substituting for
$B(k)$
from (2.38) into (2.36), and then multiplying both sides by
\begin{align} \frac { { \tilde {\rho } } }{\sqrt {1 - { \tilde {\rho } } ^2}} \, {_2F_1}( - n , n + 1 / 2 ; 1 ; { \tilde {\rho } } ^2 ) ,\end{align}
and integrating both sides with respect to
$ { \tilde {\rho } } = \rho / \varepsilon$
from
$0$
to
$1$
(see also Erdogan & Bahar Reference Erdogan and Bahar1964), we reach
\begin{align} \sum _{m = 0}^\infty \beta _m \int _0^\infty \left [ { \tilde {k} } \left (1 - \frac { { \tilde {k} } }{ { \tilde \xi } }\right ) \right ]^{-1} J_{2m+1/2}( { \tilde {k} } ) \, J_{2n+1/2}( { \tilde {k} } ) \; \textrm { d} { \tilde {k} } = \delta _{0n} \quad \text{for} \quad n = 0, 1, {\cdots} , \end{align}
where
$ { \tilde {k} } = k \varepsilon$
,
$ { \tilde \xi } ^2 = { \tilde {k} } ^2 + i { \tilde \eta } ^2$
and
$ { \tilde \eta } ^2 = \eta ^2 \varepsilon ^2$
. Equation (2.40) constitutes a system of linear algebraic questions for
$\beta _m$
that can be numerically solved in a manner similar to (2.26). Given (2.16), (2.34) and (2.38), it is straightforward to show that the leading-order contribution to
$F_1$
is given by
\begin{align} F_1^{(1)} &= - 4 \pi \eta ^2 p^\star \, \int _0^ \varepsilon \big({\boldsymbol{u}}_1^{(0)} \big |_{z = 0} \boldsymbol{\cdot }{ \boldsymbol{e }}_z\big) \rho \; \textrm { d}\rho \notag \\ &= 4 i \eta ^2 \sqrt {2 \varepsilon \pi } {p^\star }^2 \sum _{m = 0}^\infty \beta _m \int _0^ \varepsilon \rho \int _0^\infty \sqrt {k} \, J_{2m+1/2}(k \varepsilon ) \, J_0 (k \rho ) \; \textrm { d}k \textrm { d}\rho \notag \\ &= - 8 i \eta ^2 {p^\star }^2 \beta _0 \, \varepsilon . \end{align}
2.4. Limit of low oscillatory Reynolds number
Williams (Reference Williams1966) elegantly demonstrated that the force exerted on an arbitrarily shaped body undergoing slow oscillations can be expressed asymptotically in the limit of small
$\eta$
as
\begin{align} F = F_s \left (1 - \frac {F_s}{6 \pi } \, e^{i \pi /4}\, \eta \right ) + \mathcal{O}(\eta ^2), \end{align}
where
$F_s$
represents the Stokes (
$\eta = 0$
) drag on the body. Below, we derive
$F_s$
for our annular disk, continuing to assume that
$ \varepsilon \ll 1$
. By fixing
$\eta = 0$
, the equations for the pressure distribution and force on the solid disk in the first sub-problem, and the velocity distribution in the
$z$
direction and the force in the second problem, simplify to the following (refer to Happel & Brenner Reference Happel and Brenner1983; Davis Reference Davis1991; Tanzosh & Stone Reference Tanzosh and Stone1996
a, for more detail):
\begin{align} p_{0,s} & = \frac {4}{\pi } \, \int _0^\infty \, \sin k \, e^{-kz} J_0(k \rho ) \; {\rm d}k, \end{align}
\begin{align} F_{0,s} & = - 4 \pi \int _0^1 p_0 \, \rho \; \textrm { d}\rho = -16, \end{align}
\begin{align} {\boldsymbol{u}}_{1,s}^{(1)} \boldsymbol{\cdot }{ \boldsymbol{e }}_z & = \frac {4}{\pi ^2} \, \int _0^\infty \, \frac {d}{{\rm d}k} \left [ \frac {\sin (k \varepsilon )}{k} \right ] \left ( 1 + k z\right ) \, e^{-kz} J_0(kr) \; {\rm d}k, \end{align}
\begin{align} F_{1,s}^{(1)} & = - 4 \pi \, \hat {p}_{0,s} \big |_{\rho = 0} \int _0^ \varepsilon \big( {\boldsymbol{u}}_{1,s}^{(0)}\big |_{z = 0} \boldsymbol{\cdot }{ \boldsymbol{e }}_z \big) \rho \; \textrm { d}\rho = \frac { 64 \varepsilon ^3 } { 3 \pi ^2 }, \end{align}
where
$\hat {p}_{0,s} \big |_{\rho = 0} = 4/\pi$
. In this setting, closed-form expressions for the pressure and velocity fields are also obtainable in oblate spheroidal coordinates (for additional details, see again Happel & Brenner Reference Happel and Brenner1983; Davis Reference Davis1991; Tanzosh & Stone Reference Tanzosh and Stone1996a
). Given (2.44) and (2.46), the Stokes drag on the annular disk in the limit of small pore radius can be written asymptotically as
\begin{align} F_s = -16 \left ( 1 - \frac { 4 \varepsilon ^3 } { 3 \pi ^2 } \right ) + {\cdots} , \end{align}
which is precisely what is arduously derived by Davis (Reference Davis1991b ).
2.5. Limit of high oscillatory Reynolds number
When the oscillatory Reynolds number becomes very large (
$\eta \to \infty$
), we consider two sub-limits based on the magnitude of
$ { \tilde \eta } = \eta \varepsilon$
as discussed below.
2.5.1. Sub-limit of
$ { \tilde \eta } \gg 1$
In cases where both
$\eta$
and
$ { \tilde \eta }$
are much greater than one, it can be shown that the original problem, to the leading order, reduces to the following mixed boundary-value problem for the velocity potential
$\phi$
:
\begin{align} & {\nabla} ^2 \phi = 0, \end{align}
with
\begin{align} { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol{\nabla }} \phi &= 1 \quad \text{for} \quad z = 0 \quad \text{and} \quad \varepsilon \leqslant \rho \leqslant 1, \notag \\ \phi &= 0 \quad \text{for} \quad z = 0 \quad \text{when} \quad 0 \leqslant \rho \lt \varepsilon \quad \text{or} \quad \rho \gt 1, \notag \\ \phi &\to 0 \quad \text{as} \quad \sqrt {\rho ^2 + z^2} \to \infty , \end{align}
where the velocity and pressure fields in this limit are calculated from
${\boldsymbol{u}}_\infty = { \boldsymbol{\nabla }} \phi$
and
$p_\infty = - i \eta ^2 \phi$
, respectively. The potential fields and forces associated with the sub-problems of the above problem are given by (see, e.g. Lamb Reference Lamb1945; Sneddon Reference Sneddon1966; Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007)
\begin{align} \phi _0 &= \frac {2}{\pi } \int _0^\infty \, \frac {d}{{\rm d}k} \left [ \frac {\sin k}{k} \right ] e^{-kz} J_0(kr) \; \textrm { d}k, \end{align}
\begin{align} F_{0,\infty } &= 4 \pi i \eta ^2 \int _0^1 \phi _0 \, \rho \; \textrm { d}\rho = - \frac {8}{3} i \eta ^2, \end{align}
\begin{align} \phi _1^{(1)} &= - \frac {2}{\pi } \hat {\phi }_0 \big |_{r=0} \int _0^\infty \, \frac {\sin (k \varepsilon )}{k} \, e^{-kz} J_0(kr) \; \textrm { d}k, \end{align}
\begin{align} F_{1,\infty }^{(1)} &= 4 \pi i \eta ^2 \, \phi _0 \big |_{r=0} \int _0^ \varepsilon { \boldsymbol n } \boldsymbol{\cdot }{ \boldsymbol{\nabla }} \phi _1^{^{(0)}} \rho \; \textrm { d}\rho = \frac { 32 \varepsilon }{\pi ^2} i \eta ^2, \end{align}
where
\begin{align} \hat {\phi }_0 \big |_{r=0} & = - 2 / \pi, \\[-6pt] \nonumber \end{align}
\begin{align} {\boldsymbol{u}}_{0,\infty } & = { \boldsymbol{\nabla }} \phi _0, \\[-6pt] \nonumber \end{align}
\begin{align} {\boldsymbol{u}}^{(1)}_{1,\infty } & = { \boldsymbol{\nabla }} \phi ^{(1)}_1, \\[-6pt] \nonumber \end{align}
\begin{align} p_{0,\infty } & = - i \eta ^2 \phi _0 , \\[-6pt] \nonumber \end{align}
\begin{align} p^{(1)}_{1,\infty } & = - i \eta ^2 \phi ^{(1)}_1 . \\[6pt] \nonumber \end{align}
Hence, in the potential flow limit, the hydrodynamic force on an annular disk asymptotically takes the form
\begin{align} F_\infty = - \frac {8}{3} i \eta ^2 \left ( 1 - \frac {12 \varepsilon }{\pi ^2} \right ) + {\cdots} . \end{align}
It is worth noting that the velocity potentials in (2.50) and (2.52) can be alternatively expressed in closed form in oblate spheroidal coordinates (see, e.g. Lamb Reference Lamb1945).
2.5.2. Sub-limit of
$ { \tilde \eta } \ll 1$
In the regime in which
$ \varepsilon \ll 1 / \eta$
while
$\eta \gg 1$
, the leading-order solution to the first problem remains the same as the one given by (2.50) and (2.51), i.e.
\begin{align} F_0 = F_{0,\infty } + \mathcal{O}( \eta \ln \eta ). \end{align}
The origin of the
$\mathcal{O}( \eta \ln \eta )$
correction in the above equation will be discussed in § 3 (see also table 1).
Table 1. Asymptotic expressions describing the variations of
$a_0$
,
$b_0$
and
$p^\star$
as functions of
$\eta$
, along with the dependence of
$\beta _0$
on
$ { \tilde \eta }$
. The coefficients
$\mathcal{A}$
,
$\mathcal{B}$
and
$\mathcal{C}$
are numerically computed and their values are reported in table 2.
Figure 2. Plots of dimensionless (a) added mass and (b) viscous damping coefficients of a solid disk (denoted by
$a_0$
and
$b_0$
, respectively), along with the (c) real and (d) imaginary components of the normalised pressure at the centre of the solid disk (defined as
$p^\star = \hat {p}_0 |_{\rho = 0} / \eta ^2$
), all versus the square root of the oscillatory Reynolds number
$\eta$
. Red circular and blue triangular symbols correspond to the asymptotic formulas for the low- and high-frequency limits, respectively, as listed in table 1.
Since
$ { \tilde \eta } \ll 1$
, the solution to the second problem follows from (2.45) and (2.46) for Stokes flow through a circular hole in an infinite wall. Given that the driving pressure here is
$ \hat {p}_{0,\infty } |_{\rho = 0}$
rather than
$\hat {p}_{0,s} |_{\rho = 0}$
, we then obtain
\begin{align} F_{1,\infty ,s} = \left ( \frac { \hat {p}_{0,\infty } \big |_{\rho = 0} }{ \hat {p}_{0,s} \big |_{\rho = 0} } \right )^2 F_{1,s}^{(1)} + \mathcal{O}\big( \varepsilon ^3 \eta ^3 \ln \eta \big) = - \frac { 16 \varepsilon ^3 \eta ^4}{ 3 \pi ^2 } + \mathcal{O}\big( \varepsilon ^3 \eta ^3 \ln \eta \big). \end{align}
3. Results and discussion
3.1. Baseline: a solid disk
For the baseline case of a solid disk, we solve a linear system of equations based on (2.26) to obtain
$\alpha _m$
for
$m = 0, 1, \ldots , 10$
across a range of
$\eta$
from
$10^{-2}$
to
$10^{3}$
. Accordingly, we retain the first ten terms in the series and evaluate the improper integrals using a global adaptive quadrature with a strict absolute tolerance to maintain accuracy. Using the computed values of
$\alpha _m$
, we employ (2.5), (2.28) and (2.30) to determine the added mass and damping coefficients, denoted by
$a_0$
and
$b_0$
, respectively, and to compute the normalised pressure at the centre of the disk, defined as
$p^\star = \hat {p}_0 \big |_{\rho = 0} / \eta ^2$
. For the pressure, repeated Shanks transformation is used to accelerate the convergence of the series.
The results of these calculations are presented as solid black lines in figure 2, where red circular and blue triangular symbols represent, respectively, the asymptotic expressions for
$a_0$
,
$b_0$
and
$p^\star$
in the limits of low and high frequencies, as outlined in table 1. A key observation is that the asymptotic curves closely overlap in an intermediate range of
$\eta$
, demonstrating their combined accuracy across the entire frequency spectrum. The low-frequency asymptotic formulas result directly from the derivation presented in § 2.4. These expressions are found to be in full agreement with the results reported by Zhang & Stone (Reference Zhang and Stone1998). However, for the high-frequency asymptotics, only the
$\mathcal{O}(1)$
terms, resulting from the potential flow theory calculations of § 2.5.1, match with previous studies. Additionally, the coefficients of the
$\mathcal{O}(\eta ^{-1} \ln \eta )$
and
$\mathcal{O}(\eta ^{-1})$
terms are computed numerically and are listed in table 2 for reference.
Table 2. Numerically computed values of the coefficients appearing in table 1.
Of particular significance is the presence of
$\mathcal{O}(\eta ^{-1} \ln \eta )$
terms in the high-frequency expansions. These terms were previously overlooked in the literature, yet they play a critical role in accurately capturing the behaviour of the oscillating disk in the high-frequency regime. While prior studies have hinted at a singular limit when transitioning from an oscillating oblate spheroid to a zero-thickness disk (Loewenberg Reference Loewenberg1993; Zhang & Stone Reference Zhang and Stone1998), our analysis rigorously confirms the necessity of incorporating logarithmic terms alongside the inverse frequency terms. To further support the inclusion of
$\mathcal{O}(\eta ^{-1} \ln \eta )$
terms, below, we present an asymptotic analysis of the integrals involved in (2.26).
3.2. Evaluating integrals appearing in (
2.26
) in the limit of
$\eta \to \infty$
For clarity, we first establish our argument for a representative case before extending it to more general scenarios. Specifically, we consider the case where
$ m = 0$
and
$ n = 1$
, in which the corresponding integral can be split into three parts as follows:
\begin{align} \int _0^\infty k \left ( 1 - \frac {k}{\xi } \right ) J_{1/2}(k) \, J_{5/2}(k) \textrm { d}k &= \int _0^1 k \left ( 1 - \frac {k}{\xi } \right ) J_{1/2}(k) \, J_{5/2}(k) \textrm { d}k \notag \\ &\quad + \int _1^\eta k \left ( 1 - \frac {k}{\xi } \right ) J_{1/2}(k) \, J_{5/2}(k) \textrm { d}k \notag \\ &\quad + \int _\eta ^\infty k \left ( 1 - \frac {k}{\xi } \right ) J_{1/2}(k) \, J_{5/2}(k) \textrm{ d}k, \end{align}
where the relevant Bessel functions and their product are given by
\begin{align} J_{1/2}(k) & = \sqrt {\frac {2}{\pi k}} \sin k, \end{align}
\begin{align} J_{5/2}(k) & = \sqrt {\frac {2}{\pi k}} \left [ \left ( \frac {3}{k^2} - 1 \right ) \sin k - \frac {3 \cos k}{k} \right ]\!, \end{align}
\begin{align} J_{1/2}(k) J_{5/2}(k) & = \frac {1}{\pi k} \left [ \left ( \frac {3}{k^2} - 1 \right ) + \left ( 1 - \frac {3}{k^2} \right ) \cos 2k - \frac {3}{k} \sin 2k \right ]\!. \\[12pt] \nonumber \end{align}
The term in the parenthesis in (3.1) can be expanded as a Taylor series in the region
$k \leqslant \eta$
as
\begin{align} 1 - \frac {k}{\xi } = 1 - \frac {\left ( k / \eta \right )}{\sqrt {2}} \left [ 1 + \frac {\left (k / \eta \right )^2}{2} \right ] + {\cdots} + i \frac {\left ( k / \eta \right )}{\sqrt {2}} \left [ 1 - \frac {\left ( k / \eta \right )^2}{2} \right ] + {\cdots} . \end{align}
Similarly, in the region
$k \geqslant \eta$
, it takes the form
\begin{align} 1 - \frac {k}{\xi } = \frac {3 \left ( \eta / k \right )^4}{8} \left [ 1 - \frac {35 \left ( \eta / k \right )^4 }{48} \right ] + {\cdots} + i \frac { \left ( \eta / k \right )^2 }{2} \left [ 1 - \frac { 5 \left ( \eta / k \right )^4 }{8} \right ] + {\cdots} . \end{align}
Substituting (3.4) and these expansions in the integrals on the right-hand side of (3.1), we find that a term proportional to
$ (k \eta )^{-1}$
emerges in the middle integral, covering the range
$1 \leqslant k \leqslant \eta$
. The integral of this term contributes a term proportional to
$\ln \eta$
, thereby justifying the appearance of
$\mathcal{O}(\eta ^{-1}\ln \eta )$
in the asymptotic expressions reported in table 1 for the limit
$\eta \to \infty$
.
This result extends naturally to other integrals, except for the special case of
$ m = n = 0$
. In particular, the systematic emergence of an
$\mathcal{O} [ (k \eta )^{-1} ]$
term can be demonstrated by repeating the steps outlined above and utilising the following general formula for half-integer order Bessel functions in (2.26):
\begin{align} J_{2n + 1/2}(k) = \sqrt {\frac {2}{\pi k}} \left [ \mathcal{R}_{2n,1/2}(k) \sin k - \mathcal{R}_{2n-1,3/2}(k) \cos k \right ]\!. \end{align}
Here, the Lommel polynomial is defined as
\begin{align} \mathcal{R}_{\ell ,\nu }(k) =\sum _{l=0}^{[\ell /2]} \frac {\left (-1\right )^l \left (\ell - l\right )! \varGamma (\ell - l + \nu )}{l! \left (\ell - 2 l\right )! \varGamma (l + \nu )} \left (\frac {2}{k}\right )^{\ell - 2 l}, \end{align}
where the square brackets denote the floor function, which rounds down to the nearest integer.
Figure 3. Plots of the (a) real and (b) imaginary components of
$\beta$
versus
$ { \tilde \eta } = \eta \varepsilon$
. Red circular and blue triangular symbols correspond to the asymptotic formulas for the limits of
$ { \tilde \eta } \to 0$
and
$ { \tilde \eta } \to \infty$
, respectively, as listed in table 1.
3.3. Added mass and damping coefficients of an annular disk
As the final step, we solve an additional linear system of equations based on (2.40) to determine
$\beta _m$
over the range
$10^{-2} \leqslant { \tilde \eta } \leqslant 10^3$
. For small values of
$ { \tilde \eta }$
, we retain the first ten terms in the series, progressively increasing this count to a maximum of forty terms for the largest
$ { \tilde \eta }$
considered. The results for
$\beta _0$
are presented in figure 3, where symbols denote the asymptotic expressions reported in table 1. As previously indicated, the coefficients in this table that are not directly derived from the analyses of §§ 2.4 and 2.5 are computed numerically (see table 2). Once again, we observe the remarkable accuracy of these simple formulas in capturing the variations of both the real and imaginary components of
$\beta _0$
as functions of
$ { \tilde \eta }$
. We note that the inclusion of
$\mathcal{O}( { \tilde \eta } ^{-1} \ln { \tilde \eta } )$
terms in the asymptotic limit of
$ { \tilde \eta } \to \infty$
follows from the same reasoning used in § 3.2. Naturally, for the integrals in (2.40), different Taylor series expansions must be employed for the term that is multiplied by the product of Bessel functions.
Figure 4. Plots of
$a / a_0$
(left column) and
$b / b_0$
(right column) versus
$ \varepsilon$
. The top, middle and bottom rows correspond to
$\eta = 10^{-2}, 10^{-1} \, \text{and} \, 10^0$
, respectively. The solid black lines represent our theoretical calculations, whereas the dashed green curves show the results from finite-element numerical simulations. The inset plots illustrate the deviations of the normalised coefficients from unity. The red circles and purple squares denote the asymptotic formulas from table 3 in the limiting cases of
$\eta \to 0$
, and
$\eta \to \infty$
and
$ { \tilde \eta } \to 0$
, respectively.
Figure 5. Plots of
$a / a_0$
(left column) and
$b / b_0$
(right column) versus
$ \varepsilon$
. The top, middle and bottom rows correspond to
$\eta = 10^1, 10^2 \, \text{and} \, 10^3$
, respectively. The solid black lines represent our theoretical calculations, whereas the dashed green curves show the results from finite-element numerical simulations. The inset plots illustrate the deviations of the normalised coefficients from unity. The purple squares and blue triangles denote the asymptotic formulas from table 3 in the limiting cases of
$\eta \to \infty$
and
$ { \tilde \eta } \to 0$
, and
$\eta \to \infty$
and
$ { \tilde \eta } \to \infty$
, respectively.
With
$\beta _0$
determined, alongside the previously computed
$p^\star$
and
$\alpha _0$
from the baseline problem, we now analyse the added mass and damping coefficients of the annular disk (see (2.5), (2.17), (2.28) and (2.41)). The solid black lines in figures 4 and 5 illustrate the variations of the normalised added mass and damping coefficients (
$a / a_0$
and
$b / b_0$
, respectively) as functions of the pore size (
$\varepsilon$
) across a broad range of oscillatory Reynolds numbers (
$\eta$
). From top to bottom, the rows in figure 4 correspond to
$\eta = 10^{-2}, 10^{-1} \, \text{and} \, 10^0$
, while those in figure 5 represent
$\eta = 10^1, 10^2 \, \text{and} \, 10^3$
. The red circles, purple squares, and blue triangles denote, respectively, the asymptotic formulas listed in table 3 (derived from table 1) in the limits of
$\eta \to 0$
,
$\eta \to \infty$
and
$ { \tilde \eta } \to 0$
and
$\eta \to \infty$
and
$ { \tilde \eta } \to \infty$
. Insets highlight deviations of the normalised coefficients from unity, emphasising how porosity-induced corrections to the added mass and damping coefficients scale with
$ \varepsilon$
in different regimes (see table 3).
The green dashed curves in figures 4 and 5 represent the results of a numerical solution to the original problem (i.e. (2.3) and (2.4)) in the frequency domain, obtained using a second-order finite-element approach as implemented in COMSOL Multiphysics (Zimmerman Reference Zimmerman2006; Pepper & Heinrich Reference Pepper and Heinrich2017). These results serve as a benchmark for assessing the accuracy of our small-pore assumption, which led us to retain only the leading term in the Taylor series expansion of
$\hat {p}_0$
in (2.29). To model the unbounded domain in the simulations, the outer boundary was approximated as a large hemisphere centred at the disk’s midpoint. The radius of this hemisphere was varied between
$100 R$
and
$2000 R$
depending on the value of
$\eta$
. The two-dimensional axisymmetric computational domain was discretised using triangular elements with adaptive mesh refinement to ensure accuracy and efficiency. The validity of this numerical approach was confirmed through comparison with the calculations in § 3.1 for a solid disk, as well as through grid and domain-size independence studies, which verified that the computed added mass and damping coefficients changed by less than 1 % upon further refinement.
A preliminary examination of both figures reveals that our asymptotic formulas accurately capture the behaviour of normalised coefficients over a broad range of
$ \varepsilon$
, demonstrating remarkable precision except for the high-
$\eta$
asymptote’s subpar performance in figure 5(b). Furthermore, our theoretical predictions closely align with finite-element simulations for
$ \varepsilon$
up to 1/2, and in many cases, even beyond.
A closer inspection indicates that both normalised coefficients trend downward with increasing pore size for
$\eta \lesssim 1$
. As highlighted in the insets, their departures from unity follow a
$ \varepsilon ^3$
scaling (see figure 4). For
$\eta \lesssim 10^{-1}$
, the asymptotic expressions for
$\eta \to 0$
provide excellent accuracy but become less reliable as
$\eta$
grows. At larger
$\eta$
, the behaviour of
$a / a_0$
is well described by a combination of asymptotic formulas for
$\eta \to 0$
and
$ { \tilde \eta } \to 0$
at small pore sizes, and
$\eta \to \infty$
and
$ { \tilde \eta } \to \infty$
at bigger pore sizes. In this regime, the initial
$ \varepsilon ^3$
scaling for the decrease in
$a / a_0$
transitions to a linear trend as
$ \varepsilon$
increases (see the left column of figure 5).
Table 3. Asymptotic expressions for the normalised added mass and damping coefficients (
$a / a_0$
and
$b / b_0$
, respectively) in the limiting cases of
$\eta \to 0$
,
$\eta \to \infty$
and
$ { \tilde \eta } \to 0$
and
$\eta \to \infty$
and
$ { \tilde \eta } \to \infty$
. These formulas are derived based on those reported in table 1. The value of the coefficients
$\mathcal{A}$
and
$\mathcal{B}$
are given in table 2.
A key finding of this study is that, for
$\eta \gtrsim 1$
, while the added mass coefficient continues to decrease with increasing the pore size, the damping coefficient initially rises with
$ \varepsilon$
, reaches a peak and then declines as the pore size increases further. This initial increase is predicted by our asymptotic analysis, which anticipates a rise proportional to
$ \varepsilon ^3$
(see § 2.5.2 and table 3). Notably, despite the oscillatory Reynolds number based on the disk’s outer radius (
$\eta$
) being high, the Reynolds number based on the pore radius,
$ { \tilde \eta } = \eta \varepsilon$
, remains well below unity. The contrast in regimes between the first and second problems gives rise to the initial enhancement of the damping coefficient due to porosity.
The peak of
$b / b_0$
and its corresponding pore size (
$ \varepsilon _{max}$
) are nicely predicted by our asymptotic formula in the limit of
$\eta \to \infty$
and
$ { \tilde \eta } \to \infty$
for
$\eta \gtrsim 10^2$
(see the last two rows of the right column of figure 5). By setting the derivative of the asymptotic formula for
$b / b_0$
with respect to
$ \varepsilon$
to zero, we obtain
\begin{align} \varepsilon _{max} = \frac {{\mathcal{A}}^\beta _i}{\pi \left ( {\mathcal{A}}^p_r \ln \eta + {\mathcal{B}}^p_r \right )}. \end{align}
Figure 6 illustrates the variations of
$ \varepsilon _{max}$
and the corresponding
$b / b_0$
and
$a / a_0$
as functions of
$\eta$
, with asymptotic predictions denoted by blue triangles. As evident from (3.9),
$ \varepsilon _{max}$
decreases logarithmically with increasing
$\eta$
. In the asymptotic limit
$\eta \to \infty$
,
$b / b_0$
and
$a / a_0$
approach very slowly to 1.62 and 1, respectively. Although achieving about 60 % damping improvement may be challenging due to the requirement of very high
$\eta$
, reaching a 30 % enhancement is feasible in many engineering applications.
Figure 6. Plots of (a) the pore size corresponding to the maximum normalised damping coefficient (
$ \varepsilon _{max}$
), (b) the peak value of
$b / b_0$
, and (c)
$a /a_0$
at
$ \varepsilon _{max}$
, all as functions of
$\eta$
. The solid black lines represent our theoretical calculations, whereas the dashed green curves show the results from finite-element numerical simulations. Also, the blue triangles represent the predictions from (3.9) and its substitution into the asymptotic formulas for
$b / b_0$
and
$a / a_0$
in the limit
$\eta \to \infty$
and
$ { \tilde \eta } \to \infty$
, as reported in table 3.
Figure 7. (a)–( f) plots of the normalised total force
$|F| / |F_0|$
versus
$ \varepsilon$
for
$\eta$
, with the insets showing the corresponding variations of the phase angle
$\phi$
as a function of
$ \varepsilon$
. The panels correspond to
$\eta = 10^{-2}$
(a),
$10^{-1}$
(b),
$1$
(c),
$10$
(d),
$10^{2}$
(e) and
$10^{3}$
( f). Also, the green dashed lines represent the results of finite-element numerical simulations.
To conclude this section, we analyse the supplementary plots in figure 7, which depict the variations of the following two quantities with
$ \varepsilon$
for the considered values of
$\eta$
:
\begin{align}& \frac {|F|}{|F_0|} = \frac {\sqrt {a^2 + b^2}}{\sqrt {a_0^2 + b_0^2}}, \end{align}
\begin{align}&\,\,\, \phi = \tan ^{-1}\left (\frac {a}{b}\right )\!. \\[12pt] \nonumber \end{align}
As before, dashed green lines indicate simulation results. A notable observation is that, despite the increase in the damping coefficient for sufficiently high
$\eta$
, the normalised total force always trends downward due to porosity, regardless of the oscillatory Reynolds number (see the main panels of figure 7). Furthermore, the phase angle
$\phi$
increases from
$0$
to
$\pi /2$
as
$\eta$
approaches
$\infty$
(see the insets of figure 7). To reconcile the contrasting behaviours of the normalised total force with that of the normalised damping coefficient, we consider the relation
\begin{align} b = \frac {|F|}{\eta ^2} \cos \phi . \end{align}
At sufficiently high
$\eta$
, the reduction in the total force due to the porosity of the disk is more than offset by the decrease in the phase angle. Since
$\phi$
approaches
$\pi /2$
in this regime, even a small decrease in
$\phi$
leads to a relatively large change in the magnitude of
$\cos \phi$
, thereby influencing the damping behaviour significantly. This counterintuitive increase in the damping coefficient is thus purely a result of the phase shift effect and not due to additional flow separation or vortex shedding.
4. Summary
We developed a semi-analytical framework to investigate the hydrodynamics of an annular disk undergoing small-amplitude, broadside sinusoidal oscillations in an unbounded fluid. Utilising the superposition principle, the reciprocal theorem and asymptotic analysis, we systematically explored the influence of porosity and oscillation frequency on the disk’s added mass and damping coefficients. By decomposing the problem into two sub-problems – one describing a solid disk in heave motion and the other modelling flow through a circular orifice in an infinite wall – we derived asymptotic expressions that accurately capture the variations of the hydrodynamic coefficients across a wide range of oscillatory Reynolds numbers and pore sizes. Numerical simulations confirmed the validity of our asymptotic formulations, demonstrating their accuracy for inner-to-outer radius ratios up to at least
$1/2$
.
A key finding of our analysis is the identification of a previously overlooked logarithmic term in the high-frequency asymptotic expansions, which plays a crucial role in accurately characterising the behaviour of the added mass and damping coefficients. More importantly, we uncovered a non-trivial dependence of the damping coefficient on pore size. While the added mass coefficient decreases monotonically with increasing porosity, the damping coefficient initially rises, reaching a maximum before eventually declining. Remarkably, in the extreme inertial limit, the peak damping coefficient surpasses that of a solid disk by up to 62 %, potentially challenging conventional assumptions about how porosity affects oscillatory fluid–structure interactions. The transition from a viscosity-dominated regime to an inertial regime reveals intricate interplay between permeability effects and oscillatory flow behaviour, highlighting the complex role of porosity in the fluid–structure dynamics.
These insights have considerable implications for engineering applications such as designing offshore platforms and wave energy converters, where porosity could be strategically optimised to enhance stability and efficiency. Our findings suggest that tuning porosity offers a promising means of controlling hydrodynamic forces, particularly in high-frequency oscillatory environments. Future work may extend this framework to explore nonlinear effects or larger pore geometries; while this would require retaining higher-order terms and modifying parts of the derivation, leading-order nonlinear corrections may be accessible via the reciprocal theorem, and the large-pore limit could be studied by modelling the annular disk as a ring, following Davis (Reference Davis1991b ).
Funding
Financial support from the Water Power Technologies Office of the U.S. Department of Energy through a Seedling grant is gratefully acknowledged.
Declaration of interests
The authors declare no conflict of interest.
Open access
