fsolve() only returns the trivial solution

Question

0 个投票

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 个评论
显示无隐藏无

dpb 2023-8-9

在 MATLAB Online 中打开

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

Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.

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)

0.2699 0.2817

disp(fun(L))

1.0e-05 * 0.8560 -0.2957

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?

João Lazari 2023-8-9

在 MATLAB Online 中打开

I'm looking to construct the equilibrium curve from the adjusted parameters A and B. The curve shuld look something like the figure, as it was used to adjust the parameters.

The bounds for x are 0 and 1, as it is a mole fraction. The should some 1 if in the same phase, but here we have fractions for one component in two different phases, that are in equilibrium. For each temperature, in the figure, we have a liquid mixture, in blue (alpha phase) in equilibrium with another liquid mixture, in red (beta phase).