multicomponentFluid: Correct boundary conditions of non-solved species

Whilst the cell values of non-solved species do not change, the boundary
values might, and correcting them is necessary for certain
post-processing operations to produce sensible results.

applications/modules/multicomponentFluid/thermophysicalPredictor.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,10 +45,10 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
 
     forAll(Y, i)
     {
+        volScalarField& Yi = Y_[i];
+
         if (thermo_.solveSpecie(i))
         {
-            volScalarField& Yi = Y_[i];
-
             fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
@@ -67,6 +67,10 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
 
             fvConstraints().constrain(Yi);
         }
+        else
+        {
+            Yi.correctBoundaryConditions();
+        }
     }
 
     thermo_.normaliseY();

applications/modules/multiphaseEuler/phaseSystem/phaseModel/phaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -296,11 +296,9 @@ public:
             //- Access the species mass fractions
             virtual PtrList<volScalarField>& YRef() = 0;
 
-            //- Return the active species mass fractions
+            //- Should the given specie be solved for? I.e., is it active and
-            virtual const UPtrList<volScalarField>& YActive() const = 0;
+            //  not the default specie?
-
+            virtual bool solveSpecie(const label speciei) const = 0;
-            //- Access the active species mass fractions
-            virtual UPtrList<volScalarField>& YActiveRef() = 0;
 
             //- Return the fuel consumption rate matrix
             virtual tmp<fvScalarMatrix> R(volScalarField& Yi) const = 0;

applications/modules/multiphaseEuler/phaseSystem/phaseModels/MulticomponentPhaseModel/MulticomponentPhaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,19 +47,7 @@ Foam::MulticomponentPhaseModel<BasePhaseModel>::MulticomponentPhaseModel
 )
 :
     BasePhaseModel(fluid, phaseName, referencePhase, index)
-{
+{}
-    PtrList<volScalarField>& Y = this->thermo_->Y();
-
-    forAll(Y, i)
-    {
-        if (this->thermo_->solveSpecie(i))
-        {
-            const label j = YActive_.size();
-            YActive_.resize(j + 1);
-            YActive_.set(j, &Y[i]);
-        }
-    }
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -142,18 +130,12 @@ Foam::MulticomponentPhaseModel<BasePhaseModel>::YRef()
 
 
 template<class BasePhaseModel>
-const Foam::UPtrList<Foam::volScalarField>&
+bool Foam::MulticomponentPhaseModel<BasePhaseModel>::solveSpecie
-Foam::MulticomponentPhaseModel<BasePhaseModel>::YActive() const
+(
+    const label speciei
+) const
 {
-    return YActive_;
+    return this->thermo_->solveSpecie(speciei);
-}
-
-
-template<class BasePhaseModel>
-Foam::UPtrList<Foam::volScalarField>&
-Foam::MulticomponentPhaseModel<BasePhaseModel>::YActiveRef()
-{
-    return YActive_;
 }
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModels/MulticomponentPhaseModel/MulticomponentPhaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,14 +52,6 @@ class MulticomponentPhaseModel
 :
     public BasePhaseModel
 {
-protected:
-
-    // Protected data
-
-        //- Pointer list to active species
-        UPtrList<volScalarField> YActive_;
-
-
 public:
 
     // Constructors
@@ -82,6 +74,7 @@ public:
         //- Correct the species fractions
         virtual void correctSpecies();
 
+
         // Species
 
             //- Return whether the phase is pure (i.e., not multi-component)
@@ -99,11 +92,9 @@ public:
             //- Access the species mass fractions
             virtual PtrList<volScalarField>& YRef();
 
-            //- Return the active species mass fractions
+            //- Should the given specie be solved for? I.e., is it active and
-            virtual const UPtrList<volScalarField>& YActive() const;
+            //  not the default specie?
-
+            virtual bool solveSpecie(const label speciei) const;
-            //- Access the active species mass fractions
-            virtual UPtrList<volScalarField>& YActiveRef();
 };
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModels/PurePhaseModel/PurePhaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,24 +104,12 @@ Foam::PurePhaseModel<BasePhaseModel>::YRef()
 
 
 template<class BasePhaseModel>
-const Foam::UPtrList<Foam::volScalarField>&
+bool Foam::PurePhaseModel<BasePhaseModel>::solveSpecie
-Foam::PurePhaseModel<BasePhaseModel>::YActive() const
+(
+    const label speciei
+) const
 {
-    // Y_ has never been set, so we are returning an empty list
+    return false;
-
-    return Y_;
-}
-
-
-template<class BasePhaseModel>
-Foam::UPtrList<Foam::volScalarField>&
-Foam::PurePhaseModel<BasePhaseModel>::YActiveRef()
-{
-    FatalErrorInFunction
-        << "Cannot access the species fractions of for a pure phase"
-        << exit(FatalError);
-
-    return Y_;
 }
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModels/PurePhaseModel/PurePhaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,11 +96,9 @@ public:
             //- Access the species mass fractions
             virtual PtrList<volScalarField>& YRef();
 
-            //- Return the active species mass fractions
+            //- Should the given specie be solved for? I.e., is it active and
-            virtual const UPtrList<volScalarField>& YActive() const;
+            //  not the default specie?
-
+            virtual bool solveSpecie(const label speciei) const;
-            //- Access the active species mass fractions
-            virtual UPtrList<volScalarField>& YActiveRef();
 };
 
 

applications/modules/multiphaseEuler/thermophysicalPredictor.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,24 +47,34 @@ void Foam::solvers::multiphaseEuler::compositionPredictor()
     {
         phaseModel& phase = fluid_.multicomponentPhases()[multicomponentPhasei];
 
-        UPtrList<volScalarField>& Y = phase.YActiveRef();
+        UPtrList<volScalarField>& Y = phase.YRef();
         const volScalarField& alpha = phase;
         const volScalarField& rho = phase.rho();
 
         forAll(Y, i)
         {
-            fvScalarMatrix YiEqn
+            if (phase.solveSpecie(i))
-            (
+            {
-                phase.YiEqn(Y[i])
+                fvScalarMatrix YiEqn
-             ==
+                (
-               *specieTransfer[Y[i].name()]
+                    phase.YiEqn(Y[i])
-              + fvModels().source(alpha, rho, Y[i])
+                 ==
-            );
+                   *specieTransfer[Y[i].name()]
+                  + fvModels().source(alpha, rho, Y[i])
+                );
 
-            YiEqn.relax();
+                YiEqn.relax();
-            fvConstraints().constrain(YiEqn);
+
-            YiEqn.solve("Yi");
+                fvConstraints().constrain(YiEqn);
-            fvConstraints().constrain(Y[i]);
+
+                YiEqn.solve("Yi");
+
+                fvConstraints().constrain(Y[i]);
+            }
+            else
+            {
+                Y[i].correctBoundaryConditions();
+            }
         }
     }