2.1. Porous boundary condition

Along a porous aerofoil, the perturbation flow velocity on the aerofoil surface, $w$, is related to the local seepage flow rate directed along the unit normal to the wing surface, $w_{s}$, by

(2.1)\begin{equation} w(x,t)=w_{s}(x,t)+\frac{\partial y_{a}}{\partial x}(x,t) + \frac{\partial y_{a}}{\partial t}(x,t), \end{equation}

where the function $y_{a}(x,t) = y_{c}(x,t) - \alpha (t)$ defines the mean surface of the aerofoil relative to the angle of attack. This linearised boundary condition assumes that the deformation of the aerofoil mean line is small: $|y_{a}|\ll 1$. Additionally, the aerofoil thickness is assumed to be small. Aerofoils with small but non-zero thickness can typically be handled with thin aerofoil theory (see § 4.2 of Van Dyke (Reference Van Dyke1964) or Baddoo & Ayton (Reference Baddoo and Ayton2018)) but in this problem aerofoil thickness does not contribute to the pressure jump across the chord. Thus, aerofoil thickness effects decouple from the hydrodynamic effects of porosity. We therefore do not consider the effects of thickness further here. The seepage velocity depends on the pressure jump across the aerofoil and on the local porous structure. Porosity imparts non-dimensional hydrodynamic inductance $\rho _{e}$ and resistance $\varPhi$ due to the viscous fluid–solid interactions of the flow within the porous medium. Combining these effects yields the unsteady porosity boundary condition (Morse & Ingard Reference Morse and Ingard1986, p. 254)

(2.2)\begin{equation} 2 \rho_{e}(x) \frac{\partial w_{s}}{\partial t}(x,t) + \varPhi(x) w_{s}(x,t) = - {\rm \Delta} p(x,t), \end{equation}

where ${\rm \Delta} p$ is the local difference of the surface pressure (upper minus lower). The dimensional scalings used in this work identify the so-called effective density $\rho _{e} = m_{e}/V > 1$, where $m_{e} \ge 1$ is the effective mass or structure factor of the porous medium, and $V$ is the fluid-to-solid volumetric fraction of the porous matrix (Bliss Reference Bliss1982; Attenborough Reference Attenborough1983). Similarly, the dimensionless flow resistance $\varPhi$ is the pressure drop per unit length per unit mean flow velocity through the porous aerofoil scaled by $\rho U/(2 l)$. Strictly speaking, $\rho _{e}$ and $\varPhi$ are also functions of frequency, although we do not consider that detail here.

The linearised Bernoulli equation for unsteady flow enables the pressure to be expressed as a function of the velocity potential $\phi$ as

(2.3)\begin{equation} p(x,y,t) = -\left(\frac{\partial \phi}{\partial x} + \frac{\partial \phi}{\partial t}\right).\end{equation}

Applying (2.3) to $y = 0^\pm$, $x>-1$, and taking the difference from the upper and lower sides yields expressions for the pressure jump along the aerofoil and the wake

(2.4a)\begin{gather} {\rm \Delta} p(x,t) = -2 \left(\gamma_{a}(x,t) + \frac{\partial }{\partial t} \int_{-1}^x \gamma_{a}(\xi,t)\,\textrm{d} \xi \right), \quad -1 < x < 1, \end{gather}

(2.4b)\begin{gather}{\rm \Delta} p(x,t) = -2 \left(\gamma_{w}(x,t) + \frac{\partial }{\partial t} \int_{1}^x \gamma_{w}(\xi,t)\,\textrm{d} \xi + \frac{\textrm{d} \varGamma }{\textrm{d} t}(t)\right), \quad 1 < x, \end{gather}

where $\varGamma$ represents the circulation around the aerofoil. On the other hand, application of the Plemelj formulae (Ablowitz & Fokas Reference Ablowitz and Fokas2003) to the Biot–Savart law shows that the fluid normal velocity on the wing, $w$, is related to the vorticity distributions on the aerofoil ($\gamma _{a}$) and in the wake ($\gamma _{w}$) through the following singular integral equation (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, (5-313a))

(2.5)\begin{equation} w(x,t)= \frac{1}{2{\rm \pi}}\int\kern-10pt_{-1}^{1}\frac{\gamma_{a}(\xi,t)}{\xi-x}\,\textrm{d}\xi +\frac{1}{2{\rm \pi}}\int_{1}^{\infty}\frac{\gamma_{w}(\xi,t)}{\xi-x} \,\textrm{d}\xi, \quad -1 < x < 1, \end{equation}

where the bar on the integral denotes the Cauchy principal value. Note that once $\gamma _a$ and $\gamma _w$ have been determined, the full complex velocity field can be computed from the Biot–Savart law. Combining (2.1), (2.2), (2.4a) and (2.5) yields the integral equation

(2.6)\begin{align} & \left[ 2 \rho_{e}(x) \frac{\partial}{\partial t} + \varPhi(x) \right] \left\{ \frac{1}{2{\rm \pi}}\int\kern-10pt_{-1}^{1}\frac{\gamma_{a}(\xi,t)}{\xi-x}\,\textrm{d}\xi + \frac{1}{2{\rm \pi}} \int_{1}^{\infty}\frac{\gamma_{w}(\xi,t)}{\xi-x}\, \textrm{d}\xi - \frac{\partial y_{a}}{\partial t} - \frac{\partial y_{ a}}{\partial x} \right\} \nonumber\\ &\quad = 2 \left(\gamma_{a}(x,t) + \frac{\partial }{\partial t} \int_{-1}^x \gamma_{a}(\xi,t)\,\textrm{d} \xi \right), \end{align}

for $-1 < x < 1$.

We further assume harmonic motions of reduced frequency $k$, so that we may write the vorticity distributions and mean camber line as

(2.7a–c)\begin{equation} {\gamma}_{{a}}(x,t) = \hat{\gamma}_{{a}}(x)\, \textrm{e}^{\textrm{i} k t}, \quad {\gamma}_{{w}}(x,t) = \hat{\gamma}_{{w}}(x)\, \textrm{e}^{\textrm{i} k t} \quad {\rm and} \quad y_{a}(x,t) = \hat{y}_{a}(x)\, \textrm{e}^{\textrm{i} k t},\end{equation}

where the real part is assumed. Since we do not allow any pressure jump across the wake, we require the expression on the right-hand side of (2.4b) to vanish. Solving the associated integral equation yields

(2.8)\begin{equation} \hat{\gamma}_{{w}}(x) = - \textrm{i} k \,\textrm{e}^{\textrm{i} k (1-x)} \int_{-1}^1 \hat{\gamma}_{a}(\xi)\,\textrm{d} \xi = -\textrm{i} k \hat{\varGamma}\,\textrm{e}^{\textrm{i} k (1-x)}, \end{equation}

where Kelvin's circulation theorem has been applied to enforce that the net circulation of the aerofoil and its wake vanishes, and $\hat {\varGamma }$ is the circulation around the aerofoil with the time dependence factored out. We also enforce the Kutta condition, namely that the pressure jump vanishes at the trailing edge

(2.9)\begin{equation} {\rm \Delta} p(1) = 0. \end{equation}

Substitution of (2.7a–c) and (2.8) into (2.6) yields

(2.10)\begin{equation} \mathcal{L}\hat{\gamma}_{{a}} = f_{a}(x) + \hat{\varGamma} f_{w}(x),\end{equation}

for $-1 < x < 1$, where we use the notation $\mathcal {L}$ to represent the singular linear operator

(2.11)\begin{equation} \mathcal{L}f \triangleq - {\textrm{i} k}\psi(x,k) \int_{-1}^x f(\xi)\,\textrm{d} \xi + \frac{1}{\rm \pi}\int\kern-10pt_{-1}^{1}\frac{f(\xi)}{\xi-x}\,\textrm{d} \xi - \psi(x,k) f(x), \end{equation}

and

(2.12)\begin{equation} \psi (x,k) \triangleq \frac{4}{2 \textrm{i} k \rho_{e}(x) + \varPhi(x)}.\end{equation}

The forcing functions $f_{a}$ and $f_{w}$ are defined as

(2.13)\begin{gather} f_{a}(x) \triangleq 2 \left( \frac{\textrm{d} \hat{y}_{a}}{\textrm{d} x}(x) + \textrm{i} k \hat{y}_{a}(x) \right), \end{gather}

(2.14)\begin{gather}f_{w}(x) \triangleq \frac{\textrm{i} k}{\rm \pi} \int_{1}^\infty \frac{\textrm{e}^{\textrm{i} k (1-\xi)}}{\xi - x}\, \textrm{d} \xi. \end{gather}

The subscript notations ‘${a}$’ and ‘${w}$’ are again employed here to symbolise that $f_{a}$ corresponds to contributions from the mean aerofoil profile and its motions, whereas $f_{w}$ corresponds to contributions from the unsteady wake. The problem is now to determine the function $\hat {\gamma }_{a}$, from which the constant $\hat {\varGamma }$ can be found via integration per (2.8).

The operator $\mathcal {L}$ consists of two parts: a Volterra part (the first term in (2.11)), and a singular Fredholm integral part (the second and third terms in (2.11)). Accordingly, we refer to (2.10) as a singular Fredholm–Volterra integral equation (SF–VIE). The literature on these types of integral equations is apparently non-existent: the closest comparisons that could be found by the authors considered only the case where the kernel is weakly singular (Darwish Reference Darwish1999; Abdou Reference Abdou2003) and not the Cauchy principal value considered in the present work. In particular, it is the presence of the Volterra part of $\mathcal {L}$ that precludes the possibility of a solution using the classical singular integral equation methods of Muskhelishvili (Reference Muskhelishvili1946). Consequently, we now seek a numerical solution by expanding $\hat {\gamma }_{{a}}$ into an appropriate series of basis functions.

3.1. Motivation – the steady case

We first consider the case where the field is steady ($k=0$) and the wake vanishes, as illustrated in figure 1(a). The SF–VIE (2.10) for the bound vorticity distribution becomes

(3.1)\begin{equation} \frac{1}{\rm \pi}\int\kern-10pt_{-1}^{1}\frac{\hat{\gamma}_{a}(\xi)}{\xi-x}\,\textrm{d} \xi -\psi(x,0) \hat{\gamma}_{a}(x) = 2\frac{\textrm{d} \hat{y}_{a}}{\textrm{d} x}(x), \quad -1 < x < 1. \end{equation}

In the impermeable case ($\psi \equiv 0$), the typical solution approach is to expand $\hat {\gamma }_{a}$ in terms of weighted Chebyshev polynomials (Rienstra Reference Rienstra1992). However, this approach dictates the behaviour of the vorticity distribution at the endpoints. In particular, the vorticity distribution is usually written as

(3.2)\begin{equation} \hat{\gamma}_{a}(x) = \hat{\gamma}_{0} \sqrt{\frac{1-x}{1+x}} + \sqrt{1-x^2} \sum_{n=1}^N \hat{\gamma}_{n} U_{n-1} (x), \end{equation}

where $U_n$ are the Chebyshev polynomials of the second kind and $\hat {\gamma }_{n}$ are coefficients to be determined. Consequently, $\hat {\gamma }_{a}$ possesses a square-root singularity at the leading edge and a square-root zero at the trailing edge. This series necessarily satisfies the steady Kutta condition at the trailing edge. However, as we will now show, this choice of basis expansion leads to invalid results at the endpoints when the aerofoil is permeable.

By sending $x \rightarrow -1$, we obtain the asymptotic behaviours

(3.3a)\begin{gather} \psi(x,0)\hat{\gamma}_{a}(x) \sim \psi(-1,0) \hat{\gamma}_{0} \sqrt{\frac{2}{1+x}}, \end{gather}

(3.3b)\begin{gather}\frac{1}{\rm \pi}\int\kern-10pt_{-1}^{1}\frac{\hat{\gamma}_{a}(\xi)}{\xi-x}\,\textrm{d} \xi\sim \varPhi^\ast(x), \end{gather}

(3.3c)\begin{gather}2 \frac{\textrm{d} y_{a}}{\textrm{d} x}(x) \sim 2 \frac{\textrm{d} y_{a}}{\textrm{d} x}(-1), \end{gather}

where $\varPhi ^\ast (x)= o((1+x)^{-1/2})$ according to Muskhelishvili (Reference Muskhelishvili1946, (29.8)). Substitution of these limits into (3.1) results in a contradictory equation where the left-hand side scales like $(1+x)^{-1/2}$ whereas the right-hand side tends to a constant as $x \rightarrow -1$. Asymptotic analysis at the trailing edge generates similar contradictions. Consequently, the Chebyshev expansion generates spurious results at both endpoints, and the $\hat {\gamma }_{n}$ coefficients for $n>1$ must account for the contradiction, resulting in a slowly converging series. The modification of the square-root behaviour at the endpoints due to porosity is embedded in the partially porous aerofoil solution by Iosilevskii (Reference Iosilevskii2011, Reference Iosilevskii2013) and is detailed in the generalised porous aerofoil solution by Hajian & Jaworski (Reference Hajian and Jaworski2017).

Suppose we do not explicitly enforce the square-root behaviour of $\hat {\gamma }_{a}$ at the endpoints and instead express $\hat {\gamma }_{a}$ in the form

(3.4)\begin{equation} \hat{\gamma}_{a}(x) =w^{\alpha,-\beta}(x) \hat{\gamma}_{a}^\ast(x), \end{equation}

where $w^{a,b}$ represents the weight function

(3.5)\begin{equation} w^{a,b}(x) \triangleq (1-x)^a (1+x)^b.\end{equation}

The function $\hat {\gamma }_{a}^\ast (x)$ is Hölder continuous on $x \in [-1,1]$ and is finite and non-zero at $x=\pm 1$. The constants $\alpha$ and $\beta$ in (3.4) are unknown and must be found as a part of the solution. Using the new expansion (3.4), the limits (3.3a) and (3.3b) instead become

(3.6a)\begin{gather} \psi(x,0)\hat{\gamma}_{a}(x) \sim \psi(-1,0) \hat{\gamma}_{a}^\ast(-1) \frac{2^\alpha}{(1+x)^{\beta}}, \end{gather}

(3.6b)\begin{gather}\frac{1}{\rm \pi}\int\kern-10pt_{-1}^{1}\frac{\hat{\gamma}_{a}(\xi)}{\xi-x}\,\textrm{d} \xi\sim \hat{\gamma}_{a}^\ast(-1) \frac{2^\alpha \cot({\rm \pi} \beta)}{(1+x)^{\beta}}+ \varPhi^\ast(x), \end{gather}

where now $\varPhi ^*(x) = o((1+x)^{-\beta })$. Accordingly, by matching the singularities in the above two terms through (3.1), we obtain the following expression for $\beta$:

(3.7)\begin{equation} \beta = \frac{1}{\rm \pi} \cot^{-1} \left(\psi(-1,0)\right) + n_\beta, \end{equation}

for $n_\beta \in \mathbb {Z}$. A similar procedure at the endpoint $x = 1$ yields a similar expression for $\alpha$

(3.8)\begin{equation} \alpha = \frac{1}{\rm \pi} \cot^{-1} \left(\psi(1,0)\right) + n_\alpha, \end{equation}

for $n_\alpha \in \mathbb {Z}$. Physically speaking, we require a finite force when integrating the surface pressure round the leading edge, and that the Kutta condition holds at the trailing edge. Accordingly, we restrict $n_\alpha \geq 0$ and $n_\beta \geq -1$. Consequently, we seek an expansion of the vorticity distribution as a sequence of weighted Jacobi polynomials of the form

(3.9)\begin{equation} \hat{\gamma}_{a}(x) = \hat{\gamma}_{0} w^{\alpha,-\beta}(x)+ w^{\alpha,1 - \beta}(x) \sum_{n=1}^\infty \hat{\gamma}_{n} P_{n-1}^{\alpha,1-\beta} (x), \end{equation}

where $P_n^{a,b}$ represents the $n$th Jacobi polynomial with parameters $a$ and $b$. The Jacobi polynomials are a classical family of orthogonal polynomials (Szegö Reference Szegö1939) and represent a generalisation of Chebyshev polynomials. Some important properties of the Jacobi polynomials are catalogued in appendix A. In the unsteady case, the presence of an effective density in the porous boundary condition implies that the parameters $\alpha$ and $\beta$ may be complex valued. However, this scenario is not an issue, as Kuijlaars, Martinez-Finkelshtein & Orive (Reference Kuijlaars, Martinez-Finkelshtein and Orive2005) have established a theory for Jacobi polynomials with generalised parameters.

Note that the inverse cotangent function in (3.7) and (3.8) decreases monotonically for positive arguments. Therefore, the effect of porosity is to decrease the strength of both the leading-edge singularity and the trailing-edge zero. In the large porosity limit, $\psi \rightarrow \infty$, the singularity and zero vanish and we have $\alpha = \beta = 0$. Accordingly, the pressure jump along the chord also vanishes in this limit.

We may now substitute our Jacobi polynomial expansion (3.9) into the singular integral equation (3.1) and collocate at the Jacobi nodes to determine the coefficients $\hat {\gamma }_{n}$ following the procedure of Baddoo, Hajian & Jaworski (Reference Baddoo, Hajian and Jaworski2019) to furnish a solution to the full unsteady problem. This numerical technique is an example of a spectral method (Trefethen Reference Trefethen2000), where the unknown function is expanded globally in terms of basis functions whose coefficients are chosen by collocation.

Before we apply this strategy to the unsteady problem, we point out that an alternative weight function should be used for higher-accuracy solutions. For uniformly porous aerofoils with constant $\psi$, the remainder terms in the asymptotic expansions (3.6a) and (3.6b) are regular. Accordingly, the Jacobi weight function (3.5) precisely captures the behaviour at the endpoints. However, in the general case of non-constant-porosity profiles, the weight function should be written more generally as

(3.10)\begin{equation} w^{\alpha,\beta}(x)= (1-x)^{\alpha(x)} (1+x)^{\beta(x)}, \end{equation}

where $\alpha$ and $\beta$ are regular at $x = \pm 1$ (Hajian, Jaworski & Grace Reference Hajian, Jaworski and Grace2018). Although a set of orthogonal polynomials could in principle be constructed with the general weight function (3.10) via Gram–Schmidt orthogonalisation, we find it more appropriate to use the weight function (3.5) with the Jacobi polynomials due to the availability of many useful identities (see appendix A) and practical ease of computation. Whilst our choice of weight function precludes the possibility of spectral accuracy, only a few polynomials are usually required to obtain a degree of accuracy that is finer than the size of other physical quantities that are being ignored within the assumptions of the theoretical model.

3.2. Unsteady solution

We now adapt the steady solution (3.9) to the full unsteady problem (2.10). The SF–VIE in the unsteady case is distinct from the singular integral equation of the steady case in a number of ways. Firstly, the forcing term $f_{w}$ is not regular but possesses a logarithmic singularity at $x=1$. Secondly, the coefficient $\hat {\varGamma }$ multiplying $f_{w}$ is unknown a priori because it is proportional to the aerofoil circulation. Thirdly, the integral equation (2.10) now contains a term of Volterra type. We will now adapt the solution approach in § 3.1 to address these issues simultaneously.

We first address the fact that the forcing term $f_{w}$ is not regular as $x\rightarrow 1$ in (2.10). In particular, we have the logarithmic behaviour (Abramowitz & Stegun Reference Abramowitz and Stegun1964)

(3.11)\begin{equation} f_{w}(x) \sim -\frac{\textrm{i} k}{\rm \pi} \log(1-x), \quad \textrm{as } x \rightarrow 1. \end{equation}

For the asymptotic matching procedure at the endpoints, we therefore require that the left side of (2.10) possesses a logarithmic singularity of the same strength of $f_{w}$ at $x = 1$. The only way to generate this logarithmic singularity is through the principal value part of the operator $\mathcal {L}$ (Muskhelishvili Reference Muskhelishvili1946). In particular, we require

(3.12)\begin{equation} \hat{\gamma}_{a}(1) = -\textrm{i} k \hat{\varGamma}. \end{equation}

This behaviour may be alternatively derived by enforcing the Kutta condition and requiring the pressure to vanish at the trailing edge in (2.4a). Accordingly, we adapt the expansion in (3.4) and seek a solution of the form

(3.13)\begin{equation} \hat{\gamma}_{a}(x) =w^{\alpha,-\beta}(x) \hat{\gamma}_{a}^\ast(x) - \textrm{i} k 2^{\beta-1} \hat{\varGamma} w^{0,1-\beta}(x), \end{equation}

where $\hat {\gamma }^\ast _{a}$ is a smooth function. We now note the leading-order asymptotic behaviours as $x\rightarrow -1$:

(3.14a)\begin{gather} \int_{-1}^x \hat{\gamma}_{a}(\xi)\,\textrm{d} \xi \sim \frac{(1+x)^{1-\beta}2^{\alpha} }{1-\beta} \hat{\gamma}^\ast_{a}(-1), \end{gather}

(3.14b)\begin{gather}\psi(x,k)\hat{\gamma}_{a}(x) \sim \frac{\psi(-1,k)2^{\alpha} }{(1+x)^\beta} \hat{\gamma}^\ast_{a}(-1), \end{gather}

(3.14c)\begin{gather}\frac{1}{\rm \pi}\int\kern-10pt_{-1}^1 \frac{\hat{\gamma}_{a}(\xi)}{\xi - x} \,\textrm{d} \xi \sim \frac{\cot(\beta {\rm \pi})2^{\alpha} }{(1+x)^\beta }\hat{\gamma}^\ast_{a}(-1) +\varPhi^\ast(x), \end{gather}

(3.14d)\begin{gather}f_{a}(x) + \hat{\varGamma} f_{w}(x) \sim f_{a}(-1) + \hat{\varGamma} f_{w}(-1). \end{gather}

Consequently, we see that the Volterra part of the SF–VIE does not contribute to the asymptotic behaviour at the leading edge, and the expression for $\beta$ is the same as for the steady case (3.7) with $\psi (-1,0)$ replaced with $\psi (-1,k)$.

We now inspect the behaviour as $x \rightarrow 1$ and track terms at higher orders to ensure the correct asymptotic matching. In this limit, we obtain the behaviours

(3.15a)\begin{gather} \int_{-1}^x \hat{\gamma}_{a}(\xi)\,\textrm{d} \xi \sim \hat{\varGamma}, \end{gather}

(3.15b)\begin{gather}\psi(x,k) \hat{\gamma}_{a}(x) \sim \psi(1,k) \left(- \textrm{i} k \hat{\varGamma} + \frac{(1-x)^\alpha}{2^\beta} \hat{\gamma}_{a}^\ast(1) \right), \end{gather}

(3.15c)\begin{gather}\frac{1}{\rm \pi}\int\kern-10pt_{-1}^1 \frac{\hat{\gamma}_{a}(\xi)}{\xi - x} \,\textrm{d} \xi \sim -\frac{\textrm{i} k \hat{\varGamma }}{\rm \pi} \log(1-x) + \text{const.} + \cot(\alpha {\rm \pi})2^{-\beta}(1-x)^\alpha \hat{\gamma}_{a}(1), \end{gather}

(3.15d)\begin{gather}f_{a}(x) + \hat{\varGamma} f_{w}(x) \sim f_{a}(1) - \frac{\textrm{i} k \hat{\varGamma}}{\rm \pi} \log(1-x). \end{gather}

We note that the logarithmic singularities on the third and fourth lines cancel by virtue of the Kutta condition. The unknown constant $\hat {\varGamma }$ will be found through the collocation procedure, so we must choose $\alpha$ such that the leading-order zero (proportional to $(1-x)^\alpha$) vanishes. Consequently, we again obtain the same expression for $\alpha$ as in the steady case (3.8) with $\psi (1,0)$ replaced with $\psi (1,k)$.

Varying the flow resistance $\varPhi$ can have different effects on the Jacobi parameters $\alpha$ and $\beta$, as illustrated in figure 2(a). For very large values of $\varPhi$, the Jacobi parameters are close to $1/2$ and the aerofoil is almost impermeable. As $\varPhi$ is decreased, the real part of the Jacobi parameters tends to one of three values: if $k \rho _{e}<2$ the real part tends to 0, if $k \rho _{e} = 2$ it tends to $1/4$ and if $k \rho _{e}>2$ it tends to $1/2$. Similar behaviour can be observed as the effective density $\rho _{e}$ is increased in figure 2(b). The parameters transition from a lower value of approximately 0.1 to the impermeable limit of 1/2. This behaviour is due to the branch point of $\cot ^{-1}(z)$ at $z = \textrm {i}$.

Figure 2.

The real () and imaginary () parts of the first term of the Jacobi parameters $\alpha$ and $\beta$ as functions of porous medium resistivity $\varPhi$ and effective density $\rho _{e}$ for a range of reduced frequencies. The values of $\varPhi$ and $\rho _{e}$ here are taken to be constant with respect to chordwise position so $\alpha$ and $\beta$ values can be determined by a value of $\psi$ at the aerofoil endpoints in an application. In panel (a), $\rho _{e} = 1$ and $k$ is 0 (purple), 1, 2, 3, 4 (green) and in panel (b) $\varPhi = 1$ and $k = 0, 0.25, 0.5, 0.75, 1.$

The behaviour at the endpoints motivates an expansion of $\hat {\gamma }_{a}$ of the form

(3.16)\begin{equation} \hat{\gamma}_{a}(x) = - \textrm{i} k \hat{\varGamma} w^{0,1-\beta}(x)2^{\beta-1} + w^{\alpha,-\beta}(x) \hat{\gamma}_0 + w^{\alpha,1-\beta}(x) \sum_{n=1}^\infty \hat{\gamma}_n P_{n-1}^{\alpha,1-\beta}(x). \end{equation}

Recall that we do not know the circulation $\hat {\varGamma }$ a priori because it must be determined as part of the solution. However, integrating (3.16) yields

(3.17)\begin{align} \hat{\varGamma} = -2 \,\textrm{i} k \hat{\varGamma} B(2-\beta,1)+ \hat{\gamma}_0 2^{1+\alpha-\beta} B(1-\beta,1+\alpha) + \hat{\gamma}_1 2^{2+\alpha - \beta} B(2-\beta,1+\alpha), \end{align}

where we have used (A4) to express the result in terms of the beta function, $B$. Rearranging (3.17) then yields an equation for $\hat {\varGamma }$ in terms of the first two coefficients of the Jacobi expansion

(3.18)\begin{equation} \hat{\varGamma}= \frac{\hat{\gamma}_0 2^{1+\alpha-\beta} B(1-\beta,1+\alpha) + \hat{\gamma}_1 2^{2+\alpha - \beta} B(2-\beta,1+\alpha)}{1+2\,\textrm{i} k B(2- \beta,1)}. \end{equation}

It proves convenient to express

(3.19)\begin{equation} \hat{\varGamma} = \hat{\gamma}_0 \hat{\varGamma}_0 + \hat{\gamma}_1 \hat{\varGamma}_1, \end{equation}

so that the SF–VIE (2.10) may be expressed in the new form

(3.20)\begin{equation} \mathcal{H}\hat{\gamma}_{a}= f_{a}(x), \end{equation}

where the new, regularised operator $\mathcal {H}$ is defined as

(3.21)\begin{equation} \mathcal{H}\hat{\gamma}_{a} \triangleq \mathcal{L}\hat{\gamma}_{a}(x)- (\hat{\gamma}_0 \hat{\varGamma}_0 + \hat{\gamma}_1 \hat{\varGamma}_1) f_{{w}}(x). \end{equation}

It is straightforward to show that $\mathcal {H}\hat {\gamma }_{a}(x)$ is bounded at the endpoints. This regularisation improves the conditioning of the collocation matrix below.

We may now find approximate solutions for the coefficients $\hat {\gamma }_n$ by a collocation procedure. In particular, we truncate the infinite sum in (3.16) at $N$ and collocate the resulting (3.20) at the zeros of the Jacobi polynomial $P_{N+1}^{\alpha ,-\beta }$. The result is an $(N+1) \times (N+1)$ system of linear equations for the coefficients $\hat {\gamma }_n$. When constructing the system of equations, it is useful to note that the effect of the operator $\mathcal {L}$ on each individual weighted Jacobi polynomial can be computed using standard functions. For example,

(3.22)\begin{equation} \mathcal{L} w^{a,b} P_n^{a,b}(x) = \frac{w^{a,b}(x)}{\rm \pi} Q_n^{a,b}(x) -\textrm{i} k \psi(x,k) I_n^{a,b}(x) -\psi(x,k) w^{a,b}(x)P_n^{a,b}(x), \end{equation}

where $I_n^{a,b}$ is the integral of the weighted Jacobi (see (A7)) and $Q_n^{a,b}$ is the associated Jacobi function of the second kind (see (A5)). This result allows us to rapidly evaluate the collocation matrix below. The system of equations for the coefficients $\hat {\gamma }_n$ is given by $\boldsymbol {A} \hat {\boldsymbol {\gamma }} = \boldsymbol {f}_{a}$, where

(3.23)\begin{gather} \boldsymbol{A}^T= \begin{bmatrix} L_{0,0}^{\alpha,-\beta}+\hat{\varGamma}_0 (C L_{0,0}^{0,1-\beta}-f_{w}(x_0)) & \cdots & L_{0,N}^{\alpha,-\beta}+\hat{\varGamma}_0 (C L_{0,N}^{0,1-\beta}-f_{w}(x_N))\\ L_{0,0}^{\alpha,1-\beta}+\hat{\varGamma}_1 (C L_{0,0}^{0,1-\beta}-f_{w}(x_0)) & \cdots & L_{0,N}^{\alpha,1-\beta}+\hat{\varGamma}_1 (C L_{0,N}^{0,1-\beta}-f_{w}(x_N))\\ L_{1,0}^{\alpha,1-\beta} & \cdots & L_{1,N}^{\alpha,1-\beta}\\ \vdots & \ddots & \vdots\\ L_{N-1,0}^{\alpha,1-\beta} & \cdots & L_{N-1,N}^{\alpha,1-\beta} \end{bmatrix}, \end{gather}

(3.24a,b)\begin{gather}\hat{\boldsymbol{\gamma}} = \begin{bmatrix} \hat{\gamma}_0\\ \vdots\\ \hat{\gamma}_N \end{bmatrix}, \quad \boldsymbol{f}_{a} = \begin{bmatrix} f_{a}(x_0)\\ \vdots\\ f_{a}(x_N) \end{bmatrix}, \end{gather}

and where $C =- \textrm {i} k 2^{\beta - 1}$ and $L_{i,j}^{a,b} = \mathcal {L} w^{a,b}P^{a,b}_{i}(x_j)$. We solve this system in the least-squares sense using, for example, the backslash operation in Matlab.

For low reduced frequencies that are typically of interest ($k<3$), fewer than 10 Jacobi polynomials are usually required to resolve the vorticity distribution. As the reduced frequency increases, more polynomials are required to resolve an asymptotic inner region near the trailing edge; we comment on this scenario further in appendix B.

We now present several extensions to our method, including the case where the porosity distribution is discontinuous. The numerical method verifies the solution for quasi-steady aerodynamics and establishes both the circulatory and non-circulatory vorticity distributions for generalised unsteady aerodynamics of porous aerofoils.

3.3. Solution for discontinuous porosity distributions

The case of a discontinuous porosity profile is now considered. This scenario is motivated in part by the investigation by Geyer & Sarradj (Reference Geyer and Sarradj2014), who showed that, depending on the porous material, aerofoils with porosity at the trailing-edge section only can still achieve a noticeable noise reduction, while maintaining a certain level of aerodynamic performance over a fully porous aerofoil. A schematic of a partially porous aerofoil is illustrated in figure 3. When the discontinuity is located at $x=c$, the original SF–VIE (2.10) may be partitioned into two integral equations

(3.25)\begin{equation} \frac{1}{\rm \pi}\int\kern-10pt_{-1}^{1}\frac{\hat{\gamma}_{{a}}(\xi)}{\xi-x}\,\textrm{d} \xi - \psi_{i}(x,k) \left(\hat{\gamma}_{{a}}(x) + \textrm{i} k \int_{-1}^x \hat{\gamma}_{{a}}(\xi)\,\textrm{d} \xi \right) = f_{{a}}(x) + \hat{\varGamma} f_{{w}}(x) , \end{equation}

where

(3.26)\begin{align} \psi_{i}(x,k)= \begin{cases} \psi_{l}(x,k) & \text{for } -1 < x < c,\\ \psi_{r}(x,k) & \text{for }c < x < 1.\end{cases} \end{align}

Note that the subscripts $l$ and $r$ correspond the left and right sides of the discontinuity so that $\psi _{l}(c^-,k) \neq \psi _{r}(c^+,k)$.

Figure 3.

Example of an aerofoil with a discontinuous porosity distribution, where the discontinuity occurs at dimensionless chordwise position $x =c$. (a) Illustration of an aerofoil with a discontinuous porosity distribution where the porosity variation along the chord is indicated by at colour gradient. (b) A representative discontinuous porosity distribution. Note that the porosity along the forward section need not be constant.

We require the pressure jump across the wing to vanish at $x = c$ to ensure that there is no discontinuity in the seepage velocity (2.2). In particular, asymptotic analysis close to the discontinuity (Baddoo et al. Reference Baddoo, Hajian and Jaworski2019) reveals that, to leading order,

(3.27)\begin{equation} \left|\hat{\gamma}_{{a}} (x)+ \textrm{i} k \int_{-1}^x \hat{\gamma}_{{a}}(\xi)\,\textrm{d} \xi \right| \sim |x-c|^\lambda {\hat{\gamma}}^\ast_{{a}}(x), \quad {\rm as}\;x \rightarrow c, \end{equation}

where

(3.28)\begin{equation} \lambda = \frac{1}{\rm \pi} [ \cot^{-1} (\psi_{l}(c^-,k)) - \cot^{-1} (\psi_{r}(c^+,k)) ], \end{equation}

and where ${\hat {\gamma }}^\ast _{{a}}(x)$ is regular in $-1 < x < c$ and $c < x < 1$ but may be discontinuous at $x=c$.

These results reflect some physical characteristics of partially porous aerofoils. On one hand, if the junction transitions from less permeable to more permeable ($\Re [\psi _{l}(c^-,k)] < \Re [\psi _{r}(c^+,k)]$) then the pressure vanishes at $c$ and the junction behaves as a second trailing edge. On the other hand, if the junction transitions from more permeable to less permeable ($\Re [\psi _{l}(c^-,k)] > \Re [\psi _{r}(c^+,k)]$) then $\lambda <0$ and the pressure is singular at the junction. Therefore, in this case the junction behaves like a second leading edge.

These observations motivate two separate expansions for $\hat {\gamma }_{{a}}$ in the left and right regions of the forms

(3.29a)\begin{gather} \hat{\gamma}_{l}(\tau_{l}) = \hat{\gamma}_{l,0} w^{\lambda,-\beta}(\tau_{l})+ 2^{\beta-1} \varPi w^{0,1- \beta}(\tau_{l})+ w^{\lambda,1-\beta}(\tau_{l})\sum_{n=1}^\infty \hat{\gamma}_{{l},n} P_{n-1}^{\lambda,1-\beta} ( \tau_{l} ), \end{gather}

(3.29b)\begin{gather}\hat{\gamma}_{r}(\tau_{r}) = 2^{-\alpha} \varPi w^{\alpha,0}(\tau_{r}) + \varLambda w^{0,\lambda}(\tau_{r}) + w^{\alpha,\lambda}(\tau_{r}) \sum_{n=1}^\infty \hat{\gamma}_{{r},n} P_{n-1}^{\alpha,\lambda} ( \tau_{r}), \end{gather}

where $\varPi$ and $\varLambda$ are constants, and we have introduced the rescaled variables

(3.30a)\begin{gather} \tau_{l}(x) = -1+2 \left(\frac{x+1}{1+c}\right), \quad -1 < x < c, \end{gather}

(3.30b)\begin{gather}\tau_{r}(x) = 1+ 2 \left(\frac{x-1}{1-c}\right), \quad c < x < 1, \end{gather}

so that $-1<\tau _{l},\, \tau _{r} <1$. We now seek to express the constants $\varPi$ and $\varLambda$ in (3.29a) and (3.29b) in terms of the unknown coefficients $\hat {\gamma }_{{l},n}$ and $\hat {\gamma }_{{r},n}$. Beginning with the constant $\varPi$, we note the two relations

(3.31a,b)\begin{equation} \hat{\gamma}_{{a}}(c) = \varPi \quad \textrm{and}\quad \hat{\gamma}_{{a}}(c)=-\textrm{i} k \int_{-1}^c \hat{\gamma}_{{a}}(\xi)\,\textrm{d} \xi. \end{equation}

The latter expression may be evaluated using the expansion (3.29a) and the quadrature formula (A4). A simple rearrangement then allows us to express $\varPi$ in terms of $\hat {\gamma }_{l,0}$ and $\hat {\gamma }_{l,1}$ as

(3.32)\begin{equation} \varPi = \varPi_0 \hat{\gamma}_{l,0} + \varPi_1 \hat{\gamma}_{l,1}, \end{equation}

where

(3.33a)\begin{gather} \varPi_0 = \frac{-\textrm{i} k (1+c)2^{\lambda -\beta}B(1-\beta,1+\lambda)}{1+ \textrm{i} k (1+c) B(2- \beta,1)}, \end{gather}

(3.33b)\begin{gather}\varPi_1 = \frac{-\textrm{i} k (1+c)2^{\lambda +1 - \beta}B(1-\beta,1+\lambda)}{1+ \textrm{i} k (1+c) B(2 - \beta,1)}. \end{gather}

We now seek to express $\varLambda$ in terms of the coefficients $\hat {\gamma }_{{l},n}$ and $\hat {\gamma }_{{r},n}$. By employing an approach similar to § 3.2, it is straightforward to show that the new expression for $\varLambda$ is

(3.34)\begin{equation} \varLambda = \varLambda_{l,0} \hat{\gamma}_{l,0} + \varLambda_{l,1} \hat{\gamma}_{l,1} + \varLambda_{r,1} \hat{\gamma}_{r,1}, \end{equation}

where

(3.35a)\begin{gather} \varLambda_{l,0} = \varPi_0 \mathcal{M} (\textrm{i}/k + (1-c)B(1,1+\alpha)), \end{gather}

(3.35b)\begin{gather}\varLambda_{l,1} = \varPi_1 \mathcal{M} (\textrm{i}/k + (1-c)B(1,1+\alpha)), \end{gather}

(3.35c)\begin{gather}\varLambda_{r,1} = \mathcal{M} (1-c) 2^{\alpha+\lambda} B(1+\alpha,1+\lambda), \end{gather}

(3.35d)\begin{gather}\mathcal{M} = \frac{ -\textrm{i} k 2^{-\lambda}}{1+\textrm{i} k (1-c)B(1,1+\lambda)}. \end{gather}

We may now substitute the expansions (3.29a) and (3.29b) into (3.25). By collocating at the Jacobi nodes on the forward and aft sections, we obtain a system of linear equations for the unknown coefficients. During the procedure we encounter the Cauchy integral of the weighted Jacobi polynomials without the principal value, which can be calculated using (A6).

Although only a single discontinuity was considered in this example, any finite number of discontinuities could be modelled using the same approach.

3.4. Circulatory and non-circulatory solutions

The solution to the full unsteady problem in § 3.2 may be separated into circulatory and non-circulatory parts by writing

(3.36)\begin{equation} \hat{\gamma}_{{a}}(x) = \hat{\gamma}_{{a}}^C(x) + \hat{\gamma}_{{a}}^{NC}(x), \end{equation}

where the superscripts $C$ and $NC$ denote the circulatory and non-circulatory contributions, respectively. The circulatory part is sometimes referred to the wake-induced component because it contains information about the effect of the downstream wake on the aerofoil. Conversely, the non-circulatory part is sometimes referred to as the added mass component as it represents the effects of the unsteady sloshing of the flow about the aerofoil. Recent research into the origin of added mass led Leonard & Roshko (Reference Leonard and Roshko2001) and Eldredge (Reference Eldredge2010) to postulate that its associated force may be represented solely by inviscid theory, even in viscous and separated flows. Corkery, Babinsky & Graham (Reference Corkery, Babinsky and Graham2019) confirmed experimentally the ability of inviscid theory to represent added mass effects as a non-circulatory component that depends only on body geometry and its motion, where the circulatory terms in turn measure the viscous effects associated with the bound vorticity and the wake. The circulatory and non-circulatory components combine to give the full vorticity distribution on the aerofoil.

As its name suggests, the non-circulatory component is the solution to (2.10) subject to the auxiliary requirement that its net circulation is identically zero,

(3.37)\begin{equation} \int_{-1}^1 \hat{\gamma}^{NC}_{{a}}(\xi)\,\textrm{d} \xi = 0. \end{equation}

The problem of finding $\hat {\gamma }_{a}^{NC}$ subject to the SF–VIE (2.10) and to both the non-circulatory condition (3.37) and the Kutta condition is generally ill posed. In other words, the Kutta condition (2.9) cannot be applied to the circulatory and non-circulatory solutions individually. As such, we permit singularities at the trailing edge in both $\hat {\gamma }_{{a}}^C$ and $\hat {\gamma }_{{a}}^{NC}$. These singularities are perfectly valid and appear in the impermeable case detailed in, for example, Bisplinghoff et al. (Reference Bisplinghoff, Ashley and Halfman1996, §§ 5 and 6). The Kutta condition for the full solution is then enforced by specifying that these singularities are equal and opposite and ${\rm \Delta} p(1)=0$ when the circulatory and non-circulatory solutions are combined according to (3.36), thereby connecting the bound circulation to the wake strength.

It is simpler to derive the non-circulatory solution and then use the full solution and (3.36) to determine the circulatory solution. Following the analysis of § 3.2, we seek an expansion of the form

(3.38)\begin{equation} \hat{\gamma}_{{a}}^{NC}(x) = \varTheta w^{\alpha-1,1-\beta}(x) + w^{\alpha,-\beta}(x) \hat{\gamma}_0^{NC} + w^{\alpha,1-\beta}(x) \sum_{n=1}^N \hat{\gamma}_n^{NC} P_{n-1}^{\alpha,1-\beta}(x), \end{equation}

where the constant $\varTheta$ is chosen so that the circulation vanishes

(3.39)\begin{equation} \varTheta = \frac{-\hat{\gamma}_0^{NC} B(1+\alpha,1-\beta) - \hat{\gamma}_1^{NC} 2B(1+\alpha, 2-\beta)}{B(\alpha,2-\beta)}. \end{equation}

At this point a collocation scheme similar to the procedure described in § 3.2 can determine the coefficients ${\hat {\gamma }_n^{NC}}$.

3.5. Quasi-steady solution

The quasi-steady problem is equivalent to the steady problem described in Hajian & Jaworski (Reference Hajian and Jaworski2017) with the exception that the kinematic boundary condition is replaced by the instantaneous unsteady boundary condition (2.1). Specifically, the quasi-steady assumption augments the angle of attack due to aerofoil pitch or camber with the effective angle of attack from heaving motions. With this modification, the singular integral equation (2.6) becomes

(3.40)\begin{equation} \frac{1}{\rm \pi}\int\kern-10pt_{-1}^{1}\frac{\hat{\gamma}_{{a}}^Q(\xi)}{\xi-x}\,\textrm{d}\xi - \psi(x,0) \hat{\gamma}_{{a}}^Q(x) = 2\left( \textrm{i} k \hat{y}_{{a}}(x)+ \frac{\textrm{d} \hat{y}_{{a}}}{\textrm{d} x} (x) \right), \end{equation}

where the superscript $Q$ denotes the quasi-steady solution. This singular integral equation may be solved analytically using classical techniques (Muskhelishvili Reference Muskhelishvili1946). The solution is rendered unique by enforcing the Kutta condition, and the corresponding solution is

(3.41)\begin{align} \hat{\gamma}_{{a}}^Q(x) &= \frac{-2}{1+(\psi(x,0))^2} \left\{\psi(x,0) \left(\textrm{i} k \hat{y}_{{a}}(x)+ \frac{\textrm{d} \hat{y}_{{a}}}{\textrm{d} x} (x) \right) \right.\nonumber\\ &\quad + \left.\frac{Z(x)}{\rm \pi} \int\kern-10pt_{-1}^1 \left( \textrm{i} k \hat{y}_{{a}}(\xi)+ \frac{\textrm{d} \hat{y}_{{a}}}{\textrm{d} x} (\xi) \right)\frac{\textrm{d} \xi }{Z(\xi)(\xi - x)} \right\}, \end{align}

where

(3.42)\begin{equation} Z(x) = \sqrt{1 + (\psi(x,0))^2} \exp \left[\int\kern-10pt_{-1}^1 \frac{j(\xi)}{\xi - x}\,\textrm{d} \xi\right], \end{equation}

with $j(x) = (1/{\rm \pi} ) \cot ^{-1}(\psi (x,0))$.

In special cases, such as uniformly porous aerofoils with simple camber descriptions, the singular integrals in (3.41) can be calculated analytically, and the full solution can be expressed in closed form. Otherwise, the solution (3.41) contains nested singular integral equations that can be computed numerically using Gauss–Jacobi quadrature.

4.1. Continuous porosity distributions

The numerical method in this work is amenable to any Hölder continuous porosity distribution. Using (2.4a) and (A7) to construct the aerofoil pressure jump from the numerical scheme, figure 4 compares the magnitudes of the unsteady pressure distributions for aerofoils with linear and parabolic reciprocal flow resistance distributions (cf. (2.12)) undergoing pitching or heaving motions of unit amplitude. The Kutta condition is clearly satisfied at the trailing edge for all cases presented. Figure 4 indicates that the introduction of porosity decreases the pressure jump across the aerofoil under both pitching and heaving motions. As the porosity parameter $\psi$ increases along the chord, the pressure distribution decreases except in a small region localised to the trailing edge. This reduction corresponds to a significant reduction in the unsteady lift. We also note the rapid changes in pressure at the trailing edge $x=1$, which is caused by the reduction in the strength of the trailing-edge zero by (3.8). This behaviour is associated with the reduction in vortex shedding at the porous trailing edge.

Figure 4.

Magnitude of the unsteady pressure distributions for a porous aerofoil with continuous porosity distributions undergoing harmonic motions with unit amplitude: (a) pitching about the leading edge with linear reciprocal dimensionless flow resistance for $k=0.5$; (b) heaving with parabolic reciprocal dimensionless flow resistance for $k=0.1$. $\rho _{e}=1.2$ in both cases. The impermeable limit is indicated by the thick black line.

4.2. Discontinuous porosity distributions

The numerical method in § 3.3 is now applied to aerofoils with discontinuous porosity distributions. Specifically, the numerical scheme is demonstrated for aerofoils with an impermeable leading-edge section for $x<0$ and a constant-porosity section for $x>0$. Figure 5 plots the surface pressure jump for pitching motions at $k=0.5$ about the leading edge ($x=-1$) and heaving motions at $k=0.1$. The seepage velocity is continuous across the chord, as evidenced by the vanishing jump at pressure located at the junction $x = 0$. Note that even a small amount of porosity is sufficient to enforce a zero pressure jump at the junction.

Figure 5.

Magnitude of the unsteady pressure distributions on a partially porous aerofoil undergoing pitching or heaving motions with unit amplitude: (a) pitching about the leading edge with linear reciprocal dimensionless flow resistance for $k=0.5$; (b) heaving with parabolic reciprocal dimensionless flow resistance for $k=0.1$. The porosity distributions are plotted on a logarithmic scale in the figure insets, where $\psi = 0$ for $x<0$. The impermeable case is indicated by the thick black line. $\rho _{e}=1.2$ in both cases.

For both the pitching and plunging cases and at different reduced frequencies, the pressure jump along the permeable section decreases when the porosity is large. Accordingly, the permeable section of the wing behaves effectively as an extension of the wake: even if the pressure jump vanishes, such as in the orange curve in figure 5(b), the vorticity distribution does not, by (2.4a). For these cases with an impermeable forward section, the aerofoil solution for large porosity in the aft section is the same as the solution of an impermeable aerofoil truncated at $x = c$, i.e. the solution attained as rescaling the characteristic length scale by a factor of $(c+1)/2$, or equivalently, using the mappings

(4.1a–c)\begin{equation} x \mapsto \frac{(c+1)x+(c-1)}{2}, \quad y_{a}(x) \mapsto \frac{c+1}{2}y_{a}, \quad k \mapsto \frac{c+1}{2} k. \end{equation}

6.1. Theodorsen function

The Theodorsen function (Theodorsen Reference Theodorsen1935) may be interpreted as the ratio of the wake-induced (circulatory) lift $L^C$ to the quasi-steady lift $L^Q$ (cf. Taha & Rezaei Reference Taha and Rezaei2019; Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 279),

(6.1)\begin{equation} C(k) = \frac{L^C}{L^{Q}}.\end{equation}

This ratio can be expressed in closed form for an impermeable aerofoil as

(6.2)\begin{equation} C(k) = \frac{K_1(\textrm{i} k)}{K_0(\textrm{i} k) + K_1(\textrm{i} k)} \quad \mbox{for}\ \psi(x,k)\, \equiv 0,\end{equation}

where $K_0$ and $K_1$ are modified Bessel functions of the second kind. The lack of an analytic solution to the SF–VIE (2.10) precludes the derivation of a similar expression for the permeable case in terms of standard functions. However, we may use the numerical solution presented in § 3 to compute the relevant circulatory and quasi-steady lift quantities to construct numerically a porous extension to the Theodorsen function via (6.1). Specifically, the wake-induced lift is computed by the numerical scheme in § 3.4, whereas the quasi-steady lift may be determined analytically using the solution presented in § 3.5. Results in this section are constructed using harmonically heaving motions and a particular family of porosity distributions, although other motions and porosity profiles can be examined using our numerical method.

The porosity profile (2.12) is set by two dimensionless properties of the porous medium, flow resistance $\varPhi$ and effective density $\rho _{e}$, each of which has different effects on the relationship between the wake and quasi-steady lift. Figure 7 illustrates these separate effects of $\varPhi$ and $\rho _{e}$ on the magnitude and phase of the porous extension of the Theodorsen function, whose classical result for impermeable aerofoils (6.2) is indicated by the thick black curves. The left column of figure 7 examines the effect of reducing the dimensionless flow resistance of the porous medium. For these porosity configurations, the magnitude of the Theodorsen function is larger than the traditional, impermeable Theodorsen function at low reduced frequency; however, this trend reverses at high frequency, as illustrated in figure 7. Given that aerofoil porosity cannot increase the magnitude of the quasi-steady lift (cf. § 3.5), figure 7 implies that an effect of the dimensionless flow resistance is to modify the wake-induced lift. This observation is consistent with the weakening of the trailing-edge zero due to porosity at that position, as discussed in § 3.1.

Figure 7.

Complex, magnitude and phase plots of the porous Theodorsen function for a range of flow resistance distributions (a,c,e) and effective densities (b,d,f). The dimensionless flow resistance is given by the reciprocal linear profile $\varPhi =1/(\mu (1+x))$ so that the leading edge is impermeable in every case. On the left, $\mu$ ranges from 0, 0.05, 0.10, 0.15, 0.2 (black to orange) and $\rho _{e}=1.5$; on the right, the effective density takes values $\infty$, 4, 3, 2, 1 (black to green) and $\mu =0.05$. Thus, the black curve is the classical Theodorsen function for impermeable aerofoils. In the complex plots, the points $k=0,10^{-3},10^{-2},10^{-1},1,10$ are indicated by $\square$ with $k=0$ representing the rightmost part of the curve.

The right column of figure 7 fixes the chordwise flow resistance distribution $\varPhi (x)$ and steps through a set of constant effective density values representative of the full range of physically relevant values. The influence of effective density on the magnitude of the Theodorsen function is marginal at low reduced frequencies, say, below $k \approx 1$, as would be expected of the term arising from the unsteady contribution of the aerofoil boundary condition (2.2). The role of effective density becomes pronounced at larger reduced frequencies, and the magnitude of the porous Theodorsen function decreases with decreasing effective density.

Although the magnitude of the Theodorsen function is relatively insensitive to changes in dimensionless flow resistance and effective density below $k \approx 1$, the bottom row of figure 7 demonstrates that the phase depends strongly on these parameters at lower reduced frequencies and in a complicated way. The reduction of effective density increases the maximum value of phase lag and the reduced frequency at which it occurs. Smaller values of the dimensionless flow resistance of the porous medium lead to larger maximum phase lag but lower reduced frequencies at these extrema. Viewed together, the phase lag of the porous Theodorsen function is greater and has peaks at larger reduced frequencies than the classical impermeable scenario. The phase angle of the Theodorsen function plays a critical role in the energy transfer between the fluid flow and aerofoil motions, and the above observations suggest that the porosity distribution of the aerofoil may be tuned to shift the flutter boundary of an aeroelastic aerofoil in a desired way (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, pp. 279–280).

Lastly, we note in the first row of figure 7 that the ratio between the wake-induced lift and quasi-steady lift tends to a finite, real value with a large reduced frequency that will depend on both the non-dimensional flow resistance $\varPhi$ and effective density $\rho _{e}$ of the porous material.

However, caution should be exercised in the large-frequency limit, which is beyond the realm of physical validity of the mathematical modelling employed in this paper: as discussed by Howe (Reference Howe1979) and noted more recently by Weidenfeld & Manela (Reference Weidenfeld and Manela2016), the present modelling assumptions for porosity are only valid when there is Stokes flow in the pores passing through the wing, which is rendered invalid at large frequencies. We expect that a higher-order (e.g. quadratic) porosity law would yield more physically meaningful results at high frequencies; an exploration of such porosity models is beyond the scope of the present work and is not pursued here.

6.2. Sears function

We now consider a uniform flow in the horizontal direction with a transverse sinusoidal gust, as illustrated in figure 8. The gust convects at the free-stream velocity and has unit amplitude and reduced frequency $k$. The interaction between the harmonic gust generates unsteady lift on the stationary aerofoil, as well as pressure perturbations that propagate to the acoustic far field as sound waves. The unsteady lift response is described by the Sears function (Sears Reference Sears1941), which may be written in closed form for an impermeable aerofoil as

(6.3)\begin{equation} S(k) = C(k) [ J_0(k)- \textrm{i} J_1(k) ] + \textrm{i} J_1(k) \quad \mbox{for } \psi(x,k)\, \equiv 0, \end{equation}

where $J_n$ are Bessel functions of the first kind of order $n$.

Figure 8.

The harmonic gust response problem considered by Sears (Reference Sears1941). A symmetric, stationary porous aerofoil is subjected to an unsteady vertical sinusoidal disturbance convected at the velocity of the free stream, resulting in unsteady wake generation. Acoustic waves are also generated that can be modelled using the Sears response function.

To construct a porous extension of the Sears function, the harmonic gust of unitary amplitude enters the porous boundary condition (2.1) as

(6.4)\begin{equation} w(x,t) = w_s(x,t) -\textrm{e}^{\textrm{i} k (t-x) }. \end{equation}

Therefore, the forcing function $f_{a}$ in the SF–VIE becomes

(6.5)\begin{equation} f_{a}(x) = \textrm{e}^{-\textrm{i} k x}, \end{equation}

and we may apply the numerical scheme developed in § 3 to find the resulting total unsteady lift: the Sears function. Whilst the Volterra part of the SF–VIE renders it impossible to find analytic forms for the Sears function, our numerical scheme is sufficiently fast and robust that we may produce an accurate approximation for a range of porosity gradients.

We now explore the effects of aerofoil porosity distributions illustrated in figure 9 on the Sears gust response function. Figure 10 presents the complex-valued information of the porous Sears function analogue as a function of the dimensionless flow resistance $\varPhi$ and effective density $\rho _{e}$ in the left and right columns, respectively. The porosity distribution has a fixed value of unity at the trailing edge for all cases considered, and the permeable length of the aerofoil is varied. We plot the unsteady lift normalised by the quasi-steady lift particular to the given porosity distribution, such that the Sears function asymptotes to unity for vanishing reduced frequency. This quasi-steady lift follows from (3.41) with the forcing on the right-hand side of (3.40) replaced by (6.5). Note that the phase comparison in the bottom row of figure 10 involves a factor of $\textrm {e}^{-\textrm {i} k}$ that unwinds the Argand diagram representation of the Sears function to more easily distinguish differences due to $\varPhi$ and $\rho _{e}$.

Figure 9.

The (reciprocal) piecewise linear flow resistance profiles used to compute the porous Sears function. The colours correspond to the curves in figure 10.

Figure 10.

Complex, magnitude and phase plots of the porous Sears function for a range of flow resistance distributions (left) and effective densities (right). The dimensionless flow resistance is given by a reciprocal piecewise linear profiles illustrated in figure 9. The leading edge is impermeable in every case. On the left, the colours correspond to the flow resistance profiles in figure 9 and $\rho _e=1.5$; on the right, the effective density takes values $\infty$, 4, 3, 2, 1 (black to green) and the flow resistance profile is given by the rightmost curve in figure 9. The black curves represent the impermeable Sears function. The points $k=0,0.1,1,2,5,10$ are indicated by $\square$ with $k=0$ representing the rightmost part of the curve.

The results show that the magnitude of the porous Sears function increases with decreasing flow resistance. However, the actual unsteady lift is smaller at a given reduced frequency when compared to the impermeable case, as the Sears function must be multiplied by the quasi-steady lift introduced in the normalisation. Changes in the phase relative to the impermeable aerofoil limit are most pronounced in this porosity configuration for reduced frequencies between O($10^{-1}$) and O(1) that lie in the typical range of aerospace interest. Note the crossover in the porous Sears function phase curves with the impermeable limit at a particular value of reduced frequency, which occurs for variations in both the dimensionless flow resistance and effective density in the bottom row of figure 10. Effective density variations do not lead to significant changes in the porous Sears function for the cases considered here, where the deviation between the porous and impermeable results in the right column of figure 10 depend on the dimensionless flow resistance. At the limits of very small or large reduced frequency, the curves for the porous Sears function collapse to the classical impermeable result in terms of both phase and amplitude.

6.3. Wagner function

With numerical solutions in hand for the porous Theodorsen and Sears functions, we now turn our attention to time-domain aerodynamic responses. It is sufficient in this linearised physical system to determine indicial functions, which are the aerodynamic responses to impulsive system changes. Duhamel's integral may then be used to construct the aerodynamic response to an arbitrary aerofoil motion or gust field (Fung Reference Fung1993, pp. 277–280). Wagner (Reference Wagner1925) first solved the canonical problem of the transient lift experienced by an aerofoil whose angle of attack changes instantaneously from zero to a fixed value; herein we develop the porous analogue to the so-called Wagner function.

Recalling the non-dimensional scalings chosen in § 2, the linearised normal wing velocity for an impulsive unit change in angle of attack modifies the boundary condition (2.1) to

(6.6)\begin{equation} w(x,t) = w_{s}(x,t) + H(t). \end{equation}

Note that the Heaviside function $H$ may be expressed in the frequency domain as

(6.7)\begin{equation} H(t) = \int_{-\infty}^{\infty} \frac{\textrm{e}^{\textrm{i} k t}}{\textrm{i} k} \,\textrm{d} k, \end{equation}

where the path of integration passes below the pole at $k=0$. Carrying through the analysis of § 3.2 with boundary condition (6.6) furnishes the circulatory lift history due to an impulsive change in angle of attack,

(6.8)\begin{equation} L^C(t) = L^{Q} \int_{-\infty}^\infty \frac{C(k)}{\textrm{i} k}\, \textrm{e}^{\textrm{i} k t} \,\textrm{d} k. \end{equation}

The quasi-steady lift $L^Q$ corresponding to the asymptotic lift value at long times for a given porosity profile is

(6.9)\begin{equation} L^{Q} = \int_{-1}^1 \frac{4}{1+(\psi(x,0))^2} \left\{\psi(x,0) +\frac{Z(x)}{\rm \pi} \int\kern-10pt_{-1}^1\frac{\textrm{d} \xi }{Z(\xi)(\xi - x)} \right\}{\textrm{d} x}, \end{equation}

where $Z$ is defined in (3.42). The ratio $L^C(t)/L^{Q}$ from (6.8) defines the Wagner function,

(6.10)\begin{equation} \phi(t)= \int_{-\infty}^\infty \frac{C(k)}{\textrm{i} k} \,\textrm{e}^{\textrm{i} k t} \textrm{d} k. \end{equation}

Therefore, the porous Theodorsen function analogue in the frequency domain from § 6.1 generates the associated time-domain Wagner function for impulsive aerofoil motions.

Figure 11 plots the porous Wagner function against convective time for the same ranges of dimensionless flow resistance and effective density examined in figure 7, i.e. varying the flow resistance in figure 11(a) and setting different constant effective density values in figure 11(b). In each case, the Wagner function asymptotes to unity at large time, which is anticipated from (6.10) by the fact that the Theodorsen function is unity at $k=0$ regardless of the porosity profile. Aerofoil porosity decreases the magnitude of the Wagner function at short times, whose precise value depends on the porosity parameters. Figure 11(b) makes clear that effective density controls the short-time behaviour of the Wagner function and has a marginal influence at large times. However, figure 11(a) indicates counterintuitively that a smaller flow resistance of the porous aerofoil leads to a larger initial value of the Wagner function, and a lesser flow resistance yields a faster approach to the quasi-steady lift value. Additionally, figures 11(a) and 11(b) show that the late-time behaviour is controlled by the flow resistance, whereas the early-time behaviour is controlled by the effective density. These can observations can be interpreted in light of the asymptotic solutions derived in appendix B. Therein, we show that, to leading order, the low-frequency (long-time) lift depends only on the flow resistance of the porous medium, whereas as the high-frequency (short-time) lift depends only on the effective density.

Figure 11.

Porous extension of the Wagner function $\phi (t)$ for an impulsive change in angle of attack as a function of convective time. The ranges of the dimensionless flow resistance $\varPhi$ and effective density $\rho _{e}$ of the porous aerofoil are the same as in figure 7. The black curve indicates the classical Wagner function for an impermeable aerofoil. In (a) the flow resistance ranges from high (black) to low (orange), and in (b) the effective density varies from high (black) to low (yellow).

Figure 12 clarifies the roles of the porosity parameters by considering only the unsteady circulatory lift, which is also the metric of practical interest. We show here that the dimensionless flow resistance controls the rise time and asymptotic value of the circulatory lift at long times, and the effective density governs the short-time aerodynamic response only. Therefore, the effective density of the porous medium is the key parameter to reduce the transient effect of aerofoil motion on the resulting unsteady lift, where the flow resistance distribution may be set by steady aerodynamic constraints of a particular aerial system.

Figure 12.

Circulatory lift history due to an impulsive change in angle of attack of a porous aerofoil as a function of convective time. The ranges of the dimensionless flow resistance $\varPhi$ and effective density $\rho _{e}$ of the porous aerofoil are the same as in figure 7. The black curve corresponds to an impermeable aerofoil. In (a) the flow resistance ranges from 0 (black) to high (orange), and in (b) the effective density varies from high (black) to low (yellow).

In complement to the passive unsteady aerodynamic control that tailored porosity distributions could provide, we note the possible predictive limitations that are tied to the underlying linear assumptions of the flow field. Beckwith & Babinsky (Reference Beckwith and Babinsky2009) demonstrate experimentally that the impermeable Wagner function predicts well the lift history at long times for small angles of attack, where measured early-time responses overshoot the theoretical predictions due to delayed stall effects not present in the Wagner model. Larger impulsive changes to the angle of attack lead to persistent nonlinear flow structures that must be integrated into the aerodynamic lift model to improve its predictive capacity (Pitt Ford & Babinsky Reference Pitt Ford and Babinsky2014; Li & Wu Reference Li and Wu2015; Stevens & Babinsky Reference Stevens and Babinsky2017). However, unsteady porous aerofoils may introduce important nonlinear flow features in real flows even at small angles, and an experimental campaign beyond the scope of the present work is needed to examine the relevant fluid mechanics of porous aerofoils and the predictive capability of the present model.

6.4. Küssner function

The Küssner function (Küssner Reference Küssner1936) is the indicial aerodynamic response to a stationary aerofoil encountering a sharp-edged gust. Like the relationship between the Theodorsen function and the Wagner function, the Sears harmonic gust response function is connected via the Fourier transform to the Küssner sharp-edged gust function, which can be used to predict the aerodynamic response to arbitrary linear gusts by appeal to Duhamel's integral. The sharp-edged gust problem is carried out similarly to § 6.2 with boundary condition (2.1) now given by

(6.11)\begin{equation} w(x,t) = w_s(x,t) + H(x-t+1). \end{equation}

The frequency-domain representation of the Heaviside function (6.7) may be used again to arrive at the transient total lift on the aerofoil,

(6.12)\begin{equation} L(t) = L^Q \int_{-\infty}^\infty \frac{S(k)}{\textrm{i} k} \,\textrm{e}^{\textrm{i} k (t-1)} \,\textrm{d} k. \end{equation}

Note the dependence of the unsteady lift on the Sears function $S(k)$, into which all of the information about the porosity distribution is embedded. The Küssner function is the ratio $L(t)/L^Q$ of the unsteady lift to the quasi-steady lift (6.9),

(6.13)\begin{equation} \psi(t) = \int_{-\infty}^\infty \frac{S(k)}{\textrm{i} k} \,\textrm{e}^{\textrm{i} k (t-1)} \,\textrm{d} k, \end{equation}

where the contour of integration passes below the origin.

Figures 13 and 14 illustrate the effects of the dimensionless flow resistance $\varPhi$ and effective density $\rho _{e}$ of the porous aerofoil on the Küssner function and on the total unsteady lift, respectively. The descriptions of $\varPhi$ and $\rho _{e}$ are the same as in figure 10, where the porous Sears functions computed in § 6.2 generate the porous Küssner functions in figures 13 and 14. The results in these figures differ by a factor of the quasi-steady lift that depends on the flow resistance distribution, where accelerated rise times in figure 13(a) for aerofoils with less resistance have instead slower actual rise times in the lift coefficient and a smaller asymptotic lift value in figure 12(a), as discussed previously in the context of the Wagner function. Figures 13(b) and 14(b) show that the effective density has a marginal influence on the sharp-edged gust response. Therefore, the dimensionless flow resistance distribution has the dominant influence on the porous Küssner function.

Figure 13.

Porous extension of the Küssner function $\psi (t)$ for a sharp-edged gust as a function of convective time. The ranges of dimensionless flow resistance $\varPhi$ and effective density $\rho _{e}$ of the porous aerofoil are the same as in figure 10. The black curve corresponds to the classical Küssner function for an impermeable aerofoil. In (a) the flow resistance ranges from high (black) to low (orange) and in (b) the effective density varies from high (black) to low (yellow).

Figure 14.

Circulatory lift history of a stationary porous aerofoil encountering a sharp-edged gust as a function of convective time. The ranges of dimensionless flow resistance $\varPhi$ and effective density $\rho _{e}$ of the porous aerofoil are the same as in figure 10. The black curve indicates results for an impermeable aerofoil. In (a) the flow resistance ranges from high (black) to low (orange) and in (b) the effective density varies from high (black) to low (yellow).

Comparisons between the impermeable Küssner function to both experimental and computational simulations (Biler, Badrya & Jones Reference Biler, Badrya and Jones2019; Sedky, Lagor & Jones Reference Sedky, Lagor and Jones2020) are favourable for gust amplitudes of the order of or smaller than the free-stream velocity, i.e. moderate to small gust ratios. Larger gust ratios can lead to signification deviations from the linear theory, especially at long times due to large-scale vortex structures generated by the gust encounter (Perrotta & Jones Reference Perrotta and Jones2017; Andreu-Angulo et al. Reference Andreu-Angulo, Babinsky, Biler, Sedky and Jones2020). One might anticipate that the robustness of the Küssner function for gust encounters with significant nonlinear flow features might carry over to porous aerofoils, whose response in short times is similar to the impermeable case. Experimental unsteady lift measurements of porous aerofoils that are beyond the scope of the present work would be necessary to assess this claim and the effects of porosity on gust response in real flows.