# Matrix dimensions must agree. Error Message appears.

9 views (last 30 days)
Abdullahi Geedi on 13 May 2019
Answered: Dennis on 13 May 2019
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

Jan on 13 May 2019
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 on 13 May 2019
I've uploaded it. Though im not sure If its the function as it runs fine with gas mixture of 3 components.

Dennis on 13 May 2019
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)
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.