代码可复制到https://hub.gke2.mybinder.org/user/lijil168-requirements-l6zexquh/tree运行

1、发动机悬置模态及解耦参考

2、发动机动力学激励计算参考

3、用数组和矩阵两种方式计算刚度矩阵，并对比结果，原文公式推导有点点错误。python用多维数组完成矩阵运算，很简洁而且可读性好。

4、思路：由发动机爆压计算倾覆力矩，建立状态空间动力学模型，进行模态及解耦率计算，最后进行仿真，得到发动机质心的状态，再换算到各个悬置点。

5、结果

分别为质量矩阵、刚度矩阵、自然频率、模态振型矩阵(每列对应振型向量)、解耦率(每列对应模态阶次，行从上到下对应x,y,z,theta_x,theta_y,theta_z 6个自由度)、阻尼矩阵。

5.1发动机激励计算结果

5.1悬置系统仿真结果

6、完整的仿真代码

​

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

import numpy as np

from numpy import sin

from numpy import cos

from numpy import arcsin as asin

from numpy import sqrt

from numpy import tan

from matplotlib.pyplot import*

from scipy import interpolate

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_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_pi:第i个悬置的p向刚度；

#theta_pi：p轴与系统x轴的夹角；

#phi_pi：p轴与系统y轴的夹角；

#psi_pi:p轴与系统z轴的夹角；

K_xxi=k_pi*cos(theta_pi)**2+k_qi*cos(theta_qi)**2+k_ri*cos(theta_ri)**2

K_yyi=k_pi*cos(phi_pi)**2+k_qi*cos(phi_qi)**2+k_ri*cos(phi_ri)**2

K_zzi=k_pi*cos(psi_pi)**2+k_qi*cos(psi_qi)**2+k_ri*cos(psi_ri)**2

K_xyi=k_pi*cos(theta_pi)*cos(phi_pi)+k_qi*cos(theta_qi)*cos(phi_qi)+k_ri*cos(theta_ri)*cos(phi_ri)

K_xzi=k_pi*cos(theta_pi)*cos(psi_pi)+k_qi*cos(theta_qi)*cos(psi_qi)+k_ri*cos(theta_ri)*cos(psi_ri)

K_yzi=k_pi*cos(psi_pi)*cos(phi_pi)+k_qi*cos(psi_qi)*cos(phi_qi)+k_ri*cos(psi_ri)*cos(phi_ri)

return [K_xxi,K_xyi,K_xzi,K_yyi,K_yzi,K_zzi]

def K_sys(K_xxi,K_xyi,K_xzi,K_yyi,K_yzi,K_zzi,Ai,Bi,Ci):

Ai=np.array(Ai)

Bi=np.array(Bi)

Ci=np.array(Ci)

K_xx=sum(K_xxi)

K_xy=sum(K_xyi)

K_xz=sum(K_xzi)

K_yy=sum(K_yyi)

K_yz=sum(K_yzi)

K_zz=sum(K_zzi)

K_alpha_alpha=sum(K_yyi*Ci**2)+sum(K_zzi*Bi**2)-2*sum(K_yzi*Bi*Ci) ###

K_beta_beta=sum(K_xxi*Ci**2)+sum(K_zzi*Ai**2)-2*sum(K_xzi*Ai*Ci)

K_gamma_gamma=sum(K_xxi*Bi**2)+sum(K_yyi*Ai**2)-2*sum(K_xyi*Ai*Bi)   ####原公式有误！！！

K_alpha_beta=sum(K_xzi*Bi*Ci)+sum(K_yz*Ai*Ci)-sum(K_zzi*Ai*Bi)-sum(K_xyi*Ci**2)

K_beta_gamma=sum(K_xyi*Ai*Ci)+sum(K_xz*Ai*Bi)-sum(K_xxi*Ci*Bi)-sum(K_yzi*Ai**2) ###

K_alpha_gamma=sum(K_xyi*Ci*Bi)+sum(K_yz*Ai*Bi)-sum(K_yyi*Ci*Ai)-sum(K_xzi*Bi**2)

K_x_alpha=sum(K_xzi*Bi)-sum(K_xyi*Ci)

K_x_beta=sum(K_xxi*Ci)-sum(K_xzi*Ai)

K_x_gamma=sum(K_xyi*Ai)-sum(K_xxi*Bi)

K_y_alpha=sum(K_yzi*Bi)-sum(K_yyi*Ci)

K_y_beta=sum(K_xyi*Ci)-sum(K_yzi*Ai)

K_y_gamma=sum(K_yyi*Ai)-sum(K_xyi*Bi)

K_z_alpha=sum(K_zzi*Bi)-sum(K_yzi*Ci)

K_z_beta=sum(K_xzi*Ci)-sum(K_zzi*Ai)

K_z_gamma=sum(K_yzi*Ai)-sum(K_xzi*Bi)

K=np.array([[K_xx,K_xy,K_xz,K_x_alpha,K_x_beta,K_x_gamma],\

[K_xy,K_yy,K_yz,K_y_alpha,K_y_beta,K_y_gamma],\

[K_xz,K_yz,K_zz,K_z_alpha,K_z_beta,K_z_gamma],\

[K_x_alpha,K_y_alpha,K_z_alpha,K_alpha_alpha,K_alpha_beta,K_alpha_gamma],\

[K_x_beta,K_y_beta,K_z_beta,K_alpha_beta,K_beta_beta,K_beta_gamma],\

[K_x_gamma,K_y_gamma,K_z_gamma,K_alpha_gamma,K_beta_gamma,K_gamma_gamma]])

return K

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

def PointSensor(Yxyzi,x_p,y_p,z_p):

#将刚体质心加速度Yxyzi根据传感器的的坐标转换到Pxyzi：

TransMatrix=np.array([[1,0,0,0,z_p,y_p],\

[0,1,0,z_p,0,x_p],\

[0,0,1,y_p,x_p,0]])

Pxyzi=Yxyzi@TransMatrix.T

return Pxyzi

def EngineSim(mPiston,mConnectingRod,l,lB,R,Ap,IC,offside_cylinder,pr,crank_angle,omega):

#mPiston=430/1000;               #kg

#mConnectingRod=440/1000;        #kg

#l=140/1000;                 #m,Connecting rod pin-pin length

#lB=37/1000;             #m,Connecting rod pin B-CG lengh

#R=49/1000;                  #m,Crank radius

#Ap=5800;               #mm2,Pistion area

#IC=0.0015;                  #kg.m2,Connecting rod inertia

#offside_cylinder=0/1000   #气缸与曲轴偏置距  add by lijilin 2020.12.23

#num_cylinders=6

#-------Pressure in combustion chamber-------

#pr=[18,32,32.5,32,20,15,10,8,6,5,3,1.2,0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.0,2.0,4,9,15,18];

#crank_angle=[0,20,23,26,50,60,70,80,100,110,150,190,200,220,230,250,280,300,330,360,380,400,440,460,480,500,540,600,630,660,690,710,720];

#omega=3000; #rpm speed of engine

pi=np.pi

pr=np.array(pr)

n=np.size(pr)

pr.resize(n,1)

crank_angle=np.array(crank_angle)

crank_angle.resize(n,1)

omeg=omega*pi/30;

Rl=R/l;       #define ratio R over l

lA=l-lB;

ang=crank_angle*pi/180;

sa=sin(ang);

ca=cos(ang);

s2a=sin(2*ang);

c2a=cos(2*ang);

#beta=asin((Rl*sa));

beta=asin((R*sa-offside_cylinder)/l)  #####update by lijilin 2020.12.23

ka=ca+Rl*c2a/cos(beta)+Rl**3*s2a**2/cos(beta)**3/4;

#//Piston acceleration:

aP=R*ka*omeg**2;

#//onnecting rod angular acceleration:

alpha_c=Rl*omeg**2*sa/cos(beta);

#//Connecting rod GC acceleration:

k3=lA*sa/l;

k4=ca+Rl*c2a*lB/cos(beta)/l;

agx=-R*omeg**2*k3;

agy=-R*omeg**2*k4;

ag=sqrt(agx**2+agy**2);

#//Piston pressure force:

Fp=pr*Ap/9.8;

#//Piston inertia force:

FIP=-mPiston*aP;

#//Resultant piston force:

FPt=Fp+FIP;

#//Connection rod inertia forces(N):

FIx=-mConnectingRod*agx;

FIy=-mConnectingRod*agy;

#//Conneting rod inertia torque(Nm):

TIG=IC*alpha_c;

#//Crank-pin bearing forces:

Bx=(-FIP-Fp+lB*FIy/l)*tan(beta)-lA*FIx/l-TIG/l/cos(beta);

By=Fp+FIP-FIy;

#//Engine torque:

Te=R*(By*sa-Bx*ca);

#//气缸壁侧压力：

Fw=-FIx-Bx;

#倾覆力矩：

T_w=-Fw*(R*ca+l*cos(beta))

T_offside=-By*offside_cylinder

T_all=T_w+T_offside

crank_angle_i=np.linspace(0,719,720)

#"nearest","zero"为阶梯插值

#slinear 线性插值

#"quadratic","cubic" 为2阶、3阶B样条曲线插值

# 'slinear’, 'quadratic’ and 'cubic’ refer to a spline interpolation of first, second or third order)

func_interp=interpolate.interp1d(crank_angle[:,0],T_all[:,0],"cubic")

T_all_i=func_interp(crank_angle_i).reshape(720,1)

psi_shift=int(720/num_cylinders)

T_shift=np.zeros((720,num_cylinders))

T_shift[:,0]=T_all_i[:,0]

for i in range(num_cylinders-1):

i=i+1

T_shift[:-psi_shift*i,i]=T_all_i[psi_shift*i:,0]

T_shift[-psi_shift*i+1:,i]=T_all_i[:psi_shift*i-1,0]

T_output=np.sum(T_shift,1)

return [Fp,FIP,FPt,FIx,FIy,crank_angle,TIG,Bx,By,crank_angle,Te,crank_angle,Fw,T_w,T_offside,T_all,crank_angle_i,T_shift,T_output]

###===============================================================================================================

#######################begin:计算发动机激励#################################

#-----input-----

mPiston=430/1000;               #kg

mConnectingRod=440/1000;        #kg

l=140/1000;                 #m,Connecting rod pin-pin length

lB=37/1000;             #m,Connecting rod pin B-CG lengh

R=49/1000;                  #m,Crank radius

Ap=5800;               #mm2,Pistion area

IC=0.0015;                  #kg.m2,Connecting rod inertia

offside_cylinder=0/1000   #气缸与曲轴偏置距  add by lijilin 2020.12.23

num_cylinders=6

omega=3000

#-------Pressure in combustion chamber--------

pr=[18,32,32.5,32,20,15,10,8,6,5,3,1.2,0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.0,2.0,4,9,15,18];

crank_angle=[0,20,23,26,50,60,70,80,100,110,150,190,200,220,230,250,280,300,330,360,380,400,440,460,480,500,540,600,630,660,690,710,720];

[Fp,FIP,FPt,FIx,FIy,crank_angle,TIG,Bx,By,crank_angle,Te,crank_angle,Fw,\

T_w,T_offside,T_all,crank_angle_i,T_shift,T_output]=EngineSim(mPiston,mConnectingRod,l,lB,R,Ap,IC,offside_cylinder,pr,crank_angle,omega)

#-----output------

figure(figsize=(20,20))

subplot(3,3,1)

plot(crank_angle,np.c_[Fp,FIP,FPt]);#活塞气压力 惯性力  合力

title('Piston pressure force Piston inertia force  Resultant piston force')

subplot(3,3,2)

plot(crank_angle,np.c_[FIx,FIy]);#连杆惯性力x   y

title('Connection rod inertia forces(N)x   y')

subplot(3,3,3)

plot(crank_angle,TIG);#惯性力矩

title('Conneting rod inertia torque(Nm)')

subplot(3,3,4)

plot(crank_angle,np.c_[Bx,By]);#曲柄销力x y

title('Crank-pin bearing forcesx y')

subplot(3,3,5)

plot(crank_angle,Te);#发动机力矩

title('Engine torque')

subplot(3,3,6)

plot(crank_angle,Fw);#气缸侧压

title('offside_Force')

subplot(3,3,7)

plot(crank_angle,np.c_[T_w,T_offside,T_all])

title('offside_Torque_w z  all')

subplot(3,3,8)

plot(crank_angle_i,T_shift[:,:])

title('offside_Torque all')

subplot(3,3,9)

plot(crank_angle_i,T_output,label='$mean: %f$' %np.mean(T_output)) ####

title('offside_Torque_output  all')

legend()

#######################end:计算发动机激励#################################

##==================================================================================================================

#######################begin:计算动力总成悬置模态及解耦率#################################

[m,Ixx,Iyy,Izz,Ixy,Ixz,Iyz]=[153.2,8.2,3.79,7.79,0.83,0.22,1.08]

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

Ai=[0.01953,-0.00631,0.10237]

Bi=[-0.39667,0.51314,-0.01586]

Ci=[0.05429,0.19266,-0.24009]

theta_pi=[0,0,0]

theta_qi=[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]

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

phi_qi=[0,0,0]

phi_ri=[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]

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

psi_ri=[0,0,0]

k_pi=[196000,154000,210000]

k_qi=[49000,154000,28000]

k_ri=[154000,189000,28000]

#[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))

#######################end:计算动力总成悬置模态及解耦率#################################

##==================================================================================================================

##############begin:建立状态空间方程并用发动机倾覆力矩激励进行仿真 state space function########################

import matplotlib.pyplot as plt

from scipy.signal import lti, lsim

import pandas as pd

#plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签

#plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

c_pi=[100,100,100]

c_qi=[100,100,100]

c_ri=[100,100,100]

C_sys=K_sys_byMatrix(c_pi,c_qi,c_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri,Ai,Bi,Ci)

print(C_sys)

M=M_sys

K=K_sys

C=C_sys

n=6

G=np.r_[np.c_[C,M],np.c_[M,np.zeros((n,n))]] #G=[C,M;M,O]

H=np.r_[np.c_[K,np.zeros((n,n))],np.c_[np.zeros((n,n)),-M]] #H=[K,O;O,-M]

ssA=-np.linalg.inv(G)@H;

ssB=np.linalg.inv(G)@np.row_stack([np.eye(n),np.zeros((n,n))]); #G dX +H X =E u => dX="-G\H" X +"G\E" u

#ssB=[zeros(n),inv(M);inv(M),-inv(M)^2]*[eye(n);zeros(n)];

#ssC=[-M\K,-M\C];

ssC=np.column_stack([-np.linalg.inv(M)@K,-np.linalg.inv(M)@C])

ssD=np.linalg.inv(M); #设输出Y= d^2X=[-M\K,-M\C]*[X,dX]+[M]\U

sys=lti(ssA,ssB,ssC,ssD);

#------------从文件中读入发动机倾覆力矩激励数据：---------

#data = pd.read_csv('testdata.txt', header=None, )

#data.head(10)

#Fs=data.values[0]

#y=data.values[1:]

#t=np.linspace(0,1/Fs[0]*(np.size(y)-1), num=np.size(y))  #Fs[0]才是采样频率的值！！！

#------------使用计算的发动机倾覆力矩激励数据：---------

t=crank_angle_i/360/omega*60

y=T_output

y.resize(np.size(y))

y1=(y-np.mean(y)) ###########

U=np.zeros((np.size(y),6))

U[:,3]=y1

X0=np.zeros((n*2)) #X0=zeros(2*n,1);%X0=[x1_0,...,dx1_0,...]

tout,Y,X=lsim(sys,U,t,X0);

#-------计算悬置点的加速度：

Pxyzi1=PointSensor(Y,Ai[0],Bi[0],Ci[0])

Pxyzi2=PointSensor(Y,Ai[1],Bi[1],Ci[1])

Pxyzi3=PointSensor(Y,Ai[2],Bi[2],Ci[2])

plt.figure(figsize=(20,20))

plt.subplot(3,3,1)

plt.plot(t,U) #######

plt.xlabel('Torque of engine(Nm)');

#grid on

plt.grid(alpha=0.3)

plt.subplot(3,3,2)

plt.plot(t,X)  ######

plt.xlabel('6 DOFs of position,speed of engine');

plt.subplot(3,3,3)

for i in range(3):

plt.plot(t,Y[:,i],label='$Dof:%s$' % ["x","y","z","theta_x","theta_y","theta_z"][i])

plt.xlabel('G acc. of speed/(mps^2)');

plt.legend()

plt.subplot(3,3,4)

for i in [3,4,5]:

plt.plot(t,Y[:,i],label='$Dof:%s;;rms: %f$' %(["x","y","z","theta_x","theta_y","theta_z"][i],np.std(Y[:,i])))

plt.xlabel('G acc. of ang.speed/(radps^2)');

plt.legend()

plt.subplot(3,3,5)

for i in range(3):

plt.plot(t,Pxyzi1[:,i],label='$Dof:%s;;rms: %f$' %(["x","y","z"][i],np.std(Pxyzi1[:,i])))

plt.xlabel('Mount1 acc. of speed/(mps^2)');

plt.legend()

plt.subplot(3,3,6)

for i in range(3):

plt.plot(t,Pxyzi2[:,i],label='$Dof:%s;;rms: %f$' %(["x","y","z"][i],np.std(Pxyzi2[:,i])))

plt.xlabel('Mount2 acc. of speed/(mps^2)');

plt.legend()

plt.subplot(3,3,7)

for i in range(3):

plt.plot(t,Pxyzi3[:,i],label='$Dof:%s;;rms: %f$' %(["x","y","z"][i],np.std(Pxyzi3[:,i])))

plt.xlabel('Mount3 acc. of speed/(mps^2)');

plt.legend()

plt.show()