From 793fec252113df45a6e085dfef09dae8de054c09 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Fri, 9 May 2008 11:47:33 +0100
Subject: [PATCH] Added ddtPhiCorr support for moving mesh
---
.../multiphase/interDyMFoam/interDyMFoam.C | 7 +-
.../solvers/multiphase/interDyMFoam/pEqn.H | 2 +-
.../EulerDdtScheme/EulerDdtScheme.C | 199 +++++-------
.../backwardDdtScheme/backwardDdtScheme.C | 288 ++++++++----------
4 files changed, 212 insertions(+), 284 deletions(-)
diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index b3046c7fbb9..e46cba7db9a 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
// Make the fluxes absolute
- if (mesh.changing() && correctPhi)
+ if (mesh.changing())
{
fvc::makeAbsolute(phi, U);
}
@@ -99,8 +99,11 @@ int main(int argc, char *argv[])
#include "correctPhi.H"
}
+ // Keep the absolute fluxes for use in ddtPhiCorr
+ surfaceScalarField phiAbs("phiAbs", phi);
+
// Make the fluxes relative to the mesh motion
- if (mesh.changing() && correctPhi)
+ if (mesh.changing())
{
fvc::makeRelative(phi, U);
}
diff --git a/applications/solvers/multiphase/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interDyMFoam/pEqn.H
index 6ae105b3bcf..5bd1a98d40a 100644
--- a/applications/solvers/multiphase/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interDyMFoam/pEqn.H
@@ -9,7 +9,7 @@
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
- //+ fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+ + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
);
phi = phiU +
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index 1a5d4b498dd..ee3381932e5 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------* \
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
@@ -21,7 +21,7 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
+
\*---------------------------------------------------------------------------*/
#include "EulerDdtScheme.H"
@@ -279,7 +279,7 @@ EulerDdtScheme<Type>::fvmDdt
scalar rDeltaT = 1.0/mesh().time().deltaT().value();
fvm.diag() = rDeltaT*mesh().V();
-
+
if (mesh().moving())
{
fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
@@ -314,7 +314,7 @@ EulerDdtScheme<Type>::fvmDdt
scalar rDeltaT = 1.0/mesh().time().deltaT().value();
fvm.diag() = rDeltaT*rho.value()*mesh().V();
-
+
if (mesh().moving())
{
fvm.source() = rDeltaT
@@ -375,50 +375,30 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
(
const volScalarField& rA,
const GeometricField<Type, fvPatchField, volMesh>& U,
- const fluxFieldType& phi
+ const fluxFieldType& phiAbs
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
IOobject ddtIOobject
(
- "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+ "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phiAbs.name() + ')',
mesh().time().timeName(),
mesh()
);
- if (mesh().moving())
- {
- return tmp<fluxFieldType>
- (
- new fluxFieldType
- (
- ddtIOobject,
- mesh(),
- dimensioned<typename flux<Type>::type>
- (
- "0",
- rDeltaT.dimensions()*rA.dimensions()*phi.dimensions(),
- pTraits<typename flux<Type>::type>::zero
- )
- )
- );
- }
- else
- {
- tmp<fluxFieldType> phiCorr =
- phi.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
+ tmp<fluxFieldType> phiCorr =
+ phiAbs.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
- return tmp<fluxFieldType>
+ return tmp<fluxFieldType>
+ (
+ new fluxFieldType
(
- new fluxFieldType
- (
- ddtIOobject,
- fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr())
- *fvc::interpolate(rDeltaT*rA)*phiCorr
- )
- );
- }
+ ddtIOobject,
+ fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime(), phiCorr())
+ *fvc::interpolate(rDeltaT*rA)*phiCorr
+ )
+ );
}
@@ -429,7 +409,7 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
const volScalarField& rA,
const volScalarField& rho,
const GeometricField<Type, fvPatchField, volMesh>& U,
- const fluxFieldType& phi
+ const fluxFieldType& phiAbs
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
@@ -437,111 +417,94 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
IOobject ddtIOobject
(
"ddtPhiCorr("
- + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+ + rA.name() + ','
+ + rho.name() + ','
+ + U.name() + ','
+ + phiAbs.name() + ')',
mesh().time().timeName(),
mesh()
);
- if (mesh().moving())
+ if
+ (
+ U.dimensions() == dimVelocity
+ && phiAbs.dimensions() == dimVelocity*dimArea
+ )
{
return tmp<fluxFieldType>
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned<typename flux<Type>::type>
- (
- "0",
- rDeltaT.dimensions()*rA.dimensions()
- *rho.dimensions()*phi.dimensions(),
- pTraits<typename flux<Type>::type>::zero
+ rDeltaT
+ *fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
+ *(
+ fvc::interpolate(rA*rho.oldTime())*phiAbs.oldTime()
+ - (fvc::interpolate(rA*rho.oldTime()*U.oldTime())
+ & mesh().Sf())
)
)
);
}
- else
+ else if
+ (
+ U.dimensions() == dimVelocity
+ && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
{
- if
+ return tmp<fluxFieldType>
(
- U.dimensions() == dimVelocity
- && phi.dimensions() == dimVelocity*dimArea
- )
- {
- return tmp<fluxFieldType>
+ new fluxFieldType
(
- new fluxFieldType
+ ddtIOobject,
+ rDeltaT
+ *fvcDdtPhiCoeff
(
- ddtIOobject,
- rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
- *(
- fvc::interpolate(rA*rho.oldTime())*phi.oldTime()
- - (fvc::interpolate(rA*rho.oldTime()*U.oldTime())
- & mesh().Sf())
- )
+ U.oldTime(),
+ phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
)
- );
- }
- else if
+ *(
+ fvc::interpolate(rA*rho.oldTime())
+ *phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
+ - (
+ fvc::interpolate
+ (
+ rA*rho.oldTime()*U.oldTime()
+ ) & mesh().Sf()
+ )
+ )
+ )
+ );
+ }
+ else if
+ (
+ U.dimensions() == rho.dimensions()*dimVelocity
+ && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
+ {
+ return tmp<fluxFieldType>
(
- U.dimensions() == dimVelocity
- && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
- )
- {
- return tmp<fluxFieldType>
+ new fluxFieldType
(
- new fluxFieldType
- (
- ddtIOobject,
- rDeltaT
- *fvcDdtPhiCoeff
- (
- U.oldTime(),
- phi.oldTime()/fvc::interpolate(rho.oldTime())
- )
- *(
- fvc::interpolate(rA*rho.oldTime())
- *phi.oldTime()/fvc::interpolate(rho.oldTime())
- - (
- fvc::interpolate
- (
- rA*rho.oldTime()*U.oldTime()
- ) & mesh().Sf()
- )
- )
+ ddtIOobject,
+ rDeltaT
+ *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
+ *(
+ fvc::interpolate(rA)*phiAbs.oldTime()
+ - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
)
- );
- }
- else if
+ )
+ );
+ }
+ else
+ {
+ FatalErrorIn
(
- U.dimensions() == rho.dimensions()*dimVelocity
- && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
- )
- {
- return tmp<fluxFieldType>
- (
- new fluxFieldType
- (
- ddtIOobject,
- rDeltaT
- *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
- *(
- fvc::interpolate(rA)*phi.oldTime()
- - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
- )
- )
- );
- }
- else
- {
- FatalErrorIn
- (
- "EulerDdtScheme<Type>::fvcDdtPhiCorr"
- ) << "dimensions of phi are not correct"
- << abort(FatalError);
+ "EulerDdtScheme<Type>::fvcDdtPhiCorr"
+ ) << "dimensions of phiAbs are not correct"
+ << abort(FatalError);
- return fluxFieldType::null();
- }
+ return fluxFieldType::null();
}
}
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index 0c30c641c3f..f36312b64f3 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -21,7 +21,7 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
+
\*---------------------------------------------------------------------------*/
#include "backwardDdtScheme.H"
@@ -531,58 +531,38 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
mesh()
);
- if (mesh().moving())
- {
- return tmp<fluxFieldType>
- (
- new fluxFieldType
- (
- ddtIOobject,
- mesh(),
- dimensioned<typename flux<Type>::type>
- (
- "0",
- rDeltaT.dimensions()*rA.dimensions()*phi.dimensions(),
- pTraits<typename flux<Type>::type>::zero
- )
- )
- );
- }
- else
- {
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(U);
+ scalar deltaT = deltaT_();
+ scalar deltaT0 = deltaT0_(U);
- scalar coefft = 1 + deltaT/(deltaT + deltaT0);
- scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
- scalar coefft0 = coefft + coefft00;
+ scalar coefft = 1 + deltaT/(deltaT + deltaT0);
+ scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+ scalar coefft0 = coefft + coefft00;
- return tmp<fluxFieldType>
+ return tmp<fluxFieldType>
+ (
+ new fluxFieldType
(
- new fluxFieldType
- (
- ddtIOobject,
- rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+ ddtIOobject,
+ rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+ *(
+ fvc::interpolate(rA)
*(
- fvc::interpolate(rA)
- *(
- coefft0*phi.oldTime()
- - coefft00*phi.oldTime().oldTime()
- )
- - (
- fvc::interpolate
+ coefft0*phi.oldTime()
+ - coefft00*phi.oldTime().oldTime()
+ )
+ - (
+ fvc::interpolate
+ (
+ rA*
(
- rA*
- (
- coefft0*U.oldTime()
- - coefft00*U.oldTime().oldTime()
- )
- ) & mesh().Sf()
- )
+ coefft0*U.oldTime()
+ - coefft00*U.oldTime().oldTime()
+ )
+ ) & mesh().Sf()
)
)
- );
- }
+ )
+ );
}
@@ -593,7 +573,7 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
const volScalarField& rA,
const volScalarField& rho,
const GeometricField<Type, fvPatchField, volMesh>& U,
- const fluxFieldType& phi
+ const fluxFieldType& phiAbs
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
@@ -601,152 +581,134 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
IOobject ddtIOobject
(
"ddtPhiCorr("
- + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+ + rA.name() + ','
+ + rho.name() + ','
+ + U.name() + ','
+ + phiAbs.name() + ')',
mesh().time().timeName(),
mesh()
);
- if (mesh().moving())
+ scalar deltaT = deltaT_();
+ scalar deltaT0 = deltaT0_(U);
+
+ scalar coefft = 1 + deltaT/(deltaT + deltaT0);
+ scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+ scalar coefft0 = coefft + coefft00;
+
+ if
+ (
+ U.dimensions() == dimVelocity
+ && phiAbs.dimensions() == dimVelocity*dimArea
+ )
{
return tmp<fluxFieldType>
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned<typename flux<Type>::type>
- (
- "0",
- rDeltaT.dimensions()*rA.dimensions()
- *rho.dimensions()*phi.dimensions(),
- pTraits<typename flux<Type>::type>::zero
+ rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
+ *(
+ coefft0*fvc::interpolate(rA*rho.oldTime())
+ *phiAbs.oldTime()
+ - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
+ *phiAbs.oldTime().oldTime()
+ - (
+ fvc::interpolate
+ (
+ rA*
+ (
+ coefft0*rho.oldTime()*U.oldTime()
+ - coefft00*rho.oldTime().oldTime()
+ *U.oldTime().oldTime()
+ )
+ ) & mesh().Sf()
+ )
)
)
);
}
- else
+ else if
+ (
+ U.dimensions() == dimVelocity
+ && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
{
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(U);
-
- scalar coefft = 1 + deltaT/(deltaT + deltaT0);
- scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
- scalar coefft0 = coefft + coefft00;
-
- if
+ return tmp<fluxFieldType>
(
- U.dimensions() == dimVelocity
- && phi.dimensions() == dimVelocity*dimArea
- )
- {
- return tmp<fluxFieldType>
+ new fluxFieldType
(
- new fluxFieldType
+ ddtIOobject,
+ rDeltaT
+ *fvcDdtPhiCoeff
(
- ddtIOobject,
- rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
- *(
- coefft0*fvc::interpolate(rA*rho.oldTime())
- *phi.oldTime()
- - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
- *phi.oldTime().oldTime()
- - (
- fvc::interpolate
- (
- rA*
- (
- coefft0*rho.oldTime()*U.oldTime()
- - coefft00*rho.oldTime().oldTime()
- *U.oldTime().oldTime()
- )
- ) & mesh().Sf()
- )
- )
+ U.oldTime(),
+ phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
)
- );
- }
- else if
- (
- U.dimensions() == dimVelocity
- && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
- )
- {
- return tmp<fluxFieldType>
- (
- new fluxFieldType
- (
- ddtIOobject,
- rDeltaT
- *fvcDdtPhiCoeff
- (
- U.oldTime(),
- phi.oldTime()/fvc::interpolate(rho.oldTime())
- )
+ *(
+ fvc::interpolate(rA*rho.oldTime())
*(
- fvc::interpolate(rA*rho.oldTime())
- *(
- coefft0*phi.oldTime()
- /fvc::interpolate(rho.oldTime())
- - coefft00*phi.oldTime().oldTime()
- /fvc::interpolate(rho.oldTime().oldTime())
- )
- - (
- fvc::interpolate
+ coefft0*phiAbs.oldTime()
+ /fvc::interpolate(rho.oldTime())
+ - coefft00*phiAbs.oldTime().oldTime()
+ /fvc::interpolate(rho.oldTime().oldTime())
+ )
+ - (
+ fvc::interpolate
+ (
+ rA*rho.oldTime()*
(
- rA*rho.oldTime()*
- (
- coefft0*U.oldTime()
- - coefft00*U.oldTime().oldTime()
- )
- ) & mesh().Sf()
- )
+ coefft0*U.oldTime()
+ - coefft00*U.oldTime().oldTime()
+ )
+ ) & mesh().Sf()
)
)
- );
- }
- else if
+ )
+ );
+ }
+ else if
+ (
+ U.dimensions() == rho.dimensions()*dimVelocity
+ && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
+ {
+ return tmp<fluxFieldType>
(
- U.dimensions() == rho.dimensions()*dimVelocity
- && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
- )
- {
- return tmp<fluxFieldType>
+ new fluxFieldType
(
- new fluxFieldType
- (
- ddtIOobject,
- rDeltaT
- *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
+ ddtIOobject,
+ rDeltaT
+ *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
+ *(
+ fvc::interpolate(rA)
*(
- fvc::interpolate(rA)
- *(
- coefft0*phi.oldTime()
- - coefft00*phi.oldTime().oldTime()
- )
- - (
- fvc::interpolate
+ coefft0*phiAbs.oldTime()
+ - coefft00*phiAbs.oldTime().oldTime()
+ )
+ - (
+ fvc::interpolate
+ (
+ rA*
(
- rA*
- (
- coefft0*U.oldTime()
- - coefft00*U.oldTime().oldTime()
- )
- ) & mesh().Sf()
- )
+ coefft0*U.oldTime()
+ - coefft00*U.oldTime().oldTime()
+ )
+ ) & mesh().Sf()
)
)
- );
- }
- else
- {
- FatalErrorIn
- (
- "backwardDdtScheme<Type>::fvcDdtPhiCorr"
- ) << "dimensions of phi are not correct"
- << abort(FatalError);
+ )
+ );
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "backwardDdtScheme<Type>::fvcDdtPhiCorr"
+ ) << "dimensions of phiAbs are not correct"
+ << abort(FatalError);
- return fluxFieldType::null();
- }
+ return fluxFieldType::null();
}
}
--
GitLab