Skip to content
Snippets Groups Projects
Commit 380bfb09 authored by andy's avatar andy
Browse files

ENH: Re-worked particle pressure gradient and virtual mass forces

parent a4ac852e
Branches
Tags
No related merge requests found
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "PressureGradientForce.H"
#include "fvcDdt.H"
#include "fvcGrad.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -33,12 +34,13 @@ Foam::PressureGradientForce<CloudType>::PressureGradientForce
(
CloudType& owner,
const fvMesh& mesh,
const dictionary& dict
const dictionary& dict,
const word& forceType
)
:
ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
UName_(this->coeffs().lookup("U")),
gradUPtr_(NULL)
ParticleForce<CloudType>(owner, mesh, dict, forceType, true),
UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
DUcDtInterpPtr_(NULL)
{}
......@@ -50,7 +52,7 @@ Foam::PressureGradientForce<CloudType>::PressureGradientForce
:
ParticleForce<CloudType>(pgf),
UName_(pgf.UName_),
gradUPtr_(NULL)
DUcDtInterpPtr_(NULL)
{}
......@@ -66,18 +68,48 @@ Foam::PressureGradientForce<CloudType>::~PressureGradientForce()
template<class CloudType>
void Foam::PressureGradientForce<CloudType>::cacheFields(const bool store)
{
static word fName("DUcDt");
bool fieldExists = this->mesh().template foundObject<volVectorField>(fName);
if (store)
{
const volVectorField& U = this->mesh().template
lookupObject<volVectorField>(UName_);
gradUPtr_ = fvc::grad(U).ptr();
if (!fieldExists)
{
const volVectorField& Uc = this->mesh().template
lookupObject<volVectorField>(UName_);
volVectorField* DUcDtPtr = new volVectorField
(
fName,
fvc::ddt(Uc) + (Uc & fvc::grad(Uc))
);
DUcDtPtr->store();
}
const volVectorField& DUcDt = this->mesh().template
lookupObject<volVectorField>(fName);
DUcDtInterpPtr_.reset
(
interpolation<vector>::New
(
this->owner().solution().interpolationSchemes(),
DUcDt
).ptr()
);
}
else
{
if (gradUPtr_)
DUcDtInterpPtr_.clear();
if (fieldExists)
{
delete gradUPtr_;
gradUPtr_ = NULL;
const volVectorField& DUcDt = this->mesh().template
lookupObject<volVectorField>(fName);
const_cast<volVectorField&>(DUcDt).checkOut();
}
}
}
......@@ -95,11 +127,24 @@ Foam::forceSuSp Foam::PressureGradientForce<CloudType>::calcCoupled
{
forceSuSp value(vector::zero, 0.0);
const volTensorField& gradU = *gradUPtr_;
value.Su() = mass*p.rhoc()/p.rho()*(p.U() & gradU[p.cell()]);
vector DUcDt =
DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
value.Su() = mass*p.rhoc()/p.rho()*DUcDt;
return value;
}
template<class CloudType>
Foam::scalar Foam::PressureGradientForce<CloudType>::massAdd
(
const typename CloudType::parcelType&,
const scalar
) const
{
return 0.0;
}
// ************************************************************************* //
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -38,6 +38,7 @@ SourceFiles
#include "ParticleForce.H"
#include "volFields.H"
#include "interpolation.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -53,13 +54,15 @@ class PressureGradientForce
:
public ParticleForce<CloudType>
{
// Private data
protected:
// Protected data
//- Name of velocity field
const word UName_;
//- Velocity gradient field
const volTensorField* gradUPtr_;
//- Rate of change of carrier phase velocity interpolator
autoPtr<interpolation<vector> > DUcDtInterpPtr_;
public:
......@@ -75,7 +78,8 @@ public:
(
CloudType& owner,
const fvMesh& mesh,
const dictionary& dict
const dictionary& dict,
const word& forceType = typeName
);
//- Construct copy
......@@ -99,8 +103,8 @@ public:
// Access
//- Return const access to the velocity gradient field
inline const volTensorField& gradU() const;
//- Return the rate of change of carrier phase velocity interpolator
inline const interpolation<vector>& DUcDtInterp() const;
// Evaluation
......@@ -117,6 +121,13 @@ public:
const scalar Re,
const scalar muc
) const;
//- Return the added mass
virtual scalar massAdd
(
const typename CloudType::parcelType& p,
const scalar mass
) const;
};
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -26,23 +26,20 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class CloudType>
const Foam::volTensorField& Foam::PressureGradientForce<CloudType>::gradU()
const
inline const Foam::interpolation<Foam::vector>&
Foam::PressureGradientForce<CloudType>::DUcDtInterp() const
{
if (gradUPtr_)
{
return *gradUPtr_;
}
else
if (!DUcDtInterpPtr_.valid())
{
FatalErrorIn
(
"const volTensorField& PressureGradientForce<CloudType>::gradU()"
"const"
) << "gradU field not allocated" << abort(FatalError);
return *reinterpret_cast<const volTensorField*>(0);
"inline const Foam::interpolation<Foam::vector>&"
"Foam::PressureGradientForce<CloudType>::DUcDtInterp() const"
) << "Carrier phase DUcDt interpolation object not set"
<< abort(FatalError);
}
return DUcDtInterpPtr_();
}
......
......@@ -24,8 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "VirtualMassForce.H"
#include "fvcDdt.H"
#include "fvcGrad.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -34,14 +32,12 @@ Foam::VirtualMassForce<CloudType>::VirtualMassForce
(
CloudType& owner,
const fvMesh& mesh,
const dictionary& dict
const dictionary& dict,
const word& forceType
)
:
ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
Cvm_(readScalar(this->coeffs().lookup("Cvm"))),
DUcDtPtr_(NULL),
DUcDtInterpPtr_(NULL)
PressureGradientForce<CloudType>(owner, mesh, dict, forceType),
Cvm_(readScalar(this->coeffs().lookup("Cvm")))
{}
......@@ -51,11 +47,8 @@ Foam::VirtualMassForce<CloudType>::VirtualMassForce
const VirtualMassForce& vmf
)
:
ParticleForce<CloudType>(vmf),
UName_(vmf.UName_),
Cvm_(vmf.Cvm_),
DUcDtPtr_(NULL),
DUcDtInterpPtr_(NULL)
PressureGradientForce<CloudType>(vmf),
Cvm_(vmf.Cvm_)
{}
......@@ -71,36 +64,7 @@ Foam::VirtualMassForce<CloudType>::~VirtualMassForce()
template<class CloudType>
void Foam::VirtualMassForce<CloudType>::cacheFields(const bool store)
{
if (store && !DUcDtPtr_)
{
const volVectorField& Uc = this->mesh().template
lookupObject<volVectorField>(UName_);
DUcDtPtr_ = new volVectorField
(
"DUcDt",
fvc::ddt(Uc) + (Uc & fvc::grad(Uc))
);
DUcDtInterpPtr_.reset
(
interpolation<vector>::New
(
this->owner().solution().interpolationSchemes(),
*DUcDtPtr_
).ptr()
);
}
else
{
DUcDtInterpPtr_.clear();
if (DUcDtPtr_)
{
delete DUcDtPtr_;
DUcDtPtr_ = NULL;
}
}
PressureGradientForce<CloudType>::cacheFields(store);
}
......@@ -114,12 +78,10 @@ Foam::forceSuSp Foam::VirtualMassForce<CloudType>::calcCoupled
const scalar muc
) const
{
forceSuSp value(vector::zero, 0.0);
forceSuSp value =
PressureGradientForce<CloudType>::calcCoupled(p, dt, mass, Re, muc);
vector DUcDt =
DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
value.Su() = mass*p.rhoc()/p.rho()*Cvm_*DUcDt;
value.Su() *= Cvm_;
return value;
}
......
......@@ -28,7 +28,6 @@ Description
Calculates particle virtual mass force
SourceFiles
VirtualMassForceI.H
VirtualMassForce.C
\*---------------------------------------------------------------------------*/
......@@ -36,9 +35,7 @@ SourceFiles
#ifndef VirtualMassForce_H
#define VirtualMassForce_H
#include "ParticleForce.H"
#include "volFields.H"
#include "interpolation.H"
#include "PressureGradientForce.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -52,22 +49,13 @@ namespace Foam
template<class CloudType>
class VirtualMassForce
:
public ParticleForce<CloudType>
public PressureGradientForce<CloudType>
{
// Private data
//- Name of velocity field
const word UName_;
//- Virtual mass coefficient - typically 0.5
scalar Cvm_;
//- Rate of change of carrier phase velocity
volVectorField* DUcDtPtr_;
//- Rate of change of carrier phase velocity interpolator
autoPtr<interpolation<vector> > DUcDtInterpPtr_;
public:
......@@ -82,7 +70,8 @@ public:
(
CloudType& owner,
const fvMesh& mesh,
const dictionary& dict
const dictionary& dict,
const word& forceType = typeName
);
//- Construct copy
......@@ -104,15 +93,6 @@ public:
// Member Functions
// Access
//- Return the rate of change of carrier phase velocity
inline const volVectorField& DUcDt() const;
//- Return the rate of change of carrier phase velocity interpolator
inline const interpolation<vector>& DUcDtInterp() const;
// Evaluation
//- Cache fields
......@@ -143,8 +123,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "VirtualMassForceI.H"
#ifdef NoRepository
#include "VirtualMassForce.C"
#endif
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment