diff --git a/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H b/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H
index 02375780ac..076ad5ea03 100644
--- a/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H
+++ b/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H
@@ -32,6 +32,9 @@ Description
         rho = rho0*(1 - beta*(T - T0))
     \endverbatim
 
+    Coefficient mixing is very inaccurate and not supported,
+    so this equation of state is not applicable to mixtures.
+
 SourceFiles
     BoussinesqI.H
     Boussinesq.C
@@ -172,7 +175,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H b/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H
index 74e3a581b3..eeefe7fbb9 100644
--- a/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H
+++ b/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -187,18 +187,7 @@ inline void Foam::Boussinesq<Specie>::operator+=
     const Boussinesq<Specie>& b
 )
 {
-    scalar Y1 = this->Y();
-    Specie::operator+=(b);
-
-    if (mag(this->Y()) > small)
-    {
-        Y1 /= this->Y();
-        const scalar Y2 = b.Y()/this->Y();
-
-        rho0_ = Y1*rho0_ + Y2*b.rho0_;
-        T0_ = Y1*T0_ + Y2*b.T0_;
-        beta_ = Y1*beta_ + Y2*b.beta_;
-    }
+    NotImplemented;
 }
 
 
@@ -218,31 +207,8 @@ inline Foam::Boussinesq<Specie> Foam::operator+
     const Boussinesq<Specie>& b2
 )
 {
-    Specie sp(static_cast<const Specie&>(b1) + static_cast<const Specie&>(b2));
-
-    if (mag(sp.Y()) < small)
-    {
-        return Boussinesq<Specie>
-        (
-            sp,
-            b1.rho0_,
-            b1.T0_,
-            b1.beta_
-        );
-    }
-    else
-    {
-        const scalar Y1 = b1.Y()/sp.Y();
-        const scalar Y2 = b2.Y()/sp.Y();
-
-        return Boussinesq<Specie>
-        (
-            sp,
-            Y1*b1.rho0_ + Y2*b2.rho0_,
-            Y1*b1.T0_ + Y2*b2.T0_,
-            Y1*b1.beta_ + Y2*b2.beta_
-        );
-    }
+    NotImplemented;
+    return b1;
 }
 
 
@@ -270,18 +236,8 @@ inline Foam::Boussinesq<Specie> Foam::operator==
     const Boussinesq<Specie>& b2
 )
 {
-    Specie sp(static_cast<const Specie&>(b1) == static_cast<const Specie&>(b2));
-
-    const scalar Y1 = b1.Y()/sp.Y();
-    const scalar Y2 = b2.Y()/sp.Y();
-
-    return Boussinesq<Specie>
-    (
-        sp,
-        Y2*b2.rho0_ - Y1*b1.rho0_,
-        Y2*b2.T0_   - Y1*b1.T0_,
-        Y2*b2.beta_ - Y1*b1.beta_
-    );
+    NotImplemented;
+    return b1;
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
index 29805f93fd..88a4ee1fd8 100644
--- a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
@@ -170,7 +170,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
index dc16404ba2..9576e11500 100644
--- a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -338,19 +338,19 @@ inline void Foam::PengRobinsonGas<Specie>::operator+=
     const PengRobinsonGas<Specie>& pg
 )
 {
-    scalar Y1 = this->Y();
+    scalar X1 = this->Y()/this->W();
     Specie::operator+=(pg);
 
     if (mag(this->Y()) > small)
     {
-        Y1 /= this->Y();
-        const scalar Y2 = pg.Y()/this->Y();
+        X1 *= this->W()/this->Y();
+        const scalar X2 = this->W()*pg.Y()/(pg.W()*this->Y());
 
-        Tc_ = Y1*Tc_ + Y2*pg.Tc_;
-        Vc_ = Y1*Vc_ + Y2*pg.Vc_;
-        Zc_ = Y1*Zc_ + Y2*pg.Zc_;
+        Tc_ = X1*Tc_ + X2*pg.Tc_;
+        Vc_ = X1*Vc_ + X2*pg.Vc_;
+        Zc_ = X1*Zc_ + X2*pg.Zc_;
         Pc_ = RR*Zc_*Tc_/Vc_;
-        omega_ = Y1*omega_ + Y2*pg.omega_;
+        omega_ = X1*omega_ + X2*pg.omega_;
     }
 }
 
@@ -392,12 +392,12 @@ Foam::PengRobinsonGas<Specie> Foam::operator+
     }
     else
     {
-        const scalar Y1 = pg1.Y()/sp.Y();
-        const scalar Y2 = pg2.Y()/sp.Y();
+        const scalar X1 = sp.W()*pg1.Y()/(pg1.W()*sp.Y());
+        const scalar X2 = sp.W()*pg2.Y()/(pg2.W()*sp.Y());
 
-        const scalar Tc = Y1*pg1.Tc_ + Y2*pg2.Tc_;
-        const scalar Vc = Y1*pg1.Vc_ + Y2*pg2.Vc_;
-        const scalar Zc = Y1*pg1.Zc_ + Y2*pg2.Zc_;
+        const scalar Tc = X1*pg1.Tc_ + X2*pg2.Tc_;
+        const scalar Vc = X1*pg1.Vc_ + X2*pg2.Vc_;
+        const scalar Zc = X1*pg1.Zc_ + X2*pg2.Zc_;
 
         return PengRobinsonGas<Specie>
         (
@@ -406,7 +406,7 @@ Foam::PengRobinsonGas<Specie> Foam::operator+
             Vc,
             Zc,
             RR*Zc*Tc/Vc,
-            Y1*pg1.omega_ + Y2*pg2.omega_
+            X1*pg1.omega_ + X2*pg2.omega_
         );
     }
 }
@@ -438,27 +438,14 @@ Foam::PengRobinsonGas<Specie> Foam::operator==
     const PengRobinsonGas<Specie>& pg2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(pg1)
-     == static_cast<const Specie&>(pg2)
-    );
-
-    const scalar Y1 = pg1.Y()/sp.Y();
-    const scalar Y2 = pg2.Y()/sp.Y();
-
-    const scalar Tc = Y2*pg2.Tc_ - Y1*pg1.Tc_;
-    const scalar Vc = Y2*pg2.Vc_ - Y1*pg1.Vc_;
-    const scalar Zc = Y2*pg2.Zc_ - Y1*pg1.Zc_;
-
     return PengRobinsonGas<Specie>
     (
-        sp,
-        Tc,
-        Vc,
-        Zc,
-        RR*Zc*Tc/Vc,
-        Y2*pg2.omega_ - Y1*pg1.omega_
+        static_cast<const Specie&>(pg1) == static_cast<const Specie&>(pg2),
+        NaN,
+        NaN,
+        NaN,
+        NaN,
+        NaN
     );
 }
 
diff --git a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H
index c6ba0ea890..59ce65bf6a 100644
--- a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H
+++ b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H
@@ -27,6 +27,9 @@ Class
 Description
     Adiabatic perfect fluid equation of state.
 
+    Coefficient mixing is very inaccurate and not supported,
+    so this equation of state is not applicable to mixtures.
+
 SourceFiles
     adiabaticPerfectFluidI.H
     adiabaticPerfectFluid.C
@@ -169,7 +172,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H
index 8c838244fa..593649a43d 100644
--- a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H
+++ b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -199,19 +199,7 @@ inline void Foam::adiabaticPerfectFluid<Specie>::operator+=
     const adiabaticPerfectFluid<Specie>& pf
 )
 {
-    scalar Y1 = this->Y();
-    Specie::operator+=(pf);
-
-    if (mag(this->Y()) > small)
-    {
-        Y1 /= this->Y();
-        const scalar Y2 = pf.Y()/this->Y();
-
-        p0_ = Y1*p0_ + Y2*pf.p0_;
-        rho0_ = Y1*rho0_ + Y2*pf.rho0_;
-        gamma_ = Y1*gamma_ + Y2*pf.gamma_;
-        B_ = Y1*B_ + Y2*pf.B_;
-    }
+    NotImplemented;
 }
 
 
@@ -231,37 +219,8 @@ inline Foam::adiabaticPerfectFluid<Specie> Foam::operator+
     const adiabaticPerfectFluid<Specie>& pf2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(pf1)
-      + static_cast<const Specie&>(pf2)
-    );
-
-    if (mag(sp.Y()) < small)
-    {
-        return adiabaticPerfectFluid<Specie>
-        (
-            sp,
-            pf1.p0_,
-            pf1.rho0_,
-            pf1.gamma_,
-            pf1.B_
-        );
-    }
-    else
-    {
-        const scalar Y1 = pf1.Y()/sp.Y();
-        const scalar Y2 = pf2.Y()/sp.Y();
-
-        return adiabaticPerfectFluid<Specie>
-        (
-            sp,
-            Y1*pf1.p0_ + Y2*pf2.p0_,
-            Y1*pf1.rho0_ + Y2*pf2.rho0_,
-            Y1*pf1.gamma_ + Y2*pf2.gamma_,
-            Y1*pf1.B_ + Y2*pf2.B_
-        );
-    }
+    NotImplemented;
+    return pf1;
 }
 
 
@@ -290,23 +249,8 @@ inline Foam::adiabaticPerfectFluid<Specie> Foam::operator==
     const adiabaticPerfectFluid<Specie>& pf2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(pf1)
-     == static_cast<const Specie&>(pf2)
-    );
-
-    const scalar Y1 = pf1.Y()/sp.Y();
-    const scalar Y2 = pf2.Y()/sp.Y();
-
-    return adiabaticPerfectFluid<Specie>
-    (
-        sp,
-        Y2*pf2.p0_    - Y1*pf1.p0_,
-        Y2*pf2.rho0_  - Y1*pf1.rho0_,
-        Y2*pf2.gamma_ - Y1*pf1.gamma_,
-        Y2*pf2.B_     - Y1*pf1.B_
-    );
+    NotImplemented;
+    return pf1;
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
index 41a5e887df..e910e1ea84 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
@@ -28,6 +28,9 @@ Description
     Incompressible, polynomial form of equation of state, using a polynomial
     function for density.
 
+    Coefficient mixing is very inaccurate and not supported,
+    so this equation of state is not applicable to mixtures.
+
 Usage
     \table
         Property     | Description
@@ -176,7 +179,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
index 6a4304a91f..6c3c307a86 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
@@ -184,16 +184,7 @@ inline void Foam::icoPolynomial<Specie, PolySize>::operator+=
     const icoPolynomial<Specie, PolySize>& ip
 )
 {
-    scalar Y1 = this->Y();
-    Specie::operator+=(ip);
-
-    if (mag(this->Y()) > small)
-    {
-        Y1 /= this->Y();
-        const scalar Y2 = ip.Y()/this->Y();
-
-        rhoCoeffs_ = Y1*rhoCoeffs_ + Y2*ip.rhoCoeffs_;
-    }
+    NotImplemented;
 }
 
 
@@ -213,31 +204,8 @@ Foam::icoPolynomial<Specie, PolySize> Foam::operator+
     const icoPolynomial<Specie, PolySize>& ip2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(ip1)
-      + static_cast<const Specie&>(ip2)
-    );
-
-    if (mag(sp.Y()) < small)
-    {
-        return icoPolynomial<Specie, PolySize>
-        (
-            sp,
-            ip1.rhoCoeffs_
-        );
-    }
-    else
-    {
-        const scalar Y1 = ip1.Y()/sp.Y();
-        const scalar Y2 = ip2.Y()/sp.Y();
-
-        return icoPolynomial<Specie, PolySize>
-        (
-            sp,
-            Y1*ip1.rhoCoeffs_ + Y2*ip2.rhoCoeffs_
-        );
-    }
+    NotImplemented;
+    return ip1;
 }
 
 
@@ -263,20 +231,8 @@ Foam::icoPolynomial<Specie, PolySize> Foam::operator==
     const icoPolynomial<Specie, PolySize>& ip2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(ip1)
-     == static_cast<const Specie&>(ip2)
-    );
-
-    const scalar Y1 = ip1.Y()/sp.Y();
-    const scalar Y2 = ip2.Y()/sp.Y();
-
-    return icoPolynomial<Specie, PolySize>
-    (
-        sp,
-        Y2*ip2.rhoCoeffs_ - Y1*ip1.rhoCoeffs_
-    );
+    NotImplemented;
+    return ip1;
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H
index 6b30f00961..1321317b3c 100644
--- a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H
+++ b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H
@@ -157,7 +157,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H
index 2c48cd814f..4486de6272 100644
--- a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H
+++ b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H
@@ -185,16 +185,17 @@ inline void Foam::incompressiblePerfectGas<Specie>::operator+=
     const incompressiblePerfectGas<Specie>& ipg
 )
 {
-    scalar Y1 = this->Y();
-    Specie::operator+=(ipg);
-
-    if (mag(this->Y()) > small)
+    if (notEqual(pRef_, ipg.pRef_))
     {
-        Y1 /= this->Y();
-        const scalar Y2 = ipg.Y()/this->Y();
-
-        pRef_ = Y1*pRef_ + Y2*ipg.pRef_;
+        FatalErrorInFunction
+            << "pRef " << pRef_ << " for "
+            << (this->name().size() ? this->name() : "others")
+            << " != " << ipg.pRef_ << " for "
+            << (ipg.name().size() ? ipg.name() : "others")
+            << exit(FatalError);
     }
+
+    Specie::operator+=(ipg);
 }
 
 
@@ -214,31 +215,21 @@ inline Foam::incompressiblePerfectGas<Specie> Foam::operator+
     const incompressiblePerfectGas<Specie>& ipg2
 )
 {
-    Specie sp
+    if (notEqual(ipg1.pRef_, ipg2.pRef_))
+    {
+        FatalErrorInFunction
+            << "pRef " << ipg1.pRef_ << " for "
+            << (ipg1.name().size() ? ipg1.name() : "others")
+            << " != " << ipg2.pRef_ << " for "
+            << (ipg2.name().size() ? ipg2.name() : "others")
+            << exit(FatalError);
+    }
+
+    return incompressiblePerfectGas<Specie>
     (
-        static_cast<const Specie&>(ipg1)
-      + static_cast<const Specie&>(ipg2)
+        static_cast<const Specie&>(ipg1) + static_cast<const Specie&>(ipg2),
+        ipg1.pRef_
     );
-
-    if (mag(sp.Y()) < small)
-    {
-        return incompressiblePerfectGas<Specie>
-        (
-            sp,
-            ipg1.pRef_
-        );
-    }
-    else
-    {
-        const scalar Y1 = ipg1.Y()/sp.Y();
-        const scalar Y2 = ipg2.Y()/sp.Y();
-
-        return incompressiblePerfectGas<Specie>
-        (
-            sp,
-            Y1*ipg1.pRef_ + Y2*ipg2.pRef_
-        );
-    }
 }
 
 
@@ -264,19 +255,10 @@ inline Foam::incompressiblePerfectGas<Specie> Foam::operator==
     const incompressiblePerfectGas<Specie>& ipg2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(ipg1)
-     == static_cast<const Specie&>(ipg2)
-    );
-
-    const scalar Y1 = ipg1.Y()/sp.Y();
-    const scalar Y2 = ipg2.Y()/sp.Y();
-
     return incompressiblePerfectGas<Specie>
     (
-        sp,
-        Y2*ipg2.pRef_ - Y1*ipg1.pRef_
+        static_cast<const Specie&>(ipg1) == static_cast<const Specie&>(ipg2),
+        NaN
     );
 }
 
diff --git a/src/thermophysicalModels/specie/equationOfState/linear/linear.H b/src/thermophysicalModels/specie/equationOfState/linear/linear.H
index d49103c5b8..c8aac35688 100644
--- a/src/thermophysicalModels/specie/equationOfState/linear/linear.H
+++ b/src/thermophysicalModels/specie/equationOfState/linear/linear.H
@@ -31,6 +31,9 @@ Description
         rho = rho0 + psi*p
     \endverbatim
 
+    Coefficient mixing is very inaccurate and not supported,
+    so this equation of state is not applicable to mixtures.
+
 SourceFiles
     linearI.H
     linear.C
@@ -158,7 +161,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/linear/linearI.H b/src/thermophysicalModels/specie/equationOfState/linear/linearI.H
index 6ec30b4793..ced54ac47f 100644
--- a/src/thermophysicalModels/specie/equationOfState/linear/linearI.H
+++ b/src/thermophysicalModels/specie/equationOfState/linear/linearI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -150,17 +150,7 @@ inline void Foam::linear<Specie>::operator+=
     const linear<Specie>& pf
 )
 {
-    scalar Y1 = this->Y();
-    Specie::operator+=(pf);
-
-    if (mag(this->Y()) > small)
-    {
-        Y1 /= this->Y();
-        const scalar Y2 = pf.Y()/this->Y();
-
-        psi_ = Y1*psi_ + Y2*pf.psi_;
-        rho0_ = Y1*rho0_ + Y2*pf.rho0_;
-    }
+    NotImplemented;
 }
 
 
@@ -180,33 +170,8 @@ inline Foam::linear<Specie> Foam::operator+
     const linear<Specie>& pf2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(pf1)
-      + static_cast<const Specie&>(pf2)
-    );
-
-    if (mag(sp.Y()) < small)
-    {
-        return linear<Specie>
-        (
-            sp,
-            pf1.psi_,
-            pf1.rho0_
-        );
-    }
-    else
-    {
-        const scalar Y1 = pf1.Y()/sp.Y();
-        const scalar Y2 = pf2.Y()/sp.Y();
-
-        return linear<Specie>
-        (
-            sp,
-            Y1*pf1.psi_ + Y2*pf2.psi_,
-            Y1*pf1.rho0_ + Y2*pf2.rho0_
-        );
-    }
+    NotImplemented;
+    return pf1;
 }
 
 
@@ -233,21 +198,8 @@ inline Foam::linear<Specie> Foam::operator==
     const linear<Specie>& pf2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(pf1)
-     == static_cast<const Specie&>(pf2)
-    );
-
-    const scalar Y1 = pf1.Y()/sp.Y();
-    const scalar Y2 = pf2.Y()/sp.Y();
-
-    return linear<Specie>
-    (
-        sp,
-        Y2*pf2.psi_  - Y1*pf1.psi_,
-        Y2*pf2.rho0_ - Y1*pf1.rho0_
-    );
+    NotImplemented;
+    return pf1;
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H
index 30eaf30ca4..495fd5f79c 100644
--- a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H
+++ b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H
@@ -25,7 +25,14 @@ Class
     Foam::perfectFluid
 
 Description
-    Perfect gas equation of state.
+    Simple extension of the perfect gas equation of state to liquids by
+    the addition of a constant density off-set.
+
+    Coefficient mixing is very inaccurate and not supported,
+    so this equation of state is not applicable to mixtures.
+
+    This equation of state is rather inaccurate and has been superseded by
+    rPolynomial which is much more accurate and supports mixtures.
 
 SourceFiles
     perfectFluidI.H
@@ -157,7 +164,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H
index 9204872f40..933908c86b 100644
--- a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H
+++ b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H
@@ -160,17 +160,7 @@ inline void Foam::perfectFluid<Specie>::operator+=
     const perfectFluid<Specie>& pf
 )
 {
-    scalar Y1 = this->Y();
-    Specie::operator+=(pf);
-
-    if (mag(this->Y()) > small)
-    {
-        Y1 /= this->Y();
-        const scalar Y2 = pf.Y()/this->Y();
-
-        R_ = 1.0/(Y1/R_ + Y2/pf.R_);
-        rho0_ = Y1*rho0_ + Y2*pf.rho0_;
-    }
+    NotImplemented;
 }
 
 
@@ -190,33 +180,8 @@ inline Foam::perfectFluid<Specie> Foam::operator+
     const perfectFluid<Specie>& pf2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(pf1)
-      + static_cast<const Specie&>(pf2)
-    );
-
-    if (mag(sp.Y()) < small)
-    {
-        return perfectFluid<Specie>
-        (
-            sp,
-            pf1.R_,
-            pf1.rho0_
-        );
-    }
-    else
-    {
-        const scalar Y1 = pf1.Y()/sp.Y();
-        const scalar Y2 = pf2.Y()/sp.Y();
-
-        return perfectFluid<Specie>
-        (
-            sp,
-            1.0/(Y1/pf1.R_ + Y2/pf2.R_),
-            Y1*pf1.rho0_ + Y2*pf2.rho0_
-        );
-    }
+    NotImplemented;
+    return pf1;
 }
 
 
@@ -243,34 +208,8 @@ inline Foam::perfectFluid<Specie> Foam::operator==
     const perfectFluid<Specie>& pf2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(pf1)
-     == static_cast<const Specie&>(pf2)
-    );
-
-    if (mag(sp.Y()) < small)
-    {
-        return perfectFluid<Specie>
-        (
-            sp,
-            pf1.R_,
-            pf1.rho0_
-        );
-    }
-    else
-    {
-        const scalar Y1 = pf1.Y()/sp.Y();
-        const scalar Y2 = pf2.Y()/sp.Y();
-        const scalar oneByR = Y2/pf2.R_ - Y1/pf1.R_;
-
-        return perfectFluid<Specie>
-        (
-            sp,
-            mag(oneByR) < small ? great : 1/oneByR,
-            Y2*pf2.rho0_ - Y1*pf1.rho0_
-        );
-    }
+    NotImplemented;
+    return pf1;
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H
index dbbc8103bf..5cdd38f380 100644
--- a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H
+++ b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H
@@ -141,7 +141,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C
index 85b8937640..62d3f71691 100644
--- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C
+++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,9 +53,9 @@ void Foam::rhoConst<Specie>::write(Ostream& os) const
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Specie>
-Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst<Specie>& ico)
+Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst<Specie>& rc)
 {
-    ico.write(os);
+    rc.write(os);
     return os;
 }
 
diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H
index 0c8af21904..3725802afd 100644
--- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H
+++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H
@@ -143,7 +143,7 @@ public:
             //- Return entropy [J/kg/K]
             inline scalar S(const scalar p, const scalar T) const;
 
-            //- Return compressibility rho/p [s^2/m^2]
+            //- Return compressibility [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
             //- Return compression factor []
diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H
index 3373614b2f..d32d08cdb4 100644
--- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H
+++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,11 +45,11 @@ template<class Specie>
 inline Foam::rhoConst<Specie>::rhoConst
 (
     const word& name,
-    const rhoConst<Specie>& ico
+    const rhoConst<Specie>& rc
 )
 :
-    Specie(name, ico),
-    rho_(ico.rho_)
+    Specie(name, rc),
+    rho_(rc.rho_)
 {}
 
 
@@ -129,17 +129,14 @@ inline Foam::scalar Foam::rhoConst<Specie>::CpMCv(scalar p, scalar T) const
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Specie>
-inline void Foam::rhoConst<Specie>::operator+=(const rhoConst<Specie>& ico)
+inline void Foam::rhoConst<Specie>::operator+=(const rhoConst<Specie>& rc)
 {
-    scalar Y1 = this->Y();
-    Specie::operator+=(ico);
+    const scalar Y1 = this->Y();
+    Specie::operator+=(rc);
 
     if (mag(this->Y()) > small)
     {
-        Y1 /= this->Y();
-        const scalar Y2 = ico.Y()/this->Y();
-
-        rho_ = Y1*rho_ + Y2*ico.rho_;
+        rho_ = this->Y()/(Y1/rho_ + rc.Y()/rc.rho_);
     }
 }
 
@@ -156,14 +153,14 @@ inline void Foam::rhoConst<Specie>::operator*=(const scalar s)
 template<class Specie>
 inline Foam::rhoConst<Specie> Foam::operator+
 (
-    const rhoConst<Specie>& ico1,
-    const rhoConst<Specie>& ico2
+    const rhoConst<Specie>& rc1,
+    const rhoConst<Specie>& rc2
 )
 {
     Specie sp
     (
-        static_cast<const Specie&>(ico1)
-      + static_cast<const Specie&>(ico2)
+        static_cast<const Specie&>(rc1)
+      + static_cast<const Specie&>(rc2)
     );
 
     if (mag(sp.Y()) < small)
@@ -171,18 +168,15 @@ inline Foam::rhoConst<Specie> Foam::operator+
         return rhoConst<Specie>
         (
             sp,
-            ico1.rho_
+            rc1.rho_
         );
     }
     else
     {
-        const scalar Y1 = ico1.Y()/sp.Y();
-        const scalar Y2 = ico2.Y()/sp.Y();
-
         return rhoConst<Specie>
         (
             sp,
-            Y1*ico1.rho_ + Y2*ico2.rho_
+            sp.Y()/(rc1.Y()/rc1.rho_ + rc2.Y()/rc2.rho_)
         );
     }
 }
@@ -192,33 +186,24 @@ template<class Specie>
 inline Foam::rhoConst<Specie> Foam::operator*
 (
     const scalar s,
-    const rhoConst<Specie>& ico
+    const rhoConst<Specie>& rc
 )
 {
-    return rhoConst<Specie>(s*static_cast<const Specie&>(ico), ico.rho_);
+    return rhoConst<Specie>(s*static_cast<const Specie&>(rc), rc.rho_);
 }
 
 
 template<class Specie>
 inline Foam::rhoConst<Specie> Foam::operator==
 (
-    const rhoConst<Specie>& ico1,
-    const rhoConst<Specie>& ico2
+    const rhoConst<Specie>& rc1,
+    const rhoConst<Specie>& rc2
 )
 {
-    Specie sp
-    (
-        static_cast<const Specie&>(ico1)
-     == static_cast<const Specie&>(ico2)
-    );
-
-    const scalar Y1 = ico1.Y()/sp.Y();
-    const scalar Y2 = ico2.Y()/sp.Y();
-
     return rhoConst<Specie>
     (
-        sp,
-        Y2*ico2.rho_ - Y1*ico1.rho_
+        static_cast<const Specie&>(rc1) == static_cast<const Specie&>(rc2),
+        NaN
     );
 }
 
diff --git a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H
index f26e226883..016054d3ae 100644
--- a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H
+++ b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H
@@ -202,7 +202,7 @@ public:
 
         // Fundamental equation of state properties
 
-            //- Liquid compressibility rho/p [s^2/m^2]
+            //- Liquid compressibility [s^2/m^2]
             //  Note: currently it is assumed the liquid is incompressible
             inline scalar psi(scalar p, scalar T) const;
 
diff --git a/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H b/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H
index 57912a6bb6..741cef0819 100644
--- a/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H
+++ b/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H
@@ -124,7 +124,7 @@ public:
             //- Liquid density [kg/m^3]
             virtual scalar rho(scalar p, scalar T) const = 0;
 
-            //- Liquid compressibility rho/p [s^2/m^2]
+            //- Liquid compressibility [s^2/m^2]
             //  Note: currently it is assumed the liquid is incompressible
             virtual scalar psi(scalar p, scalar T) const = 0;
 
diff --git a/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H b/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H
index 1842ca0e39..c606b6e121 100644
--- a/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H
+++ b/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H
@@ -107,7 +107,7 @@ public:
             //- Liquid density [kg/m^3]
             inline scalar rho(scalar p, scalar T) const;
 
-            //- Liquid compressibility rho/p [s^2/m^2]
+            //- Liquid compressibility [s^2/m^2]
             //  Note: currently it is assumed the liquid is incompressible
             inline scalar psi(scalar p, scalar T) const;