| 3 | 1/1 | 返回列表 |
| 查看: 1342 | 回復(fù): 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 ] |
» 搶金幣啦!回帖就可以得到:
+1/189
+1/181
+5/120
+1/90
+1/58
+1/48
+1/45
+1/41
+1/35
+1/29
+1/14
+2/12
+1/8
+1/8
+1/6
+1/5
+1/3
+1/2
+1/1
+1/1
| 3 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 求調(diào)劑 +5 | Zhangbod 2026-03-21 | 7/350 |
|
|---|---|---|---|---|
|
[考研] 085600材料與化工306 +4 | z1z2z3879 2026-03-21 | 4/200 |
|
|
求調(diào)劑 +3 | 13341 2026-03-20 | 3/150 |
|
|
[考研] 313求調(diào)劑 +4 | 肆叁貳壹22 2026-03-19 | 4/200 |
|
|
[考研] 306求0703調(diào)劑一志愿華中師范 +5 | 紙魚ly 2026-03-21 | 5/250 |
|
|
[考研] 268求調(diào)劑 +9 | 簡單點0 2026-03-17 | 9/450 |
|
|
[考研] 求調(diào)劑 +3 | 白QF 2026-03-21 | 3/150 |
|
|
[考研] 南昌大學材料專碩311分求調(diào)劑 +6 | 77chaselx 2026-03-20 | 6/300 |
|
|
[考研] 化學求調(diào)劑 +4 | 臨澤境llllll 2026-03-17 | 5/250 |
|
|
[考研] 一志愿南昌大學,327分,材料與化工085600 +9 | Ncdx123456 2026-03-19 | 9/450 |
|
|
[考研] 304求調(diào)劑 +7 | 司空. 2026-03-18 | 7/350 |
|
|
[考研] 289求調(diào)劑 +6 | 懷瑾握瑜l 2026-03-20 | 6/300 |
|
|
[考研] 求調(diào)劑 +3 | 暗涌afhb 2026-03-16 | 3/150 |
|
|
[考研] 085600材料與化工求調(diào)劑 +6 | 緒幸與子 2026-03-17 | 6/300 |
|
|
[考研] 材料工程專碩調(diào)劑 +5 | 204818@lcx 2026-03-17 | 6/300 |
|
|
[考博] 26博士申請 +3 | 1042136743 2026-03-17 | 3/150 |
|
|
[考研] 326求調(diào)劑 +5 | 上岸的小葡 2026-03-15 | 6/300 |
|
|
[考研] 有沒有道鐵/土木的想調(diào)劑南林,給自己招師弟中~ +3 | TqlXswl 2026-03-16 | 7/350 |
|
|
[考研] 一志愿蘇州大學材料工程(085601)專碩有科研經(jīng)歷三項國獎兩個實用型專利一項省級立項 +6 | 大火山小火山 2026-03-16 | 8/400 |
|
|
[考研] 一志愿南京大學,080500材料科學與工程,調(diào)劑 +4 | Jy? 2026-03-16 | 4/200 |
|