stabilityBlendingFactor.H 14.5 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
8
    Copyright (C) 2018-2020 OpenCFD Ltd.
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
-------------------------------------------------------------------------------
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::stabilityBlendingFactor

Group
    grpFieldFunctionObjects

Description
33
    Computes the \c stabilityBlendingFactor to be used by the
34
    local blended convection scheme. The output is a surface field weight
35
    between 0-1.
36

37 38
    The weight of a blended scheme, i.e. \c w, is given by a function of
    the blending factor, \c f:
39 40

    \f[
41
        w = f_{scheme_1} + (1 - f_{scheme_2})
42 43
    \f]

44
    The factor is calculated based on six criteria:
45 46 47 48 49 50 51 52
    \verbatim
      1. mesh non-orthogonality field
      2. magnitude of cell centres gradient
      3. convergence rate of residuals
      4. faceWeight
      5. skewness
      6. Courant number
    \endverbatim
53 54 55

    The user can enable them individually.

56 57
    For option 1, the following relation is used, where \f$\phi_1\f$ is
    the non-orthogonality:
58
    \f[
59 60
        fNon =
            min
61 62 63 64
            (
                max
                (
                    0.0,
65 66
                    (\phi_1 - max(\phi_1))
                    /(min(\phi_1) - max(\phi_1))
67 68 69
                ),
                1.0
            )
70
    \f]
71

72 73 74
    For option 2, the following relation is used, where \f$\phi_2\f$ is
    the magnitude of cell centres gradient (Note that \f$\phi_2 = 3\f$
    for orthogonal meshes):
75

76 77
    \f[
        fMagGradCc =
78 79 80 81 82
            min
            (
                max
                (
                    0.0,
83 84
                    (\phi_2 - max(\phi_2))
                    / (min(\phi_2) - max(\phi_2))
85 86 87
                ),
                1.0
            )
88 89 90 91
    \f]

    For option 3, a PID control is used in order to control residual
    unbounded fluctuations for individual cells.
92

93 94 95 96 97 98
    \f[
        factor =
            P*residual
            + I*residualIntegral
            + D*residualDifferential
    \f]
99

100
    where \c P, \c I and \c D are user inputs.
101

102
    The following relation is used:
103
    \f[
104
        fRes = (factor - meanRes)/(maxRes*meanRes);
105
    \f]
106

107
    where
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
    \vartable
        meanRes | Average(residual)
        maxRes  | User input
    \endvartable

    Note that \f$f_{Res}\f$ will blend more towards one as
    the cell residual is larger then the domain mean residuals.


    For option 4, the following relation is used, where \f$\phi_4\f$ is
    the face weight (Note that \f$\phi_4 = 0.5\f$ for orthogonal meshes):

    \f[
        ffaceWeight = min
        (
            max
            (
                0.0,
                (min(\phi_4) - \phi_4)
                / (min(\phi_4) - max(\phi_4))
            ),
            1.0
        )
    \f]


    For option 5, the following relation is used, where \f$\phi_5\f$ is
    the cell skewness:

    \f[
        fskewness =
        min
        (
            max
            (
                0.0,
                (\phi_5    - max(\phi_5))
                / (min(\phi_5) - max(\phi_5))
            ),
            1.0
        )
    \f]


    For option 6, the following relation is used:
153

154 155 156 157 158
    \f[
        fCoWeight = min(max((Co - Co1)/(Co2 - Co1), 0), 1)
    \f]

    where
159 160 161 162 163 164 165 166 167 168 169
    \vartable
        Co1 | Courant number below which scheme2 is used
        Co2 | Courant number above which scheme1 is used
    \endvartable

    The final factor is determined by:

    \f[
        f = max(fNon, fMagGradCc, fRes, ffaceWeight, fskewness, fCoWeight)
    \f]

170 171
    An indicator (volume) field, named \c blendedIndicator
    is generated if the log flag is on:
172 173 174 175 176 177
    - 1 represent scheme1 as active,
    - 0 represent scheme2 as active.

    Additional reporting is written to the standard output, providing
    statistics as to the number of cells used by each scheme.

178 179 180 181 182 183 184 185
    Operands:
    \table
      Operand        | Type           | Location
      input          | -              | -
      output file    | dat | $FOAM_CASE/postProcessing/\<FO\>/\<time\>/\<file\>
      output field   | volScalarField | $FOAM_CASE/\<time\>/\<outField\>
    \endtable

186
Usage
187
    Minimal example by using \c system/controlDict.functions:
188 189 190
    \verbatim
    stabilityBlendingFactor1
    {
191
        // Mandatory entries (unmodifiable)
192
        type                stabilityBlendingFactor;
193
        libs                (fieldFunctionObjects);
194

195 196 197 198 199 200
        // Mandatory entries (unmodifiable)
        field               <field>;    // U;
        result              <outField>; // UBlendingFactor;

        // Optional entries (runtime modifiable)
        tolerance           0.001;
201

202
        // Any of the options can be chosen in combinations
203

204 205 206
        // Option-1
        switchNonOrtho      true;
        nonOrthogonality    nonOrthoAngle;
207 208 209
        maxNonOrthogonality 20;
        minNonOrthogonality 60;

210 211
        // Option-2
        switchGradCc        true;
212 213 214
        maxGradCc           3;
        minGradCc           4;

215 216
        // Option-3
        switchResiduals     true;
217 218 219 220 221 222
        maxResidual         10;
        residual            initialResidual:p;
        P                   1.5;
        I                   0;
        D                   0.5;

223 224 225 226 227 228 229 230 231 232 233 234 235
        // Option-4
        switchFaceWeight    true;
        maxFaceWeight       0.3;
        minFaceWeight       0.2;

        // Option-5
        switchSkewness      true;
        maxSkewness         2;
        minSkewness         3;

        // Option-6
        switchCo            true;
        U                   U;
236 237
        Co1                 1;
        Co2                 10;
238 239

        // Optional (inherited) entries
240 241
        ...
    }
242
    \endverbatim
243

244 245
    Example of function object specification to calculate the \c residuals used
    by \c stabilityBlendingFactor. The following writes 'initialResidual:p'
246
    field
247 248 249 250 251 252 253 254 255
    \verbatim
    residuals
    {
        type            residuals;
        libs            (utilityFunctionObjects);
        writeFields     true;
        writeControl    writeTime;
        fields          (p);
    }
256 257
    \endverbatim

258
    where the entries mean:
259
    \table
260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
      Property     | Description                           | Type | Req'd | Dflt
      type         | Type name: stabilityBlendingFactor    | word |  yes  | -
      libs         | Library name: fieldFunctionObjects    | word |  yes  | -
      field        | Name of operand field                 | word |  yes  | -
      result  | Name of surface field to be used in the localBlended scheme <!--
          --> | word | yes
      switchNonOrtho | Select non-orthogonal method        | bool | no | false
      nonOrthogonality | Name of the non-orthogonal field  <!--
                   --> | word | no | nonOrthoAngle
      maxNonOrthogonality| Maximum non-orthogonal for scheme2 | scalar | no | 20
      minNonOrthogonality| Minimum non-orthogonal for scheme1 | scalar | no | 60
      switchGradCc | Select cell centre gradient method    | bool | no | false
      maxGradCc| Maximum gradient for scheme2              | scalar | no | 2
      minGradCc| Minimum gradient for scheme1              | scalar | no | 4
      switchResiduals | Select residual evolution method   | bool | no | false
      residual    | Name of the residual field | word | no | initialResidual:p
      maxResidual| Maximum residual-mean ratio for scheme1 | scalar | no | 10
      P       | Proportional factor for PID                | scalar | no | 3
      I       | Integral factor for PID                    | scalar | no | 0
      D       | Differential factor for PID                | scalar | no | 0.25
      switchFaceWeight | Select face weight method         | bool | no | false
      faceWeight | Name of the faceWeight field       | word | no | faceWeight
      maxFaceWeight | Maximum face weight for scheme1      | scalar | no | 0.2
      minFaceWeight | Minimum face weight for scheme2      | scalar | no | 0.3
      switchSkewness   | Select skewness method            | bool | no | false
      skewness | Name of the skewness field           | word | no | skewness
      maxSkewness | Maximum skewness for scheme2           | scalar | no | 2
      minSkewness | Minimum skewness for scheme1           | scalar | no | 3
      switchCo         | Select Co blended method          | bool | no | false
      U   | Name of the flux field for Co blended          | word | no | U
      Co1 | Courant number below which scheme2 is used     | scalar | no | 1
      Co2 | Courant number above which scheme1 is used     | scalar | no | 10
      tolerance    | Tolerance for number of blended cells | scalar | no | 0.001
293 294
    \endtable

295 296 297
    The \c result entry is the field which is read by the \c localBlended scheme
    specified in \c fvSchemes. This name is determined by the \c localBlended
    class.
298

299 300 301 302
    The inherited entries are elaborated in:
     - \link functionObject.H \endlink
     - \link fieldExpression.H \endlink
     - \link writeFile.H \endlink
303

304
    Usage by the \c postProcess utility is not available.
305

306
See also
307 308 309 310 311
    - Foam::functionObject
    - Foam::functionObjects::fieldExpression
    - Foam::functionObjects::fvMeshFunctionObject
    - Foam::functionObjects::writeFile
    - ExtendedCodeGuide::functionObjects::field::stabilityBlendingFactor
312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340

SourceFiles
    stabilityBlendingFactor.C

\*---------------------------------------------------------------------------*/

#ifndef functionObjects_stabilityBlendingFactor_H
#define functionObjects_stabilityBlendingFactor_H

#include "fieldExpression.H"
#include "writeFile.H"
#include "volFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace functionObjects
{

/*---------------------------------------------------------------------------*\
                       Class stabilityBlendingFactor Declaration
\*---------------------------------------------------------------------------*/

class stabilityBlendingFactor
:
    public fieldExpression,
    public writeFile
{
341
    // Private Member Data
342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404

        //- Cell based blended indicator
        volScalarField indicator_;

        // Switches

            //- Switch for non-orthogonality
            Switch nonOrthogonality_;

            //- Switch for grad of cell centres
            Switch gradCc_;

            //- Switch for residuals
            Switch residuals_;

            //- Switch for face weight
            Switch faceWeight_;

            //- Switch for skewness
            Switch skewness_;

            //- Switch for Co
            Switch Co_;


        // Lower and upper limits

            //- Maximum non-orthogonality for fully scheme 2
            scalar maxNonOrthogonality_;

            //- Minimum non-orthogonality for fully scheme 1
            scalar minNonOrthogonality_;

            //- Maximum  gradcc for fully scheme 2
            scalar maxGradCc_;

            //- Minimum  gradcc for fully scheme 1
            scalar minGradCc_;

            //- Maximum ratio to average residual for scheme 2
            scalar maxResidual_;

            //- Minimum face weight for fully scheme 2
            scalar minFaceWeight_;

            //- Maximum face weight for fully scheme 1
            scalar maxFaceWeight_;

            //- Maximum skewness for fully scheme 2
            scalar maxSkewness_;

            //- Minimum skewness for fully scheme 1
            scalar minSkewness_;

            //- Maximum Co for fully scheme 2
            scalar Co1_;

            //- Minimum Co for fully scheme 1
            scalar Co2_;


        // File names

405
            //- Name of the non-orthogonality field
406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445
            word nonOrthogonalityName_;

            //- Name of the face weight field
            word faceWeightName_;

            //- Name of the skewnes field
            word skewnessName_;

            //- Name of the residual field
            word residualName_;

            //- Name of the U used for Co based blended
            word UName_;


        //- Tolerance used when calculating the number of blended cells
        scalar tolerance_;


        //- Error fields
        scalarField error_;
        scalarField errorIntegral_;
        scalarField oldError_;
        scalarField oldErrorIntegral_;

        //- Proportional gain
        scalar P_;

        //- Integral gain
        scalar I_;

        //- Derivative gain
        scalar D_;


    // Private Member Functions

        //- Init fields
        bool init(bool first);

446 447 448
        //- Calculate statistics
        void calcStats(label&, label&, label&) const ;

449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476
        //- Calculate the blending factor field and return true if successful
        virtual bool calc();


protected:

    // Protected Member Functions

        //- Write the file header
        virtual void writeFileHeader(Ostream& os) const;


public:

    //- Runtime type information
    TypeName("stabilityBlendingFactor");


    // Constructors

        //- Construct from Time and dictionary
        stabilityBlendingFactor
        (
            const word& name,
            const Time& runTime,
            const dictionary& dict
        );

477 478 479 480 481 482
        //- No copy construct
        stabilityBlendingFactor(const stabilityBlendingFactor&) = delete;

        //- No copy assignment
        void operator=(const stabilityBlendingFactor&) = delete;

483 484

    //- Destructor
485
    virtual ~stabilityBlendingFactor() = default;
486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507


    // Member Functions

        //- Read the stabilityBlendingFactor data
        virtual bool read(const dictionary&);

        //- Write the stabilityBlendingFactor
        virtual bool write();
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace functionObjects
} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //