Matrix dimensions must agree. Error Message appears.
显示 更早的评论
Hi
I'm running a function to calculate fugacity and other ouput. The function works some of the time. Like if i run it with the input below:
Which is a gas mixture of 3 components-.
Pc=[4.59E6 3.36E6 2.47E6]
Tc=[190.6 469.7 568.7]
w=[0.011 0.251 0.396]
kappa=zeros(3)
eta=zeros(3)
P=0.1E6
T=300
z=[0.85 0.09 0.06]
Ko=[332 0.75 0.03]
But when i try to run a gas mixture of 9 components, with this input:
Pc=[4.599E6 4.872E6 4.248E6 3.796E6 3.37E6 3.025E6 2.49E6 1.817E6 1.401E6];
Tc=[190.6 305.3 369.8 425.1 469.7 507.6 568.7 658.1 722];
T=320;
Omega=[0.012 0.1 0.152 0.2 0.252 0.301 0.346 0.577 0.721];
P=0.2E6;
kappa=zeros(9);
eta=zeros(9);
ki=(Pc/P.*exp(5.37.*(1+Omega).*(1-Tc/T)));
Ko=ki;
z= [0.7 0.08 0.07 0.03 0.01 0.02 0.04 0.02 0.03];
this error message Appears :
Matrix dimensions must agree.
Error in PTFLASH_VLE_PRVDW (line 28)
while max(abs((fV-fL)./fV))>0.00001 && iter<1000 && psi_0*psi_1<
end
The complete Code is below with different subprograms.
The function:
function [x,y,alphaV,fL,fV] = PTFLASH_VLE_PRVDW(Pc, Tc, Omega, kappa, eta, P, T, z, Ko)
%The funtion PTFLASH_VLE_PRVDW makes a PT-FLASH calculation for a mixture
%at given P, T and overall composition (z) using PR-CEoS and VDW MIXING RULES.
%This function requires MIXRULES_VDW and FUGACITY_INMIX_PRVDW funtions to run.
%All input/output data are expressed in SI.
%P[Pa], T[K], w[dimensionless], V[m^3/mol], Z[dimensionless]
%giving a firt value for c
c=length(z);
%firt verification for alphaV between 0 and 1
psi_0=sum(z.*(Ko-1));
psi_1=sum(z.*(Ko-1)./Ko);
if psi_0*psi_1>=0
x=0;
y=0;
alphaV=0;
fL=0;
fV=0;
else
iter=0;
fL=zeros(1,c);
fV=[1,1,1];
K=Ko;
while max(abs((fV-fL)./fV))>0.00001 && iter<1000 && psi_0*psi_1<0
[ alphaV ] = RACHFORDRICE_BISECT( K,z );
x=z./((1-alphaV)+alphaV*K);
y=K.*x;
[fiL,~,fL,~] = fugacity_INMIX_PRVDW( Pc,Tc,Omega,kappa,eta,P,T,x );
[~,fiV,~,fV] = fugacity_INMIX_PRVDW( Pc,Tc,Omega,kappa,eta,P,T,x );
K=fiL./fiV;
psi_0=sum(z.*(K-1));
psi_1=sum(z.*(K-1)./K);
iter=iter+1;
end
end
end
The subprograms:
function [fiL,fiV,fL,fV] = FUGACITY_INMIX_PRVDW(Pc, Tc, w, kappa, eta, P, T, x)
%The funtion FUGACITY_INMIX_PRVDW calculates the fugacity (f) and fugacity
%coefficient (fi) for a fluid IN A MIXTURE at given pressure P, temperature T
%and composition using PR-CEoS and VDW_MIXING RULES. This function requires
%ZMIX_PRVDW funtion to run.
%All input/output data are expressed in SI.
%P[Pa], T[K], f[Pa], fi[dimensionless]
%universal gas constant
R=8.314;
%PR attractive constant and covolume at critical point
b=0.07780*R*Tc./Pc;
k=0.37464+1.54226*w-0.26992*w.^2;
alpha=(1+k.*(1-sqrt(T./Tc))).^2;
a=0.45724*alpha.*(R.^2*(Tc.^2))./Pc;
%calling MIXRULES_VDW function
[am, bm, aa] = MIXRULES_VDW(a, b, kappa, eta, x);
%calculation of Am and Bm
Am=am*P/(R*T)^2;
Bm=bm*P/(R*T);
%calculation of A and B
A=aa*P/(R*T)^2;
B=b*P/(R*T);
%calculation of ZL and ZV using ZMIX_PRVDW function
[Z,V,ZL,ZV,VL,VV] = ZMIX_PRVDW(Pc, Tc, w, kappa, eta, P, T, x);
%giving a firt value for c
c=length(x);
%calculation of fugacity coeffiecient
for i=1:c
c(i)=x*(A(i,:))';
lnfiL(i)=B(i)/Bm*(ZL-1)-log(ZL-Bm)-Am/(2*sqrt(2)*Bm)*(2*c(i)/Am-B(i)/Bm)*log((ZL+(1+sqrt(2))*Bm)/(ZL+(1-sqrt(2))*Bm));
lnfiV(i)=B(i)/Bm*(ZV-1)-log(ZV-Bm)-Am/(2*sqrt(2)*Bm)*(2*c(i)/Am-B(i)/Bm)*log((ZV+(1+sqrt(2))*Bm)/(ZV+(1-sqrt(2))*Bm));
fiL(i)=exp(lnfiL(i));
fiV(i)=exp(lnfiV(i));
end
%calculation of fugacity
fL=P*fiL.*x;
fV=P*fiV.*x;
end
Another subprogram:
function [Z,V,ZL,ZV,VL,VV] = ZMIX_PRVDW(Pc, Tc, w, kappa, eta, P, T, x)
%The funtion ZMIX_PRVDW calculates the compressibility factor z and the specific
%molar volume V for a mixture at given pressure P, temperature T and concentration
%using PR-CEoS and VDW MIXING RULES. This function requires REALROOTS_CARDANO and
% MIXRULES_VDW funtions to run.
%All input/output data are expressed in SI.
%P[Pa], T[K], w[dimensionless], V[m^3/mol], Z[dimensionless]
%universal gas constant
R=8.314;
%PR attractive constant and covolume at critical point
b=0.07780*(R*Tc)./Pc;
k=0.37464+1.54226*w-0.26992*w.^2;
alpha=(1+k.*(1-sqrt(T./Tc))).^2;
a=0.45724*alpha.*(R^2.*(Tc.^2))./Pc;
%calling MIXRULES_VDW function
[am, bm, aa] = MIXRULES_VDW(a, b, kappa, eta, x);
%calculation of Am and Bm
Am=am*P/(R*T)^2;
Bm=bm*P/(R*T);
%calculation of A and B
a2=-1+Bm;
a1=Am-3*Bm^2-2*Bm;
a0=-Am*Bm+Bm^2+Bm^3;
%calculation of Z and V using REALROOTS_CARDANO function
Z=REALROOTS_CARDANO(a2,a1,a0);
V=Z*R*T/P;
%selection of Z and V values for liquid and vapor
if min(V)<b
ZL=max(Z);
ZV=max(Z);
VL=max(V);
VV=max(V);
else
ZL=min(Z);
ZV=max(Z);
VL=min(V);
VV=max(V);
end
end
The last Subprogram:
function [RR]=REALROOTS_CARDANO(a2,a1,a0)
%This function calculates the real roots of a polynomial of degree 3 using
%an algorithm based on the Cardano's formula.
%The polynomial must be written as: x^3 + a2*x^2 + a1*x + a0.
%Input data: a2, a1, a0.
%Output data: array containing the real roots.
%The output array can be 1x1 or 1x3.
Q=(a2^2-3*a1)/9;
R=(2*a2^3-9*a1*a2+27*a0)/54;
if Q==0 && R==0 %Case of single real root with multiplicity 3
RR=-a2/3;
else
if Q^3-R^2>=0
theta=acos(R/(Q^(3/2)));
RR(1)=-2*sqrt(Q)*cos(theta/3)-a2/3;
RR(2)=-2*sqrt(Q)*cos((theta+2*pi)/3)-a2/3;
RR(3)=-2*sqrt(Q)*cos((theta+4*pi)/3)-a2/3;
else
RR=-sign(R)*((sqrt(R^2-Q^3)+abs(R))^(1/3)+Q/(sqrt(R^2-Q^3)+abs(R))^(1/3))-a2/3;
end
end
end
2 个评论
Jan
2019-5-13
It is not enough to see two different inputs. Without knowing the code it is impossible to guess the reason of the problem or to suggest a solution. So please post details about the failing function.
Abdullahi Geedi
2019-5-13
采纳的回答
更多回答(0 个)
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!