配气机构的动力学分析

配气机构动力学分析课程设计
目录
一、配气机构的机构简图 ........................................ 错误!未定义书签。

二、配气机构运动学计算分析 (1)
1)配气机构中间参数法的代数分析 (1)
2）运初始值的设定及简化计算 (3)
三、配气机构动力学计算分析 (8)
1)受力分析及微分方程的建立 (8)
2）配气机构质量的换算及方程参数的计算 (10)
3)动力学微分方程的求解 (12)
四、配气机构动力学优化比较 (16)
参考文献： (23)
附件： (24)
配气机构的运动学和动力学分析
一、配气机构的机构简图
其自由度为5432352621F n p p =--=⨯-⨯-= 主动件为凸轮轴，输出件为气门。

二、配气机构的运动学计算分析
1、配气机构中间参数法的代数分析
由上面的机构简图可以得到，摇臂轴与凸轮轴的竖直位移为： 000cos cos cos cos T T T T y l l h l l h H αγαγ++=++=
化简得到：
000(cos cos )(cos cos )T T T l l h h ααγγ-+-=- （1）
摇臂轴与凸轮轴的水平位移：
00sin sin sin sin T T x l l l l H αγαγ+=+=
化简得到：
00(sin sin )(sin sin )0T l l ααγγ-+-= （2）
上面（1）（2）两式对时间求导得到
sin sin cos cos 0
T T T T dh dh l l dt d l l α
γα
γωαωγωϕωαωγ⎧
+==⋅⎪⎨
⎪--=⎩ 解得cos sin()T T h l αωγωαγ'=
- cos sin()
T h l γωα
ωαγ'=--
其中αω，γω分别为摇臂和推杆的角速度，两式对时间求导得到摇臂和推杆的角加速度为：
22
22
(cos sin )sin()cos()()cos [sin()]cos sin []sin()
cos sin()sin()
[sin()]cos cos cos()[]sin()sin()
T T T T T T T T T T T T T T T T h h l l h l h h l h l l l h h l l l γαγαωγωγωαγαγωωωγ
εαγωα
ωγαγωγαγαγαγωγωααγαγαγ''''-⋅----=
-''-
-''-=---''-+---
2
22223cos [sin()]cos cos cos()cos ()sin()sin ()T T T T T T h l h h l l ωγ
αγωγωγαγλααγαγ'
-'''-+=---
同理，得到推杆的角加速度为
22223
cos cos cos cos()()sin()sin ()
T T T h h l l γωαωγλααγελαγαγ'''+-=-+-- 其中T
l l
λ=
即为挺柱和推杆长度比 根据机构简图上的几何关系，00ββαα-=- 0(cos cos )V V l h ββ-=对时间求导可以得到
sin sin V
V V dh l l dt
βαβωβω=⋅=⋅ 22
2
(cos sin )V V d h l dt
ααβωβε=⋅+⋅ 将摇臂的角速度，角加速度带入可以得到：
cos cos sin sin sin()sin()
V V T V T T T dh l h l h dt l l ωγ
γββωαγαγ''=⋅=--
22222223
22223
cos cos cos cos()cos {cos []sin [()]}sin()sin()sin ()
cos sin ()[cos sin()cos sin ]sin()sin ()V T T T V T T T V V T T T T d h h h h l dt l l l l l h h l l ωγωγωγαγλα
ββαγαγαγωγβωγαγβλαβαγαγ''''-+=⋅+⋅----'''=+-----气门传动机构的传动比
00sin sin 1
sin()sin()V V
V V T T T T T T dh dh l l dt dt i h dh l l h h dt
ββωαγαγωω'==≈=
--'' 对中间参数进行线性近似可以得到
00000020000000000020000sin cos sin cos()
()[]
sin()sin()sin ()sin sin()
()sin()sin ()
V V T T V V
T T l l i l l l l l l βββαγαααγαγαγβαγβαααγαγ-=+-------≈
+---
2、运动初始值的设定及运动学计算的简化计算
初始参数的设定：
凸轮轴转速：1000r/min 故2104.72/60
n
rad s πω== 运动开始时推杆与竖直位置成5度角，摇臂水平且摇臂轴两端摇臂成一条直线（即机构简图中所示1OO 和2OO 在一条直线上），故05γ= 0090αβ==，
αβ=。

由此可得气门传动机构的传动比计算可简化为V
T
l i l =
气门传动机构的传动比：1.385 挺柱最大升程7.42mm 初始时刻挺柱升程00T h =
气门最大升程可计算得为10.28mm
根据摇臂的空间尺寸我们计算出摇臂两端的长度分别为：
46
47.86cos12.86cos9.67
V l mm mm ==⨯
33.2
34.55cos12.86cos9.67
T l mm mm ==⨯
由于上文的计算分析采用了中间参数（即摇臂轴和数值方向的夹角ɑ）的分析方法，因而应首先近似计算该参数值。

（1）式左侧可以近似为
00000(cos cos )(cos cos )[cos()cos()]T T l l l ααγγαγαγ-+-=---
故0
000cos()cos()T T T
h h l αγαγ--=--
根据初始设定值并画图调试得到较精确的中间参数计算值的简化计算式为
arccos(cos85)34.55
T h
α=- （3）
上式即为中间参数ɑ的计算式
下面进行配气机构各部件运动参数的简化计算，由于推杆的偏转角度γ一般比较小，因而可以近似认为00γγ==。

因T
l l ，所以λ为一个小量，又α接近于直
角，所以由推杆角速度的表达式知道γω很小（表达式中存在高阶小量cos λα） 根据高次方凸轮型线的设计得到在凸轮上升段、下降段挺柱的运动参数分别为： A 、上升缓冲段（024.5ϕ≤<）
0.25(1cos3.67)mm T h ϕ=-
0.9175sin3.67T h ϕ'= 3.367cos3.67T h ϕ''=
上升基本工作段（24.587ϕ≤<）
261422
878787877.6713.780()8.257() 2.343()0.447()62.562.562.562.5
T h mm ϕϕϕϕ----=-+-+513218787878725.265(
)45.417()30.071()9.015()62.562.562.562.5
T h ϕϕϕϕ----'=-+- 41220
87878723.161208.18(
)358.372()173.551()62.562.562.5
T h ϕϕϕ---''=-+-+ B 、下降基本工作段（87149.5ϕ≤<）
26142287
87
87
87
7.6713.780()8.257(
) 2.343(
)0.447(
)62.5
62.5
62.5
62.5
T h mm
ϕϕϕϕ----=-+-+5132187
87
87
87
25.265(
)45.417(
)30.071(
)9.015(
)62.5
62.5
62.5
62.5T h ϕϕϕϕ----'=-+-+
4122087
87
87
23.161208.18(
)358.372(
)173.551(
)62.5
62.5
62.5
T h ϕϕϕ---''=-+-+
下降缓冲段（149.5174ϕ≤<）
0.25[1cos3.67(174)]mm T h ϕ=--
0.9175sin[3.67(174)]T h ϕ'=-- 3.367cos[3.67(174)]T h ϕ''=-
下面进行各运动件的运动学简化计算，这里均取00γγ==以简化计算 A 、推杆角速度和角加速度的简化计算
cos sin()248tan[arccos(cos85)]34.55T T T h h h l γωαω
ωαγ''=-≈-
-- 222232232
23cos cos cos cos()()
sin()sin ()1()
tan sin 10.14()
34.55248tan[arccos(cos85)]sin [arccos(cos85)]34.5534.55T T T T T T T T T T h h l l h h l l h h h h γωαωγλααγελαγαγωωλαα
ωω'''+-=-+--'''≈-+'''≈-+--
代入上面的挺柱升程与凸轮轴转角关系，并用Matlab 计算绘图得到
可见摆杆的角速度很小，几乎为零，因而在复杂计算中可以忽略推杆的角速度从而简化计算。

推杆的角加速度也不大的，但要考虑其往复摆动的惯性力并进行校核计算。

B 、摇臂角速度和角加速度的简化计算
cos sin()sin 34.55sin[arccos(cos85)]34.55
T T T T T T h h h h l l αωγωω
ωαγα'''=
≈≈
-- 222232222
22cos cos cos()cos ()
sin()sin ()1()sin sin tan 34.55sin[arccos(cos85)]34.55
1
(
)
34.55sin [arccos(cos85)]tan[arccos(cos8534.55T T T T T T T T T T
T T T h h l l h h l l h h h h h αωγωγαγλαεαγαγωωαααωω'''-+=---'''≈-''≈
-'---)]
34.55
摇臂由于一直做高速往复摆动，角速度和角加速都比较高。

其加速度峰值出现在工作段起始和终止部分，需要做强度校核。

C 、气门的加速度和角加速度的简化计算
cos sin 1.381sin()V V V T T T T T
dh l l
h h h dt l l γβωωωαγ'''=≈≈- 2222223222222cos sin ()[cos sin()cos sin ]sin()sin ()
()
tan 6.668
1.381()
34.55tan [arccos(cos85)]
34.55
V V V T T T T V V T T T T T T T d h l l h h dt l l l l h h l l h h h ωγβωγαγβλαβαγαγλωωαωω'''=+-----'''≈-'''≈--
气门的加速度峰值同样出现在工作段附近，在设计气门弹簧时需要考虑其最大负惯性力，即最大负加速度。

三、配气机构的动力学计算分析
配气机构是一个多质量的振动系统，由于只需着重研究气门的运动情况，故考虑阻尼、燃气作用力、气门座刚度、气门间隙等因素建立单质量系统。

根据作用在质量M 上的力平衡条件，即可建立气门运动的微分方程。

其动力学计算模型如下图所示：
上图中整个配气机构换算到气门端的质量为M ，整个气门驱动机构的刚度为0K ，两端分别连接凸轮轴和气门。

气门弹簧刚度为K ，气门座刚度为Z K ，01 D D 是
阻尼系数。

1、受力分析及微分方程的建立
作用在质量M 上的力描述如下：
1）、气门弹簧力：1()V P h f '+ 其中1f 为气门弹簧的预紧变形量
2）、配气机构弹性变形力0()T V K ih h h δ-- 3）、气门落座时气门做的反作用力()Z V K h ∆-
其中Z K 为气门座刚度，∆为气门座在气门弹簧预紧力及燃气作用力作用下的初变形量。

1g
Z
P f F K '+∆=
4）、燃气作用力g F ，即由于气缸中气体压力与排气管中气体压力差而作用在气门头部面积上的力。

进气门可取为0，排气门可根据排气门开启前后气缸中气压与排气管气压的变化曲线求得。

5）、气门座摩擦力Z F （忽略） 6）、配气机构外阻尼力1V
dh D dt
，表现为配气机构中运动件与固定件之间的摩擦力
7）、配气机构内阻尼力0()V
T dh dh D i
dt dt
-，表现为配气机构运动件之间的摩擦 根据以上各力的平衡条件，得到以下微分方程：
20110
2()()()()V V V
T T V Z V V g Z d h dh dh dh M K ih h h K h P h f F F D D i dt dt dt dt
δ'=--+∆--+---+-
以凸轮轴转角ϕ为自变量
2V V V dh dh dh n dt d d ωπϕϕ
== 222
22222
24V V V
d h d h d h n dt d d ωπϕϕ== 代入上式得到
201222
101
[()()()]41
[()]
2V T V Z V V g Z V V T d h M K ih h h K h P h f F F d n
dh dh dh D D i n d d d δϕππϕϕϕ
'=--+∆--+-----
对上面的微分方程补充初始条件：
气门刚刚打开的瞬间凸轮轴转角0ϕϕ=，0
0V V
h h ϕϕϕϕϕ
==∂=
=∂
引入函数()
()V h u ϕϕϕ
∂=
∂带入上面的微分方程得到一阶微分方程组 01021[()()()]Z V V u D D u K K P h M h u ωϕϕωϕ
∂⎧'=-+-+++Φ⎪∂⎪
⎨
∂⎪=⎪∂⎩ （*） 式中001()()()T
T Z dh D i
K ih h K P f F d ϕωδϕ
'Φ=+-+∆-+，是关于凸轮轴转角的已知函数。

（*）式可运用经典四阶Runge-Kutta 方法进行数值求解。

初始条件为：
V V
T
h h ih u ϕϕϕϕ
ϕϕϕϕϕ
ϕ
=====∂∂==∂∂
2、配气机构质量的换算及方程参数的计算 a 、质量换算 2
22000321121
(12)[()()]K K K
M M i M i M i K K K =+++ 由于凸轮轴刚度远大于配气机构的整体刚度10K K ，故上式可简化为：
20322
()K M M M i
K =+ 挺柱的刚度较大，可以忽略挺柱对系统刚度的影响。

假设推杆在靠近挺柱一端是
不动的，则推杆端质量21
3
P M M =，P M 为推杆质量。

气门端质量321
3R V s V
I M M M l =++，V M 为气门质量，忽略弹簧上座及锁夹质量，
s M 为弹簧质量，R I 为摇臂的转动惯量。

b 、系统刚度的计算
1）、凸轮轴刚度 122
3()c b b E JL
K l L l =
-
其中凸轮轴材料弹性模数52110 1.110/c E GPa N mm ==⨯
凸轮轴横断面直径惯性矩4
4
444307.951032
32
D J mm mm ππ⨯=
=
=⨯
凸轮轴支撑轴颈内端面之间距离253mm
挺柱中心线到支撑轴颈间的距离：进气凸轮48.5ib l mm =， 排气凸轮82.5ob l mm =
故547
12222
33 1.1107.9510253/ 6.7410/()48.5(25348.5)c i b b E JL K N mm N m l L l ⨯⨯⨯⨯⨯===⨯-⨯- 54712222
33 1.1107.9510253
/ 3.3510/()82.5(25382.5)c o b b E JL K N mm N m l L l ⨯⨯⨯⨯⨯===⨯-⨯- 2）、推杆刚度 2EA
K l
=
其中推杆材料的弹性模数52100110/E GPa N mm ==⨯ 推杆横断面积222450.265A mm mm π=⨯= 推杆长度244l mm =
故57211050.265
/ 2.0610/244
EA K N mm N m l ⨯⨯===⨯ 3）、摇臂、摇臂轴及摇臂座刚度 3K 预估值定为7510/N m ⨯
4）、单质量系统刚度
进气系统换算的刚度为： 6
02222
777
32111/7.0610/11 1.385 1.385510 6.7410 2.0610i i K N m N m i i K K K =
==⨯++++
⨯⨯⨯ 排气系统换算刚度为： 602222
777
32111
/ 5.8710/11 1.385 1.385510 3.3510 2.0610o o K N m N m i i K K K =
==⨯++++
⨯⨯⨯ 5）、气门座刚度 7510/Z K N m =⨯
小组成员计算的配气机构各零件惯量参数如下：
挺柱质量135g 推杆质量90.6P M g = 摇臂质量125.8g
摇臂转动惯量换算到气门端的当量质量2275613.8
33.0147.86
R V I g g l == 进气门质量124.88iV M g = 排气门质量112.32oV M g = 气门弹簧质量106.8s M g = 所以进气系统的当量质量：
220032222
62
7
11
()()33117.0610124.88106.0833.0190.6(1.385)33 2.0610200.05i R i iV s P V K K I M M M i
M M M i K l K g g
=+=+++⨯=+⨯++⨯⨯⨯⨯=
排气系统的当量质量：
220032222
62
7
11
()()3311 5.8710112.32106.0833.0190.6(1.385)33 2.0610185.39i R o iV s P V K K I M M M i
M M M i K l K g g
=+=+++⨯=+⨯++⨯⨯⨯⨯=
c 、阻尼系数的确定
内阻尼一般很小，可取00D = 外阻尼可用阻尼比来进行估算：
Ω=
故进气系统外阻尼：
1220.1/237.68/i D s kg s ==⨯= 排气系统外阻尼：
1220.1/208.64/o D s kg s ==⨯= d 、其他参数的确定
气门弹簧预紧力定为200P N '= 气门弹簧预压缩量定为110f mm =
排气门开启前后缸内和排气管压差约为4bar
所以燃气作用力3251.2610410504g o F S p m Pa N -=∆=⨯⨯⨯= 气门座初变形量：
14
200504
0.014510
g
Z
P f F mm mm K '++∆==
=⨯ 弹簧刚度：320020000/1010N
P N m m
-'=
=⨯
忽略气门座的摩擦力Z F
气门间隙进气门0.2h mm δ=，排气门0.3h mm δ= 3、动力学微分方程的求解 初始值的设定： 气门刚刚开启时刻，气门端的位移应刚刚抵消气门间隙和气门座在弹簧预紧力及燃气作用力作用下的初变形量，即00.20.0140.214iV h h mm mm δ=+∆=+= 排气门00.30.0140.314oV h h mm mm δ=+∆=+= 对于进气系统的初始值计算如下： 此时的凸轮轴转角018.4ϕ=
V
T
h ih
u ϕϕϕϕ
ϕϕϕ
====∂==∂
通过Matlab 求解得到气门升程、气门速度和气门加速度随曲轴转角的变化关系如下：
凸轮轴转角018.4时摇臂末端并未接触气门顶部，此时仅是消除气门间隙和气
门座的预变形。

上图可以看出，进气门在气门开启和气门落座过程中均有速度波动，在凸轮型线中缓冲段和工作段的交界处也有速度波动，但波动量几乎不可见。

下图为进气系统各类曲线的总览
排气门系统的初始值计算如下： 凸轮轴转角023.07ϕ
=
V
T
h ih u ϕϕϕϕ
ϕϕϕ
====∂==∂
通过Matlab 求解得到气门升程、气门速度和气门加速度随曲轴转角的变化关系如下：
同理，凸轮轴转角023.07时排气系统的摇臂未接触排气门顶部，此时仅是消
除气门间隙和气门座的预变形。

上图可以看出，同样在排气门起座和落座的时，气门出现速度波动，由于燃气作用力的影响等，排气门的速度波动明显大于进气门。

在凸轮型线缓冲段和工作段的交界处，也存在速度波动，但相对而言波动很小。

下图为排气系统各类曲线：
四、配气机构动力学优化比较
进气门速度的理想运动学曲线和动力学分析的速度曲线比较
进气门最大升程10.3496mm
进气门最大速度1.693m/s
进气门最大加速度830.6730m/s^2
动力学升程的响应慢于运动学，且由于气门间隙的影响，运动学峰值高出约
0.2mm
排气门最大升程10.1789mm
排气门最大速度1.7m/s
排气门最大加速度835.9052m/s^2
同样，动力学升程的响应慢于运动学，且由于气门间隙的影响，运动学峰值高出约0.3mm
下面为运用2100发动机原有凸轮型线得到的进排气系统升程曲线、速度曲线和加速度曲线。

进气系统：
排气系统：
通过对比原2100配气系统动力学循环曲线和运用高次方凸轮型线后的2100配气系统动力学循环曲线可以得到以下结论：
1、在进排气升程变化不大的情况下，采用多项式高次方凸轮型线可以显著的减少进排气系统的速度和加速度幅值，速度由2.5m/s左右降至1.7m/s左右，加速
度幅值由1900m/s^2左右降至800m/s^2左右。

2、采用高次方凸轮型线使得进排气系统在开启和落座时刻的速度、加速度波动幅度明显减少，速度几乎没有波动。

这样在使用材料和加工工艺相同的情况下保证了配气系统具备较高的使用寿命，较低是噪声，同时可以降低了成本。

参考文献
[1]李佳, 刘震涛, 刘忠民, 等. 内燃机配气机构系统刚度的虚拟设计[J]. 农业工程学报, 2012, 28(6): 44-49.
[2]赵冬青, 苏铁熊, 赵振锋, 等. 顶置凸轮轴配气机构运动学和动力学计算[J]. 车用发动机, 2004 (6): 1-3.
[3]史绍熙. 柴油机设计手册[J]. 北京: 中国农业机械出版社, 1984.
[4]高文志, 李明海, 袁文华, 等.内燃机课程设计[J]. 北京: 中国水利水电出版社, 2010.8
[5]内燃机设计
附件
运动学计算的Matlab代码
x=linspace(0,174,1740);
y=0.25*(1-cosd(3.67*x)).*(x>=0&x<24.5)+(7.67-13.780*((87-x)/62.5).^2+8.257 *((87-x)/62.5).^6-2.343*((87-x)/62.5).^14+0.447*((87-x)/62.5).^22).*(x>=24.5& x<87)+(7.67-13.780*((x-87)/62.5).^2+8.257*((x-87)/62.5).^6-2.343*((x-87)/62.
5).^14+0.447*((x-87)/62.5).^22).*(x>=87&x<149.5)+0.25*(1-cosd(3.67*(174-x) )).*(x>=149.5&x<174);
y1=0.9175*sind(3.67*x).*(x>=0&x<24.5)+(25.265*((87-x)/62.5)-45.417*((87-x)/ 62.5).^5+30.071*((87-x)/62.5).^13-9.015*((87-x)/62.5).^21).*(x>=24.5&x<87)+ (-25.265*((x-87)/62.5)+45.417*((x-87)/62.5).^5-30.071*((x-87)/62.5).^13+9.01 5*((x-87)/62.5).^21).*(x>=87&x<149.5)+(-0.9175*sind(3.67*(174-x))).*(x>=149 .5&x<174);
y2=3.367*cosd(3.67*x).*(x>=0&x<24.5)+(-23.161+208.18*((87-x)/62.5).^4-358 .372*((87-x)/62.5).^12+173.551*((87-x)/62.5).^20).*(x>=24.5&x<87)+(-23.161 +208.18*((x-87)/62.5).^4-358.372*((x-87)/62.5).^12+173.551*((x-87)/62.5).^20 ).*(x>=87&x<149.5)+(3.367*cosd(3.67*(174-x))).*(x>=149.5&x<174);
wt=104.72*y1./(248*tan(acos(cosd(85)-y/34.55)));
plot(x,wt,'b')
et=-y2*104.72^2./(248*tan(acos(cosd(85)-y/34.55)))+0.14*(104.72*y1/34.55).^ 2./(sin(acos(cosd(85)-y/34.55))).^3;
plot(x,et,'r')
wy=y1*104.72./(34.55*sin(acos(cosd(85)-y/34.55)));
plot(x,wy,'b')
ey=y2*104.72^2./(34.55*sin(acos(cosd(85)-y/34.55)))-(104.72*y1/34.55).^2./(( sin(acos(cosd(85)-y/34.55))).^2.*tan(acos(cosd(85)-y/34.55)));
plot(x,ey,'r')
wq=1.381*104.72*y1;
plot(x,wq,'b')
aq=1.381*y2*104.72^2-(104.72*y1/34.55).^2*6.668./((tan(acos(cosd(85)-y/34.
55))).^2);
plot(x,aq,'r')。