ガスが流出する圧力容器を囲う箱の中の圧力上昇の計算

まだうまくいっていない

clear; close all; clc;

R = 287; %Air. J/kg·K.
rho = 1.293;%Air. f 1.293 kg m−3

P1_init = 8e6;%MPa
pa = 0.1013e6;%MPa
P2_init = pa;

dt = 0.01;%second
t_max = 10;%second
ts = 0:dt:t_max;%second

LWP = 85.73;%mm
SWP = 40.64;%mm
Width = 3.18;%mm
total_area = 0.5*(LWP+Width)*(SWP+Width);
void_area = 0.5*LWP*SWP;
metal_part = total_area - void_area;
phi = metal_part / total_area; %about 0.1

Cd = 1.0;

D = 0.092;%m
V1 = pi*D^2/4 * 0.725;%m^3
V2 = 1.22*1.22*1.6;%m^3
T = 300;%K

Atotal = 1.22*1.22;%m^2

P1 = P1_init;
P2 = P2_init;
P1s = [];
P2s = [];
for idx=1:length(ts)
      
    dp1dt = - phi * Atotal * Cd * P1 / V1 * sqrt(2*(P1-P2)/R/T);
    dp2dt = rho * phi * Atotal * Cd / V2 * sqrt(2*(P1-P2)/rho);

    P1 = P1 + dp1dt*dt;
    P2 = P2 + dp2dt*dt;

    P1s(end+1) = P1;
    P2s(end+1) = P2;
end

plot(ts,P1s./pa);
figure();
plot(ts,P2s./pa);


big;