matlab 相平衡续贴 lidaxue1987兄请留意

diyifan001 · 发表于 1970-1-1 08:00:00

本帖最后由 diyifan001 于 2011-4-20 18:54 编辑

开这个帖子是续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 · 发表于 1970-1-1 08:00:00

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

小飞侠 · 发表于 1970-1-1 08:00:00

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

lidaxue1987 · 发表于 1970-1-1 08:00:00

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

lidaxue1987 · 发表于 1970-1-1 08:00:00

本帖最后由 lidaxue1987 于 2011-4-19 10:18 编辑

楼主你好
我把你的两个程序合并为一个
最后程序返回的物系里面组分的气液压缩因子
请楼主帮忙测试，我测试没问题
=========================
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 · 发表于 1970-1-1 08:00:00

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

lidaxue1987 · 发表于 1970-1-1 08:00:00

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

diyifan001 · 发表于 1970-1-1 08:00:00

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

lidaxue1987 · 发表于 1970-1-1 08:00:00

本帖最后由 lidaxue1987 于 2011-4-19 10:46 编辑

回复 diyifan001 的帖子

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

diyifan001 · 发表于 1970-1-1 08:00:00

回复小飞侠的帖子

姑娘
我看好你吆

diyifan001 · 发表于 1970-1-1 08:00:00

回复 lidaxue1987 的帖子

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

小飞侠 · 发表于 1970-1-1 08:00:00

回复 diyifan001 的帖子

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

高马跑 · 发表于 1970-1-1 08:00:00

是不是可以换用aspen来算？

lidaxue1987 · 发表于 1970-1-1 08:00:00

回复高马跑的帖子

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

diyifan001 · 发表于 1970-1-1 08:00:00

回复高马跑的帖子

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

高马跑 · 发表于 1970-1-1 08:00:00

回复 lidaxue1987 的帖子

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

高马跑 · 发表于 1970-1-1 08:00:00

回复 diyifan001 的帖子

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

lidaxue1987 · 发表于 1970-1-1 08:00:00

回复高马跑的帖子

我可不是什么大侠呵呵

diyifan001 · 发表于 1970-1-1 08:00:00

回复高马跑的帖子

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

lidaxue1987 · 发表于 1970-1-1 08:00:00

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

matlab 相平衡续贴 lidaxue1987兄请留意

评分

相关帖子

matlab 相平衡 续贴 lidaxue1987兄请留意

评分

相关帖子

matlab 相平衡续贴 lidaxue1987兄请留意