From 0b321e3eea07dfb3ab8c5ee4f3f4024f8be5b58c Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Fri, 15 Sep 2023 10:35:59 +0100 Subject: [PATCH] isothermalFilm: Corrected impingement pressure transfer The previous implementation was dimensionally inconsistent and was missing a factor of the VbyA field. This change will, in most cases, reduce the total impingement pressure contribution. --- .../CloudFilmTransfer/CloudFilmTransfer.C | 12 +---- .../filmCloudTransfer/filmCloudTransfer.C | 52 +++++-------------- .../filmCloudTransfer/filmCloudTransfer.H | 9 ---- .../isothermalFilm/momentumPredictor.C | 11 ++-- 4 files changed, 22 insertions(+), 62 deletions(-) diff --git a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C index 9c52b1e837..19bfa85a7d 100644 --- a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C +++ b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C @@ -83,21 +83,12 @@ void Foam::CloudFilmTransfer::absorbInteraction const parcelThermo& thermo = static_cast&>(this->owner()).thermo(); - // Patch face normal - const vector& nf = pp.faceNormals()[facei]; - // Patch velocity const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; // Relative parcel velocity const vector Urel = p.U() - Up; - // Parcel normal velocity - const vector Un = nf*(Urel & nf); - - // Parcel tangential velocity - const vector Ut = Urel - Un; - const liquidProperties& liq = thermo.liquids().properties()[0]; // Local pressure @@ -107,8 +98,7 @@ void Foam::CloudFilmTransfer::absorbInteraction ( facei, mass, // mass - mass*Ut, // tangential momentum - mass*mag(Un), // impingement pressure + mass*Urel, // momentum mass*liq.Hs(pc, p.T()) // energy ); diff --git a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C index 5e5b847c86..7af2188488 100644 --- a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C +++ b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C @@ -75,7 +75,6 @@ Foam::wordList Foam::fv::filmCloudTransfer::addSupFields() const { return wordList { - "pi", film_.alpha.name(), film_.thermo.he().name(), film_.U.name() @@ -136,36 +135,6 @@ inline Foam::fv::filmCloudTransfer::CloudToFilmTransferRate } -void Foam::fv::filmCloudTransfer::addSup -( - fvMatrix& eqn, - const word& fieldName -) const -{ - if (debug) - { - Info<< type() << ": applying source to " << eqn.psi().name() << endl; - } - - // Droplet impingement pressure - if (fieldName == "pi") - { - eqn += - CloudToFilmTransferRate - ( - pressureFromCloud_, - dimPressure*dimVolume - ); - } - else - { - FatalErrorInFunction - << "Support for field " << fieldName << " is not implemented" - << exit(FatalError); - } -} - - void Foam::fv::filmCloudTransfer::addSup ( const volScalarField& rho, @@ -240,11 +209,20 @@ void Foam::fv::filmCloudTransfer::addSup Info<< type() << ": applying source to " << eqn.psi().name() << endl; } - eqn += CloudToFilmTransferRate(momentumFromCloud_, dimMomentum); - - if (ejection_.valid()) + if (fieldName == film_.U.name()) { - eqn -= fvm::Sp(alpha()*rho()*ejection_->rate(), eqn.psi()); + eqn += CloudToFilmTransferRate(momentumFromCloud_, dimMomentum); + + if (ejection_.valid()) + { + eqn -= fvm::Sp(alpha()*rho()*ejection_->rate(), eqn.psi()); + } + } + else + { + FatalErrorInFunction + << "Support for field " << fieldName << " is not implemented" + << exit(FatalError); } } @@ -260,13 +238,11 @@ void Foam::fv::filmCloudTransfer::resetFromCloudFields() { massFromCloud_.setSize(nCloudPatchFaces); momentumFromCloud_.setSize(nCloudPatchFaces); - pressureFromCloud_.setSize(nCloudPatchFaces); energyFromCloud_.setSize(nCloudPatchFaces); } massFromCloud_ = 0; momentumFromCloud_ = Zero; - pressureFromCloud_ = 0; energyFromCloud_ = 0; cloudFieldsTransferred_ = true; @@ -281,13 +257,11 @@ void Foam::fv::filmCloudTransfer::parcelFromCloud const label facei, const scalar mass, const vector& momentum, - const scalar pressure, const scalar energy ) { massFromCloud_[facei] += mass; momentumFromCloud_[facei] += momentum; - pressureFromCloud_[facei] += pressure; energyFromCloud_[facei] += energy; } diff --git a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H index e7afbaf6b1..a1e58e61b8 100644 --- a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H +++ b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H @@ -75,7 +75,6 @@ class filmCloudTransfer scalarField massFromCloud_; vectorField momentumFromCloud_; - scalarField pressureFromCloud_; scalarField energyFromCloud_; //- Switch to ensure the ejection rate is not updated until @@ -147,13 +146,6 @@ public: // Add explicit and implicit contributions to compressible equation - //- Add explicit droplet impingement contribution to pressure field - virtual void addSup - ( - fvMatrix& eqn, - const word& fieldName - ) const; - //- Add explicit contribution to phase continuity virtual void addSup ( @@ -193,7 +185,6 @@ public: const label facei, const scalar mass, const vector& momentum, - const scalar pressure, const scalar energy ); diff --git a/applications/modules/isothermalFilm/momentumPredictor.C b/applications/modules/isothermalFilm/momentumPredictor.C index 0af820d6dd..db3cdbd708 100644 --- a/applications/modules/isothermalFilm/momentumPredictor.C +++ b/applications/modules/isothermalFilm/momentumPredictor.C @@ -76,8 +76,9 @@ Foam::solvers::isothermalFilm::pe() const // Update the pressure, mapping from the fluid region as required p.correctBoundaryConditions(); - // Add the droplet impingement pressure - p.ref() += mesh.time().deltaT()*fvModels().source(p, "pi")().Su(); + // Add the pressure caused normal momentum sources (e.g., parcels impinging + // with a normal velocity) + p.ref() += VbyA*(nHat & (fvModels().source(alpha, rho, U) & U)); return p; } @@ -90,6 +91,10 @@ void Foam::solvers::isothermalFilm::momentumPredictor() // Calculate the surface tension coefficient const volScalarField sigma(this->sigma()); + // Get the momentum source and remove any normal components + fvVectorMatrix alphaRhoUsource(fvModels().source(alpha, rho, U)); + alphaRhoUsource.source() -= nHat*(nHat & alphaRhoUsource.source()); + tUEqn = ( fvm::ddt(alpha, rho, U) + fvm::div(alphaRhoPhi, U) @@ -97,7 +102,7 @@ void Foam::solvers::isothermalFilm::momentumPredictor() + momentumTransport->divDevTau(U) == contactForce(sigma) - + fvModels().source(alpha, rho, U) + + alphaRhoUsource ); fvVectorMatrix& UEqn = tUEqn.ref();