simulkade / PVTtool

A Matlab toolbox for the PVT calculation using the cubic equations of state (EOS)
24 stars 17 forks source link

Stability Testing #2

Open qingqi-code opened 2 years ago

qingqi-code commented 2 years ago

Dear Ali,

I am very grateful to study PVTtool. Here I have a question. When I run code of the 'testcase2':

% test case
% liquid-two-phase phase transition
clc; clear;
% Define the components and load pure physsical properties
[component, comp_flag] = addComponents({'CH4', 'C2H6', 'C4H10', ...
'C6H14', 'C10H22', 'C15H32'});
% Define the thermodynamic models
T0 = 300; % [K]
p0 = 100e5; % [Pa]
thermo1 = addThermo();
thermo1.EOS = @PREOS;
mixture1 = addMixture(component, T0, p0);
% Define flash options
options.accuracy = 1e-7;
options.iteration = 100;
% Negative flash
[vapor_y, liquid_x, vapor_frac] = ...
        vleflashnegative(mixture1, thermo1, options)
%mixture1.temperature=T_range(i);
[liquid_z, vapor_z, fugacity, HR] = PREOS(mixture1, thermo1)    

[s1, sl, sv] = stabilityTest(mixture1, thermo1)

I found that the result is:

vapor_y =

    0.8653    0.1097    0.0198    0.0048    0.0004    0.0000

liquid_x =

    0.4053    0.1472    0.1165    0.1114    0.1099    0.1097

vapor_frac =

   -0.5188

liquid_z =

    0.5924

vapor_z =

    0.5924

fugacity =

    1.9419    0.3795    0.0352    0.0037    0.0001    0.0000

HR =

  -3.0503e+04

s1 =

     1     1

sl =

    0.9910

sv =

    1.0000

The problem is

s1 =

     1     1

that means it converged to trivial solution. If so, the system is single-phase but the simulation gives me two phase data. I think this code did not implement stability analysis to flash calculation directly. If it is problem, I would be happy to help modify this model.

Thank you so much.

simulkade commented 2 years ago

Hi @qingqi-code apologies for the late reply. You are right; i have not incorporated the stability test in the flash calculation function. I've been looking for some time to pay some attentions to this package, which is becoming more and more hard to find :-) You are more than welcome to make changes and improve the code.