## 3 Finite Element Model
### Mathematical Model  
Small deformations of the silicone VFs can be modeled using a displacement field \(\mathbf{u}\), which describes the displacement of two linear transversely isotropic elastic bodies [27]. The VFs were oriented such that the medial-lateral direction defines the \(x\)-axis, the \(y\)-axis is in the anterior-posterior direction, and the \(z\)-axis is positive in the superior direction. As a simplification, the dynamics of the two VFs are treated as symmetric and displacement is treated as uniform along the length of the VFs; as such, the displacement field does not vary along the \(y\)-axis, i.e.,  
\[\mathbf{u}=\mathbf{u}(\mathbf{x},t)=u_{x}(\mathbf{x},t)\hat{i}+u_{z}(\mathbf{ x},t)\hat{j}, \tag{1}\]  
where \(\mathbf{x}=[x,\ z]^{\mathrm{T}}\) is a vector of planar spatial coordinates. This reduction of physical dimension greatly reduces the computational cost of the FE model.  
Given a deformation with displacement field \(\mathbf{u}\) and an arbitrary variation \(\mathbf{u}^{\delta}\), the virtual work principle gives [28]  
\[\int_{\Omega}\left(\mathbf{u}^{\delta}\right)^{\mathrm{T}}\rho(\mathbf{x}) \hat{\mathbf{u}}\ \mathrm{d}V+\int_{\Omega}\boldsymbol{\sigma}(\mathbf{u}): \boldsymbol{\varepsilon}(\mathbf{u}^{\delta})\ \mathrm{d}V=\int_{\Gamma}\left(\mathbf{u}^{\delta}\right)^{ \mathrm{T}}\mathbf{f}_{\mathrm{s}}\ \mathrm{d}S, \tag{2}\]  
where \(\Omega\) and \(\Gamma\) are the 2D VF domain and its boundary, respectively, \(\rho\) is the material density, \(\boldsymbol{\sigma}\) and \(\boldsymbol{\varepsilon}\) are the stress and strain tensors, and \(\mathbf{f}_{\mathrm{s}}\) is the surface force, in this case arising from aerodynamics.  
Since the vocal folds are modeled as isotropic in the \(xz\)-plane, [PERSON]'s law gives [28]  
\[\boldsymbol{\sigma}(\mathbf{u}):\boldsymbol{\varepsilon}(u^{\delta})= \boldsymbol{\tilde{\sigma}}(\mathbf{u})^{\mathrm{T}}\boldsymbol{\tilde{ \varepsilon}}(u^{\delta})=\boldsymbol{\tilde{\varepsilon}}(u^{\delta})^{ \mathrm{T}}\boldsymbol{\tilde{\sigma}}(\mathbf{u})=\boldsymbol{\tilde{ \varepsilon}}(u^{\delta})^{\mathrm{T}}\mathbf{C}\boldsymbol{\tilde{\varepsilon }}(\mathbf{u}), \tag{3}\]  
where \(\boldsymbol{\tilde{\varepsilon}}(\mathbf{u})=\left[\frac{\partial u_{x}}{ \partial x},\ \frac{\partial u_{z}}{\partial z},\ \frac{\partial u_{z}}{\partial x}+\frac{\partial u_{x}}{ \partial z}\right]^{\mathrm{T}}\) and \(\mathbf{C}\) is a positive definite matrix given by  
\[\mathbf{C}(E,\
u)=\begin{pmatrix}\lambda+2\mu&\lambda&0\\ \lambda&\lambda+2\mu&0\\ 0&0&\mu\end{pmatrix}=\mu(E,\
u)\begin{pmatrix}\frac{2(1-\
u)}{1-2\
u}&\frac{2 \
u}{(1-2\
u)}&0\\ \frac{2\
u}{(1-2\
u)}&\frac{2(1-\
u)}{1-2\
u}&0\\ 0&0&1\end{pmatrix}. \tag{4}\]Here \(\lambda\) and \(\mu\) are the Lame parameters that are related to the [PERSON]'s modulus \(E\) and [PERSON]'s ratio \(\
u\) via  
\[\mu(E,\
u)=\frac{E}{2(1+\
u)},\ \ \ \lambda(E,\
u)=\frac{\
u E}{(1+\
u)(1-2\
u)}. \tag{5}\]  
Note that for \(\mathbf{C}(E,\
u)\) to be positive definite requires \(E>0\) and \(-1<\
u<1/2\).  
Finally, following [PERSON] et al. [27], the viscous damping experienced by the VFs is modeled by replacing the shear modulus, \(\mu(E,\
u)\), with \(\mu+\eta\mathrm{d}/\mathrm{d}t\), where \(\eta\) is viscosity.