[MATLAB毕业设计] matlab有限差分法求解偏微分方程出错～～～

[lkbiandou]

在求解过程中，模型方程也给出了，就是方程中一个代数的计算结果有问题，请教给为帮我看一下，我想让每一个T(i,j,k)都对应一个反应速率rc，但是现在求解出来的结果却是每一页上的T（i，j，k）都只能对应得到一个rc，而不是rc(i,j,k).麻烦大家帮我看看吧。为了这个我已经耗费太多的时间进去了！
CODE：
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=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];
%------------------------------------------------------------------------------



感谢各位了！！！谢谢！！！