View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002418OpenFOAMContributionpublic2017-01-02 08:242017-01-13 14:26
ReporterSvensen Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status closedResolutionsuspended 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Product Versiondev 
Summary0002418: Time-averaged Wall Shear Stress
DescriptionIn the research of blood flow, commonly it is required to know not only the Wall Shear Stress, but an Time-average Wall Shear Stress (TAWSS) over the certain period of time (1 second for example).

The proposed contribution is a functionObject which is based on wallShearStress object. It implements this functionality. For more details about TAWSS please look on equation (14) at http://utopia.duth.gr/soulis/Publications%20(Selected%20Items)/Relative%20Residence%20Time%20and%20Oscillatory%20Shear%20Index%20of%20Non-.pdf

This functionality is often requested on cfd-online (for example, https://www.cfd-online.com/Forums/cfx/97793-time-averaged-shear-stress.html), but still was not implemented in OpenFOAM.

If I remember it right, Merge Request on GitHub is not the primary way to contribute the code now, that's why I uploaded it here.

Any comments, remarks?
Tagscontribution, functionObject

Activities

Svensen

2017-01-02 08:24

reporter  

time-averaged-WSS.patch (31,721 bytes)    
From 708a1eb93df69a689bba061c6476c1c1f2794275 Mon Sep 17 00:00:00 2001
From: Svensen <ssindeev@yandex.ru>
Date: Sat, 3 Dec 2016 11:45:43 +0100
Subject: [PATCH 1/7] Initial commit. no errors in compilation

---
 src/functionObjects/field/Make/files               |   1 +
 .../timeAverageWallShearStress.C                   | 267 +++++++++++++++++++++
 .../timeAverageWallShearStress.H                   | 178 ++++++++++++++
 3 files changed, 446 insertions(+)
 create mode 100644 src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
 create mode 100644 src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H

diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index de5be18..7de6df7 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -53,6 +53,7 @@ MachNo/MachNo.C
 turbulenceFields/turbulenceFields.C
 yPlus/yPlus.C
 wallShearStress/wallShearStress.C
+timeAverageWallShearStress/timeAverageWallShearStress.C
 wallHeatFlux/wallHeatFlux.C
 
 writeCellCentres/writeCellCentres.C
diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
new file mode 100644
index 0000000..8ed8fc4
--- /dev/null
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
@@ -0,0 +1,267 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "timeAverageWallShearStress.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
+#include "wallPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(timeAverageWallShearStress, 0);
+    addToRunTimeSelectionTable(functionObject, timeAverageWallShearStress, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::functionObjects::timeAverageWallShearStress::writeFileHeader(const label i)
+{
+    // Add headers to output data
+    writeHeader(file(), "Wall shear stress");
+    writeCommented(file(), "Time");
+    writeTabbed(file(), "patch");
+    writeTabbed(file(), "min");
+    writeTabbed(file(), "max");
+    file() << endl;
+}
+
+
+void Foam::functionObjects::timeAverageWallShearStress::calcShearStress
+(
+    const volSymmTensorField& Reff,
+    volVectorField& shearStress
+)
+{
+    shearStress.dimensions().reset(Reff.dimensions());
+
+    forAllConstIter(labelHashSet, patchSet_, iter)
+    {
+        label patchi = iter.key();
+
+        vectorField& ssp = shearStress.boundaryFieldRef()[patchi];
+        const vectorField& Sfp = mesh_.Sf().boundaryField()[patchi];
+        const scalarField& magSfp = mesh_.magSf().boundaryField()[patchi];
+        const symmTensorField& Reffp = Reff.boundaryField()[patchi];
+
+        ssp = (-Sfp/magSfp) & Reffp;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::timeAverageWallShearStress::timeAverageWallShearStress
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    logFiles(obr_, name),
+    patchSet_()
+{
+    volVectorField* wallShearStressPtr
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                type(),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector
+            (
+                "0",
+                sqr(dimLength)/sqr(dimTime),
+                Zero
+            )
+        )
+    );
+
+    mesh_.objectRegistry::store(wallShearStressPtr);
+
+    read(dict);
+    resetName(typeName);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::timeAverageWallShearStress::~timeAverageWallShearStress()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::timeAverageWallShearStress::read(const dictionary& dict)
+{
+    fvMeshFunctionObject::read(dict);
+
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+    patchSet_ =
+        mesh_.boundaryMesh().patchSet
+        (
+            wordReList(dict.lookupOrDefault("patches", wordReList()))
+        );
+
+    Info<< type() << " " << name() << ":" << nl;
+
+    if (patchSet_.empty())
+    {
+        forAll(pbm, patchi)
+        {
+            if (isA<wallPolyPatch>(pbm[patchi]))
+            {
+                patchSet_.insert(patchi);
+            }
+        }
+
+        Info<< "    processing all wall patches" << nl << endl;
+    }
+    else
+    {
+        Info<< "    processing wall patches: " << nl;
+        labelHashSet filteredPatchSet;
+        forAllConstIter(labelHashSet, patchSet_, iter)
+        {
+            label patchi = iter.key();
+            if (isA<wallPolyPatch>(pbm[patchi]))
+            {
+                filteredPatchSet.insert(patchi);
+                Info<< "        " << pbm[patchi].name() << endl;
+            }
+            else
+            {
+                WarningInFunction
+                    << "Requested wall shear stress on non-wall boundary "
+                    << "type patch: " << pbm[patchi].name() << endl;
+            }
+        }
+
+        Info<< endl;
+
+        patchSet_ = filteredPatchSet;
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::timeAverageWallShearStress::execute()
+{
+    typedef compressible::turbulenceModel cmpModel;
+    typedef incompressible::turbulenceModel icoModel;
+
+    volVectorField& wallShearStress =
+        const_cast<volVectorField&>
+        (
+            mesh_.lookupObject<volVectorField>(type())
+        );
+
+    tmp<volSymmTensorField> Reff;
+    if (mesh_.foundObject<cmpModel>(turbulenceModel::propertiesName))
+    {
+        const cmpModel& model =
+            mesh_.lookupObject<cmpModel>(turbulenceModel::propertiesName);
+
+        Reff = model.devRhoReff();
+    }
+    else if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
+    {
+        const icoModel& model =
+            mesh_.lookupObject<icoModel>(turbulenceModel::propertiesName);
+
+        Reff = model.devReff();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to find turbulence model in the "
+            << "database" << exit(FatalError);
+    }
+
+    calcShearStress(Reff(), wallShearStress);
+
+    return true;
+}
+
+
+bool Foam::functionObjects::timeAverageWallShearStress::write()
+{
+    logFiles::write();
+
+    const volVectorField& wallShearStress =
+        obr_.lookupObject<volVectorField>(type());
+
+    Log << type() << " " << name() << " write:" << nl
+        << "    writing field " << wallShearStress.name() << endl;
+
+    wallShearStress.write();
+
+    const fvPatchList& patches = mesh_.boundary();
+
+    forAllConstIter(labelHashSet, patchSet_, iter)
+    {
+        label patchi = iter.key();
+        const fvPatch& pp = patches[patchi];
+
+        const vectorField& ssp = wallShearStress.boundaryField()[patchi];
+
+        vector minSsp = gMin(ssp);
+        vector maxSsp = gMax(ssp);
+
+        if (Pstream::master())
+        {
+            file() << mesh_.time().value()
+                << token::TAB << pp.name()
+                << token::TAB << minSsp
+                << token::TAB << maxSsp
+                << endl;
+        }
+
+        Log << "    min/max(" << pp.name() << ") = "
+            << minSsp << ", " << maxSsp << endl;
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
new file mode 100644
index 0000000..50c070e
--- /dev/null
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::timeAverageWallShearStress
+
+Group
+    grpForcesFunctionObjects
+
+Description
+    Calculates and write the time average wall shear stress (TAWSS) at wall patches as
+    the volScalarField field 'TAWSS'.
+
+        \f[
+            Stress = R \dot n
+        \f]
+
+    where
+    \vartable
+        R       | stress tensor
+        n       | patch normal vector (into the domain)
+    \endvartable
+
+    The shear-stress symmetric tensor field is retrieved from the turbulence
+    model.  All wall patches are included by default; to restrict the
+    calculation to certain patches, use the optional 'patches' entry.
+
+    Example of function object specification:
+    \verbatim
+    timeAverageWallShearStress1
+    {
+        type        timeAverageWallShearStress;
+        libs        ("libfieldFunctionObjects.so");
+        ...
+        patches     (".*Wall");
+    }
+    \endverbatim
+
+Usage
+    \table
+        Property | Description               | Required    | Default value
+        type     | type name: timeAverageWallShearStress | yes        |
+        patches  | list of patches to process | no         | all wall patches
+    \endtable
+
+See also
+    Foam::functionObject
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObjects::logFiles
+    Foam::functionObjects::pressureTools
+    Foam::functionObjects::timeControl
+
+SourceFiles
+    timeAverageWallShearStress.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_timeAverageWallShearStress_H
+#define functionObjects_timeAverageWallShearStress_H
+
+#include "fvMeshFunctionObject.H"
+#include "logFiles.H"
+#include "volFieldsFwd.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class timeAverageWallShearStress Declaration
+\*---------------------------------------------------------------------------*/
+
+class timeAverageWallShearStress
+:
+    public fvMeshFunctionObject,
+    public logFiles
+{
+
+protected:
+
+    // Protected data
+
+        //- Optional list of patches to process
+        labelHashSet patchSet_;
+
+
+    // Protected Member Functions
+
+        //- File header information
+        virtual void writeFileHeader(const label i);
+
+        //- Calculate the shear-stress
+        void calcShearStress
+        (
+            const volSymmTensorField& Reff,
+            volVectorField& shearStress
+        );
+
+
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        timeAverageWallShearStress(const timeAverageWallShearStress&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const timeAverageWallShearStress&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("timeAverageWallShearStress");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        timeAverageWallShearStress
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~timeAverageWallShearStress();
+
+
+    // Member Functions
+
+        //- Read the wallShearStress data
+        virtual bool read(const dictionary&);
+
+        //- Calculate the wall shear-stress
+        virtual bool execute();
+
+        //- Write the wall shear-stress
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
-- 
1.9.1


From b6acb04076a7733e3faa7d2bff2a51fbce77c212 Mon Sep 17 00:00:00 2001
From: Svensen <ssindeev@yandex.ru>
Date: Sat, 3 Dec 2016 13:58:17 +0100
Subject: [PATCH 2/7] TAWSS is computed and written only in [start; stop]
 period

---
 .../timeAverageWallShearStress.C                   | 122 ++++++++++-----------
 .../timeAverageWallShearStress.H                   |   6 +
 2 files changed, 67 insertions(+), 61 deletions(-)

diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
index 8ed8fc4..96134a3 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
@@ -90,7 +90,9 @@ Foam::functionObjects::timeAverageWallShearStress::timeAverageWallShearStress
 :
     fvMeshFunctionObject(name, runTime, dict),
     logFiles(obr_, name),
-    patchSet_()
+    patchSet_(),
+    startTime_(GREAT),
+    stopTime_(GREAT)
 {
     volVectorField* wallShearStressPtr
     (
@@ -179,6 +181,9 @@ bool Foam::functionObjects::timeAverageWallShearStress::read(const dictionary& d
 
         patchSet_ = filteredPatchSet;
     }
+    
+    dict.lookup("startTime") >> startTime_;
+    dict.lookup("stopTime") >> stopTime_;
 
     return true;
 }
@@ -186,78 +191,73 @@ bool Foam::functionObjects::timeAverageWallShearStress::read(const dictionary& d
 
 bool Foam::functionObjects::timeAverageWallShearStress::execute()
 {
-    typedef compressible::turbulenceModel cmpModel;
-    typedef incompressible::turbulenceModel icoModel;
-
-    volVectorField& wallShearStress =
-        const_cast<volVectorField&>
-        (
-            mesh_.lookupObject<volVectorField>(type())
-        );
-
-    tmp<volSymmTensorField> Reff;
-    if (mesh_.foundObject<cmpModel>(turbulenceModel::propertiesName))
-    {
-        const cmpModel& model =
-            mesh_.lookupObject<cmpModel>(turbulenceModel::propertiesName);
-
-        Reff = model.devRhoReff();
-    }
-    else if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
-    {
-        const icoModel& model =
-            mesh_.lookupObject<icoModel>(turbulenceModel::propertiesName);
-
-        Reff = model.devReff();
+	const scalar currentTime = mesh_.time().value();
+	Info << currentTime << endl;
+	
+	if (currentTime >= startTime_ && currentTime <= stopTime_)
+	{
+		Info << "TAWSS computing" << endl;
+		typedef compressible::turbulenceModel cmpModel;
+		typedef incompressible::turbulenceModel icoModel;
+
+		volVectorField& wallShearStress =
+		    const_cast<volVectorField&>
+		    (
+		        mesh_.lookupObject<volVectorField>(type())
+		    );
+
+		tmp<volSymmTensorField> Reff;
+		if (mesh_.foundObject<cmpModel>(turbulenceModel::propertiesName))
+		{
+		    const cmpModel& model =
+		        mesh_.lookupObject<cmpModel>(turbulenceModel::propertiesName);
+
+		    Reff = model.devRhoReff();
+		}
+		else if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
+		{
+		    const icoModel& model =
+		        mesh_.lookupObject<icoModel>(turbulenceModel::propertiesName);
+
+		    Reff = model.devReff();
+		}
+		else
+		{
+		    FatalErrorInFunction
+		        << "Unable to find turbulence model in the "
+		        << "database" << exit(FatalError);
+		}
+
+		calcShearStress(Reff(), wallShearStress);
     }
     else
     {
-        FatalErrorInFunction
-            << "Unable to find turbulence model in the "
-            << "database" << exit(FatalError);
+    	Info << "TAWSS out of bounds" << endl;
     }
 
-    calcShearStress(Reff(), wallShearStress);
-
     return true;
 }
 
 
 bool Foam::functionObjects::timeAverageWallShearStress::write()
 {
-    logFiles::write();
-
-    const volVectorField& wallShearStress =
-        obr_.lookupObject<volVectorField>(type());
-
-    Log << type() << " " << name() << " write:" << nl
-        << "    writing field " << wallShearStress.name() << endl;
-
-    wallShearStress.write();
-
-    const fvPatchList& patches = mesh_.boundary();
-
-    forAllConstIter(labelHashSet, patchSet_, iter)
+	const volVectorField& wallShearStress =
+    obr_.lookupObject<volVectorField>(type());
+		    
+	const scalar currentTime = mesh_.time().value();
+	
+	if (currentTime == stopTime_)
+	{
+		logFiles::write();
+
+		Log << type() << " " << name() << " write:" << nl
+		    << "    writing field " << wallShearStress.name() << endl;
+
+		wallShearStress.write();  
+    }
+    else
     {
-        label patchi = iter.key();
-        const fvPatch& pp = patches[patchi];
-
-        const vectorField& ssp = wallShearStress.boundaryField()[patchi];
-
-        vector minSsp = gMin(ssp);
-        vector maxSsp = gMax(ssp);
-
-        if (Pstream::master())
-        {
-            file() << mesh_.time().value()
-                << token::TAB << pp.name()
-                << token::TAB << minSsp
-                << token::TAB << maxSsp
-                << endl;
-        }
-
-        Log << "    min/max(" << pp.name() << ") = "
-            << minSsp << ", " << maxSsp << endl;
+    	Info << "TAWSS is not written" << endl;
     }
 
     return true;
diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
index 50c070e..1ddaa18 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
@@ -106,6 +106,12 @@ protected:
 
         //- Optional list of patches to process
         labelHashSet patchSet_;
+        
+        //- Start time of WSS collecting
+        scalar startTime_;
+        
+        //- Stop time of WSS collecting
+        scalar stopTime_;
 
 
     // Protected Member Functions
-- 
1.9.1


From d266c0c3ab0be7f896b15b919185c752fadc7554 Mon Sep 17 00:00:00 2001
From: Svensen <ssindeev@yandex.ru>
Date: Sun, 4 Dec 2016 01:22:14 +0100
Subject: [PATCH 3/7] bugfix for previous commit

---
 .../timeAverageWallShearStress.C                   | 49 ++++++++++++++--------
 .../timeAverageWallShearStress.H                   |  5 +++
 2 files changed, 36 insertions(+), 18 deletions(-)

diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
index 96134a3..352eaa7 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
@@ -92,7 +92,25 @@ Foam::functionObjects::timeAverageWallShearStress::timeAverageWallShearStress
     logFiles(obr_, name),
     patchSet_(),
     startTime_(GREAT),
-    stopTime_(GREAT)
+    stopTime_(GREAT),
+    TAWSS_(new volScalarField
+        (
+            IOobject
+            (
+                "TAWSS",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "TAWSS-init",
+                sqr(dimLength)/sqr(dimTime),
+                0
+            )
+        ))
 {
     volVectorField* wallShearStressPtr
     (
@@ -191,10 +209,10 @@ bool Foam::functionObjects::timeAverageWallShearStress::read(const dictionary& d
 
 bool Foam::functionObjects::timeAverageWallShearStress::execute()
 {
-	const scalar currentTime = mesh_.time().value();
-	Info << currentTime << endl;
+	const scalar currentTime_ = mesh_.time().value();
+	Info << currentTime_ << endl;
 	
-	if (currentTime >= startTime_ && currentTime <= stopTime_)
+	if (currentTime_ >= (startTime_-SMALL) && currentTime_ <= (stopTime_+SMALL))
 	{
 		Info << "TAWSS computing" << endl;
 		typedef compressible::turbulenceModel cmpModel;
@@ -229,10 +247,12 @@ bool Foam::functionObjects::timeAverageWallShearStress::execute()
 		}
 
 		calcShearStress(Reff(), wallShearStress);
+		*TAWSS_ += mag(wallShearStress);
     }
     else
     {
-    	Info << "TAWSS out of bounds" << endl;
+    	Info << "TAWSS at " << currentTime_ << " out of bounds" << endl;
+    	Info << "startTime_ = " << startTime_ << " stopTime_ = " << stopTime_ << endl;
     }
 
     return true;
@@ -240,20 +260,13 @@ bool Foam::functionObjects::timeAverageWallShearStress::execute()
 
 
 bool Foam::functionObjects::timeAverageWallShearStress::write()
-{
-	const volVectorField& wallShearStress =
-    obr_.lookupObject<volVectorField>(type());
-		    
-	const scalar currentTime = mesh_.time().value();
-	
-	if (currentTime == stopTime_)
-	{
-		logFiles::write();
+{  
+	const scalar currentTime = obr_.time().value();
 
-		Log << type() << " " << name() << " write:" << nl
-		    << "    writing field " << wallShearStress.name() << endl;
-
-		wallShearStress.write();  
+	if (mag(currentTime-stopTime_) < SMALL)
+	{
+		Info << "TAWSS is written !!!" << endl;
+		TAWSS_->write();
     }
     else
     {
diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
index 1ddaa18..d9390b6 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
@@ -80,6 +80,8 @@ SourceFiles
 
 #include "fvMeshFunctionObject.H"
 #include "logFiles.H"
+#include "volFields.H"
+#include "surfaceFields.H"
 #include "volFieldsFwd.H"
 #include "HashSet.H"
 
@@ -112,6 +114,9 @@ protected:
         
         //- Stop time of WSS collecting
         scalar stopTime_;
+        
+        //- Time average Wall Shear Stress
+        volScalarField* TAWSS_;
 
 
     // Protected Member Functions
-- 
1.9.1


From 0266b5e188aa7dfad1c20cd283aad19021bb8afb Mon Sep 17 00:00:00 2001
From: Svensen <ssindeev@yandex.ru>
Date: Sun, 4 Dec 2016 11:26:33 +0100
Subject: [PATCH 4/7] TAWSS computation according to integral formula

---
 .../timeAverageWallShearStress.C                   | 32 +++++++++++-----------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
index 352eaa7..de4897b 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
@@ -48,7 +48,7 @@ namespace functionObjects
 void Foam::functionObjects::timeAverageWallShearStress::writeFileHeader(const label i)
 {
     // Add headers to output data
-    writeHeader(file(), "Wall shear stress");
+    writeHeader(file(), "Time average wall shear stress (TAWSS)");
     writeCommented(file(), "Time");
     writeTabbed(file(), "patch");
     writeTabbed(file(), "min");
@@ -190,7 +190,7 @@ bool Foam::functionObjects::timeAverageWallShearStress::read(const dictionary& d
             else
             {
                 WarningInFunction
-                    << "Requested wall shear stress on non-wall boundary "
+                    << "Requested time average wall shear stress on non-wall boundary "
                     << "type patch: " << pbm[patchi].name() << endl;
             }
         }
@@ -214,7 +214,11 @@ bool Foam::functionObjects::timeAverageWallShearStress::execute()
 	
 	if (currentTime_ >= (startTime_-SMALL) && currentTime_ <= (stopTime_+SMALL))
 	{
-		Info << "TAWSS computing" << endl;
+		if (debug)
+		{
+			Info << "TAWSS computing" << endl;
+		}
+		
 		typedef compressible::turbulenceModel cmpModel;
 		typedef incompressible::turbulenceModel icoModel;
 
@@ -247,14 +251,11 @@ bool Foam::functionObjects::timeAverageWallShearStress::execute()
 		}
 
 		calcShearStress(Reff(), wallShearStress);
-		*TAWSS_ += mag(wallShearStress);
-    }
-    else
-    {
-    	Info << "TAWSS at " << currentTime_ << " out of bounds" << endl;
-    	Info << "startTime_ = " << startTime_ << " stopTime_ = " << stopTime_ << endl;
-    }
 
+		// TAWSS = TAWSS + WSS*dt;
+		*TAWSS_ += (mag(wallShearStress)*obr_.time().deltaTValue());
+    }
+    
     return true;
 }
 
@@ -265,14 +266,13 @@ bool Foam::functionObjects::timeAverageWallShearStress::write()
 
 	if (mag(currentTime-stopTime_) < SMALL)
 	{
-		Info << "TAWSS is written !!!" << endl;
+		scalar T = stopTime_ - startTime_;
+		*TAWSS_ /= T;
+		
 		TAWSS_->write();
+		Info << "TAWSS is written" << endl;
     }
-    else
-    {
-    	Info << "TAWSS is not written" << endl;
-    }
-
+    
     return true;
 }
 
-- 
1.9.1


From 9f5de26788d4852c68f5f38bb84388f147a63737 Mon Sep 17 00:00:00 2001
From: Svensen <ssindeev@yandex.ru>
Date: Sun, 4 Dec 2016 13:19:54 +0100
Subject: [PATCH 5/7] TAWSS computation works right

---
 .../field/timeAverageWallShearStress/timeAverageWallShearStress.C  | 7 +++++--
 .../field/timeAverageWallShearStress/timeAverageWallShearStress.H  | 2 ++
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
index de4897b..d1de447 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
@@ -253,7 +253,7 @@ bool Foam::functionObjects::timeAverageWallShearStress::execute()
 		calcShearStress(Reff(), wallShearStress);
 
 		// TAWSS = TAWSS + WSS*dt;
-		*TAWSS_ += (mag(wallShearStress)*obr_.time().deltaTValue());
+		*TAWSS_ += mag(wallShearStress);//*obr_.time().deltaTValue());
     }
     
     return true;
@@ -266,7 +266,10 @@ bool Foam::functionObjects::timeAverageWallShearStress::write()
 
 	if (mag(currentTime-stopTime_) < SMALL)
 	{
-		scalar T = stopTime_ - startTime_;
+		scalar dt = obr_.time().deltaTValue();
+		scalar T = (stopTime_ - startTime_) + dt;
+		T /= dt;
+		
 		*TAWSS_ /= T;
 		
 		TAWSS_->write();
diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
index d9390b6..10f554e 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.H
@@ -61,6 +61,8 @@ Usage
         Property | Description               | Required    | Default value
         type     | type name: timeAverageWallShearStress | yes        |
         patches  | list of patches to process | no         | all wall patches
+        startTime | start computing TAWSS from startTime | no		|	yes
+        stopTime | stop computing TAWSS	at stopTime	| no		|	yes
     \endtable
 
 See also
-- 
1.9.1


From 72d85d3f5db1d4ee315da06439a2c0f9b21f6971 Mon Sep 17 00:00:00 2001
From: Svensen <ssindeev@yandex.ru>
Date: Sun, 4 Dec 2016 13:47:12 +0100
Subject: [PATCH 6/7] startTime and stopTime have the default falues

---
 .../field/timeAverageWallShearStress/timeAverageWallShearStress.C     | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
index d1de447..98da5c7 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
@@ -200,8 +200,8 @@ bool Foam::functionObjects::timeAverageWallShearStress::read(const dictionary& d
         patchSet_ = filteredPatchSet;
     }
     
-    dict.lookup("startTime") >> startTime_;
-    dict.lookup("stopTime") >> stopTime_;
+	startTime_ = dict.lookupOrDefault("startTime", obr_.time().startTime().value());
+    stopTime_ = dict.lookupOrDefault("stopTime", obr_.time().endTime().value());
 
     return true;
 }
-- 
1.9.1


From d1d8de48224f5cda9fc1ff4a16bc7460fb782561 Mon Sep 17 00:00:00 2001
From: Svensen <ssindeev@yandex.ru>
Date: Mon, 2 Jan 2017 09:07:37 +0100
Subject: [PATCH 7/7] works OK

---
 .../field/timeAverageWallShearStress/timeAverageWallShearStress.C       | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
index 98da5c7..1f1e552 100644
--- a/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
+++ b/src/functionObjects/field/timeAverageWallShearStress/timeAverageWallShearStress.C
@@ -253,7 +253,7 @@ bool Foam::functionObjects::timeAverageWallShearStress::execute()
 		calcShearStress(Reff(), wallShearStress);
 
 		// TAWSS = TAWSS + WSS*dt;
-		*TAWSS_ += mag(wallShearStress);//*obr_.time().deltaTValue());
+		*TAWSS_ += mag(wallShearStress);
     }
     
     return true;
-- 
1.9.1
time-averaged-WSS.patch (31,721 bytes)   

henry

2017-01-05 12:27

manager   ~0007594

Time-averaging fields is generally useful and not limited to 'wallShearStress' so it would be better if the averaging methods are implemented in a more general way so they are applicable to as wide a range of applications as possible.

The 'fieldAverage' 'functionObject' already provided with OpenFOAM is fairly general and includes support for averaging windows and appears to fit your requirements. Have you already tried this averaging method and found it lacking is some manner? If so it would be better to extend this 'functionObject' to cover your needs rather than write a completely specialized one only applicable to averaging 'wallShearStress'.

Svensen

2017-01-06 05:53

reporter   ~0007601

'fieldAverage' will sum a 'wallShearStress' VECTOR and divide it by a number of time steps.

In my case I need to sum the 'wallShearStress' MAGNITUDE and divide it by a number of time steps. As I see from the documentation for 'fieldAverage' it is not possible to operate with a vector MAGNITUDE.

The second point is that my 'Time-averaged Wall Shear Stress' object is only the beginning. I've just wanted to get some reaction on this idea. My plan is to implement all other formulas (15)-(17) from the attached document. And I think that it is better to do with the separate functionObjects for each equation.
However, you, as architect, may suggest another solution to this case. It maybe a single functionObject with some flags like "TAWSS yes; TAWSSV no; OSI no; RRT no; ... " and so on.

However, as I said before, these parameters are useful for analysis of blood flow and are of particular interest for the community in this field.

Any suggestions ?

henry

2017-01-06 07:39

manager   ~0007602

You could use the 'mag' 'functionObject' to take the magnitude of the 'wallShearStress' before averaging or it might be generally useful to add a 'mag' option to the 'wallShearStress' 'functionObject'.

As architect and maintainer of OpenFOAM I look for the most general implementations of functionality which introduce least code and code duplication to increase the usefulness of OpenFOAM to the largest number of users with the least increase in maintainence cost and overhead.

henry

2017-01-13 14:26

manager   ~0007633

Waiting for feedback.

Issue History

Date Modified Username Field Change
2017-01-02 08:24 Svensen New Issue
2017-01-02 08:24 Svensen File Added: time-averaged-WSS.patch
2017-01-02 08:24 Svensen Tag Attached: contribution
2017-01-02 08:24 Svensen Tag Attached: functionObject
2017-01-05 12:27 henry Note Added: 0007594
2017-01-06 05:53 Svensen Note Added: 0007601
2017-01-06 07:39 henry Note Added: 0007602
2017-01-13 14:26 henry Assigned To => henry
2017-01-13 14:26 henry Status new => closed
2017-01-13 14:26 henry Resolution open => suspended
2017-01-13 14:26 henry Note Added: 0007633