2020年iypt第六题下沉的碗（二）仿真动画代码
代码允许转载，但是放在附录时一定要记得把出处和制作人标清楚！！！！！！！！！！！！！！！！！！！！！！！！！！！
常见问题：
1.代码复制一定要复制到编辑器中，而不是命令行
2.一定保证符号全英文
3.对代码修改后尤其是方程部分修改后，求解错误是正常现象，请一定要首先检查方程是否正确
4.这个代码求解方程是龙格库塔法求解的方程，如果不了解数值计算的可以先了解一下数值计算！！
软件：matlab
代码：
%--------------
%2020iypt下沉的碗仿真动画
%使用ode45求解动力学方程后进行动画模拟
%这里的模型使用的是无限长的碗，图画的为有限长所以会有没过边界后还在下降情况
%作者：多一份问候 时间：2020/1/22 仅在b战发布
%引用或转载请标明出处
%----------------
global M g p R r hd N N2 L
M=0.13581;%千克
g=9.8;
p=1000;%水的密度
R=0.05;
r=0.003;
hd=0.002;%碗的厚度
N=0.0001;%阻尼系数小孔的
N2=1;%碗的阻力系数
L=0.1;%碗的高度
t1=[0,14];
H1=[0;0;0;0];
[t,H] = ode45(@dH,t1,H1);%求解动力学方程
% hold on
% plot(t,H(:,1:3));%绘制的曲线随时间的变化图
% plot(t,H(:,1));%碗的底端高度
% plot(t,H(:,2));%碗的下落速度
% plot(t,H(:,3));%碗内液面高度
%动画绘制
%在上述运行完毕后得到
A=size(H);
loops = A(1,1);
F(loops) = struct('cdata',[],'colormap',[]);
v = VideoWriter('peaks135.31kg0.5_0.003.avi');%创建动画
open(v);
for j=1:loops
Hz{j}=num2str(H(j));
end
for i=1:A(1,1)
figure(i)
hold on
axis([-0.25 0.25 -0.4 0.1])%设置坐标轴范围x：-0.25，0.25 y：-0.4 0.1
%碗的图形绘制
%这里的图形绘制认为是随着页面的高度变画而变化的
plot([R,R],[L-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%右侧的碗
plot([-R,-R],[L-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%左侧的碗
plot([-R,-r],[-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%碗底部的小孔左侧
plot([r,R],[-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%碗底部的小孔右侧
plot([-R,R],[-H(i,1)+H(i,3),-H(i,1)+H(i,3)],'LineWidth',2,'Color','b');%碗内液面
plot([-0.25,R],[-H(i,1)+H(i,3),-H(i,1)+H(i,3)],'LineWidth',0.5,'Color','b','LineStyle','--');%碗内液面指示线
text(-0.25,-H(i,1)+H(i,3)-0.006,Hz{1,i});%碗内液面指示线数字
plot([-0.25,R],[L-H(i,1),0.1-H(i,1)],'LineWidth',0.5,'Color','b','LineStyle','--');%碗内液面指示线
text(-0.25,L-H(i,1)-0.006,Hz{1,i});%碗内液面指示线数字
%水面高度认为是不变的
plot([-0.25,-R],[0,0],'LineWidth',2,'Color','c');%液面高度左侧
plot([R,0.25],[0,0],'LineWidth',2,'Color','c');%页面高度右侧
F(i) = getframe(gcf);%保存帧
writeVideo(v,F(i));%写入动画
close
end
close(v);
function dHdt=dH(t,H)
global M g p R r hd N N2
dHdt=zeros(4,1);
dHdt(1)=H(2);
dHdt(2)=g-p*g*pi*(H(1)-H(3))*(R^2-r^2)/M-N2*H(2);
dHdt(3)=sqrt(-2*g*r^4/R^4*(H(3)-H(1))-N*H(4)*hd);
dHdt(4)=H(3);
end
