求教matlab dde23计算铣削动力学时滞微分方程的数值解
本帖最后由 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,,,opts1,state);%sol = dde23(ddefun,lags,history,tspan) 时间区间应该很大为了缩短计算时间设置为
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,,opts2,state);
direction=1;
end
if direction==1 %表示未变形切屑厚度由小于0增大到正值
state = - state;
sol = dde23(@M_DDE,7.3e-4,sol,,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 = 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 表示未变形切屑厚度由正值减小到0function = events_increase(t,y,Z,state)%state的作用
%EXER7EThe 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
页:
[1]