[com1] Replace viscosity by consistency factor for Herschel-Bulkley rheology

avaframe/com1DFA/DFAfunctionsCython.pyx

@@ -95,10 +95,9 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double alpha2TauyObrienAndJulien = cfg.getfloat('alpha2TauyObrienAndJulien')
   cdef double beta2TauyObrienAndJulien = cfg.getfloat('beta2TauyObrienAndJulien')
   cdef double alphaObrienAndJulien = cfg.getfloat('alphaObrienAndJulien')
-  cdef double alpha1EtaHerschelAndBulkley = cfg.getfloat('alpha1EtaHerschelAndBulkley')
-  cdef double beta1EtaHerschelAndBulkley = cfg.getfloat('beta1EtaHerschelAndBulkley')
   cdef double alpha2TauyHerschelAndBulkley = cfg.getfloat('alpha2TauyHerschelAndBulkley')
   cdef double beta2TauyHerschelAndBulkley = cfg.getfloat('beta2TauyHerschelAndBulkley')
+  cdef double kHerschelAndBulkley = cfg.getfloat('kHerschelAndBulkley')
   cdef double nHerschelAndBulkley = cfg.getfloat('nHerschelAndBulkley')
   cdef double alpha1EtaBingham = cfg.getfloat('alpha1EtaBingham')
   cdef double beta1EtaBingham = cfg.getfloat('beta1EtaBingham')
@@ -170,7 +169,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double gravAccNorm, accNormCurv, effAccNorm, gravAccTangX, gravAccTangY, gravAccTangZ, forceBotTang, sigmaB, tau
   cdef double muVoellmyRaster, xsiVoellmyRaster
   cdef double shearRate, etaObrienAndJulien, tauyObrienAndJulien, lmObrienAndJulien
-  cdef double lambdaBagnold, cObrienAndJulien, etaHerschelAndBulkley, tauyHerschelAndBulkley, etaBingham, tauyBingham                                                                          
+  cdef double lambdaBagnold, cObrienAndJulien, tauyHerschelAndBulkley, etaBingham, tauyBingham                                                                          
   # variables for interpolation
   cdef int Lx0, Ly0, LxEnd0, LyEnd0, iCell, iCellEnd
   cdef double w[4]
@@ -340,12 +339,10 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
             tau = tauyObrienAndJulien + etaObrienAndJulien * shearRate + cObrienAndJulien * (shearRate * shearRate)
           elif frictType == 11:
             ## Herschel and Bulkley
-            # viscosity
-            etaHerschelAndBulkley = alpha1EtaHerschelAndBulkley * math.exp(beta1EtaHerschelAndBulkley * cvSediment)
             # yield shear stress
             tauyHerschelAndBulkley = alpha2TauyHerschelAndBulkley * math.exp(beta2TauyHerschelAndBulkley * cvSediment)
             # shear stress
-            tau = tauyHerschelAndBulkley + etaHerschelAndBulkley * math.pow(shearRate, nHerschelAndBulkley)
+            tau = tauyHerschelAndBulkley + kHerschelAndBulkley * math.pow(shearRate, nHerschelAndBulkley)
           elif frictType == 12:
             ## Bingham
             # viscosity

avaframe/com1DFA/checkCfg.py

@@ -135,10 +135,9 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
         "alpha2TauyObrienAndJulien",
         "beta2TauyObrienAndJulien",
         "alphaObrienAndJulien",
-        "alpha1EtaHerschelAndBulkley",
-        "beta1EtaHerschelAndBulkley",
         "alpha2TauyHerschelAndBulkley",
         "beta2TauyHerschelAndBulkley",
+        "kHerschelAndBulkley",
         "nHerschelAndBulkley",
         "alpha1EtaBingham",
         "beta1EtaBingham",

avaframe/com1DFA/com1DFA.py

@@ -932,8 +932,7 @@ def createReportDict(avaDir, logName, relName, inputSimLines, cfg, reportAreaInf
             "type": "columns",
             "model": "Herschel and Bulkley",
             "n": cfgGen["nHerschelAndBulkley"],
-            "alpha1": cfgGen["alpha1EtaHerschelAndBulkley"],
-            "beta1": cfgGen["beta1EtaHerschelAndBulkley"],
+            "K": cfgGen["kHerschelAndBulkley"],
             "alpha2": cfgGen["alpha2TauyHerschelAndBulkley"],
             "beta2": cfgGen["beta2TauyHerschelAndBulkley"],
         }

avaframe/com1DFA/com1DFACfg.ini

@@ -401,16 +401,14 @@ beta2TauyObrienAndJulien = 13.11
 # parameter, default value Takahashi (1978)
 alphaObrienAndJulien = 0.01
 #++++++++++++Herschel and Bulkley
-# viscosity parameter
-## TODO adjusting, default value O´Brien and Julien (1985)
-alpha1EtaHerschelAndBulkley = 0.650
- ## TODO adjusting, default value O´Brien and Julien (1985)
-beta1EtaHerschelAndBulkley = 16.81
 # yield shear stress parameter
  ## TODO adjusting, default value O´Brien and Julien (1985)
 alpha2TauyHerschelAndBulkley = 0.00886
  ## TODO adjusting, default value O´Brien and Julien (1985)
 beta2TauyHerschelAndBulkley = 13.11
+# consistency factor K
+## TODO adjusting
+kHerschelAndBulkley = 50
 # flow exponent
 ## TODO adjusting
 nHerschelAndBulkley = 1.5

avaframe/tests/test_DFAfunctionsCython.py

@@ -839,10 +839,9 @@ def test_computeForceC():
         "mucoulomb": "0.155",
         "mucoulombminshear": "0.155",
         "tau0coulombminshear": "70",
-        "alpha1EtaHerschelAndBulkley": "0.650",
-        "beta1EtaHerschelAndBulkley": "16.81",
         "alpha2TauyHerschelAndBulkley": "0.00886",
         "beta2TauyHerschelAndBulkley": "13.11",
+        "kHerschelAndBulkley": "50",
         "nHerschelAndBulkley": "1.5",
         "alpha1EtaBingham": "0.650",
         "beta1EtaBingham": "16.81",
@@ -985,19 +984,15 @@ def test_computeForceC():
 
     ## Herschel and Bulkley
     # read input parameters
-    test_alpha1EtaHerschelAndBulkley = float(cfg["GENERAL"]["alpha1EtaHerschelAndBulkley"])
-    test_beta1EtaHerschelAndBulkley = float(cfg["GENERAL"]["beta1EtaHerschelAndBulkley"])
     test_alpha2TauyHerschelAndBulkley = float(cfg["GENERAL"]["alpha2TauyHerschelAndBulkley"])
     test_beta2TauyHerschelAndBulkley = float(cfg["GENERAL"]["beta2TauyHerschelAndBulkley"])
+    test_kHerschelAndBulkley = float(cfg["GENERAL"]["kHerschelAndBulkley"])
     test_nHerschelAndBulkley = float(cfg["GENERAL"]["nHerschelAndBulkley"])
     # compute tau
-    test_etaHerschelandBulkley = test_alpha1EtaHerschelAndBulkley * math.exp(
-        test_beta1EtaHerschelAndBulkley * test_cvSediment
-    )
     test_tauyHerschelandBulkley = test_alpha2TauyHerschelAndBulkley * math.exp(
         test_beta2TauyHerschelAndBulkley * test_cvSediment
     )
-    test_tauHerschelAndBulkley = test_tauyHerschelandBulkley + test_etaHerschelandBulkley * math.pow(
+    test_tauHerschelAndBulkley = test_tauyHerschelandBulkley + test_kHerschelAndBulkley * math.pow(
         test_shearRate[0], test_nHerschelAndBulkley
     )
     test_tauArray[1] = test_tauHerschelAndBulkley

docs/theoryCom1DFA.rst

@@ -480,17 +480,18 @@ With
    &\mathbf{x} \qquad &\text{position vector}\end{aligned}
 
 For solid-fluid mixtures like debris flows, where the properties of the fluid phase are dominating the flow process,
-the shear stress tensor is, among other things, a function of dynamic viscosity :math:`\eta_m` and shear rate :math:`\dot\gamma` (:ref:`theoryCom1DFA:Rheological models`).
+the shear stress tensor is, among other things, a function of a consistency factor :math:`K` 
+and the shear rate :math:`\dot\gamma` (:ref:`theoryCom1DFA:Rheological models`).
 
 .. math::
-    \tau^{(b)}_i = f(\eta_m,\dot\gamma,\overline{h},\rho_m,t,\mathbf{x})
+    \tau^{(b)}_i = f(K,\dot\gamma,\overline{h},\rho_m,t,\mathbf{x})
     :label: rheological model
 
 With
 
 .. math::
    \begin{aligned}
-   &\eta_m \qquad &\text{dynamic viscosity of the solid-fluid mixture}\\
+   &K \qquad &\text{consistency factor of the solid-fluid mixture}\\
    &\dot\gamma \qquad &\text{shear rate}\\
    &\overline{h} \qquad &\text{average flow thickness}\\
    &\rho_m \qquad &\text{density of the solid-fluid mixture}\\
@@ -662,12 +663,12 @@ and non-Newtonian fluids, in which a yield shear stress has to be exceeded or sh
 The rheological models are incorporated into a single, general form, which can be expressed as follows:
 
 .. math::
-    \tau^{(b)} = \tau_y + \eta_m \cdot \dot\gamma^n + C \cdot \dot\gamma^2
+    \tau^{(b)} = \tau_y + K \cdot \dot\gamma^n + C \cdot \dot\gamma^2
     :label: rheology-general
 
 The yield shear stress :math:`\tau_y` :math:`[Pa]` defines a lower limit below which no flow takes place (see :math:`\tau_0` at :ref:`samosatfrict`). 
-The dynamic viscosity :math:`\eta_m` :math:`[Pa \cdot s]` quantifies the internal frictional force between two neighbouring layers of the mixture in relative motion. 
-:math:`\dot\gamma` is the flow velocity gradient, or shear rate, :math:`\frac{du}{dz}` :math:`[s^{-1}]` along the axis normal to the slope (flow thickness).
+:math:`K` is a consistency factor :math:`[Pa \cdot s^{n}]`. :math:`\dot\gamma` is the flow velocity gradient,
+or shear rate, :math:`\frac{du}{dz}` :math:`[s^{-1}]` along the axis normal to the slope (flow thickness).
 The flow index :math:`n` describes the rheological behaviour of the mixture as :cite:`KaZhHaHe2025`:
 
 - :math:`n = 0` a rate-independent solid-like behaviour,
@@ -679,22 +680,24 @@ The flow index :math:`n` describes the rheological behaviour of the mixture as :
 :math:`C` incorporates the turbulent and dispersive shear stresses :math:`[kg \cdot m^{-1}]`,
 which considers the inertial impact between the mixture particles as well :cite:`ObJu1985`.
 
-Depending on how the parameters are selected, you can choose between three different models. :numref:`Overview-rheological-models-table` gives an overview
-about the relation between the implemented rheological models and the used parameters:
+Three models that vary in complexity use different parameter combinations. For the :ref:`bingham` and the :ref:`obrienandjulien` rheologies (:math:`n = 1`)
+the consistency factor corresponds to the bulk dynamic viscosity :math:`\eta_m` :math:`[Pa \cdot s]` of the which quantifies the internal frictional force
+between two neighbouring layers of the solid-fluid mixture in relative motion.
+:numref:`Overview-rheological-models-table` gives an overview about the relation between the implemented rheological models and the used parameters:
 
 .. _Overview-rheological-models-table:
 
 .. table:: Overview of the implemented rheological models and their parameters according to :eq:`rheology-general`
 
-    +-------------------------------------+-----------------------------+------------------------+------------------------+
-    | :math:`\boldsymbol{Model}`          | :math:`\boldsymbol{\tau_y}` | :math:`\boldsymbol{n}` | :math:`\boldsymbol{C}` |
-    +=====================================+=============================+========================+========================+
-    | *O´Brien and Julien*                | :math:`> 0`                 | :math:`= 1`            | :math:`> 0`            |
-    +-------------------------------------+-----------------------------+------------------------+------------------------+
-    | *Herschel and Bulkley*              | :math:`> 0`                 | :math:`\neq 1`         | :math:`= 0`            |
-    +-------------------------------------+-----------------------------+------------------------+------------------------+
-    | *Bingham*                           | :math:`> 0`                 | :math:`= 1`            | :math:`= 0`            |
-    +-------------------------------------+-----------------------------+------------------------+------------------------+
+    +-------------------------------------+-----------------------------+------------------------+------------------------+------------------------+
+    | :math:`\boldsymbol{Model}`          | :math:`\boldsymbol{\tau_y}` | :math:`\boldsymbol{K}` | :math:`\boldsymbol{n}` | :math:`\boldsymbol{C}` |
+    +=====================================+=============================+========================+========================+========================+
+    | *O´Brien and Julien*                | :math:`> 0`                 | :math:`\eta_m`         | :math:`= 1`            | :math:`> 0`            |
+    +-------------------------------------+-----------------------------+------------------------+------------------------+------------------------+
+    | *Herschel and Bulkley*              | :math:`> 0`                 | :math:`K`              | :math:`\neq 1`         | :math:`= 0`            |
+    +-------------------------------------+-----------------------------+------------------------+------------------------+------------------------+
+    | *Bingham*                           | :math:`> 0`                 | :math:`\eta_m`         | :math:`= 1`            | :math:`= 0`            |
+    +-------------------------------------+-----------------------------+------------------------+------------------------+------------------------+
 
 
 :numref:`Overview-rheological-models-fig` illustrates the behaviour of the graphs of the rheological models due to shear rate.
@@ -731,6 +734,8 @@ in following substitution :cite:`Iv1997, DeIv2001, GiFlSaHe2020`:
 
 where :math:`\overline{u}` is the depth-averaged flow velocity and :math:`\overline{h}` is the flow thickness.
 
+.. _obrienandjulien:
+
 O´Brien and Julien
 +++++++++++++++++++
 The quadratic rheological model proposed by O´Brien and Julien :cite:`ObJu1985` accounts for various shear stress components, including
@@ -761,17 +766,21 @@ of sediment :math:`[kg \cdot m^{-3}]`, :math:`d_s` is the sediment size :math:`[
 The quadratic rheological model is particularly suitable for simulating debris flows with high sediment concentrations, in which energy dissipation and
 resistance due to turbulence and particle collisions are significant.
 
+.. _herschelandbulkley:
+
 Herschel and Bulkley
 +++++++++++++++++++++++++
 The model by Herschel and Bulkley :cite:`HeBu1926` is expressed by an empirical power-law equation:
 
 .. math::
-    \tau^{(b)} = \tau_y + \eta_m \cdot \dot{\gamma}^n
+    \tau^{(b)} = \tau_y + K \cdot \dot{\gamma}^n
     :label: herschelAndBulkley
 
 where the flow index :math:`n` describes the rheological behaviour of the mixture (see above). The factor in front of the shear rate
-was originally introduced as a consistency factor :math:`K`. In :py:mod:`com1DFA` :math:`K` equals to the bulk dynamic viscosity :math:`\eta_m`.
-This rheology applies to fine-grained soil-water mixtures that exhibit shear-thinning or shear-thickening behavior, respectively, with increasing shear rates.
+was originally introduced as a consistency factor :math:`K`. This rheology applies to fine-grained soil-water mixtures that
+exhibit shear-thinning or shear-thickening behavior, respectively, with increasing shear rates.
+
+.. _bingham:
 
 Bingham
 ++++++++++++++++