Skip to content

Human Phonation

In this example, we will simulate the acoustic wave propagation inside a 3D model of the human vocal tract. This example is based on the simulation workflow that was used in the simVoice joint project of FAU Erlangen, TU Wien and TU Graz 1, 2

Thereby the acoustic sources for the used PCWE are computed from the incompressible pressure field that was obtained by a CFD simulation of the flow 3, 4, 5. Due to the size of the input data, the source files for this example are not provided in the documentation. Please contact Andreas Wurzinger instead.


Mesh

The acoustic Mesh was created with the commercial software Ansys ICEM. It consists 4 regions: The larynx region LARYNX, the vocal tract VT, the propagation region PR, and the perfectly matched layer region PML.

CAA Setup


XML flow pressure Interpolation

In order to perform a acoustic simulation on the presented mesh, the acoustic source data has to be interpolated onto it. Therefore, CFSdat (a subpackage included in openCFS) is used.

As openCFS does, CFSdat uses an xml file to configure the workflow.

  • Basic structure including XML header
<?xml version="1.0" encoding="UTF-8"?>
<cfsdat xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.cfs++.org/simulation">
    <pipeline>

        USERINPUT

    </pipeline>
</cfsdat>
  • Define data length

First we define the length and timestep size of the CFD data. We simulate one full cycle of the vocal folds' motion, which are osscilation at f_0=148 \mathrm{Hz}. Therefore, we define numsteps= 675 timesteps with a stepsize of \Delta t=1\cdot 10^{-5} s.

<stepValueDefinition>
    <startStop>
        <startStep value="0"/>              <!-- same as the first Ensight step number, but ATTENTION consider delta t !!! -->
        <numSteps value="675"/>           <!-- step numbers to calculate -->
        <startTime value="1e-05"/>          <!-- start from the time that is equal to the 7501th step from the Ensight -->
        <delta value="1e-05"/>              <!-- delta t -->
        <deleteOffset value="no"/>          <!-- if NO then the time will be the same as in the CFD (Ensight) -->
    </startStop>
</stepValueDefinition>
  • Define mesh input

The CFD data was exported in Ensight-Gold format as the file: 'CFD_results.case'.

<meshInput id="input" gridType="fullGrid">
    <inputFile>
        <ensight fileName="cfd_data_path/CFD_results.case">
            <variableList>
                <variable CFSVarName="fluidMechPressure" EnsightVarName="Pressure"/>
            </variableList>
        </ensight>
    </inputFile>
</meshInput>
  • Conservative field interpolation

For the field interpolation the conservative cut-cell volume approach is used to interpolate the cell data onto the acoustic mesh saved in the hdf5 file 'CAA_mesh.h5'.

<interpolation type="FieldInterpolation_Conservative_CutCell" id="interp" inputFilterIds="input">
    <targetMesh>
        <hdf5 fileName="CAA_mesh.h5"/>
    </targetMesh>
    <singleResult>
        <inputQuantity resultName="fluidMechPressure"/>
        <outputQuantity resultName="fluidMechPressure_interpolated"/>
    </singleResult>
    <regions>
        <sourceRegions>
            <region name="Background"/>
        </sourceRegions>
        <targetRegions>
            <region name="LARYNX"/>
            <region name="VT"/>
        </targetRegions>
    </regions>
</interpolation>
  • Define output

Finally we define the output filename by meshOutput id to save the interpolated pressure field in for both reagions LARYNX and VT.

<meshOutput id="InterpolatedPressureField" inputFilterIds="interp">
    <outputFile>
        <hdf5 extension="cfs" compressionLevel="6" externalFiles="no"/> <!-- compression level (greater number = bigger compression; default=1) -->
    </outputFile>
    <saveResults>
        <result resultName="fluidMechPressure_interpolated">
            <allRegions/>
        </result>
    </saveResults>
</meshOutput>
  • Run CFSdat
cfsdat interpolatePressure

XML Source term computation

As we now have obtained the conservatively interpolated incompressible flow pressure on the acoustic grid, we again use CFSdat to compute the RHS source of the PCWE \partial p / \partial t. The convective part of the RHS source is neglected in this example due to its negligible impact on the resulting sound spectrum in the far field (see references).

The xml file is structured as before with the inputs:

  • Define data length

We define the length and timestep size of the CFD data as done for the flow pressure Interpolation.

<stepValueDefinition>
    <startStop>
        <startStep value="0"/>              <!-- same as the first Ensight step number -->
        <numSteps value="675"/>             <!-- step numbers to calculate -->
        <startTime value="1e-05"/>          <!-- start from the time that is equal to the *startStep*-th step from the Ensight -->
        <delta value="1e-05"/>              <!-- delta t -->
        <deleteOffset value="no"/>          <!-- if NO then the time will be the same as in the CFD (Ensight) -->
    </startStop>
</stepValueDefinition>
  • Define input

This time we define the previously computed interpolated flow pressure data.

<meshInput id="input" gridType="fullGrid">
    <inputFile>
        <hdf5 fileName="results_hdf5/InterpolatedPressureField.cfs"/>
    </inputFile>
</meshInput>
<timeDeriv1 inputFilterIds="input" id="pressure_time_derivative">
    <singleResult>
        <inputQuantity resultName="fluidMechPressure_interpolated"/>
        <outputQuantity resultName="acouRhsLoad"/>
    </singleResult>
</timeDeriv1>
  • Define output
<meshOutput id="source_dpdt" inputFilterIds="pressure_time_derivative">
    <outputFile>
        <hdf5 extension="cfs" compressionLevel="6" externalFiles="no"/><!-- compression level (greater number = bigger compression; default=1) -->
    </outputFile>
    <saveResults>
        <result resultName="acouRhsLoad">
            <allRegions/>
        </result>
    </saveResults>
</meshOutput>

XML Simulation Setup

For the propagation simulation

  • Basic structure including XML header
<?xml version="1.0"?>
<cfsSimulation xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xmlns="http://www.cfs++.org/simulation">

    FILE_FORMATS

    DOMAIN

    <sequenceStep>

        ANALYSIS_TYPE

        <pdeList>

            PDE

        </pdeList>

        SOLVER_SETTINGS


    </sequenceStep>
  • File formats
<fileFormats>
    <input>
        <hdf5 fileName="source_dpdt.cfs"/>
    </input>
    <output>
        <hdf5 id="h5"/>
        <text id="txt"/>
    </output>
    <materialData file="mat.xml" format="xml"/>
</fileFormats>
  • Computational domain

We define the material for all regions aswell as the non-conforming interface. Additionally, we define a microphone point at which we want to obtain the acoustic pressure for postprocessing.

<domain geometryType="3d">
    <regionList>
        <region name="LARYNX" material="air"/>
        <region name="VT" material="air"/>
        <region name="PR" material="air"/>
        <region name="PML" material="air"/>
    </regionList>
    <ncInterfaceList>
        <ncInterface name="IF1" masterSide="IF_VT" slaveSide="IF_PR"/>
    </ncInterfaceList>
    <nodeList>
        <nodes name="mic">
            <coord x="0.24657" y="0.009" z="-0.049069"/>
        </nodes>
    </nodeList>
</domain>
  • Analysis Type

We compute a transient simulation of one full cycle of the vocal folds' motion, which are osscilation at f_0=148 \mathrm{Hz}. Therefore, we define numsteps= 675 timesteps with a stepsize of \Delta t=1\cdot 10^{-5} s.

<analysis>
    <transient>
        <numSteps>675</numSteps> 
        <deltaT>1e-5</deltaT>
    </transient>
</analysis>

Including definition of PML regions, absorbing boundary conditions, non-conforming interfaces and Temporal Blending of the RHS source with the blending function f(t)

\begin{equation} f(t) = \left\{ \begin{array}{lll} -\frac{1}{\rho_{0} c^2}\left[1-cos\left(\frac{0.5 \pi}{8\cdot 10^{-4}} t\right)\right] & \mathrm{if} & t < 8\cdot 10^{-4} s\\ -\frac{1}{\rho_{0} c^2} & \mathrm{else} \end{array}\right. \end{equation}

with the density \rho_{0} = 1.204 \frac{kg}{m^3} and the speed of sound c = 343.4 \frac{m}{s}.

As output, we define the acouPotentialD1 \frac{\partial \psi_a}{\partial t }, saved both for the whole domain as .cfs output and at the microphone point as .txt ouput.

<acoustic formulation="acouPotential" timeStepAlpha="-0.3">
    <regionList>
        <region name="LARYNX"/>
        <region name="VT"/>
        <region name="PR"/>
        <region name="PML" dampingId="dampPML2"/>
    </regionList>
    <ncInterfaceList>
        <ncInterface name="IF1" formulation="Nitsche" nitscheFactor="50"/>
    </ncInterfaceList>
    <dampingList>
        <pml id="dampPML2">
            <propRegion>
                <direction comp="x"     min="0.19"      max="0.275"/>
                <direction comp="y"     min="-0.051"    max="0.069"/>
                <direction comp="z"     min="-0.0525"   max="0.0675"/>
            </propRegion>
            <type>inverseDist</type>
            <dampFactor>1.0</dampFactor>
        </pml>
    </dampingList>
    <bcsAndLoads>
        <absorbingBCs volumeRegion="LARYNX" name="IF_ABC"/>
        <rhsValues name="LARYNX">
            <grid>
                <defaultGrid quantity="acouRhsLoad" dependtype="GENERAL" >
                    <globalFactor> ((t lt 8e-4)? (-1.0)/(1.204*343.4*343.4)*(1*(1-(cos(0.5*pi/8e-4*t))^2)) : (-1.0)/(1.204*343.4*343.4))</globalFactor>
                </defaultGrid>
            </grid>
        </rhsValues>
        <rhsValues name="VT">
            <grid>
                <defaultGrid quantity="acouRhsLoad" dependtype="GENERAL">
                    <globalFactor> ((t lt 8e-4)? (-1.0)/(1.204*343.4*343.4)*(1*(1-(cos(0.5*pi/8e-4*t))^2)) : (-1.0)/(1.204*343.4*343.4))</globalFactor>
                </defaultGrid>
            </grid>
        </rhsValues>
    </bcsAndLoads>
    <storeResults>
        <nodeResult type="acouPotentialD1">
            <allRegions/>                   
            <nodeList>
                <nodes name="mic" outputIds="txt"/>
            </nodeList>           
        </nodeResult>                
    </storeResults>
</acoustic>
  • Solver Settings
<linearSystems>
    <system>
        <solverList>
            <pardiso id="default">
            </pardiso>
        </solverList>
    </system>
</linearSystems>
  • Run the simulation
cfs propagation

Results

As a simulation result we obtain the acoustic pressure field by multiplying the time derivative of the scalar acoustic potential with the density

p_a = \rho_0 \frac{\partial \psi_a}{\partial t }

which can be visualized with Paraview

Acoustic Pressure

\frac{\partial \psi_a}{\partial t } in the (x,y) plane.

Acoustic Pressure2

Contour planes at \frac{\partial \psi_a}{\partial t } = 5 \frac{m^2}{s^2}.

Further interpretations can be done by calculating the Amplitude Spectral Density of the acoustic pressure at the microphone point and plot the spectral distribution of the sound pressure level (SPL). The following figure compares the simulation data of 20 vocal fold oscillation cycles with a measurement signal.

SPL


References


  1. Stefan Schoder, Michael Weitz, Paul Maurerlehner, Alexander Hauser, Sebastian Falk, Stefan Kniesburges, Michael Döllinger, and Manfred Kaltenbacher. Hybrid aeroacoustic approach for the efficient numerical simulation of human phonation. The Journal of the Acoustical Society of America, 147(2):1179–1194, 2020. 

  2. Stefan Schoder, Paul Maurerlehner, Andreas Wurzinger, Alexander Hauser, Sebastian Falk, Michael Döllinger, and Manfred Kaltenbacher. Aeroacoustic sound source characterization of the human voice production - perturbed convective wave equation. Applied Sciences, 2021. 

  3. Hossein Sadeghi, Stefan Kniesburges, Manfred Kaltenbacher, Anne Schützenberger, and Michael Döllinger. Computational models of laryngeal aerodynamics: potentials and numerical costs. Journal of Voice, 33(4):385–400, 2019. 

  4. Hossein Sadeghi, Michael Döllinger, Manfred Kaltenbacher, and Stefan Kniesburges. Aerodynamic impact of the ventricular folds in computational larynx models. The Journal of the Acoustical Society of America, 145(4):2376–2387, 2019. 

  5. Hossein Sadeghi, Stefan Kniesburges, Sebastian Falk, Manfred Kaltenbacher, Anne Schützenberger, and Michael Döllinger. Towards a clinically applicable computational larynx model. Applied Sciences, 9(11):2288, 2019.