Honeycomb Structures and Resolution in FEM

[ShashankanB]

TL;DR Summary

I struggle to implement the perturbation of a wave in a Honeycomb structure : i can't have the results which can be find on internet.

Hi,
I am a engineering student and I am in apprenticeship with a mechanics lab : I am studying about Honeycomb propagation and periodic structures. I hadn't have any courses in FEM or in Bloch's theory, I've done all by myself : but there is the problem, i didn't practice or use viable courses. I tired to implement the propagation in 2D of a wave in these honeycombs, with the following code in Python :

[CODE lang="python" title="Honeycomb wave propagation (function names are in french)"]from Matrices import mass_matrix,stiffness_matrix
from Maillage import Elements_decoupe_unique
from Maillage import Affiche_elements_decoupe as AED
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eig

rayon = 0.01
element = Elements_decoupe_unique(2, 1, rayon, opt = False, offset = True)
periode = np.array([element[5]] + [element[7]] + [element[3]] + [element[6]] + [element[4]])
nb = len(periode) #Nombre d'elements
nn = nb + 1 #Nombre de noeuds

#Initialisation
h = rayon*np.array([0.5,0.5,0.5,0.5,1])
K = np.zeros((2*len(h)+2,2*len(h)+2)) #Matrice de Raideur
F = np.zeros(2*len(h)+2) #Matrice des forces
M = np.zeros((2*len(h)+2,2*len(h)+2)) #Matrice des masses

L = np.sqrt((periode[0,1] - periode[1,1])[0] ** 2 + (periode[0,1] - periode[1,1])[1] ** 2)
#Valeurs tabulées
E = 71e9
I_zz = 1e-6
r = 2700
s = 1e-4

for i in range(len(h)):
K[i*2 : 2*i + 4,2*i : 2*i + 4] += stiffness_matrix(E,I_zz,h)
M[i*2 : 2*i + 4,2*i : 2*i + 4] += mass_matrix(r,s,h)

def Porpagation2D():
"""
Effectue la méthode de Bloch-Floquet sur un maillage hexagonal avec une perturbation bidimensionnelle

Argument :

Retoune:
- plot : affichage de la courbe de dispersion

"""
a = h[-1]
t0,L0 = 0.005, rayon
w0 = np.pi**2*t0/L0**2*np.sqrt(E/(12*r))

# Reciprocal lattice basis vectors
b1 = (2 * np.pi / a) * np.array([1, -1 / np.sqrt(3)])
b2 = (2 * np.pi / a) * np.array([0, 2 / np.sqrt(3)])

# Define high-symmetry points in reduced coordinates
P_G = np.array([0, 0])
P_K = np.array([2/3, 1/3])
P_M = np.array([1/2, 0])

# Chemin de la zone de Brillouin
num_points = 100
path = np.vstack([
np.linspace(P_G, P_K, num_points),
np.linspace(P_K, P_M, num_points),
np.linspace(P_M, P_G, num_points)
])


kx_ky = path @ np.array([b1, b2])
omega_values = np.zeros((len(kx_ky), 2 * nn), dtype=complex)
for i,(k1, k2) in enumerate(kx_ky):
# Construct the transformation matrix T
T = np.eye(2 * nn, dtype = complex)

T[0,2] = np.exp(1j*(k1*L))
T[1,3] = np.exp(1j*(k1*L))
T[2,0] = np.exp(-1j*(k1*L))
T[3,1] = np.exp(-1j*(k1*L))

T[4,6] = np.exp(1j*(k2*L*np.sqrt(3)/2 + k1*L/2))
T[5,7] = np.exp(1j*(k2*L*np.sqrt(3)/2 + k1*L/2))
T[6,4] = np.exp(-1j*(k2*L*np.sqrt(3)/2 + k1*L/2))
T[7,5] = np.exp(-1j*(k2*L*np.sqrt(3)/2 + k1*L/2))

K_transformed = np.dot(np.conj(T).T, np.dot(K, T))
M_transformed = np.dot(np.conj(T).T, np.dot(M, T))

eigenvalues, _ = eig(K_transformed, M_transformed)
omega_values[i, :] = np.sqrt(eigenvalues)/w0

distances = np.sqrt(np.sum(np.diff(kx_ky, axis=0) ** 2, axis=1))
k_path = np.insert(np.cumsum(distances), 0, 0)
for i in range(omega_values.shape[1]):
plt.plot(k_path, np.real(omega_values[:, i]),'.', lw = 0, ms = 5, c = 'black')

plt.grid()
plt.xlabel(r"Vecteur d'onde réduit (k) $m^{-1}$")
plt.xticks([0, k_path[num_points-1], k_path[2*num_points-1], k_path[-1]], ["Γ", "K", "M", "Γ"])
plt.ylabel(r"Pulsation temporelle normalisée ($\Omega$) rad/s ")
plt.title("Courbe de dispersion")
plt.show()

Porpagation2D()[/CODE]

But the result is absurd (left, mine and right from a research paper):

VS

If you have any suggestion or help, I am thanking you in advance,

P.S. : Please don't mind any English errors.