10.7 多個桶子水位問題

多桶子並列成一排，底下以旁通管接通，若其中一個水桶之水位為Ui，則其與相鄰之水桶流水量Ｖ應為左右兩流量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)

程式結構

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

同樣，此程式若無參數，則可使用預設資料，若打入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

執行後之水位變化圖