From da1e6ea8d036dee5ec84c9f434c93737e41baf02 Mon Sep 17 00:00:00 2001 From: sergio Date: Mon, 1 Mar 2021 19:25:21 -0800 Subject: [PATCH 1/2] ENH: Improvement to overset that allows multiple motion solvers to operate with overset 1) Adding zoneMotion to rigidBodyMotion 2) Introducing PID to prescribedRotation restraint 3) Making drivenLinearMotion read total displacement 4) When drivenLinearMotion is used sixDof and rigid-body solvers write total displacement --- .../overInterDyMFoam/overInterDyMFoam.C | 11 ++-- .../points0/points0MotionSolver.C | 14 +++- .../drivenLinearMotion/drivenLinearMotion.C | 29 ++------- .../drivenLinearMotion/drivenLinearMotion.H | 6 -- src/overset/include/setCellMask.H | 3 - .../prescribedRotation/prescribedRotation.C | 58 ++++++++++++----- .../prescribedRotation/prescribedRotation.H | 28 +++++++- .../rigidBodyMotion/rigidBodyMotion.C | 17 +++++ .../rigidBodyMotion/rigidBodyMotion.H | 6 ++ .../rigidBodyMeshMotion/rigidBodyMeshMotion.C | 64 +++++++++++++++++-- .../rigidBodyMeshMotion/rigidBodyMeshMotion.H | 9 ++- .../sixDoFRigidBodyMotionSolver.C | 7 +- .../background/constant/dynamicMeshDict | 3 +- 13 files changed, 188 insertions(+), 67 deletions(-) diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C b/applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C index 4e8c7b2891..1dd554aba7 100644 --- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C @@ -97,9 +97,15 @@ int main(int argc, char *argv[]) dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0) ); - #include "correctPhi.H" + if (correctPhi) + { + #include "correctPhi.H" + } #include "createUf.H" + #include "setCellMask.H" + #include "setInterpolatedCells.H" + turbulence->validate(); if (!LTS) @@ -108,9 +114,6 @@ int main(int argc, char *argv[]) #include "setInitialDeltaT.H" } - #include "setCellMask.H" - #include "setInterpolatedCells.H" - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; diff --git a/src/dynamicMesh/motionSolvers/displacement/points0/points0MotionSolver.C b/src/dynamicMesh/motionSolvers/displacement/points0/points0MotionSolver.C index 9f57740a8e..e0c04e7cfe 100644 --- a/src/dynamicMesh/motionSolvers/displacement/points0/points0MotionSolver.C +++ b/src/dynamicMesh/motionSolvers/displacement/points0/points0MotionSolver.C @@ -48,7 +48,6 @@ Foam::IOobject Foam::points0MotionSolver::points0IO(const polyMesh& mesh) "points0", IOobject::READ_IF_PRESENT ); - IOobject io ( "points0", @@ -71,6 +70,19 @@ Foam::IOobject Foam::points0MotionSolver::points0IO(const polyMesh& mesh) io.rename("points"); } + // Choose the latest points if points0 is not present +// if (!io.typeHeaderOk()) +// { +// const word instance = +// mesh.time().findInstance +// ( +// mesh.meshDir(), +// "points", +// IOobject::MUST_READ +// ); +// io.instance() = instance; +// io.rename("points"); +// } return io; } diff --git a/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.C b/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.C index 5eef22353a..d2da6f71cb 100644 --- a/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.C +++ b/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2019-2020 OpenCFD Ltd. + Copyright (C) 2019-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -56,26 +56,21 @@ Foam::solidBodyMotionFunctions::drivenLinearMotion::drivenLinearMotion : solidBodyMotionFunction(SBMFCoeffs, runTime), CofGvelocity_(SBMFCoeffs.get("CofGvelocity")), - normal_(SBMFCoeffs.getOrDefault("normal", Zero)), CofGvel_ ( IOobject ( CofGvelocity_, time_.timeName(), + "uniform", time_, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), dimensionedVector(dimless, Zero) - ), - displacement_(Zero) + ) { read(SBMFCoeffs); - if (mag(normal_) > SMALL) - { - normal_ /= (mag(normal_) + SMALL); - } } @@ -84,24 +79,10 @@ Foam::solidBodyMotionFunctions::drivenLinearMotion::drivenLinearMotion Foam::septernion Foam::solidBodyMotionFunctions::drivenLinearMotion::transformation() const { - scalar deltaT = time_.deltaT().value(); - - vector velocity = CofGvel_.value(); - - // Take out normal component - if (mag(normal_) > SMALL) - { - velocity = CofGvel_.value() - ((CofGvel_.value() & normal_)*normal_); - } - - DebugInFunction << "Vel on plane :" << velocity << endl; - - // Translation of centre of gravity with constant velocity - //const vector displacement = velocity*t; - displacement_ += velocity*deltaT; + DebugInFunction << "displacement :" << CofGvel_.value() << endl; quaternion R(1); - septernion TR(septernion(-displacement_)*R); + septernion TR(septernion(-CofGvel_.value())*R); DebugInFunction << "Time = " << time_.value() << " transformation: " << TR << endl; diff --git a/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.H b/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.H index 46096124ab..324ba140ad 100644 --- a/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.H +++ b/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.H @@ -66,15 +66,9 @@ class drivenLinearMotion //- Name of the meshObject to dum CofG velocity word CofGvelocity_; - //- Normal plane direction to restrict movement on a plane - vector normal_; - //- Uniform vector to follow uniformDimensionedVectorField CofGvel_; - //- Last displacement - mutable vector displacement_; - // Private Member Functions diff --git a/src/overset/include/setCellMask.H b/src/overset/include/setCellMask.H index d29e7005f1..b4af6fa7dd 100644 --- a/src/overset/include/setCellMask.H +++ b/src/overset/include/setCellMask.H @@ -32,15 +32,12 @@ Description const cellCellStencilObject& overlap = Stencil::New(mesh); const labelList& cellTypes = overlap.cellTypes(); - //label nHoles = 0; - cellMask.primitiveFieldRef() = 1.0; forAll(cellMask, cellI) { if (cellTypes[cellI] == cellCellStencil::HOLE) { cellMask[cellI] = 0.0; - //nHoles++; } } cellMask.correctBoundaryConditions(); diff --git a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C index 87ba49ca39..fd6503fef1 100644 --- a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C +++ b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2018-2020 OpenCFD Ltd. + Copyright (C) 2018-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -61,7 +61,11 @@ Foam::RBD::restraints::prescribedRotation::prescribedRotation ) : restraint(name, dict, model), - omegaSet_(model_.time(), "omega") + omegaSet_(model_.time(), "omega"), + omega_(Zero), + prevMom_(Zero), + error0_(Zero), + integral0_(Zero) { read(dict); } @@ -123,23 +127,31 @@ void Foam::RBD::restraints::prescribedRotation::restrain // Adding rotation to a given body vector omega = model_.v(model_.master(bodyID_)).w(); + scalar Inertia = mag(model_.I(model_.master(bodyID_)).Ic()); // from the definition of the angular momentum: - // moment = Inertia*ddt(omega) - const scalar relax = 0.5; + // moment = Inertia*ddt(omega) - vector moment - ( - ( - relax - * Inertia - * (omegaSet_.value(model_.time().value()) - omega) - / model_.time().deltaTValue() - & a - ) - * a - ); + vector error = omegaSet_.value(model_.time().value()) - omega; + vector integral = integral0_ + error; + vector derivative = (error - error0_); + + vector moment = ((p_*error + i_*integral + d_*derivative) & a ) * a; + moment = moment*Inertia/model_.time().deltaTValue(); + +// vector moment +// ( +// ( +// Inertia +// * (omegaSet_.value(model_.time().value()) - omega) +// / model_.time().deltaTValue()/relax_ +// & a +// ) +// * a +// ); + + moment = relax_*moment + (1- relax_)*prevMom_; if (model_.debug) { @@ -149,12 +161,20 @@ void Foam::RBD::restraints::prescribedRotation::restrain << " moment " << moment << endl << " oldDir " << oldDir << endl << " newDir " << newDir << endl + << " inertia " << Inertia << endl + << " error " << error << endl + << " integral " << integral << endl + << " derivative " << derivative << endl << " refDir " << refDir << endl; } // Accumulate the force for the restrained body - fx[bodyIndex_] += spatialVector(moment, Zero); + fx[bodyIndex_] += model_.X0(bodyID_).T() & spatialVector(moment, Zero); + + prevMom_ = moment; + error0_ = error; + integral0_ = integral; } @@ -180,6 +200,12 @@ bool Foam::RBD::restraints::prescribedRotation::read const scalar magAxis(mag(axis_)); + coeffs_.readEntry("relax", relax_); + + coeffs_.readEntry("p", p_); + coeffs_.readEntry("i", i_); + coeffs_.readEntry("d", d_); + if (magAxis > VSMALL) { axis_ /= magAxis; diff --git a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H index 7533e46b8a..c49a9a113e 100644 --- a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H +++ b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2018 OpenCFD Ltd. + Copyright (C) 2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -44,6 +44,10 @@ Usage referenceOrientation | Orientation | no | I axis | Rotation axis (in reference) | yes | omega | Angular velocity (rad/s) | yes | + relax | Relax moment with prevoius iter | yes + p | Propoptional corrector for PDI | yes + d | Differencial corrector for PDI | yes + i | Integral corrector for PDI | yes \endtable SourceFiles @@ -85,6 +89,26 @@ class prescribedRotation //- Rotational velocity [rad/sec] TimeFunction1 omegaSet_; + //- Cache omega + mutable vector omega_; + + //- Cache previous momentum + mutable vector prevMom_; + + //- Relax momentum + scalar relax_; + + //- PID constants + + mutable vector error0_; + + mutable vector integral0_; + + mutable scalar p_; + + mutable scalar i_; + + mutable scalar d_; public: @@ -137,7 +161,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -} // End namespace RBD +} // End namespace restraints } // End namespace RBD } // End namespace Foam diff --git a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C index fe05c117b0..f0f2b59cc7 100644 --- a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C +++ b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C @@ -186,6 +186,23 @@ void Foam::RBD::rigidBodyMotion::status(const label bodyID) const } +const Foam::vector Foam::RBD::rigidBodyMotion::vCofR(const label bodyID) const +{ + const spatialVector velCofR(v(bodyID, Zero)); + + return velCofR.l(); +} + + +const Foam::vector Foam::RBD::rigidBodyMotion::cCofR(const label bodyID) const +{ + const spatialTransform CofR(X0(bodyID)); + + return CofR.r(); +} + + + Foam::tmp Foam::RBD::rigidBodyMotion::transformPoints ( const label bodyID, diff --git a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H index dd9038aef3..284d5b12b9 100644 --- a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H +++ b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H @@ -189,6 +189,12 @@ public: //- Report the status of the motion of the given body void status(const label bodyID) const; + //- Report linear velocity of the given body + const vector vCofR(const label bodyID) const; + + //- Report CofR of the given body + const vector cCofR(const label bodyID) const; + // Transformations diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C index e27d8d1c80..06291f52a2 100644 --- a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C +++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2016-2017 OpenFOAM Foundation - Copyright (C) 2016-2020 OpenCFD Ltd. + Copyright (C) 2016-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -120,7 +120,10 @@ Foam::rigidBodyMeshMotion::rigidBodyMeshMotion rhoInf_(1.0), rhoName_(coeffDict().getOrDefault("rho", "rho")), ramp_(Function1::NewIfPresent("ramp", coeffDict())), - curTimeIndex_(-1) + curTimeIndex_(-1), + CofGvelocity_(coeffDict().getOrDefault("CofGvelocity", "none")), + bodyIdCofG_(coeffDict().getOrDefault