| 查看: 1507 | 回复: 2 | |||
[交流]
【求助】求高手帮看一下程序（泡点温度计算）
|
|
使用SRK方程和wilson方程计算泡点温度和组成 m文件中，y（1）是组分1 的气相组成，y（2）是温度 m文件如下 function f=myfun(y,x1) p=101325; T1c=579.35;V1c=219E-06;Z1c=0.259;p1c=56.9E+05;w1=0.193; ANTA1=36.6016;ANTB1=-2979.4;ANTC1=-10.104;ANTD1=1.1445E-09;ANTE1=3.2472E-06; T2c=591.79;V2c=315.8e-06;Z2c=0.264;p2c=41.09e+05;w2=0.264; ANTA2=34.0755;ANTB2=-3.0379e+03;ANTC2=-9.1635;ANTD2=1.0289E-11;ANTE2=2.7035E-06; b=[4198.6 -3001.6]; x2=1-x1; p1s=10^(ANTA1+ANTB1/y(2)+ANTC1*log10(y(2))+ANTD1*y(2)+ANTE1*y(2)* y(2))*133.322; p2s=10^(ANTA2+ANTB2/y(2)+ANTC2*log10(y(2))+ANTD2*y(2)+ANTE2*y(2)* y(2))*133.322; SRKa1=0.42748*8.31441*8.31441*T1c*T1c/p1c*(1+(0.480+1.574*w1-0.176*w1^2)*(1-(y(2)/T1c)^0.5))^2; SRKb1=0.08664*8.31441*T1c/p1c; SRKk1=-8.31441*y(2)/p1s; SRKm1=SRKa1/p1s-SRKb1*8.31441*y(2)/p1s-(SRKb1)^2; SRKn1=-SRKa1*SRKb1/p1s; SRKp1=[1 SRKk1 SRKm1 SRKn1]; SRKV1=roots(SRKp1); V1=max(SRKV1); Z1=p1s*V1/8.31441/y(2); fai1s=exp((Z1-1)-log(p1s*(V1-SRKb1)/8.31441/y(2))-SRKa1*log(1+SRKb1/SRKV1)/SRKb1/8.31441/y(2)); SRKa2=0.42748*8.31441*8.31441*T2c*T2c/p2c*(1+(0.480+1.574*w2-0.176*w2^2)*(1-(y(2)/T2c)^0.5))^2; SRKb2=0.08664*8.31441*T2c/p2c; SRKk2=-8.31441*y(2)/p2s; SRKm2=SRKa2/p2s-SRKb2*8.31441*y(2)/p2s-(SRKb2)^2; SRKn2=-SRKa2*SRKb2/p2s; SRKp2=[1 SRKk2 SRKm2 SRKn2]; SRKV2=roots(SRKp2); V2=max(SRKV2); Z2=p2s*V2/8.31441/y(2); fai2s=exp((Z2-1)-log(p2s*(V2-SRKb2)/8.31441/y(2))-SRKa2*log(1+SRKb2/V2)/SRKb2/8.31441/y(2)); Tc12=(T1c*T2c)^0.5; Vc12=((nthroot(V1c,3)+nthroot(V2c,3))/2)^3; Zc12=(Z1c+Z2c)/2; pc12=Zc12*8.31441*Tc12/Vc12; w12=(w1+w2)/2; SRKa12=(SRKa1*SRKa2)^0.5; SRKam=y(1)*y(1)*SRKa1+2*y(1)*(1-y(1))*SRKa12+(1-y(1))*(1-y(1))*SRKa2; SRKbm=y(1)*SRKb1+(1-y(1))*SRKb2; SRKkm=-8.31441*y(2)/p; SRKmm=SRKam/p-SRKbm*8.31441*y(2)/p-(SRKbm)^2; SRKnm=-SRKam*SRKbm/p; SRKpm=[1 SRKkm SRKmm SRKnm]; SRKVm=roots(SRKpm); Vm=max(SRKVm); Zm=p*Vm/8.31441/y(2); fai1=exp(SRKb1/SRKbm*(Zm-1)-log(p*(Vm-SRKbm)/8.31441/y(2))+SRKam*(SRKb1/SRKbm-2*(y(1)*SRKa1+(1-y(1))*SRKa12)/SRKam)*log(1+SRKbm/SRKVm)/SRKbm/8.31441/y(2)); fai2=exp(SRKb2/SRKbm*(Zm-1)-log(p*(Vm-SRKbm)/8.31441/y(2))+SRKam*(SRKb2/SRKbm-2*(y(1)*SRKa12+(1-y(1))*SRKa2)/SRKam)*log(1+SRKbm/SRKVm)/SRKbm/8.31441/y(2)); V1L=V1c*Z1c^((1-y(2)/T1c)^(2/7)); V2L=V2c*Z2c^((1-y(2)/T2c)^(2/7)); A12=(V2L/V1L)*exp(-b(1)/8.314/y(2)); A21=(V1L/V2L)*exp(-b(2)/8.314/y(2)); r1cal=exp(-log(x1+A12*x2)+x2*(A12/(x1+A12*x2)-A21/(x2+A21*x1))); r2cal=exp(-log(x2+A21*x1)-x1*(A12/(x1+A12*x2)-A21/(x2+A21*x1))); f=abs(y(1)*p*fai1-r1cal*p1s*fai1s*x1*exp(V1L*(p-p1s)/8.31441/y(2)))+abs((1-y(1))*p*fai2-r2cal*p2s*fai2s*x2*exp(V2L*(p-p2s)/8.31441/y(2))); main文件如下 x1=[0.9643 0.9165 0.8322 0.7520 0.6696 0.6362 0.5576 0.4416 0.3864 0.3407 0.2988 0.2628 0.2271 0.1877 0.1482 0.0933 0.0586 0.0375]; n=18; for i=1:n x0=[1;355]; [x fval]=fminsearch(@(x) myfun(x,x1(i)),x0) x1(i)=x(1); x2(i)=x(2); end; x1; x2; 运行提示Subscripted assignment dimension mismatch. Error in ==> fminsearch at 205 fv(:,1) = funfcn(x,varargin{:}); [ Last edited by flying_fay on 2010-12-1 at 11:26 ] |
回复此楼» 猜你喜欢
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急
已经有5人回复
-
有多少人是今天查系统知道结果的？
已经有18人回复
-
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有5人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有4人回复
-
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急
已经有5人回复
-
我就是申请一个面上项目而已，这评审意见是按照杰青的条件评的吧？
已经有4人回复
-
国自然评审意见
已经有15人回复
-
为什么资助数各大高校都创新高，自己申请怎么就这么难
已经有13人回复
-
国自然面上复盘~欢迎讨论
已经有14人回复
-
系统查不到
已经有11人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 求高手帮忙看一下我的结构式画的有没有问题？
已经有16人回复
- 求matlab编程代码，求帮助，谢谢！
已经有11人回复
- 跪求高手帮忙看看这组数据如何进行拟合
已经有7人回复
- 很简单的Fortran赋值出错啦，求大神们帮忙找bug...
已经有11人回复
- 求高手帮忙看一段英文，涉及到数据分析问题，如能给出说得通的全文翻译另有50BB相赠
已经有5人回复
- 高手 帮忙看看这两数据
已经有5人回复
- 求高手分析一下DSC 曲线
已经有10人回复
- 请matlab高手看一下这个问题怎样编程才好，附程序代码，谢谢！
已经有9人回复
- 麻烦高手帮忙解一下三元二次方程组！谢谢谢谢！
已经有4人回复
- 高手帮我看看，这种示意图用什么软件画比较方便？
已经有7人回复
- matlab高手帮忙做个自定义复杂函数拟合，含有复数，
已经有25人回复
- 硫酸铵沉淀的一些问题，求高手帮忙看一下！
已经有12人回复
- 求SWMM软件高手指导
已经有15人回复
- 帮忙计算半衰期
已经有13人回复
- 急！！！！液质图谱超诡异~~~不知道是不是被污染了 ~~~求高手帮忙看一下
已经有21人回复
- 求有机合成高人帮在下看一下这个嘧啶成环反应能发生吗？？？？？？？
已经有18人回复
- 求高人帮忙看一下我这个Fortran程序的错误出在哪里？
已经有12人回复
- 求助 哪位傅里叶红外高人帮忙看一下我的红外图谱
已经有18人回复
- 请各位高人看一下这够不够JBC的要求？
已经有4人回复
- 运行一个python程序，出现问题，请高手赐教
已经有9人回复
- 刚投就被拒稿，请高手帮忙分析一下
已经有18人回复
- LABVIEW四阶龙格库塔法求解一阶微分方程组
已经有12人回复
2楼2010-12-07 09:18:15
3楼2010-12-07 10:40:43
