diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
index f043febe81..1a71c63e05 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,14 +36,13 @@ Foam::scalar Foam::DistortedSphereDragForce<CloudType>::CdRe
     const scalar Re
 ) const
 {
+    // (AOB:Eq. 35; LMR:Eq. 9)
     if (Re > 1000.0)
     {
         return 0.424*Re;
     }
-    else
-    {
-        return 24.0*(1.0 + (1.0/6.0)*pow(Re, 2.0/3.0));
-    }
+
+    return 24.0*(1.0 + (1.0/6.0)*pow(Re, 2.0/3.0));
 }
 
 
@@ -70,13 +70,6 @@ Foam::DistortedSphereDragForce<CloudType>::DistortedSphereDragForce
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::DistortedSphereDragForce<CloudType>::~DistortedSphereDragForce()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -90,14 +83,16 @@ Foam::forceSuSp Foam::DistortedSphereDragForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    forceSuSp value(Zero);
-
     // Limit the drop distortion to y=0 (sphere) and y=1 (disk)
-    scalar y = min(max(p.y(), 0), 1);
+    const scalar y = min(max(p.y(), scalar(0)), scalar(1));
 
-    value.Sp() = mass*0.75*muc*CdRe(Re)*(1 + 2.632*y)/(p.rho()*sqr(p.d()));
-
-    return value;
+    // (LMR:Eq. 10)
+    return
+        forceSuSp
+        (
+            Zero,
+            mass*0.75*muc*CdRe(Re)*(1.0 + 2.632*y)/(p.rho()*sqr(p.d()))
+        );
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
index 7cdfba6116..386ed8ba15 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
@@ -126,7 +126,7 @@ SourceFiles
 namespace Foam
 {
 /*---------------------------------------------------------------------------*\
-                       Class DistortedSphereDragForce Declaration
+                   Class DistortedSphereDragForce Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
@@ -168,9 +168,12 @@ public:
             );
         }
 
+        //- No copy assignment
+        void operator=(const DistortedSphereDragForce<CloudType>&) = delete;
+
 
     //- Destructor
-    virtual ~DistortedSphereDragForce();
+    virtual ~DistortedSphereDragForce() = default;
 
 
     // Member Functions
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C
index d943cc0e79..7ac32aaddc 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C
@@ -31,19 +31,16 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::scalar Foam::ErgunWenYuDragForce<CloudType>::CdRe
-(
-    const scalar Re
-) const
+Foam::scalar Foam::ErgunWenYuDragForce<CloudType>::CdRe(const scalar alphacRe)
+const
 {
-    if (Re > 1000.0)
+    // (ZZB:Eq. 14, GLSLR:Table 3)
+    if (alphacRe < 1000.0)
     {
-        return 0.44*Re;
-    }
-    else
-    {
-        return 24.0*(1.0 + 0.15*pow(Re, 0.687));
+        return 24.0*(1.0 + 0.15*pow(alphacRe, 0.687));
     }
+
+    return 0.44*alphacRe;
 }
 
 
@@ -85,13 +82,6 @@ Foam::ErgunWenYuDragForce<CloudType>::ErgunWenYuDragForce
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::ErgunWenYuDragForce<CloudType>::~ErgunWenYuDragForce()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -105,7 +95,7 @@ Foam::forceSuSp Foam::ErgunWenYuDragForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    scalar alphac(alphac_[p.cell()]);
+    const scalar alphac = alphac_[p.cell()];
 
     if (alphac < 0.8)
     {
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H
index 41cc8eb897..e50ce32a03 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H
@@ -168,8 +168,8 @@ class ErgunWenYuDragForce
 
     // Private Member Functions
 
-        //- Drag coefficient multiplied by Reynolds number
-        scalar CdRe(const scalar Re) const;
+        //- Drag coefficient multiplied by volume fraction and Reynolds number
+        scalar CdRe(const scalar alphacRe) const;
 
 
 public:
@@ -200,9 +200,12 @@ public:
             );
         }
 
+        //- No copy assignment
+        void operator=(const ErgunWenYuDragForce<CloudType>&) = delete;
+
 
     //- Destructor
-    virtual ~ErgunWenYuDragForce();
+    virtual ~ErgunWenYuDragForce() = default;
 
 
     // Member Functions
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C
index 04237685ea..fccb825f1b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +33,7 @@ License
 template<class CloudType>
 Foam::scalar Foam::NonSphereDragForce<CloudType>::CdRe(const scalar Re) const
 {
+    // (HL:Eq. 10-11)
     return 24.0*(1.0 + a_*pow(Re, b_)) + Re*c_/(1 + d_/(Re + ROOTVSMALL));
 }
 
@@ -78,13 +80,6 @@ Foam::NonSphereDragForce<CloudType>::NonSphereDragForce
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::NonSphereDragForce<CloudType>::~NonSphereDragForce()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -98,11 +93,7 @@ Foam::forceSuSp Foam::NonSphereDragForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    forceSuSp value(Zero);
-
-    value.Sp() = mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d()));
-
-    return value;
+    return forceSuSp(Zero, mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d())));
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H
index 2d5a791232..b87ef89d4f 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H
@@ -165,7 +165,7 @@ protected:
     // Protected Data
 
         //- Ratio of surface of sphere having same volume as particle to
-        //  actual surface area of particle (0 < phi <= 1)
+        //- actual surface area of particle (0 < phi <= 1)
         scalar phi_;
 
 
@@ -214,9 +214,12 @@ public:
             );
         }
 
+        //- No copy assignment
+        void operator=(const NonSphereDragForce<CloudType>&) = delete;
+
 
     //- Destructor
-    virtual ~NonSphereDragForce();
+    virtual ~NonSphereDragForce() = default;
 
 
     // Member Functions
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C
index 4f5d5ca8d5..14d96101d2 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C
@@ -66,13 +66,6 @@ Foam::PlessisMasliyahDragForce<CloudType>::PlessisMasliyahDragForce
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::PlessisMasliyahDragForce<CloudType>::~PlessisMasliyahDragForce()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -86,11 +79,11 @@ Foam::forceSuSp Foam::PlessisMasliyahDragForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    scalar alphac(alphac_[p.cell()]);
+    const scalar alphac = alphac_[p.cell()];
 
-    scalar cbrtAlphap(cbrt(1.0 - alphac));
+    const scalar cbrtAlphap = cbrt(1.0 - alphac);
 
-    scalar A =
+    const scalar A =
         26.8*pow3(alphac)
        /(
             sqr(cbrtAlphap)
@@ -99,9 +92,8 @@ Foam::forceSuSp Foam::PlessisMasliyahDragForce<CloudType>::calcCoupled
           + SMALL
         );
 
-    scalar B =
-        sqr(alphac)
-       /sqr(1.0 - sqr(cbrtAlphap));
+    // (P:Eq. 36)
+    const scalar B = sqr(alphac)/sqr(1.0 - sqr(cbrtAlphap));
 
     return forceSuSp
     (
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
index da24387aea..89254663a9 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
@@ -168,9 +168,12 @@ public:
             );
         }
 
+        //- No copy assignment
+        void operator=(const PlessisMasliyahDragForce<CloudType>&) = delete;
+
 
     //- Destructor
-    virtual ~PlessisMasliyahDragForce();
+    virtual ~PlessisMasliyahDragForce() = default;
 
 
     // Member Functions
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C
index 10fd75711f..4b14d1b18b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,14 +33,13 @@ License
 template<class CloudType>
 Foam::scalar Foam::SphereDragForce<CloudType>::CdRe(const scalar Re) const
 {
+    // (AOB:Eq. 35)
     if (Re > 1000.0)
     {
         return 0.424*Re;
     }
-    else
-    {
-        return 24.0*(1.0 + 1.0/6.0*pow(Re, 2.0/3.0));
-    }
+
+    return 24.0*(1.0 + (1.0/6.0)*pow(Re, 2.0/3.0));
 }
 
 
@@ -67,13 +67,6 @@ Foam::SphereDragForce<CloudType>::SphereDragForce
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::SphereDragForce<CloudType>::~SphereDragForce()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -87,11 +80,8 @@ Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    forceSuSp value(Zero);
-
-    value.Sp() = mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d()));
-
-    return value;
+    // (AOB:Eq. 34)
+    return forceSuSp(Zero, mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d())));
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H
index 1fd0bb0772..cb71dfd0f2 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H
@@ -169,9 +169,12 @@ public:
             );
         }
 
+        //- No copy assignment
+        void operator=(const SphereDragForce<CloudType>&) = delete;
+
 
     //- Destructor
-    virtual ~SphereDragForce();
+    virtual ~SphereDragForce() = default;
 
 
     // Member Functions
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C
index 70f5cbb6e0..06bf9fbe85 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C
@@ -31,16 +31,15 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::scalar Foam::WenYuDragForce<CloudType>::CdRe(const scalar Re) const
+Foam::scalar Foam::WenYuDragForce<CloudType>::CdRe(const scalar alphacRe) const
 {
-    if (Re > 1000.0)
+    // (ZZB:Eq. 14, GLSLR:Table 3)
+    if (alphacRe < 1000.0)
     {
-        return 0.44*Re;
-    }
-    else
-    {
-        return 24.0*(1.0 + 0.15*pow(Re, 0.687));
+        return 24.0*(1.0 + 0.15*pow(alphacRe, 0.687));
     }
+
+    return 0.44*alphacRe;
 }
 
 
@@ -82,13 +81,6 @@ Foam::WenYuDragForce<CloudType>::WenYuDragForce
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::WenYuDragForce<CloudType>::~WenYuDragForce()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -102,7 +94,7 @@ Foam::forceSuSp Foam::WenYuDragForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    scalar alphac(alphac_[p.cell()]);
+    const scalar alphac = alphac_[p.cell()];
 
     return forceSuSp
     (
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
index dd3dd88575..cc930b14b6 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
@@ -158,8 +158,8 @@ class WenYuDragForce
 
     // Private Member Functions
 
-        //- Drag coefficient multiplied by Reynolds number
-        scalar CdRe(const scalar Re) const;
+        //- Drag coefficient multiplied by volume fraction and Reynolds number
+        scalar CdRe(const scalar alphacRe) const;
 
 
 public:
@@ -190,9 +190,12 @@ public:
             );
         }
 
+        //- No copy assignment
+        void operator=(const WenYuDragForce<CloudType>&) = delete;
+
 
     //- Destructor
-    virtual ~WenYuDragForce();
+    virtual ~WenYuDragForce() = default;
 
 
     // Member Functions