From 6247f0630b325eb17e5a28098a2cc7983bb9314a Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Thu, 21 Mar 2013 17:05:27 +0000
Subject: [PATCH] compressibleInterFoam: Improve handling of compressibility
and instate support for transonic flow
---
.../multiphase/compressibleInterFoam/TEqn.H | 7 ++-
.../compressibleInterFoam.C | 2 -
.../compressibleInterFoam/createFields.H | 17 +-----
.../multiphase/compressibleInterFoam/pEqn.H | 55 ++++++++++-------
.../twoPhaseMixtureThermo.H | 10 ++++
.../laminar/depthCharge2D/system/fvSchemes | 1 +
.../constant/transportProperties | 59 -------------------
.../laminar/depthCharge3D/system/fvSchemes | 1 +
8 files changed, 54 insertions(+), 98 deletions(-)
delete mode 100644 tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties
diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
index 92a2bd64d8f..d97e8b2a354 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
@@ -1,11 +1,11 @@
{
- solve
+ fvScalarMatrix TEqn
(
fvm::ddt(rho, T)
+ fvm::div(rhoPhi, T)
- fvm::laplacian(twoPhaseProperties.alphaEff(turbulence->mut()), T)
+ (
- p*fvc::div(phi)
+ fvc::div(fvc::absolute(phi, U), p)
+ fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
)
*(
@@ -14,5 +14,8 @@
)
);
+ TEqn.relax();
+ TEqn.solve();
+
twoPhaseProperties.correct();
}
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 3148382ad2e..a42d9119804 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -104,8 +104,6 @@ int main(int argc, char *argv[])
}
}
- rho = alpha1*rho1 + alpha2*rho2;
-
runTime.write();
Info<< "ExecutionTime = "
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 102d4d57a99..c7289b23f90 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -38,9 +38,9 @@
volScalarField& p = twoPhaseProperties.p();
volScalarField& T = twoPhaseProperties.T();
- const volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
+ volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
const volScalarField& psi1 = twoPhaseProperties.thermo1().psi();
- const volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
+ volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
const volScalarField& psi2 = twoPhaseProperties.thermo2().psi();
volScalarField rho
@@ -93,18 +93,5 @@
compressible::turbulenceModel::New(rho, U, rhoPhi, twoPhaseProperties)
);
- Info<< "Creating field dpdt\n" << endl;
- volScalarField dpdt
- (
- IOobject
- (
- "dpdt",
- runTime.timeName(),
- mesh
- ),
- mesh,
- dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
- );
-
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 71745c0168b..73babb08f04 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -26,28 +26,44 @@
tmp<fvScalarMatrix> p_rghEqnComp1;
tmp<fvScalarMatrix> p_rghEqnComp2;
- //if (transonic)
- //{
- //}
- //else
+ if (pimple.transonic())
{
surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
+ p_rghEqnComp1 =
+ fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
+ + correction
+ (
+ psi1*fvm::ddt(p_rgh)
+ + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+ );
+ deleteDemandDrivenData(p_rghEqnComp1().faceFluxCorrectionPtr());
+ p_rghEqnComp1().relax();
+
+ p_rghEqnComp2 =
+ fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+ + correction
+ (
+ psi2*fvm::ddt(p_rgh)
+ + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+ );
+ deleteDemandDrivenData(p_rghEqnComp2().faceFluxCorrectionPtr());
+ p_rghEqnComp2().relax();
+ }
+ else
+ {
p_rghEqnComp1 =
fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
- + fvc::div(phid1, p_rgh)
- - fvc::Sp(fvc::div(phid1), p_rgh);
+ + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
p_rghEqnComp2 =
fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
- + fvc::div(phid2, p_rgh)
- - fvc::Sp(fvc::div(phid2), p_rgh);
+ + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
}
- // Thermodynamic density needs to be updated by psi*d(p) after the
- // pressure solution - done in 2 parts. Part 1:
- //twoPhaseProperties.rho() -= psi*p_rgh;
+ // Cache p_rgh prior to solve for density update
+ volScalarField p_rgh_0(p_rgh);
while (pimple.correctNonOrthogonal())
{
@@ -69,8 +85,8 @@
if (pimple.finalNonOrthogonalIter())
{
- // Second part of thermodynamic density update
- //twoPhaseProperties.rho() += psi*p_rgh;
+ //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+ //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
dgdt =
(
@@ -88,15 +104,14 @@
p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
- // twoPhaseProperties.correct();
+ // Update densities from change in p_rgh
+ rho1 += psi1*(p_rgh - p_rgh_0);
+ rho2 += psi2*(p_rgh - p_rgh_0);
- Info<< "max(U) " << max(mag(U)).value() << endl;
- Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
+ rho = alpha1*rho1 + alpha2*rho2;
K = 0.5*magSqr(U);
- if (twoPhaseProperties.dpdt())
- {
- dpdt = fvc::ddt(p);
- }
+ Info<< "max(U) " << max(mag(U)).value() << endl;
+ Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
}
diff --git a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H
index b2e636a4c47..608dc91a9b1 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H
@@ -94,6 +94,16 @@ public:
return thermo2_();
}
+ rhoThermo& thermo1()
+ {
+ return thermo1_();
+ }
+
+ rhoThermo& thermo2()
+ {
+ return thermo2_();
+ }
+
//- Update properties
virtual void correct();
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes
index 469090b5e0c..ff41a7c06ee 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes
@@ -34,6 +34,7 @@ divSchemes
div(phid2,p_rgh) Gauss upwind;
div(rho*phi,T) Gauss upwind;
div(rho*phi,K) Gauss upwind;
+ div(phi,p) Gauss upwind;
div(phi,k) Gauss vanLeer;
div((muEff*dev2(T(grad(U))))) Gauss linear;
}
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties
deleted file mode 100644
index 564df56f208..00000000000
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties
+++ /dev/null
@@ -1,59 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| ========= | |
-| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
-| \\ / O peration | Version: dev |
-| \\ / A nd | Web: www.OpenFOAM.org |
-| \\/ M anipulation | |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
- version 2.0;
- format ascii;
- class dictionary;
- location "constant";
- object transportProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-phases (water air);
-
-water
-{
- transportModel Newtonian;
- nu 1e-06;
- rho 1000;
- k 0; // 0.613;
- Cv 4179;
-
- equationOfState
- {
- type perfectFluid;
-
- rho0 1000;
- R 3000;
- }
-}
-
-air
-{
- transportModel Newtonian;
- nu 1.589e-05;
- rho 1;
- k 0; // 2.63e-2;
- Cv 721;
-
- equationOfState
- {
- type perfectFluid;
-
- rho0 0;
- R 287;
- }
-}
-
-pMin pMin [ 1 -1 -2 0 0 0 0 ] 10000;
-
-sigma sigma [ 1 0 -2 0 0 0 0 ] 0.07;
-
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes
index 469090b5e0c..ff41a7c06ee 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes
@@ -34,6 +34,7 @@ divSchemes
div(phid2,p_rgh) Gauss upwind;
div(rho*phi,T) Gauss upwind;
div(rho*phi,K) Gauss upwind;
+ div(phi,p) Gauss upwind;
div(phi,k) Gauss vanLeer;
div((muEff*dev2(T(grad(U))))) Gauss linear;
}
--
GitLab