momentumTransportModels::linearViscousStress::divDevTau: Refactored for flux()

The divDevTau matrix now caches the explicit correction flux if fluxRequired is
true for U so that calling flux() on the matrix returns the flux of the complete
stress, not just the implicit part.

src/MomentumTransportModels/compressible/compressibleMomentumTransportModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -138,7 +138,7 @@ public:
         //- Return the effective stress tensor including the laminar stress
         virtual tmp<volSymmTensorField> devTau() const = 0;
 
-        //- Return the source term for the momentum equation
+        //- Return the stress matrix for the momentum equation
         virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U) const = 0;
 };
 

src/MomentumTransportModels/momentumTransportModels/linearViscousStress/linearViscousStress.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,7 +70,7 @@ Foam::linearViscousStress<BasicMomentumTransportModel>::devTau() const
 {
     return volSymmTensorField::New
     (
-        IOobject::groupName("devTau", this->alphaRhoPhi_.group()),
+        this->groupName("devTau"),
         (-(this->alpha_*this->rho_*this->nuEff()))
        *dev(twoSymm(fvc::grad(this->U_)))
     );
@@ -86,8 +86,12 @@ Foam::linearViscousStress<BasicMomentumTransportModel>::divDevTau
 {
     return
     (
-      - fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
       - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
+      + this->divDevTauCorr
+        (
+            -(this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))),
+            U
+        )
     );
 }
 
@@ -102,8 +106,12 @@ Foam::linearViscousStress<BasicMomentumTransportModel>::divDevTau
 {
     return
     (
-      - fvc::div((this->alpha_*rho*this->nuEff())*dev2(T(fvc::grad(U))))
       - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
+      + this->divDevTauCorr
+        (
+            -(this->alpha_*rho*this->nuEff())*dev2(T(fvc::grad(U))),
+            U
+        )
     );
 }
 

src/MomentumTransportModels/momentumTransportModels/momentumTransportModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,8 @@ License
 #include "surfaceFields.H"
 #include "wallFvPatch.H"
 #include "nearWallDist.H"
+#include "fvcFlux.H"
+#include "fvmDiv.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -130,6 +132,28 @@ bool Foam::momentumTransportModel::read()
 }
 
 
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::momentumTransportModel::divDevTauCorr
+(
+    const tmp<volTensorField>& devTauCorr,
+    volVectorField& U
+) const
+{
+    return fvm::divc
+    (
+        tmp<surfaceVectorField>
+        (
+            new surfaceVectorField
+            (
+                groupName("devTauCorrFlux"),
+                fvc::flux(devTauCorr)
+            )
+        ),
+        U
+    );
+}
+
+
 void Foam::momentumTransportModel::validate()
 {}
 

src/MomentumTransportModels/momentumTransportModels/momentumTransportModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,7 @@ SourceFiles
 #include "IOdictionary.H"
 #include "primitiveFieldsFwd.H"
 #include "volFields.H"
-#include "surfaceFieldsFwd.H"
+#include "surfaceFields.H"
 #include "fvMatricesFwd.H"
 #include "geometricOneField.H"
 #include "viscosity.H"
@@ -139,6 +139,11 @@ public:
         //- Const access to the coefficients dictionary
         virtual const dictionary& coeffDict() const = 0;
 
+        inline word groupName(const word& name) const
+        {
+            return IOobject::groupName(name, alphaRhoPhi_.group());
+        }
+
         //- Helper function to return the name of the turbulence G field
         inline word GName() const
         {
@@ -205,6 +210,17 @@ public:
         //- Return the stress tensor [m^2/s^2]
         virtual tmp<volSymmTensorField> sigma() const = 0;
 
+        //- Return the explicit stress correction matrix
+        //  for the momentum equation.
+        //  Combined with the implicit part by divDevTau
+        //  Also creates and registers the devTauCorrFlux field
+        //  which is cached in the matrix if fluxRequired.
+        virtual tmp<fvVectorMatrix> divDevTauCorr
+        (
+            const tmp<volTensorField>& devTauCorr,
+            volVectorField& U
+        ) const;
+
         //- Validate the fields after construction
         //  Update derived fields as required
         virtual void validate();