My code appears to be skipping iterations?

%Input initial conditions
Pi=300000; %Air pressure applied
Po=101325; %Pressure discharged
g=9.81; %Acceleration due to gravity
RHO=1400; %density
Dt=0.05; %Diameter of throat
At=pi*(Dt.^2)/4 %Area of throat
Df=6.096; %Diameter of feed tank
Af=pi*(Df.^2)/4 %Area of feed tank
Dpc=Df/5; %Diameter of pumping chamber
Apc=pi*(Dpc.^2)/4 %Area of pumping chamber
Cd=0.975; %Coefficient of discharge
hpc=1.71; %Height of pumping chamber
hinitial=1.71; %Initial height in pumping chamber
hfinitial=6.14; %initial height in feed tank
Cp=0.7 %diffusivity coefficient

%All details up to this point must be input by user
%Calculation of further initial conditions
Ptinitial=RHO*g*hfinitial %Initial pressure static head
Qiinitial=Cd*At*sqrt((2*(Pi+(RHO*g*hinitial)-Ptinitial)/RHO))%Initial flowrate into nozzle
Qoinitial=(At/(sqrt(1-Cp)))*sqrt(2*(Pi-Po)/RHO) %Initial flowrate out of RFD
QHinitial=Qoinitial-Qiinitial %Flowrate of feed tank
tp=sqrt(2/g)*(Apc/(Cd*At))*((sqrt(((Pi-Ptinitial)/(RHO*g))+hpc))-sqrt((Pi-Ptinitial)/(RHO*g)))%Pumping time
n=round(tp) %Rounded pumping time
%Assigning an empty matrix for all values to be calculated
RFDoutput=zeros(n,6) %Output matrix
RFDoutput(:,1)=0:1:n-1 %Column one for time

%Assigning initial conditions to first row
RFDoutput(1,1)=0 %Initial time
RFDoutput(1,2)=hfinitial 
RFDoutput(1,3)=hinitial 
RFDoutput(1,4)=Ptinitial 
RFDoutput(1,5)=Qiinitial 
RFDoutput(1,6)=QHinitial

t=1

while hinitial>0
    hdecf(t)=QHinitial*t/Af %Height decrease

      hdecpc(t)=Qiinitial*t/Apc

      hnewH(t)=hfinitial-hdecf(t)
      hfinitial=hnewH(t)
      RFDoutput(t+1,2)=hfinitial

      hnewpc(t)=hinitial-hdecpc(t)
      hinitial=hnewpc(t)
      RFDoutput(t+1,3)=hinitial
      % CODE IS NOT WORKING FOR THE HEIGHT AND Qi

      Pt(t)=RHO*g*hnewH(t)
      RFDoutput(t+1,4)=Pt(t)

      Qi(t)=Cd*At*sqrt(2*((Pi+(RHO*g*hnewpc(t))-Pt(t))/RHO))
      RFDoutput(t+1,5)=Qi(t)
      Qiinitial=Qi(t)

      QH(t)=Qoinitial-Qi(t)
      RFDoutput(t+1,6)=QH(t)
      QHinitial=QH(t)
      t=t+1
  end

if hfinitial>Dpc
    h = warndlg('Warning: H empty when h=0')
    %Warning message for empty H
end

%Begin plotting various parameters
x=RFDoutput(:,1)
y1=RFDoutput(:,2)
subplot(2,2,1)
plot(x,y1)
title('L level in H as a function of time')
xlabel('Time increment (s)')
ylabel('H level (m)')

y2=RFDoutput(:,3)
subplot(2,2,2)
plot(x,y2)
title('H level in feed tank as a function of time')
xlabel('Time increment (s)')
ylabel('Feed tank level (m)')

y3=RFDoutput(:,4)
subplot(2,2,3)
plot(x,y3)
title('Static pressure at nozzle as a function of time')
xlabel('Time increment (s)')
ylabel('Nozzle Pressure (Pa)')

y4=RFDoutput(:,5)
subplot(2,2,4)
plot(x,y4)
title('Flowrate through nozzle as a function of time')
xlabel('Time increment (s)')
ylabel('Nozzle Flowrate, Qi (m^3 /s)')

Just to explain what I'm trying to do. My while loop is supposed to use an initial flowrate (Q) and at a time increment of 1 second, calculate the decrease in height in the level in a tank. For each result, I want it to put various parameters into the RFDoutput matrix. It seems to be skipping loops/finishing early? It should be filling the full matrix for RFDoutput, so it should have the same number of rows as the pumping time (tp). I want it to terminate when the height in the pumping chamber (hpc) becomes zero. Thanks in advance.

Answer 1

[Walter Roberson]

Your hdecpc is not commented.

You have

hdecpc(t)=Qiinitial*t/Apc

which is something that increases with time.

You have

        hnewpc(t)=hinitial-hdecpc(t)
        hinitial=hnewpc(t)

so your height is decreasing by this value whose magnitude increases with time. Effectively you have an accelerating rate of decrease in height. I suspect you want a linear increase (or perhaps a decrease that varies with pressure of overlaying water.)