WFEM 2D

Authors: Maveryck Andrés Garzon Espejo , Unbekannter Benutzer (ge32qad) , Roman Maier 

Supervisors:  Sebastian Schopper, Tom Hicks

Introduction

The Wave Finite Element Method (WFEM) and Bloch's Theorem applied to the Boundary Conditions are powerful tools for analyzing vibrations in repetitive 2D cell structures. The WFEM is a numerical technique used to solve partial differential equations that describe wave propagation. It involves dividing the domain into smaller, manageable elements and solving for the unknown solution within each element. On the other hand, Bloch's Theorem provides a framework for applying periodic boundary conditions to wave propagation problems. When applied to the Boundary Conditions, Bloch's Theorem requires that the solution to the wave equation be periodic at the boundary, making it well suited for modeling periodic structures such as those found in 2D cells.

By combining the Wave FEM with Bloch's Boundary Conditions, it is possible to analyze the vibrations in 2D cells and solve an eigenvalue problem. The eigenvalue problem involves finding the frequencies and modes of vibration for a given system. By solving this problem, it is possible to gain a deeper understanding of the behavior and properties of the 2D cells, including its resonant frequencies, energy distribution, and stability. Lately, these kind of analysis have been applied to design and study Acoustic Metamaterials as these materials typically make use of resonators that can interact with wavelengths much larger than the unit-cell length scale [1].

So far, this WFEM has only been implemented for 1D- and 2D-Periodic structures with rectangular unit cells in bm-fem. In the scope of this project, the WFEM (inverse approach) should be implemented for 2D-periodic structures with unit cells of various shapes. 

Main Objectives

• Extend the WFEM in bm-fem code for 2D-periodic structures to rhombic and hexagonal lattices.
• Validate the implementation with analytic solutions and metamaterial designs from literature.

Secondary Objectives

• Implement a function to process automatically meshes to be used to solve a 2D WFEM problem.
• Generate a function that returns the type of lattice imported into the solver.Generate a function to draw the Brillouin Zone associated to the lattice vectors of the cell.

Theory

Mathematical Model

In this project, undamped two-dimensional structures are considered (yet approaches of Mode Synthesis, as pointed out in [1] can include damping formulations) . In general, a dynamical mechanical system without damping can be described by the following equation:

where  $//<![CDATA[ \begin{array}{l}\mathbf{K}\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\mathbf{M}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mathbf{f}\end{array} //]]>$ are the mass matrix, stiffness matrix and load vector repectively, whereas $//<![CDATA[ \begin{array}{l}\mathbf{q}\end{array} //]]>$ denotes the nodal displacement vector. The matrices and load vectors are discrete representations of the continuous problem, discretized according to any simple Finite Element formulation (e.g. 3-Node Plane-Strain). By assuming that the discretized structure has  $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$ degrees of freedom (DOF), then the dimension of $//<![CDATA[ \begin{array}{l}\mathbf{K}\end{array} //]]>$ and  $//<![CDATA[ \begin{array}{l}\mathbf{M}\end{array} //]]>$ is $//<![CDATA[ \begin{array}{l}n \times n\end{array} //]]>$, while $//<![CDATA[ \begin{array}{l}\mathbf{f}\end{array} //]]>$and$//<![CDATA[ \begin{array}{l}\mathbf{q}\end{array} //]]>$ are of dimension $//<![CDATA[ \begin{array}{l}n \times 1\end{array} //]]>$.

For free vibration, the system oscillates with natural frequencies without external loads, hence the load vector vanishes.

The latter equation (2) represents the free vibration of a mechanical system. The Fourier Ansatz, whereby the displacement can be represented by a time-dependent harmonic function, is used:

Then the system of equations can be transformed into an eigenvalue problem.

The system has $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$ eigenvalues  with eigenfrequencies  $//<![CDATA[ \begin{array}{l}\omega^{2}\end{array} //]]>$ and corresponding eigenmodes or eigenshapes $//<![CDATA[ \begin{array}{l}\mathbf{\phi}\end{array} //]]>$. 

As the cells studied are distributed periodically (i.e. are repetitive in space) and supposing there are infinitely many repetitions, then properties of the cell must behave periodically as well. Bloch's Theorem, as a generalization to Floquet's Theorem, relates periodicity of any property in space [2]. By supposing that the waveforms are periodically distributed  and by applying a Fourier Transform in space to(3) the following relation is obtained:

In (5) the term $//<![CDATA[ \begin{array}{l}\mathbf{\tilde{\phi}}\end{array} //]]>$ represents a periodic waveform function and, $//<![CDATA[ \begin{array}{l}\mathbf{x}=\left[x, y \right]^{T}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mathbf{k} = \left[k_{1}, k_{2} \right]^{T}\end{array} //]]>$ represent the position vector and wave vector respectively. As $//<![CDATA[ \begin{array}{l}\mathbf{\tilde{\phi}}\end{array} //]]>$ is assumed to be periodic, then the following periodicity constraints must be satisfied:

In (6), the vector $//<![CDATA[ \begin{array}{l}\mathbf{r} = a_{1} \mathbf{e_{1}} + a_2 \mathbf{e_{2}}\end{array} //]]>$ is called the lattice vector, which is the vector that represents the primary translation vectors of the lattice. The constants $//<![CDATA[ \begin{array}{l}a_1\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}a_2\end{array} //]]>$ are integers that are related to the periodicity as this constants define translation among cells. The covariant basis vectors $//<![CDATA[ \begin{array}{l}\mathbf{e_{1}}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mathbf{e_{2}}\end{array} //]]>$ define a basis to identify positions inside the cell. Often, these vectors are referred as direct lattice basis vectors [3] as these vectors span the edges of the unit cell and their respective magnitude denotes the lengths of the parallel pairs.

By combining Eqs. (5) and (6), Bloch's Theorem can be expressed by a shift which relates boundary conditions of an unit cell [1]. 

The components of the wave vector are complex values: $//<![CDATA[ \begin{array}{l}k_{1} = \varepsilon_{1} + i \beta_{1}\end{array} //]]>$,   $//<![CDATA[ \begin{array}{l}k_{2} = \varepsilon_{2} + i \beta_{2}\end{array} //]]>$ in which the complex part  $//<![CDATA[ \begin{array}{l}\varepsilon_{1,2}\end{array} //]]>$ represents actually a periodic phase shift, whereas the real part $//<![CDATA[ \begin{array}{l}i \beta_{1,2}\end{array} //]]>$ represents the change in amplitude. By assuming that the real component of the wave vector is zero $//<![CDATA[ \begin{array}{l}i\beta_1=i\beta_2=0\end{array} //]]>$    (i.e. the waves do not change in amplitude by a periodic shift in position), then the frequency dependence of the waves. By this assumption, the wave vector can be represented as:

wherein $//<![CDATA[ \begin{array}{l}\eta \in \mathbb{N}^{2}\end{array} //]]>$. Then, by this definition, (7) can be represented as follows (this kind of representation is found in most references such as [3], [8] and [9]):

In this context the coefficients  $//<![CDATA[ \begin{array}{l}n_{1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}n_2\end{array} //]]>$ are arbitrary and real valued.

Reciprocal Space, Brillouin Zone and Lattices

In solid-state physics, the first Brillouin zone is a uniquely defined primitive cell in reciprocal space. Let  $//<![CDATA[ \begin{array}{l}\mathbf{e^{1}}\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\mathbf{e^{2}}\end{array} //]]>$ be two linear independent vectors which can be used as basis of the space in $//<![CDATA[ \begin{array}{l}\mathbb{C}^{2}\end{array} //]]>$ Mathematically, the definition of the reciprocal space depends on the definition of the lattice vector space. This is shown in (10):

In the latter, $//<![CDATA[ \begin{array}{l}\delta^{i}_{j}\end{array} //]]>$ represents the Kronecker delta. The wave vector  $//<![CDATA[ \begin{array}{l}\mathbf{k}\end{array} //]]>$ can be represented by a linear combination of the reciprocal basis. According to Gonella et. al. [9], the main advantage of the representation of the wave vector using the reciprocal space is that it describes the periodicity of the frequency/wavenumber relation, whereas the direct lattice vectors represent the spatial periodicity of the lattice. To show this statement, suppose that the representation of the wave vector is:

Then, by using (11) and (10), then the coefficients $//<![CDATA[ \begin{array}{l}k_{i}\end{array} //]]>$ can be calculated as:

Now suppose regard the repeatability of (8), then by supposing there exists some  $//<![CDATA[ \begin{array}{l}\mathbf{\eta '}\end{array} //]]>$ such that:

Then, a solution by using this definition would be expressed as:

By this, it can be shown that the coefficients are: $//<![CDATA[ \begin{array}{l}k_{i}' =2\pi \mathbf{\eta'} \cdot e_{i} = k_{i} + 2 \pi m_{i}\end{array} //]]>$. This shows, in fact the advantage shown by [9], but pose the difficulty of defining uniqueness due to the periodicity. Yet, by searching possible combinations of the wavevectors wherein just a period of the waves is investigated. This is known as the Brillouin Zone. The period corresponds to a region in the reciprocal lattice whose area equals the area of the reciprocal lattice’s unit cell [9]. By different kind of symmetries of the unit cell, the Brillouin zone can be further reduced to an irreducible Brillouin zone [2].

In Figure 1, the different kinds of 2D lattices are shown with their respective Brillouin Zones and Irreducible Brillouin Zone (IBZ) representations. The grey areas represent the (IBZ) for each of the lattices. Thereby stepping all over such areas with combinations of wavevectors, uniquely defined eigenmodes and eigenfrequencies can be found. Obviously, depending on the symmetry (Maurin et. al. [3] points relations to symmetry and point groups), the grey areas could be reduced further than in the traditional sketch as shown in Figure 1.

Figure 1. Representation of the different main types of 2D lattices and their respective sketches of the Irreducible Brillouin Zone (IBZ). From left to right: oblique, rectangular, rhombic, square and hexagonal lattices. Extracted from [3].  

One relevant aspect of these 2D lattices is that each of them is represented by a parallelogram [3]. This poses and important approach on how to identify the type of lattice given a set of nodes of a FEM mesh. 

Application of Bloch's Boundary Conditions into the Discrete Model

From Eqs. (3) and (7) a relationship can be stated between the displacements of the nodes on the boundary of the lattice.

Representation of Boundaries of an unit cell

Figure 2. Representation of Boundaries of typical 2D lattice. Extracted from [1].

For illustrative purpose, consider a discretized cell as shown in Figure 2. The nodes lying on the boundaries can be segmented into groups depending on the location. By using the shift stated in (7), the displacement on the left boundary $//<![CDATA[ \begin{array}{l}\mathbf{q}_{L}\end{array} //]]>$ is related to the displacement of the right boundary $//<![CDATA[ \begin{array}{l}\mathbf{q}_{R}\end{array} //]]>$ as follows:

Analogously, the displacement bottom $//<![CDATA[ \begin{array}{l}\mathbf{q}_{B}\end{array} //]]>$ and top $//<![CDATA[ \begin{array}{l}\mathbf{q}_{T}\end{array} //]]>$ nodes relationship can be expressed as:

In the same fashion, by taking into account the set of corner nodes $//<![CDATA[ \begin{array}{l}\left[ \mathbf{q}_{BL}, \mathbf{q}_{BR}, \mathbf{q}_{TL}, \mathbf{q}_{TR} \right]^{T}\end{array} //]]>$, the following relations can be obtained:

The relationships expressed in (15), (16) and (17) represent a node synthesis as some nodes of the discrete cell can be expressed in terms of others, thus reducing the DOFs to solve of the system. The synthesis can be expressed in matrix form as follows:

The latter equation shows that the matrix $//<![CDATA[ \begin{array}{l}\mathbf{R}\end{array} //]]>$ is composed of different identity $//<![CDATA[ \begin{array}{l}\mathbf{I}\end{array} //]]>$ submatrices, which are not equivalent in size. By this method of synthesis, the related pairs of edges (top-bottom) and (left-right) should have the same number of nodes. Furthermore, the nodes should be aligned such that the shift described in (7). For 2D cells, the corner nodes are one element sets. 

From (18), the eigenvalue problem described in (4) is modified as follows:

In the latter the operator $//<![CDATA[ \begin{array}{l}\mathbf{{}^H}\end{array} //]]>$ denotes the Hermitian Transpose or Complex Conjugate of a matrix and $//<![CDATA[ \begin{array}{l}\mathbf{\bar{\phi}}\end{array} //]]>$ denotes a synthesized eigenmode. On the other hand the pre- and postmultiplication of the mass and stiffness matrices by $//<![CDATA[ \begin{array}{l}\mathbf{R^{H}}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mathbf{R}\end{array} //]]>$, respectively, performs a reduction of the size of those matrices as well. 

The modified eigenvalue problem, shown in (19), is solved by identifying the frequencies of wave propagation by stepping around the Brillouin Zone of the lattice. In [1], this approach is called as  the "$//<![CDATA[ \begin{array}{l}\mathbf{w(k)}\end{array} //]]>$ solution". This approach is considered simple as combinations of different wavevectors $//<![CDATA[ \begin{array}{l}\mathbf{k}\end{array} //]]>$ are used to calculate respectively the propagation coefficients $//<![CDATA[ \begin{array}{l}\mu_{1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mu_{2}\end{array} //]]>$ to compute matrix $//<![CDATA[ \begin{array}{l}\mathbf{R}\end{array} //]]>$. After this calculation , the eigenvalue problem can be solved. 

Algorithms and methods

For this project, the existing "BlochInverseSolvingApproach" solver was modified to handle the processing of non-rectangular repetitive structures. First, the Usual Simulation Procedure is introduced.

Usual Simulation Procedure

The workflow, which is described in Figure 3, was not modified for this project. But you can see the process how the user gets to a solution with the bm_fem tool. After the geometry and the corresponding finite element mesh have been created in GiD, the model is imported into Matlab. Then the material parameters are defined and the solver can be started. After that the dispersion curves are plotted. If the user wants to get additional results, e.g. the eigenmodes, a postprocessing in GID is necessary.

Figure 3. Flowchart for usual simulation procedure.

The procedure is mainly simillar to the way a Finite Element Code works, wherein a mesh is required. In the process, this mesh has to be imported by using bm-fem designated functions such that the mesh is managed as an object within such framework. Afterwards, the material properties are set on the nodes and, then, the solver can be called. After running the solver, the solver have two combinations of outputs: one comprises a set of dispersion curves encapsulated within a matrix

where $//<![CDATA[ \begin{array}{l}s\end{array} //]]>$ represents the number of steps defined to loop over the Brillouin Zone and $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$ the number of modes ; the other combination comprises the same set of dispersion curves as well as a combination of dispersion coefficients 

calculated within the solver and a set of eigenmodes associated to certain branch of the dispersion curves  ,

which has this representation as a matrix of vectors. If the user were to choose the second combination of outputs, then with other functions coupled within the bm-fem framwork, there is the possibility to export these modes for visualization in GiD. A sample code on how to use the solver is the following:

io=MdpaInput('Plate20Shell4.mdpa'); % Some input file 
model = io.readModel(); %read the model

model.getAllNodes.addDof(["DISPLACEMENT_X","DISPLACEMENT_Y"]); % dofs for shell as example here


%Element properties: here aluminium
E = 7.31e10;
nu = 0.35;
rho = 2770;
t = 1e-7;

model.getAllElements.setPropertyValue('YOUNGS_MODULUS',E);
model.getAllElements.setPropertyValue('POISSON_RATIO',nu);
model.getAllElements.setPropertyValue('NUMBER_GAUSS_POINT',3);
model.getAllElements.setPropertyValue('DENSITY',rho);
model.getAllElements.setPropertyValue('THICKNESS',t);

% Create Bloch Solver
solver = BlochInverse2DSolvingStrategy(model);

% define phases and number of bands
numPhases = 50;             % 'resolution' of the dispersion curve

% Set a starting o reference zero (to avoid singularities in the eigenvalue problem)
zero = 1e-60;

% Create two vectors to define a contour within the Brillouin Zone (in this case Gamma - X)
phasesX = linspace(zero, zero, numPhases);
phasesY = linspace(zero, pi, numPhases);

% compute first 2 frequency bands/modes
numBands = 3;         

% solve --- Use the second set of outputs
[counter, frequencies, preEigVector, prePropConst] = solver.solve(phasesX, phasesY, numBands);

% Get the lattice vectors from the solver
[r1, r2] = solver.calculateLatticeVectors();

% Band of visualization
kk = 3;

% Visualization Steps
femModel = appendUnitCells2DInLattice(model,3,3,[r1;r2],prePropConst, preEigVector{kk});
filename = convertStringsToChars("examples/hexagonal_sample_mode_" + num2str(kk) +"_");
out = GiDOutput(filename);
out.setPreference('writeModelParts',false);
out.writeMesh(femModel);
out.setPreference('nodalResults',["DISPLACEMENT"]);
out.writeResults(femModel);

Solver Routine

For convenience, in the latter excerpt of code, the function "appendUnitCells2DInLattice" was implemented. This function enables the replication of the cells in the directions of the lattice vectors. This is an extension of the older function "appendUnitCells2D" as this function just considered the lattice to be replicated in a rectangular or regular manner.

On the other hand, the solver was part of the core of the project, as the extension of the previous code to handle all cases of 2D Periodic structures required extra steps. Some of the steps are encapsulated in the Initialization routine, which is shown in Figure 6. However, in Figure 4, some steps are added; namely two. The first one is that after the initialization, the solver can automatically draw the Brillouin zone of the identified lattice. This is helpful for the purpose of visualization, as the user can actually interact and check the extension and shape of the Brillouin Zone, and also check the computation of the lattice vectors and the reciprocal/contravariant basis is computed properly. 

Figure 4. Flowchart of the new Solver Implementation.

After that, the solver has to reformulate the way the wavevectors are given. The user needs to input the wave vectors as interpreted with the reciprocal basis (with $//<![CDATA[ \begin{array}{l}k_1\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}k_2\end{array} //]]>$ bounded by $//<![CDATA[ \begin{array}{l}[-\pi,\pi]\end{array} //]]>$) . Yet, due to the fact that the lattice vectors are used and saved with respect to the Euclidean basis  $//<![CDATA[ \begin{array}{l}[1,0]^{T}\end{array} //]]>$ and   $//<![CDATA[ \begin{array}{l}[0,1]^{T}\end{array} //]]>$, and the matrices use the same basis, then a transformation should be performed to handle rotations and scalings to accomodate that. Figure 5 shows an example of that, by considering that in the reciprocal space, the magnitudes of the reciprocal space of any reciprocal basis is the same. In this example a square lattice is rotated, yet in the reciprocal space the problem is equivalent, but in the Euclidean Space the definition is different. That step accomodates the transformation to enable that the computation of (9) is performed. 

Figure 5. Rotation of square lattice and redefinition of the vectors. Left: Unrotated square lattice. Right: Rotated square lattice.

Initialization Routine

Finally, the initialization routine has been the most crucial part of the project, as it had to be generalized for any 2D lattice and not just rectangular cells aligned with the $//<![CDATA[ \begin{array}{l}xy\end{array} //]]>$ axis. In this case the approach was to find a bounding parallelogram (as mentioned in the section "Reciprocal Space, Brillouin Zone and Lattices") of the FEM mesh. In Figure 6, the process of finding the convex hull is separated from the function which finds the bounding parallelogram in the code. Yet is it worth mentioning that the main requirement of the bounding-parallelogram-finding algorithm starts by searching the convex hull. In MATLAB the function convhull is called which performs such computation. This function takes as inputs the positions (in space) of all the nodes of the mesh to compute the convex hull and returns all the positions of the nodes which makes up the convex hull, i.e. the nodes located along the boundaries. By having this information, then the bounding parallelogram could be found. The main idea of the algorithm is based on the rotating calipers algorithm. The code of this algorithm was written by D'Errico [4] and was carefully saved as an external code in the framework of bm-fem (the definition and steps of the implementation by D'Errico [4] follows the ones described in this link). 

Figure 6. Flowchart of the new Initialization Routine.

The outputs of the bounding parallelogram finding function are the positions of the corners of the parallelogram. By having this information, the covariant basis vectors can be defined. As mentioned in [1], the definition is technically arbitrary as there are possible 8 combinations based on which of the corners of the parallelogram is chosen as an "origin" and how to order the lattice vectors as 1 or 2. To set the algorithm to work automatically, the origin corner was set as the possible leftmost one. In case the parallelogram is found to be a rectangle, then the origin corner is the lowest of all the leftmost corners. In case of the ordering, the definition on which one to select as  $//<![CDATA[ \begin{array}{l}\mathbf{e_{1}}\end{array} //]]>$ is by setting the other pointing corner as the next one in the Counter-Clockwise direction and  $//<![CDATA[ \begin{array}{l}\mathbf{e_{2}}\end{array} //]]>$ as the next one in Clockwise direction.

Detection of the Brillouin zone

With this in mind, a representation of the direct lattice basis, defined by a parallelogram, is shown in Figure 7.

Figure 7. Parallelogram and definition of the direct lattice vectors.

By having the definition of the lattice vectors, the lattice can be classified as shown in Table 1.

Boundary nodes

After the determination of the lattice type, the nodes, which intersect each of the sides of the parallelogram, are found and classified. The nodes of the mesh are classified in different arrays by checking which side it intersects (i.e. bottom, right, top, left). If the node of the cell corresponds to one of the corners of the parallelogram, it is assigned correspondingly as bottom-left, bottom-right, top-right or top-left. Also some lattices could not have "corner" nodes and just nodes intersecting the edges of the parallelogram. The nodes intersecting the edges are sorted with the idea that the associativity of the nodes, as shown in eqs (15) and (16), is fulfilled. After the sorting is performed, finally the quality checks are carried out. Firstly, it is checked that the number of nodes of pair of edges (top-bottom and left-right) are the same, Secondly it is checked the associativity of the nodes on the edges. The associativity is checked by setting a vector pointing from corresponding bottom to top (or left to right) node and ensure that the vector described is the same as one of the corresponding lattice vectors. The visual representation of the checks is shown in Figure 8. 

Figure 8. Visual representation of checks to ensure the nodes are sorted correctly. Left: Equal number of nodes. Right: Node associativity

Validation

For the validation, three different cases are considered, which are taken from two different publications 3 and 7. There are different assumptions regarding the kinematic relations. Both in the validation cases and in the following example, the kinematic model plain strain is used.

Material properties 

The cells are made of isotropic material aluminum with the following material properties:

Geometry

The unit cells possessed a constant thickness.

Analytical Solution 

For validation, analytical cases are inevitable, therefore, according to [7], a homogeneous  plate is considered. Whereby the following wavenumber can be calculated: 

$//<![CDATA[ \begin{array}{l}\mathrm{k_x}\end{array} //]]>$ and$//<![CDATA[ \begin{array}{l}\mathrm{k_y}\end{array} //]]>$ describes the wavenumber in x- and y-direction, $//<![CDATA[ \begin{array}{l}\mathrm{d}\end{array} //]]>$ is the dimension of the unit cell and $//<![CDATA[ \begin{array}{l}\mathrm{\theta}\end{array} //]]>$ discribes the angle the betweeen the lattice vectors $//<![CDATA[ \begin{array}{l}\mathrm{e_1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mathrm{e_2}\end{array} //]]>$.  The phase velocity describes the propagation speed of a wave with the same phase. The phase velocity of a transverse or shear wave in a homogeneous plate can be calculated as follows:

$//<![CDATA[ \begin{array}{l}\mathrm{\mu}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mathrm{\lambda}\end{array} //]]>$ are the Lamè Parameters. These can be calculated as follows: 

So that the eigenfrequency in 7 can be calculated as follows:

The dimensonless frequency is calculated in: 

Normalization

The normalization of the dispersions is done with the following equation: 

$//<![CDATA[ \begin{array}{l}\mathrm{r_1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\mathrm{r_2}\end{array} //]]>$ describes the lattice vectors. 

Results

Table 2 gives an overview of the cases used, the used finite element mesh is presented, as well as an animation for the calculation of the dispersion curves by following an specific path in the Brillouin zone, and some examples from the literature compared to the results obtained from the bm_fem tool. Three different unit cells are now considered in Table 2. The first one is a Square unit cell, the second a oblique unit cell and the thrid a sierpinski cell (named after the creator of the cell).  For the first two cases, the path in the brillouin zone is defined from $//<![CDATA[ \begin{array}{l}\mathrm{\Gamma}\end{array} //]]>$ to $//<![CDATA[ \begin{array}{l}\mathrm{X}\end{array} //]]>$ . For the Sierpinksi cell, the path as in the brillouin zone of 3 is used (red path in column and row fourth of Table 2). In order to validate the implementation of the code, the dispersion curves from the literature are compared with the results from the bm_fem tool. There is a very good agreement with the cases from the literature (see Table 2, columns 4, 5). For the first two cases, exists an analytical solution which are shown as solid lines and  the numerical solution is shown as a dotted line (see. 7  Page 18). If we compare the analytical and the numerical solution of  7 to the results of the square unit cell and the oblique unit cell of the bm_fem tool, we see that the dispersion curves are almost identical. 

Table 2: Validation cases from 3 and 7

	Mesh generated with GiD (blue) and mesh from Literature (grey)	Brillouin zone and the Dispersion Curve by the bm_fem tool	Dispersion Curve by Literature	Dispersion Curve by bm_fem tool	Sources
Square unit cell					7
Oblique unit cell					7
Sierpinski cell					3

The Sierpinski cell represents a special feature. There is a band gap in between, which means that for certain frequencies there is no wave propagation. This band gap is in 3 a thick gray line. In the result of the bm_fem tool the horizontal gap between the red curves is the so-called band gap. This band gap represents a range of values for which there is no solution to the eigenvalue problem, furthermore, the position and the width of the band gap is influenced by the distribution of the mass in the cell and the distribution of the stiffness. This is further explained in [10] where a cell is modeled as a lumped-mass configuration. From the validation cases it is concluded that the bm_fem tool provides correct results and different examples can be analyzed, which will be explained in the next section Examples. 

Examples

Other examples can be calculated. A special case is the rhombic cell. Furthermore, there are cases in the literature that serve as an example for the application of the code. A rhombic cell and Oblique cell with different holes is shown in Table 3. 

Table 3: rhombic cell and Oblique cell with different holes

You can find more examples in the appendix.

A wide variety of geometries can be solved with the bm_fem tool according to the natural frequencies and the eigenmodes in two dimensions. Thus, e.g. engineering problems in acoustics can be solved for small repeating structures. 

References

1. D. Krattiger and M. I. Hussein, ‘Generalized Bloch mode synthesis for accelerated calculation of elastic band structures’, Journal of Computational Physics, vol. 357, pp. 183–205, 3 2018.

2. W. Kong and P. Löcherer, Implementation of WFEM in 2D and 3D, 20-Feb-2022. [Online]. Available: Implementation of WFEM in 2D and 3D. [Accessed: 06-Feb-2023]. 

3. F. Maurin, C. Claeys, E. Deckers, and W. Desmet, ‘Probability that a band-gap extremum is located on the irreducible Brillouin-zone contour for the 17 different plane crystallographic lattices’, International Journal of Solids and Structures, vol. 135, pp. 26–36, 3 2018.

4. J. D'Errico, “A suite of minimal bounding objects,” MathWorks, 26-Jan-2012. [Online]. Available: https://de.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects. [Accessed: 07-Feb-2023]. 

5. X. Xu, X. Wang, and Y. Mei, ‘A design method for absorption of low-frequency noise using acoustic metamaterials’, 8 2019, vol. 569.

6. N. Guarín-Zapata and J. Gómez, 'Evaluation of the Spectral Finite Element Method With the Theory of Phononic Crystals'. Departamento de Ingeniería Civil-Universidad EAFIT, Medellín, Colombia, 2015.

7. C. Valencia, J. Gómez and N. Guarín-Zapata, 'A general purpose element-based approach to compute dispersion relations in periodic materials with existing finite element codes'. Universidad EAFIT-Departamento de Ingeniería Civil, Medellín, Colombia, 2019.

8. M. Ruzzene, F. Scarpa, and F. Soranna, ‘Wave beaming effects in two-dimensional cellular structures’, Smart Materials and Structures, vol. 12, pp. 363–372, 6 2003.

9. S. Gonella and M. Ruzzene, ‘Analysis of in-plane wave propagation in hexagonal and re-entrant lattices’, Journal of Sound and Vibration, vol. 312, pp. 125–139, 4 2008.
10. Y. Li, X. Wang and G. Yan 'Analytical dispersion curves and bandgap boundaries for quadrilateral lattices', European Journal of Mechanics / A Solids, vol 97, 2023
11. T. Hicks, ‘Characterization of Dispersive Wave Propagation in Periodic Media of Different Symmetry Groups’, Technical University of Munich, 11 2022.

Appendix

A square cell with rectangular and square holes is shown in Table A1. In addition to the dispersion curve, the first two eigenmodes are shown.

Table A1: Square cell with rectangular and square holes