Matrix dimensions must agree. Error Message appears.

1 次查看(过去 30 天)
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
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
Abdullahi Geedi 2019-5-13
I've uploaded it. Though im not sure If its the function as it runs fine with gas mixture of 3 components.

请先登录,再进行评论。

采纳的回答

Dennis
Dennis 2019-5-13
In your function PTFLASH_VLE_PRVDW you define:
fL=zeros(1,c);
fV=[1,1,1];
next you try to substract fV - fL, which will only work if fL has 3 elements.
fV=ones(1,c)
might solve your problem.
You might also just remove that condition from your while loop, since if fL can only contain zeros und Fv can only contain ones it will always be true.

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Oil, Gas & Petrochemical 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by