python发动机悬置解耦计算-按重心处整车坐标系解耦

一、参考

http://www.doc88.com/p-6834780959259.html

二、python计算结果

三、python程序

# _*_ coding:UTF-8 _*_

import numpy as np

from numpy import cos

def M_sys(m,Ixx,Ixy,Ixz,Iyy,Iyz,Izz):

    return np.array([[m,0,0,0,0,0],\

                  [0,m,0,0,0,0],\

                  [0,0,m,0,0,0],\

                  [0,0,0,Ixx,-Ixy,-Ixz],\

                  [0,0,0,-Ixy,Iyy,-Iyz],\

                  [0,0,0,-Ixz,-Iyz,Izz]])

def K_sys_byMatrix(k_pi,k_qi,k_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri,Ai,Bi,Ci):

    #参数均采用数组格式

    n=np.size(Ai)

    ki=np.zeros((n,3,3))

    Ti=np.zeros_like(ki)

    BBi=np.zeros((n,3,6))

    Ki=np.zeros((n,6,6))

    #K=np.zeros((6,6))

    for i in range(n):

        ki[i]=np.array([[k_pi[i],0,0],\

                        [0,k_qi[i],0],\

                        [0,0,k_ri[i]]])

        Ti[i]=np.array([[cos(theta_pi[i]),cos(phi_pi[i]),cos(psi_pi[i])],\

                        [cos(theta_qi[i]),cos(phi_qi[i]),cos(psi_qi[i])],\

                        [cos(theta_ri[i]),cos(phi_ri[i]),cos(psi_ri[i])]])

        BBi[i]=np.array([[1,0,0,0,Ci[i],-Bi[i]],\

                        [0,1,0,-Ci[i],0,Ai[i]],\

                        [0,0,1,Bi[i],-Ai[i],0]])

        Ki[i]=BBi[i].T@Ti[i].T@ki[i]@Ti[i]@BBi[i]

    K=np.sum(Ki,axis=0)

    return K

def Energy_DistributionJ(M,Value,Vector_Matrix):

    #3维能量分布：（阶次、Dof、Dof)

    n=M.shape[0]

    KE_klj=np.zeros((n,n,n))    

    for j in range(n):

        Value_j=Value[j] #0维数组：特征值 频率的平方

        Vector_r=Vector_Matrix[:,j] #一维数组：行向量

        Vector_c=Vector_r.reshape(n,1) #改成二维数组：列向量

        KE_klj[j]=0.5*Value_j*M*Vector_c*Vector_r

    return KE_klj

def Energy_percent(KE_klj): 

    n=KE_klj.shape[0]

    #KE_klj 3维能量分布：（阶次、Dof、Dof=页、行、列）

    #将每行Dof的能量合并缩维，得到（各Dof总能量占比,阶）：

    Energy_perDof=np.sum(KE_klj,2)  #按行合并，第3维压缩掉，成2维数组：矩阵（阶次、各Dof能量）

    Energy_AllDof_r=np.sum(Energy_perDof,1) #按行合并，第2维压缩掉，成一维数组：行向量[1阶总能量、2阶总能量、...]

    Energy_AllDof_c=Energy_AllDof_r.reshape(n,1) #改成二维数组：列向量[[1阶总能量],[2阶总能量]、...]

    Energy_percent=100.0*Energy_perDof/Energy_AllDof_c #2维数组：矩阵（阶次、各Dof能量占比）

    Energy_percent=Energy_percent.T ########二维数组转置，得到矩阵（各Dof总能量占比,阶）

    return Energy_percent

###########input###########################################################################

[m,Ixx,Iyy,Izz,Ixy,Ixz,Iyz]=[1224.4,79.385,523.73,491.332,-5.866,-76.403,-4.698]

PG=[0.266249,-0.011964,0.186296]

ABC=[[-0.525,-0.413,0.110],[-0.525,0.413,0.110],[0.675,-0.380,0.151],[0.675,0.380,0.151]]

ABC=np.array(ABC)-np.array(PG)

Ai,Bi,Ci=ABC[:,0],ABC[:,1],ABC[:,2]

theta_pi=[0*np.pi/180,0*np.pi/180,0*np.pi/180,0*np.pi/180] #3

theta_qi=[90*np.pi/180,90*np.pi/180,90*np.pi/180,90*np.pi/180]

theta_ri=[90*np.pi/180,90*np.pi/180,90*np.pi/180,90*np.pi/180] #93

phi_pi=[90*np.pi/180,90*np.pi/180,90*np.pi/180,90*np.pi/180]

phi_qi=[0,0,0,0]

phi_ri=[90*np.pi/180,90*np.pi/180,90*np.pi/180,90*np.pi/180]

psi_pi=[90*np.pi/180,90*np.pi/180,90*np.pi/180,90*np.pi/180] #87

psi_qi=[90*np.pi/180,90*np.pi/180,90*np.pi/180,90*np.pi/180]

psi_ri=[0*np.pi/180,0*np.pi/180,0*np.pi/180,0*np.pi/180] #3

k_pi=[1386000,1386000,2745000,2745000]

k_qi=[354000,354000,367500,367500]

k_ri=[1068000,1068000,2025000,2025000]

###########end input###########################################################################

M_sys=M_sys(m,Ixx,Ixy,Ixz,Iyy,Iyz,Izz)

#[K_xxi,K_xyi,K_xzi,K_yyi,K_yzi,K_zzi]=K_element_i(k_pi,k_qi,k_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri)

#K=K_sys(K_xxi,K_xyi,K_xzi,K_yyi,K_yzi,K_zzi,Ai,Bi,Ci)

#print(K*(abs(K)>0.001))

K_sys=K_sys_byMatrix(k_pi,k_qi,k_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri,Ai,Bi,Ci)

Value, Vector_Matrix = np.linalg.eig(np.linalg.inv(M_sys)@K_sys)

#np.set_printoptions(formatter={'float': '{: 0.3f}'.format})

np.set_printoptions(precision=3, suppress=True)

print(M_sys)

print(K_sys)

print(np.sqrt(Value)/2/np.pi)

print(Vector_Matrix)

#print(Vector_Matrix.T@M_sys@Vector_Matrix)

#print(Vector_Matrix.T@K_sys@Vector_Matrix)

KE_klj=Energy_DistributionJ(M_sys,Value,Vector_Matrix)

print(Energy_percent(KE_klj))

四、virtual.lab motion模态计算验证

五、Excite PU模态计算验证

当地板水平放置，结果见上面，当地板逆时针转45度，模态结果见下图，自然频率不变。

virtual.lab motion （使用BUSH力单元）和PU（使用NONL）模态计算的结果完全一致，和理论计算接近。但是如果PU使用emo1，在地板偏转一个角度，自然频率会出现较大的变化，原因不明！！
见下图：