## 3 Finite Element Model
### Finite Element Approximation  
The displacement field is approximated with a linear combination of piece-wise linear basis functions \(\{\phi_{i}\}_{i=1}^{N}\), given by  
\[\mathbf{u}(\mathbf{x},t)\approx \begin{pmatrix}\sum_{i=1}^{N}\alpha_{i}(t)\phi _{i}(\mathbf{x})\\ \sum_{i=1}^{N}\beta_{i}(t)\phi_{i}(\mathbf{x})\end{pmatrix}  \tag{6}\]  
where each \(\phi_{i}(\mathbf{x})\) is defined over a set of connected elements that approximate the geometry of the VFs. The FE approximation is then found by substituting this approximation into Equation (2), giving  
\[\mathbf{M}\ddot{\mathbf{\theta}}+\mathbf{D}\dot{\mathbf{\theta}}+\mathbf{K}\dot{\mathbf{ \theta}}=\mathbf{F}, \tag{7}\]  
where  
\[\mathbf{M}=\int_{\Omega}\rho\mathbf{\Phi}^{\mathrm{T}}\mathbf{\Phi}\;\mathrm{d}V,\ \ \ \mathbf{D}=\int_{\Omega}\eta\mathbf{\Phi}_{d}^{\mathrm{T}}\mathbf{S}\mathbf{\Phi}_{d}\; \mathrm{d}V,\ \ \ \mathbf{K}=\int_{\Omega}\mathbf{\Phi}_{d}^{\mathrm{T}}\mathbf{C}\mathbf{\Phi}_{d}\; \mathrm{d}V,\ \ \ \text{and}\ \ \ \mathbf{F}=\int_{\Gamma}\mathbf{\Phi}^{\mathrm{T}}\mathbf{f}_{s}\; \mathrm{d}S. \tag{8}\]  
Here,  
\[\mathbf{\theta}(t) = \left[\alpha_{1}(t),\ldots\alpha_{N}(t),\beta_{1}(t),\ldots\beta _{N}(t)\right]^{\mathrm{T}}, \tag{9}\] \[\mathbf{\Phi}(\mathbf{x}) = \begin{pmatrix}\phi_{1}(\mathbf{x})&\ldots&\phi_{N}(\mathbf{x})&0 &\ldots&0\\ 0&\ldots&0&\phi_{1}(\mathbf{x})&\ldots&\phi_{N}(\mathbf{x})\end{pmatrix} ,\] (10) \[\mathbf{\Phi}_{d}(\mathbf{x}) = \begin{pmatrix}\partial_{x}\phi_{1}(\mathbf{x})&\ldots&\partial_{ x}\phi_{N}(\mathbf{x})&0&\ldots&0\\ 0&\ldots&0&\partial_{z}\phi_{1}(\mathbf{x})&\ldots&\partial_{z}\phi_{N}( \mathbf{x})\\ \partial_{z}\phi_{1}(\mathbf{x})&\ldots&\partial_{z}\phi_{N}(\mathbf{x})& \partial_{x}\phi_{1}(\mathbf{x})&\ldots&\partial_{x}\phi_{N}(\mathbf{x}) \end{pmatrix} , \tag{11}\]  
and integration is interpreted element wise. The computation of \(\mathbf{M}\), \(\mathbf{D}\), and \(\mathbf{K}\) now follow the standard finite element construction [28]. The boundary of the VFs was split into two sections: (_i_) a section where external forces are applied, denoted \(\Gamma_{\text{free}}\), and (_ii_) a section where the nodes do not move, which represents the walls of the trachea, denoted \(\Gamma_{\text{fixed}}\). The vector \(\mathbf{F}\) is the nodal force vector and is the result of the transverse surface force. In this work aerodynamic pressure is treated as the only stress leading to an external force on the surface nodes. The pressure at each node along the free surface is modeled using a one-dimensional Bernoulli flow model, which gives  
\[p(s,t)=\left\{ \begin{array}{rl}p_{\text{sub}}-(p_{\text{sub}}-p_{\text{sup }})\left(\frac{A_{\text{sep}}(t)}{A(s,t)}\right)^{2},&A(s,t)<A_{\text{sep}}\\ p_{\text{sup}},&A(s,t)\geq A_{\text{sep}}\end{array} \right. \tag{12}\]where \(s\in\Gamma_{\text{free}}\), \(p_{\text{sub}}\) and \(p_{\text{sup}}\) are the subglottal and supraglottal pressures, respectively, \(A_{\text{min}}(t)\) is the minimum glottal area at time \(t\), \(A(s,t)\) is the area at location \(s\) at time \(t\), and \(A_{\text{sep}}=1.3A_{\text{min}}(t)\) is the glottal area at the location of flow separation [29; 30]. Vocal fold collision is modeled by restricting the \(x\)-coordinate of each node from crossing the midline. Should a time step result in a node crossing the midline, its \(x\)-coordinate is forced to the midline. The force required to restrict the node to the midline is directly the contact force. Contact pressure can be obtained by dividing by the edge area (distance between the the node and its neighbor). Finally, time integration was achieved through the finite differencing scheme set out by [PERSON] et al. [27]. The time integration gives the displacement field over time, yielding a position time series for each node in the finite element mesh.