Implementation of WFEM in 2D and 3D

Authors: Weiyi Kong, Paula Löcherer

Supervisors: Sebastian Schopper, Lukas Kleine-Wächter

Introduction

Many structures such as honeycomb can be considered to be periodic. One property of periodic structures is that they can show stop band behavior. The identification of stop bands, i.e. zones with no free wave propagation in structures is a major concern in the study of structural behavior, since the noise and vibration levels can be reduced by these frequency zones with attenuated structural response. 

Wave finite element method combines the conventional finite element method with the periodic boundary conditions. Instead of investigating the whole structure, computational cost can be largely reduced by limiting calculation to a unit cell due to the periodicity. In the inverse approach, the stop bands are identified by solving eigenvalue problems for a set of propagation constants, which is related to the wave numbers.

The inverse approach of one-dimension and two-dimension has already been implemented in the Matlab package bm-fem developed by the chair of structural mechanics at TUM. The goal of this project is to extend the scripts to the three-dimensional case and validate the codes by a benchmark case. 

Theory

In this project, we consider undamped three-dimensional structures. In general, a dynamic mechanical system without damping can be described by the equation

where $//<![CDATA[ \begin{array}{l}\bm{M}\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\bm{K}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{F}\end{array} //]]>$ are the mass matrix, stiffness matrix and load vector repectively, while $//<![CDATA[ \begin{array}{l}\bm{q}\end{array} //]]>$ denotes the nodal displacement vector. For a continuum, the discrete matrices and vectors can be obtained from the finite element method. Assume that the discretized structure has $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$ degree of freedom, then the dimension of $//<![CDATA[ \begin{array}{l}\bm{M}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{K}\end{array} //]]>$ is $//<![CDATA[ \begin{array}{l}n\times n\end{array} //]]>$, while $//<![CDATA[ \begin{array}{l}\bm{F}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{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 vector of degrees of freedom can be expressed by an exponential approach.

Then the system equation can be transformed into an eigenvalue problem.

The system has $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$ eigenvalues $//<![CDATA[ \begin{array}{l}\omega^2_r\end{array} //]]>$ and corresponding eigenfrequencies $//<![CDATA[ \begin{array}{l}\bm{\Phi}_r\end{array} //]]>$, which is called eigenmode or mode shape.

In the following sections, the analysis of the infinite periodic structure is reduced to a unit cell by Bloch's theorem, which generalizes Floquet's theorem for the second and third dimension. The three-dimensional cube is modelled by the finite element approach. By imposing compatible displacements and force balance at the boundaries, a reduced system is deduced.  To solve the free vibration problem of the reduced system, the inverse approach is applied, where the propagation vector is purely imaginary. The eigenvalues represent frequencies of free wave propagation, which means that wave can propagate without attenuation.

Three-dimensional wave finite element method

In this section, the WFEM is extended to three-dimensional structures. We assume that the unit structure is a cube with length $//<![CDATA[ \begin{array}{l}L_x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}L_y\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}L_z\end{array} //]]>$ in $//<![CDATA[ \begin{array}{l}x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}y\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}z\end{array} //]]>$ direction respectively. We utilize conventional finite element methods to calculate the stiffness matrix and mass matrix of the discretized unit-cell and then apply periodic boundary conditions to obtain the reduced matrices.

The boundary degrees of freedom can be categorized into three groups corresponding to the corners, edges and faces.

where

Similarly, the vector of forces can be written as:

where

Bloch's theorem and wave propagation

Bloch's theorem relates the nodal displacements and forces of a structure to the neighboring structure, which enables us to reduce the system of the periodic problem to a unit cell. Assume that we have a reference point with position vector $//<![CDATA[ \begin{array}{l}\bm{p}_{ref}\end{array} //]]>$ with respect to the local coordinate system of the regarding unit cell. Then the position vector of any point inside the analyzed unit-cell with respect to a reference unit-cell can be written as:

where $//<![CDATA[ \begin{array}{l}\bm{d}_1\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\bm{d}_2\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{d}_3\end{array} //]]>$ are vectors in $//<![CDATA[ \begin{array}{l}x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}y\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}z\end{array} //]]>$ direction with length $//<![CDATA[ \begin{array}{l}L_x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}L_y\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}L_z\end{array} //]]>$ respectively. The natural numbers $//<![CDATA[ \begin{array}{l}n_1\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}n_2\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}n_3\end{array} //]]>$ represents the number of unit-cells moved in $//<![CDATA[ \begin{array}{l}x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}y\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}z\end{array} //]]>$ direction with respect to the reference cell.

Bloch's theorem states that the relative change of the state variables from any location in the unit cell to the reference point can be expressed by an exponential term. When a wave travels from one unit-cell to the other at frequency $//<![CDATA[ \begin{array}{l}\omega\end{array} //]]>$, the wave amplitude $//<![CDATA[ \begin{array}{l}V\end{array} //]]>$ at position $//<![CDATA[ \begin{array}{l}\bm{p}\end{array} //]]>$ can be formulated. It can be further reduced as following, because the relative change of $//<![CDATA[ \begin{array}{l}V(\bm{p})\end{array} //]]>$ is independent of the position of the unit-cell within the periodic structure: 

where $//<![CDATA[ \begin{array}{l}V(\bm{p}_{ref})\end{array} //]]>$ is a certain plane wave disturbance perpendicular to the directions of propagation and $//<![CDATA[ \begin{array}{l}\bm{\kappa}\end{array} //]]>$ is the complex wave vector, where the real part describes the change in phase while the imaginary part describes the change in amplitude in $//<![CDATA[ \begin{array}{l}x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}y\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}z\end{array} //]]>$ direction. Alternatively, the wave vectors can be expressed by propagation constants which are defined as:

Thus, the phase change of the wave is represented by the imagary part of propagation constants, while the decay of amplitude is described by the real part.

Three-dimensional inverse approach

In the inverse approach, we assume free wave motion without external loads and attenuation. The propagation constants are hence strictly imaginary. For the three-dimensional unit cell considered, the degrees of freedom $//<![CDATA[ \begin{array}{l}\bm{q}_{K1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{q}_{K7}\end{array} //]]>$ can be related to each other as

Similarly, the nodal loads $//<![CDATA[ \begin{array}{l}\bm{f}_{K1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{f}_{K7}\end{array} //]]>$ satisfies

Analogously, we can obtain the relationships of other boundary nodal displacements and forces.

From the displacement compatibility and force balance boundary conditions, a reduction matrix $//<![CDATA[ \begin{array}{l}\bm{R}\end{array} //]]>$ can be derived as

The degree of freedom of the system can be reduced by the reduction matrix. The reduced displacement vector is defined as

The complete vector and the reduced degree of freedom vector are related by

Since the internal forces are equal to zero, the equilibrium of force can be expressed as

where $//<![CDATA[ \begin{array}{l}\bm{R}^{*T}\end{array} //]]>$is the complex conjugate of the matrix $//<![CDATA[ \begin{array}{l}\bm{R}\end{array} //]]>$.

The equation of motion can thus be formulated as

Therefore, the equation of motion is reduced to an eigenvalue problem:

where $//<![CDATA[ \begin{array}{l}\bm{M}^{red}=\bm{R}^{*T}\bm{M}\bm{R}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{K}^{red}=\bm{R}^{*T}\bm{K}\bm{R}\end{array} //]]>$are the reduced mass and stiffness matrices. The frequencies of free wave propagation can hence be calculated by solving the eigenvalue problem for $//<![CDATA[ \begin{array}{l}(-\omega^2\bm{M}^{red}+\bm{K}^{red})\end{array} //]]>$. In the case of strictly imaginary propagation constants, the matrix $//<![CDATA[ \begin{array}{l}(-\omega^2\bm{M}^{red}+\bm{K}^{red})\end{array} //]]>$ is Hermitian so the solutions for $//<![CDATA[ \begin{array}{l}\omega\end{array} //]]>$ are always real.

Irreducible Brillouin contour and band structure 

To obtain the stop bands, we must solve the eigenvalue problem for every possible combination of propagation constants $//<![CDATA[ \begin{array}{l}k_x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}k_y\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}k_z\end{array} //]]>$. The corresponding solution $//<![CDATA[ \begin{array}{l}\omega\end{array} //]]>$, which are the frequencies for free propagation make up the dispersion surface in the three-dimensional wave domain.

In mathematics and solid-state physics, the first Brillouin zone is a uniquely defined primitive cell in reciprocal space. The solutions of periodic structures can be completely characterized by their behavior in a single Brillouin zone. Therefore, we do not need to solve the problem for the entire wave domain. By the symmetry of the unit cell, the Brillouin zone can be further reduced to an irreducible Brillouin zone, which is the smallest unit to capture all information to obtain the band gaps.

To further reduce the computational cost, the problem is only investigated along the contour of the irreducible Brillouin zone, as the maxima and minima of each dispersion surface can be found on the boundary. For the three-dimensional unit cell, the corner points of the irreducible Brillouin zone are $//<![CDATA[ \begin{array}{l}O=(0,0,0)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}A=(0,\pi,0)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}B=(0,\pi,\pi)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}C=(\pi,\pi,\pi)\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}O=(0,0,0)\end{array} //]]>$. We only consider the wave numbers $//<![CDATA[ \begin{array}{l}(\kappa_x, \kappa_y, \kappa_z)\end{array} //]]>$ respectively the propagation constants $//<![CDATA[ \begin{array}{l}(k_x, k_y, k_z)\end{array} //]]>$ varying along the contour (O-A-B-C).

Algorithm and Implementation

In this project, we utilize the finite element tool bm-mfem, which is a Matlab package developed by the chair of structural mechanics at TUM. The package is capable to read a mesh file with various element types (e.g. tetrahedron element with 4 nodes, tetrahedron element with 10 nodes and hexahedron element with 8 nodes) and assemble the stiffness and mass matrices of the discrete system.

Firstly, we generate a .mdpa mesh file in GID, which is a universal pre and post processor for numerical simulations in science and engineering. In wave finite element, it suffices to generate the mesh of a single unit-cell. The generated mesh file contains information about nodes, elements and boundary conditions. The second step is to load the model, add specified degrees of freedom to each node and set material properties. For a three-dimensional solid structure, we have three degrees of freedom, i.e., displacements in x, y and z direction at each node.

After the model is built, a solver for the three-dimensional inverse approach is created. The wave numbers along the contour of Brillouin zone and number of bands to compute are specified and passed to the solver. The corresponding eigenfrequencies, eigenvectors and propagation constants are returned. Finally, a post-processing function is used to append the repeating structures to the existing model of unit-cell. The displacements are calculated according to Bloch's theorem. In addition, a GID model is exported for visualization.

Initialization

To apply the three-dimensional inverse approach, the first class function implemented in the solver is initialization. Generally, the information of boundary node indices are not included in the .mdpa generated by GID. Therefore, it is necessary to identify the boundary nodes and degrees of freedom, which is a prerequisite for imposing periodic boundary conditions.

The boundary node indices are found by their coordinates. The algorithm is illustrated by the search of K1 nodes, which is implemented by the function findVertex1Nodes.

function [nodeIdsK1] = findVertex1Nodes(obj)

First, the size of K1 node index vector is counted by the following codes:

for i=1:l_xCoords
    if nodeXcoords(i) == minX && nodeYcoords(i) == minY && nodeZcoords(i) == minZ
         n = n+1;
    end
end

Here minX, minY and minZ represent the smallest coordinate in x, y and z directions respectively. Then nodes with x coordinate equal to minX, y coordinate equal to minY and z coordinate equal to minZ are identified as K1 nodes and the returned K1 index vector is stored in the variable nodeIdsK1.

for i=1:l_xCoords
    if nodeXcoords(i) == minX && nodeYcoords(i) == minY && nodeZcoords(i) == minZ
        m = m+1;
        nodeIdsK1(m) = nodeIds(i);
    end
end

After K1 nodes have been found, we continue to search for other vertex nodes. Note that the number of K nodes are all identical due to symmetry, the tedious procedure for the identification of vector size n can be avoided by using nodeIdsK1 as input variable. For example, nodeIdsK2 can be detected by

nodeIdsK2 = obj.findMaxXNodes(nodeIdsK1);

The function findMaxXNodes finds nodes whose y and z coordinates are the same as the input variable, while x coordinate is equal to maxX, i.e., the largest x coordinate.

function [nodeIdsMaxX] = findMaxXNodes(obj, nodeIdsN)
...
for i = 1:length(nodeIdsN)
    for j = 1:length(nodeXcoords)
        if nodeYcoords(j) == NNodeArray(i).getY && ...
            nodeZcoords(j) == NNodeArray(i).getZ && ...
            nodeXcoords(j) == maxX && ...
            any(nodeIdsMaxX == j) == false
            nodeIdsMaxX(i) = nodeIds(j);
        end
    end
end
...

Analogously, other K vertices can be found by findMaxXNodes, findMaxYNodes and findMaxZNodes functions. The identification of edge node indices and face node indices are implemented in the same manner. The obtained edge or face index vectors also contain vertices or edge nodes; thus, the deletion of these repeated nodes is followed. For instance, the code below deletes E9 nodes from F1. 

[~, posF1, ~] = intersect(nodeIdsF1, nodeIdsE9);
nodeIdsF1(posF1) = [];

The node index vectors are then stored in class properties KNodeIds, ENodeIds and FNodeIds, which are cells with size $//<![CDATA[ \begin{array}{l}(1\times 8)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}(1\times 12)\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}(1\times 6)\end{array} //]]>$ respectively.

The next step is to search boundary degrees of freedom according the node indices. 

function [KdofIds, EdofIds, FdofIds, IdofIds] = getAllDofIds(obj)

The fixed boundary nodes are then deleted from the degree of freedom index vectors. The returned arrays are stored in class variables Kdofs, Edofs, Fdofs and Idofs.

Finally, the stiffness matrix and mass matrix of the discretized model is assembled by functions implemented in bm-mfem and assgined to class properties massMatrix and stiffnessMatrix.

obj.massMatrix = obj.assembler.assembleGlobalMassMatrix(obj.femModel);
obj.stiffnessMatrix = obj.assembler.assembleGlobalStiffnessMatrix(obj.femModel);

3D inverse solving strategy

Solve function is one method function inside BlochInverse3DSolvingStrategy class. The inputs of this function are the phase-vectors in x, y and z-direction, depicting the variation of propagation constants $//<![CDATA[ \begin{array}{l}(k_x, k_y, k_z)\end{array} //]]>$ and the number of bands, which is the number of eigenfrequencies to be calculated.

function [counter, frequencies, preEigVector, prePropConst] = solve(obj, phasesX, phasesY, phasesZ,  numberOfBands)

It returns the eigenfrequencies and a related counter to numerize them, as well as the affiliated eigenvectors. 

In the following sections, some of the important secondary class method functions are briefly explained.

Implementation of the reduction matrix

For easier assembling of the reduction matrix, the input phase vectors are first transformed into exponential terms which define the relationship between degrees of freedom at boundary nodes. Three vectors $//<![CDATA[ \begin{array}{l}\bm{muX}=e^{i*\bm{phasesX}}\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\bm{muY}=e^{i*\bm{phasesY}}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\bm{muZ}=e^{i*\bm{phasesZ}}\end{array} //]]>$ are calculated and used to initialize the reduction matrix.

For saving storage space, especially for fine meshes, the R-matrix is already initialized as a sparse matrix with the dimensions [length($//<![CDATA[ \begin{array}{l}\bm q\end{array} //]]>$) x length($//<![CDATA[ \begin{array}{l}\bm q^{red}\end{array} //]]>$)].

j = [1:K1,1:K1,1:K1,1:K1,1:K1,1:K1,1:K1,1:K1,K1+1:E1,...
     E1+1:E4,K1+1:E1,E1+1:E4,K1+1:E1,E1+1:E4,K1+1:E1,E1+1:E4,E4+1:E9,E4+1:E9,E4+1:E9,E4+1:E9,...
     E9+1:F1,F1+1:F4,E9+1:F1,F1+1:F4,F4+1:F5,F4+1:F5,...
     F5+1:I];
i = 1:length(j);
v = [ones(1,numK1),muX*ones(1,numK1),muX*muZ*ones(1,numK1),muZ*ones(1,numK1),muY*ones(1,numK1),muX*muY*ones(1,numK1),muX*muY*muZ*ones(1,numK1),muY*muZ*ones(1,numK1),...               
	 ones(1,numE1),muY*ones(1,numE4),muX*ones(1,numE1),ones(1,numE4),muZ*ones(1,numE1),muY*muZ*ones(1,numE4),muX*muZ*ones(1,numE1),muZ*ones(1,numE4),ones(1,numE9),muX*ones(1,numE9),muY*ones(1,numE9),muX*muY*ones(1,numE9),...
     ones(1,numF1),muX*ones(1,numF4),muY*ones(1,numF1),ones(1,numF4),ones(1,numF5),muZ*ones(1,numF5),...
     ones(1,numI)];

R = sparse(i,j,v);

where K1, E1, E4, E9, F1, F4, F5 and I are the endindex of the corresponding nodes inside the reduced displacement vector, while numK1, numE1, numE4, numE9, numF1, numF4, numF5 and numI are the number of nodes accordingly.

Computation of the reduced matrices and vectors

For the reduction of the matrices, already implemented functions are used to obtain beforehand the sorted stiffness matrices for the object. They are both of size [length($//<![CDATA[ \begin{array}{l}\bm q\end{array} //]]>$) x length($//<![CDATA[ \begin{array}{l}\bm q\end{array} //]]>$)] and can be reduced to [length($//<![CDATA[ \begin{array}{l}\bm q^{red}\end{array} //]]>$) x length($//<![CDATA[ \begin{array}{l}\bm q^{red}\end{array} //]]>$)] by premultiplying the complex-conjugate of R and postmultiplying R ($//<![CDATA[ \begin{array}{l}\bm{M}^{red}=\bm{R}^{*T}\bm{M}\bm{R}\end{array} //]]>$,$//<![CDATA[ \begin{array}{l}\bm{K}^{red}=\bm{R}^{*T}\bm{K}\bm{R}\end{array} //]]>$) as described before. The sorted matrices are subsequently also saved as sparse matrices.

Solving of the eigenvalue problem

The reduced eigenvalue problem (22) itself is solved inside a loop for all phases. After applying the beforehand explained reduction to the matrices, the eigenvalue problem can be solved by using the in MATLAB implemented function eigs() to calculate the first numberOfBands eigenfrequencies.

omega = eigs(Kred,Mred,numberOfBands)

Then the eigenvalues are transformed to frequencies by using the relation $//<![CDATA[ \begin{array}{l}f = \sqrt{ \frac{| \omega |}{2\pi}}\end{array} //]]>$ and stored with the corresponding propagation constants, before starting the next iteration of the loop.

The same procedure is utilized inside another loop to calculate the eigenvectors, if they are also of interest. The reduced displacements are regained as explained in the following paragraph.

Adding reduced displacements

Afterwards, the readding of the displacements reduced by the Bloch theorem could be implemented by the matrix multiplication $//<![CDATA[ \begin{array}{l}\bm q = \bm R * \bm q^{red}\end{array} //]]>$. However, a direct approach is implemented alternatively in the solver in order to optimize the calculation time.

The reduced displacements are directly regained by the underlying relations of the R-matrix, e.g. $//<![CDATA[ \begin{array}{l}\bm{q}_{K7}=e^{i(\left | k_x \right | +\left | k_y \right |+\left | k_z \right |)}\bm{q}_{K1}\end{array} //]]>$ (16).

...
modeK7 = mode(1:K1, :) * muX * muY * muZ;
...

After calculating all the displacements and appending them to the reduced vector, the complete vector has to be sorted by another function. In contrast for the matrix multiplication they would have already been in the right order.

It is noteworthy that the calculation of the reduction matrices could be conducted in a similar manner to optimize the calculation time. However, the relations for the mass and stiffness matrices are much more complex than for the displacements. Therefore, an easier and clearer approach is implemented instead to avoid programming mistakes.

Computation and visualization of the stop bands

For validation purposes, a small Matlab script was created, where varying propagation constants along the Brillouin Contour are handed to the solve function of an already initialized model. The phases vectors ranges with a given resolution (numPhases) from the respective coordinates of the corner points in the order $//<![CDATA[ \begin{array}{l}O,A,B,C,O\end{array} //]]>$. They are passed to the solver together with the number of bands (numBands) that need to be calculated.

% tetrahedral Brilloin zone:
phasesX = [linspace(zero, zero, numPhases), linspace(zero, zero, numPhases), linspace(zero, pi, numPhases), linspace(pi, zero, numPhases)];
phasesY = [linspace(zero, pi,   numPhases), linspace(pi,   pi,   numPhases), linspace(pi,   pi, numPhases), linspace(pi, zero, numPhases)];
phasesZ = [linspace(zero, zero, numPhases), linspace(zero, pi,   numPhases), linspace(pi,   pi, numPhases), linspace(pi, zero, numPhases)];

% solve:
[counter, frequencies] = solver.solve(phasesX, phasesY, phasesZ, numBands);

The returned frequencies and affiliated counters are then visualized and compared with reference stop bands from literature or other programs. For a clearer understanding, vertical lines are displayed at the corner points of the Brillouin Zone.

Visualization of the mode shapes

The model of the system consists of $//<![CDATA[ \begin{array}{l}(n+1\times m+1 \times p+1)\end{array} //]]>$ cells is generated by a auxiliary function appendUnitCells3D. 

function [femModel] = appendUnitCells3D(femModel, n, m, p, eigenvalue, eigenvector)

This function appends n cells in x direction, m cells in y direction and p cells in z direction to the original model in the order of z to y to x direction. The input variable eigenvalue is the propagation constant matrix of one band with size $//<![CDATA[ \begin{array}{l}[3\times N_{phase}]\end{array} //]]>$, where $//<![CDATA[ \begin{array}{l}N_{phase}\end{array} //]]>$ denotes the size of phasesX, phasesY and phasesZ, while the input variable eigenvector is the nodal displacement matrix of a band with size $//<![CDATA[ \begin{array}{l}[N_{dof} \times N_{phase}]\end{array} //]]>$, with $//<![CDATA[ \begin{array}{l}N_{dof}\end{array} //]]>$ representing the degrees of freedom. 

The nodal displacement vector corresponding to the phases of all appended cells are also calculated. The variable used to assemble these eigenvectors of is a cell of size $//<![CDATA[ \begin{array}{l}[(n+1)\times (m+1) \times (p+1)]\end{array} //]]>$, each entry storing the displacement matrix of a unit cell. The calculation is based on Bloch's theorem.

...
tmp{(m+1)*(p+1)*i + 1} = eigenvalue(1, k)*tmp{(m+1)*(p+1)*(i-1) + 1};
...
tmp{(m+1)*(p+1)*i + (p+1)*j + 1} = eigenvalue(2, k)*tmp{(m+1)*(p+1)*i + (p+1)*(j-1) + 1};
...
tmp{(m+1)*(p+1)*i + (p+1)*j + q + 1} = eigenvalue(3, k)*tmp{(m+1)*(p+1)*i + (p+1)*j + q};
...

The displacements of the current cell are related to the neighboring ones by $//<![CDATA[ \begin{array}{l}k_x\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}k_y\end{array} //]]>$ or $//<![CDATA[ \begin{array}{l}k_z\end{array} //]]>$.

These modes are then stored in the model as displacements at different "time" steps, the number of time steps is hence equal to $//<![CDATA[ \begin{array}{l}N_{phase}\end{array} //]]>$. The appended femModel can be exported to GID post-processing files by the following commands.

out = GiDOutput('outputfile');
out.setPreference('writeModelParts',false);
out.writeMesh(femModel);
out.setPreference('nodalResults',["DISPLACEMENT"]);
out.writeResults(femModel);

Validation and results

Benchmark problem for three-dimensional unit-cell

To validate our codes, the vibrational behavior of a three-dimensional unit-cell inside an infinite undamped cube is used. The material characteristics which are of a type of steel and dimensions of the unit-cell in this example are summarized as follows:

The unit-cell can be further reduced using the irreducible Brillouin contour along the corner points $//<![CDATA[ \begin{array}{l}O=(0,0,0)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}A=(0,\pi,0)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}B=(0,\pi,\pi)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}C=(\pi,\pi,\pi)\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}O=(0,0,0)\end{array} //]]>$ as defined above.

The rule of thumb for appropriate mesh size

The element size in element–based acoustic computations should be related to the wavelength. Most commonly the rule of thumb "six to ten nodes per investigated wavelength" is used as a guideline for mesh size.

The wavelength is computed as $//<![CDATA[ \begin{array}{l}\lambda=\frac{c}{f}\end{array} //]]>$. Thus the sound velocity of shear and pressure waves for the given material properties are calculated as follows.

Therefore, the two wavelengths can be approximated by using an estimate of $//<![CDATA[ \begin{array}{l}f_{est}=2750Hz\end{array} //]]>$ for the maximum frequency in the first band.

Hence a mesh size of about $//<![CDATA[ \begin{array}{l}0.195 m\end{array} //]]>$ to $//<![CDATA[ \begin{array}{l}0.117 m\end{array} //]]>$ is fine enough according to the rule of thumb. With the unit-cells dimensions of $//<![CDATA[ \begin{array}{l}1m\end{array} //]]>$ in each direction, a mesh size of roughly 5-10 elements per edge seems to be an appropriate number.

Results for the band plots

The first step for validating the code is to compare the frequencies of the first band along the irreducible Brillouin contour with the results from [Ramirez] as shown below:

For a first validation, two unit-cells with a tetrahedron mesh were inspected. Regarding the rule of the first one with five elements per edge the second with 10. The coarser mesh created in GID is depicted below, along with the first band for both unit cells.

The plot for the coarser mesh is depicted on the left-hand side and the plot for the finer one is depicted on the right-hand side. Both are a good approximation of the given band, only the highest eigenfrequency $//<![CDATA[ \begin{array}{l}f\approx2750Hz\end{array} //]]>$ is slightly overestimated by the coarse mesh, which is in harmony with the rule of thumb, where a slightly smaller mesh size was calculated (min 0.195m).

Furthermore, different mesh sizes and also different element types are compared to a plot for the first six bands gained by the computation of an unit-cell with identical properties with the program Comsol Multiphysics provided by our supervisors.

Firstly, the results for the already discussed two unit-cells with a tetrahedron mesh as well a finer mesh with 25 elements per edge are depicted below. (From left to right the mesh gets finer).

It seems that a coarser mesh yields a better approximation of the results from Comsol in higher frequencies. This contradicts the expected results regarding the rule of thumb, that higher frequencies are depicted better by a finer mesh. To save computation time, for the finest discretization a smaller number of phases along the Brillouin Contour was used. Therefore, the shape of the contour cannot be fully compared to the other plots, because the smaller number of phases leads to edgier curves in the plot.

Then the influence of different types of mesh elements is depicted. Instead of tetrahedrons, hexahedrons were used for the mesh. Again, the first one with 5 elements per edge the second with 10. The coarser mesh created in GID is depicted below, as well as the first 7 bands for both unit cells.

Again, the coarser mesh gives the better approximation. Furthermore, it seems better than the tetrahedron mesh of the same size, but the results for a finer mesh clearly show mistakes between $//<![CDATA[ \begin{array}{l}A=(0,\pi,0)\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}B=(0,\pi,\pi)\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}C=(\pi,\pi,\pi)\end{array} //]]>$.

Lastly, a second-order tetrahedron element was used, again with a mesh size of 5 and 10 elements per edge. The results are depicted below, the finer mesh on the right.

The results of the 2nd order mesh with 5 elements per side seems to be between the results from first-order tetrahedron elements with mesh size of 10 and 25. The number of nodes of the second-order element lies between these two, therefore this result is as expected. On the other hand, the result for the finer mesh has more in common with the 10-element hexahedron plot.

Results for the Visualization of Mode Shapes

One cell in each direction is appended to the model and exported to .msh and .res files, which can be loaded and processed in GID. In this section, we only visualize the real part of the eigenmodes since the imaginary part representing the phase change is not of concern. Our main interest is in the mode shapes of the first two bands at turning points of the Brillouin zone contour, namely point O, A, B and C. Since the validation is conducted qualitatively, a coarse tetrahedron mesh with 5 elements on each edge is employed in this case.

The mode shapes at point O, A, B and C of the first band is scaled and depicted as below.

Note that the mode at point O represents the rigid body translation in x direction.

The following pictures illustrate the nodal displacements of the second band. The mode shape at point O is the rigid body translation in y direction.

Although the largest amplitude of the displacements are close to each other in the first two bands, the different movements can be clearly observed from the figures.

Conclusion and outlook

To summarize the deviations of the band plots, it seems there is a problem for larger models, concerning the number of DOFs. In [Waki] is stated that the eigenvalue problem can be ill-conditioned because it contains a wide range of eigenvalues. This problem becomes worse when a finer discretization is used, because then the system matrices become larger and comprise more eigenvalues.

The results of only the first band are in a fairly good agreement with the curves as presented in [Ramirez]. For coarser discretization and a limited number of eigenvalues the inverse approach can therefore be used as implemented for further research.

References

Ramirez, Jose Daniel Perez: Investigation of Periodic Structures for the Design of Acoustic Metamaterials Using Unit-Cell Modeling,Technical University of Munich, master thesis, 2017

Waki, Y.: Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides” (2009)

Claus C. Claeys, Karel Vergote, Paul Sas, Wim Desmet: On the potential of tuned resonators to obtain low-frequency vibrational stop bands in periodic panels, Journal of Sound and Vibration, Volume 332, Issue 6, 2013, Pages 1418-1436, ISSN 0022-460X, https://doi.org/10.1016/j.jsv.2012.09.047

Matlab package bm-fem developed by the Chair of Structural Mechanics at Technical University of Munich

Supervisors Sebastian Schopper and Lukas Kleine-Wächter work with bm-fem and Comsol Multiphysiks

Speed of Sound for three dimensional Solids from: https://en.wikipedia.org/wiki/Speed_of_sound#Three-dimensional_solids (14.02.2022)