View Issue Details

IDProjectCategoryView StatusLast Update
0003148OpenFOAMBugpublic2019-01-04 13:56
Reporteroertel59Assigned To 
PrioritynormalSeverityminorReproducibilityalways
Status newResolutionopen 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Product Versiondev 
Fixed in Version 
Summary0003148: reactingMultiphaseEulerFoam: Dictionary entry order for the phases involved affects results
DescriptionAt least for reactingMultiphaseEulerFoam, it doesn't feel unreasonable to expect that the results are not affected by the order in which the phases involved are specified. I.e. in the phaseProperties dictionary, a change from

phases (air water)

to

phases (water air)

should yield identical results. This is only the case though if the simulation is incompressible, as evident from tests with a laminar bubbleColumn case using different entry orders and combinations of equationOfState. The original solver gives

original/rhoConst/log.diff.alpha.airMean is empty
original/rhoConstPerfectFluid/log.diff.alpha.airMean contains entries
original/perfectGasRhoConst/log.diff.alpha.airMean contains entries
original/perfectGasPerfectFluid/log.diff.alpha.airMean contains entries

The attached patch separates the source term assembly and the solution of the alphaEqns. If only the changes to multiphaseSystem are considered, I get:

multiphaseSystem/rhoConst/log.diff.alpha.airMean is empty
multiphaseSystem/rhoConstPerfectFluid/log.diff.alpha.airMean is empty
multiphaseSystem/perfectGasRhoConst/log.diff.alpha.airMean is empty
multiphaseSystem/perfectGasPerfectFluid/log.diff.alpha.airMean contains entries

Apparently, the order how the pEqnComps are added to the pEqn plays a role as well (round-off behaviour?) If the pEqnComps are first summed up and then added to the pEqn as a whole, this yields:

multiphaseSystemAndpEqn/rhoConst/log.diff.alpha.airMean is empty
multiphaseSystemAndpEqn/rhoConstPerfectFluid/log.diff.alpha.airMean is empty
multiphaseSystemAndpEqn/perfectGasRhoConst/log.diff.alpha.airMean is empty
multiphaseSystemAndpEqn/perfectGasPerfectFluid/log.diff.alpha.airMean is empty.

Admittedly, I am not so sure though whether the change to pEqn is convincing enough to be considered.
TagsNo tags attached.

Activities

oertel59

2019-01-04 12:49

reporter  

original.tar.gz (52,001 bytes)
multiphaseSystemAndpEqn.tar.gz (17,890 bytes)
multiphaseSystem.tar.gz (29,785 bytes)
0001-reactingMultiphaseEulerFoam-Fixes-a-bug-where-input-.patch (6,719 bytes)
From d00c51e3b73636fae3a85c9f81942098f6e07629 Mon Sep 17 00:00:00 2001
From: "Lehnigk, Ronald (FWDC) - 10659" <r.lehnigk@hzdr.de>
Date: Fri, 4 Jan 2019 12:54:34 +0100
Subject: [PATCH] reactingMultiphaseEulerFoam: Fixes a bug where input order
 affects results

Prior to this commit, the order in which the phases involved are specified
through the entry "phases" in the phaseProperties dictionary affected the
results in case of an compressible simulation.

Patch contributed by Institute of Fluid Dynamics,
Helmholtz-Zentrum Dresden - Rossendorf (HZDR)
---
 .../multiphaseSystem/multiphaseSystem.C       | 67 ++++++++++++-------
 .../reactingMultiphaseEulerFoam/pU/pEqn.H     | 17 +++--
 2 files changed, 55 insertions(+), 29 deletions(-)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 62729c9b0..63f679db6 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,32 +202,38 @@ void Foam::multiphaseSystem::solveAlphas()
     }
     MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs);
 
-    // Solve for the moving phase alphas
+    // Compute alpha equation source terms for the moving phases
+    PtrList<volScalarField::Internal> Sps(movingPhases().size());
+    PtrList<volScalarField::Internal> Sus(movingPhases().size());
     forAll(movingPhases(), movingPhasei)
     {
         phaseModel& phase = movingPhases()[movingPhasei];
         volScalarField& alpha = phase;
 
-        surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];
-        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
-        phase.correctInflowOutflow(alphaPhi);
-
-        volScalarField::Internal Sp
+        Sps.set
         (
-            IOobject
+            phase.index(),
+            new volScalarField
             (
-                "Sp",
-                mesh_.time().timeName(),
-                mesh_
-            ),
-            mesh_,
-            dimensionedScalar(dimless/dimTime, 0)
+                IOobject
+                (
+                    "Sp",
+                    mesh_.time().timeName(),
+                    mesh_
+                ),
+                mesh_,
+                dimensionedScalar(dimless/dimTime, 0)
+            )
         );
 
-        volScalarField::Internal Su
+        Sus.set
         (
-            "Su",
-            min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U()))
+            phase.index(),
+            new volScalarField
+            (
+                "Su",
+                min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U()))
+            )
         );
 
         if (phase.divU().valid())
@@ -238,12 +244,12 @@ void Foam::multiphaseSystem::solveAlphas()
             {
                 if (dgdt[celli] > 0.0)
                 {
-                    Sp[celli] -= dgdt[celli];
-                    Su[celli] += dgdt[celli];
+                    Sps[phase.index()][celli] -= dgdt[celli];
+                    Sus[phase.index()][celli] += dgdt[celli];
                 }
                 else if (dgdt[celli] < 0.0)
                 {
-                    Sp[celli] +=
+                    Sps[phase.index()][celli] +=
                         dgdt[celli]
                        *(1 - alpha[celli])/max(alpha[celli], 1e-4);
                 }
@@ -265,29 +271,40 @@ void Foam::multiphaseSystem::solveAlphas()
                 {
                     if (dgdt2[celli] < 0.0)
                     {
-                        Sp[celli] +=
+                        Sps[phase.index()][celli] +=
                             dgdt2[celli]
                            *(1 - alpha2[celli])/max(alpha2[celli], 1e-4);
 
-                        Su[celli] -=
+                        Sus[phase.index()][celli] -=
                             dgdt2[celli]
                            *alpha[celli]/max(alpha2[celli], 1e-4);
                     }
                     else if (dgdt2[celli] > 0.0)
                     {
-                        Sp[celli] -= dgdt2[celli];
+                        Sps[phase.index()][celli] -= dgdt2[celli];
                     }
                 }
             }
         }
+    }
+
+    // Solve for the moving phase alphas
+    forAll(movingPhases(), movingPhasei)
+    {
+        phaseModel& phase = movingPhases()[movingPhasei];
+        volScalarField& alpha = phase;
+
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
+        phase.correctInflowOutflow(alphaPhi);
 
         MULES::explicitSolve
         (
             geometricOneField(),
             alpha,
             alphaPhi,
-            Sp,
-            Su
+            Sps[phase.index()],
+            Sus[phase.index()]
         );
 
         phase.alphaPhiRef() = alphaPhi;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index f50de4b7b..a92bb2811 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -361,13 +361,22 @@ while (pimple.correct())
 
         {
             fvScalarMatrix pEqn(pEqnIncomp);
-
-            forAll(phases, phasei)
             {
-                if (pEqnComps.set(phasei))
+                fvScalarMatrix pEqnComp
+                (
+                    p_rgh,
+                    pEqn.dimensions()
+                );
+
+                forAll(phases, phasei)
                 {
-                    pEqn += pEqnComps[phasei];
+                    if (pEqnComps.set(phasei))
+                    {
+                        pEqnComp += pEqnComps[phasei];
+                    }
                 }
+
+                pEqn += pEqnComp;
             }
 
             pEqn.solve();
-- 
2.20.1

henry

2019-01-04 13:56

manager   ~0010252

Currently the phase fractions are solved in the order specified in the "phases" list and the latest solution of each is used so the order is likely to affect the solution unless an outer iteration is converged.
It is possible to cache the intermediate results to ensure the order does not matter but it will be slightly slower (less convergent?), use more memory and add significant complexity, I guess this is what you have done in the patch you provide, but is it worth it? If you loop to convergence do you still get differences? How large are the difference? Do they matter if they are less than the numerical error of the solution?

Issue History

Date Modified Username Field Change
2019-01-04 12:49 oertel59 New Issue
2019-01-04 12:49 oertel59 File Added: original.tar.gz
2019-01-04 12:49 oertel59 File Added: multiphaseSystemAndpEqn.tar.gz
2019-01-04 12:49 oertel59 File Added: multiphaseSystem.tar.gz
2019-01-04 12:49 oertel59 File Added: 0001-reactingMultiphaseEulerFoam-Fixes-a-bug-where-input-.patch
2019-01-04 13:56 henry Note Added: 0010252