12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Numerical simulation of the sibilant /s/ sound articulation process
HsuehJui Lu
Sakuya Sugimoto
Department of Computational Science
Department of Computational Science
Kobe University
Kobe University
1-1 Rokkodaicho, Nada Ward, Kobe, Hyogo
1-1 Rokkodaicho, Nada Ward, Kobe, Hyogo
sugimotodot@stu.kobe-u.ac.jp
Lsray@stu.kobe-u.ac.jp
Tsukasa Yoshinaga
ChungGang Li
Department of Mechanical Engineering
Department of Computational Science
Toyohashi University of Technology
Kobe University
1-1 Hibarigaoka, Tempaku-cho, Toyohashi, Aichi
1-1 Rokkodaicho, Nada Ward, Kobe, Hyogo
yoshinaga@me.tut.ac.jp
cgli@aquamarine.kobe-u.ac.jp
Kazunori Nozaki
Akiyoshi Iida
Dental Hospital
Department of Mechanical Engineering
Osaka University
Toyohashi University of Technology
1-8 Yamadaoka, Suita, Osaka
1-1 Hibarigaoka, Tempaku-cho, Toyohashi, Aichi
knozaki@dent.osaka-u.ac.jp
iida@me.tut.ac.jp
Makoto Tsubokura
Department of Computational Science
Kobe University
1-1 Rokkodaicho, Nada Ward, Kobe, Hyogo
tsubo@tiger.kobe-u.ac.jp

ABSTRACTThe numerical simulation of the articulation process of thesibilant /s/ sound was conducted with the moving realistic
vocal tract geometry of the Japanese word /usui/ using animplicit compressible flow solver. The moving immersedboundary method with the hierarchical structure grid was
adopted to approach the complex geometry of the humanspeech organs. The sound emitted from the vocal tract was
analyzed by the spectrogram to clarify the effect ofvelopharyngeal closure. The sibilant /s/ sound was gradually
generated by reducing the constriction size. However, while thevelopharyngeal valve was not completely closed, the soundwas affected by the nasal emission. After the closure of thevalve, a typical /s/ sound was produced by the turbulent airflow
in the vocal tract. These results provide the underlying insights
necessary for clinical treatment of speech difficulty for sibilant
sound.

escapes through the nasal cavity during the sibilant soundproduction, the nasal emission occurs. The nasal emission is
inaudible with the large velopharyngeal opening. However,when the velopharyngeal opening becomes small, the undesired
sound might be generated.
The effect of the velopharyngeal opening and themechanism of the nasal emission have been investigated by themeasurements and numerical simulations. Bunton and Story(2012) examined the relationship among perceptual ratings ofthe hypernasality, and noted the differences in nasalance, sizeof the nasal valve opening and perceptual ratings ofhypernasality among the three English vowels, /i/, /u/, and /ɑ/.
The results indicated that the listeners could detect thehypernasality for the high and low vowels with nasal port areasfrom 0.01 to 0.15 cm2, respectively. Moreover, the further work
(Bunton, 2015) indicated that the perceptual ratings of thehypernasality and the nasalance increased until nasal port areas
reached 0.16 cm2.
A series of studies (Sundström et al., 2020) was conductedby utilizing the numerical flow simulations to investigate theeffect of the velopharyngeal insufficiency and indicated that theratio of channel areas at the velopharyngeal valve and the oral
constriction can be considered as a factor to determine theairflow configuration during the production of sibilant /s/sound. In addition, the results showed that the acoustic energy
in the nasal cavity and the far-field sound were directly related
to the size of the velopharyngeal opening. However, the almost
all studies of the flow simulation used a fixed vocal tract
geometry which only provides the stationary result instead of
the time-varying sound production in the actual human speechprocess.
In the present study, therefore, the numerical simulation ofsibilant /s/ production with the realistically moving vocal tract
was conducted to investigate the articulation process ofvelopharyngeal closure and tongue movement. We simulate thesibilant /s/ sound production process in the part of the Japanese

INTRODUCTION
The sibilant sound is one of the most difficult articulationsin the consonant production, and it is known to acquire in thelatest leaning process. Therefore, it is often reported inarticulation disorders due to the insufficient tongue formation.
The sibilant /s/ is generated by the airflow passing through aconstriction formed between the alveolar ridge and the tongue
tip. When the airflow is forced to go through the narrow
constriction (sibilant groove) of the tongue tip, the turbulent jet
flow is generated, and the velocity fluctuation downstream
from the constriction is considered as the source of the sound.
The sibilant fricatives are characterized as broadband noise in
the frequency range above the characteristic peaks,
approximately 4–7 kHz for /s/.
The velopharyngeal valve, which consists of the velum(soft palate), the lateral and posterior pharyngeal walls, controls
the amount of airflow going to the nasal cavity. If the flow

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
word /usui/ which is extracted after the vowel production of /u/to the sibilant sound /s/, including the tongue elevation and thevelopharyngeal valve closure. The acoustic field was directlypredicted with the flow field by the compressible Navier–Stokes equations (Li and Tsubokura, 2017). To simulate theairflow in the complex geometry of the human vocal tract, themoving immersed boundary method with the hierarchicalstructure grid system was applied. By analyzing thespectrogram of the generated sound, current study clarifies theeffect of tongue movement and velopharyngeal closure during
the articulation process of the sibilant /s/ production.

with a normal dentition of angle class I. The CT scan was
conducted by a 320-row Area Detector CT (20 fps; 320 slices
of 512 x 512 pixels). The subject pronounced the Japaneseword /usui/ for 1.6 s while the subject was sitting down. The
vocal tract geometry was collected from the scans of thephonation process and extracted from the images based on thebrightness values of scan, as shown in Figure 1, including theoral cavity, nasal cavity, tongue constriction andvelopharyngeal valve. The geometries used in the current study
were the duration from the end of vowel /u/ (Figure 1(a) and1(b)) to the sibilant sound /s/ (Figure 1(c) and 1(d)), which
includes the movement of the tongue elevation and the closureof the velopharyngeal valve. The cross-sectional areas of the
tongue and velum constrictions are reduced from 121 mm2 and55.6 mm2 to 12.9 mm2 and 0.01 mm2, respectively.

GOVERNING EQUATION
The governing equations are the Navier-Stokes equations,

U F1 F2 F30t x1 x2 x3(1)

where t is the time, xi indicates three directions in Cartesian
coordinate system (i = 1, 2, 3), and the conservative vector U is(2)
where the flux vectors Fi are

 ui




 ui u1  P i1   Ai1 i  1, 2,3
 ui u2  P i 2   Ai 2
Fi  
 ui u3  P i 3   Ai 3



T 
   e  P  ui   Aij u j  k

xi 


(3)

where ρ is the density, ui indicates the velocity components (i =1, 2, 3), P is the pressure, δij is the Kronecker delta, µAij is thestress term, e is the total energy, k is the thermal conductivity,
and the T is the temperature.
The current framework mainly based on the previous study(Lu et al., 2021). The LUSGS implicit algorithm and theadaptively switched time stepping scheme were applied for thetime advancement. The Roe scheme with 5th order MUSCL
was used to approximates the convective terms. And the second
order central difference method was adopted to calculate themagnitudes of the viscous terms. The immersed boundary
method was adopted for the geometry of the realistic humanvocal tract and the tongue and velum movement during thephonation process. To reduce the grid generation time for thegeometry and to provide better load balancing and higher
performance on parallelization, the hierarchical structure grid
was applied in the simulations. Using the hierarchical structure
grid system can make the working time required to build thecomputational grids shorter and simultaneously provide betterload balancing and higher performance for parallel
computations.

Figure 1: MRI images and vocal tract geometry of (a)(b) the
end of vowel /u/ to (c)(d) the sibilant sound /s/ in the phonation
process of /usui/.

COMPUTATIONAL CONDITION
After testing the time and space resolutions in the
preliminary simulation, the time step was set as 2 × 10 −6 s andthe minimum grid size was chosen to 0.05 mm near theconstriction to keep the accuracy of the immersed boundaryaround the turbulent region. The total grid number wasapproximately 1.82 × 108. The grid distribution on the midsagittal plane is shown in Figure 2(a). Although thephysiologically flow rate varies during the real phonation
process, the constant flow rate was set in this simulation to
focus on examining the effect of the movement of the speechorgans. The uniform inlet velocity was set to 1.5 m/s on thepharynx inlet (the red line shown in Figure 2(a)), which
resulted in a physiological flowrate of 450 cm3/s. Since thesibilant /s/ sound will be generated by the geometry of the
vocal tract (Yoshinaga et al., 2019), therefore, there is no voicesource set in this simulation.

VOCAL TRACT MODEL
To simulate the pronouncing process, the geometry of a
vocal tract was extracted from 4D-CT images of a subject
pronouncing the Japanese word /usui/. The subject was a 42year-old Japanese male who self-reported no speech disorders

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
frequency was 50 kHz with 256 points. The overlap ratio was
0.9 for the spectrogram calculation.

RESULT AND DISCUSSION
The flow and the pressure fluctuation field after the /u/ tothe /s/ during the phonation process at the sagittal plane areshown in Figure 3. At t = 0.01 s (Figure 3(a) and 3(d)), thevelopharyngeal valve was opened with the cross-sectional areaof 55.6 mm2 which let the airflow go to both nasal and oralcavities. Specifically, there was 35% airflow go through thevelopharyngeal valve and 65% airflow goes to the oral cavity.
Therefore, the velocity of the mainstream at the oral cavity was
only around 10 m/s. In addition, because the tongue was at thelower position, the constriction of oral cavity was not narrow
enough to generate the /s/ sound. At t = 0.05 s (Figure 3(b) and3(e)), due to the closing velopharyngeal valve, the higher flow
resistance to the upper chamber decreased the flow rate from
35% to 13% in the nasal cavity and increased the flow ratefrom 65% to 87% in the oral cavity. Furthermore, the reductionof the cross-sectional area of the oral cavity caused by theelevated tongue position leaded the velocity of the mainstream
in the oral cavity was accelerated to around 20 m/s and theperiodic pressure wave was observed near the lips in Figure3(e). Besides, because of the narrow velopharyngeal valve, thevelocity became turbulent at the nasal cavity which resulted in
the pressure wave at the exit of the nose in Figure 3(e). At t =0.10 s (Figure 3(c) and 3(f)), since the velopharyngeal valve
was completely closed, there was no nasal sound emission
observed in pressure fluctuation field. At this time, since the tipof the tongue moved posterior direction and elevated towards
the hard palate, the velocity of the mainstream in the oral cavity
reached around 25 m/s and the amplitude of the sound wasincreased at the far-field (Figure 3(f)).

Figure 2: (a) The grid distribution (every 16th gridline isshown for clarity) (b) The computational domain.

The computational domain is shown in Figure 2(b). The x1
is defined as the anterior-posterior direction; the x2 is defined as
the inferior–superior direction; the x3 is defined as thetransverse direction. To keep the flow in the computational
domain from being polluted by reflecting pressure waves at theouter computational boundary, an absorbing boundary
condition was used as the outlet condition. Fast Fouriertransform (FFT) using the Hann window was applied to thepressure waveforms sampled at 80 mm from the lip outlet toanalyze the far-field sound spectrogram. The FFT sampling

Figure 3: (a-c) Instantaneous flow field and (d-f) pressure fluctuation along the sagittal plane at (a)(d) t = 0.01 s, (b)(e) t = 0.05 s and(c)(f) t = 0.10 s.

312th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
The spectrogram of the propagating sound collected at 80mm from the lips is shown in Figure 4(a). For t = 0.01 s to t =0.03 s, the amplitude of the broadband noise around 5–9 kHz
was observed with relatively small amplitudes due to the
opening velopharyngeal valve which leaded to part of the
airflow going through the nasal cavity. Besides, since the crosssectional area of the oral tract was deceased with the tongue
elevation, the amplitude peak at 6 kHz gradually appeared aftert = 0.02 s. For t = 0.03 s to t = 0.07 s, because of the closing
velopharyngeal valve, the broadband noise was amplified in theextended frequency range of 3–11 kHz. After t = 0.08 s,
because the size of the velopharyngeal valve became smaller
than that of the tongue constriction, the airflow tended to go to
the oral cavity rather than the nasal cavity. In addition, since
the tongue constriction at this moment was narrow enough forthe production of the sibilant sound, the typical /s/ sound was
generated with the characteristics peak after 4k Hz, and thebroadband noise was observed around 4–12k Hz.
Figure 4(b) shows the measurement result of the phonation
process by the real subject (Yoshinaga et al., 2019). Although
the frequency range of the /s/ sound is different, thespectrogram during the co-articulation by the current study (t =0.01 s to t = 0.08 s) is similar to that shown in Figure 4(b).
Furthermore, the relation between the size of thevelopharyngeal valve and the sound amplitude observed at the
far-field was consistent with the previous studies (Kummer,
1992; Sundström, 2020).

incisor angle on production of sibilant/s. Scientific Reports,
11(1), 1-11.
Sundström, E., Boyce, S., and Oren, L., 2020, Effects ofvelopharyngeal openings on flow characteristics of nasal
emission. Biomechanics and modeling in mechanobiology, 113.
Yoshinaga, T., Nozaki, K., and Wada, S., 2019, A
simplified vocal tract model for articulation of [s]: The effect of
tongue tip elevation on [s]. PloS one, 14(10), e0223382.

CONCLUSION
In the current study, the sibilant /s/ phonation process was
simulated with the realistical vocal tract movement of theJapanese word /usui/. By analyzing the spectrogram, theinteraction between the sound generated by the tongue
movement and the nasal emission caused by thevelopharyngeal opening were observed. The sibilant /s/ soundwas generated gradually by the reduction of the constriction
size with the tongue elevation. However, while thevelopharyngeal valve was not completely closed, the soundwas influenced by the nasal emission. At the end of thesimulation, the typical sibilant /s/ sound was produced by the
tongue constriction when the velopharyngeal was completely
closed. These results provide clear description of the /s/ soundarticulation process which contribute to the clinical treatment
for the sibilant sound production.

REFERENCES
Bunton, K., and Story, B. H., 2012, The relation of nasality
and nasalance to nasal port area based on a computationalmodel. The Cleft palate-craniofacial journal, 49(6), 741-749.
Bunton, K., 2015, Effects of nasal port area on perceptionof nasality and measures of nasalance based on computationalmodeling.
Kummer, A. W., Curtis, C., Wiggs, M., Lee, L., and Strife,
J. L., 1992, Comparison of velopharyngeal gap size in patientswith hypernasality, hypernasality and nasal emission, or nasal
turbulence (rustle) as the primary speech characteristic. The
Cleft palate-craniofacial journal, 29(2), 152-156.
Li, C. G., and Tsubokura, M., 2017, An implicit turbulencemodel for low-Mach Roe scheme using truncated Navier–Stokes equations. Journal of Computational Physics, 345, 462474.
Lu, H., Yoshinaga, T., Li, C., Nozaki, K., Iida, A., andTsubokura, M., 2021, Numerical investigation of effects of

Figure 4. The spectrogram during co-articulation and sibilant/s/ sound (a) by the current study (b) measurement for the
subject (Yoshinaga et al., 2019).

4