diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C index ae3a2e8829a6acdc4effe27f543505072e62ff50..3e8e5473638747b839ac8770fac12f8ee0d649cc 100644 --- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C +++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -148,7 +148,7 @@ int main(int argc, char *argv[]) volScalarField rAB(1.0/BEqn.A()); surfaceScalarField rABf("rABf", fvc::interpolate(rAB)); - phiB = fvc::flux(B) + rABf*fvc::ddtCorr(B, phiB); + phiB = fvc::flux(B); while (bpiso.correctNonOrthogonal()) { diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C index 28530b882f60ec7c76b660bbc9d9118524a8d227..6ac069ed5f26cb124686c1ef98bf1b0ae00e6938 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -649,14 +649,14 @@ CoEulerDdtScheme<Type>::fvcDdtUfCorr const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf ) { + const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT())); + if ( U.dimensions() == dimVelocity && Uf.dimensions() == dimDensity*dimVelocity ) { - const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT())); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -676,7 +676,8 @@ CoEulerDdtScheme<Type>::fvcDdtUfCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr + this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime()) + *rDeltaT*phiCorr ) ); } @@ -686,7 +687,32 @@ CoEulerDdtScheme<Type>::fvcDdtUfCorr && Uf.dimensions() == dimDensity*dimVelocity ) { - return fvcDdtUfCorr(U, Uf); + fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime()); + fluxFieldType phiCorr + ( + phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + Uf.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phiUf0, + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { @@ -708,14 +734,14 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr const fluxFieldType& phi ) { + dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); + if ( U.dimensions() == dimVelocity && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -737,8 +763,13 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr) - *rDeltaT*phiCorr + this->fvcDdtPhiCoeff + ( + rhoU0, + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr ) ); } @@ -748,7 +779,31 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - return fvcDdtPhiCorr(U, phi); + fluxFieldType phiCorr + ( + phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + phi.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C index f84e5d5401a0527cb41b96f46c8624843ad65b16..6d6f2061b593464b8017491e0a56d3284e2a4fdf 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -1358,7 +1358,12 @@ CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, mesh().Sf() & Uf.oldTime()) + this->fvcDdtPhiCoeff + ( + rhoU0, + mesh().Sf() & Uf.oldTime(), + rho.oldTime() + ) *( mesh().Sf() & ( @@ -1377,7 +1382,64 @@ CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr && Uf.dimensions() == rho.dimensions()*dimVelocity ) { - return fvcDdtUfCorr(U, Uf); + DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 = + ddt0_<GeometricField<Type, fvPatchField, volMesh>> + ( + "ddtCorrDdt0(" + U.name() + ')', + U.dimensions() + ); + + DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 = + ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>> + ( + "ddtCorrDdt0(" + Uf.name() + ')', + Uf.dimensions() + ); + + dimensionedScalar rDtCoef = rDtCoef_(ddt0); + + if (evaluate(ddt0)) + { + ddt0 = + rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime()) + - offCentre_(ddt0()); + } + + if (evaluate(dUfdt0)) + { + dUfdt0 = + rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime()) + - offCentre_(dUfdt0()); + } + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + U.name() + ',' + Uf.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + mesh().Sf() & Uf.oldTime(), + rho.oldTime() + ) + *( + mesh().Sf() + & ( + (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0())) + - fvc::interpolate + ( + rDtCoef*U.oldTime() + offCentre_(ddt0()) + ) + ) + ) + ) + ); } else { @@ -1453,7 +1515,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phi.oldTime()) + this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime()) *( (rDtCoef*phi.oldTime() + offCentre_(dphidt0())) - fvc::dotInterpolate @@ -1473,7 +1535,57 @@ CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - return fvcDdtPhiCorr(U, phi); + DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 = + ddt0_<GeometricField<Type, fvPatchField, volMesh>> + ( + "ddtCorrDdt0(" + U.name() + ')', + U.dimensions() + ); + + DDt0Field<fluxFieldType>& dphidt0 = + ddt0_<fluxFieldType> + ( + "ddtCorrDdt0(" + phi.name() + ')', + phi.dimensions() + ); + + dimensionedScalar rDtCoef = rDtCoef_(ddt0); + + if (evaluate(ddt0)) + { + ddt0 = + rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime()) + - offCentre_(ddt0()); + } + + if (evaluate(dphidt0)) + { + dphidt0 = + rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime()) + - offCentre_(dphidt0()); + } + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + U.name() + ',' + phi.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime()) + *( + (rDtCoef*phi.oldTime() + offCentre_(dphidt0())) + - fvc::dotInterpolate + ( + mesh().Sf(), + rDtCoef*U.oldTime() + offCentre_(ddt0()) + ) + ) + ) + ); } else { diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C index 49705716a1bd3c4387a3149c854d9a78cddc6bd1..1c328d3f47bf0c64f4ee32dee96da7b58d3f6c7a 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -557,14 +557,14 @@ EulerDdtScheme<Type>::fvcDdtUfCorr const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf ) { + dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); + if ( U.dimensions() == dimVelocity && Uf.dimensions() == rho.dimensions()*dimVelocity ) { - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -584,7 +584,8 @@ EulerDdtScheme<Type>::fvcDdtUfCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr + this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime()) + *rDeltaT*phiCorr ) ); } @@ -594,7 +595,32 @@ EulerDdtScheme<Type>::fvcDdtUfCorr && Uf.dimensions() == rho.dimensions()*dimVelocity ) { - return fvcDdtUfCorr(U, Uf); + fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime()); + fluxFieldType phiCorr + ( + phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + Uf.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phiUf0, + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { @@ -616,14 +642,14 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr const fluxFieldType& phi ) { + dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); + if ( U.dimensions() == dimVelocity && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -645,8 +671,13 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr) - *rDeltaT*phiCorr + this->fvcDdtPhiCoeff + ( + rhoU0, + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr ) ); } @@ -656,7 +687,31 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - return fvcDdtPhiCorr(U, phi); + fluxFieldType phiCorr + ( + phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + phi.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C index 8594425be322d002eeeb93c27a2a9fb63f3b98c4..73bda9c99d26b5d5bbd328d0c9e96e98337c46cc 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -654,14 +654,14 @@ SLTSDdtScheme<Type>::fvcDdtUfCorr const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf ) { + const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT())); + if ( U.dimensions() == dimVelocity && Uf.dimensions() == dimDensity*dimVelocity ) { - const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT())); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -681,7 +681,8 @@ SLTSDdtScheme<Type>::fvcDdtUfCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr + this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime()) + *rDeltaT*phiCorr ) ); } @@ -691,7 +692,32 @@ SLTSDdtScheme<Type>::fvcDdtUfCorr && Uf.dimensions() == dimDensity*dimVelocity ) { - return fvcDdtUfCorr(U, Uf); + fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime()); + fluxFieldType phiCorr + ( + phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + Uf.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phiUf0, + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { @@ -713,14 +739,14 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr const fluxFieldType& phi ) { + dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); + if ( U.dimensions() == dimVelocity && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -742,8 +768,13 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr) - *rDeltaT*phiCorr + this->fvcDdtPhiCoeff + ( + rhoU0, + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr ) ); } @@ -753,7 +784,31 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - return fvcDdtPhiCorr(U, phi); + fluxFieldType phiCorr + ( + phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + phi.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C index 8e2d5c70fea1d5413627f86cbe4eeef8eaa36456..29e704619767a61f68940d0e547b275c9e7bb79b 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License @@ -753,21 +753,21 @@ backwardDdtScheme<Type>::fvcDdtUfCorr const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf ) { + dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); + + 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 && Uf.dimensions() == rho.dimensions()*dimVelocity ) { - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - - scalar deltaT = deltaT_(); - scalar deltaT0 = deltaT0_(U); - - scalar coefft = 1 + deltaT/(deltaT + deltaT0); - scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); - scalar coefft0 = coefft + coefft00; - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -789,7 +789,12 @@ backwardDdtScheme<Type>::fvcDdtUfCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, mesh().Sf() & Uf.oldTime()) + this->fvcDdtPhiCoeff + ( + rhoU0, + mesh().Sf() & Uf.oldTime(), + rho.oldTime() + ) *rDeltaT *( mesh().Sf() @@ -807,7 +812,37 @@ backwardDdtScheme<Type>::fvcDdtUfCorr && Uf.dimensions() == rho.dimensions()*dimVelocity ) { - return fvcDdtUfCorr(U, Uf); + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + Uf.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + mesh().Sf() & Uf.oldTime(), + rho.oldTime() + ) + *rDeltaT + *( + mesh().Sf() + & ( + (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime()) + - fvc::interpolate + ( + coefft0*U.oldTime() + - coefft00*U.oldTime().oldTime() + ) + ) + ) + ) + ); } else { @@ -829,21 +864,21 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr const fluxFieldType& phi ) { + dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); + + 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 && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - - scalar deltaT = deltaT_(); - scalar deltaT0 = deltaT0_(U); - - scalar coefft = 1 + deltaT/(deltaT + deltaT0); - scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); - scalar coefft0 = coefft + coefft00; - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -865,7 +900,7 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phi.oldTime()) + this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime()) *rDeltaT *( (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime()) @@ -884,7 +919,29 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - return fvcDdtPhiCorr(U, phi); + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + phi.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime()) + *rDeltaT + *( + (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime()) + - fvc::dotInterpolate + ( + mesh().Sf(), + coefft0*U.oldTime() - coefft00*U.oldTime().oldTime() + ) + ) + ) + ); } else { diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C index fd4f502ed991ea0fa282ebee81da93fac88470fc..45bc97047eb3dbe74e020b959a555dcc398ed7c8 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License @@ -178,7 +178,7 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff ddtCouplingCoeff -= min ( mag(phiCorr) - /(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)), + *mesh().time().deltaT()*mesh().deltaCoeffs()/mesh().magSf(), scalar(1) ); } @@ -216,6 +216,19 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff } +template<class Type> +tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff +( + const GeometricField<Type, fvPatchField, volMesh>& U, + const fluxFieldType& phi, + const fluxFieldType& phiCorr, + const volScalarField& rho +) +{ + return fvcDdtPhiCoeff(U, phi, phiCorr/fvc::interpolate(rho)); +} + + template<class Type> tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff ( @@ -227,6 +240,23 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff } +template<class Type> +tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff +( + const GeometricField<Type, fvPatchField, volMesh>& U, + const fluxFieldType& phi, + const volScalarField& rho +) +{ + return fvcDdtPhiCoeff + ( + U, + phi, + (phi - fvc::dotInterpolate(mesh().Sf(), U))/fvc::interpolate(rho) + ); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H index 8a1dfbaa8de439fd766febe6aee2123ffef64c29..3d293e0bdd14bf4b30c3a9863edf3973ad369c2b 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License @@ -219,12 +219,27 @@ public: const fluxFieldType& phiCorr ); + tmp<surfaceScalarField> fvcDdtPhiCoeff + ( + const GeometricField<Type, fvPatchField, volMesh>& U, + const fluxFieldType& phi, + const fluxFieldType& phiCorr, + const volScalarField& rho + ); + tmp<surfaceScalarField> fvcDdtPhiCoeff ( const GeometricField<Type, fvPatchField, volMesh>& U, const fluxFieldType& phi ); + tmp<surfaceScalarField> fvcDdtPhiCoeff + ( + const GeometricField<Type, fvPatchField, volMesh>& U, + const fluxFieldType& phi, + const volScalarField& rho + ); + virtual tmp<fluxFieldType> fvcDdtUfCorr ( const GeometricField<Type, fvPatchField, volMesh>& U, diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C index cfa891be09cf1e2da11b6e3203159137bf77b599..94fabc81b8f45033107e67cfc516a4b43abf4977 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -413,14 +413,14 @@ localEulerDdtScheme<Type>::fvcDdtUfCorr const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf ) { + const surfaceScalarField rDeltaT(fvc::interpolate(localRDeltaT())); + if ( U.dimensions() == dimVelocity && Uf.dimensions() == dimDensity*dimVelocity ) { - const surfaceScalarField rDeltaT(fvc::interpolate(localRDeltaT())); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -440,7 +440,8 @@ localEulerDdtScheme<Type>::fvcDdtUfCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr + this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime()) + *rDeltaT*phiCorr ) ); } @@ -450,7 +451,32 @@ localEulerDdtScheme<Type>::fvcDdtUfCorr && Uf.dimensions() == dimDensity*dimVelocity ) { - return fvcDdtUfCorr(U, Uf); + fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime()); + fluxFieldType phiCorr + ( + phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + Uf.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phiUf0, + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { @@ -472,14 +498,14 @@ localEulerDdtScheme<Type>::fvcDdtPhiCorr const fluxFieldType& phi ) { + const surfaceScalarField rDeltaT(fvc::interpolate(localRDeltaT())); + if ( U.dimensions() == dimVelocity && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - GeometricField<Type, fvPatchField, volMesh> rhoU0 ( rho.oldTime()*U.oldTime() @@ -501,8 +527,13 @@ localEulerDdtScheme<Type>::fvcDdtPhiCorr mesh().time().timeName(), mesh() ), - this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr) - *rDeltaT*phiCorr + this->fvcDdtPhiCoeff + ( + rhoU0, + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr ) ); } @@ -512,7 +543,31 @@ localEulerDdtScheme<Type>::fvcDdtPhiCorr && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea ) { - return fvcDdtPhiCorr(U, phi); + fluxFieldType phiCorr + ( + phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) + ); + + return tmp<fluxFieldType> + ( + new fluxFieldType + ( + IOobject + ( + "ddtCorr(" + + rho.name() + ',' + U.name() + ',' + phi.name() + ')', + mesh().time().timeName(), + mesh() + ), + this->fvcDdtPhiCoeff + ( + U.oldTime(), + phi.oldTime(), + phiCorr, + rho.oldTime() + )*rDeltaT*phiCorr + ) + ); } else { diff --git a/tutorials/incompressible/pimpleFoam/LES/channel395/constant/fvOptions b/tutorials/incompressible/pimpleFoam/LES/channel395/constant/fvOptions index d9146368ed7e8f3d7382c4d8259c5c5ad4376d70..af32b36e5a225125300acb58c9f1479aa330abbb 100644 --- a/tutorials/incompressible/pimpleFoam/LES/channel395/constant/fvOptions +++ b/tutorials/incompressible/pimpleFoam/LES/channel395/constant/fvOptions @@ -21,7 +21,7 @@ momentumSource selectionMode all; - fields (U); + fields (U); Ubar (0.1335 0 0); }