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
29
30
31
32
33
34
35
36
37
38
39
40
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
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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::DiscreteFractureExtensiveQuantities
*/
#ifndef EWOMS_DISCRETE_FRACTURE_EXTENSIVE_QUANTITIES_HH
#define EWOMS_DISCRETE_FRACTURE_EXTENSIVE_QUANTITIES_HH
#include <opm/models/immiscible/immiscibleextensivequantities.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
* \ingroup ExtensiveQuantities
*
* \brief This class expresses all intensive quantities of the discrete fracture model.
*/
template <class TypeTag>
class DiscreteFractureExtensiveQuantities : public ImmiscibleExtensiveQuantities<TypeTag>
{
using ParentType = ImmiscibleExtensiveQuantities<TypeTag>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
enum { dimWorld = GridView::dimensionworld };
enum { numPhases = FluidSystem::numPhases };
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using DimVector = Dune::FieldVector<Scalar, dimWorld>;
public:
/*!
* \copydoc MultiPhaseBaseExtensiveQuantities::update()
*/
void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{
ParentType::update(elemCtx, scvfIdx, timeIdx);
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.interiorFace(scvfIdx);
unsigned insideScvIdx = scvf.interiorIndex();
unsigned outsideScvIdx = scvf.exteriorIndex();
unsigned globalI = elemCtx.globalSpaceIndex(insideScvIdx, timeIdx);
unsigned globalJ = elemCtx.globalSpaceIndex(outsideScvIdx, timeIdx);
const auto& fractureMapper = elemCtx.problem().fractureMapper();
if (!fractureMapper.isFractureEdge(globalI, globalJ))
// do nothing if no fracture goes though the current edge
return;
// average the intrinsic permeability of the fracture
elemCtx.problem().fractureFaceIntrinsicPermeability(fractureIntrinsicPermeability_,
elemCtx, scvfIdx, timeIdx);
auto distDirection = elemCtx.pos(outsideScvIdx, timeIdx);
distDirection -= elemCtx.pos(insideScvIdx, timeIdx);
distDirection /= distDirection.two_norm();
const auto& problem = elemCtx.problem();
fractureWidth_ = problem.fractureWidth(elemCtx, insideScvIdx,
outsideScvIdx, timeIdx);
assert(fractureWidth_ < scvf.area());
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto& pGrad = extQuants.potentialGrad(phaseIdx);
unsigned upstreamIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
// multiply with the fracture mobility of the upstream vertex
fractureIntrinsicPermeability_.mv(pGrad,
fractureFilterVelocity_[phaseIdx]);
fractureFilterVelocity_[phaseIdx] *= -up.fractureMobility(phaseIdx);
// divide the volume flux by two. This is required because
// a fracture is always shared by two sub-control-volume
// faces.
fractureVolumeFlux_[phaseIdx] = 0;
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
fractureVolumeFlux_[phaseIdx] +=
(fractureFilterVelocity_[phaseIdx][dimIdx] * distDirection[dimIdx])
* (fractureWidth_ / 2.0) / scvf.area();
}
}
public:
const DimMatrix& fractureIntrinsicPermeability() const
{ return fractureIntrinsicPermeability_; }
Scalar fractureVolumeFlux(unsigned phaseIdx) const
{ return fractureVolumeFlux_[phaseIdx]; }
Scalar fractureWidth() const
{ return fractureWidth_; }
const DimVector& fractureFilterVelocity(unsigned phaseIdx) const
{ return fractureFilterVelocity_[phaseIdx]; }
private:
DimMatrix fractureIntrinsicPermeability_;
DimVector fractureFilterVelocity_[numPhases];
Scalar fractureVolumeFlux_[numPhases];
Scalar fractureWidth_;
};
} // namespace Opm
#endif
|
