A stochastic estimation framework for interpretable force modelling in flapping-wing aerodynamics

Abstract

Unsteady aerodynamic forces in flapping wings arise from complex, nonlinear flow structures that challenge predictive modelling. In this work, we introduce a data-driven framework that links experimentally observed flow structures to sectional pressure loads on physical grounds. The methodology combines proper orthogonal decomposition and quadratic stochastic estimation (QSE) to model and interpret these forces using phase-resolved velocity fields from particle image velocimetry measurements. The velocity data are decomposed in a wing-fixed frame to isolate dominant flow features, and pressure fields are reconstructed by solving the Poisson equation for incompressible flows. The relationship between velocity and pressure modes is captured through QSE, which accounts for nonlinear interactions and higher-order dynamics. We introduce an uncertainty-based convergence criterion to ensure model robustness. Applied to a flapping airfoil, the method predicts normal and axial forces with less than 6 % average error using only two velocity modes. The resulting model reveals an interpretable underlying mechanism: linear terms in the QSE model the circulatory force linked to the formation of vortices on the wing, while quadratic terms capture the nonlinear component due to added-mass effects and flow–vorticity interactions. This data-driven yet physically grounded approach offers a compact tool for modelling the unsteady aerodynamics in flapping systems with potential to generalise to other problems.

1. Introduction

1.1. Background and motivation

Natural flyers, and swimmers, rely on flapping wings, or fins, to generate large unsteady fluid forces, enabling them to achieve efficient propulsion and high manoeuvrability even at low Reynolds numbers (Lighthill 1969; Ellington 1984). Inspired by these capabilities, researchers have increasingly focused on unsteady aerodynamics to design bio-inspired micro air vehicles (MAVs) that operate in the low-Reynolds-number regime, ( $Re \leq O(10^4)$ , Shyy et al. 2007). Scientific literature on the topic agrees that the increased performance of flapping wings is due to the formation of high vorticity structures on the leading and trailing edge of the wing (McCroskey 1982). The flapping motion delays stall on the wing, increasing the chord-normal force. The projection of this force over the aerodynamic axes would, in turn, contribute to the increased lift and thrust on the wing, depending on the specific kinematics.

Traditionally, flapping-wing aerodynamics has been modelled using linear potential-flow theory in the high-Reynolds-number limit, assuming small oscillations of a two-dimensional (2-D) airfoil. The pioneering works of Theodorsen (1935) and Garrick (1937) fall within this category. Theodorsen’s lift model includes both an added-mass and a circulatory contribution incorporating both a quasi-steady and wake-induced effects. Garrick’s model extends it for thrust prediction by projecting lift and accounting for leading-edge suction due to the formation of a leading-edge vortex (LEV). While elegant and analytically tractable, these models are limited to small-amplitude oscillations and inviscid high-Reynolds flows, limiting their applicability to MAVs. This mismatch between theory and experimentally observed flow physics was shown by Pitt Ford & Babinsky (2013), illustrating that classical circulation-based models fail to predict the lift enhancement caused by LEVs in highly unsteady, separated flows.

To overcome these limitations and capture the inherently 3-D and nonlinear nature of flapping flight, more recent models have built upon blade element theory (see Branlard 2017). For example, Sane & Dickinson (2002) decompose the chord-normal force into an added-mass force, a quasi-steady circulatory force, a force depending on the wing rotation and a force due to the wake capture. Pesavento & Wang (2004) further refined this model by adding a viscous force term. However, the circulatory components in both models required empirical fitting based on experimental (Sane & Dickinson 2002) or numerical (Pesavento & Wang 2004) data to account for dynamic stall effects.

In contrast to these physics-based approaches, Brunton, Rowley & Williams (2013) derived a purely data-driven force model for real-time force estimation and control purposes. They used balanced proper orthogonal decomposition to derive a reduced-order surrogate model from aerodynamic force data, using the potential linear models as a skeleton. While valuable in control-oriented applications, these models are inherently system specific, and do not offer new insight into the mechanisms of aerodynamic force generation. As shown also by Pitt Ford & Babinsky (2013), the force is highly dependent on the specific flow structures generated on the wing during the flapping motion, and thus this relation should be exploited if the final objective is to provide interpretable force model.

This motivates a broader question: Is it possible to construct a data-driven model that connects flow structures to aerodynamic forces in a way that is both accurate and physically interpretable? Rather than focusing solely on force estimation, the present work addresses this question by aiming to extract dominant flow structures and relate them to aerodynamic forces through interpretable physical mechanisms.

1.2. Flow-based data-driven approaches to unsteady aerodynamics

Flow-based decompositions can provide a clearer link between the generation of aerodynamic forces and the velocity-field structures using a data-driven approach, thus unveiling the physics of the phenomenon. Proper orthogonal decomposition (POD, Berkooz, Holmes & Lumley 1993) offers a basis for reduced-order modelling that can isolate dominant structures in complex flows. For instance, Liang & Dong (2015) applied POD to the flow over a 3-D flapping wing, using a POD-Galerkin projection of the impulse equation (see Noca, Shiels & Jeon 1999) to determine the modal contribution to the force. However, because the decomposition was performed in an inertial frame, in which the wing is moving, it primarily captured the wake dynamics, not allowing to identify the effect of the flow features generated on the wing’s surface. A more accurate decomposition of the flow over the wing requires the flow to be mapped in a reference frame in which the wing is fixed. The choice of the decomposition reference frame has been addressed by several authors, e.g. Lewin & Haj-Hariri (2005) propose performing POD in a wing-fixed reference frame while Troshin & Seifert (2018) propose remapping the flow domain in a steady domain using a volume-preserving transformation.

While these decompositions can provide a compact representation of the flow behaviour over the wing, a further level of complexity is introduced when their relation to the force is sought, especially when this relation must be identified by solving partial differential equations in non-inertial or transformed domains. To address this, Raiola, Discetti & Ianiro (2021) performed POD in a wing-fixed reference frame and applied linear stochastic estimation (LSE) to relate force measurements to the velocity-field modes. This approach, referred to in here as LSE-POD or as extended POD in the literature (see Borée 2003), allowed modelling the unsteady aerodynamic force, linking it to a limited set of velocity modes representative of flow features around the wing. While LSE-POD provides an interpretable link between force and flow structures, the relationship uncovered is purely correlational and lacks grounding in the governing physics, the Navier–Stokes equations. This highlights a key limitation of existing data-driven approaches and motivates the development of a model that remains interpretable by more faithfully adhering to the underlying physics.

The relationship between kinematic quantities (e.g. the velocity field) and dynamic quantities (e.g. the pressure field or the force) has been explored, for other fluid mechanics problems, by means of high-order stochastic estimation methods. For example, Naguib, Wark & Juckenhöfel (2001) estimated the relation between hot-wire measurements and wall-mounted microphones in a boundary layer, finding it necessary to include a quadratic term in the estimation in order to provide a proper modelling. Similarly, Ausseur et al. (2006) used a time-domain quadratic stochastic estimation to estimate the POD time coefficients of the measured velocity field in an unsampled instant in the separated flow over an airfoil, while Baars & Tinney (2014) developed a POD-based spectral higher-order stochastic estimation and tested it by relating a simulated supersonic-jet pressure wavepacket to its acoustic pressure signature in the far field. These studies suggest the question: Is a nonlinear (e.g. quadratic) relation between velocity and dynamic quantities better suited to capturing the physics dictated by the Navier–Stokes equations and to improving the interpretability of the results?

In response to this question, nonlinear modelling strategies have been proposed. Notable examples with applications in fluid dynamics are the sparse identification of nonlinear dynamical systems (SINDy Loiseau, Noack & Brunton 2018), neural networks for order reduction and dynamic prediction (Linot & Graham 2020; Fukami et al. 2021) or Gaussian-process-based system identification (Raissi & Karniadakis 2018). While these approaches have achieved compact and accurate reconstructions of the flow dynamics, they prioritise predictive power over physical interpretability, and typically incorporate governing equations weakly, if at all (see Amico, Matteucci & Cafiero 2025 for a review on different strategies to enhance interpretability in Artificial Intelligence-based fluid models). Exceptions are physic-informed neural networks (Cai et al. 2021), although these embed physical constraints within black-box architectures, offering limited insight into the underlying flow mechanisms. Furthermore, many of these approaches rely on latent embeddings whose relation to measurable physical quantities is unclear, complicating both validation and control-oriented applications.

Beyond these limitations, a more fundamental difficulty emerges: constructing interpretable, low-order force models grounded in the physics of the flow requires the choice of an appropriate low-dimensional representation. If the dimensionality is too low the model may omit essential flow features; if too high the model becomes unnecessarily complex, prone to overfitting and computationally demanding. In linear frameworks based on POD, truncation is typically based on energy content or cumulative variance criteria. On the other side of the spectrum, fully nonlinear model-free approaches either embed this reduced basis selection implicitly, through regularisation (Raissi & Karniadakis 2018) and sparsity promoting architectures (Brunton, Proctor & Kutz 2016) or select it by maximising downstream performance. When nonlinear bases are built these heuristics are no longer applicable, creating the challenge of selecting a modal basis expressive enough to support nonlinear coupling while maintaining interpretability and robustness.

Alternatively to these data-driven approaches, the connection between flow features and aerodynamic forces can be retrieved through a rigorous and interpretable physical modelling. A classical example is force-element analysis (Chang 1992), which decomposes pressure forces into different contributions associated with flow, derived directly from the Navier–Stokes equations. Despite its interpretability, this framework relies on a weak formulation that requires both the actual flow and a potential-flow solution satisfying specific boundary conditions, a requirement that can become restrictive in realistic flow scenarios. Our objective is to pursue an intermediate path: a data-driven decomposition that retains the interpretability of physically derived models, while avoiding their stringent requirements, and that remains mindful of the critical issue of reduced-order representation.

1.3. Objectives and scope of the present contribution

We hypothesise that a quadratic, data-driven model grounded in the physics of the flow, specifically the quadratic stochastic estimation (QSE) of pressure modes from velocity POD modes, can recover pressure and force estimates that are both accurate and physically interpretable. This model should enable the identification of dominant flow structures responsible for aerodynamic force generation, while maintaining consistency with the governing equations and allowing for modal truncation without compromising physical insight. Building on top of the work by Raiola et al. (2021), we test this hypothesis by applying POD to velocity fields measured via planar particle image velocimetry (PIV) in the flow around a flapping airfoil. The decomposition is performed in a wing-fixed reference frame to capture the flow structures forming near the wing’s surface. The pressure field around the wing is estimated by solving the Poisson equation for pressure in incompressible flows, using the measured velocity field. Recognising that the Galerkin projection of the Poisson equation leads to quadratic modal interactions, we use QSE to model the quadratic contribution of the velocity field on the pressure field. The QSE-POD framework is complemented by an uncertainty quantification to assess convergence and robustness. Finally, the pressure distribution is integrated over the wing’s surface to quantify the modal contribution to the aerodynamic force. This approach enables the extraction of physically interpretable flow structures and their connection to force generation via nonlinear interactions. A central question we investigate is whether QSE-POD, despite being a data-driven framework, can yield pressure modes that are compatible with those obtained via a Galerkin projection of the Poisson equation. As we will demonstrate later, the estimated pressure modes satisfy this equation on a mode-by-mode basis, reinforcing the physical credibility of the model.

The remainder of this manuscript is structured as follows. In § 2, we describe the methodology used in this study. This includes the experimental set-up and velocity-field acquisition (§ 2.1), the pressure estimation via the Poisson equation integration (§ 2.2), the flow field decomposition using POD (§ 2.3), the formulation of the QSE-POD method (§ 2.4) and the uncertainty quantification approach used to assess its convergence (§ 2.5). In § 3, we present and discuss the results: we evaluate the effect of mode truncation on the model convergence (§ 3.1), interpret the dominant flow and pressure structures linked to force generation (§§ 3.2 and 3.3) and compare QSE-POD with the LSE-POD approach (§ 3.4). Finally, conclusions are drawn in § 4, highlighting the main contributions of this work.

2. Methodology

2.1. Experimental set-up and velocity-field acquisition

The experimental set-up used in this work is the same as that used by Raiola et al. (2021) and it will be briefly reported here. The flow field around a flapping airfoil in sinusoidal heaving/pitching motion has been measured in the water tunnel facility of the Universidad Carlos III de Madrid, featuring a $2.5\,\text{m} \times 0.5\,\text{m} \times 0.55\,\text{m}$ glass test section. An aluminium 2-D wing with a NACA 0012 section and aspect ratio of 16.3 has been tested at a chord-based Reynolds number of 3600 (within the range of interest of MAVs). Motion has been provided through a 4-bar linkage driven by 2 linear actuators, providing both heaving $h$ and pitching $\theta$ motion as described in (2.1), in which $c$ is the wing chord and $f$ is the flapping frequency corresponding to a Strouhal number $St=2cf/U_\infty =0.2$ , where $U_{\infty}$ is the freestream flow speed,  which has been shown to be a typical value for natural flyers over a range of size and flight speed (Shyy et al. 2007). The analysed mean and amplitude values of the pitch, $\theta _m$ and $\theta _0$ , respectively, are collected in table 1. Wing kinematic parameters were chosen as values expected to maximise the propulsive efficiency, as reported in the publications where measurements have been originally reported (Moriche et al. 2020; Raiola et al. 2021)

(2.1a) \begin{align} &h(t)=c \,\sin (2\pi f t), \end{align}

(2.1b) \begin{align} &\theta (t)=\theta _m+\theta _0\,\sin (2\pi f t+\pi /2). \end{align}

Table 1. Pitching characteristics for each experimental case.

Planar PIV measurements are carried out in a x–y plane placed at a distance of 5.5c from the tip. The flow is seeded with neutrally buoyant polyamide particles with $56\,\mu \text{m}$ diameter. The PIV system is composed by a dual-cavity Nd:Yag Quantel Evergreen laser ( $200\,\text{mJ pulse}^{-1}$ at $15\,\text{Hz}$ ) and a 5.5 M pixel Andor sCMOS camera. The PIV images resolution is approximately $8.5\,\text{pixels mm}^{-1}$ for the present experiment. The PIV system is synchronised with the movement system, providing phase-averaged measurements of 55 snapshots each for 80 different phases of the flapping motion. The PIV image background is removed via an eigenbackground removal procedure (Mendez et al. 2017). Flow fields are extracted through an iterative multi-grid/multi-pass image deformation algorithm (Scarano 2001) with final interrogation windows of $32\times 32$ pixels and $75\,\%$ overlap. Following the typical figure of merit for PIV uncertainty provided by Westerweel (1997), the measurement error is estimated to be 0.8 % of the free-stream velocity. Two separate experiments with illumination on different sides of the wing have been carried out and their phase-averaged fields have been merged in order to provide a complete map of the velocity field around the wing. Some velocity-field snapshots are represented in figure 1, illustrating the baseline flow considered in this study at different phases of the flapping motion in terms of $t/T$ , where $ T = 1/\!f$ is the flapping period. The black arrows represent the velocity field obtained from PIV and colours represent the vorticity in the plane-normal direction $\omega _z = (\boldsymbol{\nabla }\times \boldsymbol{u}) \,\boldsymbol{\cdot }\, \boldsymbol{k}$ , where $\boldsymbol{u}$ is the velocity field (as defined below) and $\boldsymbol{k}$ is the plane-normal unitary vector.

The aerodynamic force is measured with an ATI Industrial Automations Nano-17 load cell with IP68 waterproofing connected between the wing and the linkage. According to the manufacturer’s specification and the test conditions, the resolution of the measured force coefficients is estimated to be $0.03$ . Signal detrending is applied to remove the load cell long-period thermal drift and a low-pass filter with threshold at $St=5$ is applied to remove electromagnetic noise. Force measurements are phase averaged and corrected for the wing’s inertia force and for the waterproofing-induced height-dependent bias, both determined experimentally. The final uncertainty on the force coefficients is estimated experimentally to be lower than $0.1$ .

Figure 1. (a) Case with $\theta _m=0^\circ$ and $\theta _0=10^\circ$ . (b) Case with $\theta _m=10^\circ$ and $\theta _0=10^\circ$ . Baseline PIV flow data time series for different cases in the wing-fixed reference frame. The black arrows represent the velocity field and colours show the plane-normal vorticity field.

Figure 2. (a) Case with $\theta _m=0^\circ$ and $\theta _0=10^\circ$ . (b) Case with $\theta _m=10^\circ$ and $\theta _0=10^\circ$ . Normal (left) and axial (right) force coefficients as estimated from the Poisson equation (red) and as measured from the load cell (black).

2.2. Pressure estimation via the Poisson equation

Pressure fields are estimated from the measured velocity fields by means of the Poisson equation for pressure (see (2.2), Van Oudheusden 2013), where $p$ is the pressure field, $\boldsymbol{u}$ is the velocity field and $\rho$ is the fluid density

(2.2) \begin{equation} {\nabla}^2{\!p}=-\rho \boldsymbol{\nabla \,\boldsymbol{\cdot }\,}(\boldsymbol{u}\boldsymbol{\cdot \boldsymbol{\nabla }})\boldsymbol{u}. \end{equation}

A Matlab-based finite-difference scheme has been implemented, based on the work of Reimer & Cheviakov (2013), supporting mixed Dirichlet- and Neumann-type boundary conditions. Dirichlet conditions with pressure equivalent to null have been imposed on the upwind, upper and lower boundaries of the domain. Non-homogeneous Neumann conditions have been applied on the downwind boundary as well as on the wing’s surface. The Neumann conditions have been derived directly from the Navier–Stokes momentum (see (2.3)), in which $\mu$ stands for the fluid viscosity

(2.3) \begin{equation} \boldsymbol{\nabla \!}p=-\rho \frac {\partial \boldsymbol{u}}{\partial t} -\rho (\boldsymbol{u\,\boldsymbol{\cdot }\,\boldsymbol{\nabla }})\boldsymbol{u}+\mu {\nabla} ^2\boldsymbol{u}. \end{equation}

To reduce discretisation errors in the pressure reconstruction, the original PIV mesh has been interpolated on a mesh with double the resolution. The increased resolution allows for more accurate estimation of wall-normal velocity gradients, which influence the specification of boundary conditions and the accuracy of the Poisson solver. During the interpolation the flow has been assumed to move with the wing speed at the wing’s surface to improve the near-wall estimate. It is important to note that, with a known velocity field, the solution to the non-homogeneous Poisson equation, (2.2), is a Green’s function integral. This implies that pressure will be inherently smoother than velocity, but also that small-scale contributions of the source term will have limited impacted on pressure.

The sectional loads on the wing have been computed by interpolating the pressure field over the airfoil’s surface and then computing the line integral according to the local normal direction of the surface (see (2.4)), i.e.

(2.4) \begin{equation} \boldsymbol{F}_{\textit{Poisson}}(t) = -\oint _{{S}} {p}(\boldsymbol{\boldsymbol{x}},{t}) \boldsymbol{n} (\boldsymbol{{x}}){\rm d}\boldsymbol{{x}}, \end{equation}

where $S$ stands for the airfoil’s surface. The results have been compared with the measured loads in order to validate the pressure integration procedure. The results of the comparison are reported in figure 2 for the cases with $\theta _m=0^\circ$ and $\theta _0=10^\circ$ and $\theta _m=10^\circ$ and $\theta _0=10^\circ$ , showing that the pressure reconstruction captures the dominant pressure structures, while necessarily lacking the full fidelity of the experimental measurements. To quantify the agreement, root mean square error (RMSE) and coefficient of determination $(R^2)$ are computed over the full flapping cycle. For the chord-normal component, we find $R^2 = 0.8246$ and $\textit{RMSE} = 0.2232$ N m⁻¹, while for the axial force, $R^2 = 0.3638$ and $\textit{RMSE} = 0.5718$ N m⁻¹ on average across all the cases. These discrepancies stem primarily from 3-D effects, as for the solution of Poisson equation we are assuming 2-D flow while experiments measure the whole forces experienced by the wing. These observations are supported by the study of Moriche et al. (2020) where 2-D Direct Numerical Simulation flapping-wing simulations were compared with this study’s experimental configuration, displaying uncertainties similar in shape and magnitude. Similar results are obtained for all the experimental cases, but are not reported here for briefness sake.

2.3. Proper orthogonal decomposition of the flow field

Proper orthogonal decomposition (Berkooz et al. 1993) is a mathematical procedure which identifies a set of orthonormal eigenfunctions providing the most compact representation, in a least-square sense, of a set of observations. For the purposes of this work, let the observations correspond to the phase-averaged velocity fields measured with PIV at different flapping phases, which leads to the bi-orthogonal POD definition (Aubry 1991). The velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ is approximated as

(2.5) \begin{equation} \boldsymbol{u}(\boldsymbol{x},t)=\sum _{i=0}^{n_m} \psi ^{(i)} (t) \sigma ^{(i)} \boldsymbol{{\phi }}^{(i)}(\boldsymbol{x}), \end{equation}

with $\boldsymbol{x}$ the spatial coordinate and $t$ representing the flapping phase (non-dimensionalised time over the flapping cycle). The set of vectorial functions $\boldsymbol{\phi }^{(i)}$ constitute the basis of the spatial decomposition of the fluctuating velocity field, $\psi ^{(i)}$ constitutes the basis of the temporal decomposition and $\sigma ^{(i)}$ are the magnitudes associated with each mode; the symbol $n_m$ indicates the number of modes we account for in the approximation.

The functions $\boldsymbol{\phi }^{(i)}, \ \psi ^{(i)}$ and magnitudes $\sigma ^{(i)}$ are computed through the snapshot POD method (Sirovich 1987). Briefly, this method finds a set of orthogonal temporal functions that optimises the decomposition of the two-point temporal auto-correlation function of the flow

(2.6) \begin{equation} C(t,t^{\prime}) = \int [{u}(\boldsymbol{x},t){u}(\boldsymbol{x},t^{\prime})+ {v}(\boldsymbol{x},t){v}(\boldsymbol{x},t^{\prime})]\, {\rm d}\boldsymbol{x}, \end{equation}

where $u$ and $v$ are the components in the $x$ and $y$ directions, respectively, of the velocity field $\boldsymbol{u}$ . The solution to this problem is given by the Fredholm equation, i.e. the solutions of the problem are the eigenfunctions of the temporal correlation function

(2.7) \begin{equation} \int C(t,t^{\prime})\psi ^{(i)}(t^{\prime})\,{\rm d}t^{\prime} = \sigma ^{{(i)}^2}\psi ^{(i)}(t). \end{equation}

It can be proven that the spatial modes, $\boldsymbol{\phi }^{(i)}$ , instead are the eigenfunctions of the temporally averaged spatial correlation function. To obtain the complete vectorial basis $\boldsymbol{\phi }^{(i)} = [\phi _u^{(i)}, \phi _v^{(i)}]$ , the velocity field is projected onto the temporal eigenfunctions according to $\phi _u^{(i)} = \sigma ^{{(i)}^{-1}}\langle u,\psi ^{(i)}\rangle _T$ and $\phi _v^{(i)} = \sigma ^{{(i)}^{-1}}\langle v,\psi ^{(i)}\rangle _T$ , where $\langle \,\,,\,\rangle _T$ indicates the ${L}^2$ inner product in the temporal domain.

In the present work, it will be assumed that the time-averaged velocity field corresponds to $\psi ^{(0)}\sigma ^{(0)} \boldsymbol{\phi }^{(0)}$ by construction, with $\psi ^{(0)}$ being constant in time. In practice, POD will be applied in the wing-fixed reference frame, as was done for the LSE-POD model (Raiola et al. 2021), to obtain velocity structures in the near-wing region. In order to restrict the POD to the fluid domain, the values inside the solid wing have been set to zero prior to performing the POD. As a result, the inner region contributes no kinetic energy even if it is not explicitly excluded from the computation of spatial norms and inner products. This approach, commonly used in PIV-based analyses, allows us to compute a decomposition that is solely driven by spatial correlations in the surrounding fluid region.

2.4. Quadratic stochastic estimation for pressure modal decomposition

The relationship of the velocity modes to the pressure field can be determined from the Poisson equation (see (2.2)). Performing a Galerkin projection of such an equation onto the space spanned by the velocity-based POD modes (see (2.8)), it can be observed that the pressure depends quadratically from the velocity POD modes

(2.8a) \begin{align} &{\nabla} ^2\!p(\boldsymbol{x},t)=-\rho \sum _{i=0}^{n_m}\sum _{j=0}^{n_m}\psi ^{(i)}(t)\psi ^{(j)}(t)\sigma ^{(i)}\sigma ^{(j)} Q^{(i,j)}(\boldsymbol{x}), \end{align}

(2.8b) \begin{align} &Q^{(i,j)}(\boldsymbol{x})=\boldsymbol{\nabla \,\boldsymbol{\cdot }\,}(\boldsymbol{\phi ^{(i)}}(\boldsymbol{x})\boldsymbol{\cdot \boldsymbol{\nabla }})\boldsymbol{\phi ^{(j)}}(\boldsymbol{x}). \end{align}

It is worth noting that in (2.8) the terms having either $i=0$ or $j=0$ are necessarily linear in nature, since one of the temporal functions involved is constant in time (as $\psi ^{(0)}$ stands for the time average). While the POD-Galerkin projection can be used to assess the quadratic contribution of the velocity modes onto the pressure, this approach is complicated by the necessity of restating the boundary conditions in a matching quadratic form (see Noack, Papas & Monkewitz 2005). More so, the Poisson equation should be solved in the non-inertial wing-fixed reference frame, complicating even more the solution to the problem.

To overcome these complications, the physical relation can be captured using QSE. Following the work by Ausseur et al. (2006), let us consider the solution to (2.8) that takes the form

(2.9) \begin{equation} p(\boldsymbol{x},t)=\sum _{i=0}^{n_m}\sum _{j=0}^{n_m}\psi ^{(i)}(t)\psi ^{(j)}(t) \sigma _p^{(i,j)} P^{(i,j)}(\boldsymbol{x}), \end{equation}

with $\sigma _p^{(i,j)}$ and $P^{(i,j)}$ being, respectively, the magnitude and the spatial distribution of the pressure field contribution correlated to the velocity-field mode pair $(i,j)$ . Note that (2.9) does not represent a POD of pressure: neither the temporal modes $\psi ^{(i)}\psi ^{(j)}$ nor the spatial modes $P^{(i,j)}$ are necessarily orthogonal and there is no energy ranking imposed on $\sigma _p^{(i,j)}$ . This quadratic model finds a direct relation between velocity-field modes and their contribution to pressure. The QSE computes the $P^{(i,j)}$ mode by minimising the ${L}^2$ error, i.e. the mean square error for a discrete database. The solution to this problem is provided by the solution of the matrix problem

(2.10) \begin{equation} [\langle \psi \psi ,\psi \psi \rangle ][\sigma P] = [\langle p,\psi \psi \rangle ], \end{equation}

in which $P^{(i,j)}$ is the term to be solved, $\langle \,\,,\,\rangle _T$ indicates the ${L}^2$ inner product in the temporal domain and where the definition of the matrices of the problem is given by

(2.11a) \begin{align} [\langle \psi \psi ,\psi \psi \rangle ] &= \begin{bmatrix} \langle \psi ^{(1)}\psi ^{(1)},\psi ^{(1)}\psi ^{(1)}\rangle & \ldots & \langle \psi ^{(n_m)}\psi ^{(n_m)},\psi ^{(1)}\psi ^{(1)}\rangle \\ \vdots & \ddots & \vdots \\ \langle \psi ^{(n_m)}\psi ^{(n_m)},\psi ^{(1)}\psi ^{(1)}\rangle & \ldots & \langle \psi ^{(n_m)}\psi ^{(n_m)},\psi ^{(n_m)}\psi ^{(n_m)}\rangle \end{bmatrix}\!, \end{align}

(2.11b) \begin{align} [\sigma P] & = \begin{bmatrix} \sigma _p^{(1,1)}P^{(1,1)}\\ \vdots \\ \sigma _p ^{(n_m,n_m)}P^{(n_m,n_m)} \end{bmatrix}\!, \end{align}

(2.11c) \begin{align} [\langle p,\psi \psi \rangle ] &= \begin{bmatrix} \langle p,\psi ^{(1)}\psi ^{(1)}\rangle \\ \vdots \\ \langle p,\psi ^{(n_m)}\psi ^{(n_m)}\rangle \end{bmatrix}\!. \end{align}

As a consequence of this approach, there is no guarantee that the resulting modes $(P)$ will be orthogonal and as such they may present redundant pressure features. In the present work the pressure contributions $P^{(i,j)}$ and $P^{(j,i)}$ , correlated to the mode pairs $(i,j)$ and $(j,i)$ , respectively, have been summed up and considered once, i.e. the contributions with $i\lt j$ have been set to zero. The solution to (2.10) provides a pressure decomposition once the number of modes $n_m$ to be retained is selected. This pressure decomposition can be used to find the contribution each mode has to aerodynamic force generation according to

(2.12a) \begin{align} \boldsymbol{\sigma _{F_{\textit{QSE}}}}^{(i,j)} &= -\oint _{S}\sigma _p^{(i,j)}P^{(i,j)}(\boldsymbol{x})\boldsymbol{n}(\boldsymbol{x})\,{\rm d}\boldsymbol{x}, \end{align}

(2.12b) \begin{align} \boldsymbol{F}_{\!\textit{QSE}} (t) &= \sum _{i=0}^{n_m}\sum _{j=0}^{n_m} \boldsymbol{\sigma _{F_{\textit{QSE}}}}^{(i,j)}\psi ^{(i)}(t)\psi ^{(j)}(t), \end{align}

which, after non-dimensionalising, can be expressed as a force coefficient $\boldsymbol{C}_{F_{\textit{QSE}}} = \boldsymbol{F}_{\!\textit{QSE}}/\rho U_\infty S$ . Note that, given the stochastic nature of QSE, the obtained results will depend heavily on the number of modes considered, it being possible to obtain non-converged modes if $n_m$ is not selected appropriately for our data quality.

2.5. Uncertainty quantification

Uncertainty in the QSE-POD framework arises from both measurement noise (aleatoric) and modelling assumptions (epistemic). Pressure estimation relies on phase-averaged velocity fields and the numerical solution of the Poisson equation, both of which contribute to aleatoric and epistemic uncertainty, respectively. Importantly, these uncertainty sources are not separable in the final model: because QSE estimates pressure from second-order velocity correlations via regression, both types of uncertainty become embedded in the estimated modal coefficients.

Note that the epistemic uncertainty associated with the effect of modal truncation to QSE cannot be evaluated analytically due to the regressor nature of the QSE model. To address this, we adopt an empirical approach based on robustness to noise: we study how the sensitivity of the model to added noise varies with the number of retained modes $n_m$ . While this does not quantify uncertainty in the physical sense, it evaluates the statistical convergence of our estimator, revealing truncation orders beyond which the model becomes unstable or prone to overfitting.

Two metrics can be built to assess the performance of the model. First, for a general measured magnitude $q$ and its prediction $q_{\textit{QSE}}$ , it is possible to define the residual of the prediction $r_q = \|q-q_{\textit{QSE}}\|/\|q\|$ , where $\|\,\boldsymbol{\cdot }\,\|$ refers to the $L^2$ norm in the whole spatio-temporal domain of $q$ . Analogously, the signal level of the reconstruction can be defined as $S_q = 1-r_q$ . Secondly, it is possible to define the noisiness of said prediction $N_q = \|\Delta q_{\textit{QSE}}\|/\|q\|$ , where $\Delta q_{\textit{QSE}}$ is the uncertainty map of $q_{\textit{QSE}}$ . A QSE-POD model capable of making accurate predictions will report signals close to 1, while a precise method will yield reduced noise values. For the sake of interpretability, we will consider a model convergent if the difference between the pressure signal and noise is close to one, $S_p-N_p \approx 1$ .

Simultaneously, this analysis takes into account an important characteristic of the present dataset: only a single realisation is available for each flapping parameter combination. This prevents the use of conventional train/test data splits to evaluate the generalisation of the model. To address this, our truncation strategy is designed not only to capture dominant flow features, but also to act as a form of regularisation, limiting the inclusion of potentially redundant or collinear terms in the quadratic basis. By constraining the modal order, we aim to mitigate overfitting and ensure robust reconstructions, even in the absence of independent validation data.

To assess the model’s reliability in practice, it is necessary to evaluate the uncertainty associated with the QSE-POD predictions $\Delta q_{\textit{QSE}}$ . This is achieved through an uncertainty quantification procedure that combines two complementary techniques: statistical uncertainty propagation, based on local linearisation, and Monte Carlo simulation, based on stochastic sampling.

Statistical uncertainty propagation (Coleman & Steele 2009) is based on the first-order truncation of a Taylor expansion. As shown in (2.13), the uncertainty of $\boldsymbol{y}$ , denoted as $\boldsymbol{\Delta y}$ , can be approximated from the Jacobian $\boldsymbol{\nabla \!}\boldsymbol{f}(\boldsymbol{x})$ and the uncertainty in $\boldsymbol{x}$ , i.e. $\boldsymbol{\Delta x}$

(2.13a) \begin{align} \boldsymbol{y} + \boldsymbol{\Delta y} = \boldsymbol{f}(\boldsymbol{x} + \boldsymbol{\Delta x}) &= \boldsymbol{f}(\boldsymbol{x}) + |\boldsymbol{\nabla \!}\boldsymbol{f}(\boldsymbol{x})|\boldsymbol{\Delta x} + O(||\boldsymbol{\Delta x}||_2), \end{align}

(2.13b) \begin{align} \boldsymbol{\Delta y} &\approx |\boldsymbol{\nabla \!}\boldsymbol{f}(\boldsymbol{x})|\boldsymbol{\Delta x}. \end{align}

Which is to say, the uncertainty of the output is the uncertainty of the input times the sensitivity of the output to the input. This has been used to estimate the singular value uncertainty $\Delta \sigma , \ \Delta \sigma _p$ computations, which act as weights in the projection of the velocity field $\boldsymbol{u}$ onto the temporal modes $\psi$ . Since these singular values are essentially Euclidean norms in the discretised domain, their uncertainty can be propagated analytically as in (2.14)

(2.14a) \begin{align} \Delta \sigma ^{(j)} &= \frac {2}{\sigma ^{(j)} }\sum _k^{N_p}\phi ^{(j)}(x_k)\Delta \phi ^{(j)}(x_k), \end{align}

(2.14b) \begin{align} \Delta \sigma _p^{(i,j)} &= \frac {2}{\sigma _p^{(i,j)} }\sum _k^{N_p}P^{(i,j)} (x_k)\Delta P^{(i,j)} (x_k). \end{align}

When an analytical Jacobian is not available, we apply Monte Carlo simulation. As described by Papadopoulos & Yeung (2001), Monte Carlo is a stochastic method that estimates uncertainty by sampling from the input space, converging towards uncertainty propagation results. A set of $M =1000$ samples is drawn from a uniform distribution ( $\mathcal{U}$ ) within the measurement bounds, i.e.

(2.15) \begin{equation} X = \{\boldsymbol{x}_1,\boldsymbol{x}_2,\ldots ,\boldsymbol{x}_M\}\sim \mathcal{U}_{[\bar {x}-\Delta x, \bar {x}+\Delta x]}. \end{equation}

The assumption of uniformly distributed data within the uncertainty ensures a conservative estimate of the mean and the width of the confidence interval, as the uniform distribution imposes greater variance than the normal distribution, another commonly used alternative. Next, the output set is computed as $F = \{\boldsymbol{f}(\boldsymbol{x}_1),\ldots ,\boldsymbol{f}(\boldsymbol{x}_M)\}$ , from which both the expected value $\mathbb{E}[F]$ and a $95\,\%$ confidence interval are extracted. To check convergence regarding the number of samples $(M=1000)$ , Monte Carlo simulation was repeated three times, for different random number seeds. This experiment resulted in three distinct mean and confidence intervals showing differences below $0.1\,\%$ in force coefficient estimates.

Monte Carlo is particularly useful for mode and force uncertainty estimation, where the Jacobian is either unknown or intractable. The combination of these two techniques enables the full quantification of the QSE-POD predicted force $\boldsymbol{F}_{\!\textit{QSE}}$ and its uncertainty $\boldsymbol{\Delta F}_{\textit{QSE}}$ . A summary of the uncertainty quantification procedure is collected in table 2.

Table 2. Uncertainty quantification summary.

Note that $\Delta p$ is set to the numerical integration error. As $p$ is obtained from the integration of (2.2), the samples obtained from Monte Carlo will have implicit in them some numerical error $p = p +O(\Delta x^2)$ . After performing Monte Carlo simulation, it is observed that the numerical error dominates over measurement uncertainty propagation, i.e. $\Delta x^2\gt \gt \Delta p_\textit{Monte-Carlo}$ . Then, to be conservative in our computations, the final pressure uncertainty is set to numerical integration error $\Delta p =\Delta x^2$ .

3. Results and discussion

3.1. Model convergence and truncation order selection for QSE-POD

As addressed in § 2.4, the QSE model for pressure depends on the number of modes $n_m$ retained in the low-order reconstruction of the flow. It is expected that the greater $n_m$ is, the smaller the residual is, increasing the value of $S_p$ . On the other hand, increasing the number of modes considered is related to an increase in the noisiness of the data $N_p$ . The selection of $n_m$ , therefore, will be performed in order to maximise $S_p-N_p$ . In the following, only the modes whose kinetic energy content is above 5 % of the total ( $\sigma ^{(k)^2}/\sum _{i}\sigma ^{(i)^2}\gt \epsilon =0.05$ ) will be considered for the computation of $S_p-N_p$ . Following this criterion, only low-order reconstructions with $n_m\leq 6$ , i.e. involving modes 0–6 (see figure 3), will be explored for the maximisation. Note that the 5 % threshold has been set arbitrarily and should be expanded to consider more modes if no clear optimum is found within the range of interest. Although this threshold is empirical, the resulting range coincides with the physical interpretation given in prior studies: Raiola et al. (2021) showed that modes 1–2 predominantly capture the periodic formation and shedding of leading- and trailing-edge vortices, while modes 3–6 are associated with smaller-scale vorticity features. In this sense, our selection reflects not only a statistical criterion but also aligns with known features of the flapping flow’s coherent structures and their dominance over unsteady aerodynamic loads.

Figure 3. On the left: Singular values kinetic energy (black) and cumulative kinetic energy content (red) with the selected cutoff value $\epsilon = 0.05$ . On the right: Pressure signal (black) and noisiness (red) for different $n_m$ values; solid lines show the mean value and the shaded region is between the maximum and minimum values for the different experimental cases.

Within this range, $S_p-N_p$ is maximal for $n_m = 2$ on average between all the different experimental cases (see figure 3). This maximum is driven mainly by the noisiness: the pressure noise level is highly sensitive to the cutoff mode number $n_m$ changes, sharply increasing for $n_m\gt 2$ , unlike the pressure signal level which starts at $S_p = 0.75$ for $n_m=1$ and asymptotically increases towards $S_p = 1$ . The high sensitivity of reconstruction noise $N_p$ to $n_m$ makes highly relevant the quantification of uncertainty for any outputs obtained using this framework. This sharp rise in noise can be interpreted as a symptom of overfitting, potentially arising from the emergence of near-collinear elements in the quadratic basis used to reconstruct the pressure field. Such collinearities do not prevent us solving the regression problem; however, they can amplify the model’s sensitivity to noise and degrade its generalisability. This observation motivates the adoption of a conservative truncation threshold that balances physical fidelity with robustness. In the following, the velocity, pressure and force decomposition will be obtained for a value $n_m=2$ and their mutual relation will be analysed.

3.2. Structure of velocity modes

In the following, the results of the velocity-field decomposition are shown for the flapping-wing cases B (with $\theta _m=0^\circ$ and $\theta _0=10^\circ$ ) and E (with $\theta _m=10^\circ$ and $\theta _0=10^\circ$ ). The trends reported for these representative cases are, however, applicable to all experimental conditions that share the same mean pitch setting (i.e. observations for case B extend to cases A and C, while those for case E extend to cases D and F). The POD of the velocity field is reported in figure 4 for the 0th mode (corresponding to the time average), the 1st mode and the 2nd mode. The contribution of these three modes to the full field is reported in table 3, showing that these three modes recover more than 70 % of the kinetic energy of the flow.

Figure 4. Proper orthogonal decomposition of the flow for the cases with $\theta _m = 0^\circ$ and $\theta _0 = 10^\circ$ (left) and $\theta _m = 10^\circ$ and $\theta _0 = 10^\circ$ (right): (1st row) 0th POD spatial mode; (2nd row) 1st POD spatial mode; (3rd row) 2nd POD spatial mode. The spatial modes are represented in terms of vorticity (contour plot) and velocity (quiver plot). The temporal modes in the insets are represented for a complete oscillation $t/T \in [0,1]$ over each spatial mode.

The 0th mode represents a free stream angled by the mean pitch and two counter-rotating vortices on the two sides of the airfoil. The 1st and 2nd modes represent vortical structures aligned along the mean flow direction at different positions. Temporally, they show a quasi-sinusoidal behaviour with period equal to the flapping period and a time shift of $\pi /2$ between them, suggesting that they model the same spatio-temporally coherent feature. When summed together accordingly to the POD decomposition in (2.5), they model the evolution of the quasi-steady bound vorticity on the wing, as well as the one due to LEVs and trailing-edge vortices (TEVs), during the flapping cycle. This vorticity is shown to grow on the wing up to being released in the wake along the mean flow direction. This vorticity can be assimilated to a trail of vortices forming on the wing and being advected with the flow after separation. The evolution of this vorticity aligns with the heaving motion for cases with $\theta _m=0$ , while it shows signatures of delayed stall for $\theta _m=10$ , such as the persistence of attached vorticity structures during the upstroke and a later detachment of the shear layer. The spatial average of the velocity, instead, can be assimilated to the heave velocity of the wing. These results are coherent with the ones observed by Raiola et al. (2021).

Table 3. Kinetic energy content of the velocity-field modes normalised by the total kinetic energy ( ${\sigma ^{(k)^2}}/{{\sum} _i\sigma ^{(i)^2}}$ ) for the cases B and E.

Figure 5. Frequency analysis of POD modes (not including the mean mode). The first (black squares) and second (red circles) mode contributions are reported separately. Data for the cases with $\theta _m = 0^\circ$ , $\theta _0 = 10^\circ$ and $\theta _m = 10^\circ$ , $\theta _0 = 10^\circ$ are presented on the left and rght plot, respectively.

In addition to their spatial interpretation, the temporal dynamics of the POD coefficients can be analysed through its power spectral density (PSD). Figure 5 reports the PSD of the non-mean modes, showing that their energy is concentrated at the fundamental flapping frequency, with weaker peaks at flapping frequency harmonics. This dominant peak is consistent with the periodic growth and shedding of vortical structures from the airfoil. Remarkably, the second mode exhibits a sparser distribution of energy, containing contributions at higher harmonics that are more typical of higher-order modes, especially for the cases with $\theta _m = 10^\circ$ . This indicates that part of the higher-frequency dynamics is embedded within the second POD mode, rather than being isolated in separate modes. Such behaviour is consistent with the nonlinear nature of the flow: while the main vortical features are locked to the flapping frequency, weaker oscillations at higher frequencies naturally arise from nonlinear interactions and are redistributed across the low-order modes by the POD.

3.3. Pressure modal decomposition and force generation

Using the first 3 modes of the velocity field, a QSE is performed to reconstruct the pressure field and identify the pressure structures associated with them. Based on their modal indices, QSE-POD modes can be grouped into three contributions: a constant term $(0,0)$ which corresponds to the mean flow and is time invariant; linear term (i.e. $(0,j)$ modes with $j\neq 0$ ) which mirror the temporal behaviour of the original velocity POD modes 1 and 2 and describe their interaction with the mean flow; quadratic term (i.e. $(i,j)$ modes with $i,j\neq 0$ ) which introduce nonlinear modal interactions between modes 1 and 2 and new dynamics into the model.

To clarify this decomposition, recall that the temporal coefficient $\psi ^{(0)}$ is constant, since it corresponds to the temporal average. Consequently, the mode $(0,0)$ yields a steady pressure field. Additionally, pressure modes involving repeated indices $(i,i)$ generate positive-definite temporal products $\psi ^{(i)} \psi ^{(i)} \geq 0$ , which result in net force contributions that persist over time.

Before proceeding with further analysis, it is worth clarifying the interpretive scope of our discussion. The QSE is a statistical regression method and, as such, it does not guarantee a causal or orthogonal decomposition. All the observations are based on correlations, which do not necessarily imply a direct physical mechanism even if they might suggest the presence of it. Any claim can be speculative without the support of a much bigger set of data which can help generalise the observations. Aware of these limitations, we adopt a physically motivated lens to interpret the QSE outputs. In particular, we follow the spirit of the force-element approach introduced by Chang (1992), in which aerodynamic forces are understood as arising from the spatial distribution of flow structures. In this framework, pressure-induced forces are expressed as

(3.1) \begin{align} \int _Sp(\boldsymbol{\alpha \,\boldsymbol{\cdot}\,n}){\rm d}A& = \rho \int _S\phi \frac {\partial \boldsymbol{u}}{\partial t}\boldsymbol{\cdot\, n} {\rm d}A-\frac {\rho }{2}\int _S|\boldsymbol{u}^2|(\boldsymbol{\alpha \,\boldsymbol{\cdot }\,n}){\rm d}A \nonumber \\ &\quad + \mu \int _S\boldsymbol{n}\times \boldsymbol{\omega \,\boldsymbol{\cdot }\,}\boldsymbol{\nabla }\phi {\rm d}A - \rho \int _V \boldsymbol{u}\times \boldsymbol{\omega \,\boldsymbol{\cdot }\,}\boldsymbol{\nabla }\phi {\rm d}V, \end{align}

where $\boldsymbol{\alpha }$ is the unitary vector indicating the direction of the component of interest of the force, $\phi$ is a general potential function fulfilling that ${\nabla} ^2\phi = 0$ and, at the boundary, $\boldsymbol{\nabla }\phi \boldsymbol{\cdot n} = (\boldsymbol{\alpha \boldsymbol{\cdot } n})$ and $\boldsymbol{\omega } = \boldsymbol{\nabla }\times \boldsymbol{u}$ is the vorticity of the field. The physical interpretation separates between the no-vorticity, or added-mass, forces represented in the first line of the right-hand side of the equation  and the vorticity-related forces which are included in the bottom part of the right-hand side. The present study follows the same spirit as the force-element analysis, seeking to simplify the interpretation of force generation by associating pressure patterns with combinations of dominant velocity features, but introduces a data-driven formulation that infers these relationships directly from the measured velocity modes. In this sense, our analysis aims to provide spatial and temporal support to the physical ideas behind force-element theory, using experimentally derived modal pressure fields to interpret how unsteady aerodynamic forces arise.

In the following, the individual pressure modes derived from QSE-POD are examined, focusing on their spatial structure and associated force contributions. While the modes are not guaranteed to be orthogonal or uniquely linked to isolated flow features, their structure reveals consistent patterns that can be related to dominant aerodynamic mechanisms. By analysing these contributions, we identify which combinations of velocity modes are most influential in shaping the pressure field and generating aerodynamic forces.

3.3.1. Mode-by-mode contributions

The spatial distribution $P^{(i,j)}$ , representing the pressure contribution of each velocity mode pair $(i,j)$ to the pressure, is shown in figures 6 and 7 for the experimental cases B and E, respectively. Each panel includes the associated temporal pattern $\psi ^{(i)}\psi ^{(j)}$ , illustrating the evolution of pressure throughout the oscillation period. The relevance – in terms of their variance – that these patterns have for the pressure field ( $\sigma _p^{(i,j)}$ ) and on the force generated ( $\boldsymbol{\sigma^{(i,j)} _{F_{\textit{QSE}}}}$ ) is summarised in figure 8. The force contribution is represented in terms of chordwise ( $\sigma ^{(i,j)}_{F_x}$ ) and chord-normal ( $\sigma ^{(i,j)}_{F_y}$ ) components .

Figure 6. Spatial distribution $P^{(i,j)}$ of the contribution of the mode pair $(\psi ^{(i)},\psi ^{(j)})$ to pressure for the modes from 0th to 2nd along an oscillation. The temporal modes are represented for $t/T\in [0,1]$ over each spatial mode. The pressure distribution is reported in the colour map. Green and black quiver maps indicate the velocity field contained in the spatial POD modes $\phi ^{(i)}$ and $\phi ^{(j)}$ , respectively. The red vector shows the normalised force generated by the pressure distribution. Represented here is the case with $\theta _m = 0^\circ$ and $\theta _0 = 10^\circ$ .

Figure 7. Spatial distribution $P^{(i,j)}$ of the contribution of the mode pair $(\psi ^{(i)},\psi ^{(j)})$ to pressure for the modes from 0th to 2nd. The pressure distribution is reported in the colour map. Green and black quiver maps indicate the velocity field contained in the spatial POD modes $\phi ^{(i)}$ and $\phi ^{(j)}$ , respectively. The red vector shows the direction of the force generated by the pressure distribution. Represented here is the case with $\theta _m = 10^\circ$ and $\theta _0 = 10^\circ$ .

Figure 8. (a) Case with $\theta _m=0^\circ$ and $\theta _0=10^\circ$ . (b) Case with $\theta _m=10^\circ$ and $\theta _0=10^\circ$ . Squared magnitude of the contribution of the mode pair $(\psi ^{(i)},\psi ^{(j)})$ to: (left) pressure $\sigma _p$ ; (middle) chordwise force $\sigma _{F_x}$ ; (right) chord-normal force $\sigma _{F_y}$ . Results are indicated in terms of percentage of the total contribution of the first $3$ modes.

Interestingly, while all the mode pairs $(i,j)$ contribute significantly to the pressure, only few pairs contribute to the forces. In particular, for the cases with $\theta _m = 0^\circ$ , the $(0,0)$ mode dominates the generation of axial forces, with a less important contribution from quadratic modes. However, in the cases with $\theta _m = 10^\circ$ , this dominance shifts mainly toward quadratic modes, which represent quadratic interactions between the first two non-constant velocity modes, with an important contribution also from linear modes. This shift can be attributed to the misalignment between the mean flow and the chordwise direction, which introduces additional asymmetries affecting the normal force component.

These observations are supported by the spatial pressure distributions: for $\theta _m = 0^\circ$ , the $(0,0)$ mode is almost symmetric with respect to the chord, with strong pressure variations (and thus force) in the chord direction; for $\theta _m = 10^\circ$ this symmetry is lost and the pressure gradient turns mostly in the normal direction. An opposite behaviour occurs for linear modes: for $\theta _m = 0^\circ$ they produce an anti-symmetric pressure gradient with respect to the chord (and thus almost no chordwise force contribution); for $\theta _m = 10^\circ$ this symmetry is lost and part of the pressure gradient is directed along the chordwise direction. Consistently, chordwise asymmetries in the quadratic modes are present for both cases, but become more pronounced when $\theta _m = 10^\circ$ .

In contrast, normal force generation is primarily driven by linear modes. For $\theta _m = 0^\circ$ , the most influential contributor is the $(0,1)$ mode, whose relative impact increases with greater $\theta _0$ . Although this trend persists for $\theta _m = 10^\circ$ , it becomes less dominant mostly in favour of the constant contribution $(0,0)$ , even if a slightly increased contribution from quadratic modes is observed. Again, this behaviour is reflected in figures 6 and 7: for $\theta _m = 0^\circ$ , the $(0,1)$ mode presents a strong pressure gradient aligned in the chord-normal direction, leading to normal force generation. When $\theta _m = 10^\circ$ , this directional pattern becomes less exclusive of the $(0,1)$ mode, as both the constant and quadratic contributions show more important pressure gradients in the normal direction.

To summarise the force contribution by different modes, while quadratic modes play a key role in axial force production, normal forces are largely governed by linear terms. However, when a mean angle is present, symmetries of the modes are broken and contributions start to mix up between force components.

As an empirical verification of the physical link between velocity and pressure modes, we assess how closely the QSE-POD pressure modes satisfy the Poisson equation on a mode-by-mode basis. In practice, the residuals of the approximation ${\nabla} ^2P^{(i,j)}\approx \boldsymbol{\nabla \,\boldsymbol{\cdot }\,} [(\boldsymbol{\phi }^{(i)}\,\boldsymbol{\cdot\, \boldsymbol{\nabla }})\boldsymbol{\phi }^{(j)}]$ are evaluated, yielding mean normalised residuals below $1\,\%$ . Specifically, the maximum mean residual observed across all cases and truncations remains below 0.20 % globally, below 0.25 % in the near-field region ( $x \in [-0.5, 2.5]$ , $y \in [-1.5, 1.5]$ ) and below 0.61 % in the region adjacent to the airfoil, where forces are integrated. Peak local errors can be higher, with maximum values of 37.4 % globally, 30.6 % in the near field and 8.5 % around the airfoil. These higher values are spatially localised and occur primarily in regions far from the wing where non-physical values due to frame rotation appear. Note that these results are only true for the selected low-order representation of the pressure decomposition, which may be accurate in the present dataset but does not imply that QSE-POD pressure modes generally satisfy the Poisson equation. They should therefore be interpreted as an empirical property of the analysed cases, rather than a guaranteed feature of the method.

These results validate QSE-POD as an effective surrogate for recovering the Galerkin-projected solution of the pressure Poisson (see (2.8)), particularly in the regions most relevant to force generation. Directly obtaining this decomposition via Galerkin projection remains significantly more challenging, as the nonlinear nature of the Poisson equation and the presence of moving boundaries in a wing-fixed frame require careful treatment of boundary conditions and mode interactions. In contrast, QSE-POD enables a data-driven decomposition that, when applied with an appropriate truncation order and to a well-posed quadratic problem, can approximate the pressure contributions of kinematic features on a stochastic basis, yielding results that are consistent with physically sound behaviour.

3.3.2. Combined modal contributions

While the pressure patterns reported in figures 6 and 7 have a direct relation to the corresponding velocity-field modes, their behaviour is difficult to comprehend since, as with the velocity-field POD modes, they do not represent, individually, spatio-temporally coherent features. A more clear way to analyse the flow patterns with QSE-POD analysis is by combining different modes to construct spatio-temporally coherent behaviours, as reported in figures 9 and 10. Pressure contributions are separated into linear (2nd row) and quadratic (3rd row) contributions of the spatio-temporal coherent features represented by velocity-field modes 1 and 2. The non-constant fluctuating pressure field obtained from the solution of the Poisson equation (1st row) is reported for reference. The corresponding contribution to the aerodynamic force through time is reported on the top. The quiver plot represents the evolution of the spatio-temporal coherent features contained in modes 1 and 2.

Figure 9. Spatial distribution of the contribution of the linear (middle) and quadratic (bottom) modes to pressure and force for the case with $\theta _m = 0^\circ$ and $\theta _0 = 10^\circ$ . At the top we show the Poisson equation pressure and forces for reference. The body axial and perpendicular forces generated by each pressure distribution are represented in red and blue, respectively. The pressure distribution is reported in the colour map. The quiver map indicates the non-mean velocity field as reconstructed from modes 1 and 2. The mean velocity and pressure are not accounted for.

Figure 10. Spatial distribution of the contribution of the linear (middle) and quadratic (bottom) modes to pressure and force for the case with $\theta _m = 10^\circ$ and $\theta _0 = 10^\circ$ . At the top we show the Poisson equation pressure and forces for reference. The body axial and perpendicular forces generated by each pressure distribution are represented in red and blue, respectively. The pressure distribution is reported in the colour map. The quiver map indicates the non-mean velocity field as reconstructed from modes 1 and 2. The mean velocity and pressure is not accounted for.

From this visualisation we can infer how the shown pressure distribution generates the corresponding forces. The linear modes represent pressure oscillations at the same frequency as the flapping motion. For the cases where $\theta _m = 0^\circ$ a chord-symmetric pressure distribution appears to be shed downstream, inducing mainly chord-normal forces. When $\theta _m = 10^\circ$ , this symmetry is broken. The resulting asymmetry in the vortex propagation direction leads to pressure imbalances along the chord, which cause the force vector to rotate clockwise. This reorientation reflects the influence of asymmetries in unsteady flow features, particularly the formation and shedding of LEVs, which are now shed deflected with respect to the chord line.

The quadratic modes, instead, are associated with pressure fluctuations that exhibit two distinct peaks per flapping cycle, resulting in a more complex temporal structure. In all the cases, they produce pulsating pressure patterns along the mean flow direction, contributing more to the chord-wise forces than their linear counterparts. This effect is accentuated with mean pitch, as can be noticed in cases with $\theta _m = 10^\circ$ , where flow causes these pulsating mean patterns to induce stronger chord-normal forces compared with the $\theta _m = 0^\circ$ case. Moreover, across all configurations, the quadratic modes are associated with net-negative chord-normal forces, whereas positive axial forces emerge primarily in cases with non-zero-mean pitch. This force signatures mimic the formation and shedding of vortices which, in asymmetric unsteady flow conditions, tend to reorient aerodynamic loading, reducing the chord-normal force in exchange for positive axial contributions. The contrasting behaviour between zero and non-zero-mean pitch cases, which emerges in the present dataset, points to the fundamental role of vortex asymmetries and stall dynamics in reorienting the aerodynamic forces.

In zero-mean pitch cases ( $\theta _m = 0^\circ$ ), both components of the linear forces display oscillations at the flapping frequency and roughly in phase with heave velocity (see figure 9). This behaviour is lost in the mean pitch cases (see figure 10), particularly at the onset of the upstroke ( $t/T = 0.75$ ) where pressure symmetry also breaks down. Furthermore, quadratic modes display low-pressure oscillations for cases with $\theta _m = 0^\circ$ , whereas mean pitch cases exhibit pressure bursts, particularly at the initiation of the upstroke. This phenomenon can be interpreted in the light of force-element analysis (see (3.1)); in particular, focusing on the $\boldsymbol{u}\times \boldsymbol{\omega }$ term, associated with the Kutta–Joukowski lift, which can be regarded as the main contributor to the aerodynamic forces. It is possible to expand that term to each modal contribution

(3.2) \begin{equation} \int _V \boldsymbol{u}\times \boldsymbol{\omega \,\boldsymbol{\cdot }\,}\boldsymbol{\nabla }\phi {\rm d}V = \int _V (\boldsymbol{u}_0\times \boldsymbol{\omega }_0 + \boldsymbol{u}_{1+2}\times \boldsymbol{\omega }_0 + \boldsymbol{u}_0\times \boldsymbol{\omega }_{1+2} + \boldsymbol{u}_{1+2}\times \boldsymbol{\omega }_{1+2}) \,\boldsymbol{\cdot }\,\boldsymbol{\nabla }\phi {\rm d}V, \end{equation}

where $\boldsymbol{u}_0$ and $\boldsymbol{\omega }_0$ are the constant velocity and vorticity fields, respectively, and $\boldsymbol{u}_{1+2}$ and $\boldsymbol{\omega }_{1+2}$ are the velocity and vorticity fields associated with the first two POD modes.

As discussed in § 3.2, the vorticity $\boldsymbol{\omega }_0$ contained in the mode 0 can be interpreted as two counter-rotating vortices on the two sides of the airfoil, while the velocity field $\boldsymbol{u}_0$ can be, on average, assimilated to a constant velocity aligned along the average flow direction. Accordingly, the force related to the $\boldsymbol{u}_0\times \boldsymbol{\omega }_0$ term of (3.2), the constant contribution, is aligned with the normal to the average flow direction. While $\theta _m=0^\circ$ , the intensity of these vortical structures is mostly balanced, producing little to no force in the chord-normal direction. When a mean pitch is included, an unbalance exists, producing a force affecting both, chord-normal and chordwise components.

The vorticity $\boldsymbol{\omega }_{1+2}$ contained in the modes 1 and 2 can be associated with a trail of vortices forming on the wing and being advected with the flow after separation. The evolution of this vorticity aligns with the heaving motion for cases with $\theta _m=0^\circ$ peaking at $t/T = 0.5$ , while it is affected by delayed stall for $\theta _m=10^\circ$ . The spatial average of the velocity $\boldsymbol{u}_{1+2}$ close to the airfoil, instead, can be assimilated to the heave velocity of the wing. The $\boldsymbol{u}_{1+2}\times \boldsymbol{\omega }_0$ term of (3.2), therefore, represents the interaction between the heave velocity and the vortical structures of mode 0, since $\boldsymbol{\omega }_0$ retains significant values near the airfoil, where $\boldsymbol{u}_{1+2}$ approaches the heave velocity. The $\boldsymbol{u}_0\times \boldsymbol{\omega }_{1+2}$ term instead, accounts for the interaction between the trail of vortices in modes 1 and 2 with the average velocity. Both terms contribute to QSE linear forces. For zero-mean pitch cases, the $\boldsymbol{u}_{1+2}\times \boldsymbol{\omega }_0$ term would produce nearly no forces due to the orientation of $\boldsymbol{u}_{1+2}$ and the quasi-symmetry of the vorticity $\boldsymbol{\omega }_0$ . The $\boldsymbol{u}_0\times \boldsymbol{\omega }_{1+2}$ term, instead, would produce a force directed along the chord-normal direction. When the symmetry is broken through the introduction of a mean pitch, the second term can be expected to produce a relevant contribution, while the third term would follow the temporal pattern of delayed vortex shedding produced by the deep stall. The combination of these effects leads to the observed shift in phase of the peak of the linear force contributions.

Finally, the $\boldsymbol{u}_{1+2}\times \boldsymbol{\omega }_{1+2}$ term includes the interaction between the oscillating velocity and the vorticity shedding, which would contribute to the quadratic forces. For zero-mean pitch cases, this term is expected to produce forces along the chordwise direction due to vortex generation and shedding close to the airfoil, while the contribution to chord-normal forces comes from the interaction of the wake vortices and the velocity induced by such vortices. For $\theta _m = 10^\circ$ cases, this term is expected to be strongly affected by the deep stall behaviour and by the delayed vortex shedding. In particular, the expected delay in shedding is consistent with the bursts observed during the upstroke, when unsteady circulation is expected to peak over the airfoil.

Figure 11. (a) Case with $\theta _m=0^\circ$ and $\theta _0=10^\circ$ . (b) Case with $\theta _m=10^\circ$ and $\theta _0=10^\circ$ . Frequency analysis of linear and quadratic modes (without including the mean mode). The left plots show the frequency content of normal forces and in the right column the same is shown for the axial forces. The linear mode (black squares) and quadratic mode (red circles) contributions are reported separately.

Crucially, such interpretation can be made here because the QSE framework includes quadratic temporal modes $\psi ^{(i)}\psi ^{(j)}$ providing the decomposition with the ability to represent higher-frequency dynamics. This capability appears particularly relevant in the axial force (see figure 11). In this sense, quadratic modes do not merely complement the linear ones, but may play an essential role in accounting for the broadband fluctuations observed in the forces. In figure 11 the power spectrum of the quadratic and linear forces $|\mathcal{F}[F_i](f)|^2$ is plotted against the harmonics of the flapping frequency $f/f_0$ , where $f_0$ represents the flapping frequency. The spectra are analysed up to the tenth harmonic. For the case with $\theta _m = 0^\circ$ , the frequency content exhibits sharp, well-defined peaks: linear modes concentrate their energy at odd harmonics of the flapping frequency, while quadratic modes concentrate on even harmonics. The resulting peaks are sharp and well separated, reflecting the symmetry and periodicity of the motion. In contrast, when $\theta _m = 10^\circ$ , the spectral content becomes more broadband. Although linear and quadratic modes still peak near $f/f_0 = 1$ and $f/f_0 = 2$ , respectively, the harmonic separation becomes less distinct and the energy spreads across neighbouring frequencies. This peak broadening suggests that the presence of mean pitch introduces asymmetries and flow nonlinearities that break the ideal harmonic structure, leading to more complex, less periodic force signals. Moreover, quadratic modes are less compact in frequency, as their power spectra decay more slowly compared with those of the linear modes. This suggests the emergence of higher-frequency harmonics as quadratic terms contain frequencies that originally did not appear in the linear modes (see also figure 5).

Both behaviours, the emergence of harmonics and the spread in frequency content, are well-documented phenomena in wake-dominated flows (Elgar, Van Atta & Gharib 1990). Quadratic interactions generate triadic couplings among velocity features, through which self-interaction triads of the fundamental shedding frequency produce higher harmonics (Vandiver & Jong 1987). These triadic interactions arise from the nonlinear terms and correspond to energy exchanges among frequency triplets satisfying the resonance condition $f_1 + f_2 = f_3$ (Waleffe 1992). Figure 5 shows a strict dominance of the flapping frequency, the first and the second harmonics. While the linear QSE modes follow this behaviour too, quadratic modes present a much different picture, with the first harmonic dominant in their spectra and presenting higher frequencies. This supports the idea that quadratic modes capture the effects that triadic interactions have on forces, as most of the high-frequency behaviour is contained within these modes, which must arise from the quadratic terms $\boldsymbol{u}\times \boldsymbol{\omega }$ and $|\boldsymbol{u}|^2$ of force-element analysis.

Moreover, triadic interactions are linked to energy transfer across frequencies and to cascade processes in turbulent flows (Kinjangi & Foti 2023). Their influence becomes more pronounced in cases with mean pitch, where energy is distributed over a broader range of frequencies. Furthermore, the QSE decomposition recovers high-frequency dynamics even though the retained POD modes are dominated by low-frequency content (see figure 5). This suggests that the quadratic interactions among POD modes may generate triadic couplings capable of reconstructing higher-frequency behaviour.

Interestingly, some time-average behaviour also appears in the quadratic modes. This is due to the fact that quadratic terms like $\psi ^{(i)}\psi ^{(i)}$ are always non-negative, and thus they will produce a net contribution over an oscillation cycle. As a result, these modes can carry steady components that, potentially, could also be attributed to the mean mode $(0,0)$ . Similar behaviour is present also for other harmonic components in higher-order modes, which could potentially lead to information mixing between the time-constant and time-varying parts of the reconstruction in the QSE framework, as a result of the quadratic temporal modes not being orthogonal. To address this problem, we applied an a priori truncation criterion based on uncertainty quantification, ensuring this overlap between the dynamics does not lead to unacceptable information mixing in the quadratic modes. This further demonstrate the idea that proper truncation before performing QSE is of paramount importance.

Taken together, these observations support a physically grounded interpretation of the QSE-POD modes. The linear modes appear to capture the effects of vorticity generation and shedding synchronised with the wing’s motion, producing pressure forces that arise from the interaction of these structures with the mean velocity field, as dictated by the structure of the Galerkin-projected Poisson equation. This includes contributions from both LEVs and TEVs, as well as bound vorticity around the wing. This is supported by the fact that, in flapping-wing aerodynamics, LEVs are closely tied to the primary motion (Dickinson, Lehmann & Sane 1999), and TEVs also shed once per half-stroke in typical cases (Rival, Schönweitz & Tropea 2011). As such, pressure fluctuations that mirror the flapping frequency are consistent with these near-field, motion-synchronised structures. Furthermore, the pressure gradients in the linear modes predominantly align with the chord-normal direction, in agreement with force contributions linked to bound vorticity and vortex formation at the wing’s surface (Sane 2003), and this behaviour is altered when mean flow and chord directions do not align.

In contrast, the quadratic pressure modes likely represent nonlinear interactions between flow features, which would arise from flow–vorticity interactions and added-mass terms as derived from force-element analysis (Chang 1992). These modes exhibit higher-frequency content and spatial patterns aligned with the wake, suggesting they may encode a complex dynamics like LEV-TEV mutual and self-induction interactions in the $\boldsymbol{u}\times \boldsymbol{\omega }$ term. Although these interpretations are not directly enforced in the data-driven decomposition, they provide a plausible aerodynamic explanation for the observed spatio-temporal structures and their influence on unsteady force production.

3.4. Comparison between QSE-POD and LSE-POD

So far we have discussed the possibilities that QSE-POD offers for obtaining physical insights, but to study its impact it is key to verify its force modelling capabilities. To do so, we will compare QSE-POD with other stochastic estimation alternatives for flow-based modal decomposition. The LSE-POD method (Raiola et al. 2021) has demonstrated good performance in modelling forces in flapping-wing problems. Briefly, LSE-POD projects the measured forces onto the basis formed by the velocity-field POD modes, yielding the estimated contribution that these modes make to the forces. This way, LSE-POD uses these contributions to predict the future expected forces. As discussed by Raiola et al. (2021), 3 modes (i.e. a LOR with threshold $n_m = 2$ ) are necessary to recreate the most relevant flow features, representing the circulation evolution around the airfoil. These modes contribute to the bulk of aerodynamic forces, in particular in the chord-normal direction. Additional modes (i.e. a LOR with threshold $n_m = 6$ ) allow us to include higher-order shedding phenomena in the wake of the airfoil, which refines force models, especially in the chordwise direction. As the QSE-POD framework will increase the number of modes through quadratic modes, QSE-POD with $n_m = 2$ will be compared against LSE-POD with $n_m = 6$ . This will ensure that the differences in quality are not due to significant differences in the size of the modal basis used. Although smaller in magnitude, chordwise pressure forces play a key role in capturing oscillations in axial loading. Their accurate reconstruction is critical for resolving unsteady thrust dynamics, particularly in the presence of wake asymmetries and vortex shedding. This is especially relevant in data-driven models, where capturing these fluctuations can significantly improve predictions of time-resolved aerodynamic forces.

This is where QSE offers a clear advantage. Because it incorporates nonlinear (quadratic) interactions, QSE is better equipped to capture the physical mechanisms underlying chordwise force generation. As shown in figure 12, QSE-POD demonstrates notably improved accuracy in reconstructing axial forces, while LSE-POD performs slightly better for normal forces , within the bounds of measurement uncertainty. This trend is consistent across all experimental cases: QSE yields a mean axial residual of $r_{C_x} = \ 5.5\,\%$ and a normal force residual of $r_{C_y}= \ 5.6\,\%$ compared with $11.3\,\%$ and $3.0\,\%$ , respectively, for LSE, where $r_{C_x} = \lVert C_{x,\textit{model}}-C_{x,\textit{poisson}}\rVert ^2/\lVert C_{x,\textit{poisson}}\rVert ^2$ . These results reflect the underlying structure of the problem: while the chord-normal force is dominated by near-linear dynamics that LSE can approximate, axial forces arises more strongly from nonlinear flow interactions that QSE is designed to capture.

Figure 12. (a) Case with $\theta _m=0^\circ$ and $\theta _0=10^\circ$ . (b) Case with $\theta _m=10^\circ$ and $\theta _0=10^\circ$ . Normal (left) and axial (right) force coefficients as estimated from the Poisson equation (black), QSE-POD with $n_m = 2$ (red) and LSE-POD with $n_m = 6$ (blue).

Although QSE-POD achieves accurate force reconstruction using only the first few velocity modes, this does not imply that higher-order modes are irrelevant. In QSE, higher-frequency force fluctuations can be captured through quadratic interactions between low-order modes, effectively encoding some of the dynamics that would otherwise require higher-order modes in a linear framework. However, this does not exclude the possibility that higher-order velocity modes also contribute to the force through their own quadratic interactions or linear correlations with pressure features. This is probably the reason why LSE with $n_m = 6$ can outperform the reconstruction QSE with $n_m=2$ even though they share the same number of modes, as LSE is likely to encode higher-order dynamics in its higher-order modes. This is supported by Mohammadi, Morton & Martinuzzi’s (2025) comments on how POD modes can contain high-order dynamics. It remains an open question to what extent these higher-order modes are dynamically independent or instead reflect nonlinear couplings with the dominant flow structures. This could also explain why linear models like LSE-POD still capture certain higher-frequency effects up to a point, despite relying solely on low-dimensional projections.

This is supported by figure 3 (right), where the cumulative gain in reconstructed force variance $S_p$ begins to plateau around $0.9$ , even as additional modes are included. This suggests that, although QSE captures a large portion of the dynamics with a compact basis, a residual contribution remains that could be linked to higher-order modal interactions. As such, the success of QSE should not be seen as a reason to dismiss the role of higher-order velocity features altogether, but rather as evidence that much of the dynamic content attributed to such modes can be efficiently captured through nonlinear (i.e. quadratic) interactions among a reduced number of dominant modes, leaving only a small residual potentially associated with genuinely higher-order effects.

4. Conclusion

This work introduces a data-driven framework that links experimentally observed flow structures to sectional pressure loads on physical grounds. The method involves extracting the POD modes from a measured velocity-field dataset. Leveraging the well-known quadratic relationship between pressure and velocity fields, the pressure-Poisson equation ${\nabla} ^2p = \boldsymbol{\nabla }\,\boldsymbol{\cdot }\,[(u\,\boldsymbol{\cdot }\,\boldsymbol{\nabla })u]$ ), the approach applies QSE to identify the quadratic modal decomposition $P^{(i,j)}$ of the pressure distribution. This framework is here referred to as QSE-POD.

Applied to flapping-wing experiments, QSE-POD produces compact and physically meaningful reduced-order models. With only two velocity modes, it reconstructs dominant flow features, distinguishes linear and quadratic contributions to force generation and achieves accurate predictions of both normal and axial forces. Compared against LSE-POD (i.e. EPOD), QSE-POD provides more accurate reconstructions of axial forces, while requiring fewer modes and preserving comparable accuracy for normal forces. This highlights the efficiency of the quadratic formulation in capturing quadratic dynamics, which becomes increasingly relevant in cases with non-zero-mean pitch. This balance of accuracy and compactness highlights its value as a reduced-order alternative to more data demanding approaches.

Beyond the flapping-wing aerodynamics, the proposed framework provides a general tool for modelling nonlinear velocity–pressure couplings in unsteady flows. Its flexibility suggests promising applications in aeroacoustics, where velocity-based measurements are often more accessible than pressure data. By exploiting the quadratic structure of the pressure-Poisson equation, QSE-POD can be applied to problems ranging from airfoil noise to turbulent jet noise, offering a reduced-order means to connect coherent velocity structures with their acoustic signature. More broadly, the framework could be integrated into flow-control strategies, where accurate yet compact estimation of unmeasured flow quantities is critical for real-time decision making.