co2model.hh 11.3 KB
Newer Older
1
2
3
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
Bernd Flemisch's avatar
Bernd Flemisch committed
4
 *   See the file COPYING for full copying permissions.                      *
5
6
7
8
9
10
11
12
 *                                                                           *
 *   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          *
Bernd Flemisch's avatar
Bernd Flemisch committed
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
18
19
20
21
 *   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
 *
22
 * \brief Adaption of the fully implicit scheme to the CO2Model model.
23
24
25
26
 */
#ifndef DUMUX_CO2_MODEL_HH
#define DUMUX_CO2_MODEL_HH

27
#include <dumux/implicit/2p2c/2p2cmodel.hh>
28
29
30
31
32

namespace Dumux
{
/*!
 * \ingroup CO2Model
33
 * \brief Adaption of the BOX or CC scheme to the non-isothermal two-phase two-component flow model.
Thomas Fetzer's avatar
Thomas Fetzer committed
34
 *
35
36
37
38
39
40
41
42
43
 *   See TwoPTwoCModel for reference to the equations used.
 *   The CO2 model is derived from the 2p2c model. In the CO2 model the phase switch criterion
 *   is different from the 2p2c model.
 *   The phase switch occurs when the equilibrium concentration
 *   of a component in a phase is exceeded, instead of the sum of the components in the virtual phase
 *   (the phase which is not present) being greater that unity as done in the 2p2c model.
 *   The CO2VolumeVariables do not use a constraint solver for calculating the mole fractions as is the
 *   case in the 2p2c model. Instead mole fractions are calculated in the FluidSystem with a given
 *   temperature, pressurem and salinity.
44
45
 *   The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
 * 	 problem file. Make sure that the according units are used in the problem setup. useMoles is set to false by default.
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
 *
 */

template<class TypeTag>
class CO2Model: public TwoPTwoCModel<TypeTag>
{
    typedef TwoPTwoCModel<TypeTag> ParentType;

     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;

     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
     enum {
         switchIdx = Indices::switchIdx,

         wPhaseIdx = Indices::wPhaseIdx,
         nPhaseIdx = Indices::nPhaseIdx,
         wCompIdx = Indices::wCompIdx,
         nCompIdx = Indices::nCompIdx,

         wPhaseOnly = Indices::wPhaseOnly,
         nPhaseOnly = Indices::nPhaseOnly,
         bothPhases = Indices::bothPhases,

72
73
         pwsn = TwoPTwoCFormulation::pwsn,
         pnsw = TwoPTwoCFormulation::pnsw,
74
75
76
77
78
79
80
81
82
83
84
85
         formulation = GET_PROP_VALUE(TypeTag, Formulation)
     };

     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
     typedef typename GridView::template Codim<0>::Iterator ElementIterator;
     enum {
         dim = GridView::dimension,
         dimWorld = GridView::dimensionworld
     };

     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
86
     static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
87
88
     enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
     enum { dofCodim = isBox ? dim : 0 };
89

90
public:
91
92
93
94
95
96
97
98
99
100
101
102
103


     /*!
      * \brief Update the static data of all vertices in the grid.
      *
      * \param curGlobalSol The current global solution
      * \param oldGlobalSol The previous global solution
      */
     void updateStaticData(SolutionVector &curGlobalSol,
                           const SolutionVector &oldGlobalSol)
     {
         bool wasSwitched = false;

104
105
         for (unsigned i = 0; i < ParentType::staticDat_.size(); ++i)
             ParentType::staticDat_[i].visited = false;
106
107
108

         FVElementGeometry fvGeometry;
         static VolumeVariables volVars;
109
110
111
         ElementIterator eIt = this->gridView_().template begin<0> ();
         const ElementIterator &eEndIt = this->gridView_().template end<0> ();
         for (; eIt != eEndIt; ++eIt)
112
         {
113
             fvGeometry.update(this->gridView_(), *eIt);
114
             for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
115
             {
116
                 int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
117

118
                 if (ParentType::staticDat_[globalIdx].visited)
119
120
                     continue;

121
                 ParentType::staticDat_[globalIdx].visited = true;
122
123
                 volVars.update(curGlobalSol[globalIdx],
                                this->problem_(),
124
                                *eIt,
125
126
127
                                fvGeometry,
                                scvIdx,
                                false);
128
                 const GlobalPosition &globalPos = eIt->geometry().corner(scvIdx);
129
130
131
132
133
                 if (primaryVarSwitch_(curGlobalSol,
                                       volVars,
                                       globalIdx,
                                       globalPos))
                 {
134
                     this->jacobianAssembler().markDofRed(globalIdx);
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
                     wasSwitched = true;
                 }
             }
         }

         // make sure that if there was a variable switch in an
         // other partition we will also set the switch flag
         // for our partition.
         if (this->gridView_().comm().size() > 1)
             wasSwitched = this->gridView_().comm().max(wasSwitched);

         ParentType::setSwitched_(wasSwitched);
     }

 protected:


     /*!
153
154
      * \brief Performs variable switch at a vertex, returns true if a
      *        variable switch was performed.
155
156
157
158
159
160
161
162
      */
     bool primaryVarSwitch_(SolutionVector &globalSol,
                              const VolumeVariables &volVars, int globalIdx,
                              const GlobalPosition &globalPos)
       {
         typename FluidSystem::ParameterCache paramCache;
           // evaluate primary variable switch
           bool wouldSwitch = false;
163
           int phasePresence = ParentType::staticDat_[globalIdx].phasePresence;
164
165
166
167
168
169
           int newPhasePresence = phasePresence;

           // check if a primary var switch is necessary
           if (phasePresence == nPhaseOnly)
           {

170
               Scalar xnw = volVars.moleFraction(nPhaseIdx, wCompIdx);
171
172
173
174
175
               Scalar xnwMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, nPhaseIdx);

               if(xnw > xnwMax)
                   wouldSwitch = true;

176
               if (ParentType::staticDat_[globalIdx].wasSwitched)
177
178
179
180
181
182
183
184
185
186
                   xnwMax *= 1.02;

               //If mole fraction is higher than the equilibrium mole fraction make a phase switch
               if(xnw > xnwMax)
               {
                   // wetting phase appears
                   std::cout << "wetting phase appears at vertex " << globalIdx
                             << ", coordinates: " << globalPos << ", xnw > xnwMax: "
                             << xnw << " > "<< xnwMax << std::endl;
                   newPhasePresence = bothPhases;
187
                   if (formulation == pnsw)
188
                       globalSol[globalIdx][switchIdx] = 0.0;
189
                   else if (formulation == pwsn)
190
191
192
193
194
195
                       globalSol[globalIdx][switchIdx] = 1.0;
               }
           }
           else if (phasePresence == wPhaseOnly)
           {

196
               Scalar xwn = volVars.moleFraction(wPhaseIdx, nCompIdx);
197
198
199
200
201
               Scalar xwnMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, wPhaseIdx);

               //If mole fraction is higher than the equilibrium mole fraction make a phase switch
               if(xwn > xwnMax)
                   wouldSwitch = true;
202
               if (ParentType::staticDat_[globalIdx].wasSwitched)
203
204
205
206
207
208
209
210
211
212
213
                   xwnMax *= 1.02;


               if(xwn > xwnMax)
               {
                   // non-wetting phase appears
                   std::cout << "non-wetting phase appears at vertex " << globalIdx
                             << ", coordinates: " << globalPos << ", xwn > xwnMax: "
                             << xwn << " > "<< xwnMax << std::endl;

                   newPhasePresence = bothPhases;
214
                   if (formulation == pnsw)
215
                       globalSol[globalIdx][switchIdx] = 0.999;
216
                   else if (formulation == pwsn)
217
218
219
220
221
222
                       globalSol[globalIdx][switchIdx] = 0.001;
               }
           }
           else if (phasePresence == bothPhases)
           {
               Scalar Smin = 0.0;
223
               if (ParentType::staticDat_[globalIdx].wasSwitched)
224
225
226
227
228
229
230
                   Smin = -0.01;

               if (volVars.saturation(nPhaseIdx) <= Smin)
               {
                   wouldSwitch = true;
                   // nonwetting phase disappears
                   std::cout << "Nonwetting phase disappears at vertex " << globalIdx
231
                             << ", coordinates: " << globalPos << ", sn: "
232
233
234
                             << volVars.saturation(nPhaseIdx) << std::endl;
                   newPhasePresence = wPhaseOnly;

235
236
237
                   if(!useMoles) //mass-fraction formulation
                   {
					   globalSol[globalIdx][switchIdx]
238
						   = volVars.massFraction(wPhaseIdx, nCompIdx);
239
240
241
242
                   }
                   else //mole-fraction formulation
                   {
					   globalSol[globalIdx][switchIdx]
243
					       = volVars.moleFraction(wPhaseIdx, nCompIdx);
244
                   }
245
246
247
248
249
250
               }
               else if (volVars.saturation(wPhaseIdx) <= Smin)
               {
                   wouldSwitch = true;
                   // wetting phase disappears
                   std::cout << "Wetting phase disappears at vertex " << globalIdx
251
                             << ", coordinates: " << globalPos << ", sw: "
252
253
254
                             << volVars.saturation(wPhaseIdx) << std::endl;
                   newPhasePresence = nPhaseOnly;

255
256
257
                   if(!useMoles) //mass-fraction formulation
                   {
					   globalSol[globalIdx][switchIdx]
258
						   = volVars.massFraction(nPhaseIdx, wCompIdx);
259
260
261
262
                   }
                   else //mole-fraction formulation
                   {
						globalSol[globalIdx][switchIdx]
263
						= volVars.moleFraction(nPhaseIdx, wCompIdx);
264
                   }
265
266
267
               }
           }

268
269
           ParentType::staticDat_[globalIdx].phasePresence = newPhasePresence;
           ParentType::staticDat_[globalIdx].wasSwitched = wouldSwitch;
270
271
272
273
274
275
276
277
278
           return phasePresence != newPhasePresence;
       }


};

}

#endif