The suppressed sound of inhomogeneous ﬂapping airfoils
Michael WEIDENFELD(1) , Eran ARAD(2)(1) RAFAEL, Israel, wmichael@rafael.co.il(2) RAFAEL, Israel

Abstract
We study the aeroacoustic problem of a ﬂapping airfoil in a high-Reynolds and low-Mach regime, and propose
novel means to reduce the far-ﬁeld sound originated by the ﬂuid-structure interaction. The present model consists of a two-dimensional, elastic airfoil, periodically actuated at the leading edge. The near-ﬁeld is analyzed
using thin airfoil theory, and a discrete-vortex wake-model. The derived airfoil motion and wake dynamics arethen introduced as a source term to Powell-Howe’s acoustic analogy, to project the sound towards a distant observer. Considering the elastic ﬂapping conﬁguration, we seek optimal material properties and a linear thickness
distribution to reduce the sound of an otherwise rigid airfoil while retaining the lift amplitude of the rigid conﬁguration. To this end, the aeroacoustic model is introduced into an optimization scheme, yielding an effective
Pareto front, reﬂecting the trade-off between the competing objectives of lower sound levels and aerodynamic
efﬁciency. Compared with the rigid heaving airfoil, over 10 [dB] sound reduction was obtained for the optimal ﬂexible conﬁguration with linear thickness distributions, producing the same lift amplitude. The optimal
conﬁgurations are obtained by time shifting the motion- and wake-dipoles to antiphase mode, thus mitigating
the total radiated sound, while keeping the lift-amplitude value. The opposing sound dipoles follow from thesynchronized disposition of the airfoils’ motion and circulation signals.
Keywords: Vortex sound, Optimization, Airfoil noise, Flapping noise, Vibro-acoustics

Introduction
The canonical ﬂuid-structure interaction problem involving a ﬂapping ﬂexible sheet (ﬂag, airfoil, sail, etc.) immersed in a uniform ﬂow, is analyzed in the present work in an attempt to study novel means to mitigate theassociated far-ﬁeld noise. In general, the unsteady ﬁlament motion, and the interaction between the surface andwake vorticity, give rise to acoustic pressure-waves which then radiate as sound (1). In line with the ongoinginterest in controlling noise generated aerodynamically by airborne vehicles (2), for both civil and military applications, the study of far-ﬁeld sound originating from unsteady surface motion, plays a key role.
The associated noise problem of ﬂapping elastic sheets was recently studied in theoretical works which examined the sound radiation of rigid and ﬂexible airfoils in unsteady motion and ﬂow conditions (3, 4). Inspired byowls “silent ﬂight” technology (5), incorporating non-uniform geometry and elastic properties as a methodology
for noise reduction, Weidenfeld and Manela (6) studied the role of airoil’s permeability in reducing ﬂapping
noise. Their results reveal the negating effect of the opposing seepage- and motion-dipoles, on the total soundradiation.
In the context of ﬂapping-wings noise, variable geometry and mass distribution is often considered, and introduced into various optimization schemes (7, 8), where a Pareto front of best possible designs is sought. The
objective of the present work is to investigate means to control and reduce the sound of a ﬂapping elasticﬁlament, with regards to the aerodynamic performance, by varying geometry and structural features along thechord. To this end, we consider a two-dimensional setup of a thin elastic airfoil, periodically actuated at the
upstream leading-edge, and free at the downstream trailing-edge. The near-ﬁeld description includes a discretewake model, and the far-ﬁeld acoustics is obtained via Powell-Howe acoustic analogy (1), valid for low-Machand high-Reynolds number ﬂow regime as is considered in the present work. Considering this setup, an optimization problem is formulated and a Pareto-optimal solution is obtained.

596

Problem formulation

Figure 1. Problem schematics: An elastic airfoil immersed in a uniform ﬂow and actuated at the leading edgeby small amplitude heaving motion. The model includes trailing edge wake shedding in the form of discrete
point vortices with negative (red) and positive (blue) circulations.
The problem schematics is shown in ﬁgure 1. An inextensible, thin ﬂexible airfoil extending over the interval−a ≤ x1 ≤ a is immersed in a uniform ﬂow of mean density ρ0 with velocity U in the x1 -direction. The
structure has a span L in the x3 -direction and thickness distribution σ (x1 ) with max(σ )  a  L such that the
airfoil motion may be regarded two-dimensional in the x1 x2 -plane. At t ≥ 0 leading-edge sinusoidal heaving
actuation is applied,(1)ξ (x1 = −a) = ε h a sin(ωht).
Here, ξ (x1 ,t) marks the ﬁlament displacement in the x2 direction, ε h  1 the scaled heaving amplitude (hereafter, overbars marks non-dimensional quantities), and ωh denotes the heaving frequency. The displacement
ξ (x1 ,t) is governed by the linear Euler-Bernoulli’s beam equation


∂2
∂ 2ξ
∂ 2ξ
(2)
ρs (x1 ) 2 + 2 B(x1 ) 2 = −Δp(x1 ,t).
∂t
∂ x1
∂ x1
In (2), ρs (x1 ) and B(x1 ) mark the mass per unit area and bending rigidity chord-wise distributions, in units[kg/m2 ] and [N·m], respectively. For a uniform structure in the x3 direction, these properties are given byρs (x1 ) = ρm σ (x1 ),(3)

B(x1 ) = Eσ (x1 )/12.

(4)

3

In (3)-(4), ρm is the material density, and E is Young’s modulus. Assuming high Reynolds and low Machconditions, and a small amplitude motion induced by the leading-edge heaving, the ﬂow is regarded potential,
and thin airfoil theory approximation is applied. At the small amplitude motion and induced angles of attack,the ﬂow is considered attached throughout the airfoil and so vorticity is released continuously from the trailingedge. The present wake model consists of discrete line vortices, where at each time step t = nΔt, a new vortex
is formed and placed adjacent to the instantaneous trailing edge position (x1 , x2 ) = (a +UΔt, ξ (a,t)). The vortex
strength Γk is ﬁxed by Kelvin’s theorem and the instantaneous airfoil circulation, and its trajectory follows fromthe potential ﬂow computation.

Aeroacoustic scheme
The entire aeroacoustic problem is described in details by Weidenfeld et al. (9) and a short summary is given
here, for brevity. As follows from thin airfoil theory, the airfoil is represented by its circulation distribution

597

along the chord γ(x1 ,t). Turning to complex notation, z = x1 + ix2 represents location on the complex z-plane
and W (z) marks the conjugate velocity ﬂow ﬁeld derived from the potential theory. The airfoil impermeability
condition is then given by
∂ξ
∂ξ
+U
(5)
= − Im {W (z)|−a≤x1 ≤a } .
∂t∂ x1
The pressure jump Δp = p+ − p− across the airfoil, appearing on the right-hand side of (2), marks the ﬂuid
loading owing to the lower (p− ) and upper (p+ ) surface pressure difference, and is related to the airfoil circulation γ(x1 ,t) via the unsteady Bernoulli’s equation∂ x1
γ(s,t) ds.
(6)
∂t −a
Wake vortices dynamics is coupled to the system through the right-hand side of the impermeability condition(5). In addition, Kelvin’s theorem is imposedΔp(x1 ,t) = ρ0Uγ(x1 ,t) + ρ0

n−1

 a

k=1

−a

ΓN = − ∑ Γk −γ(x1 ,t) dx1 ,(7)ensuring the system total circulation vanishes at all times, and ﬁxing the strength of the recently shed vortex
ΓN . To complete the formulation of the problem and set regular ﬂow along the airfoil surface and ﬁnite velocity
at the trailing edge, we apply the unsteady Kutta condition, ﬁxing the vanishing trailing edge circulation. Inline with the equation of motion (2), and free-end condition at x1 = a, the well imposed problem consists ofthe four boundary condition

 2 
 3 

∂ξ
∂ ξ
∂ ξ
ξ (−a,t) = ε h a sin(ωht),
= 0,
=
0,
= 0,
(8)
∂ x1 (−a,t)∂ x12 (a,t)∂ x13 (a,t)
and accompanied by the homogeneous initial state of the airfoil at t = 0.
Having formulated the near-ﬁeld problem, we now turn to the acoustic radiation scheme in which the former
serves as an effective “source term” to be radiated to the far-ﬁeld. In the present low-Mach and high-Reynolds
ﬂow regime, the acoustic pressure is governed by the linear inhomogeneous wave equation (1)


1 ∂2
∂ 2ξ
2
−
∇
δ ( fa ) + ρ0 ∇ · (Ω × v) ,
(9)
p
=
ρ
0
∂t 2
c20 ∂t 2
where c0 marks the mean speed of sound, v is the ﬂuid velocity, δ is the Dirac delta function and the function
fa = fa (x,t) vanishes on the airfoil surface. The vorticity ﬁeld, denoted by Ω is given by the sum of all wakevortices circulation
n


(10)
Ω = x̂3 ∑ Γk δ x − xΓk (t) .
k=1

In (10), xΓk (t) marks the instantaneous k-th vortex location, and the hat denotes a unit vector. The total acousticpressure is decomposed into the sum of motion (pξ ) and wake (pw ) dipole contributions,ptot (x,t) = pξ (x,t) + pw (x,t),
where

 a
dτ
∂ 2ξ
√
a2 − y21 2 (y1 , τ) dy1 ,∂τ−∞ tr − τ −a⎛
 tr V (2) (τ) dτ
n
∂ ⎝
ρ0 Γk
Γk
√
sin θ
pw (x,t) = ∑ ∂ttr − τ
−∞
|x|
2π
2c
0
k=1

pξ (x,t) =ρ0 cos θ ∂
2π 2c0 |x| ∂t− cos θ

598

(11)

 tr

 tr

dτ
√
−∞ tr − τ



∂Y2
∂Y2
(1)
(1)
VΓk (τ).
−VΓk (τ)∂ y2∂ y1 xΓ (τ)
k

(12)

(13)

( j)

In (12)-(13), VΓk marks the k-th vortex velocity component in the x j -direction. The angle θ speciﬁes the farﬁeld observer’s inclination, measured clockwise from the x2 -axis. Accordingly, the lift and suction dipoles are
identiﬁed as the terms involving the cos θ = x2 / |x|, and sin θ = x1 / |x|, respectively.
We focus on acoustic wavelengths λ0 far larger than the airfoil chord λ0  a, which equivalently speciﬁes
ω 0 M  1, where ω 0 = ω0 a/U is the reduced acoustic frequency. In line with the low Mach number working
frame, this conﬁnes the analysis to low-moderate acoustic frequencies satisfying ω 0  1/M, and allowing theapplication of compact scheme using the two-dimensional compact Green’s function (1), introduced into equations (12)-(13).
Before proceeding and for future reference we introduce the total acoustic lift dipole (θ = 0)
pLift
tot (x2 ,t) = ptot (x1 = 0, x2 ,t) = pξ (x1 = 0, x2 ,t) + pw (x1 = 0, x2 ,t)(14)

Optimization scheme
Consider a rigid heaving airfoil, producing unsteady lift with amplitude Lrigid , and radiating sound to a distantobserver in the x2 -direction with amplitude ptot rigid , where the norm sign x denotes the maximum normof the variable x, throughout the paper. We seek an alternate elastic, non-homogeneous airfoil, with optimal
thickness distribution along the chord, to minimize the amplitude of the far-ﬁeld acoustic lift dipole withoutcompromising on the rigid airfoil’s lift amplitude. Accordingly, The optimization problem is given byminimizewith respect tosubject topLifttot ,(15a)

σ (x1 ), ρm , E,

(15b)

L ≥ Lrigid (±33%),(15c)

where the cost function is the lift dipole amplitude, the input variables are the airfoil thickness distribution, material density, and the Young modulus, and the minimization is subjected to a minimum lift amplitude constraint
with a 33% tolerance for faster convergence and to obtain a sufﬁciently wide Pareto Front. The optimization
R software (10), and involves a mixed-strategies practice, inprocess is executed using the modeFRONTIER
which three methods are combined consecutively: the evolutionary MOGA-II (11), the heuristic Nelder & Mead
Simplex (12), and the gradient based NLPQLP method (13). The genetic search is the process through which
the a Pareto front materializes. The Simplex and gradient based search are then used to ﬁne-tune a design
around a speciﬁc constraint value.

Scaling and numerical procedure
The problem outlined hitherto, consists of the near- and far-ﬁeld descriptions and the optimization framework.
We now proceed and devise the non-dimensional form of the problem using a, U, a/U, ρ0U 2 , and 2πaU
as reference quantities for the length, velocity, time, pressure and vortices circulation, respectively. The nondimensional quantities, marked by overbars, governing the aeroacoustic problem are
ρ=

ρm
,
ρ0

E=

E
,
12ρ0U 2σ (x1 ) =σ (x1 ),a

L(t) =

L(t),ρ0U 2 a(16)

denoting the material to air density ratio, normalized elasticity, thickness ratio, and the normalized lift force,
respectively. Applying the same scaling to the acoustic problem (equations (12), (13), and (11)), the nondimensional acoustic pressure takes the form
ptot (x,t) =M
πtot =8 |x|M π (t r ) + πw (t r ) ,8 |x| ξ(17)599

where πξ , and πw represent the acoustic motion- and wake- “kernels”, respectively. In addition, the coupled
dynamical-acoustic problem is further governed by x2aε h , ω h = ωh , θ = arccos
(18)
|x|
U
marking the scaled heaving amplitude and angular frequency, and the observer far-ﬁeld directivity, respectively.
The scaled optimization problem is then speciﬁed byminimizewith respect tosubject to

Lift
,
πtot

(19a)

σ (x1 ), ρ, E,(19b)

L ≥ Lrigid (±33%),(19c)

In all cases considered, the heaving amplitude and angular frequency were ﬁxed by ε h = 0.01, and ω h = 1 inaccordance with the small amplitude motion and attached ﬂow throughout the airfoil assumptions, set in Sec.
3. In the present setup, only combinations leading to the system’s periodic state are considered. The numerical
procedure requires discretization in both space and time. Thus the airfoil circulation distribution γ(x1 ,t) is
decomposed into Nγ = 100 discrete vortices, and the time step is set to π/65, which proved sufﬁcient for
convergence. The system of equations was integrated using Euler’s explicit method for 8 heaving cycles, unlessdivergence was detected ﬁrst, originating from unstable parameters combination. We focus on linear thickness
distribution of following form
σ (x1 ) = β1 x1 + β0 ,(20)
and begin by seeking ﬁrst an optimal elastic structure whose thickness is uniform throughout the chord (β1 = 0).
We then move on to a structure whose thickness varies linearly along the chord.

Results
6.1 Uniform thickness distribution
In the present section we seek a heaving elastic ﬁlament with uniform thickness which produces the same lift
amplitude as it’s rigid counterpart and radiates lower sound levels to a far-ﬁeld observer. All far-ﬁeld results
are given in the lift dipole direction and so the π Lift superscript is omitted from all dipoles for brevity. Figure(2) presents a scatter map relating the sound reduction with respect to the rigid case (vertical axis), and thelift amplitude ratio (horizontal axis). The blue markers represent all feasible conﬁgurations (satisfying the lift
amplitude constraints), and the red markers represent the discrete Pareto front. The horizontal dashed line marks
zero sound reduction, and the vertical dashed line marks the rigid setup lift amplitude. The uniform thickness
distribution results, identiﬁed by the square markers, indicate that no stable elastic conﬁguration of uniform
thickness produces lower sound levels without an aerodynamic compromise. A minor sound reduction of about∼ 2 [dB] is obtained for the substantial cost of a 33% in lift amplitude. We therefore proceed to examine
variable thickness conﬁgurations, starting with a linear distribution.
6.2 Linear thickness distribution
By linearly varying the airfoil thickness along the chord, numerous conﬁgurations are found to both reducesound levels and satisfy the lift amplitude constraint, as illustrated by the dot markers in ﬁgure (2). In ﬁgure(2), the red dot markers represent the actual discrete Pareto front obtained from the optimization process, andthe solid curve is a suggested theoretical continuous front in the form of linear relation between the soundreduction and lift amplitude ratio. In the uniform thickness case, there was no indication for such linear relation. Additionally, the results illustrated in ﬁgure (2), visualize the differences between the two Pareto fronts,
highlighting the radical improvement in sound reduction, simply by allowing the the airfoil thickness to varylinearly along the chord.

600

10

uniform distribution

0
-10
-20
linear distribution
-30

0.8

1

1.2

1.4

1.6

1.8

Figure 2. Scatter map of noise reduction (in [dB]) vs. lift ratio depicting real conﬁgurations (blue markers)
and the discrete Pareto front (red markers). Square markers represent uniform thickness distribution results, dots
represent linear thickness distribution results. The solid black line marks a theoretical Pareto front. The horizontal dashed line marks zero sound reduction, and the vertical dashed line marks the rigid setup lift amplitude.
The dashed lines intersections mark the rigid conﬁguration setup.
The optimal results are rationalized through the mechanism suggested in ﬁgures (3a)-(3c). All three ﬁgures
present the motion and wake lift-dipole contributions, along with the total lift kernel, obtained for the rigid
setup (3a), the optimal conﬁguration satisfying L = Lrigid (3b), and the optimal conﬁguration satisfying
L = 0.67Lrigid (3c). The optimal conﬁgurations clearly present similar trends, in ﬁgures (3b)-(3c), whichsigniﬁcantly differ from that of the rigid conﬁguration in (3a). The results indicate that optimal conﬁgurations
are obtained when two conditions are met:
1. The motion dipole is similar in magnitude to that of the rigid case,
2. The motion and wake dipoles are in antiphase.
The motion sound reﬂects the unsteady forces exerted on the ﬂuid by the airfoil motion, consequently, keeping
it’s magnitude similar to that of the rigid case, maintains the lift constraint. In addition, setting the motion andwake dipoles in antiphase acts to mitigate the overall sound of the ﬂuid-structure interaction.
Turning to the near-ﬁeld dynamics, under the non-dimensional actuation frequency of ω h = 1, all Pareto designs
exhibit an airfoil motion resembling the ﬁrst mode of a cantilever beam, as suggested by the deﬂection obtained
for the L = Lrigid optimal results presented in ﬁgure (4a). Each point along the chord features a phase lagwith respect to the leading edge motion. To represent the entire beam, we may deﬁne an average beam motion
 1 1ξ =ξ (x1 ,t) dx1 .
2 −1

(21)

Recalling Kelvin’s theorem from eq. (7), the airfoil and wake circulations amount to zero, so that Γairfoil (t) =−Γwake (t). Figure (4b) compares between the average beam motion and both circulation time-variations over a
period of motion for the L = Lrigid optimal design. The trends provide a complementary dynamical viewof the optimal results, suggesting that synchronizing the airfoil motion and circulation is required for ultimately
reducing far-ﬁeld sound. It follows that the motion and wake-circulation are set to be in antiphase, whichthen promotes the counterpart antiphased sound-dipoles leading to the desired “silent” outcome. Figure (4c)

601

(a)

30

(b)

30

15

15

15

0

0

0

-15

-15

-15

T/4

T/2

3T/4

T

T/4

T/2

3T/4

(c)

30

T

T/4

T/2

3T/4

T

Figure 3. Sound components for three problem setups. (a) Rigid setup, (b) optimal conﬁguration satisfying
L = Lrigid , (c) optimal conﬁguration satisfying L = 0.67Lrigid . Three different curves mark the motion), wake dipole (), and total sound dipole ().
dipole (presents the absolute phase-difference between average beam motion and wake-circulation (blue), and airfoilcirculation (red), for all Pareto optimal results. The results reinforce the previous observation and indicatethat all Pareto results exhibit locked-in motion and airfoil-circulation signals, and antiphased motion and wakecirculation signals.
Examining the geometry features, the results indicate that all Pareto-optimal designs share a negative slope(a)
Downwards

5

0

0

-0.01

-5

0

0.5

1

3 /4
/2
/4

-10

Upwards
-0.5

(c)

5 /4

10

0.01

-0.02

(b)

15

phase [rad]

0.02

T/4

T/2

3T/4

T

60

80

100

120

Figure 4. Near ﬁeld dynamics. (a) Airfoil instantaneous deﬂection at ξ (−1,t) = 0. The blue curve corresponds
to a downward accelerating airfoil, and the red curve corresponds to upward acceleration. (b) The solid red
curve represents the average airfoil motion ξ , the solid red curve represents the wake circulation, and thedotted red curve represents the airfoil circulation. Results are plotted over a motion period for the L = Lrigid
optimal design. (c) Blue curve marks the absolute phase-difference between average beam motion and wakecirculation, and the red curve marks the absolute phase-difference between average beam motion and airfoilcirculation. Both curves represent the Pareto-optimal results for linear thickness distribution and plotted forvarious lift-constraint values.
resulting with a thicker leading edge and a thinner trailing edge, where the higher-lift designs feature smaller
slope values and the lower-lift designs are of slightly larger slopes.

Conclusion
In the present study we considered the effect of non-uniform thickness distribution on the sound radiation of anelastic, leading-edge actuated, two-dimensional thin airfoil, in a low-Mach potential uniform ﬂow. Using thin602

airfoil theory along with a discrete wake model, the near-ﬁeld dynamics was obtained and introduced as aneffective dipole-type source term to the Powell-Howe acoustic analogy, to obtain the far-ﬁeld sound.
The aeroacoustic setup, was then introduced into an optimization scheme where optimal material properties andthickness distribution was sought. The objective was to minimize ﬂapping sound with respect to a rigid heaving
conﬁguration while retaining the rigid setup lift-amplitude value. It was found that uniform ﬂags inefﬁciently
reduce sound in about ∼2 [dB], compared with the rigid airfoil, at the cost of 33% of the lift-amplitude.
In contrast, the optimal linear thickness conﬁguration was found to reduce ﬂapping sound in about ∼10 [dB]compared with the rigid setup, while retaining the same lift amplitude. The Pareto front was found to exhibit
an approximated linear relation between sound reduction and lift-amplitude ratio, exhibiting a further substantial
sound reduction of up to ∼30 [dB] for smaller lift amplitude values.
The mechanism through which sound reduction is obtained for a speciﬁc lift-amplitude value is two-fold: (1)The motion and wake dipoles are shifted to counteract each other and reduce the total sound emitted, (2) the
motion dipole, reﬂecting the unsteady forces excreted on the ﬂuid by the airfoil motion, is ﬁxed in magnitude
to the rigid setup value, thus maintaining the desired lift amplitude value. The antiphase merit between the twosound dipoles was found to follow from phase-locking the airfoil motion and circulation.

REFERENCES
1 Howe MS. Theory of Vortex Sound. Cambridge Texts in Applied Mathematics. Cambridge University Press;
2003.
2 Zaporozhets O, Tokarev V, Attenborough K. Aircraft Noise. New York: Taylor and Francis; 2011.
3 Manela A. Vibration and sound of an elastic wing actuated at its leading edge. Journal of Sound andVibration. 2011;331:638–650. doi:10.1016/j.jsv.2011.09.020.
4 Manela A. On the acoustic radiation of a pitching airfoil.
doi:10.1063/1.4816295.

Physics of Fluids. 2013;25(7):071906.

5 Kroeger RA, Grushka HD, Helvey TC. Low speed aerodynamics for ultra-quiet ﬂight. Tech Rep AFFDLTR-71-75. 1972;doi:10.1016/0022-460X(71)90105-2.
6 Weidenfeld M, Manela A. On the attenuating effect of permeability on the low frequency sound of an
airfoil. Journal of Sound and Vibration. 2016;375:275–288. doi:10.1016/j.jsv.2016.04.002.
7 Stanford BK, Beran PS. Analytical Sensitivity Analysis of an Unsteady Vortex-Lattice Method for FlappingWing Optimization. Journal of Aircraft. 2010;47(2):647–662. doi:10.2514/1.46259.
8 Stewart EC, Patil MJ, Canﬁeld RA, Snyder RD. Aeroelastic Shape Optimization of a Flapping Wing. Journalof Aircraft. 2014;53(3):13–17. doi:10.2514/1.C033278.
9 Weidenfeld M, Arad E. Mitigating the sound of a ﬂapping airfoil using optimal structural properties distributions. Journal of Sound and Vibration. 2018 oct;432:235–248. doi:10.1016/j.jsv.2018.06.016.
10 ESTECO. modeFRONTIER, [Software Package]; 2014.
11 Poles S. MOGA-II an improved multi-objective genetic algorithm. ESTECO - Technical Report. 2003;.
12 Poles S. The simplex method; 2003. 2003-005.
13 Schittkowski K. NLPQLP : A Fortran Implementation of a Sequential Quadratic Programming Algorithmwith Distributed and Non-Monotone Line Search; 2010.

603