凸轮机构的弹性动力学

凸轮机构的弹性动力学分析（附MATLAB 代码）

【问题】已知一凸轮系统，欲使其考虑弹性因素后从动件的真实运动规律按照余弦加速度运动规律运动，建立该凸轮系统的弹性动力学模型，分析其未考虑弹性因素时从动件的运动规律，并绘制出从动件的理论运动规律及考虑弹性因素后的真实运动规律。凸轮系统的运动及动力参数自定。程序代码需提供电子版，并说明运行环境。 【解答】

一、建立动力学模型

取图1所示的凸轮机构为研究对象，图2为其所对应的动力学模型。

图1：凸轮机构运动简图 图2：凸轮机构的动力学模型

为使得问题简化，力学模型中忽略了凸轮轴的扭转变形、弯曲变形以及回位弹簧的阻尼作用。图2中k 为系统等效弹簧的刚度，c 为凸轮机构从动组件的阻尼系数，k h 为回位弹簧的刚度，F 0为回位弹簧的预紧力，M 为凸轮机构在从动件侧的当量质量，x 为与凸轮廓线有关的等效凸轮升程（图中所示的凸轮并非真正的凸轮，其廓线对应的升程与真实凸轮廓线对应的升程x 0具备关系x rx 0，其中r 为摇臂比。因为x 与x 0仅相差一个比例系数r ，为了便，y 为于叙述，后文将只注重分析x 与从动件输出的关系，而不再专门区别x 与x 0的差异）从动件的实际升程。

二、建立动力学方程

该机构的自由度为1，利用牛顿第二定律建立运动微分方程：

d 2y dy dx

M 2=-k (y -x ) -c (-) -k h y -F 0 （式1）

dt dt dt

设凸轮转动的角速度为ω，它与时间微分dt 、凸轮转角微分d ϕ具有关系：

dt =

d ϕ

ω

（式2）

将（式2）代入（式1）并整理可得：

d 2y dy dx M ω+c ω+(k +k ) y =c ω+kx -F 0 （式3） h 2

d ϕd ϕd ϕ

2

微分方程（式3）有两层含义：①若已知从动件的真实运动规律，可求解出凸轮在高速运转条件下考虑弹性变形影响的理论轮廓；②若已知凸轮廓线，可求解考虑弹性变形的从动件的动力学响应。

三、运动方程的求解 （一）凸轮轮廓的设计

kg ，凸轮的转速n =1200r /min ，c =55. 7148N ⋅s /m ，已知条件如下：M =0. 08459

为避免余弦加速度运动产生的冲击，k h =10400N /m ，k =3194800N /m ，F 0=400N ；取凸轮的推程运动角和回程运动角均为180，远休止角和近休止角均为0，从动件的最大升程h =6. 2mm 。

( 1 - cos ϕ ) 根据已知条件，可以确定从动件的位移方程 y =将上式代入（式3）可得：

h

2

h h 1dx [M ω2-(k +k h )]cos ϕ+c ωsin ϕ+(k +k h ) h +F 0=c ω+kx （式4） 222d ϕ

由于（式4）对应的常微分方程难以求出解析解，这里利用MATLAB 求解出其数值解并与位移方程比较如下图：

-3

升程/m

凸轮的转角/rad

若需要考虑机构的弹性变形，凸轮的轮廓应按照上图中的红色曲线进行设计。

（二）由已知廓线求解从动件的真实运动规律

由于系统的弹性变形，从动件的真实运动y 与等效凸轮升程x 不再相等，当然，从动件的真实运动速度、加速度与对应的理论值也不相等。

由于弹簧的预紧力F 0为常数，它只影响系统振动的初始平衡位置，故在分析从动件的运动规律时不再考虑，从而（式1）被简化为：

d 2y dy dx

M 2+c +(k +k h ) y =c +kx （式5）

dt dt dt

根据振动理论，系统自由振动的固有频率ωn =

k +k h

(1-ξ2) M

其中阻尼比

ξ=

c

2M (k +k h )

代入相关数值计算可得

ωn =6147rad /s ，ξ=0. 0535

如果从动件按照前述的余弦加速度规律运动，则 x =将上式代入（式5）可得：

h

(1-cos ϕ) 2

d 2y dy dx M ω+c ω+(k +k ) y =c ω+kx （式6）

h

d ϕ2d ϕd ϕ

2

利用MATLAB 求解（式6）的数值解，结果如下（图中的n =

-3

ωn

）： ω

位移/m

凸轮的转角/rad

-3

速度/（m /s ）

凸轮的转角/rad

-3

加速度/（m /s 2）

凸轮的转角/rad

从上图中的位移、速度、加速度分析看，当考虑到系统的弹性和阻尼后，工作端的运动规律发生改变，x ≠y 。只有当n 很大时，也就是说当系统的固有频率ωn 很大时（刚度大），而

且凸轮的角速度很小时，y 才接近x 值。

程序附录：

% filename: tulun.m

% function: cooperate with jisuanlunkuo.m to calculate the curve of the cam

function dx=tulun(a,x) %% 凸轮机构的结构参数

h=6.2*10^(-3); % 升程

m=0.08459; % 凸轮机构的当量质量 n=1200;w=n/60*2*pi; % 凸轮转速及角速度 c=55.7148; % 阻尼 k1=10400; % 回位弹簧刚度 k=3194800; % 系统等效弹簧刚度 F=400; % 回位弹簧预紧力

%% 理论轮廓的微分方程

dx=(h/2*(m*w^2-k-k1)*cos(a)+h/2*c*w*sin(a)+0.5*(k+k1)*h+F-k*x)/(c*w);

% filename: jisuanlunkuo.m

% function: cooperate with tulun.m to calculate the curve of the cam %% 计算考虑弹性变形的凸轮的理论轮廓曲线 [a,x]=ode45('tulun' ,[0 2*pi],0); plot(a,x,'r' ); hold on ;

%% 作出未考虑弹性变形的凸轮的理论轮廓曲线 a1=linspace(0,2*pi); h=6.2*10^(-3); y=h/2*(1-cos(a1)); plot(a1,y,'-b' ); %% 添加标注 grid on ;

legend(' 考虑弹性变形后凸轮的理论轮廓' , ' 未考虑弹性变形的凸轮理论轮廓' ); xlabel(' 凸轮的转角/rad'); ylabel(' 升程/m');

% filename: yundongfenxi.m

% function: cooperate with fenxi.m to calculate the displacement,velocity % and acceleration of the cam mechanism when the elastic deformation of % the cam is not considered in design function dy=yundongfenxi(a,y) %% 凸轮机构的结构参数 h=6.2*10^(-3); m=0.08459;

c=55.7148; k1=10400; w=1537;

% 不同情况下将凸轮的角速度值分别设为1537,2049,3074 k=3194800; %% 运动分析微分方程 dy=zeros(2,1); dy(1)=y(2);

dy(2)=(c*w*h/2*sin(a)+k*h/2*(1-cos(a))-(k+k1)*y(1)-c*w*y(1))/(m*w*w);

% filename: fenxi.m

% function: cooperate with yundongfenxi.m to calculate the displacement,velocity

% and acceleration of the cam mechanism when the elastic deformation of % the cam is not considered in design %% 绘制从动件的位移曲线

[a,y]=ode45('yundongfenxi' ,[0 2*pi],[0 0]); n=length(a); figure(1);

plot(a,y(:,1),'--b' );

% 不同情况下将线型分别设为--b ，--g ，--k hold on ;

%% 绘制从动件的速度曲线 figure(2); a1=zeros(n-1,1); for i=1:(n-1) a1(i)=a(i); end

y1=diff(y(:,1))./diff(a); plot(a1,y1,'--b' );

% 不同情况下将线型分别设为--b ，--g ，--k hold on ;

%% 绘制从动件的加速度曲线 figure(3); a2=zeros(n-2,1); for i=2:(n-1) a2(i-1)=a1(i); end

y2=diff(y1(:,1))./diff(a1); plot(a2,y2,'--b' );

% 不同情况下将线型分别设为--b ，--g ，--k hold on ;

% fliename: zuihou.m

% function: give the ideal displacement,velocity and acceleration of the % cam mechanism %% 绘制理想位移曲线 a0=linspace(0,2*pi); n=length(a0);

y0=3.1*10^(-3)*(1-cos(a0)); figure(1); plot(a0,y0,'r' ); hold on ;

%% 绘制理想速度曲线 figure(2); a01=zeros(1,n-1); for i=1:(n-1) a01(i)=a0(i); end

y01=diff(y0)./diff(a0); plot(a01,y01,'r' ); hold on ;

%% 绘制理想加速度曲线 figure(3); a02=zeros(1,n-2); for i=2:(n-1) a02(i-1)=a01(i); end

y02=diff(y01)./diff(a01); plot(a02,y02,'r' ); hold on ; %% 添加标注 figure(1);

xlabel(' 凸轮的转角/rad'); ylabel(' 位移/m');

legend('n=4', 'n=3', 'n=2', ' 标准输出' ); figure(2);

xlabel(' 凸轮的转角/rad'); ylabel(' 速度/（m/s）' );

legend('n=4', 'n=3', 'n=2', ' 标准输出' ); figure(3);

xlabel(' 凸轮的转角/rad'); ylabel(' 加速度/（m/s^2）' );

legend('n=4', 'n=3', 'n=2', ' 标准输出' );

程序运行说明：在MATLAB 主窗口中运行jisuanlunkuo.m 文件可以得到考虑弹性变形后凸轮的理论轮廓曲线和未考虑弹性变形的凸轮理论轮廓曲线对比图；初次运行fenxi.m 可以得到

凸轮角速度

ω=1537rad /s

情况下从动件的位移、速度和加速度图，之后需要按程序说明

中的位置将yundongfenxi.m 文件中的角速度分别改为

2049rad /s

和

3074rad /s

并相应地

将fenxi.m 文件中的线型改为--g 和--k ，然后再次运行fenxi.m 文件可以得到以及

ω=2049rad /s

ω=3074rad /s

情况下从动件的位移、速度和加速度图（每次修改后注意保存），最后

运行zuihou.m 文件得到最终的从动件运动量的对比分析图。

凸轮机构的弹性动力学分析（附MATLAB 代码）

【问题】已知一凸轮系统，欲使其考虑弹性因素后从动件的真实运动规律按照余弦加速度运动规律运动，建立该凸轮系统的弹性动力学模型，分析其未考虑弹性因素时从动件的运动规律，并绘制出从动件的理论运动规律及考虑弹性因素后的真实运动规律。凸轮系统的运动及动力参数自定。程序代码需提供电子版，并说明运行环境。 【解答】

一、建立动力学模型

取图1所示的凸轮机构为研究对象，图2为其所对应的动力学模型。

图1：凸轮机构运动简图 图2：凸轮机构的动力学模型

为使得问题简化，力学模型中忽略了凸轮轴的扭转变形、弯曲变形以及回位弹簧的阻尼作用。图2中k 为系统等效弹簧的刚度，c 为凸轮机构从动组件的阻尼系数，k h 为回位弹簧的刚度，F 0为回位弹簧的预紧力，M 为凸轮机构在从动件侧的当量质量，x 为与凸轮廓线有关的等效凸轮升程（图中所示的凸轮并非真正的凸轮，其廓线对应的升程与真实凸轮廓线对应的升程x 0具备关系x rx 0，其中r 为摇臂比。因为x 与x 0仅相差一个比例系数r ，为了便，y 为于叙述，后文将只注重分析x 与从动件输出的关系，而不再专门区别x 与x 0的差异）从动件的实际升程。

二、建立动力学方程

该机构的自由度为1，利用牛顿第二定律建立运动微分方程：

d 2y dy dx

M 2=-k (y -x ) -c (-) -k h y -F 0 （式1）

dt dt dt

设凸轮转动的角速度为ω，它与时间微分dt 、凸轮转角微分d ϕ具有关系：

dt =

d ϕ

ω

（式2）

将（式2）代入（式1）并整理可得：

d 2y dy dx M ω+c ω+(k +k ) y =c ω+kx -F 0 （式3） h 2

d ϕd ϕd ϕ

2

微分方程（式3）有两层含义：①若已知从动件的真实运动规律，可求解出凸轮在高速运转条件下考虑弹性变形影响的理论轮廓；②若已知凸轮廓线，可求解考虑弹性变形的从动件的动力学响应。

三、运动方程的求解 （一）凸轮轮廓的设计

kg ，凸轮的转速n =1200r /min ，c =55. 7148N ⋅s /m ，已知条件如下：M =0. 08459

为避免余弦加速度运动产生的冲击，k h =10400N /m ，k =3194800N /m ，F 0=400N ；取凸轮的推程运动角和回程运动角均为180，远休止角和近休止角均为0，从动件的最大升程h =6. 2mm 。

( 1 - cos ϕ ) 根据已知条件，可以确定从动件的位移方程 y =将上式代入（式3）可得：

h

2

h h 1dx [M ω2-(k +k h )]cos ϕ+c ωsin ϕ+(k +k h ) h +F 0=c ω+kx （式4） 222d ϕ

由于（式4）对应的常微分方程难以求出解析解，这里利用MATLAB 求解出其数值解并与位移方程比较如下图：

-3

升程/m

凸轮的转角/rad

若需要考虑机构的弹性变形，凸轮的轮廓应按照上图中的红色曲线进行设计。

（二）由已知廓线求解从动件的真实运动规律

由于系统的弹性变形，从动件的真实运动y 与等效凸轮升程x 不再相等，当然，从动件的真实运动速度、加速度与对应的理论值也不相等。

由于弹簧的预紧力F 0为常数，它只影响系统振动的初始平衡位置，故在分析从动件的运动规律时不再考虑，从而（式1）被简化为：

d 2y dy dx

M 2+c +(k +k h ) y =c +kx （式5）

dt dt dt

根据振动理论，系统自由振动的固有频率ωn =

k +k h

(1-ξ2) M

其中阻尼比

ξ=

c

2M (k +k h )

代入相关数值计算可得

ωn =6147rad /s ，ξ=0. 0535

如果从动件按照前述的余弦加速度规律运动，则 x =将上式代入（式5）可得：

h

(1-cos ϕ) 2

d 2y dy dx M ω+c ω+(k +k ) y =c ω+kx （式6）

h

d ϕ2d ϕd ϕ

2

利用MATLAB 求解（式6）的数值解，结果如下（图中的n =

-3

ωn

）： ω

位移/m

凸轮的转角/rad

-3

速度/（m /s ）

凸轮的转角/rad

-3

加速度/（m /s 2）

凸轮的转角/rad

从上图中的位移、速度、加速度分析看，当考虑到系统的弹性和阻尼后，工作端的运动规律发生改变，x ≠y 。只有当n 很大时，也就是说当系统的固有频率ωn 很大时（刚度大），而

且凸轮的角速度很小时，y 才接近x 值。

程序附录：

% filename: tulun.m

% function: cooperate with jisuanlunkuo.m to calculate the curve of the cam

function dx=tulun(a,x) %% 凸轮机构的结构参数

h=6.2*10^(-3); % 升程

m=0.08459; % 凸轮机构的当量质量 n=1200;w=n/60*2*pi; % 凸轮转速及角速度 c=55.7148; % 阻尼 k1=10400; % 回位弹簧刚度 k=3194800; % 系统等效弹簧刚度 F=400; % 回位弹簧预紧力

%% 理论轮廓的微分方程

dx=(h/2*(m*w^2-k-k1)*cos(a)+h/2*c*w*sin(a)+0.5*(k+k1)*h+F-k*x)/(c*w);

% filename: jisuanlunkuo.m

% function: cooperate with tulun.m to calculate the curve of the cam %% 计算考虑弹性变形的凸轮的理论轮廓曲线 [a,x]=ode45('tulun' ,[0 2*pi],0); plot(a,x,'r' ); hold on ;

%% 作出未考虑弹性变形的凸轮的理论轮廓曲线 a1=linspace(0,2*pi); h=6.2*10^(-3); y=h/2*(1-cos(a1)); plot(a1,y,'-b' ); %% 添加标注 grid on ;

legend(' 考虑弹性变形后凸轮的理论轮廓' , ' 未考虑弹性变形的凸轮理论轮廓' ); xlabel(' 凸轮的转角/rad'); ylabel(' 升程/m');

% filename: yundongfenxi.m

% function: cooperate with fenxi.m to calculate the displacement,velocity % and acceleration of the cam mechanism when the elastic deformation of % the cam is not considered in design function dy=yundongfenxi(a,y) %% 凸轮机构的结构参数 h=6.2*10^(-3); m=0.08459;

c=55.7148; k1=10400; w=1537;

% 不同情况下将凸轮的角速度值分别设为1537,2049,3074 k=3194800; %% 运动分析微分方程 dy=zeros(2,1); dy(1)=y(2);

dy(2)=(c*w*h/2*sin(a)+k*h/2*(1-cos(a))-(k+k1)*y(1)-c*w*y(1))/(m*w*w);

% filename: fenxi.m

% function: cooperate with yundongfenxi.m to calculate the displacement,velocity

% and acceleration of the cam mechanism when the elastic deformation of % the cam is not considered in design %% 绘制从动件的位移曲线

[a,y]=ode45('yundongfenxi' ,[0 2*pi],[0 0]); n=length(a); figure(1);

plot(a,y(:,1),'--b' );

% 不同情况下将线型分别设为--b ，--g ，--k hold on ;

%% 绘制从动件的速度曲线 figure(2); a1=zeros(n-1,1); for i=1:(n-1) a1(i)=a(i); end

y1=diff(y(:,1))./diff(a); plot(a1,y1,'--b' );

% 不同情况下将线型分别设为--b ，--g ，--k hold on ;

%% 绘制从动件的加速度曲线 figure(3); a2=zeros(n-2,1); for i=2:(n-1) a2(i-1)=a1(i); end

y2=diff(y1(:,1))./diff(a1); plot(a2,y2,'--b' );

% 不同情况下将线型分别设为--b ，--g ，--k hold on ;

% fliename: zuihou.m

% function: give the ideal displacement,velocity and acceleration of the % cam mechanism %% 绘制理想位移曲线 a0=linspace(0,2*pi); n=length(a0);

y0=3.1*10^(-3)*(1-cos(a0)); figure(1); plot(a0,y0,'r' ); hold on ;

%% 绘制理想速度曲线 figure(2); a01=zeros(1,n-1); for i=1:(n-1) a01(i)=a0(i); end

y01=diff(y0)./diff(a0); plot(a01,y01,'r' ); hold on ;

%% 绘制理想加速度曲线 figure(3); a02=zeros(1,n-2); for i=2:(n-1) a02(i-1)=a01(i); end

y02=diff(y01)./diff(a01); plot(a02,y02,'r' ); hold on ; %% 添加标注 figure(1);

xlabel(' 凸轮的转角/rad'); ylabel(' 位移/m');

legend('n=4', 'n=3', 'n=2', ' 标准输出' ); figure(2);

xlabel(' 凸轮的转角/rad'); ylabel(' 速度/（m/s）' );

legend('n=4', 'n=3', 'n=2', ' 标准输出' ); figure(3);

xlabel(' 凸轮的转角/rad'); ylabel(' 加速度/（m/s^2）' );

legend('n=4', 'n=3', 'n=2', ' 标准输出' );

程序运行说明：在MATLAB 主窗口中运行jisuanlunkuo.m 文件可以得到考虑弹性变形后凸轮的理论轮廓曲线和未考虑弹性变形的凸轮理论轮廓曲线对比图；初次运行fenxi.m 可以得到

凸轮角速度

ω=1537rad /s

情况下从动件的位移、速度和加速度图，之后需要按程序说明

中的位置将yundongfenxi.m 文件中的角速度分别改为

2049rad /s

和

3074rad /s

并相应地

将fenxi.m 文件中的线型改为--g 和--k ，然后再次运行fenxi.m 文件可以得到以及

ω=2049rad /s

ω=3074rad /s

情况下从动件的位移、速度和加速度图（每次修改后注意保存），最后

运行zuihou.m 文件得到最终的从动件运动量的对比分析图。