From 70bcb2bc7ac3eba3a1d061a5577bd3f0e44e7952 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Wed, 5 Oct 2022 15:38:29 +0100 Subject: [PATCH] fvModels::waveForcing: Added a continuity correction to momentum equation to compensate for phase-change generated by the phase-fraction forcing. This removes spurious velocity generated around the wave. --- src/waves/fvModels/forcing/forcing.C | 64 +++++++++++--------- src/waves/fvModels/forcing/forcing.H | 9 +++ src/waves/fvModels/waveForcing/waveForcing.C | 34 ++++++++--- src/waves/fvModels/waveForcing/waveForcing.H | 4 +- 4 files changed, 71 insertions(+), 40 deletions(-) diff --git a/src/waves/fvModels/forcing/forcing.C b/src/waves/fvModels/forcing/forcing.C index 7c85681aef..579480637b 100644 --- a/src/waves/fvModels/forcing/forcing.C +++ b/src/waves/fvModels/forcing/forcing.C @@ -116,57 +116,63 @@ void Foam::fv::forcing::readCoeffs() } -Foam::tmp Foam::fv::forcing::forceCoeff() const +Foam::tmp Foam::fv::forcing::scale() const { - tmp tforceCoeff + tmp tscale ( - new volScalarField::Internal + volScalarField::Internal::New ( - IOobject - ( - typedName("forceCoeff"), - mesh().time().timeName(), - mesh() - ), + typedName("scale"), mesh(), - dimensionedScalar(lambda_.dimensions(), scale_.valid() ? 0 : 1) + dimensionedScalar(dimless, scale_.valid() ? 0 : 1) ) ); - scalarField& forceCoeff = tforceCoeff.ref(); + + scalarField& scale = tscale.ref(); forAll(origins_, i) { const vectorField& c = mesh().cellCentres(); const scalarField x((c - origins_[i]) & directions_[i]); - forceCoeff = max(forceCoeff, scale_->value(x)); + scale = max(scale, scale_->value(x)); } - forceCoeff *= lambda_.value(); - // Write out the force coefficient for debugging if (debug && mesh().time().writeTime()) { - volScalarField vForceCoeff - ( - IOobject - ( - typedName("forceCoeff"), - mesh().time().timeName(), - mesh() - ), - mesh(), - lambda_.dimensions(), - zeroGradientFvPatchField::typeName - ); - vForceCoeff.primitiveFieldRef() = forceCoeff; - vForceCoeff.correctBoundaryConditions(); - vForceCoeff.write(); + tscale->write(); + } + + return tscale; +} + + +Foam::tmp Foam::fv::forcing::forceCoeff +( + const volScalarField::Internal& scale +) const +{ + tmp tforceCoeff + ( + volScalarField::Internal::New(typedName("forceCoeff"), lambda_*scale) + ); + + // Write out the force coefficient for debugging + if (debug && mesh().time().writeTime()) + { + tforceCoeff->write(); } return tforceCoeff; } +Foam::tmp Foam::fv::forcing::forceCoeff() const +{ + return forceCoeff(scale()); +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::fv::forcing::forcing diff --git a/src/waves/fvModels/forcing/forcing.H b/src/waves/fvModels/forcing/forcing.H index 9f47498db9..f81a4d6291 100644 --- a/src/waves/fvModels/forcing/forcing.H +++ b/src/waves/fvModels/forcing/forcing.H @@ -80,6 +80,15 @@ protected: //- Non-virtual read void readCoeffs(); + //- Return the scale distribution + tmp scale() const; + + //- Return the force coefficient given the scale distribution + tmp forceCoeff + ( + const volScalarField::Internal& scale + ) const; + //- Return the force coefficient tmp forceCoeff() const; diff --git a/src/waves/fvModels/waveForcing/waveForcing.C b/src/waves/fvModels/waveForcing/waveForcing.C index fbdc045475..d3b9d5bbe6 100644 --- a/src/waves/fvModels/waveForcing/waveForcing.C +++ b/src/waves/fvModels/waveForcing/waveForcing.C @@ -29,6 +29,9 @@ License #include "fvmSup.H" #include "addToRunTimeSelectionTable.H" +#include "fvcDdt.H" +#include "fvcDiv.H" + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam @@ -56,7 +59,7 @@ Foam::fv::waveForcing::waveForcing liquidPhaseName_(coeffs().lookup("liquidPhase")), alphaName_(IOobject::groupName("alpha", liquidPhaseName_)), UName_(coeffs().lookupOrDefault("U", "U")), - forceCoeff_(this->forceCoeff()) + scale_(this->scale()) {} @@ -76,8 +79,10 @@ void Foam::fv::waveForcing::addSup { if (fieldName == alphaName_) { - eqn -= fvm::Sp(forceCoeff_(), eqn.psi()); - eqn += forceCoeff_()*alphaWaves_(); + const volScalarField::Internal forceCoeff(this->forceCoeff(scale_)); + + eqn -= fvm::Sp(forceCoeff, eqn.psi()); + eqn += forceCoeff*alphaWaves_(); } } @@ -91,34 +96,45 @@ void Foam::fv::waveForcing::addSup { if (fieldName == UName_) { - eqn -= fvm::Sp(rho*forceCoeff_(), eqn.psi()); - eqn += rho*forceCoeff_()*Uwaves_(); + const volScalarField::Internal forceCoeff(rho*this->forceCoeff(scale_)); + + eqn -= fvm::Sp(forceCoeff, eqn.psi()); + eqn += forceCoeff*Uwaves_(); + + const surfaceScalarField& rhoPhi = + mesh().lookupObject("rhoPhi"); + + eqn += fvm::Sp + ( + scale()*(fvc::ddt(rho)()() + fvc::div(rhoPhi)()()), + eqn.psi() + ); } } bool Foam::fv::waveForcing::movePoints() { - forceCoeff_ = this->forceCoeff(); + scale_ = this->scale(); return true; } void Foam::fv::waveForcing::topoChange(const polyTopoChangeMap&) { - forceCoeff_ = this->forceCoeff(); + scale_ = this->scale(); } void Foam::fv::waveForcing::mapMesh(const polyMeshMap& map) { - forceCoeff_ = this->forceCoeff(); + scale_ = this->scale(); } void Foam::fv::waveForcing::distribute(const polyDistributionMap&) { - forceCoeff_ = this->forceCoeff(); + scale_ = this->scale(); } diff --git a/src/waves/fvModels/waveForcing/waveForcing.H b/src/waves/fvModels/waveForcing/waveForcing.H index 7f10b80937..5b01a595e3 100644 --- a/src/waves/fvModels/waveForcing/waveForcing.H +++ b/src/waves/fvModels/waveForcing/waveForcing.H @@ -115,8 +115,8 @@ class waveForcing //- Velocity field calculated from the current wave form tmp Uwaves_; - //- Forcing coefficient field - tmp forceCoeff_; + //- Forcing coefficient scale field + tmp scale_; public: