problem.hh 16 KB
Newer Older
Thomas Fetzer's avatar
Thomas Fetzer committed
1
2
3
4
5
6
7
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
Thomas Fetzer's avatar
Thomas Fetzer committed
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 * \ingroup RANSTests
 * \brief Pipe flow test for the staggered grid RANS model
 *
24
25
 * This test simulates pipe flow experiments performed by John Laufer in 1954
 * \cite Laufer1954a.
Thomas Fetzer's avatar
Thomas Fetzer committed
26
27
28
29
 */
#ifndef DUMUX_PIPE_LAUFER_PROBLEM_HH
#define DUMUX_PIPE_LAUFER_PROBLEM_HH

30
#include <dumux/common/properties.hh>
31
32
#include <dumux/common/parameters.hh>
#include <dumux/common/timeloop.hh>
33
#include <dumux/common/numeqvector.hh>
34

35
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
36
#include <dumux/freeflow/rans/problem.hh>
37
#include <dumux/freeflow/turbulencemodel.hh>
38
#include <dumux/freeflow/turbulenceproperties.hh>
Thomas Fetzer's avatar
Thomas Fetzer committed
39

40
namespace Dumux {
41

Thomas Fetzer's avatar
Thomas Fetzer committed
42
43
44
45
46
47
48
49
/*!
 * \ingroup NavierStokesTests
 * \brief  Test problem for the one-phase (Navier-) Stokes problem in a channel.
 *
 * This test simulates is based on pipe flow experiments by
 * John Laufers experiments in 1954 \cite Laufer1954a.
 */
template <class TypeTag>
50
class PipeLauferProblem : public RANSProblem<TypeTag>
51
{
52
    using ParentType = RANSProblem<TypeTag>;
Thomas Fetzer's avatar
Thomas Fetzer committed
53

54
    using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
55
56
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using FluidState = GetPropType<TypeTag, Properties::FluidState>;
57
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
58
59
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
60
    using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
61
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
Thomas Fetzer's avatar
Thomas Fetzer committed
62

63
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
64
    using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
65
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
66
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
67
68
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
Thomas Fetzer's avatar
Thomas Fetzer committed
69

Thomas Fetzer's avatar
Thomas Fetzer committed
70
71
    using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;

72
    static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
73

Thomas Fetzer's avatar
Thomas Fetzer committed
74
public:
75
76
    PipeLauferProblem(std::shared_ptr<const GridGeometry> gridGeometry)
    : ParentType(gridGeometry), eps_(1e-6)
Thomas Fetzer's avatar
Thomas Fetzer committed
77
78
    {
        inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
79
        inletTemperature_ = getParam<Scalar>("Problem.InletTemperature", 283.15);
80
        wallTemperature_ = getParam<Scalar>("Problem.WallTemperature", 303.15);
81
82
83
84

        FluidSystem::init();
        Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
        FluidState fluidState;
Timo Koch's avatar
Timo Koch committed
85
        fluidState.setPressure(0, 1e5);
86
        fluidState.setTemperature(temperature());
Timo Koch's avatar
Timo Koch committed
87
88
        Scalar density = FluidSystem::density(fluidState, 0);
        Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, 0) / density;
89
        Scalar diameter = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
90

91
        // ideally the viscosityTilde parameter as inflow for the Spalart-Allmaras model should be zero
92
93
94
        viscosityTilde_ = 1e-3 * turbulenceProperties.viscosityTilde(inletVelocity_, diameter, kinematicViscosity);
        turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(inletVelocity_, diameter, kinematicViscosity);

95
        if (ModelTraits::turbulenceModel() == TurbulenceModel::komega || ModelTraits::turbulenceModel() == TurbulenceModel::sst)
96
97
98
99
100
101
102
103
            dissipation_ = turbulenceProperties.dissipationRate(inletVelocity_, diameter, kinematicViscosity);
        else
            dissipation_ = turbulenceProperties.dissipation(inletVelocity_, diameter, kinematicViscosity);

        if (ModelTraits::turbulenceModel() == TurbulenceModel::oneeq)
            initializationTime_ = getParam<Scalar>("TimeLoop.Initialization", 1.0);
        else
            initializationTime_ = getParam<Scalar>("TimeLoop.Initialization", -1.0);
104
105
106

        turbulenceModelName_ = turbulenceModelToString(ModelTraits::turbulenceModel());
        std::cout << "Using the "<< turbulenceModelName_ << " Turbulence Model. \n";
107
        std::cout << std::endl;
Thomas Fetzer's avatar
Thomas Fetzer committed
108
109
110
111
112
113
114
    }

   /*!
     * \name Problem parameters
     */
    // \{

115
    bool isOnWallAtPos(const GlobalPosition &globalPos) const
116
    {
117
118
        return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_
               || globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_;
119
    }
Thomas Fetzer's avatar
Thomas Fetzer committed
120
121

   /*!
122
     * \brief Returns the temperature [K] within the domain for the isothermal model.
Thomas Fetzer's avatar
Thomas Fetzer committed
123
124
     */
    Scalar temperature() const
125
    { return inletTemperature_; }
Thomas Fetzer's avatar
Thomas Fetzer committed
126
127

    // \}
128

Thomas Fetzer's avatar
Thomas Fetzer committed
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
   /*!
     * \name Boundary conditions
     */
    // \{

   /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary control volume.
     *
     * \param globalPos The position of the center of the finite volume
     */
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
        BoundaryTypes values;

144
        // common boundary types for all turbulence models
145
        if(isOutlet_(globalPos))
146
            values.setDirichlet(Indices::pressureIdx);
147
        else
148
        {
149
            // walls and inflow
150
151
            values.setDirichlet(Indices::velocityXIdx);
            values.setDirichlet(Indices::velocityYIdx);
152
        }
153
154

#if NONISOTHERMAL
155
156
157
        if(isOutlet_(globalPos))
            values.setOutflow(Indices::energyEqIdx);
        else
158
            values.setDirichlet(Indices::temperatureIdx);
159
#endif
160
        // turbulence model-specific boundary types
161
        setBcTypes_(values, globalPos);
162

Thomas Fetzer's avatar
Thomas Fetzer committed
163
164
165
        return values;
    }

166
167
168
169
170
171
    /*!
     * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
     *
     * \param element The finite element
     * \param fvGeometry The finite-volume geometry
     * \param scv The sub control volume
172
     * \param pvIdx The primary variable index in the solution vector
173
174
175
176
177
178
     */
    template<class Element, class FVElementGeometry, class SubControlVolume>
    bool isDirichletCell(const Element& element,
                         const FVElementGeometry& fvGeometry,
                         const SubControlVolume& scv,
                         int pvIdx) const
179
    { return isDirichletCell_(element, fvGeometry, scv, pvIdx); }
180

181
     /*!
182
      * \brief Evaluate the boundary conditions for a dirichlet values at the boundary.
183
184
185
186
187
      *
      * \param element The finite element
      * \param scvf the sub control volume face
      * \note used for cell-centered discretization schemes
      */
188
    PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
Thomas Fetzer's avatar
Thomas Fetzer committed
189
    {
190
        const auto globalPos = scvf.ipGlobal();
191
        PrimaryVariables values(initialAtPos(globalPos));
192

193
#if NONISOTHERMAL
194
        values[Indices::temperatureIdx] = isOnWallAtPos(globalPos) ? wallTemperature_ : inletTemperature_;
195
#endif
196

197
        return values;
Thomas Fetzer's avatar
Thomas Fetzer committed
198
199
    }

200
201
202
203
204
205
206
    /*!
     * \brief Evaluate the boundary conditions for fixed values at cell centers
     *
     * \param element The finite element
     * \param scv the sub control volume
     * \note used for cell-centered discretization schemes
     */
207
    PrimaryVariables dirichlet([[maybe_unused]] const Element& element, const SubControlVolume& scv) const
208
    {
209
        if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::kepsilon
210
211
                   || ModelTraits::turbulenceModel() == TurbulenceModel::komega
                   || ModelTraits::turbulenceModel() == TurbulenceModel::sst)
212
213
214
215
216
217
218
            return dirichletTurbulentTwoEq_(element, scv);
        else
        {
            const auto globalPos = scv.center();
            PrimaryVariables values(initialAtPos(globalPos));
            return values;
        }
219
    }
Thomas Fetzer's avatar
Thomas Fetzer committed
220
221

   /*!
222
     * \brief Evaluate the initial value for a control volume.
Thomas Fetzer's avatar
Thomas Fetzer committed
223
224
225
     *
     * \param globalPos The global position
     */
226
    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
Thomas Fetzer's avatar
Thomas Fetzer committed
227
    {
228
        // common initial conditions for all turbulence models
229
        PrimaryVariables values(0.0);
230
        values[Indices::pressureIdx] = 1.0e+5;
231
232
233
        values[Indices::velocityXIdx] = time() > initializationTime_
                                        ? inletVelocity_
                                        : time() / initializationTime_ * inletVelocity_;
234
        if (isOnWallAtPos(globalPos))
235
            values[Indices::velocityXIdx] = 0.0;
236
#if NONISOTHERMAL
237
        values[Indices::temperatureIdx] = isOnWallAtPos(globalPos) ? wallTemperature_ : inletTemperature_;
238
#endif
239
        // turbulence model-specific initial conditions
240
        setInitialAtPos_(values, globalPos);
Thomas Fetzer's avatar
Thomas Fetzer committed
241
242
243
244
        return values;
    }

    // \}
245

Thomas Fetzer's avatar
Thomas Fetzer committed
246
    void setTimeLoop(TimeLoopPtr timeLoop)
247
    { timeLoop_ = timeLoop; }
Thomas Fetzer's avatar
Thomas Fetzer committed
248
249

    Scalar time() const
250
    { return timeLoop_->time(); }
Thomas Fetzer's avatar
Thomas Fetzer committed
251
252

private:
253
    bool isInlet_(const GlobalPosition& globalPos) const
254
    { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
255

256
    bool isOutlet_(const GlobalPosition& globalPos) const
257
    { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
Thomas Fetzer's avatar
Thomas Fetzer committed
258

259
    //! Initial conditions for the komega, kepsilon and lowrekepsilon turbulence models
260
261
    void setInitialAtPos_([[maybe_unused]] PrimaryVariables& values,
                          [[maybe_unused]] const GlobalPosition &globalPos) const
262
    {
263
264
265
        if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 0) // zero equation models
            return;
        else if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 1)  // one equation models
266
        {
267
268
269
270
271
272
273
274
275
276
277
278
279
280
            values[Indices::viscosityTildeIdx] = viscosityTilde_;
            if (isOnWallAtPos(globalPos))
                values[Indices::viscosityTildeIdx] = 0.0;
        }
        else // two equation models
        {
            static_assert(numTurbulenceEq(ModelTraits::turbulenceModel()) == 2, "Only reached by 2eq models");
            values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_;
            values[Indices::dissipationIdx] = dissipation_;
            if (isOnWallAtPos(globalPos))
            {
                values[Indices::turbulentKineticEnergyIdx] = 0.0;
                values[Indices::dissipationIdx] = 0.0;
            }
281
282
283
284
        }
    }

    //! Boundary condition types for the one-eq turbulence model
285
286
    void setBcTypes_([[maybe_unused]] BoundaryTypes& values,
                     [[maybe_unused]] const GlobalPosition& pos) const
287
    {
288
289
290
        if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 0) // zero equation models
            return;
        else if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 1)  // one equation models
291
        {
292
293
294
295
            if(isOutlet_(pos))
                values.setOutflow(Indices::viscosityTildeIdx);
            else // walls and inflow
                values.setDirichlet(Indices::viscosityTildeIdx);
296
        }
297
        else // two equation models
298
        {
299
300
301
302
303
304
305
306
307
308
309
310
            static_assert(numTurbulenceEq(ModelTraits::turbulenceModel()) == 2, "Only reached by 2eq models");
            if(isOutlet_(pos))
            {
                values.setOutflow(Indices::turbulentKineticEnergyEqIdx);
                values.setOutflow(Indices::dissipationEqIdx);
            }
            else
            {
                // walls and inflow
                values.setDirichlet(Indices::turbulentKineticEnergyIdx);
                values.setDirichlet(Indices::dissipationIdx);
            }
311
312
313
314
        }
    }

    template<class Element, class FVElementGeometry, class SubControlVolume>
315
    bool isDirichletCell_([[maybe_unused]] const Element& element,
316
                          const FVElementGeometry& fvGeometry,
317
                          [[maybe_unused]] const SubControlVolume& scv,
318
                          const int& pvIdx) const
319
    {
320
321
        if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::kepsilon)
        {
322
            const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element);
323
324
325
326
327
328
329
            // For the kepsilon model we set fixed values within the matching point and at the wall
            if (this->inNearWallRegion(eIdx))
                return pvIdx == Indices::turbulentKineticEnergyEqIdx || pvIdx == Indices::dissipationEqIdx;
            if (this->isMatchingPoint(eIdx))
                return pvIdx == Indices::dissipationEqIdx;
            return false;
        }
330
331
        else if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::komega ||
                           ModelTraits::turbulenceModel() == TurbulenceModel::sst)
332
333
334
335
336
337
338
339
340
        {
            // For the komega model we set a fixed dissipation (omega) for all cells at the wall
            for (const auto& scvf : scvfs(fvGeometry))
                if (isOnWallAtPos(scvf.center()) && pvIdx == Indices::dissipationIdx)
                    return true;
            return false;
        }
        else
            return ParentType::isDirichletCell(element, fvGeometry, scv, pvIdx);
341
342
    }

343
    //! Specialization for the kepsilon and komega
344
345
    template<class Element, class SubControlVolume>
    PrimaryVariables dirichletTurbulentTwoEq_(const Element& element,
346
                                              const SubControlVolume& scv) const
347
348
349
    {
        const auto globalPos = scv.center();
        PrimaryVariables values(initialAtPos(globalPos));
350
        unsigned int  elementIdx = this->gridGeometry().elementMapper().index(element);
351

352
353
354
355
356
357
358
359
360
        if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::kepsilon)
        {
            // For the kepsilon model we set a fixed value for the turbulent kinetic energy and the dissipation
            values[Indices::turbulentKineticEnergyEqIdx] = this->turbulentKineticEnergyWallFunction(elementIdx);
            values[Indices::dissipationEqIdx] = this->dissipationWallFunction(elementIdx);
            return values;
        }
        else
        {
361
362
            static_assert(ModelTraits::turbulenceModel() == TurbulenceModel::komega
                       || ModelTraits::turbulenceModel() == TurbulenceModel::sst, "Only valid for Komega");
363
            // For the komega model we set a fixed value for the dissipation
364
            const auto wallDistance = ParentType::wallDistance(elementIdx);
365
            using std::pow;
366
            values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
367
368
369
                                                    / (ParentType::betaOmega() * wallDistance * wallDistance);
            return values;
        }
370
371
    }

Thomas Fetzer's avatar
Thomas Fetzer committed
372
373
    Scalar eps_;
    Scalar inletVelocity_;
374
375
    Scalar inletTemperature_;
    Scalar wallTemperature_;
376
    Scalar initializationTime_;
377
378
379
    Scalar viscosityTilde_;
    Scalar turbulentKineticEnergy_;
    Scalar dissipation_;
380
    std::string turbulenceModelName_;
Thomas Fetzer's avatar
Thomas Fetzer committed
381
382
    TimeLoopPtr timeLoop_;
};
383
} // end namespace Dumux
Thomas Fetzer's avatar
Thomas Fetzer committed
384
385

#endif