View Issue Details

IDProjectCategoryView StatusLast Update
0000993OpenFOAMBugpublic2018-07-10 11:21
Reporterdkxls Assigned Touser2 
PrioritynormalSeverityminorReproducibilityN/A
Status resolvedResolutionfixed 
PlatformLinux x86_64OSopenSUSEOS Version12.2
Summary0000993: [BreakupModel]: Inconsistent or no reading of TAB/ETAB coefficients leading to wrong sphere deviation
DescriptionIn the present implementation, the TAB coefficients in BreakupModel are initialized to 0.0 and only read/set to proper values if "solveOscillationEq" is true. In case the coefficients are read, then only ones in the BreakupModel constructor.
This lead to the following issues:

1. The TAB model sets it's coefficients to the ones from the BreakupModel. If "solveOscillationEq" is false, the TAB constructor sets it to true, but coefficients are still 0.0 as the reading takes place in the BreakupModel constructor prior to setting solveOscillationEq to true.

2. The ETAB model reads a own set of coefficients, which could lead to different coefficients in the ETAB and the BreakupModel (and hence in the solveTABEq function). Besides, also here it is possible that the BreakupModel coefficients are not read and left at 0.0 as in the TAB model.

These issues are not exactly critical as they can be circumvented by setting the appropriate switches correctly, but this require the user to exactly know which switches to set where.
I only understood it after examining the source code - the tutorial was rather confusing!
Additional InformationThe oscillation equation should probably be moved from the SprayParcel to the BreakupModel, as it is an integral part of it and requires also the coefficients from there.
Although, the DistortedSphereDragForce drag model (code provided in bug #991) requires the oscillation equation to be solved.
Tagsbreakup, spray

Activities

dkxls

2013-09-04 19:37

reporter  

0001-Fix-inconsistent-breakup-coefficient-reading.patch (5,035 bytes)   
diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
index f714be3..9439fc4 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
@@ -64,27 +64,42 @@ Foam::BreakupModel<CloudType>::BreakupModel
 (
     const dictionary& dict,
     CloudType& owner,
-    const word& type
+    const word& type,
+    const Switch solveOscillationEq = false
 )
 :
     SubModelBase<CloudType>(owner, dict, typeName, type),
-    solveOscillationEq_(this->coeffDict().lookup("solveOscillationEq")),
+    solveOscillationEq_(solveOscillationEq),
     y0_(0.0),
     yDot0_(0.0),
-    TABComega_(0.0),
-    TABCmu_(0.0),
-    TABWeCrit_(0.0)
+    TABComega_(8.0),
+    TABCmu_(10.0),
+    TABWeCrit_(12.0)
 {
+    // Give the option to solve the oscillation equation, even if not required 
+    // by the breakup model
+    if (!solveOscillationEq_)
+    {
+        this->coeffDict().readIfPresent
+        (
+            "solveOscillationEq",
+            solveOscillationEq_
+        );
+    }
+
+    // Check for non-default TAB coefficients
     if (solveOscillationEq_)
     {
         const dictionary TABcoeffsDict(dict.subDict("TABCoeffs"));
-        y0_ = TABcoeffsDict.template lookupOrDefault<scalar>("y0", 0.0);
-        yDot0_ = TABcoeffsDict.template lookupOrDefault<scalar>("yDot0", 0.0);
-        TABComega_ =
-            TABcoeffsDict.template lookupOrDefault<scalar>("Comega", 8.0);
-        TABCmu_ = TABcoeffsDict.template lookupOrDefault<scalar>("Cmu", 10.0);
-        TABWeCrit_ =
-            TABcoeffsDict.template lookupOrDefault<scalar>("WeCrit", 12.0);
+        Switch defaultCoeffs = TABcoeffsDict.lookup("defaultCoeffs");
+        if (!defaultCoeffs)
+        {
+            TABcoeffsDict.lookup("y0") >> y0_;
+            TABcoeffsDict.lookup("yDot0") >> yDot0_;
+            TABcoeffsDict.lookup("Comega") >> TABComega_;
+            TABcoeffsDict.lookup("Cmu") >> TABCmu_;
+            TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
+        }
     }
 }
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
index e0672ca..09cbe36 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
@@ -97,7 +97,8 @@ public:
         (
             const dictionary& dict,
             CloudType& owner,
-            const word& type
+            const word& type,
+            const Switch solveOscillationEq
         );
 
         //- Construct copy
diff --git a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
index fb48040..7f08e47 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
@@ -34,7 +34,7 @@ Foam::ETAB<CloudType>::ETAB
     CloudType& owner
 )
 :
-    BreakupModel<CloudType>(dict, owner, typeName),
+    BreakupModel<CloudType>(dict, owner, typeName, true),
     Cmu_(10.0),
     Comega_(8.0),
     k1_(0.2),
@@ -55,15 +55,6 @@ Foam::ETAB<CloudType>::ETAB
 
     scalar k21 = k2_/k1_;
     AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow4(WeTransition_);
-
-    if (!BreakupModel<CloudType>::solveOscillationEq_)
-    {
-        Info<< "Warning: solveOscillationEq is set to "
-            << BreakupModel<CloudType>::solveOscillationEq_ << nl
-            << " Setting it to true in order for the ETAB model to work."
-            << endl;
-        BreakupModel<CloudType>::solveOscillationEq_ = true;
-    }
 }
 
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
index 0fe325f..53f416f 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
@@ -34,7 +34,7 @@ Foam::TAB<CloudType>::TAB
     CloudType& owner
 )
 :
-    BreakupModel<CloudType>(dict, owner, typeName),
+    BreakupModel<CloudType>(dict, owner, typeName, true),
     Cmu_(BreakupModel<CloudType>::TABCmu_),
     Comega_(BreakupModel<CloudType>::TABComega_),
     WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
@@ -52,16 +52,6 @@ Foam::TAB<CloudType>::TAB
             (1.0 - exp(-xx)*(1.0 + xx + sqr(xx)/2.0 + pow3(xx)/6.0))*rrd100;
     }
 
-    if (!BreakupModel<CloudType>::solveOscillationEq_)
-    {
-        WarningIn("Foam::TAB<CloudType>::TAB(const dictionary&, CloudType&)")
-            << "solveOscillationEq is set to "
-            << BreakupModel<CloudType>::solveOscillationEq_ << nl
-            << " Setting it to true in order for the TAB model to work."
-            << endl;
-        BreakupModel<CloudType>::solveOscillationEq_ = true;
-    }
-
     if (SMDCalcMethod_ == "method1")
     {
         SMDMethod_ = method1;

dkxls

2013-09-04 19:37

reporter  

0001-Fix-inconsistent-ETAB-coefficient-reading.patch (1,410 bytes)   
From c0e01d022e1923d080872c812408e5f78165ebe2 Mon Sep 17 00:00:00 2001
From: Armin Wehrfritz <armin.wehrfritz@aalto.fi>
Date: Wed, 4 Sep 2013 21:32:39 +0300
Subject: [PATCH] Fix inconsistent ETAB coefficient reading

---
 src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C |    9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
index 7f08e47..24a33e1 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
@@ -35,21 +35,18 @@ Foam::ETAB<CloudType>::ETAB
 )
 :
     BreakupModel<CloudType>(dict, owner, typeName, true),
-    Cmu_(10.0),
-    Comega_(8.0),
+    Cmu_(BreakupModel<CloudType>::TABCmu_),
+    Comega_(BreakupModel<CloudType>::TABComega_),
     k1_(0.2),
     k2_(0.2),
-    WeCrit_(12.0),
+    WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
     WeTransition_(100.0),
     AWe_(0.0)
 {
     if (!this->defaultCoeffs(true))
     {
-        this->coeffDict().lookup("Cmu") >> Cmu_;
-        this->coeffDict().lookup("Comega") >> Comega_;
         this->coeffDict().lookup("k1") >> k1_;
         this->coeffDict().lookup("k2") >> k2_;
-        this->coeffDict().lookup("WeCrit") >> WeCrit_;
         this->coeffDict().lookup("WeTransition") >> WeTransition_;
     }
 
-- 
1.7.10.4

dkxls

2013-09-04 19:46

reporter   ~0002461

Last edited: 2013-09-05 00:24

The patches I uploaded fix the two issues and should guarantee consistent coefficients for the breakup models.

dkxls

2013-09-04 20:38

reporter  

0001-Enhance-TAB-coefficient-reading-and-add-enableOscill.patch (4,043 bytes)   
From 48872ed9df56504a8b3a60eff752ff758cc75049 Mon Sep 17 00:00:00 2001
From: Armin Wehrfritz <armin.wehrfritz@aalto.fi>
Date: Wed, 4 Sep 2013 22:26:42 +0300
Subject: [PATCH] Enhance TAB coefficient reading and add
 enableOscillationEq()

---
 .../BreakupModel/BreakupModel/BreakupModel.C       |   54 ++++++++++++++++----
 .../BreakupModel/BreakupModel/BreakupModel.H       |   10 ++++
 2 files changed, 54 insertions(+), 10 deletions(-)

diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
index 9439fc4..f5fa020 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
@@ -25,6 +25,30 @@ License
 
 #include "BreakupModel.H"
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::BreakupModel<CloudType>::readTABCoeffs(const bool printMsg)
+{
+    const dictionary TABcoeffsDict(this->dict().subDict("TABCoeffs"));
+    bool def = TABcoeffsDict.lookupOrDefault<bool>("defaultCoeffs", false);
+    if (!def)
+    {
+        TABcoeffsDict.lookup("y0") >> y0_;
+        TABcoeffsDict.lookup("yDot0") >> yDot0_;
+        TABcoeffsDict.lookup("Comega") >> TABComega_;
+        TABcoeffsDict.lookup("Cmu") >> TABCmu_;
+        TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
+    }
+    if (printMsg && def)
+    {
+        Info<< incrIndent;
+        Info<< indent << "Employing default TAB coefficients" << endl;
+        Info<< decrIndent;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -90,16 +114,7 @@ Foam::BreakupModel<CloudType>::BreakupModel
     // Check for non-default TAB coefficients
     if (solveOscillationEq_)
     {
-        const dictionary TABcoeffsDict(dict.subDict("TABCoeffs"));
-        Switch defaultCoeffs = TABcoeffsDict.lookup("defaultCoeffs");
-        if (!defaultCoeffs)
-        {
-            TABcoeffsDict.lookup("y0") >> y0_;
-            TABcoeffsDict.lookup("yDot0") >> yDot0_;
-            TABcoeffsDict.lookup("Comega") >> TABComega_;
-            TABcoeffsDict.lookup("Cmu") >> TABCmu_;
-            TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
-        }
+        readTABCoeffs(true);
     }
 }
 
@@ -114,6 +129,25 @@ Foam::BreakupModel<CloudType>::~BreakupModel()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
+void Foam::BreakupModel<CloudType>::enableOscillationEq(const bool printMsg)
+{
+    if (!solveOscillationEq_)
+    {
+        solveOscillationEq_ = true;
+
+        if (printMsg)
+        {
+            Info<< incrIndent;
+            Info<< indent << "Enable oscillation equation solution" << endl;
+            Info<< decrIndent;
+        }
+
+        readTABCoeffs(printMsg);
+    }
+}
+
+
+template<class CloudType>
 bool Foam::BreakupModel<CloudType>::update
 (
     const scalar dt,
diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
index 09cbe36..ac54b73 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
@@ -55,6 +55,12 @@ class BreakupModel
     public SubModelBase<CloudType>
 {
 
+    // Private Member Functions
+
+        //- Read TAB coefficients
+        void readTABCoeffs(const bool printMsg);
+
+
 protected:
 
     // Protected data
@@ -161,6 +167,10 @@ public:
 
     // Member Functions
 
+        //- Enable the solution of the oscillation equation and check that the 
+        //  TAB coefficients are set correctly
+        void enableOscillationEq(const bool printMsg);
+
         //- Update the parcel properties and return true if a child parcel
         //  should be added
         virtual bool update
-- 
1.7.10.4

dkxls

2013-09-04 20:44

reporter   ~0002462

The last patch enhances the TAB coefficient reading and adds a function 'enableOscillationEq()' that makes it possible to enable the solution of the oscillation equation and checks that the coefficients are set correctly.

The last patch needs to be applied atop the other two patches.
Sorry for the bit messy patch arrangement, I therefor also uploaded all changed source files. Hope this improves the readability.

The code is tested and works as expected.

dkxls

2013-09-04 20:45

reporter  

BreakupModel.H (6,983 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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::BreakupModel

Description
    Templated break-up model class

SourceFiles
    BreakupModel.C
    BreakupModelNew.C

\*---------------------------------------------------------------------------*/

#ifndef BreakupModel_H
#define BreakupModel_H

#include "IOdictionary.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                       Class BreakupModel Declaration
\*---------------------------------------------------------------------------*/

template<class CloudType>
class BreakupModel
:
    public SubModelBase<CloudType>
{

    // Private Member Functions

        //- Read TAB coefficients
        void readTABCoeffs(const bool printMsg);


protected:

    // Protected data

        Switch solveOscillationEq_;

        scalar y0_;
        scalar yDot0_;
        scalar TABComega_;
        scalar TABCmu_;
        scalar TABWeCrit_;


public:

    //- Runtime type information
    TypeName("breakupModel");

    //- Declare runtime constructor selection table
    declareRunTimeSelectionTable
    (
        autoPtr,
        BreakupModel,
        dictionary,
        (
            const dictionary& dict,
            CloudType& owner
        ),
        (dict, owner)
    );


    // Constructors

        //- Construct null from owner
        BreakupModel(CloudType& owner);

        //- Construct from dictionary
        BreakupModel
        (
            const dictionary& dict,
            CloudType& owner,
            const word& type,
            const Switch solveOscillationEq
        );

        //- Construct copy
        BreakupModel(const BreakupModel<CloudType>& bum);

        //- Construct and return a clone
        virtual autoPtr<BreakupModel<CloudType> > clone() const
        {
            return autoPtr<BreakupModel<CloudType> >
            (
                new BreakupModel<CloudType>(*this)
            );
        }


    //- Destructor
    virtual ~BreakupModel();


    //- Selector
    static autoPtr<BreakupModel<CloudType> > New
    (
        const dictionary& dict,
        CloudType& owner
    );


    // Access

        inline const Switch& solveOscillationEq() const
        {
            return solveOscillationEq_;
        }

        inline const scalar& y0() const
        {
            return y0_;
        }

        inline const scalar& yDot0() const
        {
            return yDot0_;
        }

        inline const scalar& TABComega() const
        {
            return TABComega_;
        }

        inline const scalar& TABCmu() const
        {
            return TABCmu_;
        }

        inline const scalar& TABWeCrit() const
        {
            return TABWeCrit_;
        }


    // Member Functions

        //- Enable the solution of the oscillation equation and check that the 
        //  TAB coefficients are set correctly
        void enableOscillationEq(const bool printMsg);

        //- Update the parcel properties and return true if a child parcel
        //  should be added
        virtual bool update
        (
            const scalar dt,
            const vector& g,
            scalar& d,
            scalar& tc,
            scalar& ms,
            scalar& nParticle,
            scalar& KHindex,
            scalar& y,
            scalar& yDot,
            const scalar d0,
            const scalar rho,
            const scalar mu,
            const scalar sigma,
            const vector& U,
            const scalar rhoc,
            const scalar muc,
            const vector& Urel,
            const scalar Urmag,
            const scalar tMom,
            scalar& dChild,
            scalar& massChild
        );
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#define makeBreakupModel(CloudType)                                           \
                                                                              \
    typedef CloudType::sprayCloudType sprayCloudType;                         \
    defineNamedTemplateTypeNameAndDebug                                       \
    (                                                                         \
        BreakupModel<sprayCloudType>,                                         \
        0                                                                     \
    );                                                                        \
    defineTemplateRunTimeSelectionTable                                       \
    (                                                                         \
        BreakupModel<sprayCloudType>,                                         \
        dictionary                                                            \
    );


#define makeBreakupModelType(SS, CloudType)                                   \
                                                                              \
    typedef CloudType::sprayCloudType sprayCloudType;                         \
    defineNamedTemplateTypeNameAndDebug(SS<sprayCloudType>, 0);               \
                                                                              \
    BreakupModel<sprayCloudType>::                                            \
        adddictionaryConstructorToTable<SS<sprayCloudType> >                  \
            add##SS##CloudType##sprayCloudType##ConstructorToTable_;


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "BreakupModel.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
BreakupModel.H (6,983 bytes)   

dkxls

2013-09-04 20:45

reporter  

BreakupModel.C (5,451 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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 "BreakupModel.H"

// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //

template<class CloudType>
void Foam::BreakupModel<CloudType>::readTABCoeffs(const bool printMsg)
{
    const dictionary TABcoeffsDict(this->dict().subDict("TABCoeffs"));
    bool def = TABcoeffsDict.lookupOrDefault<bool>("defaultCoeffs", false);
    if (!def)
    {
        TABcoeffsDict.lookup("y0") >> y0_;
        TABcoeffsDict.lookup("yDot0") >> yDot0_;
        TABcoeffsDict.lookup("Comega") >> TABComega_;
        TABcoeffsDict.lookup("Cmu") >> TABCmu_;
        TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
    }
    if (printMsg && def)
    {
        Info<< incrIndent;
        Info<< indent << "Employing default TAB coefficients" << endl;
        Info<< decrIndent;
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CloudType>
Foam::BreakupModel<CloudType>::BreakupModel
(
    CloudType& owner
)
:
    SubModelBase<CloudType>(owner),
    solveOscillationEq_(false),
    y0_(0.0),
    yDot0_(0.0),
    TABComega_(0.0),
    TABCmu_(0.0),
    TABWeCrit_(0.0)
{}


template<class CloudType>
Foam::BreakupModel<CloudType>::BreakupModel
(
    const BreakupModel<CloudType>& bum
)
:
    SubModelBase<CloudType>(bum),
    solveOscillationEq_(bum.solveOscillationEq_),
    y0_(bum.y0_),
    yDot0_(bum.yDot0_),
    TABComega_(bum.TABComega_),
    TABCmu_(bum.TABCmu_),
    TABWeCrit_(bum.TABWeCrit_)
{}


template<class CloudType>
Foam::BreakupModel<CloudType>::BreakupModel
(
    const dictionary& dict,
    CloudType& owner,
    const word& type,
    const Switch solveOscillationEq = false
)
:
    SubModelBase<CloudType>(owner, dict, typeName, type),
    solveOscillationEq_(solveOscillationEq),
    y0_(0.0),
    yDot0_(0.0),
    TABComega_(8.0),
    TABCmu_(10.0),
    TABWeCrit_(12.0)
{
    // Give the option to solve the oscillation equation, even if not required 
    // by the breakup model
    if (!solveOscillationEq_)
    {
        this->coeffDict().readIfPresent
        (
            "solveOscillationEq",
            solveOscillationEq_
        );
    }

    // Check for non-default TAB coefficients
    if (solveOscillationEq_)
    {
        readTABCoeffs(true);
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CloudType>
Foam::BreakupModel<CloudType>::~BreakupModel()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class CloudType>
void Foam::BreakupModel<CloudType>::enableOscillationEq(const bool printMsg)
{
    if (!solveOscillationEq_)
    {
        solveOscillationEq_ = true;

        if (printMsg)
        {
            Info<< incrIndent;
            Info<< indent << "Enable oscillation equation solution" << endl;
            Info<< decrIndent;
        }

        readTABCoeffs(printMsg);
    }
}


template<class CloudType>
bool Foam::BreakupModel<CloudType>::update
(
    const scalar dt,
    const vector& g,
    scalar& d,
    scalar& tc,
    scalar& ms,
    scalar& nParticle,
    scalar& KHindex,
    scalar& y,
    scalar& yDot,
    const scalar d0,
    const scalar rho,
    const scalar mu,
    const scalar sigma,
    const vector& U,
    const scalar rhoc,
    const scalar muc,
    const vector& Urel,
    const scalar Urmag,
    const scalar tMom,
    scalar& dChild,
    scalar& massChild
)
{
    notImplemented
    (
        "bool Foam::BreakupModel<CloudType>::update"
        "("
            "const scalar, "
            "const vector&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "const scalar, "
            "const scalar, "
            "const scalar, "
            "const scalar, "
            "const vector&, "
            "const scalar, "
            "const scalar, "
            "const vector&, "
            "const scalar, "
            "const scalar, "
            "scalar&, "
            "scalar&"
        ");"
    );

    return false;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "BreakupModelNew.C"

// ************************************************************************* //

BreakupModel.C (5,451 bytes)   

dkxls

2013-09-04 20:45

reporter  

ETAB.C (5,699 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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 "ETAB.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CloudType>
Foam::ETAB<CloudType>::ETAB
(
    const dictionary& dict,
    CloudType& owner
)
:
    BreakupModel<CloudType>(dict, owner, typeName, true),
    Cmu_(BreakupModel<CloudType>::TABCmu_),
    Comega_(BreakupModel<CloudType>::TABComega_),
    k1_(0.2),
    k2_(0.2),
    WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
    WeTransition_(100.0),
    AWe_(0.0)
{
    if (!this->defaultCoeffs(true))
    {
        this->coeffDict().lookup("k1") >> k1_;
        this->coeffDict().lookup("k2") >> k2_;
        this->coeffDict().lookup("WeTransition") >> WeTransition_;
    }

    scalar k21 = k2_/k1_;
    AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow4(WeTransition_);
}


template<class CloudType>
Foam::ETAB<CloudType>::ETAB(const ETAB<CloudType>& bum)
:
    BreakupModel<CloudType>(bum),
    Cmu_(bum.Cmu_),
    Comega_(bum.Comega_),
    k1_(bum.k1_),
    k2_(bum.k2_),
    WeCrit_(bum.WeCrit_),
    WeTransition_(bum.WeTransition_),
    AWe_(bum.AWe_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CloudType>
Foam::ETAB<CloudType>::~ETAB()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class CloudType>
bool Foam::ETAB<CloudType>::update
(
    const scalar dt,
    const vector& g,
    scalar& d,
    scalar& tc,
    scalar& ms,
    scalar& nParticle,
    scalar& KHindex,
    scalar& y,
    scalar& yDot,
    const scalar d0,
    const scalar rho,
    const scalar mu,
    const scalar sigma,
    const vector& U,
    const scalar rhoc,
    const scalar muc,
    const vector& Urel,
    const scalar Urmag,
    const scalar tMom,
    scalar& dChild,
    scalar& massChild
)
{
    scalar r  = 0.5*d;
    scalar r2 = r*r;
    scalar r3 = r*r2;

    scalar semiMass = nParticle*pow3(d);

    // inverse of characteristic viscous damping time
    scalar rtd = 0.5*Cmu_*mu/(rho*r2);

    // oscillation frequency (squared)
    scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;

    if (omega2 > 0)
    {
        scalar omega = sqrt(omega2);
        scalar romega = 1.0/omega;

        scalar We = rhoc*sqr(Urmag)*r/sigma;
        scalar Wetmp = We/WeCrit_;

        scalar y1 = y - Wetmp;
        scalar y2 = yDot*romega;

        scalar a = sqrt(y1*y1 + y2*y2);

        // scotty we may have break-up
        if (a + Wetmp > 1.0)
        {
            scalar phic = y1/a;

            // constrain phic within -1 to 1
            phic = max(min(phic, 1), -1);

            scalar phit = acos(phic);
            scalar phi = phit;
            scalar quad = -y2/a;
            if (quad < 0)
            {
                phi = constant::mathematical::twoPi - phit;
            }

            scalar tb = 0;

            if (mag(y) < 1.0)
            {
                scalar theta = acos((1.0 - Wetmp)/a);

                if (theta < phi)
                {
                    if (constant::mathematical::twoPi - theta >= phi)
                    {
                        theta = -theta;
                    }
                    theta += constant::mathematical::twoPi;
                }
                tb = (theta - phi)*romega;

                // breakup occurs
                if (dt > tb)
                {
                    y = 1.0;
                    yDot = -a*omega*sin(omega*tb + phi);
                }
            }

            // update droplet size
            if (dt > tb)
            {
                scalar sqrtWe = AWe_*pow4(We) + 1.0;
                scalar Kbr = k1_*omega*sqrtWe;

                if (We > WeTransition_)
                {
                    sqrtWe = sqrt(We);
                    Kbr =k2_*omega*sqrtWe;
                }

                scalar rWetmp = 1.0/Wetmp;
                scalar cosdtbu = max(-1.0, min(1.0, 1.0 - rWetmp));
                scalar dtbu = romega*acos(cosdtbu);
                scalar decay = exp(-Kbr*dtbu);

                scalar rNew = decay*r;
                if (rNew < r)
                {
                    d = 2.0*rNew;
                    y = 0.0;
                    yDot = 0.0;
                }
            }
        }
    }
    else
    {
        // reset droplet distortion parameters
        y = 0;
        yDot = 0;
    }

    // update the nParticle count to conserve mass
    nParticle = semiMass/pow3(d);

    // Do not add child parcel
    return false;
}


// ************************************************************************* //
ETAB.C (5,699 bytes)   

dkxls

2013-09-04 20:45

reporter  

TAB.C (6,350 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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 "TAB.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CloudType>
Foam::TAB<CloudType>::TAB
(
    const dictionary& dict,
    CloudType& owner
)
:
    BreakupModel<CloudType>(dict, owner, typeName, true),
    Cmu_(BreakupModel<CloudType>::TABCmu_),
    Comega_(BreakupModel<CloudType>::TABComega_),
    WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
    SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
{
    // calculate the inverse function of the Rossin-Rammler Distribution
    const scalar xx0 = 12.0;
    const scalar rrd100 =
        1.0/(1.0 - exp(-xx0)*(1.0 + xx0 + sqr(xx0)/2.0 + pow3(xx0)/6.0));

    forAll(rrd_, n)
    {
        scalar xx = 0.12*(n + 1);
        rrd_[n] =
            (1.0 - exp(-xx)*(1.0 + xx + sqr(xx)/2.0 + pow3(xx)/6.0))*rrd100;
    }

    if (SMDCalcMethod_ == "method1")
    {
        SMDMethod_ = method1;
    }
    else if (SMDCalcMethod_ == "method2")
    {
        SMDMethod_ = method2;
    }
    else
    {
        SMDMethod_ = method2;
        WarningIn("Foam::TAB<CloudType>::TAB(const dictionary&, CloudType&)")
            << "Unknown SMDCalculationMethod. Valid options are "
            << "(method1 | method2). Using method2" << endl;
    }
}


template<class CloudType>
Foam::TAB<CloudType>::TAB(const TAB<CloudType>& bum)
:
    BreakupModel<CloudType>(bum),
    Cmu_(bum.Cmu_),
    Comega_(bum.Comega_),
    WeCrit_(bum.WeCrit_),
    SMDCalcMethod_(bum.SMDCalcMethod_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CloudType>
Foam::TAB<CloudType>::~TAB()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class CloudType>
bool Foam::TAB<CloudType>::update
(
    const scalar dt,
    const vector& g,
    scalar& d,
    scalar& tc,
    scalar& ms,
    scalar& nParticle,
    scalar& KHindex,
    scalar& y,
    scalar& yDot,
    const scalar d0,
    const scalar rho,
    const scalar mu,
    const scalar sigma,
    const vector& U,
    const scalar rhoc,
    const scalar muc,
    const vector& Urel,
    const scalar Urmag,
    const scalar tMom,
    scalar& dChild,
    scalar& massChild
)
{
    cachedRandom& rndGen = this->owner().rndGen();

    scalar r = 0.5*d;
    scalar r2 = r*r;
    scalar r3 = r*r2;

    scalar semiMass = nParticle*pow3(d);

    // inverse of characteristic viscous damping time
    scalar rtd = 0.5*Cmu_*mu/(rho*r2);

    // oscillation frequency (squared)
    scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;

    if (omega2 > 0)
    {
        scalar omega = sqrt(omega2);
        scalar We = rhoc*sqr(Urmag)*r/sigma;
        scalar Wetmp = We/WeCrit_;

        scalar y1 = y - Wetmp;
        scalar y2 = yDot/omega;

        scalar a = sqrt(y1*y1 + y2*y2);

        // scotty we may have break-up
        if (a+Wetmp > 1.0)
        {
            scalar phic = y1/a;

            // constrain phic within -1 to 1
            phic = max(min(phic, 1), -1);

            scalar phit = acos(phic);
            scalar phi = phit;
            scalar quad = -y2/a;
            if (quad < 0)
            {
                phi = constant::mathematical::twoPi - phit;
            }

            scalar tb = 0;

            if (mag(y) < 1.0)
            {
                scalar coste = 1.0;
                if ((Wetmp - a < -1) && (yDot < 0))
                {
                    coste = -1.0;
                }

                scalar theta = acos((coste-Wetmp)/a);

                if (theta < phi)
                {
                    if (constant::mathematical::twoPi - theta >= phi)
                    {
                        theta = -theta;
                    }
                    theta += constant::mathematical::twoPi;
                }
                tb = (theta-phi)/omega;

                // breakup occurs
                if (dt > tb)
                {
                    y = 1.0;
                    yDot = -a*omega*sin(omega*tb + phi);
                }

            }

            // update droplet size
            if (dt > tb)
            {
                scalar rs =
                    r/(1.0 + (4.0/3.0)*sqr(y) + rho*r3/(8*sigma)*sqr(yDot));

                label n = 0;
                scalar rNew = 0.0;
                switch (SMDMethod_)
                {
                    case method1:
                    {
                        #include "TABSMDCalcMethod1.H"
                        break;
                    }
                    case method2:
                    {
                        #include "TABSMDCalcMethod2.H"
                        break;
                    }
                }

                if (rNew < r)
                {
                    d = 2*rNew;
                    y = 0;
                    yDot = 0;
                }
            }
        }
    }
    else
    {
        // reset droplet distortion parameters
        y = 0;
        yDot = 0;
    }

    // update the nParticle count to conserve mass
    nParticle = semiMass/pow3(d);

    // Do not add child parcel
    return false;
}


// ************************************************************************* //
TAB.C (6,350 bytes)   

user2

2013-09-05 18:19

  ~0002467

Thanks for the report - fixed by commit 1064355

Issue History

Date Modified Username Field Change
2013-09-04 19:16 dkxls New Issue
2013-09-04 19:16 dkxls Tag Attached: breakup
2013-09-04 19:16 dkxls Tag Attached: spray
2013-09-04 19:16 dkxls Tag Attached: ETAB
2013-09-04 19:16 dkxls Tag Attached: TAB
2013-09-04 19:37 dkxls File Added: 0001-Fix-inconsistent-breakup-coefficient-reading.patch
2013-09-04 19:37 dkxls File Added: 0001-Fix-inconsistent-ETAB-coefficient-reading.patch
2013-09-04 19:46 dkxls Note Added: 0002461
2013-09-04 20:38 dkxls File Added: 0001-Enhance-TAB-coefficient-reading-and-add-enableOscill.patch
2013-09-04 20:44 dkxls Note Added: 0002462
2013-09-04 20:45 dkxls File Added: BreakupModel.H
2013-09-04 20:45 dkxls File Added: BreakupModel.C
2013-09-04 20:45 dkxls File Added: ETAB.C
2013-09-04 20:45 dkxls File Added: TAB.C
2013-09-05 00:24 dkxls Note Edited: 0002461
2013-09-05 18:19 user2 Note Added: 0002467
2013-09-05 18:19 user2 Status new => resolved
2013-09-05 18:19 user2 Fixed in Version => 2.2.x
2013-09-05 18:19 user2 Resolution open => fixed
2013-09-05 18:19 user2 Assigned To => user2
2018-07-10 11:18 administrator Tag Detached: ETAB
2018-07-10 11:21 administrator Tag Detached: TAB