fsolve() only returns the trivial solution
2 views (last 30 days)
Show older comments
Hello everyone,
I'm trying to solve a phase equilibrium problem, using the acctivity coefficients through the van Laar model.
I've adjusted the parameters from experimental data and now I'm trying to create the equilibrium curve, but I only get that the molar fraction of component 1 in the alpha phase is equal to the its own molar fraction in beta phase. Wich only should occur in the upper critical point.
I've already tryed to vary the initial guess, but with no results. Does anyone know how to solve this?
Tc = 285.9631; %system critical temperature
T = linspace(270, Tc, 100); %temperature range
A(1) = 175.10385; %A parameter of van Laar
A(2) = -0.591981;
B = 1.1819; %Parameter B of van Laar
x0 = [0.2 0.9]; %composition of component 1 in alpha and beta phase initial guess
for i= 1:length(T)
Ti = T(i);
fun = @(x)eql_Laar(Ti, A, B, x);
L(i,:) = fsolve(fun, x0);
end
function F = eql_Laar(T, a, B, x)
A = a(1) + a(2)*T;
x1_alfa = x(1);
x2_alfa = 1 - x1_alfa;
x1_beta = x(2);
x2_beta = 1 - x1_beta;
gama_1_alfa = exp(A/((1 + ((A*x1_alfa)/(B*x2_alfa)))^2));
gama_2_alfa = exp(B/((1 + ((B*x2_alfa)/(A*x1_alfa)))^2));
gama_1_beta = exp(A/((1 + ((A*x1_beta)/(B*x2_beta)))^2));
gama_2_beta = exp(B/((1 + ((B*x2_beta)/(A*x1_beta)))^2));
F(1) = x1_alfa*gama_1_alfa - x1_beta*gama_1_beta;
F(2) = x2_alfa*gama_2_alfa - x2_beta*gama_2_beta;
2 Comments
dpb
on 9 Aug 2023
Tc = 285.9631; %system critical temperature
T = linspace(270, Tc, 100); %temperature range
A(1) = 175.10385; %A parameter of van Laar
A(2) = -0.591981;
B = 1.1819; %Parameter B of van Laar
x0 = [0.2 0.9]; %composition of component 1 in alpha and beta phase initial guess
for i= 1 %:length(T)
Ti = T(i);
fun = @(x)eql_Laar(Ti, A, B, x);
L(i,:) = fsolve(fun, x0);
end
xx=x0;
for i=1:5
F(i,:)=fun(x0);
x0=x0+[0.1 -0.1];
xx=[xx;x0];
end
xx(end,:)=[];
plot(xx(:,1),F(:,1),xx(:,2),F(:,2))
legend('F1','F2','location','northwest')
xlabel('X'), ylabel('F')
disp(L)
disp(fun(L))
function F = eql_Laar(T, a, B, x)
A = a(1) + a(2)*T;
x1_alfa = x(1);
x2_alfa = 1 - x1_alfa;
x1_beta = x(2);
x2_beta = 1 - x1_beta;
gama_1_alfa = exp(A/((1 + ((A*x1_alfa)/(B*x2_alfa)))^2));
gama_2_alfa = exp(B/((1 + ((B*x2_alfa)/(A*x1_alfa)))^2));
gama_1_beta = exp(A/((1 + ((A*x1_beta)/(B*x2_beta)))^2));
gama_2_beta = exp(B/((1 + ((B*x2_beta)/(A*x1_beta)))^2));
F(1) = x1_alfa*gama_1_alfa - x1_beta*gama_1_beta;
F(2) = x2_alfa*gama_2_alfa - x2_beta*gama_2_beta;
end
So, what are you looking for here? The intersection point of F1, F2?
What are bounds for x? If those are fractions, shouldn't they sum to 1.0, or is this some other measure?
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!