zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

rigidBodyMeshMotion: Added support for dynamic mesh refinement/unrefinement

Browse Source 
The floatingObject tutorial has been update to demonstrate this functionality by
adding the following topoChanger entry to dynamicMeshDict:

topoChanger
{
    type    refiner;

    libs    ("libfvMeshTopoChangers.so");

    // How often to refine
    refineInterval  1;

    // Field to be refinement on
    field           alpha.water;

    // Refine field in between lower..upper
    lowerRefineLevel 0.001;
    upperRefineLevel 0.999;

    // Have slower than 2:1 refinement
    nBufferLayers   1;

    // Refine cells only up to maxRefinement levels
    maxRefinement   1;

    // Stop refinement if maxCells reached
    maxCells        200000;

    // Flux field and corresponding velocity field. Fluxes on changed
    // faces get recalculated by interpolating the velocity. Use 'none'
    // on surfaceScalarFields that do not need to be reinterpolated.
    correctFluxes
    (
        (phi none)
        (nHatf none)
        (rhoPhi none)
        (alphaPhi.water none)
        (meshPhi none)
        (ghf none)
    );

    // Write the refinement level as a volScalarField
    dumpLevel       true;
}

Note that currently only single rigid body motion is supported (but multi-body
support will be added shortly) and the Crank-Nicolson scheme is not supported.
This commit is contained in:
Henry Weller
2021-11-02 14:11:52 +00:00
parent 465b9382c0
commit 37c7d6b9ac
12 changed files with 356 additions and 209 deletions
Download Patch File Download Diff File Expand all files Collapse all files

45
src/OpenFOAM/fields/Fields/transformField/transformField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,6 +29,45 @@ License
// * * * * * * * * * * * * * * * global functions * * * * * * * * * * * * * //
void Foam::transformPoints
(
vectorField& rtf,
const spatialTransform& tr,
const vectorField& tf
)
{
forAll(rtf, i)
{
rtf[i] = tr.transformPoint(tf[i]);
}
}
Foam::tmp<Foam::vectorField> Foam::transformPoints
(
const spatialTransform& tr,
const vectorField& tf
)
{
tmp<vectorField > tranf(new vectorField(tf.size()));
transformPoints(tranf.ref(), tr, tf);
return tranf;
}
Foam::tmp<Foam::vectorField> Foam::transformPoints
(
const spatialTransform& tr,
const tmp<vectorField>& ttf
)
{
tmp<vectorField > tranf = New(ttf);
transformPoints(tranf.ref(), tr, ttf());
ttf.clear();
return tranf;
}
void Foam::transform
(
vectorField& rtf,
@ -36,7 +75,7 @@ void Foam::transform
const vectorField& tf
)
{
tensor t = q.R();
const tensor t = q.R();
TFOR_ALL_F_OP_FUNC_S_F(vector, rtf, =, transform, tensor, t, vector, tf)
}
@ -73,7 +112,7 @@ void Foam::transformPoints
const vectorField& tf
)
{
vector T = tr.t();
const vector T = tr.t();
// Check if any translation
if (mag(T) > vSmall)

18
src/OpenFOAM/fields/Fields/transformField/transformField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -36,6 +36,7 @@ SourceFiles
#define transformField_H
#include "transform.H"
#include "spatialTransform.H"
#include "quaternion.H"
#include "septernion.H"
#include "vectorField.H"
@ -111,6 +112,21 @@ tmp<Field<sphericalTensor>> transformFieldMask<sphericalTensor>
);
//- Transform given vectorField of coordinates with the given spatialTransform
void transformPoints(vectorField&, const spatialTransform&, const vectorField&);
//- Transform given vectorField of coordinates with the given spatialTransform
tmp<vectorField> transformPoints(const spatialTransform&, const vectorField&);
//- Transform given tmp<vectorField> of coordinates
// with the given spatialTransform
tmp<vectorField> transformPoints
(
const spatialTransform&,
const tmp<vectorField>&
);
//- Rotate given vectorField with the given quaternion
void transform(vectorField&, const quaternion&, const vectorField&);

6
src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionSolver/solidBodyMotionSolver.C
View File

@ -186,10 +186,6 @@ void Foam::solidBodyMotionSolver::updateMesh(const mapPolyMesh& mpm)
motionSolver::updateMesh(mpm);
// Map points0_. Bit special since we somehow have to come up with
// a sensible points0 position for introduced points.
// Find out scaling between points0 and current points
// Get the new points either from the map or the mesh
const pointField& points =
(
@ -206,7 +202,7 @@ void Foam::solidBodyMotionSolver::updateMesh(const mapPolyMesh& mpm)
if (oldPointi >= 0)
{
label masterPointi = mpm.reversePointMap()[oldPointi];
const label masterPointi = mpm.reversePointMap()[oldPointi];
if (masterPointi == pointi)
{

13
src/fvMeshTopoChangers/refiner/fvMeshTopoChangersRefiner.C
View File

@ -1072,7 +1072,8 @@ Foam::fvMeshTopoChangers::refiner::refiner(fvMesh& mesh)
meshCutter_(mesh),
dumpLevel_(false),
nRefinementIterations_(0),
protectedCells_(mesh.nCells(), 0)
protectedCells_(mesh.nCells(), 0),
timeIndex_(-1)
{
// Read static part of dictionary
readDict();
@ -1266,6 +1267,16 @@ Foam::fvMeshTopoChangers::refiner::~refiner()
bool Foam::fvMeshTopoChangers::refiner::update()
{
// Only refine on the first call it a time-step
if (timeIndex_ != mesh().time().timeIndex())
{
timeIndex_ = mesh().time().timeIndex();
}
else
{
return false;
}
// Re-read dictionary. Chosen since usually -small so trivial amount
// of time compared to actual refinement. Also very useful to be able
// to modify on-the-fly.

3
src/fvMeshTopoChangers/refiner/fvMeshTopoChangersRefiner.H
View File

@ -186,6 +186,9 @@ class refiner
//- Protected cells (usually since not hexes)
PackedBoolList protectedCells_;
//- The time index used for updating
label timeIndex_;
// Private Member Functions

120
src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C
View File

@ -181,126 +181,12 @@ void Foam::RBD::rigidBodyMotion::status(const label bodyID) const
}
Foam::tmp<Foam::pointField> Foam::RBD::rigidBodyMotion::transformPoints
Foam::spatialTransform Foam::RBD::rigidBodyMotion::transform0
(
const label bodyID,
const pointField& initialPoints
const label bodyID
) const
{
// Calculate the transform from the initial state in the global frame
// to the current state in the global frame
spatialTransform X(X0(bodyID).inv() & X00(bodyID));
tmp<pointField> tpoints(new pointField(initialPoints.size()));
pointField& points = tpoints.ref();
forAll(points, i)
{
points[i] = X.transformPoint(initialPoints[i]);
}
return tpoints;
}
Foam::tmp<Foam::pointField> Foam::RBD::rigidBodyMotion::transformPoints
(
const label bodyID,
const scalarField& weight,
const pointField& initialPoints
) const
{
// Calculate the transform from the initial state in the global frame
// to the current state in the global frame
spatialTransform X(X0(bodyID).inv() & X00(bodyID));
// Calculate the septernion equivalent of the transformation for 'slerp'
// interpolation
septernion s(X);
tmp<pointField> tpoints(new pointField(initialPoints));
pointField& points = tpoints.ref();
forAll(points, i)
{
// Move non-stationary points
if (weight[i] > small)
{
// Use solid-body motion where weight = 1
if (weight[i] > 1 - small)
{
points[i] = X.transformPoint(initialPoints[i]);
}
// Slerp septernion interpolation
else
{
points[i] =
slerp(septernion::I, s, weight[i])
.transformPoint(initialPoints[i]);
}
}
}
return tpoints;
}
Foam::tmp<Foam::pointField> Foam::RBD::rigidBodyMotion::transformPoints
(
const labelList& bodyIDs,
const List<const scalarField*>& weights,
const pointField& initialPoints
) const
{
List<septernion> ss(bodyIDs.size() + 1);
ss[bodyIDs.size()] = septernion::I;
forAll(bodyIDs, bi)
{
const label bodyID = bodyIDs[bi];
// Calculate the transform from the initial state in the global frame
// to the current state in the global frame
spatialTransform X(X0(bodyID).inv() & X00(bodyID));
// Calculate the septernion equivalent of the transformation
ss[bi] = septernion(X);
}
tmp<pointField> tpoints(new pointField(initialPoints));
pointField& points = tpoints.ref();
List<scalar> w(ss.size());
forAll(points, i)
{
// Initialise to 1 for the far-field weight
scalar sum1mw = 1;
forAll(bodyIDs, bi)
{
w[bi] = (*(weights[bi]))[i];
sum1mw += w[bi]/(1 + small - w[bi]);
}
// Calculate the limiter for wi/(1 - wi) to ensure the sum(wi) = 1
scalar lambda = 1/sum1mw;
// Limit wi/(1 - wi) and sum the resulting wi
scalar sumw = 0;
forAll(bodyIDs, bi)
{
w[bi] = lambda*w[bi]/(1 + small - w[bi]);
sumw += w[bi];
}
// Calculate the weight for the stationary far-field
w[bodyIDs.size()] = 1 - sumw;
points[i] = average(ss, w).transformPoint(initialPoints[i]);
}
return tpoints;
return X0(bodyID).inv() & X00(bodyID);
}

29
src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H
View File

@ -181,33 +181,8 @@ public:
// Transformations
//- Transform the given initial pointField of the specified body
// to correspond to the current motion state
tmp<pointField> transformPoints
(
const label bodyID,
const pointField& initialPoints
) const;
//- Transform the given initial pointField of the specified body
// to correspond to the current motion state scaled using
// 'slerp' interpolation
tmp<pointField> transformPoints
(
const label bodyID,
const scalarField& weight,
const pointField& initialPoints
) const;
//- Transform the given initial pointField of the specified body
// to correspond to the current motion state scaled using
// 'slerp' interpolation
tmp<pointField> transformPoints
(
const labelList& bodyIDs,
const List<const scalarField*>& weights,
const pointField& initialPoints
) const;
//- Return the transformation of bodyID relative to the initial time
spatialTransform transform0(const label bodyID) const;
//- Write

267
src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
View File

@ -25,6 +25,7 @@ License
#include "rigidBodyMeshMotion.H"
#include "polyMesh.H"
#include "mapPolyMesh.H"
#include "pointPatchDist.H"
#include "pointConstraints.H"
#include "timeIOdictionary.H"
@ -73,8 +74,7 @@ Foam::rigidBodyMeshMotion::bodyMesh::bodyMesh
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
IOobject::NO_WRITE
),
pointMesh::New(mesh),
dimensionedScalar(dimless, 0)
@ -166,45 +166,17 @@ Foam::rigidBodyMeshMotion::rigidBodyMeshMotion
}
}
const pointMesh& pMesh = pointMesh::New(mesh);
// Calculate scaling factor everywhere for each meshed body
forAll(bodyMeshes_, bi)
{
const pointMesh& pMesh = pointMesh::New(mesh);
const pointPatchDist pDist(pMesh, bodyMeshes_[bi].patchSet_, points0());
pointPatchDist pDist(pMesh, bodyMeshes_[bi].patchSet_, points0());
bodyMeshes_[bi].weight_.primitiveFieldRef() =
bodyMeshes_[bi].weight(pDist.primitiveField());
pointScalarField& scale = bodyMeshes_[bi].weight_;
// Scaling: 1 up to di then linear down to 0 at do away from patches
scale.primitiveFieldRef() =
min
(
max
(
(bodyMeshes_[bi].do_ - pDist.primitiveField())
/(bodyMeshes_[bi].do_ - bodyMeshes_[bi].di_),
scalar(0)
),
scalar(1)
);
// Convert the scale function to a cosine
scale.primitiveFieldRef() =
min
(
max
(
0.5
- 0.5
*cos(scale.primitiveField()
*Foam::constant::mathematical::pi),
scalar(0)
),
scalar(1)
);
pointConstraints::New(pMesh).constrain(scale);
// scale.write();
pointConstraints::New(pMesh).constrain(bodyMeshes_[bi].weight_);
}
}
@ -217,6 +189,31 @@ Foam::rigidBodyMeshMotion::~rigidBodyMeshMotion()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Type Foam::rigidBodyMeshMotion::bodyMesh::weight
(
const Type& pDist
) const
{
// Scaling: 1 up to di then linear down to 0 at do away from patches
Type weight = min(max((do_ - pDist)/(do_ - di_), scalar(0)), scalar(1));
// Convert the weight function to a cosine
weight =
min
(
max
(
0.5 - 0.5*cos(weight*Foam::constant::mathematical::pi),
scalar(0)
),
scalar(1)
);
return weight;
}
Foam::tmp<Foam::pointField>
Foam::rigidBodyMeshMotion::curPoints() const
{
@ -309,25 +306,84 @@ void Foam::rigidBodyMeshMotion::solve()
// Update the displacements
if (bodyMeshes_.size() == 1)
{
pointDisplacement_.primitiveFieldRef() = transformPoints
const septernion transform0
(
bodyMeshes_[0].bodyID_,
bodyMeshes_[0].weight_,
points0()
) - points0();
this->transform0(bodyMeshes_[0].bodyID_)
);
vectorField& pointDisplacement = pointDisplacement_.primitiveFieldRef();
const pointField& points0 = this->points0();
const scalarField& weight = bodyMeshes_[0].weight_;
forAll(points0, pointi)
{
// Move non-stationary points
if (weight[pointi] > small)
{
// Use solid-body motion where weight = 1
if (weight[pointi] > 1 - small)
{
pointDisplacement[pointi] =
transform0.transformPoint(points0[pointi])
- points0[pointi];
}
// Slerp septernion interpolation
else
{
pointDisplacement[pointi] =
slerp(septernion::I, transform0, weight[pointi])
.transformPoint(points0[pointi])
- points0[pointi];
}
}
}
}
else
{
labelList bodyIDs(bodyMeshes_.size());
List<const scalarField*> weights(bodyMeshes_.size());
forAll(bodyIDs, bi)
List<septernion> transforms0(bodyMeshes_.size() + 1);
forAll(bodyMeshes_, bi)
{
bodyIDs[bi] = bodyMeshes_[bi].bodyID_;
weights[bi] = &bodyMeshes_[bi].weight_;
// Calculate the septernion equivalent of the transformation
transforms0[bi] = septernion(transform0(bodyMeshes_[bi].bodyID_));
}
pointDisplacement_.primitiveFieldRef() =
transformPoints(bodyIDs, weights, points0()) - points0();
transforms0[bodyMeshes_.size()] = septernion::I;
vectorField& pointDisplacement = pointDisplacement_.primitiveFieldRef();
const pointField& points0 = this->points0();
List<scalar> w(transforms0.size());
forAll(points0, pointi)
{
// Initialise to 1 for the far-field weight
scalar sum1mw = 1;
forAll(bodyMeshes_, bi)
{
w[bi] = bodyMeshes_[bi].weight_[pointi];
sum1mw += w[bi]/(1 + small - w[bi]);
}
// Calculate the limiter for wi/(1 - wi) to ensure the sum(wi) = 1
scalar lambda = 1/sum1mw;
// Limit wi/(1 - wi) and sum the resulting wi
scalar sumw = 0;
forAll(bodyMeshes_, bi)
{
w[bi] = lambda*w[bi]/(1 + small - w[bi]);
sumw += w[bi];
}
// Calculate the weight for the stationary far-field
w[bodyMeshes_.size()] = 1 - sumw;
pointDisplacement[pointi] =
average(transforms0, w).transformPoint(points0[pointi])
- points0[pointi];
}
}
// Displacement has changed. Update boundary conditions
@ -338,6 +394,121 @@ void Foam::rigidBodyMeshMotion::solve()
}
void Foam::rigidBodyMeshMotion::updateMesh(const mapPolyMesh& mpm)
{
// pointMesh already updates pointFields
motionSolver::updateMesh(mpm);
// Get the new points either from the map or the mesh
const pointField& points =
(
mpm.hasMotionPoints()
? mpm.preMotionPoints()
: mesh().points()
);
const pointMesh& pMesh = pointMesh::New(mesh());
pointField points0(points);
for (int iter=0; iter<3; iter++)
{
// Calculate scaling factor everywhere for each meshed body
forAll(bodyMeshes_, bi)
{
const pointPatchDist pDist(pMesh, bodyMeshes_[bi].patchSet_, points0);
pointScalarField& weight = bodyMeshes_[bi].weight_;
forAll(points0, pointi)
{
const label oldPointi = mpm.pointMap()[pointi];
if (oldPointi >= 0)
{
if (mpm.reversePointMap()[oldPointi] != pointi)
{
weight[pointi] = bodyMeshes_[bi].weight(pDist[pointi]);
}
}
else
{
FatalErrorInFunction
<< "Cannot determine co-ordinates of introduced vertices."
<< " New vertex " << pointi << " at co-ordinate "
<< points[pointi] << exit(FatalError);
}
}
pointConstraints::New(pMesh).constrain(weight);
}
forAll(points0, pointi)
{
const label oldPointi = mpm.pointMap()[pointi];
if (oldPointi >= 0)
{
if (mpm.reversePointMap()[oldPointi] == pointi)
{
points0[pointi] = points0_[oldPointi];
}
else
{
if (bodyMeshes_.size() == 1)
{
// Use solid-body motion where weight = 1
if (bodyMeshes_[0].weight_[pointi] > 1 - small)
{
points0[pointi] =
transform0(bodyMeshes_[0].bodyID_).inv()
.transformPoint(points[pointi]);
}
// Slerp septernion interpolation
else
{
points0[pointi] =
slerp
(
septernion::I,
septernion(transform0(bodyMeshes_[0].bodyID_)),
bodyMeshes_[0].weight_[pointi]
).invTransformPoint(points[pointi]);
}
}
else
{
NotImplemented;
// labelList bodyIDs(bodyMeshes_.size());
// List<const scalarField*> weights(bodyMeshes_.size());
// forAll(bodyIDs, bi)
// {
// bodyIDs[bi] = bodyMeshes_[bi].bodyID_;
// weights[bi] = &bodyMeshes_[bi].weight_;
// }
// points0[pointi] = invTransformPoints
// (
// bodyIDs,
// weights,
// points[pointi]
// );
}
}
}
}
}
points0_.transfer(points0);
// points0 changed - set to write and check-in to database
points0_.rename("points0");
points0_.writeOpt() = IOobject::AUTO_WRITE;
points0_.instance() = mesh().time().timeName();
points0_.checkIn();
}
bool Foam::rigidBodyMeshMotion::write() const
{
timeIOdictionary dict

8
src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -93,6 +93,9 @@ class rigidBodyMeshMotion
const label bodyID,
const dictionary& dict
);
template<class Type>
inline Type weight(const Type& pDist) const;
};
@ -159,6 +162,9 @@ public:
//- Solve for motion
virtual void solve();
//- Update local data for topology changes
virtual void updateMesh(const mapPolyMesh&);
//- Write motion state information for restart
virtual bool write() const;

9
src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C
View File

@ -30,6 +30,7 @@ License
#include "timeIOdictionary.H"
#include "uniformDimensionedFields.H"
#include "forces.H"
#include "transformField.H"
#include "OneConstant.H"
#include "mathematicalConstants.H"
#include "addToRunTimeSelectionTable.H"
@ -272,9 +273,9 @@ void Foam::rigidBodyMeshMotionSolver::solve()
{
forAllConstIter(labelHashSet, bodyMeshes_[bi].patchSet_, iter)
{
label patchi = iter.key();
const label patchi = iter.key();
pointField patchPoints0
const pointField patchPoints0
(
meshSolver_.pointDisplacement().boundaryField()[patchi]
.patchInternalField(meshSolver_.points0())
@ -282,9 +283,9 @@ void Foam::rigidBodyMeshMotionSolver::solve()
meshSolver_.pointDisplacement().boundaryFieldRef()[patchi] ==
(
transformPoints
Foam::transformPoints
(
bodyMeshes_[bi].bodyID_,
transform0(bodyMeshes_[bi].bodyID_),
patchPoints0
) - patchPoints0
)();

43
tutorials/multiphase/interFoam/RAS/floatingObject/constant/dynamicMeshDict
View File

@ -76,4 +76,47 @@ mover
}
topoChanger
{
type refiner;
libs ("libfvMeshTopoChangers.so");
// How often to refine
refineInterval 1;
// Field to be refinement on
field alpha.water;
// Refine field in between lower..upper
lowerRefineLevel 0.001;
upperRefineLevel 0.999;
// Have slower than 2:1 refinement
nBufferLayers 1;
// Refine cells only up to maxRefinement levels
maxRefinement 1;
// Stop refinement if maxCells reached
maxCells 200000;
// Flux field and corresponding velocity field. Fluxes on changed
// faces get recalculated by interpolating the velocity. Use 'none'
// on surfaceScalarFields that do not need to be reinterpolated.
correctFluxes
(
(phi none)
(nHatf none)
(rhoPhi none)
(alphaPhi.water none)
(meshPhi none)
(ghf none)
);
// Write the refinement level as a volScalarField
dumpLevel true;
}
// ************************************************************************* //

4
tutorials/multiphase/interFoam/RAS/floatingObject/system/fvSchemes
View File

@ -16,7 +16,7 @@ FoamFile
ddtSchemes
{
default CrankNicolson 0.9;
default Euler;
}
gradSchemes
@ -26,7 +26,7 @@ gradSchemes
divSchemes
{
div(rhoPhi,U) Gauss vanLeerV;
div(rhoPhi,U) Gauss vanLeerV;
div(phi,alpha) Gauss interfaceCompression vanLeer 1;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;