diff --git a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C index ba5818bb96..883d52b002 100644 --- a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C +++ b/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::fv::multiphaseMangrovesSource::dragCoeff(const volVectorField& U) const +{ + tmp tdragCoeff = tmp::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::fv::multiphaseMangrovesSource::inertiaCoeff() const +{ + tmp tinertiaCoeff = tmp::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& eqn, + const label fieldi +) +{ + const volVectorField& U = eqn.psi(); + + fvMatrix 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 mangrovesEqn ( - IOobject - ( - 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 MangrovesEqn - ( - - fvm::Sp(rho*DragForceMangroves, U) - - rho*InertiaForceMangroves*fvm::ddt(rho, U) + - fvm::Sp(rho*dragCoeff(U), U) + - rho*inertiaCoeff()*fvm::ddt(U) ); // Contributions are added to RHS of momentum equation - eqn += MangrovesEqn; + eqn += mangrovesEqn; } diff --git a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.H b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.H index 0e8ec5b03a..217eb3f16d 100644 --- a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.H +++ b/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 dragCoeff(const volVectorField& U) const; - //- Disallow default bitwise copy construct - multiphaseMangrovesSource(const multiphaseMangrovesSource&); - - //- Disallow default bitwise assignment - void operator=(const multiphaseMangrovesSource&); + //- Return the inertia force coefficient + tmp inertiaCoeff() const; public: @@ -126,6 +132,13 @@ public: // Add explicit and implicit contributions + //- Add implicit contribution to momentum equation + virtual void addSup + ( + fvMatrix& eqn, + const label fieldi + ); + //- Add implicit contribution to phase momentum equation virtual void addSup ( diff --git a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C index ddbfce77fe..59f50d3a30 100644 --- a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C +++ b/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::fv::multiphaseMangrovesTurbulenceModel::kCoeff +( + const volVectorField& U +) const +{ + tmp tkCoeff = tmp::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::fv::multiphaseMangrovesTurbulenceModel::epsilonCoeff +( + const volVectorField& U +) const +{ + tmp tepsilonCoeff = tmp::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,82 +175,48 @@ void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup { const volVectorField& U = mesh_.lookupObject(UName_); - // Set alpha as zero-gradient - const volScalarField& epsilon = - mesh_.lookupObject(epsilonName_); - const volScalarField& k = mesh_.lookupObject(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 epsilonEqn ( - - fvm::Sp(epsilonMangroves, epsilon) - ); - eqn += epsilonEqn; + - fvm::Sp(epsilonCoeff(U), eqn.psi()) + ); + eqn += epsilonEqn; } else if (eqn.psi().name() == kName_) { - kMangroves.correctBoundaryConditions(); fvMatrix kEqn ( - - fvm::Sp(kMangroves, k) - ); + - fvm::Sp(kCoeff(U), eqn.psi()) + ); + eqn += kEqn; + } +} + + +void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldi +) +{ + const volVectorField& U = mesh_.lookupObject(UName_); + + if (eqn.psi().name() == epsilonName_) + { + fvMatrix epsilonEqn + ( + - fvm::Sp(rho*epsilonCoeff(U), eqn.psi()) + ); + eqn += epsilonEqn; + } + else if (eqn.psi().name() == kName_) + { + fvMatrix 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& eqn, - const label fieldi -) -{ - const volVectorField& U = mesh_.lookupObject(UName_); - - // Set alpha as zero-gradient - const volScalarField& epsilon = - mesh_.lookupObject(epsilonName_); - const volScalarField& k = mesh_.lookupObject(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 epsilonEqn - ( - - fvm::Sp(rho*epsilonMangroves, epsilon) - ); - eqn += epsilonEqn; - } - else if (eqn.psi().name() == "k") - { - Info<< " applying source to k stuff.... " << endl; - kMangroves.correctBoundaryConditions(); - fvMatrix kEqn - ( - - fvm::Sp(rho*kMangroves, k) - ); - eqn += kEqn; - } - -} - -void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup -( - const volScalarField& alpha, - const volScalarField& rho, - fvMatrix& eqn, - const label fieldi -) -{ - NotImplemented; -} // ************************************************************************* // diff --git a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.H b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.H index 7131062183..00f3b54c47 100644 --- a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.H +++ b/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 kCoeff(const volVectorField& U) const; - //- Disallow default bitwise copy construct - multiphaseMangrovesTurbulenceModel - ( - const multiphaseMangrovesTurbulenceModel& - ); - - //- Disallow default bitwise assignment - void operator=(const multiphaseMangrovesTurbulenceModel&); + //- Return the epsilon coefficient + tmp epsilonCoeff(const volVectorField& U) const; 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& eqn, - const label fieldi - ); // IO