Plotting Using a function with varying input : Error Matrix dimensions must agree.

Hi

I've function to calculate vapour fraction( alphaV) which uses Pressure as an input. I would like to see how the vapour fraction varies with changes in pressure particularly in the range of

P1=4.83E2
P2=2.963E7

My input is :

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=4.83E2:500:29620000;
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];

[x,y,alphaV,fL,fV] = PTFLASH_VLE_PRVDW(Pc, Tc, Omega, kappa, eta, P, T, z, Ko)

But the problem i'm currently having is i get an error message:

Error using  / 
:Matrix dimensions must agree.

Error in Plotting (line 10)
ki=(Pc/P.*exp(5.37.*(1+Omega).*(1-Tc/T)));

My 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=ones(1,c);
    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

My subprogram for fugacity:

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

My Subprogram for Compressibility

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

Subprogram for ALpha V

function [ alphaV ] = RACHFORDRICE_BISECT( K,z )
%This function applies the bisection method to the Rachford-Rice equation
%for alphaV in (0,1).
%Input/output are on molar basis. 

a=0;
b=1;
%psi_0=sum(z.*(K-1))
%psi_1=sum((z.*(K-1))./K)

while b-a>0.000001
    alphaV=(a+b)/2;
    psi_alphaV=sum((z.*(K-1))./((1-alphaV)+alphaV*K));
    psi_b=sum((z.*(K-1))./((1-b)+b*K));
    if psi_alphaV*psi_b<0
        a=alphaV;
    else
        b=alphaV;
    end
end

end

Last subprogram for real roots:

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

[madhan ravi]

Didn't verify your code but try putting ./ at that line and see what happens.

[Abdullahi Geedi]

Same except no the dot is also part of the error message.

Error using  ./ 
Matrix dimensions must agree.

Error in Plotting (line 10)
ki=(Pc./P.*exp(5.37.*(1+Omega).*(1-Tc./T)));

[Dennis]

You need to rethink that.

First

P=29620000E7:500:4.83E2;

is empty because 29620000E7 is bigger than 4.83E2.

You can fix this by changing 500 to -500. You will then notice that the array needs more than 4 Terabyte of RAM. Most likely you mean 29620000 or 2.9620000E7.

Last Pc only has 10ish (did not count) elements, so P only should have 10 aswell and i do not see how.

[Abdullahi Geedi]

The Range and data. Was a mistake i saw and corrected it now.

P=4.83E2:500:29620000;

That is more the attempt i wanted to make. But like you said i only have ten components, there can only have 10 Pc.

I was thinking of using Ko as constant even though i change the pressure. Is that possible to plot it then.

[Adam Danz]

What do you expect to be the output of this: Pc/P ?

A vector? A scalar? A matrix?

Pc is a vector with 9 values. P is a vector with 59240 values.

Pc'./P  %this produces a [9-by-59240] matrix

Is that what you're looking for?

[Dennis]

I think you could use a loop:

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];
kappa=zeros(9);
eta=zeros(9);
z= [0.7 0.08 0.07 0.03 0.01 0.02 0.04 0.02 0.03];
alphaV=zeros(59240,1); % <- might need to change this!
for P=4.83E2:500:29620000
ki=(Pc/P.*exp(5.37.*(1+Omega).*(1-Tc/T)));
[x,y,alphaV(P),fL,fV] = PTFLASH_VLE_PRVDW(Pc, Tc, Omega, kappa, eta, P, T, z, ki);
end

I did not check if alphaV is a scalar/vector/matrix, you might need to change the preallocation of alphaV. The version above overwrites x,y,fL,fV in every iteration - if you need those values you need to index them aswell.

[Abdullahi Geedi]

I want to know how pressure affects vapour fraction. I've done it the long way of using one P at a time, using 22 Pressure, thereby generating 22 AlphaV. But i need more dots on my plot.

I fgured i could do thing by using a range with step-size.

If possible i would like the matlab script to generate 9 Ki values for every pressure, which allows me then to get vapour fraction value as well. Then i can plot AlphaV against P

[Abdullahi Geedi]

Didn't work. tried the revised method.