T. Oshima, M.

Imano, A full finite-volume time-domain approach towards general-purpose code

development for sound propagation prediction with unstructured mesh, Proceedings of Inter-Noise
2008 (Shanghai), No. 287, 15 pages in CD-ROM, 2008.
The code is available at http://openfoamwiki.net/index.php/Contrib_potentialWaveFoam

A Full Finite-Volume Time-Domain Approach towards

General-purpose Code Development for Sound Propagation
Prediction with Unstructured Mesh

Takuya Oshima1 , Masashi Imano2

1
Faculty of Engineering, Niigata University,
8050 Igarashi-Ninocho, Niigata City, 9502181, Japan
2
School of Engineering, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo, 1138656, Japan

ABSTRACT
While finite-difference time-domain (FDTD) approach has widely been accepted as a
simple, fast and proven measure for numerical sound propagation prediction, it has
suffered enforcement of orthogonal mesh usage and lack of general-purpose solver code.
Due to the weaknesses we even now have to write solvers and pre/postprocessing codes
on case-by-case basis. This has made real industrial applications of the technique with
complex geometries difficult. The issue was addressed here through introduction of
a full finite-volume time-domain (FVTD) approach meant as a replacement for the
FDTD technique. The main strength of the FVTD approach in principle is a great
flexibility in mesh handlings which allows full unstructured meshes containing arbi-
trary shapes of polyhedra. Thus the strength opens possibility of using a vast variety
of general-purpose pre/postprocessors designed for finite volume or finite element
meshes. The proposed FVTD technique, along with an acoustic impedance boundary
condition specifically developed for use with the technique, was formulated, imple-
mented and tested using solutions obtained by the FDTD technique as benchmarks.
Both techniques were confirmed to produce identical results under identical geome-
try, mesh and computational conditions. The demanded processor times and memory
usages for FVTD calculations were more than ten times of FDTD calculations, which
still was thought to be allowable up to medium-sized problems with recent advance-
ments in processor performance taken into account. The results obtained under full
unstructured tetrahedral meshes, however, showed numerical dispersions and diffu-
siveness, which indicated necessity of further works.

1 INTRODUCTION
While finite-difference time-domain (FDTD) approach has widely been accepted as
a simple, fast and proven measure for numerical sound propagation prediction, it has
suffered enforcement of orthogonal mesh usage and lack of general-purpose solver code.

1
Email address: oshima@eng.niigata-u.ac.jp
2
Email address: imano@arch.t.u-tokyo.ac.jp
Thus we even now often have to write solvers and pre- and postprocessing codes on case-
by-case basis. The resultant negative effects arising from the shortcomings are mainly
twofold:

1. With recent advancements in processor performance taken into account, the short-
comings has made relative cost of human powers with regard to pre- and post-
processings much higher than the cost devoted for numerical simulation itself, in
particular for small- and medium-sized cases. The situation has made real indus-
trial applications of the FDTD technique for complex geometries difficult.

2. Despite the wide use of the FDTD technique in many literatures, we still do not,
and will not be able to, have a common case sharing framework for accumulat-
ing case examples which would stimulate communications between computational
acousticians unless the circumstance changes.

On the other hand, if we were to propose a new approach as an alternative to the

FDTD technique, we also have to keep in mind that the FDTD technique is so widely
used since it does have its own strengths:

1. The FDTD technique is mathematically simple and thus easy to implement.

2. Thanks to the simplicity, its computational cost per cell is extremely low.

The possible new approach is expected to maintain these features at least to some extent.
All the issue was addressed here through introduction of a full finite-volume time-
domain (FVTD) approach meant as a replacement for the FDTD technique, as opposed to
a mixed finite-difference and finite-volume approach meant as a complement for handling
body-fitted cells in FDTD computation [1]. The main strength of the full FVTD approach
in principle is a great flexibility in mesh handlings which allows using unstructured meshes
containing arbitrary shapes of polyhedra, while maintaining relative simplicity compared
to other more advanced approaches such as BEM or FEM. Hence the strength opens
possibility of using a vast variety of general-purpose pre- and postprocessors designed for
finite volume or finite element meshes, keeping computational costs per cell relatively low
at the same time.
To fully exploit the inherent feature of finite volume approach, the proposed technique
has been implemented as a user application on top of an open-source finite volume based
toolkit, OpenFOAM [2]. With this approach, not only the developers can make maximum
use of its tried and proved finite volume operators and I/O libraries, but also a user can
get full access to the included mesh format converters and postprocessing exporters.
The implementation was tested for contrasting its accuracy and computational costs
against the FDTD approach. Furthermore, several formulations of the acoustic impedance
boundary condition were carried out comparative tests under orthogonal and unstructured
meshes to choose the accompanying implementation of the boundary condition with the
main FVTD implementation.

2 FORMULATION AND IMPLEMENTATION

2.1 Finite-Volume Formulation
The wave propagation equation represented in velocity potential φ is written as the
following equation.
∂2φ
= c20 ∇2 φ (1)
∂t2
 Inside of
computational domain
Boundary
cell face b
φP φb
Internal
cell face f Cell of
φN interest P

Figure 1: Unstructured mesh system.

where t, c0 , ∇2 are time, propagation speed of the wave and Laplacian operator respec-
tively. Using φ, pressure p and particle velocity u are written as follows.
∂φ
p = ρ , (2)
∂t
u = −∇φ. (3)

Eq. (1) is discretized under unstructured grid system as shown in Fig. 1 where the
definition point of physical quantities are taken at the barycenter of each control volume
(CV). For the left hand side of Eq. (1), by integrating over the CV with time-invariant
volume V and applying central time-differential scheme, we get
∂2 ∫ φn+1 − 2φn + φn−1
φ dV ≈ V
∂t2 V ∆t2
where φn−1 , φn , φn+1 denote the values of φ at the (n − 1)-th, n-th, (n + 1)-th steps of
time step ∆t. For the right hand side, by integrating Eq. (1) within a CV and applying
divergence theorem, we get
∫ ∫
c20 ∇2 φdV = c20 dS · ∇φ
V S
∑
≈ c20 S f · (∇φ)f (4)
f

where S f denotes the face area vector of the f-th face that constitutes polyhedral CV in
question as follows.

S f = Sf nf (5)

where Sf , nf are the area and the unit outward normal vector of the face f respectively.
If a vector connecting the centers of the CV P and its adjacent CV N, dP N , is parallel
to S f , S f (∇φ)f is written in terms of ∂φ/∂nf , the surface-normal gradient of φ. Thus the
term within the summation in the rightmost hand side of Eq. (4) is discretized as follows.
∂φ
S f · (∇φ)f = Sf
∂nf
φN − φP
≈ Sf (6)
|dPN |
However, if dP N is nonorthogonal to S f , S f has to be decomposed into its orthogonal part
∆f and nonorthogonal part kf .

S f · (∇φ)f = ∆f · (∇φ)f + kf · (∇φ)f

 a) Underrelaxed correction b) Orthogonal correction c) Over-relaxed correction
Figure 2: Nonorthogonal mesh treatment vectors.

The first term of the right hand side of the equation above, the orthogonal part, is
discretized similarly to Eq. (6) as follows.
φN − φP
∆f · (∇φ)f ≈ |∆f |
|dPN |
The nonorthogonal part, (∇φ)f in the second term, is given by interpolating the gradient
of φ at the centers of CVs P and N.

(∇φ)f = fx (∇φ)P + (1 − fx )(∇φ)N (7)

Here, the interpolation coefficient fx and the gradient (∇φ)P are given as follows.

fN
fx = ,
|dPN |
1 ∫
(∇φ)P = dS φ
V S
1 ∑
≈ Sφf
V f

where φf is a face-interpolated value of φ at the center of CVs.

2.2 Nonorthogonal Corrections for Laplacian Term

The orthogonal and nonorthogonal component vectors ∆f and kf can be calculated
arbitrarily. In the present paper we are comparing three techniques proposed by Jasak
[3].

Underrelaxed correction:
dPN · S f
∆f = dPN (8)
|dPN |2

Orthogonal correction:
dPN
∆f = Sf (9)
|dPN |

Overrelaxed correction:
dPN
∆f = S2 (10)
dPN · S f f
Jasak concludes overrelaxed correction to be the best correction technique as a result of
testing the techniques with solutions of the Laplacian equation under 10◦ –60◦ skewed two
dimensional quadrilateral meshes. However, typical skewnesses for tetrahedral meshes,
which most of meshing softwares support for three-dimensional geometries, typically fits
in a relatively small range of under 20◦ . Furthermore, behaviors for the wave equation is
not known. Hence the effectiveness of the techniques will be reinvestigated in the following
section.

2.3 Rigid Boundary Conditions

On acoustically rigid boundaries b, normal component of particle velocity ub is fixed
to zero.

ub = n b · u b = 0 (11)

Substituting the relationship above to Eqs. (3) and (5) leads to the equation below which
represents the surface normal gradient of φ being zero.

S b · (∇φ)b = 0

2.4 Normal Incidence Acoustic Impedance Boundary Conditions

Acoustic impedance at boundaries under normal incidence condition z is given as
p
z = ,
ub
where ub is the normal component of particle velocity. Substituting Eqs. (2) and (3) to
the equation above, we get an advection equation of velocity potential φ.
∂φ 1 ∂φ
= − , (12)
∂nb cb ∂t
where
z
cb = .
ρ
The advection equation Eq. (12) has to be discretized to calculate the surface normal
gradient of φ at the time step n + 1 at the boundary face barycenter b shown in Fig. 1,
namely ∂φ/∂nb |n+1 . In the deriving process of the discretized forms of Eq. (12), three
discretization schemes are applied for each of time and spacial directions, which leads to
nine combinations of schemes in total. The discretized equations are shown below, with
φ subscripted by P and b being the values of φ at barycenters of boundary-internal CV
and boundary face respectively and ∆n being the distance between the barycenters. The
abbreviated type name that appears in the subtitle of each formulation is referred to later
in Section 4.

2.4.1 Upwind Types

Upwind types are derived regarding time derivatives at P as approximations of those
at b.

Central differencing formulation (Type U-C)The formulation is obtained by applying a

second order central differencing scheme to the time derivative of Eq. (12).
¯n+1
∂φ ¯¯ 1 φn+1 − φnP
¯ = − P
(13)
∂nb ¯ cb ∆t
Second-order backward formulation (Type U-B) The formulation is obtained by applying
second order backward differencing scheme [4] to the time derivative.
¯n+1
∂φ ¯¯ 1 3φn+1 − 4φnP + φn−1
¯ = − P P
(14)
∂nb ¯ cb 2∆t

Crank-Nicholson type formulation (Type U-CN) Crank-Nicholson interpolation scheme is

applied to obtain the spatial derivative of Eq. (12) at time step n + 1/2, while the time
derivative is discretized with a second order differencing scheme.
 ¯n+1 ¯n 
¯
1  ∂φ ¯ ∂φ ¯¯  1 φn+1 − φnP
¯ + ¯ = − P
(15)
2 ∂nb ¯ ∂nb ¯ cb ∆t

By reducing the equation above, we have a recurring formula with regard to the spatial
derivative.
¯n+1 ¯n
∂φ ¯¯ 2 φn+1 − φnP ∂φ ¯¯
¯ = − P
− ¯ (16)
∂nb ¯ cb ∆t ∂nb ¯

2.4.2 Predictor-Corrector Types

As noted above, the spatial derivative values ∂φ/∂nb |n+1 obtained by Eqs. (13), (14)
and (16) are in fact those at the CV barycenter P, ∂φ/∂nb |n+1
P . To extrapolate the value
to the boundary b, we first calculate a predictor with
¯n+1
∂φ ¯¯
φn+1
∗ = φn+1 + ¯ ∆n, (17)
P
∂nb ¯P

and then apply a predictor-corrector scheme [4] to the surface-normal direction to obtain
the surface gradient ∂φ/∂nb |n+1 .

Central differencing formulation (Type PC-C) First we calculate the predictor similarly to
the upwind formulation, Eq. (13),
¯n+1
∂φ ¯¯ 1 φn+1 − φnP
¯ = − P
∂nb ¯P cb ∆t

and then we obtain a corrected value, by way of Eq. (17), as follows.

¯n+1 ( ¯n+1 
∂φ ¯¯ 1 ∂φ ¯¯ φn+1 − φnb 
¯ = ¯ − ∗
∂nb ¯ 2 ∂nb ¯P cb ∆t

φnb in the equation above is calculated by the following equation (the same applies here-
after).
¯n
∂φ ¯¯
φnb = φnP + ¯ ∆n (18)
∂nb ¯
Second-order backward formulation (Type PC-B) Similarly to Eq. (14), we obtain a pre-
dictor
¯n+1
∂φ ¯¯ 3φn+1 − 4φnP + φPn−1
¯ = − P
,
∂nb ¯P 2cb ∆t

followed by a corrected value by way of Eq. (17),

¯n+1  ¯n+1 
∂φ ¯¯ 1  ∂φ ¯¯ 3φn+1 − 4φnb + φn−1
¯ = ¯ − ∗ b .
∂nb ¯ 2 ∂nb ¯P 2cb ∆t

Crank-Nicholson type formulation (Type PC-CN) By replacing ∂φ/∂nb |n in the correspond-

ing upwind formulation to ∂φ/∂nb |nP , we obtain
¯n+1 ¯n
∂φ ¯¯ 2 φn+1 − φnP ∂φ ¯¯
¯ = − P
− ¯
∂nb ¯P cb ∆t ∂nb ¯P

followed by a corrected value by way of Eq. (17),

¯n+1 ( ¯n+1 ¯n 
∂φ ¯¯ 1 ∂φ ¯¯ 2 φn+1
∗ − φnb ∂φ ¯¯ 
¯ = ¯ − − ¯ .
∂nb ¯ 2 ∂nb ¯P cb ∆t ∂nb ¯

2.4.3 Algebraic Types

The left hand side in Eq. (12) is discretized as
¯n+1
∂φ ¯¯ φn+1
b − φn+1
P
¯ =
∂nb ¯ ∆n

while( the right hand

) side is discretized at the boundary surface b. Then we can ob-
tain φn+1
b − φn+1
P /∆n, namely the second order differentiated form of ∂φ/∂nb |n+1 , by
algebraic reduction.

Central differencing formulation (Type A-C) By discretizing both sides of Eq. (12), we
obtain
φn+1 − φn+1 1 φn+1 − φnb
b P
= − b
.
∆n cb ∆t
( )
After approximating φn+1
b − φn+1
P /∆n to ∂φ/∂nb |n+1 , the equation above is reduced to
¯n+1
∂φ ¯¯ φnb + φn+1
P
¯ =
∂nb ¯ cb ∆t + ∆n

Second-order backward formulation (Type A-B) Similarly to Type A-C, we obtain

¯n+1
∂φ ¯¯ 4φnb − φn−1
b − 3φnP
¯ = .
∂nb ¯ 2cb ∆t + 3∆n
Crank-Nicholson type formulation (Type A-CN) Similarly to Type A-C, we obtain
( ¯n )
1 φn+1 − φn+1 ∂φ ¯¯ 1 φn+1 − φnb
b P
+ ¯ =− b
2 ∆n ∂nb ¯ cb ∆t

and its final reduced form

¯n+1 ¯n
∂φ ¯¯ 2 ( ) cb ∆t ∂φ ¯¯
¯ = φnb − φn+1 − ¯ . (19)
∂nb ¯ cb ∆t + 2∆n P
cb ∆t + 2∆n ∂nb ¯

3 TESTS FOR VALIDATION AND NONORTHOGONAL CORRECTION TECHNIQUES

3.1 Computational Setups
To validate the proposed FVTD technique and to test the correction techniques under
unstructured meshes, a comparative test using a sound propagation problem in a closed
cube of 1 m × 1 m × 1 m, one of the AIJ-BPCA (Benchmark Platform on Computational
Methods for Architectural and Environmental Acoustics) [5] problems, was carried out.
The detail of the tested cases are shown in Tab. 1.

Case 1 The problem was solved using a conventional FDTD code written in Fortran
77 employing a pressure-particle velocity leapfrog scheme. The case is meant to be the
benchmark case to which the results obtained by the proposed technique is compared for
validation. Each edge of the cube was divided to 81 subedges to create a mesh of cell
width ∆x = 0.0123 m and the number of cells 531 441 (Fig. 3(b)). The time step ∆t and
the Courant number Co were set to 0.02 ms and 0.96 respectively.

Case 2 The problem was solved with the proposed technique under a hexahedral orthog-
onal mesh and setup both identical to the ones for Case 1.

Case 3 The problem was solved with the proposed technique under a nonuniform tetra-
hedral unstructured mesh automatically generated by an open-source mesher, Gmsh [6].
The characteristic length lc (the length with which each edge of the cube is divided) is set
to 0.025 m, to make a mesh with the number of CVs 531 333 (roughly the same as Cases
1 and 2). The ratio of maximum and minimum CV edge lengths of the generated mesh
was 6.32. The time step ∆t was set to 0.0049 ms to keep the maximum Courant number
to 0.99. In this case no nonorthogonal techniques were applied.

Cases 4, 5 and 6 The setups are same as Case 3, except that the underrelaxed, orthogonal
and overrelaxed nonorthogonal correction techniques were applied to Cases 4, 5 and 6
respectively.

Common Conditions For all cases the initial values of φ were set to represent the pressure
and particle velocity conditions of

 cos 8πr + 1
(r < 0.125)
p−1/2 (r) = 2 [Pa] (20)

0 (otherwise)
u0 (r) = 0 (21)

where r [m] is the distance from the center of the cube. All cases were run up to t = 0.04s.
 z
0.5

R4
O 0.2 0.3 R3
0.4 0.5 y
0.5
x
R2
Dimension in [m]
(a) Problem geometry (b) Surface (Hexahedral)

(c) Surface (Tetrahedral) (d) Internal (Tetrahedral)

Figure 3: (a) Problem geometry of the benchmark problem AIJ-BPCA B0-1T, (b) Surface
mesh for Cases 1 and 2 (coarsened by factor of 2 for visibility), (c) Surface mesh for Cases 3–6
(coarsened by factor of 2) and (d) Surface and internal mesh (coarsened by factor of 8) for Cases
3–6.

Table 1: Computational setups.

Case 1 2 3 4 5 6
Approach FDTD FVTD
Type of mesh — Hexahedral Unstructured tetrahedral
Number of cells/CVs 813 = 531 441 531 333
∆x [m] 0.0123 —
lc [m] — 0.025 (40 elements per edge)
∆t [ms] 0.02 0.0049
c0 [m/s] 343.7
Courant number 0.96 0.99 (max)
Nonorth. correction — Uncorrected Underrelaxed Orthogonal Overrelaxed
Initial condition A single wave of offset cosine (Eqs. (20), (21))

3.2 Results and Discussions

The transient sound pressure waveforms at the receiving point R2 shown in Fig. 3(a)
are plotted in Fig. 4, using the result of Case 1 as the benchmark case for comparison
with other cases. From Fig. 4(a), one can see that the results of FDTD and FVTD
 0.4
Pressure [Pa] a) Cases 1, 2 Case 1 (FDTD)
0.2 Case 2 (FVTD-Hexahedral)
0
-0.2
-0.4
0.4
b) Cases 1, 3 Case 1 (FDTD)
Pressure [Pa]

0.2 Case 3 (FVTD-Tetrahedral-Uncorrected)

0
-0.2
-0.4

c) Cases 1, 4 Case 1 (FDTD)

Case 4 (FVTD-Tetrahedral-Underrelaxed)

d) Cases 1, 5 Case 1 (FDTD)

Case 5 (FVTD-Tetrahedral-Orthogonal)

0.4
e) Cases 1, 6 Case 1 (FDTD)
Pressure [Pa]

0.2 Case 6 (FVTD-Tetrahedral-Overrelaxed)

0
-0.2
-0.4
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
Time [s]

Figure 4: Transient sound pressure waveforms at the receiving point R2.

techniques agrees so precisely that they can be regarded as virtually identical results.
From the results one can verify the proposed FVTD technique has the same accuracy as
a conventional FDTD under identical geometry, mesh and computational setups.
On the other hand, from the comparison of Cases 1 and 3 in Fig. 4(b), the waveform
obtained by the FVTD technique under the tetrahedral mesh is phasing forward in about
1.5 % and the overall waveform is gradually dispersing over time. In addition, the results of
Cases 4–6, which were meant to confirm the effects of the correction techniques to correct
the unwanted behavior observed in Case 3, are shown in Fig. 4(c)–(e) respectively. Despite
the employment of the correction techniques, the drifts of the phase did not improve, or
even worse, waveforms started to oscillate and diverge eventually. Although the starting
times of the oscillations become slightly later for relaxed cases, the overall characteristics
do not differ much.
From the results one can conclude that while one can expect identical results between
FDTD and FVTD techniques under identical setups, there remains works for the FVTD
technique in reducing the phase error coming from nonorthogonalities of unstructured
grids.
 Table 2: Processor and memory usages.
Case 1 2 3 4
Processor [s] 28.0 343 865 2 013
Per time step [s] 0.0140 0.172 0.106 0.247
Memory [MB] 18 301 260 260

Figure 5: Geometry of one-dimensional tube.

3.3 Computational Loads

To compare the proposed technique with the conventional FDTD from the standpoint
of computational loads, processor times and memory usages were instrumented for Cases
1–4, as shown in Tab. 2. The instrumentations were carried out on an Opteron 2.4 GHz 64-
bit Linux platform. The FVTD computations turned out to require more than ten times
for processor and memory usages. Hence it should be noted that, from the computational
load of view, the proposed FVTD technique is not meant to completely replace FDTD
especially in large cases, but should rather be used for small to medium cases where rapid
preprocessing (case setup) and postprocessing have of particular importance.
Although Cases 2 and 3 have roughly the same number of cells, one may notice that
Case 3 requires smaller amount of computational time per time step. This is because the
computational load required in calculating Laplacian is determined mostly by the number
of faces per CV, as shown in Eq. (4). It is also shown that, from Cases 3 and 4, applying
nonorthogonal correction technique more than doubles the processor usage.

4 COMPARATIVE TESTS FOR THE ACOUSTIC IMPEDANCE BOUNDARY CONDITION

The nine different formulations of the normal acoustic impedance boundary condition
proposed in Section 2.4 were tested under an acoustic tube as a one dimensional problem,
square-shaped domain as a two-dimensional problem, and an eighth of a sphere as a
three-dimensional problem.

4.1 One-dimensional Tests

4.1.1 Computational Setups
To test the formulations under normal incidence conditions, an acoustic tube shown
in Fig. 5 was uniformly meshed with ∆x being 0.1 m. The right end of the tube was
employed an impedance boundary condition with the characteristic impedance of the air,
ρc0 . All other boundaries were set to have been rigid. The time step ∆t was set to 0.291
ms which corresponds to the Courant number of 1.0. The half wave of a cosine pulse with
the wavelength 5∆x was given at the left end of the tube as an initial condition.

4.1.2 Results
Total acoustic energy levels in the tube over time were shown in Fig. 6. The for-
mulation that corresponds to each Type is defined in Section 2.4. Types U-C, PC-CN
and A-CN similarly show good attenuations after t = 0.01 s, when the absorption of
the wavefronts that reached the right end of the tube starts. On the other hand, Types
 Type Type
Type Type
Type Type

Type Type
Type Type
Type Type

Type Type
Type Type
Type Type

Figure 6: Attenuation of total acoustic energy Figure 7: Attenuation of total acoustic energy
in the tube over time. in the square domain over time (∆t = 0.05
ms).

PC-C, PC-B, A-C and A-B only show smaller attenuation than the former three types.
Furthermore Types U-B and U-CN quickly diverged as soon as the wavefronts reached
the right end.

4.2 Two-dimensional Tests

4.2.1 Setup
Following the reference [7], an acoustic wave propagation problem in a flat square
domain was solved. The geometry shown in Fig. 10 was orthogonally meshed with
the cell width 0.05 m. Acoustic impedance boundary condition with the characteristic
impedance of the air, ρc0 , was given for all edges. The time step ∆t was set to two values
of 0.05 ms and 0.1 ms, which correspond to the Courant numbers of 0.486 and 0.972. A
cosine pulse with its radius being 10 times the cell width was given at the center of the
square.

4.2.2 Results
Total acoustic energy in the field over time are plotted in Figs. 7 and 8 for ∆t = 0.05
ms and 0.1 ms respectively. When ∆t = 0.05 ms, the energies diverged for Types U-
CN, PC-B and PC-CN. It can be seen that types that employed higher order schemes in
 Type Type
Type Type
Type Type

Type Type
Type Type
Type Type

Type Type
Type Type
Type Type

Figure 8: Attenuation of total acoustic energy Figure 9: Attenuation of total acoustic energy
in the square domain over time (∆t = 0.1 ms). in the eighth of a sphere over time.

Figure 10: Geometry of computational domain.

a given formulation type have a tendency of divergence, except for the algebraic types
which did not diverge (Types A-C, A-B, A-CN). On the contrary, when ∆t = 0.1 ms, all
types of the upwind type (Types U-C, U-B, U-CN) and Type PC-CN diverged. Again, no
divergence was observed for the algebraic type. Especially Type A-CN showed the best
attenuation of −35 dB at t = 0.04 s, and the variation between the results of ∆t = 0.05
ms and 0.1 ms was smallest among the all types.
 Figure 11: Geometry of eighth-spherical computational domain.

4.3 Three-dimensional Tests

4.3.1 Setup
To test the boundary conditions under a three-dimensional problem, an eighth of
a sphere shown in Fig. 11 was meshed with tetrahedral cells of characteristic length
lc = 0.0125 m. The spherical surface was given an acoustic impedance boundary condition
with characteristic impedance of the air, ρc0 . Rigid boundary conditions was imposed for
other boundary surfaces. Time step ∆t was set to 2.36 × 10−5 s which corresponds to
the maximum Courant number of 0.99. The calculations were run up to 0.01 s. A cosine
pulse of radius 0.125 m was given at the center of the sphere (point S shown in Fig. 11)
as the initial condition.

4.3.2 Results
The results are shown in Fig. 9 as total energy attenuations in the sound fields
over time. All of the predictor-corrector types (Types PC-C, PC-B, PC-CN) and Type
U-CN resulted in divergence. In the remaining converged types, the algebraic types
showed better attenuations in 2–3 dB than the upwind types at t = 0.002 s where the
spherical wave reaches the spherical boundary surface. Within the algebraic types Type
A-C indicates slightly better attenuation in about 0.5 dB.

4.4 Discussion
Looking through all the tests under one-, two- and three- dimensional geometries,
in overall the algebraic types have good characteristics in that in most tests they show
good attenuation characteristics, and in that even in worst tests they do not diverge.
Further, of the three algebraic types, Type A-CN was among the best types in one- and
three-dimensional tests, and the unarguable best under two dimensional tests. Thus one
can conclude Type A-CN as the best performing formulation for the normal acoustic
impedance boundary condition.

5 CONCLUSIONS
To overcome several inherent shortcomings in the currently widely used FDTD-based
sound propagation technique, such as enforced usage of orthogonal meshes and lack of
general-purpose solver code, the authors presented a fully finite-volume time-domain
(FVTD) approach. The proposed FVTD technique, along with an acoustic impedance
boundary condition specifically developed for use with the technique, was formulated and
implemented on top of an open-source finite volume based toolkit, OpenFOAM. The im-
plementation was tested using solutions obtained by the FDTD technique as benchmarks.
Both techniques were confirmed to produce identical results under identical geometry,
mesh and computational conditions. The demanded processor times and memory us-
ages for FVTD calculations were more than ten times of FDTD calculations, which still
was thought to be allowable up to medium-sized problems with recent advancements
in processor performance taken into account. The boundary condition proved to show
good attenuations on one-, two- and three-dimensional meshes including an unstructured
mesh. The overall results obtained under full unstructured tetrahedral meshes, however,
showed numerical dispersions and diffusiveness, despite of nonorthogoanl corrections. The
nonorthogoanl correction problems indicated necessity of further works.

ACKNOWLEDGEMENTS
Parts of the work were supported by JSPS Grant-in-Aid for Scientific Research (A)
19206062, (B) 19360264 and by MEXT Grant-in-Aid for Young Scientists (B) 19760402.

REFERENCES
[1] Botteldooren, D. Acoustical finite-difference time-domain simulation in a quasi-
cartesian grid. J. Acoust. Soc. Am., Vol. 95, No. 5, pp. 2313–2319, 1994.5.

[2] Weller, H. G., Tabor, G., Jasak, H., and Fureby, C. A tensorial approach to computa-
tional continuum mechanics using object-oriented techniques. Computers in Physics,
Vol. 12, No. 6, pp. 620–631, December 1998.

[3] Jasak, H. Error analysis and estimation for the finite volume method with applications
to fluid flows. PhD thesis, Imperial College, June 1996.

[4] Ferziger, J. H. and Perić, M. Computational Methods for Fluid Dynamics. Springer-
Verlag (Berlin), 1996. (Japanese version).

[5] Architectural Institute of Japan. Benchmark platform on computational methods for

architectural / environmental acoustics. http://gacoust.hwe.oita-u.ac.jp/AIJ-BPCA/.

[6] Geuzaine, C. and Remacle, J.-F. Gmsh: a three-dimensional finite element mesh gen-
erator with built-in pre- and post-processing facilities. http://www.geuz.org/gmsh/.

[7] Naito, Y., Yokota, T., Sakamoto, S., and Tachibana, H. The study of complete
absorption boundary for open region calculation by FDM. Proc. The 2000 autumn
meeting of the Acoustical Society of Japan, Vol. II, pp. 751–752, September 2000.