###### Abstract
The high variation in the speed of sound is based on the density variation of the fluid mixture and lowers the convergence rate of numerical solutions, because it depends on the lowest speed of sound. In order to solve the problem, the preconditioning method first proposed by [PERSON] [10] is frequently employed. In the method, the eigenvalues of the flux matrix are modified to reduce the variation in the speed of sound over the entire computational domain [11]. In the preconditioning schemes for steady solutions, the time terms play the role of an iteration index. With respect to unsteady simulations, dual-time stepping methods that use pseudo time in the governing equations are employed [12]. When the residual of the dependent variables is lower than the prescribed criteria in the pseudo time step, the numerical solution simulation proceeds to the next step with respect to physical time. The numerical methods were used in preceding studies, including [PERSON] et al. [13; 14], [PERSON] et al. [15; 16], and [PERSON] et al. [17; 18]. The studies developed numerical methods based on dual-time stepping preconditioning techniques to solve incompressible Navier-Stokes equations, and validated the same by comparing their numerical solutions with experimental results [19]. [PERSON] et al. [20] extended the methods to solve fluid-structure-interaction problems by combining six degree of freedom (DOF) models with the incompressible code and predicting the motion of high speed supercavitating vehicles. [PERSON] et al. [21; 22], [PERSON] et al. [23], and [PERSON] and [PERSON] [24] added isothermal models to improve the accuracy of a numerical solution for cavitation flow. [PERSON] et al. [25; 26], [PERSON] et al. [27; 28; 29], and [PERSON] et al. [30] extended incompressible numerical methods to consider compressibility by solving the compressible Navier-Stokes equations, and attempted to validate the compressible codes by applying the same to various cavitation flows. [PERSON] et al. [31; 32] also developed numerical methods to more accurately solve compressible cavitation flows by investigating the effects of turbulence models and viscous flux vectors. Essentially, the numerical method is employed in the current study. As described above, several numerical techniques were developed to solve physical issues related to the accurate and efficient simulation of cavitation flow. However, the numerical simulation of cavitation flows requires significantly higher computational times due to increases in the number of equations and increases in the variation in the speed of sound when compared with those of single-phase fluid flow. Therefore, several studies have focused on resolving the issue from a different viewpoint that involved adopting parallel computing techniques [33; 34]. Domain decomposition methods are generally used to apply the parallel computing techniques for the simulation of cavitation flow, and this divides the complete computational domain into subdomains that are separately treated by each processor. This implies that increases in the number of processors increase the number of subdomains and decrease the computation time, because the size of the domain for each processor decreases. However, the simulation time does not linearly decrease with increases in the number of processors, because the interface boundary between subdomains increases. Specifically, the treatment of the interface boundary between the subdomains is critical for the accuracy and efficiency of the parallel computation. The interface boundary is categorized into non-overlapping and overlapping methods. In the non-overlapping method, a simple linear interpolation method is applied for information exchange at the interface boundary between subdomains. The method is fast and simple. However, error due to low-order interpolation can adversely affect the accuracy of the numerical result if the interface is subject to a high gradient flow that cannot be avoided when the number of subdomains increases.