## 2 Numerical Methods  
Elastic wave propagation through use of a phononic crystal was studied using a two-dimensional finite element method implemented in Solid Mechanics COMSOL modules [36]. The governing equation is given by Equation (1):  
\[\rho\frac{\partial^{2}u}{\partial t^{2}}=f_{\
u}+\
abla\cdot\sigma. \tag{1}\]  
where \(\rho\) and \(u\) are the density and displacement of the elastic materials, the volumetric force \(f_{\
u}\) is the force per deformed volume, and \(\sigma\) is the stress tensor. In this work, the proposed phononic crystals consist of high-stiffness inclusions placed in a low-stiffness matrix. The low-stiffness matrix in this study is epoxy, whereas the high-stiffness inclusions are composed of tungsten carbide [34]. The material properties are summarized in Table 1.  
The dispersion relations describing the wave propagation behavior in a periodic phononic crystal was obtained by solving the eigenvalue problem with the periodic bound ary condition. The Floquet-Bloch periodic boundary conditions were applied to the four edges of the unit cell to ensure a periodically repeating nature, as indicated by Equation (2).  
\[u(\overset{\rightarrow}{r})=u_{k}(r)e^{i\overset{\rightarrow}{k}\cdot\overset{ \rightarrow}{r}}. \tag{2}\]  
where \(k\) is the periodicity vector with a length \(\mid k\mid=2\pi/\Lambda\), in which \(\Lambda\) is the period in that direction and \(r\) is the position vector. Modal analysis was performed to obtain the eigenfrequency corresponding to a given Bloch wave vector \(k=(k_{x},\,k_{y})\). The study was performed on a range of wave vectors covering the irreducible Brillouin zone (IBZ), which spans from \(\Gamma\) to X, X to M, and M to \(\Gamma\), as shown in Figure 1a. The PnC is infinite and periodic in the \(x,\,y\) directions with a lattice constant \(a_{o}=1\) m. For the elastic wave transmission simulation, the proposed epoxy/tungsten carbide composite unit cells were centered in the simulation domain with a 1 m distance away from the absorption boundary. An incident force of 1 N was applied in front of the composite structure in the propagating direction, and edges of the simulation domain were applied with the perfectly matched layers [37], which absorb all outgoing wave energy. The transmission loss as a function of frequency is given by  
\[TL(dB)=10 log(P_{in}/P_{out}). \tag{3}\]  
where \(P\) is the power derived from mechanical energy flux, and \(P_{in}\) and \(P_{out}\) are the power of the incident and transmitted elastic waves.  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Material/Property** & **Elastic Modulus E (GPa)** & **Density \(\rho\) (kg/m\({}^{3}\))** & **Poisson’s Ratio \(\
u\)** \\ \hline Epoxy & 4.344 & 1180 & 0.3679 \\ Tungsten Carbide & 387.6 & 13,800 & 0.3459 \\ \hline \hline \end{tabular}
\end{table}
Table 1: The material properties of the matrix and inclusion used in the study.  
Figure 1: (**a**) The irreducible Brillouin zone (IBZ) in a 2D phononic crystal. (**b**) The phononic crystal consisting of one-square inclusion in the soft host, where the filling fraction is controlled by varying the parameters \(a_{i}\) and \(a_{o}\). (**c**) The dispersion curves of a phononic crystal with a square inclusion, where \(a_{o}=1\) m, \(a_{i}=0.25\) m, FF = 0.5, and Young’s modulus ratio (YMR) = 90. (**d**) For filling fraction (FF) = 0.5, variations in bandwidth with filling fraction at YMR = 9, 90, and 300, respectively.