search

Cantera / cantera

Chemical kinetics, thermodynamics, and transport tool suite
https://cantera.org
Other
580 stars 341 forks source link

Multi-Phase Equilibria #160

Open GoogleCodeExporter opened 9 years ago

GoogleCodeExporter commented 9 years ago 
What steps will reproduce the problem?
1. Run a multi-phase equilibrium run with all constituents initially in the gas 
phase. (used Pb, O, C from nasa gas and PbO(yw), PbO(rd), PbCO3(s), Pb(s), 
Pb(l), PbO(l), and graphite from nasa condensed)
2. Re-run with constituents in a random mix of condensed and gas phases
3. Compare results

What is the expected output? What do you see instead?
They should yield the same result. Instead the equilibrium run where I start 
with a random mixture of phases causes the condensed phase species not to 
equilibrate with the gas. Equilibrating in this manner seems to always add to 
the condensed phases, but never takes away. 

If I instead start from the gas phase I always get more realistic results. How 
come the equilibrate does not re-arrange all phases?

What version of the product are you using? On what operating system?
Cantera 2.0.2, and python 2.7 on windows 7 x64.

Please provide any additional information below.

Original issue reported on code.google.com by jdemut...@gmail.com on 15 Apr 2013 at 8:47

GoogleCodeExporter commented 9 years ago 
Can you provide an input file and script that demonstrates the observed 
behavior?

Original comment by yarmond on 15 Apr 2013 at 9:43

GoogleCodeExporter commented 9 years ago 
Sure. Here is a simplified file with the two cases (you'll need to comment out 
the cases one at a time. See lines 64-78)

Original comment by jdemut...@gmail.com on 15 Apr 2013 at 11:11

Attachments:

testProblemCantera.m:

close all
clear all
clc

dbstop if error
%% Constants
ag          = 6.02214E+23;      %atoms/mole Avogadro's Number
ec          = 1.60217646e-19;  %J/eV       Electron Charge
kb          = 1.3806503e-23;   %J/K        Boltzman Constant
hp          = 6.626068e-34;    %J-s        Plank's Constant
cl          = 2.99792458e8;    %m/s        Speed of Light
Ru          = kb*ag;           %J/mol-K    Universal Gas constant

%% Initialize Chemistry Data
RootFile = 'Pbcompoundtest.cti';
gas     = IdealGasMix(RootFile,'Pbcompoundtest');

H2OL    = importPhase(RootFile,'H2O_L');
PbL     = importPhase(RootFile,'Pb_L');
PbOL    = importPhase(RootFile,'PbO_L');
H2Os    = importPhase(RootFile,'H2O_s');
Cgr     = importPhase(RootFile,'C_g_r');
Pbs     = importPhase(RootFile,'Pb_s');
PbOrd   = importPhase(RootFile,'PbO_r_d');
PbOyw   = importPhase(RootFile,'PbO_y_w');
PbO2s   = importPhase(RootFile,'PbO2_s');
Pb3O4s  = importPhase(RootFile,'Pb3O4_s');
Pb2CO4s = importPhase(RootFile,'Pb2CO4_s');
PbC2O4s = importPhase(RootFile,'PbC2O4_s');
PbCO3s  = importPhase(RootFile,'PbCO3_s');
liquid  = {H2OL PbL PbOL};
solid   = {H2Os Cgr Pbs PbOrd PbOyw PbO2s Pb3O4s Pb2CO4s PbC2O4s PbCO3s};

% phase   = {gas};
phase   = [{gas} liquid solid];

% Define the extensive number of moles per phase in the initial mixture
% (This is only valid if there is more than one phase)
CondenseArray(1)                = {1.0};    %Assumes first phase is the gas
CondenseArray(2:numel(phase),1) = {0.0};    %All other phases start with no initial constituents
mix                             = Mixture([phase' CondenseArray]);  %Create Mixture

MM = zeros(nSpecies(phase{1})+2,1);
Ntot = 0;
for iP = 1:numel(phase)
    NPsp = nSpecies(phase{iP}); %Number of species phase{iP}
    if NPsp == 1        
        MM(Ntot+1:Ntot+NPsp,:)  = meanMolarMass(phase{iP});
    else
        MM(Ntot+1:Ntot+NPsp,:)  = molarMasses(phase{iP});
    end
    name(Ntot+1:Ntot+NPsp) = speciesNames(phase{iP});
    for iName = Ntot+1:Ntot+NPsp
        assignin('caller',name{iName},iName)
    end 
    Ntot = Ntot+NPsp;
    NSP(iP) = NPsp;
end
Nsp = sum(Ntot);

M       = zeros(numel(MM),1);           %Initialize Array for element Molar Masses

%% Here is where the big difference is:
%Case 1 (all in the gas phase)... notice how the output of phaseMoles(mix)
%shows that everything is as lead carbonate.
M(O) = 3;
M(Pb) = 1;
M(C) = 1;

%Case 2 (mixed gas and solid phase)... same stoichiometry, but now I have
%tied up some of the lead and oxygen as lead oxide yw (not stable at this
%temperature). Run the equilibrium run and it does not yeild the same
%results. Notice in the output of phaseMoles(mix) we have 3 moles of gas,
%and 1 mole of lead oxide yw... it never equilibrated...
% M(O) = 2;
% M(PbO_y_w) = 1;
% M(C) = 1;

Comp = [name{1} ':' num2str(M(1))];
for iN = 2:Nsp
    Comp = [Comp ', ' name{iN} ':' num2str(M(iN),'%0.15e')];
end

Temp = 400;
Press = 101325;
setSpeciesMoles(mix, Comp);
setTemperature(mix,Temp);
setPressure(mix, Press);

equilibrate(mix,'TP', 1e-20, 1e5);
% equilibrate(mix,'TP');
mix
phaseMoles(mix)
GoogleCodeExporter commented 9 years ago 
Attached is a version of the provided test case adapted for the Python 
interface, which demonstrates most of the same behavior. There are a few 
differences, which I think start to explain what's going on here.

Running the Python version, both of the 'equilibrate' calls result in 
exceptions. In the Matlab version, these exceptions are incorrectly being 
ignored, and the Matlab interface needs to be updated to show these errors.

For the first case, there is a convergence failure, although the resulting 
state of the mixture is apparently ok.

In the second case, the error message is:

    "condensed-phase speciesPbO_y_w is excluded since its thermo properties are not valid at this temperature, but it has non-zero moles in the initial state."

So the solver is aborting without modifying the mixture state. I don't know if 
there is a good fix for the MultiPhaseEquil solver in this case.

In the Python interface, there is another equilibrium solver available, which 
is referred to as the "VCSnonideal" solver in the C++ interface, and can be 
used with the vcs_equilibrate function. Unfortunately, I don't think there's 
any way to access this solver from Matlab.

Original comment by yarmond on 18 Apr 2013 at 10:04

Attachments: testProblem.py:

from Cantera import *

RootFile = 'Pbcompoundtest.cti'
gas = IdealGasMix(RootFile,'Pbcompoundtest')

H2Os    = importPhase(RootFile,'H2O_s')
Cgr     = importPhase(RootFile,'C_g_r')
Pbs     = importPhase(RootFile,'Pb_s')
PbOrd   = importPhase(RootFile,'PbO_r_d')
PbOyw   = importPhase(RootFile,'PbO_y_w')
PbO2s   = importPhase(RootFile,'PbO2_s')
Pb3O4s  = importPhase(RootFile,'Pb3O4_s')
Pb2CO4s = importPhase(RootFile,'Pb2CO4_s')
PbC2O4s = importPhase(RootFile,'PbC2O4_s')
PbCO3s  = importPhase(RootFile,'PbCO3_s')

solid   = [H2Os, Cgr, Pbs, PbOrd, PbOyw, PbO2s, Pb3O4s, Pb2CO4s, PbC2O4s, PbCO3s]

mix = Mixture([gas] + solid)

# (Mostly) good case:
print '+++ Case 1 using MultiPhaseEquil +++'
mix.set(Moles='O:3, Pb:1, C:1', T=400, P=101325)
try:
    mix.equilibrate('TP', err=1e-20, maxsteps=100000, maxiter=100000)
except Exception as e:
    print e
print mix.phaseMoles()

# Bad case:
print '+++ Case 2 using MultiPhaseEquil +++'
mix.set(Moles='O:2, PbO_y_w:1, C:1', T=400, P=101325)
try:
    mix.equilibrate('TP', err=1e-20, maxsteps=100000, maxiter=100000)
except Exception as e:
    print e
print mix.phaseMoles()

# But the VCS solver appears to work:
print '+++ Case 2 using VCS_nonideal +++'
mix.set(Moles='O:2, PbO_y_w:1, C:1', T=400, P=101325)
mix.vcs_equilibrate('TP', rtol=1e-20, maxsteps=10000, maxiter=10000)
print mix.phaseMoles()
speth commented 2 months ago

I think there is no longer a Matlab-specific issue here. In the "experimental" Matlab toolbox, we now use the clib mix_equilibrate which uses the "auto" solver option and will start by trying the VCS solver and then fall back on the MultiPhaseEquil/Gibbs solver if that fails.

speth commented 2 months ago

The example here seems to have bit-rotted a bit (still written for the "old" Python toolbox and Python 2.7) and I think we've lost the original input file (the link is to a file that was hosted on cantera.org at some point but I can no longer find it on the server). I was able to partially recreate it -- some of the species aren't in nasa_condensed.yaml, even though the original post claims they were. This new version demonstrates the same failure mode as before for the Gibbs solver, and exhibits an error I've never seen before for the VCS solver.

Using this input file: pbcompoundtest.yaml.txt and Cantera 3.1.0a2 (d37a76bf6):

import cantera as ct

RootFile = 'pbcompoundtest.yaml'
gas = ct.Solution(RootFile,'Pbcompoundtest')

H2Os    = ct.Solution(RootFile,'H2O_s')
Cgr     = ct.Solution(RootFile,'C_g_r')
Pbs     = ct.Solution(RootFile,'Pb_s')
PbOrd   = ct.Solution(RootFile,'PbO_r_d')
PbOyw   = ct.Solution(RootFile,'PbO_y_w')
PbO2s   = ct.Solution(RootFile,'PbO2_s')
Pb3O4s  = ct.Solution(RootFile,'Pb3O4_s')

solid   = [H2Os, Cgr, Pbs, PbOrd, PbOyw, PbO2s, Pb3O4s]

mix = ct.Mixture([gas] + solid)

# (Mostly) good case:
print('+++ Case 1 using MultiPhaseEquil +++')
mix.T = 400
mix.P = ct.one_atm
mix.species_moles = 'O:3, Pb:1, C:1'
try:
    mix.equilibrate('TP', solver='gibbs', rtol=1e-10, max_steps=100000, max_iter=100000)
except Exception as e:
    print(e)
print(mix.phase_moles())

# Bad case:
print('+++ Case 2 using MultiPhaseEquil +++')
mix.T = 400
mix.P = ct.one_atm
mix.species_moles = 'O:2, PbO(yw):1, C:1'
try:
    mix.equilibrate('TP', solver='gibbs', rtol=1e-10, max_steps=100000, max_iter=100000)
except Exception as e:
    print(e)
print(mix.phase_moles())

# But the VCS solver appears to work:
print('+++ Case 2 using VCS_nonideal +++')
mix.T = 400
mix.P = ct.one_atm
mix.species_moles = 'O:2, PbO(yw):1, C:1'
try:
    mix.equilibrate('TP', solver='vcs', rtol=1e-10, max_steps=10000, max_iter=10000)
except Exception as e:
    print(e)
print(mix.phase_moles())

produces the output:

+++ Case 1 using MultiPhaseEquil +++
[0.999999999999979, 0.0, 1.0000000000307008, 0.0, 0.0, 0.0, 1.0000000000000104, 0.0]
+++ Case 2 using MultiPhaseEquil +++

*******************************************************************************
CanteraError thrown by MultiPhaseEquil::MultiPhaseEquil:
condensed-phase speciesPbO(yw) is excluded since its thermo properties are 
not valid at this temperature, but it has non-zero moles in the initial state.
*******************************************************************************

[3.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]
+++ Case 2 using VCS_nonideal +++

*******************************************************************************
CanteraError thrown by VCS_SOLVE::vcs_popPhaseID:
Index out of bounds due to logic error.
*******************************************************************************

[3.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]