co2model.hh 12.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


     /*!
      * \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;
103
104
         int succeeded;
         try {
105

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

109
110
111
112
113
             FVElementGeometry fvGeometry;
             static VolumeVariables volVars;
             ElementIterator eIt = this->gridView_().template begin<0> ();
             const ElementIterator &eEndIt = this->gridView_().template end<0> ();
             for (; eIt != eEndIt; ++eIt)
114
             {
115
116
117

                 fvGeometry.update(this->gridView_(), *eIt);
                 for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
118
                 {
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
                     int dofIdxGlobal = this->dofMapper().map(*eIt, scvIdx, dofCodim);

                     if (ParentType::staticDat_[dofIdxGlobal].visited)
                         continue;

                     ParentType::staticDat_[dofIdxGlobal].visited = true;
                     volVars.update(curGlobalSol[dofIdxGlobal],
                             this->problem_(),
                             *eIt,
                             fvGeometry,
                             scvIdx,
                             false);
                     const GlobalPosition &globalPos = eIt->geometry().corner(scvIdx);
                     if (primaryVarSwitch_(curGlobalSol,
                             volVars,
                             dofIdxGlobal,
                             globalPos))
                     {
                         this->jacobianAssembler().markDofRed(dofIdxGlobal);
                         wasSwitched = true;
                     }
140
141
                 }
             }
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
             succeeded = 1;
         }
         catch (Dumux::NumericalProblem &e)
         {
             std::cout << "\n"
                     << "Rank " << this->problem_().gridView().comm().rank()
                     << " caught an exception while updating the static data." << e.what()
                     << "\n";
             succeeded = 0;
         }
         //make sure that all processes succeeded. If not throw a NumericalProblem to decrease the time step size.
         if (this->gridView_().comm().size() > 1)
             succeeded = this->gridView_().comm().min(succeeded);

         if (!succeeded) {
             if(this->problem_().gridView().comm().rank() == 0)
                 DUNE_THROW(NumericalProblem,
                         "A process did not succeed in updating the static data.");
             return;
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
         }

         // 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:


     /*!
176
177
      * \brief Performs variable switch at a vertex, returns true if a
      *        variable switch was performed.
178
179
      */
     bool primaryVarSwitch_(SolutionVector &globalSol,
180
                              const VolumeVariables &volVars, int dofIdxGlobal,
181
182
183
184
185
                              const GlobalPosition &globalPos)
       {
         typename FluidSystem::ParameterCache paramCache;
           // evaluate primary variable switch
           bool wouldSwitch = false;
186
           int phasePresence = ParentType::staticDat_[dofIdxGlobal].phasePresence;
187
188
189
190
191
192
           int newPhasePresence = phasePresence;

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

193
               Scalar xnw = volVars.moleFraction(nPhaseIdx, wCompIdx);
194
195
196
197
198
               Scalar xnwMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, nPhaseIdx);

               if(xnw > xnwMax)
                   wouldSwitch = true;

199
               if (ParentType::staticDat_[dofIdxGlobal].wasSwitched)
200
201
202
203
204
205
                   xnwMax *= 1.02;

               //If mole fraction is higher than the equilibrium mole fraction make a phase switch
               if(xnw > xnwMax)
               {
                   // wetting phase appears
206
                   std::cout << "wetting phase appears at vertex " << dofIdxGlobal
207
208
209
                             << ", coordinates: " << globalPos << ", xnw > xnwMax: "
                             << xnw << " > "<< xnwMax << std::endl;
                   newPhasePresence = bothPhases;
210
                   if (formulation == pnsw)
211
                       globalSol[dofIdxGlobal][switchIdx] = 0.0;
212
                   else if (formulation == pwsn)
213
                       globalSol[dofIdxGlobal][switchIdx] = 1.0;
214
215
216
217
218
               }
           }
           else if (phasePresence == wPhaseOnly)
           {

219
               Scalar xwn = volVars.moleFraction(wPhaseIdx, nCompIdx);
220
221
222
223
224
               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;
225
               if (ParentType::staticDat_[dofIdxGlobal].wasSwitched)
226
227
228
229
230
231
                   xwnMax *= 1.02;


               if(xwn > xwnMax)
               {
                   // non-wetting phase appears
232
                   std::cout << "non-wetting phase appears at vertex " << dofIdxGlobal
233
234
235
236
                             << ", coordinates: " << globalPos << ", xwn > xwnMax: "
                             << xwn << " > "<< xwnMax << std::endl;

                   newPhasePresence = bothPhases;
237
                   if (formulation == pnsw)
238
                       globalSol[dofIdxGlobal][switchIdx] = 0.999;
239
                   else if (formulation == pwsn)
240
                       globalSol[dofIdxGlobal][switchIdx] = 0.001;
241
242
243
244
245
               }
           }
           else if (phasePresence == bothPhases)
           {
               Scalar Smin = 0.0;
246
               if (ParentType::staticDat_[dofIdxGlobal].wasSwitched)
247
248
249
250
251
252
                   Smin = -0.01;

               if (volVars.saturation(nPhaseIdx) <= Smin)
               {
                   wouldSwitch = true;
                   // nonwetting phase disappears
253
                   std::cout << "Nonwetting phase disappears at vertex " << dofIdxGlobal
254
                             << ", coordinates: " << globalPos << ", sn: "
255
256
257
                             << volVars.saturation(nPhaseIdx) << std::endl;
                   newPhasePresence = wPhaseOnly;

258
259
                   if(!useMoles) //mass-fraction formulation
                   {
260
					   globalSol[dofIdxGlobal][switchIdx]
261
						   = volVars.massFraction(wPhaseIdx, nCompIdx);
262
263
264
                   }
                   else //mole-fraction formulation
                   {
265
					   globalSol[dofIdxGlobal][switchIdx]
266
					       = volVars.moleFraction(wPhaseIdx, nCompIdx);
267
                   }
268
269
270
271
272
               }
               else if (volVars.saturation(wPhaseIdx) <= Smin)
               {
                   wouldSwitch = true;
                   // wetting phase disappears
273
                   std::cout << "Wetting phase disappears at vertex " << dofIdxGlobal
274
                             << ", coordinates: " << globalPos << ", sw: "
275
276
277
                             << volVars.saturation(wPhaseIdx) << std::endl;
                   newPhasePresence = nPhaseOnly;

278
279
                   if(!useMoles) //mass-fraction formulation
                   {
280
					   globalSol[dofIdxGlobal][switchIdx]
281
						   = volVars.massFraction(nPhaseIdx, wCompIdx);
282
283
284
                   }
                   else //mole-fraction formulation
                   {
285
						globalSol[dofIdxGlobal][switchIdx]
286
						= volVars.moleFraction(nPhaseIdx, wCompIdx);
287
                   }
288
289
290
               }
           }

291
292
           ParentType::staticDat_[dofIdxGlobal].phasePresence = newPhasePresence;
           ParentType::staticDat_[dofIdxGlobal].wasSwitched = wouldSwitch;
293
294
295
296
297
298
299
300
301
           return phasePresence != newPhasePresence;
       }


};

}

#endif