ENH: support time-dependent actuationDiskSource disk direction (#3099)

CONTRIBUTORS.md

@@ -9,6 +9,7 @@ It is likely incomplete...
 - Yu Ankun
 - Tetsuo Aoyagi
 - Akira Azami
+- Pete Bachant
 - William Bainbridge
 - Gabriel Barajas
 - Kutalmis Bercin

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C

@@ -90,6 +90,8 @@ void Foam::fv::actuationDiskSource::writeFileHeader(Ostream& os)
         writeFile::writeCommented(os, "P");
     }
 
+    writeFile::writeCommented(os, "diskDir");
+
     os  << endl;
 }
 
@@ -208,14 +210,7 @@ Foam::fv::actuationDiskSource::actuationDiskSource
             scalarMinMax::ge(VSMALL)
         )
     ),
-    diskDir_
+    diskDir_(Function1<vector>::New("diskDir", coeffs_, &mesh)),
-    (
-        coeffs_.getCheck<vector>
-        (
-            "diskDir",
-            [&](const vector& vec){ return mag(vec) > VSMALL; }
-        ).normalise()
-    ),
     UvsCpPtr_(Function1<scalar>::New("Cp", coeffs_, &mesh)),
     UvsCtPtr_(Function1<scalar>::New("Ct", coeffs_, &mesh)),
     monitorCells_()
@@ -237,6 +232,25 @@ Foam::fv::actuationDiskSource::actuationDiskSource
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::vector Foam::fv::actuationDiskSource::diskDir() const
+{
+    const scalar t = mesh_.time().timeOutputValue();
+    const vector dir(diskDir_->value(t));
+    const scalar magDir = Foam::mag(dir);
+
+    if (magDir < SMALL)
+    {
+        FatalErrorInFunction
+            << "Actuator disk surface-normal vector is zero: " << dir
+            << " at time=" << t
+            << exit(FatalError);
+    }
+
+    // normalised:
+    return dir/magDir;
+}
+
+
 void Foam::fv::actuationDiskSource::addSup
 (
     fvMatrix<vector>& eqn,
@@ -295,19 +309,10 @@ bool Foam::fv::actuationDiskSource::read(const dictionary& dict)
                 << exit(FatalIOError);
         }
 
-        dict.readIfPresent("diskDir", diskDir_);
+        // TBD: runTime re-reading of "diskDir" ?
-        diskDir_.normalise();
-        if (mag(diskDir_) < VSMALL)
-        {
-            FatalIOErrorInFunction(dict)
-                << "Actuator disk surface-normal vector is zero: "
-                << "diskDir = " << diskDir_
-                << exit(FatalIOError);
-        }
 
         return true;
     }
-
     return false;
 }
 

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H

@@ -170,7 +170,7 @@ Usage
       type       | Type name: actuationDiskSource      | word | yes      | -
       diskArea   | Actuator disk planar surface area   | scalar | yes    | -
       diskDir    | Surface-normal vector of the actuator disk <!--
-               --> pointing downstream                 | vector    | yes | -
+               --> pointing downstream                 | Function1 | yes | -
       Cp         | Power coefficient                   | Function1 | yes | -
       Ct         | Thrust coefficient                  | Function1 | yes | -
       monitorMethod | Type of incoming velocity monitoring method   <!--
@@ -207,8 +207,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef actuationDiskSource_H
+#ifndef Foam_actuationDiskSource_H
-#define actuationDiskSource_H
+#define Foam_actuationDiskSource_H
 
 #include "cellSetOption.H"
 #include "writeFile.H"
@@ -276,7 +276,8 @@ protected:
         scalar diskArea_;
 
         //- Surface-normal vector of the actuator disk pointing downstream
-        vector diskDir_;
+        //- as a function of time
+        autoPtr<Function1<vector>> diskDir_;
 
         //- Velocity vs power coefficients
         autoPtr<Function1<scalar>> UvsCpPtr_;
@@ -364,14 +365,11 @@ public:
 
         // Access
 
-            //- Normal disk direction
+            //- Normal disk direction, evaluated at timeOutputValue
-            const vector& diskDir() const
+            vector diskDir() const;
-            {
-                return diskDir_;
-            }
 
             //- Disk area
-            scalar diskArea() const
+            scalar diskArea() const noexcept
             {
                 return diskArea_;
             }

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSourceTemplates.C

@@ -110,12 +110,13 @@ void Foam::fv::actuationDiskSource::calcFroudeMethod
     }
 
     // (BJSB:Eq. 3.9)
+    const vector diskDir = this->diskDir();
     const scalar a = 1.0 - Cp/Ct;
-    const scalar T = 2.0*rhoRef*diskArea_*magSqr(Uref & diskDir_)*a*(1 - a);
+    const scalar T = 2.0*rhoRef*diskArea_*magSqr(Uref & diskDir)*a*(1 - a);
 
     for (const label celli : cells_)
     {
-        Usource[celli] += ((cellsV[celli]/V())*T)*diskDir_;
+        Usource[celli] += ((cellsV[celli]/V())*T)*diskDir;
     }
 
     if
@@ -128,7 +129,7 @@ void Foam::fv::actuationDiskSource::calcFroudeMethod
         writeCurrentTime(os);
 
         os  << Uref << tab << Cp << tab << Ct << tab << a << tab << T
-            << endl;
+            << tab << diskDir << endl;
     }
 }
 
@@ -221,12 +222,13 @@ void Foam::fv::actuationDiskSource::calcVariableScalingMethod
     const scalar CpStar = Cp*pow3(magUref/magUdisk);
 
     // Compute calibrated thrust/power (LSRMTK:Eq. 5)
-    const scalar T = 0.5*rhoRef*diskArea_*magSqr(Udisk & diskDir_)*CtStar;
+    const vector diskDir = this->diskDir();
-    const scalar P = 0.5*rhoRef*diskArea_*pow3(mag(Udisk & diskDir_))*CpStar;
+    const scalar T = 0.5*rhoRef*diskArea_*magSqr(Udisk & diskDir)*CtStar;
+    const scalar P = 0.5*rhoRef*diskArea_*pow3(mag(Udisk & diskDir))*CpStar;
 
     for (const label celli : cells_)
     {
-        Usource[celli] += (cellsV[celli]/totalV*T)*diskDir_;
+        Usource[celli] += (cellsV[celli]/totalV*T)*diskDir;
     }
 
     if
@@ -240,7 +242,7 @@ void Foam::fv::actuationDiskSource::calcVariableScalingMethod
 
         os  << Uref << tab << Cp << tab << Ct << tab
             << Udisk << tab << CpStar << tab << CtStar << tab << T << tab << P
-            << endl;
+            << tab << diskDir << endl;
     }
 }
 

src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSourceTemplates.C

@@ -48,9 +48,8 @@ addRadialActuationDiskAxialInertialResistance
     scalarField Tr(cells.size());
 
     tensor E(Zero);
-    E.xx() = diskDir_.x();
+    const vector diskDir = this->diskDir();
-    E.yy() = diskDir_.y();
+    E.diag(diskDir);
-    E.zz() = diskDir_.z();
 
     const Field<vector> zoneCellCentres(mesh().cellCentres(), cells);
     const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);

tutorials/incompressible/simpleFoam/turbineSiting/constant/fvOptions

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2312                                 |
+|  \\    /   O peration     | Version:  v2406                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -21,7 +21,7 @@ disk1
     selectionMode   cellSet;
     cellSet         actuationDisk1;
     diskArea        40;
-    diskDir         (1 0 0);
+    diskDir         constant (1 0 0);  // Or other Function1
     writeToFile     true;
     sink            true;
     Cp              0.386;
@@ -45,7 +45,7 @@ disk2
     Cp              0.53;
     Ct              0.58;
     diskArea        40;
-    diskDir         (1 0 0);
+    diskDir         (1 0 0);  // Or other Function1
     upstreamPoint   (581753 4785663 1070);
 }