problem.hh 8.86 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
21
22
23
24
25
26
27
28
29
 *   (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 SolidEnergyTests
 * \brief A piece of stone heated by a network of channels
 */
#ifndef DUMUX_TEST_SOLIDENERGY_PROBLEM_HH
#define DUMUX_TEST_SOLIDENERGY_PROBLEM_HH

#include <string>
#include <cmath>

30
#include <dumux/common/boundarytypes.hh>
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/problem.hh>

namespace Dumux {

/*!
 * \ingroup SolidEnergyTests
 * \brief A piece of stone heated by a network of channels
 */
template <class TypeTag>
class SolidEnergyProblem : public PorousMediumFlowProblem<TypeTag>
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;
45
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
46
47
48
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using NeumannValues = GetPropType<TypeTag, Properties::NumEqVector>;
49
    using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
50
51
52
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
    using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
53
54
55

    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
56
57
58
59
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
    using SubControlVolume = typename GridGeometry::SubControlVolume;
    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
60
61
62
    using PointSource = GetPropType<TypeTag, Properties::PointSource>;

public:
63
64
    SolidEnergyProblem(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
    : ParentType(gridGeometry, paramGroup)
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
    {
        name_ = getParam<std::string>("Problem.Name");
        temperatureHigh_ = getParam<Scalar>("Problem.TemperatureHigh");
        temperatureInit_ = getParam<Scalar>("Problem.TemperatureInit");
        lambdaSolid_ = getParam<Scalar>("Component.SolidThermalConductivity");
    }

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

    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     */
    const std::string& name() const
    {
        return name_;
    }
    // \}

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

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
     * \param globalPos The position for which the bc type should be evaluated
     */
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
        BoundaryTypes bcTypes;
        bcTypes.setAllNeumann();
        return bcTypes;
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
     * This is the method for the case where the Neumann condition is
     * potentially solution dependent
     *
     * \param element The finite element
     * \param fvGeometry The finite-volume geometry
     * \param elemVolVars All volume variables for the element
116
     * \param elemFluxVarsCache Flux variables caches for all faces in stencil
117
118
119
120
121
122
123
124
     * \param scvf The sub control volume face
     *
     * Negative values mean influx.
     * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
     */
    NeumannValues neumann(const Element& element,
                          const FVElementGeometry& fvGeometry,
                          const ElementVolumeVariables& elemVolVars,
125
                          const ElementFluxVariablesCache& elemFluxVarsCache,
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
158
159
160
161
162
163
164
165
166
167
168
169
170
                          const SubControlVolumeFace& scvf) const
    {
        const auto& volVars = elemVolVars[scvf.insideScvIdx()];
        const Scalar value = 0.02*(volVars.temperatureSolid() - temperatureInit_)/0.002;
        return NeumannValues(value);
    }



    // \}

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

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
     * \param globalPos The position for which the initial condition should be evaluated
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
    {
        return PrimaryVariables(temperatureInit_);
    }

    /*!
     * \brief Applies a vector of point sources. The point sources
     *        are possibly solution dependent.
     *
     * \param pointSources A vector of PointSource s that contain
              source values for all phases and space positions.
     *
     * For this method, the values method of the point source
     * has to return the absolute rate values in units
     * \f$ [ \textnormal{unit of conserved quantity} / s ] \f$.
     * Positive values mean that the conserved quantity is created, negative ones mean that it vanishes.
     * E.g. for the mass balance that would be a mass rate in \f$ [ kg / s ] \f$.
     */
    void addPointSources(std::vector<PointSource>& pointSources) const
    {
        // add point sources along a sinus
171
        const auto& gg = this->gridGeometry();
172
173
174
175
176
177
178
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
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
        const auto ampl = (gg.bBoxMax()[1] - gg.bBoxMin()[1])/3.0;
        const auto len = (gg.bBoxMax()[0] - gg.bBoxMin()[0]);
        const auto freq = 1.0/len;

        for (int i = 0; i < 100; ++i)
        {
            const auto posx = double(i)/99.0*len;
            const auto posy = ampl*std::sin(2*M_PI*freq*posx) + ampl*3.0/2.0;
            pointSources.emplace_back(GlobalPosition({posx, posy}), typename PointSource::Values({len/100.0}));
        }
    }

    /*!
     * \brief Evaluate the point sources (added by addPointSources)
     *        for all phases within a given sub-control-volume.
     *
     * This is the method for the case where the point source is
     * solution dependent
     *
     * \param source A single point source
     * \param element The finite element
     * \param fvGeometry The finite-volume geometry
     * \param elemVolVars All volume variables for the element
     * \param scv The sub control volume
     *
     * For this method, the values() method of the point sources returns
     * the absolute conserved quantity rate generated or annihilate in
     * units \f$ [ \textnormal{unit of conserved quantity} / s ] \f$.
     * Positive values mean that the conserved quantity is created, negative ones mean that it vanishes.
     * E.g. for the mass balance that would be a mass rate in \f$ [ kg / s ] \f$.
     */
    void pointSource(PointSource& source,
                     const Element &element,
                     const FVElementGeometry& fvGeometry,
                     const ElementVolumeVariables& elemVolVars,
                     const SubControlVolume &scv) const
    {
        const auto& volVars = elemVolVars[scv];
        const Scalar charLen = 0.1;
        static const Scalar nusselt = getParam<Scalar>("Problem.NusseltNumber");
        const Scalar interfacialArea = 2*M_PI*0.001;
        const Scalar lambdaSolidFluidTransfer = lambdaSolid_;
        source *= -nusselt*interfacialArea*lambdaSolidFluidTransfer*(volVars.temperatureSolid() - temperatureHigh_)/charLen;
    }

    // \}

private:
    Scalar temperatureHigh_, temperatureInit_, lambdaSolid_;
    static constexpr Scalar eps_ = 1e-6;
    std::string name_;
};

} //end namespace Dumux

#endif