localresidual.hh 19.6 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
22
 * \ingroup NavierStokesModel
 * \copydoc Dumux::NavierStokesResidualImpl
23
24
25
26
 */
#ifndef DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH
#define DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH

27
28
#include <dune/common/hybridutilities.hh>

29
#include <dumux/common/properties.hh>
30
#include <dumux/discretization/method.hh>
31
#include <dumux/discretization/extrusion.hh>
32
#include <dumux/assembly/staggeredlocalresidual.hh>
33
#include <dumux/freeflow/nonisothermal/localresidual.hh>
34

35
namespace Dumux {
36

37
38
39
40
41
42
43
44
namespace Impl {
template<class T>
static constexpr bool isRotationalExtrusion = false;

template<int radialAxis>
static constexpr bool isRotationalExtrusion<RotationalExtrusion<radialAxis>> = true;
} // end namespace Impl

45
46
47
48
/*!
 * \ingroup NavierStokesModel
 * \brief Element-wise calculation of the Navier-Stokes residual for models using the staggered discretization
 */
49
50
51
52
53

// forward declaration
template<class TypeTag, DiscretizationMethod discMethod>
class NavierStokesResidualImpl;

54
template<class TypeTag>
55
class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
56
: public StaggeredLocalResidual<TypeTag>
57
58
59
{
    using ParentType = StaggeredLocalResidual<TypeTag>;
    friend class StaggeredLocalResidual<TypeTag>;
60

61
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
62
63
64
65
66
67
68
69
70
71

    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
    using VolumeVariables = typename GridVolumeVariables::VolumeVariables;

    using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
    using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;

    using GridFaceVariables = typename GridVariables::GridFaceVariables;
    using ElementFaceVariables = typename GridFaceVariables::LocalView;
72

73
74
75
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
76
77
78
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
    using GridView = typename GridGeometry::GridView;
79
    using Element = typename GridView::template Codim<0>::Entity;
80
81
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
82
    using Extrusion = Extrusion_t<GridGeometry>;
83
    using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
84
85
86
87
    using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
    using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
    using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
88

89
90
    static constexpr bool normalizePressure = getPropValue<TypeTag, Properties::NormalizePressure>();

91
92
    using CellCenterResidual = CellCenterPrimaryVariables;
    using FaceResidual = FacePrimaryVariables;
93

94
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
95

96
public:
97
    using EnergyLocalResidual = FreeFlowEnergyLocalResidual<GridGeometry, FluxVariables, ModelTraits::enableEnergyBalance(), (ModelTraits::numFluidComponents() > 1)>;
98

99
100
101
102
    // account for the offset of the cell center privars within the PrimaryVariables container
    static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;
    static_assert(cellCenterOffset == ModelTraits::dim(), "cellCenterOffset must equal dim for staggered NavierStokes");

103
    //! Use the parent type's constructor
104
105
    using ParentType::ParentType;

106
    //! Evaluate fluxes entering or leaving the cell center control volume.
107
108
109
110
111
112
113
114
115
    CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
                                                        const Element &element,
                                                        const FVElementGeometry& fvGeometry,
                                                        const ElementVolumeVariables& elemVolVars,
                                                        const ElementFaceVariables& elemFaceVars,
                                                        const SubControlVolumeFace &scvf,
                                                        const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        FluxVariables fluxVars;
116
117
        CellCenterPrimaryVariables flux = fluxVars.computeMassFlux(problem, element, fvGeometry, elemVolVars,
                                                                   elemFaceVars, scvf, elemFluxVarsCache[scvf]);
118

119
        EnergyLocalResidual::heatFlux(flux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
120
121
122
123

        return flux;
    }

124
    //! Evaluate the source term for the cell center control volume.
125
126
127
128
129
130
131
    CellCenterPrimaryVariables computeSourceForCellCenter(const Problem& problem,
                                                          const Element &element,
                                                          const FVElementGeometry& fvGeometry,
                                                          const ElementVolumeVariables& elemVolVars,
                                                          const ElementFaceVariables& elemFaceVars,
                                                          const SubControlVolume &scv) const
    {
132
133
134
        CellCenterPrimaryVariables result(0.0);

        // get the values from the problem
135
        const auto sourceValues = problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scv);
136
137

        // copy the respective cell center related values to the result
138
139
        for (int i = 0; i < result.size(); ++i)
            result[i] = sourceValues[i + cellCenterOffset];
140
141

        return result;
142
143
144
    }


145
    //! Evaluate the storage term for the cell center control volume.
146
147
148
149
150
    CellCenterPrimaryVariables computeStorageForCellCenter(const Problem& problem,
                                                           const SubControlVolume& scv,
                                                           const VolumeVariables& volVars) const
    {
        CellCenterPrimaryVariables storage;
151
        storage[Indices::conti0EqIdx - ModelTraits::dim()] = volVars.density();
152

153
        EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
154

155
156
157
        return storage;
    }

158
    //! Evaluate the storage term for the face control volume.
159
160
161
    FacePrimaryVariables computeStorageForFace(const Problem& problem,
                                               const SubControlVolumeFace& scvf,
                                               const VolumeVariables& volVars,
162
                                               const ElementFaceVariables& elemFaceVars) const
163
164
    {
        FacePrimaryVariables storage(0.0);
165
        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
166
167
168
169
        storage[0] = volVars.density() * velocity;
        return storage;
    }

170
    //! Evaluate the source term for the face control volume.
171
    FacePrimaryVariables computeSourceForFace(const Problem& problem,
172
173
                                              const Element& element,
                                              const FVElementGeometry& fvGeometry,
174
175
                                              const SubControlVolumeFace& scvf,
                                              const ElementVolumeVariables& elemVolVars,
176
                                              const ElementFaceVariables& elemFaceVars) const
177
178
    {
        FacePrimaryVariables source(0.0);
179
        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
180
        source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
181
        source += problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
182

183
184
185
        // Axisymmetric problems in 2D feature an extra source terms arising from the transformation to cylindrical coordinates.
        // See Ferziger/Peric: Computational methods for fluid dynamics chapter 8.
        // https://doi.org/10.1007/978-3-540-68228-8 (page 301)
186
        if constexpr (ModelTraits::dim() == 2 && Impl::isRotationalExtrusion<Extrusion>)
187
188
189
190
191
192
193
194
195
196
        {
            if (scvf.directionIndex() == Extrusion::radialAxis)
            {
                // Velocity term
                const auto r = scvf.center()[scvf.directionIndex()];
                source -= -2.0*insideVolVars.effectiveViscosity() * elemFaceVars[scvf].velocitySelf() / (r*r);

                // Pressure term (needed because we incorporate pressure in terms of a surface integral).
                // grad(p) becomes div(pI) + (p/r)*n_r in cylindrical coordinates. The second term
                // is new with respect to Cartesian coordinates and handled below as a source term.
197
198
199
200
201
                const Scalar pressure =
                    normalizePressure ? insideVolVars.pressure() - problem.initial(scvf)[Indices::pressureIdx]
                                      : insideVolVars.pressure();

                source += pressure/r;
202
203
204
            }
        }

205
206
207
        return source;
    }

208
    //! Evaluate the momentum flux for the face control volume.
209
210
211
212
213
    FacePrimaryVariables computeFluxForFace(const Problem& problem,
                                            const Element& element,
                                            const SubControlVolumeFace& scvf,
                                            const FVElementGeometry& fvGeometry,
                                            const ElementVolumeVariables& elemVolVars,
214
                                            const ElementFaceVariables& elemFaceVars,
215
216
217
                                            const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        FluxVariables fluxVars;
218
        return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
219
220
    }

221
    /*!
222
223
     * \brief Evaluate boundary conditions for a cell center dof
     */
224
225
226
227
228
229
230
231
    CellCenterResidual computeBoundaryFluxForCellCenter(const Problem& problem,
                                                        const Element& element,
                                                        const FVElementGeometry& fvGeometry,
                                                        const SubControlVolumeFace& scvf,
                                                        const ElementVolumeVariables& elemVolVars,
                                                        const ElementFaceVariables& elemFaceVars,
                                                        const ElementBoundaryTypes& elemBcTypes,
                                                        const ElementFluxVariablesCache& elemFluxVarsCache) const
232
    {
233
        CellCenterResidual result(0.0);
234

235
236
237
        if (scvf.boundary())
        {
            const auto bcTypes = problem.boundaryTypes(element, scvf);
238

239
            // no fluxes occur over symmetry boundaries
240
241
242
243
244
            if (bcTypes.isSymmetry())
                return result;

            // treat Dirichlet and outflow BCs
            result = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
245

246
247
248
249
            // treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
            static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
            if (bcTypes.hasNeumann())
            {
250
251
252
253
                const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
                const auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf);

                for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
254
                {
255
                    if (bcTypes.isNeumann(eqIdx + cellCenterOffset))
256
                        result[eqIdx] = neumannFluxes[eqIdx + cellCenterOffset] * extrusionFactor * Extrusion::area(scvf);
257
                }
258
            }
259
260
261

            // account for wall functions, if used
            incorporateWallFunction_(result, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars);
262
        }
263
        return result;
264
265
    }

266
    /*!
267
     * \brief Evaluate Dirichlet (fixed value) boundary conditions for a face dof
268
     */
269
270
271
272
273
274
275
276
277
    void evalDirichletBoundariesForFace(FaceResidual& residual,
                                        const Problem& problem,
                                        const Element& element,
                                        const FVElementGeometry& fvGeometry,
                                        const SubControlVolumeFace& scvf,
                                        const ElementVolumeVariables& elemVolVars,
                                        const ElementFaceVariables& elemFaceVars,
                                        const ElementBoundaryTypes& elemBcTypes,
                                        const ElementFluxVariablesCache& elemFluxVarsCache) const
278
279
280
281
282
283
    {
        if (scvf.boundary())
        {
            // handle the actual boundary conditions:
            const auto bcTypes = problem.boundaryTypes(element, scvf);

284
            if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
285
            {
286
                // set a fixed value for the velocity for Dirichlet boundary conditions
287
                const Scalar velocity = elemFaceVars[scvf].velocitySelf();
288
                const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
289
290
                residual = velocity - dirichletValue;
            }
291
            else if(bcTypes.isSymmetry())
292
            {
293
                // for symmetry boundary conditions, there is no flow accross the boundary and
294
                // we therefore treat it like a Dirichlet boundary conditions with zero velocity
295
                const Scalar velocity = elemFaceVars[scvf].velocitySelf();
296
297
298
                const Scalar fixedValue = 0.0;
                residual = velocity - fixedValue;
            }
299
300
301
302
        }
    }

    /*!
303
     * \brief Evaluate boundary fluxes for a face dof
304
305
306
307
308
309
310
311
312
313
314
315
316
317
     */
    FaceResidual computeBoundaryFluxForFace(const Problem& problem,
                                            const Element& element,
                                            const FVElementGeometry& fvGeometry,
                                            const SubControlVolumeFace& scvf,
                                            const ElementVolumeVariables& elemVolVars,
                                            const ElementFaceVariables& elemFaceVars,
                                            const ElementBoundaryTypes& elemBcTypes,
                                            const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        FaceResidual result(0.0);

        if (scvf.boundary())
        {
318
319
            FluxVariables fluxVars;

320
321
322
323
324
325
326
327
            // handle the actual boundary conditions:
            const auto bcTypes = problem.boundaryTypes(element, scvf);
            if (bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
            {
                // the source term has already been accounted for, here we
                // add a given Neumann flux for the face on the boundary itself ...
                const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
                result = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
328
                                         * extrusionFactor * Extrusion::area(scvf);
329
330

                // ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume
331
                result += fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
332
            }
333
            else if(bcTypes.isDirichlet(Indices::pressureIdx))
334
            {
335
                // we are at an "fixed pressure" boundary for which the resdiual of the momentum balance needs to be assembled
336
                // as if it where inside the domain and not on the boundary (source term has already been acounted for)
337
338
339
340
                result = fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());

                // incorporate the inflow or outflow contribution
                result += fluxVars.inflowOutflowBoundaryFlux(problem, element, scvf, elemVolVars, elemFaceVars);
341
342
            }
        }
343
        return result;
344
345
346
347
    }

private:

348
349
350
    //! do nothing if no turbulence model is used
    template<class ...Args, bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<!turbulenceModel, int> = 0>
    void incorporateWallFunction_(Args&&... args) const
351
    {}
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372

    //! if a turbulence model is used, ask the problem is a wall function shall be employed and get the flux accordingly
    template<bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<turbulenceModel, int> = 0>
    void incorporateWallFunction_(CellCenterResidual& boundaryFlux,
                                  const Problem& problem,
                                  const Element& element,
                                  const FVElementGeometry& fvGeometry,
                                  const SubControlVolumeFace& scvf,
                                  const ElementVolumeVariables& elemVolVars,
                                  const ElementFaceVariables& elemFaceVars) const
    {
        static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
        const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();

        // account for wall functions, if used
        for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
        {
            // use a wall function
            if(problem.useWallFunction(element, scvf, eqIdx + cellCenterOffset))
            {
                boundaryFlux[eqIdx] = problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx]
373
                                                           * extrusionFactor * Extrusion::area(scvf);
374
375
            }
        }
376
    }
377

378
379
380
381
382
383
384
385
    //! Returns the implementation of the problem (i.e. static polymorphism)
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }

    //! \copydoc asImp_()
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }
};
386
387

} // end namespace Dumux
388

389
#endif