11/22/2006

10.7 多個桶子水位問題

多桶子並列成一排,底下以旁通管接通,若其中一個水桶之水位為Ui,則其與相鄰之水桶流水量V應為左右兩流量VL、VR之和,即:


  VL(i)=[U(i-1)-U(i)]*C*DT
  VR(i)=[U(i+1)-U(i)]*C*DT
V(i)=VL(i)+VR(i)
=[U(i-1)-2*U(i)+U(i+1)]*C*DT/A(i)

上述所討論為通式,但實際上前後兩側之水桶均僅有一個旁通管,故要特別處理,即:

V(1)=[U(2)-U(1)]*C*DT/A(1)
V(end)=[U(end-1)-U(end)]*C*DT/A(1)

經過DT的時間後,各桶內水位可依下式調整之:

U(i)=U(i)+V(i)


程式結構


程式結構仍然與前述之二桶水位相同,只是現在U與D均改為向量,而且其大小應相同。本程式增加一個最終時間T,使程式超過此時間時可以停止運算。

同樣,此程式若無參數,則可使用預設資料,若打入pails(1),則會逐項輸入,若直接餵入參數,則會依設定之參數執行,但需注意其大小應一致。

執行結果,各桶之水位會印出,而且會繪出各桶之水位變化圖。

程式內容



function [ff]=pails(U,D,C,delt,T)
% pails.m
% program to calculate waterflows between pails.
% U:water level of pails, cm
% D:diameter of pails, cm
% T:time elapsed through the experiment.
% C:flow coefficent, in cm^3/cm of water level difference.
% Designed by D.S. Fon, Bime, NTU, Date:Nov. 26,2006
while 1
if nargin==0,
UU=[20 0 30 80 35];DD=[5 5 5 5 5];C=1;delt=1;T=100;
elseif U==1,
UUx=[20 0 30 80 35];DDx=[5 5 5 5 5];Cx=1;deltx=1;Tx=100;
disp('Enter following values in vector==> ')
U1=input('Initial water level of pails,cm==> ','s');
if isempty(U1), UU=UUx; else UU=str2num(U1);end
D1=input('Diameter of pails,cm==> ','s');
if isempty(D1), DD=DDx; else DD=str2num(D1);end
T=input('Input the total elapsed time, seconds==> ');
if isempty(T), T=Tx;end
delt=input('Input a time interval, seconds==> ');
if isempty(delt), delt=deltx;end
C=input('Input flow coefficient, cm^3/cm ==>');
if isempty(C), C=Cx;end
else
UU=U;DD=D;
end
AA=pi*DD.^2/4;
npail=length(UU);
tt(1)=0;i=0;nn=ceil(T/delt);ff=zeros(nn,npail+1);
clf;
ff(1,:)=[0, UU];
for i=2:nn;
tt(i)=delt*(i-1);
for j=2:npail-1
dU(j)=(UU(j-1)-2*UU(j)+UU(j+1))*C*delt/AA(j);
end
dU(1)=(UU(2)-UU(1))*C*delt/AA(1);
dU(npail)=(UU(npail-1)-UU(npail))*C*delt/AA(npail);
UU=UU+dU;
ff(i,:)=[tt(i) UU];
line(tt(i),UU,'Marker','.');
end
m=5;if m>npail, m=1;end
for i=2:npail+1
text(tt(m),ff(m,i),['P-',num2str(i)]);
end
xlabel('Time, seconds');
ylabel('Water level at each pail');
grid on;figure(gcf);
yn=input('Continue?(Y/N)[Y]==> ','s');
if ~isempty(yn) & upper(yn)~='Y', break; end
end



執行結果:

>> pails

ans =

0 20.0000 0 30.0000 80.0000 35.0000
1.0000 18.9814 2.5465 31.0186 75.1617 37.2918
2.0000 18.1444 4.8336 31.8167 70.9848 39.2205
3.0000 17.4665 6.8857 32.4373 67.3722 40.8383
4.0000 16.9276 8.7259 32.9152 64.2417 42.1896
5.0000 16.5099 10.3756 33.2787 61.5231 43.3127
6.0000 16.1975 11.8544 33.5507 59.1572 44.2402
7.0000 15.9763 13.1806 33.7499 57.0933 44.9999
8.0000 15.8339 14.3706 33.8911 55.2886 45.6158
----------------------------------------------------------
95.0000 29.3899 30.7720 33.0040 35.2304 36.6036
96.0000 29.4603 30.8153 33.0037 35.1870 36.5337
97.0000 29.5293 30.8578 33.0034 35.1444 36.4651
98.0000 29.5970 30.8994 33.0032 35.1026 36.3978
99.0000 29.6633 30.9402 33.0030 35.0616 36.3319