A20:
VIBRO-ACOUSTICS

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0014

The effect of particle surface roughness on granular energy
dissipation performance
Furkan TERZIOGLU1; Jem A. RONGONG2; Charles E. LORD3
1

Department of Mechanical Engineering, University of Sheffield, United Kingdom

2

Department of Mechanical Engineering, University of Sheffield, United Kingdom

3

Department of Mechanical Engineering, University of Sheffield, United Kingdom

ABSTRACT
Granular materials can be effectively utilised to reduce structure-borne noise and structural vibrations.
However, since granular damping shows non-linear characteristics and is affected by several different
parameters such as particle shape, a granular damper should be designed considering the factors that may
change granular energy dissipation behaviour. The particles used in a granular damper are generally assumed
as having a certain geometric shape. However, particle surfaces can deviate from the pre-defined shape due
to manufacturing limitations or wear during operation. This experimental study investigates the extent ofparticle surface distortion on the energy dissipation performance of a partially filled granular damper exposed
to harmonic excitation parallel to the direction of gravity. The effect of the particle surface is evaluated overa range of amplitude and frequency conditions. Moulded and 3D printed spherical particles are used toindicate smooth and non-smooth particle shapes, respectively. In order to quantify the particle surface
smoothness, a set of surface roughness measurements is carried out. A dissipated power measurement test rig
is used to assess dissipative performances of smooth and non-smooth spherical particle types in a cylindrical
granular damper enclosure.
Keywords: granular damping, surface roughness, 3D printed particle

1. INTRODUCTION
Relatively small particles placed in a cavity on a vibrating structure produce energy dissipation via
inter-particle interactions. This passive damping technique is called as “granular damping” or “particle
damping” (1). Because of its effectiveness even under challenging environmental conditions, it has
been employed in a wide range of engineering applications to attenuate noise and reduce vibrations(2–5). Nevertheless, it shows significant non-linear characteristics, and the design of granular dampers
depends on many different factors such as clearance, excitation direction, particle shape.
A granular damper is practically a combination of impact damper and friction damper. Thus, it
exhibits the typical dynamic properties of both dampers (6) by dissipating the vibrational energy
through inelastic collisions and frictional interactions within the granular medium. But, the sensitivityof granular damper performance to coefficient of restitution and friction coefficient of individual
particles is somewhat found to be relatively low compared to the other damper parameters (7–9) as
there is an apparent trade-off between the two principal dissipation sources.
However, these findings are based on the assumption of absolute smooth particle surfaces (inparticle-level simulations) although the surface roughness is one of the primary factors that directly
affects the individual particle coefficient of restitution and coefficient of friction (10,11). As the wear
(due to long operational duration) or the manufacturing inconsistencies can cause a level of surface
roughness on the employed particles in the damper, this therefore can lead the possible performancechanges in granular damper designs.
This paper focuses determining the particle surface roughness effect on the overall granular damper
fterzioglu1@sheffield.ac.uk
j.a.rongong@sheffield.ac.uk
3c.lord@sheffield.ac.uk
12

effectiveness conducting experimental studies. Two different particle surface types, whose surface
characteristics are quantified using surface roughness measurements, are used for this investigation.
The individual contact properties of particles are measured to verify the mentioned numerical findingsof the previous studies considering the contact property changes along with the surface roughness.

2. SMOOTH AND NON-SMOOTH PARTICLESIn order to examine the effect of particle surface roughness on the granular energy dissipation, a
moulded smooth particle and a 3D printed non-smooth particle were employed throughout this study.
These particles are illustrated in Figure 1. The material properties of these particles are based onacrylic-based polymer, and the nominal particle diameters are 5 mm.

(a)

(b)

Figure 1 – Investigated spherical particles: (a) smooth and (b) non-smooth.
As can be seen in Figure 1, the smooth particle surface seems perfectly spherical whilst the non smooth particle has apparent manufacturing marks. To determine the difference between the two
surfaces, a set of surface roughness measurements was carried out using the Alicona InfiniteFocus
optical surface measurement machine. The example surface pictures of particle types captured by themachine are demonstrated in Figure 2.
(a)

(b)

Figure 2 – Particle surfaces: (a) smooth and (b) non-smooth.
To quantify each particle surface roughness, the average mean roughness parameter (Ra) was
computed integrating the depth profile of a selected path line on a particle surface, z(x), as:
Ra =

1
Lpath

Lpath

 z ( x ) dx

(1)0

where L path is the total length of the selected path line. The measurement results are summarised inTable 1.
Table 1 – Measured average mean surface roughness values.
Particle

Smooth particle,

Non-smooth particle –Non-smooth particle –noRa [µm]perpendicular to fibres, Ra [µm]on fibre, Ra [µm]10.15 ± 0.02

11.22 ± 0.58

0.34 ± 0.04

2

0.16 ± 0.02

12.47 ± 0.24

0.34 ± 0.03

3

0.17 ± 0.04

11.66 ± 1.04

0.38 ± 0.02

4

0.15 ± 0.02

11.56 ± 0.80

0.35 ± 0.05

5

0.17 ± 0.01

12.96 ± 0.69

0.33 ± 0.04

Average

0.16

11.97

0.35

As can be seen from Table 1, 5 different specimens were considered for each particle type in themeasurements to obtain an average surface characteristic for the particle types. For each specimen, 4different path lines were constructed to determine the variation in the results depending on theselection of path position on the specimen surfaces. Although the directional orientation of this pathselection does not matter for the smooth particle, two different path direct ions were particularly used
for the non-smooth particle to establish its surface roughness because of the effect of printed fibres.
The selection of path direction can be clearly seen in Figure 3 with the exemplified path depth
functions.
z(x)Periodic deep depths
between fibres

The path
perpendicular to
fibre directionality
x

z(x)

x

The path on fibre

Figure 3 – Demonstration of different path directions on the non-smooth particle.
The surface roughness difference between the two different path directions of non -smooth particle
is apparent as can be seen in Table 1. The smooth particle has very low surface roughness comparedto the non-smooth particle.
It should be noted that the sphericity of particles is maintained. Although the non -smooth particle
has relatively large surface roughness values when the fibre-perpendicular paths are considered, theratio of its average mean surface roughness to the diameter is about 0.2%. Hence, for the granular
energy dissipation characteristics of non-spherical particle shapes, the authors address the previous
studies (12–14).

3. COLLISIONAL PROPERTIES OF INDIVIDUAL PARTICLESTo show the effect of surface roughness on the normal collisional behaviour of each introduced
particle, a simple experimental drop test rig was proposed as shown in Figure 4. In this test rig, thetested particle is released from an unknown height to the base plate, and it should be allowed to makeat least 3 successive impacts with the flat plate surface within the measurement sequence. The force
history of each impact is captured by the force transducer mounted between the base plate and a fixed
ground.

Cylindrical base plate
made from acrylicbased polymer

Force transducer

Figure 4 – Simple drop test rig.

A sample measurement which shows 3 subsequent impacts of a smooth particle with the base isillustrated in Figure 5. A basic algorithm was used to identify the impacts. It also calculates the times
between the impacts, T1-2 and T 2-3 , and the contact durations of impacts, ΔT 1, ΔT 2 and ΔT3 .
Impact 1

Impact 2

Impact 3

T1-2

ΔT1

T2-3

ΔT2

ΔT3

Figure 5 – A sample drop test measurement.
In this experimental rig, it was aimed to measure and compare the impact-based dissipation level,i.e., coefficient of restitution (COR), and the contact stiffness of the smooth and non-smooth particles.
The determination of these properties was carried out utilising the theoretical approach presented by
Nagurka and Huang (15). Nagurka and Huang (15) have basically modelled the bouncing behaviourof a spherical ball (i.e., particle) on a flat surface as a single-degree-freedom system as shown in
Figure 6. In this model, m particle is the individual particle mass; kcontact is the contact stiffness; c contact
is the contact damping coefficient determined by the contact damping ratio, ζcontact which is a functionof COR (16).

mparticle

kcontact

ccontact =ζcontact{2(kcontactmparticle)0.5}

Figure 6 – Single-degree-of-freedom model of an individual impact.
Employing this model, COR, e contact is obtained using the time gaps between the captured impacts
in the proposed test rig as:
econtact = T2 −3 / T1− 2(2)The contact stiffness is determined using the measured COR as:
mparticle
 2 + ln 2 ( econtact ) 
kcontact =(3)2 ( Tcontact )where the average contact duration can be found by,Tcontact = ( T1 + T2 + T3 ) / 3(4)
The measured properties are shown in Figure 7 depending on the impact velocity obtained as 2T 12 /g where g is the gravity. Note that the experienced impact velocities are much smaller than the yield
velocity. Therefore, there is no certain effect of impact velocity on the measured COR values.
It is apparent from the presented results that the smooth and the non-smooth particles exhibit very
different collisional properties. The mean COR is about 0.91 and the mean contact stiffness isapproximately 3.9×10 5 N/m for the smooth particle whilst they respectively are 0.83 and 2.1×10 5 N/mfor the non-smooth particle. This shows that increase in particle surface roughness reduces COR andcontact stiffness. It should be noted that this result is consistent with literature (10).
It should be noted that sampling rate, base plate vibration after the loss of contact and particle
releasing are the most significant factors that affect the accuracy of impact capturing in this test rig.
It was observed that the detected contact durations were very sensitive to th ese issues, and, therefore
the contact stiffness values had large uncertainties. However, as the time gaps between subsequent
impacts are insensitive to mis-determination of contact durations, it can be said that the estimation of
COR with this method is reliable despite its simplicity.

Figure 7 – Measured collisional properties: (a) COR and (b) contact stiffness.

4. GRANULAR ENERGY DISSIPATION EFFECTIVENESS OF PARTICLESThe experimental test rig used in this study to measure the dissipated energy of a vibrated granularmedium is shown in Figure 8. It was proposed as analogue to the pioneering Yang’s work (17). As canbe realised from Figure 8, no main structure is present in this test rig. In this way, the dissipated energy
effectiveness can be obtained as independent from the host structure dynamics, and, therefore
consistent comparisons can be conducted at any frequency and amplitude.
The used cylindrical enclosure has the inside diameter of 40 mm and the void height of 40 mm.
This provided 0.47 volume fill ratio as 358 particles were employed in the damper enclosure for bothsmooth and non-smooth particle types. The enclosure material is the same as the particles. It shouldbe also noted that the first natural frequency of damper enclosure is higher than 1200 Hz.
Accelerometer
Enclosure lid
Securing ring

Smoothparticles

Main damper
enclosureNon-smooth
particles
Forcetransducer
ShakerFigure 8 – Granular energy dissipation measurement rig with smooth and non-smooth particles.
The theoretical equations and detailed explanations to calculate the dissipated energy using thistest rig can be found elsewhere (7). The computation procedure is basically based on simultaneous
measurement of developing acceleration and force while the shaker applies harmonic excitation to thedamper enclosure. Thus, the acceleration and force signals were captured for each excitation condition(i.e., amplitude and frequency) considered here using the mounted transducers on the damper
enclosure outer surface as shown in Figure 8.
For each measurement, the granular medium was allowed to reach the steady-state condition while
it was driven by vibration. The acceleration and force signals were then captured for 5 seconds toobtain the average dissipated energy achieved by the damper. This averaged measurement provided toeliminate possible uncertainties that can occur in a single vibration cycle measurement. The usedsampling rate for the signal acquisition was 40 kHz for accurate measurement of the corresponding
signals.
To allow consistent comparisons, the granular damping efficiency parameter (18) was used:2(5)
Damping efficiency = Edissipated 2 / ( 2g ) M particle
where E dissipated is the average dissipated energy in a vibration cycle by the damper; ω is the excitation
frequency in rad/s; Г is the non-dimensional acceleration amplitude determined dividing the excitation
acceleration amplitude by the standard gravity, g; and M particle is the total particle mass employed.
The granular damping efficiency levels of smooth and non-smooth particles are shown in Figure 9
for different excitation cases in order to determine the effect of surface roughness on the damping
under a variety of dynamic conditions of granular medium. For example, the peak in Figure 9a andthe largest peak in Figure 9b are the results of optimum collective collision condition where thegranular medium exhibits the onset of bouncing bed phase. Besides, the smaller peak in Figure 9b andthe peaks in Figure 9c and Figure 9d are responsible to the optimum condition of fluidisation motional
behaviour in the granular medium. More detailed explanations on the relation between the dynamic
granular motion and the achieved granular damping effectiveness can be found in elsewhere (19,20).
(a)

(b)

(c)

(d)

(e)

(f)

Figure 9 – Damping efficiency comparisons for varying vibration amplitudes at (a) 20 Hz, (b) 40 Hz, (c)125 Hz, (d) 160 Hz, (e) 320 Hz and (f) 625 Hz.
As can be seen in Figure 9, the granular damping efficiency levels of smooth and non-smooth
particles are very similar to each other along the whole vibration amplitude axis at each excitation
frequency investigated although there are significant differences between them in terms of surfac eroughness (see Table 1) and individual collisional properties (see Figure 7). Both particle types

captured the same dynamic characteristic behaviours for each excitation case as the same volume fillratio was used for both. This also showed that the obtained ratio of the average mean surface roughness
to the diameter for the non-smooth particle type (i.e., 0.2%) does not affect the overall sphericalparticle shape.
It is somewhat surprising that the smooth particle exhibits more effective granular damping formost of the excitation cases despite its higher COR value than the non -smooth particle. It can be
explained that the smooth particles generate more intensive inter-particle interactions (i.e., larger
relative velocity and larger deformation) within the shaken granular medium which results more
efficient energy dissipation by the particles. This finding also confirms the previous results obtained
from the numerical parameter sensitivity analyses of granular dampers (7,9). It should be noted that
although it was not measured in this study it has been reported elsewhere that the coefficient of friction
grows with increasing surface roughness (11). Therefore, this can be regarded as another factor that
reduces intensity of generated interactions in the granular medium, and, therefore the granular
damping effectiveness decreases as also shown in another parametric study (8).

5. CONCLUSIONSThe presented experimental study has investigated the effect of particle surface roughness in a
granular damper subjected to axial vertical harmonic vibrations. Two apparently different sphericalparticle surface types were considered for this comparative study: smooth and non-smooth. The
quantifying properties of investigated particle types were identified by conducting surface roughness
measurements and simple drop tests. This study has revealed that the smooth particle generally
provides more effective energy dissipation in granular dampers. However, the associated difference
between the employed particle types has been found very small even though the particle types has
exhibited significantly different individual collisional properties. The results have also indicated thatthe amplitude and frequency dependent granular damping characteristic, w hich most-likely depends
on the dynamic motional behaviour of granular medium, is insensitive to the particle surface
roughness as the used surface roughness level has not considerably affected the overall particle shapeand the volume fill ratio. Therefore, there is no major effect of surface roughness on granular damping
performance if the particle surface distortion does not create a certain different geometry from theoriginal particle shape.

REFERENCES
1. Panossian H V. Structural damping enhancement via non-obstructive particle damping technique. J Vib
Acoust. 1992;114(1):101–5.
2. Panossian H V. Non-obstructive particle damping: New experiences and capabilities. In: 49th
AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. Schaumburg,
Illinois; 2008.
3. Simonian SS. Particle damping applications. In: 45th AIAA/ASME/ASCE/AHS/ASC Structures,
Structural Dynamics and Materials Conference. Palm springs, California; 2004. p. 4145–61.
4. Hu L, Shi Y, Yang Q, Song G. Sound reduction at a target point inside an enclosed cavity using particle
dampers. J Sound Vib. 2016;384:45–55.
5. Xu Z, Yu Wang M, Chen T. A particle damper for vibration and noise reduction. J Sound Vib. 2004;270(4–5):1033–40.
6. Mao K, Wang MY, Xu ZZ, Chen T. Simulation and characterization of particle damping in transient
vibrations. J Vib Acoust. 2004;126(2):202–11.
7. Wong CX, Daniel MC, Rongong JA. Energy dissipation prediction of particle dampers. J Sound Vib.
2009;319(1–2):91–118.
8. Lu Z, Masri SF, Lu X. Parametric studies of the performance of particle dampers under harmonic
excitation. Struct Control Heal Monit. 2011;18(1):79–98.
9. Wang Y, Liu B, Tian A, Tang W. Experimental and numerical investigations on the performance of particle
dampers attached to a primary structure undergoing free vibration in the horizontal and vertical directions.
J Sound Vib. 2016;371:35–55.
10. Sandeep CS, Luo L, Senetakis K. Effect of grain size and surface roughness on the normal coefficient ofrestitution of single grains. Materials (Basel). 2020;13(4).
11. Li X, Dong M, Jiang D, Li S, Shang Y. The effect of surface roughness on normal restitution coefficient,adhesion force and friction coefficient of the particle-wall collision. Powder Technol. 2020;362:17–25.

12. Terzioglu F, Rongong JA, Lord CE. The dissipative characteristics of oblate particles in granular dampers.
In: Proceedings of the International Conference on Structural Dynamic , EURODYN. 2020. p. 4851–66.
13. Pourtavakoli H, Parteli EJR, Pöschel T. Granular dampers: Does particle shape matter? New J Phys.
2016;18(7).
14. Sánchez M, Carlevaro CM, Pugnaloni LA. Effect of particle shape and fragmentation on the response ofparticle dampers. J Vib Control. 2014;20(12):1846–54.
15. Nagurka M, Huang S. A mass-spring-damper model of a bouncing ball. In: Proceedings of the American
Control Conference. Boston, Massachusetts; 2004. p. 499–504.
16. Schwager T, Pöschel T. Coefficient of restitution for viscoelastic spheres: The effect of delayed recovery.
Phys Rev E - Stat Nonlinear, Soft Matter Phys. 2008;78(5):1–12.
17. Yang MY. Development of master design curves for particle impact dampers. The Pennsylvania State
University; 2003.
18. Sack A, Heckel M, Kollmer JE, Zimber F, Pöschel T. Energy dissipation in driven granular matter in theabsence of gravity. Phys Rev Lett. 2013;111(1):1–5.
19. Terzioglu F, Rongong JA, Lord CE. Motional phase maps for estimating the effectiveness of granular
dampers. Mech Syst Signal Process. 2022;Submitted.
20. Terzioglu F, Rongong JA, Lord CE. Construction of motional phase maps for granular dampers. In:
INTER-NOISE 2022 - 51st International Congress and Exposition on Noise Control Engineering: NoiseControl in a More Sustainable Future. Glasgow; 2022.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022, in Gyeongju, Korea
Extended abstractABS-0068

Study on the vibration attenuation properties of the S-type
phononic crystals
Tongtong Huo1; Zhijing Wu2; Fengming Li3
1

College of Aerospace and Civil Engineering, Harbin Engineering University, China

2

College of Aerospace and Civil Engineering, Harbin Engineering University, China

3

College of Aerospace and Civil Engineering, Harbin Engineering University, China

ABSTRACTThe center wavelength of the first band gap for a Bragg phononic crystal corresponds to the elastic wave
length approximately twice of the lattice period length. It is difficult for Bragg phononic crystals to obtain
the band gap in a low frequency range with a small cell size. How to design a unit cell to obtain wider andlower band gaps is a research topic worth exploring. This paper proposed a kind of S-type phononic crystal
and investigated a new kind of arrangement of multiple unit cells to widen the band-gap frequency range.
The S-type phononic crystal folds the unit cells to shorten the thickness and achieves the vibration isolation
effect. The S-type structure also could connect different unit cells to extend the band gaps. The numericalresults of the proposed structure are calculated by the softwar COMSOL Multiphysics. The band gaps andfrequency responses of the structure are calculated based on the finite element method. Furthermore, the Stype phononic crystal is printed and the vibration experiment is launched. The experimental and simulation
results are compared to verify the vibration isolation performance of the structure. This thin phononic crystal
structure provides a technical way for the application of phononic crystals in low-frequency vibration andnoise reduction.
Keywords: Phononic crystals; S-type; COMSOL Multiphysics; Vibration experiment

1. INTRODUCTION
Phononic crystal [1,2] is the periodic structure composed of one or more materials withacoustic/elastic waves control. In some frequency ranges, the vibration could propagate in the
structure, these frequency ranges are called pass bands. But the vibration and noise attenuation inother frequency ranges and called band gaps. Phononic crystals has attracted the attention of many
researchers due to the band gap characteristics. These structure have a wide application prospect in
the field of engineering [3,4] .
The research on the phononic crystals was gradually deepening [5,7]. The numerous phononic
crystals proposed mainly focus on the design of unit cell to obtain a wider and lower band gap. This
paper investigates the arrangement and combined with multiple unit cells to widen the band gaps
frequency range.

2. ModelsIn this section, Fig.1(a) displays the arrangement of traditional 1D phononic crystal, and the unit
cells are arranged periodically in one direction. The corresponding unit cell is displayed in Fig.1(c).
It contains two rectangles of l_5*l_1 and l_6*l_4 . A hole is digged in the left one. The length of thehole is l_7,and the height is l_3. As shown in Fig.1(b), a different arrangement of phononic crystal is
proposed,which called “s” structure. Each layer is connected through a rectangle with the width of a.
The thickness of the structure is b. The periodic structure can produce Bragg band-gaps.

1
2

httsoar@163.com
wuzhijing@hrbeu.edu.cn

(a)

(c)

(b)

Fig. 1. (a) Traditional 1D phononic crystal, (b) “S” structure, (c) the unit cell

3. ResultCompared the frequency response curves of the different structure. As shown in Fig.2, the “S”
structure retains the band gaps of the original structure. The “S” structure could shorten the length,
but it still achieves the effect of vibration reduction.

normal structure"S" structure0-2

Acceleration (m/s2)

-4
-6
-8
-10
-12
-14
-16
-18
0

1000

2000

3000

4000

5000

6000

Frequency(Hz)Fig. 2. The “S” structure and normal structure frequency response curves.

ACKNOWLEDGEMENTS
This research is supported by the National Natural Science Foundation of China (Nos. 12072086
and 11761131006) and the Research Team Project of Heilongjiang Natural Science Foundation (No.
TD2020A001).

REFERENCES
1.

Wang Y F, Wang Y Z, Wu B, Chen WQ, Wang YS. Tunable and active phononic crystals andmetamaterials. Applied Mechanics Reviews, 2020, 72: 040801.
2. Liu Z Q, Zhang X, M ao Y, Zhu Y Y, Yang Z, Chan C T, Sheng P. Locally resonant sonic materials.
Science, 2000, 289: 1734–1736.
3. He J X, Li H L, Tian Y H, Zhang Q Z, Lu Z X, Lan J X. Numerical analysis of viscous dissipation in
microchannel sensor based on phononic crystal. Micromachines, 2021, 12(8): 994.
4. Li L J, Gang X Y, Sun Z Y, Zhang X X, Zhang F. Design of phononic crystals plate and application in
vehicle sound insulation. Advance in Engineering Software, 2018, 125: 19–26.
5. Wu Z J, Li F M, Zhang C Z. Band-gap analysis of a novel lattice with a hierarchical periodicity usingthe spectral element method. Journal of Sound and Vibration, 2018, 421: 246–260.
6. Wen S, Xiong Y, Hao S, Li F, Zhang C. Enhanced band-gap properties of an acoustic metamaterial beamwith periodically variable cross-sections. International Journal of Mechanical Sciences, 2019,
166(13):105229.
7. Hao S, Wu Z, Li F, Zhang C. Numerical and experimental investigations on the band-gap characteristics
of metamaterial multi-span beams. Physics Letters A. 2019, 383: 126029.

ABS-0174

Investigation on the vibration behavior of elastic porous materials forminiaturized systems
Zenong CAI(1) , Vicente Cutanda HENRIQUEZ(1) , Aminul ISLAM(2) , Frieder LUCKLUM(1)(1) Centre for Acoustic-Mechanical Microsystems (CAMM), Department of Electrical and Photonics Engineering, Technical

University of Denmark (DTU), DK-2800, Kongens Lyngby, Denmark, zenca@dtu.dk(2) Centre for Acoustic-Mechanical Microsystems, Department of Civil and Mechanical Engineering, Technical University ofDenmark, DK-2800, Kongens Lyngby, Denmark

ABSTRACTIn this study, we investigate eight kinds of porous materials with different geometric parameters, densities, andporosity, and the influence of these variables on their mechanical behavior is discussed. The vibration behaviorof beams made of different porous materials is dynamically characterized using a laser vibrometer, and themobility data of these porous materials are calculated. This study can be beneficial for finding suitable elastic
porous materials for different vibroacoustic applications in miniaturized systems.

Keywords: Vibration, Vibroacoustics, Porous material, Mechanical mobility

1

INTRODUCTION

Well-designed elastic porous materials can be used for vibration control in miniaturized systems due to theirdifferent and more complex behavior on elastic wave attenuation than nonporous materials.
Former investigation [1] finds that the dynamic moduli of polyurethane foams are frequency dependent, because of their constituent materials’ viscoelasticity [2], so dynamic stiffness and damping should be used in theapplication of sound enclosure structures with foams rather than static stiffness. Many different methods canbe used for obtaining dynamic mechanical material properties over a wide range of frequency [3], however,
the research about the relation between micro-geometric structure of porous materials and their vibroacoustic
dynamics behavior is not sufficient.
In this study, eight materials with different geometric parameters are selected and experimentally studied dynamically. In the second section, the parameters of the selected materials are described, followed by the introductionof the measurement system. The third section discusses the initial experimentally measured mobility results.

2

MATERIAL AND METHOD

2.1 Materials for experiments
This section describes the basic parameters of the materials used for the experiments, as shown in Table 1.
The geometric parameters are measured from microscopic images. From the different structures, the eight materials can be classified into three types: Bubble-like pores, reticulated structure, and closed pores. Among them,the microscopic images of types A, B, and C are shown in Figure 1(a) to Figure 1(c), and the photo of all thetested beams is Figure 1(d). All the skeletons of eight samples are polymer materials.

Table 1. Parameters of Materials
Length of the beam(mm)

Thickness(mm)

Density(kg/m3 )

0.5

196

9

27.6

0.5

0.5

197

18

27.6

Reticulated structure No.1

1

0.1

201

18

8.8

Reticulated structure No.2

1

1

199

9

24.5

Reticulated structure No.3

1

0.7

199

24

87.9

Closed pores No.1

0

0.3

201

9

40.6

Closed pores No.2

0

0.5

201

24

51.5

Closed pores No.3

0

0.2

199

6

102.4

Type

Open Porosity Ratio

Bubble-like pores No.1

0.5

Bubble-like pores No.2

(a) Type A Bubble-like pores

Average poresize (mm)(b) Type B Reticulated structure(c) Type C Closed pores(d) photo of tested beams

Figure 1. Photos of samples for the tests. (a) (b) (c) are micrographs of pores type A, B, and C, respectively.

2.2 Test system
In this study, as shown in Figure 2, we use the Pulse system (3560, Brüel & Kjær A/S) to collect data. The
porous materials are cut into the shape of beams and screwed at the center on a shaker (4810, Brüel & KjærA/S), which provides external excitation (white noise, 500m Vrms). A force transducer (8200, Brüel & KjærA/S) is between shaker and beam to measure force data. The force transducer needs an amplifier (NEXUS,
Brüel & Kjær A/S) to link with Pulse. A laser vibrometer (PDV-100, Polytec GmbH) is used to measure thevelocity at the center of the beam, in other words, measures the velocity of the wave amplitude of the bending
wave which comes from the shaker.

Figure 2. Sketch of the test system.
The shaker is clamped to the support so that the beam can be measured in a horizontal position to the ground,as shown in Figure 3.

Figure 3. The position of shaker, beam, force transducer, and laser measurement point.

3

RESULT AND DISCUSSION

The mechanical mobility of eight beams is calculated by dividing the force signal by the velocity signal. The
results are shown in Figure 4. Firstly, From Figure 4, the beams with the same kind of structure share theFigure 4. The position of shaker, beam, force transducer, and laser measurement point.

similar shape of the curve, except for reticulated structure No.3 and closed pores No.3, which have too large
thicknesses, so the bending waves don’t really propagate through those beams. It means that the different pore
structures influence the dynamic mechanic behavior.
For reticulated structure, it’s interesting that there is an anti-resonance with wide bandwidth around 1500 Hz,
which cannot be seen in the curves from other structures. It might be because the solid phase of the reticulated
structure is constructed by interconnecting thin strips, and others are constructed by interconnecting faces, nomatter if the pore is open or closed. And this difference causes the wave propagation in the solid phasedifferent. However, the details of this anti-resonance should be investigated in the future. Besides, between the
curves of closed pores No.1 and No.3, because of larger density, No.3 has more resonance frequencies than
No.1.
Comparing bubble-like pores and closed pores, the significant difference is after 1500 Hz, the curves of bubblelike pores cannot show more resonance, but the curves of closed pores are still clear. So initially we can saythat the ratio of open-cell and closed-cell also affects the dynamic mechanical behavior of porous materials.
Notably the results include the influence of force transducer, which introduces extra mass to the system andneeds to be subtracted in the future. After the subtraction, the loss factor and Young’s modulus will be calculated, according to the resonance peak values of the curves, which will provide more information.

4

CONCLUSION

In this study, according to the initial vibration tests and the results of mobility, we can conclude that the
microstructure of the porous material has a significant influence on its dynamic mechanical properties, which
need to be investigated in the future. More dynamic mechanical parameters will be calculated according to theresults of this paper and further experiments.
In addition to the study shown in this paper, we also investigate the acoustic performance of the same group of
porous materials. Combined with another related research on the mechanical performance of porous material in
acoustic measurements, this research provides the basis for our work on the design and optimization of porous
materials for microacoustic and vibroacoustic applications on millimeter-scale devices, for example, hearing aids.

REFERENCES[1] Kingsbury HB, Cho KH, Powers WR. Dynamic complex modulus of polyurethane foams. Journal of Cellular Plastics. 1978 Mar;14(2):113-7.
[2] Lee WE. Cellular solids, structure and properties. Materials Science and Technology. 2000 Feb 1;16(2):233.
[3] Jaouen L, Renault A, Deverge M. Elastic and damping characterizations of acoustical porous materials: Available experimental methods and applications to a melamine foam. Applied acoustics. 2008 Dec
1;69(12):1129-40.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0400

Waveguide absorbers based on spiral acoustic black holes
for damping vibration in a plate structure
Seongmin PARK1; Jae Yeon LEE2; Wonju JEON1
1

Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Republic of
Korea
2

1st R&D Center, Agency for Defense Development, Republic of Korea

ABSTRACT
We propose a waveguide absorber (WGA) based on a spiral acoustic black hole (ABH) to reduce vibrationin plates. The geometrical dimensions of the WGA are determined so that it has a lower cut-on frequency
than a target frequency by using the impedance method. We also determine the attachment location anddirection of the WGA on a plate to maximize its damping performance based on the structural intensity
analysis. A couple of WGAs attached to a plate could reduce the mobility of the plate by more than 10 dB.
In addition, the WGAs had far less weight but could achieve larger peak reductions compared to thickdamping material or spiral ABHs attached to edges.
Keywords: Acoustic black hole (ABH), Vibration damping, Elastic wave

1. INTRODUCTION
The acoustic black hole (ABH) is a wedge structure from which an elastic wave does not escape.
When the thickness of the wedge structure is smoothly decreased to zero by the power-law ℎ(𝑥) =𝜀𝑥 𝑚 (𝑚 ≥ 2), the group velocity and phase velocity of incident waves approache to zero and elastic
wave energy is highly localized near the wedge tip (1). Many researchers have conducted theoretical,
numerical, and experimental studies on the ABH to use it in vibration damping, noise re duction,energy harvesting. (2-6). Recently, a concept of spiral ABH was proposed by Lee and Jeon (7) toreduce the space occupied by the ABH, and its damping performance was validated experimentally
by Park et al. (8). In addition, the cut-on frequency of the spiral ABH has been thoroughly investigatedby developing an impedance method, which provided the possibility to design the spiral ABH wit h adesired cut-on frequency (9).
Although the spiral ABH is a compact solution for damping vibration in thin beams, it is bulky andheavy when the host structure is a plate. There have been several attempts such as using ABH-based
periodic structures (10, 11) for vibration suppression using lightweight devices, but these attempts
decreases the structural integrity of host structures, so it is difficult to use the structures in practice.
Therefore, in this study, we propose another type of spiral ABH that has far less volume and weightcompared to the original spiral ABH. The proposed spiral ABH is mounted locally on the plate surface,
and it is utilized as a waveguide absorber (WGA) for vibration damping of thin plates.

2. RESULTS AND DISCUSSIONS
Figure 1(a) describes four WGAs based on spiral ABHs attached to a simply supported host platewith a uniform thickness of 9 mm and edge lengths of 700 mm and 6 mm. Each WGA is composed ofa spiral ABH with a gap distance of 12 mm and an arc length of 800 mm, a damping layer with lengt hof 80 mm and a thickness of 3 mm, and a connecting rod with length of 30 mm and height of 20 mm.
The widths of the spiral ABH, the damping layer, and the connecting rod are 30 mm. The materials
for the spiral ABH and the damping layer are Aluminum 6061 T6 and butyl rubber, respectively.
1
2

kaistian13@kaist.ac.kr
jaeyeonlee@add.re.kr

Detailed explanations on determining the values of geometrical parameters, location and direction of
attachment are in reference (12).
The surface-averaged mobilities of the host structure with and without the four WGAs are shownin Fig. 1(b). Large peak reductions about 12.2−21.0 dB could be achieved by using the four WGAs,
compared to the single WGA. The peak reductions for the WGAs are compared to those of damping
material treatment and spiral ABHs attached to plate edges, as shown in Fig. 1(b). For damping
treatment, damping material with the thickness of 3 mm is covered on the top surface of the plate. The
spiral ABHs attached to edges have the same geometrical dimensions as the WGA. When comparedto the case of damping material treatment (1.2kg), vibration in the host plate is more dampened byusing the spiral ABHs (11.5kg) or WGAs (0.6kg). The WGAs are locally attached on the plate, andthose are much lighter than the spiral ABHs attached to the edges.

Figure 1 – (a) Four WGAs based on spiral ABHs attached to a host plate and (b) the surface-averaged
mobilities of the host plate and the plate with four WGAs, damping material, and two spiral ABHs

REFERENCES
1. M. A. Mironov, Propagation of a flexural wave in a plate whose thickness decrease smoothly to zero in a
finite interval, Sov. Phys. Acoust., 34 (1988) 318-319.
2. V. V. Krylov and F. J. B. S. Tilman, Acoustic ‘black holes’ for flexural waves as effective vibration dampers,
J. Sound Vib., 274 (2004) 605-619.
3. E. P. Bowyer, V. V. Krylov, Experimental study of sound radiation by plates containing circular indentationsof power-law profile, Appl. Acoust., 88 (2015) 30-37.
4. L. Zhao, S. C. Conlon, F. Semperlotti, An experimental study of vibration based energy harvesting in
dynamically tailored structures with embedded acoustic black holes, Smart Mater. Struct., 24 (2015)065039.
5. J. Y. Lee and W. Jeon, Exact solution of Euler-Bernoulli equation for acoustic black holes via generalized
hypergeometric differential equation, J. Sound Vib., 452 (2019) 191-204.
6. S.-W. Seo, J. Y. Lee, K.-J. Kim, and W. Jeon, Wave-based analysis of dual acoustic black holes for anechoic
termination of shock testing devices, Wave Motion, 95 (2020) 102468.
7. J. Y. Lee and W. Jeon, Vibration damping using a spiral acoustic black hole, J. Acoust. Soc. Am. 141 (2017)
1437-1445.
8. S. Park, M. Kim, and W. Jeon, Experimental validation of vibration damping using an Archimedean spiralacoustic black hole, J. Sound Vib., 459 (2019) 114838.
9. J. Y. Lee and W. Jeon, Wave-based analysis of the cut-on frequency of curved acoustic black holes, J.
Sound Vib., 492 (2021) 115731.
10. S. Park, W. Jeon, Ultra-wide low-frequency band gap in a tapered phononic beam, J. Sound Vib. 499(2021) 115977.
11. L. Tang, L. Cheng, K. Chen. Complete sub-wavelength flexural wave band gaps in plates with periodic
acoustic black holes. J. Sound Vib., 502 (2021) 116102.
12. S. Park, J. Y. Lee, and W. Jeon, Vibration damping of plates using waveguide absorbers based on spiralacoustic black holes J. Sound Vib., 521 (2022) 116685.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022, in Gyeongju, Korea
Extended abstractABS-0417

A study on intuitive connection between structural mode andradiated sound power
Dongjoon Kim1; Jung-Han Woo2; Pyung-Sik Ma2; Yun-Ho Seo2
1

Korea Institute of Machinery & Materials, Korea

2

Korea Institute of Machinery & Materials, Korea

ABSTRACTThe radiated sound power is one of the frequently used indicators for evaluating of sound radiation
characteristics of structure. For precise measurement of sound power, both sound pressure and velocity shallbe measured. But when the structure is simple such as plate-like structure, the radiated sound power can becalculated by finely discretizing the plate as elemental sound radiators and measuring structural vibration
based on discrete calculation method (DCM). Then, the relationship between the calculated sound power andmeasured structural vibration can be expressed with some indicators: radiation efficiency, radiation mode,
structural radiation mode. These indicators help to understand the sound radiation characteristics of structure,
but there is a limit to use when we try to reduce radiation noise through structural modification as thoseindicators also change. Also, those indicators have less accuracy of calculation and limit of the available
frequency range. In this paper, an indicator called Q-factor is introduced which intuitively connect structuralmodes and radiated sound power. Q-factor can correspond to various input components and the size of the
Q-factor represents the contribution of each structural mode to the sound radiation, thereby helping to moreeasily identify the sound radiation characteristics of the structure.
Keywords: Sound radiation, structural mode, discrete calculation method

1. INTRODUCTION
Sound radiation from plate-like structures such as ship hull and deck, frame of washing machine
has been studied for a long time. To reduce the unintended sound radiation, excitation sources andcomposing plates’ vibration characteristics should be analyzed. Changing excitation sources’ location
and working frequency is an easiest way, but in most cases, excitation sources such as engine, turbine,
motor can’t be moved to other location. Many engineers have been studied to calculate far field soundpower using the object’s surface velocity. DCM uses the discretized plate elements’ velocity byassuming as baffled pistons and radiation impedance to calculated sound powe r (1). Also, there weremany studies to calculate sound power of a rectangular plate having simple boundary conditions
numerically (2–5). Calculated sound power can be an indicator to analyze the vibroacoustic
characteristics of the structure, but it does not tell us how to reduce the acoustic noise or how tomodify the structures to reduce the sound. Aforementioned studies (1–5) also calculated the radiation
efficiency dividing the sound power by spatially averaged squared velocity of the structure . The
radiation efficiency not only represents the structure’s sound radiation efficiency, but also represents
the relationship between the structure’s vibration and radiated sound power. This can be a huge hint
for finding structural modes contribute highly to sound radiation, so many studies used the radiation
efficiency and have been conducted on faster calculation of radiation efficiency using asymptotic andempirical formulas (4–7). Similarly, radiation mode calculated by performing eigen-analysis ofradiation resistance matrix (8,9) informs the radiation characteristics at specific frequency. The loworder radiation modes are dominant at low frequencies, so adopting a structural modification that can
suppress the low-order mode is effective in reducing noise.
1

djkim@kimm.re.kr

However, the radiation efficiency and radiation mode have limits in high frequency range. Since theradiation efficiency of the radiation modes is close to one, all radiation modes contribute equally tosound radiation. Also, the radiation efficiency varies depending on not only structure’s dimension but
also location of excitation source, which means that using derived empirical formulas of averaged
radiation efficiency results in a large error. In this point of view, aforementioned indicators are notfriendly indicator for structural modification to reduce noise. In this paper, Q-factor is introduced tocover the disadvantages of existing methods. Q-factor is an indicator that can directly evaluate thestructural mode’s contribution on sound radiation, which utilizes actual structural characteristics
rather than virtual indicators like radiation modes and does not change depending on the lo cation of
excitation source unlike radiation efficiency. Derivation of Q-factors are presented and comparison
with existing methods are also described in below.

2. Derivation of radiation efficiency, radiation mode, and Q-factor
Normal velocity 𝑣 of a thin plate with uniform thickness can be expressed as
∞

∞

𝑛=1

𝑛=1

𝑗𝜔𝐹𝜙𝑛 (𝑥, 𝑦)𝜙𝑛 (𝑥0 , 𝑦0 )𝑣(𝑥, 𝑦) = ∑= ∑ 𝑎𝑛 𝜙𝑛 (𝑥0 , 𝑦0 )2 (1 + 𝑗𝜂) − 𝜔 2 ]𝑀𝑛 [𝜔𝑛,𝑛(1)

where 𝜔 is frequency, 𝜙 is mode shape vector, 𝑀𝑛 is modal mass of nth mode, 𝜂 is damping loss
factor, (𝑥0, 𝑦0 ) is input location, 𝜔𝑛 is the natural frequency of the structure, and 𝑎𝑛 is modal
velocity amplitude. Using eq. (1), force mobility 𝑌𝑣𝐹 can be expressed as∞𝑣(𝑥, 𝑦)𝑗𝜔𝜙𝑛 (𝑥, 𝑦)𝜙𝑛 (𝑥0 , 𝑦0 )𝑌𝑣𝐹 ==∑2 (1 + 𝑗𝜂) − 𝜔 2 ]𝐹(𝑥0 , 𝑦0 )𝑀𝑛 [𝜔𝑛,𝑛(2)𝑛=1Also, Sound radiation from a baffled plate can be expressed as𝑠𝑠
𝑊 = Re{𝑣 𝐻 𝑝} = Re{𝑣 𝐻 𝑍𝑣} = 𝑣 𝐻 𝑅𝑣
2
2

(3)

where W represents the sound power, 𝑠 represents elemental radiator area, Re represents the realpart, 𝑝 represents pressure matrix, 𝐻 represents the Hermitian operator, 𝑍 represents the acoustic
impedance matrix, and 𝑅 represents the radiation resistance matrix. Diagonal and off-diagonalcomponents of 𝑅 matrix are defined as
𝑅𝑖𝑖 (𝑘) = 𝜌𝑐𝑠𝑖 (1 −𝑅𝑖𝑗 (𝑘) =𝐽1 (2𝑘𝑟))
𝑘𝑟

2𝜌𝑐𝑘 2 𝑠𝑖 𝑠𝑗 𝐽1 (𝑘𝑟𝑖 ) 𝐽1 (𝑘𝑟𝑗 ) sin 𝑘𝑑𝑖𝑗
[
]
𝜋
𝑘𝑟𝑖
𝑘𝑟𝑗
𝑘𝑑𝑖𝑗

(4)
(5)

where 𝜌 is density of a plate, 𝑐 is sound speed in air, 𝐽1 is first-order Bessel function, 𝑘 iswavenumber, 𝑑𝑖𝑗 is the distance between element 𝑖 and 𝑗 , 𝑟 is equivalent radius of elemental
radiator.
Using eq. (1) and (3), radiation efficiency 𝜎 can be derived as𝜎=𝑊̅̅̅2 〉
𝜌𝑐𝑆𝑇 〈𝑣(6)̅̅̅2 〉 is spatially averaged squared velocity of a plate. If 〈𝑣̅̅̅2 〉
where 𝑆𝑇 is total area of a plate and 〈𝑣is calculated using only a single mode, calculated 𝜎 becomes a structural mode’s radiation efficiency.
By performing eigen-analysis of 𝑅 matrix, radiation mode 𝜑 can be derived and eq. (3) can berearranged as
𝑊 = 𝑣 𝐻 𝑅𝑣 = 𝑣 𝐻 (𝜑𝐻 Λ𝜑)𝑣

(7)

where Λ represents eigenvalue and 𝜑 represents eigenvector of 𝑅 matrix. Similar to structuralmode, radiation efficiency of each radiation mode can be expressed as
𝜎𝑖 =

2
Λ
𝜌𝑐𝑆𝑇 𝑖𝑖(8)

where 𝑖 represents the order of radiation mode.
By combining eq. (1) and (3), 𝑊 can be expressed in a new form as (10)𝑊 = 𝑣 𝐻 𝑅𝑣 = (𝑎𝐻 𝜙 𝐻 )𝑅(𝜙𝑎) = 𝑎𝐻 𝑃𝑎 = (Φ𝐻 Ω𝐻 )𝑃(ΩΦ) = Φ𝐻 𝑄Φ

(9)𝑄 = Ω𝐻 𝜙 𝐻 𝑅𝜙Ω(10)

where Ω and Φ areΩ𝑛 (ω) =𝑗𝜔
2 (1 + 𝑗𝜂) − 𝜔 2 ]𝑀𝑛 [𝜔𝑛,𝑛Φ𝑛 (𝑥0 , 𝑦0 ) = 𝜙𝑛 (𝑥0 , 𝑦0 )𝐹

(11)
(12)

From eq. (9), 𝑄 represents Q-factor of a plate. Diagonal term of 𝑄 represents each structural mode’scontribution on sound radiation.

3. Numerical evaluation
In this section, numerical investigation on a simple rectangular finite plate is described.
Specification of plate model is shown in table 1.
Table 1 – Specification of plate
𝐿

60 mm

𝑊

40 mm

𝑇

1 mm

𝑚𝑠

1700 kg/m3

𝐸

13.47 GPa

𝜈

0.384

𝜂

0.04

The origin of coordinate system was set as the center of the plate and the grid spacing was set as 1
mm. Then, radiation efficiency of the plate can be calculated as Figure 1. (x, y) in the legend represents
the location of input force. When the location of input force varies, the radiation efficiency also varies
significantly as shown in Figure 1. Also, the calculated radiation efficiency using empirical formula
shows a significant difference in most frequency range. It can lead to an overestimation orunderestimation of the radiated sound power. In addition, when plate is excited by multiple sources,
above expressions are not applicable as the sum between radiation efficiency is not possible. In thisFigure 1 – Radiation efficiency variation according to the location of excitation source ( 𝐹=1)
perspective, radiation efficiency is a good indicator to predict sound power, but it can be said that it
is not effective in terms of accuracy and utilization.
Figure 2 shows the normalized radiation mode shapes of a plate. In case of low frequency (as shownin Figure 2. (a)), the radiation mode shapes are similar to low-order structural mode shapes such as
rocking motion and piston-like motion. Also, their radiation efficiencies as shown in Figure 3
represents the effectiveness on sound radiation. In low frequency, low-order radiation mode dominates
the sound radiation and one can get a hint on how to modify the structure for sound reduction. However,

Figure 2 – Normalized radiation mode shapes of a plate at (a) 𝑘𝑙=1.36 and (b) 𝑘𝑙=15.6

Figure 3 – Radiation efficiencies of radiation modes of a plate

Figure 4 – Q-factor of a plate at (a) 𝑘𝑙=1.36 and (b) 𝑘𝑙=15.6
in case of high frequency (as shown in Figure 2. (b)), the radiation mode shapes become more complex
and similar to high-order structural mode shapes. Furthermore, similar to Figure 1, the radiation
efficiency of each radiation modes becomes close to one. It means that every radiation mode
contributes equally to sound radiation, so it becomes difficult to judge intuitively which part needs
modification. Also, radiation mode calculation requires a large amount of memory and calculation
time to perform eigen-analysis for each frequency.
Figure 4 represents the Q-factor as 3d bar graph. In low frequency range as shown in Figure 4. (a),only a low-order structural mode has meaningful Q value. This means that 1st structural mode
contributes largely to sound radiation. In case of high frequency range as shown in Figure 4. (b), morebars with different values appeared. This means that many high-order structural modes have involved
in sound radiation and each mode’s degree of participation is expressed by the Q value. To sum up,

Q-factor intuitively connects mode shape and radiated sound power, which can be used as an indicator
to find the structural modes participating in sound radiation and to determine the degree ofparticipation of those modes.

4. Conclusion--

Radiation efficiency and radiation modes can be an indicator for pretending the relationshipbetween the structural vibration and radiated sound power.
However, those indicators can only be utilized effectively in low frequency range , and there isa limit to intuitively understanding the relationship between structural modes and soundradiation
Q-factor can be calculated multiplying resistance matrix by structural mode shape vectors.
Q-factor inform the participating structural modes and their participation level in soundradiation.
Q-factor is more intuitive than radiation efficiency and radiation modes.

ACKNOWLEDGEMENTS
This work was supported by Korea Institute of Machinery & Materials (NK238B).

REFERENCES
1.

Hashimoto N. Measurement of sound radiation efficiency by the discrete calculation method. Appl
Acoust. 2001;62(4):429–46.
2. G. Maidanik. Response of ribbed panels to reverberant acoustic fields. J Acoust Soc Am.
1962;34(1962):809–26.
3. Leppington FG, Broadbent EG, Heron KH. The Acoustic Radiation Efficiency of Rectangular Panels.
Proc R Soc Lond A Math Phys Sci. 1982;382(1783):245–71.
4. Xie G, Thompson DJ, Jones CJC. The radiation efficiency of baffled plates and strips. J Sound Vib.
2005;280(1–2):181–209.
5. Wallace CE. Radiation Resistance of a Rectangular Panel. J Acoust Soc Am. 1972;51(3B):946–52.
6. Putra A, Thompson DJ. Radiation efficiency of unbaffled and perforated plates near a rigid reflecting
surface. J Sound Vib [Internet]. 2011;330(22):5443–59. Available from:http://dx.doi.org/10.1016/j.jsv.2011.05.033
7. Leppington FG, Broadbent EG, Heron KH. Acoustic Radiation From Rectangular Panels With
Constrained Edges. Proc R Soc London, Ser A Math Phys Sci. 1984;393(1804):67–84.
8. Elliott SJ, Johnson ME. Radiation modes and the active control of sound power. J Acoust Soc Am.
1993;94(4):2194–204.
9. Ji L, Stuart Bolton J. Coupling mechanism analysis of structural modes and sound radiations of a tire
tread band based on the S-mode technique. Appl Acoust [Internet]. 2015;99:161–70. Available from:http://dx.doi.org/10.1016/j.apacoust.2015.06.013
10. Kim D, Park N. Calculation and reduction of sound radiation from a thin plate structure excited bycomplex inputs. J Sound Vib [Internet]. 2020;484:115517. Available from:https://doi.org/10.1016/j.jsv.2020.115517

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0511

Equivalent comfort contours for foot-transmitted vibration
generated from uniaxial accelerator pedal simulators
Junsun YOO1; Seonbin LIM1; Dongju YOON2; No-Cheol PARK1
1
2

Department of Mechanical Engineering, Yonsei University, South Korea

Research & Development Division, Hyundai Motor Co., Ltd, South Korea

ABSTRACT
Vehicle vibration not only causes discomfort in the driver and passengers, but also increases mental as well
as physical fatigue. Therefore, so-called vibrotactile sensation test is being performed to develop vehicle
components based on the driver’s reaction and improve riding comfort by reducing the driver’s discomfort
due to the vibration and noise occurred in vehicles. The vehicle accelerator pedal is constantly contacted by
the driver while driving, and it frequently generates various vibration. In this study, a total of three uniaxial
accelerator pedal simulators were built to conduct a vibrotactile sensation test with 28 participants. From thetest, the sensitivity of the foot sole of each participant to the vibration magnitudes at frequencies generated
from the uniaxial accelerator pedal simulators was extracted. And equivalent comfort contours, representing
the magnitudes of vibration required to give the identical vibrotactile sensation to the participants in thefrequency range, were derived through the quantification of the sensitivity of the foot sole. The equivalent
comfort contours include the perception thresholds of the participants for foot-transmitted vibration. It is
expected that the equivalent comfort contours can be applied to the vehicle accelerator pedal vibration topredict whether drivers feel discomfort.
Keywords: Foot-transmitted vibration, Vibration perception threshold, Equivalent Comfort Contours

1. INTRODUCTION
Vehicle vibrations not only have an adverse effect on the driving performance and safety byincreasing the physical and mental fatigue of drivers, but also lower the vibrotactile quality of vehicles
by causing discomforts [1]. Therefore, automotive industries have been focused on improving thevibrotactile quality of vehicle components. They, recently, have performed vibrotactile sensation testwith respect to the body part being stimulated by vibrations generated from vehicle componen t [2, 3].
The results derived by quantifying and analyzing the subjective data obtained from the participantsfor vibrotactile sensation test can be reflected in the design of vehicle component [4, 5]. The vibrationstransmitted from the accelerator pedal, brake pedal, steering wheel, and seat that are in physicalcontact with the driver in the vehicle greatly affect the riding comfort. Among them, there are manystudies on the steering wheel because it influences the greatest effect on the driver’s riding comfort[6, 7]. On the other hand, there has not much study on the accelerator pedal vibration even though it
constantly transmits vibration to the foot during driving.
This paper presents equivalent comfort contours, representing vibrotactile sensation perceived bydrivers in regard to vibrations transmitted from the vehicle accelerator pedal to the sole of the foot.
To produce the equivalent comfort contours, vibrotactile sensation test were carried out using threeuniaxial pedal vibration simulators reflecting vehicle accelerator pedal. Vibration magnitudes atvarious frequencies occurring in the vehicle were applied to the foot sole of the participants for thevibrotactile sensation test through the pedal vibration simula tors. By quantifying the subjective
responses of the participants for vibrotactile sensation test to vibrotactile stimuli, equivalent comfort
contours considering the sensitivity of the foot sole are derived. The equivalent comfort contours
indicate the actual vibration magnitude required for drivers to have the same vibrotactile sensation ateach frequency. The equivalent comfort contours can contribute to the vehicle accelerator pedal to be
12

pnch@yonsei.ac.kr
junsunyoo0118@naver.com

developed in the future in consideration of the discomfort perceived by drivers according to themagnitude of the foot-transmitted vibration for each frequency.

2. PERCEPTION THRESHOLD
2.1 Procedure
As shown in Fig. 1, three uniaxial accelerator pedal vibration simulators were made to perform thevibrotactile sensation test for perception threshold with an up-down (staircase) algorithm [22]. The
participants for the vibrotactile sensation test were classified into three groups: expert, nonexpert, andoverall. Expert group is made up of people who conduct vibration-related research, and nonexpert
group is made up of random people who have not experienced for vibration -related research. All
participants for the vibrotactile sensation test were instructed to wear sneakers and apply a consistentlevel of pedal effort to the accelerator pedal vibration simulator. In addition, they were not exposedto the accelerator pedal vibration for a long time in order not to feel tiredness or numb in the feet.

(a)

(b)

(c)

Figure 1 – Uniaxial accelerator pedal vibration simulators for (a) X-axis (b) Y-axis (c) Z-axis
2.2 Results and Discussion
Figure 2 shows the perception thresholds of the participants for the foot-transmitted vibration from theaccelerator pedal vibration simulator were determined depending on experimental group. The perception
thresholds were produced based on the medians of the data obtained from the vibrotactile sensation test to
exclude outliers. The perception thresholds allow to figure out the frequency range perceived by participants
as relatively sensitive or insensitive to the foot-transmitted vibration. It can be seen that the perception
thresholds are located in the lowest in the frequency range between 63 and 125 Hz, meaning that theparticipants perceive most sensitive at that frequency range.

(a)
(b)
(c)
Figure 2 – Perception thresholds of (a) expert, (b) nonexpert, and (c) overall for foot-transmitted vibrationfrom accelerator pedal vibration simulator

3. EQUIVALENT COMFORT CONTOURS
3.1 Procedure
The vibrotactile sensation test for equivalent comfort contours was performed with the same participants
participating in the experiment for perception threshold. In this experiment, two sinusoidal vibrations in Table
1 were applied separately with an interval of 2 seconds to the foot sole of the participants through theaccelerator pedal vibration simulator. For the experiment, the vibrotactile sensation due to the magnitude of
the test vibration was set to exceed the vibrotactile sensation of the perception threshold. By using Eq. (1)based on the Stevens’ power law, the objectivity of the subjective data acquired from the experiment was
secured and the equivalent comfort contours for foot-transmitted vibration were derived.
Table 1 – Vibrations used in the vibrotactile sensation test for equivalent comfort contours
Type of vibration
Frequency [Hz]Magnitude [dB]
FeatureReference vibration

100

108

20, 25, 31.5, 40, 50,
Test vibration

63, 80, 100, 125,

90 ~ 120

160, 200

Sinusoidal vibration
including cosinetapered ends

1

𝜓 𝑛
(1)
𝜑 − 𝜑0 = ( )𝑘where 𝜑 is the vibration magnitude and 𝜑0 is the perception threshold. On the right side of theequation, 𝜓 is the vibrotactile sensation magnitude. n and k are vibrotactile sensation growth and aconstant, respectively.
3.2 Results and Discussion
Equivalent comfort contours, as shown in Fig. 3, demonstrate the quantified vibrotactile sensation
perceived by the participants depending on the magnitude of foot-transmitted vibration for eachfrequency based on the sensitivity of the foot sole. In other words, equivalent comfort contours
represent the vibration magnitude actually needed in order to give the same vibrotactile sensation to
the drivers.
A low magnitude in equivalent comfort contours means a great sensitivity to the foot-transmitted
vibration at that frequency and magnitude. In the frequency range between 63 and 125 Hz, participants
perceived the foot-transmitted vibrations most sensitively. Thus, the participants have the greatest
discomfort in those frequency range. On the other hand, they become insensitive to the foottransmitted vibrations as towards the lower or higher frequencies.

(a)
(b)
(c)
Figure 3 – Equivalent comfort contours for foot-transmitted vibration through a Z-axis accelerator
pedal vibration simulator depending on experimental groups: (a) expert, (b) nonexpert, and (c) overa ll

4. CONCLUSIONSIn this study, three uniaxial accelerator pedal vibration simulators were made to investigate thesensitivity of the foot sole from the vibrotactile sensation test with respect to the vibration magnitudes
at frequencies. Then, equivalent comfort contours were derived by quantifying the sensitivity of thefoot sole. The equivalent comfort contours indicate the actual vibration magnitude needed for drivers
to have the same vibrotactile sensation at frequencies, and therefore eligible to discover the frequency
range where the drivers perceive as relatively sensitive or insensitive to the foot-transmitted vibration.
Moreover, it can be known a specific frequency and the corresponding vibration magnitude at which
the drivers begin to perceive the vibration being transmitted to their foot. Based on the equivalent
comfort contours, the vehicle accelerator pedal to be developed in the future can be improved
considering the discomfort perceived by drivers depending on the magnitude of the foot -transmitted
vibration at each frequency.

ACKNOWLEDGEMENTS
This work was supported by the Hyundai Motor Company, South Korea, and by the Materials andComponents Technology Development Program (20011013), funded by the Ministry of Trade,Industry & Energy (MOTIE, South Korea).

REFERENCES
1. Zhou Y, Chen S. Vehicle ride comfort analysis with whole-body vibration on long-span bridges subjected
to crosswind. Journal of Wind Engineering and Industrial Aerodynamics. 2016;155:126-140.
2. Morioka M, Griffsrfin M J. Equivalent comfort contours for vertical vibration of steering wheels: effectof vibration magnitude, grip force, and hand position. Applied ergonomics. 2009;40(5):817-825.
3. Morioka M, Griffin M J. Magnitude-dependence of equivalent comfort contours for fore-and-aft, lateral,and vertical vibration at the foot for seated persons. Journal of Sound and Vibration. 2010;329(14): 29392952.
4. Mansfield N, Sammonds G, Nguyen L. Driver discomfort in vehicle seats–Effect of changing roadconditions and seat foam composition. Applied ergonomics. 2015;50:153-159.
5. Beard G F, Griffin M J. Discomfort of seated persons exposed to low frequency lateral and roll oscillation:
Effect of backrest height. Applied Ergonomics. 2016;54:51-61.
6. Ye S, Hou L, Zhang P, Bu X, Xiang J, Tang H, Lin J. Transfer path analysis and its application in
low-frequency vibration reduction of steering wheel of a passenger vehicle. Applied Acoustics.
2020;157:107021
7. Yoo W S, Na S D, Kim M S. Relationship between subjective and objective evaluations of steering wheel
vibration. Journal of mechanical science and technology. 2011;25(7):1695-1701.
8. Morioka M, Griffin M J. Magnitude-dependence of equivalent comfort contours for fore-and-aft, lateraland vertical hand-transmitted vibration. Journal of sound and vibration. 2006;295(3-5):633-648.

ABS-0585

Vibroacoustic characterization of foamed and fiber reinforced plastic parts
Arne ROTERMUND(1) ; Steffen HOFFMANN(1) ; Jörg HAIN(2) ; Sabine C. LANGER(1)(1) Institut für Akustik, Technische Universität Braunschweig, Germany, a.rotermund@tu-braunschweig.de(2) Bio-Materials and Polymers, Volkswagen AG, Group Innovation, Wolfsburg, Germany

ABSTRACT
Thermoplastics are widely used materials in the automotive industry and various industrial applications. The
injection molding process enables fast and cost-effective production of products or parts made of thermoplastics. Reinforcing the plastics with glass-fibers can enhance mechanical properties like the specific stiffness andtherefore lead to a potential lightweight construction. Foaming can be used in the injection molding process tocreate integral foams with solid skin layers and a porous core which can further increase the specific stiffness
and additionally leads to less warpage and better thermal insulation. The vibroacoustic behavior of a part madefrom integral foam also differs from that of a compact injection molded part. Besides the inhomogeneous density distribution, foaming leads to different fiber orientations and consequently to a change of the anisotropic
viscoelastic properties. Increasing requirements for acoustic comfort lead to a demand for reliable numericalmodels and input parameters to predict the vibroacoustic behavior of these parts. In this paper, the characterization of the dynamic mechanical properties of beam specimens of integral foams from polyvinyl chloride andglass-fiber reinforced polypropylene is presented.

Keywords: Foam, Fiber reinforced plastics

1

INTRODUCTION

Injection molding is a technology widely used. Pure thermoplastics as well as reinforced thermoplastics (e.g.
with glass fibers) are cost-effective and fast to process for many products, especially in the automotive industry.
With an increasing demand for a lightweight design in order to reduce emissions foaming of parts becomes
more attractive. Plastic foams can be found in various products like thermal insulations or sound absorbers
for decades. The microcellular injection molding (MIM) technology has already shown its potential in some
applications. The MuCell R process is a MIM technology that has already been used for interior parts for carsof different manufacturers with a resulting weight reduction of up to 20% [1, 2].
Foaming of fiber-reinforced thermoplastic parts can lead to a reduction of weight and material used with theadditional potential to increase dimensional accuracy, thermal and acoustic insulation as well as an improved
impact strength [3]. As found in [1] MIM leads to a displacement of fibers from the parts center towards theparts surfaces. Fiber concentrations and global fiber orientation therefore differ from compact parts. This change
also influences the mechanical properties [4] and is likely to lead to a varying vibroacoustic behaviour.
The process investigated in [1], as well as in this paper, is the MIM process Ku-FizzTM , developed at
Volkswagen AG. The process is described in detail in [1] and a weight reduction between 5 and 16% compared
to compact material is reported there for plate specimens made of glass fiber reinforced polypropylene with
30 wt% glass fibers.
While some MIM processes, like MuCell R , have already been studied for a while, Ku-FizzTM differs from
MuCell R in some major aspects which is why the Ku-FizzTM process and its influence on the produced parts
need to be studied further. One major difference is the different pressure operating window for the physicalblowing agent. Ku-FizzTM employs a much lower gas pressure of 5 to 30 bar compared to MuCell R with80 to 200 bar [5]. Depending on the actual gas pressure, Ku-FizzTM can create a comparable or a different
foam and fiber microstructure compared to MuCell R [1]. Another major difference between the two MIM

processes is that in the Ku-FizzTM process, gas is added directly into the feed hopper. This enables a good
homogeneity of the polymer melt with the diffused gas controlled by the screw speed and temperature controlof the tools. These differences enable the Ku-FizzTM process to achieve different skin layers compared to other
MIM processes, as was shown in [4]. These differences of Ku-FizzTM from other MIM processes lead to thequestion how much the dynamic mechanical properties are influenced by them.
In this paper the dynamic mechanical properties are investigated through the analysis of the vibroacoustic
behaviour of specimens produced with the Ku-FizzTM process. Therefore beam specimens are cut out of plate
specimens (geometry as shown in [1]) and investigated with an extended Resonance-curve method based on
ISO 6721-3 [6]. The material investigated in this work is a foamed glass fiber reinforced polypropylene with
20 wt% glass fibers. The rigid integral foam (IF) specimens can be considered as a sandwich-structured material
with dense surface layers and a porous core layer. As a reference, another polymeric structural foam thatshows typical sandwich behaviour is compared with the material produced with Ku-FizzTM . Therefore polyvinyl
chloride (PVC) foam specimens with a varying foam structure are measured.

2

MATERIAL PARAMETER ESTIMATION

The method for the determination of dynamic mechanical properties of plastics under flexural vibration used inthis paper is based on the Resonance-curve method in ISO 6721-3 [6] and is extended in the context of thispaper. The aim of the method is to determine the flexural storage modulus Efi0 of a vibrating beam specimen atseveral resonance frequencies fri . Efi0 is calculated with Eq. 1:
Efi0 =

2  2√fri
4π 3ρ l 2
t
ki2

(1)

In addition to the density ρ, length l and thickness t, a factor ki is needed that depends on the mode. The
calculation of ki , the flexural loss modulus Efi00 and the loss factor ηfi is further described in [6].
In order to excite the specimens with an electromagnetic field, a light metal plate with a mass of 0.02 g isapplied to the specimens. The beams are mounted on textile threads and their position in relation to the end of
the beam lm is chosen so that the influence on the deflection shapes is as small as possible. After measuring
Laser vibrometer

Electric magnet

Textile threads
Metal plate
Foam beam sample
z

z
y

lm
x

l

lm

x

w

t

y

Figure 1. Experimental setup for material property determination according to ISO 6721-3 [6]the surface velocity Lv (Resonance-curve) at one end of the beam as shown in Fig. 2, the areas with occurringresonances are selected. To determine both the loss factor ηfi and the resonance frequencies fri , the peak fitting
method is applied, which is also well suited for samples with high structural damping. Therefore a surrogate
model (Eq. 2) [7, p. 48]) for the single mass oscillator is fitted to the selected data using the method of leastsquares.
ni
H( jΩ) = 2≈ Lv(2)
ωi − Ω2 + j2ηi ωi Ω
Since the resonance frequencies of the light weight foam specimens are shifted by the mass of the metal plates
applied at the beams ends, these are to be considered in the material characterization. Especially for the thickspecimens, the shear deformation must not be neglected neither. In order to consider both aspects in the theory
for material characterization, a model based on the Timoshenko-Ehrenfest beam theory presented in [8] is used.
In the mentioned paper the frequency equation (Eq. 3) for a Timoshenko beam with two masses placed at theends is derived.
√2  2
r + Λβ1 β2 − r2 + s2 sin λ1 sinh λ22ξ Λ (cos λ1 cosh λ2 − 1) − ξ 2 4 + ξ 2 s2 − r2
r





− (β1 + β2 ) −Λ 4 + ξ 2 (s2 − r2 )2 λ1 λ12 − s2 ξ 2 sin λ1 cos λ2 − λ2 λ2 + s2 ξ 2 cos λ1 sinh λ2 = 0
To calculate the frequency equation, dimensionless quantities
r
r
r
ρA
Ef I
1
1 I
2
ξ = ωl
, s=
, r=
,
Ef I
l κGA
l A

(3)

Λ = 1 − Ω2 r2 s2

and the two parameters λ1 and λ2 are required for the case ξ > 1/rs
s
s
p
p
2
2
2
2
2
2
2
ξ (r + s ) + ξ 4 + ξ (s − r )ξ 2 (r2 + s2 ) − ξ 4 + ξ 2 (s2 − r2 )2, λ2 =
λ1 =22where A is the cross-section of the beam, I the moment of inertia, ρ the density and l the length of thespecimen. The shear modulus G and shear correction factor κ are assumed to be as follows:κ=

10(1 + ν),(12 + 11ν)G=

Ef
2 (1 + ν)

Throughout the paper a Poisson’s ratio ν of 0.3 is used. Although proper material parameters depend on theposition along the beam height, length and width, only homogenized values for the beam specimens are used inthis paper. The square roots of the frequency equation represent the natural angular frequencies. The material
parameters are determined with an inverse minimization problem using least squares. The error between themeasured fri and calculated natural frequencies fTi is used as residuum.
Nfri − fTi
frii=1

r( f ) = ∑The starting value (flexural storage modulus Ef0 ) for the fitting is calculated using Eq. 1. In order to savecomputational costs, the frequency vector ω is preconditioned to the bands around the measured resonancefrequencies. Where the width of the frequency bands into higher and lower frequencies are a quarter of the
distance to the adjacent peak. For the root determination in the pre-configured frequency bands, Brent’s method[9] is used. Tab. 1 shows Ef0 calculated according to Eq. 1 and with the presented method (Eq. 3) based on thefirst eight natural frequencies of a FEM simulation and their percentage deviation from FEM input parameters(Ef0 : 1300 MPa, l: 0.5 m, t: 10 mm, w: 20 mm, ρ: 510 kg/m3 ). The FE model is a beam geometry modeled
in Abaqus.The elements are 20-node brick elements with quadratic ansatz functions and reduced integration(C3D20R). The mounting is not considered.
Table 1. Method comparison to FEM simulation (Eq. 1 and Eq. 3)
Ef0 , ISO 6721-3 (Eq. 1)
err in %Ef0 , Timoshenko (Eq. 3)
err in %

1297
-0.2
1300
0.0

1294
-0.47
1300
0.0

1286
-1.08
1298
-0.18

1283
-1.29
1302
0.16

1274
-2.02
1302
0.12

1265
-2.68
1304
0.28

1255
-3.5
1305
0.42

1242
-4.43
1307
0.54

101

19

76

1
0

0

1

2

Loss factor ηf
0.05

1.3

0.04

1.2

0.03

3

0

f in kHz1
2
f in kHz

ηf

59

La

125

Ef0 in GPa

Storage modulus Ef0
1.4

La calculated98

fr [%]Lv

Lv measured

3

Figure 2. Measured frequency response of a PVC foam specimen mounted in mode 3
Tab. 1 shows a underestimation of Ef0 up to 4.5% between the ISO 6721-3 (Eq. 1) and the FEM input
parameters. The method based on Timoshenko–Ehrenfest beam theory shows an overestimation of 0.5% for thesame mode order. All the following analyses are therefore only performed with the presented method.
As an example for the measurement of a beam sample, the measured surface velocity is shown in Fig. 2
over the frequency. On the right, the flexural storage modulus Ef0 and the flexural loss factor ηfi are shown
determined by the presented method. Based on these material parameters an FE simulation is carried out. In
the model the mass of the metal plate is considered as a point mass. The deviation between the measured andcalculated (FEM) resonance frequencies ∆ fr is used as an evaluation criterion. This is shown in Fig. 2 in thebottom on the left. In the following, only flexural storage modules Ef0 and associated flexural loss factors ηfi
are presented that show a resonant frequency deviation of less than 2% in a comparison of simulations andmeasurements.

3

EXPERIMENTAL CHARACTERIZATION

The materials under research in this paper are foamed glass fiber reinforced polypropylene (POLYfill PPH
GF5020 HC - PP-GF20 - Polykemi AB) and different foams of PVC. The PP-GF20 is produced with the
Ku-FizzTM process with different process parameters. The gas pressure of the physical blowing agent varied
between 7.5, 10 and 20 bar. The melt temperature varied between 220 ◦C and 230 ◦C. PVC is compared inform of free foam (FF) (Forex® by 3A Composites GmbH), which is a foam with a homogeneous structure,and IF (Vekaplan S by the VEKA AG), which is a foam with an inhomogeneous structure consisting of solid
skin layers and a porous core. The different foam structures of PVC FF and IF specimens are shown in
Fig. 3. Microscopic images of specimens produced with Ku-FizzTM were shown in [1] and [4]. The PVC FF
is measured with different thicknesses (Tab. 2) and compared with a PVC IF. The PVC FF and IF specimens
all have a much higher foaming rate (Approx. 64-71%) than the PP-GF20 IF specimens (6-8%).
Table 2. PVC & PP-GF20 beam specimens
Material
PVC
PVC
PVC
PVC
PVC
PP-GF20
PP-GF20

type
FF
FF
FF
FF
IF
IF (H)IF (V)

Thickness, mm
3
4
5
10
10
3.4
3.4

Length, mm
275
275
275
500
500
380
96.7

Width, mm
10
10
10
20
20
10
10

Density, kg/m3
517
417
441
450
510
968
968

1 mm
(a) 3 mm PVC FF(b) 5 mm PVC FF(c) 10 mm PVC FF(d) 10 mm PVC IF

Figure 3. Microscopic images of the foamed PVC samples of different thicknesses and types of foaming (allpictures are taken with the same microscopic settings and share the scale shown in (a))
All beam specimens are cut out of plate specimens. In case of PP-GF20 the beams are cut out at different
locations of the plate with an orientation parallel to the melt flow (in the following described as horizontal andabbreviated by H) and an orientation perpendicular to the melt flow (in the following described as vertical andabbreviated by V), which is shown in Fig. 7. As presented in [1], in the Ku-FizzTM process fibers are displaced
through the cell growth from the center to the surface of the plate. Together with the typically varying fiber
volume fractions of fiber reinforced injection molded parts this leads to a slightly varying density from thespecimens on the outside to the ones of the inside of the plate. The density variations are in both horizontal
and vertical direction. The dimensions and densities (apparent densities) of the specimens are shown in Tab. 2
as mean values over all specimens. The length of a few of the PP-GF20 specimens had to be slightly reduced
from their initial length in order to eliminate defects in the specimens. The lengths of the PVC specimens has
been chosen that the length to thickness ratio is at least 50 for all specimens as is recommended in [6]. For
the PP-GF20 specimens the length is limited by the plates geometry. A length to thickness ratio of 50 cantherefore only be achieved with the horizontal but not withe the vertical vertical specimens which is why an
implementation of the consideration of shear has been chosen as described in section 2. The experiments are
conducted at the ambient temperature of the laboratory which is 19 ◦C.
Fig. 4 shows the flexural storage moduli Ef0 and flexural loss factors ηf calculated with the measured frequency responses from left to right for PVC FF with an increasing thickness of 3, 4, 5 and 10 mm as well as
PVC IF with a thickness of 10 mm on the right. The parameters of the specimens are shown in Tab. 2. All
measured foam samples do not show a logarithmic increase of storage modulus over frequency typical for PVC

Ef0 in GPa

1.5
1.3
1.1

0.9
0.05

3 mm FF

4 mm FF5 mm FF

10 mm FF

10 mm FF

ηf

0.04
0.03
0.02

0

1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3f in kHzf in kHzf in kHzf in kHzf in kHz

Figure 4. Flexural storage modulus Ef0 and loss factors ηf of PVC free-foam (FF) and integral-foam (IF) specimens different thicknesses (see Tab. 2).

H7
6.5

5.4
0.05

7.5 bar 230 ◦C

10 bar 230 ◦C

20 bar 230 ◦C

Specimens

Ef0 in GPa

7.6

20 bar 220 ◦C

ηf

0.04
H1

0.03
0.02

0

1
2
f in kHz30

1
2f in kHz30

1
2f in kHz30

1
2f in kHz3

Figure 5. Flexural storage modulus Ef0 and flexural loss factors ηf of PP-GF20 IF beam specimens (in melt
flow direction) (see Tab. 2).
and other thermoplastics, as shown in [10]. With increasing material thickness, there is a linear decrease instorage modulus over frequency for the FF specimens. This effect becomes even clearer when looking at theright plot in Fig. 4. The IF also shows a linear decrease of the flexural storage modulus versus frequency witha significantly higher gradient. This linear decrease of bending stiffness over frequency is often found in sandwich structures [11]. In this case the linear decrease in the flexural storage modulus of the sandwich structure
is superimposed on the logarithmically increasing material stiffness of the compact material.
In the following it is investigated if the PP-GF20 IF specimens with a much lower foaming rate show a similar behaviour. In order to additionally investigate the location-dependent material behavior of the microcellular
injection molded PP-GF20 plate structure, horizontal and vertical specimens are produced, as shown in Fig. 7.
Fig. 5 shows the determined flexural storage moduli and loss factors of the horizontally cut out specimens fordifferent injection pressures and temperatures with the colorbar indicating the position of the specimen. Regardless of pressure and temperature, Fig. 5 shows similar curves for the loss factors for all specimens. The material
stiffness, however, shows significantly higher values for the specimens located at the edge (H1 and H7). Here, a

V30
3.5

ηf

2.5
0.08

7.5 bar 230 ◦C

10 bar 230 ◦C

20 bar 230 ◦C

Specimens

Ef0 in GPa

4.5

20 bar 220 ◦C

0.06
0.04

V1
0

1

2 3 4
f in kHz

5

0

1

2 3 4
f in kHz

5

0

1

2 3 4
f in kHz

5

0

1

2 3 4
f in kHz

5

Figure 6. Flexural storage modulus Ef0 and flexural loss factors ηf of PP-GF20 IF beam specimens (perpendicular to melt flow direction) (see Tab. 2).

stronger orientation of the fibers in the injection direction presumably takes place in the manufacturing process,which contributes to an increased bending stiffness in the specimens. Also, as the density of the H1 and H7
specimens is slightly higher, the fiber concentration is assumed to be higher based on the findings of [1] wherea fiber displacement towards the mold surfaces is reported.
A comparison of the right plot with the other plots in Fig. 5 indicates that this effect depends less on the
foaming pressure and more on the processing temperature. The effect of increasing stiffness towards the edgeof the flat structures due to fiber orientation is also shown for the different process executions on the left handside in Fig. 7.
H1
H2

H6
H7
7.3

V29
V30

Ef0 in GPa

Ef0 in GPa

4.2

6.4

3.4

5.5
0.04

2.6
0.07

30
V

25
V

20
V

15
V

10
V

1
H
2
H
3
H
4
H
5
H
6
H
7

H

V

0.05

0.02

5

0.06

1

0.03

V

ηf

ηf

V1
V2

Flow direction
7.5 bar 230 ◦C
10 bar 230 ◦C
20 bar 230 ◦C
20 bar 220 ◦C

Figure 7. Averaged values of Ef0 and ηf of the first 3 modes of Fig. 5 and Fig. 6
The results of the specimens vertically cut out of the plate structure are shown in Fig. 6. The specimens also
exhibit location-dependent material parameters. Additionally, there is a correlation between the material stiffness and the loss factors determined. The effect of increasing stiffness towards the edge of the flat structures(especially for the V1-V4 specimens near the gate) can also be seen for the four different process parameter combinations. If the fiber concentration or orientation (and thus also the bending stiffness) increases, thehomogenized loss factors of the composite material decrease. This relationship can also be seen in averaged
quantities for the first 3 modes on the right side in Fig. 7. All four process parameter combinations lead to thesame effect that the highest flexural moduli were measured in the specimens close to the gate of the injection
molding process. Based on the findings of [1], it is assumed that is due to the varying fiber orientations, thatare slightly influenced by the process parameters and differ significantly at different positions of the plate.

4

CONCLUSIONSIn this paper, the Resonance-curve method extended with the consideration of shear deformation and additional
mass is used for the determination of dynamic flexural properties of foamed materials. The error comparedto the formula for for flexural modulus of ISO 6721-3 can be significantly reduced by using the Timoshenko
theory presented in [8].
The flexural storage modulus and the flexural loss factor of foamed plastics were considered for different
PVC foams. It is shown that the homogenized flexural modulus is dependent on the structure thickness in theinvestigated specimens as it influences the structures density and most likely the proportion of the skin layer
thickness to the core layer thickness. Comparing the 3 mm FF specimens with the 10 mm specimens, which havea similar density, the samples of the FF show a typical logarithmic increase of the flexural storage modulus overthe frequency while the linear decrease of the flexural storage modulus of the IF seems to indicate sandwich
effects. The linear decrease in the flexural storage modulus of the sandwich structure is superimposed on thelogarithmically increasing material stiffness of the starting material.
The IF specimens produced with the Ku-FizzTM process do not exhibit a significant typical sandwich behaviour as is seen with PVC foam specimens. The rather low density reduction compared to the compact
specimens might be a reason for this. However, the varying fiber orientations and concentrations throughout
the plate reported in [1] have a significant influence on the dynamic mechanical properties of the material. The
decreasing fiber orientation perpendicular to the flow with growing distance to the injection gate is agrees wellwith the decreasing flexural storage modulus.

ACKNOLEDGEMENTS
This research was part of the ACTion project funded by the Federal Ministry of Education and Research of
Germany (BMBF). The authors would like to thank the BMBF for funding the research on integral foams
produced with the Ku-FizzTM process.

REFERENCES[1] Simon SA, Hain J, Osswald T. Effect o gas pressure on the microstructure of parts foamed with the novel
microcellular injection molding technology Ku-FizzTM . SPE Polymers, Vol. 2 (4), 2021, p. 311-324.
[2] Heitkamp H, Betsche M. More than Just Foam Bubbles, Kunststoffe international, 12/2014, p. 14-19.
[3] Jeon B, Han JW, Lee KS, Cha SW. Improvement of the Mechanical Properties of Biodegradable Polymers
Using a Microcelluar Foaming Process and Natural By-Products, Polymer-Plastics Technology and Engineering, Vol. 51 (4), 2012, p. 401-406.
[4] Gómez-Monterde J, Hain J, Sánchez-Soto M, Maspoch MLl. Microcellular injection molding: A comparison
between MuCell process and the novel micro-foaming technology IQ foam. Journal of Materials Processing
Technology, Vol. 268, 2019, p. 162-170.
[5] Anderson JR, Cha SW, Laing DE, Pierrick DE, Stevenson JE. Injection molding of microcellular material.
EP0952908, 1998.
[6] ISO 6721-3:2021, Plastics — Determination of dynamic mechanical properties — Part 3: Flexural vibration— Resonance-curve method, 2021.
[7] VDI 3830-5: Werkstoff- und Bauteildämpfung - Blatt 5, Verein Deutscher Ingenieure, Beuth Verlag, Berlin,
2005.
[8] Shi W, Shean ZB, Peng XL, Li XF. Frequency equation and resonant frequencies of free–free Timoshenko
beams with unequal end masses, International Journal of Mechanical Sciences, Vol. 115-116, 2016, p. 406415.
[9] Brent RP. Algorithms for Minimization without Derivatives, Prentice-Hall, Englewood Cliffs, New Jersey,
1973.
[10] Lagakos N, Jarzynski J, Cole JH, Buraro JA. Frequency and temperature dependance of elastic moduli ofpolymers, Journal of Applied Physics, Vol. 59 (12), 1986, p. 4017-4031.
[11] Nilsson AC, Baro S, Piana EA. Vibro-acoustic properties of sandwich structures, Applied Acoustics, Vol.
139, 2018, p. 259-266.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022, in Gyeongju, Korea
Extended abstractABS-0682

A Study on the Wave Characteristics of the Periodic Structures for
the Vibration Attenuation in a Rotating Ring
Kyu-Beom Kim 1; Jeong-Guon Ih1
1

Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Korea

ABSTRACTIn this work, the wave characteristics are studied on the periodic structure attenuating the vibration in arotating ring, which is often the main power transmission element of many machines. The rotating ring
structure is modeled as an arc-shaped Euler-Bernoulli beam, on which the periodic structures in the form ofstiffeners are distributed. The propagation constant of the wave is obtained by solving the eigenvalue problem,
utilizing the transfer matrix method and Floquet’s theorem. When the stiffeners are periodically attachedalong the ring, the stopbands of flexural or longitudinal wave propagation appear at the resonance frequencies
of a unit cell in a periodic arrangement. In the case of circular rings, flexural and longitudinal waves arecoupled at some narrow frequency bands, which produces additional stopbands. When it comes to a rotating
ring, the wave propagation characteristics depend on the forward or backward directions because of the
Coriolis force. This effect causes asymmetric stopbands for the forward and backwards waves. Such acharacteristic helps increase the bandwidth of the stopbands of the propagating structural waves due to theexcitation on the rotating ring with periodic structures.
Keywords: Periodic structure, Rotating ring, Structural wave, Wave coupling, Stiffener, Stopband

1. INTRODUCTION

This study is concerned with the elastic wave, such as flexural and longitudinal ones, propagation
on a ring and the resulting vibration response of it. When a traveling wave propagates in thecircumferential direction of a circular ring structure, it can be interpreted as an infinite periodic
structure when a periodic substructure is attached along the circumferential direction. This study
investigates the effect of such a periodic structure on its vibration attenuation effect. The rotating ring
structure of a machine is simplified into a thin circular beam, and the periodic substructure issimplified into a periodically attached stiffener.
The mechanical impedance changes periodically along the propagation direction in the periodic
structure, as described by Bragg scattering. As well known, Bragg scattering results in wave
attenuation whenever the offset interference of waves happens due to the repeated structural pattern
for a frequency band, called the stopbands. Such an offset interference occurs when the length of aunit cell corresponds to an integer multiple of the half wavelength. In stopband regions, resonanceoccurs within a propagation cycle, causing waves to be trapped near the discontinuity. While a periodic
structure has passbands, in which waves can propagate without obstruction [1,2]. Such a phenomenon
was also described in beams and plates with periodic discontinuities [3-5]. Various methods have beenused for analyzing the wave propagation in a periodic structure, for example, the transfer matrix
method, plane wave expansion method, wave finite element method, etc. [6-8]. Floquet's theorem isalso employed, assuming an infinitely repeated structure.
Wave propagation in the ring propagates indefinitely along the ring as far as the medium is lossless.
Simultaneously, as a locally finite structure between discontinuities, resonance occurs when the waves
of equal phases repeatedly meet at a point. In addition, the flexural wave and the longitudinal wave
interact with each other in the circular ring. The wave propagation characteristics in a rotating ringdiffer with propagation direction due to the presence of the Coriolis force [9]. Two different uniquenatural frequencies exist for each mode. A difference occurs due to the Doppler effect in reference toa fixed coordinate system because the waves propagate in a rotating coordinate system [10]. When1

kbkim759@kaist.ac.kr, jgih21@gmail.com

flexural and longitudinal waves interact within a periodic structure, wave coupling occurs, resulting
in an additional stopband [11-13]. Unlike resonance, the phase difference of a stopband caused bywave coupling is not a multiple of π between one periodic interval.
In this work, we analyze the characteristics of a structure presented as a simplified thin curved
beam-shaped ring and periodic stiffeners. The effect of the asymmetric stiffeners on a ring is studied,which enables the interaction of flexural and longitudinal waves. If the phase difference does not
happen in the stopband due to asymmetrical stiffener and curved beam characteristics, the resonancedoes not occur, resulting in a significant vibration reduction. When such a circular structure also
rotates, it is exciting and practical to see the effects of wave propagation in both directions and whether
the stopband, due to the periodic structure, would become asymmetrical.

2. ANALYSIS METHOD
2.1

Transfer Matrix Considering the Periodic Stiffeners

Figure 1-Example geometry of a beam and stiffener set constituting a unit cell. Symmetric and asymmetricstiffeners can be seen, and the spatial periodicity of stiffeners is L.
Figure 1 illustrates an element including a stiffener, which is a unitary part of the periodic structureof concern. The transfer matrix T of each element between two adjacent stiffeners represents therelationship between the wave variables v L on the left-hand side and v R on the right-hand side, of
which their relationship is given as v R = Tv L . Each wave variable represents the displacement andforce associated with a specific propagating wave. The transfer matrix Tcell of a unit cell can beobtained by multiplying the transfer matrices of the elements occupying a basic periodic system asfollows:v ≡ [ξ η θ

Q M ] ≡ qTTN=
v ( x0 + L) Tcell v=( x0 ), Tcell

T

FT  ,

n

∏T .
i =1

(1)
(2a,b)i

Here, x0 is the origin point where a unit cell begins, L denotes the unit cell length, and Ti is thetransfer matrix of an element of the whole structure. Because such a unit cell repeats infinitely, wave
propagation characteristics in the periodic structure can be easily obtained using the matrix Tcell .
The stiffener is assumed to be a point mass with roational inertia. Also, it is assumed that there isno deformation within an element. Based on such assumptions, the transfer matrix can be obtained byreferring to the boundary conditions of an element as
 1
 0

 0
Tstiff =  2
 −ω ms 0
 2
 ω ms r01
00−ω 2 ms
000
12
ω ms rs0−ω 2 I zs

0
0
0
1
0
0

0
0
0
0
1
0

0
0 
0
.
0
0

1 

(3)

Here, ms and I zs denote the mass and moment of inertia of the stiffener, respectively, and rs is thecenter distance between the adjacent stiffeners. Note that T34 and T61 are nonzero in the asymmetricstiffener, which means that the interaction between flexural and longitudinal waves exists.
The transfer matrix of a curved beam element can be determined from governing equations
following the Kirchhoff-Love plate assumptions. To satisfy the assumptions, the thickness of the beam
must be less than λ 6 and R 10 , where R means the radius of curvature of the curved element [10].
The governing equations can be converted into a form of the state equation as follows:
1

 0
R

1
 −
0

R

 0
0
dv
= Λv
=
,Λ 
dx

 −ω 2 ρ A
0−ω 2 ρ A
 0

 0
0

0

1
EA

1

0

0

0

0

0
1
R
0

0 −
0


0 

0
0 


1 
0
EI  .

1
0 
R


0
0 

−1 0 
0(4a,b)The terms including 1 / R indicate that the interaction between flexural and longitudinal waves exists.
The transfer matrix can be written as Tcurve = eΛl , where l is the length of the curved beam.
For a rotating structure, additional terms are involved with the rotating speed Ω in the transfer
matrix as well as in the coefficient matrix of the state equations. The terms proportional to Ω are2
the Coriolis force terms, and the terms proportional to Ω are the centrifugal force terms. Coriolis
force terms are imaginary, which occur in asymmetric wave propagation.
When the periodic structure is configured in the axial direction of a one-dimensional structure withan interval L , one can express a wave vector v as follows using the Floquet’s theorem:v ( x + L) =e− µ L v( x).

(5)

Propagation constant µ is given by =µ µ R + j µ I , of which µ R denotes the amount of attenuation,and µ I means the phase difference. The sign of the propagation constant indicates the direction ofwave propagation. In the frequency range of a stopband, µ is a real number, and in a passband, µis a purely imaginary number in the periodic structures. If the transfer matrix method is employed,
one can write the following expressions:v( x +=L) Tcell v=( x) e− µ L v( x),(Tcell − λi I=)ψ i 0.

(6a,b)

Here, Tcell is the transfer matrix of one periodic cell, λi = e− µi L is an eigenvalue of Tcell , called the
Floquet exponent, and ψ i = qTi

fiT 

T

denotes an eigenvector, indicating the waveform.

2.2 Periodically Stiffened Stationary Circular Ring
For the periodically stiffened circular beam ring, the transfer matrix of a unit cell Tcell is amultiplication of Tcurve and Tstiff . In the calculation, the following ring parametric values, as
indicated in Fig.1, are used: Young's modulus, E = 70 GPa ; density, ρ = 2700 kg/m3 ; thickness,
h = 3 mm ; radius, R = 0.2 m ; length, L = 2π R / 18 ; number of separate elements, N = 18 ; width,b = 10 mm . For the stiffener, they are: density, ρ s = 7800 kg/m3 ; length, ls = 9.5 mm ; thickness,
hs = 3 mm ; width, bs , sym = 5 mm ; bs , asym = 10 mm .

The propagation constants of this test structure are calculated in the frequency domain, as depicted in
Figure 2. Near the 4200 Hz, one can find a stopband due to wave coupling. In the case of a circular ring withperiodic stiffeners of asymmetrical type, the amount of attenuation is 40% larger, and the width of the
stopband region is 36% wider than the circular ring with a symmetrical stiffener type. The ring frequency ofthe circular ring without stiffener is 4052 Hz. but the longitudinal wave appears at a lower frequency than
the ring frequency. It means that the attenuation caused by the wave coupling effect can occur at a lower
frequency.

Figure 2- Imaginary part of the propagation constant of the waves propagating on a circular beam withsymmetric/asymmetric stiffener:
symmetric stiffener;, flexural wave with asymmetric stiffener;, longitudianl wave with asymmetric stiffener;, flexural wave with, longitudinal withsymmetric stiffener.
2.3 Periodically Stiffened Rotating Circular Ring
For the simulation of a periodically-stiffened rotating circular ring, the rotational speed is set as
255 Hz. Other computation conditions are the same as the previous one. The imaginary part of the
propagation constant in a circular ring with and without rotation is shown in Figure 3. The stopbands
caused by the wave coupling appear differently in frequency regions depending on the wave
propagation direction. It happens that the wave in one direction can propagate while the other isblocked.

Figure 3-Imaginary part of the propagation constant of the waves propagating on a rotating ring withperiodic stiffeners:, longitudinal wave with rotation;, flexural wave with rotation;, longitudinal wave without rotation;, flexural wave without rotation.

3. CONCLUSION

In this work, the transfer matrix and the propagation constants of the periodically stiffened straight
beam, stationary circular ring, and rotating circular ring are estimated using the transfer matrix method
and eigenvalue problems based on Floquet’s theorem. The physical effect of a stiffener on the main
structure is considered a concentrated mass and moment of inertia in obtaining the transfer matrix.
In the structure with asymmetric stiffeners, the wave coupling caused by the interaction between
flexural and longitudinal waves occurs well. In general, the phase difference between unit cells is
multiple of π in such a case. However, the phase difference is not multiple of π in the stopband
due to wave coupling. Therefore, the wave in the stopband due to wave coupling does not cause
resonance in the entire ring structure. In the rotating ring, wave propagation characteristics aredifferent depending on their propagation direction because of the Coriolis force. Accordingly, due tothe wave coupling effect, the stopbands appear at other frequency regions for each propagation
direction. It is thought that such a physical phenomenon can find some interesting applications.
Future research tasks include the use of meta-materials, which interact more effectively withbending and longitudinal waves, and the practical application to the rotating machines by using theintentional wave coupling effect.

ACKNOWLEDGEMENTS

This work was partially supported by the NRF grant (No. 2020R1I1A206675112) in Korea.

REFERENCES
1.

2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.

Brillouin L. Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices. New York:
Dover Publications; 1953.
Mead DJ. Free wave propagation in periodically supported infinite beams. J. Sound Vib. 1970; 11: 181197.
Heckl M. Investigations on the vibrations of grillages and other simple beam structures. J. Acoust. Soc.
Am. 1964; 36: 1335.
Brunskog J. A wave approach to structural transmission loss in periodic structures: Thin beam case.
Acta Acust. Acust. 2005; 91: 91-102.
Jeon J, Ih JG, Brunskog J. Wave transmission and vibration response in periodically stiffened plates
using a free wave approach. J. Acoust. Soc. Am. 2021; 149: 3694-3702.
Thomson WT. Matrix solution for the vibration of nonuniform beams. J. Appl. Mech. 17, 337–339, 1950.
Li L. Use of Fourier series in the analysis of discontinuous periodic structures. J. Opt. Soc. Am. 1996;
A13: 1870–1876.
Mencik JM, M. N. Ichchou MN. Wave finite elements in guided elastodynamics with internal fluid. Int.
J. Solids Struc. 2007; 44: 2148-2167.
Huang SC, Soedel W. Effects of Coriolis acceleration on the free and forced in-plane vibrations ofrotating rings on elastic foundation. J. Sound Vib. 1987; 115: 253-274.
Endo M, Hatamura K, Sakata M, O. Taniguchi O. Flexural vibration of a thin rotating ring. J. Sound
Vib. 1984; 92: 261-272.
Mead DJ, Markus S. Coupled flexural–longitudinal wave motion in a periodic beam. J. Sound Vib.
1983; 90: 1-24.
Roy AK, Plunkett R. Wave attenuation in periodic structures. J. Sound Vib. 1986; 104: 395-410.
Friis L, Ohlrich M. Coupling of flexural and longitudinal wave motion in a periodic structure withasymmetrically arranged transverse beams. J. Acoust. Soc. Am. 2005; 118: 3010.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022, in Gyeongju, Korea
Extended abstractABS-1006

Experimental Study based on Vibration Analysis for Source
Localization in the Waterworks pipeline

Sun-Ho Lee 1,2; Choon-Su Park1,2; Dong-Jin Yoon1,2,*
1

Korea Research Institute of Standards and Science, 267 Gajeong-Ro, Yuseong-Gu, Daejeon, 34113,

Republic of Korea
2

University of Science and Technology, 217 Gajeong-Ro, Yuseong-Gu, Daejeon, 34113, Republic of Korea

ABSTRACT
Waterworks pipeline damage can cause not only economic loss due to leaks but also settlement. In particular,
more than 14 % of waterworks accidents in domestic are derived from third-party damage. Thus, apreemptive response for precaution of the accident is essential to solve the safety and loss problem at the site.
There are several reports of pipeline monitoring via elastic wave measurements. Particularly, fluid-structural
coupled vibration is effective for source localization in a long distance. In this study, damage detection using
accelerometers was performed in large diameter buried waterworks pipelines with a diameter and distanceof 1,200 mm and 136 m. The optimal frequency band for calculating an accurate damage detection wasdetermined by observing frequency characteristics. Our studies yielded an empirical damage detection in thepipeline that was in good agreement with the artificial damage localization.
Keywords: waterworks, pipeline, large diameter, source localization

1. INTRODUCTION
Modern cities need infrastructure for normal operation. And most infrastructure is installed on the ground,but some are buried underground due to space utilization problems. In particular, under the urban cities exist
buried pipelines for water supply. When damage occurs on the buried waterworks pipelines, various problems
will arise, such as leakage, settlement, supply interruption, incident and etc. Therefore, early detection ofpipeline damage is essential for the function of the city.
In order to solve the problem of buried pipelines, in the early days, damage assessment was dependent on
skilled workers' hearing, like the doctor's stethoscope. However, errors by subjective decisions were always
accompanied. In the 1970s, correlation analysis of acoustics was proposed to solve the error of subjective
decisions[1]. However, the detectable distance was limited to a few meters[2]. Then, it was known that awave speed analysis was needed to obtain accurate results from the correlation analysis[3]. Accordingly, theanalysis of fluid-filled pipe was performed, and as a result, waves that contribute to energy convey wereidentified with the vibration mode of fluid-filled pipe[4,5].
In addition, theoretically, wave propagation characteristic and the variables that determine the speed was
described[6]. Especially Lee and Yoon carried out the experiment for long distance abnormal detection in thein-service pipeline with a diameter of 100 mm[7].
However, most studies have been carried out on small-diameter and short-distance pipelines, the damage
monitoring for large diameter and long distance studies are insufficient.
In this study, the damage detection, which occurred by thirty-party damage, was performed on a longdistance large-diameter waterworks pipeline with a diameter and distance of 1,200 mm and 136 m. The
optimal frequency band for calculating an accurate damage detection was determined by observing
frequency characteristics. Our studies yielded an empirical damage detection in the pipeline that was*corresponding author : djyoon@kriss.re.kr

in good agreement with the artificial damage localization.

2. Fundamentals of source localization in the fluid-filled pipe
As shown in Figure1, assuming the buried pipeline was damaged, the vibration caused by damage
such as third-party damage will be propagated toward both sides. And also the vibration mode will
have appeared in various mode shapes. At the time, the structural characteristics of the fluid-filledpipe generate unique vibration characteristics.
The motion of a fluid-filled pipe can express this wave. The longitudinal waves, which propagate
through compression-expansion motion in the fluid, are constrained by the elastic pipe. In addition,these motions will cause the pipe's contraction and expansion. In other words, it propagates to thelongitudinal direction of the pipe through the inner fluid and will generate radial displacement. Thisvibration is very effective in long-range detection because of its low attenuation and ease of measuring
using uniaxial accelerometers. Additionally, this wave is called a quasi-longitudinal wave.
As well as mentioned above, the wave speed was needed to obtain accurate results of source
localization. It can be expressed by
2𝐵𝐵𝑓𝑓

𝑐𝑐 = 𝑐𝑐𝑓𝑓 �1 + 𝐸𝐸ℎ 𝑎𝑎−𝜌𝜌ℎ𝜔𝜔
𝑎𝑎2

2

−

�

1
2

(1)

where Cf is the longitudinal wave speed of the fluid in free space, Bf is the bulk modulus of fluid, a is pipe
diameter, E is young’s modulus, ρ is the density of pipe, h is the thickness of the pipe and ω is the angularfrequency.
Equation(1) represents the dispersive characteristics of the quasi-longitudinal wave, which are
related to the ω. This means that the accurate result of source localization depends on select of thefrequency band with the higher signal-to-noise ratio.

Figure 1. The schematics of third-party damage on the buried pipe

3. Expriments
In this research, the experiments were conducted at a real underground waterworks pipeline in Jeonju-Jucity, Korea. As shown in Figure 2, the 1,100 mm large diameter steel pipe is buried about 3 m underground.
The total length of the under-test pipeline is about 578 m. And the accelerometers (Sensitivity 10 V/g) wereattached to both ends of the pipe. And the artificial damage using an impact hammer (Sensitivity 0.225 mV/N)
was reproduced in the simulation of TPI.

Figure 2. Sensor position and the schematics of artificial damage

4. Result and discussion
The source localization using cross-correlation techniques can be calculated from the time delay
between two sensors. And also, it is necessary to improve the signal-to-noise ratio for an accurate
result. As a matter of fact, the main frequency band for the bandpass filter was chosen from observing
the power spectral density and coherence. Then, the time delay is obtained from the cross-correlation
function as shown in Figure 3.

Figure 3. Algorithm for the time delay estimation

Figure 4 is the raw signals obtained from each sensor with the force of 20 kN artificial damage. It
propagated 341 m and 237 m, respectively. As mentioned above, in order to improve the signal-tonoise ratio, the main frequency band was selected by observing the power spectral density andcoherence.
-5

Signal(t) of Sensor 1

10]

Signal(t) of Sensor 2

1

0
-1
-2
0

0.2

0.4

0.6

0.8

1

Acceleration [m/s1

Acceleration [m/s

-4

2

2

2

]

10
2

0

-1

-2
0

0.2

Time [sec]

0.4

0.6

0.8

1

Time [sec]Figure 4. Output signal(t) of each sensor

In order to verify the performance of the algorithm, artificial damage about 30 times was applied,
and the results are shown in Figure 5. Consequently, the result of damage detection in buried pipes ispossible with a very high probability in the optimal frequency band 30-80Hz band. However, thisresult may not be the correct answer at the current experimental stage. so it was thought that more
additional experiments to be needed in various environments.
.

Figure 5 Result of artificial damage test

5. Conclusions
This paper discussed the theoretical and experimental study on the detection of damages such as
third party damage for buried pipelines. The theoretical study describes the propagation of wave in
fluid-filled pipe and the background of source localization theory. For experimental verification, theartificial damage using an impact hammer was reproduced in the simulation of TPI.
Theoretical analysis of wave propagation speed in fluid-filled pipes is also described. A new
algorithm for the cross-correlation function is suggested and the results showed a promising damage
detection possibility in real underground water pipelines.
ACKNOWLEDGEMENTS
This research was supported by Development of safety measurement technology for Infrastructureindustry funded by Korea Research Institute of Standards and Science (KRISS – 2022 – GP20220010) and Korea Environment Industry & Technology Institute(KEITI) through Advanced Water
Management Research Project funded by Korea Ministry of Environment(MOE) (127587)

References
1. BRAINARD JR, Frank S. Leakage problems and the benefits off leak detection programs. JournalAmerican Water Works Association, 1979, 71.2: 64-65.
2. COLE, E. Shaw. Methods of leak detection: an overview. Journal‐ American Water Works Association,
1979, 71.2: 73-75.
3. FUCHS, Helmut V.; RIEHLE, Rainer. Ten years of experience with leak detection by acoustic signal
analysis. Applied acoustics, 1991, 33.1: 1-19.
4. HUNAIDI, Osama; CHU, Wing T. Acoustical characteristics of leak signals in plastic water distribution
pipes. Applied Acoustics, 1999, 58.3: 235-254.
5. FULLER, C. R.; FAHY, Frank J. Characteristics of wave propagation and energy distributions incylindrical elastic shells filled with fluid. Journal of sound and vibration, 1982, 81.4: 501-518.
.6. MUGGLETON, J. M.; BRENNAN, M. J.; LINFORD, P. W. Axisymmetric wave propagation in fluidfilled pipes: wavenumber measurements in in vacuo and buried pipes. Journal of Sound and Vibration,
2004, 270.1-2: 171-190.
7. LEE, Young-Sup; YOON, Dong-Jin; JEONG, Jung-Chae. Leak location detection of underground water
pipes using acoustic emission and acceleration signals. Journal of the Korean Society for Nondestructive
Testing, 2003, 23.3: 227-236.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0026

Control of Traveling Bending Waves in A Ring Strip Speaker
Ki-Ho LEE1; Jeong-Guon IH2
1, 2 MechanicalEngineering , Korea Advanced Institute of Science and Technology ,Korea 1 Currently at Samsung Electronics Co., Korea

ABSTRACTA strip speaker adopts a thin beam as the sound radiator controlled by the traveling-wave control method (TCM) using a small number of array actuators. The previous works investigate the strip speaker consistingof a straight beam and three actuators : the primary actuator acting as a wave source at the middle of thebeam and two control actuators near the beam edges to eliminate the reflected bending waves. In this work,
the TCM -controlled strip speaker adopts a closed -ring strip without the structural boundary to further
reduce the actuators . The bending wave continuously propagating along the ring periphery exhibits
similar characteristics to an infinite beam . The ring strip is considered the straight beam of which bothedges are connected , thus requiring a single control actuator only . The control result reveals that two
actuators are enough to suppress the structural resonances , radiating a smooth sound spectrum . The
directivity of radiated sound field depends on actuator configuration due to the axisymmetry of bending
wavefield on a ring. Based on the present work, one can design a ring strip structure to be a novel acoustic
transducer satisfying the desired performance.
Keywords: Active control, Bending wave, Inverse problem, Strip speaker, Traveling wave control, Ring

1. INTRODUCTION
The strip speaker has been suggested as a simplified panel speaker, which adopts a thin beam andsmall actuators (1, 2). It was shown that the traveling-wave control method (TCM) effectively avoids
the strip speaker's structural resonances. Three actuators are employed to this end: one for the primary
excitation and two for the wave control at the beam edges. In this work, a circular ring strip controlled by
the TCM is considered an extended concept of a straight strip speaker. While the finite straight beam has
wave reflections at the edges , waves on the circular ring strip propagate indefinitely like an infinite
beam in both directions along the ring. Such an axisymmetric wave field can radiate the sound with abroad directivity , which bears a significant potential as a speaker . As shown in Fig. 1, two different
actuator arrays can be conceived of as possible layouts of the TCM -controlled ring strip as the soundradiator. The first layout configures the primary and control actuators with a short spacing in a nearlycollocated manner. The second layout positions the two actuators oppositely to each other. The vibroacoustic characteristics of the ring strip controlled by each layout are predicted from the simulation.
(a)

(b)

Figure 1 – Configuration of actuators to control the circular ring strip:
to indicate the angular coordinate θ;1
2

kiho7.lee@gmail.com
jgih21@gmail.com, primary actuator;, center line crossing the origin, control actuator. (a) Layout 1, (b) layout 2.

2. CONTROL OF TRAVELING WAVES ON A RING STRIPThe governing equation of motion for the bending wavefield normal to the surface of a thinhomogeneous ring strip with a large radius of curvature is given by (3, 4)
D  4 v EA
2v
+ 2 v + m ' 2 = F ( 0 ) .
4
4
a 
a


(1)

Here, a denotes the ring radius from ring center to neutral surface, A the cross-sectional area definedas bh, b the ring width, h the thickness of the ring, m’ the mass density (= h), D the dynamic stiffness(=Eh3 /12(1- 2)), E the Young’s modulus,  the Poisson’s ratio, and F the excitation force on the ring
surface at θ0 position. Similar to the straight beam, the bending wavefield consists of the four wavecomponent as (5)v ( ) = C+ e− jkb + C− e jkb + C+n e − kb + C−n e kb ,(2)

where C indicates the complex coefficient of each wave component, kb the bending wavenumber. The
index n denotes the evanescent component, and the subscripts + and - represent the propagating
direction of the propagating waves. For eliminating the propagating waves by using the TCM, theinput gain of the control actuator E is estimated by applying the inverse technique to the matrix
equation as follows (1, 6):
†
H
H
ˆ
E
  =  K   Cd  = ( K K ) K Cd .
−1

(3)

Here, K denotes the transfer matrix between the input gain of the control actuator and the target wave
component, Cd the desired coefficient of the target wave component, † the pseudo-inverse operator,
and the superscript H means the Hermitian transpose of a matrix. In practical situations, the measured
transfer function is employed to secure the effect of the unknown complicated boundary condition of
the actual beam structure. The objective of Control layout 1 is to eliminate the wave propagating along
the anti-clockwise direction, which originates from the primary actuator. Control layout 2 is employed
to maintain the outgoing waves from the primary actuator only, similar to the straight beam control
case using three actuators.

3. VIBRATION FIELD AFTER CONTROLIn the simulation, a test ring strip made of steel is selected with a = 100 mm, h = 1 mm, and b =30 mm. For Control layout 1, the primary and the control actuators are located at 180° and
190°positions, respectively. Control layout 2 configures the angular position of the primary actuatorat 180° and the control actuator at 0° for symmetry. The bending wave control using the TCM isconducted with a force of 1 N driven by the primary actuator.
As a control example, Fig. 2 demonstrates the controlled vibration field of the ring strip at the fifth
resonance frequency of 823 Hz, which is represented as a function of angular coordinate θ along thering. In the case of Control layout 1, one can find that the phase is linear to the angle, which meansthat the wave propagates along the ring periphery indefinitely.
(a)

(b)

Figure 2 – Vibration field of the ring strip at the resonance frequency of 823 Hz:, Control layout 1;, Control layout 2. (a) Amplitude, (b) phase.

, uncontrolled;On the other hand, the controlled field of Control layout 2 reveals the resultant phase decreasing
linearly from the primary actuator position at θ= 180°. One can find that the controlled field by Control
layout 1 has a significant amplitude even similar to, although a little bit less than, the resonant
response of the uncontrolled field. It can be explained that the single wave propagating along the ring
strip invokes constructive interference by itself indefinitely. In the case of Control layout 2, thevibrational energy generated by the primary actuator is absorbed by the control actuator after thecontrol, in principle. The driving-point mobility of the primary actuator shown in Fig. 3 exhibits thatControl layout 1 has prominent sharp peaks at the resonances due to the self-interaction of the
propagating wave, as mentioned above.
In contrast, the structural resonances are smoothly suppressed by using Control layout 2. It should
be recalled that the fundamental objective of vibration control is to suppress the spiky resonant
response of the structure. Thus, one can conclude that Control layout 2 is better effective in converting
the ring strip into a sound radiator by moderating the spectral fluctuation caused by the structural
resonances for a broad frequency range.
(a)

(b)

Figure 3 – Driving-point mobility of the primary actuator on the test ring strip:, Control layout 1;, Control layout 2. (a) Amplitude, (b) phase.

, uncontrolled;

4. CHARACTERISTICS OF THE RADIATED SOUND
The far-field sound pressure radiated from the outer surface of the ring strip is calculated as (7)j0 e jk0 R 
n +1 bsinc ( 0.5k0 b cos  )p ( R,  ,  ) =vn ( − j )cos n ,(4)
 k0 R sin  n =0
H n ' ( k0 a sin  )where R is the distance of the sound observation point from the center of the ring, φ the elevation
angle,  0 the density of air medium, k0 the free-field acoustic wave number, and Hn ’ the first derivativeof the Hankel function of the first kind of the n th order. The radiated sound pressure level at a position
1 m apart from the ring center and normal to the primary actuator with θ= 180°, φ= 90° is depicted inFig. 4.

Figure 4 – Sound pressure level radiated from the circular ring 1-m apart from the ring center andnormal to the primary actuator:, uncontrolled;, Control layout 1;, Control layout 2.

The sound spectrum radiated from the ring strip controlled by Control layout 1 has multiple sharp
fluctuations at the resonance frequencies. On the other hand, the controlled field by Control layout 2
radiates the smooth spectrum by suppressing the resonant responses. However, one can find thespectral troughs at harmonics of 314 Hz in the case of Control layout 2, which are caused by theminimum volume velocity of the controlled vibration field (2). The directive sound radiation
characteristic contributes to the observed sound at a field point in addition to the contribution fromthe vibration field. The axisymmetry of the radiating surface affects the sound pressure amplitude
radiated from the ring strip. It is thought that the broad spectral troughs of the sound field radiatedfrom Control layout 2 can be compensated by adjusting the electric input gain of the primary actuatorto obtain the flat spectral shape.

5. CONCLUSION
In this work, the vibration control of the ring strip using only two actuators is carried out toimplement it as a loudspeaker. When the primary and control actuators are positioned in the very
vicinity, it eliminates a propagating wave component originating from the primary actuator, dependingon their arranging order. The controlled vibration field is composed of a traveling bending wave
almost indefinitely propagating on the ring. The radiated sound spectrum reveals that such aconfiguration for active control of structural wave effectively smoothes the anti-resonant responses
owing to the repeated constructive interaction of the propagating wave by itself. However, theremaining very spiky peaks at resonances can be problematic as a loudspeaker. When the primary andcontrol actuators are arranged at the far opposite positions over the ring, similar to the straight strip
speaker, the wave energy generated by the primary actuator is absorbed by the control actuator in thecontrolled field. The control result shows that the structural resonances of the test ring strip aresuppressed, thus radiating the smooth sound spectrum from the controlled vibration field. Some dull
troughs can be observed in the radiated sound spectrum, which may be compensated by simple gain
control of the input signal. Based on the findings in this study, the ring strip speaker can be designed
as a broadband acoustic radiator by using only a pair of actuators. It is thought that the sound radiatedfrom the rear surface of the ring-type beam can be isolated by sealing the ring in the sideway as acavity filled with porous materials.

ACKNOWLEDGEMENTS
This work was partially supported by the NRF grant (No. 2020R1I1A206675112) in Korea.

REFERENCES
1. Lee K, Ih J, Jung D. A strip speaker using the traveling bending wave on a beam controlled by threeactuators. J Sound Vib. 2021; 116136.
2. Lee K, Ih J. Use of the Normalized Design Parameters for Designing a Strip Speaker Operated by theTraveling-Wave Control Method. IEEE Access. 2022; 10: 22941-22951.
3. Soedel W. Vibration of Shells and Plates, New York, USA: Marcel Dekker; 1981.
4. Leissa AW. Vibration of Shells, Washington D.C., USA: NASA; 1973.
5. Cremer L, Heckl M, Unger EE. Structural-borne Sound, Berlin: Springer; 1973.
6. Bai MR, Ih J, Benesty J. Acoustic Array Systems: Theory, Implementation, and Application, Singapore:
Wiley; 2013.
7. Junger MC, Feit D. Sound, Structures, and Their Interaction, New York: MIT Press; 1985.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0159

A study on the underwater radiated noise prediction of ship
structures using the empirical radiation efficiency
Seong-Hyun LEE1; Yun-Ho SEO1; Pyung-Sik Ma1, Bong-Ki Kim1 and Sock-Kyu Lee2
1

Korea Institute of Machinery & Materials, Korea
2

LIG Nex1, Korea

ABSTRACT
In order to manage the underwater radiated noise of a ship, it is helpful to monitor vibration signals using
sensors attached to the hull. The underwater radiated noise at a speed lower than the cavitation inception
speed can be predicted using the acceleration level and acoustic radiation efficiency of a submerged plate. It
is practically impossible to measure the acceleration level for the entire plate immersed in water, and theacceleration level shows a large difference depending on the location of the sensor. The distribution ofacceleration levels of ship’s submerged plates is measured and investigated. Acceleration levels nearstiffeners show lower than those in areas without stiffeners. Although many studies have been conducted on
the radiation efficiency of small plates, that of submerged plates of large structures such as ships is not clearly
defined. The empirical radiation efficiency is investigated using measured acceleration levels and underwaterradiated noise. To improve the accuracy of underwater radiated noise prediction, the optimal location of theaccelerometer is studied using the correlation with the measured underwater radiated noise.
Keywords: Underwater radiated noise, Submerged plate, Radiation efficiency

1. INTRODUCTION
Recently, there has been an increasing demand for the reduction of underwater radiat ed noise toprotect the aquatic ecosystem. Underwater radiated noise of ships can be divided into noise generated
by machinery operating for propulsion, noise generated by propellers and flow noise. The main noise
source before cavitation is machinery noise, which is the noise generated by the machine transmitted
to the hull and radiated into the water. After cavitation occurs, the noise generated by propellers is themost important. Below CIS (Cavitation Inception Speed), noise can be predicted by the vibration level
and the radiation efficiency of the submerged plate. The radiation efficiency of the submerged plate
has been presented theoretically, numerically, and experimentally by various researchers, but most of
it is a study on a small-sized panel, so there is a limit to its application to a large structure such as aship. Measuring a ship's underwater radiated noise is not easy because it has to be measured in quiet
sea area using a measuring device such as a wireless buoy, which is time consuming and expensive.
Instead of measuring the underwater radiated noise, it can be predicted by installing acceleration
sensors on the hull. In this study, the number and location of accelerometers were investigated to
predict the noise, and the empirical radiation efficiency was presented using the measurement results.

2. VIBRATION DISTRIBUTION MEASUREMENT
In order to predict underwater radiated noise, the vibration distribution of the submerged plate ofan operating ship was measured. Since it was impossible to measure the vibration level of all plates
in contact with water, accelerometers were installed at specific locations in the four engine
compartments where propulsion equipment was installed. A total of 97 sensors were attached, andFigure 1 and Table 1 show the number and location of the accelerometers installed on each pla te.
Accelerometers were installed in various positions in the longitudinal direction, and sensors were
installed in only one position in the vertical direction. Because the tank structure was located in thelower part and the sound absorbing material was installed in the upper part, so the space that can beinstalled is limited.
PCB 352C33 was used for the accelerometer, and measurement was performed using LAN -XI from

B&K. Figure 2 shows the acceleration distribution measured in each engine room. The x -axis shows
the frame number and the y-axis shows the frequency. The vertical line marked with a red dotted line
indicates the part stiffened by the web-frame, and it can be seen that the acceleration level is measured
low near the stiffener. The distribution of the acceleration level is different for each frequency/position,and it is not easy to measure the vibration representative value of the entire plate using a limited
number of sensors.
Table 1 – Accelerometer quantity and position
Position

Quantity

No.1 ME Room[Port] 22 pts. [STBD] 10 pts.

No.2 ME Room[Port] 16 pts. [STBD] 10 pts.

No.3 ME Room[Port] 10 pts. [STBD] 9 pts.

No.4 ME Room[Port] 10 pts. [STBD] 9 pts.

FWD

Tank

Tank
#14

#13
#2

#14
#15
#16
#17
#18
#19
#20
#21
#22

#3

#17
#1

#28

#29
#30
#31

#6

WEB Frame section
Local stiffener
Accelerometer on stiffener
Accelerometer on non-stiffener

#7
#8

#10

#19

#9

#21

#22

#11

탈출구

#23

Tank

Tank

탈출구

#14

#2

#15
#16
#18

#20

#9

#32

#13

#8

#19

#10

#33
#34

#17

#5
#6

#18

#5

#24

#26
#27

#3
#4

No. 3 ME Room

#16

#4

No. 1 ME Room

#23
#25

#2

#15

#3
#22
#21
#20
#19

#11

No. 2 ME Room

#4
#5

#20
#21

#6

#22

#7

#23

#8

#24
#25

#9

#26

#10

#27

#11

No. 4 ME Room

#9
#8
#7
#6
#5
#4

#18
#17
16
#15
#14

#3

#13

#2

탈출구

#1

AFT

#28

Figure 1 – Positions of accelerometersNo.1 ME Room (STBD)

No.2 ME Room (STBD)

No.1 ME Room (PORT)

No.2 ME Room (PORT)

No.3 ME Room (STBD)

No.4 ME Room (STBD)

No.3 ME Room (PORT)

No.4 ME Room (PORT)

Figure 2 – Measured acceleration distribution: (x-axis) frame number, (y-axis) frequency.

3. RADIATION EFFICIENCY
The radiation efficiency of a finite plate considering the fluid-loading effect have been investigatedtheoretically and experimentally by various researchers [1~5], but there is a limit to its application tolarge structures such as the submerged plate of a ship. The underwater radiated noise emitted by thevibration of the ship’s plate can be explained by the following expression:
W = srad r c A < v 2 >(1)
where, W is the acoustic power (Watt), s rad is the acoustic radiation efficiency, r is the density ofseawater, c is the speed of sound in the sea, A is the area of vibrating plate, and < v 2 > is a squared
averaged velocity of the plate. Radiation efficiency can be estimated by using the vibration level ofthe submerged plate and the underwater radiated noise measured during the operation of the ship.
When measuring the acceleration of the plate described in Chapter 2, the underwater radiated noise
was measured for each speed by using a wireless buoy. However, the measured underwater radiatednoise includes the propeller noise and fluid noise, and it is impossible to extract only the noise fromradiation from the plate. The radiation efficiency was estimated using the URN measurement result,which may contain fluid noise in the vessel before CIS, and the vibration level of the submerged plate.
Figure 3 shows the estimated radiation efficiency for various speed conditions below CIS andaveraged one. When measuring underwater radiated noise, it is difficult to accurately measureradiation efficiency in a band below 63Hz because background noise is larger than noise generated byships. In the band above 1000Hz, the influence of propeller noise and flow noise is domina nt, so thecalculated radiation efficiency can include the effect of propeller noise and flow noise.
In order to compare the radiation efficiency obtained through the measurement, a precedingradiation efficiency model was calculated. It is important to se t the area as the radiation efficiency
depends on the area of the immersed plate. Assume that the area of a single plate fixed by thelongitudinal stiffener and the transverse stiffener is 2.4 (m)  0.6 (m), and the size of the plate tocalculate the radiation efficiency is (2.4(m)  2)  (0.6(m)  5) [4]. The plate thickness was assumed
to be 0.01 (m) based on the measured thickness of the ship. Figure 4 compares the radiation efficiency
model of previous researchers with the empirical radiation efficiency calculated with measurements.
It can be seen that it has a large value compared to the existing models in the frequency band of 315
Hz or higher, and additional examination is needed to determine whether the actual radiation
efficiency has a large value or is a result of the propeller and flow noise.
0

Radiation efficiency

-10

-20

A knots
B knots
C knots
D knots
E knots
F knots
G knots
H knots
I knots
Averaged

-30

-40

-50

-60
31.5

63

125

250

500

1000

2000

4000

1/3 Oct. Band Center Frequency (Hz)Figure 3 – Estimated radiation efficiency.

8000

0

Maidanik's model [1]Hattori et al.'s model [2]Cheng et al.'s model [3]Han's model [4]Kim's model [5]Estimated by measurements

Radiation efficiency

-10

-20

-30

-40

-50

-60
31.5

63

125

250

500

1000

2000

4000

8000

1/3 Oct. Band Center Frequency (Hz)

4. CONCLUSION
In order to predict the underwater radiated noise using the sensor attached to the hull, the vibratio ndistribution of the receiving plate was measured using a number of sensors and the empirical radiation
efficiency was presented using the vibration and underwater radiated noise measurement results. Since
the vibration of the submerged plate is distributed differently depending on the position a nd frequency,
sufficient consideration is required when measuring the representative vibration of the immersed plateusing a limited number of sensors. The empirical radiation efficiency proposed using the measured
results shows a larger value in the band above 315Hz compared to the previous models, so it isnecessary to study this difference further.

ACKNOWLEDGEMENTS
This research was supported by the Defense Acquisition Program Administration and the Defense
Industry Technology Center (project name: CIS monitoring technology for ship propeller using hull
attached sensors, contract number: UC200001D), and joint research with LIG Nex1 (contract number:
LIGNEX1-2020-0949 (00)) was carried out.

REFERENCES
1. G. Maidanik, Response of ribbed panels to reverberant acoustic fields, Journal of the Acoustical Societyof America, 1962; 34: 809-826
2. K. Hattori, K. Nakamachi, and M. Sanada, Prediction of underwater sound radiated from ship’s hull byusing statistical energy analysis, Inter-Noise 85, 1985; Vol. II: 645.
3. Z. Cheng, J. Fan, B. Wang, and W. Tang, Radiation efficiency of submerged rectangular plates, Applied
Acoustics, 2012; 73: 150-157.
4. H. Han and K Lee, Specification of the hull vibration to control underwater radiated noise by estimation
with modified experimental sound radiation efficiency, Journal of Mechanical Science and Technology,
2014; 28: 3425-3432.
5. H.-S. Kim, J.-S. Kim, B.-K. Kim and S.-R. Kim, A Study on the approximate formula for radiation
efficiency of a simply supported rectangular plate in water, Transactions of the Korean Society for Noiseand Vibration Engineering, 2014; 24(1): 21~27.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0239

BIW Sensitivity Analysis by Using ERP and Verification.
Junmin Park1;
1

Heejung Kwon2

Vehicle Development /Chassis/NVH Group, Hyundai Motor Europe Technical Center GmbH, Germany
2

NVH CAE Team/R&D Division, Hyundai Motor Group, Korea

ABSTRACTOne of the significant factors of interior noise in automobiles is the structure born noise caused byvibration of vehicle body transferred from noise sources (powertrain; road excitation; ancillaries). Many
vehicle body panels have drum like membrane shapes, so they can amplify the vibration, thereby increasing
interior noise levels. In this research, to solve this problem in the early stage of vehicle development, theinterior noise is predicted using CAE including ERP(Equivalent Radiated Power) analysis method, focusing
on the panels of the BIW. In addition, the validity of ERP from CAE is verified by conducting Laser
Vibrometer scanning measurements to precisely investigate the dynamic characteristics of the vehicle body(BIP). At the same time, road noise TPA with blocked force is conducted, in order to acquire accuracy andefficiency of prediction for interior noise generated by panels. Finally, effective solutions to structural noiseproblems that could occur during development are presented and verified.
Keywords: ERP (Equivalent Radiated Power), TPA (Transfer Path Analysis), BIP(Body In Prime)
NTF (Noise Transfer Function), Vibro-Acoustic Sensitivity. Laser Scanning Vibrometer.

1. INTRODUCTION
Structure-borne noise, one of main factors in cabin noise of automobile, is usually generated byvibration, transmitted from powertrain, road surface, tire and ancillaries. In particular, the panels of
the vehicle body vibrate like a drum, so it is a major cause of structure -borne noise. Recently, withgrowth of electric vehicles, the structure-borne noise from rough road, where the masking effect ofengine noise has disappeared, sounds louder and more present to customers, even though the roadnoise level has not increased. Accordingly, automobile manufacturers are trying to reduce road noise
through various numerical and experimental approaches.
The purpose of this research is to present a design guideline for reducing structural vibration andoptimizing BIW(Body In White) structures, that are excited by rough road vibration, by using ERP
analysis. In general, techniques such as mode analysis, ODS(Operational Defection Shape) analysis,
and panel contribution are used for optimal vehicle body design. Among them, the panel design
guideline by ERP is presented, which can best represent the characteristics of structural vibration.
panel

ERP  C  Vn2 S
surf

Whereas,V : normal velocity of panel
S : area of element
1
C :   Radiation Loss Factor     c
2
 : fluid density
c : speed of sound in fluid

1
2

jmpark@hyundai-europe.com
kwonheejung@hyundai.com(1)

The ERP as shown in equation (1) can be obtained numerically, e.g with a commercial program
like Nastran. However, due to the high measurement effort, cases of experimental verification of the
ERP method are rare. In this study, the ERP values from CAE simulation are verified experimentallyusing a laser scanning vibrometer. The aim of this investigation is to reduce structural vibration andsound radiation from the vehicle body by optimizing the ERP value.

2. Verification of ERP.
2.1 Transfer Path Analysis
ERP has the advantage of being able to conveniently and quickly presented the vibration
characteristics of the panel with respect to the frequency and giving a conservative estimation of theradiated noise.
To select an excitation point for the experimental and numerical verification of the ERP, a TPA(Transfer Path Analysis) was performed first in CAE and Test. With this TPA, the main contributing
paths to interior noise were identified.

Figure 1 – TPA result by CAE

Figure 2 – TPA result by Test
𝑃 = ∑𝑠𝑛𝑖=𝑠1(𝑁𝑇𝐹𝑖 × 𝐹𝑖 )(2)

Whereas,P = Interior noise
sn = number of path
NTF = Noise Transfer Function
F = input force

As shown in Figs 1 and 2, the main paths identified by the TPA
are 'Front Left Shock Absorber x direction', 'Rear right shock
Absorber z direction' and 'Front rear left subframe z-dir. According
to CAE simulation and test, these transfer paths have the greatest
influence on interior noise. In this study, the front left shock
absorber x-dir and the rear right shock absorber z-dir were chosen
as main excitation points.
Fig.3 is a graph obtained by normalizing the NTF (NoiseTransfer Function) with input force to the body. The NTF is
measured or calculated and the input force to the body is calculatedusing Transfer Path Analysis based on TPA equation (2). Throughthis, it was possible to identify the main paths that have the greatest

influence on the interior noise. Calculating the ERP values by excitation of these paths is the most
reasonable to relate the ERP value and the noise.
We used Blocked Force In-Situ Method and classical TPA method to verify the spindle load used
for TPA calculation and road noise simulation. Blocked force is very useful for CAE simulation
because it is independent of the FRF changes of the receiver.
2.2 Laser Scanning Vibrometer Test.
To verify the calculated ERP by CAE, the surface velocity of the panel vibration is measured witha laser scanning vibrometer. The body is excited with a shaker at the input points identified in section
2.1. In order to obtain a good spatial discretization of the ERP values in the test, the vibration signal
for about over 2,800 nodes is measured. The validity of ERP is verified by comparing and correlating
the measured and calculated values. In addition, a parameter study is done by adding mass at certain
points of the body and by measuring a body with and without damping sheets.

Figure 4 – Laser Scanning Vibrometer with BIP and Shaker Positions
As shown in Fig. 4, to measure the vibration of the panel, the BIP (BIW + Windshield) is mounte don air springs, and the Front Left Shock Absorber (Fig. 5) and the Rear Right Shock Absorber, which
are the main paths of the interior noise, are excited. The vibration of each BIP panel is measured witha 3D laser scanning vibrometer positioned by two robot arms. In order to measure ERP values in roadnoise frequency are up to 500Hz, it is necessary to set a small distance between scanning points for
sufficient spatial discretization. As mentioned above, about 2,800 ~3,000 points on the entire BIP
panel were automatically measured over night as the laser scanning vibrometer with positioned
automatically by two robot arms.

Figure 5 – Scanning points of BIP

Figure 6 – Normal vector velocity of vibration.

Depending on the main excitation path and the direction of the excitation, the sensitivity of thereacting panel is also different. In the case of the Front Right Shock Absorber x -dir, mainly dash board
and windshield reacted sensitively. In the case of the Rear Left Shock Absorber z-dir, the sensitivityof the floor and roof was high. Table 1 summarizes the configuration of each experiment.

First, the BIP composed only of steel andglass (without damping sheets) was measured
and compared with the CAE analysis results.
The purpose of this is to exclude the non-linear
component of the damping sheet. Then, the BIP
with the damping sheet attached was measured,and finally, a mass study was carried out by
installing a mass damper on the most sensitive
part. CAE simulations were also calculated and verified under the same boundary conditions.
One of the problems when using a 3D laser
scanner is, that each X, Y, and Z velocity value
is measured in a global coordinate system, but
the velocity value normal to the surface is not
automatically accessible. To solve this issue,
considering the fact that the main vibration
direction of the panel is out of plane, thedifference between the calculated value
assuming the Vector Norm value (equation (3))
as Normal Velocity and the Panel ERP value
calculated in Nastran is exactly same. And eventhough the big difference between number ofnodes in the FE model and scanning points, thecalculated ERP is within 5%, as shown in figure
7. Therefore, the Vector Norm value can be used to verify ERP consistency.
𝑉𝑛 = √𝑉𝑥2 + 𝑉𝑦2 + 𝑉𝑧2
2.3 Verification.
2.3.1 Mode shapes of the main peaks.
Table 2 – Mode shapes of main peaks of test and CAE(3)

In table 2, the ODS (Operating Deformation Shape) of the main peaks from figure 7 are compared,
showing only minor differences in the natural frequencies and deformation shapes between CAE
simulation and laser scanning vibrometer test.
2.3.2 Frequency band for verification
The measurement of the laser 3D scanner
vibrometer in this experiment guarantees over 80%
coherence reliability from 40Hz to 340Hz as shownin figure 8.. Since it mainly covers the frequency
range of structural noise sufficiently, this frequency
band was set as the ERP verification area.

2.3.3 ERP comparison between test and CAE

Figure 9 – ERP comparison graph.
When shaking at the front shock absorber x-dir, which is the first excitation point, mainly the dashboard and windshield react the most sensitively, so the ERP test and simulati on of these panels werecompared. Exciting the body at the second point, Rear Shock Absorber z-dir, the main response is at
the floor and roof, so these panels were compared. As shown in Fig. 9, the test and simulation was
well matching in the main peaks of the booming (70 ~ 180Hz), tire cavity (180 ~ 220Hz), and rumble
band (220 ~ 340Hz), which greatly affect cabin noise.
However, the predefined road noise frequency band (see table 3) of the standardized road noise
analysis too wide to be appropriated for ERP analysis. Therefore, a 1/3 Octave band with A weighting
is used for the analsys. This allows a better link of the noise analysis to human perception.
Table 3 – ERP frequency band
Frequency band
Road noise band
1/3 octave band

Figure 10 – 1/3 Octave ERP comparison graph.
The test and the simulation results were compared in 1/3 octave as shown in Fig.10. Due to only
minor differences between Test and CAE reuslts, the ERP values in 1/3 octaves are used as a mainreference for body design in the early development of the vehicle.

3. Parameter study of ERP.
3.1 Damping sheet effects.
Two BIPs fromt the same lot with (figure 11) and without damping sheet were prepared and thevelocity of the vibration of the panels were measured under the same conditions.
On the dash board, the panel vibration was reduced due to thedamping sheet effect in a wide frequency band from 100Hz to
350Hz(see figure 12, first line) On the center floor and rear floor,
the effect of the damping sheet mainly occurs between 150Hz and
350Hz, as shown in the lower section of figure 12.

Figure 11 – Damping sheet on the BIP

Figure 12 – Test results of damping sheet influence

Figure 13 – CAE results of damping sheet influence

The CAE simulation also confirms, that the panel vibrations were reduced due to the damping
sheet effect from 100Hz, as shown in figure 13.
3.2 Mass damper effects.
Through tests and simulations, the ODS of the frequency with the largest ERP value was identified
and a mass damper 1.5kg was installed on the part with the largest amplitude.

Figure 14 – Mass damper(1.5kg) on dash board and floor.

Figure 15 – Test results of mass(1.5kg) influence

Figure 16 – CAE results of mass(1.5kg) influence
As seen in figure 15, the ERP value of the dash board was improved at 160 Hz, 202Hz and 275Hz,
which also affects the interior noise. On the other side, there is no perticular influence on the ERP
value of the center floor. The CAE results shows similar results as seen below in figure 16.

4. Road noise
4.1 Road noise prediction
After verifying the ERP method on BIW/BIP level, road noise was simulated using a full vehicle
FE model with an added mass damper (section 3.2), utilizing the ERP analysis. The interior noise intire cavity and booming frequency area were significantly reduced by 5 dB(A) ~ 6 dB(A) comparedto the baseline model as shown in figure 17.

Figure 17 – Road noise prediction with mass damper by CAE

5. CONCLUSIONSERP analysis is used to analyze structural born noise and predict the sound radiation from panelsof the vehicle body. First, a TPA analysis was performed to identify the main contributing paths tointerior noise. Based on the main transfer paths, the ERP values are calculated to approximate theinterior noise.
Using a laser 3D scanning vibrometer, the velocity of panel vibration was measured. A shaker wasused to excite the body at front right shock absorber x-dir and rear left shock absorber z-dir, which
are the main paths identified by the TPA. The ERP value from the measurements is calculated using
the Vector Norm and is used to verify the ERP results from CAE simulation. In addition, the sensitivityof the panel was analyzed more precisely by utilizing the 1/3 octave band rather than the wide roadnoise frequency band when evaluating.
In order to understand the effect of applying damping sheets and a mass damper to the vehicle
body, the ERP was analyzed experimentally and by CAE simulation. The results of this comparison
show, that ERP results obtained from simulation are reliable and therefore suitable to optimize thebody structure in an early development stage. Depend on this verification, also road noise can bepredicted on vehicle level by CAE simulation.
Finally, based on this research project, road noise and body related vibration issues can be solved
by CAE in early vehicle design stage by setting and managing ERP based target values.

REFERENCES
1. Marinus Luegmair and Hannes Münch , Advanced Equivalent Radiated Power (ERP) Calculation for
Early Vibroacoustic Product Optimization. The 22nd International Congress on Sound and Vibration.
2. Seungchan Choi, Hyukjin Oho, Kwanghyun Hwang, Kyoung-Jin Chang. Noise Reduction Design
Methodology for the Passenger Vehicle Considering Acoustic Characteristics of the Interior Cavity.
Transactions of KSAE, Vol. 29, No. 3, pp.263-271 (March, 2021).
3. Matthias Klaerner, Steffen Marburg, Lothar Kroll3 FE based measures for structure borne sound radiation.
Inter Noise 2014. p. 216
4. Reger, J., Schneider, T., Ehlert, C. Optimization of Car Body Parts Regarding ERP Using a Genetic Algorithm and Morphing, Proceedings of the 3rd European Conference on Computational Mechanics,
Lisbon,borne, Portugal, 5 – 8 June, (2006). sound radiation. Inter Noise 2014. p. 216
5. Dominik Grauvogl , Jeffrey Krabinski, NVH Comparison of a Novel Wound Field Synchronous
Machine with Magnetic Asymmetry with a PSM for a HV Electric Drive, JSAE 2021 Grauvogl.
6. Aya Karim & Jesper Brindhag, Structure–Acoustic Interaction between Vehicle Floor Panels and Carpets
ISRN LUTVDG/TVSM--21/5251--SE (1-69) | ISSN 0281-6679

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0083

A study on sound absorption of double-elastic plates in impedance
tubes with consideration of near-field effects
Hyun-Sil KIM1; Pyung-Sik MA; Seong-Hyun LEE
1

Korea Institute of Machinery and Materials, Korea

ABSTRACTIn this study, sound absorption of the double thin plates in impedance tubes are considered, where the plates
are elastic and have micro perforations. Cross-section of the impedance tubes are circular or rectangular. The
sound absorption coefficient is calculated after solving structure-acoustic coupled equations. When the plates(with or without micro-perforations) are flexible, there occur interactions between the plates and the acoustic
fields. Accordingly, the waves inside the cavities have higher order terms in addition to the plane wave,
where these higher order terms are near-field waves and exponentially decay near the plates. The soundabsorption coefficient of the micro-perforated plates (MPPs) are also computed by the finite element method(FEM). The comparison of the analytical solutions and the results by the FEM shows good agreement, which
confirms the accuracy of the proposed method. It is shown that as more near-field waves or higher orderterms are included, the results show convergence and the peak frequency moves toward a lower frequency.
Keywords: Sound absorption, Impedance tube, double plates

1. INTRODUCTION
Micro-perforated panels or plates (MPPs) are considered as next-generation sound absorbers, as
they are free of fibrous or porous materials, resistant to humidity, and lightweight. MPPs are usually
made of thin metal sheets or plastics, with the perforation diameter typically being in the submillimeter range. MPPs are widely used as sound-absorbing materials in many engineering
applications such as silencers [1,2], small-scale enclosures [3], noise barriers [4], etc.
When one measures sound absorption coefficient using an impedance tube, there exists an air cavity
between the specimen and the backing wall. When the MPPs are rigid, all of the waves inside the tube areplane waves. However, when the MPPs (with or without micro-perforations) are flexible, there occur
interactions between the MPPs and the acoustic fields. Accordingly, the wave inside the cavity has higher
order terms in addition to the plane wave, and the reflected wave includes a near-field wave or an
exponentially decaying wave near the specimen.
Toyoda et al. [5] studied the sound absorption of an MPP with a rigid backing wall installed in a circularduct. They solved structure-acoustic coupled equations that accounted for near-field waves. When a finite
elastic plate with a backing air cavity is surrounded by an infinite rigid baffle and exposed to a plane incident
wave in a semi-infinite acoustic medium, one must solve a structure-acoustic coupled problem that takes into
account semi-infinite acoustic media, a finite elastic plate and an air cavity surrounded by a rigid baffle. This
involves complicated mathematical formulations, and analytic solutions are available only for simple
geometries [6,7].
Kim et al. [8] investigated the near-field effects on the sound transmission and absorption of elastic MPPs
in impedance tubes. For low-frequency sound absorption of a circular plate, they studied how the peakfrequency is affected by the MPP parameters such as the cavity depth, radius and surface mass density of theplate.
In this study, the effects of the near-field waves on the sound absorption of flexible double MPPs in
impedance tubes are investigated considering circular and rectangular cross-sections. The interactions
between the flexible MPPs and acoustic fields including near-field waves are formulated. Predictions of thehskim@kimm.re.kr

1

sound absorption coefficient from the proposed method are compared with the results by the finite element
method (FEM) to confirm the accuracy of the proposed method.

2. FORMULATION

Fig. 1. Sound absorption of double-elastic MPPs in an impedance tube. For a circular tube, the radius is a ,and for a rectangular tube, the cross-section is L1  L2 .
2.1 Circular impedance tubei (  t − kz )In a circular impedance tube with radius a , a plane wave Ae
is incident onto double-elastic

MPPs, as shown in Fig. 1, where i = −1 ,  is the angular frequency, k =  / c , and c is the speed ofsound in air. A plane incident wave induces only axisymmetric vibration modes in circular MPPs.
it

Accordingly, the reflected wave with the omission of the common term e

is given bypr ( r , z ) = Ae +ikz

 R e J ( r / a ) ,ik z z

s

0

s

(1)

s =0

where (r , z ) is the circular cylindrical coordinate, Rs is the coefficient of the reflected wave andJ 0 ( s r / a ) is the Bessel function of order zero. It is assumed that the impedance tube is rigid. Accordingly,
the parameter  s is determined from the boundary condition, pr / r = 0 or dJ 0 ( s r / a) / dr = 0 at

r = a , where  s is a root of the equation J1 ( s ) = 0 and  0 = 0 , 1 = 3.8317 ,  2 = 7.0156 , etc.
The wavenumber k z in thez direction satisfies the following relationship:( r / a )2 + kz2 = k 2 .
s(2)

In Eq. (2), k z becomes an imaginary number when k   s / ak z = −ik z = −i ( s / a ) 2 − k 2 .

(3)Inside the first cavity (0  z  q1 ) , the wave is generated aspm = [ Bs e−ik z + Cs eik z ]J 0 (  s r / a ) .
z

z

(4)

s =0

The wave in the second cavity (q1  z  (q1 + q2 )) is generated aspt =  [ Ds e−ik z + Es eik z ]J 0 (  s r / a ) .
z

s =0

The boundary condition at the jth MPP ( j = 1, 2 ) is

z

(5)

−

1 p

= vj
at z = z j ,(6)i zin which z1 = 0 and z2 = q . At z = q1 + q2 , velocity at the rigid block becomes zero. In Eq. (6),  isthe density of air and v j is the average velocity of the jth MPP, which is related to the plate velocity v j , p
and the fluid velocity v j , f through the hole according to the following equation [9]:v j = v j , p (1 −  j ) + v j , f  j .

(7)

In Eq. (7),  j is the perforation ratio defined as  j =  d j / 4l j , where d j represents the diameter2

2of the hole and l j is the distance between the holes of the jth MPP. The impedance of the jth MPP isrelated to the pressure difference according to [9]Z j , resist (v j , f − v j , p ) + Z j , react v j , f = p j .

(8)

The impedance of the hole, Z j = Z j , resist + Z j , react , is given in Ref. [10]. Eliminating v j , f from Eqs.
(7) and (8) yields

vj =  jv j, p + j p j

,

Zj

(9)

in which  j = 1 −  j ( Z j , react / Z j ) . The pressure difference p j is given byp1 = 2 A +  ( Rs − Bs − Cs ) J 0 (  s r / a )at z = 0,

(10)

s =0



p2 = [ Bs e−ik q + Cs eik q − Ds e −ik q − Es eik q ] J 0 (  s r / a )
z 1

z 1

z 1

z 1

at z = q1 .

(11)

s =0

The velocity continuity of the MPP in Eq. (6) yields

Rs = − Bs + Cs
− Bs e− ik z q1+ Cs eik z q1= − Ds e

at z = 0,− ik z q1+ Es eik z q1(12)
at z = q1 .

(13)

At z = q1 + q2 , zero velocity at the rigid block requires− Ds e − ik ( q + q ) + Es eik ( q + q ) = 0.
z

1

2

z

1

2

(14)

If the displacement of the jth MPP is  j =  j e ( v j , p = i j ), the governing equation becomesit

D j 4 j − M j 2 j = p j ,(15)

where D j = E j h j / 12(1 − j ) , M j =  pj h j , and E j ,  j ,  pj are the Young’s modulus, Poisson’s ratio
32

and density of the jth MPP, respectively. In Eq. (15), structural damping is considered by introducing acomplex Young’s modulus E (1 + i ) , where  is the loss factor.
The axisymmetric displacement of a circular plate is expressed as j (r ) =  a j ,mm (r ) .

(16)

m =1

In this study, plates are all assumed to have clamped boundary conditions. The mode m ( r ) satisfying theclamped boundary condition m = 0 at r = a is given in Ref. [11].
Using Eq. (16), Eq. (15) can be rewritten as

M j  a j , m ( 2j ,m −  2 )m (r ) = p j .

(17)

m =1

After taking some mathematical manipulation (see the details in Ref. [8]), one has matrix equations for {Rs }and {Ds }C1R + C2 D = F,

(18)

C3 R + C4 D = 0.

(19)

In Eqs. (18) and (19), C j ( j = 1, 2, 3, 4) are the coefficient matrices, R and D are the unknown vectors

RT = {R0Rs

RN } and DT = {D0s

DN } , where N s is the summation limit for s , and F

Dssdenotes a column matrix including A . The size of the matrix C j is NT  NT , in which NT = N s + 1 .
For numerical examples, we consider sound absorption of MPPs made of steel or aluminum plates. The
material properties of aluminum are  p = 2700 kg/m3, E = 7.2 1010 Pa, and  = 0.34 while those ofsteel are  p = 7800 kg/m3, E = 2.11011 Pa, and  = 0.31. The loss factors for steel and aluminum areassumed as  = 0.01. The properties of air are  = 1.2 kg/m3 and c = 343 m/s, and the viscosity is 0 = 17.9  10 −6 kg/ms. The diameter of the cylindrical impedance tube is 100 mm, and the cut-off frequency
above which the wave decays exponentially is given by 2 f / c = 1.8412 / a , or f = 2010 Hz. The
frequency of interest in this study is below 1600 Hz, for which k z in Eq. (2) is real when s = 0 andbecomes imaginary when s  1 .
After solving the matrix equations (18) and (19), the sound absorption coefficient  is defined as

 = 1− | 1 + R0 / A |2 .

(20)

2.2 Numerical examples for circular MPPs
Figure 2(a) shows the absorption coefficient  of double aluminum MPPs, where h1 = h2 = 0.3 mm,

d1 = d 2 = 0.3 mm, q1 = q2 = 30 mm and the perforation ratios are  1 =  2 = 0.05%. In Fig. 2, thepredictions by Eqs. (18) and (19) when N s = 0
COMSOL [12], showing good agreement.

Absorption Coefficient

1.0

to 4 are compared with the result by the FEM using

FEM
A

0.8

0.6

0.4

0.2

0.0
200

400

600

800

Hz

(a)

1000

1200

1400

1600

12.1 Hz

Absorption Coefficient

1.0

12.1 HzNr = 00.8

Nr = 10.6

0.4

Nr = 2,3,4
0.2

FEM
0.0
1180

1200

1220

1240

1260

Hz

(b)
Fig. 2. Sound absorption of circular aluminum MPPs when h1 = d1 = 0.3 mm, h2 = d 2 = 0.3 mm,

 1 =  2 = 0.05 % and q1 = q2 = 30 mm: (a) at a frequency range of (0, 1600 Hz), and (b) enlargement of
the circle ‘A’ in (a).

Fig. 2(b) is an enlargement of the circle ‘A’ in Fig. 2(a). It is observed that as N s increases from zero
to 4,  shows convergence and the predictions when N s = 2, 3 and 4 are indistinguishable. Recall thatN s = 0 represents the plane wave solution, while N s  1 denotes the solutions including the higher orderterms. As discussed in the earlier work [8], the peaks move toward a lower frequency as higher order terms
are included.
2.3 Rectangular impedance tube
When the cross-section of an impedance tube is rectangular, L1  L2 , the reflected wave with theomission of the common term ei tis expressed as

pr = Aeikz +   Rrs eik z cos( rx / L1 ) cos( sy / L2 ).
z

r =0

(21)

s =0

The wavenumber k z in the z direction satisfies the following relationship: 2 (r / L1 ) 2 + ( s / L2 ) 2  + k z2 = k 2 .

(22)

The waves in the first cavity (0  z  q1 ) and in the second cavity (q1  z  (q1 + q2 )) are defined in thesame manner similar to Eq. (21) using the cosine functions as the mode shape.
When k  22

(r / L ) + (s / L )  in Eq. (22), k becomes an imaginary number
2

1

2

2

z





k z = −ik z = −i  2 (r / L1 )2 + (s / L2 ) 2 − k 2 .

(23)

Taking similar manipulation to the equations for the circular impedance tube in Eqs. (1) to (17) results in thesimilar matrix equations as Eqs. (18) and (19). The cross-section of the impedance tube is 0.26 m x 0.21 m,of which the corresponding cut-off frequency is 660 Hz. The frequency of interest for the rectangular cross-section in this study is below 600 Hz, for which k z in Eq. (22) is real only when r = s = 0 , andbecomes an imaginary number for other values of r and s .

kz

2.4 Numerical examples for rectangular MPPs
Figure 3(a) shows the absorption coefficient of double steel MPPs, where h1 = d1 = 0.2 mm, h2 = d 2 =0.3 mm, q1 = q2 = 30 mm and the perforation ratios are  1 =  2 = 0.08%. Fig. 3(b) shows anenlargement of circle ‘A’ in Fig. 3(a) around 348 Hz. The absorption coefficient  shows convergence as
N r increases from zero to 6. Note that N r = N s in all numerical examples for rectangular MPPs, and theresults with N r − 1 are identical to those with N r , because matrix equations vanish when ( m + r ) is aneven number. Accordingly, N r increases from 0, 2, 4, 6, etc.
1.0

Nr=6
FEM

Absorption Coefficient

0.8

A
0.6

0.4

0.2

0.0
100

200

300

400

500

600

Hz

(a)
1.0

3.9 HzFEM

Nr = 0Absorption Coefficient

0.8

Nr = 2
Nr = 40.6

Nr = 6,8

0.4

0.2

0.0
320

330

340

350

360

370

Hz

(b)
Fig. 3. Sound absorption of steel MPPs when h1 = d1 = 0.2 mm, h2 = d 2 = 0.3 mm,  1 =  2 = 0.08 % and

q1 = q2 = 30 mm: (a) at a frequency range of (0, 600 Hz), and (b) enlargement of the circle ‘A’ in (a).

Figure 4(a) shows the absorption coefficient of double aluminum MPPs, where N r increases from 0 to
6. The MPP parameters are h1 = d1 = 0.7 mm, h2 = d 2 = 1.0 mm, q1 = q2 = 30 mm and the perforation
ratios are  1 =  2 = 0.05%. Fig. 4(b) shows an enlargement of circle ‘A’ in Fig. 4(a) around 510 Hz and
535 Hz. The absorption coefficient  shows convergence as N r increases from zero to 6. Figs. 2, 3 and4 show good agreement between the predictions and the results by the FEM.

1.0

Absorption Coefficient

0.8

A
0.6

0.4

0.2

0.0
100

200

300

400

500

600

Hz

(a)

8.1 Hz6.0 Hz1.0

Nr = 4,6

Nr = 0Absorption Coefficient

0.8

Nr = 2
0.6

0.4

0.2

FEM
0.0
500

510

520

530

540

550

Hz

(b)
Fig. 4. Sound absorption of rectangular aluminum MPPs when h1 = d1 = 0.7 mm, h2 = d 2 = 1.0 mm,

 1 =  2 = 0.05 % and q1 = q2 = 30 mm: (a) at a frequency range of (0, 600 Hz), and (b) enlargement of
the circle ‘A’ in (a).

3. CONCLUSIONSThe sound absorption coefficient of double elastic MPPs in impedance tubes was computed by solving
structure-acoustic coupled equations, and compared with the results by the FEM, which confirmed theaccuracy of the proposed method. As additional higher order terms are included, the results showed
convergence. It was shown that as more near-field waves are included, the peak frequency moves toward alower frequency.

ACKNOWLEDGEMENTS
The present study was supported by the project “Artificial intelligence-based diagnosis and damage
management technology” at the Korea Institute of Machinery and Materials.

REFERENCES
1. Bravo T, Maury C, Pinhède C. Optimisation of micro-perforated cylindrical silencers in linear andnonlinear regimes. J. Sound Vib. 2016;363:359-379.
2. Yu X, Cheng L, You X. Hybrid silencers with micro-perforated panels and internal partitions. J. Acoust.
Soc. Am. 2015;137(2);951-962.
3. Yang C, Cheng L. Sound absorption of microperforated panels inside compact acoustic enclosures. J.
Sound Vib. 2016;360:140-155.
4. Asdrubali F, Pispola G. Properties of transparent sound-absorbing panels for use in noise barriers. J. Acoust.
Soc. Am. 2007;121(1):214-221.
5. Toyoda M, Mu RL, Takahashi D. Relationship between Helmholtz-resonance absorption and panel-type
absorption in finite flexible microperforated-panel absorbers. Appl. Acoust. 2010;71:315-320.
6. Mechel FP. Panel absorber, J. Sound Vib. 2001;248:43-70.
7. Kim B-K, Kang H-J, Kim J-S, Kim H-S, Kim S-R. Tunneling effect in sound transmission loss
determination: Theoretical approach. J. Acoust. Soc. Am. 2004;115(5):2100-2109.
8. Kim H-S, Ma P-S, Kim B-K, Kim S-R, Lee S-Y. Near-field effects on the sound transmission andabsorption of elastic micro-perforated plates in impedance tubes. J. Sound Vib. 2012;499:116001.
9. Takahashi D, Tanaka M. Flexural vibration of perforated plates and porous elastic materials under acoustic
loading. J. Acoust. Soc. Am. 2002;112(4):1456-1464.
10. Maa D-Y. Microperforated-panel wideband absorbers. Noise Control Eng. J. 1987;29(3):77-84.
11. Leissa AW. Vibration of Plates, Acoustical Society of America; 1993 (Chapter 4).
12. COMSOL Multiphysics, COMSOL Multiphysics Modeling Guide: Version 4.4, COMSOL AB,
Stockholm, 2013.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0336

Airborne Interior Noise Contribution Analysisof an Electric Vehicle Based on SEA
Ki-Sang CHAE1; Eun-Soo JO2; Dong-Chul PARK3; Hyun Woo KIM3;
Woo-Keun SONG4; Andreas SCHUHMACHER5; Moon-Ju HWANG5
123
45

Hyundai Motor Company, South Korea

Brüel & Kjær Sound and Vibration Measurement A/S, Denmark

ABSTRACTSEA (Statistical Energy Analysis) is widely applied to analyze the airborne contribution of vehicle interior
noise at high frequencies such as above 400 Hz, usually in 1/3 octave band without phase information. For
electric vehicles, the motor-related or power electric (PE)-related airborne contribution is very serious to thepassengers, and it is very important to identify the mainly contributing paths and regions. This study ispresenting an SEA model build-up process for an electric vehicle from the acoustic source estimation, theintermediate exterior acoustic transferring, the sound transmission through the body panels, and the interior
acoustic transferring. The acoustic sources of the vehicle are experimentally estimated by applying SPC(Source-Path-Contribution) for several operating conditions with an electric vehicle. The exterior
propagations from the acoustic sources to the vehicle body panels are measured in terms of acoustic transfer
functions. The body-side SEA panels and interior cavities are modeled by using the geometry of the vehicle.
Some of the SEA panels’ sound transmission loss (STL) are set by using SEA simulations and some is set as
the one of the database measured in the previous vehicles. Finally, SEA-simulated source contributions in
the interior noise are compared to the SPC-estimated contributions.
Keywords: SEA, Electric Vehicle, Source, Path, Contribution

1. INTRODUCTION
The market share of EVs is continuously increasing. Automobile manufacturers need to prepare
prediction tools to secure the indoor noise quality of EVs. In this study, the vehicle noise source
quantification method was applied to an electric vehicle, and the SEA-based interior noise analysis (1~ 4) was performed by using the quantified noise sources. There is an another paper by the authors
about this study (5), which is presenting the study in the viewpoint of binaural auralization by
combining the SPC (6, 7) and the SEA (8, 9). This paper is presenting the study in the viewpoint of
SEA model build-up and the airborne contribution analysis.

2. SEA MODEL BUILD-UP PROCESSThe outline of SEA modeling using the quantified noise sources in this study is shown in Figure 1.
The source-path-control (SPC) technique was applied to experimentally quantify the vehicle noise
sources in driving conditions as the volume velocities (Q ith) as shown in Figure 2 (6, 7). The exteriorside transfer functions from the noise sources to the vehicle SEA panels, where the SEA panel division
is as shown in Figure 3, were measured by using a monopole source (<|p ext, panel_jth |>/ |Q ith|) (8) asshown in Figure 4 and Figure 5, and used to calculate the sound pressure loading on the panels (<|pext,22
panel_jth | > operating ) by using the noise sources’ volume velocity in operating conditions (|Q ith | operating ).
Hear, < > denotes the average over the panel area. The sound transmission coefficients or sound1
2
4kschae@hyundai.com
goodwater@hyundai.com
Wookeun.Song@bksv.com

transmission loss (STL) of the vehicle SEA panels were prepared by applying SEA simulations, or byacoustic measurements, or selected from the STL database of previous vehicles. The interior-side
transfer functions (<pear2>/Π cavity near panel_jth) from the interior SEA cavity near the panel to the ear
position were calculated by using the SEA model for the vehicle interior space. Using the STL andthe the interior-side transfer function, the pressure ratios or pressure transmissibility (|<p ear 2>/<pext,2
panel_jth >|SEA ) of the ear position sound pressure over the exterior sound pressure loading s were
calculated for all the SEA panels as shown in equation (1) and Figure 6 and Figure 7. The ear position
sound pressure was calculated by multiplying <p ext, panel_jth2> estimated and |<pear2>/<p ext, panel_jth2>|SEA for
all the SEA panels and summing the multiplications (9).
The applied SEA model build-up process is different with traditional automotive SEA applications(2 ~ 4), at which all the input data are put to SEA simulation tools. In this study, the pressure ratiosfrom the panel exterior sound pressure loading to the ear position sound pressure are calculated fromthe body-side SEA models and summarized in a text-based database (Excel file), and utilized to
calculate the ear position binaural sound signal using Matlab (5) with the estimated source volume
velocities in time signal and their transferring functions to SEA panels. For the purpose of generating
binaural sound signal at the ear position, the time-based data processing had been studied and applied
to mix the time-based source volume velocities and the 1/3 octave-based transfer function from thesource to the receiver (5).

Figure 1 – SEA modeling process

10 dB

Figure 2 – SPC-estimated airborne sources

Figure 3 – SEA panel division (right-side view)

Microphones at the source positions(measuring |p|2)

A monopole source over a SEA panel(measuring |Q|2 during area scanning)

Figure 4 – Measuring the exterior transfer functions
Time (secs)

Time (secs)(b)
Exerior Transfer Funtion (dB)

Exerior Transfer Funtion (dB)

(a)

20 dB

20 dB

Figure 5 – Measured exterior transfer functions. (a) From the sources to the windshield area, (b)from the sources to the front-foot floor area2pear

① Panel SPL Loading② SEA panelpTransmitted power Π③ Power-transmitted (incoming) SEA cavity

2exterior, panel jth

Figure 6 – The concept of SEA-based pressure ratio
2
pear
2
pext,
panel jthApanel jth
40c0 panel jth 2preceiving
cavity power -incoming cavity

(1)

SEA

1.E+00
1.E-01
1.E-02
1.E-03Pressure Ratio

1.E-04
1.E-05

1.E-06
1.E-07
1.E-08
1.E-0910 dB

1.E-101.E-11
1.E-12
1.E-13
1.E-14
1.E-15

400

500

630

800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000

Frequency (Hz)Figure 7 – SEA-based pressure ratio of the vehicle (y axis in log)

3. AIRBORNE INTERIOR NOISE CONTRIBUTION ANALYSISThe calculated SEA-based airborne interior noise spectrum at the driver’s ear position is shown in
Figure 8, and this spectrum is resembling the SPC-estimated spectrum at the driver’s right ear. At thedriving condition of 40 mph constant speed, the airborne contribution by SEA or SPC is similar to themeasured between 500 Hz and 2 kHz. Below 500 Hz, the structure-borne contribution was identified
as the main contribution by the SPC analysis, and above 2 kHz, the difference between the measured
and the airborne contribution may be caused by wind noise. For the WOT condition, the similarities
between the measured, the SEA-calculated airborne, the SPC-estimated airborne seems higher, andthe gap between the measured and the airborne seems to be smaller compared to the constant -speed
driving condition. Figure 9 shows the comparison of SPC-estimated spectrum and SEA-based
spectrum along with the time for the case of WOT condition. In the case of WOT conditions where
the order components of the motor noise are important, the main order and the higher orders are wellvisible on the color map in the SPC-estimated spectrogram, while the main order is well visible in the
SEA-based spectrogram, but the higher orders are not clearly visible in the SEA-based spectrogram.
Figure 10 shows the calculated SEA panel contribution to the interior noise at the driver’s ear position,and the mainly contributing areas could be identified. Compared to the 40 mph case, it can be seenthat the high frequency component is larger under the WOT condition.
SEA
Measured-Driver Left
Measured-Driver Right
SPC-Airborne-Driver Left
SPC-Airborne-Driver Right

10 dB
(a)

SEA
Measured-Driver Left
Measured-Driver Right
SPC-Airborne-Driver Left
SPC-Airborne-Driver Right

10 dB
(b)

Figure 8 – Calculated airborne interior noise at the driver’s ear position. (a) Running at 40 mph
constant speed on a rough road, (b) running under WOT condition on a smooth road (RPM
averaged)Time (secs)

Time (secs)

Frequency (kHz)Frequency (kHz)Frequency (kHz)Frequency (kHz)Frequency (kHz)

(b)

(a)

10 dB
Time (secs)

10 dB

Time (secs)

Figure 9 – Calculated airborne interior noise contribution under WOT running condition on asmooth road. (a) SPC-estimated, (b) SEA-based.

20 dBSound Pressure Level (dB)

Sound Pressure Level (dB)

20 dB

(a)

(b)

Figure 10 – Calculated airborne interior noise SEA-panel contribution at the driver’s ear position.
(a) Running at 40 mph constant speed on a rough road, (b) running under WOT condition on asmooth road (RPM averaged)

4. CONCLUSIONSIn this study, an SEA-based simulation technique linking the time signal of the vehicle noise source
by SPC and the body-side SEA model was studied and applied to an electric vehicle. The SPCestimated airborne contribution and the SEA-based airborne contribution to the interior noise showedsimilar spectral characteristics. In the SEA-based airborne contribution, the main order of motor noise
was well implemented, but the higher orders were not well implemented. Further improvements areneeded in this regard. The mainly contributing areas to the receiving position could be identified fromthe calculated SEA-panel contribution, and the interior noise countermeasurements could be studied
by considering the mainly contributing areas.

REFERENCES
1. RH Lyon, Statistical energy analysis of dynamical systems - Theory and Applications, MIT Press, 1975.
2. Z. Yu, S.T. Raveendra, H. Mo, C. Mollo, “Development of General SEA Template for Four-Door
Passenger Cars,” SAE 2003-01-1558, 2003.
3. K. S. Chae, J. H. Jeong, C. M. Park, and J. W. Yoo, “Prediction and Improvement of High Frequency
Road Noise of a Mid-Size Sedan”, SAE Technical Paper 2007-01-2307, 2007.
4. Thibault Lafont, Claudio Bertolini, Francesca Ronzio, Theophane Courtois, Davide Caprioli, "Application of Statistical Energy Analysis on a car: from the vehicle modeling to parts targeting", InterNoise 2016, 2016.
5. Eun Soo Jo, Ki-Sang Chae, Dong Chul Park, Wookeun Song, Andreas Schuhmacher, Moonju Hwang,“Auralization of electric vehicle's interior noise SEA simulation”, Proceedings of the InterNoise 2021,
2021
6. A. Schuhmacher and D. Tcherniak, “Engine Contribution Analysis Using a Noise and Vibration
Simulator”, International Conference on Noise and Vibration Engineering (ISMA2006) Vol. 43, 2006.
7. D. Tcherniak and A. Schuhmacher, “Application of decomposition-based technique in NVH source
contribution analysis”, Proceedings of the IMAC-XXVII (Florida, USA), 2009.
8. K. S. Chae and Y. H. Park, “An Efficient Exterior SPL Loading Measurment for SEA Airborne Prediction”,
Hyundai-Kia Motor Group R&D Internal Report, 2018.
9. K. S. Chae, B. H. Kim, J. H. Kim, and S. W. Lee, “Development and Application of a Simplified
Simulation Program of Vehicle Insulation Performance”, Hyundai-Kia Motor Group Conference, 2016.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0409

Prediction of sound radiation and transmission of corrugated panel
using wavenumber domain approach
Jungsoo RYUE1
1

University of Ulsan, Korea

ABSTRACT
Corrugated plates have a wide range of applications in various engineering fields. Corrugated plates areusually approximated into equivalent orthotropic plates, but this assumption is not appropriate at highfrequencies in which the local deformation of corrugation works. In this study, vibro-acoustic features ofcorrugated panels are investigated by using wavenumber domain approach, considering the corrugate panels
as waveguide structures. The results are validated by comparing them with those from the orthotropic plates.
Keywords: Corrugated panel, Sound radiation, Sound transmission

1. INTRODUCTION
In various branches of engineering, corrugated panels are widely used due to its increased bending
stiffness. Due to its one directional corrugation, corrugated panels can be simplified as orthotropic
panels (1~4), but this assumption does not work properly at high frequencies where the local
vibrations from each corrugated strip occur. For this high frequency range, numerical approaches suchas a finite element and boundary element method have been applied (5, 6). However, this conventional
FE/BE method is not much applicable for high frequency analysis because it takes heavy computing
burden due to the huge model size.
In this study, vibro-acoustic features of corrugated panels are investigated by using wavenumber
domain approach, named wavenumber finite element and boundary element (WFE/BE) method . Thismethod is for waveguide structures which have a uniform cross-section along its length (7, 8). So, in
this work, corrugate panels are considered as waveguide structures which are infinite in length. Using
this method, radiation ratio and sound transmission loss are predicted and the compared with thosefrom the approximate orthotropic panels.

2. PANEL MODELINGThe corrugated steel panel used in this study is shown in Figure 1(a), brought from a panel inreference (2). Since the panel is treated as a waveguide structure, the panel is assumed to be infinite
along the 𝑥 direction. At both ends in 𝑦 direction, the simply supported boundary condition is

(a)

(b)

Figure 1 – (a) Trapezoidally corrugated panel, (b) design parameters for corrugation
1

jsryue@ulsan.ac.kr

Table 1 – Geometry of the corrugation
𝑏2 , 𝑏4 , mm

𝑏3 (= 2𝑏1 = 2𝑏5), mm

𝜙5 , 𝜙5 , degree

ℎ, mm

𝐻, mm

𝐿𝑐 , mm

95

46

50

0.65

55

250

applied. The geometrical property of the corrugation is shown in Figure 1(b) and Table 1. Material
properties are set to nominal values of mild steel. That is, 𝐸 = 200 GPa, 𝜌 = 7800 kg/m3 , 𝜈 =0.3 and damping factor is set to 𝜂 = 0.05.
2.1 Orthotropic Panel Model
The corrugated panel has different bending stiffnesses in 𝑥 and 𝑦 directions. For the vibration
at low frequencies in which the bending wavelength in 𝑦 direction is much larger than the period of
the corrugation profile, the panel’s bending stiffnesses in the two orthogonal directions can be written
as (1~4)𝐸ℎ

𝐵𝑥 = (1−𝜈2 )𝐿 ∑𝑛 𝑏𝑛 (𝑧𝑛2 +𝑦ℎ2 +𝑏𝑛224+ℎ2 −𝑏𝑛224

cos 2𝜙𝑛 ) , 𝐵𝑦 =

𝐸ℎ3

𝐿𝑦

12(1−𝜈2 ) ∑𝑛 𝑏𝑛(1)

The governing equation of the orthotropic panel is
𝐵𝑥

𝜕4𝑤
𝜕𝑥 4

+ 2𝐵𝑥𝑦

𝜕4𝑤
𝜕𝑥 4

+ 𝐵𝑦

𝜕4 𝑤
𝜕𝑥 4

+𝜇

𝜕4𝑤
𝜕𝑥 4

=𝐹

(2)

where 𝑤 denotes the vertical displacement, 𝐵𝑥𝑦 = √𝐵𝑥 𝐵𝑦 and 𝜇 denotes mass per unit area. For
the panel shown in Figure 1 and Table 1, two bending stiffnesses are 𝐵𝑥 = 38.6 kNm, 𝐵𝑥 = 4.5 Nm
and the mass per unit area is 𝜇 = 8692 kg/m2 .
2.2 Wavenumber Finite Element and Boundary Element (WFE/BE) Model
The schematic diagram of WFE/BE modelling for the corrugated panel is shown in Figure 2. The
governing equation for the wavenumber finite element which models the structural part is given bỹ = 𝐅̃
[𝐊 4 (−𝑖𝑘𝑥 )4 + 𝐊 2 (−𝑖𝑘𝑥 )2 + 𝐊1 (−𝑖𝑘𝑥 ) + 𝐊 0 − 𝜔2 𝐌] 𝚽
(3)where 𝐊 𝑖 are the stiffness matrices, 𝐌 is the mass matrix, 𝜔 is the angular frequency, 𝑘𝑥 is thẽ and 𝐅̃ denote displacement and external force vectors in wavenumber
𝑥-directional wavenumber, 𝚽domain, respectively. The wavenumber boundary element for the fluid (air) is expressed by
̃
𝜕𝚿

̃
𝐏

̃ − 𝐆(𝜔, 𝑘𝑥 ) = in
𝐇(𝜔, 𝑘𝑥 )𝚿
𝜕𝐧
𝑗𝜔

(4)

̃ is the velocity
where 𝐇 and 𝐆 are the matrices for the Green’s function, stiffness matrices, 𝚿̃ /𝜕𝐧 is the fluid particle velocity into the normal direction, 𝐧 , and 𝐏
̃in is thepotential and 𝜕𝚿
incident pressure vector. In this modelling, two acoustic domains are defined: Upper semi -infinite
space and lower semi-infinite space with baffled condition. Since the WFE and WBE are coupled on
the upper and lower surfaces, the coupled WFE/BE equations arẽ = 𝐅̃ + 𝑖𝜔𝜌0 𝐂1 𝚿̃[𝐊(𝑘𝑥 ) − 𝜔2 𝐌] 𝚽
(5)
̃

̃

(6)
where
and 𝐂2 are matrices allocating thecoupled WFEs and WBEs, 𝐈 denotes the identity matrix, 𝜌0 is the fluid density. By solving Eqs
̃ (𝜔, 𝑘𝑥 ) , fluid pressure 𝒑̃ and fluid velocitỹ (𝑘𝑥 ) = 𝑖𝜔𝜌0 𝚿(4)~(6), the panel displacement 𝚽̃̃(𝜔, 𝑘𝑥 ) = −𝜕𝚿/𝜕𝐧 can be found.
𝒗
𝐈 𝜕𝚿/𝜕𝐧 − 𝑖𝜔𝐂2 𝚽 = 𝟎
𝐊 = 𝐊 4 (−𝑖𝑘𝑥 )4 + 𝐊 2 (−𝑖𝑘𝑥 )2 + 𝐊1 (−𝑖𝑘𝑥 ) + 𝐊 0 , 𝐂1

Figure 2 – Schematic diagram of the WFE/BE modelling for the corrugated panel

3. PREDICTION OF VIBRO-ACOUSTIC PERFORMANCE
3.1 Dispersion Relations
In order to identify the waves which exist in the corrugated panel, dispersion relations are obtained
as shown in Figure 3 by solving Eq.(3) for 𝐅̃ = 0. For comparison, dispersion curves obtained fromthe equivalent orthotropic panel (green lines) and a strip for flange in the corrugation (red lines) areplotted in Figure 3. For the flange strip, a reduced width with a factor of 0.9 for the 1 st wave and 0.8
for the 2 nd wave are adopted (9). Also the acoustic wave is illustrated as a red dashed line in Figure 3.
From the dispersion curves, it is obvious that the orthotropic approximation is only valid up to about
200 Hz. From about 300 Hz, local waves guided within flange strips are generated. Hence, the local
vibrations of flange and web should be included for the high frequency analysis.
From the wave curves of orthotropic panel and acoustic wave, it can be found that the firstcoincidence frequency is around 250 Hz. It is expected that the radiation feature of the panel will bechanged abruptly around this frequency.

Figure 3 – Dispersion relations of the corrugated panel. (WFE, blue lines; Orthotropic panel,green lines; flange strip, red lines)
3.2 Radiation Ratio
As shown in Figure 2, a point force is applied on the upper side of the panel and the sound radiation
takes place from the upper side of the panel. In this situation, the panel is assumed to be baffled for
𝑦 < 0 and 𝑦 > 𝐿𝑦 . The radiation ratio is given by2

𝜎 = 𝛱𝑟𝑎𝑑 /𝜌0 𝑐𝐿𝑦 ⟨𝑣 ⟩inf(7)
2

where 𝛱𝑟𝑎𝑑 is the radiated power, 𝑐 is the speed of sound, ⟨𝑣 ⟩inf is the time and space averaged
squared velocity of the panel.
In this calculation, the location where the vertical (𝑧-direction) force is applied on is divided into
three categories: on the flange, on the web and on the junction where the flange and web are connected.
This division is introduced based on the fact that the local bending stiffness at these three locations
are quite different. Within each group of excitation locations, the point force is applied all the possible
locations by shifting it. The radiation ratios predicted by WFE/BE and the equivalent orthotropic panel
are illustrated in Figure 4. It can be seen from Figure 4 that the WFE/BE results are quite dependingon the excitation locations, especially at high frequencies where the local vibration from corrugation
takes place. In addition, it is sure that the orthotropic panel model is not useful at high frequencies,
particularly for the flange excitation. In case of junction excitation, Figure 4 shows that the corrugatedpanel tend to behave more like the orthotropic panel even at high frequencies.
To include the local vibration of the corrugation into the approximate model, two additional
simplified models are proposed here and combined with the existing orthotropic panel. First one is anequivalent orthotropic panel corresponding to a single corrugation of the panel (This is for middle
frequency range); Second one is a single plate strip corresponding to the flange which is the widest
strip (This is for high frequency range). The radiation ratio predicted from this revised approximate
model is shown in Figure 4 and it is found that the approximate model gives very good agreement

with the WFE/BE result for flange excitation.

Figure 4 –Radiation ratios predicted for the corrugated panel
3.3 Sound Transmission Loss
The sound transmission loss through the panel is defined by
𝑇𝐿 = −10 log10 𝜏𝑓(8)where 𝜏𝑓 denotes the power transmission coefficient. For an acoustic incidence into the orthotropic
panel as shown in Figure 5(a), the power transmission coefficient can be obtained by (1,2)
2𝜋/2𝜏𝑓 = 𝜋 ∫0

𝜃

−2

𝑍

0
𝑟
∫0 |1 + 2𝜌 𝑐 cos 𝜃|0𝑑(sin2 𝜃)𝑑𝜙(9)where 𝑍𝑟 is the surface impedance of the panel(1,2) and 𝜃0 denotes upper limit of the incident angle
for field incidence which is 78 degrees. On the other hand, the power transmission coefficient in
WFE/BE analysis is expressed as (8)
𝜋

𝜋

𝜋

0

0

𝜋

𝜏𝑓 = ∫𝛽2 ∫𝛼2 𝜏(𝛼, 𝛽) sin2 𝛽 sin 𝛼 𝑑𝛼𝑑𝛽 / ∫𝛽2 ∫𝛼2 sin2 𝛽 sin 𝛼 𝑑𝛼𝑑𝛽
0

(10)

0

where 𝛼0 and 𝛽0 correspond to the lower limit of the incident angle, which is 12 degrees in theconfiguration in Figure 5(b).
The sound transmission losses predicted by WFE/BE method and orthotropic panel model areillustrated in Figure 6(a), together with the measured one which is in reference (2). Figure 6(a) showsthat the TL from the WFE/BE method has very good agreement with the measured result. To
investigate the transmission paths through the corrugation, the TLs through the flange and web are
separated and plotted in Figure 6(a). From this separation, it is found that input sound is transmitted
mainly through the flange at low frequencies but as frequency increases sound transmission throughthe web rises.

(b)

(a)

Figure 5 – Coordinates and parameters for acoustic excitation. (a) Orthotropic panel model, (b) WFE/BEmodel.

(a)(b)

Figure 6 – (a) Comparison of sound transmission losses, (b) Image plot of transmitted power overlaid on
dispersion curves in frequency-wavenumber domain.
To reveal the contribution of flange and web more clearly, the transmitted power is illustrated infrequency-wavenumber domain as an image plot in Figure 6(b) overlaid on the dispersion diagram. In
Figure 6(b), the bright region implies the strong transmission, which means TL becomes lower in thebright region. By comparing the results in Figure 6(a) and (b), it is clear that local vibrations of flange
and web drops transmission loss.

4. CONCLUSIONS
Corrugated panels are commonly approximated as orthotropic panels. However, this model is not
valid at high frequencies where local modes of the flange and web occur. Furthermore, theconventional finite element and boundary element method has a limitation in high frequency analysisof the corrugated panel. In this study, vibro-acoustic features of corrugated panels were predicted bymeans of the WFE/BE method and confirmed that WFE/BE approach is a useful tool for this purpose.
Also, in the radiation ratio prediction, a revised approximate model is proposed in this study, which
consists of three plate models: orthotropic global plate (for low frequency), orthotropic plate
corresponding to a single corrugation (for mid-frequency) and local plate strip for flange (for highfrequency). For sound transmission, the WFE/BE method showed a good agreement with experiments
and the contribution of the local vibrations in flange and web was well described. Further investigation
will be carried out to find a simple prediction model on sound transmission. Also, statistical energy
analysis (SEA) for corrugated panels will be attempted as a future work.

ACKNOWLEDGEMENTS
This research was supported by the Basic Science Research Program through the National Research
Foundation of Korea (NRF), funded by the Ministry of Education (2019R1F1A1046725).

REFERENCES
1. Hansen CH. Sound transmission loss of corrugated panels. Noise Con. Eng. 1993; 40(2); 187-197.
2. Windle RM. The Prediction of Industrial Noise and its Transmission through Metal Cladding Systems.
PhD Thesis, University of Salford, UK; 1995.
3. Hopkins C. Sound Insulation. New York, USA: Routledge; 2007.
4. Bies DA, Hansen CH. Engineering Noise Control: Theory and Practice. London, UK: Spon Press; 2003.
5. Samanta A. Mukhopadhyay M. Finite element static and dynamic analysis of folded plates. Engineering
Structures 1999;21:277-287.
6. Piscoya R., Ochmann M. Transmission loss of corrugated plates. Forum Acusticum 2011; 27 June-1 July
2011; Aalborg, Denmark 2011. p. 211-215.
7. Nilsson CM, Jones CJC, Thompson DJ, Ryue J. A waveguide finite element and boundary element

approach to calculating the sound radiated by railway and tram rails. J. Sound Vib. 2009; 321;813-836.
8. Kim H, Ryue J, Thompson DJ, Müller AD. Application of a wavenumber domain numerical method tothe prediction of the radiation efficiency and sound transmission of complex extruded panels. J. Sound
Vib. 2019; 449; 98-120 .
9. Li H, Squicciarini G, Thompson DJ, Ryue J, Xiao X, Yao D, Chen J. A modelling approach for noisetransmission through extruded panels in railway vehicles. J. Sound Vib. 2021; 502; 116095.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0118

Improved identification of sound sources at low frequencies using the twostep steepest descent method based on the equivalent source method
Laixu Jiang 1; Guangping Zou 1; Jeong-Guon Ih 2
1
2

College of Aerospace and Civil Engineering, Harbin Engineering University, China

Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Korea

ABSTRACT
Source identification is challenging when the underdetermined systems consist of relatively fewer
measurement data than the sources to reconstruct. Wideband acoustic holography (WBH) is a sparse
regularization-based source identification method with satisfactory reconstruction accuracy. In the previous
work, it is found that the initial source strength distribution significantly affects the results, so it is challenging
to identify low-frequency coherent sources using WBH alone. In this work, the improved WBH is proposed,
motivated to refine its deficiencies at low-frequency reconstruction. The initial source strength of WBH ismodified by introducing the results of the fastest descent method. Combined with the a priori assumption of
the source strength, label the equivalent source of the actual source and adjust the threshold filtering
conditions to ensure that these equivalent source strengths are not involved in the filtering process during theiteration. The present work shows that the reconstruction accuracy of the modified WBH is similar to that of
the WBH at mid-high frequencies. Meanwhile, it can bear better reconstruction and localization accuracy at
low frequencies without degrading the reconstruction performance at the mid to high-frequency sound field.
Keywords: Acoustical holography, Underdetermined situation, Iterative methods, Data compression

1. INTRODUCTION
The reconstruction of acoustic parameters and localization for sound sources in the sound field canhelp engineers to investigate the vibro-acoustic mechanism and noise control. Near-field acoustic
holography (NAH) has been widely used for sound field reconstruction and localization [1 , 2].
Recently, the wideband acoustical holography method (WBH) has been proposed [3, 4], which canhelp reduce the number of measuring points further for source localization and extend the frequency
range. Unlike the conventional ESM based on the L 2-norm regularization, the WBH method uses thegradient iteration algorithm that includes a dynamic filtering process to enforce sparsity in themonopole source model for finding the sparse solution of the equivalent source amplitude. The
detailed solving process can be referred to in Ref. [3]. In this way, the reconstruction accuracy of
WBH is much better than conventional ESM in the middle to high frequency. However, due to theinherent restrictions of the dedicated iterative solvers in the WBH method, the reconstruction accuracy
at low frequencies is not satisfactory as that at mid and high frequencies, especially when dealing withcoherent sources.
This paper proposes an improved WBH method, called the MWBH, which aims to refine the WBH
method to improve the reconstruction accuracy at low frequencies. Firstly, the influence of the initialsource strengths on the reconstruction results in the WBH method is investigated. Next, based on thiseffect, the steepest descent method is introduced to first reconstruct the sound field as the priorassumption of the source strength distribution. Finally, the calculated results of the source strength
replace the zero-source strength as the initial value to start the iteration process. Also, the threshold
filter condition is adjusted to ensure that these equivalent source strengths of the actual source strength
are not involved in the threshold filtering process during the iteration. A simulation example with twocoherent pulsating spheres and two loudspeaker experiments are used to evaluate this method. Since
the WBH method is also developed based on the steepest descent method, the suggested method iscalled the equivalent source method based on the two-step steepest descent method.
1
2

lxjiang92@gmail.com
jgih21@gmail.com

2. MWBH description
In the iterative solving process of the WBH method, the setting of two parameters is very
influential to the results: one is the initial source strength, 𝐪0 , which is the input for the beginning of
the iteration, and the WBH method sets the initial source strength to zero, 𝐪0 = 0. The second is thethreshold filter condition, which sets the equivalent source strength below the threshold to zero, andthose above the threshold are preserved while the sidelobe is removed. The algorithm first sets ahigher threshold and gradually decreases the threshold in the iterative process to obtain a sparse
equivalent source strength solution with suppressed sidelobe. The threshold setting in the WBH
method is related to the highest source strength obtained at each iteration, and the threshold function
is defined as [3]Tk = 10− Dk
qk +1 , max
20

(1)

,

where 𝐷𝑘 denotes the dynamic filtering range with the suggested as 60 and |𝐪k+1, max | is theamplitude of the largest element in |𝐪𝑘+1 |. Accordingly, the approximate solution |𝐪𝑘+1 | after (𝑘 +1)-th iteration step can be written as [3]q ,qk +1, i =  k +1, i 0,qk +1, i  Tk

if

(2)

.

othewiseThis work focuses on the above two parameters, initial source strength and threshold filter condition,
to refine them accordingly.
2.1 Effect of initial source strength
Here, three kinds of source strength values are employed as the initial source strengths in the WBH
method to examine the effect of initial source strengths on the WBH results: (1) zero vector (original
WBH method); (2) result of Tikhonov regularization method (TR-WBH); (3) result of steepest descentmethod (SDM-WBH). Two coherent pulsating spheres are taken as the test example. The radius of thepulsating sphere is 𝑎 = 0.01 𝑚 , and its velocity amplitude is 𝑢𝑎 = 0.01 𝑚/𝑠 . The sound pressureradiated from a pulsating sphere is expressed asP=

0c0 ka 2

R(1 + ka )2ua (kR + j )e− jkR ,(3)

where 𝜌0 is the air density, 𝑐0 the speed of sound, k the wavenumber, and R the distance from thepulsating sphere to the field point. Locations of two coherent pulsating spheres are (-0.16, 0, 0) m and(0.16, 0, 0) m, respectively. The hologram, reconstruction, and equivalent source planes are locatedat 0.3 m and 0.05 m above the plate and 0.1 cm below the pulsating sphere. The spacing of hologram
points is 0.06 m on a 0.3 m × 0.3 m plane. The reconstruction point spacing is 0.02 m over a 0.40 m
× 0.40 m reconstruction plane. The equivalent source plane grid is the same as the reconstruction
plane. Uniform random noise with a signal-to-noise ratio of 30 dB is added to the pressure data forthe simulation, reflecting an actual measurement situation with noise contamination.

Reconstruction error, /%

70
60
50
40
30
20
10
0

500

1000

1500

2000

2500

Frequency, Hz

Figure 1 – A comparison of the reconstruction errors with different initial values :TR-WBH;, SDM-WBH, WBH;,

The reconstruction error is presented in the root-mean-squared error norm (𝑒𝑟𝑚𝑠 , in )) given by𝑒𝑟𝑚𝑠 =̂𝑠 ‖‖𝐏𝑠 −𝐏2
̂𝑠 ‖‖𝐏2

̂𝑠 denote the theoretical and the reconstructed pressure,× 100% , where 𝐏𝑠 and 𝐏

respectively. The reconstruction errors of the WBH with different initial values are shown in F ig. 1.
Figure 1 reveals that the reconstruction errors of the TR-WBH and SDM-WBH methods are
significantly lower than the WBH method at lower frequencies. From here, it can be concluded thatthe initial source strength greatly influences the calculation results, and a reasonable a priori
assumption of the sound field distribution is important for improving the reconstruction performancein the WBH method. Therefore, the MWBH algorithm uses the preliminary calculation by the steepest
descent method at the start of the iteration, 𝐪0 = 𝐪𝑠𝑑𝑚 , instead of the zero initial value, 𝐪0 = 0.
2.2 Refinement of the threshold filter condition
In the WBH method, the source strength below the threshold is set to zero in each iteration by
setting the threshold filter function to suppress the sidelobe. It also introduces uncertainty: theequivalent source strengths of the actual sources should be set to zero in the iterations. Although thesesource strengths are recovered in subsequent iterations, it creates a bias in the results. Therefore, theidea of the MWBH method on the threshold filtering condition is to label the equivalent source
strengths of the actual source. These source strengths are not involved in the threshold filtering processduring the iterative process. The MWBH algorithm is implemented as follows: 1) Assume there are
M hologram points and N equivalent sources. The relationship between the sound pressure at thehologram point and the equivalent source strength is an underdetermined matrix, i.e., M < N; 2)
Employ the steepest descent method to obtain the equivalent source strength vector 𝐪𝑠𝑑𝑚 , thenarranging the equivalent source strength vectors in descending to obtain the subscript index 𝐼𝑚 forthe first m values, retain them, where m is one-fourth of the number of equivalent sources [5]; 3)
Replace the initial value 𝐪0 = 0 with 𝐪0 = 𝐪𝑠𝑑𝑚 , and then solving the equivalent source strengths
vector by the WBH method. During the iteration process, the threshold filtering condition in Eq. (2)is changed toq ,qk +1, i =  k +1, i 0,

if qk +1, i  Tk or i  I mothewise(4);4) Obtain the output equivalent source strength vector, q.

3. Numerical Simulations with coherent pulsating spheres
The reconstruction performance of MWBH is tested to reconstruct the coherent source sound field.
The simulation conditions are the same as in Section 3.1. Figure 2 compares the reconstruction results
using the WBH and the MWBH methods at different measurement distances (d = 0.1 m, 0.2 m, 0.3 m,
and 0.4 m).

Figure 2 – A comparison of reconstruction errors using the WBH method and MWBH method. The
measurement distance is varied as (a) d = 0.1 m, (b) d=0.2 m, (c) d=0.3 m, (d) d=0.4 m:, WBH
method;, MWBH.

It can be seen from Fig. 2 that, as expected, the reconstruction accuracy of the WBH method at highfrequencies is still acceptable. Still, the reconstruction accuracy at low frequencies is not good.
Besides, the WBH method is influenced very much by the measurement distance, and the sound fieldreconstruction performance is lower when the measurement distance is far. Conversely, thereconstruction error of MWBH is less than 15) at the four-measurement distance for the analyzed
frequencies. Especially at low frequencies, the MWBH method dramatically improves thereconstruction accuracy with the WBH method and is less affected by the measurement distance.

4. Experimental verification with coherent sound sources
A dual loudspeaker with a spacing of 0.30 m is used as the sound source in an anechoic chamber
to validate the proposed method. Two loudspeakers produce the same phase, and a microphone arrayof 0.3 m × 0.3 m collects the sound pressure signal for a span of 500-2500 Hz. The microphone arrayof 36 measuring points is 0.30 m from the speaker plane. The sound pressure on the reconstruction
plane at 0.02 m in front of the loudspeaker is reconstructed using WBH and MWBH, respectively.
The reconstruction results are shown in Fig. 3. One can see that the two sources can be identified
using both WBH and MWBH methods at 2.5 kHz. However, when the analysis frequency is low, e.g.,
500 Hz, the WBH method identifies a strong source in the middle of the two sources. In contrast, the
MWBH method can still identify two sources separately. The result reveals that, compared with the
WBH method, the MWBH method can significantly improve the reconstruction accuracy of coherent
sound sources at low frequencies without degrading the reconstruction performance at the mid to highfrequency sound field.

Figure 3 –The reconstructed sound pressure distribution on the reconstruction plane using WBH and
MWBH methods at 500 Hz and 2500 Hz.

5. CONCLUSIONSIn this paper, a new method called the MWBH method is proposed to solve the low-frequency
problem of the WBH method, which is to improve the reconstruction accuracy. Based on the fact thatthe initial source strength setting in the WBH method has a significant impact on the reconstruction
accuracy, the fastest descent method is introduced to obtain the source strength. The filtering thresholdcondition in the WBH method is modified to ensure that the actual source strength is not to be zero
during the iterative process. Simulation with two coherent pulsation spheres shows that the MWBH
can significantly improve the reconstruction accuracy at low frequencies compared to the WBH
method for the model and frequencies tested in this paper. Moreover, comparing the reconstruction
results under four different measurement distances reveals that the reconstruction accuracy of the
MWBH method is relatively less affected by the increase of the measurement distance compared to
the WBH method. The experiment with two speakers also validates the findings in the simulation

without degrading the reconstruction performance of the mid-high frequency compared to the WBH
method.

ACKNOWLEDGEMENTS
This work was partially funded by the China Scholarship Council (No. 202006680090) and the
NRF of Korea (No. 2020R1I1A2066751).

REFERENCES
1. Williams EG. Fourier Acoustics: Sound Radiation and Near-Field Acoustical Holography. Academic
Press, 1999.
2. Jeon IY, Ih JG. On the holographic reconstruction of vibroacoustic fields using equivalent sources andinverse boundary element method. J. Acoust. Soc. Am. 2005; 118 (6): 3473-3482.
3. Hald J. Fast wideband acoustical holography. J. Acoust. Soc. Am. 2016; 139(4): 1508–1517.
4. Hald J. A comparison of iterative sparse equivalent source methods for near-ﬁeld acoustical holography. J.
Acoust. Soc. Am. 2018; 143(6): 3758–3769.
5. Blumensath T, Davies M E. Iterative hard thresholding for compressed sensing. Applied andComputational Harmonic Analysis. 2009; 27(3): 265-274.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0176

Identification of acoustical modes inside a flanged cylindric duct
based on out-duct measurement
Kang GAO1; Weikang JIANGThe State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, China

ABSTRACT
Mode identification is widely applied for analyzing the noise inside ducts. Conventional out-duct
measurement methods could only identify radiation modes, which contain the mixture of both incident modes
and reflective modes. A methodology based on an out-duct placed microphone array is proposed to identify
and decompose incident and reflective modes in-duct, by which acoustic modes in the finite cylindric ductwith an infinite flange are studied to obtain forward and backward modes solutions in closed-forms.
According to the generalized radiation impedance theory, complex reflection coefficients between incident
modes and reflective modes could be represented as an infinite matrix equation, which is unable to solve
practically. A simplified and solvable equation is developed to solve the complex reflection coefficients byonly considering self-reflection for a single incident mode, after which the propagation models of bothforward and backward modes are established based on Tyler and Sofrin’s model. Through the measurementof an out-duct microphone array, in-duct forward modes could be identified by the inverse method. In-ductsound sources are furtherly reconstructed through superpositions of incident modes identification results. The
effectiveness of the presented method is verified through comparisons of in-duct measurements from both
numerical simulations and tests.
Keywords: Duct modes, Mode identification, Complex reflection coefficients, Out-duct measurement

1. INTRODUCTION
Rapidly developed air transportation industry has brought grievous noise pollution and damaged
the physical and psychological comfort of cabin passengers. Noise reduction in civil aircraft has beena never-ending goal for achievement to designers and manufacturers. Among the noise sources of theairliner, the aero-engine fan is always one of the dominant sources [1].
Aero-engine fans are installed inside the nacelle, which is usually simplified as a cylindrical duct.
According to duct acoustics, any sound field inside a duct could be represented by the superpositionof duct modes [2]. That means the physical mechanism, propagation, and radiation properties of the
fan noise could be analyzed through duct modes. Through the measurement of microphones, ductmode amplitudes could be calculated by solving a matrix inverse problem, this process is known as
duct mode identification, which is a widely used experimental method for analy zing inner duct noises.
Mode identification (MI) includes Azimuthal Mode Identification (AMI) [3] and Radial Mode
Identification (RMI) [4]. AMI is often carried out through a ring microphone array [5] and calculated
easily by spatial Fourier transformation into the mode domain [6]. For RMI, the sound information
along the duct axis or radial direction is measured sufficiently to decompose the radial modes, whichmeans the microphones need to be distributed along the duct axial [7] or radial direction [8] .
Limited by the duct geometry, microphones for MI could only be mounted on the duct w all, theduct cross-section, or be placed outside the duct. In contrast to in-duct measurement, measurementwith an out-of-duct array is more attractive, which causes no intrusion to the in-duct physical fields
as well as the duct wall [1] and is more flexible on the spatial layout. The first attempt of applying outduct measurement to MI was performed by Thomas et al. [9] , who adopted a rotatable and axially
movable ring microphone array placed near the outlet of the duct. Lan et al. [10] firstly applied
microphones located on the turbulence control screen (TCS) to measure the radiation sound field andaccomplished both AMI and RMI. That method was then applied to MI of a full-scale engine inlet [11].
1gk_sjtu@sjtu.edu.cn

Castres and Joseph [12] proposed an inverse technique for MI using the sound pressure measured by
microphones on TCS. The robustness and accuracy of this inverse process were analyzed after which
the method was investigated through experiments on a fan rig and the effects of different arraygeometries on the modal recognition results are compared [13] . Lewy [14] developed a method to deduce
the amplitude of primary modes through an arc array measured in the far field. Yu and Huang [15] usedordinary out-duct linear arrays to detect modes, and the performance of the inverse reconstruction was
improved through the compressive sensing method.
However, current MI methods based on out-duct measurement actually could only identify
radiation modes that contain the mixture of both incident modes and reflective modes. When modes
generated by in-duct sources propagate to the duct outlet, reflective modes would arise, and bothincident modes and reflective modes on the cross-section at the duct outlet would radiate to external
space. To analyze the physical mechanism and propagation of in -duct sources, the incident modes
should be identified accurately.
The objective of this study is to develop a methodology to identify and decompose incident andreflective modes inside the duct utilizing an exterior placed microphone array. A monopole inside a
hard-wall cylindrical duct with an infinite flange is considered to derive a closed-form solution of thetransfer function from the source to the near-field microphones. The generalized radiation impedancetheory [16,17] is applied to calculate complex reflection coefficients between the incident and reflectivemodes. The radiation modes are identified firstly by the inverse techniques, after which the incident
modes are extracted. Finally, in-duct sound sources are reconstructed through superpositions ofincident modes. The effectiveness of the presented method is verified through numerical simulationsin Section 3 and experiments in Section 4.

2. THEORY ANALYSIS
2.1 Sound Field in A Flanged Duct of Semi-infinite Length
Considering a semi-infinite circular duct with a flange shown in Figure 1, both of the duct walls
are perfectly rigid. The radius of the duct is R. The sound field is stimulated by a monopole locatedat (rs ,θs ,zs ) represented in cylindrical coordinates with its volume velocity amplitude Q0. According
to the Helmholtz equation, the complex acoustic pressure in-duct is given as follows [18]∇2 𝑝 + 𝑘 2 𝑝 = −𝑗𝜌𝜔Q0 𝛿(𝒙 − 𝒙𝒔 )(1)
Where ∇ is the gradient operator.𝑘 and 𝜔 is the wavenumber and angular frequency respectively, 𝜌means the density of the medium, 𝑗 is the imaginary unit and 𝛿 is the Dirac function. 𝒙𝑠 = (𝑟𝑠 , 𝜃𝑠 , 𝑧𝑠 )represents the vector form of the source’s location as well as 𝒙 = (𝑟, 𝜃, 𝑧) represents the vector formof the location to be solved. Eq. (1) could be solved by modes decomposition as follows+∞

𝑝(𝑟, 𝜃, 𝑧) = ∑

+∞

∑
𝑚=−∞

(𝐴𝑚𝑛 𝑒 −𝑗𝑘𝑧 𝑧 + 𝐵𝑚𝑛 𝑒 𝑗𝑘𝑧 𝑧 ) 𝐽𝑚 (𝑘𝑚𝑛 𝑟)𝑒 𝑗𝑚𝜃

(2)

𝑛=0

herein, 𝑚 and 𝑛 represent the azimuthal and radial order of mode, respectively. Under the coordinate
system shown in Figure 1, 𝐴𝑚𝑛 and 𝐵𝑚𝑛 represents the modal amplitude of the incident mode andthe reflective mode, respectively. 𝐽𝑚 denotes the first-kind Bessel function at order m. 𝑘𝑚𝑛 is′ (𝑘decided by the boundary condition, and for the rigid duct wall, 𝑘𝑚𝑛 satisfies 𝐽𝑚
𝑚𝑛 𝑅) = 0, in which“′” means derivation. The sound field is stimulated by a monopole located at (𝑟𝑠 , 𝜃𝑠 , 𝑧𝑠 ) in-duct, andthe mode amplitude of incident mode becomes as follows(3)
𝐴𝑚𝑛 = 𝜌𝜔𝑄0 𝐽𝑚 (𝑘𝑚𝑛 𝑟𝑠 )𝑒 −𝑗𝑚𝜃𝑠 𝑒 𝑗𝑘𝑧 𝑧𝑠 /(2𝑘𝑧 𝛬𝑚𝑛 )The term 𝛬𝑚𝑛 is introduced as a factor to satisfy the mode normalization condition. Only propagativemodes are considered. According to the wave equation and Eq. (1), the acoustic velocity along theaxial direction at the cross-section of the outlet could be expressed as follows𝑚+𝑛0
1𝑣𝑧 (𝑟, 𝜃) = ∑∑ 𝛺𝑚𝑛 (𝐴𝑚𝑛 − 𝐵𝑚𝑛 ) 𝐽𝑚 (𝑘𝑚𝑛 𝑟)𝑒 𝑗𝑚𝜃
(4)
𝜌𝑐
𝑚=𝑚−
𝑛=0
Where 𝑐 is the sound velocity and 𝛺𝑚𝑛 = 𝑘𝑧 /𝑘.
Herein, the wavelength of the sound is considered to be much smaller than the dimensions of theflange, which means the flange could be regarded as an infinite one. Zorumski [16] derived thetheoretical relationship between the mode amplitude 𝐴𝑚𝑛 from the incident mode (𝑚, 𝑛) and mode
amplitudes from the reflective modes are as follows

𝐵𝑚𝑙 = 𝑅𝑚𝑙𝑛 𝐴𝑚𝑛 (𝑙 = 1,2,3, ⋯ , +∞)(5)
Where 𝐵𝑚𝑙 is the mode amplitude of reflective mode (𝑚, 𝑙), and 𝑅𝑚𝑙𝑛 is the complex reflective
coefficient which decays as the increase of 𝑙.

Figure 1– Schematic diagram of a semi-infinite length flanged duct
For one incident mode propagating through the open end of the duct, multiple modes with the sameazimuthal order will be excited according to Eq. (5), in which 𝑙 must be truncated at some finite
value and solved. In an ideal situation, as the growth of 𝑙, 𝑅𝑚𝑙𝑛 decays rapidly to zero. However,according to the calculation results of 𝑅𝑚𝑙𝑛 by Norris et al. [19] and Wendoloski et al. [20], theconvergence of this calculation is fairly slow. Only self-reflection is considered here for a single
incident mode, followed by
𝐵𝑚𝑛 = 𝑅𝑚𝑛𝑛 𝐴𝑚𝑛(6)Though Eq. (6) is significantly simplified compared to Eq. (5), the computational accuracy of 𝐵𝑚𝑛
is still assured owing to the calculation results of Shao et al. [17] which indicated that the self-reflectionof an incident mode is predominant in most instances. According to Zorumski [16] , 𝑅𝑚𝑛𝑛 solves asfollows
𝑅𝑚𝑛𝑛 = (1 + 𝑍𝑚𝑛𝑛 𝛺𝑚𝑛 )/(𝑍𝑚𝑛𝑛 𝛺𝑚𝑛 − 1)
(7)[17]
Where 𝑍𝑚𝑛𝑛 is the generalized radiation impedance. According to Shao et al.
, 𝑍𝑚𝑛𝑛 could be
calculated as follows
𝑍𝑚𝑛𝑛 = −1𝜋⁄2′2 (𝑘𝑅∞′2 (𝑘𝑅1
sin 𝜑 3 𝐽𝑚
sin 𝜑)
cosh 𝜉 3 𝐽𝑚cosh 𝜉)
{∫
𝑑𝜑
+
𝑗
𝑑𝜉}
∫
2
2
2
2
2
2
2 2
𝑊 0 [sin 𝜑 − (𝑘𝑚𝑛 ) ⁄𝑘 ]0 [cosh 𝜉 − (𝑘𝑚𝑛 ) ⁄𝑘 ]Where 𝑊 = 2 (1 − (𝑘

𝑚2

2
𝑚𝑛 𝑅)

(8)

).

2.2 Radiation Field Out of The Duct
According to Eq. (4), the acoustic velocity along the axial direction at the cross -section of the
open-end for one mode (𝑚, 𝑛) could be denoted as(9)𝑣𝑚𝑛 = 𝑉mn 𝐽𝑚 (𝑘𝑚𝑛 𝑟)𝑒 𝑗𝑚𝜃
Where 𝑉mn = 𝛺𝑚𝑛 (𝐴𝑚𝑛 − 𝐵𝑚𝑛 )/𝜌𝑐 is the axial velocity mode component. Consider one measurement
position 𝑃 locates at (ℎ𝑝 , 𝜃𝑝 , 𝜑𝑝 ) denoted by spherical coordinates shown in Figure 1.The soundpressure radiated by a single-mode is calculated by the Rayleigh integral [21] as follows𝑗𝜌𝑐𝑘
𝑝𝑚𝑛 (ℎ𝑝 , 𝜃𝑝 , 𝜑𝑝 ) =𝑉 𝑒 −𝑗𝑘ℎ𝑝
2𝜋ℎ𝑝 𝑚𝑛(10)
2𝜋𝑅× ∫ 𝑑𝜃 ∫ 𝐽𝑚 (𝑘𝑚𝑛 𝑟)𝑒 𝑗[𝑚𝜃+𝑘𝑟(𝑐𝑜𝑠 𝜃 𝑐𝑜𝑠 𝜃𝑝 +𝑠𝑖𝑛 𝜃 𝑠𝑖𝑛 𝜃𝑝 ) 𝑠𝑖𝑛 𝜑𝑝 ] 𝑟𝑑𝑟
00Using the properties of the first-kind Bessel function, Eq. (10) could be furtherly written following
𝑠𝑖𝑛 𝜑𝑝 𝐽𝑚 ′ (𝑘𝑅 𝑠𝑖𝑛 𝜑𝑝 ) 𝑒 −𝑗𝑘ℎ𝑝 𝑗𝑚(𝜃𝑝 −𝜋)2𝑝𝑚𝑛 (ℎ𝑝 , 𝜃𝑝 , 𝜑𝑝 ) = 𝑗 2𝑚+1 𝑘 2 𝑅𝑉𝑚𝑛 𝐽𝑚 (𝑘𝑚𝑛 𝑅)
𝑒
(11)
ℎ𝑝
𝑘𝑚𝑛 2 − 𝑘 2 𝑠𝑖𝑛 𝜑𝑝 2
That is, the sound pressure in the exterior space outside the duct could also be decomposed into modes,
which would be expressed as follows
𝑝(ℎ𝑝 , 𝜃𝑝 , 𝜑𝑝 ) = ∑

𝑚+
𝑚=𝑚−

∑

𝑛0
𝑛=0

𝑝𝑚𝑛 (ℎ𝑝 , 𝜃𝑝 , 𝜑𝑝 )(12)

2.3 Mode Identification Based on The Inverse Technique
Assume that a microphone array is placed in the exterior space of the duct. The number ofmicrophones is denoted as 𝑁mic, the number of modes radiating from the opening is denoted as 𝑁mode.
According to Eq. (10) and Eq. (11), the complex acoustic pressure at the microphone array could be
represented in form of vector and matrix as[𝒑𝒎𝒊𝒄 ] = [𝑮𝒎𝒏−𝟏 ][𝑽𝒎𝒏 ](13)
Where [𝒑𝒎𝒊𝒄 ] is a 𝑁mic dimensional column vector, and [𝑽𝒎𝒏 ] is a 𝑁mode dimensional column
vector which contains modes coefficients 𝑉mn . [𝑮𝒎𝒏−𝟏 ] is the transfer function matrix which is
independent of the sources. When 𝑁mic > 𝑁mode, the least-squares solution [𝑽𝒎𝒏 ] is given by[𝑽𝒎𝒏 ] = [𝑮𝒎𝒏−𝟏 ]† [𝒑𝒎𝒊𝒄 ]
(14)
†
𝑯
−𝟏
𝑯
Here [𝑮𝒎𝒏−𝟏 ] = ([𝑮𝒎𝒏−𝟏 ] [𝑮𝒎𝒏−𝟏 ]) [𝑮𝒎𝒏−𝟏 ] is the Moore-Penrose inverse of [𝑮𝒎𝒏−𝟏 ] ; thesuperscript H represents complex conjugate transposition.
The calculation process of Eq. (14) is the identification process of radiation modes but isrepresented by the axial velocity modes. Once the modal coefficients of the normal velocity areascertained, the acoustic pressure and axial acoustic velocity on the cross-section of the duct open end
could be reconstructed, meanwhile, the acoustic field in the far-field out-duct could also be calculated.
Substituting Eq. (6) into Eq. (11) and Eq. (12), a similar formula would be derived as[𝒑𝒎𝒊𝒄 ] = [𝑮𝒎𝒏−𝟐 ][𝑨𝒎𝒏 ](15)This means amplitudes of incident modes could also be solved by the inverse techniques as Eq. (1 4).
Denote the element of the transfer function matrix [𝑮𝒎𝒏−𝟏 ] and [𝑮𝒎𝒏−𝟐 ] as 𝑔𝑚𝑛−1 and 𝑔𝑚𝑛−2
respectively. The relationship between 𝑔𝑚𝑛−1 and 𝑔𝑚𝑛−2 could be represented as
𝑔𝑚𝑛−2 = 𝑔𝑚𝑛−1 × 𝛺𝑚𝑛 (1 − 𝑅𝑚𝑛𝑛 )/𝜌𝑐
(16)

3. NUMERICAL SIMULATIONS
3.1

Generation of Simulation Data

Figure 2– Schematic diagram of a semi-infinite length flanged duct
As shown in Figure 2, a flanged duct with a duct radius 0.3𝑚 is considered, the radius of theflange is 1.2𝑚; the sound source in-duct is a monopole locates at (0.2,0, −1.2)𝑚 in coordinates of
Figure 2; the sound-field in-duct and in the near field are calculated by the finite element method(FEM), perfectly matched layers (PML) are adopted at the left of the duct and outer space of the flange
to simulate an infinite environment with no reflections; a spiral array of 56 microphones is settled
1.2𝑚 away from the flange plane whose center locates at (0,0,1.2)𝑚 in coordinates of Figure 2. The
analysis frequency of the sound source is from 200Hz to 1700Hz with an interval 100Hz.
3.2 Results of MI and Source Reconstruction
In each analysis frequency, the inverse technique in section 2.3 is performed to identify theradiation modes and incident modes. Then the axial acoustic velocity at the duct opening and the innersource could be reconstructed. Figure 3(a) demonstrates the result of radiation modes identification
at a frequency of 1700 Hz in which the number of propagable modes is 27. The reconstruction resultof axial acoustic velocity at the duct opening is shown in Figure 3(b).

Once the acoustic velocity along the z-direction is reconstructed, the sound pressure of an arbitrary
position out-duct could be recalculated according to Eq. (11) and Eq. (12). Assume an observation
point whose position is (0.3,0.3√3, 1.3) in the coordinate of Figure 2, the sound pressure at theobservation under the analysis frequency band is calculated and shown in Figure 4. Figure 4(a) shows
the calculation results of sound pressure level, while Figure 4(b) demonstrates the calculation
results of phase angle. Both the calculation results in Figure 4(a) and Figure 4(b) are comparedto the calculation results by FEM simulation and have been coincident well with FEM results,
which verify the effectiveness of the radiation modes identification method.

(a)(b)
Figure 3–Results of MI and axial velocity reconstruction at 1700 Hz(a) mode identification; (b) reconstruction at the duct open end(a)(b)
Figure 4–The acoustic pressure at the calculation position(a) Sound pressure level; (b) Phase angle
The generalized radiation impedance 𝑍𝑚𝑛𝑛 could be solved by Eq. (8), in this paper, the complex
Simpson's formula is adopted to calculate the integral numerically, after which the self -reflection
coefficient could also be solved by Eq. (7). Ultimately, the incident modes be identified through Eq.
(15) and Eq. (16) based on out-duct measurements.

(a)(b)
Figure 5–Calculation results of 𝑅𝑚𝑛𝑛 and incident modes identification at 1700 Hz(a) 𝑅𝑚𝑛𝑛 ; (b) Incident modes identification
Figure 5(a) demonstrates the calculation results of self-reflection coefficients for all propagativemodes under 1700 Hz. Based on the results of Figure 5(a), amplitudes of incident modes could
furtherly be calculated and are shown in Figure 5(b). At each analysis frequency, the amplitude ofradiation modes and incident modes are calculated through Eq. (15) and Eq. (18), after which the self reflective modes could be calculated by Eq. (8). Calculation results of incident modes and reflectivemodes at all analysis frequencies are shown in Figure 6(a) and Figure 6(b), respectively.

(a)(b)
Figure 6–Calculation results of modes identification at all analysis frequencies(a) Incident modes identification; (b) Reflective modes identification
As shown in section 2.1, incident modes are generated by the decomposition of the sound source,which means the inner sound source could also be reconstructed by the superposition of identified
incident modes. The sound pressure amplitude at the source plane in 1700 Hz is calculated by theidentified incident modes and shown in Figure 7(a), where the cross mark represents the position of
the source. Calculation results are represented in dB with reference sound pressure equal to
2 × 10−5 Pa. As shown in Figure 7(a), a clear source is reconstructed at the source plane, whose
position matches the actual source well. The deviation of radial position is mainly due to the limited
propagative modes and would decline with a higher frequency. Compared to the theoretical source
map is shown in Figure 7(b), the source strength calculated by the proposed method is 3 dB smaller.
The deviation is caused by the simplification of duct end reflection and is acceptable.

(a)(b)
Figure 7– Results of source reconstruction at 1700 Hz(a) Calculated source map; (b) Theoretical source map

4. EXPERIMENTAL VERIFICATION
4.1 Experimental Facilities and Setup
An experimental platform was established as shown in Figure 8, which is mainly composed of asteel circular duct whose inner radius was 0.3m , brackets to keep the duct level, two small
loudspeakers whose rated power was both 1W, and a wooden baffle at the opening whose parameters
were 2m×2m×0.01m. A spiral microphone array of 56 channels was mounted 1.25 m away from thebaffle. A microphone is placed about 2.5m away from the center of the duct opening and measured at
the same time with the microphone array. The center of the microphone array was coincident with theduct center for the convenience of calculation. Glass wool was filled at the left of the duct to simulate
a semi-infinite propagation environment. Some sound-absorbing sponges were placed on the groundto decrease the reflection of the ground. Chirp signals ranging from 200Hz~2300Hz were input inboth two loudspeakers. The input signals were generated by a signal generator and the frequency
sweep time was 2 seconds. The sampling frequency was 16384Hz, and the sampling time was 20

seconds.

Figure 8–Experimental setup
4.2

Experimental Results

(a)

(b)

(c)
(d)
Figure 9– Results of source reconstruction based on incident modes identification(a) 600 Hz; (b) 1000 Hz(c) 1700 Hz (d) 2300 Hz
The incident modes were first identified through Eq. (15). Based on the calculation results ofincident modes, the inner duct sources were reconstructed by calculating the sound pressure amplitude
at the source plane. Results of source reconstruction were partly shown in Figure 9 and represented
in dB, where the cross mark represents the actual position of sound sources. Both two sources could
be identified in Figure 9(a), Figure 9(c), and Figure 9(d), in which the positions of identified sources
have good coincidence with the actual positions. Only one source was reconstructed in Figure 9(b),this was mainly caused by the difference in frequency response between the two loudspeakers, whichcaused the strength at 1000 Hz of one source is much larger than another. The reconstruction resultsof inner-duct sound sources verified the effectiveness of the proposed method to identify incident
modes based on out-duct measurement.

5. CONCLUSIONS
A methodology based on out-duct measurement was proposed to identify duct modes. Comparedto mode identification methods based on in-duct measurements, it avoids the installation limitation.

Compared to current mode identification methods based on out-duct measurements, not only radiation
modes but also incident and reflective modes are identified, which could furtherly accomplish soundsources reconstruction with good resolution and accuracy.
The proposed method is theoretically limited to a situation without flow. Further research could be
conducted on flow duct situations with spatial flow field distribution.

ACKNOWLEDGEMENTS
This work is supported by the National Natural Science Foundation of China (Grant No.51835008).

REFERENCES
1. Bu H, Huang X, Zhang X. An overview of testing methods for aeroengine fan noise. Prog Aerosp Sci.
2021 Jul;124:100722.
2. Rienstra SW, Hirschberg A. An Introduction to Acoustics. Technische Universiteit Eindhoven; 2001.
3. Sijtsma P, Brouwer H. Deconvolution of azimuthal mode detection measurements. J Sound Vib. 2018
May;422:1–14.
4. Tapken U, Pardowitz B, Behn M. Radial mode analysis of fan broadband noise.23rd AIAA/CEAS
Aeroacoustics Conf; 5-9 June 2017; Denver, United states 2017.
5. Caldas LC. In-duct beamforming and mode detection using a circular microphone array for thecharacterisation of broadband aeroengine fan noise. Universidade de São Paulo; 2016.
6. Sijtsma P, Zillmann J. In-Duct and Far-Field Mode Detection Techniques. 13th AIAA/CEAS
Aeroacoustics Conf;21-23 May 2007; Rome, Italy 2007.
7. Enghardt L, Holewa A, Tapken U. Comparison of Different Analysis Techniques to Decompose a Broadband Ducted Sound Field in its Mode Constituents. 13th AIAA/CEAS Aeroacoustics Conf;21-23 May
2007; Rome, Italy 2007.
8. Sutliff DL, Dahl MD. Techniques for analyzing rotating rake mode measurements over passive treatment.
Int J Aeroacoustics. 2016;15(4–5):430–61.
9. Thomas R, Farassat F, Clark L, Gerhold C, Kelly J, Becker L. A mode detection method using theazimuthal directivity of a turbofan model. 5th AIAA/CEAS Aeroacoustics Conf Exhibit;10-12 May 1999;
Bellevue, WA, United states 1999.
10. Lan J, Premo J, Sutliff D. Inlet Mode Measurements with an Inflow Control Device Microphone Array.
8th AIAA/CEAS Aeroacoustics Conf Exhibit;17-19 June 2002; Breckenridge, CO, United states 2002.
11. Brusniak L, Underbrink J, Nesbitt E, Lynch D, Martinez M. Phased Array Measurements of Full-Scale
Engine Exhaust Noise. 13th AIAA/CEAS Aeroacoustics Conf;21-23 May 2007; Rome, Italy 2007.
12. Castres FO, Joseph PF. Mode detection in turbofan inlets from near field sensor arrays. J Acoust Soc Am.
2007 Feb;121(2):796–807.
13. Castres FO, Joseph PF. Experimental investigation of an inversion technique for the determination ofbroadband duct mode amplitudes by the use of near-field sensor arrays. J Acoust Soc Am. 2007
Aug;122(2):848–59.
14.Lewy S. Inverse method predicting spinning modes radiated by a ducted fan from free-field measurements.
J Acoust Soc Am. 2005 Feb;117(2):744–50.
15. Yu W, Huang X. Reconstruction of Aircraft Engine Noise Source using Beamforming and Compressive
Sensing. IEEE Access. 2018; 6:11716–26.
16. Zorumski WE. Generalized radiation impedances and reflection coefficients of circular and annular ducts.
J Acoust Soc Am. 1973 Dec;54(6):1667–73.
17. Shao W, Mechefske CK. Analyses of radiation impedances of finite cylindrical ducts. J Sound Vib. 2005
Aug;286(1–2):363–81.
18. Zhong D, Yang D, Zhu M. Improvement of sound source localization in a finite duct using beamforming
methods. Appl Acoust. 2016 Feb 1;103:37–46.
19. Norris AN, Sheng IC. Acoustic radiation from a circular pipe with an infinite flange. J Sound Vib. 1989
Nov 22;135(1):85–93.
20. Wendoloski JC, Fricke FR, McPhedran RC. Boundary Conditions Of A Flanged Cylindrical Pipe. J Sound
Vib. 1993 Mar 22;162(1):89–96.
21. McAlpine A, Daymond-King AP, Kempton AJ. Sound radiation from a flanged inclined duct. J Acoust
Soc Am. 2012 Dec;132(6):3637–46.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022 in Gyeongju, Korea

ABS-0195

Door vibration contribution separation of refrigerator equippedwith top-mounted compressor using transfer path analysis
Chisato SHIBAHARA1; Junji YOSHIDA1
1

Osaka Institute of Technology, Japan

ABSTRACTIn these years, home refrigerator with a compressor installed at the top area has been developed to improvethe usability. On the other hand, the rotational speed of the compressor is required to be operated at lower
speed for energy saving. However, this transition sometimes makes a large vibration and the vibration makes
rattling noise of inner bottles. Hence, effective vibration countermeasure is essential to compatible the energy
saving with the usability. In this study, we then focused on the technique for finding out factors increasing
the door vibration using transfer path analysis (TPA). Through the door vibration measurement, the vibration
was large at a low speed range. The factors were investigated using TPA. In the analysis, the door vibration
was separated into factors such as input force and transfer function of body in each direction. The result
showed the main factors of the vibration was the body vibration transmission characteristics along left-right
direction and the compressor force. As the verification, we increased stiffness of a rubber bush under thecompressor to change the high contributing input force. As the result, the vibration at the low speed range
could be changed and applicability of TPA to the refrigerator was confirmed.
Keywords: Refrigerator, Compressor, TPA, Force identification, Contribution

1. INTRODUCTION
In recent years, due to increase of the global warming concern, home appliances are required to
improve the electricity consumption performance. Accordingly, the rotational speed of the compressor
had better to be set at low speed to reduce the electrical power consumption. On the other hand,
refrigerators are required not only the low electricity consumption performance but also thefunctionality. Hence, some refrigerators are equipped with the compressors at the upper area, where
the users are hard to touch. However, in case a compressor is installed on the upper part, the vibrationof the refrigerator generally becomes large because of the upper input point.
In this study, we then measured and analyzed the vibration of the door at various compressor
rotational speed to find out the main factor increasing the vibration through transfer path analysis(TPA) and attempted the vibration reduction considering the contribution.

2. VIBRATION MEASUREMENT OF REFRIGERATOR DOORIn this section, we measured vibration of a refrigerator in which the compressor was installed at
top area. As the measurement point, door upper point was selected because the door vibration
sometime shakes bottles and makes noise. In the measurement, the rotational speed of the upper
mounted compressor was changed from low to high speed to find the large vibration speed and thefrequency band.
2.1 Experimental Conditions
Large size refrigerator (1.8 m height and 100 kg weight) in which the compressor was installed at
top side was used for the measurement as shown in Figure 1. For the measurement, an acceleration
sensor was attached to the top of the door to measure vibrations along front -back, left-right, and updown directions, respectively. In the measurement, the door was opened and the rotational speed of
the compressor was changed from 16 to 25 rps in each 1 rps step. In each rotational speed, the vibration
1

junji.yoshida@oit.ac.jp

was recorded for about 10 s at each speed.

Figure 1 – Vibration measurement condition and acceleration attachment position
2.2 Vibration Displacement of the Door
Figure 2 shows the first order vibration displacement at each compressor rotational speed
calculated from the obtained vibration acceleration.

Figure 2 – First order door vibration displacement in each direction
Horizontal and vertical axes show the rotational speed and the displacement, respectively. Blue,
red and green solid lines show the vibration displacement along front-back, left-right and up-down
directions. As a result, it was found that the vibration was observed to be large at around 18 rps band
along front-back direction. Accordingly, the factors causing the large vibration was investigated usingtransfer path analysis in the following sections.

3. TRANSFER PATH ANALYSIS FOR ANALYSING THE CONTRIBUTION OF
DOOR VIBRATION FROM THE COMPRESSOR (1-7)
In this section, transfer path analysis (TPA) was applied to the refrigerator to find the main factor
increasing the door vibration. In the TPA, the influence of output force of compressor, transfer
characteristic of compressor rubber bush, transfer characteristic of refrigerator body in each direction
on the door vibration are analyzed.
3.1

Estimation of the Output Force from the Compressor and Input Force to the

Refrigerator Body Using Inverse Matrix Method
Here, we estimated output force of the compressor and transferred input force to the refrigerator
body in each direction by using inverse matrix method. We then conducted a hammering test to obtainpoint inertance of the compressor in three directions and transfer functions from compressor side to
the refrigerator body where the compressor was installed to estimate the output force and obtaintransfer characteristics of compressor rubber bush. As the measurement point of compressor, thecompressor was supported by four rubber bushes as shown in Figure 3. Hence, measuring point
inertance and transfer functions at all rubber bush attachment points are considered to be the best. On
the other hand, there was not enough space for the measurement in general. Then, the representative
single point at both compressor side and refrigerate side were selected for the measurement as shownin Figure 3 in this study.

Figure 3 – Vibration position of compressor and refrigerator body
Figure 4 shows the point inertance of the compressor in each direction. By using this point
inertance and the vibration at the point in the operational condition, the output force from thecompressor in each direction was estimated.

Figure 4 – Point inertance of the compressor in each direction
Horizontal and vertical axes show the frequency and the transfer function, respectively. Blue, redand green solid lines show the transfer function along front-back, left-right and up-down directions.
In the measurement of compressor vibration at the operational condition, the identical point of thecompressor was used and the acceleration along each direction was measured at each rotational speed.
After then, the acceleration signals were divided by the point inertance to estimate the output force
from the compressor by using inverse matrix of the point inertance. Figure 5 shows the calculated
output force from the compressor in each direction.

Figure 5 – The calculated output force from the compressor in each direction

Horizontal and vertical axes show the rotational speed and the force, respectively. Blue, red andgreen solid lines show the calculated force along front-back, left-right and up-down directions. As
shown in Figure 5, the output force was found to be large at low speed condition. Subsequently, therefrigerator body vibration at the compressor attachment point was calculated in each direction bymultiplying the estimated compressor output force with the transfer function. Figure 6 shows theestimated vibration displacement of the refrigerator body in each direction.

Figure 6 – Estimated vibration displacement at around the input force point to the refrigerator bodyHere, a hammering test was again conducted to estimate the input force of the refrigerator body inthe same manner as the compressor output force estimation except for the measurement point. In th iscase, the point inertance of the input force point of the refrigerator body was measured in the three
directions as shown in Figure 7 (a). In addition, the vibration at the refrigerator body along the three
directions was also measured at the operational condition. And the input force to the refrigerator bodywas estimated by dividing the operational vibration by the point inertance. Figure 7 (b) shows theestimated input force level at the refrigerator body side.

(b) Input force to the refrigerator body(a) Point inertance
Figure 7 – Measured point inertance and calculated input force to the refrigerator body
3.2 Vibration Transfer Characteristics of Refrigerator Body
In this part, we evaluated the vibration transfer characteristics of the refrigerator. Hammering test
was then carried out again between the input force point of the refrigerator body to the door as
indicated in Figure 1. Figure 8 shows the transfer function from the input point of the refrigerator
body along three directions to the door top side along front-back direction.

Figure 8 – Transfer function from the input force point of the refrigerator body to the doorAs the result, it was found that the transfer function from the body along left-right direction wasparticularly large. This indicates that the refrigerator body had natural frequency at around 18 Hz andthe door vibration has a possibility to be increased by the transfer function at the operational condition .
Here, we conducted vibration mode analysis to the refrigerator body (8). Figure 9 shows thevibration mode of the refrigerator body at 18 Hz where the transfer function had the peak.

Figure 9 – Vibration mode at 18 Hz
The result shows that the vibration mode was twist mode at the upper part of the refrigerator body.
And this mode increased the transfer function from the input force along left -side to the door.
3.3 Contribution Analysis to the Refrigerator Door Vibration
The contribution from input point of the refrigerator body to the door can be obtained in eachdirection by multiplying the estimated input force to the refrigerator body obtained in Sec. 3.1 and thetransfer function of the refrigerator obtained in Sec. 3.2. In addition, the total door vibration along thetarget direction can be obtained by summarizing the contributions in three directions. Figure 10 shows
the contribution elements of the door vibration along front-back direction. In the elements, compressor
output force, input force to refrigerator body, transfer function and contribution w ere illustrated ineach direction. In addition, force transmission rate from the compressor to the refrigerator body wascalculated by using the input and output force in each direction.

Figure 10 – Contribution elements to door vibration along front-back direction
From the upper contribution element, the contribution along left-right direction was observed to
be dominant to the overall door vibration at the target rotational speed of 18 rps. And the large leftright contribution was found to be made by the transfer function (twisted mode) along left -right
direction. In addition, the transferred force along the direction was also large in comparison to theother directions. This indicates that the door vibration can be decreased effectively by the transfer
function reduction at the 18 Hz. However, a major structural modification is required to reduce thispeak because the vibration mode is twist vibration behavior at the upper part of the refrigerator body.
Then, we attempted to decrease the transfer force to the refrigerator body from the compressor in t henext section.

4. VIBRATION REDUCTION CONSIDERING THE CONTRIBUTION
For the reduction of the transferred force to the body at the target speed, we attempted to ex change
the compressor rubber mount bush stiffness. In the consideration of the mount rubber stiffness, theoutput force along left-right direction became larger at lower frequency as shown in Figure 11. And
the force transmission rate had peak at around 21 Hz. Therefore, the force transmission rate and thecontribution were expected to decrease at the target frequency of 18 Hz if the peak frequency becomes
higher by the increase of the rubber bush stiffness as shown in arrow in Figure 11.

Figure 11 – Left-right contribution elements to door front-back vibration(Pick up elements from Figure 10)
4.1 Rubber Bush Stiffness Increase
Here, we actually increased the rubber bush stiffness around 2 times harder and the force
transmission rate was again calculated through transfer path analysis. Figure 12 shows the comparison
before and after modification of rubber bush.

(a) Front - Back(b) Left - Right(c) Up - Down
Figure 12 – Force transmission rate comparison before and after modification
Horizontal and vertical axes show the rotational speed and the force transmission rate, respectively.
Solid and dotted curves show the force transmission rate of original (low stiffness) and modified (highstiffness) rubber bush. As shown in the upper figures, the force transmission rate along the highcontributing left-right direction was decreased at the target frequency band. Then, the door vibration
at the operational condition was measured before and after modification. Figure 13 shows the door
vibration comparison before and after modification in each direction.

(a) Front - Back(c) Up - Down(b) Left - Right
Figure 13 – Comparison of door vibration displacement due to rubber bush hardness change

Horizontal and vertical axes show the rotational speed and the displacement, respectively. Solidand dotted lines show the original and modified vibration displacement. As shown in this figure, thedoor vibration along front-back direction, where the vibration was large at the original condition,
could be decreased at the target (large) rotational speed by increasing the stiffness. In addition, themodification was also confirmed not to increase significantly to the door vibration along the otherdirections.

5. CONCLUSIONSIn this study, we attempted to analyze the factor of the refrigerator’s door vibration in detail by
employing TPA for the effective countermeasure. Through the door vibration measurement, thevibration was large at a low speed range. In the application of TPA, the door vibration was separated
into compressor output force, force transmission rate by rubber mount, transfer function of body ineach direction. The result showed the main factors of the vibration was the body vibration transmission
characteristics along left-right direction and the compressor force. Then, the stiffness of a rubber bushunder the compressor was increased to decrease the input force at the large vibration range. As theresult, the vibration at the low speed range was decreased effectively and applicability of TPA to therefrigerator was confirmed.

ACKNOWLEDGEMENTS
We would like to express our gratitude to Panasonic Corporation Appliances Co., Ltd. for providing
the refrigerator and each part in advancing this study.

REFERENCES
1. Van der Auweraer H, Wyckaert W, Van Der Linden P, Noise and Vibration Transfer Path Analysis, Lecture
series - van Kareman Institute for fluid dynamics, van Kareman Institute, Rhode-Saint-Genèse, Belgium
1979.
2. Van der Auweraer H, Mas P, Dom S, Vecchio A, Janssens K, Van de Ponseele P, Transfer Path Analysis
in the Critical Path of Vehicle Refinement: The Role of Fast, Hybrid and Operational Path Analysis, SAE
2007-01-2352, Proceedings SEA Noise and Vibration Conference, St.Charles (IL), USA 2007. p. 15-17.
3. Noumura K, Yoshida J, Method of Transfer Path Analysis for Vehicle Interior Sound with No Exci-tation
Experiment, FISITA2006 Proceedings F2006D183, Japan 2006. p. 1-10.
4. Yoshida J, Kawagoe N, Kawamura T, Improvement of Bicycle Riding Comfort by Reduction of Seat
Vibration, Journal of System Design and Dynamics, Japan 7 (3), 2013. p. 293-303.
5. Yamashita R, Factor analysis of noise for drum type washing machine using operational TPA, Inter-Noise
2015, San Francisco, USA 2015. p. 1-10.
6. Bermayer W, Brandl F, Hoeldrich R, Sontacci A, Brandl S, Fliesser W, Sound Engineering based on
Source Contributions and Transfer Path Results, JSAE Paper 318/2007/5399, 2007. p. 15–17.
7. Yoshida J, Tanaka K, Contribution analysis of vibration mode utilizing operational TPA, Mechani-cal
Engineering Journal 3 (1), 2016. p. 1–15.
8. Inoue Y, et al, Shindo no Kangaekata Toraekata, (How to think and grasp vibration), Ohmsha, 1998.
p. 139–146. (In Japanese)

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022, in Gyeongju, Korea

ABS-0292

Application of operational TPA to handle vibration of electric
powered lawn mower
Senta UEGAKI1; Junji YOSHIDA1
1

Osaka Institute of Technology, Japan

ABSTRACTIn this study, applicability of operational transfer path analysis (OTPA) to noise and vibration, in which theoperational condition is constant, was considered. As an instance, handle vibration of an electric poweredlawn mower was investigated for understanding the contribution from motor vibration along each direction
to the handle. OTPA was originally proposed for contribution separation of vehicle noise and vibration, andthe contribution separation accuracy is reported to be improved by using signals measured at various
conditions. However, the lawn mower was operated at a constant speed and could not be operated undervarious conditions. In the vibration measurement, vibration signals at four motor attachment points along
front-back, left-right and up-down directions were measured as the reference signals and the handle frontback vibration was recorded as the response signal. In the transfer function calculation from the obtained
signals, a noise reduction method was applied to prevent accuracy deterioration by high correlation among
signals. As the result, the contribution from the motor vertical vibration was indicated to be the main
contributor. As the verification, motor weight was increased to reduce the motor vertical vibration and the
response vibration was observed to decrease and the accuracy of the method was confirmed.
Keywords: Operational transfer path analysis, Lawn mower, Noise reduction

1. BACKGROUND AND PURPOSE
Lawnmower is a kind of machine used to cut grass in gardens and other areas. They are used ataround the world due to easy operation and ability to cut wide grass area. In recent years, electricmotor powered lawnmowers have become popular for the low environmental load as shown in
Figure 1. However, even though electric lawnmowers have good environmental performance, if thevibration is large, this increases the risk of hand-arm vibration syndrome. Therefore, small vibrationis still very important performance for operator protection(1-5).

Figure 1 – Lawnmower and the measurement vibration direction

1

junji.yoshida@oit.ac.jp

On the other hand, vibration reduction often requires structural modification or additional weight.
Hence, it is important to understand the main factors increasing the handle vibration in order toeffectively reduce vibration. In this study, we then applied operational transfer path analysis (OTPA)
as an analysis method for the contribution separation. This method has been proposed for contribution
separation technique of vehicle interior noise using only operational signals(6) and reported to benecessary to measure various operational data set to obtain accurate contribution(7). However, theelectric lawnmower employed in this study is generally used at a constant rotational speed and wecannot obtain signals under various operating conditions. Then, we consider how to obtain accurate
contribution by applying OTPA to a lawnmower in several ways through the actual vibration reduction
verification.

2. LAWNMOWER VIBRATION MEASUREMENT AND OTPA

160
150
140
130
120
110
100
90
80

Front-Back
Right-Left
Up-Down

0

20

40
60
Frequency [Hz]

80

100

VAL [dB]

VAL [dB]

2.1 Measurement of Lawnmower Vibration at Operational Condition
Handle vibration of the lawnmower was firstly measured at the constant motor rotational speed at2800 rpm to identify the vibration reduction frequency band and grip position as shown in Figure 1.
In the measurement, the vibration at the upper handle position and lower handle position along frontback, left-right, and up-down directions were measured for 10 s. The vibration acceleration levels at
the lower and upper handle positions are shown in Figure 2 (a) and (b), respectively.
160
150
140
130
120
110
100
90
80

Front-Back
Right-Left
Up-Down

0

20

40
60
Frequency [Hz]

80

100

(a) Vibration at the lower handle position(b) Vibration at the upper handle position
Figure 2 – Vibration acceleration levels at handle along front-back, left-right, and up-down direction
Vertical and horizontal axes in both figures indicate the vibration acceleration level (VAL) and thefrequency. Red, blue, and green bars indicate the VALs along front-back, left-right, and up-down
directions, respectively. These results show that the vibration at upper point along front-back direction
was the largest, and the vibration at 47 Hz, which is the frequency of the first -order vibration of the
motor indicated by the arrow in Figure 2, was the largest. Accordingly, we attempt to identify the main
contribution of the large vibration by applying OTPA to reduce the front-back handle vibration at
upper point at 47 Hz.
Procedure of OTPA
In this section, procedure of OTPA is described as follows.
In OTPA, the contribution of each reference point to the response point is obtained by multiplying
each reference signal with the transfer function. The transfer function in this method is calculatedusing principal component regression method. The calculation procedure is shown in Figure 3(8-11).
2.2

Reference
signalPrincipal
component

Ain1

T1

Ain2

T2

Ain3

T3

Ain4

T4

Response
signalAout

Multiple regression
Principal component
analysisanalysis
Figure 3 – Calculation image of transfer function of original OTPA using principal component regression method

Firstly, principal component analysis is applied to the reference signal matrix [ Ain] by singularvalue decomposition to remove correlation among reference signals as shown in Eqs. (1) and (2). The
calculated uncorrelated signals are principal component [T] .

Ain  = U S V T
T  = Ain V  = U S 

(1)
(2)

After eliminating the noise component having very low level, multiple regression analysis isapplied between the remained (signal) principal components [T] and the response signal [Aout] toobtain the influence [B] of each principal component to the response signal as shown in Eqs. (3) and(4) .

Aout  = T B(3)

B = (T T T ) T T Aout 

(4)

−1

Transfer function from reference signal to response signal [H] is calculated by multiplying thecoefficient [V], that connects reference signal to principal component, and the regression coefficient[B], that connects principal component and response signal as shown in Eq. (5)H  = V ( T T T ) T T Aout −1(5)Finally, the reference point contribution and principal component contribution to the response
point are calculated as shown in Eqs. (6) and (7), respectively.

Acont  = Ain H (6)

Tcont  = T B (7)This is the outline for obtaining contributions by OTPA.

3. VIBRATION CONTRIBUTION SEPARATION TO HANDLE FROM MOTOR
3.1
Application of OTPA to Lawnmower
In OTPA application, the main input force is the inertia from rotating motor and blade. Hence, we
calculated the contribution of input forces along three directions to the handle vibration by using themethod. For the application, the handle vibration at the upper point along the front -back direction was
set as the response point as shown in Figure 4 (a). About the reference point, acceleration signals ataround force input point are generally used to estimate the high contributing part. In this product,
there are four motor attachment point as shown in Figure 4 (b). However, measuring all attachment
point vibration along three directions (total 12 points) was time-consuming and we supposed thevibrations had very high correlation each other at the low frequency band due to close point as shownin Figure 4 (b). Then, in the first OTPA, we set a single attachment point (No.1 in Figure 4 (b)) as thereference point. Accordingly, the vibration of the response point and the reference point along threedirections were measured simultaneously for obtaining the contribution from the motor vibration
along each direction to the handle vibration.

Up-Down
Back

Front-Back

4

1

3

2

Left

Right
Front

(a) Handle as response point(b) Motor attachment point as reference points
Figure 4 – Measurement points for OTPA

VAL [dB]In this experiment, the lawnmower was operated at a constant motor speed of 2800 rpm for 10 sunder same conditions. Figure 5 shows the contributions of each motor vibrations to the handle
vibration along front-back direction.
160
150
140
130
120
110
100
90
80

Sum
Front-Back
Right-Left
Up-Down

40

45
50
Frequency [Hz]55

Figure 5 – Calculated contribution to handle vibration by OTPA
Red, blue, and green bars show the contribution of the motor contribution along front-back, leftright and up-down directions, respectively. In addition, black bar shows summarized contributions.
This result shows that the contributions from motor vibration along front-back and left-right directions
are larger than the up-down direction at the target frequency of 47 Hz. In the next section, we verified
the first OTPA application through a simple structure modification.
3.2 Verification of Contribution Separation Accuracy
In the verification test, 17 kg weight was placed on the top of the motor to decrease the reference
point vibration along vertical direction where the contribution was lowest. If the calculated
contribution result was correct, the response point vibration will not decrease significantly by themodification because only low contribution reference point was measured. In this condition, thevibration at reference and response points were again measured. Figure 6 shows the comparison before
and after modification at each point.

Original
Modified
VAL [dB]

VAL [dB]

160
150
140
130
120
110
100
90
80
40

45
50
Frequency [Hz]

160
150
140
130
120
110
100
90
80
40

55

Original
Modified
VAL [dB]

VAL [dB](a) Vibration at the response point
160
150
140
130
120
110
100
90
80
40

45
50
Frequency [Hz]55

Original
Modified

160
150
140
130
120
110
100
90
80

45
50
Frequency [Hz]55(b) Vibration at the attachment point
along front-back direction
Original
Modified

40

45
50
Frequency [Hz]55(c) Vibration at the attachment point(d) Vibration at the attachment point
along left-right direction
along up-down direction
Figure 6 – Vibration acceleration comparison before and after modification
Black and red bars show the vibration acceleration levels before and after modification. As theresult, not only the vibration of reference point along up-down direction (low contribution point) was
decreased largely, but the response point (handle vibration) was also decreased so much. This meansthat the contribution of the up-down vibration at the motor attachment point had high contribution in
actual and OTPA in the first case did not separate contribution accurately. As the reason why thecontribution accuracy was low in this OTPA application is considered as follows. The first possibility
is the inadequate operational condition for obtaining robust transfer function. Another possibility is
the number of reference point. In this application, we set single reference point for carrying out theanalysis easily but there were four input points in actual. This reduction may decrease the accuracy.
Then, in the next section, we discuss how to improve the accuracy by modifying the experimental andanalytical conditions.

4. ACCURACY IMPROVEMENT OF CONTRIBUTION SEPARATION
4.1 Improvement by Increasing Number of Reference Points
In the previous section, OTPA using single reference point was not able to calculate accurate
contribution. And as the reason, constant operational condition and the number of reference points
were considered. However, the operational condition of this product was not changed easily. Therefore,
in this section, the number of reference points was increased from single point to the original fourattachment points along three directions (total 12 points), and the contribution was calculated again
by OTPA as second application. In addition, to identify the high contribution direction of motor
vibration, vibration contributions in the same direction at four reference points were summarized
considering the phase. Figure 7 shows the contribution of motor vibration along each direction to thehandle vibration along front-back direction.

Sum
Front-Back
Right-Left
Up-Down

40

45
50
Frequency [Hz]VAL [dB]

VAL [dB]

160
150
140
130
120
110
100
90
80

160
150
140
130
120
110
100
90
80

55

Sum
Front-Back
Right-Left
Up-Down

40

45
50
Frequency [Hz]55(a) OTPA result in first application(b) OTPA result in second application(Single reference point)(Original four reference points)
Figure 7 – Contribution separation comparison according to the reference point number
Red, blue, and green bars show the contribution of the motor contribution along front-back, leftright and up-down directions, respectively. In addition, black bar shows the summarized contributions.
By increasing the number of reference points to the original four attachment points, the contributionof motor vibration along up-down direction was observed to be the highest as shown in in Figure 7(b). However, the contribution difference was very small. Furthermore, all calculated contributions
were much higher than the summarized contribution. This result makes us misunderstand that allcontributions are high and reduction of all contributions are necessary. This inconvenient situation is
known to occur when the correlation between reference point signals is high as same as
multicollinearity in multiple regression analysis(12, 13). In this case, the influence of measurement
error is increased very largely in the calculated contribution. We then consider how to remove thenoise component for improvement the accuracy more in the following section.
4.2 Noise Component Reduction in Contribution Separation
In the calculation of transfer function of OTPA, Eq. (5) was employed.

H  = V ( T T T ) T T Aout −1(5)Contribution Ratio [%]In Eq. (5), if the correlation between reference points is high, the inverse matrix may become
unstable (multicollinearity), and the transfer function increase due to the influence of very small noise
components in the measurement signal. Therefore, appropriate noise r eduction can be considered as
one of the method to improve the contribution accuracy
Figure 8 shows the cumulative principal component (PC) contribution ratio of the reference point
signal at 47 Hz, which was used for contribution separation in the previou s section. The contribution
ratio shows the summarized percentage of each PC contribution ([T] matrix in Eq. (5)) from largest
PC1 to smallest PC12. Noting that the number of PC equals to the number of reference point
theoretically.
100
98
96
94

92
90
0

1

2

3

4

5

6

7

8

9 10 11 12

Principal Component

Figure 8 – Cumulative contribution ratio of principal components at 47 Hz

VAL [dB]

Vertical and horizontal axes indicate the cumulative contribution ratio, and the PC number. And
black and red bars show the cumulative contribution at each PC number and only PC1 contribution ,respectively. The results show that the PC 1 contribution was dominant and the other PC did not havesignificant contribution at the target frequency. This suggests that reference point signals had very
high correlation and the noise influence may be increased in the reference point contribution.
We then reduce the principal components which contribution was 1% or less at all frequency bandas the noise components and contribution to the response point was again calculated as shown in
Figure 9.
160
Sum
150
Front-Back140
Right-Left
130
Up-Down
120
110
100
90
80
40
45
50
55
Frequency [Hz]Figure 9 – Contribution separation results with noise reduction
Red, blue, and green bars show the contribution of the motor contribution along front-back, leftright and up-down directions, respectively. In addition, black bar shows the summarized contributions.
As the result, the contribution along up-down direction was calculated as much larger than the othercontribution. This tendency is very similar with the verification test in Sec. 3.2. In addition, the highest
contribution is almost same as the summarized contribution. From these results, the contribution
separation accuracy seems to be improved so much by the noise reduction. Finally, the contributionof the electric lawnmower which operational condition is constant could be obtained very accurately
by increasing the number of reference points to the original four attachment points and removing thenoise component.

5. CONCLUSION
In this study, we considered a contribution separation technique to identify the high contributing
vibration direction of a motor attachment point by applying OTPA for the effective vibration reductionof an electric lawnmower. In the vibration measurement at the handle, the vibration at the upper point
along front-back direction at the motor first-order frequency was found to be large. Then, we applied
OTPA to separate contribution easily by setting single motor attachment point along three directions
as a reference point and the handle vibration along front-back as the response point. As the result, theobtained contribution was observed to be incorrect through a simple verification test. We thenincreased the number of reference points from one to original four motor attachment points andapplied OTPA again. The result indicated a certain accuracy improvement but not enough.
Subsequently, we applied noise component reduction procedure in addition to the reference point
increase. As the result, the contribution accuracy was increased so much and the contribution of theelectric lawnmower which operational condition was constant could be obtained correctly byincreasing the number of reference points to the original four attachment points and removing thenoise component.

REFERENCES
1. Lindell H, Lawnmower with Vibration Reduction, US 2006, US20060156705A1.
2. Bacon J R, Mower with Shock Mounts,US 1997.US4015407A.
3. Ryan A C, Balanced Reciprocating, Brush Cutting, Mower Attachment, US 2001.US6314707B1.
4. Foster R, Independent four Wheel Vibration Damping System for Riding Mowers, US 2001,
US7708292B2.

5. Tint P, Tarmas G, Koppel T, Reinhold K, Kalle S, Vibration and Noise caused by Lawn Maintenance
Machines in Association with Risk to Health, Agronomy Research Biosystem, Estonian 2012.10(1). p.
251-260.
6. Noumura K, Theory and example of operational transfer path analysis, Acoustical Society of Japan, Japan
2011. 67(4). p.163-168.
7. Yoshida J, Noumura K, Basic Consideration for Improvement of TPA Method Employing Running Data,
Japan Society of Mechanical Engineers, Japan 2011. 77(777). p. 1720-1728.
8. Van der Auweraer H, Wyckaert W, Van Der Linden P, Noise and Vibration Transfer Path Analysis, Lecture
series - van Kareman Institute for fluid dynamics, van Kareman Institute, Rhode-Saint-Genèse, Belgium
1979.
9. Van der Auweraer H, Mas P, Dom S, Vecchio A, Janssens K, Van de Ponseele P, Transfer Path Analysis
in the Critical Path of Vehicle Refinement: The Role of Fast， Hybrid and Operational Path Analysis,
Proc. SEA Noise and Vibration Conference, St.Charles (IL), USA 2007.p. 15-17.
10. Noumura K, Yoshida J, Method of Transfer Path Analysis for Vehicle Interior Sound with No Excitation
Experiment, Proc. FISITA2006 World Automotive Congress, Japan 2006.p. 1-10.11.
Noumura K,
Theory and example of operational transfer path analysis, Acoustical Society of Japan, Japan 2011. 67(4).
p.163-168.
11. Yoshida J, Kawagoe N, Kawamura T, Improvement of Bicycle Riding Comfort by Reduction of Seat
Vibration, Journal of System Design and Dynamics, Japan 2013.7(3).p. 293-303.
12. Yutaka S, Suppressor Variable and Multi-collinearity in Multiple Linear Regression Analysis-Correlationand Explanation Power-, Bulletin of the Shizuoka Institute of Science and Technology, Japan 2010.18.p.
179-184.
13. Wakui Y, Wakui S, Illustrated Regression Analysis - A Practical Prediction Method for Revealing
Complex Statistical Data, Nippon Jitsugyo Publishing Co., Ltd. Japan 2002.p. 104-107. 147-211.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022, in Gyeongju, Korea
Extended abstractABS-0670

An excitation technique in a vibro-acoustic inverse method
for the efficient radiation of panel speaker
Jung-Han Woo1; Dongjoon Kim1; Yun-Ho Seo1; Pyung-Sik Ma1
1

Korea Institute of Machinery and Materials, Korea

In the case of directly excitation of a plate to radiate sound using a thin panel structure, a sound quality
deteriorates due to the influence of many resonance modes of the rectangular plate, and sound power andradiation efficiency are also reduced. To overcome this problem, an indirect vibration-acoustic rendering
technique was applied to control the vibration field composed of multiple virtual speakers and baffles using
an array of vibrators placed at the corners1. The direct rendering method using the relationship between theforce of the vibrator array and the sound pressure at the measurement point in space was also applied without
any restrictions on the vibration field of the plate surface2. Normal forces excited from actuators array on theplate are only considered to generate the desired sound field in previous studies3. However, not only normal
forces but also moments are provoked by adjacent actuators4 and they are used in the inverse vibro-acoustic
problem for an efficient input power. Angular velocity can be approximated by using finite difference method
as show in Fig. 1 and moment mobility can be indirectly calculated by using power equation and reciprocity
with force mobilities5. To compare the performance characteristics in generating the rendered sound radiationby using the two rendering methods as explained above, a thin, flat, homogeneous, isotropic, and simply
supported aluminum plate is used for the numerical simulation. The dimensions of the plate are 0.7 m inlength, 0.4 m in width, and 1 mm in thickness. The material properties of the panel are as follows: density,
2700 kg/m3; Young's modulus, 60 GPa; Poisson's ratio, 0.33; loss factor, 0.01. The boundary condition in thepractical situation is usually not easy to precisely define; in this work, a simply supported boundary conditionis assumed. In the edges of the plate, 196 actuators are located, which are 30-mm apart from the edge andwith 10-mm uniform spacing to avoid the spatial aliasing. Inside the panel, a total of 2079 vibration
observation points are distributed uniformly with 10 mm spacing. The generated radiation field within thefrequency range of interest is investigated by setting LW= 98 dB or Lp= 90 dB as the desired sound strength.
The required power with the suggested method is reduced as 5 dB. Although the calculated energy with has
no difference with employing regularization technique, actuators positions which needs a high power andcan be regarded as important point are moved adequately from the corner to the center. It is observed that the
input power of the array actuators used in both indirect and direct rendering methods can be reasonably
reduced by adopting suggested technique to fulfill the generation of an omnidirectional acoustic field.
Keywords: Vibro-acoustic inverse problem, Panel speaker, Vibration field control, Actuator array

Figure 1. Schematic diagram of the forces and moments for the finite difference method.
1

jhwoo@kimm.re.kr

140
120
100
80
60
40
IRM: F
DRM: F
IRM: F+M
DRM: F+M

20
0

0

200

400

600

800

1000

Frequency (Hz)Figure 2. Calculated input energies of the actuators for generating a uniform hemispherical sound radiation.

Figure 3. Reconstruction error of the acoustic field and required input energy of the actuators for generatinga uniform hemispherical sound radiation.

ACKNOWLEDGEMENTS
This work was supported by the research project ‘‘Artificial Intelligence-Based Diagnosis andDamage Management Technology’’ at Korea Institute of Machinery and Materials (KIMM) underGrant NK238B.

REFERENCES
1. Woo J.-H., Ih J.-G. Generation of virtual speaker and baffle on a thin plate controlled by an actuator arrayat the boundary. IEEE/ASME Trans. Mechatronics. 2019;24:1197-1207.
2. Woo J.-H., Ih J.-G., Park Y. Comparison of two vibro-acoustic inverse methods to radiate a uniform soundfield from a plate. J. Sound Vib. 2019;458:445-457.
3. Woo J.-H., Ih J.-G. Generation of the directional sound field from a panel speaker by suing the rendering
method. Proc.of the Spring Meeting of the ASK; 6-7 May 2021; Jeju, Korea.
4. Sanderson, Michael A., and C. R. Fredö. Direct measurement of moment mobility: Part I: a theoretical
study. J. Sound Vib. 1995;179.4:669-684.
5. Kim, D., Park N.-C. Calculation and reduction of sound radiation from a thin plate structure excited bycomplex inputs. J. Sound Vib. 2020;484:115517.

PROCEEDINGS of the24th International Congress on Acoustics
October 24 to 28, 2022, in Gyeongju, Korea
Extended abstractABS-0722

Production line diagnostics using a microphone array
and artificial intelligence processing
Abdelhakim DARRAZ1,2; Jean-Hugh THOMAS1; Anthony LARCHER3;
Jérôme ANTONI4; Pierre MOLLON2; Gabriel KIRIE2;
1

Laboratoire d’Acoustique de l’Université du Mans, Institut d’Acoustique-Graduate School (IA-GS), France
2
3

VALEO Equipements Electriques Moteur, France

Laboratoire d’Informatique de l’Université du Mans, LIUM EA4023, France
4

Laboratoire Vibration Acoustique, France

ABSTRACTThe objective of the study is to lead to a diagnosis on a production line. Specifically, the aim is to classify
automotive starters by recording the sound emitted by these systems when they are powered. A microphonearray is used for this purpose in the vicinity of the starter. Each starter must be classified into one of two
categories: healthy or faulty. The difficulty of the problem is the environment of the production line, whichcontains both background noise and unwanted noise from any noisy operation in the vicinity of the area.
Acoustic imaging provides spatial information on the origin of the sound. This study is interested in therelevance of the descriptors extracted from the acquired signals, which will then be presented to a deep neuralnetwork for decision making. The descriptors are based either on acoustic mappings or on time-frequency
representations of an acoustic signal which is the resultant of beamforming focused towards a precise pointof the starter. The approach is applied to a database collected on the production line. The results obtained on
this database are presented and compared as a function of the input parameters.
Keywords: Acoustic diagnosis, Microphone array, Deep neural networks

1. INTRODUCTION
The object manufactured and tested on the production line is an automotive starter. A starter is a
DC electric motor used to run an internal combustion engine until the engine's self-combustion takes
over. The aim is to be able to distinguish between a healthy starter that can be fitted to a vehicle anda starter with a manufacturing defect, which must be removed from the line. The starter is thus fed
into the line leading to an engine speed of around 7200 rpm. During this operation, the starter emits asound that allows an expert operator to recognize its good condition. In the presence of a defect, thesound perceived is different according to the expert. The idea of the study is therefore to set up asystem for acquiring acoustic signals on the production line. Automatic signal processing should thenallow the acoustic diagnosis to be made, i.e. either accepting the starter as compliant or rejecting it asdefective. One problem encountered is the presence of unwanted and disruptive noise on the line,
which is why the system implemented will also use directional sound monitoring.

2. MEASUREMENT SYSTEM
A microphone array was built for this study. It consists of 16 digital MEMS arranged in two parallel
planes very close to each other. The arrangement of the microphones allows a satisfactory directivity
to be obtained in the [1000, 6000] Hz frequency range. The antenna is placed on the production line
35 cm from the starter under test. A Mini DSP data acquisition card was used to record the data andtransmit it via a USB cable to a computer for processing. Figure 1 highlights an example of the sound1

Abdelhakim.darraz@valeo.com

recorded by a microphone of the array during the test of a starter.

Figure 1 – Time-frequency representation of an acoustic signal acquired by the array

3. DIAGNOSIS PROCESSThe objective is to set up a diagnostic chain capable of making a binary decision based on theacquired time signals: healthy starter or defective starter. The decision organ will be a neural network.
Several approaches are envisaged and the results will then be compared. The first processing step
responds to the need to limit the sensitivity of the method to surrounding noise. This involves the useof beamforming processing (1) that exploits the number of microphones and thus allows an increasein the signal-to-noise ratio. The signals acquired by the antenna are thus focused on a point of thestarter, leading to the obtaining of a single temporal signal. It is this signal that carries the information
to judge the state of the starter.
In the present study, the choice was made not to extract the relevant parameters from this temporalsignal. Previous work has already been carried out in this direction, exploiting descriptors such asenergy in one-third octave band, Mel Frequency Cepstrum Coefficients … (2). The idea here is totrust neural networks to make the best use of the available information to make the right decision.
However, two types of neural networks will be used: LSTM recurrent networks (3) which seem well
adapted to time series and convolutional networks (4) which have provided interesting results in image
processing applications. The structure of the convolutional network is the Imagenet network (5).
The first case concerns the use of the recurrent network on the focused time signal. In the second
case, the convolutional network is used, with a time-frequency representation of the focused signal
using beamforming. Thus it is the temporal evolution of the frequency content that is highlighted fora point on the starter. In the third case, the spatial information will be further exploited. The signals
acquired by the microphone antenna are focused at several points of the starter on a regular grid.
Several temporal signals are obtained. On each of these signals, four descriptors are sequentiallycalculated: the Root Mean Squared, the Crest Factor, the Peak to Peak and the kurtosis values (6).
This leads to four mappings of the starter area representing averaged spatio-temporal information.
In order to also take into account the frequency information contained in the noise emitted by thestarter, the operation is repeated using signals that are focused and then filtered in well-chosenfrequency bands. There are three such bands: [0 1200] Hz, [2000 6500] Hz, [0 7000] Hz. These bands
were selected after analysis of the acoustic signal produced by a starter test. Finally 12 mappings will
be presented to the convolutional network.

4. DATA BASE
The three diagnostic processes described will be evaluated on a database of acoustic signals
recorded during the testing of 116 starters, 55 of which are judged by experts to be healthy and theremaining 61 to be defective. The samples have a duration of 3s and are taken at a rate of 48 kHz,consisting of 16 signals for each starter from the microphone antenna placed 35 cm away. 80% of thedatabase is used for neural network training. 10% is used for the validation of the neural model during
training, trying to avoid overlearning, and the last 10% tests the diagnostic process once the network
is fixed.
The performance of the three diagnostic processes is assessed by the indicators Accuracy, Recall
and Precision defined according to:
+

=

=

(1)

+

=

(2)

+

(3)

where Positive and Negative respectively correspond to defective and compliant starters.

5. Results
The results from the data base are presented in Tables 1 and 2.
Table 1 – Performance indicators for the three diagnosis processes for the training set
Diagnostics

process/

Accuracy

RecallPrecision

Process 1

74.8 %

76.7 %

79.7 %

Process 2

94.6 %

92.1 %

97.7 %

Process 3

79.9 %

72.1 %

87.5 %

Indicators

Table 2 – Performance indicators for the three diagnosis processes for the test set
Diagnostics

process/

Accuracy

RecallPrecision

Process 1

54.7 %

55.6 %

56.4 %

Process 2

79.4 %

75.6 %

84.2 %

Process 3

71.2 %

66.7 %

77.3 %

Indicators

6. ANALYSISThe tables clearly show that the diagnostic process working from the spectrogram of the focusedsignal gives the best results both in the learning phase and in the test. It also seems that using temporalsignals without looking for a more adequate representation is not the right strategy. Process 2, forwhich the spatial dimension of the problem was clearly considered, does not seem to provide any
improvement over the spectrogram of a single signal focused at a point on the starter. It may be that
averaging the kurtosis or Crest Factor indicators is not optimal. It is indeed probable that the structureof the signal evolves over the time period considered, which is information carried by a spectrogram
but is no longer so when calculating an indicator leading to a single value over this time period.
For all the diagnostic processes, better results can also be observed in learning than in testing. The

Precision indicator is also better than the Recall indicator showing a better potential of the processes
to detect faulty starters than healthy ones.

CONCLUSIONS
Diagnostic processes have been proposed for health monitoring on the production line ofautomotive starters. These processes use acoustic data for decision making using a MEMS microphoneantenna. They exploit the possibility of beamforming to focus the sound towards a direction, which isthe direction of a point on the starter motor and rely on artificial intelligence techniques for decision
making namely convolutional or recurrent neural networks. The results for the case study showed that
it was quite possible to obtain interesting results from an existing neural network and that the use ofa spectrogram as input was a relevant choice. As a follow-up to this work, further studies will be
carried out to see whether spatial information could not be better taken into account in the diagnostic
process.

REFERENCES
1. Van Veen BD, and Buckley KM, Beamforming: A versatile approach to spatial filtering. IEEE ASSP Mag.
1988; 5(2), 4–24.
2. Darraz A, Antoni J, Mollon P, Kirie G, Thomas JH, Larcher A, Rotating machine diagnosis using acoustic
imaging and artificial intelligence. Surveillance, Vishno and AVE conferences, INSA-Lyon, Universitéde Lyon, Jul 2019, Lyon, France.
3. Hochreiter S, and Schmidhuber J. Long short-term memory. Neural computation. 1997:1735-80.
4. Lecun Y, and Bengio Y. Convolutional networks for images, speech, and time-series. 1995.
5. Krizhevsky A, Sutskever I, and Hinton GE. Imagenet classification with deep convolutional neuralnetworks, in Advances in Neural Information Processing Systems 25, edited by Pereira F, Burges JC,
Bottou L, and Weinberger KQ (Curran Associates, Inc.); 2012. p. 1097-1105.
6. Cardenas Cabada E, Leclere Q, Antoni J, and Hamzaoui N, Fault detection in rotating machines with
beamforming: Spatial visualization of diagnosis features. Mechanical Systems and Signal Processing.
2017; 97, 33-43.

ABS-0928

Bayesian recovery of sinusoidal rough surface from scattered acousticpressure on a linear microphone array
Michael-David JOHNSON(1) , Timo LÄHIVAARA(2) , Jacques CUENCA(3) , Mansour ALKMIM(3,6) , Laurent DE RYCK(3) ,Giulio DOLCETTI(4) , Anton KRYNKIN(5)(1) Department of Mechanical Engineering, the University of Sheffield, Sheffield, United Kingdom mjohnson9@sheffield.ac.uk(2) Department of Applied Physics, University of Eastern Finland, Kuopio, Finland(3) Siemens Industry Software NV, Leuven, Belgium(4) Department of Civil, Environmental and Mechanical Engineering, University of Trento, Trento, Italy(5) Department of Mechanical Engineering, the University of Sheffield, Sheffield, United Kingdom(6) Department of Mechanical Engineering, KU Leuven, Leuven, Belgium

ABSTRACTThe recovery of the properties or geometry of a rough surface from scattered sound is of interest in many applications including medicine and river flow or structural health monitoring. Both classical inversion algorithms
and more recent data-driven machine learning techniques have shown great potential but are often sensitive to
uncertainties. In addition, the information about the confidence in the prediction is often limited. Aiming tofurther the developments of the inverse techniques and account for uncertainties, a Bayesian framework is usedin the present work to infer the properties of a rough surface, parameterised as a superposition of sinusoidal
components. The forward model used for the inversion relies on the Kirchhoff approximation of the scatteredsound above the surface of interest. The Kirchhoff approximation implies a validity region in the parameter
space which is incorporated in the Bayesian formulation. This accelerates convergence while yielding solutions
that have physical meaning. An implementation based on the Metropolis-Hastings algorithm is tested on bothsynthetic surface realisations and experimental samples produced in a laboratory setting. The robustness of theestimation in the presence of noise is evaluated.

Keywords: Inverse scattering, Kirchhoff approximation, surface reconstruction, Bayesian framework

1

INTRODUCTION

Non-intrusive acoustic inversion to recover the elevation or other properties of rough surfaces is highly importantand an active research field [1, 2, 3, 4], especially in the setting of river monitoring, where intruding the flow
could change the properties significantly. Recent work has attempted to utilise the Kirchhoff approximation invarious settings to approximate the elevation of a rough surface, first with a single frequency matrix inversion
method [5], continuing to a multiple frequency approach [6], and a machine learning approach, utilising random
forests [7]. However, although the performance of these methods have been promising, the statistical measureof the uncertainty is an open question. This report re-frames the Kirchhoff approximation into a Bayesian
viewpoint, allowing for well-developed sampling schemes such as the Adaptive Metropolis [8] algorithm tobe used to infer the uncertainties of the recovery. Further, the rough surface is parameterised into a sum ofsinusoidal functions, thus potentially reducing the number of parameters to recover, as well as allowing the useof the validity criterion of the Kirchhoff approximation as a constraint.

Figure 1. Pictoral respresentation of the experiment setup, as well as key variables that are described in the
document.

2

THEORYThe Kirchhoff approximation was used as the underlying model, obtained through the approximation of a scattered signal via an assumption based on reflections from a tangent plane. The suitability of this approximation
for the conditions analysed in this paper was presented in Krynkin et al. [9] [10], and developed further in
Dolcetti et al. [6].
Instead of directly attempting to recover the surface elevation, a sinusoidal parameterisation of the surface
is defined by


N
2πx
ζ (x, θ ) = ∑ θ1,i cos+ θ3,i ,
(1)
θ2,i
i=1
where θ , θlk ∈ R(3,N) , a 3 x N matrix, is a collection of amplitudes, wavelengths and phases that form the roughsurface, with N the number of sinusoids. For the derivation of the scattering equation, θ will be temporarily
omitted from the argument of ζ . Let the rough surface be defined by a function ζ (x) as shown in Equation (1)which satisfies the following condition [11]:sin(φ ) >1,
(kh)1/3(2)

where h is the radius of curvature of the surface, k is the acoustic wavenumber and φ is the incident angle of
the acoustic wave. When Equation (2) is satisfied, Kirchhoff approximation can be used to model reflectionsof an acoustic wave from a rough surface. In 2D, for a source with a given directivity pattern, the scattered
acoustic pressure ps can be expressed as [10]:ps (R) =12kπi

A(x, 0) ik(R1 +R2 )√(qy − qx γ)dx,
e
R1 R2
−∞

Z ∞

(3)

where, as shown in Figure 1, the values R1 and R2 are the Euclidean distance from the source at (x1 , y1 ) andreceiver at (x2 , y2 ) to a given point (x, ζ (x)) on the surface, respectivelyqR1 = (x1 − x)2 + (y1 − ζ (x))2 ,
(4)
R2 =

q

(x2 − x)2 + (y2 − ζ (x))2 .

(5)

In Equation (3), R = (x2 , y2 ), γ = dζ (x)/dx; qx and qy are the x and y components of q = −k∇S (R1 + R2 ) withthe gradient ∇S = (∂ /∂ x, ∂ /∂ y). The directivity term A(r) is far-field radiation from a baffled piston, which isgiven by [12],2J1 (ka sin(φ (r) − (−φ0 + π/2))),(6)
ka sin(φ (r) − (−φ0 + π/2))
where a is the aperture, J1 is the Bessel function of the first kind, φ0 is the angle of inclination of the source
main axis to the Ox-axis, and φ (r) is the angle between the vector produced from the location of the source
and the point r with the Oy-axis.
This can be generalised in matrix form, when there are M receivers as
A(r) =

p = {p(R j ), j = 1, .., M}.

(7)

In practical applications, noise is present in measured data as a sum of many noise factors such as backgroundnoise and electrical interference. In the case of Gaussian additive noise assumption, the equation becomesp̃ = p + ε,(8)where ε ∼ N (0, Iσ 2 ) is drawn from a normal distribution and there is no covariance between receivers. Anotherapproach to Equation (8) is to treat the overall equation in the perspective of probability. In the Bayesianframework, Equation (8) can be interpreted as
P(θ |p̃) ∝ P(p̃|θ )P(θ ).

(9)

The posterior distribution is proportional to the likelihood function P(p̃|θ ) multiplied by the prior belief P(θ ),which can include further constraints on the physicality, by incorporating Equation (2). The Metropolis-Hastings
algorithm [13] was used as the model to sample from the posterior distribution, with a further modification
named the Adaptive Metropolis scheme [8]. A uniform prior was selected for all parameters, only allowing
surfaces which satisfied the condition (2).

3

EXPERIMENT SETUP

Experimental data used in this paper for validation purposes was collected with 34 1/4′′ microphones arranged
in a linear orientation, and a loudspeaker (Visaton G 25 FFL), arranged in a linear orientation. The specifics ofthe setup can be found in [7]. A sinusoidal surface (with amplitude ζ1 = 0.0015 m and wavelength ζ2 = 0.05m) was machined from an aluminium block with a length of 0.55 m in the x-direction, and a width of 0.35
m. A broadband acoustic signal was recovered, and 14,000 Hz was selected as the frequency to be used in
the simulation. In Figure 4, a rough surface composed of multiple sinusoids was created from fibreboard, and18,650 Hz was selected.

4

DISCUSSION

Figure 2 presents the distributions and densities of each parameter recovered for the sinusoidal surface using
the Adaptive Metropolis scheme. The mode of the distributions are close to the true parameter values (errors
of approximately 15% for amplitude, 2% for wavelength, and close to 0 for phase). The amplitude has rather awide spread, which could be because of the frequency selected for the simulation being larger than the surface
amplitude. The densities show that this recovery is not multi-modal. Figure 3 presents the actual surface
recovery of the mean function, as well as a shaded region containing the parameters that form a 68% credibleinterval (also reported in Figure 2, above the plots). It can be observed that the first half of the surface isrecovered better than the latter. This could be because of the finite receiver setup, as well the distribution of
the energy given from the source location and angle. To showcase an extension past three-parameters, a fortyparameter surface was also recovered (N=20, wavelengths are fixed). Here, similar behaviour of the method
can be observed, i.e. fitting better in the first half of the surface. However, the uncertainty is rather wide, byobservation of the samples in Figure 4. This could be improved by incorporating a multiple frequency approach,as explored in Ref. [6].
This algorithm can be used as a method of acquired model certainty, and could be used as a final step in asurface recovery pipeline.

Figure 2. Corner plot produced with [14]. The leading diagonal contains the distributions of the amplitude,
wavelength, and phase respectively. And inside presents the density between a given 2 parameters.

Figure 3. The true 3 parameter surface, plotted against the mean surface provided from taking the mean predictions from the algorithm. The filled region contributes to surfaces created from the bounds shown in Figure 2

Figure 4. The true experimental 40 parameter surface, plotted against the mean surface provided by the algorithm. The black line are random samples from the algorithm.

REFERENCES[1] G. Dolcetti, A. Krynkin, K. V. Horoshenkov, and S. J. Tait, “An acoustic technique to measure the velocityof shallow turbulent flows remotely,” in Free Surface Flows and Transport Processes, pp. 181–194, Springer,
2018.
[2] A. Nichols, S. Tait, K. Horoshenkov, and S. Shepherd, “A non-invasive airborne wave monitor,” Flow
Measurement and Instrumentation, vol. 34, pp. 118–126, 2013.
[3] S. P. Walstead and G. B. Deane, “Reconstructing surface wave profiles from reflected acoustic pulses using
multiple receivers,” The Journal of the Acoustical Society of America, vol. 136, no. 2, pp. 604–613, 2014.
[4] L. Hellequin, J.-M. Boucher, and X. Lurton, “Processing of high-frequency multibeam echo sounder data
for seafloor characterization,” IEEE Journal of Oceanic Engineering, vol. 28, no. 1, pp. 78–89, 2003.
[5] A. Krynkin, K. Horoshenkov, A. Nichols, and S. Tait, “A non-invasive acoustical method to measure
the mean roughness height of the free surface of a turbulent shallow water flow,” Review of Scientific
Instruments, vol. 85, no. 11, p. 114902, 2014.
[6] G. Dolcetti, M. Alkmim, J. Cuenca, L. De Ryck, and A. Krynkin, “Robust reconstruction of scattering
surfaces using a linear microphone array,” Journal of Sound and Vibration, vol. 494, p. 115902, 2021.
[7] M.-D. Johnson, A. Krynkin, G. Dolcetti, M. Alkmim, J. Cuenca, and L. D. Ryck, “Surface shape reconstruction from phaseless scattered acoustic data using a random forest algorithm,” Journal of the AcousticalSociety of America, vol. Accepted, 2022.
[8] H. Haario, E. Saksman, and J. Tamminen, “An adaptive Metropolis algorithm,” Bernoulli, vol. 7, pp. 223–
242, 2001.
[9] A. Krynkin, K. V. Horoshenkov, and T. Van Renterghem, “An airborne acoustic method to reconstruct a
dynamically rough flow surface,” The Journal of the Acoustical Society of America, vol. 140, no. 3, pp. 2064–
2073, 2016.

[10] A. Krynkin, G. Dolcetti, and S. Hunting, “Acoustic imaging in application to reconstruction of rough rigid
surface with airborne ultrasound waves,” Review of Scientific Instruments, vol. 88, no. 2, p. 024901, 2017.
[11] F. G. Bass and I. M. Fuks, Wave Scattering from Statistically Rough Surfaces: International Series in Natural
Philosophy, vol. 93. Elsevier, 2013.
[12] P. M. Morse and K. U. Ingard, Theoretical acoustics. Princeton university press, 1986.
[13] W. K. Hastings, “Monte carlo sampling methods using markov chains and their applications,” 1970.
[14] D. Foreman-Mackey et al., “corner. py: Scatterplot matrices in python.,” J. Open Source Softw., vol. 1,
no. 2, p. 24, 2016.

ABS-0557

Numerical analysis on defect detection using structural intensity in solid
bodiesAxel Helmut SCHUHMANN(1,2) ; Nikolai KLEINFELLER(1) ; Shahrokh SEPEHRIRAHNAMA(2) ; Sebastian OBERST(2) ;Christian ADAMS(1) ; Tobias MELZ(1)(1) Technical University Darmstadt, Mechanical Engineering, Department of System Reliability, Adaptronics and Machine

Acoustics, Germany, axel_schuhmann@outlook.de(2) University of Technology Sydney, Faculty of Engineering and Information Technology, Centre for Audio, Acoustics andVibration (CAAV), School of Mechanical and Mechatronic Engineering Australia, Sebastian.Oberst@uts.edu.au

ABSTRACT
Analysing structural intensity (SI) offers the possibility to assess the transmission of wave energy within astructure. Measurement of SI has been mainly focused on thin shells and beams. In this work a measurement
method is presented to evaluate SI within solid, homogeneous, and isotropic bodies. The method is based on
the reciprocity principle, a fundamental assumption in linear vibroacoustics. It allows the reconstruction of thestructural intensity field within the bulk of a solid body from the measured surface velocities on the exterior
boundaries. From the preliminary results, we demonstrate the capability of this method in approximating the
spatial variation of the reconstructed stress and velocity fields using finite element simulation results. Inspiredby the reciprocity-based method, we also demonstrate a cavity detection technique using the structural intensity
measured along a closed path on a surface of a solid block. Despite some discrepancies in the estimate of themagnitude, the method works well in principle for the benchmark problem of a rectangular cube and it can beverified using our recently set up experimental test. Our proposed method provides an alternative energy-based
SI detection technique that may perform as well as those exploiting velocity/acceleration or strain/stress.

Keywords: Structural Intensity, Measurement, Damage Detection

1

INTRODUCTION

A structural intensity’s (SI) vector field describes the propagation of structure-borne wave energy within anelastic body [1]. Researchers aim to use the information of the SI vector field to better understand structuresand gain new insights into the propagation and influence of structure-borne wave propagation. Some recent
applications include the analysis of mechanical structures, e. g., for the noise reduction in airplanes [2] or theinvestigation of transmission paths within the front of a car [3], as well as bio-mechanical applications, e. g.,
to bones fractures due to head injuries [4] or studies about the energy propagation of the human tympanic
membrane [5]. A damaged beam or shell structure has an altered SI vector field. Several studies [6, 7, 8, 9, 10]
demonstrate that it is still challenging to extract the necessary information from the SI vector field to analyseor localise damages in these structures.
In the 1990s, the first formulation for the prediction of the SI based on numerical models was developed, thus
giving researchers the opportunity to investigate more geometrically complex structures or volume structures[11, 12]. Measurement techniques for SI use shell [13, 14] or plate theories [15, 17] to estimate the stressterms from displacement measurements. A measurement method that is applicable for volume structures does
not exist, yet. This raises the question whether it might be possible to develop such a measurement method forthe estimation of the SI vector field in solid volume bodies. More precisely, how can surface measurements offer
the necessary information to determine the SI vector field accurately enough within a solid volume structure?
This leads to the next question: How can the information based on these surface measurements be used as

irregularity detection indicator?
This paper is a first step to answer these two questions and structured as in the following. A brief introduction is
given before the theoretical approach for the estimation of full-field internal data based on surface information
published by H OLLAND et al. [18] is outlined. Verification of the method is conducted using a benchmark
model which is presented in a subsequent section. The numerical implementation in the software A BAQUS isapplied and first promising results are calculated. Finally, the work concludes with an analysis covering the ideaof surface measurements as irregularity detection indicators based on SI magnitudes at the surface.

2

THEORYThis section provides the theoretical background for the description of energy propagation in volume structures
based on SI. Furthermore, an approach is presented on the basis of which dynamic full-field internal stresses
can be predicted through 3D surface laser D OPPLER vibrometry (3D-SLDV) .
2.1 Structural intensity
A fundamental idea of analysing energy flow is the intensity concept. Analogously to sound intensity, SI describes the energy flow per unit area. The time averaged SI ⟨Ii (t)⟩ can be calculated by using the complex
stress σ̃i j (ω) tensor and the complex conjugate velocity vector ṽ∗j (ω) [1], thus giving SI for mono-frequent,stationary signals as1 (1)⟨Ii (t)⟩ = − Re σ̃i j (ω) · ṽ∗j (ω) .
2
In Equation (1) the tilde symbol ( ˜ ) denotes complex values and the asterisk (∗ ) indicates complex conjugate
values with ω being the angular frequency.
Equation (1) is valid for any elastic solid body. Usually, SI is mostly analysed for the cases of thin shells
and beams because due to mathematical assumptions like K IRCHHOFF-Hypothesis the calculation of SI can bedone for the shell’s mid-plane or the beam’s neutral axis [17, 13, 14]. The main reason of this circumstanceis the lack of availability of the stresses σ̃i j (ω) and velocities ṽ∗j (ω) within the structure under investigation.
For a successful SI analysis of a volume structure in the experimental context, the knowledge of the full-fieldinternal stresses σ̃i j (ω) and the corresponding velocities ṽ j (ω) is required as explained in the following.
2.2 Prediction of full-field internal stresses of volume structures
The approach to estimate the full-field internal stresses σ̃i j (ω) is based on the work of H OLLAND et al. [18].
The fundamental idea is to utilise the principles of the structural boundary element method which is experimentally validated to be capable of determining the internal dynamic stresses σ̃i j (ω) and strains ε̃i j (ω) in a solid,
homogeneous, isotropic, and vibrating body [18]. For this purpose, all components of the velocity vector ṽ j on
the entire surface S of the structure under investigation have to be measured to estimate the full-field internal
stresses σ̃i j (ω). First, a reciprocity integral formulation [19]

Z 
Z
(1)
(1) (2)(2)(2) (1)(1) (2)(1) (1)( fi − ρ üi )ui − ( fi − ρ üi )ui dV = (σi j ui − σi j ui ) n j dS.
(2)
V

S

is applied to connect the behavior of two different states (1) and (2) of the structure under investigation. The
applicability of the reciprocity integral formulation assumes the linearity of the structure. Figure 1 shows the twostates (1) and (2). State (1) describes the measured state. The full-field internal stresses σ̃i j (ω) and velocitiesṽ j required for SI analysis are estimated on the basis of surface velocities ṽ j . State (2) is initially unknown,
G:yand to describe it, the G REEN’s functions G of displacements ũG:yj:k and internal stresses σ̃i j:k,i which are thesolutions of the equations of motion
G:y˜σ̃iG:y(3)
j:k,i + f i = ρ ü j:k
to a time-harmonic body forcef˜i = δ (x − y)e−iωt(4)

(1): Measured

(1)

Ti

(2): Calculated

(2)

=0

σ̃i j ni
x
r

(1)

ũi

z

y
y

(2)
f˜i

x
(1)

fi

(2)
at x = y : f˜i = δ (x − y)e−iωt δik=0

Figure 1. Methods of calculating structural intensity in a solid body using (a) the velocity and stresses acrossthe volume, and (b) the surface velocity and mechanical traction in conjunction with the reciprocity principle.
in the kth direction at a driving point x = y within the structure, are applied. The superscript G : y denotesthat the variable is based on a fundamental solution due to the excitation at the position y. The indices j : krepresent the jth component of the displacement vector ũ based on an excitation of the time-harmonic bodyforce f˜i in the kth direction. Here, the G REEN’s function of the displacements [20]1(U1 δ jk −U2 r, j r,k )4π µ ik r  2 
 ik reikT r
i
1
e T
kL
i
1
e L
U1 =
+
−
−
−
,
2
2
r
kT r (kT r)
r
kT
kL r (kL r)r ik r  2 
 ik r
3i
3
e T
kL
3i
3
e L
U2 = 1 +
−
−
1+
−
.
2
2
kT r (kT r)
r
kT
kL r (kL r)
r
ũG:y
j:k =

(5)

for the general three-dimensional case are solved in the spatial domain with the wave numbers
kT =

ω
ω
and kL = .
qT
qL

(6)

and their corresponding phase velocities qT , qL . The subscript T denotes shear waves while L represents longitudinal waves. The distancer = |y − x|(7)changes with the position x for an arbitrary point within the structure with respect to a fixed driving point y.
Next, H OOK’s law


G:y
G:y
G:y
σ̃iG:y
=
λ
ũ
δ
+
µ
ũ
+
ũ
(8)
i
j
j:k, j
mk,m
ik, j
jk,i
with L AMÉ-Parameters λ , µ and the first order spatial derivatives of the displacements ũG:y in the m, j-directions
offer the G REEN’s function of the internal stresses σ̃iG:yj:k, j [20]. The estimation of the displacements ũk on thebasis of ACHENBACH’s derivation [19]ũk y =

1

Z

e−iωt

S

(−σ̃iG:yj:k u j ) ni dS(9)is given by a surface integral based on the knowledge of displacements u j at the surface S of the structure
under investigation. The derivation of Equation (9) combines the reciprocity integral formulation (Equation (2))

G:yand the G REEN’s function for the displacement ũG:yj:k (Equation (5)) as well as the stresses σ̃i j:k, j (Equation (8)).
The spatial derivatives of the displacements ũkε̃kl y =ũk,l + ũl,k−1=22iωe−iωt

Z h

i
G:y
(σ̃iG:y
+
σ̃
)v
ni dS
j
j:k,l
i j:l,kS(10)

approximate the strain tensor ε̃kl . The displacements u j on the surface S can either be estimated by measurements or simulations. In Equation (10), the displacements u j have been replaced by the velocities v j including
the factor 1/iω since the measured quantity using the 3D-SLDV is the surface velocity. The final step is toestimate the internal stresses σ˜i j
σ̃i j = 2µ ε̃i j + λ ε̃kk δi j
(11)
by using H OOKE’s law. Since the partial derivative of displacements uk with respect to time tṽk y =∂ ũk y
∂t

=

1

Z

e−iωt

S

(−σ̃iG:y
j:k v j ) ni dS(12)provide the velocities ṽk , the components of the Equation (1) required to calculate SI can be obtained. In conclusion, Equations (11) and (12) offer an option to estimate the complex stress tensor σ̃i j (ω) and the complex
velocity vector v˜j ∗ (ω) at any internal point of the structure under investigation, if the velocities v j across theentire surface are available.
Consequently, SI can be estimated at any arbitrary position y within the volume structure by Equation(1). As stated before, the method assumes a solid, homogeneous, isotropic, vibrating, volume structure withaccessible faces for measuring the velocities v j across the entire surface. Additionally, it implies a free boundary
condition. It states that under the conditions mentioned above, surface measurements are theoretically sufficientto estimate the SI within a volume structure. In the next two Sections, the method is implemented and verified
for the example of a benchmark structure.

3

BENCHMARK STRUCTURE AND ITS NUMERICAL MODELLING

Fig. 3 presents the benchmark model that is used in this work. The model’s characteristics are presented inTable 1. In Fig. 3 a) the dynamic excitation force F acts at the top of the model, while in b) a discrete damper
element is attached at the bottom of the model. This enables energy to be transferred from the excitation
location (source) to the damper element (sink). The numerical model has been simulated using the FE software
A BAQUS 2021.HF3. In a steady-state-dynamic analysis the complex stress tensors at the integration points aswell as their complex velocity vectors are obtained.
Table 1. Material properties benchmark model
Dimensions:

100 mm · 100 mm · 200 mm

Modulus of Elasticity:

70000 N/mm2

Density:

2.7 g/mm3

Poison’s Ratio:

0.35

2000 linear C3D8R elements have been used for the simulation and apart from a discrete viscous damping
element no other dissipating mechanism such as structural (material) or R AYLEIGH damping has been applied.
The number of elements ensures the fast determination of numerical results and limits the number of surface
elements to be approximated. The C3D8R linear element is chosen because it offers an easy implementationof the method due to the limited numbers of nodes. For a more accurate but more computational expensive
representation quadratic elements or fully integrated, incompatible modes element formulations can be chosen[16]. However, the numerical model has been replicated in an experimental setup and it’s characterization has

a)

c)

b)

Frame

5N

SpringRubberSprings

bands

200 mmShakerImpedancehead
100 mm

100 mmDamperTest object

Figure 2. A finite element model of a solid rectangular cube made of Aluminium with (a) a harmonic point
force F applied on the top, and (b) four springs attached at between the bottom face and the ground. A dashpot
damper is placed underneath to allow the flow of structural intensity. Panel (c) shows the recently preparedexperimental set up based on the standard shaker test and corresponding to the proposed numerical model. The
experimental results are yet to be collected.
been conducted and found to be sufficient at the time of investigation. This investigation serves as a preparation
for the application of the method on a real test object. The experimental set-up is presented in Fig. 2 c). The
exciting force is applied by a shaker. The shaker is decoupled from the test object using a soft spring andrubber bands. After verification of the method using numerical data as input, the experimental set-up can beused to measure the surface velocities required for the method.

4

MEASUREMENT METHOD

4.1 Numerical implementation
Section 2 offers the theoretical approach to estimate the internal stress tensors as well as the velocity vectors
at an arbitrary point y. Using this information the SI vector field can be calculated. From the continuous andprecise formulation of the surface integrals in Equations (11) and (12) an approximation must be derived toapply it on a real world problem.
In this work the method is tested using a numerical model. Therefore, input data for the method is extracted
from an Abaqus Output Database file (.odb). An approximation of the surface integrals from Equation (10)and (12) is applied using the mid-point rule which can be found in [21]. The calculation is performed using a
Matlab script. Each integrand is summarised into a function eskl (x, y, ω, n, v) or osk (x, y, ω, n, v). For the strain
tensor at an internal point y the approximation is then given asε̃kl y ≈−1 · ∑Ns=1 eskl (x, y, ω, n, v)∆A iωt
·e .
2iω

(13)

Fig. 3 is illustrated to gain a better understanding of the method. It shows a scheme of the numerical volume
model with its local coordinate system in the bottom left corner. For the calculations of internal stress components at a specific point, as shown in Equation (13), it is necessary to calculate the integrand eskl (x, y, ω, n, v)
at the centre of each partial surface and multiply it with the corresponding area ∆A. The point where thestress tensor or velocity vector is calculated is defined by the vector y. The input value y for the function

F∆A

osk or eskl

v

n

x

r

z
x

y
y

Figure 3. Schematic representation of the measurement method of SI using the surface velocity v and discretization into finite boundary element. The principle of reciprocity is being used to estimate the SI at the point ywithin the volume of the solid.
eskl (x, y, ω, n, v) is consistent for one specific point. In contrast to that, the x vector changes for each surface
area element. According to the face of the model the unit normal vector n changes, too. Since velocities v of
the numerical model are given at the nodes of the surface elements, four nodes are averaged to calculate thevelocity vector in the centre of it. This is valid for regular meshes. In case of using a different element type,
implying an irregular mesh, weighting of the velocity components at the nodes is required. The excitation force
F acts at a constant frequency. Afterwards the products of all partial surfaces need to be summed up. The
surface integral is thus approximated. It needs to be mentioned that the whole model is discretised. For thesake of clarity only the top surface shows the discretisation in Fig. 3.
4.2 Results
Testing the method is performed by comparing the method’s results with the results generated by the numericalsimulation. Results for a specific point exactly in the middle of the numerical model offered promising findings
for the velocity vector and and the stress tensor. The analysis was therefore extended to a qualitative comparisonof the middle plane of the model. In Fig. 4 the plane is illustrated.
x-z-plane at y = 0.05z
x y

Figure 4. Visualization of the middle plane where results are calculated forFig. 5 shows the distribution of the real normal stress magnitudes in the excitation direction for a section of
the middle plane. The magnitudes in Fig. 5 have been normalised using the maximum value from the method’s
results. In Fig. 5 a) the distribution on the basis of the method is shown while in b) the distribution based

Figure 5. Contour plots of the real normalised stress field Re {σ̃zz } at the middle section parallel to the z-x
plane. Panel (a) shows the reconstructed stress field using the reciprocity-based method, while panel (b) shows
the stresses directly from the finite element simulations.
on the FE simulation results are presented. The contours show almost identical spread of Re{σ̃zz } with the
middle cross section. The magnitude of the reconstructed stress field, cf. Fig. 5 a), is two orders of magnitude
smaller than the finite elements results, cf. Fig. 5 b). Despite this discrepancy in the magnitude of the stresses
due to numerical errors based on mesh characteristics, single-precision computations, and the use of linear, but
reduced integrated solid elements (which could lead to hourglassing [16]), the high-quality match between thefield variations indicate the potential of using the surface-measured velocity to qualitatively estimate the stressfield within the middle section. Considering that this is a preliminary analysis, it implies the feasibility of
improving this method as an alternative energy-based method for estimation of SI in a solid body.
Fig. 6 shows the qualitative comparison of the real part of the velocity component Re{ṽz } in excitation direction
z for the same section of the middle plane as in Fig. 5. The magnitudes of Re{ṽz } have been normalised by thecorresponding maximum magnitude from the reconstructed field and the FE results. The colour coding is scaled
to include maximum and minimum of both results. It is noted that the reconstructed field shows approximately
the same velocity throughout the plane-section implying a uniform distribution. The FE results also show asimilar uniform velocity distribution across the plane-section compared to the new method, except for somenodal points where the numerical errors are large and causes variation in the contour values. Similar to thestress results, the comparison of the velocity fields ṽ j show a discrepancy in the estimated magnitude of thereconstructed field in Fig. 6 a), which can be due to the aforementioned reasons. Nevertheless, the qualitative
comparison of the contours shows the capability of the reconstruction method to estimate the uniform velocity
field within the middle section.
4.3 Discussion
For the real normal stress component in excitation direction, which is the dominant component, a good qualitative agreement could be found between the method’s results and the FE simulation results. Apart from somenumerical artefacts, the real velocity component in excitation direction offers a uniform distribution in the FE
simulation results. The uniform distribution can also be found in the reconstructed field. Both times, for stressas well as velocity, only the qualitative distributions are matching while quantitative confirmation is still pending. A possible source of error for the inconsistent magnitudes has already been investigated. The missing
sub-tessellation algorithm for the boundaries of the model, mentioned in [18], is not found to dominate results and serve as a key problem. Nevertheless, there are other problems that might cause the inconsistency
between the results’ magnitudes. As derived in Section 2, the method requires a free boundary condition forFigure 6. Contour plots of the normalised velocity results Re{ṽz } at the middle section parallel to the z-x plane.
Panel (a) shows the reconstructed velocity field using the reciprocity-based method, while panel (b) shows thevelocity results directly from the finite element simulations.
the measured state (1). Such conditions can only be approximated in practice, especially, at excitation points
and positions of damping. An adapted integral formulation for those two points represents a possible solutionthat will be investigated in the near future. Further, results in [18] are validated using the flexural wave theory which is based on the E ULER -B ERNOULLI beam theory, including the ’normal hypothesis’ which does notallow for shear deformations. As a result, certain types of waves such as longitudinal strain waves remain disregarded. This might be another reason for the inconsistent magnitudes. The force in this work rather excites
longitudinal strain waves instead of transversal flexural waves. A way of examining this problem would be tochange the exciting force so that mainly flexural waves are excited and compare the results again. Apart fromthat the most rudimentary numerical element C3D8R has been used. This is another possible reason for thelack of matching magnitudes. Since the results are preliminary further investigations are required to evaluate
the methods potential and verify it before the method shall be tested on the experimental set-up presented in
Section 3.

5

SI AS DAMAGE DETECTION INDICATOR IN VOLUME STRUCTURES

SI is intended to be used to either detect cracks in plates [6], diagnose flaws in beams [7], or to identify
damaged plates [10]. In these papers SI deals as the quantity that yields information on flaws, cracks or
damages generally. In this work it is investigated whether SI measurements on the surface might be utilised fordamage detection such as a cavity inside a bulk structure. The analysis is based on numerical results.
5.1 Analysis of the global SI distribution
As a first step the global SI distributions are analysed. To do so, a post-processing code to calculate SI vector
fields is created and applied. It is based on Equation 1. Two models are compared which differ only ina cylindrical cavity within one of the models. The cylindrical cavity models a damage within the volume
structure. The cavity is 60 mm long and its diameter measures 20 mm. The models are similar to the one
analysed in Section 4. Fig. 7 shows the global SI vector field distribution of the two models using an in-housevisualization method [22]. The cylindrical cavity seems to act as a disruption for the energy flow. The SI
vector field in a) appears to be more ordered than in b), which is due to the different mesh. From a global
view the influence of the cavity is hardly visible. Only when zoomed in the effect of the cavity is noticeable.
Therefore, in figure 7 c) a zoomed section of b) is presented. The cavity is visible in this figure. It shows the
SI vectors spreading around the cavity. Since there is less space available to transport energy from excitation

Figure 7. The SI field for (a) an isotropic and flawless solid block, compared to (b) the SI field in the same
block with a cylindrical cavity along the z-axis at the middle. Panel (c) shows the zoomed-in view of the SI
field around the cavity. In both cases, a point force is applied at the middle of the the top and a dashpot isplaced at the same location on the opposite face (bottom face), and the block is supported by springs on the
ground, as shown in Fig. 2
point towards the damper element, it is assumed that magnitudes of the vectors around the cavity in a) arehigher than compared to the same vectors in b). This assumption should also be valid for SI vectors at thesurface. It is therefore reasonable to analyse the two models’ SI surface vectors, velocity vectors, stress tensors
and its components.
5.2 Analysis of SI distribution along a path
The analysis of the surface quantities is conducted by the use of path plots in A BAQUS. Paths plots enable a
comparison of numerically calculated values along a specified path. In this work the paths are placed on a faceof the numerical model. Fig. 8 shows the position of the two paths. In Fig. 8 a), the path is created directly
over the cavity, aiming for the surface quantities to contain the information about the cavity. In b) it is found at
the same place but in the model without the cavity. Figure 9 shows the SI vector components along the path of
a)

path

b)

Figure 8. Rectangular paths on the side face to plot the structural intensity along for (a) the model with acylindrical cavity, and (b) the reference flaw-less model.
each model. Comparing the different values offers the following. First, the z-component of the SI vectors yield
the largest magnitude. The y-components are close to zero along these paths. Additionally, the x-components
include both positive and negative SI values and are also close to zero. Moreover, for each component x, y, orFigure 9. Comparison of SI vector components along the path, shown in Fig. 8, to study the effects of thepresence of the cavity and demonstrate the potential for cavity detection using path-based measurement of SI.
z, the graphs are qualitatively very similar to each other. The fact that the magnitude of the z-component in themodel with the cylindrical cavity is larger than without the cavity is again supporting the assumption that the
reduced space leads to higher SI magnitudes, which could be used for damage detection.
The analysis conducted in this section offers first promising results. The idea that damage detection could beachieved by analysing magnitudes of SI at the surface seems to be supported. The result that due to the reduced
space SI magnitude in z-direction for the model with cavity are larger is supporting the general applicabilityof the idea. SI analysis at the surface might indeed offer enough information to detect flaws within a volume
structure.
Nevertheless, the analysis in this work is far from complete. There are different limiting factors that need tobe considered. The findings in this work are only valid for one specific path. A more complete understandingof the magnitude differences requires the analysis of more paths along the surface. Moreover, the comparisonof parallel sides of the paths from different models could offer a better understanding. Another limiting factor
is that the cavity in this work represents no real damage, flaw, or crack. A more accurate model might offerdifferent insights. It needs to be shown statistically that there is a significant difference found for different cavity
sizes or even cracks when compared with flawless models. A benchmark cavity could help to investigate therelative influence of the cavity’s size. Besides, further investigations are required to evaluate if the differences inmagnitude are even measurable. Depending on the measurement system that is used, this might not be possible
or too expensive to achieve the required resolution.

6

CONCLUSIONSThis work is motivated by two questions raised in the introduction in Section 1. Recapturing the first question
from the introduction: Might it be possible to develop a measurement method for the estimation of the SI vector
field in solid volume bodies. More precisely, how can surface measurements offer the necessary information todetermine the entire SI vector field within a solid volume structure? In this work, a theoretical approach fora measurement technique of SI within a solid, homogeneous, and isotropic body with accessible faces on thebasis of the findings in [18] is presented. The theoretical background is explained. A brief introduction of theimplementation is outlined and the numerical model for the test of the method is shown. Preliminary results
are presented and offer promising results. Further investigations need to be conducted. In subsequent projects
the method shall be tested on an experimental set-up that has already been prepared.
The second question from the introduction was: How can the information based on these surface measurements
be used as irregularity detection indicator? In this work the potential of SI as irregularity detection indicator
has been shown by SI surface magnitude analysis using paths. However, only one specific path is analysed andno statement can be made about general application. Nevertheless, the first promising findings should not bedisregarded. On the contrary, they should be considered a possible starting point for further studies in this field.
Additional paths and planes need to be analysed to possibly verify the first promising findings. Furthermore,different shapes and sizes of the cavity and their influence on SI could be investigated. A realistic cavity or
crack model might offer new insights. Finally, experimental studies need to show that the numerical results
match reality. The measurement of SI on a test object should therefore be approached.

ACKNOWLEDGEMENTS
The results in this paper were obtained as part of a master’s thesis carried out in cooperation with the Universityof Technology Sydney. The first author thanks all co-authors for input and help during this time.

REFERENCES[1] Pavić G. Structure–borne energy flow. Handbook of noise and vibration control, John Wiley & Sons, Inc,
Hoboken, New Jersey. USA 2007.
[2] Winter R, Heyen S, Biedermann J B. Experimental Structure Intensity Analysis of an Airbus A400M fuselage structure using high-resolution vibration measurements. Proceedings of the 23rd International Congresson Acoustics, Aachen 2019. p.8118.
[3] Stoewer T, Ebert J, Melz T. Opportunities and limitations of structural intensity calculation regarding uncertainties in the NVH design of complex vehicle body structures. Uncertainty in Mechanical Engineering II,
807, 2015. p.78–86.
[4] Lee H P. Assessment of head injury of children due to golf ball impact. Computer Methods in Biomechanics
and Biomedical Engineering 13(5), 2010. p.1476–8259.
[5] Pires F. Structural intensity assessment & material identification of a human tympanic membrane. PhD thesis, Universität von Antwerpen. Antwerpen 2020.
[6] Lee H P, Lim S P, Khun M S. Diversion of energy flow near crack tips of a vibrating plate using thestructural intensity technique. Journal of Sound and Vibration, 296(3), 2006, p.602–622.
[7] Li T Y, Zhang W H, Liu T G. Vibrational power flow analysis of damaged beam structures. Journal ofSound and Vibration, 242(1), 2001. p.59–68.
[8] Li T Y, Zhu X, Zhao Y, Hu X F. The wave propagation and vibrational energy flow characteristics of a platewith a part-through surface crack. International Journal of Engineering Science, 47(10), 2009. p.1025–1037.
[9] Tian X, Liu G, Gao Z, Chen P, Mu W. Crack detection in offshore platform structure based on structuralintensity approach. Journal of Sound and Vibration, 389, 2017. p.236–249.
[10] Wong W O, Wang X Q, Cheng L. Modal power flow analysis of a damaged plate. Journal of Sound andVibration, 320(1-2), 2009. p.84–100.
[11] Gavrić L, Pavić G. A finite element method for computation of structural intensity by the normal mode
approach. Journal of Sound and Vibration, 164(1), 1993. p.29–43.

[12] Hambric S. Power flow and mechanical intensity calculations in structural finite element analysis. Journalof Vibration and Acoustics, 112(4), 1990. p.1048–9002.
[13] Pires, F. and Vanlanduit, S. and Dirckx, J. J. J.. Structural intensity analysis on irregular shells. Journal of
Vibration and Acoustics, 141(1), 2019, 011011-1–011011-12.
[14] Pires F, Avril S, Vanlanduit S, Dirckx J J J. Structural intensity assessment on shells via a finite element
approximation. The Journal of the Acoustical Society of America, 145(1), 2019. p.312–326.
[15] Eck T, Walsh S J. Measurement of vibrational energy flow in a plate with high energy flow boundary
crossing using electronic speckle pattern interferometry. Applied Acoustics, 73(9), 2012. p.936–951.
[16] Oberst S, Lai JCS, Marburg S. Guidelines for numerical vibration and acoustic analysis of disc brakesqueal using simple models of brake systems. Journal of Sound and Vibration, 332(9), 2013. p.2284 2284
– 2299.
[17] Pires F, Muyshondt P G G, Keustermans W, Vanlanduit S, Roozen N B, Dirckx J J J. Structural intensity analysis of flat plates based on digital stroboscopic holography measurements. Journal of Sound andVibration, 428, 2018. p.168–178.
[18] Holland S D, Renshaw J, Roberts R. Measurement of dynamic full-field internal stresses through surface
laser Doppler vibrometry. Applied Physics Letters, 91 (13), 2007. p.134101.
[19] Achenbach J D. Reciprocity in elastodynamics. Cambridge monographs on mechanics, 2003.
[20] Beskos D E. Boundary element methods in mechanics /ed. by Dimitri E. Beskos. Mechanics and mathematical methods Series 3, Computational methods in mechanics, 1987.
[21] Schäfer M. Computational Engineering - introduction to numerical methods. Springer, 2006.
[22] Zimmermann R. Entwicklung einer Visualisierungsmethode zur Auswertung und Analyse der Strukturintensität an Volumenstrukturen (Development of a visualization method for the evaluation and analysis of thestructure intensity on volume structures). Bachelor thesis, TU-Darmstadt, Darmstadt, 2019.