An immersed boundary method for fluid--structure--acoustics interactions involving large deformations and complex geometries

TL;DR

This paper introduces an immersed boundary method that effectively models fluid-structure-acoustics interactions involving large deformations and complex geometries, validated through benchmark cases and applications to bio-inspired flapping foils.

Contribution

The paper develops a novel IB method combining high-order schemes and finite element modeling for complex fluid-structure-acoustics interactions.

Findings

The IB--WENO scheme outperforms implicit IB--lattice Boltzmann in accuracy.

IB--TENO scheme shows lower dissipation in acoustic wave simulations.

Validated applications include insect flight and energy harvesting with flexible foils.

Abstract

This paper presents an immersed boundary (IB) method for fluid--structure--acoustics interactions involving large deformations and complex geometries. In this method, the fluid dynamics is solved by a finite difference method where the temporal, viscous and convective terms are respectively discretized by the third-order Runge-Kutta scheme, the fourth-order central difference scheme and a fifth-order W/TENO (Weighted/Targeted Essentially Non-oscillation) scheme. Without loss of generality, a nonlinear flexible plate is considered here, and is solved by a finite element method based on the absolute nodal coordinate formulation. The no-slip boundary condition at the fluid--structure interface is achieved by using a diffusion-interface penalty IB method. With the above proposed method, the aeroacoustics field generated by the moving boundaries and the associated flows are inherently solved. In order to validate and verify the current method, several benchmark cases are conducted: acoustic waves scattered from a stationary cylinder in a quiescent flow, sound generation by a stationary and a rotating cylinder in a uniform flow, sound generation by an insect in hovering flight, deformation of a red blood cell induced by acoustic waves and acoustic waves scattered by a stationary sphere. The comparison of the sound scattered by a cylinder shows that the present IB--WENO scheme, a simple approach, has an excellent performance which is even better than the implicit IB--lattice Boltzmann method. For the sound scattered by a sphere, the IB--TENO scheme has a lower dissipation compared with the IB--WENO scheme. Applications of this technique to model fluid-structure-acoustics interactions of flapping foils mimicking an insect wing section during forward flight and flapping foil energy harvester are also presented, considering the effects of foil shape and flexibility.

Figures

Equations

$${\rho }_s\frac{{\partial }^2\mathbf{X}}{\partial t^2}+\frac{\partial }{\partial s}\left[(K_S|\frac{\partial \mathbf{X}}{\partial s}|-1)\frac{\partial \mathbf{X}}{\partial s}\right]+K_B\frac{{\partial }^4\mathbf{X}}{\partial s^4}={\mathbf{F}}_{\mathbf{f}},$$

$$\frac{\partial Q}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}+\frac{\partial H}{\partial z}-\frac{1}{\Re }(\frac{\partial F_u}{\partial x}+\frac{\partial G_v}{\partial y}+\frac{\partial H_v}{\partial z})=S,$$

$$Q=[\rho ,\rho u,\rho v,\rho w,E{]}^T,F=[\rho u,\rho u^2+P,\rho uv,\rho uw,(E+P)u{]}^T,$$

$$G=[\rho v,\rho uv,\rho v^2+P,\rho vw,(E+P)v{]}^T,H=[\rho w,\rho uw,\rho vw,\rho w^2+P,(E+P)w{]}^T,$$

$$F_u=[0,{\tau }_{xx},{\tau }_{xy},{\tau }_{xz},b_x{]}^T,G_v=[0,{\tau }_{xy},{\tau }_{yy},{\tau }_{yz},b_y{]}^T,H_v=[0,{\tau }_{xz},{\tau }_{yz},{\tau }_{zz},b_z{]}^T,$$

$$b_x=u{\tau }_{xx}+v{\tau }_{xy}+w{\tau }_{xz},b_y=u{\tau }_{xy}+v{\tau }_{yy}+w{\tau }_{yz},b_z=u{\tau }_{xz}+v{\tau }_{yz}+w{\tau }_{zz},$$

$${\mathbf{F}}_f=\alpha {\int }_0^t({\mathbf{U}}_{ib}-\mathbf{U})dt+\beta ({\mathbf{U}}_{ib}-\mathbf{U}),$$

$${\mathbf{U}}_{ib}(s,t)={\int }_V\mathbf{u}(x,t){\delta }_h(\mathbf{X}(s,t)-\mathbf{x})d\mathbf{x},$$

$$\mathbf{f}(\mathbf{x},t)=-{\int }_{\Gamma }{\mathbf{F}}_f(s,t){\delta }_h(\mathbf{X}(s,t)-\mathbf{x})ds,$$

$${\delta }_h(x,y,z)=\frac{1}{h^3}\lambda (\frac{x}{h})\lambda (\frac{y}{h})\lambda (\frac{z}{h}).$$

$$\lambda (r)=\{\begin{array}{ll}\frac{1}{8}(3-2|r|+\sqrt{1+4|r|-4|r{|}^2}),& 0\le |r|<1\\ \frac{1}{8}(5-2|r|-\sqrt{-7+12|r|-4|r{|}^2}),& 1\le |r|<2\\ 0,& 2\le |r|.\end{array}$$

$$P_a=ϵ\exp [-\ln (2)\frac{(x-4{)}^2+y^2}{0.04}]$$

$$\mathrm{Re}={\rho }_f{U}_{0}D/\mu ,M=U_0/c,$$

$$\alpha (t)={\alpha }_0[1+sin(2\pi ft)],$$

$$P_a=A^{\ast }\exp [-\ln (2)\frac{(y-1{)}^2}{0.04}]+A^{\ast }\exp [-\ln (2)\frac{(y+1{)}^2}{0.04}],$$

$$A_p=-A\exp \{-B\log (2)[(x-2{)}^2+y^2+z^2]\}\cos (\omega t),$$

$${\mathbf{X}}_0(t)=\frac{A_0}{2}\cos (2\pi ft)[\cos \beta ,\sin \beta ],\alpha (t)=\frac{{\alpha }_m}{2}\sin (2\pi ft+\varphi ).$$

$$\frac{A_0}{L},\mathrm{Re}=\frac{{\rho }_{f}UL}{\mu },U_r=\frac{U_0}{U},M=\frac{\pi U}{2c},m^{\ast }=\frac{{\rho }_s}{{\rho }_{f}L},{\omega }^{\ast }=\frac{2\pi f}{{\omega }_n},$$

$$C_T=\frac{2F_T}{{\rho }_f{U}^{2}L},C_L=\frac{2F_L}{{\rho }_f{U}^{2}L},C_p=\frac{-2{\int }_0^L\mathbf{f}\cdot \mathbf{v}dl}{{\rho }_f{U}^{3}L},$$

Full text

K****eywords Immersed boundary method  $\cdot$ large deformations  $\cdot$ complex geometries  $\cdot$ fluid–structure–acoustics interaction

1 Introduction

Flow-induced pressure fluctuations extensively exist in engineering applications, such as aviation technologies [1], ventilation [2] and biomechanics [3]. In many of these applications, the major task is to identify the noise source and then to reduce or eliminate it (see e.g. rotorcrafts [1] and ventilations [2]). Flow-induced sound also plays a crucial role in the biomechanics, such as phonation (flow-induced sound in the larynx [4]) and heart murmurs (blood-flow-induced noise [5]). Accurate prediction of the flow-induced noise in these systems can be utilized to optimize low noise engineering design. The acoustic prediction in biomechanics also has potential in medical diagnosis [3]. However, in most of the engineering applications, fluid-structure-acoustics interactions involving complex geometries are involved, making the accurate prediction of the sound generation extremely challenging.

Numerical methods for computational aeroacoustics (CAA) can be categorized into three groups [6]. The first group of methods calculate the near field flow dynamics by using computational fluid dynamics (CFD) techniques, and predict the far-field sound generation by employing acoustic analogies, such as the pioneering analogy proposed by Lighthill [7] and the extension of this analogy to include the influence of the solid boundaries in the sound field [8]. Since Lighthill’s work, great efforts have been made to improve the method [9, 10, 11]. The second group of methods, also called acoustic/viscous splitting methods, decompose the compressible viscous equations into incompressible and perturbed compressible parts. The perturbation in the far-field represents the acoustic quantities [12]. The third group of methods employ direct numerical simulation (DNS) to solve the compressible Navier-Stokes equations. By using DNS, the generation and propagation processes of the sound in the near and intermediate fields can be calculated, without suffering from the limitations such as the low Mach number and compactness of the source region [6, 13].

In order to resolve the sound pressure which is much smaller compared to the ambient pressure [1], high-order finite difference methods (e.g. compact schemes [14] and Weighted/Targeted Essentially Non-oscillation schemes [15, 16]) on body-conformal meshes are widely used. In these methods, it is difficult to generate high-quality structured meshes around the complex and moving boundaries. The finite volume method has the ability to handle complex boundaries. However, it is still challenging to handle moving boundaries. In addition, it suffers from the low-order accuracy and introduces extra dissipation and dispersion [17]. The high-order overset grid method [18] has been successfully applied for CAA in complex geometries. But, the applications of this type of methods are limited due to their complexity of implementation and low efficiency at low Mach numbers [3]. As an efficient method for fluid–structure interaction (FSI), the immersed boundary (IB) method, first developed by Peskin [19, 20], has been extended to acoustics [3, 17]. In this method, structured mesh is generated initially and fixed during the computation. Therefore, the regeneration of mesh is avoided. Because of its simple boundary treatment, IB method has gained popularity for a wide range of applications [21, 22, 23, 24]. In the hybrid method proposed by Seo and Mittal [3], the IB method and the hydrodynamic/acoustic splitting technique are combined to handle the acoustic problems involving complex geometries in low-Mach number flow. It encounters challenges when the Mach number rises and the complex geometries move. Penalty IB (pIB) is a typical IB method [25], where the IB is conceptually split into two Lagrangian components: one component is massless and interacts with the fluid exactly as the traditional IB method, and the other component carrying mass is connected to the massless component by virtual stiff springs. However, most of the previous studies based on pIB method focus on FSI problems without acoustics [21, 26, 27, 28, 29, 30]. Sun et al. presented an IB method which considers the linear Euler equations for acoustic scattering modeling [17]. In their work, the moving boundary and fluid viscosity are not explored. The ability of the pIB method for fluid–structure–acoustics interactions involving complex geometries has not been well explored. This is the motivation for this work.

Flapping foil has drawn growing attention recently, due to its high performance in micro aerial vehicles (MAVs) [31, 32, 33, 34, 35] and power generators [36, 37, 38, 39]. The aerodynamic characteristics of the flapping foil including its thrust, lift and power efficiency have been extensively studied by researchers [31, 40]. However, the aeroacoustics induced by the flapping foil has not been well understood [41]. Although the high performance of insect flight can be achieved with low noises, the flapping sound generated by the foil may have potential effects on biological functions, such as the communication using aposematic signals with the locomotion-induced sound [41]. Moreover, the study on the sound induced by flapping foil may also have potential applications in medical inspections, such as the phonation of human which performs like a flapping foil induced sound [4]. Due to its importance in many areas, the sound generation by flapping foils is numerically studied by the present method. The effects of both the elasticity and the geometrical shape of the foil on the force and sound generation process are analyzed. The numerical examples presented here can also enrich the limited database of fluid–structure–acoustics interactions.

In this paper, an immersed boundary method introduced in our previous work [30] is extended to fluid–structure–acoustics interactions involving large deformations and complex geometries. The organization of the paper is as follows. The numerical approach is briefly introduced in Section 2. Several validations including acoustic waves scattered by a stationary cylinder, sound generation by a stationary and a rotating cylinder, sound generation by an insect in hovering flight, deformation of a red blood cell induced by acoustic waves and acoustic waves scattered by a stationary sphere are presented in Section 3. Application of the current numerical method in modelling flapping foil induced acoustics is presented in Section 4. Finally, conclusions are given in Section 5.

2 Numerical method

The current numerical method includes three important components: the structure, compressible fluid and fluid–structure interaction. Without loss of generality, a flexible plate immersed in the two-dimensional fluid is used as an example to introduce the structure dynamics. The plate is assumed to be elastic and its dynamics is governed by the following nonlinear equation [22, 42, 30]

[TABLE]

where $\boldsymbol{X}$ is the Lagrangian coordinates of the flexible beam, $\rho_{s}$ is linear density, $K_{S}$ and $K_{B}$ are respectively the stretching and bending rigidity, $s$ is the arc coordinate, and $\boldsymbol{F}_{f}$ is the external force acting on the beam. The absolute nodal coordinate formulation (ANCF) proposed by Shabana [43, 44, 45] is adopted to solve Eq. 1. This method was combined with the IB method by Wang et al. [30].

The flow dynamics considered here are governed by the compressible viscous Navier–Stokes equations

[TABLE]

where $\rho$ is the fluid density, $u$, $v$ and $w$ are respectively the three velocity components, $P$ is the pressure, $E$ is the total energy, $S$ is a general source term including the IB-imposed Eulerian force and other body forces, Re is the Reynolds number, and $\tau_{ij}$ is the shear stress.

In the fluid solver, the fifth-order Weighted Essentially Non-oscillation (WENO) scheme proposed by Liu et al. [15] is used for the spatial discretization of the convective term. We also introduce a recently developed scheme, the Target Essentially Non-oscillation (TENO) [16], to discretize the convective term, to compare the performance of the IB–WENO with that of the IB–TENO. For the viscous terms, a fourth-order central difference scheme is used to discretize the spatial derivatives. For all unsteady equations involved in flow solver, the third-order TVD Runge-Kutta method is used for temporal discretization [46]. The dynamics of the fluid and flexible structures are solved independently. The interaction force is calculated explicitly using a feedback law [47] based on the pIB method [25]. The interaction force between the fluid and the structure can be determined by the feedback law [25]

[TABLE]

where $\boldsymbol{U}_{ib}$ is the boundary velocity obtained by interpolation at the IB, $\boldsymbol{U}$ is the structure velocity, and $\alpha$ and $\beta$ are large positive constants. It is noted that this method does not transform the Lagrangian density of structure into the Eulerian density, as did in Zhu and Peskin [48]. Therefore, it is not necessary to apply the incompressible limitation or to modify the continuity equation [48, 30]. Eq. 7 is the Lagrangian force acting on the structure. The Lagrangian force acting on the fluid by the immersed boundary is -$\boldsymbol{F}_{f}$, which is spread onto fluid nodes to achieve the boundary condition. The interpolation of the velocity and the spreading of the Lagrange force to the adjacent grid points are expressed as

[TABLE]

[TABLE]

where $\boldsymbol{u}$ is the fluid velocity, $\boldsymbol{X}$ is the coordinates of structural nodes, $\boldsymbol{x}$ is the coordinates of fluid, $s$ is the arc coordinate for a two dimensional domain, $V$ is the fluid domain, $\Gamma$ is the structure domain and $\delta_{h}$ is the smoothed Dirac delta function [20], expressed as

[TABLE]

In this paper, the four-point delta function introduced by Peskin [20] is used

[TABLE]

Instead of using uniform mesh which is time consuming, a non-uniform mesh is used here to improve the computational efficiency [23]. The mesh is uniform in both directions within a small inner box containing the solid for achieving good accuracy and interpolation in the IB method, and it is stretched in the remaining of the computational domain. Non-reflecting boundary conditions are applied on the boundaries [49, 50]. In addition, the solver is parallelized using hybrid OpenMP and MPI [51].

3 Validations

Validations of the present solver for fluid–structure interactions including flow over a stationary cylinder, structure dynamics, deformation of a flexible panel induced by shock waves in a shock tube, an inclined flexible plate in a hypersonic flow, and shock-induced collapse of a cylindrical helium cavity in the air have been conducted in our previous work [30]. Here, we focus on validations of acoustics modelling including acoustic waves scattered from a stationary cylinder in a quiescent flow, sound generation by a stationary and a rotating cylinder in a uniform flow, sound generation by an insect in hovering flight, deformation of a red blood cell induced by acoustic waves and acoustic waves scattered by a stationary sphere.

3.1 Acoustic waves scattered by a stationary cylinder

In order to verify the accuracy of the present numerical method on uniform Cartesian meshes in simulating the acoustic problem, acoustic waves scattered by a stationary cylinder is validated against the published data [52, 53]. In this simulation, a rigid cylinder with a diameter of $D=1.0$ is fixed with its center at (0, 0). The initial conditions are localized pressure perturbations of Gaussian distribution [52]

[TABLE]

where $\epsilon=10^{-3}$. The initial density and pressure are respectively 1.0 and $1/\gamma$ ($\gamma=1.4$), and the initial velocity of the fluid is zero.

We used a uniform Cartesian mesh over a computational domain of ($20D\times 20D$) with three different mesh spacings: $D/100$, $D/50$ and $D/25$. The time histories of the fluctuating pressure, $\Delta p=P-p_{\infty}$ where $P$ and $p_{\infty}$ are respectively the current and ambient pressure, at points (2, 0) and (2, 2) are plotted in Fig. 1 along with the DNS results reported in Ref. [52]. According to the comparison in Fig. 1, all the three mesh spacings are able to capture the fluctuating pressure. The small discrepancies depicted in the scaled up figures (see insets in Fig. 1) can be diminished when the mesh size is $D/100$.

Fig. 2 presents the comparison of present results with those calculated by the immersed boundary–lattice Boltzmann method (IB-LBM) in Ref. [54]. The mesh spacing for the fluid is $D/40$, same as that used in Ref. [54]. The results show that the present IB–WENO scheme has an excellent performance on capturing the scattered acoustic pressure. Specifically, the sound pressure predicted by the current explicit IB is significantly better than that predicted by the explicit IB–LBM, as reported in Ref. [54].

We further present the fluctuating pressure contour for mesh spacing of $D/50$ at three instants in Fig. 3 along with figures from Ref. [53] to illustrate the propagation of the acoustic wave. As shown in Fig. 3(a), a principal pulse is generated due to the initial pressure perturbation. When the acoustic wave impacts on the cylinder, a reflected wave off from the cylinder surface is generated and yields to a secondary acoustic wave, as demonstrated by Fig. 3(b). As reported by Liu and Vasilyes [52], two parts of the principal wave front split by the cylinder traverse its span, collide, and merge, thereby generating a third acoustic wave front (see Fig. 3(c)). These three acoustic wave fronts are all well captured according to the time history of fluctuating pressure presented in Fig. 1. This benchmark validation shows that the present method has a good accuracy in capturing acoustic waves.

3.2 Sound generation by a stationary cylinder in a uniform flow

Flow around a stationary cylinder has been extensively studied theoretically, experimentally and numerically [55, 56, 57], since the work of Strouhal on aeolian tones. Most of previous studies on the sound generation due to flow past a circular cylinder was done by using hybrid or acoustic/viscous splitting methods to reduce the computational expense [6]. However, direct numerical simulation (DNS) is an effective way to identify both the aerodynamics and characteristic features of the sound accurately. Inoue and Hatakeyama studied the sound generation by a stationary [6] and rotating [58] cylinder in a uniform flow using DNS, and clarified the relation between the vortex/flow dynamics and the sound pressure.

In this section, a stationary cylinder in a uniform flow is considered. The cylinder is located at the origin. In order to discuss the sound pressure, we define $r$ (nondimensionalized by the diameter of the cylinder) as the distance from the origin and $\theta$ as the circumferential angle. The non-dimensional parameters that govern this problem are defined as follows

[TABLE]

where $D$ is the diameter of the cylinder, $U_{0}$ is the inlet velocity, and $\rho_{f}$, $\mu$ and $c$ are respectively the density, viscosity and sound speed of the fluid in the far field. In the current stationary case, Re=150 and $M=0.2$. Three mesh regions are used in the simulations to improve the computational efficiency, including a cylinder occupied region, a sound region and a sponge region (similar as that in Ref. [6]) around the cylinder. The cylinder occupied region extends from ($-1.25D,-1.25D$) to ($1.25D,1.25D$). The sound region is $212D$ in width and length. The sponge region is as large as $400D$ in both width and length to diminish the reflections from the boundary. Non-reflecting boundary are applied on the external boundaries of the sponge region. In the cylinder occupied region, the mesh spacing is $D/40$. The mesh spacing in the sound region increases to $D/5$. There is a sinusoidal transition between the cylinder occupied region and sound region. The maximum mesh spacing in the sponge region is $12.5D$. In order to achieve a fast transition to the asymmetric Kármán vortex street, an initial velocity perturbation is applied in the near wake.

The time histories of the drag ($C_{D}$) and lift ($C_{L}$) coefficients scaled by the $0.5\rho_{f}U_{0}^{2}D$ are presented in Fig. 4, with data reported in Ref. [6] for comparison. The amplitudes of $C_{L}$ and mean value of $C_{D}$ are $0.525$ and $1.40$, respectively, which agree well with the results ($C_{L}=0.520$ and $C_{D,m}=1.39$) from Ref. [6]. The Strouhal number calculated by the present method is about 0.178, which is also close to 0.175 in Ref. [59] and 0.183 in Ref. [6].

The fluctuating pressure $\Delta\tilde{p}$ is defined by $\Delta\tilde{p}(x,y,t)=\Delta p(x,y,t)-\Delta\bar{p}(x,y)$, where $\Delta p$ denotes the total fluctuating pressure and $\bar{p}$ is the time average pressure. $\Delta p=p-p_{\infty}$ with $p_{\infty}$ being the ambient pressure [6]. In present paper, the fluid density $\rho_{f}$ and sound speed $c$ are used to nondimensionalize the pressure. Fig. 5 shows decay of the pressure peaks (include both positive and negative peaks) measured at $\theta=90^{o}$, which agrees well with results from Ref. [6]. The pressure peaks of the cylindric pressure waves generated by the unsteady flow tend to decay in proportion to $r^{-\frac{1}{2}}$, which agrees well with the theoretical prediction by Landau and Lifshitz [60], as indicated by the dashed line in Fig. 5.

We further increase the Reynolds number to 1000 (Ma=0.2), and assume it is still in laminar flow regime. The same problem has also been examined by Brentner et al. [61], who used the Lighthill acoustic analogy to separate flow dynamics and acoustics calculations. Fig. 6 presents the comparison of the time histories of $C_{D}$ and $C_{L}$ with the data from Ref. [61]. The present drag coefficient and St number are respectively 1.60 and 0.215, which agree with the numerical results (1.56 and 0.238) in Ref. [61]. The discrepancy (2.6% in the drag coefficient and 9.7% in St) is probably due to the compressibility of the fluid neglected in the reference, as indicated in Ref. [61] where the compressible solver predicts lower frequency than the incompressible one.

Fig. 7 shows the SPL (reference pressure $20uPa$) in the frequency space obtained by the Fast Fourier transform (FFT). For comparison, data from in Ref [61, 62] are shown in the same figure. The present results agree well with those from the references. The maximum discrepancies of the peak values are about 5%, which is probably induced by the neglected quadrupole source in FW-H. It should be noted that the results from Ref. [61] were calculated in two steps: the incompressible N-S equations were first solved, then the acoustics was computed by solving Ffowcs Williams–Hawkings (FW-H) equation using the resolved unsteady flow as the input with the quadrupole source neglected, while, the present method obtains the sound pressure from DNS. The experimental results from Ref. [62] show much lower peaks at high Reynolds number, but the frequency space still coincides with that at low Reynolds number. According to Ref. [61], it is still challenging to accurately predict the SPL in the frequency space by using turbulence models such as $k-\omega$ and SST. In the future, our effort will be made to incorporate large eddy simulation and wall models into the current solver to address this challenge. Fig. 8 is a snapshot of the fluctuating pressure contours for future validation of newly developed method.

3.3 Sound generation by a rotating cylinder in a uniform flow

Here, we consider the sound generated by a rotating cylinder in a uniform flow, which was studied by Inoue et al. [58]. The non-dimensional parameters governing this problem are Re, $M$ (defined by Eq. 13) and the velocity ratio, defined by $\alpha=M_{\theta}/M$, where $M_{\theta}$ is the anticlockwise angular velocity of the rotating cylinder. In the current rotating case, Re=160, $M=0.2$ and $\alpha=1.5$. The mesh strategy is same as that in Section 3.2.

The decay of the fluctuating pressure peaks measured at $\theta=90^{o}$ is plotted in Fig. 9. The pressure peaks also tend to decay in proportion to $r^{-\frac{1}{2}}$, similar to that discussed in the stationary case. The polar plots in Fig. 10 are the root-mean-square of the fluctuating pressure $\Delta\tilde{p}$ measured at $r=r^{\prime}(1-M{\rm cos}\theta)$ with $r^{\prime}=75$, taking the Doppler effect into account. The data from Ref. [58] are also presented in Fig. 10 for comparison. The results show that the pressure measured at $\theta=90^{o}$ is larger than that measured at $\theta=0^{o}$, indicating the lift dipole dominates the sound generation. The profile is not symmetric with respect to the line $\theta=0^{o}$ and $180^{o}$ as that in the stationary condition, due to the effects of rotation. The good agreement between the present results and the computational results from Ref. [58] shows that the present method has an excellent ability to handle the acoustic simulations involving moving boundaries.

Fig. 11 shows the time histories of fluctuating pressure measured at $r=75$, $\theta=90^{o}$ and $270^{o}$. The results show that the pressure peaks measured at $\theta=90^{o}$ are higher than those at $\theta=270^{o}$, due to the fact that the anticlockwise rotation of the cylinder leads to asymmetrical lift and drag. This is also indicated by the root-mean-square of $\Delta\tilde{p}$ plotted in Fig. 10.

3.4 Sound generation by an insect in hovering flight

In this section, an insect in hovering flight is considered to validate the current solver in handling moving body with relatively complex geometry configuration. The schematic of the problem is shown in Fig. 12. A circular cylinder (body) with a diameter of $0.5L$, and two elliptic cylinders (wings) with the dimensions of $L$ and $0.4L$ are used to model the insect. The two wings flap symmetrically, and the flapping motion of the wings is prescribed by

[TABLE]

where the amplitude of the flapping angle $\alpha_{0}$ is $25.3^{o}$ and $f$ is the flapping frequency. The Reynolds number defined by $\rho_{f}U_{max}L/\mu$ is 200, where $U_{max}$ is the maximum wing tip velocity. The Strouhal number defined by $fL/U_{max}$ is 0.25, and the Mach number based on $U_{max}$ is 0.1. The computational domain extends from ($-100L$, $-100L$) to ($100L$, $100L$), with 50 points to resolve the wing length.

A direct comparison of the lifts generated by the body and wing is shown in Fig. 13. It is found that the present results agree well with the data from Ref. [63]. Fig. 14 shows that the fluctuating pressure measured at (0, $60L$) and (0, $-60L$) from present computation is also in a good agreement with the data from Ref. [63]. Good agreements of the root-mean-square of fluctuating pressure measured at a distance of $50L$ shown in Fig. 15 confirm the good performance of the current method in handling acoustics involving complex moving geometries. Additionally, a qualitatively comparison of the fluctuating pressure is presented in Fig. 16.

3.5 Deformation of a red blood cell induced by acoustic waves

In this section, we consider the deformation of a red blood cell (RBC) induced by acoustic waves to validate the current method in handling fluid–structure–acoustics interactions. This problem was experimentally studied by Mishra et al. [64]. In the simulation, a localized pressure perturbation of Gaussian distribution is applied on the fluid domain, which is expressed as

[TABLE]

where $A^{*}=A/(\rho_{f}c^{2})$ is the non-dimensional perturbation amplitude, and $c$ is the sound speed of the unperturbed fluid. The structure-to-fluid mass ratio is $m^{*}=\rho_{s}/(\rho_{f}D)=0.2$, where $D$ is the diameter of the cylindrical cell. The RBC is assumed to be elastic, and its non-dimensional stretching and bending rigidity are $K_{S}^{*}=K_{S}/(\rho_{f}c^{2}D)=0.1$ and $E_{B}^{*}=E_{B}/(\rho_{f}c^{2}D^{3})=1.0\times 10^{-6}$, respectively. The prestress is also applied on the RBC to hold it in a cylindrical shape initially. The pressure perturbation amplitudes ranging from 0.2 to 1.0 are examined. The instantaneous sound pressure contours at three instants for $A^{*}=0.4$ are shown in Fig. 17. Fig. 18 shows a qualitative comparison of the instantaneous deformation of the RBC. The results show that the shape of the RBC changes from a cylinder to an ellipse under the load of the symmetrically distributed pressure perturbations. As shown in Fig. 18, the deformation increases with the pressure amplitude, which qualitatively agrees with the experimental results. In addition, the wavy deformations of the RBCs are observed for large $A^{*}$ (e.g. $A^{*}=1.0$). This phenomenon is the dynamic response of the RBC to the external load, and is consistent with the experimental observation.

3.6 Acoustic waves scattered by a stationary sphere

In this section, acoustic waves scattered by a stationary sphere is conducted to validate the present solver in capturing acoustic waves in three-dimensional space. In this case, a rigid sphere with a radius of $R$ is fixed in the fluid with its center at (0, 0). A periodic Gaussian pressure source in the fluid domain is applied

[TABLE]

where $A=0.01$, $B=16$ and $\omega=2\pi$. All quantities are non-dimensionalized by the radius of the sphere, the ambient fluid density and sound speed.

Three mesh regions are used to improve the computational efficiency and conserve the accuracy as used in Section 3.2. Three mesh spacings of $R/20$, $R/40$ and $R/80$ within a domain of $2R\times 2R\times 2R$ covering the sphere are used. The far field mesh spacing is $R/10$. The sphere is discretized by triangular elements with a mesh spacing approximately equal to that of the fluid domain. The fluctuating pressure along $x$-direction on the line of $y=0,z=0$ and along $y$-direction on the line of $x=0,z=0$ are presented in Fig. 19 with data from Refs. [65, 66]. As shown in this figure, the fluctuating pressure from the present simulation approaches to the results from the references when the mesh is refined. Good agreement is achieved when the finest mesh spacing of $R/80$ is used, showing that the present numerical method is able to predict acoustic waves in three-dimensional space. Fig. 20 shows three views of the fluctuating pressure field at the mesh spacing of $R/80$ when the periodic solution is obtained.

To compare the performance of the IB–WENO and the IB–TENO, Fig. 21 shows the pressure fluctuation along $x$-direction. It is found that the numerical dissipation of the IB-TENO is lower compared to the IB-WENO, as demonstrated by the peak values near $x/R=3.1$ and 4.7.

4 Sound generation by flapping foils in a uniform flow

Having conducted rigorous validations, the sound generated by flapping foils in forward flight and flapping foil energy harvester is numerically studied in the two-dimensional domain by considering the geometrical shape and the flexibility of the foil. Both the force and sound generations are compared and analyzed.

4.1 Flapping foils in forward flight

4.1.1 Physical problem

In this section, the acoustic perturbations induced by flapping foils in forward flight are considered, as shown in Fig. 22 [31]. Here, the foil is clamped at the leading edge, and the clamping device undergoes a combined translational and rotational motions, given by [31]

[TABLE]

where $A_{0}$ is the translational amplitude, $f$ is the flapping frequency, $\alpha_{m}$ is the rotation amplitude, and $\beta$ is the angle between the stroke plane and the horizontal plane. The phase difference ($\phi$) is 0 unless it is mentioned.

The non-dimensional parameters including the flapping amplitude, Reynolds number, inlet velocity, Mach number, mass ratio and frequency ratio that control this problem are given by

[TABLE]

where $L$ is the chord length, $U=2fA_{0}$ is the average translational velocity of the leading edge, $f$ is the flapping frequency of the foil, $c$ is the sound speed of the fluid, $\rho_{f}$ is the fluid density, $\rho_{s}$ is the linear density of the foil, and $\omega_{n}=k_{n}^{2}/L^{2}\sqrt{E_{B}/\rho_{s}}$ with $k_{n}=1.8751$ (the frequency of the first natural vibration mode of the wing with fixed leading edge [67, 31]) and $E_{B}$ being the bending rigidity. Here, Re=100, $m*=5.0$ and $U_{r}=0.4$. The thrust, lift and power coefficients are defined as

[TABLE]

where $F_{T}$ and $F_{L}$ are respectively the thrust and lift force acting on the foil by the ambient fluid, $\boldsymbol{f}$ is the hydrodynamic traction on the foil, and $\boldsymbol{v}$ is the velocity of the foil. The fluctuating pressure $\Delta\tilde{p}$ is defined in the same way in Section 3.2. All the pressures presented afterwards are scaled by the fluid density $\rho_{f}$ and sound speed $c$.

4.1.2 Mesh convergence study

First, mesh convergence is conducted to guarantee the reliability of the solver in modeling the acoustic perturbations induced by flapping foils. Non-uniform mesh scheme is adopted to improve the computational efficiency, where a uniform refined mesh region around the foil is used to enhance the accuracy. In the simulation, the computational domain extends from ($-36.25L$, $-36.25L$) to ($36.25L$, $36.25L$). Three cases are considered for the mesh convergence study: $\Delta x=L/80$, $L/40$ and $L/20$, where $\Delta x$ is the finest mesh spacing around the foil. Other non-dimensional parameters for these three cases are: $M=0.1$, $A_{0}/L=1.25$, $\beta=45^{o}$ and $\alpha_{m}=0^{o}$. In Fig. 23, the time histories of thrust, lift and power coefficients are presented. The results show that reliable results can be achieved with $\Delta x=L/40$.

Similar to Section 3.2, we define $r$ as the distance from the origin ( nondimensionalized by $L$) and $\theta$ as the circumferential angle. The histories of the fluctuating pressure measured at $r=10,\theta=90^{o}$ are presented in Fig. 24, where it is found that the effects of the mesh difference are negligible for mesh size of $L/80$ and $L/40$. Therefore, $\Delta x=L/40$ is adopted in the later simulations. This figure also shows that the positive pressure wave peaks are significantly larger than the negative ones, which will be discussed in detail later.

4.1.3 Propagation and decay of pressure waves

In Fig. 25, the instantaneous contours of the fluctuating pressure in a flapping period are presented. It is observed that during downstroke, negative pressure pulses are generated from the upper side of the foil, and positive pressure pulses are generated from the lower side of the foil, and vice versa. During a foil stroke, positive pressure is formed on the loading face due to fluid being displaced and negative pressure is formed on the opposite face mainly due to the leading-edge vortices. The positive and negative pressure exchange sides during the switch between the downstroke and the upstroke [41]. A clear illustration of the foil-loading mechanism can be found in the work by Inada et al. [68].

The circumferential distribution of fluctuating pressure peaks are presented in Fig. 26. It shows that the pressure peaks decay with the distance. In addition, the positive peaks are much larger than the negative ones, because as mentioned above, the positive and negative pressure fluctuations are generated by different mechanisms. Fig. 26 shows that the pressure peaks decay with the increasing distance, and the positive peaks are much larger than the negative peaks. In addition, Fig. 26 shows that the pressure peaks distribution is symmetrical about the stroke plane, indicated by the dashed line in the figure. It is also noted that some markers plotted in Fig. 26 are not smoothly distributed along the circumference. A reasonable explanation is that the oscillation is introduced by the vortices. When the vortex approaches to the probes, the recorded pressure will be perturbed. As the distance increases, the fluctuating pressure along the circumference tend to be smoother, as shown in Fig. 26. This is due to the fact that the vortex effects reduce with the distance. The fluctuating pressure peaks generated by the flapping wing on the windward side ($180^{o}-270^{o}$) are slightly larger than those on the leeward side ($0^{o}-90^{o}$), as shown in Fig. 26, due to the presence of the free stream. In order to figure out the influence of the phase difference ($\phi$) between the translational and rotational motions on the fluctuating pressure, combined translational and rotational motions of the flapping wing at $\phi=0^{o},45^{o}$ and $90^{o}$ with $\alpha_{m}=\pi/4$ are simulated. Polar diagrams of the fluctuating pressure peaks are presented in Fig. 27. It is found that the fluctuating pressure peaks on the windward side decrease slightly with the phase difference. Similar to the case without rotational motion, the positive fluctuating pressure peaks on the windward side are larger than the negative ones, due to the presence of the free stream. On the leeward side, the difference between the positive and negative fluctuating pressure peaks are not significant. However, the comparison between the flapping with and without translational motion, as shown in Figs. 26 and 27, shows that the combined translational and rotational motions can reduce the fluctuating pressure peaks, especially those on the leeward side.

Shown in Fig. 28 are the polar plots of the fluctuating pressures $\Delta\tilde{p}$ at different instants from $0/6T$ to $3T/6$ in the 31st cycle. As the upstroke process is a inverse of the downstroke process, only a half period is presented. It is found that the directivity of the pressure waves fluctuates around the stroke angle $45^{o}$ at the flapping frequency. Fig. 28 also describes the whole sound generation process due to the flapping wing loading and the vortex shedding.

After the analysis of the sound generation mechanism and amplitude, we would like to further analyze the frequency of the flapping foil-induced sound by using FFT to analyze the fluctuating pressure. The fluctuating pressure of ten cycles from $10T$ to $20T$ measured at the probes along $r=18$ are used in the frequency analysis. Fig. 29 shows the polar distribution of the fluctuating pressures at the frequency of $f$ (the flapping frequency), $2f$ and $3f$. It indicates that the sound pressure is dominated by the frequency of $f$. The peaks of the fluctuating pressures at the frequencies of $2f$ and $3f$ are much lower than those at the frequency of $f$. The fluctuating pressure at the frequencies over $3f$ are negligible.

4.1.4 Comparison of the sound generated by flat plate and NACA0015 foil

Here, two flat plates (a rigid and a flexible) and a rigid NACA0015 foil undergoing combined translational and rotational motion in a uniform flow, which widely exists in the insects wings [31], are numerically simulated to study the influences of the foil flexibility and geometry.

The motion of the foil is controlled by Eq. 17. All the non-dimensional parameters are defined in Eq. 18. The amplitude of the rotating angle $\alpha_{m}=\pi/4$, other parameters for the flexible foils are: frequency ratio $\omega^{*}=0.6$ and $m^{*}=5.0$. The thrust, lift and power coefficients defined in Eq. 19 are presented in Fig. 30, along with the numerical data from Ref. [31]. It is found the results from the current compressible solver agree well with the data obtained by Tian et al. [31] using an incompressible solver. Because the Mach number used in the current solver is low ($M=0.1$), the compressibility of the fluid is not significant. As seen from Fig 30, the discrepancies between the current results and those from Ref. [31] for the flexible foil seem to be much more significant than those for the rigid foil. A plausible explanation is that the flexible foil induces higher fluid velocity around it, which makes the compressibility more significant. The histories of the thrust, lift and power coefficients in a period for the three cases are plotted in Fig. 31 for comparison. The results show that the flexibility of the plate has a remarkable influence on the force generation. The thrust and lift are enhanced significantly by the flexibility. However, the comparison between the rigid plate and NACA0015 foil indicates that the geometrical shape does not have significant effects on the aerodynamic characteristics at the conditions considered.

In order to compare the sound generation of rigid and flexible foils, the polar diagrams of the root-mean-square of $\Delta\tilde{p}$ measured at $r=30$ are presented in Fig. 32. As shown in Fig. 32 (a) and (c), the fluctuating pressure generated by the rigid plate and NACA0015 foil distributes symmetrically about the stroke plane (indicated by the dashed line in Fig. 32 (a) and (c)). However, Fig. 32 (b) shows that the distribution of the fluctuating pressure generated by the flexible foil is asymmetrical. An evident shift (about $15^{o}$ anticlockwise) is observed, as shown in Fig. 32 (b). As analyzed in Ref. [31], the deformation of the flexible plate during upstroke is higher than that during downstroke, due to the presence of free stream. Obviously, this sound shift is introduced by the asymmetrical deformation (as shown in Fig. 33) of the foil. Fig. 32 also shows that the flexible plate generates the largest fluctuating pressure, which is significantly larger than that generated by the rigid plate and NACA0015 foil. The results indicate that the flexibility influences the sound generation significantly. The effects of geometrical shape on the sound amplitudes are much smaller. In addition, it is found that the fluctuating pressure on the lower section is larger than that on the upper section. A reasonable explanation is that the lower section is located on the windward side and the upper section is on the leeward side.

The FFT is further used to study the frequency of the fluctuating pressure. According to the frequency analysis of the fluctuating pressure peaks, three main frequencies ($f$, $2f$ and $3f$) are captured in these three flapping foils. The fluctuating pressures at higher frequencies are negligible. The polar diagrams of the fluctuating pressure peaks at the three main frequencies are presented in Fig. 34. It is found that the frequency of $f$ dominates in the rigid flapping foil. However, Fig. 34 (b) shows that both the frequency of $f$ and $2f$ dominate in the flexible flapping foil. As shown in Fig. 34 (a) and (c), the fluctuating pressure for rigid foils distributes symmetrically to the stroke plane as mentioned before. Fig. 34 (b) shows that the fluctuating pressure at the frequency of $f$ generated by the flexible foil shift about $15^{o}$ anticlockwise compared with that generated by the rigid foil. However, the shift of the fluctuating pressure at the frequency of $2f$ approaches about $30^{o}$. This indicates that the propagation direction of the fluctuating pressure from a flexible foil at the frequency of $f$ and $2f$ exhibit remarkable differences, which can be used to identify the characteristics of the sound source. Additionally, the comparison between the rigid plate and NACA0015 foil also shows that the geometrical effects on the sound generation is not significant as mentioned in Section 4.1.4.

4.2 Flapping foil energy harvester

In this section, we further adopt the current numerical method to study the sound generated by a thin plate (without considering the thickness of the plate) and a NACA0015 foil (stream-lined shape) flapping as energy harvesters.

A foil undergoing prescribed plunge and pitch governed by Eq. 17 in a uniform flow is considered. The characteristic parameters defined in Eq. 18 are given as: $A_{0}/L=2.0$, $M=0.1$, $U_{r}=1.78$, $\alpha_{m}=152.6^{o}$ and $\beta=90^{o}$. The Reynolds number, defined by Re=$\rho_{f}U_{0}L/\mu$ is 1100, and the distance from leading edge to the pitching axis is $L/3$. The details of the setup can be found in Refs. [69, 70]. The computational domain based on a non-uniform Cartesian mesh has a size of $96L\times 96L$, with a uniform area extending from ($-2.5L$, $-2.5L$) to ($2.5L$, $2.5L$), where the mesh spacing is $L/100$. Extensive preliminary tests have been conducted to make sure the results are independent of the mesh size and computational domain. The time histories of non-dimensional drag ($C_{D}$, defined as $-C_{T}$), lift ($C_{L}$) and power ($C_{p}$) extracted from the fluid obtained from the present simulations agree well with those available in Ref. [69] as shown in Fig. 35.

Rigid and flexible plates without considering the thickness effect are also numerically simulated for comparison. All the non-dimensional parameters for the rigid plate are same as that for the NACA0015 foil. For the flexible plate, the mass ratio defined in Eq. 18 is 5.0, the non-dimensional bending rigidity of the flexible plate defined as $E^{*}=E_{B}/(\rho_{f}U_{0}^{2}L^{3})$ is $1.0$. Fig. 36 presents the comparison of the drag, lift and power coefficients in a period obtained from the NACA0015 foil and the two plates. It is found that the tendencies of $C_{D}$, $C_{L}$ and $C_{p}$ are consistent, though small discrepancies are observed. Specially, Fig. 36 (a) shows that the plates induce larger drag than NACA0015 foil, and the flexibility of the plate also increases the drag coefficient. Fig. 36 (b) shows that the lift generation of the rigid plate and NACA0015 is very similar, the flexibility tends to generate lower lift amplitude. The mean power coefficient (defined by $-\bar{C}_{p}$) is generally used to evaluate the efficiency of the foil in energy harvester. Here, the power coefficients in five periods are used to calculate $-\bar{C}_{p}$. The mean power coefficient of the NACA0015 foil and rigid plate are 0.85 and 0.82, respectively, compared with that of only 0.58 for the flexible plate. The results show that NACA0015 foil is the most efficient in energy harvesting among the three foils. The flexibility of the foil in this case deteriorates the energy harvesting efficiency. The instantaneous vorticity contours in a period are presented in Fig. 37. Similar vortex fields can be observed for the three cases, without remarkable differences.

Fig. 38 presents the decay of the fluctuating pressure (defined in Section 4.1) measured at $\theta=90^{o}$. It shows that the fluctuating pressure decays in proportion to the $r^{-\frac{1}{2}}$ in the intermediate and far fields, which agrees well with the theoretical result [60]. In all cases, the positive fluctuating pressure is slightly larger than the negative one, which may indicate the loading process is stronger than the vortices shedding effects on the fluctuating pressure generation at the measured point. The comparison of Fig. 38 (b) and (c) shows that the flexibility of the plate increases the positive pressure at the measured point, which can be reasonably explained by the larger tip-displacement of the flexible plate.

To illustrate the frequency characteristics of the fluctuating pressure generated by the foil, the FFT is used to analyze the fluctuating pressure. Polar diagrams of the fluctuating pressure peaks measured at $r=20$ with three frequencies $f$, $2f$ and $3f$ are shown in Fig. 39. It is found that the sound at frequency of $f$ dominates in the vertical direction for all cases, while the fluctuating pressure at frequency of $2f$ dominates in the horizontal direction, and the fluctuating pressures at higher frequency are negligible. The fluctuating pressure in the vertical direction is dominated by the frequency of the lift ($f$, see Fig. 39 (b)). Similarly, the fluctuating pressure in the horizontal direction is dominated by the frequency of the drag ($2f$, see Fig. 39 (a)). Comparison of the fluctuating pressure generated by NACA0015 foil and rigid plate shows negligible differences indicating that the geometrical shape of the foil does not significantly affect the sound generation. It is also found that the fluctuating pressure at the frequency of $2f$ induced by the flexible plate is larger than that induced by NACA0015 foil and rigid plate, which agrees with the larger amplitudes of drag induced by the flexible plates, see Fig. 36 (a). The smaller fluctuating pressure at the frequency of $f$ and $3f$ induced by the flexible plate also agrees with the smaller lift amplitude generated by the flexible plate, see Fig. 36 (b). The instantaneous contours of the fluctuating pressure $\Delta p$ in a period are further presented in Fig. 40. It is found that the difference of the fluctuating pressure field generated by the NACA0015 foil and rigid plate is not significant. The one of the flexible plate tends to be different at $0T/4$ and $2T/4$, when the deformation of the plate is significant.

Here, three additional cases at $E_{b}^{*}=0.25,$ 0.5 and 5.0 are considered to illustrate the effects of the flexibility on the sound generation. Fig. 41 shows the polar diagram of the root-mean-square of the fluctuating pressure. It shows that the fluctuating pressure generated by the flexible plates increase slightly with the flexibility. However, the flexible plate at $E_{b}^{*}=0.5$ generates a significantly larger fluctuating pressure, as the plate flaps near its natural frequency ($\omega^{*}$ defined in Eq. 18 is 0.8).

5 Conclusions

In this paper, an immersed boundary method for fluid–structure–acoustics interactions involving large deformations and complex geometries is presented. The validations including acoustic waves scattered from a stationary cylinder, sound generation by a stationary and a rotating cylinder in a uniform flow, sound generation by an insect in hovering flight, deformation of a red blood cell induced by acoustic waves and acoustic waves scattered by a stationary sphere have been conducted. Results show that the current solver has a good performance in modelling fluid–structure–acoustics interactions involving large deformations and complex geometries. It indicates that the immersed boundary method handled by delta function is accurate enough for predicting the dilatation and acoustics. We further adopt the present solver to model fluid–structure–acoustics interactions of flapping foils in forward flight and flapping foil energy harvester. The numerical examples presented here can also enrich the limited database of fluid–structure–acoustics interactions.

In the computation of flapping foils in forward flight, the present simulation captures the propagation and the decay of the sound pressure accurately. Two main mechanisms to generate positive and negative pressure are observed. The present simulation shows that the positive pressure formed on the loading face is much larger than the negative pressure generated by the vortex shedding. The directivity of the pressure wave also presents fluctuations around the stroke plane. The flexibility of the foil generates larger thrust with higher fluctuating pressure, but the geometrical shape does not have significant influences on the force and sound generation. Based on the FFT analysis of the fluctuating pressures, the fluctuating pressures are dominated by the flapping frequency $f$ of the foil for the rigid and NACA0015 flapping foil. However, both $f$ and $2f$ components are significant for the flexible foil.

In the computation of flapping foil energy harvester, it is found that the fluctuating pressure generated by NACA0015 foil, rigid plate and flexible plate are similar in terms of the frequency. The lift dominates the pressure fluctuation in the vertical direction and the drag dominates the pressure fluctuation in the horizontal direction. Some differences of the fluctuating pressure are also observed and analyzed. The results show that the geometrical shape does not have significant effects on the force and sound generation, while the flexibility of the plate tends to deteriorate the power extraction. The current flexible plate also induces larger sound at the frequency of $2f$ and weaker sound at the frequencies of $f$ and $3f$.

While many merits have been demonstrated by the cases presented, it should be pointed out for fluid–structure–acoustics interactions at large Reynolds numbers where turbulence flows exist, large eddy simulation and wall models should be considered. We will address this challenge in the near future.

Acknowledgements

Dr. F.-B. Tian is the recipient of an Australian Research Council Discovery Early Career Researcher Award (project number DE160101098). This work was conducted with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian Government.