24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

lkbiandou

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1230.5
帖子: 208
在线: 94.6小时
虫号: 945069
注册: 2010-01-19
专业: 化学反应工程

[求助] 大家帮忙看一个matlab的小程序，看看哪里出错了

在下自己写了一个程序，但是水平差，写的达不到自己的求解。希望各位帮忙看看问题出在哪里了；代码如下：有没有看明白的地方，我们一起交流，谢谢各位了。

function asshole
clear; clc
i=1:10; j=1:10; k=1:10;
T(i,j,k)=0;
Tin(1:10,1:10,1)=483.15;
x0=Tin(1:10,1:10,1);
S=fsolve(@TxEqs,x0);

function f = TxEqs(T)
global m n K dz dx  i j k
m=10; n=10;K=10;
Ramda = 0.45;    % W/(m K)
Kw=1.2;
Tw=483.15;
rho=1440;Cor=0.9;
R=4.2;
x0=0.3959; dx=0.06041;df=pi/6/10;
dz=0.1;
%Tin=483.15;
P=5;w=5; R0=0.0125;L=10;Dr=1;G=500;Cp=1;M=100;
  %反应器入口
a5=L*M*Ramda/R^2/G/Cp; a6=rho*Cor*L*M/G/Cp*(1+1);
d33=2; d44=3;
%  T(:,:,

;
for k=1:K;
for j=1:n;
      for i=1:m;
         [X,Y]=Reaction(T(i,j,k));
         [U,V]=HRTX(T(i,j,k));
         a5=L*M*Ramda/R0^2/G/Cp; a6=rho*Cor*L*M/G/Cp*(-U*X-Y*V);
      end
end
end
%--------------------------------------------------------------------------
%模型中一些代数式
for k=1:K;
for j=2:n-1;
      for i=2:m-1;
         %一些代数式
         d33(i,j,k)=1/(x0+i*dx)^2*(T(i,j+1,k)-2*T(i,j,k)+T(i,j-1,k))/df^2;
         d44(i,j,k)=(T(i+1,j,k)-2*T(i,j,k)+T(i-1,j,k))/dx^2+1/(x0+i*dx)*(T(i+1,j,k)-T(i-1,j,k))/(2*dx);
         T(i,j,k+1)=T(i,j,k)+dz*(a5*(d33(i,j,k)+d44(i,j,k))+a6);    %当前层计算下一层
      end
end
end
%--------------------------------------------------------------------------
%模型方程的表达式
for k=1:K;
for j=1:n;
      for i=1:m;
         %k=1，入口位置
         T(i,j,1)=483.15;
         %初始条件AD段的处理,i=1
         %T(2,j,k)=T(1,j,k)-dx*Kw*R/Ramda*(T(1,j,k)-Tw);  %变更k层第2行数据，第一层不变
         T(1,j,k)=1/(3+2*dx*Kw*R/Ramda)*(4*T(2,j,k)-T(3,j,k)+2*dx*Kw*R*Tw/Ramda);%刘玉兰的，i=1
         %%这里还应该加入i=m的取值
         %AB段，j=1的处理
         %T(i+1,1,k)=T(i,1,k)+0;          %变更k层第1列数据
         T(i,1,k)=(4*T(i,2,k)-T(i,3,k))/3;%刘玉兰的，j=1
         %CD段，j=n
         %T(i+1,n,k)=T(i,n,k)+0;          %变更k层第n列数据
         T(i,n+1,k)=(4*T(i,n,k)-T(i,n-1,k))/3;%刘玉兰的，j=n
         %除边界以外浓度和温度的表达式
         %d33=2; d44=3;
      end
   end
end
T

O=0;
for k=1:K;
for j=1:n;
      for i=1:m;
         [X,Y]=Reaction(T(i,j,1));
         O=O+1;
      end
end
end
X
Y
for k=1:K;
for j=1:n;
      for i=1:m;
         [U,V]=HRTX(T(i,j,1));
         O=O+1;
      end
end
end
%[U,V]=HRTX(T);
U
V

function [HCO,HCO2]=HRTX(T)
%%计算CO和CO2的焓变
HCO=[-18288.6-11.7808*T-0.700294*10^-2*T^2+0.406434*10^-4*T^3-0.458711*10^-7*T^4+0.189644*10^-10*T^5]*4.184;
HCO2=[-9048.33-5.4173*T-0.028276*T^2+0.683828*10^-4*T^3-0.659525*10^-7*T^4+2.541207*10^-11*T^5]*4.184;
[HCO;HCO2];
%--------------------------------------------------------------------------

function [r1,r2]=Reaction(T)
% Calculate the reaction rate（计算反应速度）
%k = 0.027*exp(0.021*(T-773));
%f = 15100*exp(-11000./T).*((1-x)./(11+x)-1.2*x.^2./k./(11+x).^2);
yMin=0.0182;yCOin=0.1128; yH2in=0.7717; yCO2in=0.0209; yN2in=0.0735; yH2Oin=0.0001;
y_CO=0.08;y_CO2=0.02;P=5;%%%%%%%%%%%%%%%%%%这里的反应物浓度需要跟实际反应配合起来
TBA=490.15;
BB=1-2*yCOin-2*yCO2in;
BA=1-2*y_CO-2*y_CO2;
yM=BA/BB*(yCOin+yMin+yCO2in)-y_CO-y_CO2;
yH2=BA/BB*(-2*yCOin-3*yCO2in+yH2in)+2*y_CO+3*y_CO2;
yH2O=BA/BB*(yCO2in+yH2Oin)-y_CO2;
yN2=BA/BB*yN2in;
P_Ca=P/0.101325;
KF1=exp(13.1652+9203.26/T-5.92839*log(T)-0.00352404*T+0.0000102264*T*T-0.00000000769446*T^3+2.38583E-12*T^4)*(0.101325^(-2));
KF2=exp(1.6654+4553.34/T-2.72613*log(T)-0.01422914*T+0.000017206*T*T-0.00000001106294*T^3+3.19698E-12*T^4)*(0.101325^(-2));
fH2=P*yH2*exp((0.110785+35.3324/T-5005.47/(T*T*T)-19.6109*yH2/T-20.9799*yH2*yH2/T)*P_Ca/T);
fCO2=P*y_CO2*exp((-0.343605+428.452/T-69217700/(T*T*T)-327.402*y_CO2/T-374.954*y_CO2*y_CO2/T)*P_Ca/T);
fCO=P*y_CO*exp((-0.093261+189.156/T-399940/(T*T*T)-181.527*y_CO/T+140.001*y_CO*y_CO/T)*P_Ca/T);
fM=P*yM*exp((-1.49696+997.85/T-100000000/(T*T*T)-792.109*yM/T-803.4*yM*yM/T)*P_Ca/T);
fH2O=P*yH2O*exp((-1.78527+1408.49/T-183959000/(T*T*T)-3648.32*yH2O/T-3116.5*yH2O*yH2O/T)*P_Ca/T);
%--------------------------------------------------------------------------
KCO=exp(0.22349-7.6694*10^3*(1/T-1/TBA));    %此处需要输入参数···暂时采用这个
KCO2=exp(-4.8272-8.5623*10^3*(1/T-1/TBA));
KH2=exp(-0.15458+1.2853*10^3*(1/T-1/TBA));
KT1=1.7299*10^3*exp(-36178.22/8.314/T)/3600;
KT2=1.1426*10^4*exp(-50484.09/8.314/T)/3600;
%--------------------------------------------------------------------------
Beta1=fM/(KF1*fCO*fH2^2);
Beta2=fM*fH2O/(KF2*fCO2*fH2^3);
%-------------------------------------------------------------------------
r1=KT1*fCO*(fH2^2)*(1-Beta1)/((1+KCO*fCO+KCO2*fCO2+KH2*fH2)^3);    %CO反应速率
r2=KT2*fCO2*(fH2^3)*(1-Beta2)/((1+KCO*fCO+KCO2*fCO2+KH2*fH2)^4); %CO2反应速率
[r1;r2];

结果总是出不来。大家帮我看看怎么办。

回复此楼