发动机悬置python仿真计算
发动机悬置python仿真计算
- 代码可复制到https://hub.gke2.mybinder.org/user/lijil168-requirements-l6zexquh/tree运行
- 1、发动机悬置模态及解耦参考
- 2、发动机动力学激励计算参考
- 3、用数组和矩阵两种方式计算刚度矩阵,并对比结果,原文公式推导有点点错误。python用多维数组完成矩阵运算,很简洁而且可读性好。
- 4、思路:由发动机爆压计算倾覆力矩,建立状态空间动力学模型,进行模态及解耦率计算,最后进行仿真,得到发动机质心的状态,再换算到各个悬置点。
- 5、结果
- 5.1发动机激励计算结果
- 5.1悬置系统仿真结果
- 6、完整的仿真代码
代码可复制到https://hub.gke2.mybinder.org/user/lijil168-requirements-l6zexquh/tree运行
1、发动机悬置模态及解耦参考
http://www.doc88.com/p-6364178164728.html
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()