setRDeltaT.H 5.52 KB
Newer Older
Henry Weller's avatar
Henry Weller committed
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
Henry Weller's avatar
Henry Weller committed
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2013-2016 OpenFOAM Foundation
9
    Copyright (C) 2020 OpenCFD Ltd.
Henry Weller's avatar
Henry Weller committed
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
-------------------------------------------------------------------------------
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/>.

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

{
30
    volScalarField& rDeltaT = trDeltaT.ref();
31

Henry Weller's avatar
Henry Weller committed
32 33 34
    const dictionary& pimpleDict = pimple.dict();

    // Maximum flow Courant number
35
    scalar maxCo(pimpleDict.get<scalar>("maxCo"));
Henry Weller's avatar
Henry Weller committed
36 37

    // Maximum time scale
38
    scalar maxDeltaT(pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT));
Henry Weller's avatar
Henry Weller committed
39 40 41 42

    // Smoothing parameter (0-1) when smoothing iterations > 0
    scalar rDeltaTSmoothingCoeff
    (
43
        pimpleDict.getOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
Henry Weller's avatar
Henry Weller committed
44 45 46 47 48
    );

    // Damping coefficient (1-0)
    scalar rDeltaTDampingCoeff
    (
49
        pimpleDict.getOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
Henry Weller's avatar
Henry Weller committed
50 51 52 53
    );

    // Maximum change in cell temperature per iteration
    // (relative to previous value)
54
    scalar alphaTemp(pimpleDict.getOrDefault<scalar>("alphaTemp", 0.05));
Henry Weller's avatar
Henry Weller committed
55

56 57
    // Maximum change in cell concentration per iteration
    // (relative to reference value)
58
    scalar alphaY(pimpleDict.getOrDefault<scalar>("alphaY", 1.0));
Henry Weller's avatar
Henry Weller committed
59 60 61 62 63 64 65 66

    Info<< "Time scales min/max:" << endl;

    // Cache old reciprocal time scale field
    volScalarField rDeltaT0("rDeltaT0", rDeltaT);

    // Flow time scale
    {
67
        rDeltaT.ref() =
Henry Weller's avatar
Henry Weller committed
68
        (
69 70
            fvc::surfaceSum(mag(phi))()()
           /((2*maxCo)*mesh.V()*rho())
Henry Weller's avatar
Henry Weller committed
71 72 73 74 75 76
        );

        // Limit the largest time scale
        rDeltaT.max(1/maxDeltaT);

        Info<< "    Flow        = "
77 78
            << 1/gMax(rDeltaT.primitiveField()) << ", "
            << 1/gMin(rDeltaT.primitiveField()) << endl;
Henry Weller's avatar
Henry Weller committed
79 80
    }

81 82
    // Heat release rate time scale
    if (alphaTemp < 1)
Henry Weller's avatar
Henry Weller committed
83
    {
84
        volScalarField::Internal rDeltaTT
Henry Weller's avatar
Henry Weller committed
85
        (
86
            mag(Qdot)/(alphaTemp*rho*thermo.Cp()*T)
Henry Weller's avatar
Henry Weller committed
87 88 89
        );

        Info<< "    Temperature = "
90 91
            << 1/(gMax(rDeltaTT.field()) + VSMALL) << ", "
            << 1/(gMin(rDeltaTT.field()) + VSMALL) << endl;
Henry Weller's avatar
Henry Weller committed
92

93 94 95 96 97 98 99 100 101
        rDeltaT.ref() = max(rDeltaT(), rDeltaTT);
    }

    // Reaction rate time scale
    if (alphaY < 1)
    {
        dictionary Yref(pimpleDict.subDict("Yref"));

        volScalarField::Internal rDeltaTY
Henry Weller's avatar
Henry Weller committed
102
        (
103 104 105 106 107 108 109
            IOobject
            (
                "rDeltaTY",
                runTime.timeName(),
                mesh
            ),
            mesh,
110
            dimensionedScalar(rDeltaT.dimensions(), Zero)
Henry Weller's avatar
Henry Weller committed
111
        );
112 113 114 115 116 117 118 119 120 121 122 123

        bool foundY = false;

        forAll(Y, i)
        {
            if (i != inertIndex && composition.active(i))
            {
                volScalarField& Yi = Y[i];

                if (Yref.found(Yi.name()))
                {
                    foundY = true;
124
                    const scalar Yrefi = Yref.get<scalar>(Yi.name());
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

                    rDeltaTY.field() = max
                    (
                        mag
                        (
                            reaction->R(Yi)().source()
                           /((Yrefi*alphaY)*(rho*mesh.V()))
                        ),
                        rDeltaTY
                    );
                }
            }
        }

        if (foundY)
        {
            Info<< "    Composition = "
                << 1/(gMax(rDeltaTY.field()) + VSMALL) << ", "
                << 1/(gMin(rDeltaTY.field()) + VSMALL) << endl;

            rDeltaT.ref() = max(rDeltaT(), rDeltaTY);
        }
        else
        {
            IOWarningIn(args.executable().c_str(), Yref)
                << "Cannot find any active species in Yref " << Yref
                << endl;
        }
Henry Weller's avatar
Henry Weller committed
153 154 155 156 157 158
    }

    // Update tho boundary values of the reciprocal time-step
    rDeltaT.correctBoundaryConditions();

    // Spatially smooth the time scale field
159
    if (rDeltaTSmoothingCoeff < 1)
Henry Weller's avatar
Henry Weller committed
160 161 162 163 164 165 166 167 168
    {
        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
    }

    // Limit rate of change of time scale
    // - reduce as much as required
    // - only increase at a fraction of old time scale
    if
    (
169
        rDeltaTDampingCoeff < 1
Henry Weller's avatar
Henry Weller committed
170 171 172 173 174 175
     && runTime.timeIndex() > runTime.startTimeIndex() + 1
    )
    {
        rDeltaT = max
        (
            rDeltaT,
176
            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
Henry Weller's avatar
Henry Weller committed
177 178 179
        );
    }

180 181 182
    // Update tho boundary values of the reciprocal time-step
    rDeltaT.correctBoundaryConditions();

Henry Weller's avatar
Henry Weller committed
183
    Info<< "    Overall     = "
184 185
        << 1/gMax(rDeltaT.primitiveField())
        << ", " << 1/gMin(rDeltaT.primitiveField()) << endl;
Henry Weller's avatar
Henry Weller committed
186 187 188 189
}


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