From 9812d957c9ae91e40f1d207ce54daea66b7db50e Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Mon, 30 Jan 2017 16:37:27 +0000
Subject: [PATCH] ThermalPhaseChangePhaseSystem: Improved robustness
Patch contributed by Juho Peltola, VTT.
Resolves patch request https://bugs.openfoam.org/view.php?id=2443
---
.../ThermalPhaseChangePhaseSystem.C | 59 ++++++++++---------
1 file changed, 30 insertions(+), 29 deletions(-)
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index b33289a0b73..4335bace8a9 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -375,6 +375,18 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
+ volScalarField L
+ (
+ min
+ (
+ (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
+ - (neg(iDmdt)*he1 + pos(iDmdt)*hef1),
+ 0.3*mag(hef2 - hef1)
+ )
+ );
+
+ volScalarField iDmdtNew(iDmdt);
+
if (massTransfer_ )
{
volScalarField H1
@@ -389,28 +401,13 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
Tf = saturationModel_->Tsat(phase1.thermo().p());
- scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
+ iDmdtNew =
+ (H1*(Tf - T1) + H2*(Tf - T2))/L;
- iDmdt =
- (1 - iDmdtRelax)*iDmdt
- + iDmdtRelax*(H1*(Tf - T1) + H2*(Tf - T2))
- /min
- (
- (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
- - (neg(iDmdt)*he1 + pos(iDmdt)*hef1),
- 0.3*mag(hef2 - hef1)
- );
-
- Info<< "iDmdt." << pair.name()
- << ": min = " << min(iDmdt.primitiveField())
- << ", mean = " << average(iDmdt.primitiveField())
- << ", max = " << max(iDmdt.primitiveField())
- << ", integral = " << fvc::domainIntegrate(iDmdt).value()
- << endl;
}
else
{
- iDmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
+ iDmdtNew == dimensionedScalar("0",dmdt.dimensions(), 0);
}
volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
@@ -423,16 +420,7 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
H2.boundaryFieldRef() =
max(H2.boundaryField(), phase2.boundaryField()*HLimit);
- volScalarField mDotL
- (
- iDmdt*
- (
- (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
- - (neg(iDmdt)*he1 + pos(iDmdt)*hef1)
- )
- );
-
- Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
+ Tf = (H1*T1 + H2*T2 + iDmdtNew*L)/(H1 + H2);
Info<< "Tf." << pair.name()
<< ": min = " << min(Tf.primitiveField())
@@ -440,6 +428,19 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
<< ", max = " << max(Tf.primitiveField())
<< endl;
+ scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
+ iDmdt = (1 - iDmdtRelax)*iDmdt + iDmdtRelax*iDmdtNew;
+
+ if (massTransfer_ )
+ {
+ Info<< "iDmdt." << pair.name()
+ << ": min = " << min(iDmdt.primitiveField())
+ << ", mean = " << average(iDmdt.primitiveField())
+ << ", max = " << max(iDmdt.primitiveField())
+ << ", integral = " << fvc::domainIntegrate(iDmdt).value()
+ << endl;
+ }
+
// Accumulate dmdt contributions from boundaries
volScalarField wDmdt
(
--
GitLab