Update to disengagement functionObject

* Correct computation of holdup (trough phase center of mass)
* Added switches for direction and phase names

Author: fmuni

OFsolvers/applications/birdmultiphaseEulerFoam/functionObjects/disengagement/disengagement.C

@@ -69,15 +69,17 @@ Foam::functionObjects::disengagement::disengagement
     logFiles(obr_,name),
     phases_(mesh_.lookupObject<phaseSystem>("phaseProperties").phases()),
     phaseName_(dict.lookup<word>("phase")),
+    inletPhaseName_(dict.lookup<word>("inletPhase")),
     inletPatch_(dict.lookup<word>("inlet")),
     tolerance_(dict.lookup<scalar>("tolerance")),
     nsamples_(dict.lookup<label>("nsamples")),
-    disengage_(dict.lookup("disengage")),
-    holdup_(2*nsamples_,{0.,-1.}),
+    direction_(dict.lookup<vector>("direction")),
+    disengage_(dict.lookup<bool>("disengage")),
+    phase_com_(2*nsamples_,{0.,-1.}),
     disengaged_(false)
 {
     //- Check that U has fixedValue
-    volVectorField& U = const_cast<volVectorField&>(mesh_.lookupObject<volVectorField>("U." + phaseName_));
+    volVectorField& U = const_cast<volVectorField&>(mesh_.lookupObject<volVectorField>("U." + inletPhaseName_));
     label patchI = mesh_.boundaryMesh().findPatchID(inletPatch_);
     volVectorField::Boundary& UBf = const_cast<volVectorField::Boundary&>(U.boundaryFieldRef());
 
@@ -106,9 +108,9 @@ void Foam::functionObjects::disengagement::writeFileHeader(const label i)
 {
     if (Pstream::master())
     {
-        writeHeader(file(), "holdup");
+        writeHeader(file(), "phase_com");
         writeCommented(file(), "time");
-        writeTabbed(file(), "holdup");
+        writeTabbed(file(), "phase_com");
         writeTabbed(file(), "disengaged");
 
         file() << endl;
@@ -119,13 +121,14 @@ bool Foam::functionObjects::disengagement::execute()
 {
     const volScalarField& alpha = mesh_.lookupObject<volScalarField>("alpha." + phaseName_);
 
-    //- Compute holdup
-    scalar volume = gSum(mesh_.V());
-    scalar holdup = gSum(fvc::volumeIntegrate(alpha))/volume;
-    Pair<scalar> holddata(mesh_.time().value(), holdup);
+    //- Compute phase_com
+    scalar volume = gSum(fvc::volumeIntegrate(alpha));
+    volScalarField hcoord = mesh_.C()&direction_;
+    scalar phase_com = gSum(fvc::volumeIntegrate(hcoord*alpha))/volume;
+    Pair<scalar> holddata(mesh_.time().value(), phase_com);
 
     //- See function at the beginning of file
-    add_to_list_like_static_queue(holdup_,holddata);
+    add_to_list_like_static_queue(phase_com_,holddata);
 
     //- Skip rest if no disengagement is required
     if (!disengage_) return true;
@@ -134,35 +137,35 @@ bool Foam::functionObjects::disengagement::execute()
     if (disengaged_ ) return true;
 
     //- Do the check only if the list is complete
-    if (holdup_[0].second() > 0.)
+    if (phase_com_[0].second() > 0.)
     {
-        scalar holdup_mean_long(0.);
-        scalar holdup_mean_short(0.);
-        scalar t0_long(holdup_[0].first());
+        scalar phase_com_mean_long(0.);
+        scalar phase_com_mean_short(0.);
+        scalar t0_long(phase_com_[0].first());
         scalar deltat_long(0.);
         scalar deltat_short(0.);
 
         //- The first (oldest) sample is skipped to have a well-defined dt
         for (int i = 1; i < 2*nsamples_; i++)
         {
-            scalar dt = holdup_[i].first() - t0_long;
-            t0_long = holdup_[i].first();
-            holdup_mean_long += dt*holdup_[i].second();
+            scalar dt = phase_com_[i].first() - t0_long;
+            t0_long = phase_com_[i].first();
+            phase_com_mean_long += dt*phase_com_[i].second();
             deltat_long += dt;
 
             //- Compute average on the most recent samples
             if (i >= nsamples_)
             {
-                holdup_mean_short +=  dt*holdup_[i].second();
+                phase_com_mean_short +=  dt*phase_com_[i].second();
                 deltat_short += dt;
             }
 
         }
 
-        holdup_mean_long /= deltat_long;
-        holdup_mean_short /= deltat_short;
+        phase_com_mean_long /= deltat_long;
+        phase_com_mean_short /= deltat_short;
 
-        if (mag(holdup_mean_long - holdup_mean_short) < tolerance_)
+        if (mag(phase_com_mean_long - phase_com_mean_short) < tolerance_)
         {
             if(!disengaged_)
             {
@@ -172,7 +175,7 @@ bool Foam::functionObjects::disengagement::execute()
             disengaged_ = true;
 
             //- Get boundary condition
-            volVectorField& U = const_cast<volVectorField&>(mesh_.lookupObject<volVectorField>("U." + phaseName_));
+            volVectorField& U = const_cast<volVectorField&>(mesh_.lookupObject<volVectorField>("U." + inletPhaseName_));
 
             label patchI = mesh_.boundaryMesh().findPatchID(inletPatch_);
 
@@ -199,9 +202,9 @@ bool Foam::functionObjects::disengagement::write()
 
     if (Pstream::master())
     {
-            file() << holdup_[2*nsamples_ -1].first();
+            file() << phase_com_[2*nsamples_ -1].first();
             file() << tab;
-            file() << holdup_[2*nsamples_ -1].second();
+            file() << phase_com_[2*nsamples_ -1].second();
             file() << tab;
             file() << dis;
             file() << endl;

OFsolvers/applications/birdmultiphaseEulerFoam/functionObjects/disengagement/disengagement.H

@@ -76,6 +76,9 @@ protected:
         //- Name of phase
         word phaseName_;
 
+        //- Name of phase
+        word inletPhaseName_;
+
         //- Name of inlet patch
         word inletPatch_;
 
@@ -85,11 +88,14 @@ protected:
         //- Number of samples for averaging
         label nsamples_;
 
+        //- Direction for computing the holdup
+        vector direction_;
+
         //- Optional for not disengaging
         bool disengage_;
 
-        //- Vector with holdups and times
-        List<Pair<scalar>> holdup_;
+        //- Vector with phase ceneter of mass and times
+        List<Pair<scalar>> phase_com_;
 
         //- Check if disengaged
         bool disengaged_;