View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0000999OpenFOAMBugpublic2013-09-08 14:492013-09-18 14:13
Reporterdkxls Assigned Touser2 
PriorityurgentSeveritymajorReproducibilityN/A
Status resolvedResolutionfixed 
PlatformLinux x86_64OSopenSUSEOS Version12.2
Summary0000999: [Lagrangian]: Incorrect treatment of coupled and uncoupled solutions
DescriptionFrom cloudSolution.H the definition of the 'coupled' switch is:
"Flag to indicate whether parcels are coupled to the carrier phase, i.e. whether or not to generate source terms for carrier phase".
This definition is inconsistent, as the first part says that the phases are not coupled (i.e. 0-way coupling) and the second part suggests a 1-way coupling.
A detailed description and justification of the different interaction stages is given in the paper by Elghobashi (1994).

The actual implementation leads for uncoupled simulation to a 1-way coupling, with several flaws/inconsistencies, especially with respect to spray simulation.

1. The coupling in the present implementation is applied in the *Parcel templates to all source terms (momentum, heat and mass), which leads in the case of a uncoupled calculation to a loss in mass and energy.
Where the loss in mass and energy might be justified in certain cases, it is NOT for spray simulation!
Hence the decision, if source terms are generated or not should not be made in the *Parcel templates, but only in the *Cloud templates!

2. In spray simulation it is often assumed that an intact liquid core exists. This liquid core is typically modeled by injecting large "blobs" in the size of the nozzle diameter and the extend of the liquid core region is typically determined by a breakup length. For the "blobs" in the liquid core zero-drag is assumed, but heat and mass transfer is still calculated and also contributed to the source term.
In the present SprayParcel implementation, the calculation of the droplet properties inside the liquid core is switched to uncoupled leading to wrong droplet velocities and momentum, mass and energy losses. This error results from the incorrect coupling implementation, which applies a 1-way coupling instead of a no-drag assumption (0-way coupling for the momentum transfer).

As the infrastructure in the source code is largely given with the division in coupled and uncoupled forces, the changes are rather small.
It has to be distinguished between 0, 1, and 2-way coupling, not only between coupled and uncoupled:
0-way = no momentum transfer
1-way = momentum transfer: continues --> dispersed
2-way = momentum transfer: continues <-> dispersed

This division can easily be implemented with a scalar variable indicating the coupling regime, which leads to the following calculations:
0-way coupling: Only non-coupled forces are calculated
1-way coupling: Coupled and non-coupled forces are calculated, but no contribution to the source term
2-way coupling: Coupled and non-coupled forces are calculated and changes are contributed to the source term

I uploaded a patch that fixes this issue without breaking the compatibility for non-spray solvers.
To retain the old functionality for non-spray solver a new scalar 'couplingRegime' is introduced, where couplingRegime=1 corresponds to coupled=false and couplingRegime=2 corresponds to coupled=true. The no drag assumption equals to couplingRegime=0 (i.e. 0-way coupling).
In a future release the coupled switch can also be removed entirely as it has no extra function.
Since the patch includes a lot white-space changes in the *Parcel templates, it is a bit messy. Hence I also uploaded a patch that ignores these with otherwise the same changes (git diff -b).

References:
Elghobashi, S., "On predicting particle-laden turbulent flows", Applied Scientific Research, Volume 52, Issue 4, pp 309-329, June 1994, http://dx.doi.org/10.1007/BF00936835
Additional InformationThe dieselSpray implementation in 2.0.x handles this issue correctly (i.e. no drag assumption for parcels in the liquid core).
https://github.com/OpenFOAM/OpenFOAM-2.0.x/blob/master/src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C#L137
TagsNo tags attached.

Activities

dkxls

2013-09-08 14:49

reporter  

Correct-coupling-of-parcels-and-carrier-phase_no-whitespace-changes.patch (14,918 bytes)    
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index f3a8f16..936b50c 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -88,7 +88,7 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
-        if (solution_.coupled())
+        if (solution_.couplingRegime() > 1.5)
         {
             td.cloud().relaxSources(td.cloud().cloudCopy());
         }
@@ -99,7 +99,7 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
-        if (solution_.coupled())
+        if (solution_.couplingRegime() > 1.5)
         {
             td.cloud().scaleSources();
         }
@@ -165,7 +165,7 @@ template<class CloudType>
 template<class TrackData>
 void Foam::KinematicCloud<CloudType>::evolveCloud(TrackData& td)
 {
-    if (solution_.coupled())
+    if (solution_.couplingRegime() > 1.5)
     {
         td.cloud().resetSourceTerms();
     }
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 72cad26..2dbe962 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -555,7 +555,7 @@ Foam::KinematicCloud<CloudType>::SU(volVectorField& U) const
             << max(UCoeff()).value() << endl;
     }
 
-    if (solution_.coupled())
+    if (solution_.couplingRegime() > 1.5)
     {
         if (solution_.semiImplicit("U"))
         {
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
index fd7bf31..7fab1df 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
@@ -43,6 +43,7 @@ Foam::cloudSolution::cloudSolution
     iter_(1),
     trackTime_(0.0),
     coupled_(false),
+    couplingRegime_(1.0),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
     resetSourcesOnStartup_(true),
@@ -69,6 +70,7 @@ Foam::cloudSolution::cloudSolution
     iter_(cs.iter_),
     trackTime_(cs.trackTime_),
     coupled_(cs.coupled_),
+    couplingRegime_(cs.couplingRegime_),
     cellValueSourceCorrection_(cs.cellValueSourceCorrection_),
     maxTrackTime_(cs.maxTrackTime_),
     resetSourcesOnStartup_(cs.resetSourcesOnStartup_),
@@ -90,6 +92,7 @@ Foam::cloudSolution::cloudSolution
     iter_(0),
     trackTime_(0.0),
     coupled_(false),
+    couplingRegime_(0.0),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
     resetSourcesOnStartup_(false),
@@ -126,6 +129,8 @@ void Foam::cloudSolution::read()
 
     if (coupled_)
     {
+        couplingRegime_ = 2.0;
+
         const dictionary&
             schemesDict(dict_.subDict("sourceTerms").subDict("schemes"));
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
index e155abc..5423258 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
@@ -88,6 +88,14 @@ class cloudSolution
             //  carrier phase
             Switch coupled_;
 
+            //- Coupling regime of parcels and carrier phase:
+            //  0-way coupling: Only non-coupled forces are calculated
+            //  1-way coupling: Coupled and non-coupled forces are calculated,
+            //                  but no contribution to the source term
+            //  2-way coupling: Coupled and non-coupled forces are calculated
+            //                  and changes are contributed to the source term
+            scalar couplingRegime_;
+
             //- Flag to correct cell values with latest transfer information
             //  during the lagrangian timestep
             Switch cellValueSourceCorrection_;
@@ -174,8 +182,11 @@ public:
             //- Return const access to the coupled flag
             inline const Switch coupled() const;
 
-            //- Return non-const access to the coupled flag
-            inline Switch& coupled();
+            //- Return const access to the coupling regime
+            inline scalar couplingRegime() const;
+
+            //- Return non-const access to the coupling regime
+            inline scalar& couplingRegime();
 
             //- Return const access to the cell value correction flag
             inline const Switch cellValueSourceCorrection() const;
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
index 6280b97..f2f5d1b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
@@ -101,15 +101,21 @@ inline Foam::scalar Foam::cloudSolution::trackTime() const
 }
 
 
-inline Foam::Switch& Foam::cloudSolution::coupled()
+inline const Foam::Switch Foam::cloudSolution::coupled() const
 {
     return coupled_;
 }
 
 
-inline const Foam::Switch Foam::cloudSolution::coupled() const
+inline Foam::scalar& Foam::cloudSolution::couplingRegime()
 {
-    return coupled_;
+    return couplingRegime_;
+}
+
+
+inline Foam::scalar Foam::cloudSolution::couplingRegime() const
+{
+    return couplingRegime_;
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
index 7abbc9f..6e399ed 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
@@ -105,7 +105,7 @@ inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
     volScalarField& Yi
 ) const
 {
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         if (this->solution().semiImplicit("Yi"))
         {
@@ -180,7 +180,7 @@ Foam::ReactingCloud<CloudType>::Srho(const label i) const
         )
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         scalarField& rhoi = tRhoi();
         rhoi = rhoTrans_[i]/(this->db().time().deltaTValue()*this->mesh().V());
@@ -217,7 +217,7 @@ Foam::ReactingCloud<CloudType>::Srho() const
         )
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         scalarField& sourceField = trhoTrans();
         forAll(rhoTrans_, i)
@@ -236,7 +236,7 @@ template<class CloudType>
 inline Foam::tmp<Foam::fvScalarMatrix>
 Foam::ReactingCloud<CloudType>::Srho(volScalarField& rho) const
 {
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         tmp<volScalarField> trhoTrans
         (
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index e97b6cf..9f1600b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -51,7 +51,7 @@ void Foam::ThermoCloud<CloudType>::setModels()
         ).ptr()
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         this->subModelProperties().lookup("radiation") >> radiation_;
     }
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index ea1b79f..0b9a45a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -249,7 +249,7 @@ Foam::ThermoCloud<CloudType>::Sh(volScalarField& hs) const
             << max(hsCoeff()).value() << endl;
     }
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         if (this->solution().semiImplicit("h"))
         {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 5d9e4ef..d64b433 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -132,14 +132,11 @@ void Foam::KinematicParcel<ParcelType>::calc
 
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (td.cloud().solution().coupled())
-    {
     // Update momentum transfer
     td.cloud().UTrans()[cellI] += np0*dUTrans;
 
     // Update momentum transfer coefficient
     td.cloud().UCoeff()[cellI] += np0*Spu;
-    }
 }
 
 
@@ -166,7 +163,11 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     // Momentum source due to particle forces
     const parcelType& p = static_cast<const parcelType&>(*this);
-    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
+    forceSuSp Fcp(vector::zero, 0.0);
+    if (td.cloud().solution().couplingRegime() > 0.5)
+    {
+        Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
+    }
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
     const forceSuSp Feff = Fcp + Fncp;
     const scalar massEff = forces.massEff(p, mass);
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 72ad743..78e9548 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -345,8 +345,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     {
         td.keepParticle = false;
 
-        if (td.cloud().solution().coupled())
-        {
         scalar dm = np0*mass0;
 
         // Absorb parcel into carrier phase
@@ -373,7 +371,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T0, idG, idL, idS);
 
         td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
-        }
 
         return;
     }
@@ -420,8 +417,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // 4. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (td.cloud().solution().coupled())
-    {
     // Transfer mass lost to carrier mass, momentum and enthalpy sources
     forAll(YGas_, i)
     {
@@ -479,7 +474,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td.cloud().radT4()[cellI] += dt*np0*T4;
         td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
-    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 47a1310..96ad1ba 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -375,8 +375,6 @@ void Foam::ReactingParcel<ParcelType>::calc
     {
         td.keepParticle = false;
 
-        if (td.cloud().solution().coupled())
-        {
         scalar dm = np0*mass0;
 
         // Absorb parcel into carrier phase
@@ -392,7 +390,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         td.cloud().UTrans()[cellI] += dm*U0;
 
         td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
-        }
 
         return;
     }
@@ -438,8 +435,6 @@ void Foam::ReactingParcel<ParcelType>::calc
     // 4. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (td.cloud().solution().coupled())
-    {
     // Transfer mass lost to carrier mass, momentum and enthalpy sources
     forAll(dMass, i)
     {
@@ -469,7 +464,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         td.cloud().radT4()[cellI] += dt*np0*T4;
         td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
-    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 74ece70..e98f928 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -237,8 +237,6 @@ void Foam::ThermoParcel<ParcelType>::calc
 
     //  Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (td.cloud().solution().coupled())
-    {
     // Update momentum transfer
     td.cloud().UTrans()[cellI] += np0*dUTrans;
 
@@ -260,7 +258,6 @@ void Foam::ThermoParcel<ParcelType>::calc
         td.cloud().radT4()[cellI] += dt*np0*T4;
         td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
-    }
 }
 
 
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
index 22d9811..ed96db0 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
@@ -68,13 +68,13 @@ void Foam::SprayParcel<ParcelType>::calc
     const CompositionModel<reactingCloudType>& composition =
         td.cloud().composition();
 
-    const bool coupled = td.cloud().solution().coupled();
+    const scalar coupled = td.cloud().solution().couplingRegime();
 
     // check if parcel belongs to liquid core
     if (liquidCore() > 0.5)
     {
-        // liquid core parcels should not interact with the gas
-        td.cloud().solution().coupled() = false;
+        // no drag calculation for liquid core parcels
+        td.cloud().solution().couplingRegime() = 0.0;
     }
 
     // get old mixture composition
@@ -131,7 +131,7 @@ void Foam::SprayParcel<ParcelType>::calc
     }
 
     // restore coupled
-    td.cloud().solution().coupled() = coupled;
+    td.cloud().solution().couplingRegime() = coupled;
 }
Correct-coupling-of-parcels-and-carrier-phase_no-whitespace-changes.patch (14,918 bytes)   

dkxls

2013-09-08 14:49

reporter  

Correct-coupling-of-parcels-and-carrier-phase.patch (24,466 bytes)    
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index f3a8f16..936b50c 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -88,7 +88,7 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
-        if (solution_.coupled())
+        if (solution_.couplingRegime() > 1.5)
         {
             td.cloud().relaxSources(td.cloud().cloudCopy());
         }
@@ -99,7 +99,7 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
-        if (solution_.coupled())
+        if (solution_.couplingRegime() > 1.5)
         {
             td.cloud().scaleSources();
         }
@@ -165,7 +165,7 @@ template<class CloudType>
 template<class TrackData>
 void Foam::KinematicCloud<CloudType>::evolveCloud(TrackData& td)
 {
-    if (solution_.coupled())
+    if (solution_.couplingRegime() > 1.5)
     {
         td.cloud().resetSourceTerms();
     }
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 72cad26..2dbe962 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -555,7 +555,7 @@ Foam::KinematicCloud<CloudType>::SU(volVectorField& U) const
             << max(UCoeff()).value() << endl;
     }
 
-    if (solution_.coupled())
+    if (solution_.couplingRegime() > 1.5)
     {
         if (solution_.semiImplicit("U"))
         {
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
index fd7bf31..7fab1df 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
@@ -43,6 +43,7 @@ Foam::cloudSolution::cloudSolution
     iter_(1),
     trackTime_(0.0),
     coupled_(false),
+    couplingRegime_(1.0),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
     resetSourcesOnStartup_(true),
@@ -69,6 +70,7 @@ Foam::cloudSolution::cloudSolution
     iter_(cs.iter_),
     trackTime_(cs.trackTime_),
     coupled_(cs.coupled_),
+    couplingRegime_(cs.couplingRegime_),
     cellValueSourceCorrection_(cs.cellValueSourceCorrection_),
     maxTrackTime_(cs.maxTrackTime_),
     resetSourcesOnStartup_(cs.resetSourcesOnStartup_),
@@ -90,6 +92,7 @@ Foam::cloudSolution::cloudSolution
     iter_(0),
     trackTime_(0.0),
     coupled_(false),
+    couplingRegime_(0.0),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
     resetSourcesOnStartup_(false),
@@ -126,6 +129,8 @@ void Foam::cloudSolution::read()
 
     if (coupled_)
     {
+        couplingRegime_ = 2.0;
+
         const dictionary&
             schemesDict(dict_.subDict("sourceTerms").subDict("schemes"));
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
index e155abc..5423258 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
@@ -88,6 +88,14 @@ class cloudSolution
             //  carrier phase
             Switch coupled_;
 
+            //- Coupling regime of parcels and carrier phase:
+            //  0-way coupling: Only non-coupled forces are calculated
+            //  1-way coupling: Coupled and non-coupled forces are calculated,
+            //                  but no contribution to the source term
+            //  2-way coupling: Coupled and non-coupled forces are calculated
+            //                  and changes are contributed to the source term
+            scalar couplingRegime_;
+
             //- Flag to correct cell values with latest transfer information
             //  during the lagrangian timestep
             Switch cellValueSourceCorrection_;
@@ -174,8 +182,11 @@ public:
             //- Return const access to the coupled flag
             inline const Switch coupled() const;
 
-            //- Return non-const access to the coupled flag
-            inline Switch& coupled();
+            //- Return const access to the coupling regime
+            inline scalar couplingRegime() const;
+
+            //- Return non-const access to the coupling regime
+            inline scalar& couplingRegime();
 
             //- Return const access to the cell value correction flag
             inline const Switch cellValueSourceCorrection() const;
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
index 6280b97..f2f5d1b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
@@ -101,15 +101,21 @@ inline Foam::scalar Foam::cloudSolution::trackTime() const
 }
 
 
-inline Foam::Switch& Foam::cloudSolution::coupled()
+inline const Foam::Switch Foam::cloudSolution::coupled() const
 {
     return coupled_;
 }
 
 
-inline const Foam::Switch Foam::cloudSolution::coupled() const
+inline Foam::scalar& Foam::cloudSolution::couplingRegime()
 {
-    return coupled_;
+    return couplingRegime_;
+}
+
+
+inline Foam::scalar Foam::cloudSolution::couplingRegime() const
+{
+    return couplingRegime_;
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
index 7abbc9f..6e399ed 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
@@ -105,7 +105,7 @@ inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
     volScalarField& Yi
 ) const
 {
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         if (this->solution().semiImplicit("Yi"))
         {
@@ -180,7 +180,7 @@ Foam::ReactingCloud<CloudType>::Srho(const label i) const
         )
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         scalarField& rhoi = tRhoi();
         rhoi = rhoTrans_[i]/(this->db().time().deltaTValue()*this->mesh().V());
@@ -217,7 +217,7 @@ Foam::ReactingCloud<CloudType>::Srho() const
         )
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         scalarField& sourceField = trhoTrans();
         forAll(rhoTrans_, i)
@@ -236,7 +236,7 @@ template<class CloudType>
 inline Foam::tmp<Foam::fvScalarMatrix>
 Foam::ReactingCloud<CloudType>::Srho(volScalarField& rho) const
 {
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         tmp<volScalarField> trhoTrans
         (
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index e97b6cf..9f1600b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -51,7 +51,7 @@ void Foam::ThermoCloud<CloudType>::setModels()
         ).ptr()
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         this->subModelProperties().lookup("radiation") >> radiation_;
     }
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index ea1b79f..0b9a45a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -249,7 +249,7 @@ Foam::ThermoCloud<CloudType>::Sh(volScalarField& hs) const
             << max(hsCoeff()).value() << endl;
     }
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         if (this->solution().semiImplicit("h"))
         {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 5d9e4ef..d64b433 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -132,14 +132,11 @@ void Foam::KinematicParcel<ParcelType>::calc
 
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (td.cloud().solution().coupled())
-    {
-        // Update momentum transfer
-        td.cloud().UTrans()[cellI] += np0*dUTrans;
+    // Update momentum transfer
+    td.cloud().UTrans()[cellI] += np0*dUTrans;
 
-        // Update momentum transfer coefficient
-        td.cloud().UCoeff()[cellI] += np0*Spu;
-    }
+    // Update momentum transfer coefficient
+    td.cloud().UCoeff()[cellI] += np0*Spu;
 }
 
 
@@ -166,7 +163,11 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     // Momentum source due to particle forces
     const parcelType& p = static_cast<const parcelType&>(*this);
-    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
+    forceSuSp Fcp(vector::zero, 0.0);
+    if (td.cloud().solution().couplingRegime() > 0.5)
+    {
+        Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
+    }
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
     const forceSuSp Feff = Fcp + Fncp;
     const scalar massEff = forces.massEff(p, mass);
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 72ad743..78e9548 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -345,35 +345,32 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     {
         td.keepParticle = false;
 
-        if (td.cloud().solution().coupled())
+        scalar dm = np0*mass0;
+
+        // Absorb parcel into carrier phase
+        forAll(YGas_, i)
+        {
+            label gid = composition.localToGlobalCarrierId(GAS, i);
+            td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
+        }
+        forAll(YLiquid_, i)
         {
-            scalar dm = np0*mass0;
-
-            // Absorb parcel into carrier phase
-            forAll(YGas_, i)
-            {
-                label gid = composition.localToGlobalCarrierId(GAS, i);
-                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
-            }
-            forAll(YLiquid_, i)
-            {
-                label gid = composition.localToGlobalCarrierId(LIQ, i);
-                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
-            }
+            label gid = composition.localToGlobalCarrierId(LIQ, i);
+            td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
+        }
 /*
-            // No mapping between solid components and carrier phase
-            forAll(YSolid_, i)
-            {
-                label gid = composition.localToGlobalCarrierId(SLD, i);
-                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
-            }
+        // No mapping between solid components and carrier phase
+        forAll(YSolid_, i)
+        {
+            label gid = composition.localToGlobalCarrierId(SLD, i);
+            td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
+        }
 */
-            td.cloud().UTrans()[cellI] += dm*U0;
+        td.cloud().UTrans()[cellI] += dm*U0;
 
-            td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T0, idG, idL, idS);
+        td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T0, idG, idL, idS);
 
-            td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
-        }
+        td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
 
         return;
     }
@@ -420,65 +417,62 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // 4. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (td.cloud().solution().coupled())
+    // Transfer mass lost to carrier mass, momentum and enthalpy sources
+    forAll(YGas_, i)
     {
-        // Transfer mass lost to carrier mass, momentum and enthalpy sources
-        forAll(YGas_, i)
-        {
-            scalar dm = np0*dMassGas[i];
-            label gid = composition.localToGlobalCarrierId(GAS, i);
-            scalar hs = composition.carrier().Hs(gid, pc, T0);
-            td.cloud().rhoTrans(gid)[cellI] += dm;
-            td.cloud().UTrans()[cellI] += dm*U0;
-            td.cloud().hsTrans()[cellI] += dm*hs;
-        }
-        forAll(YLiquid_, i)
-        {
-            scalar dm = np0*dMassLiquid[i];
-            label gid = composition.localToGlobalCarrierId(LIQ, i);
-            scalar hs = composition.carrier().Hs(gid, pc, T0);
-            td.cloud().rhoTrans(gid)[cellI] += dm;
-            td.cloud().UTrans()[cellI] += dm*U0;
-            td.cloud().hsTrans()[cellI] += dm*hs;
-        }
+        scalar dm = np0*dMassGas[i];
+        label gid = composition.localToGlobalCarrierId(GAS, i);
+        scalar hs = composition.carrier().Hs(gid, pc, T0);
+        td.cloud().rhoTrans(gid)[cellI] += dm;
+        td.cloud().UTrans()[cellI] += dm*U0;
+        td.cloud().hsTrans()[cellI] += dm*hs;
+    }
+    forAll(YLiquid_, i)
+    {
+        scalar dm = np0*dMassLiquid[i];
+        label gid = composition.localToGlobalCarrierId(LIQ, i);
+        scalar hs = composition.carrier().Hs(gid, pc, T0);
+        td.cloud().rhoTrans(gid)[cellI] += dm;
+        td.cloud().UTrans()[cellI] += dm*U0;
+        td.cloud().hsTrans()[cellI] += dm*hs;
+    }
 /*
-        // No mapping between solid components and carrier phase
-        forAll(YSolid_, i)
-        {
-            scalar dm = np0*dMassSolid[i];
-            label gid = composition.localToGlobalCarrierId(SLD, i);
-            scalar hs = composition.carrier().Hs(gid, pc, T0);
-            td.cloud().rhoTrans(gid)[cellI] += dm;
-            td.cloud().UTrans()[cellI] += dm*U0;
-            td.cloud().hsTrans()[cellI] += dm*hs;
-        }
+    // No mapping between solid components and carrier phase
+    forAll(YSolid_, i)
+    {
+        scalar dm = np0*dMassSolid[i];
+        label gid = composition.localToGlobalCarrierId(SLD, i);
+        scalar hs = composition.carrier().Hs(gid, pc, T0);
+        td.cloud().rhoTrans(gid)[cellI] += dm;
+        td.cloud().UTrans()[cellI] += dm*U0;
+        td.cloud().hsTrans()[cellI] += dm*hs;
+    }
 */
-        forAll(dMassSRCarrier, i)
-        {
-            scalar dm = np0*dMassSRCarrier[i];
-            scalar hs = composition.carrier().Hs(i, pc, T0);
-            td.cloud().rhoTrans(i)[cellI] += dm;
-            td.cloud().UTrans()[cellI] += dm*U0;
-            td.cloud().hsTrans()[cellI] += dm*hs;
-        }
+    forAll(dMassSRCarrier, i)
+    {
+        scalar dm = np0*dMassSRCarrier[i];
+        scalar hs = composition.carrier().Hs(i, pc, T0);
+        td.cloud().rhoTrans(i)[cellI] += dm;
+        td.cloud().UTrans()[cellI] += dm*U0;
+        td.cloud().hsTrans()[cellI] += dm*hs;
+    }
 
-        // Update momentum transfer
-        td.cloud().UTrans()[cellI] += np0*dUTrans;
-        td.cloud().UCoeff()[cellI] += np0*Spu;
+    // Update momentum transfer
+    td.cloud().UTrans()[cellI] += np0*dUTrans;
+    td.cloud().UCoeff()[cellI] += np0*Spu;
 
-        // Update sensible enthalpy transfer
-        td.cloud().hsTrans()[cellI] += np0*dhsTrans;
-        td.cloud().hsCoeff()[cellI] += np0*Sph;
+    // Update sensible enthalpy transfer
+    td.cloud().hsTrans()[cellI] += np0*dhsTrans;
+    td.cloud().hsCoeff()[cellI] += np0*Sph;
 
-        // Update radiation fields
-        if (td.cloud().radiation())
-        {
-            const scalar ap = this->areaP();
-            const scalar T4 = pow4(this->T_);
-            td.cloud().radAreaP()[cellI] += dt*np0*ap;
-            td.cloud().radT4()[cellI] += dt*np0*T4;
-            td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
-        }
+    // Update radiation fields
+    if (td.cloud().radiation())
+    {
+        const scalar ap = this->areaP();
+        const scalar T4 = pow4(this->T_);
+        td.cloud().radAreaP()[cellI] += dt*np0*ap;
+        td.cloud().radT4()[cellI] += dt*np0*T4;
+        td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
 }
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 47a1310..96ad1ba 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -375,24 +375,21 @@ void Foam::ReactingParcel<ParcelType>::calc
     {
         td.keepParticle = false;
 
-        if (td.cloud().solution().coupled())
-        {
-            scalar dm = np0*mass0;
-
-            // Absorb parcel into carrier phase
-            forAll(Y_, i)
-            {
-                scalar dmi = dm*Y_[i];
-                label gid = composition.localToGlobalCarrierId(0, i);
-                scalar hs = composition.carrier().Hs(gid, pc_, T0);
+        scalar dm = np0*mass0;
 
-                td.cloud().rhoTrans(gid)[cellI] += dmi;
-                td.cloud().hsTrans()[cellI] += dmi*hs;
-            }
-            td.cloud().UTrans()[cellI] += dm*U0;
+        // Absorb parcel into carrier phase
+        forAll(Y_, i)
+        {
+            scalar dmi = dm*Y_[i];
+            label gid = composition.localToGlobalCarrierId(0, i);
+            scalar hs = composition.carrier().Hs(gid, pc_, T0);
 
-            td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
+            td.cloud().rhoTrans(gid)[cellI] += dmi;
+            td.cloud().hsTrans()[cellI] += dmi*hs;
         }
+        td.cloud().UTrans()[cellI] += dm*U0;
+
+        td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
 
         return;
     }
@@ -438,37 +435,34 @@ void Foam::ReactingParcel<ParcelType>::calc
     // 4. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (td.cloud().solution().coupled())
+    // Transfer mass lost to carrier mass, momentum and enthalpy sources
+    forAll(dMass, i)
     {
-        // Transfer mass lost to carrier mass, momentum and enthalpy sources
-        forAll(dMass, i)
-        {
-            scalar dm = np0*dMass[i];
-            label gid = composition.localToGlobalCarrierId(0, i);
-            scalar hs = composition.carrier().Hs(gid, pc_, T0);
+        scalar dm = np0*dMass[i];
+        label gid = composition.localToGlobalCarrierId(0, i);
+        scalar hs = composition.carrier().Hs(gid, pc_, T0);
 
-            td.cloud().rhoTrans(gid)[cellI] += dm;
-            td.cloud().UTrans()[cellI] += dm*U0;
-            td.cloud().hsTrans()[cellI] += dm*hs;
-        }
+        td.cloud().rhoTrans(gid)[cellI] += dm;
+        td.cloud().UTrans()[cellI] += dm*U0;
+        td.cloud().hsTrans()[cellI] += dm*hs;
+    }
 
-        // Update momentum transfer
-        td.cloud().UTrans()[cellI] += np0*dUTrans;
-        td.cloud().UCoeff()[cellI] += np0*Spu;
+    // Update momentum transfer
+    td.cloud().UTrans()[cellI] += np0*dUTrans;
+    td.cloud().UCoeff()[cellI] += np0*Spu;
 
-        // Update sensible enthalpy transfer
-        td.cloud().hsTrans()[cellI] += np0*dhsTrans;
-        td.cloud().hsCoeff()[cellI] += np0*Sph;
+    // Update sensible enthalpy transfer
+    td.cloud().hsTrans()[cellI] += np0*dhsTrans;
+    td.cloud().hsCoeff()[cellI] += np0*Sph;
 
-        // Update radiation fields
-        if (td.cloud().radiation())
-        {
-            const scalar ap = this->areaP();
-            const scalar T4 = pow4(this->T_);
-            td.cloud().radAreaP()[cellI] += dt*np0*ap;
-            td.cloud().radT4()[cellI] += dt*np0*T4;
-            td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
-        }
+    // Update radiation fields
+    if (td.cloud().radiation())
+    {
+        const scalar ap = this->areaP();
+        const scalar T4 = pow4(this->T_);
+        td.cloud().radAreaP()[cellI] += dt*np0*ap;
+        td.cloud().radT4()[cellI] += dt*np0*T4;
+        td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
 }
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 74ece70..e98f928 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -237,29 +237,26 @@ void Foam::ThermoParcel<ParcelType>::calc
 
     //  Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (td.cloud().solution().coupled())
-    {
-        // Update momentum transfer
-        td.cloud().UTrans()[cellI] += np0*dUTrans;
+    // Update momentum transfer
+    td.cloud().UTrans()[cellI] += np0*dUTrans;
 
-        // Update momentum transfer coefficient
-        td.cloud().UCoeff()[cellI] += np0*Spu;
+    // Update momentum transfer coefficient
+    td.cloud().UCoeff()[cellI] += np0*Spu;
 
-        // Update sensible enthalpy transfer
-        td.cloud().hsTrans()[cellI] += np0*dhsTrans;
+    // Update sensible enthalpy transfer
+    td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 
-        // Update sensible enthalpy coefficient
-        td.cloud().hsCoeff()[cellI] += np0*Sph;
+    // Update sensible enthalpy coefficient
+    td.cloud().hsCoeff()[cellI] += np0*Sph;
 
-        // Update radiation fields
-        if (td.cloud().radiation())
-        {
-            const scalar ap = this->areaP();
-            const scalar T4 = pow4(this->T_);
-            td.cloud().radAreaP()[cellI] += dt*np0*ap;
-            td.cloud().radT4()[cellI] += dt*np0*T4;
-            td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
-        }
+    // Update radiation fields
+    if (td.cloud().radiation())
+    {
+        const scalar ap = this->areaP();
+        const scalar T4 = pow4(this->T_);
+        td.cloud().radAreaP()[cellI] += dt*np0*ap;
+        td.cloud().radT4()[cellI] += dt*np0*T4;
+        td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
 }
 
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
index 22d9811..ed96db0 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
@@ -68,13 +68,13 @@ void Foam::SprayParcel<ParcelType>::calc
     const CompositionModel<reactingCloudType>& composition =
         td.cloud().composition();
 
-    const bool coupled = td.cloud().solution().coupled();
+    const scalar coupled = td.cloud().solution().couplingRegime();
 
     // check if parcel belongs to liquid core
     if (liquidCore() > 0.5)
     {
-        // liquid core parcels should not interact with the gas
-        td.cloud().solution().coupled() = false;
+        // no drag calculation for liquid core parcels
+        td.cloud().solution().couplingRegime() = 0.0;
     }
 
     // get old mixture composition
@@ -131,7 +131,7 @@ void Foam::SprayParcel<ParcelType>::calc
     }
 
     // restore coupled
-    td.cloud().solution().coupled() = coupled;
+    td.cloud().solution().couplingRegime() = coupled;
 }
Correct-coupling-of-parcels-and-carrier-phase.patch (24,466 bytes)   

dkxls

2013-09-09 09:55

reporter   ~0002474

Some more reference for the no drag assumption with respect to (engine) sprays. I personally think the paper by Beale (1999) gives a good overview, but also Xin (1998) describes the topic well.

Beale, J.C.; Reitz, R.D., "Modeling Spray Atomization with the Kelvin-Helmholtz/Rayleigh-Taylor Hybrid Model," Atomization and Sprays, Vol. 9, pp. 623-650, 1999

Xin, J.; Ricart, L.; Reitz, R.D., "Computer Modeling of Diesel Spray Atomization and Combustion," Combust. Sci. and Tech. V. 137, 1-6, p. 171, 1998

user2

2013-09-18 14:13

  ~0002499

Thanks for the report - I've implemented an alternative fix in commit 949fc50

Issue History

Date Modified Username Field Change
2013-09-08 14:49 dkxls New Issue
2013-09-08 14:49 dkxls File Added: Correct-coupling-of-parcels-and-carrier-phase_no-whitespace-changes.patch
2013-09-08 14:49 dkxls File Added: Correct-coupling-of-parcels-and-carrier-phase.patch
2013-09-09 09:55 dkxls Note Added: 0002474
2013-09-18 14:13 user2 Note Added: 0002499
2013-09-18 14:13 user2 Status new => resolved
2013-09-18 14:13 user2 Fixed in Version => 2.2.x
2013-09-18 14:13 user2 Resolution open => fixed
2013-09-18 14:13 user2 Assigned To => user2