matlab 相平衡 续贴 lidaxue1987兄请留意

diyifan001

开这个帖子是续http://meng.horse/thread-28834-1-1.html这个帖子

首先，选PR方程（SRK等立方型状态方程等思路一样）编程语言matlab（最后用.net 做个界面）
第一：函数getAB，计算PR方程整理成多项式后的AB值，顺便算出a、 b 阿尔法（用al表示）K的值

function y=getAB(p,t,x)
R=8.3145;
Pr=p./x(:,1);
Tr=t./x(:,2);
a=(0.45724.*R.*R.*x(:,2).*x(:,2))./x(:,1);
b=(0.0778.*R.*x(:,2))./x(:,1);
k=0.37464+1.54226.*x(:,3)-0.26992.*x(:,3).*x(:,3);
al=(1+k.*(1-sqrt(Tr))).*(1+k.*(1-sqrt(Tr)));
A=(a.*al.*p)./(R.*R.*t.*t);
B=b.*p./(R.*t);
y=[A B a b al k];

第二：求压缩性系数方法getZ，就是通过AB的值求出三次方程的三个根

function y=getZ(A,B)
temp=size(A);
x1=ones(temp(1),1);
x2=B-1;
x3=A-2.*B-3.*B.*B;
x4=B.*B.*B-A.*B+B.*B;
X=[x1,x2,x3,x4];
y=roots(X);

第三 ：求混合物的压缩性系数，因为主要计算混合物的相平衡，因为昨晚跟个兄弟讨论hysys extension的东西，咬了会文档，所以写的有点乱，有空整理一下，慢慢来，一步一步来，免得出错
function y=HunHeZ(p,t,X,Kij)
R=8.31451;
AiBi=getAB(p,t,X(:,1:3));
% b
bi=AiBi(:,4);
tempb=bi.*X(:,4);
b=sum(tempb);
%zi*zj
HLtemp=X(:,4);
n=length(HLtemp);
HL1=repmat(HLtemp,1,n);
HL2=HL1';
Mix1=HL1.*HL2;
% a*al
a=AiBi(:,3);
al=AiBi(:,5);
aal1=a.*al;
aal2=repmat(aal1,1,n);
aal3=aal2';
Mix2=sqrt(aal2.*aal3);
%1-kij
Mix3=1-Kij;
Mixaal=Mix1.*Mix2.*Mix3;
Mixaal=sum(Mixaal);
Mixaal=sum(Mixaal);

A=(Mixaal.*p)./(R.*R.*t.*t);
B=b.*p./(R.*t);
Ztemp=getZ(A,B);
for index=1:3
if imag(Ztemp(index))~=0;
Ztemp(index)=0;
end
end
Z=Ztemp(find(Ztemp));
y=Z;


把这些文件黏贴到文本文档，改名和函数名同名的m文件，就是 "函数名".m   ,都存到一个文件夹里以备后用。

chenhaijunsjy

急等看看,学学{:1106_362:}

小飞侠

学习一下，楼主真有钻研精神啊……

lidaxue1987

没想到 楼主速度这么快～呵呵 先回复看看

lidaxue1987

楼主你好
我把你的两个程序合并为一个
最后程序返回的物系里面组分的气液压缩因子
请楼主帮忙测试，我测试没问题
=========================
function ZZZ=getroots(p,t,x)
R=8.3145;
Pr=p./x(:,1);
Tr=t./x(:,2);
a=(0.45724.*R.*R.*x(:,2).*x(:,2))./x(:,1);
b=(0.0778.*R.*x(:,2))./x(:,1);
k=0.37464+1.54226.*x(:,3)-0.26992.*x(:,3).*x(:,3);
al=(1+k.*(1-sqrt(Tr))).*(1+k.*(1-sqrt(Tr)));
A=(a.*al.*p)./(R.*R.*t.*t);
B=b.*p./(R.*t);
temp=size(A);
x1=ones(temp(1),1);
x2=B-1;
x3=A-2.*B-3.*B.*B;
x4=B.*B.*B-A.*B+B.*B;
X=[x1,x2,x3,x4];
for i=1:temp(1)
Zl(1,i)=real(min(roots(X(i,:))));
Zv(1,i)=real(max(roots(X(i,:))));
end
ZZZ=[Zl;Zv];==================================================
后续组分的逸度系数求取 同样可以添加，返回出来的是逸度系数即可，因为混合物系相平衡逸度用的比较多，呵呵

diyifan001

恩  您的我看了，异曲同工。
但是我建议写程序还是分开写好些，将来找错误好找，可以一个一个一个函数的运行，以前就是文献上的一个式子写错了，搞得半天没弄好。
另外matlab我不太懂接口的用处，.net的C#可以用接口的方式让很多状态方程共享一个算法，不知道matlab里怎么实现？避免代码繁琐。
晚上继续：混合规则算混合体系的压缩性系数，各个组分的逸度系数。

lidaxue1987

至于matlab高级应用  我就不懂了，以前也不是搞这个的  呵呵  见谅！

diyifan001

还有，只返回压缩性系数有个问题，以后迭代算法中判断Z值是否合理需要具体的Z值，进行付给迭代初值也需要Z值（不知道其他算法需不要要，有其他框图或付初值的方法可以提供一下最好）
所以函数分功能块写要好些。不知道您以为如何？

lidaxue1987

回复 diyifan001 的帖子

付给迭代初值也需要的Z值：
对于混合物液相压缩因子 可以设置成 Z0=bm;
bm是按照混合规则计算得到的
气相的可以设置成1,或者0.9
==========================
至于判断是否合理，这个就有点难了，不知道是否可以把压缩因子按照混合规则进行计算，分别得到两相的压缩因子  呵呵
============
函数分块写 没问题的，呵呵  您按照您的思路来

diyifan001

回复 小飞侠 的帖子

姑娘
我看好你吆

diyifan001

回复 lidaxue1987 的帖子

ok
有啥错误您指出来
我是业余程序员，呵呵，还有，有啥好的迭代算法还望不啬赐教

小飞侠

回复 diyifan001 的帖子

呵呵，我刚开始学习这个，本身对语言有恐惧感，觉得很复杂……慢慢学习呢

高马跑

是不是可以换用aspen来算？

lidaxue1987

回复 高马跑 的帖子

aspen算这个简直太小儿科了，问题是：我们想知道 aspen是怎么算的，以及我们自己编写的与aspen的结果又没有大的出入。呵呵

diyifan001

回复 高马跑 的帖子

aspen plus 、hysys都可以算，hysys下面的状态条从黄变绿就算好了。
呵呵

高马跑

回复 lidaxue1987 的帖子

wow 我对这两个都仅仅知道点皮毛，请各位大侠多多指教，现在在用matlab做毕业论文里面的反应器的一个模拟计算，头就大了，正学着怎么用呢

高马跑

回复 diyifan001 的帖子

Aspen还比较了解，另一个还真没接触过，我先学习一下再来玩 嘻嘻

lidaxue1987

回复 高马跑 的帖子

我可不是什么大侠  呵呵

diyifan001

回复 高马跑 的帖子

恩
一起学习，有经验了编个炮网专用的工具箱
呵呵

lidaxue1987

楼主 最近工作进展怎么样啊?呵呵 小弟还在等着好消息呢{:1106_383:}