撒克逊碗动态物理仿真

%% 撒克逊碗

clc;clear;global M g p R r hd N N2 L

g=9.8;M=0.13581;R=0.05;p=1000;%千克%水的密度

r=0.003;hd=0.002;N=0.0001;%碗的厚度%阻尼系数小孔的

N2=1;L=0.1;%碗的阻力系数%碗的高度

t1=[0,40];H1=[0;0;0;0]

[t,H] = ode45(@dH,t1,H1);%求解动力学方程

%% 动画绘制

figure('numbertitle','off','name','撒克逊碗');% 动画绘制

for i=1:10:length(H)

clf;xlim([-0.25 0.25]);ylim([-0.4 0.1]);hold on%设置坐标轴范围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');%页面高度右侧

hold off;pause(0.05);%可以在这里添加输出动图的程序

drawnow;%实时更新坐标图窗口

end

 

 

% %% 绘制的曲线随时间的变化图

% figure(2);subplot(2 ,2 ,1);plot(t,H(:,1:3));%绘制的曲线随时间的变化图

% subplot(2 ,2 ,3);plot(t,H(:,1));%碗的底端高度

% subplot(2 ,2 ,2);plot(t,H(:,2));%碗的下落速度

% subplot(2 ,2 ,4);plot(t,H(:,3));%碗内液面高度

%% 公式

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