simulkade / PVTtool

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

vleflashnegative.m seems cannot be converged. #5

Open Qiuyan0831 opened 3 months ago

Qiuyan0831 commented 3 months ago

Hi professor,

Hope this message finds you well.

I try to run your examples, but I find the function vleflashnegative.m seems not work. For example, the example 'blackoilmodel':

% test case % liquid-two-phase phase transition clc; clear; % Define the components and load pure physsical properties z=[0.1, 0.1, 0.1, 0.2, 0.3, 0.2] [component, comp_flag] = addComponents({'CH4', 'C2H6', 'C4H10', ... 'C6H14', 'C8H18', 'C10H22'}); % Define the thermodynamic models T = 100+273.15; % [K] T0= 25+273.15; % [K] surface temperature p0= 1e5; % [Pa] surface pressure thermo1 = addThermo(); thermo1.EOS = @PREOS; mixture1 = addMixture(component, T0, p0); mixture1.mole_fraction = z; % Define flash options options.accuracy = 1e-7; options.iteration = 100; % flash for surface properties [vapor_y, liquid_x, vapor_frac_surf] = ... vleflash(mixture1, thermo1, options) [liquid_z, vapor_z, fugacity, HR] = PREOS(mixture1, thermo1); R=8.314; % [J/(mol.K)] mixture1.mole_fraction=liquid_x; [liquid_z_surf, vapor_z, fugacity, HR] = PREOS(mixture1, thermo1); V_oil_surface = liquid_z_surf(1-vapor_frac_surf)RT0/p0; mixture1.mole_fraction=vapor_y; [liquid_z, vapor_z_surf, fugacity, HR] = PREOS(mixture1, thermo1); V_gas_surface = vapor_z_surfvapor_frac_surfRT0/p0;

n_point=50; p_range=linspace(p0, 300e5, n_point); for i=1:n_point p=p_range(i); % flash for reservoir properties mixture1.temperature=T; mixture1.pressure=p; mixture1.mole_fraction = z; [vapor_y, liquid_x, vapor_frac_res] = ... vleflashnegative(mixture1, thermo1, options) ; %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Here! cannot be converged. if vapor_frac_res<=0 vapor_y=z; liquid_x=z; vapor_frac_res=0; end
mixture1.mole_fraction=liquid_x; [liquid_z_res, vapor_z, fugacity, HR] = PREOS(mixture1, thermo1); V_oil_res = liquid_z_res(1-vapor_frac_res)RT/p; mixture1.mole_fraction=vapor_y; [liquid_z, vapor_z_res, fugacity, HR] = PREOS(mixture1, thermo1); V_gas_res = vapor_z_resvapor_frac_resRT/p;

%[s1, sl, sv] = stabilityTest(mixture1, thermo1) Bo(i) = V_oil_res/V_oil_surface Bg(i) = V_gas_res/V_oil_surface Rs(i) = (vapor_frac_surf-vapor_frac_res)vapor_z_surfRT0/p0/V_oil_surface rs1 = (vapor_frac_res-vapor_frac_surf)liquid_z_surfRT0/p0/V_gas_surface; if rs1<0 rs1=0; end rs(i)=rs1; end

Could you tell me the problem? AND could you please teach me what is the different between vleflash and vleflashnegative, and their applicable scene?

simulkade commented 3 months ago

Hi @Qiuyan0831 , The negative flash algorithm is described in a paper by two of the giants of fluid phase equilibria, Whitson and Michelsen. My implementation is a very basic one, used in my thesis several years ago, and can often fail. At the moment, I am not working on this code so I'm afraid you are on your own. If you come up with a solution, feel free to make a pull request.