An immersed boundary method for ﬂuid–structure–acoustics interaction atlow Reynolds numbers
Li Wang(1) , Fang-bao Tian(2) , Joseph C. S. Lai(3)(1) School of Engineering and Information Technology, University of New South Wales, Canberra ACT, 2600, Australia,l.wang@unsw.edu.au(2) School of Engineering and Information Technology, University of New South Wales, Canberra ACT, 2600, Australia,

f.tian@adfa.edu.au(3) School of Engineering and Information Technology, University of New South Wales, Canberra ACT, 2600, Australia,

j.lai@adfa.edu.au

Abstract
This paper presents an immersed boundary method for ﬂuid–structure–acoustics interactions involving complex
boundaries at low Reynolds numbers. In this method, the compressible Navier–Stokes equations are considered,
where a ﬁfth-order accuracy Weighted Essentially Non-Oscillation (WENO) scheme and a fourth-order central
difference scheme are adopted to discretize the convective and diffusion terms, respectively. The third-order
Runge-Kutta method is used for the temporal discretization. The non-linear ﬂexible structure immersed in theﬂuid is numerically solved by using a ﬁnite element method. By using a penalty immersed boundary method,the no-slip boundary between the ﬂuid and the structure is achieved. Acoustic waves scattering benchmark
problems and sound generation by a micro ﬂapping vehicle are presented to validate the present solver. Comparisons with the published data are presented to demonstrate the good performance of the present method inmodelling acoustics. The capability of using this method to study the sound generation by ﬂapping plates is also
demonstrated.
Keywords: Fluid–structure interaction, Immersed boundary method, Large deformation, Acoustics

1INTRODUCTION

Flow-induced sound extensively exists in engineering applications, such as aviation technologies [1], ventilation [2] and biomechanics [3]. Direct numerical simulation that solves the compressible Navier-Stokes equations
is an effective method for computational aeroacoustics (CAA). In the traditional computatiaonl ﬂuid dynamics(CFD), conformal mesh is generally used [4], which is time consuming especially when handling ﬂuid–structureinteraction (FSI) problems involving large deformations. The immersed boundary (IB) method developed byPeskin [5] has proven to be an efﬁcient method for this type of FSI simulations. By using this method, theconformal mesh is not needed and the regeneration of mesh is also avoided. Because of its simple boundarytreatment, IB method has gained popularity for a wide range of applications [6, 7]. Penalty IB (pIB) is a typical
IB method [8] where the IB is conceptually split into two Lagrangian components: one component is masslessand interacts with the ﬂuid exactly as the traditional IB method, and the other component carrying mass isconnected to the massless component by virtual stiff springs. However, most of the previous studies based onIB method focus on FSI problems without acoustics [9, 10, 11, 12]. Sun et al. presented an IB method which
considers the linear Euler equations for acoustic scattering modeling [13], but the viscosity of the ﬂuid is notconsidered.
Here, an immersed boundary method for ﬂuid–structure–acoustics interactions involving complex boundaries at
low Reynolds numbers are presented. The organization of the paper is as follows. The numerical approach isbrieﬂy introduced in Section 2. Two validations including acoustic waves scattered by a stationary cylinder andsound generation by a micro ﬂapping vehicle are presented in Section 3. The ﬁnal conclusions are given in2049

Section 4.

2

NUMERICAL METHOD

In this section, the in-house solver used will be introduced brieﬂy, and the details can be found in our previous
work [12, 14]. In this method, the compressible ﬂuid ﬂow, structural dynamics and ﬂuid–structure interaction
between them are considered. The structure is assumed to be elastic and its dynamics is governed by thefollowing non-linear equation
∂∂X∂ X4∂X∂X+| − 1)(1)ρs(KS |+ KB 4 = Ff ,∂t∂s∂s∂s∂swhere X is the Lagrangian coordinates of the ﬂexible beam, ρs is linear density, KS and KB are respectively
the stretching and bending rigidity, s is the arc coordinate, and F f is the external force acting on the beam.
Absolute nodal coordinates formulation proposed by Shabana [15] is adopted in the structure solver to solve
Eq. (1).
The ﬂow dynamics is governed by the two-dimensional compressible viscous Navier–Stokes equations∂Q ∂F ∂G1 ∂ Fu ∂ Gv+++ (+) = S,∂t∂x∂ y Re ∂ x∂y(2)where Q = [ρ, ρu, ρv, E]T , F = [ρu, ρu2 +P, ρuv, (E +P)u]T , G = [ρv, ρuv, ρv2 +P, (E +P)v]T , Fu = [0, τxx , τxy , bx ]T ,Gv = [0, τxy , τyy , by ]T , bx = uτxx + vτxy ,by = uτxy + vτyy . ρ is the ﬂuid density, u, v and w are respectively thethree velocity components, P is the pressure, E is the total energy, S is a general source term including theIB-imposed Eulerian force and other body forces, Re is the Reynolds number, and τi j is the shear stress.
In the ﬂuid solver, the ﬁfth-order accuracy WENO scheme proposed by Liu et al. [16] is employed for thespatial discretization of the convective term. For the viscous term, a fourth-order central difference scheme isused to discretize the spatial derivatives. For all unsteady equations in the ﬂow solver, the third-order TVD
Runge-Kutta method is used for time discretization. The dynamics of the ﬂuid and ﬂexible structures are solved
independently and the interaction force is calculated explicitly using a feedback law based on the pIB method.

3VALIDATIONS AND APPLICATIONS

3.1 Acoustics wave scattered by a stationary cylinder
In order to verify the accuracy of the method in simulating the acoustic problem, acoustic wave scattered from astationary cylinder is validated against the published data [17]. In this problem, a rigid cylinder with a diameterof D = 1.0 is ﬁxed with its center at (0, 0). The initial conditions are localized perturbations of the Gaussian
distribution of pressure P,(x − 4)2 + y2](3)
P = εexp[−ln(2)0.04
where ε = 10−3 . The initial density and pressure are respectively 1.0 and 1/γ (the adiabatic coefﬁcient γ = 1.4),and the initial velocity of the ﬂuid is zero.
We used a uniform Cartesian mesh over a computational domain of (16D × 12D) with three different meshspacings: D/100, D/50 and D/25. The time histories of the ﬂuctuation pressure, Δp = P − p0 where P and p0
are respectively the current and initial pressure, from a numerical probe located at (2, 0) are plotted in Fig. 1
along with the exact results reported in Ref. [17]. According to the comparisons in Fig. 1, all the three meshspacings can well capture the ﬂuctuating pressure, and the small discrepancies depicted in the scaled up ﬁgures
can be diminished when the mesh size is D/100.
Comparison of Δp contours obtained by Bailoor et al. [18] and this study at the time of 6.0 in Fig. 2 showsthat the collision, merge and reﬂection of the acoustic waves are well captured by our numerical method.

2050

12x 10−5Δ x= Δ y = D/100
Δ x= Δ y = D/50Δ x= Δ y = D/25
Liu and Vasilyev (2007)Δp60−6024t6810Figure 1. Time histories of the ﬂuctuating pressure measured at (2,0).

Figure 2. Fluctuating pressure at t = 6 obtained by Bailoor et al. [18] (left) and the present simulation (right).
The contour level ranges from −2.0 × 10−5 to 2.0 × 10−5 .
3.2 Sound generation by a micro ﬂapping vehicle
Here, a micro ﬂapping vehicle in hovering ﬂight is considered to validate the current solver in handling moving
body with relatively complex geometry conﬁguration. The schematic of the problem is shown in Fig. 3. A
circular cylinder (body) with a diameter of 0.5L, and two elliptic cylinders (wings) with the dimensions of Land 0.4L are used to model the micro ﬂapping vehicle. The two wings ﬂap symmetrically, and the ﬂapping
motion of the wings is prescribed byα(t) = α0 [1 + sin(2π f t)],(4)
where the amplitude of the ﬂapping angle α0 is 25.3o and f is the ﬂapping frequency. The Reynolds number
deﬁned by ρ f Umax L/μ is 200, where Umax is the maximum wing tip velocity. The Strouhal number deﬁned byf L/Umax is 0.25, and the Mach number based on Umax 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. 4. It is found that these results

2051

Figure 3. Schematic of the prescribed ﬂapping motion of a micro ﬂapping vehicle.

630-3468Figure 4. Sound generation by a micro ﬂapping vehicle in hovering ﬂight: comparison of the lift generated by
the wing and body.

90o

0.00040.00030.00020.0001180o0oPresent270o
Seo and Mittal (2011)

Figure 5. Sound generation by a micro ﬂapping vehicle in hovering ﬂight: comparison of the root-mean-squareof ﬂuctuating pressure measured at a distance of 50L.

2052

Figure 6. Schematic of a rigid (left) and ﬂexible (right) foil ﬂapping in a uniform ﬂow.
agree well with the data from Ref. [19]. The ﬂuctuating pressure Δ p̃ is deﬁned by Δ p̃(x, y,t) = Δp(x, y,t) −Δ p̄(x, y), where Δp denotes the total ﬂuctuating pressure and p̄ is the time average pressure. Δp = p − p∞ withp∞ being the ambient pressure. Good agreements of the root-mean-square of ﬂuctuating pressure measured ata distance of 50L shown in Fig. 5 conﬁrm the good performance of the current method in handling acoustics
involving complex moving geometries.
3.3 Sound generation by a ﬂexible ﬂapping plate
Here, the acoustic perturbations induced by ﬂapping foils in forward ﬂight are considered, as shown in Fig. 6 [20].
The foil is clamped at the leading edge, and the clamping device undergoes a combined translational and rotational motions, given by [20]A0X 0 (t) = cos(2π f t)[cosβ , sinβ ],2(5)αmα(t) =sin(2π f t + φ ).
2where A0 is the translational amplitude, f is the ﬂapping frequency, αm is the rotation amplitude, and β is theangle between the stroke plane and the horizontal plane. The phase difference (φ ) is 0.
The non-dimensional parameters including the ﬂapping amplitude, Reynolds number, inlet velocity, Mach number, mass ratio and frequency ratio that control this problem are given respectively byA0,
LRe =ρ f UL,μUr =U0,U

M=πU,2cm∗ =ρs,
ρf L

ω∗ =

2π f,
ωn(6)

where L is the chord length, U = 2 f A0 is the average translational velocity of the leading edge, f is the ﬂapping
frequency of the foil,c is the sound speed of the ﬂuid, ρ f is the ﬂuid density, ρs is the linear density of thefoil, and ωn = kn2 /L2 EB /ρs with kn = 1.8751 (the frequency of the ﬁrst natural vibration mode of the wingwith ﬁxed leading edge) and EB being the bending rigidity. Here, Re=100, m∗ = 5.0 and Ur = 0.4. The thrust,
lift and power coefﬁcients are deﬁned as
CT =2FT,ρ f U 2L

CL =2FL,ρ f U 2L

Cp =−2 0L F f · vdl,
ρ f U 3L

(7)

where FT and FL are respectively the thrust and lift force acting on the foil by the ambient ﬂuid, F f is thehydrodynamic traction on the foil, and v is the velocity of the foil. The ﬂuctuating pressure Δ p̃ is deﬁned in thesame way in Section 3.2. All the pressures presented afterwards are scaled by the ﬂuid density ρ f and soundspeed c.
In order to compare the sound generation of rigid and ﬂexible foils, the polar diagrams of the root-mean-squareof Δ p̃ measured at r = 30L are presented in Fig. 7. As shown in Fig. 7 (a) and (c), Δ p̃ generated by therigid plate and NACA0015 foil distributes symmetrically about the stroke plane (indicated by the dashed line2053

(a)

(b)
(c)

R90o0.0008

R




0.0006





0.0004



R

RR
 |Δp̃| = 0.0008

R0.0002

R

R |Δp̃| = 0.0008

180o0oo270 |Δp̃| = 0.0008

Figure 7. Flapping foils in forward ﬂight: root-mean-square of Δ p̃ measured at r = 30L: (a) rigid plate, (b)ﬂexible plate and (c) NACA0015 foil. Re=100, M = 0.1, Ur = 0.4, A0 /L = 1.25 and αm = π/4. The dashedlines indicate the directions of the stroke planes.
in Fig. 7 (a) and (c)). However, Fig. 7 (b) shows that the distribution of the ﬂuctuating pressure generated by
the ﬂexible foil is asymmetrical. An evident shift (about 15o anticlockwise) is observed, as shown in Fig. 7(b). As analyzed in Ref. [20], the deformation of the ﬂexible plate during upstroke is higher than that during
downstroke, due to the presence of free stream. Obviously, this sound shift is introduced by the asymmetricaldeformation of the foil. Fig. 7 also shows that the ﬂexible plate generates the largest ﬂuctuating pressure,which is signiﬁcantly larger (up to 50%) than that generated by the rigid plate and NACA0015 foil. The results
indicate that the ﬂexibility inﬂuences the sound generation signiﬁcantly. The effects of geometrical shape on thesound amplitudes are much smaller. In addition, it is found that the ﬂuctuating pressure on the lower section
is larger than that on the upper section. A reasonable explanation is that the lower section is located on thewindward side and the upper section is on the leeward side.

4CONCLUSIONSAn immersed boundary method for ﬂuid–structure–acoustics interactions involving complex boundaries at lowReynolds numbers is presented. In this method, the compressible Navier–Stokes equations and the non-linear
ﬂexible structure immersed in the ﬂuid are numerically solved by using a ﬁnite difference method and a ﬁnite
element method, respectively. By using a penalty immersed boundary method, the no-slip boundary between theﬂuid and the structure is achieved. Comparisons with the published data indicate the good performance of thismethod in modelling acoustics for acoustic waves scattering by a rigid cylinder and the sound generation by amicro ﬂapping vehicle. The capability of this method to study the sound generation by ﬂapping plates is also
demonstrated.

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.

REFERENCES[1] Colonius, Tim and Lele, Sanjiva K, 2004. “Computational aeroacoustics: progress on nonlinear problems ofsound generation”. Progress in Aerospace Sciences, 40(6), pp. 345–416.

2054

[2] Waye, K.P., Rylander, R., Benton, S. and Leventhall, H.G., 1997. “Effects on performance and work quality
due to low frequency ventilation noise”. Journal of Sound and Vibration, 205(4), pp. 467–474.
[3] Seo, Jung Hee and Mittal, Rajat, 2011. “A high-order immersed boundary method for acoustic wave scattering and low-Mach number ﬂow-induced sound in complex geometries”. Journal of Computational Physics,230(4), pp. 1000–1019.
[4] Tian, F.B., Bharti, R.P. and Xu, Y.Q., 2014. “Deforming-Spatial-Domain/Stabilized Space–Time (DSD/SST)
method in computation of non-Newtonian ﬂuid ﬂow and heat transfer with moving boundaries”. Computational Mechanics, 53(2), pp. 257–271.
[5] Peskin, C.S., 1977. “Numerical analysis of blood ﬂow in the heart”. Journal of Computational Physics,
25(3), pp. 220–252.
[6] Tian, F.B., Dai, H., Luo, H., Doyle, J.F. and Rousseau, B., 2014. “Fluid–structure interaction involving large
deformations: 3D simulations and applications to biological systems”. Journal of Computational Physics,
258, pp. 451–469.
[7] Huang, W.X., Shin, S.J. and Sung, H.J., 2007. “Simulation of ﬂexible ﬁlaments in a uniform ﬂow by theimmersed boundary method”. Journal of Computational Physics, 226(2), pp. 2206–2228.
[8] Kim, Y. and Peskin, C.S., 2007. “Penalty immersed boundary method for an elastic boundary with mass”.
Physics of Fluids, 19(5), p. 053103.
[9] Mittal, R. and Iaccarino, G., 2005. “Immersed boundary methods”. Annu. Rev. Fluid Mech., 37, pp. 239–261.
[10] Tian, Fang-Bao and Luo, Haoxiang and Zhu, Luoding and Liao, James C and Lu, Xi-Yun, 2011. “An efﬁcient immersed boundary-lattice Boltzmann method for the hydrodynamic interaction of elastic ﬁlaments”.
Journal of Computational Physics, 230(19), pp. 7266–7283.
[11] Ghias, Reza and Mittal, Rajat and Dong, Haibo, 2007. “A sharp interface immersed boundary method forcompressible viscous ﬂows”. Journal of Computational Physics, 225(1), pp. 528–553.
[12] Wang, L., Currao, G.M., Han, F., Neely, A.J., Young, J. and Tian, F.B., 2017. “An immersed boundary method for ﬂuid–structure interaction with compressible multiphase ﬂows”. Journal of ComputationalPhysics, 346, pp. 131–151.
[13] Sun, Xiaofeng and Jiang, Yongsong and Liang, An and Jing, Xiaodong, 2012. “An immersed boundary
computational model for acoustic scattering problems with complex geometries”. The Journal of the Acoustical Society of America, 132(5), pp. 3190–3199.
[14] Wang, Li and Tian, Fang-Bao and Lai, Joseph CS, 2019. “An immersed boundary method for ﬂuid–structure–acoustics interactions involving large deformations and complex geometries”. Under review.
[15] Shabana, A.A., 1997. “Flexible multibody dynamics: review of past and recent developments”. MultibodySystem Dynamics, 1(2), pp. 189–222.
[16] Liu, Xu-Dong and Osher, Stanley and Chan, Tony, 1994. “Weighted essentially non-oscillatory schemes”.
Journal of Computational Physics, 115(1), pp. 200–212.
[17] Liu, Q. and Vasilyev, O.V., 2007. “A Brinkman penalization method for compressible ﬂows in complexgeometries”. Journal of Computational Physics, 227(2), pp. 946–966.
[18] Bailoor, S., Annangi, A., Seo, J.H. and Bhardwaj, R., 2017. “Fluid–structure interaction solver for compressible ﬂows with applications to blast loading on thin elastic structures”. Applied Mathematical Modelling,
52, pp. 470–492.

2055

[19] Seo, Jung-Hee and Mittal, Rajat, 2011. “Computation of aerodynamic sound around complex stationaryand moving bodies”. 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace
Exposition, pp. 1087.
[20] Tian, Fang-Bao and Luo, Haoxiang and Song, Jialei and Lu, Xi-Yun, 2013. “Force production and asymmetric deformation of a ﬂexible ﬂapping wing in forward ﬂight”. Journal of Fluids and Structures, 36,
pp. 149–161.

2056