matlab采用四阶龙格库塔法求解问题求助

2201_75738463

已于 2022-12-08 17:11:39 修改

阅读量314

点赞数

文章标签：矩阵线性代数

于 2022-12-08 17:10:51 首次发布

本文链接：https://blog.csdn.net/2201_75738463/article/details/128240348

版权

i=1;
L=0.01;
while y3>0.0212
Nto=9707.4;
[T,y3]=ode45(@reactor,[0 L],[644.15 0.0810]);
L=L+0.1
i=i+1;
end
for j=1:i
function [p1,p2,p3,p4,p5]=fenya(y3)
p1=3.05*(0.3535-0.0810-y3);p2=3.05*0.1432;p3=3.05*y3;p4=3.05*(0.0488+0.0810-y3);p5=3.05*(0.3735-0.0810+y3);
function [CPm]=CPm(T)
CP1=(6.8712+0.03135*(T/100)+0.14138E-2*(p1/0.101325)-0.6E-6*(p1/0.101325)^2+0.1603E-3*(p1/0.101325)*(T/100))*4.184;
CP2=(4.23329-0.4145*(T/100)+0.072309*(T/100)^2-0.34116E-2*(T/100)^3+0.57726E-2*(p2/0.010325)-0.7404E-3*(p2/0.101325)*(T/100))*8.314;
CP3=(3.86771-0.23279*(T/100)+0.046135*(T/200)^2-0.2186E-2*(T/100)^3+0.42112E-2*(p3/0.101325)-0.4694E-3*(p3/0.101325)*(T/100))*8.314;
CP4=(3.18266+0.53754*(T/100)-0.020125*(T/100)^2+0.018520*(p4/0.101325)-0.22009E-2*(p4/0.101325)*(T/100))*8.314;
CP5=(0.65765-0.049712*(T/100)+0.5269E-3*(T/100)^3+0.020739*(p5/0.101325)-0.27123E-2*(p5/0.101325)*(T/100))*18.0513*4.184;
CPm=CP1*(0.3535-0.0810-y3)+CP2*0.1432+CP3*y3+CP4*(0.0488+0.0810-y3)+CP5*(0.3735-0.0810+y3);
function [HR]=rerong(T)
HR=-(10000+0.219*T-2.845E-3*T^2+0.9703E-6*T^3)*4.184;
function dy=reactor(T,y3);
f(1)=pi*16*1500/Nto*r3*0.8;
f(2)=HR/CPm*f(1);
r3=kT*p3*(p4)^0.5*(1-(p4*p1/Kp/p3/p5));
kT=1.2E-6*(0.101325)^(-0.5)*exp(-104600/Rg/T);Rg=8.314;
Kp=exp(2.3026*(2185/T-0.1102/2.3026*logT+0.6218E-3*T-1.0604E-7*T^2-2.218));
X=[i,L,f(1),f(2),CPm,HR]
y3=X[i,3];
T=X[i,4];
end