## 4 Discussion
### Comparison between the Spectral Method Algorithm and the COMSOL Simulation
In order to find out the influence of mesh grids on the COMSOL simulation results, five different meshes were used to calculate the mode wavenumber, as shown in Figure 3:3  
The wave number \(k_{z}\) of the axial propagation modes in the water-filled PMMA tube at \(f=20\) kHz was calculated using COMSOL with five different mesh files, as shown in in Figure 3:4 The differences in \(k_{z}\) at each mode between mesh 1-4 and mesh 5 are plotted in Figure 3:5 The horizontal axis is the mode index, meaning that there exist 22 modes at \(f=20\) kHz. The vertical axis is the absolute value of the difference between mesh 1-4 and mesh 5 at each mode. Figure 5 shows that the simulation result \(k_{z}\) of each mode tends to be consistent with the increase in the number of grid elements. In order to achieve the most accurate calculation results, the following simulation used mesh 5. The mode analysis in COMSOL can calculate the wave number and mode shape of a propagation mode at a specific frequency. Combined with sweep analysis, dispersion curves of the axial propagation mode of the water-filled PMMA tube can be obtained. The frequency range calculated here is 0-90 kHz and the frequency interval is 200 Hz. Figure 6 shows the comparison of the dispersion curves obtained from the COMSOL simulation and the spectral method algorithm. \begin{table}
\begin{tabular}{c c c c c c} \hline \hline
**Mesh** & **Mesh 1** & **Mesh 2** & **Mesh 3** & **Mesh 4** & **Mesh 5** \\ \hline Element Number & 2214 & 3146 & 4644 & 8698 & 27,842 \\ \hline \hline \end{tabular}
\end{table} 
Table 2: COMSOL mesh file and their numbers of elements. Figure 5: The difference of \(k_{z}\) between mesh 1–4 and mesh 5 at each mode. In Figure 6, the horizontal axis is the dimensionless parameter \(k_{1}b\), as shown in Equation (82), and the vertical axis is the normalized phase velocity \(c_{p}/c_{1}\). The vertical axis uses logarithmic coordinates. \[k_{1}b=\frac{\omega}{c_{1}}b \tag{82}\]  
The curves in Figure 6 represent the calculation results from the spectral method algorithm and the dots represent the COMSOL simulation results. It can be seen from Figure 6 that the dispersion curves of the axial propagation modes in the water-filled PMMA tube obtained from the spectral method algorithm were completely consistent with the COMSOL simulation results. Table 3 shows the time consumed by the COMSOL Multiphysics 5.6 software and the spectral method algorithm implemented in MATLAB R2022b to obtain the dispersion curves in Figure 6. The computer used here to run the COMSOL simulation and the MATLAB program of the proposed algorithm was a Lenovo desktop workstation, which has 2 Intel(r) Xeon(r) Silver 4110 CPU, and 128 GB RAM. The spectral method algorithm not only provides an accurate calculation but requires much less time and CPU resources to do so than the COMSOL simulation software. Table 4 compares the wavenumbers \(k_{z}\) of the axial propagation modes in the water-filled PMMA tube at \(f=8\) kHz provided by the COMSOL simulation and the spectral method algorithm. The error in the table was calculated using Equation (83), where \(k_{z(s)}\) represents the axial wavenumber calculated using the spectral method and \(k_{z(c)}\) represents the axial wavenumber calculated using the COMSOL simulation. \[E=\frac{\left|k_{z(s)}-k_{z(c)}\right|}{k_{z(c)}} \tag{83}\] \begin{table}
\begin{tabular}{c c} \hline \hline
**Calculation Methods** & **Calculation Time** \\ \hline COMSOL Multiphysics & Approximately 10 h \\ The spectral method algorithm & A few tens of seconds \\ \hline \hline \end{tabular}
\end{table} 
Table 3: The calculation time of the COMSOL simulation and the spectral method algorithm. Figure 6: Comparison of the dispersion curves obtained using the COMSOL simulation and the spectral method algorithm. It can be seen from Table 4 that both the COMSOL simulation and the spectral method algorithm can find 12 different modes at 8 kHz, and the axial wavenumber corresponding to the same mode is identical under two different calculation ways.