|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 feng137937 于 2019-3-29 16:33 编辑
- state = +1;
- opts1 = ddeset('Jumps',0,'Events',@events_decrease,'RelTol',1e-3,'AbsTol',1e-6);
- opts2 = ddeset('Jumps',0,'Events',@events_increase,'RelTol',1e-3,'AbsTol',1e-6);
- sol = dde23(@M_DDE,7.3e-4,[0; 0; 0;0],[0,0.1],opts1,state);%sol = dde23(ddefun,lags,history,tspan) 时间区间应该很大为了缩短计算时间设置为[0,0.1]
- direction=-1; %表示未变形切屑厚度由正值减小到0
- while sol.x(end) < 0.1
- fprintf('Restart at %5.1f.\n',sol.x(end));
- if direction==-1
- state = - state;
- sol = dde23(@M_DDE,7.3e-4,sol,[sol.x(end),0.1],opts2,state);
- direction=1;
- end
- if direction==1 %表示未变形切屑厚度由小于0增大到正值
- state = - state;
- sol = dde23(@M_DDE,7.3e-4,sol,[sol.x(end),0.1],opts1,state);
- direction=-1;
- end
- end
- yplot = sol.y;
- figure
- plot(sol.x,yplot(1,:))
- figure
- plot(sol.x,yplot(3,:))
复制代码- function yp = M_DDE(t,y,Z,state)
- m2 = 0.002075; %模态参数
- m3 = 0.001728;
- c2 = 0.1847;
- c3 = 0.1220;
- k2 = 3.4e6;
- k3 = 4.4e6;
- fz = 5e-5;
- wd=8586;
- K1=1;
- K2=-4.6e5;
- ylag1 = Z(:,1);
- phi_2=cos(10*pi*t); %第二阶模态振型
- phi_3=cos(20*pi*t); %第三阶模态振型
- if state == +1
- rh=1; % rh=1表示切屑厚度为正值
- else
- rh=0; % rh=0表示刀齿不切削工件
- end
- yp = [ y(2)
- (-c2*y(2)-k2*y(1)+ phi_2*( K1*sin(wd*t)+K2*( fz + phi_2*( ylag1(1)-y(1))+ phi_3*( ylag1(3)-y(3) )))*rh)/m2
- y(4)
- (-c3*y(4)-k3*y(3)+ phi_3*( K1*sin(wd*t)+ K2*( fz + phi_2*( ylag1(1)-y(1))+ phi_3*( ylag1(3)-y(3) )))*rh)/m3];
复制代码- function [value,isterminal,direction] = events_decrease(t,y,Z,state) %state的作用
- phi_2=cos(10*pi*t); %第二阶模态振型
- phi_3=cos(20*pi*t); %第三阶模态振型
- ylag1 = Z(:,1);
- u= phi_2*y(1)+ phi_3*y(3); %当前时间振动位移
- u_tau= phi_2*ylag1(1)+ phi_3*ylag1(3); %前一刀齿周期振动位移
- fz = 5e-5;
- value = fz+u_tau-u; %切屑厚度
- isterminal = 1;
- direction = -1; %the event function is decreasing 表示未变形切屑厚度由正值减小到0
复制代码- function [value,isterminal,direction] = events_increase(t,y,Z,state) %state的作用
- %EXER7E The event function for Exercise 7 of the DDE Tutorial.
- phi_2=cos(10*pi*t); %第二阶模态振型
- phi_3=cos(20*pi*t); %第三阶模态振型
- ylag1 = Z(:,1);
- u= phi_2*y(1)+ phi_3*y(3); %当前时间振动位移
- u_tau= phi_2*ylag1(1)+ phi_3*ylag1(3); %前一刀齿周期振动位移
- fz = 5e-5;
- value = fz+u_tau-u; %切屑厚度
- isterminal = 1;
- direction = 1; %the event function is increasing 表示未变形切屑厚度由小于0增大到正值
复制代码 利用matlab的dde23求解铣削动力学时滞方程的数值解遇到一些问题。由于振动位移的增大切削厚度会由正值减小到负值从而引起方程变化;当振动位移小到一定程度后,切屑厚度又变成正值,方程又放生变化。我是利用Events中的direction来判断切削厚度引起的方程的形式。direction=-1表示切削厚度由正值减小到负值,direction=1表示切削厚度变为正值。但是计算结果不准确。我这样利用direction来确定方程的形式是否存在问题,希望高手能给一些指导。
@牛小贱 @xiezhh @ChaChing |
-
铣削动力学方程
|