Why wont my loop populate pre-allocated matrix?

1 visualización (últimos 30 días)
Sam Jones
Sam Jones el 24 de Mayo de 2022
Comentada: Sam Jones el 24 de Mayo de 2022
I am attempting to solve the Redlich-Kwong cubic equation of state to determine the compressibility factor, Z, for 25 different values of T and P. When I run the for loop shown below the variable "Zliq" outputs as a 3x1 complex double. I need it to be a 3x25 complex double to show all iterations. Any help much appreciated.
clc
clear all
close all
syms V Z
T=[277.15 287.15 297.15 307.15 317.15 327.15 337.15 347.15 357.15 367.15 381.15 401.15 421.15 441.15 461.15 481.15 501.15 521.15 541.15 561.15 581.15 601.15 621.15 637.15 647.3];%K
P=[0.813 1.597 2.982 5.318 9.1 15 23.91 36.96 55.57 81.46 133.90 254.35 451.01 754.45 1201.0 1832.6 2696.5 3844.9 5335.5 7231.5 9603.6 12534.8 16133.1 19596.1 22120.0]*10^3; %Converted to Pa
% Gas constant
R=8.314; %m3 Pa/K mol
% Water properties
Tc=647.3; %K
Pc=22120*10^3; %Pa
% RK equation constants
Pr=P/Pc;
Tr=T/Tc;
a=(27*(R^2)*(Tc^2))/(64*Pc);
b=0.08664*(R*Tc)/Pc;
%Compressibility factor Z and specific volume
Zliq=zeros(3,25);
for i=1:25
beta(i)=0.08664*Pr(i)/Tr(i);
q(i)=(0.42748*(Tr(i).^-0.5))/(0.08664*Tr(i));
eq0=Z == beta(i) + Z*(Z+beta(i))*((1+beta(i)-Z)/(q(i)*beta(i)));
Zliq = double(solve(eq0,Z));
end

Respuesta aceptada

KSSV
KSSV el 24 de Mayo de 2022
clc
clear all
close all
syms V Z
T=[277.15 287.15 297.15 307.15 317.15 327.15 337.15 347.15 357.15 367.15 381.15 401.15 421.15 441.15 461.15 481.15 501.15 521.15 541.15 561.15 581.15 601.15 621.15 637.15 647.3];%K
P=[0.813 1.597 2.982 5.318 9.1 15 23.91 36.96 55.57 81.46 133.90 254.35 451.01 754.45 1201.0 1832.6 2696.5 3844.9 5335.5 7231.5 9603.6 12534.8 16133.1 19596.1 22120.0]*10^3; %Converted to Pa
% Gas constant
R=8.314; %m3 Pa/K mol
% Water properties
Tc=647.3; %K
Pc=22120*10^3; %Pa
% RK equation constants
Pr=P/Pc;
Tr=T/Tc;
a=(27*(R^2)*(Tc^2))/(64*Pc);
b=0.08664*(R*Tc)/Pc;
%Compressibility factor Z and specific volume
Zliq=zeros(3,25);
for i=1:25
beta=0.08664*Pr(i)/Tr(i);
q=(0.42748*(Tr(i).^-0.5))/(0.08664*Tr(i));
eq0= Z == beta + Z*(Z+beta)*((1+beta-Z)/(q*beta));
Zliq(:,i) = double(solve(eq0,Z));
end

Más respuestas (0)

Categorías

Más información sobre Function Creation en Help Center y File Exchange.

Productos


Versión

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by