From b6c5fc0ada8da47eb1464e841517de25b486e0e4 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Tue, 5 Feb 2019 14:07:59 +0000
Subject: [PATCH] ENH: Added checks for input Reynolds stresses based on Lund
coefficient constraints. Fixes #1124
---
.../turbulentDFSEMInletFvPatchVectorField.C | 74 +++++++++++++++++--
.../turbulentDFSEMInletFvPatchVectorField.H | 9 ++-
2 files changed, 75 insertions(+), 8 deletions(-)
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index 54b5bf879d3..f8b218ecadf 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
- \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd.
+ \\/ M anipulation | Copyright (C) 2016-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -854,6 +854,8 @@ turbulentDFSEMInletFvPatchVectorField
{
eddy::debug = debug;
+ checkStresses(R_);
+
// Set UMean as patch area average value
UMean_ = gSum(U_*patch().magSf())/(gSum(patch().magSf()) + ROOTVSMALL);
}
@@ -946,14 +948,74 @@ turbulentDFSEMInletFvPatchVectorField
{}
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-Foam::turbulentDFSEMInletFvPatchVectorField::~
-turbulentDFSEMInletFvPatchVectorField()
-{}
+bool Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
+(
+ const symmTensorField& Rf
+)
+{
+ // Perform checks based on constraints imposed by Lund coefficients
+ // a11..a33
+
+ for (const symmTensor& R : Rf)
+ {
+ if (R.xx() <= 0)
+ {
+ FatalErrorInFunction
+ << "Reynold stress " << R
+ << " does not obey the Lund coefficient constraint: "
+ << "R.xx() > 0"
+ << exit(FatalError);
+ }
+ scalar a11 = sqrt(R.xx());
+
+ scalar a21 = R.xy()/a11;
+
+ scalar a22_2 = R.yy() - sqr(a21);
+
+ if (a22_2 < 0)
+ {
+ FatalErrorInFunction
+ << "Reynold stress " << R
+ << " does not obey the Lund coefficient constraint: "
+ << "R.yy() - sqr(a21) >= 0"
+ << exit(FatalError);
+ }
+
+ scalar a22 = Foam::sqrt(a22_2);
+
+ scalar a31 = R.xz()/a11;
+
+ scalar a32 = (R.yz() - a21*a31)*a22;
+
+ scalar a33_2 = R.zz() - sqr(a31) - sqr(a32);
+
+ if (a33_2 < 0)
+ {
+ FatalErrorInFunction
+ << "Reynold stress " << R
+ << " does not obey the Lund coefficient constraint: "
+ << "R.zz() - sqr(a31) - sqr(a32) >= 0"
+ << exit(FatalError);
+ }
+
+ scalar a33 = Foam::sqrt(a33_2);
+
+ if (debug)
+ {
+ Pout<< "R: " << R
+ << " a11:" << a11 << " a21:" << a21 << " a31:" << a31
+ << " a21:" << a22 << " a32:" << a32
+ << " a33:" << a33
+ << endl;
+ }
+ }
+
+ return true;
+}
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::turbulentDFSEMInletFvPatchVectorField::autoMap
(
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
index 62bf5065fa9..8c96a96c591 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
- \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
+ \\/ M anipulation | Copyright (C) 2016-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -329,11 +329,16 @@ public:
);
}
- virtual ~turbulentDFSEMInletFvPatchVectorField();
+
+ //- Destructor
+ virtual ~turbulentDFSEMInletFvPatchVectorField() = default;
// Member functions
+ //- Helper function to check that Reynold stresses are valid
+ static bool checkStresses(const symmTensorField& Rf);
+
// Mapping functions
--
GitLab