## 2 Methods
### Computational Domain and Boundary Conditions  
The reconstructed model is immersed in a rectangular domain measuring 14.9 cm \(\times\) 25 cm \(\times\) 20.2 cm (Figure 3a), with the walls of the subject's face extended to the outer domain to ensure solver stability. This configuration creates a sufficiently large environment  
Figure 3: (**a**) The computational domain of the current study; (**b**) domain grid resolution; (**c**) the time history of the flow rate for a stationary vocal fold with a maximum glottal opening for current grids and finer grids measured at the mouth; (**d**) the buffer zones of the LPCE solver; (**e**) boundary conditions of the vocal fold. Figure 2: The reconstructed supraglottal tract form a MRI scan of human speaker published by [PERSON] et al. [33]. box that minimizes the influence of outer boundaries on the flow field within the airway tract and prevents reverse flow at these boundaries. Both the incompressible flow solver and the LPC acoustic solver share the same rectangular domain discretized into a non-uniform Cartesian grid of 128 \(\times\) 256 \(\times\) 128 cells, with the highest grid density concentrated around the intraglottal region (Figure 3b). To verify grid independence, a finer grid was employed for the incompressible flow, discretizing the domain into 256 \(\times\) 256 \(\times\) 256 cells with a static vocal fold at maximum glottal opening, and as illustrated in Figure 3c, the glottal volume flow rates are nearly identical between these two grid configurations. A Dirichlet pressure boundary condition of 800 Pa is applied at the airway inlet, while a Dirichlet pressure of 0 Pa is set for all outer boundaries of the environment box, as depicted in Figure 3a. To eliminate acoustic reflections from both the airway inlet and the outer boundaries of the environment box, an anechoic buffer zone [34] is enforced at the first and last ten neighboring grid points along these boundaries, as shown in Figure 3d. A non-slip, non-penetration condition is applied to all walls for the incompressible solver, whereas the LPC solver utilizes a hard wall boundary condition [23]. The density of air is set at 1.1455 kg/m\({}^{3}\), and the kinematic viscosity is established at 6.6 \(\times\) 10\({}^{-5}\) m\({}^{2}\)/s, corresponding to approximately one-fourth of the Reynolds number for normal human phonation, thereby reducing computational costs. This choice primarily impacts the turbulent flow in the supraglottal tract, as the influence of normal pressure on vocal fold vibration significantly outweighs that of viscous pressure, which is less critical in studies of normal phonation. The vocal folds are discretized into 20,643 tetrahedron elements to accurately capture the complex flow-induced motion during phonation, with a tissue density set at 1.04 \(\times\) 10\({}^{3}\) kg/m\({}^{3}\). Like other state-of-the-art FEM simulations, our model faces inherent challenges in fully representing the physical exchange of structural forces between the vocal folds during contact. To address this limitation, we implement a rigid wall contact model to prevent penetration of the vocal folds, while maintaining a small artificial gap of 0.2 mm between them during closure. While this approach is a numerical approximation of the actual contact mechanics, it enables the successful implementation of the ghost cell methodology and ensures numerical stability. This small gap does result in slight flow leakage and a non-zero glottal area during glottal closure, which are 21% and 20% of the maximum flow rate and glottal area, respectively. Traction boundary conditions are applied at the fluid-structure interaction interface on the free surface, as depicted in Figure 3e. Additionally, zero displacement boundary conditions are enforced on the fixed surfaces, ensuring stable interactions and accurate simulations of vocal fold dynamics.