Updates to the TAB modeling for distorted droplets

Resolves bug-reports:
   http://www.openfoam.org/mantisbt/view.php?id=1018
   http://www.openfoam.org/mantisbt/view.php?id=996
   http://www.openfoam.org/mantisbt/view.php?id=991

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "DistortedSphereDragForce.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::DistortedSphereDragForce<CloudType>::CdRe
+(
+    const scalar Re
+) const
+{
+    if (Re > 1000.0)
+    {
+        return 0.424*Re;
+    }
+    else
+    {
+        return 24.0*(1.0 + (1.0/6.0)*pow(Re, 2.0/3.0));
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DistortedSphereDragForce<CloudType>::DistortedSphereDragForce
+(
+    CloudType& owner,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, false)
+{}
+
+
+template<class CloudType>
+Foam::DistortedSphereDragForce<CloudType>::DistortedSphereDragForce
+(
+    const DistortedSphereDragForce<CloudType>& df
+)
+:
+    ParticleForce<CloudType>(df)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DistortedSphereDragForce<CloudType>::~DistortedSphereDragForce()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::forceSuSp Foam::DistortedSphereDragForce<CloudType>::calcCoupled
+(
+    const typename CloudType::parcelType& p,
+    const scalar dt,
+    const scalar mass,
+    const scalar Re,
+    const scalar muc
+) const
+{
+    forceSuSp value(vector::zero, 0.0);
+
+    // Limit the drop distortion to y=0 (sphere) and y=1 (disk)
+    scalar y = min(max(p.y(), 0), 1);
+
+    value.Sp() = mass*0.75*muc*CdRe(Re)*(1 + 2.632*y)/(p.rho()*sqr(p.d()));
+
+    return value;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H

@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::DistortedSphereDragForce
+
+Description
+    Drag model based on assumption of distorted spheres according to:
+
+    \verbatim
+        "Effects of Drop Drag and Breakup on Fuel Sprays"
+        Liu, A.B., Mather, D., Reitz, R.D.,
+        SAE Paper 930072,
+        SAE Transactions, Vol. 102, Section 3, Journal of Engines, 1993,
+        pp. 63-95
+    \endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DistortedSphereDragForce_H
+#define DistortedSphereDragForce_H
+
+#include "ParticleForce.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class DistortedSphereDragForce Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class DistortedSphereDragForce
+:
+    public ParticleForce<CloudType>
+{
+    // Private Member Functions
+
+        //- Drag coefficient multiplied by Reynolds number
+        scalar CdRe(const scalar Re) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("distortedSphereDrag");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        DistortedSphereDragForce
+        (
+            CloudType& owner,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct copy
+        DistortedSphereDragForce(const DistortedSphereDragForce<CloudType>& df);
+
+        //- Construct and return a clone
+        virtual autoPtr<ParticleForce<CloudType> > clone() const
+        {
+            return autoPtr<ParticleForce<CloudType> >
+            (
+                new DistortedSphereDragForce<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~DistortedSphereDragForce();
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Calculate the coupled force
+            virtual forceSuSp calcCoupled
+            (
+                const typename CloudType::parcelType& p,
+                const scalar dt,
+                const scalar mass,
+                const scalar Re,
+                const scalar muc
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DistortedSphereDragForce.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C

@@ -368,7 +368,7 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
 )
 {
     const scalar& TABCmu = td.cloud().breakup().TABCmu();
-    const scalar& TABWeCrit = td.cloud().breakup().TABWeCrit();
+    const scalar& TABtwoWeCrit = td.cloud().breakup().TABtwoWeCrit();
     const scalar& TABComega = td.cloud().breakup().TABComega();
 
     scalar r = 0.5*this->d();
@@ -385,27 +385,19 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
     {
         scalar omega = sqrt(omega2);
         scalar rhoc = this->rhoc();
-        scalar We = this->We(this->U(), r, rhoc, sigma_)/TABWeCrit;
+        scalar We = this->We(this->U(), r, rhoc, sigma_)/TABtwoWeCrit;
 
-        scalar y1 = this->y() - We;
+        // Initial values for y and yDot
-        scalar y2 = this->yDot()/omega;
+        scalar y0 = this->y() - We;
+        scalar yDot0 = this->yDot() + y0*rtd;
 
         // Update distortion parameters
         scalar c = cos(omega*dt);
         scalar s = sin(omega*dt);
         scalar e = exp(-rtd*dt);
-        y2 = (this->yDot() + y1*rtd)/omega;
 
-        this->y() = We + e*(y1*c + y2*s);
+        this->y() = We + e*(y0*c + (yDot0/omega)*s);
-        if (this->y() < 0)
+        this->yDot() = (We - this->y())*rtd + e*(yDot0*c - omega*y0*s);
-        {
-            this->y() = 0.0;
-            this->yDot() = 0.0;
-        }
-        else
-        {
-            this->yDot() = (We - this->y())*rtd + e*omega*(y2*c - y1*s);
-        }
     }
     else
     {

src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C

@@ -45,6 +45,7 @@ License
 #include "makeReactingParcelSurfaceFilmModels.H"
 
 // Spray
+#include "DistortedSphereDragForce.H"
 #include "makeSprayParcelAtomizationModels.H"
 #include "makeSprayParcelBreakupModels.H"
 
@@ -72,6 +73,7 @@ namespace Foam
     makeReactingParcelSurfaceFilmModels(basicSprayCloud);
 
     // Spray sub-models
+    makeParticleForceModelType(DistortedSphereDragForce, basicSprayCloud);
     makeSprayParcelAtomizationModels(basicSprayCloud);
     makeSprayParcelBreakupModels(basicSprayCloud);
 };

src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,7 @@ Foam::BreakupModel<CloudType>::BreakupModel
     yDot0_(0.0),
     TABComega_(0.0),
     TABCmu_(0.0),
-    TABWeCrit_(0.0)
+    TABtwoWeCrit_(0.0)
 {}
 
 
@@ -55,7 +55,7 @@ Foam::BreakupModel<CloudType>::BreakupModel
     yDot0_(bum.yDot0_),
     TABComega_(bum.TABComega_),
     TABCmu_(bum.TABCmu_),
-    TABWeCrit_(bum.TABWeCrit_)
+    TABtwoWeCrit_(bum.TABtwoWeCrit_)
 {}
 
 
@@ -70,20 +70,19 @@ Foam::BreakupModel<CloudType>::BreakupModel
 :
     CloudSubModelBase<CloudType>(owner, dict, typeName, type),
     solveOscillationEq_(solveOscillationEq),
-    y0_(0.0),
+    y0_(this->coeffDict().template lookupOrDefault<scalar>("y0", 0.0)),
-    yDot0_(0.0),
+    yDot0_(this->coeffDict().template lookupOrDefault<scalar>("yDot0", 0.0)),
-    TABComega_(0.0),
+    TABComega_(8),
-    TABCmu_(0.0),
+    TABCmu_(5),
-    TABWeCrit_(0.0)
+    TABtwoWeCrit_(12)
 {
-    if (solveOscillationEq_)
+    if (solveOscillationEq_ && dict.found("TABCoeffs"))
     {
         const dictionary coeffs(dict.subDict("TABCoeffs"));
-        y0_ = coeffs.template lookupOrDefault<scalar>("y0", 0.0);
+        coeffs.lookup("Comega") >> TABComega_;
-        yDot0_ = coeffs.template lookupOrDefault<scalar>("yDot0", 0.0);
+        coeffs.lookup("Cmu") >> TABCmu_;
-        TABComega_ = coeffs.template lookupOrDefault<scalar>("Comega", 8.0);
+        scalar WeCrit(readScalar(coeffs.lookup("WeCrit")));
-        TABCmu_ = coeffs.template lookupOrDefault<scalar>("Cmu", 10.0);
+        TABtwoWeCrit_ = 2*WeCrit;
-        TABWeCrit_ = coeffs.template lookupOrDefault<scalar>("WeCrit", 12.0);
     }
 }
 
@@ -160,4 +159,3 @@ bool Foam::BreakupModel<CloudType>::update
 #include "BreakupModelNew.C"
 
 // ************************************************************************* //
-

src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,9 +63,10 @@ protected:
 
         scalar y0_;
         scalar yDot0_;
+
         scalar TABComega_;
         scalar TABCmu_;
-        scalar TABWeCrit_;
+        scalar TABtwoWeCrit_;
 
 
 public:
@@ -153,9 +154,9 @@ public:
             return TABCmu_;
         }
 
-        inline const scalar& TABWeCrit() const
+        inline const scalar& TABtwoWeCrit() const
         {
-            return TABWeCrit_;
+            return TABtwoWeCrit_;
         }
 
 

src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,7 +116,7 @@ bool Foam::ETAB<CloudType>::update
         scalar romega = 1.0/omega;
 
         scalar We = rhoc*sqr(Urmag)*r/sigma;
-        scalar Wetmp = We/this->TABWeCrit_;
+        scalar Wetmp = We/this->TABtwoWeCrit_;
 
         scalar y1 = y - Wetmp;
         scalar y2 = yDot*romega;

src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,7 +128,7 @@ bool Foam::TAB<CloudType>::update
     {
         scalar omega = sqrt(omega2);
         scalar We = rhoc*sqr(Urmag)*r/sigma;
-        scalar Wetmp = We/this->TABWeCrit_;
+        scalar Wetmp = We/this->TABtwoWeCrit_;
 
         scalar y1 = y - Wetmp;
         scalar y2 = yDot/omega;