From 5b14a522f2ad4e4c73c0ab3cb2c7cd46b6c3bbc7 Mon Sep 17 00:00:00 2001
From: Holger Class <holger.class@iws.uni-stuttgart.de>
Date: Thu, 24 Oct 2013 14:40:15 +0000
Subject: [PATCH] adaption of 3p3c to usage of EffectiveDiffusivityModel

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11803 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/implicit/3p3c/3p3cfluxvariables.hh    | 50 ++++++++++++++-------
 dumux/implicit/3p3c/3p3cproperties.hh       |  1 +
 dumux/implicit/3p3c/3p3cpropertydefaults.hh |  9 ++++
 3 files changed, 44 insertions(+), 16 deletions(-)

diff --git a/dumux/implicit/3p3c/3p3cfluxvariables.hh b/dumux/implicit/3p3c/3p3cfluxvariables.hh
index 0a8e7b343b..578c111826 100644
--- a/dumux/implicit/3p3c/3p3cfluxvariables.hh
+++ b/dumux/implicit/3p3c/3p3cfluxvariables.hh
@@ -55,6 +55,7 @@ class ThreePThreeCFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariable
 
     typedef typename GridView::template Codim<0>::Entity Element;
     typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffectiveDiffusivityModel;
 
     enum {
         dim = GridView::dimension,
@@ -108,7 +109,7 @@ public:
         }
 
         calculateGradients_(problem, element, elemVolVars);
-        calculateporousDiffCoeff_(problem, element, elemVolVars);
+        calculatePorousDiffCoeff_(problem, element, elemVolVars);
     };
 
 private:
@@ -221,7 +222,7 @@ private:
         return sp.eval(sat);
     }
 
-    void calculateporousDiffCoeff_(const Problem &problem,
+    void calculatePorousDiffCoeff_(const Problem &problem,
                                const Element &element,
                                const ElementVolumeVariables &elemVolVars)
     {
@@ -232,6 +233,10 @@ private:
         Dune::FieldMatrix<Scalar, numPhases, numComponents> diffusionCoefficientMatrix_i = volVarsI.diffusionCoefficient();
         Dune::FieldMatrix<Scalar, numPhases, numComponents> diffusionCoefficientMatrix_j = volVarsJ.diffusionCoefficient();
 
+        // the effective diffusion coefficients at vertex i and j
+        Scalar diffCoeffI;
+        Scalar diffCoeffJ;
+
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
             // make sure to calculate only diffusion coefficents
@@ -248,25 +253,38 @@ private:
                 continue;
             }
 
-            // calculate tortuosity at the nodes i and j needed
-            // for porous media diffusion coefficient
+            // Diffusion coefficient in the porous medium
+            diffCoeffI = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsI.porosity(),
+                                                                         volVarsI.saturation(phaseIdx),
+                                                                         diffusionCoefficientMatrix_i[phaseIdx][wCompIdx]);
+            diffCoeffJ = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsJ.porosity(),
+                                                                         volVarsJ.saturation(phaseIdx),
+                                                                         diffusionCoefficientMatrix_j[phaseIdx][wCompIdx]);
+
+            // -> harmonic mean
+            porousDiffCoeff_[phaseIdx][wCompIdx] = harmonicMean(diffCoeffI, diffCoeffJ);            
 
-            Scalar tauI =
-                1.0/(volVarsI.porosity() * volVarsI.porosity()) *
-                pow(volVarsI.porosity() * volVarsI.saturation(phaseIdx), 7.0/3);
-            Scalar tauJ =
-                1.0/(volVarsJ.porosity() * volVarsJ.porosity()) *
-                pow(volVarsJ.porosity() * volVarsJ.saturation(phaseIdx), 7.0/3);
             // Diffusion coefficient in the porous medium
+            diffCoeffI = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsI.porosity(),
+                                                                         volVarsI.saturation(phaseIdx),
+                                                                         diffusionCoefficientMatrix_i[phaseIdx][nCompIdx]);
+            diffCoeffJ = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsJ.porosity(),
+                                                                         volVarsJ.saturation(phaseIdx),
+                                                                         diffusionCoefficientMatrix_j[phaseIdx][nCompIdx]);
 
             // -> harmonic mean
-            porousDiffCoeff_[phaseIdx][wCompIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][wCompIdx],
-                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][wCompIdx]);
-            porousDiffCoeff_[phaseIdx][nCompIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][nCompIdx],
-                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][nCompIdx]);
-            porousDiffCoeff_[phaseIdx][gCompIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][gCompIdx],
-                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][gCompIdx]);
+            porousDiffCoeff_[phaseIdx][nCompIdx] = harmonicMean(diffCoeffI, diffCoeffJ);
 
+            // Diffusion coefficient in the porous medium
+            diffCoeffI = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsI.porosity(),
+                                                                         volVarsI.saturation(phaseIdx),
+                                                                         diffusionCoefficientMatrix_i[phaseIdx][gCompIdx]);
+            diffCoeffJ = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsJ.porosity(),
+                                                                         volVarsJ.saturation(phaseIdx),
+                                                                         diffusionCoefficientMatrix_j[phaseIdx][gCompIdx]);
+
+            // -> harmonic mean
+            porousDiffCoeff_[phaseIdx][gCompIdx] = harmonicMean(diffCoeffI, diffCoeffJ);
         }
     }
 
diff --git a/dumux/implicit/3p3c/3p3cproperties.hh b/dumux/implicit/3p3c/3p3cproperties.hh
index 35e2249457..4cfd9f7287 100644
--- a/dumux/implicit/3p3c/3p3cproperties.hh
+++ b/dumux/implicit/3p3c/3p3cproperties.hh
@@ -59,6 +59,7 @@ NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
 
 NEW_PROP_TAG(MaterialLaw);   //!< The material law which ought to be used (extracted from the spatial parameters)
 NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
+NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity
 
 NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
 NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind parameter for the mobility
diff --git a/dumux/implicit/3p3c/3p3cpropertydefaults.hh b/dumux/implicit/3p3c/3p3cpropertydefaults.hh
index d40aab1ddd..8ea4d6fb01 100644
--- a/dumux/implicit/3p3c/3p3cpropertydefaults.hh
+++ b/dumux/implicit/3p3c/3p3cpropertydefaults.hh
@@ -38,6 +38,7 @@
 #include "3p3cnewtoncontroller.hh"
 #include "3p3clocalresidual.hh"
 
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
 #include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include <dumux/material/spatialparams/implicitspatialparams.hh>
 
@@ -126,6 +127,14 @@ SET_TYPE_PROP(ThreePThreeC, Indices, ThreePThreeCIndices<TypeTag, /*PVOffset=*/0
 //! Use ImplicitSpatialParams by default.
 SET_TYPE_PROP(ThreePThreeC, SpatialParams, ImplicitSpatialParams<TypeTag>);
 
+//! The model after Millington (1961) is used for the effective diffusivity
+SET_PROP(ThreePThreeC, EffectiveDiffusivityModel)
+{ private :
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+    typedef DiffusivityMillingtonQuirk<Scalar> type;
+};
+
 // disable velocity output by default
 SET_BOOL_PROP(ThreePThreeC, VtkAddVelocity, false);
 
-- 
GitLab