| 查看: 1453 | 回复: 5 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
大家帮忙看一个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]; 结果总是出不来。大家帮我看看怎么办。 |
回复此楼» 猜你喜欢
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有111人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
-
什么时候开奖？
已经有13人回复
-
散金币祈福
已经有94人回复
-
青基已中
已经有1人回复
-
散金币祈福
已经有185人回复
-
求助，如何提取ELK的rt-TDDFT在某一时刻的自旋密度分布和电子密度分布
已经有1人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 请matlab高手看一下这个问题怎样编程才好，附程序代码，谢谢！
已经有9人回复
- 一段matlab的小程序，不明白它的意思
已经有5人回复
- 请大侠帮忙看一下我这个matlab小程序结果为什么不返回矩阵
已经有7人回复
- MATLAB 如何调用程序代码？
已经有5人回复
- 哪一本matlab程序设计的书比较好？？
已经有12人回复
- 求助Matlab中如何把几个比特合为一个symbol
已经有4人回复
- 自己写的一个小程序，出现问题不知道怎么解决
已经有9人回复
- 大家爱帮我看看这段MATLAB程序为啥出错了
已经有5人回复
- 求助，matlab公式绘图，公式里面有两个数组，编出小程序，但不出图像
已经有4人回复
- matlab 运行出错，求帮忙啦~~·
已经有3人回复
- matlab程序问题，总是出错，麻烦大家看看问题在那里？
已经有6人回复
- 初学者伤不起啊——关于matlab...
已经有12人回复
- 用matlab编一个计算稳态概率的小程序
已经有4人回复
- matlab编一个小程序解三角形
已经有4人回复
- 请教一个fortran小程序编译出错的问题，谢谢
已经有9人回复
- 【求助】请教一个matlab程序，谢谢帮忙！！
已经有15人回复
- 【求助】刚学MATLAB有个小程序有一句看不懂 请大侠指教
已经有4人回复
- 【求助】求一个关于文件读写的c语言小程序
已经有3人回复
- 【求助】帮忙看一个程序
已经有16人回复
- 【求助】帮忙改写展宽程序
已经有11人回复
- 【求助】Weickert的各向异性扩散方程的滤波方法的matlab程序代码
已经有8人回复
3楼2013-08-06 19:32:28
regmastercn1
铁杆木虫 (正式写手)
- 应助: 10 (幼儿园)
- 金币: 5466.4
- 散金: 482
- 红花: 7
- 帖子: 983
- 在线: 739.8小时
- 虫号: 1269863
- 注册: 2011-04-19
- 专业: 水文、水资源
2楼2013-08-05 21:47:35
4楼2013-08-28 23:26:29
5楼2013-08-29 01:42:27
5
投票:
