From 745fc42dee98fcee5e3d685dc26e60a019b828ce Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Wed, 20 Apr 2022 17:22:43 +0100 Subject: [PATCH] ENH: turbulentDFSEMInlet: improve Reynolds-stress realizability checks ENH: turbulentDFSEMInlet: new realizability checking function for scalar fields --- .../turbulentDFSEMInletFvPatchVectorField.C | 105 +++++++++--------- .../turbulentDFSEMInletFvPatchVectorField.H | 12 +- .../turbulentDigitalFilterInletFvPatchField.C | 6 + 3 files changed, 71 insertions(+), 52 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C index ed9d48e091..758ecc99d2 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C @@ -811,71 +811,76 @@ turbulentDFSEMInletFvPatchVectorField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -bool Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses +void Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses ( - const symmTensorField& Rf + const symmTensorField& R ) { - // Perform checks of the stress tensor based on Cholesky decomposition - // constraints + constexpr label maxDiffs = 5; + label nDiffs = 0; - forAll(Rf, facei) + // (S:Eq. 4a-4c) + forAll(R, i) { - const symmTensor& R = Rf[facei]; + bool diff = false; - if (R.xx() <= 0) + if (maxDiffs < nDiffs) { - FatalErrorInFunction - << "Reynolds stress " << R << " at face " << facei - << " does not obey the constraint: R_xx > 0" - << exit(FatalError); + Info<< "More than " << maxDiffs << " times" + << " Reynolds-stress realizability checks failed." + << " Skipping further comparisons." << endl; + return; } - const scalar a_xx = sqrt(R.xx()); + const symmTensor& r = R[i]; - const scalar a_xy = R.xy()/a_xx; - - const scalar a_yy_2 = R.yy() - sqr(a_xy); - - if (a_yy_2 < 0) + if (r.xx() < 0) { - FatalErrorInFunction - << "Reynolds stress " << R << " at face " << facei - << " leads to an invalid Cholesky decomposition due to the " - << "constraint R_yy - sqr(a_xy) >= 0" - << exit(FatalError); - } - - const scalar a_yy = Foam::sqrt(a_yy_2); - - const scalar a_xz = R.xz()/a_xx; - - const scalar a_yz = (R.yz() - a_xy*a_xz)/a_yy; - - const scalar a_zz_2 = R.zz() - sqr(a_xz) - sqr(a_yz); - - if (a_zz_2 < 0) - { - FatalErrorInFunction - << "Reynolds stress " << R << " at face " << facei - << " leads to an invalid Cholesky decomposition due to the " - << "constraint R_zz - sqr(a_xz) - sqr(a_yz) >= 0" - << exit(FatalError); - } - - const scalar a_zz = Foam::sqrt(a_zz_2); - - if (debug) - { - Pout<< "R: " << R - << " a_xx:" << a_xx << " a_xy:" << a_xy << " a_xz:" << a_xy - << " a_yy:" << a_yy << " a_yz:" << a_yz - << " a_zz:" << a_zz + WarningInFunction + << "Reynolds stress " << r << " at index " << i + << " does not obey the constraint: Rxx >= 0" << endl; + diff = true; + } + + if ((r.xx()*r.yy() - sqr(r.xy())) < 0) + { + WarningInFunction + << "Reynolds stress " << r << " at index " << i + << " does not obey the constraint: Rxx*Ryy - sqr(Rxy) >= 0" + << endl; + diff = true; + } + + if (det(r) < 0) + { + WarningInFunction + << "Reynolds stress " << r << " at index " << i + << " does not obey the constraint: det(R) >= 0" + << endl; + diff = true; + } + + if (diff) + { + ++nDiffs; } } +} - return true; + +void Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses +( + const scalarField& R +) +{ + if (min(R) <= 0) + { + FatalErrorInFunction + << "Reynolds stresses contain at least one " + << "nonpositive element. min(R) = " << min(R) + << exit(FatalError); + } } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H index 624f3354bf..25a6c57105 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H @@ -361,8 +361,16 @@ public: // Member Functions - //- Return true if input Reynold stresses are valid - static bool checkStresses(const symmTensorField& Rf); + //- Check if input Reynolds stresses are valid + // Realizability conditions (tag:S): + // Schumann, U. (1977). + // Realizability of Reynolds‐stress turbulence models. + // The Physics of Fluids, 20(5), 721-725. + // DOI:10.1063/1.861942 + static void checkStresses(const symmTensorField& R); + + //- Check if input Reynolds stresses are valid + static void checkStresses(const scalarField& R); // Mapping diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C index b5e44b4323..80331851cb 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C @@ -28,6 +28,7 @@ License #include "turbulentDigitalFilterInletFvPatchField.H" #include "addToRunTimeSelectionTable.H" #include "faceAreaWeightAMI.H" +#include "turbulentDFSEMInletFvPatchVectorField.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -284,6 +285,11 @@ turbulentDigitalFilterInletFvPatchField << "supported by the digital filter method." << endl; } + + const scalar t = this->db().time().timeOutputValue(); + const Field R(Rptr_->value(t)); + + turbulentDFSEMInletFvPatchVectorField::checkStresses(R); }