Newton Raphson multivariate xy equilibrium

1 次查看（过去 30 天）
Juliana Quintana
Juliana Quintana 2022-3-6
编辑： Torsten 2022-3-7
Hello, I am developing a code to find the equilibrium between xy and the behavior of the temperature and each of the components. The 2 variables I have to find are y component and T the temperature. I think the code is good but obviously have some issues since it isn’t returning the right answer. If you could help me I will be really grateful.
clc
clear
close all
%Vector de x
n=100;
x=linspace(0,1,n);
%Vector para guardas las variables
yresp = zeros(n,1);
Tresp = zeros(n,1);
%Inicializacion;
for i = 1: n
iter=0;
itermax=1000;
error=1000;
tol=1e-6;
% Inicialización
m0= [0.01 ; 400]';
% Recorrido del método NR multivariado
while error>tol && iter<itermax
f0= raoult(m0,x(i)) ;
j0=jacobiano(m0,x(i)) ;
delta=-(j0\f0);
xnew=m0+delta;
% Criterio del ciclo
error=norm (delta) ;
% Reemplazo nuevo valor
m0=xnew;
end
yresp(i)= xnew(1);
Tresp(i) = xnew(2);
end
figure (1)
plot(x,Tresp, 'r',yresp,Tresp,'b')
legend ('x' ,'y')
xlabel 'x,y, Benceno'
ylabel 'T (K)'
figure (2)
plot(x,yresp)
xlabel 'x'
ylabel 'y'
legend 'Benceno'
function j = jacobiano(m,x1)
n = length(x1);
h = 1e-9;
for i=1:n
dh = zeros(n,1);
dh(i) = h;
derv = (raoult(m+dh,x1))-raoult(m,x1)/(h);
j(:,i) = derv;
end
end
CONSTANTES
function resp = Antoine(T)
%Benceno
A(1)= 6.01905; B(1) = 1204.637; C(1)=53.081;
%Tolueno
A(2)=6.08426; B(2)= 1347.2; C(2)=53.363;
resp (:,1) = 10^(A(1)-B(1)/(T+C(1)));
resp (:,2) = 10^(A(2)-B(2)/(T+C(2)));
end
function gamma= wilson (T,x)
x1 = x;
x2 = 1-x1;
%Benceno(1)
%Tolueno (2)
v251 = 90.4; v252 = 104.9 ;
vb1=96; vb2=118.2;
delta251=18.8; delta252=18.7;
tb1=80.09; tb2=110.622;
eps12=0.0851; eps21 = -0.0884;
R = 8.314;
z = 2;
beta1 =(vb1-v251)/(tb1-25);
beta2= (vb2-v252)/(tb2-25);
v1= v251+beta1*((T-273.15)-25);
v2= v252+beta2*((T-273.15)-25);
delta1= (v251/v1)*delta251;
delta2= (v252/v2)*delta252;
lambda11 = -(2/z)*v1*(delta1)^2 ;
lambda12 = -(1-eps12)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda21 = -(1-eps21)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda22 = -(2/z)*v2*(delta2)^2 ;
A21= (v1/v2)*exp(-(lambda21-lambda22)/(R*T));
A12=(v2/v1)*exp(-(lambda12-lambda11)/(R*T));
gamma(1) = exp(-log(x1+A12*x2)+(x2*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2))))));
gamma(2) = exp(-log(x2+A21*x1)-(x1*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2))))));
end
function resp = raoult(m,x1)
P = 101.325 ;
y1 = m(1);
T = m(2);
Psat = Antoine(T);
gm = wilson(x1,T);
resp(:,1) = x1.*gm(1).*Psat(1)- y1*P;
resp(:,2) = (1-x1).*gm(2).*Psat(2)- (1-y1)*P;
end

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-3-6
function main
%Vector de x
n=100;
x=linspace(0,1,n);
%Vector para guardas las variables
yresp = zeros(n,1);
Tresp = zeros(n,1);
m0= [0.01 ; 400]';
for i=1:n
fun = @(m0)raoult(m0,x(i)) ;
sol = fsolve(fun,m0);
yresp(i)= sol(1);
Tresp(i)= sol(2);
m0 = sol;
end
figure (1)
plot(x,Tresp, 'r',yresp,Tresp,'b')
legend ('x' ,'y')
xlabel 'x,y, Benceno'
ylabel 'T (K)'
figure (2)
plot(x,yresp)
xlabel 'x'
ylabel 'y'
legend 'Benceno'
end
function resp = Antoine(T)
%Benceno
A(1)= 6.01905; B(1) = 1204.637; C(1)=53.081;
%Tolueno
A(2)=6.08426; B(2)= 1347.2; C(2)=53.363;
resp (:,1) = 10^(A(1)-B(1)/(T+C(1)));
resp (:,2) = 10^(A(2)-B(2)/(T+C(2)));
end
function gamma= wilson (T,x)
x1 = x;
x2 = 1-x1;
%Benceno(1)
%Tolueno (2)
v251 = 90.4; v252 = 104.9 ;
vb1=96; vb2=118.2;
delta251=18.8; delta252=18.7;
tb1=80.09; tb2=110.622;
eps12=0.0851; eps21 = -0.0884;
R = 8.314;
z = 2;
beta1 =(vb1-v251)/(tb1-25);
beta2= (vb2-v252)/(tb2-25);
v1= v251+beta1*((T-273.15)-25);
v2= v252+beta2*((T-273.15)-25);
delta1= (v251/v1)*delta251;
delta2= (v252/v2)*delta252;
lambda11 = -(2/z)*v1*(delta1)^2 ;
lambda12 = -(1-eps12)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda21 = -(1-eps21)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda22 = -(2/z)*v2*(delta2)^2 ;
A21= (v1/v2)*exp(-(lambda21-lambda22)/(R*T));
A12=(v2/v1)*exp(-(lambda12-lambda11)/(R*T));
gamma(1) = exp(-log(x1+A12*x2)+(x2*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2))))));
gamma(2) = exp(-log(x2+A21*x1)-(x1*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2))))));
end
function resp = raoult(m,x1)
P = 101.325 ;
y1 = m(1);
T = m(2);
Psat = Antoine(T);
%gm = wilson(x1,T)
gm(1:2) = 1.0;
resp(:,1) = x1.*gm(1).*Psat(1)- y1*P;
resp(:,2) = (1-x1).*gm(2).*Psat(2)- (1-y1)*P;
end
I think you should check the calculation of the activity coefficients.
  5 个评论
显示 3更早的评论
Torsten
Torsten 2022-3-6
编辑：Torsten 2022-3-7
%Vector de x
n=100;
x=linspace(0,1,n);
%Vector para guardas las variables
yresp = zeros(n,1);
Tresp = zeros(n,1);
%Inicializacion;
for i = 1: n
iter=0;
itermax=1000;
error=1000;
tol=1e-6;
% Inicialización
m0= [0.01 ; 400];
% Recorrido del método NR multivariado
while error>tol && iter<itermax
f0= raoult(m0,x(i)) ;
j0=jacobiano(m0,x(i)) ;
delta=-(j0\f0);
xnew=m0+delta;
% Criterio del ciclo
error=norm (delta) ;
% Reemplazo nuevo valor
m0=xnew;
end
yresp(i) = xnew(1);
Tresp(i) = xnew(2);
end
figure (1)
plot(x,Tresp, 'r',yresp,Tresp,'b')
legend ('x' ,'y')
xlabel 'x,y, Benceno'
ylabel 'T (K)'
figure (2)
plot(x,yresp)
xlabel 'x'
ylabel 'y'
legend 'Benceno'
end
function j = jacobiano(m,x1)
y = raoult(m,x1);
h = 1e-8;
for i=1:numel(m)
m(i) = m(i) + h;
yh = raoult(m,x1);
j(:,i) = (yh-y)/h;
m(i) = m(i) - h;
end
end
CONSTANTES
function resp = Antoine(T)
%Benceno
A(1)= 6.01905; B(1) = 1204.637; C(1)=53.081;
%Tolueno
A(2)=6.08426; B(2)= 1347.2; C(2)=53.363;
resp(1,:) = 10^(A(1)-B(1)/(T+C(1)));
resp(2,:) = 10^(A(2)-B(2)/(T+C(2)));
end
function gamma= wilson (T,x)
x1 = x;
x2 = 1-x1;
%Benceno(1)
%Tolueno (2)
v251 = 90.4; v252 = 104.9 ;
vb1=96; vb2=118.2;
delta251=18.8; delta252=18.7;
tb1=80.09; tb2=110.622;
eps12=0.0851; eps21 = -0.0884;
R = 8.314;
z = 2;
beta1 =(vb1-v251)/(tb1-25);
beta2= (vb2-v252)/(tb2-25);
v1= v251+beta1*((T-273.15)-25);
v2= v252+beta2*((T-273.15)-25);
delta1= (v251/v1)*delta251;
delta2= (v252/v2)*delta252;
lambda11 = -(2/z)*v1*(delta1)^2 ;
lambda12 = -(1-eps12)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda21 = -(1-eps21)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda22 = -(2/z)*v2*(delta2)^2 ;
A21= (v1/v2)*exp(-(lambda21-lambda22)/(R*T));
A12=(v2/v1)*exp(-(lambda12-lambda11)/(R*T));
gamma(1) = exp(-log(x1+A12*x2)+(x2*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2))))));
gamma(2) = exp(-log(x2+A21*x1)-(x1*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2))))));
end
function resp = raoult(m,x1)
P = 101.325 ;
y1 = m(1);
T = m(2);
Psat = Antoine(T);
%gm = wilson(x1,T);
gm(1:2) = 1.0;
resp(1,:) = x1.*gm(1).*Psat(1)- y1*P;
resp(2,:) = (1-x1).*gm(2).*Psat(2)- (1-y1)*P;
end
I used your code for the version with activity coefficients equal to 1.
You will have to replace "wilson" with your corrected version.
Juliana Quintana
Juliana Quintana 2022-3-6
I just made al the changes and it works perfecty. Thank you for your help and time. Amazing work :)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Newton-Raphson Method 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by