Problem with my gamma_water values not giving me the correct pressure using fsolve

Hi,

I'm having trouble with my "gamma_water" variable not giving me the correct P (PressureNonIdeal) value. I keep getting this message : No solution found.

fsolve stopped because the

last step was ineffective. However, the vector of function

values is not near zero, as measured by the value of the

function tolerance.

Here is my current code :

DelH = -1* (-917.475 - (-600.351 + -247.184)) * 1000;
DelS = -1*(172.10 - (76.87 + 232.70));
DelV = -1*(((24.630 - (11.248))*1E-6) *1E5);
Grxn = 0; %because at equilibrium
R = 8.314;


T=948.15; %675 C
x0 = 2000;
fun = @DeltaGcalc;

%Ideal gas
gamma_water = 1
PressureIdeal=fsolve(@(P) fun(DelH,DelS,DelV,T,P,R,gamma_water),x0)

%Non-ideal gas
Ta = 873.15 ; ga = 0.5203; %For a constant pressure of 2000bar
Tb = 1273.15 ; gb = 0.7939;
gamma_water = ((gb-ga)/(Tb-Ta))*T;
PressureNonIdeal=fsolve(@(P) fun(DelH,DelS,DelV,T,P,R,gamma_water),x0)


function Grxn = DeltaGcalc(DelH,DelS,DelV,T,P,R,gamma_water)
Grxn = DelH - T*DelS + DelV*(P-1) + R*T*log(gamma_water*P); %because both the solids are pure, their activity is 1. the activity of the fluid only counts.
end

According to the question, «The pressure (PressureNonIdeal) for the value calculated with the fugacity coefficient (gamma_water) will be between 2 and 10 kbar (200 to 1 000 MPa). Note that you will have to interpolate between the values in the tables, both in terms of pressure and temperature; between any two values in the table at either constant pressure or at constant temperature you can assume that the variations are linear.»

I have attached the table to this description

Answer 1

[Torsten]

Look at the curves and see what's happening. Your first equation has two solutions while your second has none:

DelH = -1* (-917.475 - (-600.351 + -247.184)) * 1000;
DelS = -1*(172.10 - (76.87 + 232.70));
DelV = -1*(((24.630 - (11.248))*1E-6) *1E5);
Grxn = 0; %because at equilibrium
R = 8.314;


T=948.15; %675 C
x0 = 2000;
fun = @DeltaGcalc;

%Ideal gas
gamma_water = 1;
PressureIdeal=fsolve(@(P) fun(DelH,DelS,DelV,T,P,R,gamma_water),x0)

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.

PressureIdeal = 4.8188e+03

x=4e3:10:8e3;
figure()
plot(x,fun(DelH,DelS,DelV,T,x,R,gamma_water))

%Non-ideal gas
Ta = 873.15 ; ga = 0.5203; %For a constant pressure of 2000bar
Tb = 1273.15 ; gb = 0.7939;
gamma_water = ((gb-ga)/(Tb-Ta))*T;
PressureNonIdeal=fsolve(@(P) fun(DelH,DelS,DelV,T,P,R,gamma_water),x0)

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

PressureNonIdeal = 5.8907e+03

x=0.1:100:1e4;
figure()
plot(x,fun(DelH,DelS,DelV,T,x,R,gamma_water))

function Grxn = DeltaGcalc(DelH,DelS,DelV,T,P,R,gamma_water)
Grxn = DelH - T*DelS + DelV*(P-1) + R*T*log(gamma_water*P); %because both the solids are pure, their activity is 1. the activity of the fluid only counts.
end

[Sam Chak]

Hi @Mathys

The gamma water value needs to be greater or equal to 0.981292151253147 to have solutions.

format long

DelH = -1* (-917.475 - (-600.351 + -247.184)) * 1000;
DelS = -1*(172.10 - (76.87 + 232.70));
DelV = -1*(((24.630 - (11.248))*1E-6) *1E5);
Grxn = 0; %because at equilibrium
R = 8.314;

T=948.15; %675 C
x0 = 2000;
fun = @DeltaGcalc;

%% Non-ideal gas
Ta  = 873.15; 
ga  = 0.5203; %For a constant pressure of 2000bar
Tb  = 1273.15; 
gb  = 0.7939;
gamma_water = ((gb-ga)/(Tb-Ta))*T + 0.332757551253147

gamma_water = 
   0.981292151253147

% PressureNonIdeal=fsolve(@(P) fun(DelH,DelS,DelV,T,P,R,gamma_water),x0)
x   = 5800:0.00001:5900;
figure()
plot(x, fun(DelH, DelS, DelV, T, x, R, gamma_water)), xlim([5890.687 5890.6895])
grid on, grid minor
yline(0, '--')

function Grxn = DeltaGcalc(DelH,DelS,DelV,T,P,R,gamma_water)
    Grxn = DelH - T*DelS + DelV*(P-1) + R*T*log(gamma_water*P); 
    % because both the solids are pure, their activity is 1. the activity of the fluid only counts.
end

[Mathys]

Good day @Sam Chak,

I have tried sumitting answers on MatLab Grader with every ga and gb on my chart that produce a gamma_water greater or equal to 0.981292151253147, so I'm assuming that there might be an issue with my approach on the problem at that point. Thanks for the help, though!

[Mathys]

Also, @Torsten, do you mind explaining to me why the curves are so different, even though it's the exact same equation twice, just with a different gamma_water?

[Sam Chak]

Hi @Mathys, for the non-ideal gas, with the rest of the parameters unchanged, you just need to adjust gb so that .

[Torsten]

Also, @Torsten, do you mind explaining to me why the curves are so different, even though it's the exact same equation twice, just with a different gamma_water?

Not that much different - only translated by R*T*log(gamma_water) in y-direction.

DelH = -1* (-917.475 - (-600.351 + -247.184)) * 1000;
DelS = -1*(172.10 - (76.87 + 232.70));
DelV = -1*(((24.630 - (11.248))*1E-6) *1E5);
Grxn = 0; %because at equilibrium
R = 8.314;


T=948.15; %675 C
x0 = 2000;
fun = @DeltaGcalc;

%Ideal gas
gamma_water = 1;
PressureIdeal=fsolve(@(P) fun(DelH,DelS,DelV,T,P,R,gamma_water),x0)

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.

PressureIdeal = 4.8188e+03

x=0.1:100:1e4;
hold on
plot(x,fun(DelH,DelS,DelV,T,x,R,gamma_water),'b')

%Non-ideal gas
Ta = 873.15 ; ga = 0.5203; %For a constant pressure of 2000bar
Tb = 1273.15 ; gb = 0.7939;
gamma_water = ((gb-ga)/(Tb-Ta))*T;
PressureNonIdeal=fsolve(@(P) fun(DelH,DelS,DelV,T,P,R,gamma_water),x0)

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

PressureNonIdeal = 5.8907e+03

x=0.1:100:1e4;
plot(x,fun(DelH,DelS,DelV,T,x,R,gamma_water),'r')
hold off

function Grxn = DeltaGcalc(DelH,DelS,DelV,T,P,R,gamma_water)
  Grxn = DelH - T*DelS + DelV*(P-1) + R*T*log(gamma_water*P); %because both the solids are pure, their activity is 1. the activity of the fluid only counts.
end

[Mathys]

@Sam Chak Well I did just that, except for one thing : since I'm trying to find the fugacity coefficient (gamma_water) using linear interpolation, I have to adjust both gb and ga, as the pressure stays constant. So if I pick a different value of gb, for example according to P=7000 bar, then I also have to adjust ga according to P=7000 bar.

[Torsten]

@Mathys

May it be the case that your equation uses a normalized pressure instead of absolute pressure ? The term DelV*(P-1) for P in the order of 2000 bar looks quite strange to me...

[Mathys]

@Torsten I do not know the pressure, so don't think this is it. In fact, the pressure is what I am looking to obtain. I just know that the PressureNonIdeal will be between 2000 bar and 10000 bar, hence my initial guess of 2000 bar in fsolve. My PressureIdeal value is correct, according to MatLab Grader