ENH: helical: new particle force model to follow helical trajectory

src/Pstream/mpi/UPstreamRequest.C

@@ -855,7 +855,7 @@ bool Foam::UPstream::finishedRequest(const label i)
     // This allows MPI to progress behind the scenes if it wishes.
 
     int flag = 0;
-    if (i < 0 || i >= PstreamGlobals::outstandingRequests_.size())
+    if (i >= 0 && i < PstreamGlobals::outstandingRequests_.size())
     {
         auto& request = PstreamGlobals::outstandingRequests_[i];
 

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C

@@ -218,9 +218,11 @@ tmp<vectorField> atmBoundaryLayer::U(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Table 1, RH:Eq. 6, HW:Eq. 5)
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
     scalarField Un
     (
-        (Ustar(z0)/kappa_)*log(((zDir() & pCf) - groundMin - d + z0)/z0)
+        (Ustar(z0)/kappa_)*log(zEff/z0)
     );
 
     return flowDir()*Un;
@@ -235,9 +237,9 @@ tmp<scalarField> atmBoundaryLayer::k(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Eq. 21; RH:Eq. 7, HW:Eq. 6 when C1=0 and C2=1)
-    return
-        sqr(Ustar(z0))/sqrt(Cmu_)
-       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
+    return sqr(Ustar(z0))/sqrt(Cmu_)*sqrt(C1_*log(zEff/z0) + C2_);
 }
 
 
@@ -249,9 +251,9 @@ tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Eq. 22; RH:Eq. 8, HW:Eq. 7 when C1=0 and C2=1)
-    return
-        pow3(Ustar(z0))/(kappa_*((zDir() & pCf) - groundMin - d + z0))
-       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
+    return pow3(Ustar(z0))/(kappa_*zEff)*sqrt(C1_*log(zEff/z0) + C2_);
 }
 
 
@@ -263,7 +265,9 @@ tmp<scalarField> atmBoundaryLayer::omega(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGJ:Eq. 13)
-    return Ustar(z0)/(kappa_*sqrt(Cmu_)*((zDir() & pCf) - groundMin - d + z0));
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
+    return Ustar(z0)/(kappa_*sqrt(Cmu_)*zEff);
 }
 
 

src/lagrangian/intermediate/parcels/include/makeParcelForces.H

@@ -49,6 +49,7 @@ License
 #include "VirtualMassForce.H"
 #include "InterfaceForce.H"
 #include "CoulombForce.H"
+#include "HelicalForce.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,7 +71,8 @@ License
     makeParticleForceModelType(SRFForce, CloudType);                           \
     makeParticleForceModelType(VirtualMassForce, CloudType);                   \
     makeParticleForceModelType(InterfaceForce, CloudType);                     \
-    makeParticleForceModelType(CoulombForce, CloudType);
+    makeParticleForceModelType(CoulombForce, CloudType);                       \
+    makeParticleForceModelType(HelicalForce, CloudType);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/intermediate/parcels/include/makeThermoParcelForces.H

@@ -45,6 +45,7 @@ License
 #include "SRFForce.H"
 #include "VirtualMassForce.H"
 #include "CoulombForce.H"
+#include "HelicalForce.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,7 +63,8 @@ License
     makeParticleForceModelType(PressureGradientForce, CloudType);              \
     makeParticleForceModelType(SRFForce, CloudType);                           \
     makeParticleForceModelType(VirtualMassForce, CloudType);                   \
-    makeParticleForceModelType(CoulombForce, CloudType);
+    makeParticleForceModelType(CoulombForce, CloudType);                       \
+    makeParticleForceModelType(HelicalForce, CloudType);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Helical/HelicalForce.C

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "HelicalForce.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::HelicalForce<CloudType>::HelicalForce
+(
+    CloudType& owner,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
+    axisPointPtr_(Function1<vector>::New("axisPoint", this->coeffs(), &mesh)),
+    axisDirPtr_(Function1<vector>::New("axisDir", this->coeffs(), &mesh)),
+    omegaPtr_(Function1<scalar>::New("omega", this->coeffs(), &mesh)),
+    UaxialPtr_
+    (
+        Function1<scalar>::NewIfPresent("Uaxial", this->coeffs(), &mesh)
+    ),
+    tauRotationalPtr_
+    (
+        Function1<scalar>::NewIfPresent("tauRotational", this->coeffs(), &mesh)
+    ),
+    tauAxialPtr_
+    (
+        Function1<scalar>::NewIfPresent("tauAxial", this->coeffs(), &mesh)
+    ),
+    radialDampingPtr_
+    (
+        Function1<scalar>::NewIfPresent("radialDamping", this->coeffs(), &mesh)
+    )
+{
+    if (tauAxialPtr_ && !UaxialPtr_)
+    {
+        FatalIOErrorInFunction(this->coeffs())
+            << "If tauAxial is specified, Uaxial must also be specified." << nl
+            << error(FatalIOError);
+    }
+}
+
+
+template<class CloudType>
+Foam::HelicalForce<CloudType>::HelicalForce
+(
+    const HelicalForce<CloudType>& df
+)
+:
+    ParticleForce<CloudType>(df),
+    axisPointPtr_(df.axisPointPtr_.clone()),
+    axisDirPtr_(df.axisDirPtr_.clone()),
+    omegaPtr_(df.omegaPtr_.clone()),
+    UaxialPtr_(df.UaxialPtr_.clone()),
+    tauRotationalPtr_(df.tauRotationalPtr_.clone()),
+    tauAxialPtr_(df.tauAxialPtr_.clone()),
+    radialDampingPtr_(df.radialDampingPtr_.clone())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::HelicalForce<CloudType>::~HelicalForce()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::forceSuSp Foam::HelicalForce<CloudType>::calcCoupled
+(
+    const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
+    const scalar dt,
+    const scalar mass,
+    const scalar Re,
+    const scalar muc
+) const
+{
+    const scalar t = this->mesh().time().timeOutputValue();
+
+    const vector e(normalised(axisDirPtr_->value(t)));
+
+    if (mag(e) < VSMALL || !std::isfinite(mag(e)))
+    {
+        FatalIOErrorInFunction(this->coeffs())
+            << "The axisDir vector cannot be zero" << nl
+            << "axisDir: " << axisDirPtr_->value(t) << nl
+            << "time: " << t << nl
+            << error(FatalIOError);
+    }
+
+    const vector& X = p.position();
+    const vector& U = p.U();
+
+
+    // Calculate the radial distance between the parcel and rotational axis [m]
+    const vector axisPoint(axisPointPtr_->value(t));
+    const scalar s = e & (X - axisPoint);  // axial projection, s=e . (x-r0)
+    const vector Pax = axisPoint + s*e;    // closest point on axis P(x)
+    vector R = X - Pax;                    // radial vector
+    scalar r = mag(R);                     // radius, r =|R|
+
+
+    // Calculate the unit vector of the radial vector, Rhat
+    vector Rhat(vector::zero);
+    if (r > VSMALL)  // if parcel is farther than VSMALL from the rotation axis
+    {
+        Rhat = R/r;
+    }
+    else  // if parcel is effectively on the rotation axis, pick something
+    {
+        const vector n
+        (
+            (mag(e.x()) < 0.9) ? vector(1,0,0) : vector(0,1,0)
+        );
+        // When the parcel is on the axis, use the cross product to generate a
+        // perpendicular direction, ensuring a valid radial vector for further
+        // calculations
+        R = e^n;
+        r = mag(R);
+        Rhat = R/r;
+    }
+
+
+    // Calculate the velocity components of the parcel
+    const scalar Uaxial = (U & e);
+    const vector Utransverse(U - Uaxial*e);
+    const scalar Uradial = (Utransverse & Rhat);
+    const vector Utangent(Utransverse - Uradial*Rhat);
+
+
+    // Accumulate the acceleration components for the parcel
+    const scalar omega = omegaPtr_->value(t);
+    // Centripetal acceleration: directed inward toward the axis of rotation
+    // The negative sign ensures the acceleration points from the particle's
+    // position toward the axis.
+    vector a(-sqr(omega)*R);
+
+
+    if (tauRotationalPtr_)  // try to adjust for the desired tangential velocity
+    {
+        const scalar tau = tauRotationalPtr_->value(t);
+
+        if (tau < SMALL || !std::isfinite(tau))
+        {
+            FatalIOErrorInFunction(this->coeffs())
+                << "tauRotational cannot be zero or negative" << nl
+                << "tauRotational: " << tau << nl
+                << "time: " << t << nl
+                << error(FatalIOError);
+        }
+
+        const vector UtangentDesired(omega*(e^R));
+        a += (UtangentDesired - Utangent)/tau;
+    }
+
+
+    if (tauAxialPtr_ )  // try to adjust for desired axial velocity
+    {
+        if (!UaxialPtr_)
+        {
+            FatalIOErrorInFunction(this->coeffs())
+                << "If tauAxial is specified, Uaxial must also be specified."
+                << exit(FatalIOError);
+        }
+
+        const scalar tau = tauAxialPtr_->value(t);
+
+        if (tau < SMALL || !std::isfinite(tau))
+        {
+            FatalIOErrorInFunction(this->coeffs())
+                << "tauAxial cannot be zero or negative" << nl
+                << "tauAxial: " << tau << nl
+                << "time: " << t << nl
+                << error(FatalIOError);
+        }
+
+        const scalar UaxialDesired = UaxialPtr_->value(t);
+        a += ((UaxialDesired - Uaxial)/tau)*e;
+    }
+
+
+    // Apply radial damping to suppress radial oscillations
+    if (radialDampingPtr_)
+    {
+        const scalar damper = radialDampingPtr_->value(t);
+
+        if (damper < SMALL || !std::isfinite(damper))
+        {
+            FatalIOErrorInFunction(this->coeffs())
+                << "radialDamping cannot be zero or negative" << nl
+                << "radialDamping: " << damper << nl
+                << "time: " << t << nl
+                << error(FatalIOError);
+        }
+
+        a += (-damper*Uradial)*Rhat;
+    }
+
+
+    return forceSuSp(mass*a, Zero);  // (Sp, Su)
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Helical/HelicalForce.H

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::HelicalForce
+
+Group
+    grpLagrangianIntermediateForceSubModels
+
+Description
+    This class implements a helical force model for Lagrangian particles.
+    The force model applies a combination of centripetal, tangential, axial,
+    and radial damping forces to the particles to induce helical motion around
+    a specified axis.
+
+Usage
+    Minimal example by using \c constant/<cloud>/subModels/particleForces:
+    \verbatim
+    helical
+    {
+        // Mandatory entries
+        axisPoint      <Function1<vector>>;  // (0 0 0);
+        axisDir        <Function1<vector>>;  // (0 0 1);
+        omega          <Function1<scalar>>;  // 10;   [rad/s]
+
+        // Optional entries
+        Uaxial         <Function1<scalar>>;  // 1;    [m/s]
+        tauRotational  <Function1<scalar>>;  // 0.1;  [s]
+        tauAxial       <Function1<scalar>>;  // 0.1;  [s]
+        radialDamping  <Function1<scalar>>;  // 50;   [1/s]
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                  | Type | Reqd  | Deflt
+      axisPoint | A point on the rotation axis | F1\<vector\> | yes | -
+      axisDir   | Rotation axis (normalised)   | F1\<vector\> | yes | -
+      omega | Rotational speed of the parcel   | F1\<scalar\> | yes | -
+      Uaxial | Parcel axial speed in the rotation axis | F1\<scalar\> | no  | -
+      tauRotational| Time-scale constant for tangential-speed feedback <!--
+                --> | F1\<scalar\> | no | -
+      tauAxial | Time-scale constant for axial-speed feedback  <!--
+                --> | F1\<scalar\> | no | -
+      radialDamping | Radial damping coefficient | F1\<scalar\> | no | -
+    \endtable
+
+Note
+    - The \c axisDir vector must not be zero.
+    - If \c tauRotational is specified, it must be positive.
+    - If \c tauAxial is specified, it must be positive.
+    - If \c tauAxial is specified, \c Uaxial must be specified.
+    - If \c radialDamping is specified, it must be positive.
+
+SourceFiles
+    HelicalForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_HelicalForce_H
+#define Foam_HelicalForce_H
+
+#include "ParticleForce.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+template<class Type> class Function1;
+
+/*---------------------------------------------------------------------------*\
+                       Class HelicalForce Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class HelicalForce
+:
+    public ParticleForce<CloudType>
+{
+    // Private Data
+
+        //- Point on the rotation axis, r0
+        autoPtr<Function1<vector>> axisPointPtr_;
+
+        //- Rotation axis (normalised), e
+        autoPtr<Function1<vector>> axisDirPtr_;
+
+        //- Rotational speed of the parcel about the rotation axis [rad/s]
+        autoPtr<Function1<scalar>> omegaPtr_;
+
+        //- Axial speed of the parcel in the rotation axis direction [m/s]
+        autoPtr<Function1<scalar>> UaxialPtr_;
+
+        //- Time-scale constant for the tangential feedback acceleration to
+        //- keep the parcel's tangential velocity equal to the prescribed
+        //- rotational speed around the axis, if need be [s]
+        autoPtr<Function1<scalar>> tauRotationalPtr_;
+
+        //- Time-scale constant for the axial feedback acceleration to keep
+        //- the parcel's velocity along the axis equal to the prescribed axial
+        //- speed, if need be [s]
+        autoPtr<Function1<scalar>> tauAxialPtr_;
+
+        //- Radial damping coefficient to suppress radial oscillations [1/s]
+        //  e.g. 50 for strong damping
+        autoPtr<Function1<scalar>> radialDampingPtr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("helical");
+
+
+    // Constructors
+
+        //- Construct from cloud, mesh and dictionary
+        HelicalForce
+        (
+            CloudType& owner,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Copy construct
+        HelicalForce(const HelicalForce<CloudType>& df);
+
+        //- Construct and return a clone
+        virtual autoPtr<ParticleForce<CloudType>> clone() const
+        {
+            return autoPtr<ParticleForce<CloudType>>
+            (
+                new HelicalForce<CloudType>(*this)
+            );
+        }
+
+        //- No copy assignment
+        void operator=(const HelicalForce<CloudType>&) = delete;
+
+
+    //- Destructor
+    virtual ~HelicalForce();
+
+
+    // Member Functions
+
+    // Evaluation
+
+        //- Calculate the coupled force
+        virtual forceSuSp calcCoupled
+        (
+            const typename CloudType::parcelType& p,
+            const typename CloudType::parcelType::trackingData& td,
+            const scalar dt,
+            const scalar mass,
+            const scalar Re,
+            const scalar muc
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "HelicalForce.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //