I think the issue with the mode shapes is the calculation of the matrix that is being used for the eigenanalysis. The OP gives mat = np.dot(K_6, np.linalg.inv(M_6)) where it should be mat = np.dot(np.linalg.inv(M_6), K_6). However, I would recommend using the scipy.linalg.eig function which supports specifying the A and B matrices to solve the generalized eigenvalue problem (e.g. scipy.linalg.eig(K, M)). The following is the complete code to calculate and plot the mode shapes for the BeamDyn based NREL 5MW blade as generated by the BeamDyn driver, followed by the plot it produces:
import numpy as np
import yaml
import matplotlib.pylab as plt
data = yaml.load(open("bd_driver.BD.sum.yaml"), Loader=yaml.Loader)
x = np.array(data["Init_Nodes_E1"])[1:, 2]
K_BD = np.array(data["K_IEC"])
M_BD = np.array(data["M_IEC"])
# eigval, eigvec = scipy.linalg.eig(K_BD[6:,6:], M_BD[6:,6:])
eigval, eigvec = np.linalg.eig(np.dot(np.linalg.inv(M_BD[6:, 6:]), K_BD[6:, 6:]))
idx = np.argsort(eigval)
freq = np.sqrt(eigval[idx]) / 2 / np.pi
shapes = eigvec[:, idx]
fig, ax = plt.subplots(6, figsize=(5, 9), sharex=True)
for i in range(3):
ax[0].plot(x, shapes[0::6, i], label=f"Mode {i+1}")
ax[1].plot(x, shapes[1::6, i], label=f"Mode {i+1}")
ax[2].plot(x, shapes[2::6, i], label=f"Mode {i+1}")
ax[3].plot(x, shapes[3::6, i], label=f"Mode {i+1}")
ax[4].plot(x, shapes[4::6, i], label=f"Mode {i+1}")
ax[5].plot(x, shapes[5::6, i], label=f"Mode {i+1}")
ax[0].set_ylabel("UX")
ax[1].set_ylabel("UY")
ax[2].set_ylabel("UZ")
ax[3].set_ylabel("RX")
ax[4].set_ylabel("RY")
ax[5].set_ylabel("RZ")
ax[-1].set_xlabel('blade loc (m)')
plt.tight_layout()
fig.align_ylabels(ax)
Please give this a try and let us know if you get reasonable mode shapes. Thanks!
