problem_stokes.hh 14 KB
Newer Older
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       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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
21
 * \ingroup BoundaryTests
22
23
 * \brief A simple Navier-Stokes test problem for the staggered grid (Navier-)Stokes model.
 */
24

25
26
27
28
29
30
#ifndef DUMUX_STOKES_SUBPROBLEM_HH
#define DUMUX_STOKES_SUBPROBLEM_HH

#include <dune/grid/yaspgrid.hh>

#include <dumux/discretization/staggered/freeflow/properties.hh>
31

32
#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
33
34
35
36
37
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/navierstokes/problem.hh>

#include <dumux/material/fluidsystems/1padapter.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
38

39
namespace Dumux {
40
41
42
template <class TypeTag>
class StokesSubProblem;

43
namespace Properties {
44
45
46
47
// Create new type tags
namespace TTag {
struct StokesOnePTwoC { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
} // end namespace TTag
48

49
// The fluid system
50
51
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::StokesOnePTwoC>
52
{
53
  using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
54
55
56
  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
};
57
58

// Set the grid type
59
60
template<class TypeTag>
struct Grid<TypeTag, TTag::StokesOnePTwoC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
61
62

// Set the problem property
63
64
template<class TypeTag>
struct Problem<TypeTag, TTag::StokesOnePTwoC> { using type = Dumux::StokesSubProblem<TypeTag> ; };
65

66
template<class TypeTag>
67
struct EnableGridGeometryCache<TypeTag, TTag::StokesOnePTwoC> { static constexpr bool value = true; };
68
69
70
71
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOnePTwoC> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOnePTwoC> { static constexpr bool value = true; };
72
73

// Use moles
74
75
template<class TypeTag>
struct UseMoles<TypeTag, TTag::StokesOnePTwoC> { static constexpr bool value = true; };
76
77

// Do not replace one equation with a total mass balance
78
79
template<class TypeTag>
struct ReplaceCompEqIdx<TypeTag, TTag::StokesOnePTwoC> { static constexpr int value = 3; };
80
} // end namespace Properties
81
82

/*!
83
84
 * \ingroup BoundaryTests
 * \brief Test problem for the one-phase (Navier-) Stokes problem.
85
86
87
88
89
90
91
 *
 * Horizontal flow from left to right with a parabolic velocity profile.
 */
template <class TypeTag>
class StokesSubProblem : public NavierStokesProblem<TypeTag>
{
    using ParentType = NavierStokesProblem<TypeTag>;
92
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
93
94
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
95
    using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
96
97
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
98
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
99
100
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
101
102
103
104

    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;

105
    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
106
    using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
107
108

public:
109
110
    StokesSubProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
    : ParentType(gridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
111
    {
112
        problemName_  =  getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
113
114
115

        // determine whether to simulate a vertical or horizontal flow configuration
        verticalFlow_ = problemName_.find("vertical") != std::string::npos;
116
117
118
119
120
121
122
123
    }

    /*!
     * \brief The problem name.
     */
    const std::string& name() const
    {
        return problemName_;
124
125
126
127
128
129
130
131
    }

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

   /*!
132
     * \brief Returns the temperature within the domain in [K].
133
134
135
136
137
138
139
     *
     * This problem assumes a temperature of 10 degrees Celsius.
     */
    Scalar temperature() const
    { return 273.15 + 10; } // 10°C

   /*!
140
     * \brief Returns the sources within the domain.
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
     *
     * \param globalPos The global position
     */
    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
    { return NumEqVector(0.0); }
    // \}

   /*!
     * \name Boundary conditions
     */
    // \{

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
     * \param element The finite element
     * \param scvf The sub control volume face
     */
    BoundaryTypes boundaryTypes(const Element& element,
                                const SubControlVolumeFace& scvf) const
    {
        BoundaryTypes values;

        const auto& globalPos = scvf.dofPosition();

167
        if (verticalFlow_)
168
        {
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
            // inflow
            if(onUpperBoundary_(globalPos))
            {
                values.setDirichlet(Indices::velocityXIdx);
                values.setDirichlet(Indices::velocityYIdx);
                values.setDirichlet(Indices::conti0EqIdx + 1);
            }

            // left/right wall
            if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
            {
                values.setDirichlet(Indices::velocityXIdx);
                values.setDirichlet(Indices::velocityYIdx);
                values.setNeumann(Indices::conti0EqIdx);
                values.setNeumann(Indices::conti0EqIdx + 1);
            }

186
        }
187
        else // horizontal flow
188
        {
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
            if (onLeftBoundary_(globalPos))
            {
                values.setDirichlet(Indices::conti0EqIdx + 1);
                values.setDirichlet(Indices::velocityXIdx);
                values.setDirichlet(Indices::velocityYIdx);
            }
            else if (onRightBoundary_(globalPos))
            {
                values.setDirichlet(Indices::pressureIdx);
                values.setOutflow(Indices::conti0EqIdx + 1);
            }
            else
            {
                values.setDirichlet(Indices::velocityXIdx);
                values.setDirichlet(Indices::velocityYIdx);
                values.setNeumann(Indices::conti0EqIdx);
                values.setNeumann(Indices::conti0EqIdx + 1);
            }
207
208
209
210
211
212
213
        }

        if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
        {
            values.setCouplingNeumann(Indices::conti0EqIdx);
            values.setCouplingNeumann(Indices::conti0EqIdx + 1);
            values.setCouplingNeumann(Indices::momentumYBalanceIdx);
214
            values.setBeaversJoseph(Indices::momentumXBalanceIdx);
215
216
217
218
219
220
        }

        return values;
    }

    /*!
221
     * \brief Evaluates the boundary conditions for a Dirichlet control volume.
222
223
224
225
226
227
     */
    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
    {
        PrimaryVariables values(0.0);
        values = initialAtPos(globalPos);

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
        if (verticalFlow_)
        {
            // Check if this a pure diffusion problem.
            static const bool isDiffusionProblem = problemName_.find("diffusion") != std::string::npos;

            Scalar topMoleFraction = 1e-3;

            if (isDiffusionProblem)
            {
                // For the diffusion problem, change the top mole fraction after some time
                // in order to revert the concentration gradient.
                if (time() >= 1e10)
                    topMoleFraction = 0.0;
            }
            else // advection problem
            {
                // reverse the flow direction after some time for the advection problem
                if (time() >= 3e5)
                    values[Indices::velocityYIdx] *= -1.0;
            }

249
            if(globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_)
250
251
252
253
254
                values[Indices::conti0EqIdx + 1] = topMoleFraction;
        }
        else // horizontal flow
        {
            static const Scalar inletMoleFraction = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InletMoleFraction");
255
            if(globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_)
256
257
258
                values[Indices::conti0EqIdx + 1] = inletMoleFraction;
        }

259
260
261
262
263

        return values;
    }

    /*!
264
     * \brief Evaluates the boundary conditions for a Neumann control volume.
265
266
     *
     * \param element The element for which the Neumann boundary condition is set
267
     * \param fvGeometry The fvGeometry
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
     * \param elemVolVars The element volume variables
     * \param elemFaceVars The element face variables
     * \param scvf The boundary sub control volume face
     */
    template<class ElementVolumeVariables, class ElementFaceVariables>
    NumEqVector neumann(const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const ElementFaceVariables& elemFaceVars,
                        const SubControlVolumeFace& scvf) const
    {
        NumEqVector values(0.0);

        if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
        {
283
            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
284

285
            const auto tmp = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
            values[Indices::conti0EqIdx] = tmp[0];
            values[Indices::conti0EqIdx + 1] = tmp[1];
        }
        return values;
    }

    // \}

    //! Get the coupling manager
    const CouplingManager& couplingManager() const
    { return *couplingManager_; }

   /*!
     * \name Volume terms
     */
    // \{

   /*!
304
     * \brief Evaluates the initial value for a control volume.
305
306
307
     *
     * \param globalPos The global position
     */
308
    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
309
310
    {
        PrimaryVariables values(0.0);
311
312
313
314
315
316
        values[Indices::pressureIdx] = 1e5;

        static const Scalar vMax = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity", 0.0);

        auto parabolicProfile = [&](const GlobalPosition& globalPos, int coord)
        {
317
318
319
320
            return vMax * (globalPos[coord] - this->gridGeometry().bBoxMin()[coord])
                        * (this->gridGeometry().bBoxMax()[coord] - globalPos[coord])
                        / (0.25 * (this->gridGeometry().bBoxMax()[coord] - this->gridGeometry().bBoxMin()[coord])
                        * (this->gridGeometry().bBoxMax()[coord] - this->gridGeometry().bBoxMin()[coord]));
321
322
323
324
325
326
        };

        if (verticalFlow_)
            values[Indices::velocityYIdx] = parabolicProfile(globalPos, 0);
        else // horizontal flow
            values[Indices::velocityXIdx] = parabolicProfile(globalPos, 1);
327
328
329
330
331

        return values;
    }

    /*!
332
333
     * \brief Returns the intrinsic permeability of required as input parameter
     *        for the Beavers-Joseph-Saffman boundary condition.
334
     */
335
    Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
336
    {
337
        return couplingManager().couplingData().darcyPermeability(element, scvf);
338
339
340
    }

    /*!
341
342
     * \brief Returns the alpha value required as input parameter for the
     *        Beavers-Joseph-Saffman boundary condition.
343
344
345
346
347
348
     */
    Scalar alphaBJ(const SubControlVolumeFace& scvf) const
    {
        return 1.0;
    }

349
    /*!
350
     * \brief Sets the time loop pointer.
351
352
353
     */
    void setTimeLoop(TimeLoopPtr timeLoop)
    { timeLoop_ = timeLoop; }
354

355
    /*!
356
     * \brief Returns the time.
357
358
359
     */
    Scalar time() const
    { return timeLoop_->time(); }
360
361
362
363
364

    // \}

private:
    bool onLeftBoundary_(const GlobalPosition &globalPos) const
365
    { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
366
367

    bool onRightBoundary_(const GlobalPosition &globalPos) const
368
    { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
369
370

    bool onLowerBoundary_(const GlobalPosition &globalPos) const
371
    { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
372
373

    bool onUpperBoundary_(const GlobalPosition &globalPos) const
374
    { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
375
376

    Scalar eps_;
377
    bool verticalFlow_;
378
    std::string problemName_;
379
    std::shared_ptr<CouplingManager> couplingManager_;
380
    TimeLoopPtr timeLoop_;
381
};
382
} // end namespace Dumux
383
384

#endif // DUMUX_STOKES_SUBPROBLEM_HH