声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2652|回复: 1

[编程技巧] 求教matlab dde23计算铣削动力学时滞微分方程的数值解

[复制链接]
发表于 2019-3-29 16:21 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 feng137937 于 2019-3-29 16:33 编辑
  1. state = +1;
  2. opts1 = ddeset('Jumps',0,'Events',@events_decrease,'RelTol',1e-3,'AbsTol',1e-6);
  3. opts2 = ddeset('Jumps',0,'Events',@events_increase,'RelTol',1e-3,'AbsTol',1e-6);
  4. 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]
  5. direction=-1;  %表示未变形切屑厚度由正值减小到0
  6. while sol.x(end) < 0.1
  7.     fprintf('Restart at %5.1f.\n',sol.x(end));
  8.     if direction==-1
  9.         state = - state;
  10.         sol = dde23(@M_DDE,7.3e-4,sol,[sol.x(end),0.1],opts2,state);
  11.         direction=1;
  12.     end
  13.     if direction==1 %表示未变形切屑厚度由小于0增大到正值
  14.         state = - state;
  15.         sol = dde23(@M_DDE,7.3e-4,sol,[sol.x(end),0.1],opts1,state);
  16.         direction=-1;
  17.     end
  18. end
  19. yplot = sol.y;
  20. figure
  21. plot(sol.x,yplot(1,:))
  22. figure
  23. plot(sol.x,yplot(3,:))
复制代码
  1. function yp = M_DDE(t,y,Z,state)

  2. m2  = 0.002075;                  %模态参数
  3. m3  = 0.001728;
  4. c2  = 0.1847;
  5. c3  = 0.1220;
  6. k2  = 3.4e6;
  7. k3  = 4.4e6;
  8. fz  = 5e-5;
  9. wd=8586;
  10. K1=1;
  11. K2=-4.6e5;
  12. ylag1 = Z(:,1);
  13. phi_2=cos(10*pi*t);              %第二阶模态振型
  14. phi_3=cos(20*pi*t);              %第三阶模态振型
  15. if state == +1
  16.     rh=1;                        % rh=1表示切屑厚度为正值      
  17. else
  18.     rh=0;                        % rh=0表示刀齿不切削工件
  19. end
  20. yp = [ y(2)
  21.     (-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
  22.     y(4)
  23.     (-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];
复制代码
  1. function [value,isterminal,direction] = events_decrease(t,y,Z,state)  %state的作用

  2. phi_2=cos(10*pi*t);                                   %第二阶模态振型
  3. phi_3=cos(20*pi*t);                                   %第三阶模态振型
  4. ylag1 = Z(:,1);
  5. u=  phi_2*y(1)+ phi_3*y(3);                           %当前时间振动位移
  6. u_tau= phi_2*ylag1(1)+ phi_3*ylag1(3);                %前一刀齿周期振动位移
  7. fz  = 5e-5;
  8. value = fz+u_tau-u;                                   %切屑厚度   
  9. isterminal = 1;
  10. direction = -1;      %the event function is decreasing 表示未变形切屑厚度由正值减小到0
复制代码
  1. function [value,isterminal,direction] = events_increase(t,y,Z,state)  %state的作用
  2. %EXER7E  The event function for Exercise 7 of the DDE Tutorial.

  3. phi_2=cos(10*pi*t);                      %第二阶模态振型     
  4. phi_3=cos(20*pi*t);                      %第三阶模态振型  
  5. ylag1 = Z(:,1);
  6. u=  phi_2*y(1)+ phi_3*y(3);              %当前时间振动位移
  7. u_tau= phi_2*ylag1(1)+ phi_3*ylag1(3);   %前一刀齿周期振动位移
  8. fz  = 5e-5;
  9. value = fz+u_tau-u;                      %切屑厚度
  10. isterminal = 1;
  11. direction = 1;      %the event function is increasing 表示未变形切屑厚度由小于0增大到正值
复制代码
利用matlab的dde23求解铣削动力学时滞方程的数值解遇到一些问题。由于振动位移的增大切削厚度会由正值减小到负值从而引起方程变化;当振动位移小到一定程度后,切屑厚度又变成正值,方程又放生变化。我是利用Events中的direction来判断切削厚度引起的方程的形式。direction=-1表示切削厚度由正值减小到负值,direction=1表示切削厚度变为正值。但是计算结果不准确。我这样利用direction来确定方程的形式是否存在问题,希望高手能给一些指导。
@牛小贱 @xiezhh @ChaChing

铣削动力学方程

铣削动力学方程
回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-28 09:09 , Processed in 0.071149 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表