From c5027f53ba58f453e4d3aad4464e4d5aaf88e6db Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Tue, 10 Jul 2018 16:17:12 +0900
Subject: [PATCH] D, mu and k boundary fields update in correct()

---
 .../diffusivityModel/diffusivityModel.C       | 397 ++++++++++++------
 .../diffusivityModel/diffusivityModel.H       |  22 +-
 2 files changed, 286 insertions(+), 133 deletions(-)

diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index 49fbb20..658afde 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -57,11 +57,68 @@ Foam::scalar Foam::diffusivityModel::mixAvgD(const label k, const UList<scalar>&
 
     XD[k] = 0.0;
     YD[k] = 0.0;
-    
+
     return 1. / (sum(XD) + XY*sum(YD));
 }
 
 
+void
+Foam::diffusivityModel::mixAvgDi(UList<scalar>& Di, const scalarSymmetricSquareMatrix &Dij, const UList<scalar>& X, const UList<scalar>& Y)
+{
+
+    // Pure Condition Test ( Xi = 1 )
+    Switch pure = false;
+    label pureSpecieI = -1;
+    forAll (X, i)
+    {
+        if (1. - X[i] < 1e-12)
+        {
+            pureSpecieI = i;
+            pure = true;
+            break;
+        }
+    }
+
+    // Mixture average D, mu and k
+    if (pure)
+    {
+        const scalar* Dpure = Dij[pureSpecieI];
+
+        forAll (X, i)
+        {
+            Di[i] = Dpure[i];
+        }
+    }
+    else
+    {
+        forAll (X, k)
+        {
+            scalar sum1 = 0.0;
+            scalar sum2 = 0.0;
+
+            forAll (X, i)
+            {
+                if (i == k)
+                {
+                    continue;
+                }
+
+                const scalar* Dk = Dij[k];
+
+                sum1 += X[i] / Dk[i];
+                sum2 += Y[i] / Dk[i];
+            }
+
+            sum2 *= X[k] / (1. - Y[k]);
+
+            Di[k] = inv(sum1 + sum2);
+        }
+    }
+
+    return;
+}
+
+
 Foam::scalar Foam::diffusivityModel::mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W)
 {
     const label K = muI.size();
@@ -105,6 +162,122 @@ Foam::scalar Foam::diffusivityModel::mixAvgK(const UList<scalar>& k, const UList
 }
 
 
+void Foam::diffusivityModel::calculateMuD
+(
+ UList<scalar> &muI,
+ scalarSymmetricSquareMatrix &Dij,
+ const scalar rhoi,
+ const scalar Cpi,
+ const scalar pi,
+ const scalar Ti,
+ const scalar WbarI,
+ const scalar rhoQc2i,
+ const scalarField &localY,
+ const scalarField &localX
+ )
+{
+
+    const label ionStart = 1;
+
+    const label neutralStart = 1 + ions_.size();
+
+    // Electron - Electron
+    label idx = 0;
+
+    muI[0] = ecs_[idx].mu(pi,Ti,rhoQc2i);
+
+    Dij(0,0) = ecs_[idx].D(pi,Ti,rhoQc2i);
+
+    // Electron - Ions
+    idx = 1;
+    forAll (ions_, J)
+    {
+        label j = J + ionStart;
+
+        // Calculate Dij
+        Dij(0,j) = ecs_[idx].D(pi,Ti,rhoQc2i);
+        Dij(j,0) = Dij(0,j);
+        idx++;
+    }
+
+    // Electron - Neutrals
+    idx = 0;
+    forAll (neutrals_, J)
+    {
+        label j = J + neutralStart;
+
+        // Calculate Dij
+        Dij(0,j) = ens_[idx].D(pi,Ti);
+        Dij(j,0) = Dij(0,j);
+        idx++;
+    }
+
+    // Ions
+    label idx1 = 0;
+    forAll (ions_, I)
+    {
+        label i = I + ionStart;
+
+        muI[i] = ccs_[idx1].mu(pi,Ti,rhoQc2i);
+
+        Dij(i,i) = ccs_[idx1].D(pi,Ti,rhoQc2i);
+
+        idx1++;
+
+        // Ion - Ions
+        for (label J = I + 1; J < ions_.size(); J++)
+        {
+            label j = J + ionStart;
+            // Calculate Dij
+            Dij(i,j) = ccs_[idx1].D(pi,Ti,rhoQc2i);
+            Dij(j,i) = Dij(i,j);
+            idx1++;
+        }
+
+    }
+
+    label idx2 = 0;
+    forAll (ions_, I)
+    {
+        label i = I + ionStart;
+
+        // Ion - Neutrals
+        forAll (neutrals_, J)
+        {
+            label j = J + neutralStart;
+
+            // Calculate Dij
+            Dij(i,j) = cns_[idx2].D(pi,Ti);
+            Dij(j,i) = Dij(i,j);
+            idx2++;
+        }
+    }
+
+    // Neutrals
+    idx = 0;
+    forAll (neutrals_, I)
+    {
+        label i = I + neutralStart;
+
+        muI[i] = nns_[idx].mu(pi,Ti);
+        Dij(i,i) = nns_[idx].D(pi,Ti);
+        idx++;
+
+        // Neutral - Neutrals
+        for (label J = I + 1; J < neutrals_.size(); J++)
+        {
+            label j = J + neutralStart;
+
+            // Calculate Dij
+            Dij(i,j) = nns_[idx].D(pi,Ti);
+            Dij(j,i) = Dij(i,j);
+            idx++;
+        }
+    }
+
+
+}
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 
@@ -388,6 +561,7 @@ Foam::diffusivityModel::~diffusivityModel()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+
 void Foam::diffusivityModel::correct()
 {
     const volScalarField &T = thermo_.T();
@@ -417,6 +591,7 @@ void Foam::diffusivityModel::correct()
 
     scalarField localD(species_.size());
 
+    scalarField Di(species_.size());
     scalarField Dii(species_.size());
     scalarField muI(species_.size());
     scalarField kI(species_.size());
@@ -449,10 +624,7 @@ void Foam::diffusivityModel::correct()
         const scalar WbarI = Wbar[celli];
         const scalar rhoQc2i = rhoQc2[celli];
 
-        forAll (species_, i)
-        {
-            localCv[i] = thermo_.composition().Cv(i,pi,Ti);
-        }
+        thermo_.composition().Cv(localCv,pi,Ti);
 
         forAll (species_, i)
         {
@@ -460,111 +632,24 @@ void Foam::diffusivityModel::correct()
             localX[i] = Y[i][celli] * WbarI / Wpure[i];
         }
 
+        calculateMuD ( muI, Dij, rhoi, Cpi, pi, Ti, WbarI, rhoQc2i, localY, localX);
+
+        forAll (Dii, i)
+        {
+            Dii[i] = Dij(i,i);
+        }
+
+
         Zrot[0] = 0;
         forAll (ions_, i)
         {
-            // ZrotI[i] = ions_[i].Zrot(Ti);
-            ZrotI[i] = 0.0;
+            ZrotI[i] = ions_[i].Zrot(Ti);
         }
         forAll (neutrals_, i)
         {
             ZrotN[i] = neutrals_[i].Zrot(Ti);
         }
 
-        // Electron - Electron
-        label idx = 0;
-
-        muI[0] = ecs_[idx].mu(pi,Ti,rhoQc2i);
-        // muI[0] = SMALL;
-
-        Dii[0] = ecs_[idx].D(pi,Ti,rhoQc2i);
-        // Dii[0] = SMALL;
-        Dij(0,0) = Dii[0];
-
-        // Electron - Ions
-        idx = 1;
-        for (label J = 0; J < ions_.size(); J++)
-        {
-            label j = J + 1;
-            // Calculate Dij
-            Dij(0,j) = ecs_[idx].D(pi,Ti,rhoQc2i);
-            // Dij(0,j) = SMALL;
-            Dij(j,0) = Dij(0,j);
-            idx++;
-        }
-
-        // Electron - Neutrals
-        idx = 0;
-        for (label J = 0; J < neutrals_.size(); J++)
-        {
-            label j = J + ions_.size() + 1;
-            // Calculate Dij
-            Dij(0,j) = ens_[idx].D(pi,Ti);
-            // Dij(0,j) = SMALL;
-            Dij(j,0) = Dij(0,j);
-            idx++;
-        }
-
-        // Ions
-        label idx1 = 0;
-        label idx2 = 0;
-        forAll (ions_, I)
-        {
-            label i = I + 1;
-
-            muI[i] = ccs_[idx1].mu(pi,Ti,rhoQc2i);
-            // muI[i] = SMALL;
-
-            Dii[i] = ccs_[idx1].D(pi,Ti,rhoQc2i);
-            // Dii[i] = SMALL;
-            Dij(i,i) = Dii[i];
-            idx1++;
-
-            // Ion - Ions
-            for (label J = I+1; J < ions_.size(); J++)
-            {
-                label j = J + 1;
-                // Calculate Dij
-                Dij(i,j) = ccs_[idx1].D(pi,Ti,rhoQc2i);
-                // Dij(i,j) = SMALL;
-                Dij(j,i) = Dij(i,j);
-                idx1++;
-            }
-
-            // Ion - Neutrals
-            for (label j = ions_.size()+1; j < species_.size(); j++)
-            {
-                // Calculate Dij
-                Dij(i,j) = cns_[idx2].D(pi,Ti);
-                // Dij(i,j) = SMALL;
-                Dij(j,i) = Dij(i,j);
-                idx2++;
-            }
-        }
-
-        // Neutrals
-        idx = 0;
-        forAll (neutrals_, I)
-        {
-            label i = I + 1 + ions_.size();
-
-            muI[i] = nns_[idx].mu(pi,Ti);
-
-            Dii[i] = nns_[idx].D(pi,Ti);
-            Dij(i,i) = Dii[i];
-            idx++;
-
-            // Neutral - Neutrals
-            for (label j = i+1; j < species_.size(); j++)
-            {
-                // Calculate Dij
-                Dij(i,j) = nns_[idx].D(pi,Ti);
-                Dij(j,i) = Dij(i,j);
-                idx++;
-            }
-        }
-
-
         // Pure Thermal conductivity
         forAll (species_, i)
         {
@@ -587,42 +672,94 @@ void Foam::diffusivityModel::correct()
 
         }
 
-        // Pure Condition Test ( Xi = 1 )
-        Switch pure = false;
-        label pureSpecieI = -1;
-        forAll (species_, i)
+        mixAvgDi(Di, Dij, localX, localY);
+        forAll (Di, i)
         {
-            if (1. - localX[i] < 1e-12)
-            {
-                pureSpecieI = i;
-                pure = true;
-                break;
-            }
-        }
-
-        // Mixture average D, mu and k
-        if (pure)
-        {
-            UList<scalar> Dpure(Dij[pureSpecieI], Dij.n());
-
-            forAll (species_, i)
-            {
-                D_[i][celli] = Dpure[i];
-            }
-        }
-        else
-        {
-            forAll (species_, i)
-            {
-                UList<scalar> Di(Dij[i], Dij.n());
-                D_[i][celli] = mixAvgD(i, Di, localX, localY);
-            }
+            D_[i][celli] = Di[i];
         }
 
         mu_()[celli] = mixAvgMu(muI, localX, Wpure);
         k_()[celli] = mixAvgK(kI, localX) / rhoi / Cpi;
     }
 
+
+    forAll(p.boundaryField(), patchi)
+    {
+            const volScalarField::Patch &rhop = rho.boundaryField()[patchi];
+            const volScalarField::Patch &Cpp = Cp.boundaryField()[patchi];
+            const volScalarField::Patch &pp = p.boundaryField()[patchi];
+            const volScalarField::Patch &Tp = T.boundaryField()[patchi];
+            const volScalarField::Patch &Wbarp = Wbar.boundaryField()[patchi];
+            const volScalarField::Patch &rhoQc2p = rhoQc2.boundaryField()[patchi];
+
+        forAll(rhop, facei)
+        {
+            const scalar rhoi = rhop[facei];
+            const scalar Cpi = Cpp[facei];
+            const scalar pi = pp[facei];
+            const scalar Ti = Tp[facei];
+            const scalar WbarI = Wbarp[facei];
+            const scalar rhoQc2i = rhoQc2p[facei];
+
+            thermo_.composition().Cv(localCv,pi,Ti);
+
+            forAll (species_, i)
+            {
+                localY[i] = Y[i].boundaryField()[patchi][facei];
+                localX[i] = localY[i] * WbarI / Wpure[i];
+            }
+
+            calculateMuD ( muI, Dij, rhoi, Cpi, pi, Ti, WbarI, rhoQc2i, localY, localX);
+
+            forAll (Dii, i)
+            {
+                Dii[i] = Dij(i,i);
+            }
+
+
+            Zrot[0] = 0;
+            forAll (ions_, i)
+            {
+                ZrotI[i] = ions_[i].Zrot(Ti);
+            }
+            forAll (neutrals_, i)
+            {
+                ZrotN[i] = neutrals_[i].Zrot(Ti);
+            }
+
+            // Pure Thermal conductivity
+            forAll (species_, i)
+            {
+                const scalar R = Neutral::RR / Wpure[i];
+                const scalar CvTrans = (3./2.)*R;
+                const scalar CvRot = (dof[i]/2.)*R;
+                const scalar CvVib = localCv[i] - CvTrans - CvRot;
+
+                const scalar rSc = rhoi * Dii[i] / muI[i];
+
+                const scalar A = 5./2. - rSc;
+                const scalar B = Zrot[i] + (2./Neutral::pi) * ((5./3.)*CvRot + rSc);
+                const scalar AB = (2./Neutral::pi)*(A/B);
+
+                const scalar fTrans = (5./2.) * (1.0 - AB * CvRot / CvTrans);
+                const scalar fRot = rSc*(1.0+AB);
+                const scalar fVib = rSc;
+
+                kI[i] = muI[i]*(fTrans*CvTrans + fRot*CvRot + fVib*CvVib);
+
+            }
+
+            mixAvgDi(Di, Dij, localX, localY);
+            forAll (Di, i)
+            {
+                D_[i].boundaryFieldRef()[patchi][facei] = Di[i];
+            }
+
+            mu_().boundaryFieldRef()[patchi][facei] = mixAvgMu(muI, localX, Wpure);
+            k_().boundaryFieldRef()[patchi][facei] = mixAvgK(kI, localX) / rhoi / Cpi;
+        }
+    }
+
     return;
 }
 
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index d76f4aa..b866fd0 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -112,14 +112,30 @@ class diffusivityModel
         void operator=(const diffusivityModel&);
 
         //- Disallow default bitwise assignment
-        scalar mixAvgD(const label k, const UList<scalar>& Di, const UList<scalar>& X, const UList<scalar>& Y);
+        inline scalar mixAvgD(const label k, const UList<scalar>& Di, const UList<scalar>& X, const UList<scalar>& Y);
 
         //- Disallow default bitwise assignment
-        scalar mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W);
+        inline void mixAvgDi(UList<scalar>& Di, const scalarSymmetricSquareMatrix &Dij, const UList<scalar>& X, const UList<scalar>& Y);
 
         //- Disallow default bitwise assignment
-        scalar mixAvgK(const UList<scalar>& k, const UList<scalar>& X);
+        inline scalar mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W);
 
+        //- Disallow default bitwise assignment
+        inline scalar mixAvgK(const UList<scalar>& k, const UList<scalar>& X);
+
+        inline void calculateMuD
+        (
+            UList<scalar> &muI,
+            scalarSymmetricSquareMatrix &Dij,
+            const scalar rhoi,
+            const scalar Cpi,
+            const scalar pi,
+            const scalar Ti,
+            const scalar WbarI,
+            const scalar rhoQc2i,
+            const scalarField &localY,
+            const scalarField &localX
+        );
 
 public: