Numerical simulation of nonlinear viscothermal acoustic wave
propagation
Jan Albertus DE JONG(1) , Erwin Reinder KUIPERS(2)(1) ASCEE, The Netherlands, j.a.dejong@ascee.nl(2) Sonova AG, Switzerland, erwin.kuipers@sonova.com

Abstract

Behind-The-Ear (BTE) hearing instruments are typically applied to assist people with a severe toprofound hearing loss. Accordingly, very high acoustic pressure levels can be obtained in the ear
canal of the wearer. At such levels, significant distortion is present, leading to a reduced sound quality. The distortion mainly originates from the balanced armature receiver driver, operated beyond its
linear limit. However, as the acoustic pressure level in the acoustic transmission system (basicallya set of tubes) can be much higher than 140 dB SPL, it is of interest to investigate the amount of
acoustic distortion caused by nonlinear wave propagation.
As such, this paper focuses on the prediction of the distortion due to the combined effect of viscothermal and nonlinear wave propagation. To simulate these effects in the tubes, a modified form of
the model of Chester is used. In this modified model, the viscothermal effects results in a half orderderivative of the acoustic pressure to time. This fractional derivative term is approximated using the
diffusive approximation, where attention is given to the linear wave dispersion error. The correct
implementation of the model is verified using a 2D full nonlinear compressible Computational Fluid
Dynamics simulation.
Keywords: Nonlinear, viscothermal, duct acoustics

1

INTRODUCTION

To transport the amplified signal from the receiver to the inner ear, a BTE hearing aid instrument
uses several tubes with relatively small diameter, in the order of 1 mm. Due to these small channels,
the acoustic wave propagation is subject to viscothermal losses. Moreover, for high amplitudes
nonlinear effects can result in additional distortion of the wave signal. The combined effects ofhigh-amplitude nonlinearities and viscothermal losses have the potential to deteriorate the speech
intelligibility. Therefore, for hearing aid instrument manufacturers it is of importance to understandthe combined impact of these effects. One way to get more insight into the combined effects is
by utilizing a full Computational Fluid Dynamics (CFD) simulation. This often involves a too highcomputational cost. Hence, more simplified model is required. The focus in this work is on thederivation and implementation of a one-dimensional model for nonlinear and viscothermal acoustic
wave propagation in these tubes.
The first one-dimensional (1D) nonlinear model which both incorporates viscothermal and nonlinear
effects has been developed by Chester [2]. This model is based on the boundary layer approximation, i.e. the viscous boundary layer δν [11] is small compared to the hydraulic radius of the tube.
It is not our intention to provide a complete overview of existing efforts in this subject. However,as a good starting point, we would like to refer the reader to the work of Menguy and Gilbert [9]and Sugimoto and Shimizu. [10].
In these models, the viscothermal effects return either into the governing continuity equation, or in
the governing momentum equation in the form of a fractional order derivative. A fractional orderderivative is non-local in time, and requires the full time history of the computation to evaluate.

5481

This is problematic in a numerical simulation, as the required amount of memory for the solution
will grow linearly with the simulated time. As a solution to this problem, Haddar et al. [4] proposed a “diffusive approximation” to the fractional derivative. Using this diffusive approximation,
the evaluation of the fractional derivative becomes again local in time, at the cost of the evaluationof extra degrees of freedom (DOFs) per grid point. Fortunately, the amount of extra DOFs is low.
It should be noted that for boundary layers which are thick compared to the hydraulic radius, thefractional order derivative introduces a dispersion error, as it is derived by assuming a thin boundary
layer in the first place.
In this work, we will improve upon the evaluation of the diffusive approximation to minimize thisdispersion error for small tubes, especially at low frequencies. The 1D numerical model is implemented in a general framework for Finite Element Analysis (FEA) and results are compared to afull nonlinear CFD simulation.
In the next section, we will present the 1D model and its major assumptions. Using a dispersion
analysis we will show how to obtain required quadrature coefficients of the diffusive approximation.
We will shortly discuss the implementation of the 1D model. Section 3 provides a comparison of
results for this 1D model with 2D CFD results. We will finish with some conclusions and remarks
on future work in Section 4.

2

1D MODEL

Upon applying the following assumptions to the full compressible Navier Stokes equations for acalorically perfect gas [7]:• Laminar flow,• Prismatic duct,• Isothermal walls,• Large aspect ratio (length in wave propagation direction, L is much larger than transverse
length scale (hydraulic radius, rh ),• Hydraulic radius is much larger than viscous and thermal penetration depth [11, p. 7],• In the viscothermal boundary layer, viscous effects dominate the flow and nonlinear effects are
negligible, whereas in the main flow region, the viscous effects are negligible and nonlinear
effects can be of importance,
we can reduce the full system of equations to the following closed set of equations for the areaaveraged velocity u(x,t), and area-averaged pressure p(x,t) [10]:
1 


∂u
∂p
∂u
+u
+
= 0,
ρ0
∂t
∂x
∂x


1
1 ∂p
∂p
∂u
∂2p
+u
+
+C 1 = 0,γ p ∂t
∂x
∂x
∂t 2


p
p0

γ

(1)
(2)

where ρ0 is the density of the gas in quiescent conditions, p0 the atmospheric mean pressure, γ theratio of specific heats, t time, x the coordinate in the direction of the duct length and C is defined
as:

 √
γ −1
ν0
C = 1+ √,(3)
Pr rh ρ0 c20
in which Pr is the Prandtl number, ν0 the kinematic viscosity at quiescent conditions, and c0 the
speed of sound at quiescent conditions. The hydraulic radius is formally defined as rh = S/Π, where5482

S is the duct cross sectional area, and Π the wetted perimeter. For a circular duct with radius R,12

rh ≡ R/2. ∂ 1p is the half-order derivative of the pressure to the time, which is defined as
∂t 2

1

1
√
1 ≡
π
2
∂t

∂2p

Zt
0

1 ∂p
√
dτ.
t − τ ∂t

(4)
1

1

In Eq. 4, it is assumed that p(x,t < 0) = 0. To compute the current value of ∂ 2 p/∂t 2 , the complete
time history of p(t) is required. As discussed in the introduction, this is an undesired property
for numerical computations. In the next section, we will treat an implementation which solves thisproblem.
2.1 IMPLEMENTATION OF THE FRACTIONAL ORDER DERIVATIVE
Following Lombard and Matigon [8], a transformation can be made to write the fractional derivative
as
1
Z∞
∂2p
=
φ (t, θ )dθ ,(5)
1∂t 2 θ =0
where φ satisfies∂φ2 ∂p= −θ 2 φ +
,
(6)
∂t
π ∂t
with φ (t = 0, θ ) = 0. This is called the diffusive representation of the fractional derivative, due tothe similarity of Eq. 6 with, for example, Newton’s cooling law. On first sight, the main advantage
is that the evolution of φ is local in time (it contains only integer-order derivatives). And in caseof a second order finite difference scheme, only the previous time instance is required to be kept
1
1
in memory. However, to compute ∂ 2 p/∂t 2 requires the computation of an improper integral in theθ -space.
In the diffusive approximation, the improper integral is estimated using a quadrature rule, on a fixed
set of L points:Z∞

L−1

φ (t, θ )dθ ≈ ∑ µ` φ` (t),

θ =0

(7)

`=0

where µ` are the quadrature coefficients. Now, φ` evolves according to:
∂ φ`
2 ∂p
= −θ`2 φ` +
,
∂t
π ∂t

` = 0...L − 1,(8)subject to the boundary condition at t = 0:φ` (t = 0) = 0,

` = 0...L − 1.

(9)Now, the challenge is to find a proper set quadrature coefficients and points (µ` , θ` ) for the estimation of the integral of Eq. 5. Here, a trade-off has to be made between computational efficiency andaccuracy, as the higher L, the more variables need to be solved at each grid and time point in the
numerical scheme.
2.2 Diffusive approximation quadrature coefficients
Following Lombard, we could determine the quadrature coefficients by minimizing the dispersion error between the exact frequency domain fractional order derivative and the diffusive approximation.

5483

However, it is found that by doing this, for low frequencies, this results in a larger dispersion error,
when comparing to the viscothermal wave number. Notably, it is found that for tubes with a diameter of 1 mm, this leads to a factor of 2 too low imaginary part of the wave number. The imaginary
part of the wave number is related to the damping of the acoustic waves. Therefore, we propose
to minimize the difference between the exact viscothermal wave number and the viscothermal wave
number based on the diffusive approximation.
By substituting the x-derivative of the small amplitude linearized version of Eq. 2 into the timederivative of the small amplitude linearized version of Eq.1 we can derive the following (lossy) 1D
wave equation:

√
1
γ −1
ν0 ∂ ∂ 2 p
1 ∂2p ∂2p
√
−
+
1
+
(10)
1 = 0,∂ x2
c20 ∂t 2
Pr rh c20 ∂t ∂t 2
where Pr denotes the Prandtl number. Upon transforming Eq. 10 to the frequency domain using
e+iωt -convention, we can derive the following Helmholtz-like ordinary differential equation:


 r 
γ − 1 1 ν0
d2 p̂2
+k 1+ 1+ √p̂ = 0,dx2
Pr rh iω(11)√
where p̂(ω) is the Fourier transform of the acoustic pressure, k = ω/c0 and i = −1. Here, the termin front of p̂ can be recognized as an approximation to the square of the viscothermal wave number:
Γ2approx . Applying the Fourier transform to Eqs. 7-8, it is found that the diffusive approximation has
the following frequency domain representation:
"
#
L−1
L−1
2iω µ`F ∑ µ` φ` (t) ≡ χ̃(ω) p̂(ω) = ∑p̂ (ω) ,(12)2`=0`=0 π θ` + iω
where F denotes the Fourier transform operator. Hence, in the diffusive approximation, we canderive that the viscothermal wave number will be estimated as
s

√
γ −1
ν0
χ̃(ω).
(13)
Γapprox (ω) = k 1 + 1 + √Pr iωrh
For circular ducts, an exact solution (in the low reduced frequency limit) for this viscothermal wave
number is available for the viscothermal wave number, which reads [12, 3]:sΓcirc = k
where

1 + (γ − 1) fκ
.
1 − fνih
J1 (i − 1) 2rδ jh
h
i,
fj =
(i − 1) δrhj J0 (i − 1) 2rδ jh

(14)

(15)

in which Jn is the cylindrical Bessel function of the first kind nth order and δν , δκ are the viscousand thermal penetration depth’s, respectively. They are defined as:
s
s
2κ0
2µ0
;
δν =
,
(16)
δκ =
ρ0 c p ωρ0 ω
where κ0 is the thermal conductivity at quiescent conditions, c p the specific heat at constant pressure, and µ0 the dynamic viscosity at quiescent conditions.

5484

L=2
L=3
L=4

Γapprox − Γex
Γex

10−1

10−2

10−3

102

103
Frequency [Hz]104

Figure 1. Normalized error in the viscothermal wave number of the diffusive approximation with
optimized coefficients for µ` and θ` according to the objective function of Eq. 18, for different
values of L. The tube has a diameter of 1 mm (rh = 0.25 mm).
To find a set of quadrature coefficients, a set K frequencies ωk , spanning the range from f0 = 100 Hz
to fK−1 = 20 kHz is chosen, which are logarithmically spaced over this range, according to k
20 · 103 K−1
fk = 100.
100
Next, the quadrature coefficients are obtained by minimization of the objective function J :(µ` , θ` ) = (µ` , θ` ) | min J (µ` , θ` ) ,µ` ,θ`(17)

where the objective function is defined as
K−1

Γ2approx (ωk )k=0

Γ2circ (ωk )J (µ` , θ` ) = ∑

−1 ,

(18)

To solve this nonlinear least squares problem, we have used the Trust Region Reflective algorithm,
as implemented in SciPy v1.2.1 [6]. The found optimal coefficients are dependent on L and K. It isfound that for increasing K, the found coefficients µ` and θ` converge to a certain value for given
L. For K ≈ 50, the optimal values for µ` and θ` are approximately independent of K.
Figure 1 shows the error in the wave number Γ as a function of frequency for different values of
L. As visible in this result, for L = 4, the error can be kept lower than 0.3% in the frequency rangeof 100 Hz to 20 kHz.
It should be noted that the found coefficients which minimize the dispersion error are dependent on
the the hydraulic radius of the duct. This is due to the rh in the denominator in Eq. 13.
2.3 Numerical implementation of the 1D model
The 1D model has been implemented in numerical code which uses the FEniCS project [1] general
finite element framework. In this framework, the weak formulation of the 1D model is implemented

5485

p = p0 + p1 sin(2πf t)  γ−1γ
T = T0 ppL0

T = T0u=0

Figure 2. Schematic of the resonator geometry for CFD computations, including annotations of theboundary conditions. For visibility, the y-scale is multiplied by a factor of 100.
using a quadratic shape functions and second order implicit time-stepping procedure and constanttime steps, for the variables p, u, and φ0 . . . φ3 . To solve the nonlinear problem at each time step,
the Newton-Raphson algorithm is used.

3

COMPARISON WITH CFD RESULTS: λ /4 RESONATOR

In this section, we compare results of the 1D model with the full nonlinear compressible NavierStokes model (CFD). To this end, a 2D duct has been generated with a length L corresponding toa quarter wavelength: L = λ /4, where λ = c0 / f .
3.1 CFD model
We solve the full nonlinear system of equations using the rhoPimpleFoam solver, from the OpenFOAM CFD suite, version 6.2 [5]. Figure 2 shows the geometry, including boundary conditions. On
the left side, the pressure perturbation is prescribed as a sinusoidal function of time. The according
temperature is computed using the adiabatic compression/expansion relation. The velocity is computed from the model using a pressureInletOutletVelocity boundary condition. As a result,the duct is driven in its λ /4 resonance, leading to a high amplitude (at least according to linear
isentropic acoustics) on the right side of the duct. The resulting amplitude on this side of the duct
is limited by the viscothermal and nonlinear effects.
The grid typically has 30 nodes in the x-direction, and approximately 12 nodes in the y-direction.
Attention has been given to put enough grid points in the viscothermal boundary layer. A grid study
has been performed to check for near grid-independence of the result. For the CFD model, it was
found that the time step has a strong influence on the amplitude computed on the right side. This
imposed a stronger restriction on the time step than the restriction based on the CFL condition. The
time step in these simulations is set to 0.1 µs.
3.2 1D model
As we are comparing for a 2D model, for the 1D model, the quadrature coefficients are determinedby using the parallel-plate formulation of the viscothermal Rott functions:fj =tanh ((1 + i) rh /δ j ),(1 + i) rh /δ j

(19)which lead to slightly different values of µ` and θ` . We checked that the error we make with
Γapprox is in the same order as for the circular duct case. The time step has been chosen equal to
the time step used in the CFD computations, and the number of grid points equal to the number in

5486

20000

10000

10000

CFD
1D Model0

Pressure at end [Pa]Pressure at end [Pa]

20000

0

−10000−10000−20000−20000
0

25

50

75

100
Time [ms]

125

150

175

200

CFD
1D model
190

192

194

196

198

200

Time [ms]Figure 3. Comparison of the pressure perturbation response at the right side of the resonator, for
p1 = 10, 000 Pa, and f = 100 Hz. Left: result over the full simulation time. Right: zoom in on thelast period. The horizontal blue line denotes the pressure amplitude according to the Low Reduced
Frequency theory.

CFD
1D Model

80000

60000

60000

40000

40000

Pressure at end [Pa]Pressure at end [Pa]

80000

20000

0

−2000020000

0−20000−40000CFD
1D model−40000
0.0

2.5

5.0

7.5

10.0
Time [ms]

12.5

15.0

17.5

20.0

19.0

19.2

19.4

19.6

19.8

20.0

Time [ms]Figure 4. Comparison of the pressure perturbation response at the right side of the resonator, for
p1 = 10, 000 Pa, and f = 1000 Hz. Left: result over the full simulation time. Right: zoom in on thelast period. The horizontal blue line denotes the pressure amplitude according to the Low Reduced
Frequency theory.
the x-direction of the CFD model.
3.3 Results
We compare results for the computed pressure perturbation on the right side of the resonator, due topressure perturbation amplitude p1 on the left side. The number of periods has been chosen to finda sufficiently periodic steady state. Due to the relatively high viscothermal damping, approximately
20 periods are sufficient to reach periodic steady state. The amplitude of p1 is set to 1·104 Pa, which
is such high that we expect nonlinear effects to impact the result. Figure 3 compares the result for
100 Hz. As visible, the amplitudes at the right side match quite well. As a reference, the horizontal
blue line denotes the amplitude expected from Low Reduced Frequency (LRF) theory, from which
we can derive:
p1| p̂R | =,(20)
cos (ΓL)
where p̂R is the pressure phasor on the right side of the resonator. Figure 4 shows the result forthe 1 kHz frequency. As visible on the right, the amplitude is significantly lower than for the linear

5487

LRF model. Moreover, the wave profile is quite distorted from a pure sinusoid. Besides of a small
phase shift, the 1D model matches quite closely with the 2D nonlinear simulation results.

4

CONCLUSIONS AND NEXT STEPS

We have presented a 1D numerical model to compute nonlinear viscothermal wave propagation in
narrow tubes / ducts. The model is implemented in a finite element framework and has been verified
using a full nonlinear compressible Navier Stokes simulation. The next steps are to compare results
for axisymmetric ducts and ducts with jumps in cross sectional area jumps, to finally arrive at amodel of the full BTE hearing aid system.

ACKNOWLEDGEMENTS
The authors would like to acknowledge Sonova A.G. for the financial support of this research
project.

REFERENCES[1] M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring,
M. E. Rognes, and G. N. Wells. The FEniCS Project Version 1.5. Archive of Numerical Software, 3(100), 2015.
[2] W. Chester. Resonant oscillations in closed tubes. J. Fluid Mech, 18(1):44–64, 1964.
[3] J. De Jong. Numerical modeling of thermoacoustic systems. PhD thesis, Universiteit Twente,
Enschede, 2015.
[4] H. Haddar, J.-R. Li, and D. Matignon. Efficient solution of a wave equation with fractionalorder dissipative terms. Journal of Computational and Applied Mathematics, 234(6):2003–2010,
July 2010.
[5] H. Jasak, A. Jemcov, Z. Tukovic, and others. OpenFOAM: A C++ library for complex physics
simulations. In International workshop on coupled methods in numerical dynamics, volume 1000,
pages 1–20. IUC Dubrovnik Croatia, 2007.
[6] E. Jones, T. Oliphant, P. Peterson, and others. SciPy: Open source scientific tools for Python,
2001.
[7] L. Landau and E. Lifshitz. Fluid Mechanics, Vol. 6. Course of Theoretical Physics, 1987.
[8] B. Lombard and D. Matignon. Diffusive approximation of a time-fractional Burger’s equationin nonlinear acoustics. arXiv:1602.07205 [physics], Feb. 2016. arXiv: 1602.07205.
[9] L. Menguy and J. Gilbert. Weakly Nonlinear Gas Oscillations in Air-Filled Tubes; Solutions
and Experiments. Acta Acustica united with Acustica, 86(5):798–810, 2000.
[10] N. Sugimoto and D. Shimizu. Boundary-layer theory for Taconis oscillations in a helium-filled
tube. Physics of Fluids (1994-present), 20(10), Oct. 2008.
[11] G. W. Swift. Thermoacoustics: A unifying perspective for some engines and refrigerators. Acoustical Society of America, Melville, NY, USA, 2003.
[12] H. Tijdeman. On the propagation of sound waves in cylindrical tubes. Journal of Sound andVibration, 39(1):1–33, Mar. 1975.

5488