MULES limitation removed: sub-cycling time now supported on morphing meshes

Support for cell-volume interpolation during time sub-cycling now provided and
used in the MULES VoF solver to allow sub-cycling on morphing meshes.

src/OpenFOAM/db/Time/Time.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -625,21 +625,25 @@ void Foam::Time::setDeltaT(const scalar deltaT)
 Foam::TimeState Foam::Time::subCycle(const label nSubCycles)
 {
     subCycling_ = true;
+    prevTimeState_.set(new TimeState(*this));
 
-    TimeState ts = *this;
     setTime(*this - deltaT(), (timeIndex() - 1)*nSubCycles);
     deltaT_ /= nSubCycles;
     deltaT0_ /= nSubCycles;
     deltaTSave_ = deltaT0_;
 
-    return ts;
+    return prevTimeState();
 }
 
 
-void Foam::Time::endSubCycle(const TimeState& ts)
+void Foam::Time::endSubCycle()
 {
+    if (subCycling_)
+    {
         subCycling_ = false;
-    TimeState::operator=(ts);
+        TimeState::operator=(prevTimeState());
+        prevTimeState_.clear();
+    }
 }
 
 

src/OpenFOAM/db/Time/Time.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,6 +127,9 @@ protected:
         //- Is the time currently being sub-cycled?
         bool subCycling_;
 
+        //- If time is being sub-cycled this is the previous TimeState
+        autoPtr<TimeState> prevTimeState_;
+
         //- Time directory name format
         static fmtflags format_;
 
@@ -347,6 +350,18 @@ public:
                 return functionObjects_;
             }
 
+            //- Return true if time currently being sub-cycled, otherwise false
+            bool subCycling() const
+            {
+                return subCycling_;
+            }
+
+            //- Return previous TimeState if time is being sub-cycled
+            const TimeState& prevTimeState() const
+            {
+                return prevTimeState_();
+            }
+
 
         // Check
 
@@ -427,8 +442,8 @@ public:
             //- Set time to sub-cycle for the given number of steps
             virtual TimeState subCycle(const label nSubCycles);
 
-            //- Reset time after sub-cycling back to given TimeState
+            //- Reset time after sub-cycling back to previous TimeState
-            virtual void endSubCycle(const TimeState&);
+            virtual void endSubCycle();
 
             //- Return non-const access to the list of function objects
             functionObjectList& functionObjects()

src/OpenFOAM/db/Time/subCycleTime.C

@@ -32,9 +32,10 @@ Foam::subCycleTime::subCycleTime(Time& t, const label nSubCycles)
 :
     time_(t),
     nSubCycles_(nSubCycles),
-    subCycleIndex_(0),
+    subCycleIndex_(0)
-    initialTimeState_(time_.subCycle(nSubCycles_))
+{
-{}
+    time_.subCycle(nSubCycles_);
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -55,7 +56,7 @@ bool Foam::subCycleTime::end() const
 
 void Foam::subCycleTime::endSubCycle()
 {
-    time_.endSubCycle(initialTimeState_);
+    time_.endSubCycle();
 }
 
 

src/OpenFOAM/db/Time/subCycleTime.H

@@ -55,7 +55,6 @@ class subCycleTime
 
         label nSubCycles_;
         label subCycleIndex_;
-        TimeState initialTimeState_;
 
 public:
 

src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ License
 #include "fvcSurfaceIntegrate.H"
 #include "slicedSurfaceFields.H"
 #include "syncTools.H"
+
 #include "fvm.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -107,7 +108,7 @@ void Foam::MULES::explicitSolve
     {
         psiIf =
         (
-            mesh.V0()*rho.oldTime()*psi0/(deltaT*mesh.V())
+            mesh.Vsc0()*rho.oldTime()*psi0/(deltaT*mesh.Vsc())
           + Su.field()
           - psiIf
         )/(rho/deltaT - Sp.field());
@@ -325,7 +326,8 @@ void Foam::MULES::limiter
 
     const unallocLabelList& owner = mesh.owner();
     const unallocLabelList& neighb = mesh.neighbour();
-    const scalarField& V = mesh.V();
+    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
+    const scalarField& V = tVsc();
     const scalar deltaT = mesh.time().deltaTValue();
 
     const scalarField& phiBDIf = phiBD;
@@ -452,14 +454,16 @@ void Foam::MULES::limiter
 
     if (mesh.moving())
     {
+        tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
+
         psiMaxn =
             V*((rho/deltaT - Sp)*psiMaxn - Su)
-          - (mesh.V0()/deltaT)*rho.oldTime()*psi0
+          - (V0()/deltaT)*rho.oldTime()*psi0
           + sumPhiBD;
 
         psiMinn =
             V*(Su - (rho/deltaT - Sp)*psiMinn)
-          + (mesh.V0()/deltaT)*rho.oldTime()*psi0
+          + (V0/deltaT)*rho.oldTime()*psi0
           - sumPhiBD;
     }
     else

src/finiteVolume/fvMesh/fvMesh.H

@@ -269,6 +269,12 @@ public:
             //- Return old-old-time cell volumes
             const DimensionedField<scalar, volMesh>& V00() const;
 
+            //- Return sub-cycle cell volumes
+            tmp<DimensionedField<scalar, volMesh> > Vsc() const;
+
+            //- Return sub-cycl old-time cell volumes
+            tmp<DimensionedField<scalar, volMesh> > Vsc0() const;
+
             //- Return cell face area vectors
             const surfaceVectorField& Sf() const;
 

src/finiteVolume/fvMesh/fvMeshGeometry.C

@@ -285,6 +285,63 @@ const volScalarField::DimensionedInternalField& fvMesh::V00() const
 }
 
 
+tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc() const
+{
+    if (moving() && time().subCycling())
+    {
+        const TimeState& ts = time();
+        const TimeState& ts0 = time().prevTimeState();
+
+        scalar tFrac =
+        (
+            ts.value() - (ts0.value() - ts0.deltaTValue())
+        )/ts0.deltaTValue();
+
+        if (tFrac < (1 - SMALL))
+        {
+            return V0() + tFrac*(V() - V0());
+        }
+        else
+        {
+            return V();
+        }
+    }
+    else
+    {
+        return V();
+    }
+}
+
+
+tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc0() const
+{
+    if (moving() && time().subCycling())
+    {
+        const TimeState& ts = time();
+        const TimeState& ts0 = time().prevTimeState();
+
+        scalar t0Frac =
+        (
+            (ts.value() - ts.deltaTValue())
+          - (ts0.value() - ts0.deltaTValue())
+        )/ts0.deltaTValue();
+
+        if (t0Frac > SMALL)
+        {
+            return V0() + t0Frac*(V() - V0());
+        }
+        else
+        {
+            return V0();
+        }
+    }
+    else
+    {
+        return V0();
+    }
+}
+
+
 const surfaceVectorField& fvMesh::Sf() const
 {
     if (!SfPtr_)

tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement

@@ -37,7 +37,8 @@ boundaryField
         centreOfMass    (0.5 0.5 0.5);
         momentOfInertia (0.08622222 0.08622222 0.144);
         mass            9.6;
-        rhoInf          1;                      // for forces calculation
+        rhoInf          1;
+        report          on;
         value           uniform (0 0 0);
     }
 }

tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/polyMesh/boundary

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6.x                                 |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict

@@ -17,7 +17,7 @@ FoamFile
 
 application     interDyMFoam;
 
-startFrom       startTime;
+startFrom       latestTime;
 
 startTime       0;
 
@@ -29,7 +29,7 @@ deltaT          0.01;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.025;
+writeInterval   0.1;
 
 purgeWrite      0;
 
@@ -47,9 +47,9 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.2;
+maxCo           0.5;
 
-maxDeltaT       0.025;
+maxDeltaT       0.01;
 
 libs
 (

tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution

@@ -20,7 +20,7 @@ solvers
     cellDisplacement
     {
         solver          GAMG;
-        tolerance       1e-08;
+        tolerance       1e-5;
         relTol          0;
         smoother        GaussSeidel;
         cacheAgglomeration true;
@@ -35,7 +35,7 @@ solvers
         preconditioner
         {
             preconditioner  GAMG;
-            tolerance       1e-05;
+            tolerance       1e-5;
             relTol          0;
             smoother        DICGaussSeidel;
             nPreSweeps      0;
@@ -55,7 +55,7 @@ solvers
     p
     {
         solver          GAMG;
-        tolerance       1e-08;
+        tolerance       1e-8;
         relTol          0.01;
         smoother        DIC;
         nPreSweeps      0;
@@ -73,7 +73,7 @@ solvers
         preconditioner
         {
             preconditioner  GAMG;
-            tolerance       2e-09;
+            tolerance       1e-8;
             relTol          0;
             nVcycles        2;
             smoother        DICGaussSeidel;
@@ -86,7 +86,7 @@ solvers
             mergeLevels     1;
         }
 
-        tolerance       2e-09;
+        tolerance       1e-8;
         relTol          0;
         maxIter         20;
     }
@@ -95,18 +95,18 @@ solvers
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
-        tolerance       1e-06;
+        tolerance       1e-6;
         relTol          0;
         nSweeps         1;
     }
 
-    "(k|epsilon|R|nuTilda)"
+    "(k|epsilon|omega|R|nuTilda)"
     {
-        $U;
+        solver          PBiCG;
-        tolerance       1e-08;
+        preconditioner  DILU;
+        tolerance       1e-8;
         relTol          0;
     }
-
 }
 
 PISO
@@ -115,14 +115,14 @@ PISO
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
     nAlphaCorr      1;
-    nAlphaSubCycles 1;
+    nAlphaSubCycles 3;
-    cAlpha          1.5;
+    cAlpha          1;
     correctPhi      yes;
 }
 
 relaxationFactors
 {
-    U               1;
+    "(U|k|epsilon|omega|R|nuTilda)"               1;
 }