From 18200e72a61feb98d4729e2cff00b42e5ad4a616 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Fri, 1 Sep 2023 13:07:00 +0100 Subject: [PATCH] multiphaseEuler: Update the virtual-mass force implementation The central coefficient part of the virtual-mass phase acceleration matrix is now included in the phase velocity transport central coefficient + drag matrix so that the all the phase contributions to each phase momentum equation are handled implicitly and consistently without lagging contribution from the other phases in either the pressure equation or phase momentum correctors. This improves the conditioning of the pressure equation and convergence rate of bubbly-flow cases. --- .../functionObjects/phaseForces/phaseForces.C | 5 +- .../multiphaseEuler/cellPressureCorrector.C | 14 +- .../multiphaseEuler/facePressureCorrector.C | 14 +- .../MomentumTransferPhaseSystem.C | 482 ++++++++---------- .../MomentumTransferPhaseSystem.H | 18 +- .../MovingPhaseModel/MovingPhaseModel.C | 51 +- .../MovingPhaseModel/MovingPhaseModel.H | 14 +- .../StationaryPhaseModel.C | 30 +- .../StationaryPhaseModel.H | 8 +- .../phaseModel/phaseModel/phaseModel.H | 8 +- .../phaseSystems/phaseSystem/phaseSystem.H | 8 +- 11 files changed, 274 insertions(+), 378 deletions(-) diff --git a/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C b/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C index 43bafa0342..049df0d273 100644 --- a/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C +++ b/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C @@ -241,7 +241,10 @@ bool Foam::functionObjects::phaseForces::execute() *forceFields_[virtualMassModel::typeName] += fluid_.lookupInterfacialModel (interface).K() - *(otherPhase.DUDt() - phase_.DUDt()); + *( + (otherPhase.DUDt() & otherPhase.U()) + - (phase_.DUDt() & phase_.U()) + ); } if (fluid_.foundInterfacialModel(interface)) diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C index ac12ba5f77..11f48b6204 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C @@ -68,6 +68,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() rAs.setSize(phases.size()); } + PtrList HVms(movingPhases.size()); PtrList> invADs; PtrList> invADfs; { @@ -109,7 +110,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() } } - fluid.invADs(As, invADs, invADfs); + fluid.invADs(As, HVms, invADs, invADfs); } volScalarField rho("rho", fluid.rho()); @@ -168,7 +169,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() // --- Optional momentum predictor if (predictMomentum) { - InfoInFunction << endl; PtrList HbyADs; { PtrList Hs(movingPhases.size()); @@ -189,6 +189,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ) *phase.U()().oldTime() ); + + if (HVms.set(movingPhasei)) + { + Hs[movingPhasei] += HVms[movingPhasei]; + } } HbyADs = invADs & Hs; @@ -253,6 +258,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() *phase.U()().oldTime() ); + if (HVms.set(movingPhasei)) + { + Hs[movingPhasei] += HVms[movingPhasei]; + } + phiHs.set ( movingPhasei, diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C index 92960da404..05e5b07f35 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C @@ -68,9 +68,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() rAs.setSize(phases.size()); } + PtrList HVmfs(movingPhases.size()); PtrList> invADfs; { - PtrList Vmfs(fluid.Vmfs()); PtrList Afs(movingPhases.size()); forAll(fluid.movingPhases(), movingPhasei) @@ -110,14 +110,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() fvc::interpolate(A) ) ); - - if (Vmfs.set(phase.index())) - { - Afs[movingPhasei] += Vmfs[phase.index()]; - } } - invADfs = fluid.invADfs(Afs); + invADfs = fluid.invADfs(Afs, HVmfs); } volScalarField rho("rho", fluid.rho()); @@ -205,6 +200,11 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() + fvc::flux(UEqns[phase.index()].H()) ) ); + + if (HVmfs.set(movingPhasei)) + { + phiHs[movingPhasei] += HVmfs[movingPhasei]; + } } phiHbyADs = invADfs & phiHs; diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index b95b33f30e..b79499374b 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -224,87 +224,6 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() } } - // Initialise Vms table if required - if (!Vms_.size()) - { - forAllConstIter - ( - virtualMassModelTable, - virtualMassModels_, - virtualMassModelIter - ) - { - const phaseInterface& interface = - virtualMassModelIter()->interface(); - - Vms_.insert - ( - interface, - new volScalarField - ( - IOobject - ( - IOobject::groupName("Vm", interface.name()), - this->mesh().time().name(), - this->mesh() - ), - this->mesh(), - dimensionedScalar(virtualMassModel::dimK, 0) - ) - ); - } - } - - // Update the virtual mass coefficients - forAllConstIter - ( - virtualMassModelTable, - virtualMassModels_, - virtualMassModelIter - ) - { - *Vms_[virtualMassModelIter.key()] = virtualMassModelIter()->K(); - } - - // Add the virtual mass force - forAllConstIter(VmTable, Vms_, VmIter) - { - const volScalarField& Vm(*VmIter()); - const phaseInterface interface(*this, VmIter.key()); - - forAllConstIter(phaseInterface, interface, iter) - { - const phaseModel& phase = iter(); - const phaseModel& otherPhase = iter.otherPhase(); - - if (!phase.stationary()) - { - fvVectorMatrix& eqn = *eqns[phase.name()]; - - const volVectorField& U = eqn.psi(); - - const surfaceScalarField& phi = phase.phiRef(); - const tmp taphi(fvc::absolute(phi, U)); - const surfaceScalarField& aphi(taphi()); - - const volScalarField VmPhase - ( - (otherPhase/max(otherPhase, otherPhase.residualAlpha())) - *Vm - ); - - eqn -= - VmPhase - *( - fvm::ddt(U) - + fvm::div(aphi, U) - fvm::Sp(fvc::div(aphi), U) - - otherPhase.DUDt() - ) - + this->MRF_.DDt(VmPhase, U - otherPhase.U()); - } - } - } - return eqnsPtr; } @@ -373,138 +292,10 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() *Kdfs_[dragModelIter.key()] = dragModelIter()->Kf(); } - // Initialise Vms table if required - if (!Vms_.size()) - { - forAllConstIter - ( - virtualMassModelTable, - virtualMassModels_, - virtualMassModelIter - ) - { - const phaseInterface& interface = - virtualMassModelIter()->interface(); - - Vms_.insert - ( - interface, - new volScalarField - ( - IOobject - ( - IOobject::groupName("Vm", interface.name()), - this->mesh().time().name(), - this->mesh() - ), - this->mesh(), - dimensionedScalar(virtualMassModel::dimK, 0) - ) - ); - } - } - - // Update the virtual mass coefficients - forAllConstIter - ( - virtualMassModelTable, - virtualMassModels_, - virtualMassModelIter - ) - { - *Vms_[virtualMassModelIter.key()] = virtualMassModelIter()->K(); - } - - // Create U & grad(U) fields - PtrList UgradUs(this->phaseModels_.size()); - forAll(this->phaseModels_, phasei) - { - const phaseModel& phase = this->phaseModels_[phasei]; - - if (!phase.stationary()) - { - const volVectorField& U = phase.URef(); - const surfaceScalarField& phi = phase.phiRef(); - - const tmp taphi(fvc::absolute(phi, U)); - const surfaceScalarField& aphi(taphi()); - - UgradUs.set - ( - phasei, - new fvVectorMatrix - ( - fvm::div(aphi, U) - fvm::Sp(fvc::div(aphi), U) - + this->MRF().DDt(U) - ) - ); - } - } - - // Add the virtual mass force - forAllConstIter(VmTable, Vms_, VmIter) - { - const volScalarField& Vm(*VmIter()); - const phaseInterface interface(*this, VmIter.key()); - - forAllConstIter(phaseInterface, interface, iter) - { - const phaseModel& phase = iter(); - const phaseModel& otherPhase = iter.otherPhase(); - - if (!phase.stationary()) - { - const volScalarField VmPhase - ( - (otherPhase/max(otherPhase, otherPhase.residualAlpha())) - *Vm - ); - - *eqns[phase.name()] -= - VmPhase - *( - UgradUs[phase.index()] - - (UgradUs[otherPhase.index()] & otherPhase.U()) - ); - } - } - } - return eqnsPtr; } -template -Foam::PtrList -Foam::MomentumTransferPhaseSystem::Vmfs() const -{ - PtrList Vmfs(this->phaseModels_.size()); - - // Add the implicit part of the virtual mass force - forAllConstIter(VmTable, Vms_, VmIter) - { - const volScalarField& Vm(*VmIter()); - const phaseInterface interface(*this, VmIter.key()); - - forAllConstIter(phaseInterface, interface, iter) - { - const phaseModel& phase = iter(); - const phaseModel& otherPhase = iter.otherPhase(); - - const volScalarField VmPhase - ( - (otherPhase/max(otherPhase, otherPhase.residualAlpha())) - *Vm - ); - - addField(phase, "Vmf", byDt(fvc::interpolate(VmPhase)), Vmfs); - } - } - - return Vmfs; -} - - template Foam::PtrList Foam::MomentumTransferPhaseSystem::Fs() const @@ -633,44 +424,6 @@ Foam::MomentumTransferPhaseSystem::Ffs() const { PtrList Ffs(this->phaseModels_.size()); - // Add the explicit part of the virtual mass force - forAllConstIter(VmTable, Vms_, VmIter) - { - const volScalarField& Vm(*VmIter()); - const phaseInterface interface(*this, VmIter.key()); - - forAllConstIter(phaseInterface, interface, iter) - { - const phaseModel& phase = iter(); - const phaseModel& otherPhase = iter.otherPhase(); - - const volScalarField VmPhase - ( - (otherPhase/max(otherPhase, otherPhase.residualAlpha())) - *Vm - ); - - addField - ( - phase, - "Ff", - -fvc::interpolate(VmPhase) - *( - byDt - ( - fvc::absolute - ( - this->MRF().absolute(iter().phi()().oldTime()), - iter().U() - ) - ) - + otherPhase.DUDtf() - ), - Ffs - ); - } - } - // Add the lift force forAllConstIter ( @@ -871,6 +624,7 @@ template void Foam::MomentumTransferPhaseSystem::invADs ( const PtrList& As, + PtrList& HVms, PtrList>& invADs, PtrList>& invADfs ) const @@ -962,6 +716,94 @@ void Foam::MomentumTransferPhaseSystem::invADs } } + // Cache the phase acceleration As and Hs + PtrList ADUDts(movingPhases.size()); + PtrList HDUDts(movingPhases.size()); + + forAllConstIter + ( + virtualMassModelTable, + virtualMassModels_, + VmIter + ) + { + const phaseInterface& interface = VmIter()->interface(); + + forAllConstIter(phaseInterface, interface, iter) + { + const phaseModel& phase = iter(); + const label i = movingPhases[phase.index()]; + + if (i != -1 && !ADUDts.set(i)) + { + const fvVectorMatrix DUDt(phase.DUDt()); + ADUDts.set(i, DUDt.A()); + HDUDts.set(i, DUDt.H()); + } + } + } + + // Add the virtual mass contributions + forAllConstIter + ( + virtualMassModelTable, + virtualMassModels_, + VmIter + ) + { + const phaseInterface& interface = VmIter()->interface(); + const volScalarField Vm(VmIter()->K()); + + forAllConstIter(phaseInterface, interface, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + const label i = movingPhases[phase.index()]; + + if (i != -1) + { + const volScalarField VmPhase + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*Vm + ); + + { + const volScalarField AVm(VmPhase*ADUDts[i]); + + invADs[i][i] += AVm; + invADfs[i][i] += fvc::interpolate(AVm); + + addField + ( + i, + "HVm", + VmPhase*HDUDts[i], + HVms + ); + } + + const label j = movingPhases[otherPhase.index()]; + + if (j != -1) + { + const volScalarField AVm(VmPhase*ADUDts[j]); + + invADs[i][j] -= AVm; + invADfs[i][j] -= fvc::interpolate(AVm); + + addField + ( + i, + "HVm", + -VmPhase*HDUDts[j], + HVms + ); + } + } + } + } + MomentumTransferPhaseSystem::invADs(invADs); MomentumTransferPhaseSystem::invADs(invADfs); } @@ -971,17 +813,18 @@ template Foam::PtrList> Foam::MomentumTransferPhaseSystem::invADfs ( - const PtrList& As + const PtrList& Afs, + PtrList& HVmfs ) const { - const label n = As.size(); + const label n = Afs.size(); PtrList> invADfs(n); forAll(invADfs, i) { invADfs.set(i, new PtrList(n)); - invADfs[i].set(i, As[i].clone()); + invADfs[i].set(i, Afs[i].clone()); } labelList movingPhases(this->phases().size(), -1); @@ -1040,13 +883,119 @@ Foam::MomentumTransferPhaseSystem::invADfs ( "0", this->mesh(), - dimensionedScalar(As[0].dimensions(), 0) + dimensionedScalar(Afs[0].dimensions(), 0) ) ); } } } + // Cache the phase acceleration Afs and Hs + PtrList AUgradUs(movingPhases.size()); + PtrList HUgradUs(movingPhases.size()); + + forAllConstIter + ( + virtualMassModelTable, + virtualMassModels_, + VmIter + ) + { + const phaseInterface& interface = VmIter()->interface(); + + forAllConstIter(phaseInterface, interface, iter) + { + const phaseModel& phase = iter(); + const label i = movingPhases[phase.index()]; + + if (i != -1 && !AUgradUs.set(i)) + { + const fvVectorMatrix UgradU(phase.UgradU()); + AUgradUs.set(i, UgradU.A()); + HUgradUs.set(i, UgradU.H()); + } + } + } + + // Add the virtual mass contributions + forAllConstIter + ( + virtualMassModelTable, + virtualMassModels_, + VmIter + ) + { + const phaseInterface& interface = VmIter()->interface(); + const volScalarField Vm(VmIter()->K()); + + forAllConstIter(phaseInterface, interface, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + const label i = movingPhases[phase.index()]; + + if (i != -1) + { + const volScalarField VmPhase + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *Vm + ); + + const surfaceScalarField VmPhasef(fvc::interpolate(VmPhase)); + + invADfs[i][i] += + byDt(VmPhasef) + fvc::interpolate(VmPhase*AUgradUs[i]); + + addField + ( + i, + "HVmf", + VmPhasef + *byDt + ( + fvc::absolute + ( + this->MRF().absolute(phase.phi()().oldTime()), + phase.U() + ) + ) + + fvc::flux(VmPhase*HUgradUs[i]), + HVmfs + ); + + const label j = movingPhases[otherPhase.index()]; + + if (j != -1) + { + invADfs[i][j] -= + byDt(VmPhasef) + fvc::interpolate(VmPhase*AUgradUs[j]); + + addField + ( + i, + "HVmf", + -VmPhasef + *byDt + ( + fvc::absolute + ( + this->MRF().absolute + ( + otherPhase.phi()().oldTime() + ), + otherPhase.U() + ) + ) + - fvc::flux(VmPhase*HUgradUs[j]), + HVmfs + ); + } + } + } + } + invADs(invADfs); return invADfs; @@ -1194,12 +1143,6 @@ Foam::MomentumTransferPhaseSystem::ddtCorrs() const const phaseModel& phase = this->movingPhases()[movingPhasei]; const label phasei = phase.index(); - VmDdtCoeffs.set - ( - phasei, - fvm::ddt(phase.U()())().A() - ); - VmDdtCorrs.set ( phasei, @@ -1212,10 +1155,15 @@ Foam::MomentumTransferPhaseSystem::ddtCorrs() const ); } - forAllConstIter(VmTable, Vms_, VmIter) + forAllConstIter + ( + virtualMassModelTable, + virtualMassModels_, + VmIter + ) { - const volScalarField& Vm(*VmIter()); - const phaseInterface interface(*this, VmIter.key()); + const phaseInterface& interface = VmIter()->interface(); + const volScalarField Vm(VmIter()->K()); forAllConstIter(phaseInterface, interface, iter) { @@ -1235,19 +1183,7 @@ Foam::MomentumTransferPhaseSystem::ddtCorrs() const phase, "ddtCorr", fvc::interpolate(VmPhase) - *( - VmDdtCorrs[phasei] - + ( - fvc::interpolate(VmDdtCoeffs[otherPhasei]) - *( - otherPhase.Uf().valid() - ? (this->mesh_.Sf() & otherPhase.Uf()())() - : otherPhase.phi()() - ) - - fvc::flux(VmDdtCoeffs[otherPhasei]*otherPhase.U()) - ) - - VmDdtCorrs[otherPhase.index()] - ), + *(VmDdtCorrs[phasei] - VmDdtCorrs[otherPhasei]), ddtCorrs ); } diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H index 8c35852810..9333d011f3 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H @@ -79,13 +79,6 @@ class MomentumTransferPhaseSystem phaseInterfaceKey::hash > KdfTable; - typedef HashPtrTable - < - volScalarField, - phaseInterfaceKey, - phaseInterfaceKey::hash - > VmTable; - typedef HashTable < autoPtr, @@ -130,9 +123,6 @@ class MomentumTransferPhaseSystem //- Face drag coefficients KdfTable Kdfs_; - //- Virtual mass coefficients - VmTable Vms_; - // Sub Models @@ -202,10 +192,6 @@ public: //- As momentumTransfer, but for the face-based algorithm virtual autoPtr momentumTransferf(); - //- Return implicit force coefficients on the faces, for the face-based - // algorithm. - virtual PtrList Vmfs() const; - //- Return the explicit force fluxes for the cell-based algorithm, that // do not depend on phase mass/volume fluxes, and can therefore be // evaluated outside the corrector loop. This includes things like @@ -219,6 +205,7 @@ public: virtual void invADs ( const PtrList& As, + PtrList& HVms, PtrList>& invADs, PtrList>& invADfs ) const; @@ -226,7 +213,8 @@ public: //- Return the inverse of (Afs + Dfs) virtual PtrList> invADfs ( - const PtrList& Afs + const PtrList& Afs, + PtrList& HVmfs ) const; //- Returns true if the phase pressure is treated implicitly diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C index a75563de45..049c339224 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C @@ -162,8 +162,6 @@ Foam::MovingPhaseModel::MovingPhaseModel dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), 0) ), Uf_(nullptr), - DUDt_(nullptr), - DUDtf_(nullptr), divU_(nullptr), momentumTransport_ ( @@ -252,18 +250,6 @@ void Foam::MovingPhaseModel::correctKinematics() { BasePhaseModel::correctKinematics(); - if (DUDt_.valid()) - { - DUDt_.clear(); - DUDt(); - } - - if (DUDtf_.valid()) - { - DUDtf_.clear(); - DUDtf(); - } - if (K_.valid()) { K_.ref() = 0.5*magSqr(this->U()); @@ -509,40 +495,23 @@ Foam::MovingPhaseModel::alphaRhoPhiRef() const template -Foam::tmp -Foam::MovingPhaseModel::DUDt() const +Foam::tmp +Foam::MovingPhaseModel::UgradU() const { - if (!DUDt_.valid()) - { - const tmp taphi(fvc::absolute(phi_, U_)); - const surfaceScalarField& aphi(taphi()); - DUDt_ = - new volVectorField - ( - IOobject::groupName("DUDt", this->name()), - fvc::ddt(U_) + fvc::div(aphi, U_) - fvc::div(aphi)*U_ - ); - } + const tmp taphi(fvc::absolute(phi_, U_)); + const surfaceScalarField& aphi(taphi()); - return tmp(DUDt_()); + return + fvm::div(aphi, U_) - fvm::Sp(fvc::div(aphi), U_) + + this->fluid().MRF().DDt(U_); } template -Foam::tmp -Foam::MovingPhaseModel::DUDtf() const +Foam::tmp +Foam::MovingPhaseModel::DUDt() const { - if (!DUDtf_.valid()) - { - DUDtf_ = - new surfaceScalarField - ( - IOobject::groupName("DUDtf", this->name()), - byDt(phi_ - phi_.oldTime()) - ); - } - - return tmp(DUDtf_()); + return fvm::ddt(U_) + UgradU(); } diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H index 29431aceb6..b1e2463087 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H @@ -112,12 +112,6 @@ protected: //- Face velocity field autoPtr Uf_; - //- Lagrangian acceleration field (needed for virtual-mass) - mutable tmp DUDt_; - - //- Lagrangian acceleration field on the faces (needed for virtual-mass) - mutable tmp DUDtf_; - //- Dilatation rate autoPtr divU_; @@ -253,11 +247,11 @@ public: //- Access the mass flux of the phase virtual const surfaceScalarField& alphaRhoPhiRef() const; - //- Return the substantive acceleration - virtual tmp DUDt() const; + //- Return the velocity transport matrix + virtual tmp UgradU() const; - //- Return the substantive acceleration on the faces - virtual tmp DUDtf() const; + //- Return the substantive acceleration matrix + virtual tmp DUDt() const; //- Return the continuity error virtual tmp continuityError() const; diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C index f5f21a7759..fddd2054aa 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C @@ -267,28 +267,26 @@ Foam::StationaryPhaseModel::alphaRhoPhiRef() const template -Foam::tmp -Foam::StationaryPhaseModel::DUDt() const +Foam::tmp +Foam::StationaryPhaseModel::UgradU() const { - return volVectorField::New - ( - IOobject::groupName("DUDt", this->name()), - this->mesh(), - dimensionedVector(dimVelocity/dimTime, Zero) - ); + FatalErrorInFunction + << "Cannot calculate DUDt of a stationary phase" + << exit(FatalError); + + return tmp(nullptr); } template -Foam::tmp -Foam::StationaryPhaseModel::DUDtf() const +Foam::tmp +Foam::StationaryPhaseModel::DUDt() const { - return surfaceScalarField::New - ( - IOobject::groupName("DUDtf", this->name()), - this->mesh(), - dimensionedScalar(dimVolume/sqr(dimTime), 0) - ); + FatalErrorInFunction + << "Cannot calculate DUDt of a stationary phase" + << exit(FatalError); + + return tmp(nullptr); } diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H index b23aa2ce00..3d4b7ded34 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H @@ -136,11 +136,11 @@ public: //- Access the mass flux of the phase virtual const surfaceScalarField& alphaRhoPhiRef() const; - //- Return the substantive acceleration - virtual tmp DUDt() const; + //- Return the velocity transport matrix + virtual tmp UgradU() const; - //- Return the substantive acceleration on the faces - virtual tmp DUDtf() const; + //- Return the substantive acceleration matrix + virtual tmp DUDt() const; //- Return the continuity error virtual tmp continuityError() const; diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H index f42991ab34..4a355e9f95 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H @@ -359,11 +359,11 @@ public: //- Access the mass flux of the phase virtual const surfaceScalarField& alphaRhoPhiRef() const = 0; - //- Return the substantive acceleration - virtual tmp DUDt() const = 0; + //- Return the velocity transport matrix + virtual tmp UgradU() const = 0; - //- Return the substantive acceleration on the faces - virtual tmp DUDtf() const = 0; + //- Return the substantive acceleration matrix + virtual tmp DUDt() const = 0; //- Return the continuity error virtual tmp continuityError() const = 0; diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H index 4a87aaa0a4..c26213db0b 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H @@ -521,10 +521,6 @@ public: // algorithm virtual autoPtr momentumTransferf() = 0; - //- Return the implicit force coefficients for the face-based - // algorithm - virtual PtrList Vmfs() const = 0; - //- Return the force fluxes for the cell-based algorithm virtual PtrList Fs() const = 0; @@ -535,6 +531,7 @@ public: virtual void invADs ( const PtrList& As, + PtrList& HVms, PtrList>& invADs, PtrList>& invADfs ) const = 0; @@ -542,7 +539,8 @@ public: //- Return the inverse of (Afs + Dfs) virtual PtrList> invADfs ( - const PtrList& Afs + const PtrList& Afs, + PtrList& HVmfs ) const = 0; //- Returns true if the phase pressure is treated implicitly