Merge branch 'feature-snappyHexMesh' into feature-snappyHexMesh-gapRefinement

This commit is contained in:
mattijs 2015-11-03 12:55:32 +00:00
commit 4927f97c5d
45 changed files with 1853 additions and 342 deletions

View File

@ -45,6 +45,8 @@ Description
int main(int argc, char *argv[])
{
argList::noParallel();
#include "setRootCase.H"
#include "createTime.H"

View File

@ -128,7 +128,7 @@ int main(int argc, char *argv[])
}
}
#include "write.H"
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()

View File

@ -1,17 +0,0 @@
if (runTime.outputTime())
{
volVectorField Ur
(
IOobject
(
"Ur",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
U1 - U2
);
runTime.write();
}

View File

@ -493,12 +493,15 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
// Min length
scalar minDistSqr = magSqr(1e-6 * globalBb.span());
// Non-empty directions
// Geometric directions
const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
Info<< " Mesh (non-empty, non-wedge) directions " << validDirs << endl;
Info<< " Mesh has " << mesh.nGeometricD()
<< " geometric (non-empty/wedge) directions " << validDirs << endl;
// Solution directions
const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
Info<< " Mesh (non-empty) directions " << solDirs << endl;
Info<< " Mesh has " << mesh.nSolutionD()
<< " solution (non-empty) directions " << solDirs << endl;
if (mesh.nGeometricD() < 3)
{

View File

@ -133,7 +133,11 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << endl;
mesh.update();
for (int i = 0; i<2; i++)
{
mesh.update();
}
mesh.checkMesh(true);
if (checkAMI)

View File

@ -22,7 +22,7 @@ BirdCarreauCoeffs
{
nu0 [0 2 -1 0 0 0 0] 1e-03;
nuInf [0 2 -1 0 0 0 0] 1e-05;
m [0 0 1 0 0 0 0] 1;
k [0 0 1 0 0 0 0] 1;
n [0 0 0 0 0 0 0] 0.5;
}

View File

@ -22,7 +22,7 @@ BirdCarreauCoeffs
{
nu0 [0 2 -1 0 0 0 0] 1e-03;
nuInf [0 2 -1 0 0 0 0] 1e-05;
m [0 0 1 0 0 0 0] 1;
k [0 0 1 0 0 0 0] 1;
n [0 0 0 0 0 0 0] 0.5;
}

View File

@ -22,7 +22,7 @@ BirdCarreauCoeffs
{
nu0 [0 2 -1 0 0 0 0] 1e-03;
nuInf [0 2 -1 0 0 0 0] 1e-05;
m [0 0 1 0 0 0 0] 1;
k [0 0 1 0 0 0 0] 1;
n [0 0 0 0 0 0 0] 0.5;
}

View File

@ -22,7 +22,7 @@ BirdCarreauCoeffs
{
nu0 [0 2 -1 0 0 0 0] 1e-03;
nuInf [0 2 -1 0 0 0 0] 1e-05;
m [0 0 1 0 0 0 0] 1;
k [0 0 1 0 0 0 0] 1;
n [0 0 0 0 0 0 0] 0.5;
}

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -2748,7 +2748,13 @@ void Foam::globalMeshData::updateMesh()
// *** Temporary hack to avoid problems with overlapping communication
// *** between these reductions and the calculation of deltaCoeffs
label comm = UPstream::worldComm + 1;
//label comm = UPstream::worldComm + 1;
label comm = UPstream::allocateCommunicator
(
UPstream::worldComm,
identity(UPstream::nProcs(UPstream::worldComm)),
true
);
// Total number of faces.
nTotalFaces_ = returnReduce
@ -2785,6 +2791,8 @@ void Foam::globalMeshData::updateMesh()
comm
);
UPstream::freeCommunicator(comm);
if (debug)
{
Pout<< "globalMeshData : nTotalPoints_:" << nTotalPoints_ << endl;

View File

@ -1483,11 +1483,7 @@ Foam::label Foam::polyMesh::findCell
return -1;
}
if
(
Pstream::parRun()
&& (decompMode == FACE_DIAG_TRIS || decompMode == CELL_TETS)
)
if (Pstream::parRun() && decompMode == FACE_DIAG_TRIS)
{
// Force construction of face-diagonal decomposition before testing
// for zero cells.

View File

@ -199,6 +199,7 @@ $(derivedFvPatchFields)/waveTransmissive/waveTransmissiveFvPatchFields.C
$(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
$(derivedFvPatchFields)/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/prghPressure/prghPressureFvPatchScalarField.C
$(derivedFvPatchFields)/prghTotalPressure/prghTotalPressureFvPatchScalarField.C
fvsPatchFields = fields/fvsPatchFields
$(fvsPatchFields)/fvsPatchField/fvsPatchFields.C

View File

@ -0,0 +1,213 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "prghTotalPressureFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "uniformDimensionedFields.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::prghTotalPressureFvPatchScalarField::
prghTotalPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(p, iF),
UName_("U"),
phiName_("phi"),
rhoName_("rho"),
p0_(p.size(), 0.0)
{}
Foam::prghTotalPressureFvPatchScalarField::
prghTotalPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchScalarField(p, iF),
UName_(dict.lookupOrDefault<word>("U", "U")),
phiName_(dict.lookupOrDefault<word>("phi", "phi")),
rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
p0_("p0", dict, p.size())
{
if (dict.found("value"))
{
fvPatchScalarField::operator=
(
scalarField("value", dict, p.size())
);
}
else
{
fvPatchField<scalar>::operator=(p0_);
}
}
Foam::prghTotalPressureFvPatchScalarField::
prghTotalPressureFvPatchScalarField
(
const prghTotalPressureFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchScalarField(ptf, p, iF, mapper),
UName_(ptf.UName_),
phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_),
p0_(ptf.p0_, mapper)
{}
Foam::prghTotalPressureFvPatchScalarField::
prghTotalPressureFvPatchScalarField
(
const prghTotalPressureFvPatchScalarField& ptf
)
:
fixedValueFvPatchScalarField(ptf),
UName_(ptf.UName_),
phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_),
p0_(ptf.p0_)
{}
Foam::prghTotalPressureFvPatchScalarField::
prghTotalPressureFvPatchScalarField
(
const prghTotalPressureFvPatchScalarField& ptf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(ptf, iF),
UName_(ptf.UName_),
phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_),
p0_(ptf.p0_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::prghTotalPressureFvPatchScalarField::autoMap
(
const fvPatchFieldMapper& m
)
{
fixedValueFvPatchScalarField::autoMap(m);
p0_.autoMap(m);
}
void Foam::prghTotalPressureFvPatchScalarField::rmap
(
const fvPatchScalarField& ptf,
const labelList& addr
)
{
fixedValueFvPatchScalarField::rmap(ptf, addr);
const prghTotalPressureFvPatchScalarField& tiptf =
refCast<const prghTotalPressureFvPatchScalarField>(ptf);
p0_.rmap(tiptf.p0_, addr);
}
void Foam::prghTotalPressureFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const scalarField& rhop =
patch().lookupPatchField<volScalarField, scalar>(rhoName_);
const scalarField& phip =
patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
const vectorField& Up =
patch().lookupPatchField<volVectorField, vector>(UName_);
const uniformDimensionedVectorField& g =
db().lookupObject<uniformDimensionedVectorField>("g");
const uniformDimensionedScalarField& hRef =
db().lookupObject<uniformDimensionedScalarField>("hRef");
dimensionedScalar ghRef
(
mag(g.value()) > SMALL
? g & (cmptMag(g.value())/mag(g.value()))*hRef
: dimensionedScalar("ghRef", g.dimensions()*dimLength, 0)
);
operator==
(
p0_
- 0.5*rhop*(1.0 - pos(phip))*magSqr(Up)
- rhop*((g.value() & patch().Cf()) - ghRef.value())
);
fixedValueFvPatchScalarField::updateCoeffs();
}
void Foam::prghTotalPressureFvPatchScalarField::write(Ostream& os) const
{
fvPatchScalarField::write(os);
writeEntryIfDifferent<word>(os, "U", "U", UName_);
writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
p0_.writeEntry("p0", os);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchScalarField,
prghTotalPressureFvPatchScalarField
);
}
// ************************************************************************* //

View File

@ -0,0 +1,243 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::prghTotalPressureFvPatchScalarField
Group
grpGenericBoundaryConditions
Description
This boundary condition provides static pressure condition for p_rgh,
calculated as:
\f[
p_rgh = p - \rho g (h - hRef)
\f]
where
\vartable
p_rgh | Pseudo hydrostatic pressure [Pa]
p | Static pressure [Pa]
h | Height in the opposite direction to gravity
hRef | Reference height in the opposite direction to gravity
\rho | density
g | acceleration due to gravity [m/s^2]
\endtable
\heading Patch usage
\table
Property | Description | Required | Default value
rhoName | rho field name | no | rho
p | static pressure | yes |
\endtable
Example of the boundary condition specification:
\verbatim
myPatch
{
type prghTotalPressure;
rhoName rho;
p uniform 0;
value uniform 0; // optional initial value
}
\endverbatim
SeeAlso
Foam::fixedValueFvPatchScalarField
SourceFiles
prghTotalPressureFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef prghTotalPressureFvPatchScalarField_H
#define prghTotalPressureFvPatchScalarField_H
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class prghTotalPressureFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class prghTotalPressureFvPatchScalarField
:
public fixedValueFvPatchScalarField
{
protected:
// Protected data
//- Name of the velocity field
word UName_;
//- Name of the flux transporting the field
word phiName_;
//- Name of phase-fraction field
word rhoName_;
//- Total pressure
scalarField p0_;
public:
//- Runtime type information
TypeName("prghTotalPressure");
// Constructors
//- Construct from patch and internal field
prghTotalPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
prghTotalPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given
// prghTotalPressureFvPatchScalarField onto a new patch
prghTotalPressureFvPatchScalarField
(
const prghTotalPressureFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
prghTotalPressureFvPatchScalarField
(
const prghTotalPressureFvPatchScalarField&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField >
(
new prghTotalPressureFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
prghTotalPressureFvPatchScalarField
(
const prghTotalPressureFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchScalarField> clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchScalarField>
(
new prghTotalPressureFvPatchScalarField(*this, iF)
);
}
// Member functions
// Access
//- Return the rhoName
const word& rhoName() const
{
return rhoName_;
}
//- Return reference to the rhoName to allow adjustment
word& rhoName()
{
return rhoName_;
}
//- Return the total pressure
const scalarField& p0() const
{
return p0_;
}
//- Return reference to the total pressure to allow adjustment
scalarField& p0()
{
return p0_;
}
// Mapping functions
//- Map (and resize as needed) from self given a mapping object
virtual void autoMap
(
const fvPatchFieldMapper&
);
//- Reverse map the given fvPatchField onto this fvPatchField
virtual void rmap
(
const fvPatchScalarField&,
const labelList&
);
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -81,15 +81,8 @@ Foam::uniformFixedValueFvPatchField<Type>::uniformFixedValueFvPatchField
fixedValueFvPatchField<Type>(p, iF),
uniformValue_(DataEntry<Type>::New("uniformValue", dict))
{
if (dict.found("value"))
{
fvPatchField<Type>::operator==(Field<Type>("value", dict, p.size()));
}
else
{
const scalar t = this->db().time().timeOutputValue();
fvPatchField<Type>::operator==(uniformValue_->value(t));
}
const scalar t = this->db().time().timeOutputValue();
fvPatchField<Type>::operator==(uniformValue_->value(t));
}

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -141,6 +141,10 @@ void Foam::uniformInletOutletFvPatchField<Type>::updateCoeffs()
return;
}
// Update the uniform value as a function of time
const scalar t = this->db().time().timeOutputValue();
this->refValue() = uniformInletValue_->value(t);
const Field<scalar>& phip =
this->patch().template lookupPatchField<surfaceScalarField, scalar>
(

View File

@ -501,7 +501,7 @@ void Foam::fvSchemes::setFluxRequired(const word& name) const
Info<< "Setting fluxRequired for " << name << endl;
}
fluxRequired_.add(name, true);
fluxRequired_.add(name, true, true);
}

View File

@ -35,27 +35,31 @@ namespace Foam
defineTypeNameAndDebug(turbulenceFields, 0);
template<>
const char* NamedEnum<turbulenceFields::compressibleField, 6>::names[] =
const char* NamedEnum<turbulenceFields::compressibleField, 8>::names[] =
{
"R",
"devRhoReff",
"k",
"epsilon",
"mut",
"muEff",
"alphat",
"alphaEff"
"alphaEff",
"R",
"devRhoReff"
};
const NamedEnum<turbulenceFields::compressibleField, 6>
const NamedEnum<turbulenceFields::compressibleField, 8>
turbulenceFields::compressibleFieldNames_;
template<>
const char* NamedEnum<turbulenceFields::incompressibleField, 4>::names[] =
const char* NamedEnum<turbulenceFields::incompressibleField, 6>::names[] =
{
"R",
"devReff",
"k",
"epsilon",
"nut",
"nuEff"
"nuEff",
"R",
"devReff"
};
const NamedEnum<turbulenceFields::incompressibleField, 4>
const NamedEnum<turbulenceFields::incompressibleField, 6>
turbulenceFields::incompressibleFieldNames_;
const word turbulenceFields::modelName = turbulenceModel::propertiesName;
@ -174,14 +178,14 @@ void Foam::turbulenceFields::execute()
const word& f = iter.key();
switch (compressibleFieldNames_[f])
{
case cfR:
case cfK:
{
processField<symmTensor>(f, model.R());
processField<scalar>(f, model.k());
break;
}
case cfDevRhoReff:
case cfEpsilon:
{
processField<symmTensor>(f, model.devRhoReff());
processField<scalar>(f, model.epsilon());
break;
}
case cfMut:
@ -204,6 +208,16 @@ void Foam::turbulenceFields::execute()
processField<scalar>(f, model.alphaEff());
break;
}
case cfR:
{
processField<symmTensor>(f, model.R());
break;
}
case cfDevRhoReff:
{
processField<symmTensor>(f, model.devRhoReff());
break;
}
default:
{
FatalErrorIn("void Foam::turbulenceFields::execute()")
@ -222,14 +236,14 @@ void Foam::turbulenceFields::execute()
const word& f = iter.key();
switch (incompressibleFieldNames_[f])
{
case ifR:
case ifK:
{
processField<symmTensor>(f, model.R());
processField<scalar>(f, model.k());
break;
}
case ifDevReff:
case ifEpsilon:
{
processField<symmTensor>(f, model.devReff());
processField<scalar>(f, model.epsilon());
break;
}
case ifNut:
@ -242,6 +256,16 @@ void Foam::turbulenceFields::execute()
processField<scalar>(f, model.nuEff());
break;
}
case ifR:
{
processField<symmTensor>(f, model.R());
break;
}
case ifDevReff:
{
processField<symmTensor>(f, model.devReff());
break;
}
default:
{
FatalErrorIn("void Foam::turbulenceFields::execute()")
@ -263,15 +287,11 @@ void Foam::turbulenceFields::end()
void Foam::turbulenceFields::timeSet()
{
// Do nothing
}
{}
void Foam::turbulenceFields::write()
{
// Do nothing
}
{}
// ************************************************************************* //

View File

@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -53,22 +53,24 @@ Description
\heading Function object usage
\table
Property | Description | Required | Default value
type | type name: processorField | yes |
fields | fields to store (see below) | yes |
Property | Description | Required | Default value
type | type name: turbulenceFields | yes |
fields | fields to store (see below) | yes |
\endtable
Where \c fields can include:
\plaintable
R | Stress tensor
devRhoReff |
k | turbulence kinetic energy
epsilon | turbulence kinetic energy dissipation rate
nut | turbulence viscosity (incompressible)
nuEff | effective turbulence viscosity (incompressible)
mut | turbulence viscosity (compressible)
muEff | effective turbulence viscosity (compressible)
alphat | turbulence thermal diffusivity (compressible)
alphaEff | effective turbulence thermal diffusivity (compressible)
devReff |
nut | turbulence viscosity (incompressible)
nuEff | effective turbulence viscosity (incompressible)
R | Reynolds stress tensor
devReff | Deviatoric part of the effective Reynolds stress
devRhoReff | Divergence of the Reynolds stress
\endplaintable
SeeAlso
@ -109,23 +111,27 @@ public:
enum compressibleField
{
cfR,
cfDevRhoReff,
cfK,
cfEpsilon,
cfMut,
cfMuEff,
cfAlphat,
cfAlphaEff
cfAlphaEff,
cfR,
cfDevRhoReff
};
static const NamedEnum<compressibleField, 6> compressibleFieldNames_;
static const NamedEnum<compressibleField, 8> compressibleFieldNames_;
enum incompressibleField
{
ifR,
ifDevReff,
ifK,
ifEpsilon,
ifNut,
ifNuEff
ifNuEff,
ifR,
ifDevReff
};
static const NamedEnum<incompressibleField, 4> incompressibleFieldNames_;
static const NamedEnum<incompressibleField, 6> incompressibleFieldNames_;
static const word modelName;

View File

@ -31,4 +31,10 @@ sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
sixDoFRigidBodyMotionSolver/externalPointEdgePoint.C
sixDoFRigidBodyMotionSolver/pointPatchDist.C
sixDoFSolvers/sixDoFSolver/sixDoFSolver.C
sixDoFSolvers/sixDoFSolver/newSixDoFSolver.C
sixDoFSolvers/symplectic/symplectic.C
sixDoFSolvers/CrankNicolson/CrankNicolson.C
sixDoFSolvers/Newmark/Newmark.C
LIB = $(FOAM_LIBBIN)/libsixDoFRigidBodyMotion

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -204,18 +204,14 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
const pointPatch& ptPatch = this->patch();
// Store the motion state at the beginning of the time-step
bool firstIter = false;
if (curTimeIndex_ != t.timeIndex())
{
motion_.newTime();
curTimeIndex_ = t.timeIndex();
firstIter = true;
}
// Patch force data is valid for the current positions, so
// calculate the forces on the motion object from this data, then
// update the positions
motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
dictionary forcesDict;
forcesDict.add("type", forces::typeName);
@ -241,11 +237,13 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
// scalar ramp = min(max((t.value() - 5)/10, 0), 1);
scalar ramp = 1.0;
motion_.updateAcceleration
motion_.update
(
firstIter,
ramp*(f.forceEff() + motion_.mass()*g_),
ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g_)),
t.deltaTValue()
t.deltaTValue(),
t.deltaT0Value()
);
Field<vector>::operator=

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -152,14 +152,14 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
const Time& t = mesh.time();
// Store the motion state at the beginning of the time-step
bool firstIter = false;
if (curTimeIndex_ != t.timeIndex())
{
motion_.newTime();
curTimeIndex_ = t.timeIndex();
firstIter = true;
}
motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
vector gravity = vector::zero;
if (db().foundObject<uniformDimensionedVectorField>("g"))
@ -171,11 +171,13 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
}
// Do not modify the accelerations, except with gravity, where available
motion_.updateAcceleration
motion_.update
(
firstIter,
gravity*motion_.mass(),
vector::zero,
t.deltaTValue()
t.deltaTValue(),
t.deltaT0Value()
);
Field<vector>::operator=

View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "sixDoFRigidBodyMotion.H"
#include "sixDoFSolver.H"
#include "septernion.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -84,7 +85,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
momentOfInertia_(diagTensor::one*VSMALL),
aRelax_(1.0),
aDamp_(1.0),
report_(false)
report_(false),
solver_(NULL)
{}
@ -121,7 +123,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
momentOfInertia_(dict.lookup("momentOfInertia")),
aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
aDamp_(dict.lookupOrDefault<scalar>("accelerationDamping", 1.0)),
report_(dict.lookupOrDefault<Switch>("report", false))
report_(dict.lookupOrDefault<Switch>("report", false)),
solver_(sixDoFSolver::New(dict.subDict("solver"), *this))
{
addRestraints(dict);
@ -171,6 +174,12 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::sixDoFRigidBodyMotion::addRestraints
@ -257,67 +266,46 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
}
void Foam::sixDoFRigidBodyMotion::updatePosition
(
scalar deltaT,
scalar deltaT0
)
{
// First leapfrog velocity adjust and motion part, required before
// force calculation
if (Pstream::master())
{
v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
// Leapfrog move part
centreOfRotation() = centreOfRotation0() + deltaT*v();
// Leapfrog orientation adjustment
Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
Q() = Qpi.first();
pi() = rConstraints_ & Qpi.second();
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
const vector& tauGlobal
)
{
static bool first = false;
// Second leapfrog velocity adjust part, required after motion and
// acceleration calculation
// Save the previous iteration accelerations for relaxation
vector aPrevIter = a();
vector tauPrevIter = tau();
// Calculate new accelerations
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();
// Relax accelerations on all but first iteration
if (!first)
{
a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
}
first = false;
}
void Foam::sixDoFRigidBodyMotion::update
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
)
{
if (Pstream::master())
{
// Save the previous iteration accelerations for relaxation
vector aPrevIter = a();
vector tauPrevIter = tau();
// Calculate new accelerations
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();
// Relax accelerations on all but first iteration
if (!first)
{
a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
}
first = false;
// Correct velocities
v() += tConstraints_ & aDamp_*0.5*deltaT*a();
pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
solver_->solve(firstIter, fGlobal, tauGlobal, deltaT, deltaT0);
if (report_)
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,25 +25,16 @@ Class
Foam::sixDoFRigidBodyMotion
Description
Six degree of freedom motion for a rigid body. Angular momentum
stored in body fixed reference frame. Reference orientation of
the body (where Q = I) must align with the cartesian axes such
that the Inertia tensor is in principle component form.
Six degree of freedom motion for a rigid body.
Symplectic motion as per:
Angular momentum stored in body fixed reference frame. Reference
orientation of the body (where Q = I) must align with the cartesian axes
such that the Inertia tensor is in principle component form. Can add
restraints (e.g. a spring) and constraints (e.g. motion may only be on a
plane).
title = {Symplectic splitting methods for rigid body molecular dynamics},
publisher = {AIP},
year = {1997},
journal = {The Journal of Chemical Physics},
volume = {107},
number = {15},
pages = {5840-5851},
url = {http://link.aip.org/link/?JCP/107/5840/1},
doi = {10.1063/1.474310}
Can add restraints (e.g. a spring)
and constraints (e.g. motion may only be on a plane).
The time-integrator for the motion is run-time selectable with options for
symplectic (explicit), Crank-Nicolson and Newmark schemes.
SourceFiles
sixDoFRigidBodyMotionI.H
@ -66,12 +57,17 @@ SourceFiles
namespace Foam
{
// Forward declarations
class sixDoFSolver;
/*---------------------------------------------------------------------------*\
Class sixDoFRigidBodyMotion Declaration
\*---------------------------------------------------------------------------*/
class sixDoFRigidBodyMotion
{
friend class sixDoFSolver;
// Private data
//- Motion state data object
@ -117,6 +113,9 @@ class sixDoFRigidBodyMotion
//- Switch to turn reporting of motion data on and off
Switch report_;
//- Motion solver
autoPtr<sixDoFSolver> solver_;
// Private Member Functions
@ -144,6 +143,9 @@ class sixDoFRigidBodyMotion
//- Apply the restraints to the object
void applyRestraints();
//- Update and relax accelerations from the force and torque
void updateAcceleration(const vector& fGlobal, const vector& tauGlobal);
// Access functions retained as private because of the risk of
// confusion over what is a body local frame vector and what is global
@ -176,21 +178,6 @@ class sixDoFRigidBodyMotion
//- Return the current torque
inline const vector& tau() const;
//- Return the orientation at previous time-step
inline const tensor& Q0() const;
//- Return the velocity at previous time-step
inline const vector& v0() const;
//- Return the acceleration at previous time-step
inline const vector& a0() const;
//- Return the angular momentum at previous time-step
inline const vector& pi0() const;
//- Return the torque at previous time-step
inline const vector& tau0() const;
// Edit
@ -234,6 +221,10 @@ public:
sixDoFRigidBodyMotion(const sixDoFRigidBodyMotion&);
//- Destructor
~sixDoFRigidBodyMotion();
// Member Functions
// Access
@ -250,9 +241,6 @@ public:
//- Return the current centre of rotation
inline const point& centreOfRotation() const;
//- Return the centre of rotation at previous time-step
inline const point& centreOfRotation0() const;
//- Return the initial centre of mass
inline const point& initialCentreOfMass() const;
@ -298,18 +286,15 @@ public:
// Update state
//- First leapfrog velocity adjust and motion part, required
// before force calculation. Takes old timestep for variable
// timestep cases.
void updatePosition(scalar deltaT, scalar deltaT0);
//- Second leapfrog velocity adjust part
// required after motion and force calculation
void updateAcceleration
//- Symplectic integration of velocities, orientation and position.
// Changes to Crank-Nicolson integration for subsequent iterations.
void update
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
scalar deltaT,
scalar deltaT0
);
//- Report the status of the motion

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -161,36 +161,6 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
}
inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q0() const
{
return motionState0_.Q();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::v0() const
{
return motionState0_.v();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::a0() const
{
return motionState0_.a();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi0() const
{
return motionState0_.pi();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau0() const
{
return motionState0_.tau();
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::initialCentreOfRotation()
{
return initialCentreOfRotation_;
@ -261,12 +231,6 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfRotation() const
}
inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfRotation0() const
{
return motionState0_.centreOfRotation();
}
inline const Foam::point&
Foam::sixDoFRigidBodyMotion::initialCentreOfMass() const
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -85,6 +85,7 @@ Foam::sixDoFRigidBodyMotionSolver::sixDoFRigidBodyMotionSolver
patchSet_(mesh.boundaryMesh().patchSet(patches_)),
di_(readScalar(coeffDict().lookup("innerDistance"))),
do_(readScalar(coeffDict().lookup("outerDistance"))),
test_(coeffDict().lookupOrDefault<Switch>("test", false)),
rhoInf_(1.0),
rhoName_(coeffDict().lookupOrDefault<word>("rhoName", "rho")),
scale_
@ -179,31 +180,15 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
<< " points." << exit(FatalError);
}
// Store the motion state at the beginning of the time-step
// Store the motion state at the beginning of the time-stepbool
bool firstIter = false;
if (curTimeIndex_ != this->db().time().timeIndex())
{
motion_.newTime();
curTimeIndex_ = this->db().time().timeIndex();
firstIter = true;
}
// Patch force data is valid for the current positions, so
// calculate the forces on the motion object from this data, then
// update the positions
motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
dictionary forcesDict;
forcesDict.add("type", forces::typeName);
forcesDict.add("patches", patches_);
forcesDict.add("rhoInf", rhoInf_);
forcesDict.add("rhoName", rhoName_);
forcesDict.add("CofR", motion_.centreOfRotation());
forces f("forces", db(), forcesDict);
f.calcForcesMoment();
dimensionedVector g("g", dimAcceleration, vector::zero);
if (db().foundObject<uniformDimensionedVectorField>("g"))
@ -218,12 +203,44 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
// scalar ramp = min(max((this->db().time().value() - 5)/10, 0), 1);
scalar ramp = 1.0;
motion_.updateAcceleration
(
ramp*(f.forceEff() + motion_.mass()*g.value()),
ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g.value())),
t.deltaTValue()
);
if (test_)
{
motion_.update
(
firstIter,
ramp*(motion_.mass()*g.value()),
ramp*(motion_.mass()*(motion_.momentArm() ^ g.value())),
t.deltaTValue(),
t.deltaT0Value()
);
}
else
{
dictionary forcesDict;
forcesDict.add("type", forces::typeName);
forcesDict.add("patches", patches_);
forcesDict.add("rhoInf", rhoInf_);
forcesDict.add("rhoName", rhoName_);
forcesDict.add("CofR", motion_.centreOfRotation());
forces f("forces", db(), forcesDict);
f.calcForcesMoment();
motion_.update
(
firstIter,
ramp*(f.forceEff() + motion_.mass()*g.value()),
ramp
*(
f.momentEff()
+ motion_.mass()*(motion_.momentArm() ^ g.value())
),
t.deltaTValue(),
t.deltaT0Value()
);
}
// Update the displacements
pointDisplacement_.internalField() =

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -70,6 +70,10 @@ class sixDoFRigidBodyMotionSolver
//- Outer morphing distance (limit of linear interpolation region)
const scalar do_;
//- Switch for test-mode in which only the
// gravitational body-force is applied
Switch test_;
//- Reference density required by the forces object for
// incompressible calculations, required if rhoName == rhoInf
scalar rhoInf_;

View File

@ -0,0 +1,94 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "CrankNicolson.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFSolvers
{
defineTypeNameAndDebug(CrankNicolson, 0);
addToRunTimeSelectionTable(sixDoFSolver, CrankNicolson, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::CrankNicolson::CrankNicolson
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
)
:
sixDoFSolver(body),
aoc_(dict.lookupOrDefault<scalar>("aoc", 0.5)),
voc_(dict.lookupOrDefault<scalar>("voc", 0.5))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::CrankNicolson::~CrankNicolson()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::sixDoFSolvers::CrankNicolson::solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
)
{
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);
// Correct linear velocity
v() = tConstraints()
& (v0() + aDamp()*deltaT*(aoc_*a() + (1 - aoc_)*a0()));
// Correct angular momentum
pi() = rConstraints()
& (pi0() + aDamp()*deltaT*(aoc_*tau() + (1 - aoc_)*tau0()));
// Correct position
centreOfRotation() =
centreOfRotation0() + deltaT*(voc_*v() + (1 - voc_)*v0());
// Correct orientation
Tuple2<tensor, vector> Qpi =
rotate(Q0(), (voc_*pi() + (1 - voc_)*pi0()), deltaT);
Q() = Qpi.first();
}
// ************************************************************************* //

View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::sixDoFSolvers::CrankNicolson
Description
Crank-Nicolson 2nd-order time-integrator for 6DoF solid-body motion.
The off-centering coefficients for acceleration (velocity integration) and
velocity (position/orientation integration) may be specified but default
values of 0.5 for each are used if they are not specified. With the default
off-centering this scheme is equivalent to the Newmark scheme with default
coefficients.
Example specification in dynamicMeshDict:
\verbatim
solver
{
type CrankNicolson;
aoc 0.5; // Acceleration off-centering coefficient
voc 0.5; // Velocity off-centering coefficient
}
\endverbatim
SeeAlso
Foam::sixDoFSolvers::Newmark
SourceFiles
CrankNicolson.C
\*---------------------------------------------------------------------------*/
#ifndef CrankNicolson_H
#define CrankNicolson_H
#include "sixDoFSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFSolvers
{
/*---------------------------------------------------------------------------*\
Class CrankNicolson Declaration
\*---------------------------------------------------------------------------*/
class CrankNicolson
:
public sixDoFSolver
{
// Private data
//- Acceleration off-centering coefficient (default: 0.5)
const scalar aoc_;
//- Velocity off-centering coefficient (default: 0.5)
const scalar voc_;
public:
//- Runtime type information
TypeName("CrankNicolson");
// Constructors
//- Construct from a dictionary and the body
CrankNicolson
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
);
//- Destructor
virtual ~CrankNicolson();
// Member Functions
//- Drag coefficient
virtual void solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace sixDoFSolvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,115 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Newmark.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFSolvers
{
defineTypeNameAndDebug(Newmark, 0);
addToRunTimeSelectionTable(sixDoFSolver, Newmark, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::Newmark::Newmark
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
)
:
sixDoFSolver(body),
gamma_(dict.lookupOrDefault<scalar>("gamma", 0.5)),
beta_
(
max
(
0.25*sqr(gamma_ + 0.5),
dict.lookupOrDefault<scalar>("beta", 0.25)
)
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::Newmark::~Newmark()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::sixDoFSolvers::Newmark::solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
)
{
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);
// Correct linear velocity
v() =
tConstraints()
& (v0() + aDamp()*deltaT*(gamma_*a() + (1 - gamma_)*a0()));
// Correct angular momentum
pi() =
rConstraints()
& (pi0() + aDamp()*deltaT*(gamma_*tau() + (1 - gamma_)*tau0()));
// Correct position
centreOfRotation() =
centreOfRotation0()
+ (
tConstraints()
& (
deltaT*v0()
+ sqr(deltaT)*beta_*a() + sqr(deltaT)*(0.5 - beta_)*a0()
)
);
// Correct orientation
vector piDeltaT =
rConstraints()
& (
deltaT*pi0()
+ sqr(deltaT)*beta_*tau() + sqr(deltaT)*(0.5 - beta_)*tau0()
);
Tuple2<tensor, vector> Qpi = rotate(Q0(), piDeltaT, 1);
Q() = Qpi.first();
}
// ************************************************************************* //

View File

@ -0,0 +1,124 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::sixDoFSolvers::Newmark
Description
Newmark 2nd-order time-integrator for 6DoF solid-body motion.
Reference:
\verbatim
Newmark, N. M. (1959).
A method of computation for structural dynamics.
Journal of the Engineering Mechanics Division, 85(3), 67-94.
\endverbatim
Example specification in dynamicMeshDict:
\verbatim
solver
{
type Newmark;
gamma 0.5; // Velocity integration coefficient
beta 0.25; // Position integration coefficient
}
\endverbatim
SourceFiles
Newmark.C
\*---------------------------------------------------------------------------*/
#ifndef Newmark_H
#define Newmark_H
#include "sixDoFSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFSolvers
{
/*---------------------------------------------------------------------------*\
Class Newmark Declaration
\*---------------------------------------------------------------------------*/
class Newmark
:
public sixDoFSolver
{
// Private data
//- Coefficient for velocity integration (default: 0.5)
const scalar gamma_;
//- Coefficient for position and orientation integration (default: 0.25)
const scalar beta_;
public:
//- Runtime type information
TypeName("Newmark");
// Constructors
//- Construct from a dictionary and the body
Newmark
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
);
//- Destructor
virtual ~Newmark();
// Member Functions
//- Drag coefficient
virtual void solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace sixDoFSolvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,57 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "sixDoFSolver.H"
// * * * * * * * * * * * * * * * * Selector * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::sixDoFSolver> Foam::sixDoFSolver::New
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
)
{
word sixDoFSolverType(dict.lookup("type"));
Info<< "Selecting sixDoFSolver " << sixDoFSolverType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(sixDoFSolverType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn("sixDoFSolver::New")
<< "Unknown sixDoFSolverType type "
<< sixDoFSolverType << endl << endl
<< "Valid sixDoFSolver types are : " << endl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return cstrIter()(dict, body);
}
// ************************************************************************* //

View File

@ -0,0 +1,51 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "sixDoFSolver.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(sixDoFSolver, 0);
defineRunTimeSelectionTable(sixDoFSolver, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFSolver::sixDoFSolver(sixDoFRigidBodyMotion& body)
:
body_(body)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFSolver::~sixDoFSolver()
{}
// ************************************************************************* //

View File

@ -0,0 +1,190 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::sixDoFSolver
Description
SourceFiles
sixDoFSolver.C
newSixDoFSolver.C
\*---------------------------------------------------------------------------*/
#ifndef sixDoFSolver_H
#define sixDoFSolver_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "sixDoFRigidBodyMotion.H"
#include "runTimeSelectionTables.H"
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sixDoFSolver Declaration
\*---------------------------------------------------------------------------*/
class sixDoFSolver
{
protected:
// Protected data
//- The rigid body
sixDoFRigidBodyMotion& body_;
// Protected member functions
//- Return the current centre of rotation
inline point& centreOfRotation();
//- Return the orientation
inline tensor& Q();
//- Return non-const access to vector
inline vector& v();
//- Return non-const access to acceleration
inline vector& a();
//- Return non-const access to angular momentum
inline vector& pi();
//- Return non-const access to torque
inline vector& tau();
//- Return the centre of rotation at previous time-step
inline const point& centreOfRotation0() const;
//- Return the orientation at previous time-step
inline const tensor& Q0() const;
//- Return the velocity at previous time-step
inline const vector& v0() const;
//- Return the acceleration at previous time-step
inline const vector& a0() const;
//- Return the angular momentum at previous time-step
inline const vector& pi0() const;
//- Return the torque at previous time-step
inline const vector& tau0() const;
//- Acceleration damping coefficient (for steady-state simulations)
inline scalar aDamp() const;
//- Translational constraint tensor
inline tensor tConstraints() const;
//- Rotational constraint tensor
inline tensor rConstraints() const;
//- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
// and return the rotated Q and pi as a tuple
inline Tuple2<tensor, vector> rotate
(
const tensor& Q0,
const vector& pi,
const scalar deltaT
) const;
//- Update and relax accelerations from the force and torque
inline void updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal
);
public:
//- Runtime type information
TypeName("sixDoFSolver");
// Declare runtime construction
declareRunTimeSelectionTable
(
autoPtr,
sixDoFSolver,
dictionary,
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
),
(dict, body)
);
// Constructors
// Construct for given body
sixDoFSolver(sixDoFRigidBodyMotion& body);
//- Destructor
virtual ~sixDoFSolver();
// Selectors
static autoPtr<sixDoFSolver> New
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
);
// Member Functions
//- Drag coefficient
virtual void solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
) = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "sixDoFSolverI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
inline Foam::point& Foam::sixDoFSolver::centreOfRotation()
{
return body_.motionState_.centreOfRotation();
}
inline Foam::tensor& Foam::sixDoFSolver::Q()
{
return body_.motionState_.Q();
}
inline Foam::vector& Foam::sixDoFSolver::v()
{
return body_.motionState_.v();
}
inline Foam::vector& Foam::sixDoFSolver::a()
{
return body_.motionState_.a();
}
inline Foam::vector& Foam::sixDoFSolver::pi()
{
return body_.motionState_.pi();
}
inline Foam::vector& Foam::sixDoFSolver::tau()
{
return body_.motionState_.tau();
}
inline const Foam::point& Foam::sixDoFSolver::centreOfRotation0() const
{
return body_.motionState0_.centreOfRotation();
}
inline const Foam::tensor& Foam::sixDoFSolver::Q0() const
{
return body_.motionState0_.Q();
}
inline const Foam::vector& Foam::sixDoFSolver::v0() const
{
return body_.motionState0_.v();
}
inline const Foam::vector& Foam::sixDoFSolver::a0() const
{
return body_.motionState0_.a();
}
inline const Foam::vector& Foam::sixDoFSolver::pi0() const
{
return body_.motionState0_.pi();
}
inline const Foam::vector& Foam::sixDoFSolver::tau0() const
{
return body_.motionState0_.tau();
}
inline Foam::scalar Foam::sixDoFSolver::aDamp() const
{
return body_.aDamp_;
}
inline Foam::tensor Foam::sixDoFSolver::tConstraints() const
{
return body_.tConstraints_;
}
inline Foam::tensor Foam::sixDoFSolver::rConstraints() const
{
return body_.rConstraints_;
}
//- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
// and return the rotated Q and pi as a tuple
inline Foam::Tuple2<Foam::tensor, Foam::vector> Foam::sixDoFSolver::rotate
(
const tensor& Q0,
const vector& pi,
const scalar deltaT
) const
{
return body_.rotate(Q0, pi, deltaT);
}
//- Update and relax accelerations from the force and torque
inline void Foam::sixDoFSolver::updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal
)
{
body_.updateAcceleration(fGlobal, tauGlobal);
}
// ************************************************************************* //

View File

@ -0,0 +1,102 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "symplectic.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFSolvers
{
defineTypeNameAndDebug(symplectic, 0);
addToRunTimeSelectionTable(sixDoFSolver, symplectic, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::symplectic::symplectic
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
)
:
sixDoFSolver(body)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::symplectic::~symplectic()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::sixDoFSolvers::symplectic::solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
)
{
if (!firstIter)
{
FatalErrorIn("sixDoFSolvers::symplectic::solve")
<< "The symplectic integrator is explicit "
"and can only be solved once per time-step"
<< exit(FatalError);
}
// First simplectic step:
// Half-step for linear and angular velocities
// Update position and orientation
v() = tConstraints() & (v0() + aDamp()*0.5*deltaT0*a0());
pi() = rConstraints() & (pi0() + aDamp()*0.5*deltaT0*tau0());
centreOfRotation() = centreOfRotation0() + deltaT*v();
Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
Q() = Qpi.first();
pi() = rConstraints() & Qpi.second();
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);
// Second simplectic step:
// Complete update of linear and angular velocities
v() += tConstraints() & aDamp()*0.5*deltaT*a();
pi() += rConstraints() & aDamp()*0.5*deltaT*tau();
}
// ************************************************************************* //

View File

@ -0,0 +1,123 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::sixDoFSolvers::symplectic
Description
Symplectic 2nd-order explicit time-integrator for 6DoF solid-body motion.
Reference:
\verbatim
Dullweber, A., Leimkuhler, B., & McLachlan, R. (1997).
Symplectic splitting methods for rigid body molecular dynamics.
The Journal of chemical physics, 107(15), 5840-5851.
\endverbatim
Can only be used for explicit integration of the motion of the body,
i.e. may only be called once per time-step, no outer-correctors may be
applied. For implicit integration with outer-correctors choose either
CrankNicolson or Newmark schemes.
Example specification in dynamicMeshDict:
\verbatim
solver
{
type symplectic;
}
\endverbatim
SeeAlso
Foam::sixDoFSolvers::CrankNicolson
Foam::sixDoFSolvers::Newmark
SourceFiles
symplectic.C
\*---------------------------------------------------------------------------*/
#ifndef symplectic_H
#define symplectic_H
#include "sixDoFSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFSolvers
{
/*---------------------------------------------------------------------------*\
Class symplectic Declaration
\*---------------------------------------------------------------------------*/
class symplectic
:
public sixDoFSolver
{
public:
//- Runtime type information
TypeName("symplectic");
// Constructors
//- Construct from a dictionary and the body
symplectic
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
);
//- Destructor
virtual ~symplectic();
// Member Functions
//- Drag coefficient
virtual void solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace sixDoFSolvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -41,6 +41,11 @@ sixDoFRigidBodyMotionCoeffs
rhoInf 1;
report on;
solver
{
type symplectic;
}
constraints
{
yLine

View File

@ -31,9 +31,16 @@ sixDoFRigidBodyMotionCoeffs
momentOfInertia (40 921 921);
rhoInf 1;
report on;
accelerationRelaxation 0.3;
value uniform (0 0 0);
accelerationRelaxation 0.4;
solver
{
type Newmark;
}
constraints
{
zAxis

View File

@ -60,7 +60,7 @@ solvers
cacheAgglomeration true;
tolerance 5e-8;
relTol 0.001;
relTol 0;
};
p_rghFinal
@ -86,8 +86,8 @@ PIMPLE
{
momentumPredictor no;
nOuterCorrectors 1;
nCorrectors 3;
nOuterCorrectors 3;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
correctPhi yes;

View File

@ -59,7 +59,12 @@ sixDoFRigidBodyMotionCoeffs
};
report on;
accelerationRelaxation 0.3;
accelerationRelaxation 0.7;
solver
{
type Newmark;
}
constraints
{

View File

@ -117,7 +117,7 @@ solvers
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 5;
nOuterCorrectors 3;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
correctPhi yes;

View File

@ -1,67 +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;
object dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dynamicFvMesh dynamicMotionSolverFvMesh;
motionSolverLibs ("libsixDoFRigidBodyMotion.so");
solver sixDoFRigidBodyMotion;
sixDoFRigidBodyMotionCoeffs
{
patches (hull);
innerDistance 0.3;
outerDistance 1;
centreOfMass (2.929541 0 0.2);
mass 412.73;
momentOfInertia (40 921 921);
rhoInf 1;
report on;
accelerationRelaxation 0.3;
value uniform (0 0 0);
constraints
{
zAxis
{
sixDoFRigidBodyMotionConstraint line;
direction (0 0 1);
}
yPlane
{
sixDoFRigidBodyMotionConstraint axis;
axis (0 1 0);
}
}
restraints
{
translationDamper
{
sixDoFRigidBodyMotionRestraint linearDamper;
coeff 8596;
}
rotationDamper
{
sixDoFRigidBodyMotionRestraint sphericalAngularDamper;
coeff 11586;
}
}
}
// ************************************************************************* //

View File

@ -22,19 +22,22 @@ boundaryField
{
inlet
{
type fixedFluxPressure;
value $internalField;
type fixedFluxPressure;
value $internalField;
}
outlet
{
type prghPressure;
p $internalField;
value $internalField;
type prghTotalPressure;
p0 $internalField;
U U.air;
phi phi.air;
rho rho.air;
value $internalField;
}
walls
{
type fixedFluxPressure;
value $internalField;
type fixedFluxPressure;
value $internalField;
}
}

View File

@ -87,6 +87,11 @@ PIMPLE
relaxationFactors
{
fields
{
iDmdt 1;
}
equations
{
".*" 1;