From a88d0e718118f21b2c80d1fe890a86e3e723466f Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Wed, 1 Sep 2021 12:24:36 -0700
Subject: [PATCH] BUG: Fixing BC's fixed value for species in multicomponent
 system

---
 .../MultiComponentPhaseModel.C                | 25 ++++++++++++++++---
 1 file changed, 22 insertions(+), 3 deletions(-)

diff --git a/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
index 8e10cfbd39..1917b21473 100644
--- a/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
+++ b/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -90,7 +90,9 @@ MultiComponentPhaseModel
                 (
                     IOobject::groupName("X" + species_[i], phaseName),
                     fluid.mesh().time().timeName(),
-                    fluid.mesh()
+                    fluid.mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
                 ),
                 Y()[i]
             )
@@ -116,7 +118,7 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel, phaseThermo>
     {
         const dimensionedScalar Wi
         (
-            "W",
+            "Wi",
             dimMass/dimMoles,
             thermo().composition().W(i)
         );
@@ -126,6 +128,23 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel, phaseThermo>
             X_[i] = W*Y()[i]/Wi;
             Xtotal += X_[i];
             X_[i].correctBoundaryConditions();
+
+            const volScalarField::Boundary& YBf = Y()[i].boundaryField();
+
+            forAll(YBf, patchi)
+            {
+                const fvPatchScalarField& YPf = YBf[patchi];
+                if (YPf.fixesValue())
+                {
+                    scalarField& xbf = X_[i].boundaryFieldRef()[patchi];
+                    const scalarField& ybf = Y()[i].boundaryField()[patchi];
+                    const scalarField& Wbf = W.boundaryField()[patchi];
+                    forAll(xbf, facei)
+                    {
+                        xbf[facei] = Wbf[facei]*ybf[facei]/Wi.value();
+                    }
+                }
+            }
         }
     }
     X_[inertIndex_] = 1.0 - Xtotal;