Method of estimation of frequency spectrum and power of the soundgenerated by an unsteady flow through a sonic crystal
Viktor HRUŠKA(1) , Michal BEDNAŘÍK(1) , Milan ČERVENKA(1)(1) Czech Technical University in Prague, Faculty of Electrical Engineering, Technická 2, 166 27 Prague 6, Czech Republic,viktor.hruska@fel.cvut.cz

Abstract
Although the importance of sonic crystals serving as sound barriers is growing, only a little interest was given to
the sound emission from a sonic crystal due to windy weather conditions. The Curle’s aeroacoustic analogy isreviewed in brief to obtain a suitable formulation for the low Mach number flow and rigid crystal structure. The
radiation from the crystal has the dipolar characteristic and the dominating frequency corresponds to the Stouhal
number 0.3. The values slightly exceeding 50 dB[A] were found for the radiation maxima at 20 m distance. It
follows from the method design that it provides maximum estimate of the radiated power.
Keywords: sonic crystals, vortex sound, CFD simulations

INTRODUCTION
Special configurations of specifically shaped elements called sonic crystals are promising in the field of noise
reduction. It can be shown that there are regions of frequency domain (band gaps) in which the crystal iseffectively impermeable for sound waves (see e.g. [1] for further commentary). The major practical advantageof the sonic crystals is ease of their construction compared to noise barriers demanding solid foundations.
Although there are various element shapes that can form a crystal lattice cell we limit ourselves to a simple yet
common case of the crystal consisting of cylindrical rods.
In this work we investigate possible drawback of this approach. When the sonic crystal is situated outdoors thewind flowing past the rods might act as a source of sound. The topic is approached by numerical simulationsof the turbulent flow with subsequent analysis by means of the vortex sound theory. The presented method
calculates the maximum estimate of the radiated sound power.

THEORYSound generation from the airflow through a sonic crystal is inherently linked to the well-known phenomenonof flow instabilities past a cylinder. The Strouhal number St specific for the studied case is not predetermined,
but we may assume it to lie between 0.2 – 0.4. Hence, we can assess that for rod diameters in centimeters andflow speeds around 10 m/s (fresh breeze) the flow instabilities frequency should be expected in tens of hertzs.
The corresponding wavelength of the radiated sound is at least three meters, making not only the rods’ diameterbut also the height and the breadth of the crystal acoustically compact for many practical cases. It follows thatthe airflow through the crystal shall be regarded incompressible. For the same reasons we do not consider
any interaction of the flow induced sound with the crystal stop-band. Usually they are at least one order ofmagnitude apart from each other for the outdoor situated noise barriers. To sum up these considerations, we
can simulate the incompressible turbulent flow through the sonic crystal and then employ a suitable aeroacoustic
analogy.
In his well-known pair of articles [2, 3] Lighthill proposed equivalency of the sound generated by the turbulencein the real fluid and the sound produced in the ideal and stationary medium with a source governed by the socalled Lighthill stress. The resulting wave equation for density perturbations ρ 0 = ρ − ρ0 propagating at theadiabatic sound speed c0 in an unbounded Stokesian fluid reads [5]3179



∂2 0
ρ − c20 ∇2 ρ 0 = ∇ · ∇ · ρuuu + p0 − c20 ρ 0 I − σ ,(1)2∂twhere u, p, I and σ denote the fluid velocity, pressure, the unit tensor and the viscous stress tensor respectively.
The mean values are denoted with the subscript 0, the perturbation values with the prime (·)0 (e.g. p0 = p − p0 ).
In the limit of low Mach number Ma and high Reynolds number Re the pressure and density perturbations
on the right hand side cancel each other, the viscous stresses are negligible compared to the Reynolds stresses
ρ0 u u . That leaves only the double divergence ∇ · [∇ · (ρ0 u u )] on the right hand side of Eq. (1).
In the presence of stationary solid body within the flow dipole sources emerge beside the free-field quadrupolar
contributions. They are formally introduced by the so called Curle’s aeroacoustic analogy resulting from theLighthill’s theory taking into account the source imaging due to the solid boundary presence. Let H denote theHeaviside step function defined as 0 for a point inside the body and 1 everywhere else. Let us suppose the low
Mach and the high Reynolds number limit and the rigid body. The wave equation for the density perturbations
then takes the form [5, 4]



∂2
Hρ 0 − c20 ∇2 Hρ 0 = ∇ · [∇ · (Hρ0 u u )] − ∇ · p0 I + ρ0 u u · ∇H .
(2)
2∂tThe first source term corresponds to the previously investigated quadrupole source from the turbulent wake andits power is proportional to Ma5 . The second source term represents the reaction of the body surface to theunsteady flow. Its radiation characteristic is dipolar and its power proportional to Ma3 for low Mach number[5]. Hence, the qudrupole contribution is negligible compared to the dipole one when Ma → 0. For a rigid bodywith the no slip velocity condition on its surface ρ0 u u = 0 on the face of the body and therefore only the p0 I
term on the right hand side persists.
The following integral analogue of Eq. 2 is obtained by means of the three-dimensional free-space Green’sfunction:

I p0 t − |xx−yy|c0∂Hc20 ρ 0 = −ni dS(yy) ,∂ xi S 4π|xx − y |(3)where ni is the ith component of the unit normal vector on the rod surface. In the far field (Fraunhofer)
approximation (|x| → ∞, p0 = c20 ρ 0 ) it reads
p0 (xx,t) ≈xi ∂c0 |xx|2 ∂t|xx| x · yp0 t −+ni dS(yy) .
c0 c0 |xx|SI(4)It follows that the far field acoustic pressure can be calculated from the known unsteady pressure on the rods’
faces obtained by a suitable numerical simulation. Since the rods’ diameter is small compared to the typical
wavelengths we disregard any diffraction effects. The integral in Eq. (4) represents the force exerted on thesurface of the rod by the fluid. For an isolated cylinder the net force is usually predominantly perpendicular to
the free flow direction and therefore it is sometimes called lift (force). Because the relation of the complicatedflow inside the crystal to the free flow is nontrivial we don’t use this designation. However, as it will be shown
below, the radiated sound directivity suggests that the overall effect of the unsteady forces is "rather a lift thana drag".

NUMERICAL EXPERIMENTS
The numerical results presented below was simulated using Comsol Multiphysics software, version 5.4. The
fluid dynamics simulations employed the Single-phase Turbulent flow interface with the Shear Stress Transport

3180

Figure 1. Example of a monofrequency plane wave (1040 Hz) propagating through the studied sonic crystal(left), transmission loss of the crystal (right) in the direction depicted left.
model (i.e. the Unsteady Reynolds Averaged Navier Stokes Equations – URANS – approach, see e.g. [6]).
The inflow velocity was 12 m/s (corresponding to a fresh breeze), the open boundary condition is imposed on
all other boundaries. The auxilliary acoustic simulations were made within the Pressure acoustics environment.
Subsequent treatment of the simulated data was made using our own code in GNU Octave environment.
3.1 The studied caseIn order to bring the simulations close to a real example we first review the features the sonic crystal exhibits
as a sound barrier. Assuming that the outdoor placed crystals shall prevent primarily the traffic noise it isreasonable to tune the configuration to have a band-stop around 1 kHz. In the following a hexagonal latticewith the lattice constant 0.13 m and the rods’ diameter 60 mm is studied. The transmission loss (assuming
invsicid and thermally nonconducting air without mean flow) and an example of wave propagating through thecrystal are in Fig. 1.
3.2 Restriction to two dimensions
Although the turbulence is essentially a three dimensional phenomenon, we might make use of the fact that thevariation of pressures and velocities is much more pronounced across the rods then along them. A simplified
example is depicted in Fig. 2 (the symmetry boundary condition is imposed on the top and the bottom of thedomain). The isosurfaces of unsteady pressure are virtually invariant along the z direction. Naturally, we cannot
fully exclude the surface pressure perturbations along the height of the rods, especially further inside the crystal,
we just assume them to be significantly weaker than the azimuthal perturbations.
On the other hand, a complete restriction to the planar description (e.g. using the 2D instead of the 3D Green’s
function in Eq. (4)) would lead to a considerable overestimation of the radiated sound power. Therefore, wesuggest to extrude the planar flow solution to three dimensions and employ Eq. (4) on the rods’s surfaces
disregarding the bases.
Weak phase differences along the height of the rod (which would not be captured by the suggested approach)are possible. Nevertheless, their near field effects would lead to weaker far field radiation and so we are getting
the maximal radiated power estimation.

3181

Figure 2. Instantaneous isosurfaces of pressure petrurbations p0 = p − p0 , where p0 (xx) = 1/T 0T p(xx)dt with T
being a period of steady oscillations. The colorbar values are in pascals.
R3.3 Sound power and spectrum based on extrusion
A planar set-up consisting of 18 rods in 6 columns was simulated. A snapshot of the instantaneous velocity
field is depicted in Fig. 3. The von Kármán vortex street is not formed inside the crystal. However, the trailingwakes behind the cylinders are unstable and therefore pressure perturbations of the cylinders’ faces occur.
Eq. (4) is employed in a discretized form. The radiated pressure is evaluated on a circle with the center in themidpoint of the crystal and of the 20 m radius. The radiated sound spectrum is dominated by the frequency
62.5 Hz (see Fig. 4), corresponding to the Strouhal number of 0.31 (the rods’ diameter is the characteristiclength).
In Fig. (4) the sound pressure level (in dB[A]) of the dominant frequency is depicted as the height and thecolor is governed by the phase angle. It is evident that the sound radiation from the crystal resembles a dipolewith nulls at the inflow direction.

DISCUSSION
Limitations of the presented method depends on the validity of various assumptions. First, the URANS approach neglects some of the turbulent nonlinearities propagating through the crystal. The main mechanism iscaptured correctly but some subtleties regarding the perturbations along the rods’ height and fine tuning of theflow synchronization might have been left out. A logical step in the future research is therefore a comparison
between the URANS and the Large Eddy Simulation (LES) in a benchmark case. Second, the imperfectionsof the extrusion technique should be better quantified and further tested. Another point of high interest is thestability of Strouhal number with varying inflow speed and crystal geometry.
The Curle’s formulation under conditions allowing for the volume sources to be neglected is well-suited for theBoundary Element Method (BEM), which might be employed instead of the Green’s functions.
The calculations disregards the reflections from the ground. Although many of the surface types (e.g. grass)have low reflection coefficient for almost grazing incidence, which would be the case, this topic should beclosely discussed in the future.

3182

Figure 3. Instantaneous velocity field around the crystal center visualized by the Line Integral Convolution.
Phase angle

4504030−20220

Flow direction−100m

20−20−104020

10010SPL (dB[A])SPL (dB[A])60m00200400f (Hz)600800

Figure 4. Left: Sound pressure level (SPL) and phase angle of the dominant frequency evaluated on the circleof 20 m radius surrounding the sonic crystal. Right: Spectrum of the emmited sound in the main radiation
direction.
Note that for the noise hygiene purposes the effect of low frequency on the human perception has to be countedin (e.g. by conversion of decibels to phones). Generally, the low frequency noise is not percieved as loud as
the benchmarks around 1 kHz (cf. e.g. [7]).

CONCLUSIONS
The method of estimating the radiated spectrum and sound power generated from the flow through a sonic
crystal was designed, tested and discussed. The sonic crystal under windy weather conditions radiates the soundwith the dipolar characteristic and the fundamental frequency corresponding to the Strouhal number 0.3. Sincesome of the loss mechanisms were disregarded the method provides the maximum estimate of the radiatedpower (ca. 55 dB[A] in the main radiation direction).
The future goals include benchmarking of the method with the LES, study of the Strouhal number stability andcomparison of the results obtained employing the Green’s functions with the BEM simulations.

3183

ACKNOWLEDGEMENTS
This work was supported by the GACR grant No. 18-24954S and by the Grant Agency of the Czech Technical
University in Prague, grant No. SGS19/107/OHK3/2T/13.

REFERENCES[1] Ming-Hui Lu, Liang Feng, and Yan-Feng Chen. Phononic crystals and acoustic metamaterials. Materials
Today, 12(12):34–42, dec 2009.
[2] M. J. Lighthill. On sound generated aerodynamically. i. general theory. Proceedings of the Royal Society ofLondon. Series A, Mathematical and Physical Sciences, 211(1107):564–587, 1952.
[3] M. J. Lighthill. On sound generated aerodynamically. ii. turbulence as source of sound. Proceedings of theRoyal Society of London. Series A, Mathematical and Physical Sciences, 222(1148):1–32, 1954.
[4] N. Curle. The influence of solid boundaries upon aerodynamic sound. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 231(1187):505–514, sep 1955.
[5] M. S. Howe. Theory of Vortex Sound. Cambridge University Press, Cambridge, 2002.
[6] Mehrzad Tabatabaian. CFD module. Mercury Learning and Information„ Dulles, Virginia :, 2015. Mehrzad
Tabatabaian., Includes bibliographical references (pages 191-197) and index.
[7] William M. Hartmann. Signals, Sound, and Sensation (Modern Acoustics and Signal Processing). American
Institute of Physics, 2004.

3184