diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C index 11f48b6204..8b84be088d 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C @@ -69,8 +69,8 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() } PtrList HVms(movingPhases.size()); - PtrList> invADs; - PtrList> invADfs; + PtrList> invADVs; + PtrList> invADVfs; { PtrList As(movingPhases.size()); PtrList Afs(movingPhases.size()); @@ -110,7 +110,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() } } - fluid.invADs(As, HVms, invADs, invADfs); + fluid.invADVs(As, HVms, invADVs, invADVfs); } volScalarField rho("rho", fluid.rho()); @@ -157,8 +157,8 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ); } - alphaByADfs = invADfs & lalphafs; - FgByADfs = invADfs & Fgfs; + alphaByADfs = invADVfs & lalphafs; + FgByADfs = invADVfs & Fgfs; } @@ -196,7 +196,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() } } - HbyADs = invADs & Hs; + HbyADs = invADVs & Hs; } forAll(movingPhases, movingPhasei) @@ -270,8 +270,8 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ); } - HbyADs = invADs & Hs; - phiHbyADs = invADfs & phiHs; + HbyADs = invADVs & Hs; + phiHbyADs = invADVfs & phiHs; } forAll(movingPhases, movingPhasei) @@ -428,12 +428,12 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() PtrList dragCorrByADs ( - invADs & dragCorrs + invADVs & dragCorrs ); PtrList dragCorrByADfs ( - invADfs & dragCorrfs + invADVfs & dragCorrfs ); forAll(movingPhases, movingPhasei) diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C index 05e5b07f35..68c4ca1634 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C @@ -69,7 +69,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() } PtrList HVmfs(movingPhases.size()); - PtrList> invADfs; + PtrList> invADVfs; { PtrList Afs(movingPhases.size()); @@ -112,7 +112,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() ); } - invADfs = fluid.invADfs(Afs, HVmfs); + invADVfs = fluid.invADVfs(Afs, HVmfs); } volScalarField rho("rho", fluid.rho()); @@ -157,8 +157,8 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() ); } - alphaByADfs = invADfs & lalphafs; - FgByADfs = invADfs & Fgfs; + alphaByADfs = invADVfs & lalphafs; + FgByADfs = invADVfs & Fgfs; } @@ -207,7 +207,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() } } - phiHbyADs = invADfs & phiHs; + phiHbyADs = invADVfs & phiHs; } forAll(movingPhases, movingPhasei) diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index b79499374b..f982c16fdc 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -158,72 +158,6 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() ); } - // Initialise Kds table if required - if (!Kds_.size()) - { - forAllConstIter - ( - dragModelTable, - dragModels_, - dragModelIter - ) - { - const phaseInterface& interface = dragModelIter()->interface(); - - Kds_.insert - ( - dragModelIter.key(), - new volScalarField - ( - IOobject - ( - IOobject::groupName("Kd", interface.name()), - this->mesh().time().name(), - this->mesh() - ), - this->mesh(), - dimensionedScalar(dragModel::dimK, 0) - ) - ); - } - } - - // Update the drag coefficients - forAllConstIter - ( - dragModelTable, - dragModels_, - dragModelIter - ) - { - *Kds_[dragModelIter.key()] = dragModelIter()->K(); - - // Zero-gradient the drag coefficient to boundaries with fixed velocity - const phaseInterface& interface = dragModelIter()->interface(); - volScalarField& K = *Kds_[dragModelIter.key()]; - - forAll(K.boundaryField(), patchi) - { - if - ( - ( - !interface.phase1().stationary() - && interface.phase1().U()() - .boundaryField()[patchi].fixesValue() - ) - && ( - !interface.phase2().stationary() - && interface.phase2().U()() - .boundaryField()[patchi].fixesValue() - ) - ) - { - K.boundaryFieldRef()[patchi] = - K.boundaryField()[patchi].patchInternalField(); - } - } - } - return eqnsPtr; } @@ -251,47 +185,6 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() ); } - // Initialise Kdfs table if required - if (!Kdfs_.size()) - { - forAllConstIter - ( - dragModelTable, - dragModels_, - dragModelIter - ) - { - const phaseInterface& interface = dragModelIter()->interface(); - - Kdfs_.insert - ( - dragModelIter.key(), - new surfaceScalarField - ( - IOobject - ( - IOobject::groupName("Kdf", interface.name()), - this->mesh().time().name(), - this->mesh() - ), - this->mesh(), - dimensionedScalar(dragModel::dimK, 0) - ) - ); - } - } - - // Update the drag coefficients - forAllConstIter - ( - dragModelTable, - dragModels_, - dragModelIter - ) - { - *Kdfs_[dragModelIter.key()] = dragModelIter()->Kf(); - } - return eqnsPtr; } @@ -539,24 +432,24 @@ Foam::MomentumTransferPhaseSystem::Ffs() const template -void Foam::MomentumTransferPhaseSystem::invADs +void Foam::MomentumTransferPhaseSystem::invADVs ( - List>& ADs + List>& ADVs ) const { - const label n = ADs.size(); + const label n = ADVs.size(); scalarSquareMatrix AD(n); scalarField source(n); labelList pivotIndices(n); - forAll(ADs[0][0], ci) + forAll(ADVs[0][0], ci) { for (label i=0; i::invADs for (label i=0; i::invADs template template class PatchField, class GeoMesh> -void Foam::MomentumTransferPhaseSystem::invADs +void Foam::MomentumTransferPhaseSystem::invADVs ( - PtrList>>& ADs + PtrList>>& ADVs ) const { - const label n = ADs.size(); + const label n = ADVs.size(); List> ADps(n); @@ -597,50 +490,50 @@ void Foam::MomentumTransferPhaseSystem::invADs for (label j=0; j -void Foam::MomentumTransferPhaseSystem::invADs +void Foam::MomentumTransferPhaseSystem::invADVs ( const PtrList& As, PtrList& HVms, - PtrList>& invADs, - PtrList>& invADfs + PtrList>& invADVs, + PtrList>& invADVfs ) const { const label n = As.size(); - invADs.setSize(n); - invADfs.setSize(n); + invADVs.setSize(n); + invADVfs.setSize(n); - forAll(invADs, i) + forAll(invADVs, i) { - invADs.set(i, new PtrList(n)); - invADs[i].set(i, As[i].clone()); + invADVs.set(i, new PtrList(n)); + invADVs[i].set(i, As[i].clone()); - invADfs.set(i, new PtrList(n)); - invADfs[i].set(i, fvc::interpolate(As[i])); + invADVfs.set(i, new PtrList(n)); + invADVfs[i].set(i, fvc::interpolate(As[i])); } labelList movingPhases(this->phases().size(), -1); @@ -650,10 +543,43 @@ void Foam::MomentumTransferPhaseSystem::invADs movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei; } - forAllConstIter(KdTable, Kds_, KdIter) + Kds_.clear(); + + // Update the drag coefficients + forAllConstIter + ( + dragModelTable, + dragModels_, + dragModelIter + ) { - const volScalarField& K(*KdIter()); - const phaseInterface interface(*this, KdIter.key()); + const phaseInterface& interface = dragModelIter()->interface(); + + tmp tKd(dragModelIter()->K()); + volScalarField& Kd = tKd.ref(); + Kds_.insert(dragModelIter.key(), tKd.ptr()); + + // Zero-gradient the drag coefficient to boundaries with fixed velocity + forAll(Kd.boundaryField(), patchi) + { + if + ( + ( + !interface.phase1().stationary() + && interface.phase1().U()() + .boundaryField()[patchi].fixesValue() + ) + && ( + !interface.phase2().stationary() + && interface.phase2().U()() + .boundaryField()[patchi].fixesValue() + ) + ) + { + Kd.boundaryFieldRef()[patchi] = + Kd.boundaryField()[patchi].patchInternalField(); + } + } forAllConstIter(phaseInterface, interface, iter) { @@ -664,34 +590,41 @@ void Foam::MomentumTransferPhaseSystem::invADs if (i != -1) { - const volScalarField Kij + const volScalarField Kdij ( - (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*K + (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*Kd ); - const surfaceScalarField Kijf(fvc::interpolate(Kij)); + const surfaceScalarField Kdijf(fvc::interpolate(Kdij)); - invADs[i][i] += Kij; - invADfs[i][i] += Kijf; + invADVs[i][i] += Kdij; + invADVfs[i][i] += Kdijf; const label j = movingPhases[otherPhase.index()]; if (j != -1) { - invADs[i].set(j, -Kij); - invADfs[i].set(j, -Kijf); + invADVs[i].set(j, -Kdij); + invADVfs[i].set(j, -Kdijf); } } } } + // Clear the Kds_ if they are not needed for the optional dragCorrection + const pimpleNoLoopControl& pimple = this->pimple(); + if (!pimple.dict().lookupOrDefault("dragCorrection", false)) + { + Kds_.clear(); + } + for (label i=0; i::invADs ) ); - invADfs[i].set + invADVfs[i].set ( j, surfaceScalarField::New @@ -771,8 +704,8 @@ void Foam::MomentumTransferPhaseSystem::invADs { const volScalarField AVm(VmPhase*ADUDts[i]); - invADs[i][i] += AVm; - invADfs[i][i] += fvc::interpolate(AVm); + invADVs[i][i] += AVm; + invADVfs[i][i] += fvc::interpolate(AVm); addField ( @@ -789,8 +722,8 @@ void Foam::MomentumTransferPhaseSystem::invADs { const volScalarField AVm(VmPhase*ADUDts[j]); - invADs[i][j] -= AVm; - invADfs[i][j] -= fvc::interpolate(AVm); + invADVs[i][j] -= AVm; + invADVfs[i][j] -= fvc::interpolate(AVm); addField ( @@ -804,14 +737,14 @@ void Foam::MomentumTransferPhaseSystem::invADs } } - MomentumTransferPhaseSystem::invADs(invADs); - MomentumTransferPhaseSystem::invADs(invADfs); + MomentumTransferPhaseSystem::invADVs(invADVs); + MomentumTransferPhaseSystem::invADVs(invADVfs); } template Foam::PtrList> -Foam::MomentumTransferPhaseSystem::invADfs +Foam::MomentumTransferPhaseSystem::invADVfs ( const PtrList& Afs, PtrList& HVmfs @@ -819,12 +752,12 @@ Foam::MomentumTransferPhaseSystem::invADfs { const label n = Afs.size(); - PtrList> invADfs(n); + PtrList> invADVfs(n); - forAll(invADfs, i) + forAll(invADVfs, i) { - invADfs.set(i, new PtrList(n)); - invADfs[i].set(i, Afs[i].clone()); + invADVfs.set(i, new PtrList(n)); + invADVfs[i].set(i, Afs[i].clone()); } labelList movingPhases(this->phases().size(), -1); @@ -834,10 +767,15 @@ Foam::MomentumTransferPhaseSystem::invADfs movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei; } - forAllConstIter(KdfTable, Kdfs_, KdfIter) + forAllConstIter + ( + dragModelTable, + dragModels_, + dragModelIter + ) { - const surfaceScalarField& Kf(*KdfIter()); - const phaseInterface interface(*this, KdfIter.key()); + const phaseInterface& interface = dragModelIter()->interface(); + const surfaceScalarField Kdf(dragModelIter()->Kf()); forAllConstIter(phaseInterface, interface, iter) { @@ -853,18 +791,18 @@ Foam::MomentumTransferPhaseSystem::invADfs fvc::interpolate(otherPhase) ); - const surfaceScalarField Kfij + const surfaceScalarField Kdfij ( - (alphaf/max(alphaf, otherPhase.residualAlpha()))*Kf + (alphaf/max(alphaf, otherPhase.residualAlpha()))*Kdf ); - invADfs[i][i] += Kfij; + invADVfs[i][i] += Kdfij; const label j = movingPhases[otherPhase.index()]; if (j != -1) { - invADfs[i].set(j, -Kfij); + invADVfs[i].set(j, -Kdfij); } } } @@ -874,9 +812,9 @@ Foam::MomentumTransferPhaseSystem::invADfs { for (label j=0; j::invADfs const surfaceScalarField VmPhasef(fvc::interpolate(VmPhase)); - invADfs[i][i] += + invADVfs[i][i] += byDt(VmPhasef) + fvc::interpolate(VmPhase*AUgradUs[i]); addField @@ -969,7 +907,7 @@ Foam::MomentumTransferPhaseSystem::invADfs if (j != -1) { - invADfs[i][j] -= + invADVfs[i][j] -= byDt(VmPhasef) + fvc::interpolate(VmPhase*AUgradUs[j]); addField @@ -996,9 +934,9 @@ Foam::MomentumTransferPhaseSystem::invADfs } } - invADs(invADfs); + invADVs(invADVfs); - return invADfs; + return invADVfs; } diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H index 9333d011f3..bac1fffe18 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H @@ -72,13 +72,6 @@ class MomentumTransferPhaseSystem phaseInterfaceKey::hash > KdTable; - typedef HashPtrTable - < - surfaceScalarField, - phaseInterfaceKey, - phaseInterfaceKey::hash - > KdfTable; - typedef HashTable < autoPtr, @@ -117,11 +110,8 @@ class MomentumTransferPhaseSystem // Private Data - //- Drag coefficients - KdTable Kds_; - - //- Face drag coefficients - KdfTable Kdfs_; + //- Cached drag coefficients for dragCorrs + mutable KdTable Kds_; // Sub Models @@ -159,14 +149,14 @@ protected: const tmp& field ) const; - //- Invert the ADs matrix inplace - void invADs(List>& ADs) const; + //- Invert the ADVs matrix inplace + void invADVs(List>& ADVs) const; - //- Invert the ADs matrix inplace + //- Invert the ADVs matrix inplace template class PatchField, class GeoMesh> - void invADs + void invADVs ( - PtrList>>& ADs + PtrList>>& ADVs ) const; @@ -201,17 +191,19 @@ public: //- As Fs, but for the face-based algorithm virtual PtrList Ffs() const; - //- Return the inverse of (As + Ds) - virtual void invADs + //- Return the inverse of the central + drag + virtual mass + // coefficient matrix + virtual void invADVs ( const PtrList& As, PtrList& HVms, - PtrList>& invADs, - PtrList>& invADfs + PtrList>& invADVs, + PtrList>& invADVfs ) const; - //- Return the inverse of (Afs + Dfs) - virtual PtrList> invADfs + //- Return the inverse of the central + drag + virtual mass + // coefficient matrix + virtual PtrList> invADVfs ( const PtrList& Afs, PtrList& HVmfs diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H index c26213db0b..31953f14c3 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H @@ -527,17 +527,19 @@ public: //- Return the force fluxes for the face-based algorithm virtual PtrList Ffs() const = 0; - //- Return the inverse of (As + Ds) - virtual void invADs + //- Return the inverse of the central + drag + virtual mass + // coefficient matrix + virtual void invADVs ( const PtrList& As, PtrList& HVms, - PtrList>& invADs, - PtrList>& invADfs + PtrList>& invADVs, + PtrList>& invADVfs ) const = 0; - //- Return the inverse of (Afs + Dfs) - virtual PtrList> invADfs + //- Return the inverse of the face central + drag + virtual mass + // coefficient matrix + virtual PtrList> invADVfs ( const PtrList& Afs, PtrList& HVmfs