INT: Refactored waves mangrove interaction fvOptions

src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2018 OpenCFD Ltd.
      \\/     M anipulation  | Copyright (C) 2017 IH-Cantabria
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,100 @@ namespace fv
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::fv::multiphaseMangrovesSource::dragCoeff(const volVectorField& U) const
+{
+    tmp<volScalarField> tdragCoeff = tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            typeName + ":dragCoeff",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("0", dimless/dimTime, 0)
+    );
+
+    volScalarField& dragCoeff = tdragCoeff.ref();
+
+    forAll(zoneIDs_, i)
+    {
+        const scalar a = aZone_[i];
+        const scalar N = NZone_[i];
+        const scalar Cd = CdZone_[i];
+
+        const labelList& zones = zoneIDs_[i];
+
+        for (label zonei : zones)
+        {
+            const cellZone& cz = mesh_.cellZones()[zonei];
+
+            for (label celli : cz)
+            {
+                const vector& Uc = U[celli];
+
+                dragCoeff[celli] = 0.5*Cd*a*N*mag(Uc);
+            }
+        }
+    }
+
+    dragCoeff.correctBoundaryConditions();
+
+    return tdragCoeff;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::fv::multiphaseMangrovesSource::inertiaCoeff() const
+{
+    tmp<volScalarField> tinertiaCoeff = tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            typeName + ":inertiaCoeff",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("0", dimless, 0)
+    );
+
+    volScalarField& inertiaCoeff = tinertiaCoeff.ref();
+
+    const scalar pi = constant::mathematical::pi;
+
+    forAll(zoneIDs_, i)
+    {
+        const scalar a = aZone_[i];
+        const scalar N = NZone_[i];
+        const scalar Cm = CmZone_[i];
+
+        const labelList& zones = zoneIDs_[i];
+
+        for (label zonei : zones)
+        {
+            const cellZone& cz = mesh_.cellZones()[zonei];
+
+            for (label celli : cz)
+            {
+                inertiaCoeff[celli] = 0.25*(Cm+1)*pi*a*a*N;
+            }
+        }
+    }
+
+    inertiaCoeff.correctBoundaryConditions();
+
+    return tinertiaCoeff;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::multiphaseMangrovesSource::multiphaseMangrovesSource
@@ -70,6 +164,25 @@ Foam::fv::multiphaseMangrovesSource::multiphaseMangrovesSource
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::fv::multiphaseMangrovesSource::addSup
+(
+    fvMatrix<vector>& eqn,
+    const label fieldi
+)
+{
+    const volVectorField& U = eqn.psi();
+
+    fvMatrix<vector> mangrovesEqn
+    (
+      - fvm::Sp(dragCoeff(U), U)
+      - inertiaCoeff()*fvm::ddt(U)
+    );
+
+    // Contributions are added to RHS of momentum equation
+    eqn += mangrovesEqn;
+}
+
+
 void Foam::fv::multiphaseMangrovesSource::addSup
 (
     const volScalarField& rho,
@@ -79,72 +192,14 @@ void Foam::fv::multiphaseMangrovesSource::addSup
 {
     const volVectorField& U = eqn.psi();
 
-    volScalarField DragForceMangroves
+    fvMatrix<vector> mangrovesEqn
     (
-        IOobject
+      - fvm::Sp(rho*dragCoeff(U), U)
-        (
+      - rho*inertiaCoeff()*fvm::ddt(U)
-            typeName + ":DragForceMangroves",
-            mesh_.time().timeName(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("DragForceMangroves", dimDensity/dimTime, 0)
-    );
-
-    volScalarField InertiaForceMangroves
-    (
-        IOobject
-        (
-            typeName + ":InertiaForceMangroves",
-            mesh_.time().timeName(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("InertiaForceMangroves", dimless, 0)
-    );
-
-    const scalar pi = constant::mathematical::pi;
-
-    forAll(zoneIDs_, i)
-    {
-        const labelList& zones = zoneIDs_[i];
-
-        forAll(zones, j)
-        {
-            const label zonei = zones[j];
-            const cellZone& cz = mesh_.cellZones()[zonei];
-
-            forAll(cz, k)
-            {
-                const label celli = cz[k];
-                const scalar a = aZone_[i];
-                const scalar N = NZone_[i];
-                const scalar Cm = CmZone_[i];
-                const scalar Cd = CdZone_[i];
-
-                const vector& Uc = U[celli];
-
-                DragForceMangroves[celli] = 0.5*Cd*a*N*mag(Uc);
-                InertiaForceMangroves[celli] = 0.25*(Cm+1)*pi*a*a;
-            }
-        }
-    }
-
-    DragForceMangroves.correctBoundaryConditions();
-    InertiaForceMangroves.correctBoundaryConditions();
-
-    fvMatrix<vector> MangrovesEqn
-    (
-      - fvm::Sp(rho*DragForceMangroves, U)
-      - rho*InertiaForceMangroves*fvm::ddt(rho, U)
     );
 
     // Contributions are added to RHS of momentum equation
-    eqn += MangrovesEqn;
+    eqn += mangrovesEqn;
 }
 
 

src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.H

@@ -61,6 +61,15 @@ class multiphaseMangrovesSource
 :
     public option
 {
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        multiphaseMangrovesSource(const multiphaseMangrovesSource&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const multiphaseMangrovesSource&);
+
+
 
 protected:
 
@@ -86,17 +95,14 @@ protected:
             //- Zone indices
             labelListList zoneIDs_;
 
-        // Field properties
 
-private:
+    // Protected Member Functions
 
-    // Private Member Functions
+        //- Return the drag force coefficient
+        tmp<volScalarField> dragCoeff(const volVectorField& U) const;
 
-        //- Disallow default bitwise copy construct
+        //- Return the inertia force coefficient
-        multiphaseMangrovesSource(const multiphaseMangrovesSource&);
+        tmp<volScalarField> inertiaCoeff() const;
-
-        //- Disallow default bitwise assignment
-        void operator=(const multiphaseMangrovesSource&);
 
 
 public:
@@ -126,6 +132,13 @@ public:
 
         // Add explicit and implicit contributions
 
+            //- Add implicit contribution to momentum equation
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldi
+            );
+
             //- Add implicit contribution to phase momentum equation
             virtual void addSup
             (

src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2018 OpenCFD Ltd.
      \\/     M anipulation  | Copyright (C) 2017 IH-Cantabria
 -------------------------------------------------------------------------------
 License
@@ -27,13 +27,7 @@ License
 #include "mathematicalConstants.H"
 #include "fvMesh.H"
 #include "fvMatrices.H"
-#include "fvmDdt.H"
-#include "fvcGrad.H"
-#include "fvmDiv.H"
-#include "fvmLaplacian.H"
 #include "fvmSup.H"
-#include "surfaceInterpolate.H"
-#include "surfaceFields.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -53,6 +47,100 @@ namespace fv
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::fv::multiphaseMangrovesTurbulenceModel::kCoeff
+(
+    const volVectorField& U
+) const
+{
+    tmp<volScalarField> tkCoeff = tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            typeName + ":kCoeff",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("0", dimless/dimTime, 0)
+    );
+
+    volScalarField& kCoeff = tkCoeff.ref();
+
+    forAll(zoneIDs_, i)
+    {
+        const scalar a = aZone_[i];
+        const scalar N = NZone_[i];
+        const scalar Ckp = CkpZone_[i];
+        const scalar Cd = CdZone_[i];
+
+        for (label zonei : zoneIDs_[i])
+        {
+            const cellZone& cz = mesh_.cellZones()[zonei];
+
+            for (label celli : cz)
+            {
+                kCoeff[celli] = Ckp*Cd*a*N*mag(U[celli]);
+            }
+        }
+    }
+
+    kCoeff.correctBoundaryConditions();
+
+    return tkCoeff;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::fv::multiphaseMangrovesTurbulenceModel::epsilonCoeff
+(
+    const volVectorField& U
+) const
+{
+    tmp<volScalarField> tepsilonCoeff = tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            typeName + ":epsilonCoeff",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("0", dimless/dimTime, 0)
+    );
+
+    volScalarField& epsilonCoeff = tepsilonCoeff.ref();
+
+    forAll(zoneIDs_, i)
+    {
+        const scalar a = aZone_[i];
+        const scalar N = NZone_[i];
+        const scalar Cep = CepZone_[i];
+        const scalar Cd = CdZone_[i];
+
+        for (label zonei : zoneIDs_[i])
+        {
+            const cellZone& cz = mesh_.cellZones()[zonei];
+
+            for (label celli : cz)
+            {
+                epsilonCoeff[celli] = Cep*Cd*a*N*mag(U[celli]);
+            }
+        }
+    }
+
+    epsilonCoeff.correctBoundaryConditions();
+
+    return tepsilonCoeff;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::multiphaseMangrovesTurbulenceModel::multiphaseMangrovesTurbulenceModel
@@ -87,81 +175,47 @@ void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup
 {
     const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
 
-    // Set alpha as zero-gradient
-    const volScalarField& epsilon =
-        mesh_.lookupObject<volScalarField>(epsilonName_);
-    const volScalarField& k = mesh_.lookupObject<volScalarField>(kName_);
-
-    volScalarField epsilonMangroves
-    (
-        IOobject
-        (
-            typeName + ":epsilonMangroves",
-            mesh_.time().timeName(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("epsilonMangroves", dimMass/dimMass/dimTime, 0)
-    );
-
-    volScalarField kMangroves
-    (
-        IOobject
-        (
-            typeName + ":kMangroves",
-            mesh_.time().timeName(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("kMangroves", dimTime/dimTime/dimTime, 0)
-    );
-
-    forAll(zoneIDs_, i)
-    {
-        const labelList& zones = zoneIDs_[i];
-
-        forAll(zones, j)
-        {
-            const label zonei = zones[j];
-            const cellZone& cz = mesh_.cellZones()[zonei];
-
-            forAll(cz, kk)
-            {
-                const label celli = cz[kk];
-                const scalar a = aZone_[i];
-                const scalar N = NZone_[i];
-                const scalar Ckp = CkpZone_[i];
-                const scalar Cep = CepZone_[i];
-                const scalar Cd = CdZone_[i];
-
-                const vector& Uc = U[celli];
-
-                const scalar f = Cd*a*N*mag(Uc);
-                epsilonMangroves[celli] = Cep*f;
-                kMangroves[celli] = Ckp*f;
-            }
-        }
-    }
-
     if (eqn.psi().name() == epsilonName_)
     {
-    	epsilonMangroves.correctBoundaryConditions();
 	    fvMatrix<scalar> epsilonEqn
 	    (
-      	    - fvm::Sp(epsilonMangroves, epsilon)
+          - fvm::Sp(epsilonCoeff(U), eqn.psi())
         );
         eqn += epsilonEqn;
     }
     else if (eqn.psi().name() == kName_)
     {
-    	kMangroves.correctBoundaryConditions();
 	    fvMatrix<scalar> kEqn
 	    (
-      	    - fvm::Sp(kMangroves, k)
+          - fvm::Sp(kCoeff(U), eqn.psi())
+        );
+        eqn += kEqn;
+    }
+}
+
+
+void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
+
+    if (eqn.psi().name() == epsilonName_)
+    {
+	    fvMatrix<scalar> epsilonEqn
+	    (
+          - fvm::Sp(rho*epsilonCoeff(U), eqn.psi())
+        );
+        eqn += epsilonEqn;
+    }
+    else if (eqn.psi().name() == kName_)
+    {
+	    fvMatrix<scalar> kEqn
+	    (
+          - fvm::Sp(rho*kCoeff(U), eqn.psi())
         );
 	    eqn += kEqn;
     }
@@ -230,106 +284,5 @@ bool Foam::fv::multiphaseMangrovesTurbulenceModel::read(const dictionary& dict)
     }
 }
 
-void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup
-(
-    const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
-    const label fieldi
-)
-{
-    const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
-
-    // Set alpha as zero-gradient
-    const volScalarField& epsilon =
-        mesh_.lookupObject<volScalarField>(epsilonName_);
-    const volScalarField& k = mesh_.lookupObject<volScalarField>(kName_);
-
-    volScalarField epsilonMangroves
-    (
-        IOobject
-        (
-            typeName + ":epsilonMangroves",
-            mesh_.time().timeName(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("epsilonMangroves", dimMass/dimMass/dimTime, 0)
-    );
-
-    volScalarField kMangroves
-    (
-        IOobject
-        (
-            typeName + ":kMangroves",
-            mesh_.time().timeName(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("kMangroves", dimTime/dimTime/dimTime, 0)
-    );
-
-    forAll(zoneIDs_, i)
-    {
-        const labelList& zones = zoneIDs_[i];
-
-        forAll(zones, j)
-        {
-            const label zonei = zones[j];
-            const cellZone& cz = mesh_.cellZones()[zonei];
-
-            forAll(cz, kk)
-            {
-                const label celli = cz[kk];
-                const scalar a = aZone_[i];
-                const scalar N = NZone_[i];
-                const scalar Ckp = CkpZone_[i];
-                const scalar Cep = CepZone_[i];
-                const scalar Cd = CdZone_[i];
-
-                const vector& Uc = U[celli];
-
-    epsilonMangroves[celli] = Cep * Cd * a * N * sqrt ( Uc.x()*Uc.x() + Uc.y()*Uc.y() + Uc.z()*Uc.z());
-    kMangroves[celli] = Ckp * Cd * a * N * sqrt ( Uc.x()*Uc.x() + Uc.y()*Uc.y() + Uc.z()*Uc.z());
-            }
-        }
-    }
-
-    if (eqn.psi().name() == "epsilon")
-    {
-	Info<< " applying source to epsilon stuff.... " << endl;
-    	epsilonMangroves.correctBoundaryConditions();
-	fvMatrix<scalar> epsilonEqn
-	(
-      	    - fvm::Sp(rho*epsilonMangroves, epsilon)
-    	);
-	eqn += epsilonEqn;
-    }
-    else if (eqn.psi().name() == "k")
-    {
-	Info<< " applying source to k stuff.... " << endl;
-    	kMangroves.correctBoundaryConditions();
-	fvMatrix<scalar> kEqn
-	(
-      	    - fvm::Sp(rho*kMangroves, k)
-    	);
-	eqn += kEqn;
-    }
-
-}
-
-void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup
-(
-    const volScalarField& alpha,
-    const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
-    const label fieldi
-)
-{
-    NotImplemented;
-}
 
 // ************************************************************************* //

src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.H

@@ -59,6 +59,18 @@ class multiphaseMangrovesTurbulenceModel
 :
     public option
 {
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        multiphaseMangrovesTurbulenceModel
+        (
+            const multiphaseMangrovesTurbulenceModel&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const multiphaseMangrovesTurbulenceModel&);
+
+
 
 protected:
 
@@ -97,18 +109,13 @@ protected:
             word epsilonName_;
 
 
-private:
+    // Protected Member Functions
 
-    // Private Member Functions
+        //- Return the k coefficient
+        tmp<volScalarField> kCoeff(const volVectorField& U) const;
 
-        //- Disallow default bitwise copy construct
+        //- Return the epsilon coefficient
-        multiphaseMangrovesTurbulenceModel
+        tmp<volScalarField> epsilonCoeff(const volVectorField& U) const;
-        (
-            const multiphaseMangrovesTurbulenceModel&
-        );
-
-        //- Disallow default bitwise assignment
-        void operator=(const multiphaseMangrovesTurbulenceModel&);
 
 
 public:
@@ -153,14 +160,6 @@ public:
                 const label fieldi
             );
 
-            //- Add implicit contribution to phase momentum equation
-            virtual void addSup
-            (
-                const volScalarField& alpha,
-                const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const label fieldi
-            );
 
         // IO