应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】请高手帮我解答个比较疑惑的问题（UNIFAC估算）
51/1返回列表
查看: 1305  |  回复: 4
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

wzhr

新虫 (小有名气)

[交流] 【求助】请高手帮我解答个比较疑惑的问题（UNIFAC估算）

最近在做一个计算，使用基团贡献法UNIFAC-Dortmund来算液液平衡数据，不知道自己编的程序哪里出问题了，一种物质环丁砜两种不同的基团拆分法，计算出来的结果竟然一样，请高人看看我的程序哪里有问题？如果谁有这方面现成的程序，不知能否给我一份作为学习，谢谢各位！
其中各中参数为
己烷-环丁砜的液液平衡数据计算，己烷划分为2个CH3、4个CH2，环丁砜分为2个C-CH2及1个CH2SUCH2
其 MAIN-           SUB-        
GROUP          GROUP K                R(K)    Q(K)
    1 CH2          1 CH3                   0.6325  1.0608
    56 SULFONE    110 (CH2)2SU    2.6870  2.1200
   42 CY-CH2      78 CY-CH2           0.7136  0.8635
REQUIRED UNIFAC INTERACTION PARAMETERS:
  N  M   A(N,M)      A(M,N)      B(N,M)      B(M,N)      C(N,M)      C(M,N)
  1 56  1058.3      438.76     -1.1856     -1.2256      0.0000      0.0000   
  1 42 -117.10      170.90     0.54810    -0.80620    -0.98000E-03 0.12910E-02
56 42  660.47     -581.16     -2.0132      4.1194      0.0000      0.0000   

function f=fun(x,T)
q1=2*1.060+0.7081*4;
r1=0.6325*6;
q2=2.12*1+2*0.8635;
r2=2.687*1+2*0.7136;
x11=x(1);
x21=x(2);
x12=1-x11;
x22=1-x21;
V11=r1./(r1*x11+r2*x12);   
V12=r2./(r1*x11+r2*x12);
F11= q1./(q1*x11+q2*x12);
F12= q2./(q1*x11+q2*x12);
W11= r1^0.75./(r1^0.75*x11+r2^0.75*x12);
W12=r2^0.75./(r1^0.75*x11+r2^0.75*x12);
logr11c=1-W11+log(W11)-5*q1*(1-V11./F11+log(V11./F11));
logr12c=1-W12+log(W12)-5*q2*(1-V12./F12+log(V12./F12));


AF111a=1.0608*1/3/(1.0608*1/3+0.7081*2/3);
AF111b=0.7081*2/3/(1.0608*1/3+0.7081*2/3);

AF122a=1/3*2.12/(1/3*2.12+2/3*0.8635);
AF122b=2/3*0.8635/(1/3*2.12+2/3*0.8635);

I1a1b=1;
I1b1a=1;
I1a2a=exp(-(1058.31-1.1856.*T+0.*T.^2)./T);
I1b2a= exp(-(1058.31-1.1856.*T+0.*T.^2)./T);
I2b2a=exp(-(-581.16+4.1194.*T+0)./T);
I1a2b=exp(-(-117.1+0.5481.*T-0.00098*T.^2)./T);
I2b1a=exp(-(170.9-0.8062.*T+0.001291*T.^2)./T);
I1b2b= exp(-(-117.1+0.5481.*T-0.00098*T.^2)./T);
I2b1b= exp(-(170.9-0.8062.*T+0.001291*T.^2)./T);


I2a1a=exp(-(438.76-1.2256.*T+0.*T.^2)./T);
I2a1b= exp(-(438.76-1.2256.*T+0.*T.^2)./T);
I2a2b=exp(-(660.47-2.0132.*T+0)./T);
I1a1a=1; I1b1b=1; I2a2a=1; I2b2b=1;


TTa=(I1a1b.*AF111b./(I1a1b.*AF111a+I1b1b.*AF111b)+I1a1a.*AF111a./(I1a1a.*AF111a+I1b1a.*AF111b));
T111a=exp(1.0608*(1-log((I1b1a.*AF111b + I1a1a.*AF111a))-TTa));
TTb=(I1b1a.*AF111a./(I1b1a.*AF111b+I1a1a.*AF111a)+I1b1b.*AF111b./( I1a1b.*AF111a+I1b1b.*AF111b));
T111b=exp(0.7081*(1-log((I1a1b.*AF111a + I1b1b.*AF111b))-TTb));

T122a=exp(2.12*(1-log(I2b2a.*AF122b+I2a2a.*AF122a)-I2a2b.*AF122b./(I2a2b.*AF122a+I2b2b.*AF122b)-I2a2a.*AF122a./(I2b2a.*AF122b+I2a2a.*AF122a)));
T122b=exp(0.8635*(1-log(I2a2b.*AF122a+I2b2b.*AF122b)- I2b2a.*AF122a./ (I2b2a.*AF122b+I2a2a.*AF122a)-I2b2b.*AF122b./(I2a2b.*AF122a+I2b2b.*AF122b)));


X11a=2*x11./(6*x11+3*x12);
X11b=4*x11./( 6*x11+3*x12);
X12a=1*x12./( 6*x11+3*x12);
X12b=2*x12./( 6*x11+3*x12);

AF11a=1.0608.*X11a./(1.0608.*X11a+0.7081.*X11b+X12a*2.12+X12b*0.8635);
AF11b=0.7081.*X11b./(1.0608.*X11a+0.7081.*X11b+X12a*2.12+X12b*0.8635);
AF12a= X12a*2.12./( 1.0608.*X11a+0.7081.*X11b +X12a*2.12+X12b*0.8635);
AF12b=X12b*0.8635./(1.0608.*X11a+0.7081.*X11b+X12a*2.12+X12b*0.8635);

TT11a=I1a1b.*AF11b./(I1a1b.*AF11a+I2a1b.*AF12a+I2b1b.*AF12b+I1b1b.*AF11b)+I1a2a.*AF12a./(I1a2a.*AF11a+I1b2a.*AF11b+I2b2a.*AF12b+I2a2a.*AF12a)+I1a2b.*AF12b./(I1a2b.*AF11a+I2a2b.*AF12a+I1b2b.*AF11b+I2b2b.*AF12b)+I1a1a.*AF11a./(I1b1a.*AF11b+I2a1a.*AF12a+I2b1a.*AF12b+I1a1a.*AF11a);
T11a=exp(1.0608*(1-log((I1b1a.*AF11b+I2a1a.*AF12a+I2b1a.*AF12b+I1a1a.*AF11a))-TT11a));

TT11b=I1b1a.*AF11a./(I1b1a.*AF11b+I2a1a.*AF12a+I2b1a.*AF12b+I1a1a.*AF11a)+I1b2a.*AF12a./(I1a2a.*AF11a+I1b2a.*AF11b+I2b2a.*AF12b+I2a2a.*AF12a)+I1b2b.*AF12b./(I1a2b.*AF11a+I2a2b.*AF12a+I1b2b.*AF11b+I2b2b.*AF12b)+I1b1b.*AF11b./(I1a1b.*AF11a+I2a1b.*AF12a+I2b1b.*AF12b+I1b1b.*AF11b);
T11b=exp(0.7081*(1-log((I1a1b.*AF11a+I2a1b.*AF12a+I2b1b.*AF12b+I1b1b.*AF11b))-TT11b));

TT12a=I2a1a.*AF11a./(I1b1a.*AF11b+I2a1a.*AF12a+I2b1a.*AF12b+I1a1a.*AF11a)+I2a1b.*AF11b./(I1a1b.*AF11a+I2a1b.*AF12a+I2b1b.*AF12b+I1b1b.*AF11b)+I2a2b.*AF12b./(I1a2b.*AF11a+I2a2b.*AF12a+I1b2b.*AF11b+I2b2b.*AF12b)+I2a2a.*AF12a./(I1a2a.*AF11a+I2a2a.*AF12a+I2b2a.*AF12b+I1b2a.*AF11b);
T12a=exp(2.12*(1-log((I1a2a.*AF11a+I1b2a.*AF11b+I2b2a.*AF12b+I2a2a.*AF12a))-TT12a));

TT12b=I2b1a.*AF11a./(I1b1a.*AF11b+I2a1a.*AF12a+I2b1a.*AF12b+I1a1a.*AF11a)+I2b1b.*AF11b./(I1a1b.*AF11a+I2a1b.*AF12a+I2b1b.*AF12b+I1b1b.*AF11b)+I2b2a.*AF12a./(I1a2a.*AF11a+I2b2a.*AF12b+I1b2a.*AF11b+I2a2a.*AF12a)+I2b2b.*AF12b./(I1a2b.*AF11a+I1b2b.*AF11b+I2a2b.*AF12a+I2b2b.*AF12b);
T12b=exp(0.8635*(1-log((I1a2b.*AF11a+I1b2b.*AF11b+I2a2b.*AF12a+I2b2b.*AF12b))-TT12b));

logr11R=2*(log(T11a)-log(T111a))+4*(log(T11b)-log(T111b));

logr12R=(log(T12a)-log(T122a))+2*(log(T12b)-log(T122b));

logr11=logr11c+logr11R;
logr12=logr12c+logr12R;
r11=exp(logr11);
r12=exp(logr12);





V21=r1./(r1*x21+r2*x22);
V22=r2./(r1*x21+r2*x22);
F21= q1./(q1*x21+q2*x22);
F22= q2./(q1*x21+q2*x22);
W21= r1^0.75./(r1^0.75*x21+r2^0.75*x22);
W22=r2^0.75./(r1^0.75*x21+r2^0.75*x22);
logr21c=1-W21+log(W21)-5*q1*(1-V21./F21+log(V21./F21));
logr22c=1-W22+log(W22)-5*q2*(1-V22./F22+log(V22./F22));


AF211a=1.0608*1/3/(1.0608*1/3+0.7081*2/3);
AF211b=0.7081*2/3/(1.0608*1/3+0.7081*2/3);

AF222a=1/3*2.12/(1/3*2.12+2/3*0.8635);
AF222b=2/3*0.8635/(1/3*2.12+2/3*0.8635);

I1a2a=exp(-(1058.31-1.1856.*T+0.*T.^2)./T);
I1b2a= exp(-(1058.31-1.1856.*T+0.*T.^2)./T);
I2b2a=exp(-(-581.16+4.1194.*T+0)./T);
I1a2b=exp(-(-117.1+0.5481.*T-0.00098*T.^2)./T);
I2b1a=exp(-(170.9-0.8062.*T+0.001291*T.^2)./T);
I1b2b= exp(-(-117.1+0.5481.*T-0.00098*T.^2)./T);
I2b1b= exp(-(170.9-0.8062.*T+0.001291*T.^2)./T);


I2a1a=exp(-(438.76-1.2256.*T+0.*T.^2)./T);
I2a1b= exp(-(438.76-1.2256.*T+0.*T.^2)./T);
I2a2b=exp(-(660.47-2.0132.*T+0)./T);
I1a1a=1; I1b1b=1; I2a2a=1; I2b2b=1;


TTa=(I1a1b.*AF211b./(I1a1b.*AF211a+I1b1b.*AF211b)+I1a1a.*AF211a./(I1a1a.*AF211a+I1b1a.*AF211b));
T211a=exp(1.0608*(1-log((I1b1a.*AF211b + I1a1a.*AF211a))-TTa));
TTb=(I1b1a.*AF211a./(I1b1a.*AF211b+I1a1a.*AF211a)+I1b1b.*AF211b./( I1a1b.*AF211a+I1b1b.*AF211b));
T211b=exp(0.7081*(1-log((I1a1b.*AF211a + I1b1b.*AF211b))-TTb));

T222a=exp(2.12*(1-log(I2b2a.*AF222b+I2a2a.*AF222a)-I2a2b.*AF222b./(I2a2b.*AF222a+I2b2b.*AF222b)-I2a2a.*AF222a./(I2b2a.*AF222b+I2a2a.*AF222a)));
T222b=exp(0.8635*(1-log(I2a2b.*AF222a+I2b2b.*AF222b)- I2b2a.*AF222a./ (I2b2a.*AF222b+I2a2a.*AF222a)-I2b2b.*AF222b./(I2a2b.*AF222a+I2b2b.*AF222b)));


X21a=2*x21./(6*x21+3*x22);
X21b=4*x21./( 6*x21+3*x22);
X22a=1*x22./( 6*x21+3*x22);
X22b=2*x22./( 6*x21+3*x22);

AF21a=1.0608.*X21a./(1.0608.*X21a+0.7081.*X21b+X22a*2.12+X22b*0.8635);
AF21b=0.7081.*X21b./(1.0608.*X21a+0.7081.*X21b+X22a*2.12+X22b*0.8635);
AF22a= X22a*2.12./(1.0608.*X21a+0.7081.*X21b +X22a*2.12+X22b*0.8635);
AF22b=X22b*0.8635./(1.0608.*X21a+0.7081.*X21b+X22a*2.12+X22b*0.8635);

TT21a=I1a1b.*AF21b./(I1a1b.*AF21a+I2a1b.*AF22a+I2b1b.*AF22b+I1b1b.*AF21b)+I1a2a.*AF22a./(I1a2a.*AF21a+I1b2a.*AF21b+I2b2a.*AF22b+I2a2a.*AF22a)+I1a2b.*AF22b./(I1a2b.*AF21a+I2a2b.*AF22a+I1b2b.*AF21b+I2b2b.*AF22b)+I1a1a.*AF21a./(I1b1a.*AF21b+I2a1a.*AF22a+I2b1a.*AF22b+I1a1a.*AF21a);
T21a=exp(1.0608*(1-log((I1b1a.*AF21b+I2a1a.*AF22a+I2b1a.*AF22b+I1a1a.*AF21a))-TT21a));

TT21b=I1b1a.*AF21a./(I1b1a.*AF21b+I2a1a.*AF22a+I2b1a.*AF22b+I1a1a.*AF21a)+I1b2a.*AF22a./(I1a2a.*AF21a+I1b2a.*AF21b+I2b2a.*AF22b+I2a2a.*AF22a)+I1b2b.*AF22b./I1a2b.*AF21a+I2a2b.*AF22a+I1b2b.*AF21b+I2b2b.*AF22b)+I1b1b.*AF21b./(I1a1b.*AF21a+I2a1b.*AF22a+I2b1b.*AF22b+1b1b.*AF21b);
T21b=exp(0.7081*(1-log((I1a1b.*AF21a+I2a1b.*AF22a+I2b1b.*AF22b+I1b1b.*AF21b))-TT21b));

TT22a=I2a1a.*AF21a./(I1b1a.*AF21b+I2a1a.*AF22a+I2b1a.*AF22b+I1a1a.*AF21a)+I2a1b.*AF21b./(I1a1b.*AF21a+I2a1b.*AF22a+I2b1b.*AF22b+I1b1b.*AF21b)+I2a2b.*AF22b./(I1a2b.*AF21a+I2a2b.*AF22a+I1b2b.*AF21b+I2b2b.*AF22b)+I2a2a.*AF22a./(I1a2a.*AF21a+I2a2a.*AF22a+I2b2a.*AF22b+I1b2a.*AF21b);
T22a=exp(2.12*(1-log((I1a2a.*AF21a+I1b2a.*AF21b+I2b2a.*AF22b+I2a2a.*AF22a))-TT22a));

TT22b=I2b1a.*AF21a./(I1b1a.*AF21b+I2a1a.*AF22a+I2b1a.*AF22b+I1a1a.*AF21a)+I2b1b.*AF21b./(I1a1b.*AF21a+I2a1b.*AF22a+I2b1b.*AF22b+I1b1b.*AF21b)+I2b2a.*AF22a./I1a2a.*AF21a+I2b2a.*AF22b+I1b2a.*AF21b+I2a2a.*AF22a)+I2b2b.*AF22b./(I1a2b.*AF21a+I1b2b.*AF21b+I2a2b.*AF22a+I2b2b.*AF22b);
T22b=exp(0.8635*(1-log((I1a2b.*AF21a+I1b2b.*AF21b+I2a2b.*AF22a+I2b2b.*AF22b))-TT22b));

logr21R=2*(log(T21a)-log(T211a))+4*(log(T21b)-log(T211b));

logr22R=(log(T22a)-log(T222a))+2*(log(T22b)-log(T222b));

logr21=logr21c+logr21R;
logr22=logr22c+logr22R;
r21=exp(logr21);
r22=exp(logr22);

f=[x(1)*r11-x(2)*r21;(1-x(1))*r12-(1-x(2))*r22];

b=[ 300.3
307.43
323.51
348.72
363.36
372.41
377.72
384.43
393.92
403.51
412.29
422.92
429.93;
for i=1:13
x0=[0.99 0.0184];
T=b(i);
[x,fval]=fsolve(@(x) fun(x, T),x0);
x11(i)=x(1);
x21(i)=x(2);
fv(i,1:2)=fval;
end;


另一种划分只是把基团个数和基团参数换了一下，其它都一样，可是出来结果一样，不知哪里错了？

[ Last edited by wzhr on 2011-3-14 at 20:50 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:

» 抢金币啦！回帖就可以得到:

查看全部散金贴
1楼 2011-03-14 20:20:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

信彼南山

木虫 (著名写手)
★
xiegangmai(金币+1): 鼓励讨论交流！ 2011-03-15 06:48:57
wzhr(金币+5): 2011-03-15 15:44:16
程序只是算法的体现，没有算法和公式，怎么看程序对不对啊？
赞一下(1人) 回复此楼
2楼2011-03-14 23:49:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wzhr

新虫 (小有名气)
★
xiegangmai(金币+1): 鼓励讨论交流！ 2011-03-15 11:58:20
引用回帖:
Originally posted by 信彼南山 at 2011-03-14 23:49:11:
程序只是算法的体现，没有算法和公式，怎么看程序对不对啊？

UNIFAC-Dortmund公式应当比较常见，所以我就上传了，我的依据就是根据UNIFAC-Dortmund公式算活度系数，但是x11，x12，x21，x22均为未知，已知只有一系列温度T，于是还需要联立x11r11=x21r21，x12r12=x22r22，x11+x12=1，x21+x22=1来求解，相当求解八个方程八个未知数的方程组。

[ Last edited by wzhr on 2011-3-15 at 11:22 ]
赞一下(1人) 回复此楼
3楼2011-03-15 10:29:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

信彼南山

木虫 (著名写手)
引用回帖:
Originally posted by wzhr at 2011-03-15 10:29:59:
UNIFAC-Dortmund公式应当比较常见，所以我就上传了，我的依据就是根据UNIFAC-Dortmund公式算活度系数，但是x11，x12，x21，x22均为未知，已知只有一系列温度T，于是还需要联立x11r11=x21r21，x12r12=x22r22， ...

这个算法我不知道，可能是你这个专业比较通用的东西，不明白

没看见8个未知数，只有4个么
赞一下 回复此楼
4楼2011-03-16 00:28:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wzhr

新虫 (小有名气)
恩，r是x的函数，这样的话，也可以说是4个未知数
赞一下 回复此楼
5楼2011-03-17 11:28:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wzhr 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》 +8 scalable 2026-08-24 8/400 2026-08-25 12:52 by jnhyjjm
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[基金申请] 2026国自然函评费到账 +19 羊腰板 2026-08-21 22/1100 2026-08-25 10:38 by popular289
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +5 process2012 2026-08-23 7/350 2026-08-25 09:42 by huixian257
[基金申请] 我面上完蛋了 +13 且听虎啸 2026-08-20 14/700 2026-08-25 09:10 by mrkang
[基金申请] 范进中举一文的中心思想 +7 炎黄贵胄 2026-08-22 8/400 2026-08-25 08:48 by ZJTJZ
[基金申请] 只有每年这种时候来逛逛小木虫 +25 yaoyewhu2008 2026-08-20 27/1350 2026-08-25 08:01 by Equinoxhua
[基金申请] 明天应该可查了！？ +5 chengyan1220 2026-08-23 5/250 2026-08-24 23:28 by 我4大白菜
[基金申请] 能否退出参与的面上项目解除限项 +21 koalala 2026-08-24 24/1200 2026-08-24 19:25 by 家与远方
[基金申请] 让我中一个面上吧！ +13 大萍1987 2026-08-20 16/800 2026-08-24 10:23 by 太傻了
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 17/850 2026-08-24 10:20 by echo8914667
[基金申请] 科研孤儿太难了 +18 我4大白菜 2026-08-20 19/950 2026-08-24 09:47 by 凯恩广盛催化分�
[基金申请] 今天放榜吗？ +15 布布和一二 2026-08-19 16/800 2026-08-23 09:55 by 张春生
[基金申请] 时间戳今天，20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了？ +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 应该是下周三26日公布了吧？ +4 哈哈蛤？ 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
[基金申请] 重要消息，中午系统在维护 +11 yuleib84 2026-08-18 12/600 2026-08-20 11:09 by xskun
信息提示
关闭
请填处理意见
关闭