python发动机曲轴连杆动力学计算

理论参考
理论参考

#Example2.3.3-Single cylinder,4 stroke-engine torque

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*
import pandas as pd
from scipy import interpolate

#rcParams['font.sans-serif']=['STFangsong'] #用来正常显示中文标签
#rcParams['axes.unicode_minus']=False #用来正常显示负号


#-----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
pi=np.pi
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];

pr=np.array(pr)
n=np.size(pr)
pr.resize(n,1)
crank_angle=np.array(crank_angle)
crank_angle.resize(n,1)

omega=3000;
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)
    
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)
title('offside_Torque_output  all')