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