Я пытаюсь смоделировать простую ракету-бутылку с помощью ode45, но масса воздуха постоянно увеличивается, и я не могу понять, почему.Мой код всегда передает отрицательное «изменение массы воздуха» (dmair) обратно в ode45, но на следующем временном шаге моя масса будет иногда входить в функцию выше, иногда в 3-4 раза выше.Кроме того, масса должна только уменьшаться, но моя функция - создать несколько критических точек вместо одной плавной кривой.Я застрял в этой проблеме в течение нескольких дней, поэтому любая оценка будет оценена!
Вот график массы воздуха и скорости его изменения.Очевидно, что масса воздуха всегда должна уменьшаться со временем, но я в недоумении от того, как оде45 ведет себя так дико непоследовательно.Кто-нибудь может мне помочь?
Вот моя основная функция:
%% Set up constants
global p_i
global v_i
global t_i
global v_bottle
global water_density
global p_atmosphere
global area_throat
global area_bottle
global discharge_coeff
global gas_const
global mass_rocket
global mass_air
global p_end
global t_end
global g
global air_density
global drag_coeff
global vel_i
global theta_i
global x_i
global z_i
global test_stand
global Thrust_vector
global water_exhausted
water_exhausted = 0;
t_span = [0 5];
v_bottle = 0.002; %m^3,
p_i = convpres(12.1+50,'psi','Pa');
v_i = 0.001; %Total in m^3
t_i = 300; %kelvin
water_density = 998.2; %kg/m^3 from wolframalpha
p_atmosphere = convpres(12.5,'psi','Pa'); %Pascals,
area_throat = 0.000346361; %m^2, from an estimate of d = 2.1cm
area_bottle = 0.008659015; %m&2, from d = 10.5cm
gas_const = 287; %J/KgK
p_end = p_i * ((v_i/v_bottle)^1.4); %pressure when all the water has been ejected
t_end = t_i * ((v_i/v_bottle)^0.4); %temperature when all the water has been ejected
g = 9.81; %m/s^2
air_density = 0.961; %kg/m^3
drag_coeff = 0.5;
discharge_coeff = 0.8;
vel_i = 0; %m/s
theta_i = pi/4; %radians (45 degrees)
x_i = 0; %m
z_i = 0.25; %m above the ground
test_stand = 0.5; %m
%determine initial mass of the system
mass_bottle = 0.15; %kg
mass_water = water_density * (v_bottle - v_i);
mass_air = v_i * (p_i/(gas_const * t_i));
mass_rocket = mass_bottle + mass_water + mass_air;
rocket_parameters = [v_i; mass_rocket; mass_air; vel_i; theta_i; x_i; z_i;];
%% Call ODE45
[t, rocket_values] = `ode45(@rocket_ode_rework,t_span,rocket_parameters);`
А вот функция ode45 вызывает:
function time_values = rocket_ode_rework(t,rocket_par)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v_air = rocket_par(1,1);
m_rocket = rocket_par(2,1);
m_air = rocket_par(3,1);
velocity = rocket_par(4,1);
theta = rocket_par(5,1);
x = rocket_par(6,1);
z = rocket_par(7,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global p_i
global v_i
global t_i
global v_bottle
global water_density
global p_atmosphere
global area_throat
global area_bottle
global discharge_coeff
global gas_const
global mass_rocket
global mass_air
global p_end
global t_end
global g
global air_density
global drag_coeff
global vel_i
global theta_i
global x_i
global z_i
global test_stand
global Thrust_vector
global water_exhausted
gam = 1.4;
if water_exhausted == 1
v_air = v_bottle;
end
% Water Propellant Phase
if v_air < v_bottle && water_exhausted == 0
bottle_pressure = p_i*((v_i/v_air)^gam);
dvolume = discharge_coeff * area_throat * sqrt((2 / water_density) * (bottle_pressure - p_atmosphere));
dmrocket = -discharge_coeff * area_throat * sqrt(2 * water_density * (bottle_pressure - p_atmosphere));
T = 2 * discharge_coeff * area_throat * (bottle_pressure - p_atmosphere);
dmair = 0;
% Air Propellant Phase
elseif v_air >= v_bottle || water_exhausted == 1
water_exhausted = 1;
bottle_pressure = p_end * (m_air / mass_air)^gam;
bottle_density = m_air / v_bottle;
bottle_temp = bottle_pressure / (gas_const * bottle_density);
if bottle_pressure > p_atmosphere
critical_pressure = bottle_pressure * (2/(gam + 1))^(gam / (gam - 1));
if critical_pressure > p_atmosphere %choked flow
p_exit = critical_pressure;
t_exit = bottle_temp * (2 / (gam + 1));
vel_exit = sqrt(gam * gas_const * t_exit);
density_exit = (critical_pressure / (gas_const * t_exit));
else %unchoked flow
p_exit = p_atmosphere;
Mach = sqrt((-1 + (bottle_pressure / p_atmosphere)^((gam-1)/gam))*(2/(gam-1)));
t_exit = bottle_temp * (1 + ((gam-1)/2) * Mach^2);
density_exit = p_atmosphere / (gas_const * t_exit);
vel_exit = Mach * sqrt(gam * gas_const * t_exit);
end
dmair = -discharge_coeff * density_exit * area_throat * vel_exit;
dmrocket = dmair;
T = -dmair * vel_exit + (p_exit - p_atmosphere) * area_throat;
dvolume = 0; %for each non-water propellant phase, there should be no change in volume
% Ballistic Phase
else
dmrocket = 0;
dmair = 0;
T = 0;
dvolume = 0;
end
end
D = (0.5 * air_density * (velocity^2)) * drag_coeff * area_bottle;
dvel = (T - D - (m_rocket * g * sin(theta)))/m_rocket;
if velocity <= test_stand %rocket has not left the test stand
dtheta = 0;
else
dtheta = (-g * cos(theta)) / velocity;
end
dx = velocity * cos(theta);
dz = velocity * sin(theta);
Thrust_vector = [Thrust_vector; T];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
time_values(1,1) = dvolume; %volume integral
time_values(2,1) = dmrocket; %mass integral
time_values(3,1) = dmair; %air integral
time_values(4,1) = dvel;
time_values(5,1) = dtheta;
time_values(6,1) = dx;
time_values(7,1) = dz;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
Спасибо!