From c340040b40509331b9949d296d79bf1cbfd906d2 Mon Sep 17 00:00:00 2001 From: Hamidreza Date: Tue, 22 Jul 2025 13:54:23 +0330 Subject: [PATCH] Normal vector of wall is included into the wall velocity calculations --- .../grainInteractionKernels.hpp | 6 +- .../sphereInteractionKernels.hpp | 6 +- src/MotionModel/MotionModel/MotionModel.hpp | 8 +- .../entities/conveyorBelt/conveyorBelt.cpp | 32 +++- .../entities/conveyorBelt/conveyorBelt.hpp | 83 +++++----- .../multiRotatingAxis/multiRotatingAxis.hpp | 6 +- .../entities/rotatingAxis/rotatingAxis.hpp | 2 +- .../entities/rotatingAxis/rotatingAxisI.hpp | 2 +- .../entities/stationary/stationary.hpp | 2 +- .../entities/vibrating/vibrating.hpp | 142 +++++++++--------- src/phasicFlow/triSurface/triSurface.hpp | 20 ++- .../conveyorBelt/settings/geometryDict | 8 +- .../conveyorBelt/settings/settingsDict | 2 +- .../screwConveyor/settings/geometryDict | 2 +- 14 files changed, 187 insertions(+), 134 deletions(-) diff --git a/src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp b/src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp index 7ee5ef39..cfa3f543 100644 --- a/src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp +++ b/src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp @@ -250,7 +250,9 @@ struct pwInteractionFunctor real cGFj = coarseGrainFactor_[i]; realx3 xi = pos_[i]; - realx3x3 tri = triangles_(tj); + realx3x3 tri = triangles_(tj); + const realx3& normW = triangles_.normal(tj); + real ovrlp; realx3 Nij, cp; @@ -262,7 +264,7 @@ struct pwInteractionFunctor int32 mInd = wTriMotionIndex_[tj]; - auto Vw = motionModel_(mInd, cp); + auto Vw = motionModel_(mInd, cp, normW); //output<< "par-wall index "<< i<<" - "<< tj<("tangentVelocity"); + linearVelocity_ = dict.getVal("linearVelocity"); + normal_ = dict.getVal("normal"); + + if(normal_.length() > verySmallValue) + { + normal_.normalize(); + } + else + { + fatalErrorInFunction<< + " normal vector in "<< + dict.globalName() << + " cannot be zero vector "<>(iIstream& is, conveyorBelt& obj) { - - return is; + + return is; } } diff --git a/src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.hpp b/src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.hpp index d39c6dbd..ef31c7c5 100644 --- a/src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.hpp +++ b/src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.hpp @@ -120,7 +120,7 @@ public: /// Tangential velocity at point p INLINE_FUNCTION_HD - realx3 pointTangentialVel(const realx3& p)const + realx3 pointTangentialVel(const realx3& p, const realx3& wallNormal)const { realx3 parentVel(0); auto parIndex = parentAxisIndex(); @@ -128,11 +128,11 @@ public: while(parIndex != -1) { auto& ax = axisList_[parIndex]; - parentVel += ax.linVelocityPoint(p); + parentVel += ax.linVelocityPoint(p, wallNormal); parIndex = ax.parentAxisIndex(); } - return parentVel + rotatingAxis::linVelocityPoint(p); + return parentVel + rotatingAxis::linVelocityPoint(p, wallNormal); } /// Translate point p for dt seconds based on the axis information diff --git a/src/MotionModel/entities/rotatingAxis/rotatingAxis.hpp b/src/MotionModel/entities/rotatingAxis/rotatingAxis.hpp index 7d40f06e..41697824 100644 --- a/src/MotionModel/entities/rotatingAxis/rotatingAxis.hpp +++ b/src/MotionModel/entities/rotatingAxis/rotatingAxis.hpp @@ -126,7 +126,7 @@ public: /// Linear tangential velocity at point p INLINE_FUNCTION_HD - realx3 linVelocityPoint(const realx3 &p)const; + realx3 linVelocityPoint(const realx3 &p, const realx3& wallNormal)const; INLINE_FUNCTION_HD realx3 transferPoint(const realx3 p, real dt)const; diff --git a/src/MotionModel/entities/rotatingAxis/rotatingAxisI.hpp b/src/MotionModel/entities/rotatingAxis/rotatingAxisI.hpp index 10a730cc..d7b4912c 100755 --- a/src/MotionModel/entities/rotatingAxis/rotatingAxisI.hpp +++ b/src/MotionModel/entities/rotatingAxis/rotatingAxisI.hpp @@ -20,7 +20,7 @@ Licence: -----------------------------------------------------------------------------*/ INLINE_FUNCTION_HD -pFlow::realx3 pFlow::rotatingAxis::linVelocityPoint(const realx3 &p)const +pFlow::realx3 pFlow::rotatingAxis::linVelocityPoint(const realx3 &p, const realx3&)const { if(!inTimeRange()) return {0,0,0}; diff --git a/src/MotionModel/entities/stationary/stationary.hpp b/src/MotionModel/entities/stationary/stationary.hpp index 2970e197..739c77ff 100644 --- a/src/MotionModel/entities/stationary/stationary.hpp +++ b/src/MotionModel/entities/stationary/stationary.hpp @@ -60,7 +60,7 @@ public: {} INLINE_FUNCTION_HD - realx3 linVelocityPoint(const realx3 &)const + realx3 linVelocityPoint(const realx3 &, const realx3&)const { return realx3(0); } diff --git a/src/MotionModel/entities/vibrating/vibrating.hpp b/src/MotionModel/entities/vibrating/vibrating.hpp index 62241dfa..9b183bb4 100644 --- a/src/MotionModel/entities/vibrating/vibrating.hpp +++ b/src/MotionModel/entities/vibrating/vibrating.hpp @@ -39,17 +39,17 @@ class dictionary; * Creates a sinoidal virating model for a wall. The viration is defined by * frequency, amplitude and phase angle. * \f[ - \vec{v}(t) = \vec{A} sin(\vec{\omega}(t-startTime)+\vec{\phi}) + \vec{v}(t) = \vec{A} sin(\vec{\omega}(t-startTime)+\vec{\phi}) \f] \verbatim // This creates sinoidal vibration on the wall in x-direction. The viration // starts at t = 0 s and ends at t = 10 s. { - angularFreq (10 0 0); - amplitude ( 1 0 0); - phaseAngle ( 0 0 0); - startTime 0; - endTime 10; + angularFreq (10 0 0); + amplitude ( 1 0 0); + phaseAngle ( 0 0 0); + startTime 0; + endTime 10; } \endverbatim * * * @@ -64,102 +64,102 @@ class dictionary; */ class vibrating : - public timeInterval + public timeInterval { private: - - // rotation speed - realx3 angularFreq_{0,0,0}; + + // rotation speed + realx3 angularFreq_{0,0,0}; - realx3 phaseAngle_{0,0,0}; + realx3 phaseAngle_{0,0,0}; - realx3 amplitude_{0,0,0}; + realx3 amplitude_{0,0,0}; - realx3 velocity_{0,0,0}; + realx3 velocity_{0,0,0}; - realx3 velocity0_{0,0,0}; + realx3 velocity0_{0,0,0}; - INLINE_FUNCTION_HD - void calculateVelocity() - { - if(inTimeRange()) - { - velocity_ = amplitude_ * sin( angularFreq_*(time()-startTime() ) + phaseAngle_ ); - }else - { - velocity_ = {0,0,0}; - } - } + INLINE_FUNCTION_HD + void calculateVelocity() + { + if(inTimeRange()) + { + velocity_ = amplitude_ * sin( angularFreq_*(time()-startTime() ) + phaseAngle_ ); + }else + { + velocity_ = {0,0,0}; + } + } public: - TypeInfoNV("vibrating"); + TypeInfoNV("vibrating"); - FUNCTION_HD - vibrating()=default; + FUNCTION_HD + vibrating()=default; - FUNCTION_H - explicit vibrating(const dictionary& dict); - + FUNCTION_H + explicit vibrating(const dictionary& dict); + - FUNCTION_HD - vibrating(const vibrating&) = default; + FUNCTION_HD + vibrating(const vibrating&) = default; - vibrating& operator=(const vibrating&) = default; + vibrating& operator=(const vibrating&) = default; - INLINE_FUNCTION_HD - void setTime(real t) - { - if( !equal(t, time()) ) velocity0_ = velocity_; - timeInterval::setTime(t); - calculateVelocity(); - } + INLINE_FUNCTION_HD + void setTime(real t) + { + if( !equal(t, time()) ) velocity0_ = velocity_; + timeInterval::setTime(t); + calculateVelocity(); + } - INLINE_FUNCTION_HD - realx3 linVelocityPoint(const realx3 &)const - { - return velocity_; - } + INLINE_FUNCTION_HD + realx3 linVelocityPoint(const realx3 &, const realx3&)const + { + return velocity_; + } - INLINE_FUNCTION_HD - realx3 transferPoint(const realx3& p, real dt)const - { - if(!inTimeRange()) return p; - return p + static_cast(0.5*dt)*(velocity0_+velocity_); - } + INLINE_FUNCTION_HD + realx3 transferPoint(const realx3& p, real dt)const + { + if(!inTimeRange()) return p; + return p + static_cast(0.5*dt)*(velocity0_+velocity_); + } - // - IO operation - FUNCTION_H - bool read(const dictionary& dict); + // - IO operation + FUNCTION_H + bool read(const dictionary& dict); - FUNCTION_H - bool write(dictionary& dict) const; + FUNCTION_H + bool write(dictionary& dict) const; - FUNCTION_H - bool read(iIstream& is); + FUNCTION_H + bool read(iIstream& is); - FUNCTION_H - bool write(iOstream& os)const; + FUNCTION_H + bool write(iOstream& os)const; }; inline iOstream& operator <<(iOstream& os, const vibrating& obj) { - if(!obj.write(os)) - { - fatalExit; - } - return os; + if(!obj.write(os)) + { + fatalExit; + } + return os; } inline iIstream& operator >>(iIstream& is, vibrating& obj) { - if( !obj.read(is) ) - { - fatalExit; - } - return is; + if( !obj.read(is) ) + { + fatalExit; + } + return is; } } diff --git a/src/phasicFlow/triSurface/triSurface.hpp b/src/phasicFlow/triSurface/triSurface.hpp index 204e0376..d01c5841 100644 --- a/src/phasicFlow/triSurface/triSurface.hpp +++ b/src/phasicFlow/triSurface/triSurface.hpp @@ -46,6 +46,9 @@ private: deviceViewType1D dPoints_; deviceViewType1D dVectices_; + + deviceViewType1D dNormals_; + public: INLINE_FUNCTION_H @@ -53,12 +56,14 @@ public: uint32 numPoints, deviceViewType1D points, uint32 numTrianlges, - deviceViewType1D vertices ) + deviceViewType1D vertices, + deviceViewType1D normals ) : numPoints_(numPoints), numTriangles_(numTrianlges), dPoints_(points), - dVectices_(vertices) + dVectices_(vertices), + dNormals_(normals) {} INLINE_FUNCTION_HD @@ -91,11 +96,17 @@ public: INLINE_FUNCTION_HD realx3x3 operator[](uint32 i)const { return triangle(i); } + INLINE_FUNCTION_HD + const realx3& normal(uint32 i)const + { + return dNormals_[i]; + } + INLINE_FUNCTION_HD uint32 numPoints()const { return numPoints_; } INLINE_FUNCTION_HD - uint32 numTrianlges()const { return numTriangles_;} + uint32 numTriangles()const { return numTriangles_;} }; @@ -224,7 +235,8 @@ public: points_.size(), points_.deviceViewAll(), vertices_.size(), - vertices_.deviceViewAll()); + vertices_.deviceViewAll(), + normals_.deviceViewAll() ); } //// - IO operations diff --git a/tutorials/sphereGranFlow/conveyorBelt/settings/geometryDict b/tutorials/sphereGranFlow/conveyorBelt/settings/geometryDict index 251e94b7..e810bc03 100755 --- a/tutorials/sphereGranFlow/conveyorBelt/settings/geometryDict +++ b/tutorials/sphereGranFlow/conveyorBelt/settings/geometryDict @@ -7,14 +7,18 @@ objectType dictionary; fileFormat ASCII; /*---------------------------------------------------------------------------*/ -// motion model can be rotatingAxis or stationary or vibrating motionModel conveyorBelt; conveyorBeltInfo { conveyorBelt1 { - tangentVelocity (0.5 0 0); + // linear velocity of belt + linearVelocity 0.5; + + // normal vector of velocity plate + // The velocity vector is tangent to this plane + normal (0 -1 0); } } diff --git a/tutorials/sphereGranFlow/conveyorBelt/settings/settingsDict b/tutorials/sphereGranFlow/conveyorBelt/settings/settingsDict index 3d3f5658..71ed02c4 100755 --- a/tutorials/sphereGranFlow/conveyorBelt/settings/settingsDict +++ b/tutorials/sphereGranFlow/conveyorBelt/settings/settingsDict @@ -34,7 +34,7 @@ writeFormat ascii; // data writting format (ascii timersReport Yes; // report timers -timersReportInterval 0.01; // time interval for reporting timers +timersReportInterval 0.1; // time interval for reporting timers diff --git a/tutorials/sphereGranFlow/screwConveyor/settings/geometryDict b/tutorials/sphereGranFlow/screwConveyor/settings/geometryDict index 5399819e..37a7999f 100644 --- a/tutorials/sphereGranFlow/screwConveyor/settings/geometryDict +++ b/tutorials/sphereGranFlow/screwConveyor/settings/geometryDict @@ -38,7 +38,7 @@ surfaces type stlWall; // type of the wall file shell.stl; // file name in stl folder material prop1; // material name of this wall - motion none; // this surface is not moving ==> none + motion none; // this surface is not movng ==> none } }