fsolve() only returns the trivial solution

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;

[dpb]

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]

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).