# 问题一题解

## 分析

（关于求解微分方程这一步骤，我觉得直接近似求解和拟合都可以，但是拟合的话最好不要用线性，二次和指数的效果相对可能好一丢丢。）

## 求解

### 求解微分方程

d P d ρ = E ρ \frac{d P}{d \rho}=\frac{E}{\rho}

d P E = ∫ d ρ ρ + C \frac{d P}{E}=\int \frac{d \rho}{\rho}+C

ρ = 0.80853816 + 0.80853816 0.0006 P 1 + 0.0017 P ; \rho=0.80853816+0.80853816\frac{0.0006P}{1+0.0017P};

### 求解管内燃油密度变化

I ( t ) = { C A 2 [ P 1 ( t ) − P 2 ( t ) ρ 1 , t ∈ ( 0 , t 1 ) 0 , t ∈ ( t 1 , T 1 ) I(t)=\left\{\begin{array}{ll} C A \sqrt{\frac{2\left[P_{1}(t)-P_{2}(t)\right.}{\rho_{1}}} & , t \in\left(0,t_{1}\right) \\ 0 & , t \in\left(t_{1},T_1\right) \end{array}\right.

O ( t ) = { 100 t , t ∈ ( 0 , 0.2 ) 20 , t ∈ ( 0.2 , 2.2 ) − 100 t + 240 , t ∈ ( 2.2 , 2.4 ) 0 , t ∈ ( 2.4 , T 2 ) O(t)=\left\{\begin{array}{ll} 100t & , t \in\left(0, 0.2\right) \\ 20 & , t \in\left(0.2, 2.2\right) \\ -100 t+240 & , t \in\left(2.2, 2.4\right) \\ 0 & , t \in\left(2.4,T_2)\right. \end{array}\right.

m 1 ( t ) = ρ 1 I ( t ) Δ t m_{1}(t)=\rho_{1} I(t) \Delta t

m 2 ( t ) = ρ 2 O ( t ) Δ t m_{2}(t)=\rho_{2} O(t) \Delta t

m 1 ( t ) − m 2 ( t ) = ρ 1 I ( t ) Δ t − ρ 2 O ( t ) Δ t m_{1}(t)-m_{2}(t)=\rho_{1} I(t) \Delta t-\rho_{2} O(t) \Delta t

ρ 2 ( t + Δ t ) − ρ 2 ( t ) = 1 V [ ρ 1 I ( t ) − ρ 2 O ( t ) ] Δ t \rho_{2}(t+\Delta t)-\rho_{2}(t)=\frac{1}{V}\left[\rho_{1} I(t)-\rho_{2} O(t)\right] \Delta t

m i n ∑ i = 1 N ∣ P ( t k ) − P 0 ∣ N min\frac{\sum_{i=1}^{N}\left|P\left(t_{k}\right)-P_{0}\right|}{N}

# 问题二题解

## 求解

### 燃油的进入

{ x = r ⋅ cos ⁡ θ y = r ⋅ sin ⁡ θ \left\{\begin{array}{l} x=r \cdot \cos \theta \\ y=r \cdot \sin \theta \end{array}\right.

{ x m , k = x m , 0 ⋅ cos ⁡ φ k − y m , 0 ⋅ sin ⁡ φ k y m , k = x m , 0 ⋅ sin ⁡ φ k − y m , 0 ⋅ cos ⁡ φ k \left\{\begin{array}{l} x_{m, k}=x_{m, 0} \cdot \cos \varphi_{k}-y_{m, 0} \cdot \sin \varphi_{k} \\ y_{m, k}=x_{m, 0} \cdot \sin \varphi_{k}-y_{m, 0} \cdot \cos \varphi_{k} \end{array}\right.

### 燃油的喷出

m i n ∑ i = 1 N ∣ P ( t k ) − P 0 ∣ N min\frac{\sum_{i=1}^{N}\left|P\left(t_{k}\right)-P_{0}\right|}{N}

# 问题三题解

## 求解

Q d = C A 2 [ P ( t ) − P e ( t ) ] ρ Q_d = CA\sqrt{\frac{2[P(t)-P_e(t)]}{\rho}}

# 代码附录

## Q1.1

%密度函数，压力函数，出油函数
P=@(x)(x-0.80853816)/(0.001859637768-0.00170*x);
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+0.0017*x));
O=@(x)100*x.*(x<=0.2)+20.*(x>0.2&x<=2.2)-(100*x-240).*(x>2.2&x<=2.4)+0.*(x>2.4);

%规定常量
delta_t=0.01;
C=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
p0=100;
m0=density(100)*V;
ml=m0;

y_test=[]
x_test=[]
l=0;

for t=0.28752:0.28752

result=0;
t_out=0;t_in=0;flag=0;
x=[];y=[];   %画图坐标
l=l+1;

for i=0:delta_t:100000

flag=flag+1;
x(flag)=i;
y(flag)=p0;

vin=0;
if t_in>t+10
t_in=t_in-(t+10);
end
if t_in<t
vin=(C*A*(2*(160-p0)/0.8696)^(.5));
end

t_in=t_in+delta_t;
m_in=0.8696*vin*delta_t;
if t_out>100
t_out=t_out-100;
end

vout=O(t_out)*delta_t;
t_out=t_out+delta_t;
mout=vout*density(p0) ;
delta_m=m_in-mout;

rho_now=(m0+delta_m)/V;
p0=P(rho_now);

m0=m0+delta_m;

flag=flag+1;
x(flag)=i;
y(flag)=p0;
end

%搜索判定条件
y_ave=0;k=0;
for j=9950000:10000000
k=k+1;
y_ave=y_ave+y(j);
end
y_ave=y_ave/k;
y_test(l)=y_ave;
x_test(l)=t;
end

plot(x,y)
%plot(x_test,y_test);

## Q1.2

%密度函数，压力函数，出油函数
P=@(x)(x-0.80853816)/(0.001859637768-0.00170*x);
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+0.0017*x));
O=@(x)100*x.*(x<=0.2)+20.*(x>0.2&x<=2.2)-(100*x-240).*(x>2.2&x<=2.4)+0.*(x>2.4);

%规定常量
delta_t=0.01;
C=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
p0=100;
m0=density(100)*V;
ml=m0;

y_test=[]
x_test=[]
l=0;
t=0;

result=0;
t_out=0;t_in=0;flag=0;
x=[];y=[];   %画图坐标
l=l+1;

for i=0:delta_t:20000

flag=flag+1;
x(flag)=i;
y(flag)=p0;

%if i<5000
%   t=0.0001*i+0.26485;
%end
%if i>=5000
%   t=0.76485;
%end
%if i<2000
%   t=-0.0001*i+0.96685;
%end
%if i>=2000
%   t=0.76485;
%end
if i<10000
t=0.00004*i+0.36485;
end
if i>=10000
t=0.76485;
end

vin=0;
if t_in>t+10
t_in=t_in-(t+10);
end
if t_in<t
vin=(C*A*(2*(160-p0)/0.8696)^(.5));
end

t_in=t_in+delta_t;
m_in=0.8696*vin*delta_t;
if t_out>100
t_out=t_out-100;
end

vout=O(t_out)*delta_t;
t_out=t_out+delta_t;
mout=vout*density(p0) ;
delta_m=m_in-mout;

rho_now=(m0+delta_m)/V;
p0=P(rho_now);

m0=m0+delta_m;

flag=flag+1;
x(flag)=i;
y(flag)=p0;
end

%搜索判定条件
%y_ave=0;k=0;
%for j=9950000:10000000
%    k=k+1;
%    y_ave=y_ave+y(j);
%end
%y_ave=y_ave/k;
%y_test(l)=y_ave;
%x_test(l)=t;

plot(x,y)
%plot(x_test,y_test);

## Q2

%密度函数，压力函数
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+1.7*10^-3*x));
P=@(x)-(90071992547409920000*x-72826643121816530000)/(153122387330596864*x-167501279180178019);

c=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
A_oil=(5/2)^2*pi;
H=5.8446;
Vmax=H*A_oil;
M=Vmax*density(0.5);
m0=density(100)*V;
delta_t=0.01;
p_wai=0.1;
L0=0.7/tan(9/180*pi)*2.5/1.4;
A_zhen=(2.5/2)^2*pi;
A_kong=0.7^2*pi;
k=[];x=[];y=[];
bl=2*ones(156,1);
h_zhen=[al;bl;cl;zeros(9755,1)]';

for i=1:628
x(i)=sin(theta(i))*r(i);
y(i)=cos(theta(i))*r(i);
end
s=[];u=1;

for i=0:0.01:6.27
xk=[];
yk=[];
for j=1:628
xk(j)=x(j)*cos(i)-y(j)*sin(i);
yk(j)=x(j)*sin(i)-y(j)*cos(i);
end
s(u)=max(yk)-2.4130;
u=u+1;
end

s=interp1([0:0.01:6.27],s,[0:0.000001:6.27],'spline');
w=0.027249;
j=1;
u=1;
p0=100;
rho0=0.85;
m_oil=M;
fy=1;

for i=0:delta_t:100000
fy=fy+int32(w*delta_t/0.000001);
while fy>6270001
fy=fy-6270001;
m_oil=M;
end
h_now=H-s(fy);
V_now=A_oil*h_now;
rho_now=m_oil/V_now;
P_now=P(rho_now);
vin=0;
if P_now>p0
vin=c*A*(2*(P_now-p0)/rho_now)^(.5);
end
m_in=rho_now*vin*delta_t;
m_oil=m_oil-m_in;
rl=(h_zhen(u)+L0)*tan(9/180*pi);h_z=h_zhen(u);
A_out=pi*rl^2-A_zhen;
if pi*rl^2-A_zhen-A_kong>0
A_out=A_kong;
end
if h_zhen(u)==0
A_out=0;
end
u=u+1;
if u>10001
u=u-10001;
end
vout=c*A_out*(2*(p0-p_wai)/rho0)^(.5);
mout=vout*rho0*delta_t;
delta_m=m_in-mout;
rho0=(m0+delta_m)/V;
k(j)=p0;
p0=P(rho0);
j=j+1;

m0=m0+delta_m;
end

plot(k)
xlabel('时间(t/ms)')
ylabel('压力(P/MPa)')

## Q3

%密度函数，压力函数
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+1.7*10^-3*x));
P=@(x)-(90071992547409920000*x-72826643121816530000)/(153122387330596864*x-167501279180178019);

c=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
A_oil=(5/2)^2*pi;
H=5.8446;
Vmax=H*A_oil;
M=Vmax*density(0.5);
m0=density(100)*V;
delta_t=0.01;
p_wai=0.1;
L0=0.7/tan(9/180*pi)*2.5/1.4;
A_zhen=(2.5/2)^2*pi;
A_kong=0.7^2*pi;
k=[];x=[];y=[];
bl=2*ones(156,1);
h_zhen=[al;bl;cl;zeros(9755,1)]';
yanchi=5000;
for i=10001:-1:yanchi+1
h_zhenl(i)=h_zhen(i-yanchi);
end
for i=1:1:yanchi
h_zhenl(i)=0;
end

for i=1:628
x(i)=sin(theta(i))*r(i);
y(i)=cos(theta(i))*r(i);
end
s=[];u=1;

for i=0:0.01:6.27
xk=[];
yk=[];
for j=1:628
xk(j)=x(j)*cos(i)-y(j)*sin(i);
yk(j)=x(j)*sin(i)-y(j)*cos(i);
end
s(u)=max(yk)-2.4130;
u=u+1;
end

s=interp1([0:0.01:6.27],s,[0:0.000001:6.27],'spline');
w=0.0545;
j=1;
u=1;
p0=100;
rho0=0.85;
m_oil=M;
fy=1;

for i=0:delta_t:100000
fy=fy+int32(w*delta_t/0.000001);
while fy>6270001
fy=fy-6270001;
m_oil=M;
end
h_now=H-s(fy);
V_now=A_oil*h_now;
rho_now=m_oil/V_now;
P_now=P(rho_now);
vin=0;
if P_now>p0
vin=c*A*(2*(P_now-p0)/rho_now)^(.5);
end
m_in=rho_now*vin*delta_t;
m_oil=m_oil-m_in;
rl=(h_zhen(u)+L0)*tan(9/180*pi);
rl2=(h_zhenl(u)+L0)*tan(9/180*pi);

A_out=pi*rl^2-A_zhen;
if pi*rl^2-A_zhen-A_kong>0
A_out=A_kong;
end
if h_zhen(u)==0
A_out=0;
end

Al_out=pi*rl2^2-A_zhen;
if pi*rl2^2-A_zhen-A_kong>0
Al_out=A_kong;
end
if h_zhenl(u)==0
Al_out=0;
end
u=u+1;
if u>10001
u=u-10001;
end
vout=c*(A_out+Al_out)*(2*(p0-p_wai)/rho0)^(.5);
mout=vout*rho0*delta_t;
delta_m=m_in-mout;
rho0=(m0+delta_m)/V;
k(j)=p0;
p0=P(rho0);
j=j+1;

m0=m0+delta_m;
end

plot(k)
xlabel('时间(t/ms)')
ylabel('压力(P/MPa)')

09-21

09-13 2万+
09-19 6977
12-15 1万+
09-19 1万+
09-17 6631
11-16
11-11 1万+
09-17 1066