test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition: Added mass and energy conservation plots

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterAndIsopropanolEvaporation/Allrun

@@ -33,7 +33,7 @@ set rmargin at screen 0.84
 gasSpecies = '$gasSpecies'
 liquidSpecies = '$liquidSpecies'
 
-set multiplot layout 4,1
+set multiplot layout 6,1
 
 set xlabel "Time (s)"
 
@@ -61,7 +61,10 @@ plot \
     'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/T.gas' w l lc 1 notitle, \
     'postProcessing_none.gas_none.liquid/plot/0/T.liquid' w l lc 2 t 'Liquid', \
     for [g in gasSpecies] for [l in liquidSpecies] \
-    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/T.liquid' w l lc 2 notitle
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/T.liquid' w l lc 2 notitle, \
+    'postProcessing_none.gas_none.liquid/plot/0/Tf.gasAndLiquid' w l lc 3 t 'Interface', \
+    for [g in gasSpecies] for [l in liquidSpecies] \
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/Tf.gasAndLiquid' w l lc 3 notitle
 
 set ytics nomirror
 set y2tics
@@ -92,6 +95,38 @@ plot \
     for [g in gasSpecies] for [l in liquidSpecies] \
     'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/C3H8O.liquid' w l lc 2 notitle
 
+set ytics nomirror
+set y2tics
+set ylabel "Mass (kg/m^3)"
+set y2label "Energy (J/m^3)"
+
+plot \
+    'postProcessing_none.gas_none.liquid/plot/0/dMass.gas' w l lc 1 t 'Gas Mass Change', \
+    for [g in gasSpecies] for [l in liquidSpecies] \
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/dMass.gas' w l lc 1 notitle, \
+    'postProcessing_none.gas_none.liquid/plot/0/dMass.liquid' w l lc 2 t 'Liquid Mass Change', \
+    for [g in gasSpecies] for [l in liquidSpecies] \
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/dMass.liquid' w l lc 2 notitle, \
+    'postProcessing_none.gas_none.liquid/plot/0/dEnergy.gas' w l axes x1y2 lc 3 t 'Gas Energy Change', \
+    for [g in gasSpecies] for [l in liquidSpecies] \
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/dEnergy.gas' w l axes x1y2 lc 3 notitle, \
+    'postProcessing_none.gas_none.liquid/plot/0/dEnergy.liquid' w l axes x1y2 lc 4 t 'Liquid Energy Change', \
+    for [g in gasSpecies] for [l in liquidSpecies] \
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/dEnergy.liquid' w l axes x1y2 lc 4 notitle
+
+set ytics nomirror
+set y2tics
+set ylabel "Mass (kg/m^3)"
+set y2label "Energy (J/m^3)"
+
+plot \
+    'postProcessing_none.gas_none.liquid/plot/0/dMass' w l lc 1 t 'Mass Error', \
+    for [g in gasSpecies] for [l in liquidSpecies] \
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/dMass' w l lc 1 notitle, \
+    'postProcessing_none.gas_none.liquid/plot/0/dEnergy' w l axes x1y2 lc 2 t 'Energy Error', \
+    for [g in gasSpecies] for [l in liquidSpecies] \
+    'postProcessing_'.g.'.gas_'.l.'.liquid/plot/0/dEnergy' w l axes x1y2 lc 2 notitle
+
 unset multiplot
 
 EOF

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterAndIsopropanolEvaporation/system/continuity

@@ -0,0 +1,133 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+
+libs            ("libutilityFunctionObjects.so");
+
+type            coded;
+
+writeControl    timeStep;
+writeInterval   1;
+
+name            continuity;
+
+codeOptions
+#{
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+#};
+
+codeInclude
+#{
+    #include "rhoThermo.H"
+#};
+
+codeRead
+#{
+    volScalarField* dMass =
+        new volScalarField
+        (
+            IOobject("dMass", mesh().time().timeName(), mesh()),
+            mesh(),
+            dimensionedScalar(dimMass/dimVolume, 0)
+        );
+    volScalarField* dMassGas = new volScalarField("dMass.gas", *dMass);
+    volScalarField* mass0Gas = new volScalarField("mass0.gas", *dMass);
+    volScalarField* dMassLiq = new volScalarField("dMass.liquid", *dMass);
+    volScalarField* mass0Liq = new volScalarField("mass0.liquid", *dMass);
+    volScalarField* dEnergy =
+        new volScalarField
+        (
+            IOobject("dEnergy", mesh().time().timeName(), mesh()),
+            mesh(),
+            dimensionedScalar(dimEnergy/dimVolume, 0)
+        );
+    volScalarField* dEnergyLiq = new volScalarField("dEnergy.liquid", *dEnergy);
+    volScalarField* energy0Liq = new volScalarField("energy0.liquid", *dEnergy);
+    volScalarField* dEnergyGas = new volScalarField("dEnergy.gas", *dEnergy);
+    volScalarField* energy0Gas = new volScalarField("energy0.gas", *dEnergy);
+
+    const volScalarField& alphaGas =
+        mesh().lookupObject<volScalarField>("alpha.gas");
+    const volScalarField& alphaLiq =
+        mesh().lookupObject<volScalarField>("alpha.liquid");
+    const rhoThermo& thermoGas =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.gas");
+    const rhoThermo& thermoLiq =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.liquid");
+
+    *mass0Gas = alphaGas*thermoGas.rho();
+    *mass0Liq = alphaLiq*thermoLiq.rho();
+
+    // !!! Note that this is only a valid measure of energy
+    // conservation because the case is of fixed volume, has no
+    // velocity, and is set up to solve in terms of sensible internal
+    // energy
+    *energy0Gas = alphaGas*thermoGas.rho()*(thermoGas.he() + thermoGas.hc());
+    *energy0Liq = alphaLiq*thermoLiq.rho()*(thermoLiq.he() + thermoLiq.hc());
+
+    dMass->store();
+    dMassGas->store();
+    mass0Gas->store();
+    dMassLiq->store();
+    mass0Liq->store();
+    dEnergy->store();
+    dEnergyGas->store();
+    energy0Gas->store();
+    dEnergyLiq->store();
+    energy0Liq->store();
+#};
+
+codeExecute
+#{
+    volScalarField& dMass =
+        mesh().lookupObjectRef<volScalarField>("dMass");
+    volScalarField& dMassGas =
+        mesh().lookupObjectRef<volScalarField>("dMass.gas");
+    const volScalarField& mass0Gas =
+        mesh().lookupObject<volScalarField>("mass0.gas");
+    volScalarField& dMassLiq =
+        mesh().lookupObjectRef<volScalarField>("dMass.liquid");
+    const volScalarField& mass0Liq =
+        mesh().lookupObject<volScalarField>("mass0.liquid");
+
+    volScalarField& dEnergy =
+        mesh().lookupObjectRef<volScalarField>("dEnergy");
+    volScalarField& dEnergyGas =
+        mesh().lookupObjectRef<volScalarField>("dEnergy.gas");
+    const volScalarField& energy0Gas =
+        mesh().lookupObject<volScalarField>("energy0.gas");
+    volScalarField& dEnergyLiq =
+        mesh().lookupObjectRef<volScalarField>("dEnergy.liquid");
+    const volScalarField& energy0Liq =
+        mesh().lookupObject<volScalarField>("energy0.liquid");
+
+    const volScalarField& alphaGas =
+        mesh().lookupObject<volScalarField>("alpha.gas");
+    const volScalarField& alphaLiq =
+        mesh().lookupObject<volScalarField>("alpha.liquid");
+    const rhoThermo& thermoGas =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.gas");
+    const rhoThermo& thermoLiq =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.liquid");
+
+    dMassGas = alphaGas*thermoGas.rho() - mass0Gas;
+    dMassLiq = alphaLiq*thermoLiq.rho() - mass0Liq;
+    dMass = dMassGas + dMassLiq;
+
+    // !!! Note that this is only a valid measure of energy
+    // conservation because the case is of fixed volume, has no
+    // velocity, and is set up to solve in terms of sensible internal
+    // energy
+    dEnergyGas = alphaGas*thermoGas.rho()*(thermoGas.he() + thermoGas.hc()) - energy0Gas;
+    dEnergyLiq = alphaLiq*thermoLiq.rho()*(thermoLiq.he() + thermoLiq.hc()) - energy0Liq;
+    dEnergy = dEnergyGas + dEnergyLiq;
+#};
+
+// ************************************************************************* //

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterAndIsopropanolEvaporation/system/controlDict

@@ -53,28 +53,8 @@ maxDeltaT       1;
 
 functions
 {
-    plot
+    #includeFunc continuity
-    {
+    #includeFunc plot
-        lib             ("libsampling.so");
-        type            probes;
-        writeControl    timeStep;
-        writeInterval   1;
-        fields
-        (
-            air.gas
-            alpha.gas
-            alpha.liquid
-            H2O.gas
-            H2O.liquid
-            C3H8O.gas
-            C3H8O.liquid
-            p
-            T.gas
-            T.liquid
-        );
-        probeLocations  ((0 0 0));
 }
-}
-
 
 // ************************************************************************* //

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterAndIsopropanolEvaporation/system/plot

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+
+libs            ("libsampling.so");
+
+type            probes;
+
+writeControl    timeStep;
+writeInterval   1;
+
+fields
+(
+    air.gas
+    alpha.gas
+    alpha.liquid
+    H2O.gas
+    H2O.liquid
+    C3H8O.gas
+    C3H8O.liquid
+    p
+    T.gas
+    T.liquid
+    Tf.gasAndLiquid
+    dMass.gas
+    dMass.liquid
+    dMass
+    dEnergy.gas
+    dEnergy.liquid
+    dEnergy
+);
+
+probeLocations  ((0 0 0));
+
+// ************************************************************************* //

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterEvaporation/Allrun

@@ -26,7 +26,7 @@ set rmargin at screen 0.84
 
 gasSpecies = '$gasSpecies'
 
-set multiplot layout 3,1
+set multiplot layout 5,1
 
 set xlabel "Time (s)"
 
@@ -54,7 +54,10 @@ plot \
     'postProcessing_'.g.'.gas/plot/0/T.gas' w l lc 1 notitle, \
     'postProcessing_none.gas/plot/0/T.liquid' w l lc 2 t 'Liquid', \
     for [g in gasSpecies] \
-    'postProcessing_'.g.'.gas/plot/0/T.liquid' w l lc 2 notitle
+    'postProcessing_'.g.'.gas/plot/0/T.liquid' w l lc 2 notitle, \
+    'postProcessing_none.gas/plot/0/Tf.gasAndLiquid' w l lc 3 t 'Interface', \
+    for [g in gasSpecies] \
+    'postProcessing_'.g.'.gas/plot/0/Tf.gasAndLiquid' w l lc 3 notitle
 
 set ytics nomirror
 set y2tics
@@ -69,6 +72,38 @@ plot \
     for [g in gasSpecies] \
     'postProcessing_'.g.'.gas/plot/0/air.gas' w l axes x1y2 lc 2 notitle
 
+set ytics nomirror
+set y2tics
+set ylabel "Mass (kg/m^3)"
+set y2label "Energy (J/m^3)"
+
+plot \
+    'postProcessing_none.gas/plot/0/dMass.gas' w l lc 1 t 'Gas Mass Change', \
+    for [g in gasSpecies] \
+    'postProcessing_'.g.'.gas/plot/0/dMass.gas' w l lc 1 notitle, \
+    'postProcessing_none.gas/plot/0/dMass.liquid' w l lc 2 t 'Liquid Mass Change', \
+    for [g in gasSpecies] \
+    'postProcessing_'.g.'.gas/plot/0/dMass.liquid' w l lc 2 notitle, \
+    'postProcessing_none.gas/plot/0/dEnergy.gas' w l axes x1y2 lc 3 t 'Gas Energy Change', \
+    for [g in gasSpecies] \
+    'postProcessing_'.g.'.gas/plot/0/dEnergy.gas' w l axes x1y2 lc 3 notitle, \
+    'postProcessing_none.gas/plot/0/dEnergy.liquid' w l axes x1y2 lc 4 t 'Liquid Energy Change', \
+    for [g in gasSpecies] \
+    'postProcessing_'.g.'.gas/plot/0/dEnergy.liquid' w l axes x1y2 lc 4 notitle
+
+set ytics nomirror
+set y2tics
+set ylabel "Mass (kg/m^3)"
+set y2label "Energy (J/m^3)"
+
+plot \
+    'postProcessing_none.gas/plot/0/dMass' w l lc 1 t 'Mass Error', \
+    for [g in gasSpecies] \
+    'postProcessing_'.g.'.gas/plot/0/dMass' w l lc 1 notitle, \
+    'postProcessing_none.gas/plot/0/dEnergy' w l axes x1y2 lc 2 t 'Energy Error', \
+    for [g in gasSpecies] \
+    'postProcessing_'.g.'.gas/plot/0/dEnergy' w l axes x1y2 lc 2 notitle
+
 unset multiplot
 
 EOF

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterEvaporation/constant/phaseProperties

@@ -122,11 +122,6 @@ massTransfer.gas
 
 massTransfer.liquid
 (
-    (liquid in gas)
-    {
-        type spherical;
-        Le 1.0;
-    }
 );
 
 phaseTransfer

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterEvaporation/system/continuity

@@ -0,0 +1,133 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+
+libs            ("libutilityFunctionObjects.so");
+
+type            coded;
+
+writeControl    timeStep;
+writeInterval   1;
+
+name            continuity;
+
+codeOptions
+#{
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+#};
+
+codeInclude
+#{
+    #include "rhoThermo.H"
+#};
+
+codeRead
+#{
+    volScalarField* dMass =
+        new volScalarField
+        (
+            IOobject("dMass", mesh().time().timeName(), mesh()),
+            mesh(),
+            dimensionedScalar(dimMass/dimVolume, 0)
+        );
+    volScalarField* dMassGas = new volScalarField("dMass.gas", *dMass);
+    volScalarField* mass0Gas = new volScalarField("mass0.gas", *dMass);
+    volScalarField* dMassLiq = new volScalarField("dMass.liquid", *dMass);
+    volScalarField* mass0Liq = new volScalarField("mass0.liquid", *dMass);
+    volScalarField* dEnergy =
+        new volScalarField
+        (
+            IOobject("dEnergy", mesh().time().timeName(), mesh()),
+            mesh(),
+            dimensionedScalar(dimEnergy/dimVolume, 0)
+        );
+    volScalarField* dEnergyLiq = new volScalarField("dEnergy.liquid", *dEnergy);
+    volScalarField* energy0Liq = new volScalarField("energy0.liquid", *dEnergy);
+    volScalarField* dEnergyGas = new volScalarField("dEnergy.gas", *dEnergy);
+    volScalarField* energy0Gas = new volScalarField("energy0.gas", *dEnergy);
+
+    const volScalarField& alphaGas =
+        mesh().lookupObject<volScalarField>("alpha.gas");
+    const volScalarField& alphaLiq =
+        mesh().lookupObject<volScalarField>("alpha.liquid");
+    const rhoThermo& thermoGas =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.gas");
+    const rhoThermo& thermoLiq =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.liquid");
+
+    *mass0Gas = alphaGas*thermoGas.rho();
+    *mass0Liq = alphaLiq*thermoLiq.rho();
+
+    // !!! Note that this is only a valid measure of energy
+    // conservation because the case is of fixed volume, has no
+    // velocity, and is set up to solve in terms of sensible internal
+    // energy
+    *energy0Gas = alphaGas*thermoGas.rho()*(thermoGas.he() + thermoGas.hc());
+    *energy0Liq = alphaLiq*thermoLiq.rho()*(thermoLiq.he() + thermoLiq.hc());
+
+    dMass->store();
+    dMassGas->store();
+    mass0Gas->store();
+    dMassLiq->store();
+    mass0Liq->store();
+    dEnergy->store();
+    dEnergyGas->store();
+    energy0Gas->store();
+    dEnergyLiq->store();
+    energy0Liq->store();
+#};
+
+codeExecute
+#{
+    volScalarField& dMass =
+        mesh().lookupObjectRef<volScalarField>("dMass");
+    volScalarField& dMassGas =
+        mesh().lookupObjectRef<volScalarField>("dMass.gas");
+    const volScalarField& mass0Gas =
+        mesh().lookupObject<volScalarField>("mass0.gas");
+    volScalarField& dMassLiq =
+        mesh().lookupObjectRef<volScalarField>("dMass.liquid");
+    const volScalarField& mass0Liq =
+        mesh().lookupObject<volScalarField>("mass0.liquid");
+
+    volScalarField& dEnergy =
+        mesh().lookupObjectRef<volScalarField>("dEnergy");
+    volScalarField& dEnergyGas =
+        mesh().lookupObjectRef<volScalarField>("dEnergy.gas");
+    const volScalarField& energy0Gas =
+        mesh().lookupObject<volScalarField>("energy0.gas");
+    volScalarField& dEnergyLiq =
+        mesh().lookupObjectRef<volScalarField>("dEnergy.liquid");
+    const volScalarField& energy0Liq =
+        mesh().lookupObject<volScalarField>("energy0.liquid");
+
+    const volScalarField& alphaGas =
+        mesh().lookupObject<volScalarField>("alpha.gas");
+    const volScalarField& alphaLiq =
+        mesh().lookupObject<volScalarField>("alpha.liquid");
+    const rhoThermo& thermoGas =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.gas");
+    const rhoThermo& thermoLiq =
+        mesh().lookupObject<rhoThermo>("thermophysicalProperties.liquid");
+
+    dMassGas = alphaGas*thermoGas.rho() - mass0Gas;
+    dMassLiq = alphaLiq*thermoLiq.rho() - mass0Liq;
+    dMass = dMassGas + dMassLiq;
+
+    // !!! Note that this is only a valid measure of energy
+    // conservation because the case is of fixed volume, has no
+    // velocity, and is set up to solve in terms of sensible internal
+    // energy
+    dEnergyGas = alphaGas*thermoGas.rho()*(thermoGas.he() + thermoGas.hc()) - energy0Gas;
+    dEnergyLiq = alphaLiq*thermoLiq.rho()*(thermoLiq.he() + thermoLiq.hc()) - energy0Liq;
+    dEnergy = dEnergyGas + dEnergyLiq;
+#};
+
+// ************************************************************************* //

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterEvaporation/system/controlDict

@@ -53,25 +53,8 @@ maxDeltaT       1;
 
 functions
 {
-    plot
+    #includeFunc continuity
-    {
+    #includeFunc plot
-        lib             ("libsampling.so");
-        type            probes;
-        writeControl    timeStep;
-        writeInterval   1;
-        fields
-        (
-            air.gas
-            alpha.gas
-            alpha.liquid
-            H2O.gas
-            p
-            T.gas
-            T.liquid
-        );
-        probeLocations  ((0 0 0));
 }
-}
-
 
 // ************************************************************************* //

test/multiphase/reactingTwoPhaseEulerFoam/interfaceComposition/waterEvaporation/system/plot

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+
+libs            ("libsampling.so");
+
+type            probes;
+
+writeControl    timeStep;
+writeInterval   1;
+
+fields
+(
+    air.gas
+    alpha.gas
+    alpha.liquid
+    H2O.gas
+    p
+    T.gas
+    T.liquid
+    Tf.gasAndLiquid
+    dMass.gas
+    dMass.liquid
+    dMass
+    dEnergy.gas
+    dEnergy.liquid
+    dEnergy
+);
+
+probeLocations  ((0 0 0));
+
+// ************************************************************************* //