Finite element analysis on the surface characteristics of acoustic resonatorswith thermal and viscous boundary layers
Naohisa INOUE∗ , Tetsuya SAKUMA
Graduate School of Frontier Sciences, The University of Tokyo, Japan

Abstract
The thermal and viscous boundary layers (TVBL) around a resonator’s neck play important roles in the absorption, particularly when there is no porous material near the neck. This paper presents a simple implementationof a TVBL boundary condition for the three-dimensional finite element analysis of Helmholtz equation. The
implementation is verified by comparing normal incidence absorption coefficients calculated and measured forslit-type resonators. Good agreement can be observed under the appropriate calculation setting. Subsequently, adetailed investigation is conducted on surface impedances of Helmholtz resonators composed of circular holes
and square baffles. Numerical values are compared with analytical values based on the classical models such as
Maa, Allard–Ingard and so on. It is demonstrated that the presented method convincingly captures the surface
impedance on an acoustic resonator unless the TVBL thicknesses are much larger than the perforation diameter.
Keywords: Viscous and thermal boundary layer, Finite element method, Helmholtz resonator

1

INTRODUCTION

Predicting the absorption performance of resonators is of great interest for noise control engineers. The
thermal and viscous boundary layers (TVBL) around a resonator’s neck play important roles in the absorption,
particularly when there is no porous material near the neck. Several papers have been dedicated to the numerical prediction of resonator’s absorption performance. The recently proposed approaches have been based oncomputational fluid dynamics [1] or the viscothermal acoustic elements [2]. Accordingly, relatively high computation and implementation costs are required. This paper concerns another simple approach that is based on theequivalent boundary condition of TVBL in the conventional Helmholtz equation [3, 4]. Despite its effectiveness,the details of the implementation have not been presented, especially where giving the TVBL boundary condition. Furthermore, to the authors’ best knowledge, three-dimensional analysis has not been performed yet. Thismay be caused by the apparent difficulty of evaluating the tangential Lhaplacian term for the viscous effect on
the boundary.
This paper presents the finite element implementation of the TVBL boundary condition for three-dimensionalfrequency-domain problems. This implementation is applied to analyze the absorption characteristics of slit-type
and Helmholtz resonators, and the appropriate calculation setting is clarified.

2

FINITE ELEMENT IMPLEMENTATION

2.1 TVBL Boundary Condition
The time convention is e jω t throughout this paper. The boundary condition taking the effect of thermal andviscous boundary layers into account is given in [5] as
[
(
)]
∂p
1
Rt
Rv 2
−
= − jωρ0+ (1 + j)
∇
p
(1)
∂n
Zn
ρ02 c20 ρ02 ω 2 tan
where p, ρ0 , c0 , and Zn are the sound pressure, air density, sound speed, and normal-incidence surface impedance,
respectively. ∇2tan is the Laplacian operator in the tangential plane of the boundary and ∂ /∂ n is the normal∗ inoue@env-acoust.k.u-tokyo.ac.jp

2005derivative to the plane. Rt and Rv are the thermal and viscous resistances, which are related to the thicknessesof the thermal and viscous boundary layers, δ ′ and δ , as follows:
√
√
1
2µ
2µ
1
′
′
Rt = ρ0 ω (γ − 1)δ , δ =and Rv = ρ0 ωδ , δ =,(2)2
ρ0 ω Npr2ρ0 ω
where γ , µ , and Npr are the heat capacity ratio, viscous coefficient, and Prandtl number of the air, respectively.
Note that several signs in the above are modified from those in [5] owing to the time convention e jω t . To
simplify Eq. (1), it is rewritten as∂p= − jkβ p − ( jω )3/2 cT p − ( jω )−1/2 cV ∇2tan p∂n(3)where β is the normalized admittance and the frequency-independent coefficients, cT and cV , are defined as
√
√
γ −1
µ
µ
cT = 2and cV =.
(4)
ρ0 Nprρ0
c0
Obviously, the TVBL boundary condition involves a fractional derivative in the time domain expression, whichmakes it difficult to use in time domain analysis. On the other hand, this paper focuses on the treatment of theboundary condition in frequency domain analysis.
2.2 Weak Formulation of the Helmholtz Equation
The weak form of the Helmholtz equation, ∇2 p + k2 p = 0, is given as∫Ω∇δ p · ∇pdV − k∫2Ωδ p · pdV −

∫
Γ

δp

∂p
dS = 0,∂n(5)where δ p is the admissible variation and ∇ is the nabla operator. Substituting Eq. (2) into Eq. (3), thefollowing equation is obtained:∫Ω∇δ p · ∇pdV − k2∫Ω

δ p · pdV + jkβ∫Γδ p · pdS + ( jω )3/2 cT∫Γδ p · pdS + ( jω )−1/2∫
Γδ p · ∇2tan pdS = 0.

(6)

The main difficulty of the implementation is to evaluate the fifth boundary integral term. Although theinterpolation function must be differentiable twice or more to evaluate the term directly, this requirement canbe reduced by applying Green’s first identity in two-dimensions as follows:∫Γδ p · ∇2tan pdS = −∫Γ∇tan δ p · ∇tan pdS +∫∂Γδ p · (∇tan p · nc )dC(7)

where ∇tan is the gradient operator in the tangential plane and nc is the outward normal vector to the contourof the boundary surface. Accordingly, the energy dissipation along the contour should be evaluated in the latter
term, which makes the global system matrix asymmetric.
2.3 Evaluation of the Tangential Gradient Term
This section presents an algorithm to evaluate the tangential gradient term in the generalized curvilinear
local coordinate set in the boundary elements. In the first, let us consider the local coordinate (ξ , η ) thatis usually set for the numerical integration in the finite element method. Then the global coordinate values
in an element are interpolated with the shape functions Ni (ξ , η ) and the nodal coordinate values (xi , yi , zi ) as
x = ∑i Ni (ξ , η )xi and so on. The covariant basis vectors of the tangential plane, eξ and eη , are given as follows.
eξ (ξ , η ) =2006∂ x(ξ , η )∂ y(ξ , η )∂ z(ξ , η )
ex +
ey +
ez ,
∂ξ
∂ξ
∂ξ

eη (ξ , η ) =∂ x(ξ , η )∂ y(ξ , η )∂ z(ξ , η )
ex +
ey +
ez
∂η
∂η
∂η

(8)

where ex , ey , ez are the unit basis vector of the global Cartesian coordinate. On the other hand, the contravariant
basis vectors of the tangential plane, eξ and eη , are given as follows:eξ =eη × eζeζ × eξ, eη =.
eξ · (eη × eζ )eη · (eζ × eξ )(9)where eζ is the normal vector to the tangential plane defined as eζ = eξ × eη . The tangential gradient vectors
in Eq.(4) can be expressed by using the terms of the general curvilinear local coordinate as follows.
∇tan δ p =∂δ p ξ ∂δ p η∂p ξ ∂p ηe +e , ∇tan p =
e +
e
∂ξ
∂η
∂ξ
∂η

Finally, the first term of Eq.(4) can be rewritten as
)
]
(
∫
∫ ∫ [
∂δ p ∂ p ξ ξ∂δ p ∂ p ∂δ p ∂ p ξ η ∂δ p ∂ p η ηe ·e +∇tan δ p · ∇tan pdS =e ·e ++e · e |J|d ξ d η
∂ξ ∂ξ
∂ξ ∂η
∂η ∂ξ
∂η ∂η
Γ
η ξ

(10)

(11)

where |J| is the Jacobian at each point in an element and defined as |J| = |eξ × eη |. The integrand in Eq.(11)can be easily evaluated at the gauss quadrature points. Furthermore, the second term of Eq.(4) can be rewritten
as
(
)
∫
∫
∂p ξ
∂p η
δ p · (∇tan p · nc )dS =δpe · nc +e · nc |J|d ξ (, d η ).
(12)
∂ξ
∂η
∂Γ
ξ (,η )
This contour integral should be calculated for surface elements adjacent to ∂ Γ.

3

PREDICTION OF THE NORMAL-INCIDENCE ABSORPTION COEFFICIENT

3.1 Slit-type Resonators
In [3], the TVBL boundary condition was implemented in the two-dimensional boundary element method(BEM), where the term of ∇2tan p was directly evaluated using the finite differences of sound pressures between
neighboring elements. This implementation was applied to estimate the normal-incidence absorption coefficientsof slit-type resonators, and the results were confirmed to be in good agreement with measured values. In thissection, to clarify the appropriate setting of the implementation, slit-type resonators are examined by the threedimensional finite element method (FEM), and the results are compared with BEM results and the measured
values in [3].
Duct
D

CavityBD = 100 mm

Vibratory Boundary

dN
Γ

nl

N

D

∂Γ , ∂Γ

N

nl

D

Γ

N

w

C

nl ∂ΓC, ∂ΓN
C

N

Γ
Neck
Γ

D

Γ

C

BC = 50 mm

ΓDC = 50 mm

DD = 500 mm

Figure 1. Horizontal geometry of the duct with a slit-type resonator.

2007Normal−inc. Absorption Coefficient []In the FEM analysis, a 100-mm-square duct with a slit-type resonator is modeled, where all surfaces are
as-sumed to be rigid (β = 0). Figure 1 shows the horizontal geometry of the model and the notation of theboundaries. These conditions correspond to those of the measurement in [3]. As a broad setting, the TVBL
boundary condition is given to the end surface of the duct (ΓD ), the front surface of the cavity (ΓC ), andthe side, upper, and bottom surfaces of the slit neck (ΓN ). On the contours of ΓD and ΓC , ∇tan p · nc is zero
except along the edges of the neck. Thus, contour integrals are evaluated along the edges (∂ ΓD , ∂ ΓC , and ∂ ΓN ).
Second-order prism and hexahedral elements are used, and a fine mesh is applied as illustrated later in Fig.
4(a). The normal-incidence absorption coefficients are calculated by the transfer function method at 1/24-octave
band center frequencies [6].
Figure 2 shows the calculated and measured results for four slit-type resonators with different combinationsof the width and depth of the neck. The FEM results are given for three different settings, and the BEM andmeasured results are cited from [3]. The FEM results with the above broad TVBL setting are represented by
red lines. The calculated values for the 4-mm-width slit correspond relatively well with the BEM and measured
ones, whereas those for the 2-mm-width slit are slightly larger than the others. Nevertheless, this TVBL setting
for the FEM seems reasonably reliable.
Next, a narrow TVBL setting only in the neck (ΓN ) is examined. As can be seen from the results represented by blue lines, the calculated values are considerably greater than those with the broad setting. Thisoverestimation of the absorption coefficient may be caused by intensification of the edge potential around theopenings of the neck due to the lack of a tangential drag force on ΓD and ΓC .
Therefore, it is suggested that the TVBL boundary condition should be given to two surfaces along an edge.
If the asymmetric contour integral term for edges, Ic, can be neglected, the calculation becomes more efficient.
However, the absorption coefficients calculated without Ic, represented with green lines, are lower than the redlines. This demonstrates that the energy dissipation at edges around a resonator’s neck must be considered to
predict the absorption coefficient accurately.
In the above calculation, fine meshes were used because the acoustic field around the neck is expected tochange sharply. Figure 3 shows the results calculated with the fine and coarse meshes illustrated in Fig. 4.
As a fixed condition, two hexahedral elements are generated in the direction of the slit width. The maximum
lengths around the neck are 2 and 10 mm for fine and coarse meshes, respectively. It can be seen that the
1
1
0.9 w = 4 mm
0.9 w = 2 mm
0.80.8
d = 30 mm
0.7 N
0.7
0.6
0.6 dN = 30 mm
dN = 10 mm
0.50.5
dN = 10 mm
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
100 150 200 250 300 350 400 450 500 550 600 200 250 300 350 400 450 500 550 600 650 700
Frequency [Hz]Frequency [Hz]Evaluated integrals in the FEMΓΝ, ΓD, ΓC, ∂ΓΝ, ∂ΓD, ∂ΓCΓΝ, ∂ΓΝΓΝ , ΓD , ΓC

Cited from [3]BEM

Measured

Figure 2. Normal-incidence absorption coefficients of four slit-type resonators calculated by the FEM with threedifferent settings. The BEM and measured values for PVC are cited from [3].

2008

Normal−inc. Absorption Coefficient []

1
0.9
Fine
Coarse0.8
0.7
dN = 30 mm
0.60.5
dN = 10 mm
0.4
0.3
0.2
0.1
0
200 250 300 350 400 450 500 550 600 650 700
Frequency [Hz]w = 4 mm

w = 4 mmdN = 30 mm

dN = 30 mm(a) Fine mesh(b) Coarse Mesh

Normal−inc. Absorption Coefficient []Figure 4. Meshes for the examination of effects of mesh
Figure 3. Normal-incidence absorption coefficient calcu- resolution around neck.
lated with two different mesh resolutions.
(a): dd = 5 mm, dN = 5 mm, p = 25 mm, Dc = 45 mm(b): dd = 1 mm, dN = 1 mm, p = 11.2 mm, Dc = 100 mm
1
FEM
0.9
Atalla0.8
dN
pp–Sgard
Dc
0.7
Maa
dd
0.6
p
0.5
0.4
0.3
0.2
0.1
0
200 250 300 350 400 450 500 550 600 650 700200 250 300 350 400 450 500 550 600 650 700
Frequency [Hz]Frequency [Hz]Figure 5. Normal-incidence absorption coefficients calculated for two Helmholtz resonators with different sizes,
Types (a) and (b).

results for the coarse mesh are slightly greater than those for the fine mesh, but the difference is not significant
when predicting the absorption coefficient.
3.2 Helmholtz Resonators
As a typical example of three-dimensional resonators, Helmholtz resonators are analyzed under a similar setting
to that in the previous section, where the cross sections of the cavity and duct are identical and square. In the
FEM calculation, the TVBL boundary condition is given to the front and back surfaces of the perforated baffle
and the cylindrical surface of the neck. The normalized admittance β is assumed to be zero for all surfaces.
The FEM results are compared with the existing theoretical models of Atalla and Sgard [7] and Maa [8]. The
resistance and inertial correction terms are approximately derived for periodic resonators in the Atalla—Sgard
model and for a single resonator in the Maa model. Note that the resistance correction of the Maa model isbased on [7].
Figure 5 shows the calculated results for two resonators with different sizes. It can be seen that all the theoretical and FEM results are relatively correlated but slightly different in terms of the peak value and frequency.
The sound field in the duct is equivalent to the normal-incidence field over the periodic resonators. However, the2009FEM results are closer to the Atalla-–Sgard model for Type (a) but closer to the Maa model for Type (b). This
inconsistency implies that the interactive motion of the neck and cavity considerably affects the resistance andinertance corrections for the neck’s impedance. From this viewpoint, the FEM may be advantageous because
the effects of the interaction are automatically incorporated.

4

PREDICTION OF THE SURFACE IMPEDANCESThis section discusses the surface impedance of Helmholtz resonators in more detail in order to clarify theapplicable condition of the TVBL boundary condition. In this study, Maa model is employed as the theoreticalreference. The analytical effective density for the cylindrical pore, ρ ′ , shows the low- and high-frequency
asymptotic behaviors as follows.

18η 4ρ + 1 8η ,[for ω <= ωL
(13)0√ ]−1 

2
2

3
j
ω
r
ρ
0r
2 J1 (s − j′√√ρ = ρ0 1 − √≃s − j J0 (s − j
1 2ρ0 η32η
(1 − j), for ω >
(14)
= ωH
 ρ0 +rω
ρ0 r2
4.1 Numerical Models of Helmholtz Resonators
Figure 6 shows the numerical models for calculating the surface impedance by transfer function method. The
cross-section of the tube is the square with D mm side. The perforation is placed in the center of the square.
In order to better understanging, two neck models of Helmholtz resonators are employed: in model A, ρ ′ isgiven as the air density inside the resonator neck, whereas in model B, ρ0 is given and the TVBL boundary
condition is given to the peripheral surface of the neck. The TVBL boundary condition is give to the baffle
surfaces in both models. For the sake of completeness, the sound pressures averaged over the cross-section areused for the transfer function method to exclude any scattered oblique waves.
The calculation was excuted at 1/3 octave band center frequencies with the following parameters. Holeradius is set to 0.25, 0.5, 1.0 and 2.0 mm. Two combinations of neck length, L mm, and hole pitch, D mm,
are set: (a)(L, D) = (1, 10) and (b)(L, D) = (4, 20).
4.2 Results and Discussions
Figure 7 shows the surface resistance, reactance and the absorption coefficient under the normal incidence
condition. As a general tendency of the resistance, model A below ωL and model B above ωH show goodD10201 recs. per section

Vibratory

20 mm 21 mm

2r

VBC

VBC

Equivalent
Fluid

L

500 mm
50 mm
Common Setting of the Calculation

VBC

VBC

Model A

Model BTwo models around the perforation

Figure 6. Calculation settings of the surface impedances. VBC stands for viscous boundary layer boundary
condition.

2010

102
101

Neck Radius
0.25

0.5

1.0

2.0

Normalized Absolute Resistance [](a) L = 1 mm, D = 10 mm

103
102

Theoretical Values Based on Eqs. (1) to (4)1Normal−inc. Absorption Coeff. []0.9 (a) L = 1 mm,
0.8
D = 10 mm

Normalized Absolute Reactance [](a) L = 1 mm, D = 10 mm

0.7

100

101

0.6

10−1

0.5
100

10−2

0.4
0.3

10−1

10−3

0.2
0.1

10−4
62.5
104
103

125

250

500
1k
2k
Frequency [Hz]

4k

10−2
8k
62.5
104

Normalized Absolute Resistance [](b) L = 4 mm, D = 20 mm

103

102

102

101

101

100

100

10−1

10−1

125

250

500
1k
2k
Frequency [Hz]4k

08k 62.5
1

Normalized Absolute Reactance [](b) L = 4 mm, D = 20 mm

0.9

125

250

500
1k
2k
Frequency [Hz]

4k

8k

4k

8k

Normal−inc. Absorption Coeff. [](b) L = 4 mm, D = 20 mm

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1

10−2
62.5

125

250

500
1k
2k
Frequency [Hz]

4k

10−2
62.5
8k

125

250

500
1k
2k
Frequency [Hz]4k

08k 62.5

125

250

500
1k
2k
Frequency [Hz]Figure 7. Normalized absolute surface resistance and reactance and normal-incidence absorption coefficient
calculated for different combinations of neck length(L), neck radius(r) and hole pitch(D). Circle markers denotethat the value was negative. Star markers denote the lower limit of surface resistance (8η LD2 /π r4 ). Solid andchain lines are for model A and B, respectively.

correlation with the theoretical reference. The resistance of model A becomes lower than the theoretical reference in entire frequency range above ωL . On the other hand, model B does not capture the low-frequency
asymptotic behavior. Regarding the reactance, in particular for r = 0.25 mm, the model B does not reflect theapparent increase of the mass-reactance indicated as the factor 4/3 of Eq.(13), which causes the shift of theabsorption peak. These tendencies demonstrate that the reliability of the TVBL boundary condition becomes
lower below ωL for the cylindrical pores. Besides, the resistances determined by the numerical calculations
become negative in the high frequencies in several cases. A cause of this unpreferable behavior may be due tothe approximation in the derivation of the TVBL boundary condition. This point should be investigated more
carefully in future work.

5

CONCLUSION

This paper presented and verified the finite element modeling of acoustic resonators with thermal and viscous boundary layers. It was confirmed that the TVBL boundary condition should be given not only to theside surfaces in the resonator’s neck but also to the front and back surfaces of the perforated baffles. Furthermore, the energy dissipation at the neck’s edges was not negligible. Thus, an asymmetric matrix equation
should be solved. Numerical methods are advantageous for predicting the absorption performance of resonators
not only because of their capability of handling resonators with complex shapes, but also because the resistance and reactance corrections are automatically incorporated, reflecting the acoustic field around the neck.
However, it should be emphasized that the TVBL boundary condition should be applied with great care for2011some frequency ranges. In the low frequencies, the TVBL boundary condition cannot express the low-frequency
asymptotic behavior of the effective density in the pores. On the other hand, in the high frequencies, the TVBL
boundary condition employed in this paper may involve some erroneous approximation in its derivation, whichcause physically unreasonable behavior in the surface resistance. In future work, the presented results should beverified by comparison with those of more exact methods such as those in [1, 2].

REFERENCES[1] J. S. Bolton and N. Kim, “Use of CFD to calculate the dynamic resistive end and correction for microperforated materials,” Acoust. Aust., 38, 134–139 (2010).
[2] W. R. Kampinga, Y. H. Wijnant and A. de Boer, “An efficient finite element model for viscothermal acoustics,” Acta Acustica united with Acustica, 97, 618–631 (2011).
[3] K. Monjyugawa, M. Terao and H. Sekine, “Numerical analyses and high precision experiments on acoustic
resistances of slit resonators,” Proc. Spring Meet. Acoust. Soc. Jpn., 639–640 (2000) (in Japanese).
[4] S. Sakamoto, H. Mukai and H. Tachibana, “Numerical study on sound absorption characteristics of resonance-type brick/block walls,” J. Acoust. Soc. Jpn., 21, 9–15 (2000).
[5] P. M. Morse and K. U. Ingard, Theoretical Acoustics, (Princeton University Press, New Jersey, 1968), Chap.
6.
[6] N. Inoue and T. Sakuma, “Numerical investigation of effect of support conditions of poroelastic materials
in impedance tube measurement,” Acoust. Sci. & Tech., 38, 213–221 (2017).
[7] N. Atalla and F. Sgard, “Modeling of perforated plates and screens using rigid frame porous models,” J.
Sound Vib., 303, 195–208 (2007).
[8] D. Maa, “Potential of microperforated panel absorber,” J. Acoust. Soc. Am., 104, 2861–2866 (1998).

2012