[cig-commits] r1313 - in trunk/aspect: include/aspect/material_model source/material_model

[gassmoeller at dealii.org]

Author: gassmoeller
Date: 2012-10-23 05:47:02 -0600 (Tue, 23 Oct 2012)
New Revision: 1313

Modified:
   trunk/aspect/include/aspect/material_model/simple.h
   trunk/aspect/source/material_model/simple.cc
Log:
Added a compositional dependency of viscosity in simple material model. Switched off by default.

Modified: trunk/aspect/include/aspect/material_model/simple.h
===================================================================
--- trunk/aspect/include/aspect/material_model/simple.h	2012-10-23 03:41:51 UTC (rev 1312)
+++ trunk/aspect/include/aspect/material_model/simple.h	2012-10-23 11:47:02 UTC (rev 1313)
@@ -181,6 +181,7 @@
         double reference_rho;
         double reference_T;
         double eta;
+        double composition_viscosity_prefactor;
         double thermal_viscosity_exponent;
         double thermal_alpha;
         double reference_specific_heat;

Modified: trunk/aspect/source/material_model/simple.cc
===================================================================
--- trunk/aspect/source/material_model/simple.cc	2012-10-23 03:41:51 UTC (rev 1312)
+++ trunk/aspect/source/material_model/simple.cc	2012-10-23 11:47:02 UTC (rev 1313)
@@ -42,12 +42,16 @@
       const double delta_temp = temperature-reference_T;
       const double temperature_dependence = std::max(std::min(std::exp(-thermal_viscosity_exponent*delta_temp/reference_T),1e2),1e-2);
 
-      return temperature_dependence * eta;
-      /*      return (this->n_compositional_fields()>0
-                ?
-                (6.5*composition[0]+1) * eta
-                :
-                eta);*/
+      double composition_dependence = 1.0;
+      if ((composition_viscosity_prefactor != 1.0) && (composition.size() > 0))
+      {
+          // TODO: Currently using an arithmetic interpolation. Usually a geometric interpolation 
+          // is assumed more realistic for viscosity
+          composition_dependence *= (1-composition[0]);
+          composition_dependence += composition_viscosity_prefactor * composition[0];
+      }
+
+      return composition_dependence * temperature_dependence * eta;
     }
 
 
@@ -223,6 +227,9 @@
           prm.declare_entry ("Viscosity", "5e24",
                              Patterns::Double (0),
                              "The value of the constant viscosity. Units: $kg/m/s$.");
+          prm.declare_entry ("Composition viscosity prefactor", "1.0",
+                             Patterns::Double (0),
+                             "A linear dependency of viscosity on composition. Dimensionless prefactor.");
           prm.declare_entry ("Thermal viscosity exponent", "0.0",
                              Patterns::Double (0),
                              "The temperature dependence of viscosity. Dimensionless exponent.");
@@ -268,6 +275,7 @@
           reference_rho              = prm.get_double ("Reference density");
           reference_T                = prm.get_double ("Reference temperature");
           eta                        = prm.get_double ("Viscosity");
+          composition_viscosity_prefactor = prm.get_double ("Composition viscosity prefactor");
           thermal_viscosity_exponent = prm.get_double ("Thermal viscosity exponent");
           k_value                    = prm.get_double ("Thermal conductivity");
           reference_specific_heat    = prm.get_double ("Reference specific heat");