From 2e60cf789944ff2ca88731d0fa47832b6404494a Mon Sep 17 00:00:00 2001 From: JuLa96 <206707905+JuLa96@users.noreply.github.com> Date: Mon, 16 Mar 2026 19:36:59 +0100 Subject: [PATCH] Replace viscosity by consistency factor for Herschel-Bulkley rheology Update docs/theoryCom1DFA.rst The doc was updated due to changes in the formulation of the Herschel-Bulkley rheology. Co-authored-by: Paula Spannring <95042192+PaulaSp3@users.noreply.github.com> --- avaframe/com1DFA/DFAfunctionsCython.pyx | 9 ++--- avaframe/com1DFA/checkCfg.py | 3 +- avaframe/com1DFA/com1DFA.py | 3 +- avaframe/com1DFA/com1DFACfg.ini | 8 ++-- avaframe/tests/test_DFAfunctionsCython.py | 11 ++--- docs/theoryCom1DFA.rst | 49 ++++++++++++++--------- 6 files changed, 40 insertions(+), 43 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 72b859b23..3ef93f80b 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/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 diff --git a/avaframe/com1DFA/checkCfg.py b/avaframe/com1DFA/checkCfg.py index ef8289a43..83c76f780 100644 --- a/avaframe/com1DFA/checkCfg.py +++ b/avaframe/com1DFA/checkCfg.py @@ -135,10 +135,9 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""): "alpha2TauyObrienAndJulien", "beta2TauyObrienAndJulien", "alphaObrienAndJulien", - "alpha1EtaHerschelAndBulkley", - "beta1EtaHerschelAndBulkley", "alpha2TauyHerschelAndBulkley", "beta2TauyHerschelAndBulkley", + "kHerschelAndBulkley", "nHerschelAndBulkley", "alpha1EtaBingham", "beta1EtaBingham", diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index ef2993f9a..20cad75a2 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/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"], } diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 1a47466e1..dae998349 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/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 diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 816ce665b..731dec875 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/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 diff --git a/docs/theoryCom1DFA.rst b/docs/theoryCom1DFA.rst index da3fce612..8e69afed3 100644 --- a/docs/theoryCom1DFA.rst +++ b/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 ++++++++++++++++