Can someone please review my code, seems to be an issue in the for loop? Check Volw array (plus a few other), its value is not changing after 2nd instance.

% MATLAB Simulation to predict the 'Thrust' and 'Free-Flight' behaviour of
% the water rocket
clear
clc
%**********************Input section**************************************
%
%---- Inital conditions-------------
%
Volw0=0.0001;                       % Initial volume of water (m^3)
m=0.150;                            % Initial mass of empty rocket (kg)
theta_i=45*pi/180;                  % Launch angle (rad)
P0=344738;                          % Initial gauge pressure of bottle (Pa)
Volb=0.0006;                        % Bottle capacity (m^3)
Lt=0.165;                           % Length of tube (m)
r=0.001;                            % radius of tube (m)
Lr=0.288;                           % rocket length (m)
D=0.074;                            % Rocket diameter (m)
Vola0=Volb-Volw0;                   % Initial volume of air in bottle (m^3)
%
%------Constants-------------------
%
Patm=101300;                        % Atmospheric pressure (Pa)
Dw=998;                             % density of water (kg/m^3)
gamma=1.4;                          % adiabatic index for air
Rn=0.011;                           % nozzle radius (m)
Da=1.19;                            % density of air (kg/m^3)
%
%--------Other calculations-----------------
msys0=m+(Volw0*Dw);                 % Initial mass of fuel+rocket (kg)
An=(pi)*Rn^2;                       % Area of nozzle (m^2)
%-----Measured constants-------------
% 
g=9.59;                             % local gravity (m/s^2)             
C=0.0004;                           % drag constant (kg/m)
%
%================== calculation section===========================
%
% ------ Tube thrust-------
%
% Based on eq 12 in the lab guide, the velocity of the rocket at the 
% end of the tube is:
%
v_tube=sqrt(P0*(pi)*Rn^2*(1-exp(-2*C*Lt/msys0))/C);   %eq 12
%
% *****************************
% ------Water thrust-----
%
% using "v_tube" calculated above as an inital condition, loop over the
% "Phase 2-water thrust" equations
%

n=5000;
dt=0.001;                          %smaller interval for higher accuracy
t=zeros(n,1);
a=zeros(n,1);
P=zeros(n,1);
v_w=zeros(n,1);
Volw=zeros(n,1);
v_sys=zeros(n,1);
msys=zeros(n,1);
mdotw=zeros(n,1);
s=zeros(n,1);
v_sys(1)=v_tube;                     %Initial velocity of rocket
Volw(1)=Volw0;                       %Initial volume of water in the rocket
s(1)=Lt;                             %Rectilinear position at the end of the tube
%---------water thrust phase----------------

for i=1:n-1
    
    P(i)=P0*(Vola0/(Vola0+(Volw0-Volw(i))))^gamma;     %Instantaneous Pressure inside rocket -eq 19 (Pa)
    v_w(i)=-sqrt(2*(P(i)-Patm)/Dw);                    %velocity of water wrt system -eq 23 (m/s)
    mdotw(i)=An*Dw*(v_w(i))*dt;                        %mass flowrate out -eq 24 (kg/s)
    msys(i)=msys0+mdotw(i);                            %rocket mass (kg)
    a(i)=(C*v_sys(i)^2+mdotw(i)*v_w(i))/msys(i);       %instantaneous acceleration of rocket eq 18   
%----  update for the next loop (i+1)--------
    
    v_sys(i+1)=v_sys(i)+a(i)*dt;                       %new velocity of rocket
    Volw(i+1)=Volw0+(mdotw(i)*dt/Dw);                  %new Volume of water
    
%----results---------  
    s(i+1)=s(i)+v_sys(i)*dt+(a(i)*dt^2)/2;             %new position of rocket
        
    t(i+1)=t(i)+dt;
    if Volw(i+1) <=0
        break
    end

end
%

plot(t,s);

xlabel('Time (s)')

ylabel('Distance (m)')

Answer 1

[Monisha Nalluru]

The reason you are not able see the change in values in due to precision.

The differnce between the values of P(i) and P(i+1) is very small, which used in calcuatiing mdotw and Volw. So, this is reason we cannot see the change in values of Volw

Inorder to see the differnce in the array values, increase the precision, by default MATLAB uses 16 digits.

Refer the following documentation

Increase Precision of Numeric Calculation

Using digits function

Set command window display format

Using eps for floating point