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

[Jan]

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]

I've uploaded it. Though im not sure If its the function as it runs fine with gas mixture of 3 components.

Answer 1

[Dennis]

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.