zeroeq2cnisubproblem.hh 14.4 KB
Newer Older
Thomas Fetzer's avatar
Thomas Fetzer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// -*- 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    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (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
 * \brief Non-isothermal two-component ZeroEq subproblem with air flowing
 *        from the left to the right and coupling at the bottom.
 */
#ifndef DUMUX_ZEROEQTWOCNI_SUBPROBLEM_HH
#define DUMUX_ZEROEQTWOCNI_SUBPROBLEM_HH

#include <dumux/freeflow/zeroeqncni/zeroeqncnimodel.hh>
#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
29
#include <dumux/multidomain/2cnistokes2p2cni/stokesncnicouplinglocalresidual.hh>
Thomas Fetzer's avatar
Thomas Fetzer committed
30
31
32
33
34
35
36
37
38
39

namespace Dumux
{

template <class TypeTag>
class ZeroEq2cniSubProblem;

namespace Properties
{
NEW_TYPE_TAG(ZeroEq2cniSubProblem,
40
             INHERITS_FROM(BoxZeroEqncni, SubDomain));
Thomas Fetzer's avatar
Thomas Fetzer committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

// Set the problem property
SET_TYPE_PROP(ZeroEq2cniSubProblem, Problem, Dumux::ZeroEq2cniSubProblem<TypeTag>);

// Use the StokesncniCouplingLocalResidual for the computation of the local residual in the ZeroEq domain
SET_TYPE_PROP(ZeroEq2cniSubProblem, LocalResidual, StokesncniCouplingLocalResidual<TypeTag>);

// Set the property for the material parameters by extracting it from the material law.
SET_TYPE_PROP(ZeroEq2cniSubProblem,
              MaterialLawParams,
              typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);

// Use the fluid system from the coupled problem
SET_TYPE_PROP(ZeroEq2cniSubProblem,
              FluidSystem,
              typename GET_PROP_TYPE(typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag), FluidSystem));

// Disable use of mole formulation
SET_BOOL_PROP(ZeroEq2cniSubProblem, UseMoles, false);

// Disable gravity
SET_BOOL_PROP(ZeroEq2cniSubProblem, ProblemEnableGravity, false);

// Enable Navier-Stokes
SET_BOOL_PROP(ZeroEq2cniSubProblem, EnableNavierStokes, true);

// Set the properties for variable inflow BC
NEW_PROP_TAG(FreeFlowSinusVelocityAmplitude);
NEW_PROP_TAG(FreeFlowSinusVelocityPeriod);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusVelocityAmplitude, 0.0);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusVelocityPeriod, 3600.0);
NEW_PROP_TAG(FreeFlowSinusPressureAmplitude);
NEW_PROP_TAG(FreeFlowSinusPressurePeriod);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusPressureAmplitude, 0.0);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusPressurePeriod, 3600.0);
NEW_PROP_TAG(FreeFlowSinusConcentrationAmplitude);
NEW_PROP_TAG(FreeFlowSinusConcentrationPeriod);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusConcentrationAmplitude, 0.0);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusConcentrationPeriod, 3600.0);
NEW_PROP_TAG(FreeFlowSinusTemperatureAmplitude);
NEW_PROP_TAG(FreeFlowSinusTemperaturePeriod);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusTemperatureAmplitude, 0.0);
SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusTemperaturePeriod, 3600.0);
}

/*!
 * \ingroup ImplicitTestProblems
Hao Wu's avatar
Hao Wu committed
88
89
90
91
92
93
94
95
 * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
 * \brief Non-isothermal two-component ZeroEq subproblem with air flowing
 *        from the left to the right and coupling at the bottom.
 *
 * The free-flow subdomain is sized 0.75m times 0.5m. Dry and hot air is flowing from left (Dirichlet)
 * to right (outflow), at the middle third of the bottom the coupling conditions
 * are applied to all balance equations. They handle the exchange to the porous-medium
 * subdomain.
Thomas Fetzer's avatar
Thomas Fetzer committed
96
 *
Hao Wu's avatar
Hao Wu committed
97
98
 * This subproblem uses the \ref ZeroEqncniModel. It is part of a multidomain model and
 * combined with the 2p2cnisubproblem for the porous-medium domain.
Thomas Fetzer's avatar
Thomas Fetzer committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
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
153
154
155
156
157
 */
template <class TypeTag>
class ZeroEq2cniSubProblem : public ZeroEqProblem<TypeTag>
{
    typedef ZeroEq2cniSubProblem<TypeTag> ThisType;
    typedef ZeroEqProblem<TypeTag> ParentType;

    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;

    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;

    enum {
        dim = GridView::dimension
    };
    enum { // equation indices
        massBalanceIdx = Indices::massBalanceIdx,
        momentumXIdx = Indices::momentumXIdx, // Index of the x-component of the momentum balance
        momentumYIdx = Indices::momentumYIdx, // Index of the y-component of the momentum balance
        momentumZIdx = Indices::momentumZIdx, // Index of the z-component of the momentum balance
        transportEqIdx = Indices::transportEqIdx, // Index of the transport equation (massfraction)
        energyEqIdx =    Indices::energyEqIdx // Index of the energy equation (temperature)
    };
    enum { // primary variable indices
        pressureIdx = Indices::pressureIdx,
        velocityXIdx = Indices::velocityXIdx,
        velocityYIdx = Indices::velocityYIdx,
        velocityZIdx = Indices::velocityZIdx,
        massOrMoleFracIdx = Indices::massOrMoleFracIdx,
        temperatureIdx = Indices::temperatureIdx
    };
    enum {
        transportCompIdx = Indices::transportCompIdx, // water component index
        phaseCompIdx = Indices::phaseCompIdx // air component index
    };

    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::ctype CoordScalar;

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;

    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef Dune::FieldVector<CoordScalar, dim> GlobalPosition;


public:
    /*!
     * \brief The sub-problem for the ZeroEq subdomain
     *
     * \param timeManager The TimeManager which is used by the simulation
     * \param gridView The simulation's idea about physical space
     */
    ZeroEq2cniSubProblem(TimeManager &timeManager, const GridView &gridView)
158
        : ParentType(timeManager, gridView)
Thomas Fetzer's avatar
Thomas Fetzer committed
159
    {
Thomas Fetzer's avatar
[io]    
Thomas Fetzer committed
160
161
162
163
        bBoxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, LowerLeftX);
        bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, UpperRightX);
        bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
        bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, UpperRightY);
Thomas Fetzer's avatar
Thomas Fetzer committed
164
165
166
167
168
169
170
171
172
        runUpDistanceX1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX1); // first part of the interface without coupling
        runUpDistanceX2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX2); // second part of the interface without coupling

        refVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefVelocity);
        refPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefPressure);
        refMassfrac_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefMassfrac);
        refTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
    }

173
    // functions have to be overwritten, otherwise they remain uninitialized
Hao Wu's avatar
Hao Wu committed
174
    //! \copydoc Dumux::ImplicitProblem::bBoxMin()
Thomas Fetzer's avatar
Thomas Fetzer committed
175
176
177
    const GlobalPosition &bBoxMin() const
    { return bBoxMin_; }

Hao Wu's avatar
Hao Wu committed
178
    //! \copydoc Dumux::ImplicitProblem::bBoxMax()
Thomas Fetzer's avatar
Thomas Fetzer committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    const GlobalPosition &bBoxMax() const
    { return bBoxMax_; }

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

    //! \copydoc Dumux::ImplicitProblem::name()
    const std::string &name() const
    { return GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Output, NameFF); }

    // \}

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

    //! \copydoc Dumux::ImplicitProblem::boundaryTypesAtPos()
    void boundaryTypesAtPos(BoundaryTypes &values,
                            const GlobalPosition &globalPos) const
    {
        values.setAllDirichlet();


        if (onUpperBoundary_(globalPos))
        {
            values.setNeumann(transportEqIdx);
208
            values.setDirichlet(temperatureIdx);
Thomas Fetzer's avatar
Thomas Fetzer committed
209
210
211
212
213
        }

        if (onLowerBoundary_(globalPos))
        {
            values.setNeumann(transportEqIdx);
214
            values.setDirichlet(temperatureIdx);
Thomas Fetzer's avatar
Thomas Fetzer committed
215
216
217
218

            if (globalPos[0] > runUpDistanceX1_ - eps_
                && globalPos[0] < runUpDistanceX2_ + eps_)
            {
219
                values.setAllCouplingDirichlet();
220
221
                values.setCouplingNeumann(momentumXIdx);
                values.setCouplingNeumann(momentumYIdx);
Thomas Fetzer's avatar
Thomas Fetzer committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
            }
        }

        if (onRightBoundary_(globalPos))
        {
            values.setAllOutflow();

            if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos)) // corner points
                values.setAllDirichlet();
        }

        if (onLeftBoundary_(globalPos))
        {
            values.setAllDirichlet();
        }

        // the mass balance has to be of type outflow
        // it does not get a coupling condition, since pn is a condition for stokes
        values.setOutflow(massBalanceIdx);

        if (onRightBoundary_(globalPos))
        {
            values.setAllOutflow();
245
            values.setDirichlet(pressureIdx);
Thomas Fetzer's avatar
Thomas Fetzer committed
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
        }
    }

    //! \copydoc Dumux::ImplicitProblem::dirichletAtPos()
    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
    {
        values = 0.0;

        values[velocityXIdx] = xVelocity_(globalPos);
        values[velocityYIdx] = 0.0;
        values[pressureIdx] = refPressure()
                              + 1.189 * this->gravity()[1] * (globalPos[1] - bBoxMin_[1]);
        values[massOrMoleFracIdx] = refMassfrac();
        values[temperatureIdx] = refTemperature();
    }

    //! \copydoc Dumux::ImplicitProblem::neumannAtPos()
    void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
    {
        values = 0.;
    }

    //! \copydoc Dumux::ImplicitProblem::sourceAtPos()
    void sourceAtPos(PrimaryVariables &values,
                     const GlobalPosition &globalPos) const
    {
272
273
        // The source term of the mass balance has to be chosen as
        // div (q_momentum) in the problem file
Thomas Fetzer's avatar
Thomas Fetzer committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
        values = 0.0;
    }

    //! \copydoc Dumux::ImplicitProblem::initialAtPos()
    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
    {
        initial_(values, globalPos);
    }

    // \}

    //! \brief Returns the velocity at the inflow.
    const Scalar refVelocity() const
    {
        return refVelocity_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusVelocityAmplitude),
                                         GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusVelocityPeriod));
    }

    //! \brief Returns the pressure at the inflow.
    const Scalar refPressure() const
    {
        return refPressure_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusPressureAmplitude),
                                         GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusPressurePeriod));
    }

    //! \brief Returns the mass fraction at the inflow.
    const Scalar refMassfrac() const
    {
        return refMassfrac_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusConcentrationAmplitude),
                                         GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusConcentrationPeriod));
    }

    //! \brief Returns the temperature at the inflow.
    const Scalar refTemperature() const
    {
        return refTemperature_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperatureAmplitude),
                                            GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperaturePeriod));
    }

private:
    // Internal method for the initial and Dirichlet conditions
    void initial_(PrimaryVariables &values,
                  const GlobalPosition &globalPos) const
    {
        values[velocityXIdx] = xVelocity_(globalPos);
        values[velocityYIdx] = 0.;

        values[pressureIdx] = refPressure() + 1.189 * this->gravity()[1] * (globalPos[1] - bBoxMin_[1]);
        values[massOrMoleFracIdx] = refMassfrac();
        values[temperatureIdx] = refTemperature();
    }

    // returns the inflow velocity profile
    const Scalar xVelocity_(const GlobalPosition &globalPos) const
    {
        if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
            return 0.0;
        return refVelocity();
    }

    // can be used for the variation of a boundary condition
    const Scalar variation_(const Scalar amplitude, const Scalar period) const
    { return sin(2*M_PI*this->timeManager().time()/period) * amplitude; }

    bool onLeftBoundary_(const GlobalPosition &globalPos) const
    { return globalPos[0] < bBoxMin_[0] + eps_; }

    bool onRightBoundary_(const GlobalPosition &globalPos) const
    { return globalPos[0] > bBoxMax_[0] - eps_; }

    bool onLowerBoundary_(const GlobalPosition &globalPos) const
    { return globalPos[1] < bBoxMin_[1] + eps_; }

    bool onUpperBoundary_(const GlobalPosition &globalPos) const
    { return globalPos[1] > bBoxMax_[1] - eps_; }

    static constexpr Scalar eps_ = 1e-8;
    GlobalPosition bBoxMin_;
    GlobalPosition bBoxMax_;

    Scalar refVelocity_;
    Scalar refPressure_;
    Scalar refMassfrac_;
    Scalar refTemperature_;

    Scalar runUpDistanceX1_;
    Scalar runUpDistanceX2_;
};
} //end namespace

#endif // DUMUX_ZEROEQTWOCNI_SUBPROBLEM_HH