Authors: Hao Wen Yen-Lin Chen Kayhan Burak Yıldırı

Supervisor: Sebastian Resch-Schopper  Lukas Kleine-Wächter

1. Introduction

Wave propagation in periodic structures is commonly analyzed by using Wave Finite Element Method(WFEM). This method leads to a dispersion eigenvalue problem(Equation 1.1).

Eigensolutions align into frequency bands according to frequency and wave characteristics. However, solvers return eigenvalues in ascending order without considering the eigenvectors that define mode shapes. As a result, frequency bands may cross or overlap, leading to mismatched eigenmodes and their corresponding frequency bands.

Figure 1.1: Dispersion curves obtained from MATLAB, without any corrections.

This misalignment complicates the understanding of mode evolution and affects the efficiency of dispersion curve computation.

To address this issue, two potential solutions are considered:

1. Matching of eigenvalues:

This method aims to track eigenmodes accurately by establishing correct connections between eigenvalues between phase steps using cost matrices and optimization algorithm Hungarian Algorithm.

1. Mode polarization: 

This method uses the directional distribution of kinetic energy to distinguish wave modes based on their polarization characteristics.

By implementing these methods, not only WFEM can further improve the accuracy and efficiency of dispersion analysis but also provides deeper insights into wave propagation mechanisms in materials.

Figure 1.2: Infinitely extended homogeneous 2D element for analysis.

For the implementation of methods and analysis, bm-mfem MATLAB FEM code from Chair of Structural Mechanics is modified and utilized with the model shown in Figure 1.2.

2. Methodology

2.1. Matching of Eigenvalues

The approach presented in [1] demonstrates a novel technique for matching eigensolutions. This approach aims to find out the correct match of eigenvalues across the intersections of dispersion curves such that the correct eigenmodes could be restored. 

It is suggested that the information of all possible matches could be obtained through the cost matrix $//<![CDATA[ \begin{array}{l}D\end{array} //]]>$,

where $//<![CDATA[ \begin{array}{l}\lambda_j(\mu_i)\end{array} //]]>$indicates the eigenvalue $//<![CDATA[ \begin{array}{l}\lambda\end{array} //]]>$of $//<![CDATA[ \begin{array}{l}j^t^h\end{array} //]]>$mode at $//<![CDATA[ \begin{array}{l}i^t^h\end{array} //]]>$phase $//<![CDATA[ \begin{array}{l}\mu\end{array} //]]>$, and $//<![CDATA[ \begin{array}{l}u_j(\mu_i)\end{array} //]]>$represents the eigenvector $//<![CDATA[ \begin{array}{l}u\end{array} //]]>$of $//<![CDATA[ \begin{array}{l}j^t^h\end{array} //]]>$mode at $//<![CDATA[ \begin{array}{l}i^t^h\end{array} //]]>$phase $//<![CDATA[ \begin{array}{l}\mu\end{array} //]]>$.

If the first $//<![CDATA[ \begin{array}{l}m\end{array} //]]>$modes are considered and the interested range of phase is divided to $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$reference points, there will be $//<![CDATA[ \begin{array}{l}n-1\end{array} //]]>$ cost matrices of dimension $//<![CDATA[ \begin{array}{l}m \times m\end{array} //]]>$, then it could be assumed that $//<![CDATA[ \begin{array}{l}i = [1,n]\end{array} //]]>$and $//<![CDATA[ \begin{array}{l}k = i+1\end{array} //]]>$, which means that $//<![CDATA[ \begin{array}{l}i^t^h\end{array} //]]>$cost matrix represent the state between $//<![CDATA[ \begin{array}{l}\mu_i\end{array} //]]>$and $//<![CDATA[ \begin{array}{l}\mu_i_+_1\end{array} //]]>$. After these assumption, a more convenient expression of cost matrix could be expressed as 

In equation 2.2, the first term, $//<![CDATA[ \begin{array}{l}|\lambda_j(\mu_i)-\lambda_l(\mu_i_+_1)|\end{array} //]]>$, represents the difference between the two eigenvalues, and the second term, $//<![CDATA[ \begin{array}{l}\text{min}(||u_j(\mu_i)-u_l(\mu_i_+_1)||,||u_j(\mu_i)+u_l(\mu_i_+_1)||)\end{array} //]]>$, accounts for the difference between the two eigenvectors. The reason of using the $//<![CDATA[ \begin{array}{l}\text{min()}\end{array} //]]>$function is that considering the possibility of sign inversion. The weight $//<![CDATA[ \begin{array}{l}w (w>0)\end{array} //]]>$is the relative importance of the second term with respect to the first term.

For example, if the dispersion curves are far away with each other, the first term is large enough to express that the combination of $//<![CDATA[ \begin{array}{l}j^t^h\end{array} //]]>$mode at $//<![CDATA[ \begin{array}{l}i^t^h\end{array} //]]>$phase and $//<![CDATA[ \begin{array}{l}l^t^h\end{array} //]]>$mode at $//<![CDATA[ \begin{array}{l}(i+1)^t^h\end{array} //]]>$phase is not the correct match. However, when the dispersion curves are getting close to each other, the first term become too small to distinguish the difference. If the first term is smaller than a specific value of threshold, it could imply that the situation when the intersection occurs. The multiplication with $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$could increase the importance of the second term, for that the second term is usually much smaller than the first term even though the first term is already small enough to pass the threshold. 

Therefore, $//<![CDATA[ \begin{array}{l}D^i^,^i^+^1_j_,_l\end{array} //]]>$, the element at the coordinate ( $//<![CDATA[ \begin{array}{l}j,l\end{array} //]]>$) in the $//<![CDATA[ \begin{array}{l}i^t^h\end{array} //]]>$cost matrix, represents the gap between the $//<![CDATA[ \begin{array}{l}j^t^h\end{array} //]]>$mode at $//<![CDATA[ \begin{array}{l}\mu_i\end{array} //]]>$and the $//<![CDATA[ \begin{array}{l}l^t^h\end{array} //]]>$mode at $//<![CDATA[ \begin{array}{l}\mu_i_+_1\end{array} //]]>$. In other words, the $//<![CDATA[ \begin{array}{l}i^t^h\end{array} //]]>$cost matrix provides the differences of all possible combination across $//<![CDATA[ \begin{array}{l}\mu_i\end{array} //]]>$and $//<![CDATA[ \begin{array}{l}\mu_i_+_1\end{array} //]]>$. For instance, $//<![CDATA[ \begin{array}{l}D^2^2^,^2^3_2_,_3\end{array} //]]>$represents the difference between the $//<![CDATA[ \begin{array}{l}2^n^d\end{array} //]]>$eigenmode at $//<![CDATA[ \begin{array}{l}\mu_2_2\end{array} //]]>$and the $//<![CDATA[ \begin{array}{l}3^r^d\end{array} //]]>$eigenmode at $//<![CDATA[ \begin{array}{l}\mu _2_3\end{array} //]]>$.

It is required to use $//<![CDATA[ \begin{array}{l}n-1\end{array} //]]>$of cost matrices to analyze all potential combination of eigenmodes throughout $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$phases. Afterward, the correct eigenmode tracing throughout $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$phases could be achieved by applying Hungarian Algorithm [2].

Besides, it could be derived from [3] that by the approximate orthogonality of the eigenvectors of two neighboring dispersion curves, the eigenmodes could be sorted when intersection of dispersion curves occurs. Since the eigenvector of different eigenmode should be orthogonal to each other, the approximate orthogonality could offer a criteria to determine whether the eigenvectors are belong to the same pattern of motion, i.e. eigenmode. The equation could be expressed as the following:

Ideally, the inner product of eigenvectors will be either 0 or 1 if normalization of eigenvectors is correctly conducted due to orthogonality. However, since the dimension of the eigenvectors might be large, conducting normalization could deteriorate the computation efficiency. Thus in this project, the normalization process is not considered.  

2.1.1. Hungarian Algorithm

Hungarian Algorithm [3] is an optimization algorithm to deal with optimal assignment problem. The concept shown by the Bipartite graphs in Figure 2.1 and Figure 2.2 describes how Hungarian Algorithm matches eigenmodes. The optimal assignment problem involves finding a one-to-one mapping between two sets of nodes, where each edge connecting a pair of nodes which is associated with a specific cost. The goal is to determine a matching that minimizes the total cost across all selected edges. For example, there are 4 workers that need to be assigned to a task with specific cost, assigning every task with a worker by the minimum cost is the goal. 

Original	After Hungarian Algorithm
Figure 2.1: Original Bipartite graph	Figure 2.2: Bipartite graph after applying Hungarian Algorithm with minimum cost = 19

In this project, the Hungarian Algorithm is applied on every cost Matrix $//<![CDATA[ \begin{array}{l}D\end{array} //]]>$to match eigenmodes by searching the minimum gap between each eigenvalue pairs. 

2.1.2.  Result

2.1.3. Matching of Eigenvalues

2.1.4.  Revised Cost Matrix

It is observed in numerical experiment that setting the threshold and the value of the weight is difficult. Additionally, high resolution is required for finding intersection of dispersion curves in high frequency span, which consumes high computational resources. Besides, the weight $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ could be problem dependent, hence it is complicated and time-consuming to tune with. 

Therefore, an adapted method is proposed in equation(2.4) by setting weight $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ in equation(2.2) to 1, and introduce a new parameter, the penalty value $//<![CDATA[ \begin{array}{l}p\end{array} //]]>$, which depends on the orthogonality in equation(2.3) to overcome the difficulties mentioned above. Since the approximate orthogonality provides a criteria to distinguish the correct eigenmode, the penalty value $//<![CDATA[ \begin{array}{l}p\end{array} //]]>$could be added to the equation when the approximate orthogonality value is low so as to express the impossibility of the correctness of the match. 

This adapted approach offers two main advantages:
1. Reduced computational resource usage: there is no need for high resolution for the determination of the intersection, which saves considerable computational resources.

2. Easier hyper-parameter tuning: since tuning $//<![CDATA[ \begin{array}{l}p\end{array} //]]>$is a linear approach (Additive), whereas tuning $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$is a nonlinear approach (Multiplicative), which is more sensitivity and harder to predict and control.

2.1.5.  Visualization 

Figure 2.1 display the dispersion curves that needs for eigenmode tracing. The table T1 is the cost matrix computed at phase k, size of 6 x 6 indicates that we are interested in sorting the eigenmodes for the lowest 6 modes. The row number indicates the eigenmode at the current phase (k = i), and the column number represent the eigenmode at the next phase (k = i+1) . In each row, the minimum value has been colored with red, and the coordinate of this value indicates the correct match of the eigenmode from current phase (k = i) to the next phase (k = i+1). 

Figure 3.1: Original Dispersion curves

Table T1: Cost Matrix constructed by 6 bands

For example, T1(2,3) is the minimum element in the second row, and this shows that the second eigenmode should be connected to the third eigenmode. T1(3,2) is also the minimum in the third row. This result shows that, at the phase k = i, the second and third dispersion curve are crossing each other, and through this eigenmode tracing method, the second eigenmode and the third eigenmode is switched (Figure 2.2). The reason is that the values in the cost matrix shows that switching the eigenmodes at current phase(k = i) could ensure the difference within the eigenmodes to be minimum.

By visualization the dispersion curves after eigenmode tracing, the smoothness of the dispersion curves could be the proof of successful eigenmode tracing because of the fact that the dispersion curves should be continuous without kinks in classical Physics.     

Figure 3.2: Switching of eigenmodes

Finally, result like Figure 2.3 could be obtained.

2.2. Polarization

2.2.1. Method Principle

Polarization analysis is a fundamental technique for understanding the directional characteristics of wave propagation in structured materials. By examining the distribution of energy among different displacement components, this method provides essential insights into wave modes and their interactions. It plays a crucial role in identifying dispersion characteristics and analyzing wave behavior in various periodic and non-periodic media. This section discusses the principles of polarization computation, its application in dispersion curve analysis, and its significance in wave dynamics research.

Example of Polarization Method
	Figure 2.2.1: Shan Jiang et al 2018 J. Phys. D: Appl. Phys. 51 045601

Why Can the Polarization Method Distinguish Energy Bands?

The polarization method effectively distinguishes energy bands because different types of elastic waves—such as longitudinal and transverse waves—exhibit distinct vibration directions during propagation. By computing the kinetic energy distribution along different spatial directions, this approach allows for precise identification of wave characteristics and modes.

According to the study by Jiang, Hu et al. (2018), the polarization method is an effective approach for analyzing wave behavior by distinguishing different polarization characteristics. This method is particularly useful in identifying wave modes and understanding their interactions with structural periodicity or material heterogeneity, which influence band formation and wave propagation properties.

As elastic waves propagate, particle motion can be either parallel to the wave propagation direction (longitudinal waves) or perpendicular to it (transverse waves). The presence of complex structural configurations or material anisotropy affects these polarization modes differently, resulting in selective enhancement or suppression of specific wave components within certain frequency ranges.

Energy Distribution Analysis

By calculating the kinetic energy ratio (polarization ratio) along different spatial directions, wave modes can be effectively identified:

    • If kinetic energy is predominantly concentrated along the propagation direction, the wave is classified as a longitudinal wave(P wave).

P Wave: https://gpg.geosci.xyz/content/seismic/wave_basics.html

    • If kinetic energy is primarily distributed in a perpendicular direction, the wave is identified as a transverse wave(S wave).S Wave: https://gpg.geosci.xyz/content/seismic/wave_basics.html

Band Formation and Mode Differentiation

In structured materials, periodicity and anisotropy influence wave propagation by altering the energy distribution across different polarization directions, leading to the formation of energy bands and gaps. The polarization method helps determine how different wave modes are affected, whether certain modes are attenuated, allowed, or modified due to structural or material characteristics.

In summary, the polarization method serves as a powerful tool for distinguishing energy bands by analyzing wave vibration directions and energy distributions. By integrating this approach with material and structural characteristics, a comprehensive understanding of band formation and propagation properties can be achieved. This method provides valuable guidance for designing wave control strategies in various applications involving wave propagation and energy band engineering.

2.2.2. Computation Procedure

The polarization approach in eigenmode tracing is employed to identify and differentiate wave branches based on their displacement characteristics. This method systematically quantifies the directional energy distribution of each eigenmode, allowing for a precise classification of wave modes. The computation process consists of the following steps:

1. Displacement Component Calculation

    • For each eigenmode, the squared displacement components along the $//<![CDATA[ \begin{array}{l}x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}y\end{array} //]]>$and $//<![CDATA[ \begin{array}{l}z\end{array} //]]>$axes are computed.

    • This calculation provides insight into the proportion of energy stored in each directional mode, facilitating the characterization of wave polarization properties.

2. Normalization and Integration

    • The squared displacement components are integrated over the entire volume of the system to account for the global behavior of the mode.

    • A normalization factor is applied to ensure that polarization values range between 0 and 1, enabling direct comparisons across different eigenmodes.

    • The kinetic energy ratio (KER) is defined as:

    where $//<![CDATA[ \begin{array}{l}u_i\end{array} //]]>$represents the displacement component in the respective direction, and $//<![CDATA[ \begin{array}{l}V\end{array} //]]>$denotes the volume of the unit cell.

3. Polarization-Based Visualization

    • The normalized polarization values are mapped onto the dispersion curves using a color-coded representation.

    • This visualization technique provides an intuitive classification of wave branches based on their dominant polarization components:

    • Longitudinal waves → Predominantly polarized along the propagation direction $//<![CDATA[ \begin{array}{l}e_x\end{array} //]]>$.

    • Transverse waves → Strong polarization in the perpendicular directions ( $//<![CDATA[ \begin{array}{l}e_y\end{array} //]]>$or $//<![CDATA[ \begin{array}{l}e_z\end{array} //]]>$).

    • Mixed waves → A combination of multiple polarization components.

2.2.3. Results of Our Simlified Model

2.2.3.1. Key Observations

To simplify the analysis problem, we continue to use the previous 2D model, the two images depict the kinetic energy ratio in the X-direction as a function of frequency (vertical axis) and phase (horizontal axis) for 20 energy bands at different phase resolutions:

    • Left plot: Computed with 50 phase points

    • Right plot: Computed with 500 phase points

The color scale represents the kinetic energy ratio, where:

    • Red (~1.0): Strong X-direction polarization (dominantly longitudinal waves)

    • Blue (~0.0): Weak X-direction polarization (dominantly transverse or mixed modes)

    • Green/Yellow: Intermediate polarization

Initially, we considered selecting a total of 50 phases. However, after plotting the results in MATLAB, we found that the graph lacked accuracy. To improve precision, we increased the number of phases tenfold to 500 for comparison.

Comparison by Increasing the Number of Phases (Fixed 20 Bands)
Figure 2.2.3: Fixed 20 bands with 50 phases	Figure 2.2.4: Fixed 20 bands with 500 phases

Our Observations:

1. Uniform Color Distribution in Low-Frequency Bands

    • At lower frequencies, the color distribution remains relatively smooth and consistent across phases, regardless of the number of phases used (50 or 500).

    • This suggests that low-frequency modes are well-resolved even with a coarse phase resolution and exhibit stable polarization characteristics across the phase space.

    • These modes are typically governed by large-scale structural deformations, which are less sensitive to phase resolution effects.

2. Color Irregularities in High-Frequency Bands

    • In the higher frequency range, irregularities in the color distribution are observed in the 50-phase case, with abrupt changes or discontinuities.

    • These irregularities reduce significantly in the 500-phase case, where the bands become much more continuous and well-defined.

    • The observed color variations indicate that higher frequency modes exhibit more complex phase-dependent behavior, requiring finer phase resolution for proper characterization.

3. The major downside of increasing phase resolution is the significant rise in computational cost in MATLAB

    • The increase from 50 to 500 phases leads to a 10× increase in data points, resulting in longer computation time due to the need for solving additional eigenvalue problems at each phase step, and also larger memory consumption, especially when storing high-resolution energy band data.

2.2.3.2. Possible Explanation and Solution

    • In wave propagation studies, particularly in band structure analysis, the accuracy of mode identification depends on the sampling density in the phase space.

    • Higher frequency waves tend to have more oscillatory and rapidly varying behavior, meaning that insufficient sampling (e.g., only 50 phase points) can introduce numerical artifacts or “gaps” in the representation of polarization behavior.

    • By increasing the number of phase points to 500, the phase space is sampled more densely, leading to: 1.Smoother and more continuous kinetic energy distributions. 2.Better distinction between energy bands 3.Reduction of artificial discontinuities or noise.

    • Increasing the number of phases is an effective solution to improve accuracy, as demonstrated in the results, but a balance must be struck between computational efficiency and resolution.

    • If computational resources are limited, an adaptive phase resolution approach can be used (e.g., finer phase steps for higher frequencies and coarser steps for lower frequencies). Alternatively, interpolation techniques may help reconstruct missing data points in high-frequency regions without fully increasing phase density.

2.2.4.  Analysis of Mode Shapes

In this study, we selected low-frequency (Band 2), medium-frequency (Band 9), and high-frequency (Band 18)modes at a fixed phase step and used MATLAB to visualize their mode shapes, as shown in the figure.

1. Low-Frequency Mode (Band 2)

    • The deformation is minimal, and the overall structure maintains a relatively regular grid pattern, indicating that this mode is primarily characterized by global motion, with little localized deformation.

    • The displacement field is evenly distributed, suggesting that this mode corresponds to a long wavelength, where the entire region moves in sync. The energy distribution remains uniform across a large scale, rather than being concentrated in specific localized areas.

    • Low-frequency modes exhibit global motion with simple mode shapes and minimal local deformation. In this frequency range, energy propagates primarily over a large scale, with less influence from localized regions.

2. Medium-Frequency Mode (Band 9)

    • The deformation becomes more pronounced, and the overall shape is no longer uniform, indicating that local deformations are becoming more significant at this frequency range.

    • As the wavelength shortens, certain regions experience increased deformation amplitude, suggesting that energy begins to concentrate in specific areas rather than being uniformly distributed. The vibration mode exhibits a more complex shape, with non-uniform deformation amplitudes across different regions.

    • Medium-frequency modes introduce localized deformation and increasing complexity. The emergence of these modes indicates that local regions begin to influence the overall response, with shorter wavelengths and increased energy concentration.

3. High-Frequency Mode (Band 18)

    • The mode shape becomes increasingly complex, with certain regions undergoing significantly larger deformations, particularly near the edges, where high strain concentrations appear.

    • With a further reduction in wavelength, the localization effect of the deformation becomes more prominent, indicating that this mode is primarily governed by higher-order vibration patterns. The irregularity in shape increases, showing that the response becomes highly dependent on local geometric properties and boundary effects.

    • High-frequency modes display strong localization effects, complex mode shapes, and greater sensitivity to geometric variations. At this frequency range, different regions exhibit highly non-uniform responses, and localized deformation amplitudes may significantly exceed the overall average.

Why higher frequency modes are more likely to exhibit uneven color bands?

This analysis demonstrates clear differences in mode shapes across frequency bands. The transition from global motion at low frequencies to localized vibrations at high frequencies highlights the frequency-dependent nature of wave propagation. As frequency increases, mode shapes become more complex, localized deformation intensifies, wavelength shortens, and energy distribution becomes increasingly non-uniform. This phenomenon is of great significance for applications in vibration control, energy manipulation, and structural optimization.

1. Increased Mode Complexity Leads to Uneven Kinetic Energy Distribution

    • Higher frequency modes exhibit more complex deformation patterns, as seen in Band Nr.18, where localized deformations are more pronounced.

    • These localized regions concentrate kinetic energy, while other areas have significantly lower energy, leading to abrupt changes in kinetic energy ratio (KER) and causing the uneven color bands.

2. Higher Sensitivity of Kinetic Energy Ratio to Phase Variations

    • In high-order modes, even small changes in phase can drastically redistribute kinetic energy.

    • Unlike low-order modes, where energy transitions are smooth, high-frequency modes experience rapid fluctuations in KER, making color bands appear more irregular.

3. Stronger Mode Coupling Leads to Localized Energy Concentration

    • In higher frequency ranges, multiple local modes may couple together (e.g., a combination of local and global vibrations).

    • This coupling causes kinetic energy to fluctuate significantly across different phase values, where some regions may have very high energy at certain phases and much lower energy at others, leading to uneven color bands.

4. Numerical Errors and Computational Precision

    • High-frequency modes involve smaller-scale structural deformations, and insufficient mesh refinement may fail to capture the kinetic energy distribution accurately.

    • These numerical errors can cause discontinuities in kinetic energy ratio calculations, making the color bands appear more irregular in high-frequency regions.

2.2.5. Possible Optimization Strategies

Enhancing the effectiveness of polarization analysis requires targeted optimization strategies to improve visualization clarity, computational efficiency, and practical applicability. The following approaches address key challenges in polarization computation and analysis.

1. Adaptive Normalization for Multi-Band Systems

    • In higher frequency regimes, wave modes often exhibit localized deformations, resulting in uneven polarization color distributions that obscure mode differentiation.

    • Implementing an adaptive normalization scheme based on local energy distribution ensures consistent contrast in visualization, facilitating clearer identification of dominant polarization characteristics.

2. Computational Efficiency Enhancements

    • Dynamically adjusting the resolution of polarization maps based on mode complexity optimizes performance by allocating computational resources more efficiently.

    • Eigenmode filtering can eliminate weakly excited modes, reducing computational overhead while preserving the accuracy of dominant wave mode characterization.

    • Parallel processing techniques and numerical solvers tailored for high-resolution eigenmode analysis can further improve computational speed without sacrificing precision.

3. Integration with Band Gap Analysis

    • Establishing a direct correlation between polarization patterns and band gap locations enables more effective wave manipulation and control.

    • By tuning material properties and structural parameters, directional energy confinement can be optimized to enhance wave filtering and attenuation in specific frequency ranges.

    • The insights gained from polarization analysis can inform the design of waveguiding structures and energy localization mechanisms, leading to more efficient band gap engineering strategies.

3. Conclusion

In this study, we deal with a matching problem in Wave Finite Element Method (WFEM) analysis. The mismatching of eigenmodes due to the independent sorting of eigenvalues without considering the corresponding eigenvectors. To overcome this problem, matching of eigenvalues and mode polarization methods are implemented.

3.1. Matching of Eigenvalues

By constructing cost matrices that quantify the differences between eigenvalues and eigenvectors across phase steps and applying the Hungarian Algorithm, we demonstrated an effective strategy for tracing and matching eigenmodes. An adapted formulation, incorporating a penalty parameter based on approximate orthogonality, further improved the method’s robustness and simplified the tuning process. This approach ensures that even when dispersion curves cross or overlap, the correct mode connections are maintained, resulting in smooth and physically consistent dispersion curves. After the implementation, dispersion curves are sorted in a correct manner and it can be processed by computer for further use.

3.2. Polarization

The polarization method exploits the directional distribution of kinetic energy to distinguish between wave modes. By analyzing the kinetic energy ratios along different spatial directions, this method accurately classifies longitudinal, transverse, and mixed modes. The resulting color-coded visualizations of dispersion curves not only enhance our understanding of mode evolution but also confirm the effectiveness of the matching approach, especially in higher frequency regimes where mode complexity increases. Despite the fact that mode polarization method gives extensive information about the directional distribution of movements, further steps are required for the computer to track the mode.

4. Reference

[1] M. M. Alghamdi, F. Bertrand, D. Boffi, F. Bonizzoni, A. Halim, and G. Priyadarshi, "On the matching of eigensolutions to parametric partial differential equations," arXiv:2207.06145 [math.NA], Jul. 2022. [Online]. Available: https://doi.org/10.48550/arXiv.2207.06145.

[2] H. W. Kuhn, "The Hungarian method for the assignment problem," Naval Research Logistics Quarterly, vol. 2, no. 1-2, pp. 83–97, Mar. 1955. DOI: 10.1002/nav.3800020109.

[3] Y. Lu and A. Srivastava, "Level repulsion and band sorting in phononic crystals," Journal of the Mechanics and Physics of Solids, vol. 111, pp. 100–112, Feb. 2018. DOI: 10.1016/j.jmps.2017.10.021.

[4] S. Jiang, H. Hu, and V. Laude, “Low-frequency band gap in cross-like holey phononic crystal strip,” J. Phys. D: Appl. Phys., vol. 51, no. 4, p. 045601, 2018. DOI: 10.1088/1361-6463/aa9d35.

[5] “Wave picture source,” Geophysical Python Cookbook, [Online]. Available: https://gpg.geosci.xyz/content/seismic/wave_basics.html. Accessed: Feb. 13, 2025.

[6] V. Cool, E. Deckers, L. Van Belle, and C. Claeys, “A guide to numerical dispersion curve calculations: Explanation, interpretation and basic Matlab code,” Mech. Syst. Signal Process., vol. 215, p. 111393, Apr. 2024, DOI: 10.1016/j.ymssp.2024.111393.

[7] A. O. Krushynska, A. Amendola, F. Bosia, C. Daraio, N. M. Pugno, and F. Fraternali, “Accordion-like metamaterials with tunable ultra-wide low-frequency band gaps,” New J. Phys., vol. 20, no. 7, p. 073051, Jul. 2018, DOI: 10.1088/1367-2630/aad354.