diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C index 4d0927a6af..363a8fd2dc 100644 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C +++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C @@ -109,6 +109,9 @@ int main(int argc, char *argv[]) surfaceScalarField phiv_pos(U_pos & mesh.Sf()); surfaceScalarField phiv_neg(U_neg & mesh.Sf()); + fvc::makeRelative(phiv_pos, U); + fvc::makeRelative(phiv_neg, U); + volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi)); surfaceScalarField cSf_pos ( @@ -161,17 +164,9 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; mesh.movePoints(motionPtr->newPoints()); - phiv_pos = U_pos & mesh.Sf(); - phiv_neg = U_neg & mesh.Sf(); - fvc::makeRelative(phiv_pos, U); - fvc::makeRelative(phiv_neg, U); - phiv_neg -= mesh.phi(); - phiv_pos *= a_pos; - phiv_neg *= a_neg; - aphiv_pos = phiv_pos - aSf; - aphiv_neg = phiv_neg + aSf; - surfaceScalarField phi("phi", aphiv_pos*rho_pos + aphiv_neg*rho_neg); + phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; + Info<< phi.boundaryField()[0] << endl; surfaceVectorField phiUp ( @@ -183,6 +178,7 @@ int main(int argc, char *argv[]) ( aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos) + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg) + + mesh.phi()*(a_pos*p_pos + a_neg*p_neg) + aSf*p_pos - aSf*p_neg ); diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H new file mode 100644 index 0000000000..0c17353ec4 --- /dev/null +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H @@ -0,0 +1,12 @@ +{ + solve + ( + fvm::ddt(rho, e) + + fvm::div(phi, e) + - fvm::laplacian(turbulence->alphaEff(), e) + == + - p*fvc::div(phi/fvc::interpolate(rho) + mesh.phi()) + ); + + thermo.correct(); +} diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C index bac588c738..3f659d44ec 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C @@ -76,34 +76,23 @@ Foam::tmp Foam::dragModels::GidaspowErgunWenYu::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); // Wen and Yu (1966) - tmp tKWenYu(0.75*Cds*phase2_.rho()*Ur*bp/d); - volScalarField& KWenYu = tKWenYu(); - - // Ergun - forAll (beta, cellj) - { - if (beta[cellj] <= 0.8) - { - KWenYu[cellj] = - 150.0*alpha_[cellj]*phase2_.nu().value()*phase2_.rho().value() - /sqr(beta[cellj]*d[cellj]) - + 1.75*phase2_.rho().value()*Ur[cellj] - /(beta[cellj]*d[cellj]); - } - } - - return tKWenYu; + return + ( + pos(beta - 0.8) + *(0.75*Cds*phase2_.rho()*Ur*bp/d) + + neg(beta - 0.8) + *( + 150.0*alpha_*phase2_.nu()*phase2_.rho()/(sqr(beta*d)) + + 1.75*phase2_.rho()*Ur/(beta*d) + ) + ); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C index 99b5cacdf2..0a624c6b0a 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C @@ -75,15 +75,11 @@ Foam::tmp Foam::dragModels::GidaspowSchillerNaumann::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d(); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C index 79259cbe05..f68e5379d2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C @@ -72,15 +72,11 @@ Foam::tmp Foam::dragModels::SchillerNaumann::K ) const { volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phase2_.rho()*Ur/phase1_.d(); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C index 4f971b4bb7..541868cf87 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C @@ -73,15 +73,11 @@ Foam::tmp Foam::dragModels::SyamlalOBrien::K { volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6))); volScalarField A(pow(beta, 4.14)); - volScalarField B(0.8*pow(beta, 1.28)); - - forAll (beta, celli) - { - if (beta[celli] > 0.85) - { - B[celli] = pow(beta[celli], 2.65); - } - } + volScalarField B + ( + neg(beta - 0.85)*(0.8*pow(beta, 1.28)) + + pos(beta - 0.85)*(pow(beta, 2.65)) + ); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C index f1126125bb..bf279a3990 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C @@ -75,15 +75,11 @@ Foam::tmp Foam::dragModels::WenYu::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d(); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C index 4320bbad35..d458f07b75 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C @@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf } } + muff.correctBoundaryConditions(); + return tmuf; } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H index 1ee149c48e..bce95f7d5a 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H @@ -21,16 +21,27 @@ forAllIter(PtrDictionary, fluid.phases(), iter) + fvm::div(phase.phiAlpha(), U) - fvm::Sp(fvc::ddt(alpha) + fvc::div(phase.phiAlpha()), U) ) - - fvm::laplacian(alpha*nuEff, U) - - fvc::div - ( - alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/), - "div(Rc)" - ) - == - - fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U) - //- (alpha*phase.rho())*fluid.lift(phase) - + (alpha/phase.rho())*fluid.Svm(phase) + - fvm::laplacian(alpha*nuEff, U) + - fvc::div + ( + alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/), + "div(Rc)" + ) + == + - fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U) + //- (alpha*phase.rho())*fluid.lift(phase) + + (alpha/phase.rho())*fluid.Svm(phase) + - fvm::Sp + ( + slamDampCoeff + *max + ( + mag(U.dimensionedInternalField()) - maxSlamVelocity, + dimensionedScalar("U0", dimVelocity, 0) + ) + /pow(mesh.V(), 1.0/3.0), + U + ) ) ); mrfZones.addCoriolis(alpha, UEqns[phasei]); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H index 484ceae4ac..79fa72d758 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H @@ -53,6 +53,18 @@ phi += fvc::interpolate(alpha)*phase.phi(); } + scalar slamDampCoeff + ( + fluid.lookupOrDefault("slamDampCoeff", 1) + ); + + dimensionedScalar maxSlamVelocity + ( + "maxSlamVelocity", + dimVelocity, + fluid.lookupOrDefault("maxSlamVelocity", GREAT) + ); + // dimensionedScalar pMin // ( // "pMin", diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCoeffs.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCoeffs.H deleted file mode 100644 index 436b68a2f2..0000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCoeffs.H +++ /dev/null @@ -1,80 +0,0 @@ -volScalarField dragCoeff -( - IOobject - ( - "dragCoeff", - runTime.timeName(), - mesh - ), - mesh, - dimensionedScalar("dragCoeff", dimensionSet(1, -3, -1, 0, 0), 0) -); - -volVectorField liftForce -( - IOobject - ( - "liftForce", - runTime.timeName(), - mesh - ), - mesh, - dimensionedVector("liftForce", dimensionSet(1, -2, -2, 0, 0), vector::zero) -); - -volScalarField heatTransferCoeff -( - IOobject - ( - "heatTransferCoeff", - runTime.timeName(), - mesh - ), - mesh, - dimensionedScalar("heatTransferCoeff", dimensionSet(1, -1, -3, -1, 0), 0) -); - -{ - volVectorField Ur = U1 - U2; - volScalarField magUr = mag(Ur); - - if (dispersedPhase == "1") - { - dragCoeff = drag1->K(magUr); - heatTransferCoeff = heatTransfer1->K(magUr); - } - else if (dispersedPhase == "2") - { - dragCoeff = drag2->K(magUr); - heatTransferCoeff = heatTransfer2->K(magUr); - } - else if (dispersedPhase == "both") - { - dragCoeff = - ( - alpha2*drag1->K(magUr) - + alpha1*drag2->K(magUr) - ); - - heatTransferCoeff = - ( - alpha2*heatTransfer1->K(magUr) - + alpha1*heatTransfer2->K(magUr) - ); - } - else - { - FatalErrorIn(args.executable()) - << "dispersedPhase: " << dispersedPhase << " is incorrect" - << exit(FatalError); - } - - volScalarField alphaCoeff - ( - (alpha1 + minInterfaceAlpha)*(alpha2 + minInterfaceAlpha) - ); - dragCoeff *= alphaCoeff; - heatTransferCoeff *= alphaCoeff; - - liftForce = Cl*(alpha1*rho1 + alpha2*rho2)*(Ur ^ fvc::curl(U)); -} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C index e34f219ebc..b52b766de9 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C @@ -75,34 +75,23 @@ Foam::tmp Foam::dragModels::GidaspowErgunWenYu::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); // Wen and Yu (1966) - tmp tKWenYu = 0.75*Cds*phase2_.rho()*Ur*bp/d; - volScalarField& KWenYu = tKWenYu(); - - // Ergun - forAll (beta, cellj) - { - if (beta[cellj] <= 0.8) - { - KWenYu[cellj] = - 150.0*phase1_[cellj]*phase2_.nu().value()*phase2_.rho().value() - /sqr(beta[cellj]*d[cellj]) - + 1.75*phase2_.rho().value()*Ur[cellj] - /(beta[cellj]*d[cellj]); - } - } - - return tKWenYu; + return + ( + pos(beta - 0.8) + *(0.75*Cds*phase2_.rho()*Ur*bp/d) + + neg(beta - 0.8) + *( + 150.0*phase1_*phase2_.nu()*phase2_.rho()/(sqr(beta*d)) + + 1.75*phase2_.rho()*Ur/(beta*d) + ) + ); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C index ecd25eb386..25de6b1dc7 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C @@ -74,15 +74,11 @@ Foam::tmp Foam::dragModels::GidaspowSchillerNaumann::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d(); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C index 3259b5d0e0..d4d77a2bb3 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C @@ -71,15 +71,11 @@ Foam::tmp Foam::dragModels::SchillerNaumann::K ) const { volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phase2_.rho()*Ur/phase1_.d(); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C index 736be885e0..07f85b38a1 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C @@ -72,15 +72,11 @@ Foam::tmp Foam::dragModels::SyamlalOBrien::K { volScalarField beta(max(phase2_, scalar(1.0e-6))); volScalarField A(pow(beta, 4.14)); - volScalarField B(0.8*pow(beta, 1.28)); - - forAll (beta, celli) - { - if (beta[celli] > 0.85) - { - B[celli] = pow(beta[celli], 2.65); - } - } + volScalarField B + ( + neg(beta - 0.85)*(0.8*pow(beta, 1.28)) + + pos(beta - 0.85)*(pow(beta, 2.65)) + ); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C index 6cfc119c1c..6319276408 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C @@ -74,15 +74,11 @@ Foam::tmp Foam::dragModels::WenYu::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d(); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C index 4320bbad35..d458f07b75 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C @@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf } } + muff.correctBoundaryConditions(); + return tmuf; } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C index 4ca01f1a3a..c1d3126352 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C @@ -82,7 +82,6 @@ int main(int argc, char *argv[]) rho = fluid.rho(); #include "zonePhaseVolumes.H" - //#include "interfacialCoeffs.H" //#include "TEqns.H" #include "UEqns.H" diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 6ae60e42dc..2224030613 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -629,8 +629,9 @@ Foam::multiphaseSystem::dragCoeffs() const ( max ( - fvc::average(dm.phase1())*fvc::average(dm.phase2()), - //dm.phase1()*dm.phase2(), + //fvc::average(dm.phase1()*dm.phase2()), + //fvc::average(dm.phase1())*fvc::average(dm.phase2()), + dm.phase1()*dm.phase2(), dm.residualPhaseFraction() ) *dm.K diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H index 744faab968..b93ea3a020 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H @@ -37,8 +37,10 @@ U0s.set(phasei, new volVectorField(phase.U())); phi0s.set(phasei, new surfaceScalarField(phase.phi())); + mrfZones.absoluteFlux(phi0s[phasei]); phasei++; + } surfaceScalarField phi0 @@ -67,6 +69,8 @@ phase.U() = rAUs[phasei]*UEqns[phasei].H(); + phase.phi().oldTime(); + mrfZones.absoluteFlux(phase.phi().oldTime()); mrfZones.absoluteFlux(phase.phi()); phase.phi() = ( @@ -74,7 +78,8 @@ + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi()) ); mrfZones.relativeFlux(phase.phi()); - surfaceScalarField pphi0("pphi0", phase.phi()); + mrfZones.relativeFlux(phase.phi().oldTime()); + surfaceScalarField pphi0("pphi0", 1.0*phase.phi()); pphi0 += rAlphaAUfs[phasei]*(g & mesh.Sf()); multiphaseSystem::dragModelTable::const_iterator dmIter = @@ -196,7 +201,7 @@ p.relax(); mSfGradp = pEqnIncomp.flux()/Dp; - U = dimensionedVector("U", dimVelocity, vector::zero); + U = dimensionedVector("U", dimVelocity, vector::zero); phasei = 0; forAllIter(PtrDictionary, fluid.phases(), iter) @@ -254,9 +259,17 @@ phasej++; } - phase.U() += - (1.0/phase.rho()) - *rAUs[phasei]*(*dcIter())*U0s[phasej]; + // phase.U() += + // (1.0/phase.rho())*rAUs[phasei]*(*dcIter()) + // *U0s[phasej]; + + phase.U() += fvc::reconstruct + ( + fvc::interpolate + ( + (1.0/phase.rho())*rAUs[phasei]*(*dcIter()) + )*phi0s[phasej] + ); } } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseModel/phaseModel/phaseModel.C index bf65fe659c..543927e29b 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseModel/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseModel/phaseModel/phaseModel.C @@ -212,10 +212,10 @@ bool Foam::phaseModel::read(const dictionary& phaseDict) //if (nuModel_->read(phaseDict_)) { - phaseDict_.lookup("nu") >> nu_; - phaseDict_.lookup("kappa") >> kappa_; - phaseDict_.lookup("Cp") >> Cp_; - phaseDict_.lookup("rho") >> rho_; + phaseDict_.lookup("nu") >> nu_.value(); + phaseDict_.lookup("kappa") >> kappa_.value(); + phaseDict_.lookup("Cp") >> Cp_.value(); + phaseDict_.lookup("rho") >> rho_.value(); return true; } diff --git a/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C b/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C index 2dfb19118b..cfc402e14a 100644 --- a/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C +++ b/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C @@ -32,6 +32,7 @@ Description #include "volFields.H" #include "PatchEdgeFaceWave.H" #include "patchEdgeFaceInfo.H" +#include "patchDist.H" using namespace Foam; @@ -49,80 +50,115 @@ int main(int argc, char *argv[]) // Get name of patch const word patchName = args[1]; - const polyPatch& patch = patches[patchName]; - // Data on all edges and faces - List allEdgeInfo(patch.nEdges()); - List allFaceInfo(patch.size()); - - // Initial seed - DynamicList