Trouble using ODE45 for coupled non-linear differential equation

Hi,

I am attempting to model the back filling of a large vacuum vessel with Nitrogen. The flow of the Nitrogen is controlled by the valve with a certain valve flow coefficient. P2 is constant and approximately one atmosphere. I have derivived the following equations to model this.

So I wrote the following code with ODE45 to monitor the pressure over time.

% IC Vector
IC_BF = [P01, Q0, dP0, dQ0];

Unrecognized function or variable 'P01'.

t = [0 10];
options = odeset('RelTol',1e-12 );
[t_ode45,Result] = ode45(@BF_dyn, t, IC_BF, options)

function result = BF_dyn(t, x)

% Constants
R = 8.314; % J / mol·K
C_v = 0.28*(6.309e-5/1); % (gallons/min) ( 6.309e-5 (m^3/s)/ 1 (gallons/min)) = m^3/s
Tamb = 293; % K
G = 0.967;
dewar_vol = 4.4; % m^3
rho_N2 = 1.25; %kg/m^3
M_N2 = 28.02*(1/1000); % g/mol * (1 kg/1000g) = kg/mol
P02 = 6894.76; % Pa

b = (sqrt(G)/(C_v))^2;
a = rho_N2*R*Tamb/(M_N2*dewar_vol);
result = [
          -2*x(2)*x(4)*b;
          (1/a)*x(4);
          x(3);
          x(4);
    ];
end

However, when I plot the results I get the following, which goes significantly higher than P2 before attempting to go negative when in reality it should asymptote out to P2.

Am I implement this system into ODE45 wrong? Is there any way to incorporate the realtionship between P1 and P2 into ODE45 by adding additional arguement to my results vector?

Thank you for any assistance you can offer!

[Jacob]

The bottom two equations in result are storing the derivitive of pressure and flowrate as state variables so that I can reference the derivitives in the first two equations. I'm not sure if this was the correct approach to this problem

[Sam Chak]

@Jacob, For what purpose do you put these unused constants inside the BF_dyn(t, x) function?

delta_P0= 6894.76;      % Pa
P01     = 1.33322e-5;   % pa
dP01    = 0;            % Pa
P02     = 6894.76;      % Pa

[Torsten]

If you equate the two equations you posted, dQ/dt cancels out and you get

Q(t) = -R*T*M*rho/V / (2*G/Cv^2)

Does that make sense ?

Answer 1

[Torsten]

I can't recover the four differential equations you solve from the two equations you included in mathematical form.

Your equations can be solved analytically:

syms a b P01 Q0 dP0 dQ0
syms t x1(t) x2(t) x3(t) x4(t) 
odes = [diff(x1,t) == -2*b*x2*x4,diff(x2,t) ==x4/a,diff(x3,t)==x3,diff(x4,t)==x4];
conds = [x1(0)==P01,x2(0)==Q0,x3(0)==dP0,x4(0)==dQ0];
sol = dsolve(odes,conds);
sol.x1
sol.x2
sol.x3
sol.x4

[Walter Roberson]

(Unfortunately due to a glitch, the output is not rendering here on MATLAB Answers. If you run the code on your own display then it will work.)