Merge branch 'fix/porosity-deformation' into 'master'

[material][porodeformation] use correct law

See merge request !1183

dumux/material/fluidmatrixinteractions/porositydeformation.hh

@@ -51,6 +51,15 @@ public:
      * \param globalPos The global position (in the element)
      * \param refPoro The solid matrix porosity without deformation
      * \param minPoro A minimum porosity value
+     *
+     * \note Han and Dusseault (2003, doi 10.1016/S0920-4105(03)00047-0)
+     *       provide a derivation for \f$\text{d} \phi = -(1 - \phi ) \text{d} \epsilon_v\f$.
+     *       Here, \f$\epsilon_v\f$ is equal to \f$text{div} \mathbf{u}\f$.
+     *       By using an initial porosity \f$\phi_0\f$ and assuming \epsilon_{v, 0} = 0,
+     *       one obtains \f$\phi = \frac{\phi_0 - \text{div} \mathbf{u}}{1 - \text{div} \mathbf{u}}\f$,
+     *       which is the formulation for the rock mechanics sign convention. Here we are
+     *       using the continuum mechanics sign convention, thus, the final formula reads:
+     *       \f$\phi = \frac{\phi_0 + \text{div} \mathbf{u}}{1 + \text{div} \mathbf{u}}\f$.
      */
     template< class FVGridGeom, class ElemSol >
     static Scalar evaluatePorosity(const FVGridGeom& fvGridGeometry,
@@ -58,7 +67,8 @@ public:
                                    const typename FVGridGeom::GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate& globalPos,
                                    const ElemSol& elemSol,
                                    Scalar refPoro,
-                                   Scalar minPoro = 0.0)
+                                   Scalar minPoro = 0.0,
+                                   Scalar maxPoro = 1.0)
     {
         // compute divergence of diplacement at the given position
         Scalar divU = 0.0;
@@ -66,8 +76,9 @@ public:
         for (int dir = 0; dir < FVGridGeom::GridView::dimension; ++dir)
             divU += gradU[dir][dir];
 
+        using std::min;
         using std::max;
-        return max(minPoro, refPoro*(1.0+divU));
+        return min( maxPoro, max(minPoro, (refPoro+divU)/(1.0+divU)) );
     }
 
     /*!

test/references/poroelasticbox-reference.vtu

@@ -4,17 +4,17 @@
     <Piece NumberOfCells="100" NumberOfPoints="121">
       <PointData Scalars="porosity" Vectors="u">
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
-          0.609054 1.59488 0.641098 1.23626 2.0962 1.51561 1.9182 1.36868 1.12713 0.838561 0.0222568 0.113044
-          0 0 0 0 0 0 0 0 0 0 0.555164 0.507341
-          0.454823 0.398592 0.340366 0.281055 0.221316 0.161671 0.102647 0.0435787 0.0423759 0.369597 0 0
-          0 0 0.458464 1.02452 1.28526 1.11934 0.566846 0.228203 0.15421 0 0 0
-          0 0.576152 1.55383 2.03624 1.82668 0.991113 0.444493 0 0 0 0 0
-          0.588805 1.60612 2.1269 1.95409 1.15444 0.608448 0 0 0 0 0 0.491653
-          1.16136 1.52249 1.45275 0.994345 0.657352 0.0296742 0.0778257 0.130315 0.189185 0.254057 0.322164 0.389881
-          0.454241 0.514274 0.571829 0.572684 0.216191 0.891249 1.25837 1.20443 0.771554 0.145949 0 0
-          0 0.0476512 0.386535 0.437605 1.47327 2.01521 1.87109 1.10929 0.0323899 0 0 0
-          0 0.164678 0.528069 1.53229 2.04988 1.88854 1.11477 0.0277798 0 0 0 0
-          0.0726792
+          0.655203 0.868329 0.672437 0.830133 0.899819 0.861442 0.890522 0.846568 0.813685 0.749571 0 0
+          1 1 1 1 1 1 1 1 1 1 0.621734 0.586077
+          0.538282 0.473145 0.383017 0.252817 0.051129 0 0 0 0 0.431814 1 1
+          1 1 0.541949 0.795026 0.836609 0.812389 0.629529 0.0797655 0 1 1 1
+          1 0.635513 0.86485 0.896869 0.885038 0.788117 0.527552 1 1 1 1 1
+          0.643345 0.86925 0.901265 0.892533 0.818094 0.65486 1 1 1 1 1 0.57287
+          0.819178 0.862068 0.855446 0.788806 0.680536 0 0 0 0 0.173414 0.348159 0.461374
+          0.537691 0.591657 0.632757 0.633306 0.0286361 0.764376 0.833118 0.825643 0.727822 0 1 1
+          1 0 0.456711 0.520115 0.85746 0.895792 0.887766 0.81069 0 1 1 1
+          1 0 0.602325 0.86295 0.897555 0.888803 0.81162 0 1 1 1 1
+          0
         </DataArray>
         <DataArray type="Float32" Name="u" NumberOfComponents="3" format="ascii">
           0 0 0 0 0 0 0 0 0 0.0481153 0.363957 0