The evolution of vortices determines the aeroacoustics generated by a hovering wing

Abstract

The effects of the evolution of vortices on the aeroacoustics generated by a hovering wing are numerically investigated by using a hybrid method of an immersed boundary–finite difference method for the three-dimensional incompressible flows and a simplified model based on the Ffowcs Williams-Hawkings acoustic analogy. A low-aspect-ratio ($AR=1.5$) rectangular wing at low Reynolds ($Re=1000$) and Mach ($M=0.04$) numbers is investigated. Based on the simplified model, the far-field acoustics is shown to be dominated by the time derivative of the pressure on the wing surface. Results show that vortical structure evolution in the flow fields, which is described by the divergence of the convection term of the incompressible Navier–Stokes equations in a body-fixed reference frame, determines the time derivative of the surface pressure and effectively the far-field acoustics. It dominates over the centrifugal acceleration and Coriolis acceleration terms in determining the time derivative of the surface pressure. The position of the vortex is also found to affect the time derivative of the surface pressure. A scaling analysis reveals that the vortex acoustic source is scaled with the cube of the flapping frequency.

1. Introduction

The aeroacoustics generated by flapping wings is ubiquitous in aerial animal flight. Wing-flapping-generated sound (Clark 2021) is essential in courtship (Spieth 1974; Gibson & Russell 2006) and sexual recognition (Warren, Gibson & Russell 2009) for many insects. Some insects, such as Drosophila melanogaster, can respond to sounds that share similar frequency characteristics to their own species (Robert & Göpfert 2002). Apart from biological systems, aeroacoustics generated by flapping wings is of importance in micro aerial vehicles (MAVs), especially in applications like military reconnaissance and environmental monitoring where quiet operation is crucial, because applications of flapping wings (Platzer et al. 2008; Shyy et al. 2010; Ravi et al. 2016) in the MAVs (Young & Garratt 2020; Chen et al. 2024) have been growing due to their high-lift performance.

In order to explain the aeroacoustics generated by flapping wings, it is necessary to introduce their force generation and the associated mechanisms, of which the majority are about the vortical structures. Early efforts on flapping wings were focused on their force generation (Ellington 1984; Platzer et al. 2008; Shyy et al. 2010). In particular, the influences of parameters (such as kinematics, aspect ratio and Reynolds numbers) on forces have been investigated. The kinematics of simplified flapping motions on the force generation have been extensively studied, including plunging and pitching motions in forward flight (Lai & Platzer 1999; Young & Lai 2007; Visbal 2009; Kim & Gharib 2010), and translating/stroke and pitching motions in hovering flight (Birch & Dickinson 2001; Poelma, Dickson & Dickinson 2006; Garmann, Visbal & Orkwis 2013; Medina & Jones 2016). It is found that unsteady flows over the flapping wing generate vortices near the wing edges which play an important role in the force generation in addition to the added mass (Chang 1992; Sane & Dickinson 2001). These vortical structures are demonstrated in figure 1, including leading-edge vortex (LEV) (Ellington et al. 1996; Ford & Babinsky 2013; Eldredge & Jones 2019), trailing-edge vortex (TEV) (Shyy et al. 2010; Han et al. 2015), tip vortex (TV) (Ringuette, Milano & Gharib 2007; Shyy et al. 2009; Kim & Gharib 2010), root vortex (RV) (Shyy et al. 2010) and sometimes secondary vortex (SV) (Lu, Shen & Lai 2006). The LEV is the most prominent flow structure. It creates a low-pressure zone above the wing, responsible for most of the lift generation during the translation/middle stroke phase (Ellington et al. 1996; Liu et al. 1998). Garmann & Visbal (2014) found that the aerodynamic force coefficients of revolving wings are hardly changed when the aspect ratio is larger than 2, as the chordwise evolution of LEV at outboard region is limited by the trailing edge, leading to the loss of the negative pressure above the wing. Jones & Babinsky (2011) investigated a waving wing at different Reynolds numbers, finding that the LEV sheds earlier at lower Reynolds numbers, leading to the earlier disappearance of the low-pressure zone and consequently a reduction of the aerodynamic forces. Other vortices contribute to the force in a similar way, i.e. through creating a low-pressure zone. Their contributions to the force are affected by shapes, kinematics and material properties of the wing (see e.g. Tian et al. 2013; Shahzad et al. 2016; Huang et al. 2023). Therefore, their evolution and interaction are quite diverse, creating complexities in studying flapping wings. For example, the interaction between LEV and the boundary layer could generate an SV which influences the force generation (Lu et al. 2006). Overall, the variations of these parameters will lead to the change of the vortex dynamics around the flapping wing, and thus modify the aerodynamic forces which could ultimately affect the aeroacoustics generated by the wings.

Figure 1. Vortical structures around a flapping wing: LEV, TEV, TV, RV and SV.

Conservation of momentum and force partition methods can be used to uncover the relationship between the vorticity and forces (Wu, Ma & Zhou 2007). Wu (1981) proposed the vorticity moment theory based on the conservation of momentum, connecting the aerodynamic force of a body moving in a fluid to the time derivative of the total first moment of the vorticity field in an open space. Based on a derivative of moment transformations, Wu, Pan & Lu (2005) further proposed two formulae to account for the aerodynamic forces on an arbitrary body in the form of control surface integral including the viscous vortical effects. This method was later successfully adopted by Li & Lu (2012) to calculate the aerodynamic forces of flapping plates by considering the ring-like vortical structures in the wake, finding that the force on a flapping plate is dominated by the flow structures close to the plate. A force and moment partitioning method is a powerful tool for analysing complex vortex flows through delineating the contribution of vortex-induced effects (Wu et al. 2007; Menon, Kumar & Mittal 2022). For example, Seo, Menon & Mittal (2022) proposed a force partitioning method for studying the vortex evolution on the surface pressure fluctuation and acoustics of a two-dimensional (2-D) flapping foil. Further discussion on the relationship between vortices and aerodynamic forces can be found in Wu, Liu & Liu (2018). As the vortex flows dominate in flapping wings, their forces can be estimated by considering the time derivative of the vortex impulse, which is the basis of scaling laws for the aerodynamic forces of flapping wings (see e.g. Dewey et al. 2013; Kang & Shyy 2013; Lee, Choi & Kim 2015; Ehrenstein 2019; Liu, Liu & Huang 2022). In these scaling laws, the flapping frequency with its intrinsic connection with vortices’ evolution period is considered as an important factor for scaling the aerodynamic forces.

The aerodynamic forces acting on a moving body, which are determined by the vortex flows, have strong effects on the acoustics (Seo, Hedrick & Mittal 2019; Tian 2020). Therefore, the vortical flows as well as their evolution would have a strong impact on the acoustics. However, the aeroacoustic characteristics have received less attention compared with the aerodynamics and vortical flows of flapping wings but have attracted growing effort recently. The aeroacoustics of flapping wings have been investigated experimentally (see e.g. Sueur, Tuck & Robert 2005; Arthur et al. 2014; Zhou et al. 2019) and numerically (see e.g. Nagarajan & Lele 2005; Nagarajan, Hahn & Lele 2006; Bae & Moon 2008; Inada et al. 2009; Wang & Tian 2020). The noise generated by a flyer in forward and aft directions is dominated by the first harmonic (flapping frequency, $f_o$) while the noise at the sides is dominated by the second harmonic ($2f_o$) (Sueur et al. 2005; Arthur et al. 2014). This is because that the wing creates one-periodic thrust in the forward and aft directions, and two-cycle lateral forces in side directions during one flapping cycle. Vorticity features also have an impact on acoustics. For example, the vortex shedding frequency near the trailing edge of a pitching NACA0012 aerofoil at high Reynolds numbers corresponds to the central frequency of the noise hump in the spectrum (Zhou et al. 2019). For a 2-D elliptical wing undergoing hovering and forward flights, a dipole source of flapping frequency is generated from the transverse motion and higher-frequency dipole sources arise from trailing edge scattering during tangential motion (Bae & Moon 2008). The dipole sources are generated by the pressure difference between the two sides of the wings (Inada et al. 2009). A more detailed study (Seo et al. 2019) showed that the time derivative of the surface pressure is the most important noise source for the aeroacoustics of flapping wings at low Mach numbers (also see Tian 2020; Clark 2021). Flexibility of wings could enhance the mean force production and reduce the noise (Geng et al. 2017).

As noted above, the aerodynamic force fluctuations are correlated with the surface pressure variation caused by vortical structure evolution, hence new insights into the mechanism of aeroacoustics generation could be obtained by investigating the vortical structure evolution. Nedunchezian, Kang & Aono (2019) found that the kinematics of the wing play an important role in determining the noise level and that the variations of the vortical structures’ topology may account for the difference in noise level. Clark (2021) discussed the vortex–wing interaction and acknowledged that the noise generated by vortices can be promoted if the vortices are in the vicinity of a wing (also noted by Curle (1955)). By using a force partitioning method, Seo et al. (2022) investigated the influences of each vortical structure and kinematic parameters on the aeroacoustics of a 2-D pitching NACA0015 aerofoil finding that the LEV plays an important role in determining the wing surface pressure and thus the far-field acoustics. However, the influence of characteristic vortices on surface pressure fluctuations in the three-dimensional (3-D) case can be evaluated in a more straightforward way without employing the projection method (Seo et al. 2022). In addition, it remains unclear which characteristic vortex is dominant in affecting the surface pressure fluctuations in the 3-D case.

To gain insights into the flow and acoustic fields associated with flapping wings, computational fluid dynamics (CFD) is an attractive method. There are two types of CFD methods used for aeroacoustic studies: (i) direct numerical simulation (DNS) and (ii) hybrid methods. In DNS, the compressible Navier–Stokes (NS) equations are solved by high-resolution numerical schemes in which all fluid scales and sound waves are resolved (Colonius, Lele & Moin 1997; Sandberg et al. 2009). Direct numerical simulation has been successfully used to investigate the aeroacoustics of flapping wings (Wang & Tian 2019; Wang, Tian & Lai 2020). Despite its successful applications, the large computational cost has hindered DNS usage, especially for low Mach number flows ($M<0.1$) where the size of the time step may be restricted more by the time scale of acoustic wave propagation than by that of hydrodynamics. To address these challenges, hybrid methods (Hardin & Pope 1994; Lele 1997; Colonius & Lele 2004) have emerged as a promising approach. In these methods, the hydrodynamic field is acquired through direct solutions of the compressible or incompressible NS equations. Concurrently, the acoustic field is evaluated using simplified models (e.g. linearised Euler equations (Bailly & Juve 2000) and Green's function coupling with acoustic analogy (Lyrintzis 2003; Farassat 2007)). The hybrid methods (e.g. coupling the incompressible simulation with the acoustic analogy) could relax the time step restriction and the grid requirement for resolving the acoustic wave propagation in the DNS at low Mach number flows. Therefore, a hybrid method coupling an incompressible solver and Farassat formulation 1A is used in this work.

This work is to develop a connection between the pressure fluctuation and the vortex evolution for a 3-D hovering wing where rotational motions are involved, and a simple scaling law for noise source. The dominant vortex on the time derivative of the surface pressure and its evolution will be addressed. To achieve this, a zero-thickness low-aspect-ratio ($AR=1.5$) rectangular hovering wing at a Reynolds number of 1000 and a Mach number of 0.04 is numerically investigated by solving the incompressible NS equations with the diffused interface immersed boundary (IB) method. By using the wing model, the Farassat formulation 1A is first validated. A simplified acoustic model based on the Ffowcs Williams-Hawkings (FW-H) acoustic analogy is expressed in the frequency domain and validated against the Farassat formulation 1A. Furthermore, the variation of the surface pressure is analysed and the numerical results are incorporated to illustrate the underlying relationship between the surface pressure variation and vortices evolution. Finally, a scaling analysis of the vortex acoustic source is conducted.

This paper is organised as follows. The physical model and numerical method are described in § 2 with validations, mesh and computational domain convergence studies presented in the Appendix (A). The numerical results are presented and discussed in § 3, and concluding remarks are provided in § 4.

2. Physical model and numerical method

2.1. Physical model

As shown in figure 2, a zero-thickness rectangular wing with a chord length of $c$, and a wingspan of $1.5c$, is considered. A low aspect ratio ($AR=1.5$) is used due to its high power economy and low computational cost (Shahzad et al. 2016). The pivot point is 0.1c away from the wing root. For hovering flight, insect wings typically have three degrees of freedom (Ellington 1984). Here, we assume the wing undergoing stroke rotation around the $z$-axis and pitching motion around the leading edge. The stroke plane deviation is usually very small (Chen et al. 2016; Shahzad et al. 2018), and thus ignored here. The stroke and pitching motions, akin to those described in Dai, Luo & Doyle (2012), are adopted,

(2.1a,b)\begin{equation} \varPhi(t)=\frac{\rm \pi}{5}\sin\left(2{\rm \pi} f_ot+\frac{\rm \pi}{2}\right), \quad \alpha(t)=\frac{\rm \pi}{6}\sin(2{\rm \pi} f_ot), \end{equation}

where $\varPhi$ is the angle between the leading edge and the $y$-axis, $\alpha$ is the angle between the chord line and the $z$-axis, $f_o = 0.25$ is the flapping frequency and $t$ is the time ( $f_o$ and $t$ are normalised by $U_{ref}$ and $c$). Note that a small stroke angle is used to reduce the computational cost. According to the scaling laws and previous studies (see e.g. Lee et al. 2015; Wang & Tian 2020), this choice is sufficient to achieve the objective of this work.

Figure 2. The schematic shows dimensions and definitions of the stroke angle ($\varPhi$) and the pitching angle ($\alpha$) for a zero-thickness 3-D rectangular flapping wing.

Here, a Reynolds number ($Re=U_{ref}c/\nu$) of 1000 and a Mach number ($M=U_{ref}/c_o$) of 0.04 are chosen, where $U_{ref}$ is the mean tip velocity and defined as $4{\rm \pi} /5 f_o (AR+0.1)c \approx 4cf_o$, $\nu$ is the kinematic viscosity and $c_o$ is the sound speed. The chosen Reynolds number and the Mach number fall in the range of some insects (Eldredge & Jones 2019; Hightower et al. 2020) and allow us to investigate the characteristic features of both flow and acoustic fields at a relatively low computational cost.

2.2. Numerical methods

Here the incompressible NS equations solver and the FW-H acoustic analogy are used. A simplified model based on the far-field approximation and the FW-H acoustic analogy is proposed here.

2.2.1. Fluid solver

As the operation speed for bio-inspired flapping wings is low (Hightower et al. 2020), the Mach number is much less than 1. Therefore, the fluid dynamics is governed by the dimensionless 3-D incompressible NS equations,

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

(2.3)$$\begin{gather}\frac{\partial {\boldsymbol u}}{\partial t} + {\boldsymbol u}\boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol u} ={-} \boldsymbol{\nabla} p +\frac{1}{Re}\nabla^2 {\boldsymbol u} + {\boldsymbol f}, \end{gather}$$

where ${\boldsymbol u}$ is the fluid velocity, $p$ is the pressure and ${\boldsymbol f}$ is the body force. Based on previous studies (Jones & Babinsky 2010; Garmann et al. 2013; Eldredge & Jones 2019), the flow fields for flapping wings at a Reynolds number of 1000 are laminar. To verify that the flow is laminar at $Re=1000$, our results of including the dynamic Smagorinsky (DS) model (Germano et al. 1991) as the subgrid model (Smagorinsky constant $C_s$ is chosen as 0.1), not shown here, show insignificant differences (less than 1 %) in the aerodynamic force coefficients from our current simulation without the DS model.

The second-order central difference is used for the spatial derivative while the second-order Crank–Nicolson method is used for time marching. The incompressible NS equations are solved by a projection method consisting of three steps: (i) the momentum equations are solved by a Jacobian-free Newton-Krylov method to obtain the intermediate velocity; (ii) the pressure is updated by solving a Poisson equation with the generalised minimal residual method where the algebraic multigrid method is used as a preconditioner; and (iii) the velocity for the current step is updated based on the solutions of the Poisson equation and the intermediate velocity. More details can be found in Calderer et al. (2015).

The IB method (Peskin 1972; Mittal & Iaccarino 2005; Huang & Tian 2019) is used to handle the boundary condition on the wing surface. Compared with the arbitrary Lagrangian–Eulerian method (Hirt, Amsden & Cook 1974; Farhat & Lakshminarayan 2014), the fluid grid in the IB method does not have to conform to the geometry of the solid surface or be regenerated if moving boundaries are involved. Thus, the IB method is more efficient in coping with moving boundary problems, as considered here. The flow field is discretised with non-uniform Cartesian grids, and the solid surface is discretised by triangle grids. The diffused interface IB method (Peskin 1972; Goldstein, Handler & Sirovich 1993; Peskin 2002; Huang & Tian 2019), which is easy to implement and can provide reliable solutions (Huang & Tian 2019; Huang et al. 2022), is used to implement the boundary condition at the fluid–structure interface. More details of this method can be found in Ji et al. (2022, 2024).

Here a computational domain with the dimensions $-5c \le x \le 5c$, $-5c \le y \le 5c$ and $-5c \le z \le 5c$, is used. A fine grid region ($2c \times 1.7c \times 1.2c$) is discretised by uniform fine grids with its size denoted as $\varDelta _{min}$ while an outer grid region is discretised by non-uniform grids. Far-field boundary conditions are enforced at the outer boundaries while the diffused interface IB method is incorporated to achieve the no-slip boundary condition at the wing surface. The fluid solver was validated in the previous works (Ji et al. 2022, 2024). Convergence studies of the grid and computational domain have been conducted with details presented in the Appendix (A). Based on the convergence studies, $\varDelta _{min}=0.01c$ and a computational domain of $10c \times 10c \times 10c$ are used for all simulations in this work.

2.2.2. Acoustic solver

To consider the sound generated by a moving boundary, Ffowcs Williams & Hawkings (1969) used the generalised function and rewrote the Curle equation (Curle 1955) into the FW-H equation,

(2.4)\begin{equation} \left(\frac{1}{c^2_o}\frac{\partial ^2 }{\partial t^2} - \nabla^2 \right)p'(\boldsymbol{x},t)=\frac{\partial}{\partial t}[\rho_o v_n \delta(f)]-\frac{\partial}{\partial x_i}[pn_i\delta(f)]+\frac{\partial^2 H(f)T_{ij}}{\partial x_i \partial x_j}, \end{equation}

where $\rho _o$ is the ambient density, $v_n$ is the surface normal velocity, $p$ is the surface pressure, $n_i$ is the unit component of normal vector of surface, $\delta (f)$ is the Dirac function, $H(f)$ is the Heaviside function, $f=0$ denotes the moving surface and $T_{ij}$ is the Lighthill stress tensor. Note that the three source terms on the right-hand side of (2.4) are monopole, dipole and quadrupole sources, respectively. As the contributions from the quadrupole source are less significant at a low Mach number and the monopole source can be ignored due to the zero thickness (Seo et al. 2019) and incompressible restriction, the integral solution of (2.4) for the dipole source, known as the Farassat formulation 1A (Farassat 2007), is given by

(2.5)\begin{align} 4{\rm \pi} p'_L(\boldsymbol{x},t)&=\oint \left[\frac{\dot{p} \cos \theta}{c_o|\boldsymbol{r}|(1-M_r)^2}+\frac{p\dot{M}_r \cos \theta}{c_o|\boldsymbol{r}|(1-M_r)^3}+\frac{p(M_r-M^2)\cos \theta}{|\boldsymbol{r}|^2(1-M_r)^3}\right.\nonumber\\ &\quad \left.+\frac{p(\cos \theta-M_n)}{|\boldsymbol{r}|^2(1-M_r)^2}\right]_{ret}{\rm d}A, \end{align}

(2.6)\begin{equation} p'(\boldsymbol{x},t)=p'_L(\boldsymbol{x},t), \end{equation}

where $p'_L(\boldsymbol {x},t)$ corresponds to the sound pressure from the dipole source, $|\boldsymbol {r}|=|\boldsymbol {x}-\boldsymbol {y}|$ is the distance from the acoustic source ( $\boldsymbol {y}$) to the observation point ($\boldsymbol {x}$), $M_r = M_ir_i$, $M_n = M_in_i$, $r_i$ is the unit component of $\boldsymbol {r}$, $M_i$ is the $i$th component of the local Mach number vector $\boldsymbol {M}$, $\cos \theta = r_in_i$, ${\rm d}A$ is the area of a surface element, $ret$ indicates the retardation time and the second-order central difference scheme is used to evaluate the time derivative $(^{{\cdot }})$.

2.2.3. Simplified model

A simplified model can be derived for the far-field acoustics to provide a better understanding of the mechanisms of the aeroacoustic generation. Based on the above assumptions, the FW-H acoustic analogy in the far field can be approximated as follows:

(2.7)\begin{equation} \left(\frac{1}{c^2_o}\frac{\partial ^2 }{\partial t^2} - \nabla^2 \right)p'(\boldsymbol{x},t) \approx{-}\frac{\partial}{\partial x_i}[pn_i\delta(f)]. \end{equation}

Given that the acoustic wavelength up to the 10th harmonics is $\lambda =c_o/(nf_o) = 10c$ ($n=10$) which is still much larger than the size of the wing ($1.5c$), the wing can be considered as a compact source. Therefore, the far-field approximation ($\partial /\partial x_i \simeq -r_i/(|\boldsymbol {r}|c_o) \partial /\partial t$, see Howe (2003)) and the free space Green's function are incorporated to obtain the sound pressure in the time and frequency domains,

(2.8)$$\begin{gather} p'(\boldsymbol{x},t)={-}\frac{\partial}{\partial x_i} \oint_S \frac{pn_i\delta(f)}{4{\rm \pi} |\boldsymbol{x}-\boldsymbol{y}|} {\rm d}S \approx \frac{r_i}{4{\rm \pi} |\boldsymbol{r}|^2 c_o}\dot{f_{p_i}}, \quad \dot{f_{p_i}} = \oint \dot{p} n_i\, {\rm d}A, \end{gather}$$

(2.9)$$\begin{gather}p'(\boldsymbol{x},\omega)\approx \frac{{\rm e}^{{\rm i}k|\boldsymbol{r}|} r_i}{4 {\rm \pi}|\boldsymbol{r}|^2 c_o} \dot{f_{p_i}}(\boldsymbol{y},\omega), \quad k = \omega/c_o. \end{gather}$$

The above equations connect the far-field acoustic pressure with the time derivative of the pressure on the wing. Therefore, uncovering the sources of the time derivative of the pressure is of importance in understanding the aeroacoustics generated by the flapping wing. Validations of the Farassat formulation 1A and the simplified model have been conducted with details shown in the Appendix (A).

3. Results and discussions

3.1. Time derivative of surface pressure force

The dominant contributor to the sound pressure in the far-field is the time derivative of the surface pressure $\dot {p}$ and thus the time derivative of surface pressure force $\dot {f_{p_i}}$ (see (2.8)). To discuss it, $\dot {f_{p_i}}$ over a flapping cycle $t/T=0$ to 1 ($T$ is the flapping period) is shown in figure 3(a). At $t/T=0.285$, $\dot {f_{p_i}}$ in the three directions is close to zero, which will be further investigated. On the other hand, the maximum fluctuation of $\dot {f_{p_i}}$ occurs near $t/T=0.42$ when $\dot {f_{p_x}}$ achieves its minimum value while the fluctuations of other two components are also high.

Figure 3. Force properties and vortical structures: (a) time derivative of the pressure force in the three directions; (b) instantaneous dimensionless surface pressure at $t/T=0.285$; (c) instantaneous dimensionless time derivative of the surface pressure at $t/T=0.285$; (d) vortices visualised by $Q$-criterion ($Q=30$) and coloured by pressure at $t/T=0.285$; (e) instantaneous dimensionless surface pressure at $t/T=0.42$; ( f) instantaneous dimensionless time derivative of the surface pressure at $t/T=0.42$; (g) vortices visualised by $Q$-criterion ($Q=30$) and coloured by pressure at $t/T=0.42$.

The instantaneous surface pressure and its time derivative together with the vortices visualised by the $Q$-criterion (Jeong & Hussain 1995) over the suction side at $t/T=0.285$ and $t/T=0.42$ are presented in figures 3(b–d) and 3(e–g), respectively. The $Q$-criterion is used here because it can be taken as the acoustics and the surface pressure variation source, as will be discussed in § 3.3. The negative pressure over the outboard region at $t/T=0.285$ is much larger than that at $t/T=0.42$ (see figure 3b,e) but no significant difference in the vortical structures between these two time steps can be identified (see figure 3d,g). In addition, a $\boldsymbol {\varPi }$-like distribution pattern with a high negative pressure ($p<-0.7$, see figure 3e) can be identified in the inboard region at $t/T=0.42$, which is associated with the low pressure within the LEV, RV and SV (SV1). In particular, the negative pressure near the wing root at $t/T=0.42$ is larger than that at $t/T=0.285$ due to the larger RV. At $t/T=0.285$, $\dot {p}$ over the outboard region is positive which compensates for the negative value over the inboard region, yielding a small amplitude $\dot {f_{p_i}}$ (see figure 3a,c). At $t/T=0.42$, $\dot {p}$ over the suction side is generally positive, indicating that the surface pressure increases, and a region R1 with a noticeably high value is found at around a chord from the wing root (see figure 3f). Although the value of $\dot {f_{p_i}}$ can be explained by the contours of $\dot {p}$, the relationship between the vortical structure evolution and $\dot {p}$ cannot be determined from figure 3. Furthermore, it should be noted that $\dot {p}$ generally increases from the wing root to the wingtip at $t/T=0.285$. However, R1 does not occur at the wingtip at $t/T=0.42$, which may contain some important physical insights and shall be investigated. The LEV, as a prominent vortical component, may significantly affect $\dot {p}$ and result in the high $\dot {p}$ at R1 (Eldredge & Jones 2019). Therefore, the LEV circulation will be quantified and examined in the next section.

3.2. The LEV evolution

To investigate the correlation between the LEV and the high $\dot {p}$ of R1, a body-fixed reference frame $s$–$n$–$\tau$, with its origin at the corner closest to the pivot point of the wing (see figure 4a) is first established, where $s$ represents the direction from the trailing edge to the leading edge, $n$ denotes the normal to the wing surface and $\tau$ points from the wing root to the wingtip. The vorticity $\boldsymbol {\omega }=\boldsymbol {\nabla } \times \boldsymbol {u}$ is calculated in the $s$, $n$ and $\tau$ directions, respectively, and the circulation of the LEV, $\varGamma _{\tau }=\iint \omega _{\tau } \,{\rm d}A$, is defined here to quantify the vorticity strength of the LEV and its variation near R1 (see figure 3f). In figure 4(b), the integration area $A$ is defined as a rectangular region to exclude the contribution of the vorticity of SV1 (see figure 3d,g). Additionally, the LEV is bounded by $Q>0$ and $\omega _{\tau }<0$ (similar definition can be found in Chen et al. (2022)). Here, two slices S1 ($\tau =0.75c$) and S2 ($\tau =0.95c$) (see figure 4c,d), are chosen for calculating the LEV circulation. Here S1 is placed out of R1 while S2 crosses R1. It should be noted that the LEV will lose its coherence at the further outboard area where the tilting of the TV may affect the estimation of the LEV circulation. Therefore, the LEV circulation at the outboard region is not considered here. Apart from $\varGamma _{\tau }$, the vorticity fluxes $\hat {\omega }_s$ and $\hat {\omega }_n$ at S1 and S2 are obtained by integrating the other two vorticity components within the LEV region to show the strength of vorticity in the $s$ and $n$ directions.

Figure 4. The LEV structures: (a) schematic of the body-fixed reference frame $s$–$n$–$\tau$ (marked by the red colour); (b) schematic of circulation integration area and the LEV is highlighted by the contours of $\omega _\tau$; (c) schematics of S1–3, with the contours of $\dot {p}$ at $t/T=0.285$; (d) schematics of S1–3, with the contours of $\dot {p}$ at $t/T=0.42$; (e) the LEV circulation $\varGamma _\tau$, and vorticity fluxes $\dot {\hat {\omega }}_s$ and $\dot {\hat {\omega }}_n$ of S1 and S2.

Overall, it can be observed from figure 4(e) that $\varGamma _\tau$ is larger with respect to the vorticity fluxes in the other two directions ($\hat {\omega }_s$ and $\hat {\omega }_n$) for both slices, and the values at S2 are generally larger than those at S1, suggesting that the LEV is more 3-D and unsteady near R1 over the inboard region. At S1, $\hat {\omega }_s$ and $\hat {\omega }_n$ have a similar variation tendency and share a similar maximum value of around 0.35 at $t/T=0.38$. At S2, the variation of $\hat {\omega }_s$ and $\hat {\omega }_n$ is different after achieving its maximum value at $t/T=0.33$. At $t/T=0.285$, the variation trend of $\varGamma _\tau$ at S1 and S2 is similar, which may explain why $\dot {p}$ has a similar magnitude near S1 and S2 (see figure 4c). Meanwhile, $\varGamma _\tau$ of S1 and S2 at $t/T=0.42$ increases, indicating a reduction of the LEV circulation, and the increase rate of $\varGamma _\tau$ at S2 is larger with respect to that of S1, corresponding to the higher $\dot {p}$ at R1 (see figure 4d). Moreover, the difference in the variation rate at $t/T=0.42$ between S1 and S2 is larger than that at $t/T=0.285$.

Additionally, the cross-correlations between $\dot {\varGamma }_\tau$, $\dot {\hat {\omega }}_s$, and $\dot {\hat {\omega }}_n$ and $\dot {p}$ are incorporated to investigate and quantify the relationship between the time derivatives of the pressure and vortex evolution,

(3.1a,b)\begin{align} R_{\dot{\varGamma}_\tau \dot{p}} = \frac{1}{N}\sum_{n=1}^{N} \frac{(\dot{\varGamma}_\tau - \overline{\dot{\varGamma}_\tau})(\dot{p} - \bar{\dot{p}})}{\sigma_{\dot{\varGamma}_\tau} \sigma_{\dot{p}}}, \quad R_{\dot{\hat{\omega}}_i \dot{p}} = \frac{1}{N}\sum_{n=1}^{N} \frac{(\dot{\hat{\omega}}_i - \overline{\dot{\hat{\omega}}_i})(\dot{p} - \bar{\dot{p}})}{\sigma_{\dot{\hat{\omega}}_i} \sigma_{\dot{p}}}, \quad i=n,s, \end{align}

where $(^{{\cdot }})$ denotes the time derivative, $N$ is the total sample number, $\dot {\varGamma }_\tau$ and $\dot {\hat {\omega }}_i$ are the time derivatives of the LEV circulation and vorticity fluxes of the LEV in the $n$ and $s$ directions at S2, $\overline {({\cdot })}$ denotes the average values and $\sigma _{\square }$ denotes the root mean square (r.m.s.) value.

Three points (${\rm P}_{\rm S1}$, ${\rm P}_{\rm S2}$ and ${\rm P}_{\rm S3}$) located at the intersection of three slices with the midchord of the wing (see figure 4c,d) have been selected to investigate as the midchord line crosses R1. It is believed that ${\rm P}_{\rm S1}$ and ${\rm P}_{\rm S2}$ are significantly influenced by the evolution of the LEV, while ${\rm P}_{\rm S3}$ is chosen for comparison as it is less influenced by the LEV. Here $R_{\dot {\varGamma }_\tau \dot {p}}$ and $R_{\dot {\hat {\omega }}_i \dot {p}}$ over $t/T=0$ to 0.5 at the three points are summarised in table 1. The sample number $N$ is 50 with an interval of $0.01T$. Here $R_{\dot {\varGamma }_\tau \dot {p}}$ at ${\rm P}_{\rm S1}$ and ${\rm P}_{\rm S2}$ is high (with absolute values larger than 0.8), indicating that the time derivative of surface pressure at these two points is significantly related to the variation of the LEV circulation. Here $R_{\dot {\hat {\omega }}_s \dot {p}}$ and $R_{\dot {\hat {\omega }}_n \dot {p}}$ share a similar value at ${\rm P}_{\rm S1}$ while $R_{\dot {\hat {\omega }}_n \dot {p}}$ ($-$0.79) is more significant than $R_{\dot {\hat {\omega }}_s \dot {p}}$ ($-$0.67) at ${\rm P}_{\rm S2}$. One possible explanation is that the dual stage vortex tilting for the flapping wing (see Chen, Cheng & Wu (2023a); Chen et al. (2023b) for more information) have reoriented the LEV and converted the vorticity from the $\tau$ axis to the $n$ axis. Conversely, a weak cross-correlation is observed at ${\rm P}_{\rm S3}$, suggesting that the pressure variation in the outboard region is less influenced by the LEV, which loses its coherency near the wingtip.

Table 1. Cross-correlation between $\dot {\hat {\omega }}_s$, $\dot {\hat {\omega }}_n$ and $\dot {\varGamma }_\tau$ at S2 and $\dot {p}$ of ${\rm P}_{\rm S1}$, ${\rm P}_{\rm S2}$ and ${\rm P}_{\rm S3}$.

3.3. Analytical analysis on the vortices and the time derivative of the surface pressure

From the previous section, the evolution of the LEV has been identified as an important factor in determining the time derivative of the surface pressure. Here, the role of the vortices evolution on the time derivative of the surface pressure is further investigated. To do so, an analytical expression of the time derivative of the surface pressure is established with all source terms placed on the right-hand side of the expression. An analysis is conducted to investigate the relationship between the source terms and vortices evolution by incorporating the simulation results. To obtain the expression, the time derivative of the divergence of the momentum equations in conjugation with the incompressible restriction expressed in the body-fixed reference frame $s$–$n$–$\tau$ (Batchelor 1967; Garmann et al. 2013; Jardin & David 2015) is as follows:

(3.2)\begin{align} \nabla^2 \frac{\partial p}{\partial t}=\rho\left(2\frac{\partial \hat{Q}}{\partial t} -\underbrace{\frac{\partial [\boldsymbol{\nabla} \boldsymbol{\cdot} (\dot{\boldsymbol{\varOmega}} \times \boldsymbol{r_o})]}{\partial t}}_{\textit{Angular acceleration}} - \underbrace{\frac{\partial [\boldsymbol{\nabla} \boldsymbol{\cdot} ( \boldsymbol{\varOmega} \times (\boldsymbol{\varOmega} \times \boldsymbol{r_o}) )]}{\partial t}}_{\textit{Centrifugal\ acceleration}} -\underbrace{\frac{\partial [\boldsymbol{\nabla} \boldsymbol{\cdot} (2\boldsymbol{\varOmega} \times \boldsymbol{v})]}{\partial t}}_{\textit{Coriolis\ acceleration}} \right), \end{align}

where $\boldsymbol {\varOmega }$ is the rotation rate of the flapping and pitching motions, $\boldsymbol {r_o}$ is the vector from the rotation centre to a given point, $\boldsymbol {v}=\boldsymbol {u}-\boldsymbol {\varOmega } \times \boldsymbol {r_o}$ is the relative velocity and $\hat {Q} =-1/2( \boldsymbol {\nabla } \boldsymbol {v}: \boldsymbol {\nabla } \boldsymbol {v})$ which is analogous to the $Q$-criterion for incompressible flow ($Q=-1/2(\boldsymbol {\nabla } \boldsymbol {u}:\boldsymbol {\nabla } \boldsymbol {u})$, see Jeong & Hussain (1995)).

In (3.2), the divergence of the term $\dot {\boldsymbol {\varOmega }} \times \boldsymbol {r_o}$ is zero. Additionally, it should be noted that the centrifugal acceleration term $-4 \boldsymbol {\varOmega }\boldsymbol {{\cdot }} \dot {\boldsymbol {\varOmega }}$ is a constant over the flow field. As the flapping amplitude and frequency are low, it is negligible compared with the Coriolis acceleration (also reported by Garmann et al. (2013)). Therefore, (3.2) can be simplified as

(3.3)\begin{align} \nabla^2 \frac{\partial p}{\partial t}=\rho\left(2\frac{\partial \hat{Q}}{\partial t} + 4\boldsymbol{\varOmega}\boldsymbol{\cdot} \dot{\boldsymbol{\varOmega}} - \frac{\partial [\boldsymbol{\nabla} \boldsymbol{\cdot} (2\boldsymbol{\varOmega} \times \boldsymbol{v})]}{\partial t}\right)\approx 2\rho\left(\frac{\partial \hat{Q}}{\partial t} - \frac{\partial [\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\varOmega} \times \boldsymbol{v})]}{\partial t}\right). \end{align}

From (3.3), the time derivatives of $\hat {Q}$ and the divergence of the Coriolis acceleration are the major sources of the Poisson equation of $\dot {p}$. Furthermore, Green's second identity is incorporated and yields the following equation:

(3.4)\begin{align} &\iiint G(\boldsymbol{r},\boldsymbol{r}')\nabla^2 \frac{\partial p}{\partial t} - \frac{\partial p}{\partial t}\nabla^2 G(\boldsymbol{r},\boldsymbol{r}') \,{\rm d}V'\nonumber\\ &\quad =\iint G(\boldsymbol{r},\boldsymbol{r}') \frac{\partial }{\partial n}\left(\frac{\partial p}{\partial t}\right) - \frac{\partial p}{\partial t} \frac{\partial G(\boldsymbol{r},\boldsymbol{r}')}{\partial n}\,{\rm d}S', \quad \nabla^2 G(\boldsymbol{r},\boldsymbol{r}')=\delta(\boldsymbol{r}-\boldsymbol{r}'), \end{align}

where $G(\boldsymbol {r},\boldsymbol {r}')$ is Green's function, $\boldsymbol {r}$ is the observer location, $\boldsymbol {r'}$ is the source location, $\delta (\boldsymbol {r}-\boldsymbol {r}')$ is the Dirac function, ${\rm d}V'$ is the flow domain and ${\rm d}S'$ is the wing surface. In the body-fixed reference frame, the two terms on the right-hand side of (3.4), ${\partial }/{\partial n}({\partial p}/{\partial t})$ and $\partial G(\boldsymbol {r},\boldsymbol {r}')/ \partial n$, are almost zero if the volume is properly selected. The solution of (3.3) is approximated as

(3.5)\begin{equation} \frac{\partial p(\boldsymbol{r})}{\partial t} \approx 2\rho \iiint G(\boldsymbol{r},\boldsymbol{r}') \left[ \frac{\partial \hat{Q}(\boldsymbol{r}')}{\partial t} -\frac{\partial [\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\varOmega} \times \boldsymbol{v}(\boldsymbol{r}'))]}{\partial t} \right] {\rm d}V', \end{equation}

which indicates that the surface pressure derivative is related to $\partial \hat {Q}/\partial t$ and $\partial [\boldsymbol {\nabla } \boldsymbol {{\cdot }} (\boldsymbol {\varOmega }\times \boldsymbol {v})]/\partial t$ in the body-fixed reference frame. The two terms on the right-hand side of (3.5) at S1–3 (see figure 4) of $t/T=0.285$ and $t/T=0.42$ are presented in figure 5. It should be noted that S3 is included here as $\dot {p}$ at $t/T=0.285$ near S3 is much larger than that of the inboard region, implying that investigating S3 is physically meaningful.

Figure 5. Partial derivative $\partial \hat {Q}/\partial t$ ($Q=0$ marked by the black line) at (a) S1, (b) S2 and (c) S3 of $t/T=0.285$; Coriolis acceleration term at (d) S1, (e) S2 and ( f) S3 of $t/T=0.285$; $\partial \hat {Q}/\partial t$ at (g) S1, (h) S2 and (i) S3 of $t/T=0.42$; Coriolis acceleration term at ( j) S1, (k) S2 and (l) S3 of $t/T=0.42$.

At $t/T=0.285$ (see figure 5a–f), the Coriolis acceleration term ($\partial [\boldsymbol {\nabla } \boldsymbol {{\cdot }} (\boldsymbol {\varOmega }\times \boldsymbol {v})]/\partial t$) is much lower with respect to $\partial \hat {Q}/\partial t$ in the three slices. Generally, it can be found that the region with high $\partial \hat {Q}/\partial t$ is mainly enclosed by the boundary of $Q=0$ (marked by the solid black line, considered as the boundary of vortical structures), which indicates that $\partial \hat {Q}/\partial t$ arises from vortices evolution. It can be seen that the LEV grows in size from S1 to S2. Furthermore, $\partial \hat {Q}/\partial t$ in the proximity of the wing surface is negative due to the absorption of the SV (SV2). In contrast, $\partial \hat {Q}/\partial t$ at the outer side of the LEV is positive due to the vorticity fed by the shear layer emanating from the leading edge. The similar distribution pattern of $\partial \hat {Q}/\partial t$ leads to the similar $\dot {p}$ at the wing surface for S1 and S2. At S3, $\partial \hat {Q}/\partial t$ near the wing surface is less prominent but much higher in the wake region where the complex and unsteady vortical interaction (see figure 3d) may account for the higher $\partial \hat {Q}/\partial t$. The different distribution patterns and magnitudes of $\partial \hat {Q}/\partial t$ account for the difference in $\dot {p}$ between S3 and the inner region.

At $t/T=0.42$ (see figure 5g–l), the Coriolis acceleration term is also much less significant than $\partial \hat {Q}/\partial t$, suggesting that $\dot {p}$ is dominated by the latter. Here $\partial \hat {Q}/\partial t$ within the LEV for S1 and S2 is larger than that at $t/T=0.285$ and the LEV begins to detach and moves downstream. The effects of the variation of the distance between the wing surface ($\boldsymbol {r}$) and the source location ($\boldsymbol {r}'$) are represented by changing the value of $G(\boldsymbol {r},\boldsymbol {r}')$ in the integration of (3.5) and lead to the variation of $\dot {p}$, which will be discussed in the next section. At S3, the high-value region of $\partial \hat {Q}/\partial t$ moves upwards compared with that of $t/T=0.285$. Overall, $\partial \hat {Q}/\partial t$ within the LEV of S1 is less than that of S2, which correlates with the lower $\dot {p}$ at the wing surface near S1. Although the high $\partial \hat {Q}/\partial t$ in the wake region of S3 is comparable with that within the LEV of S2, its larger distance from the wing surface renders its influences less significant, which may explain the lower $\dot {p}$ near S3 than that of R1.

3.4. Position of the LEV

Despite the presence of both $\partial \hat {Q}/\partial t$ and the Coriolis acceleration terms, Green's function in (3.4) may play an important role in the formation of R1 and will be elaborated upon here. Green's function $G(\boldsymbol {r},\boldsymbol {r}')$ is related to the distance between the source $\boldsymbol {r}'$ and the wing surface $\boldsymbol {r}$, but its expression cannot be readily given here. In the body-fixed reference frame, $\boldsymbol {r}$ is set as stationary and the only variable in Green's function is $\boldsymbol {r}'$. This indicates that Green's function can be evaluated by reflecting on the source location. From previous section, it is found that $\partial \hat {Q}/\partial t$ within the LEV is the dominant source for R1 (see figure 5h). Therefore, it can be summarised that the effects of Green's function can be reflected by evaluating the location of the LEV (bounded by $Q>0$ and $\omega _{\tau }<0$). In figure 6(a), the LEV with the contours of $Q$ at four specific time steps ($t/T=0.1$, 0.2, 0.42 and 0.49) is shown. Generally, it can be found that the LEV increases its size and moves downward and backward from $t/T = 0.1$ to 0.4, while less variation can be identified between $t/T=0.4$ and 0.49. Furthermore, the location of the LEV should be quantified. Analogous to the concept of the centroid of vorticity (Huang & Green 2015; Jones et al. 2016), the centroid of $Q$ of the LEV is used to represent the source location and given as follows:

(3.6a,b)\begin{equation} Q_s=\frac{\displaystyle\iint s Q\, {\rm d}A}{\displaystyle\iint Q \,{\rm d}A}, \quad Q_n=\frac{\displaystyle\iint n Q \,{\rm d}A}{\displaystyle\iint Q \,{\rm d}A}, \end{equation}

where $Q_s$ and $Q_n$ are the coordinates of the centroid of the $Q$ within the LEV in the $s$ and $n$ directions, respectively, and $A$ is the region of the LEV (defined in § 3.2).

Figure 6. Vortex evolution: (a–d) contours of $Q$ within the LEV region at $t/T = 0.1$, 0.2, 0.42 and 0.49, respectively; (e) LEV trajectory at S2 $\tau =0.95c$, represented by the centroid of $Q$ (red arrow denotes the LEV evolution direction); ( f) time histories of the centroid of $Q$ ($Q_s$ and $Q_n$) at S2 from $t/T=0.05$ to 0.5.x.

The $Q$ centroid of the LEV at S2 ($\tau =0.95c$) from $t/T=0.05$ to 0.5 is shown in figure 6(b,c). Note that the leading edge is located at $(s,n) = (0,0)$. From figure 6(a), it can be observed that the LEV generally moves from the leading edge to the trailing edge, detaching from the wing surface to downstream. However, the opposite trend can be found near the end of the half-stroke. The time histories of $Q_s$ and $Q_n$ are shown in figure 6(c) and the extreme values occur near $t/T=0.42$ corresponding to the high $\dot {p}$ for R1, as shown in figure 3. Another interpretation of the formation of R1 is that the detachment of the LEV and the variation of $\hat {Q}$ increase the surface pressure, resulting in a high positive $\dot {p}$ at R1.

3.5. A scaling law for the vortex acoustic source

As discussed in the Section 3.3, $\partial \hat {Q}/ \partial t$ arising from the vortical structures evolution is considered as the main source for $\partial p/\partial t$ which is considered as the dominant source of the far-field noise generated by the hovering wing (see (2.9)). In other words, $\partial \hat {Q}/ \partial t$ can be regarded as the noise source from vortex evolution. It is thus worthy to scale the magnitude of $\partial \hat {Q}/ \partial t$ and $\partial p/\partial t$ with kinematic parameters (e.g. $f_o$) to provide a new insight for MAV design.

According to Lee et al. (2015), the vorticity within the LEV of a flapping wing is estimated by

(3.7)\begin{equation} {\omega} \sim \frac{U_{ref} c AR \sin(\alpha)}{A_{vor}(AR + 2)}, \end{equation}

where $A_{vor}$ is the cross area of the vortex, and $\alpha$ is the angle of attack during the middle stroke. Therefore, $\hat {Q}$ can be estimated as

(3.8)\begin{equation} \hat{Q} \simeq \frac{1}{4} \omega^2 \sim \frac{U_{ref}^2 c^2 AR^2 \sin^2(\alpha)}{4A_{vor}^2(AR + 2)^2} = \frac{4{\rm \pi}^2 c^4 AR^2(AR+0.1)^2 \sin^2(\alpha)}{25 A_{vor}^2(AR + 2)^2}f_0^2, \end{equation}

where $U_{ref} = 4{\rm \pi} (AR+0.1)c f_0/5$ is applied. The change of $\hat {Q}$ happens within a period of flapping, and thus

(3.9)\begin{equation} \frac{\partial \hat{Q}}{\partial t} \sim \frac{4{\rm \pi}^2 c^4 AR^2(AR+0.1)^2 \sin^2(\alpha)}{25 A_{vor}^2(AR + 2)^2}f_0^3. \end{equation}

The scaling law for the pressure variation with respect to time is achieved by applying (3.9) into (3.5), i.e.

(3.10)\begin{equation} \frac{\partial p}{\partial t} \sim A_{vor}\partial \hat{Q}/\partial t \sim \frac{4{\rm \pi}^2 c^4 AR^2(AR+0.1)^2 \sin^2(\alpha)}{25 A_{vor}(AR + 2)^2}f_0^3. \end{equation}

The r.m.s. values of $\dot {f_{p_i}}$ (see (2.8)) are incorporated here as an indicator to show whether the $f_o^3$ law is valid. Here, five different flapping frequencies are chosen as $f_o= \{0.05,0.1,0.125,0.2,0.25\}$, yielding five different Reynolds numbers ranging from 200 to 1000. Here $\dot {f_{p_i}}$ is calculated after twenty flapping cycles where the cycle-to-cycle variations in the pressure are negligible. The second-order central difference is adopted to calculate the time derivative of the surface pressure. The r.m.s. values of $\dot {f_{p_i}}$ are evaluated over five flapping cycles as shown in figure 7, where the r.m.s. of $\dot {f_{p_i}}$ is proportional to $f_o^3$ for all three directions. The coefficients in (3.9) and (3.10) appear as a minor shift in the vertical direction, and thus are not explicitly included. Therefore, the vortex acoustic source $\partial \hat {Q}/\partial t$ is scaled by $f_o^3$. It should be noted that the scaling law, (3.10), is a very useful tool in MAV design and other engineering applications.

Figure 7. Scaling law for the time derivatives of pressure forces.

4. Conclusions

The influences of the evolution of vortices on the aeroacoustics of a hovering rectangular wing with a low aspect ratio of 1.5 at Reynolds number of 1000 and Mach number of 0.04 is numerically investigated using a 3-D NS equation solver and diffused interface IB method. A simplified acoustic model based on the FW-H acoustic analogy suggests that the time derivative of the surface pressure ($\dot {p}$) is the dominant noise source for the far-field acoustics. A body-fixed reference frame $s$–$n$–$\tau$, with its origin placed at the corner closest to the pivot point of the flapping wing is considered. The variations of the LEV circulation and vorticity fluxes calculated at two spanwise slices S1 ($\tau =0.75c$) and S2 ($\tau =0.95c$) show that they are related to the variation of $\dot {p}$ based on the cross-correlation analysis. Particularly, a region with high $\dot {p}$ caused by the interaction (vortex tilting effects) between the LEV and the TV is identified near S2, suggesting that the vortices interaction can amplify the noise source. Moreover, the formulation for the $\dot {p}$ in the body-fixed reference frame for the flapping wing is obtained by incorporating Green's function. The time derivatives of the $Q$-value in the body-fixed reference frame ($\partial \hat {Q}/ \partial t$, representing the divergence of the convection), the angular acceleration and the Coriolis acceleration determine the time derivative of the surface pressure. Further investigation suggests that $\partial \hat {Q}/ \partial t$ arising from vortical structure evolution (especially for the LEV), is the dominant source, which overrides the influences of angular and Coriolis accelerations. The position of the LEV at S2 is calculated finding that the shedding of the LEV together with the large magnitude $\partial \hat {Q}/ \partial t$ results in the increase of the surface pressure and the high $\dot {p}$. The analytical framework established in this work can be further extended to investigate the noise sources for other rotating compact surfaces (e.g. UAV propellers). A scaling analysis shows that the vortex acoustic source $\partial \hat {Q}/ \partial t$ is scaled by $f_o^3$. Further parametric investigations of the wing, such as variations in aspect ratio, stroke angle and Reynolds number, could be conducted in future work to gain deeper insight into the relationship between vortex evolution and surface pressure fluctuations.