nutUSpaldingWallFunctionFvPatchScalarField.C 9.61 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
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2016 OpenFOAM Foundation
9
    Copyright (C) 2019-2020 OpenCFD Ltd.
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
-------------------------------------------------------------------------------
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 "nutUSpaldingWallFunctionFvPatchScalarField.H"
#include "turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"


// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

38 39
Foam::tmp<Foam::scalarField>
Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcNut() const
40 41 42 43 44 45 46 47
{
    const label patchi = patch().index();

    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
    (
        IOobject::groupName
        (
            turbulenceModel::propertiesName,
48
            internalField().group()
49 50
        )
    );
51
    const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
52 53 54 55
    const scalarField magGradU(mag(Uw.snGrad()));
    const tmp<scalarField> tnuw = turbModel.nu(patchi);
    const scalarField& nuw = tnuw();

mattijs's avatar
mattijs committed
56 57 58

    // Calculate new viscosity
    tmp<scalarField> tnutw
59
    (
mattijs's avatar
mattijs committed
60 61 62 63 64
        max
        (
            scalar(0),
            sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw
        )
65
    );
mattijs's avatar
mattijs committed
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88

    if (tolerance_ != 0.01)
    {
        // User-specified tolerance. Find out if current nut already satisfies
        // eqns.

        // Run ut for one iteration
        scalarField err;
        tmp<scalarField> UTau(calcUTau(magGradU, 1, err));

        // Preserve nutw if the initial error (err) already lower than the
        // tolerance.

        scalarField& nutw = tnutw.ref();
        forAll(err, facei)
        {
            if (err[facei] < tolerance_)
            {
                nutw[facei] = this->operator[](facei);
            }
        }
    }
    return tnutw;
89 90 91
}


92 93
Foam::tmp<Foam::scalarField>
Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
94 95 96
(
    const scalarField& magGradU
) const
mattijs's avatar
mattijs committed
97 98 99 100 101 102
{
    scalarField err;
    return calcUTau(magGradU, maxIter_, err);
}


103 104
Foam::tmp<Foam::scalarField>
Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
mattijs's avatar
mattijs committed
105 106 107 108 109
(
    const scalarField& magGradU,
    const label maxIter,
    scalarField& err
) const
110 111 112 113 114 115 116 117
{
    const label patchi = patch().index();

    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
    (
        IOobject::groupName
        (
            turbulenceModel::propertiesName,
118
            internalField().group()
119 120 121 122
        )
    );
    const scalarField& y = turbModel.y()[patchi];

123
    const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
124 125 126 127 128 129 130
    const scalarField magUp(mag(Uw.patchInternalField() - Uw));

    const tmp<scalarField> tnuw = turbModel.nu(patchi);
    const scalarField& nuw = tnuw();

    const scalarField& nutw = *this;

131
    tmp<scalarField> tuTau(new scalarField(patch().size(), Zero));
132
    scalarField& uTau = tuTau.ref();
133

mattijs's avatar
mattijs committed
134 135 136
    err.setSize(uTau.size());
    err = 0.0;

137
    forAll(uTau, facei)
138
    {
139
        scalar ut = sqrt((nutw[facei] + nuw[facei])*magGradU[facei]);
140 141
        // Note: for exact restart seed with laminar viscosity only:
        //scalar ut = sqrt(nuw[facei]*magGradU[facei]);
142

143
        if (ROOTVSMALL < ut)
144 145 146 147 148
        {
            int iter = 0;

            do
            {
149
                scalar kUu = min(kappa_*magUp[facei]/ut, 50);
150 151 152
                scalar fkUu = exp(kUu) - 1 - kUu*(1 + 0.5*kUu);

                scalar f =
153 154
                    - ut*y[facei]/nuw[facei]
                    + magUp[facei]/ut
155 156 157
                    + 1/E_*(fkUu - 1.0/6.0*kUu*sqr(kUu));

                scalar df =
158 159
                    y[facei]/nuw[facei]
                  + magUp[facei]/sqr(ut)
160 161 162
                  + 1/E_*kUu*fkUu/ut;

                scalar uTauNew = ut + f/df;
mattijs's avatar
mattijs committed
163
                err[facei] = mag((ut - uTauNew)/ut);
164 165
                ut = uTauNew;

166
                //iterations_++;
167

mattijs's avatar
mattijs committed
168 169 170 171 172 173
            } while
            (
                ut > ROOTVSMALL
             && err[facei] > tolerance_
             && ++iter < maxIter
            );
174

175
            uTau[facei] = max(0.0, ut);
176

177 178 179 180 181 182 183 184 185
            //invocations_++;
            //if (iter > 1)
            //{
            //    nontrivial_++;
            //}
            //if (iter >= maxIter_)
            //{
            //    nonconvergence_++;
            //}
186 187 188 189 190 191 192
        }
    }

    return tuTau;
}


193 194 195 196 197 198 199 200 201 202 203 204
void Foam::nutUSpaldingWallFunctionFvPatchScalarField::writeLocalEntries
(
    Ostream& os
) const
{
    nutWallFunctionFvPatchScalarField::writeLocalEntries(os);

    os.writeEntryIfDifferent<label>("maxIter", 10, maxIter_);
    os.writeEntryIfDifferent<scalar>("tolerance", 0.01, tolerance_);
}


205 206
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

207
Foam::nutUSpaldingWallFunctionFvPatchScalarField::
208 209 210 211 212 213
nutUSpaldingWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
214
    nutWallFunctionFvPatchScalarField(p, iF),
215 216 217 218 219 220
    maxIter_(10),
    tolerance_(0.01)
    //invocations_(0),
    //nontrivial_(0),
    //nonconvergence_(0),
    //iterations_(0)
221 222 223
{}


224
Foam::nutUSpaldingWallFunctionFvPatchScalarField::
225 226 227 228 229 230 231 232
nutUSpaldingWallFunctionFvPatchScalarField
(
    const nutUSpaldingWallFunctionFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
233
    nutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
234 235 236 237 238 239
    maxIter_(ptf.maxIter_),
    tolerance_(ptf.tolerance_)
    //invocations_(0),
    //nontrivial_(0),
    //nonconvergence_(0),
    //iterations_(0)
240 241 242
{}


243
Foam::nutUSpaldingWallFunctionFvPatchScalarField::
244 245 246 247 248 249 250
nutUSpaldingWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
251
    nutWallFunctionFvPatchScalarField(p, iF, dict),
252 253
    maxIter_(dict.getOrDefault<label>("maxIter", 10)),
    tolerance_(dict.getOrDefault<scalar>("tolerance", 0.01))
254 255 256 257
    //invocations_(0),
    //nontrivial_(0),
    //nonconvergence_(0),
    //iterations_(0)
258 259 260
{}


261
Foam::nutUSpaldingWallFunctionFvPatchScalarField::
262 263 264 265 266
nutUSpaldingWallFunctionFvPatchScalarField
(
    const nutUSpaldingWallFunctionFvPatchScalarField& wfpsf
)
:
267
    nutWallFunctionFvPatchScalarField(wfpsf),
268 269 270 271 272 273
    maxIter_(wfpsf.maxIter_),
    tolerance_(wfpsf.tolerance_)
    //invocations_(wfpsf.invocations_),
    //nontrivial_(wfpsf.nontrivial_),
    //nonconvergence_(wfpsf.nonconvergence_),
    //iterations_(wfpsf.iterations_)
274 275 276
{}


277
Foam::nutUSpaldingWallFunctionFvPatchScalarField::
278 279 280 281 282 283
nutUSpaldingWallFunctionFvPatchScalarField
(
    const nutUSpaldingWallFunctionFvPatchScalarField& wfpsf,
    const DimensionedField<scalar, volMesh>& iF
)
:
284
    nutWallFunctionFvPatchScalarField(wfpsf, iF),
285 286 287 288 289 290
    maxIter_(wfpsf.maxIter_),
    tolerance_(wfpsf.tolerance_)
    //invocations_(0),
    //nontrivial_(0),
    //nonconvergence_(0),
    //iterations_(0)
291 292 293
{}


294 295
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

296
Foam::nutUSpaldingWallFunctionFvPatchScalarField::
297 298
~nutUSpaldingWallFunctionFvPatchScalarField()
{
299 300 301 302 303 304 305 306 307 308 309 310 311
    //if (debug)
    //{
    //    Info<< "nutUSpaldingWallFunctionFvPatchScalarField :"
    //        << " total invocations:"
    //        << returnReduce(invocations_, sumOp<label>())
    //        << " total iterations:"
    //        << returnReduce(iterations_, sumOp<label>())
    //        << " total non-convergence:"
    //        << returnReduce(nonconvergence_, sumOp<label>())
    //        << " total non-trivial:"
    //        << returnReduce(nontrivial_, sumOp<label>())
    //        << endl;
    //}
312 313 314
}


315 316
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

317 318
Foam::tmp<Foam::scalarField>
Foam::nutUSpaldingWallFunctionFvPatchScalarField::yPlus() const
319 320 321 322 323 324 325 326
{
    const label patchi = patch().index();

    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
    (
        IOobject::groupName
        (
            turbulenceModel::propertiesName,
327
            internalField().group()
328 329 330
        )
    );
    const scalarField& y = turbModel.y()[patchi];
331
    const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
332 333 334 335 336 337 338
    const tmp<scalarField> tnuw = turbModel.nu(patchi);
    const scalarField& nuw = tnuw();

    return y*calcUTau(mag(Uw.snGrad()))/nuw;
}


339 340 341 342
void Foam::nutUSpaldingWallFunctionFvPatchScalarField::write
(
    Ostream& os
) const
343 344 345 346 347 348 349 350 351
{
    fvPatchField<scalar>::write(os);
    writeLocalEntries(os);
    writeEntry("value", os);
}


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

352 353 354 355 356 357 358 359
namespace Foam
{
    makePatchTypeField
    (
        fvPatchScalarField,
        nutUSpaldingWallFunctionFvPatchScalarField
    );
}
360 361 362


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