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.

3 views (last 30 days)
% 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)')

Answers (1)

Monisha Nalluru
Monisha Nalluru on 28 Oct 2020
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

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!