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..6a598cfce1 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", diff --git a/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.C b/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.C index 5eef22353a..177e9b6b96 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. @@ -55,27 +55,22 @@ Foam::solidBodyMotionFunctions::drivenLinearMotion::drivenLinearMotion ) : solidBodyMotionFunction(SBMFCoeffs, runTime), - CofGvelocity_(SBMFCoeffs.get("CofGvelocity")), - normal_(SBMFCoeffs.getOrDefault("normal", Zero)), - CofGvel_ + cOfGdisplacement_(SBMFCoeffs.get("cOfGdisplacement")), + CofGdisp_ ( IOobject ( - CofGvelocity_, + cOfGdisplacement_, 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 :" << CofGdisp_.value() << endl; quaternion R(1); - septernion TR(septernion(-displacement_)*R); + septernion TR(septernion(-CofGdisp_.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..977572b673 100644 --- a/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.H +++ b/src/dynamicMesh/motionSolvers/displacement/solidBody/solidBodyMotionFunctions/drivenLinearMotion/drivenLinearMotion.H @@ -63,17 +63,11 @@ class drivenLinearMotion { // Private data - //- Name of the meshObject to dum CofG velocity - word CofGvelocity_; - - //- Normal plane direction to restrict movement on a plane - vector normal_; + //- Name of the meshObject to dum CofG displacement + word cOfGdisplacement_; //- Uniform vector to follow - uniformDimensionedVectorField CofGvel_; - - //- Last displacement - mutable vector displacement_; + uniformDimensionedVectorField CofGdisp_; // 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..c8de6a89d0 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), + oldMom_(Zero), + error0_(Zero), + integral0_(Zero) { read(dict); } @@ -82,7 +86,7 @@ void Foam::RBD::restraints::prescribedRotation::restrain const rigidBodyModelState& state ) const { - vector refDir = rotationTensor(vector(1, 0, 0), axis_) & vector(0, 1, 0); + vector refDir = rotationTensor(vector(1, 0, 0), axis_)&vector(0, 1, 0); vector oldDir = refQ_ & refDir; vector newDir = model_.X0(bodyID_).E() & refDir; @@ -123,23 +127,20 @@ 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 *= Inertia/model_.time().deltaTValue(); + + moment = relax_*moment + (1- relax_)*oldMom_; if (model_.debug) { @@ -149,12 +150,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); + + oldMom_ = moment; + error0_ = error; + integral0_ = integral; } @@ -180,6 +189,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..a587da3eb5 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 oldMom_; + + //- 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..1f91e02ebf 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,12 @@ Foam::rigidBodyMeshMotion::rigidBodyMeshMotion rhoInf_(1.0), rhoName_(coeffDict().getOrDefault("rho", "rho")), ramp_(Function1::NewIfPresent("ramp", coeffDict())), - curTimeIndex_(-1) + curTimeIndex_(-1), + cOfGdisplacement_ + ( + coeffDict().getOrDefault("cOfGdisplacement", "none") + ), + bodyIdCofG_(coeffDict().getOrDefault