Plotting Using a function with varying input : Error Matrix dimensions must agree.
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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
8 Commenti
Risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!