MATLAB: Matrix dimensions must agree. Error Message appears.

errorinputmatrixmatrix array

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

Best Answer

  • 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.