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

feng137937

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. 
   
3. m2  = 0.002075;                  %模态参数
   
4. m3  = 0.001728;
   
5. c2  = 0.1847;
   
6. c3  = 0.1220;
   
7. k2  = 3.4e6;
   
8. k3  = 4.4e6;
   
9. fz  = 5e-5;
   
10. wd=8586;
    
11. K1=1;
    
12. K2=-4.6e5;
    
13. ylag1 = Z(:,1);
    
14. phi_2=cos(10*pi*t);              %第二阶模态振型
    
15. phi_3=cos(20*pi*t);              %第三阶模态振型
    
16. if state == +1
    
17.     rh=1;                        % rh=1表示切屑厚度为正值      
    
18. else
    
19.     rh=0;                        % rh=0表示刀齿不切削工件
    
20. end
    
21. yp = [ y(2)
    
22.     (-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
    
23.     y(4)
    
24.     (-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. 
   
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 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. 
   
4. phi_2=cos(10*pi*t);                      %第二阶模态振型     
   
5. phi_3=cos(20*pi*t);                      %第三阶模态振型  
   
6. ylag1 = Z(:,1);
   
7. u=  phi_2*y(1)+ phi_3*y(3);              %当前时间振动位移
   
8. u_tau= phi_2*ylag1(1)+ phi_3*ylag1(3);   %前一刀齿周期振动位移
   
9. fz  = 5e-5;
   
10. value = fz+u_tau-u;                      %切屑厚度
    
11. isterminal = 1;
    
12. direction = 1;      %the event function is increasing 表示未变形切屑厚度由小于0增大到正值

复制代码

利用matlab的dde23求解铣削动力学时滞方程的数值解遇到一些问题。由于振动位移的增大切削厚度会由正值减小到负值从而引起方程变化；当振动位移小到一定程度后，切屑厚度又变成正值，方程又放生变化。我是利用Events中的direction来判断切削厚度引起的方程的形式。direction=-1表示切削厚度由正值减小到负值，direction=1表示切削厚度变为正值。但是计算结果不准确。我这样利用direction来确定方程的形式是否存在问题，希望高手能给一些指导。
@牛小贱 @xiezhh @ChaChing