fluxvariables.hh 14.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- 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
21
22
 * \ingroup NavierStokesModel
 * \copydoc Dumux::NavierStokesFluxVariablesImpl
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 */
#ifndef DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
#define DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH

#include <dumux/common/properties.hh>
#include <dumux/discretization/fluxvariablesbase.hh>
#include <dumux/discretization/methods.hh>

namespace Dumux
{

// forward declaration
template<class TypeTag, DiscretizationMethods Method>
class NavierStokesFluxVariablesImpl;


/*!
40
41
 * \ingroup NavierStokesModel
 * \brief The flux variables class for the Navier-Stokes model using the staggered grid discretization.
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
 */
template<class TypeTag>
class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
: public FluxVariablesBase<TypeTag>
{
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using Element = typename GridView::template Codim<0>::Entity;
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
    using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
    using IndexType = typename GridView::IndexSet::IndexType;
    using Stencil = std::vector<IndexType>;

    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
    using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);

    static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);

    enum {
         // grid and world dimension
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld,

        massBalanceIdx = Indices::massBalanceIdx,
    };

    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;

    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
    typename DofTypeIndices::FaceIdx faceIdx;

public:

82
83
84
85
86
87
88
89
    /*!
    * \brief Returns the advective flux over a sub control volume face
    * \param elemVolVars All volume variables for the element
    * \param elemFaceVars The face variables
    * \param scvf The sub control volume face
    * \param upwindTerm The uwind term (i.e. the advectively transported quantity)
    * \param isOutflow Determines if we are on an outflow boundary
    */
90
91
92
93
94
95
96
97
    template<class UpwindTerm>
    static Scalar advectiveFluxForCellCenter(const ElementVolumeVariables& elemVolVars,
                                             const ElementFaceVariables& elemFaceVars,
                                             const SubControlVolumeFace &scvf,
                                             UpwindTerm upwindTerm,
                                             bool isOutflow = false)
    {
        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
98
        const bool insideIsUpstream = scvf.directionSign() == sign(velocity);
99
100
101
102
103
104
105
106
107
108
        static const Scalar upWindWeight = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");

        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
        const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];

        const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
        const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars;

        const Scalar flux = (upWindWeight * upwindTerm(upstreamVolVars) +
                            (1.0 - upWindWeight) * upwindTerm(downstreamVolVars))
109
                            * velocity * scvf.area() * scvf.directionSign();
110
111
112
113

        return flux;
    }

114
115
116
    /*!
    * \brief Computes the flux for the cell center residual.
    */
117
118
119
120
121
122
123
124
    CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
                                                        const Element &element,
                                                        const FVElementGeometry& fvGeometry,
                                                        const ElementVolumeVariables& elemVolVars,
                                                        const ElementFaceVariables& elemFaceVars,
                                                        const SubControlVolumeFace &scvf,
                                                        const FluxVariablesCache& fluxVarsCache)
    {
125
        auto upwindTerm = [](const auto& volVars) { return volVars.density(); };
126

127
        const Scalar flux = advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, false);
128

129
130
        CellCenterPrimaryVariables result(0.0);
        result[massBalanceIdx] = flux;
131

132
        return result;
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    }

    /*!
    * \brief Returns the normal part of the momentum flux
    */
   FacePrimaryVariables computeNormalMomentumFlux(const Problem& problem,
                                                  const Element& element,
                                                  const SubControlVolumeFace& scvf,
                                                  const FVElementGeometry& fvGeometry,
                                                  const ElementVolumeVariables& elemVolVars,
                                                  const ElementFaceVariables& elementFaceVars)
   {
       const auto insideScvIdx = scvf.insideScvIdx();
       const auto& insideVolVars = elemVolVars[insideScvIdx];
       const Scalar velocitySelf = elementFaceVars[scvf].velocitySelf() ;
       const Scalar velocityOpposite = elementFaceVars[scvf].velocityOpposite();
       FacePrimaryVariables normalFlux(0.0);

       if(navierStokes)
       {
           // advective part
           const Scalar vAvg = (velocitySelf + velocityOpposite) * 0.5;
155
           const Scalar vUp = (scvf.directionSign() == sign(vAvg)) ? velocityOpposite : velocitySelf;
156
157
158
159
160
161
162
163
164
165
166
167
           normalFlux += vAvg * vUp * insideVolVars.density();
       }

       // diffusive part
       const Scalar deltaV = scvf.normalInPosCoordDir() ?
                             (velocitySelf - velocityOpposite) :
                             (velocityOpposite - velocitySelf);

       const Scalar deltaX = scvf.selfToOppositeDistance();
       normalFlux -= insideVolVars.viscosity() * 2.0 * deltaV/deltaX;

       // account for the orientation of the face
168
       const Scalar sgn = -1.0 * scvf.directionSign();
169
170
171
172
173
174

       Scalar result = normalFlux * sgn * scvf.area();

       // treat outflow conditions
       if(navierStokes && scvf.boundary())
       {
175
           const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ?
176
177
                                   elemVolVars[insideScvIdx] : elemVolVars[scvf.outsideScvIdx()] ;

178
           result += velocitySelf * velocitySelf * upVolVars.density() * scvf.directionSign() * scvf.area() ;
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
228
229
230
231
232
233
234
235
236
237
238
239
240
       }
       return result;
   }

   /*!
   * \brief Returns the tangential part of the momentum flux
   */
  FacePrimaryVariables computeTangetialMomentumFlux(const Problem& problem,
                                                    const Element& element,
                                                    const SubControlVolumeFace& scvf,
                                                    const FVElementGeometry& fvGeometry,
                                                    const ElementVolumeVariables& elemVolVars,
                                                    const ElementFaceVariables& elementFaceVars)
  {
      FacePrimaryVariables tangentialFlux(0.0);
      auto& faceVars = elementFaceVars[scvf];
      const int numSubFaces = scvf.pairData().size();

      // account for all sub-faces
      for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
      {
          const auto eIdx = scvf.insideScvIdx();
          const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);

          // Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected.
          if(scvf.pairData()[localSubFaceIdx].outerParallelFaceDofIdx < 0)
          {
              // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
              auto makeGhostFace = [eIdx] (const GlobalPosition& pos)
              {
                  return SubControlVolumeFace(pos, std::vector<unsigned int>{eIdx,eIdx});
              };

              // use the ghost face to check if there is a symmetry boundary condition and skip any further steps if yes
              const auto bcTypes = problem.boundaryTypes(element, makeGhostFace(scvf.pairData()[localSubFaceIdx].virtualOuterParallelFaceDofPos));
              if(bcTypes.isSymmetry())
                continue;
          }

          // if there is no symmetry boundary condition, proceed to calculate the tangential momentum flux
          if(navierStokes)
              tangentialFlux += computeAdvectivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);

          tangentialFlux += computeDiffusivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
      }
      return tangentialFlux;
  }

private:

  FacePrimaryVariables computeAdvectivePartOfTangentialMomentumFlux_(const Problem& problem,
                                                                     const Element& element,
                                                                     const SubControlVolumeFace& scvf,
                                                                     const SubControlVolumeFace& normalFace,
                                                                     const ElementVolumeVariables& elemVolVars,
                                                                     const FaceVariables& faceVars,
                                                                     const int localSubFaceIdx)
  {
      const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
      const auto insideScvIdx = normalFace.insideScvIdx();
      const auto outsideScvIdx = normalFace.outsideScvIdx();

241
      const bool innerElementIsUpstream = ( normalFace.directionSign() == sign(transportingVelocity) );
242
243
244
245
246
247
248
249
250

      const auto& upVolVars = innerElementIsUpstream ? elemVolVars[insideScvIdx] : elemVolVars[outsideScvIdx];

      const Scalar transportedVelocity = innerElementIsUpstream ?
                                         faceVars.velocitySelf() :
                                         faceVars.velocityParallel(localSubFaceIdx);

      const Scalar momentum = upVolVars.density() * transportedVelocity;

251
      return transportingVelocity * momentum * normalFace.directionSign() * normalFace.area() * 0.5;
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
  }

  FacePrimaryVariables computeDiffusivePartOfTangentialMomentumFlux_(const Problem& problem,
                                                                     const Element& element,
                                                                     const SubControlVolumeFace& scvf,
                                                                     const SubControlVolumeFace& normalFace,
                                                                     const ElementVolumeVariables& elemVolVars,
                                                                     const FaceVariables& faceVars,
                                                                     const int localSubFaceIdx)
  {
      FacePrimaryVariables tangentialDiffusiveFlux(0.0);

      const auto insideScvIdx = normalFace.insideScvIdx();
      const auto outsideScvIdx = normalFace.outsideScvIdx();

      const auto& insideVolVars = elemVolVars[insideScvIdx];
      const auto& outsideVolVars = elemVolVars[outsideScvIdx];

      // the averaged viscosity at the face normal to our face of interest (where we assemble the face residual)
      const Scalar muAvg = (insideVolVars.viscosity() + outsideVolVars.viscosity()) * 0.5;

      // the normal derivative
      const Scalar innerNormalVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
      const Scalar outerNormalVelocity = faceVars.velocityNormalOutside(localSubFaceIdx);

      const Scalar normalDeltaV = scvf.normalInPosCoordDir() ?
                                    (outerNormalVelocity - innerNormalVelocity) :
                                    (innerNormalVelocity - outerNormalVelocity);

      const Scalar normalDerivative = normalDeltaV / scvf.pairData(localSubFaceIdx).normalDistance;
      tangentialDiffusiveFlux -= muAvg * normalDerivative;

      // the parallel derivative
      const Scalar innerParallelVelocity = faceVars.velocitySelf();

      const Scalar outerParallelVelocity = faceVars.velocityParallel(localSubFaceIdx);

      const Scalar parallelDeltaV = normalFace.normalInPosCoordDir() ?
                                   (outerParallelVelocity - innerParallelVelocity) :
                                   (innerParallelVelocity - outerParallelVelocity);

      const Scalar parallelDerivative = parallelDeltaV / scvf.pairData(localSubFaceIdx).parallelDistance;
      tangentialDiffusiveFlux -= muAvg * parallelDerivative;

296
      return tangentialDiffusiveFlux * normalFace.directionSign() * normalFace.area() * 0.5;
297
298
299
300
301
302
  }
};

} // end namespace

#endif