From 6232ee1cdc0ff3714d55b75b91eb3cab26ba4a0e Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Wed, 11 Oct 2023 23:19:35 +0200 Subject: [PATCH] TST: new test for piston theory --- tests/test_quad4r_piston_theory.py | 218 +++++++++++++++++++++++++++++ 1 file changed, 218 insertions(+) create mode 100644 tests/test_quad4r_piston_theory.py diff --git a/tests/test_quad4r_piston_theory.py b/tests/test_quad4r_piston_theory.py new file mode 100644 index 0000000..f0ff846 --- /dev/null +++ b/tests/test_quad4r_piston_theory.py @@ -0,0 +1,218 @@ +import sys +sys.path.append('..') + +import numpy as np +from numpy import isclose +from scipy.sparse.linalg import eigsh, eigs +from scipy.sparse import coo_matrix + +from pyfe3d.shellprop_utils import isotropic_plate +from pyfe3d import Quad4R, Quad4RData, Quad4RProbe, INT, DOUBLE, DOF + + +def test_quad4r_piston_theory(plot=False, refinement=1): + data = Quad4RData() + probe = Quad4RProbe() + nx = refinement*21 + ny = refinement*21 + + a = 0.8 + b = 0.5 + + E = 70.e9 # Pa + nu = 0.3 + + rho = 7.8e3 # kg/m3 + h = 0.0035 # m + + xtmp = np.linspace(0, a, nx) + ytmp = np.linspace(0, b, ny) + + dx = xtmp[1] - xtmp[0] + dy = ytmp[1] - ytmp[0] + + xmesh, ymesh = np.meshgrid(xtmp, ytmp) + ncoords = np.vstack((xmesh.T.flatten(), ymesh.T.flatten(), np.zeros_like(ymesh.T.flatten()))).T + + x = ncoords[:, 0] + y = ncoords[:, 1] + z = ncoords[:, 2] + ncoords_flatten = ncoords.flatten() + + inner = np.logical_not(isclose(x, 0) | isclose(x, a) | isclose(y, 0) | isclose(y, b)) + np.random.seed(20) + rdm = (-1 + 2*np.random.rand(x[inner].shape[0])) + np.random.seed(20) + rdm = (-1 + 2*np.random.rand(y[inner].shape[0])) + x[inner] += dx*rdm*0.4 + y[inner] += dy*rdm*0.4 + + nids = 1 + np.arange(ncoords.shape[0]) + nid_pos = dict(zip(nids, np.arange(len(nids)))) + nids_mesh = nids.reshape(nx, ny) + n1s = nids_mesh[:-1, :-1].flatten() + n2s = nids_mesh[1:, :-1].flatten() + n3s = nids_mesh[1:, 1:].flatten() + n4s = nids_mesh[:-1, 1:].flatten() + + num_elements = len(n1s) + print('num_elements', num_elements) + + KC0r = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=INT) + KC0c = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=INT) + KC0v = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=DOUBLE) + Mr = np.zeros(data.M_SPARSE_SIZE*num_elements, dtype=INT) + Mc = np.zeros(data.M_SPARSE_SIZE*num_elements, dtype=INT) + Mv = np.zeros(data.M_SPARSE_SIZE*num_elements, dtype=DOUBLE) + KA_betar = np.zeros(data.KA_BETA_SPARSE_SIZE*num_elements, dtype=INT) + KA_betac = np.zeros(data.KA_BETA_SPARSE_SIZE*num_elements, dtype=INT) + KA_betav = np.zeros(data.KA_BETA_SPARSE_SIZE*num_elements, dtype=DOUBLE) + KA_gammar = np.zeros(data.KA_GAMMA_SPARSE_SIZE*num_elements, dtype=INT) + KA_gammac = np.zeros(data.KA_GAMMA_SPARSE_SIZE*num_elements, dtype=INT) + KA_gammav = np.zeros(data.KA_GAMMA_SPARSE_SIZE*num_elements, dtype=DOUBLE) + CAr = np.zeros(data.CA_SPARSE_SIZE*num_elements, dtype=INT) + CAc = np.zeros(data.CA_SPARSE_SIZE*num_elements, dtype=INT) + CAv = np.zeros(data.CA_SPARSE_SIZE*num_elements, dtype=DOUBLE) + + N = DOF*nx*ny + + prop = isotropic_plate(thickness=h, E=E, nu=nu, calc_scf=True, rho=rho) + + quads = [] + init_k_KC0 = 0 + init_k_M = 0 + init_k_KA_beta = 0 + init_k_KA_gamma = 0 + init_k_CA = 0 + for n1, n2, n3, n4 in zip(n1s, n2s, n3s, n4s): + pos1 = nid_pos[n1] + pos2 = nid_pos[n2] + pos3 = nid_pos[n3] + pos4 = nid_pos[n4] + r1 = ncoords[pos1] + r2 = ncoords[pos2] + r3 = ncoords[pos3] + normal = np.cross(r2 - r1, r3 - r2)[2] + assert normal > 0 + quad = Quad4R(probe) + quad.n1 = n1 + quad.n2 = n2 + quad.n3 = n3 + quad.n4 = n4 + quad.c1 = DOF*nid_pos[n1] + quad.c2 = DOF*nid_pos[n2] + quad.c3 = DOF*nid_pos[n3] + quad.c4 = DOF*nid_pos[n4] + quad.init_k_KC0 = init_k_KC0 + quad.init_k_M = init_k_M + quad.init_k_KA_beta = init_k_KA_beta + quad.init_k_KA_gamma = init_k_KA_gamma + quad.init_k_CA = init_k_CA + quad.update_rotation_matrix(ncoords_flatten) + quad.update_probe_xe(ncoords_flatten) + quad.update_KC0(KC0r, KC0c, KC0v, prop) + quad.update_M(Mr, Mc, Mv, prop) + quad.update_KA_beta(KA_betar, KA_betac, KA_betav) + quad.update_KA_gamma(KA_gammar, KA_gammac, KA_gammav) + quad.update_CA(CAr, CAc, CAv) + quads.append(quad) + init_k_KC0 += data.KC0_SPARSE_SIZE + init_k_M += data.M_SPARSE_SIZE + init_k_KA_beta += data.KA_BETA_SPARSE_SIZE + init_k_KA_gamma += data.KA_GAMMA_SPARSE_SIZE + init_k_CA += data.CA_SPARSE_SIZE + + print('elements created') + + KC0 = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc() + M = coo_matrix((Mv, (Mr, Mc)), shape=(N, N)).tocsc() + KA_beta = coo_matrix((KA_betav, (KA_betar, KA_betac)), shape=(N, N)).tocsc() + KA_gamma = coo_matrix((KA_gammav, (KA_gammar, KA_gammac)), shape=(N, N)).tocsc() + CA = coo_matrix((CAv, (CAr, CAc)), shape=(N, N)).tocsc() + + print('sparse KC0 and M created') + + bk = np.zeros(N, dtype=bool) + edges = np.isclose(x, 0.) | np.isclose(x, a) | np.isclose(y, 0) | np.isclose(y, b) + bk[0::DOF][edges] = True + bk[1::DOF][edges] = True + bk[2::DOF][edges] = True + + bu = ~bk + + Kuu = KC0[bu, :][:, bu] + Muu = M[bu, :][:, bu] + KA_betauu = KA_beta[bu, :][:, bu] + KA_gammauu = KA_gamma[bu, :][:, bu] + CAuu = CA[bu, :][:, bu] + + num_eigenvalues = 7 + print('eig solver begins') + # solves Ax = lambda M x + # we have Ax - lambda M x = 0, with lambda = omegan**2 + eigvals, eigvecsu = eigsh(A=Kuu, M=Muu, sigma=-1., which='LM', + k=num_eigenvalues, tol=1e-3) + print('eig solver end') + eigvecs = np.zeros((N, eigvecsu.shape[1]), dtype=float) + eigvecs[bu, :] = eigvecsu + omegan = eigvals**0.5 + + # panel flutter analysis + def MAC(mode1, mode2): + return (mode1@mode2)**2/((mode1@mode1)*(mode2@mode2)) + + MACmatrix = np.zeros((num_eigenvalues, num_eigenvalues)) + rho_air = 1.225 # kg/m^3 + v_sound = 343 # m/s + v_air = np.linspace(1.1*v_sound, 5*v_sound, 20) + Mach = v_air/v_sound + betas = rho_air*v_air**2/np.sqrt(Mach**2 - 1) + radius = 1. + gammas = betas/(2*radius*np.sqrt(Mach**2 - 1)) + mus = betas/(Mach*v_sound)*((Mach**2 - 2)/(Mach**2 - 1)) + mus*=0 # TODO quadractic eigenvalue problem to solve including damping + + omegan_vec = [] + for i, beta in enumerate(betas): + print('analysis i', i) + # solving generalized eigenvalue problem + Keffective = Kuu + beta*KA_betauu + gammas[i]*KA_gammauu + eigvals, eigvecsu = eigs(A=Keffective, M=Muu, + k=num_eigenvalues, which='LM', sigma=-1., tol=1e-6, v0=eigvecsu[:, 0]) + # NOTE with v0=eigvecsu[:, 0]) avoids fluctuations in adjacent + # solutions, while making it a bit slower + eigvecs = np.zeros((N, num_eigenvalues), dtype=float) + eigvecs[bu, :] = eigvecsu + + if i == 0: + eigvecs_ref = eigvecs + + corresp = [] + for j in range(num_eigenvalues): + for k in range(num_eigenvalues): + MACmatrix[j, k] = MAC(eigvecs_ref[:, j], eigvecs[:, k]) + if np.isclose(np.max(MACmatrix[j, :]), 1.): + corresp.append(np.argmax(MACmatrix[j, :])) + else: + corresp.append(j) + omegan_vec.append(eigvals[corresp]**0.5) + print(np.round(MACmatrix, 2)) + + eigvecs_ref = eigvecs[:, corresp].copy() + + + omegan_vec = np.array(omegan_vec) + + if plot: + import matplotlib + matplotlib.use('TkAgg') + import matplotlib.pyplot as plt + for i in range(num_eigenvalues): + plt.plot(Mach, omegan_vec[:, i]) + plt.ylabel('$\omega_n\ [rad/s]$') + plt.xlabel('Mach') + plt.show() + + +if __name__ == '__main__': + test_quad4r_piston_theory(plot=True, refinement=1)