diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean b/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean index 5053a12956..b8fe606783 100755 --- a/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean +++ b/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean @@ -7,7 +7,6 @@ wclean libso interfacialCompositionModels wclean libso multiphaseCompressibleMomentumTransportModels wclean libso multiphaseThermophysicalTransportModels multiphaseEulerFoam/Allwclean -reactingTwoPhaseEulerFoam/Allwclean wclean libso functionObjects #------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake b/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake index 5a543f5383..db3b015506 100755 --- a/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake +++ b/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake @@ -10,7 +10,6 @@ wmake $targetType interfacialCompositionModels wmake $targetType multiphaseCompressibleMomentumTransportModels wmake $targetType multiphaseThermophysicalTransportModels multiphaseEulerFoam/Allwmake $targetType $* -reactingTwoPhaseEulerFoam/Allwmake $targetType $* wmake $targetType functionObjects #------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Allwclean b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Allwclean deleted file mode 100755 index 11f2624bb1..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Allwclean +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/sh -cd ${0%/*} || exit 1 # Run from this directory - -wclean libso twoPhaseSystem -wclean - -#------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Allwmake b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Allwmake deleted file mode 100755 index f5e9b28053..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Allwmake +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/sh -cd ${0%/*} || exit 1 # Run from this directory - -# Parse arguments for library compilation -. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments - -wmake $targetType twoPhaseSystem -wmake $targetType - -#------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/CourantNos.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/CourantNos.H deleted file mode 100644 index 13c45ad51c..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/CourantNos.H +++ /dev/null @@ -1,12 +0,0 @@ -#include "CourantNo.H" - -{ - scalar UrCoNum = 0.5*gMax - ( - fvc::surfaceSum(mag(phi1 - phi2))().primitiveField()/mesh.V().field() - )*runTime.deltaTValue(); - - Info<< "Max Ur Courant Number = " << UrCoNum << endl; - - CoNum = max(CoNum, UrCoNum); -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H deleted file mode 100644 index a132accc6a..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H +++ /dev/null @@ -1,47 +0,0 @@ -for (int Ecorr=0; Ecorr - heatTransferPtr(fluid.heatTransfer()); - - phaseSystem::heatTransferTable& - heatTransfer = heatTransferPtr(); - - if (!phase1.isothermal()) - { - fvScalarMatrix E1Eqn - ( - phase1.heEqn() - == - *heatTransfer[phase1.name()] - + alpha1*rho1*(U1&g) - + fvOptions(alpha1, rho1, phase1.thermoRef().he()) - ); - - E1Eqn.relax(); - fvOptions.constrain(E1Eqn); - E1Eqn.solve(); - fvOptions.correct(thermo1.he()); - } - - if (!phase2.isothermal()) - { - fvScalarMatrix E2Eqn - ( - phase2.heEqn() - == - *heatTransfer[phase2.name()] - + alpha2*rho2*(U2&g) - + fvOptions(alpha2, rho2, phase2.thermoRef().he()) - ); - - E2Eqn.relax(); - fvOptions.constrain(E2Eqn); - E2Eqn.solve(); - fvOptions.correct(thermo2.he()); - } - - fluid.correctThermo(); - fluid.correctContinuityError(); -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Make/files deleted file mode 100644 index 66b068b0c3..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -reactingTwoPhaseEulerFoam.C - -EXE = $(FOAM_APPBIN)/reactingTwoPhaseEulerFoam diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Make/options deleted file mode 100644 index cd2b06ffa5..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/Make/options +++ /dev/null @@ -1,26 +0,0 @@ -EXE_INC = \ - -I../include \ - -ItwoPhaseSystem/lnInclude \ - -I../phaseSystems/lnInclude \ - -I../interfacialModels/lnInclude \ - -I../interfacialCompositionModels/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \ - -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \ - -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude - -EXE_LIBS = \ - -lphaseSystem \ - -ltwoPhaseSystem \ - -leulerianInterfacialModels \ - -leulerianInterfacialCompositionModels \ - -lmultiphaseMomentumTransportModels \ - -lmultiphaseThermophysicalTransportModels \ - -lthermophysicalTransportModels \ - -lfiniteVolume \ - -lfvOptions \ - -lmeshTools \ - -lsampling diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H deleted file mode 100644 index d36742bf13..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/YEqns.H +++ /dev/null @@ -1,49 +0,0 @@ -{ - autoPtr - specieTransferPtr(fluid.specieTransfer()); - - phaseSystem::specieTransferTable& - specieTransfer(specieTransferPtr()); - - fluid.correctReactions(); - - if (!phase1.pure()) - { - UPtrList& Y1 = phase1.YActiveRef(); - - forAll(Y1, i) - { - fvScalarMatrix Y1iEqn - ( - phase1.YiEqn(Y1[i]) - == - *specieTransfer[Y1[i].name()] - + fvOptions(alpha1, rho1, Y1[i]) - ); - - Y1iEqn.relax(); - Y1iEqn.solve("Yi"); - } - } - - if (!phase2.pure()) - { - UPtrList& Y2 = phase2.YActiveRef(); - - forAll(Y2, i) - { - fvScalarMatrix Y2iEqn - ( - phase2.YiEqn(Y2[i]) - == - *specieTransfer[Y2[i].name()] - + fvOptions(alpha2, rho2, Y2[i]) - ); - - Y2iEqn.relax(); - Y2iEqn.solve("Yi"); - } - } - - fluid.correctSpecies(); -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H deleted file mode 100644 index edad54fa18..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/createFieldRefs.H +++ /dev/null @@ -1,27 +0,0 @@ -phaseModel& phase1 = fluid.phase1(); -phaseModel& phase2 = fluid.phase2(); - -const volScalarField& alpha1 = phase1; -const volScalarField& alpha2 = phase2; - -volVectorField& U1 = phase1.URef(); -surfaceScalarField& phi1 = phase1.phiRef(); -const surfaceScalarField& alphaPhi1 = phase1.alphaPhi(); - -volVectorField& U2 = phase2.URef(); -surfaceScalarField& phi2 = phase2.phiRef(); -const surfaceScalarField& alphaPhi2 = phase2.alphaPhi(); - -surfaceScalarField& phi = fluid.phi(); - -rhoThermo& thermo1 = phase1.thermoRef(); -rhoThermo& thermo2 = phase2.thermoRef(); - -volScalarField& rho1 = thermo1.rho(); -const volScalarField& psi1 = thermo1.psi(); - -volScalarField& rho2 = thermo2.rho(); -const volScalarField& psi2 = thermo2.psi(); - -const IOMRFZoneList& MRF = fluid.MRF(); -fv::options& fvOptions = fluid.fvOptions(); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/createFields.H deleted file mode 100644 index 69ad38e6ce..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/createFields.H +++ /dev/null @@ -1,67 +0,0 @@ -#include "createRDeltaT.H" -#include "readGravitationalAcceleration.H" -#include "readhRef.H" - -Info<< "Creating phaseSystem\n" << endl; - -autoPtr fluidPtr -( - twoPhaseSystem::New(mesh) -); -twoPhaseSystem& fluid = fluidPtr(); - -dimensionedScalar pMin -( - "pMin", - dimPressure, - fluid -); - -#include "gh.H" - -volScalarField& p = fluid.phase1().thermoRef().p(); - -Info<< "Reading field p_rgh\n" << endl; -volScalarField p_rgh -( - IOobject - ( - "p_rgh", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -label pRefCell = 0; -scalar pRefValue = 0.0; -if (fluid.incompressible()) -{ - p = max(p_rgh + fluid.rho()*gh, pMin); - - if (p_rgh.needReference()) - { - setRefCell - ( - p, - p_rgh, - pimple.dict(), - pRefCell, - pRefValue - ); - - p += dimensionedScalar - ( - "p", - p.dimensions(), - pRefValue - getRefCellValue(p, pRefCell) - ); - p_rgh = p - fluid.rho()*gh; - } -} -mesh.setFluxRequired(p_rgh.name()); - -PtrList rAUs; -PtrList rAUfs; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H deleted file mode 100644 index 99bc95cfd6..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H +++ /dev/null @@ -1,38 +0,0 @@ -Info<< "Constructing momentum equations" << endl; - -fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime); -fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime); - -{ - autoPtr - momentumTransferPtr(fluid.momentumTransfer()); - - phaseSystem::momentumTransferTable& - momentumTransfer(momentumTransferPtr()); - - { - U1Eqn = - ( - phase1.UEqn() - == - *momentumTransfer[phase1.name()] - + fvOptions(alpha1, rho1, U1) - ); - U1Eqn.relax(); - fvOptions.constrain(U1Eqn); - fvOptions.correct(U1); - } - - { - U2Eqn = - ( - phase2.UEqn() - == - *momentumTransfer[phase2.name()] - + fvOptions(alpha2, rho2, U2) - ); - U2Eqn.relax(); - fvOptions.constrain(U2Eqn); - fvOptions.correct(U2); - } -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H deleted file mode 100644 index 0d00e8a16e..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H +++ /dev/null @@ -1,400 +0,0 @@ -const surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1)); -const surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1); - -rAUs.clear(); - -rAUs.append -( - new volScalarField - ( - IOobject::groupName("rAU", phase1.name()), - 1.0 - /( - U1Eqn.A() - + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1) - ) - ) -); -rAUs.append -( - new volScalarField - ( - IOobject::groupName("rAU", phase2.name()), - 1.0 - /( - U2Eqn.A() - + byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2) - ) - ) -); -const volScalarField& rAU1 = rAUs[0]; -const volScalarField& rAU2 = rAUs[1]; - -const surfaceScalarField alpharAUf1 -( - fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1) -); -const surfaceScalarField alpharAUf2 -( - fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2) -); - -// Drag coefficients -const volScalarField Kd(fluid.Kd()); -const volScalarField rAUKd1(rAU1*Kd); -const volScalarField rAUKd2(rAU2*Kd); -const surfaceScalarField phiKd1(fvc::interpolate(rAUKd1)); -const surfaceScalarField phiKd2(fvc::interpolate(rAUKd2)); - -// Explicit force fluxes -PtrList phiFs(fluid.phiFs(rAUs)); -const surfaceScalarField& phiF1 = phiFs[0]; -const surfaceScalarField& phiF2 = phiFs[1]; - -// --- Pressure corrector loop -while (pimple.correct()) -{ - volScalarField rho("rho", fluid.rho()); - - // Correct p_rgh for consistency with p and the updated densities - p_rgh = p - rho*gh; - - // Correct fixed-flux BCs to be consistent with the velocity BCs - MRF.correctBoundaryFlux(U1, phi1); - MRF.correctBoundaryFlux(U2, phi2); - - // Combined buoyancy and force fluxes - const surfaceScalarField ghSnGradRho - ( - "ghSnGradRho", - ghf*fvc::snGrad(rho)*mesh.magSf() - ); - - const surfaceScalarField phigF1 - ( - alpharAUf1 - *( - ghSnGradRho - - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf()) - ) - + phiF1 - ); - - const surfaceScalarField phigF2 - ( - alpharAUf2 - *( - ghSnGradRho - - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf()) - ) - + phiF2 - ); - - // Predicted velocities - volVectorField HbyA1 - ( - IOobject::groupName("HbyA", phase1.name()), - U1 - ); - HbyA1 = - rAU1 - *( - U1Eqn.H() - + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1) - *U1.oldTime() - ); - - volVectorField HbyA2 - ( - IOobject::groupName("HbyA", phase2.name()), - U2 - ); - HbyA2 = - rAU2 - *( - U2Eqn.H() - + byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2) - *U2.oldTime() - ); - - // Correction force fluxes - PtrList ddtCorrByAs(fluid.ddtCorrByAs(rAUs)); - - // Predicted fluxes - const surfaceScalarField phiHbyA1 - ( - IOobject::groupName("phiHbyA", phase1.name()), - fvc::flux(HbyA1) - phigF1 - ddtCorrByAs[0] - ); - - const surfaceScalarField phiHbyA2 - ( - IOobject::groupName("phiHbyA", phase2.name()), - fvc::flux(HbyA2) - phigF2 - ddtCorrByAs[1] - ); - - ddtCorrByAs.clear(); - - // Drag fluxes - PtrList phiKdPhis(fluid.phiKdPhis(rAUs)); - const surfaceScalarField& phiKdPhi1 = phiKdPhis[0]; - const surfaceScalarField& phiKdPhi2 = phiKdPhis[1]; - - // Total predicted flux - surfaceScalarField phiHbyA - ( - "phiHbyA", - alphaf1*(phiHbyA1 - phiKdPhi1) + alphaf2*(phiHbyA2 - phiKdPhi2) - ); - - MRF.makeRelative(phiHbyA); - - phiKdPhis.clear(); - - // Construct pressure "diffusivity" - const surfaceScalarField rAUf - ( - "rAUf", - mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2) - ); - - // Update the fixedFluxPressure BCs to ensure flux consistency - setSnGrad - ( - p_rgh.boundaryFieldRef(), - ( - phiHbyA.boundaryField() - - ( - alphaf1.boundaryField()*phi1.boundaryField() - + alphaf2.boundaryField()*phi2.boundaryField() - ) - )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) - ); - - // Construct the compressibility parts of the pressure equation - fvScalarMatrix pEqnComp1(p_rgh, dimVolume/dimTime); - fvScalarMatrix pEqnComp2(p_rgh, dimVolume/dimTime); - { - // Density variation - if (!phase1.isochoric() || !phase1.pure()) - { - pEqnComp1 += - ( - fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi()) - - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) - )/rho1; - } - if (!phase2.isochoric() || !phase2.pure()) - { - pEqnComp2 += - ( - fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi()) - - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) - )/rho2; - } - - // Compressibility - if (!phase1.incompressible()) - { - if (pimple.transonic()) - { - const surfaceScalarField phid1 - ( - IOobject::groupName("phid", phase1.name()), - fvc::interpolate(psi1)*phi1 - ); - - pEqnComp1 += - correction - ( - (alpha1/rho1)* - ( - psi1*fvm::ddt(p_rgh) - + fvm::div(phid1, p_rgh) - - fvm::Sp(fvc::div(phid1), p_rgh) - ) - ); - - pEqnComp1.relax(); - } - else - { - pEqnComp1 += (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); - } - } - if (!phase2.incompressible()) - { - if (pimple.transonic()) - { - const surfaceScalarField phid2 - ( - IOobject::groupName("phid", phase2.name()), - fvc::interpolate(psi2)*phi2 - ); - - pEqnComp2 += - correction - ( - (alpha2/rho2)* - ( - psi2*fvm::ddt(p_rgh) - + fvm::div(phid2, p_rgh) - - fvm::Sp(fvc::div(phid2), p_rgh) - ) - ); - - pEqnComp2.relax(); - } - else - { - pEqnComp2 += (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)); - } - } - - // Option sources - if (fvOptions.appliesToField(rho1.name())) - { - pEqnComp1 -= (fvOptions(alpha1, rho1) & rho1)/rho1; - } - if (fvOptions.appliesToField(rho2.name())) - { - pEqnComp2 -= (fvOptions(alpha2, rho2) & rho2)/rho2; - } - - // Mass transfer - PtrList dmdts(fluid.dmdts()); - if (dmdts.set(0)) - { - pEqnComp1 -= dmdts[0]/rho1; - } - if (dmdts.set(1)) - { - pEqnComp2 -= dmdts[1]/rho2; - } - } - - // Cache p prior to solve for density update - const volScalarField p_rgh_0("p_rgh_0", p_rgh); - - // Iterate over the pressure equation to correct for non-orthogonality - while (pimple.correctNonOrthogonal()) - { - // Construct the transport part of the pressure equation - fvScalarMatrix pEqnIncomp - ( - fvc::div(phiHbyA) - - fvm::laplacian(rAUf, p_rgh) - ); - - if (fluid.incompressible()) - { - pEqnIncomp.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); - } - - // Solve - solve(pEqnIncomp + pEqnComp1 + pEqnComp2); - - // Correct fluxes and velocities on last non-orthogonal iteration - if (pimple.finalNonOrthogonalIter()) - { - phi = phiHbyA + pEqnIncomp.flux(); - - surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); - - // Partial-elimination phase-flux corrector - { - const surfaceScalarField phi1s - ( - phiHbyA1 + alpharAUf1*mSfGradp - ); - - const surfaceScalarField phi2s - ( - phiHbyA2 + alpharAUf2*mSfGradp - ); - - surfaceScalarField phir - ( - ((phi1s + phiKd1*phi2s) - (phi2s + phiKd2*phi1s)) - /(1 - phiKd1*phiKd2) - ); - - phi1 = phi + alphaf2*phir; - phi2 = phi - alphaf1*phir; - } - - // Set the phase dilatation rates - phase1.divU(-pEqnComp1 & p_rgh); - phase2.divU(-pEqnComp2 & p_rgh); - - // Optionally relax pressure for velocity correction - p_rgh.relax(); - - mSfGradp = pEqnIncomp.flux()/rAUf; - - // Partial-elimination phase-velocity corrector - { - const volVectorField Us1 - ( - HbyA1 - + fvc::reconstruct(alpharAUf1*mSfGradp - phigF1) - ); - - const volVectorField Us2 - ( - HbyA2 - + fvc::reconstruct(alpharAUf2*mSfGradp - phigF2) - ); - - const volVectorField U - ( - alpha1*(Us1 + rAUKd1*U2) + alpha2*(Us2 + rAUKd2*U1) - ); - - const volVectorField Ur - ( - ((1 - rAUKd2)*Us1 - (1 - rAUKd1)*Us2)/(1 - rAUKd1*rAUKd2) - ); - - U1 = U + alpha2*Ur; - U1.correctBoundaryConditions(); - fvOptions.correct(U1); - - U2 = U - alpha1*Ur; - U2.correctBoundaryConditions(); - fvOptions.correct(U2); - } - } - } - - // Update and limit the static pressure - p = max(p_rgh + rho*gh, pMin); - - // Account for static pressure reference - if (p_rgh.needReference() && fluid.incompressible()) - { - p += dimensionedScalar - ( - "p", - p.dimensions(), - pRefValue - getRefCellValue(p, pRefCell) - ); - } - - // Limit p_rgh - p_rgh = p - rho*gh; - - // Update densities from change in p_rgh - rho1 += psi1*(p_rgh - p_rgh_0); - rho2 += psi2*(p_rgh - p_rgh_0); - - // Correct p_rgh for consistency with p and the updated densities - rho = fluid.rho(); - p_rgh = p - rho*gh; - p_rgh.correctBoundaryConditions(); -} - -if (!fluid.implicitPhasePressure()) -{ - rAUs.clear(); -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H deleted file mode 100644 index 4c36b6f2c8..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H +++ /dev/null @@ -1,43 +0,0 @@ -Info<< "Constructing face momentum equations" << endl; - -fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime); -fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime); - -{ - fluid.momentumTransfer(); // !!! Update coefficients shouldn't be necessary - // This should be done on demand - - autoPtr - momentumTransferPtr(fluid.momentumTransferf()); - - phaseSystem::momentumTransferTable& - momentumTransfer(momentumTransferPtr()); - - { - U1Eqn = - ( - phase1.UfEqn() - == - *momentumTransfer[phase1.name()] - + fvOptions(alpha1, rho1, U1) - ); - U1Eqn.relax(); - fvOptions.constrain(U1Eqn); - U1.correctBoundaryConditions(); - fvOptions.correct(U1); - } - - { - U2Eqn = - ( - phase2.UfEqn() - == - *momentumTransfer[phase2.name()] - + fvOptions(alpha2, rho2, U2) - ); - U2Eqn.relax(); - fvOptions.constrain(U2Eqn); - U2.correctBoundaryConditions(); - fvOptions.correct(U2); - } -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H deleted file mode 100644 index e679cbd107..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H +++ /dev/null @@ -1,382 +0,0 @@ -surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1)); -surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1); - -surfaceScalarField alphaRhof10 -( - "alphaRhof10", - fvc::interpolate - ( - max(alpha1.oldTime(), phase1.residualAlpha()) - *rho1.oldTime() - ) -); - -surfaceScalarField alphaRhof20 -( - "alphaRhof20", - fvc::interpolate - ( - max(alpha2.oldTime(), phase2.residualAlpha()) - *rho2.oldTime() - ) -); - -// Drag coefficient -const surfaceScalarField Kdf("Kdf", fluid.Kdf()); - -// Diagonal coefficients -PtrList AFfs(fluid.AFfs()); - -rAUfs.clear(); - -rAUfs.append -( - new surfaceScalarField - ( - IOobject::groupName("rAUf", phase1.name()), - 1.0 - /( - byDt(alphaRhof10) - + fvc::interpolate(U1Eqn.A()) - + AFfs[0] - ) - ) -); -rAUfs.append -( - new surfaceScalarField - ( - IOobject::groupName("rAUf", phase2.name()), - 1.0 - /( - byDt(alphaRhof20) - + fvc::interpolate(U2Eqn.A()) - + AFfs[0] - ) - ) -); -const surfaceScalarField& rAUf1 = rAUfs[0]; -const surfaceScalarField& rAUf2 = rAUfs[1]; - -AFfs.clear(); - -// Explicit force fluxes -PtrList phiFfs(fluid.phiFfs(rAUfs)); -const surfaceScalarField& phiFf1 = phiFfs[0]; -const surfaceScalarField& phiFf2 = phiFfs[1]; - -while (pimple.correct()) -{ - volScalarField rho("rho", fluid.rho()); - - // Correct p_rgh for consistency with p and the updated densities - p_rgh = p - rho*gh; - - surfaceScalarField rhof1(fvc::interpolate(rho1)); - surfaceScalarField rhof2(fvc::interpolate(rho2)); - - // Correct fixed-flux BCs to be consistent with the velocity BCs - MRF.correctBoundaryFlux(U1, phi1); - MRF.correctBoundaryFlux(U2, phi2); - - const surfaceScalarField alpharAUf1 - ( - IOobject::groupName("alpharAUf", phase1.name()), - max(alphaf1, phase1.residualAlpha())*rAUf1 - ); - - const surfaceScalarField alpharAUf2 - ( - IOobject::groupName("alpharAUf", phase2.name()), - max(alphaf2, phase2.residualAlpha())*rAUf2 - ); - - // Combined buoyancy and force fluxes - const surfaceScalarField ghSnGradRho - ( - "ghSnGradRho", - ghf*fvc::snGrad(rho)*mesh.magSf() - ); - - const surfaceScalarField phigF1 - ( - alpharAUf1 - *( - ghSnGradRho - - alphaf2*(rhof1 - rhof2)*(g & mesh.Sf()) - ) - + phiFf1 - ); - - const surfaceScalarField phigF2 - ( - alpharAUf2 - *( - ghSnGradRho - - alphaf1*(rhof2 - rhof1)*(g & mesh.Sf()) - ) - + phiFf2 - ); - - // Predicted fluxes - surfaceScalarField phiHbyA1 - ( - IOobject::groupName("phiHbyA", phase1.name()), - phi1 - ); - - phiHbyA1 = - rAUf1 - *( - alphaRhof10*byDt(MRF.absolute(phi1.oldTime())) - + fvc::flux(U1Eqn.H()) - ) - - phigF1; - - surfaceScalarField phiHbyA2 - ( - IOobject::groupName("phiHbyA", phase2.name()), - phi2 - ); - - phiHbyA2 = - rAUf2 - *( - alphaRhof20*byDt(MRF.absolute(phi2.oldTime())) - + fvc::flux(U2Eqn.H()) - ) - - phigF2; - - // Drag fluxes - PtrList phiKdPhifs(fluid.phiKdPhifs(rAUfs)); - const surfaceScalarField& phiKdPhif1 = phiKdPhifs[0]; - const surfaceScalarField& phiKdPhif2 = phiKdPhifs[1]; - - // Total predicted flux - surfaceScalarField phiHbyA - ( - "phiHbyA", - alphaf1*(phiHbyA1 - phiKdPhif1) + alphaf2*(phiHbyA2 - phiKdPhif2) - ); - - MRF.makeRelative(phiHbyA); - - phiKdPhifs.clear(); - - // Construct pressure "diffusivity" - const surfaceScalarField rAUf - ( - "rAUf", - mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2) - ); - - // Update the fixedFluxPressure BCs to ensure flux consistency - setSnGrad - ( - p_rgh.boundaryFieldRef(), - ( - phiHbyA.boundaryField() - - ( - alphaf1.boundaryField()*phi1.boundaryField() - + alphaf2.boundaryField()*phi2.boundaryField() - ) - )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) - ); - - // Construct the compressibility parts of the pressure equation - fvScalarMatrix pEqnComp1(p_rgh, dimVolume/dimTime); - fvScalarMatrix pEqnComp2(p_rgh, dimVolume/dimTime); - { - // Density variation - if (!phase1.isochoric() || !phase1.pure()) - { - pEqnComp1 += - ( - fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi()) - - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) - )/rho1; - } - if (!phase2.isochoric() || !phase2.pure()) - { - pEqnComp2 += - ( - fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi()) - - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) - )/rho2; - } - - // Compressibility - if (!phase1.incompressible()) - { - if (pimple.transonic()) - { - const surfaceScalarField phid1 - ( - IOobject::groupName("phid", phase1.name()), - fvc::interpolate(psi1)*phi1 - ); - - pEqnComp1 += - correction - ( - (alpha1/rho1)* - ( - psi1*fvm::ddt(p_rgh) - + fvm::div(phid1, p_rgh) - - fvm::Sp(fvc::div(phid1), p_rgh) - ) - ); - - pEqnComp1.relax(); - } - else - { - pEqnComp1 += (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); - } - } - if (!phase2.incompressible()) - { - if (pimple.transonic()) - { - const surfaceScalarField phid2 - ( - IOobject::groupName("phid", phase2.name()), - fvc::interpolate(psi2)*phi2 - ); - - pEqnComp2 += - correction - ( - (alpha2/rho2)* - ( - psi2*fvm::ddt(p_rgh) - + fvm::div(phid2, p_rgh) - - fvm::Sp(fvc::div(phid2), p_rgh) - ) - ); - - pEqnComp2.relax(); - } - else - { - pEqnComp2 += (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)); - } - } - - // Option sources - if (fvOptions.appliesToField(rho1.name())) - { - pEqnComp1 -= (fvOptions(alpha1, rho1) & rho1)/rho1; - } - if (fvOptions.appliesToField(rho2.name())) - { - pEqnComp2 -= (fvOptions(alpha2, rho2) & rho2)/rho2; - } - - // Mass transfer - PtrList dmdts(fluid.dmdts()); - if (dmdts.set(0)) - { - pEqnComp1 -= dmdts[0]/rho1; - } - if (dmdts.set(1)) - { - pEqnComp2 -= dmdts[1]/rho2; - } - } - - // Cache p prior to solve for density update - const volScalarField p_rgh_0("p_rgh_0", p_rgh); - - // Iterate over the pressure equation to correct for non-orthogonality - while (pimple.correctNonOrthogonal()) - { - // Construct the transport part of the pressure equation - fvScalarMatrix pEqnIncomp - ( - fvc::div(phiHbyA) - - fvm::laplacian(rAUf, p_rgh) - ); - - if (fluid.incompressible()) - { - pEqnIncomp.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); - } - - // Solve - solve(pEqnIncomp + pEqnComp1 + pEqnComp2); - - // Correct fluxes and velocities on last non-orthogonal iteration - if (pimple.finalNonOrthogonalIter()) - { - phi = phiHbyA + pEqnIncomp.flux(); - - surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); - - // Partial-elimination phase-flux corrector - { - const surfaceScalarField phi1s - ( - phiHbyA1 + alpharAUf1*mSfGradp - ); - - const surfaceScalarField phi2s - ( - phiHbyA2 + alpharAUf2*mSfGradp - ); - - surfaceScalarField phir - ( - ((phi1s + rAUf1*Kdf*phi2s) - (phi2s + rAUf2*Kdf*phi1s)) - /(1 - rAUf1*rAUf2*sqr(Kdf)) - ); - - phi1 = phi + alphaf2*phir; - phi2 = phi - alphaf1*phir; - } - - U1 = fvc::reconstruct(MRF.absolute(phi1)); - U1.correctBoundaryConditions(); - fvOptions.correct(U1); - - U2 = fvc::reconstruct(MRF.absolute(phi2)); - U2.correctBoundaryConditions(); - fvOptions.correct(U2); - - // Set the phase dilatation rates - phase1.divU(-pEqnComp1 & p_rgh); - phase2.divU(-pEqnComp2 & p_rgh); - } - } - - // Update and limit the static pressure - p = max(p_rgh + rho*gh, pMin); - - // Account for static pressure reference - if (p_rgh.needReference() && fluid.incompressible()) - { - p += dimensionedScalar - ( - "p", - p.dimensions(), - pRefValue - getRefCellValue(p, pRefCell) - ); - } - - // Limit p_rgh - p_rgh = p - rho*gh; - - // Update densities from change in p_rgh - rho1 += psi1*(p_rgh - p_rgh_0); - rho2 += psi2*(p_rgh - p_rgh_0); - - // Correct p_rgh for consistency with p and the updated densities - rho = fluid.rho(); - p_rgh = p - rho*gh; - p_rgh.correctBoundaryConditions(); -} - -if (!fluid.implicitPhasePressure()) -{ - rAUfs.clear(); -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C deleted file mode 100644 index feab05e8b3..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C +++ /dev/null @@ -1,143 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -Application - reactingTwoPhaseEulerFoam - -Description - Solver for a system of 2 compressible fluid phases with a common pressure, - but otherwise separate properties. The type of phase model is run time - selectable and can optionally represent multiple species and in-phase - reactions. The phase system is also run time selectable and can optionally - represent different types of momentum, heat and mass transfer. - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "twoPhaseSystem.H" -#include "phaseCompressibleMomentumTransportModel.H" -#include "pimpleControl.H" -#include "localEulerDdtScheme.H" -#include "fvcSmooth.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ - #include "postProcess.H" - - #include "setRootCaseLists.H" - #include "createTime.H" - #include "createMesh.H" - #include "createControl.H" - #include "createTimeControls.H" - #include "createFields.H" - #include "createFieldRefs.H" - - if (!LTS) - { - #include "CourantNo.H" - #include "setInitialDeltaT.H" - } - - Switch faceMomentum - ( - pimple.dict().lookupOrDefault("faceMomentum", false) - ); - - #include "createRDeltaTf.H" - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (pimple.run(runTime)) - { - #include "readTimeControls.H" - - int nEnergyCorrectors - ( - pimple.dict().lookupOrDefault("nEnergyCorrectors", 1) - ); - - if (LTS) - { - #include "setRDeltaT.H" - if (faceMomentum) - { - #include "setRDeltaTf.H" - } - } - else - { - #include "CourantNos.H" - #include "setDeltaT.H" - } - - runTime++; - Info<< "Time = " << runTime.timeName() << nl << endl; - - // --- Pressure-velocity PIMPLE corrector loop - while (pimple.loop()) - { - fluid.solve(rAUs, rAUfs); - fluid.correct(); - fluid.correctContinuityError(); - - #include "YEqns.H" - - if (faceMomentum) - { - #include "pUf/UEqns.H" - #include "EEqns.H" - #include "pUf/pEqn.H" - } - else - { - #include "pU/UEqns.H" - #include "EEqns.H" - #include "pU/pEqn.H" - } - - fluid.correctKinematics(); - - if (pimple.turbCorr()) - { - fluid.correctTurbulence(); - } - } - - runTime.write(); - - Info<< "ExecutionTime = " - << runTime.elapsedCpuTime() - << " s\n\n" << endl; - } - - Info<< "End\n" << endl; - - return 0; -} - - -// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaT.H deleted file mode 100644 index 2f4fc79544..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaT.H +++ /dev/null @@ -1,37 +0,0 @@ -{ - volScalarField& rDeltaT = trDeltaT.ref(); - - const dictionary& pimpleDict = pimple.dict(); - - scalar maxCo - ( - pimpleDict.lookupOrDefault("maxCo", 0.2) - ); - - scalar maxDeltaT - ( - pimpleDict.lookupOrDefault("maxDeltaT", great) - ); - - scalar rDeltaTSmoothingCoeff - ( - pimpleDict.lookupOrDefault("rDeltaTSmoothingCoeff", 0.02) - ); - - // Set the reciprocal time-step from the local Courant number - rDeltaT.ref() = max - ( - 1/dimensionedScalar(dimTime, maxDeltaT), - fvc::surfaceSum(max(mag(phi1), mag(phi2)))()() - /((2*maxCo)*mesh.V()) - ); - - // Update the boundary values of the reciprocal time-step - rDeltaT.correctBoundaryConditions(); - - fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff); - - Info<< "Flow time scale min/max = " - << gMin(1/rDeltaT.primitiveField()) - << ", " << gMax(1/rDeltaT.primitiveField()) << endl; -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/Make/files deleted file mode 100644 index 7d143eb42e..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/Make/files +++ /dev/null @@ -1,5 +0,0 @@ -twoPhaseSystem.C -twoPhaseSystemNew.C -twoPhaseSystems.C - -LIB = $(FOAM_LIBBIN)/libtwoPhaseSystem diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/Make/options deleted file mode 100644 index e41b231125..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/Make/options +++ /dev/null @@ -1,19 +0,0 @@ -EXE_INC = \ - -I../../interfacialModels/lnInclude \ - -I../../interfacialCompositionModels/lnInclude \ - -I../../phaseSystems/lnInclude \ - -I../../multiphaseCompressibleMomentumTransportModels/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ - -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \ - -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \ - -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \ - -I$(LIB_SRC)/combustionModels/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude - -LIB_LIBS = \ - -lcombustionModels \ - -lphaseSystem diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C deleted file mode 100644 index 8949d4c9f6..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ /dev/null @@ -1,333 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2020 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "twoPhaseSystem.H" -#include "dragModel.H" -#include "virtualMassModel.H" - -#include "MULES.H" -#include "subCycle.H" -#include "UniformField.H" - -#include "fvcDdt.H" -#include "fvcDiv.H" -#include "fvcSnGrad.H" -#include "fvcFlux.H" -#include "fvcSup.H" - -#include "fvmDdt.H" -#include "fvmLaplacian.H" -#include "fvmSup.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(twoPhaseSystem, 0); - defineRunTimeSelectionTable(twoPhaseSystem, dictionary); -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::twoPhaseSystem::twoPhaseSystem -( - const fvMesh& mesh -) -: - phaseSystem(mesh), - phase1_(phaseModels_[0]), - phase2_(phaseModels_[1]) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::twoPhaseSystem::~twoPhaseSystem() -{} - - -// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // - -Foam::tmp -Foam::twoPhaseSystem::sigma() const -{ - return sigma - ( - phasePairKey(phase1().name(), phase2().name()) - ); -} - - -Foam::tmp -Foam::twoPhaseSystem::Kd() const -{ - return Kd - ( - phasePairKey(phase1().name(), phase2().name()) - ); -} - - -Foam::tmp -Foam::twoPhaseSystem::Kdf() const -{ - return Kdf - ( - phasePairKey(phase1().name(), phase2().name()) - ); -} - - -Foam::tmp -Foam::twoPhaseSystem::Vm() const -{ - return Vm - ( - phasePairKey(phase1().name(), phase2().name()) - ); -} - - -void Foam::twoPhaseSystem::solve -( - const PtrList& rAUs, - const PtrList& rAUfs -) -{ - volScalarField& alpha1 = phase1_; - volScalarField& alpha2 = phase2_; - - const dictionary& alphaControls = mesh_.solverDict(alpha1.name()); - - label nAlphaSubCycles(alphaControls.lookup