problem_stokes.hh 12.2 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
#ifndef DUMUX_STOKES_SUBPROBLEM_HH
#define DUMUX_STOKES_SUBPROBLEM_HH

28
29
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
30
#include <dumux/common/numeqvector.hh>
31

32
33
34
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/navierstokes/problem.hh>

35
namespace Dumux {
36
37

/*!
38
39
 * \ingroup BoundaryTests
 * \brief Test problem for the one-phase (Navier-) Stokes problem.
40
41
42
43
44
45
46
 *
 * Horizontal flow from left to right with a parabolic velocity profile.
 */
template <class TypeTag>
class StokesSubProblem : public NavierStokesProblem<TypeTag>
{
    using ParentType = NavierStokesProblem<TypeTag>;
47
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
48
49
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
50
    using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
51
52
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
53
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
54
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
55
    using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
56
57
58
59

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

60
    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
61
    using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
62
63

public:
64
65
    StokesSubProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
    : ParentType(gridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
66
    {
67
        problemName_  =  getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
68
69
70

        // determine whether to simulate a vertical or horizontal flow configuration
        verticalFlow_ = problemName_.find("vertical") != std::string::npos;
71
72
73
74
75
76
77
78
    }

    /*!
     * \brief The problem name.
     */
    const std::string& name() const
    {
        return problemName_;
79
80
81
82
83
84
85
86
    }

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

   /*!
87
     * \brief Returns the temperature within the domain in [K].
88
89
90
91
92
93
94
     *
     * This problem assumes a temperature of 10 degrees Celsius.
     */
    Scalar temperature() const
    { return 273.15 + 10; } // 10°C

   /*!
95
     * \brief Returns the sources within the domain.
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
     *
     * \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();

122
        if (verticalFlow_)
123
        {
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
            // 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);
            }

141
        }
142
        else // horizontal flow
143
        {
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
            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);
            }
162
163
164
165
166
167
168
        }

        if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
        {
            values.setCouplingNeumann(Indices::conti0EqIdx);
            values.setCouplingNeumann(Indices::conti0EqIdx + 1);
            values.setCouplingNeumann(Indices::momentumYBalanceIdx);
169
            values.setBeaversJoseph(Indices::momentumXBalanceIdx);
170
171
172
173
174
175
        }

        return values;
    }

    /*!
176
     * \brief Evaluates the boundary conditions for a Dirichlet control volume.
177
178
179
180
181
182
     */
    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
    {
        PrimaryVariables values(0.0);
        values = initialAtPos(globalPos);

183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
        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;
            }

204
            if(globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_)
205
206
207
208
209
                values[Indices::conti0EqIdx + 1] = topMoleFraction;
        }
        else // horizontal flow
        {
            static const Scalar inletMoleFraction = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InletMoleFraction");
210
            if(globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_)
211
212
213
                values[Indices::conti0EqIdx + 1] = inletMoleFraction;
        }

214
215
216
217
218

        return values;
    }

    /*!
219
     * \brief Evaluates the boundary conditions for a Neumann control volume.
220
221
     *
     * \param element The element for which the Neumann boundary condition is set
222
     * \param fvGeometry The fvGeometry
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
     * \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))
        {
238
            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
239

240
            const auto tmp = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
            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
     */
    // \{

   /*!
259
     * \brief Evaluates the initial value for a control volume.
260
261
262
     *
     * \param globalPos The global position
     */
263
    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
264
265
    {
        PrimaryVariables values(0.0);
266
267
268
269
270
271
        values[Indices::pressureIdx] = 1e5;

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

        auto parabolicProfile = [&](const GlobalPosition& globalPos, int coord)
        {
272
273
274
275
            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]));
276
277
278
279
280
281
        };

        if (verticalFlow_)
            values[Indices::velocityYIdx] = parabolicProfile(globalPos, 0);
        else // horizontal flow
            values[Indices::velocityXIdx] = parabolicProfile(globalPos, 1);
282
283
284
285
286

        return values;
    }

    /*!
287
288
     * \brief Returns the intrinsic permeability of required as input parameter
     *        for the Beavers-Joseph-Saffman boundary condition.
289
     */
290
    Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
291
    {
292
        return couplingManager().couplingData().darcyPermeability(element, scvf);
293
294
295
    }

    /*!
296
297
     * \brief Returns the alpha value required as input parameter for the
     *        Beavers-Joseph-Saffman boundary condition.
298
299
300
301
302
303
     */
    Scalar alphaBJ(const SubControlVolumeFace& scvf) const
    {
        return 1.0;
    }

304
    /*!
305
     * \brief Sets the time loop pointer.
306
307
308
     */
    void setTimeLoop(TimeLoopPtr timeLoop)
    { timeLoop_ = timeLoop; }
309

310
    /*!
311
     * \brief Returns the time.
312
313
314
     */
    Scalar time() const
    { return timeLoop_->time(); }
315
316
317
318
319

    // \}

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

    bool onRightBoundary_(const GlobalPosition &globalPos) const
323
    { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
324
325

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

    bool onUpperBoundary_(const GlobalPosition &globalPos) const
329
    { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
330
331

    Scalar eps_;
332
    bool verticalFlow_;
333
    std::string problemName_;
334
    std::shared_ptr<CouplingManager> couplingManager_;
335
    TimeLoopPtr timeLoop_;
336
};
337
} // end namespace Dumux
338

339
#endif