Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,17 +22,13 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    porousExplicitSourceReactingParcelFoam
+    LTSReactingParcelFoam
 
 Description
-    Transient PISO solver for compressible, laminar or turbulent flow with
-    reacting multiphase Lagrangian parcels for porous media, including explicit
-    sources for mass, momentum and energy
-
-    The solver includes:
-    - reacting multiphase parcel cloud
-    - porous media
-    - mass, momentum and energy sources
+    Local time stepping (LTS) solver for steady, compressible, laminar or
+    turbulent reacting and non-reacting flow with multiphase Lagrangian
+    parcels and porous media, including explicit sources for mass, momentum
+    and energy
 
     Note: ddtPhiCorr not used here when porous zones are active
     - not well defined for porous calculations
@@ -75,7 +71,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readSIMPLEControls.H"
+        #include "readPISOControls.H"
         #include "readChemistryProperties.H"
         #include "readAdditionalSolutionControls.H"
         #include "readTimeControls.H"
@@ -84,26 +80,24 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
+        parcels.evolve();
 
-        // --- Pressure-velocity corrector
+        #include "chemistry.H"
+        #include "timeScales.H"
+
+        #include "rhoEqn.H"
+        #include "UEqn.H"
+        #include "YEqn.H"
+        #include "hsEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
         {
-            parcels.evolve();
-
-            #include "chemistry.H"
-
-            #include "timeScales.H"
-
-            #include "rhoEqn.H"
-            #include "UEqn.H"
-            #include "YEqn.H"
-            #include "hsEqn.H"
-
             #include "pEqn.H"
-
-            turbulence->correct();
         }
 
+        turbulence->correct();
+
         if (runTime.write())
         {
             chemistry.dQ()().write();

applications/solvers/lagrangian/LTSReactingParcelFoam/Make/files

@@ -0,0 +1,3 @@
+LTSReactingParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/LTSReactingParcelFoam

applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H

@@ -73,12 +73,12 @@
 
     dimensionedScalar rhoMax
     (
-        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
+        mesh.solutionDict().subDict("PISO").lookup("rhoMax")
     );
 
     dimensionedScalar rhoMin
     (
-        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
+        mesh.solutionDict().subDict("PISO").lookup("rhoMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;

applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H

@@ -48,9 +48,6 @@
         }
     }
 
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
-
     Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 
     // Second part of thermodynamic density update

applications/solvers/lagrangian/LTSReactingParcelFoam/readTimeControls.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files

@@ -1,3 +0,0 @@
-steadyReactingParcelFoam.C
-
-EXE = $(FOAM_APPBIN)/steadyReactingParcelFoam

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C

@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Analytical.H"
-#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -59,9 +58,9 @@ template<class Type>
 typename Foam::IntegrationScheme<Type>::integrationResult
 Foam::Analytical<Type>::integrate
 (
-    const Type phi,
+    const Type& phi,
     const scalar dt,
-    const Type alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
@@ -69,9 +68,18 @@ Foam::Analytical<Type>::integrate
 
     const scalar expTerm = exp(min(50, -beta*dt));
 
-    retValue.average() =
-        alpha + (phi - alpha)*(1 - expTerm)/(beta*dt + ROOTVSMALL);
-    retValue.value() =  alpha + (phi - alpha)*expTerm;
+    if (beta > ROOTVSMALL)
+    {
+        const Type alpha = alphaBeta/beta;
+        retValue.average() = alpha + (phi - alpha)*(1 - expTerm)/(beta*dt);
+        retValue.value() =  alpha + (phi - alpha)*expTerm;
+    }
+    else
+    {
+        retValue.average() = phi;
+        retValue.value() =  phi;
+    }
+
 
     return retValue;
 }

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H

@@ -81,9 +81,9 @@ public:
         //- Perform the integration
         virtual typename IntegrationScheme<Type>::integrationResult integrate
         (
-            const Type phi,
+            const Type& phi,
             const scalar dt,
-            const Type alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C

@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Euler.H"
-#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -59,14 +58,14 @@ template<class Type>
 typename Foam::IntegrationScheme<Type>::integrationResult
 Foam::Euler<Type>::integrate
 (
-    const Type phi,
+    const Type& phi,
     const scalar dt,
-    const Type alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
     typename IntegrationScheme<Type>::integrationResult retValue;
-    retValue.value() = (phi + beta*dt*alpha)/(1.0 + beta*dt);
+    retValue.value() = (phi + alphaBeta*dt)/(1.0 + beta*dt);
     retValue.average() = 0.5*(phi + retValue.value());
 
     return retValue;

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H

@@ -78,9 +78,9 @@ public:
         //- Perform the integration
         virtual typename IntegrationScheme<Type>::integrationResult integrate
         (
-            const Type phi,
+            const Type& phi,
             const scalar dt,
-            const Type alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C

@@ -60,9 +60,9 @@ template<class Type>
 typename Foam::IntegrationScheme<Type>::integrationResult
 Foam::IntegrationScheme<Type>::integrate
 (
-    const Type phi,
+    const Type& phi,
     const scalar dt,
-    const Type alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
@@ -71,9 +71,9 @@ Foam::IntegrationScheme<Type>::integrate
         "Foam::IntegrationScheme<Type>::integrationResult"
         "Foam::IntegrationScheme<Type>::integrate"
         "("
-            "const Type, "
+            "const Type&, "
             "const scalar, "
-            "const Type, "
+            "const Type&, "
             "const scalar"
         ") const"
     );

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H

@@ -183,9 +183,9 @@ public:
         //- Perform the Integration
         virtual integrationResult integrate
         (
-            const Type phi,
+            const Type& phi,
             const scalar dt,
-            const Type alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };

src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C

@@ -114,11 +114,11 @@ Foam::CollidingCloud<CloudType>::CollidingCloud
     if (this->solution().active())
     {
         setModels();
-    }
 
-    if (readFields)
-    {
-        parcelType::readFields(*this);
+        if (readFields)
+        {
+            parcelType::readFields(*this);
+        }
     }
 }
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -360,11 +360,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     if (solution_.active())
     {
         setModels();
-    }
 
-    if (readFields)
-    {
-        parcelType::readFields(*this);
+        if (readFields)
+        {
+            parcelType::readFields(*this);
+        }
     }
 
     if (solution_.resetSourcesOnStartup())
@@ -504,7 +504,7 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
         parcel.rho() = constProps_.rho0();
     }
 
-    const scalar carrierDt = this->db().time().deltaTValue();
+    const scalar carrierDt = solution_.deltaT();
     parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
     parcel.typeId() = constProps_.parcelTypeId();
 }

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -117,6 +117,11 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     if (this->solution().active())
     {
         setModels();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this, this->composition());
+        }
     }
 
     // Set storage for mass source fields and initialise to zero
@@ -142,11 +147,6 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
         );
     }
 
-    if (readFields)
-    {
-        parcelType::readFields(*this, this->composition());
-    }
-
     if (this->solution().resetSourcesOnStartup())
     {
         resetSourceTerms();

src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C

@@ -94,11 +94,11 @@ Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
     if (this->solution().active())
     {
         setModels();
-    }
 
-    if (readFields)
-    {
-        parcelType::readFields(*this, this->composition());
+        if (readFields)
+        {
+            parcelType::readFields(*this, this->composition());
+        }
     }
 
     if (this->solution().resetSourcesOnStartup())

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -135,11 +135,11 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     if (this->solution().active())
     {
         setModels();
-    }
 
-    if (readFields)
-    {
-        parcelType::readFields(*this);
+        if (readFields)
+        {
+            parcelType::readFields(*this);
+        }
     }
 
     if (this->solution().resetSourcesOnStartup())

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -192,21 +192,20 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const parcelType& p = static_cast<const parcelType&>(*this);
     const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
-    forceSuSp Feff = Fcp + Fncp;
-    Feff.Sp() += ROOTVSMALL;
+    const forceSuSp Feff = Fcp + Fncp;
 
 
     // New particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
-    const vector ap = Uc_ + (Feff.Su() + Su)/Feff.Sp();
+    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/mass;
     const scalar bp = Feff.Sp()/mass;
 
-    Spu = Feff.Sp();
+    Spu = Feff.Sp()*dt/td.cloud().solution().deltaT();
 
     IntegrationScheme<vector>::integrationResult Ures =
-        td.cloud().UIntegrator().integrate(U, dt, ap, bp);
+        td.cloud().UIntegrator().integrate(U, dt, abp, bp);
 
     vector Unew = Ures.value();
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -333,13 +333,13 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     // Integrate to find the new parcel temperature
     IntegrationScheme<scalar>::integrationResult Tres =
-        td.cloud().TIntegrator().integrate(T, dt, ap, bp);
+        td.cloud().TIntegrator().integrate(T, dt, ap*bp, bp);
 
     scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
 
     dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
 
-    Cuh = bp;
+    Cuh = bp*dt/td.cloud().solution().deltaT();
 
     return Tnew;
 }

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -520,7 +520,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
     }
 
     const scalar time = this->owner().db().time().value();
-    const scalar carrierDt = this->owner().db().time().deltaTValue();
+    const scalar carrierDt = this->owner().solution().deltaT();
     const polyMesh& mesh = this->owner().mesh();
 
     // Prepare for next time step

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C

@@ -290,7 +290,7 @@ void Foam::PatchInteractionModel<CloudType>::patchData
         }
         else
         {
-            Up = (Cf - Cf00)/mesh.time().deltaTValue();
+            Up = (Cf - Cf00)/this->owner().solution().deltaT();
         }
 
         if (mag(dn) > SMALL)
@@ -312,7 +312,9 @@ void Foam::PatchInteractionModel<CloudType>::patchData
             // magOmega = sin(angle between unit normals)
             // Normalise omega vector by magOmega, then multiply by
             // angle/dt to give the correct angular velocity vector.
-            omega *= Foam::asin(magOmega)/(magOmega*mesh.time().deltaTValue());
+            omega *=
+                Foam::asin(magOmega)
+               /(magOmega*this->owner().solution().deltaT());
 
             // Project position onto face and calculate this position
             // relative to the face centre.

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/Allrun

@@ -8,4 +8,4 @@ cd ${0%/*} || exit 1    # run from this directory
 runApplication blockMesh
 
 # run the solver
-runApplication steadyReactingParcelFoam
+runApplication LTSReactingParcelFoam

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/controlDict

@@ -51,4 +51,8 @@ maxCo           0.5;
 
 maxDeltaT       1;
 
+alphaTauSpecie  1e10;
+
+alphaTauTemp    0.005;
+
 // ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSolution

@@ -53,8 +53,9 @@ solvers
     }
 }
 
-SIMPLE
+PISO
 {
+    nCorrectors         2;
     nNonOrthogonalCorrectors 0;
     momentumPredictor yes;
     rhoMin rhoMin [1 -3 0 0 0] 0.1;
@@ -68,10 +69,5 @@ additional
     solveSpecies    true;
 }
 
-relaxationFactors
-{
-    p               0.3;
-}
-
 
 // ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/Allrun

@@ -8,6 +8,6 @@ cd ${0%/*} || exit 1    # run from this directory
 runApplication blockMesh
 
 # run the solver
-runApplication steadyReactingParcelFoam
+runApplication LTSReactingParcelFoam
 
 # ----------------------------------------------------------------- end-of-file