PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Numerical simulation of vibration damping by granular materials
Masao TAKESHIMA1; Takumi ASAKURA2
12Tokyo University of Science, Japan

ABSTRACTIn recent years, there have been examples of efforts related to vibration reduction of plate vibration using
granular materials, impact dampers, etc. On the other hand, as a means of analyzing the influence of theapplication of granular material on the vibration characteristics of the plate-shaped structure, it is currently
said that little knowledge has been obtained. In this report, a simulation method of the damping effect ofbending oscillation of the plate by granular material combining following two kinds of methods was proposed.
One is the discrete element method (DEM) for analyzing the motion of the granular material. The other isthe finite-difference time-domain (FDTD) method for vibration analysis of the thin plate.
In this proposed method, influence of the granular material on the plate vibration was considered. Validityof the proposed method was confirmed by comparison between the results of the calculation and the modelexperiment. Then, this study confirmed the dependency of the damping effect of the plate vibration on theexternal excitation force caused by the granular material.
Keywords: Vibration damping, Wave-based numerical simulation, Granular material, FDTD.

1. INTRODUCTION
Techniques of vibration control are one of the most important factors to improve in the quality of
our current living environment, and those has been developed in various ways. Among them, thevibration control for plate-like structures such as ceilings and floors in the architectural fields has
been well studied and has been also investigated for many years from the viewpoint of a floor-impact
sound problem between dwelling units on the upper and lower floors in a residential space. The noise
control technique which realizes sufficient damping effect of the ceiling structure by laying charcoal
bag on the ceiling board has been reported [1]. In these researches, it is supposed that the damping
effect may be caused by the collision of the granular material and the plate-like oscillator. The
vibration control effect has been determined experimentally in various papers, howeveranalysis
method has not been sufficiently yet investigated. It is difficult to accurately predict the damping
effect caused by such granular materials by only considering the vibration analysis of plates, can berealized by combining that simulation method of bending-vibration of plates and particle motion of
the granular materials. While various kinds of researches related to vibration damping have beenstudied, simulation method of such phenomena has not been sufficiently developed. Therefore, it iseffective to analyze the motion of particles in addition to the vibration characteristics of the bending
vibration of plate-like structures and to efficiently couple them.
There are various kinds of methods currently available for vibration analysis of main oscillators, andthese methods are classified into energy-based, and wave-based methods. Methods based on energy
include practical methods [2, 3] using empirical formulas and statistical energy analysis (SEA) [4],which are used as predictions in aerospace engineering [5]and building fields respectively[6, 7]. While
SEA can be applied without limitation to high frequencies, it is difficult to obtain accuracy at relatively
low frequencies. To address such problems, a hybrid method combining SEA and wave finite element
method (FEM) has been applied to the problem of structure propagation noise by Cotoni et al., [8].
On the other hand, methods based the wave are FEM, BEM, FDM, etc. Since these wave numericalanalyses include phase information, they make it possible to predict vibration characteristics at low
127515070@ed.tus.ac.jp
t_asakura@rs.tus.ac.jp4233

and medium frequencies like building vibrations. On the other hand, the FDTD method is a method
firstly proposed by Yee [9], and it has numerical advantages that can easily visualize the wave
propagation mechanism, and applied to acoustics and vibration fields [10, 11].
On the other hand, to simulate such a motion of granular bodies, DEM (Discrete element method)has been used. DEM was devised by Cundall and Strack [12] to analyze the complex behavior of soil
in rock engineering, and it is a method that is still in development[13]. Moreover, analysis using DEM
is also carried out in research in a wide field such as research of dampers with particles in the field of
mechanical engineering[14].
Therefore, in this research, the simulation method which couple the DEM and FDTD method is
proposed, and the damping effect caused by the granular material were numerically studied. the detailof the proposed method and validation method were firstly described. Next, the damping effect of the
granular material was discussed through the numerical and experimental results. Time-Domain
Difference Analysis and Discrete Element Analysis

2. Time-Domain Difference Analysis and Discrete Element Analysis
2.1 Time Domain Difference Analysis
This study investigates the time domain calculation for bending-vibration analysis of thin plate. The
FDTD method is a method of simulating vibration propagation by sequentially performing the time
integration on orthogonal meshes. Vibration analysis of thin plate is performed by solving propagation
equation of bending wave for thin plate existing in x,y plane by finite–differential scheme.
2.1.1 Basic relational expression
Equation (1) shows the governing equation for bending wave propagation.
2
2
§ w2
w2 ·
w § w2w2 ·ww
w2w
 Uh 2 q
D ¨ 2  2 ¸ w  [ D ¨ 2  2 ¸ w  U hP(1)
wt © wx wy ¹
wt
wt
© wx wy ¹
Here, out-of-plane displacement of bending wave vibration w, density ߩ, thickness h, flexural rigidityDሺ‫ ܦ‬ൌ ‫݄ܧ‬ଷ Ȁͳʹሺͳ െ ߭ ଶ ሻሻ, loss factor of materialߤandߦ, external force q, respectively.
In the FDTD method, the differential term of the governing equation represented by the differential
equation is changed into a finite-differential equations, and solve them step-by-step calculation.
2.2 Discrete Element Method Analysis
DEM is a simulation that creates a numerical model consisting of many particles and simulates thedisplacement and rotation of each particle. The contact force is calculated using the contact model
between particles, and the acceleration, velocity, and displacement of the particles are calculated by
the difference method from the obtained contact force. In this study, contact force between particles
is divided into normal direction, and tangential direction as shown in Fig. 2, were three-dimensionallycalculated.
2.2.1 Basic relational expression
The acting force‫ܨ‬௡ in the normal direction and the acting force‫ܨ‬௧ in the tangential direction shownin Fig. 1 are expressed by the following equation.

Figure 1 ̿ Numerical model of the contact between particles, in(a) normal, and (b) tangential direction.
(2)
Fn K n  G n  Cn  G n

Ft

K t  G t  Ct  G t(3)Here, the displacement increment ߜ, the velocity incrementߜሶ , the spring constant K, and the viscosity4234

constant C in each direction, respectively. Since the Hertz contact model is adopted in this study, theparticles are in a non-linear contact condition, and the spring constant ‫ܭ‬௡ ǡ ‫ܭ‬௧  and viscosity
constant‫ܥ‬௡ ǡ ‫ܥ‬௧ in each direction are modeled as the following equation.
1
2G 2r 12
E 2r
2, K
Gt
Kn
G
(4, 5)tn
2 Q
3 1 Q 2
1
1(6, 7)
D m  K n  G n 4 , Ct D m  K n  G n 4Here, Young's modulus E, shear modulus G, Poisson's ratioߥ, radius r, and damping coefficientȽ ,respectively. The damping coefficient can be obtained from the coefficient of restitution. The
acceleration, velocity, and displacement for the acting force between particles expressed by theseequations can be obtained by the finite-differential equations composed based on the abovementioned
numerical conditions.

Cn

2.3 Coupling of time-domain difference analysis and discrete element analysis
In order to simulate the damping effect of granular material, considerarion of the influence ofgranular material on plate vibration, and that of plate vibration on granular material are required.
Hence, this study was performed by coupling DEM and FDTD method mentioned above. Therefore,in FDTD analysis for bending vibration of thin plate, the contact force caused by the collision withparticles and the effect of out-of-plane displacement in DEM analysis were considered. These two
considerations are mutually interfered with each other.
The plate vibration was simulated by considering the influence of out-of-plane displacement in
DEM analysis. The contact force݂௣ between the thin plate and the particles is obtained from Eqs. (2)and (3) as well as the contact forces between particles. Eqs. (8) to (11) were used for estimation of
the spring constant‫ܭ‬௡ ǡ ‫ܭ‬௧ and viscosity constant‫ܥ‬௡ ǡ ‫ܥ‬௧ .
8 rG p 1 2
4  E p Ew r1
Gt
Kn
G n 2 , Kt
(8, 9)
22
2 Q p
3 E p 1  Q w  Ew 1  Q p

Cn

1

1

D m  K n  G n 4 , Ct D m  K n  G n 4(10, 11)Here, displacement incrementߜ, displacement incrementߜሶ, spring constant K, viscosity constant C,
Young's modulus E, particle mass m, particle radius r, damping coefficientȽ , Poisson's ratio ߥ ,respectively. The subscripts p and w mean particles and wall (plate).
The increment of displacement can also be calculated from the difference between the incrementsof the out-of-plane displacement and the particle displacement, which are expressed as a weighting
coefficient from the four plate elements surrounding the colliding particle center in the x, y coordinatesof the particle center. The inclination of the plate at the contact point was also determined. Moreover,particle collision is considered in the analysis of plate vibration using the FDTD method by adding
the contact force݂௣ with the particle as the external force term q in Eq. (1). At this time,݂௣ means thez-direction component of the contact force between the particle and the plate described in the section
2.3. This contact force was also similarly distributed using four weighting coefficients around thecontact point. In this way, by considering the plate vibration and the particle collision, the force
considering the inclination of each element is given to the particle, and the z-direction component that
influences the plate vibration among the forces is added to the equation of bending wave propagation
as the external force.

3. Numerical case study
3.1 Calculation model
Figures 2 and 3 show models for simulating the damping characteristics of plate vibration by
granular materials. As shown in the figure, place a spherical particle with a radius of 4 mm on a flat
plate of 500 mm × 500 mm × 10 mm with all of the four sides and corners with free edge, and thecorners are supported by spring elements in the z-direction. A virtual box of 100 mm × 100 mm × 100
mm was set at the center of the flat plate to prevent the particles from spreading around by thevibration. The vibration characteristics of the plate were evaluated by giving a transient response fromthe negative to positive direction in the z-direction at the center of the plate as the excitation force

4235

and calculating the accelerance and acceleration of the center of the plate at that time. In this study,
the values required to the simulation are given as shown in Table 1. The physical phenomena wassimulated by using these parameters.
Table 1㸫Property value about analysis object
Particle

Plate

Young’s modulus, E

17.2

4.97

Density,ߩ[kg/m³]

11.337

1.150

Poisson's ratio,ߥ[-]

0.44

0.388

Shear modulus, G [GPa]5.6

Damping

0.912

Radius, r [m]0.004

Thickness, h [m]Figure 2㸫Simulated model.

0.799
0.01

Figure 3㸫Simulated model (overhead view).

3.2 Numerical results
To evaluate the influence of granular material on thin plate vibration, the main purpose was set toconfirm the difference between the vibrational characteristics, when granular material was applied.
Firstly, we compared time response of the acceleration with or without granular materials. Secondly,
the accelerance is compared with the case when each excitation force (20N, 50N, 80N, 120N, 200N)
is applied granular materials, and the influence by the excitation force is evaluated. The obtained timeresponses of acceleration with and without the granular materials are shown in Fig. 4. As can be seenfrom these results, it is confirmed that the peak of the acceleration was reduced by applying the
granular materials. The spike-like waveforms were also observed at about 0.015s to 0.04s. It may be
due to the jump-up phenomena of the granular materials acted by the external force of the hammering.

Figure 4㸫Timeme transient responses of the acceleration.
acceleration (Excitation fforce is 100N.)
In order to evaluate the influence of the excitation force, excitation was performed with each
excitation force, and the frequency response of the acceleration was confirmed. Firstly, Fig.5 shows
the accelerance when excitation force of 10N, 50N, and 100N was applied without particles. As can

4236

be seen from these results, it can be confirmed that the graphs match all three excitation forces. From
the fact that the accelerance was obtained by dividing acceleration by excitation force was the same,
it can be confirmed that the acceleration has linearity with respect to the excitation force in the caseof the vibration of the flat plate without granular materials.
Next, Fig.6 shows the accelerance calculated by applying granular material and giving each
excitation force (20N, 50N, 80N, 120N, 200N). From this result, it was confirmed that the acceleration
has nonlinearity with respect to the excitation force in the vibration of the flat plate with particles
from the fact that the graphs have difference about each excitation force. In addition, reduction of theresonance frequency at the resonance of the (1, 1) bending mode and reduction of the accelerance canbe confirmed by the addition of the granular material. When the excitation force is large, the resonanceof the (5, 5) and (7, 7) bending mode is seen to approach the value of the accelerance without particles.
In addition, when the excitation force is small, the resonance of the (1, 1) bending mode is seen to
approach that. Thus, there was no consistent tendency that the damping effect was larger as theexcitation force was larger or smaller.

Figureigure
ure 55㸫Frequency
Frequency responses of the accelerance
accelerance. (No Particle).
ParticP

Figure 66㸫Calculated
Calculated frequency
f quency responses of the accelera
freaccelerance.
Number of particles is 80.

4. Experimental investigation for acrylic scale model
4.1 Outline of experiment
The validity of the coupled analysis method proposed was confirmed by verification experiments
using an acrylic vibration model.

4237

The model experimental apparatus is shown in Figs. 7 and 8. A flat plate with size of 500 mm ×
500 mm × 10 mm, whose four corners were supported by hanging with rubber, was adopted. The
particles were spherical with a radius of 4 mm. In order to prevent the particles from spreading around,they were surrounded by a lightweight polyethylene buffer material. An accelerati on pickup was
attached to the center of the flat plate, and the center of the lower surface of the flat plate was excitedby an impact hammer.

Figure 7㸫Model experimenting apparatus.

Figure 8㸫Model experimenting apparatus (Size).
The excitation force was applied with the situation in which the particles are evenly distributed as
the initial condition, however it is difficult to evenly arrange the particles in the experiment. Therefore,
the following measurement method was adopted. With the particles randomly placed in the area
partitioned by the buffer material, the plate was vibrated and the excitation force and acceleration atthat time were measured. This was carried out several times, and the average value of the obtained
plurality of accelerances was calculated and used as the experimental result s.
In addition, in order to avoid the experimental result largely depending on the arrangement of theparticles, in the experiment, the plate is shaken to arrange the particles as uniformly as possible ateach vibration. A large number of excitation attempts were made to be close to the excitati on
forces(20N, 50N, 80N 120N, 200N). Ten pieces of data as close as possible to each of the five typesof excitation force were extracted from each of the obtained large number of experimental results, andthe average value of these was evaluated.

4.2 Results and Discussion
In order to confirm that the excitation force characteristics are identical between simulation andexperiment, the time-response of the excitation force is shown in Fig. 9. From this result, it can beseen that the experiment and simulation are in good agreement. Therefore, it is possible to compare
equally. Next, the time response of acceleration when applying an excita tion force of 100 N with andwithout the granular materials were confirmed in Figs. 10 and 11. The numerical results and theexperimental results were in good agreement without the addition of granular materials. On the other
hand, in the state where the granular material was applied, it agreed well up to 0.01 seconds. Howevera difference between them was occurring after that time.

4238

Fii 4 responses
7 Tii of the excitation force.
Figure 9㸫Time transient

Figure 10㸫Time transient responses of theacceleration. (Excitation force is100N. )(No Particle).

Figure 11㸫Time transient responses of
the acceleration. (Excitation force is100N. )(Number of particles is 80.).

Next, experimental results and numerical results are compared for the accelerance when applying
granular material and applying each excitation force (20N, 50N, 80N, 120N, 200N). The experimental
results are shown in Fig. 12 and the numerical results are shown in Fig. 6. As can be seen from theseresults, it was confirmed that the following tendencies were seen in the numerical results. The
excitation force dependency and the decrease of each bending mode frequency in the state where thegranular material is applied, even in the experimental results. In addition, in the experimental results,
the frequency characteristics around the (3, 3) bending mode under the condition of the excitation
force of 20 N are different from other excitation forces, however the same tendency is also seen in theanalysis results. On the other hand, in the experimental results of the excitation force of 80 N, unl ike
the analysis results, no sharp peak was found at the (1, 1) bending mode frequency.

Figure 12㸫Measured frequency responses of the accelerance.

4239

5. CONCLUSIONS
A simulation method combining DEM for the motion analysis of granular materials and FDTD forthe bending vibration analysis of plates was proposed. Through a numerical study, the influence ofthe granular material on the vibration characteristics of the plate was investigated. From the case study,
it was suggested that the damping effect by the granular materials on plate vibration can be predictedby using DEM-FDTD. The validity of the proposed numerical scheme is confirmed by comparison
with the verification experimental results by acrylic model.

REFERENCES
1. T. Yamauchi, S. Kawamata, "Reduction of Heavy-weight Floor Impact Sound by Granular Materials on
Ceiling," INCE. 2018. p. 5033-5040.
2. A. C. Nilsson, "Attenuation of structure-borne sound in superstructures on ships," J. Sound Vib55(1).
1977. p. 71–91.
3. Y. Matsuda, H. Tachibana, and K. Ishii, "Propagation characteristics of solid-borne sound in building
structures," J. Acoust. Soc. Jpn. (E), 35(11). 1979. p. 609–615.
4. R. H. Lyon, "Statistical Energy Analysis of Dynamical Systems,"MIT Press. Cambridge. 1975.
5. J. E. Manning, R. H. Lyon, and T. D. Sharton, “Transmission of sound and vibration to a shroudenclosed spacecraft,” NASA Technical Report, CR-81688. 1966.
6. M. D. McCollum and J. M. Cuschieri, “Bending and in-plane wave transmission in thick connected
plates using statistical energy analysis,” J. Acoust. Soc. Am. 88(3). 1990. p. 1480–1485.
7. H. Furukawa, K. Fujiwara, Y. Ando, and Z. Maekawa, “Analysis of the structure-borne sound in an existingbuilding by the SEA method,” Appl. Acoust. 29. 1990. p. 255–271.
8. V. Cotoni, P. Shorter, and R. Langley, “Numerical and experimental validation of a hybrid finite elementstatistical energy analysis method,” J. Acoust. Soc. Am. 122(1). 2007. p. 259–270.
9. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations inisotropic media,” IEEE Trans. Antennas Propagt. AP-14. 1966. p. 302–307.
10. D. Botteldooren, “Finite-difference time-domain simulation of low-frequency room acoustic problems,”
J. Acoust. Soc. Am. 98(6). 1995. p. 3302–3308.
11. J. LoVetri, D. Mardare, and G. Soulodre, “Modeling of the seat dip effect using the finite-difference timedomain method,” J. Acoust. Soc. Am., 100. 1996. p. 2204–2212.
12. P. A. Cundall, O. D. Strack, "A discrete numerical model for granular assemble," Geotechnique 29. 1979.
p. 47-65.
13. E. Oñate, F. Zárate, J. Miquel, M. Santasusana, M. A. Celigueta, F. Arrufat, R. Gandikota, K. Valiullin,
L. Ring, "A local constitutive model for the discrete element method. Application to geomaterials andconcrete," Computational Particle Mechanics, Volume 2, issue2; 2015. p. 139–160.
14. M. Saeki, "Impact damping with granular materials in a horizontally vibrating system," J. Sound Vib
Volume251, Issue1. 2002. p. 153-161.

4240