[com1] Rheological Models for debris-flow simulations

avaframe/com1DFA/DFAfunctionsCython.pyx

@@ -86,6 +86,24 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double muCoulomb = cfg.getfloat('mucoulomb')
   cdef double muCoulombMinShear = cfg.getfloat('mucoulombminshear')
   cdef double tau0CoulombMinShear = cfg.getfloat('tau0coulombminshear')
+  cdef double rhoSediment = cfg.getfloat('rhoSediment')
+  cdef double sizeSediment = cfg.getfloat('sizeSediment')
+  cdef double cvMaxSediment = cfg.getfloat('cvMaxSediment')
+  cdef double cvSediment = cfg.getfloat('cvSediment')
+  cdef double alpha1EtaObrienAndJulien = cfg.getfloat('alpha1EtaObrienAndJulien')
+  cdef double beta1EtaObrienAndJulien = cfg.getfloat('beta1EtaObrienAndJulien')
+  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 nHerschelAndBulkley = cfg.getfloat('nHerschelAndBulkley')
+  cdef double alpha1EtaBingham = cfg.getfloat('alpha1EtaBingham')
+  cdef double beta1EtaBingham = cfg.getfloat('beta1EtaBingham')
+  cdef double alpha2TauyBingham = cfg.getfloat('alpha2TauyBingham')
+  cdef double beta2TauyBingham = cfg.getfloat('beta2TauyBingham')
   cdef double curvAccInFriction = cfg.getfloat('curvAccInFriction')
   cdef double curvAccInTangent = cfg.getfloat('curvAccInTangent')
   cdef int curvAccInGradient = cfg.getint('curvAccInGradient')
@@ -151,6 +169,8 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double nx, ny, nz, nxEnd, nyEnd, nzEnd, nxAvg, nyAvg, nzAvg
   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                                                                          
   # variables for interpolation
   cdef int Lx0, Ly0, LxEnd0, LyEnd0, iCell, iCellEnd
   cdef double w[4]
@@ -266,6 +286,13 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
       else:
           # bottom normal stress sigmaB
           sigmaB = - effAccNorm * rho * h
+          if frictType >= 10:
+            # substitution of shear rate gamma
+            if math.fabs(h) <= 1e-9:
+               message = "h is 0 - division by 0"
+               log.error(message)
+               raise ZeroDivisionError(message)
+            shearRate = 3 * uMag / h
           if frictType == 1:
             # SamosAT friction type (bottom shear stress)
             tau = DFAtlsC.SamosATfric(rho, tau0SamosAt, Rs0SamosAt, muSamosAt, kappaSamosAt, BSamosAt, RSamosAt, uMag, sigmaB, h)
@@ -298,6 +325,36 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
             xsiVoellmyRaster = xsiRaster[indCellY, indCellX]
             # Voellmy with optional spatially variable mu and xi values provided as rasters
             tau = muVoellmyRaster * sigmaB + rho * uMag * uMag * gravAcc / xsiVoellmyRaster
+          elif frictType == 10:
+            ## O`Brien and Julien
+            # viscosity
+            etaObrienAndJulien = alpha1EtaObrienAndJulien * math.exp(beta1EtaObrienAndJulien * cvSediment)
+            # yield shear stress
+            tauyObrienAndJulien = alpha2TauyObrienAndJulien * math.exp(beta2TauyObrienAndJulien * cvSediment)
+            # Prandtl mixing length
+            lmObrienAndJulien = 0.4 * h
+            # grain concentration
+            lambdaBagnold = 1.0 / (math.pow(cvMaxSediment / cvSediment, 1.0 / 3.0) - 1.0)
+            # dispersive shear stress
+            cObrienAndJulien = rho * lmObrienAndJulien * lmObrienAndJulien + alphaObrienAndJulien * rhoSediment * lambdaBagnold * lambdaBagnold * sizeSediment * sizeSediment
+            # shear stress
+            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)
+          elif frictType == 12:
+            ## Bingham
+            # viscosity
+            etaBingham = alpha1EtaBingham * math.exp(beta1EtaBingham * cvSediment)
+            # yield shear stress
+            tauyBingham = alpha2TauyBingham * math.exp(beta2TauyBingham * cvSediment)
+            # shear stress
+            tau = tauyBingham + etaBingham * shearRate
           else:
             tau = 0.0
 

avaframe/com1DFA/checkCfg.py

@@ -130,6 +130,20 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
         "tau0coulombminshear",
         "muvoellmyminshear",
         "xsivoellmyminshear",
+        "alpha1EtaObrienAndJulien",
+        "beta1EtaObrienAndJulien",
+        "alpha2TauyObrienAndJulien",
+        "beta2TauyObrienAndJulien",
+        "alphaObrienAndJulien",
+        "alpha1EtaHerschelAndBulkley",
+        "beta1EtaHerschelAndBulkley",
+        "alpha2TauyHerschelAndBulkley",
+        "beta2TauyHerschelAndBulkley",
+        "nHerschelAndBulkley",
+        "alpha1EtaBingham",
+        "beta1EtaBingham",
+        "alpha2TauyBingham",
+        "beta2TauyBingham",
     ]
 
     # if samosATAuto check for release volume and volume class to determine which parameter setup
@@ -168,6 +182,21 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
                 frictModel.lower() in frictP and "voellmyminshear" not in frictP
                 for frictP in frictParameters
             ]
+        elif frictModel.lower() == "obrienandjulien":
+            frictParameterIn = [
+                frictModel.lower() in frictP.lower()
+                for frictP in frictParameters
+            ]
+        elif frictModel.lower() == "herschelandbulkley":
+            frictParameterIn = [
+                frictModel.lower() in frictP.lower()
+                for frictP in frictParameters
+            ]
+        elif frictModel.lower() == "bingham":
+            frictParameterIn = [
+                frictModel.lower() in frictP.lower()
+                for frictP in frictParameters
+            ]
         else:
             frictParameterIn = [frictModel.lower() in frictP for frictP in frictParameters]
 
@@ -217,6 +246,15 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
                     cfg["INPUT"]["xiFile"],
                 )
             )
+    # O´Brien and Julien friction model
+    # fetch parameters and check if cvSediment < cvMaxSediment
+    if frictModel.lower() == "obrienandjulien":
+        cvSediment = cfg['GENERAL']['cvSediment']
+        cvMaxSediment = cfg['GENERAL']['cvMaxSediment']
+        if cvSediment >= cvMaxSediment:
+            message = "cvSediment must be smaller than cvMaxSediment"
+            log.error(message)
+            raise ValueError(message)
 
     # check if explicitFriction has valid values
     if cfg["GENERAL"]["explicitFriction"] not in ["0", "1"]:

avaframe/com1DFA/com1DFA.py

@@ -808,6 +808,37 @@ def createReportDict(avaDir, logName, relName, inputSimLines, cfg, reportAreaInf
             "mu file": inputSimLines["muFile"].name,
             "xi file": inputSimLines["xiFile"].name,
         }
+    elif cfgGen["frictModel"].lower() == "obrienandjulien":
+        reportST["Friction model"] = {
+            "type": "columns",
+            "model": "O´Brien and Julien",
+            "alpha": cfgGen["alphaObrienAndJulien"],
+            "Cmax": cfgGen["cvMaxSediment"],
+            "Cv": cfgGen["cvSediment"],
+            "alpha1": cfgGen["alpha1EtaObrienAndJulien"],
+            "beta1": cfgGen["beta1EtaObrienAndJulien"],
+            "alpha2": cfgGen["alpha2TauyObrienAndJulien"],
+            "beta2": cfgGen["beta2TauyObrienAndJulien"],
+        }
+    elif cfgGen["frictModel"].lower() == "herschelandbulkley":
+        reportST["Friction model"] = {
+            "type": "columns",
+            "model": "Herschel and Bulkley",
+            "n": cfgGen["nHerschelAndBulkley"],
+            "alpha1": cfgGen["alpha1EtaHerschelAndBulkley"],
+            "beta1": cfgGen["beta1EtaHerschelAndBulkley"],
+            "alpha2": cfgGen["alpha2TauyHerschelAndBulkley"],
+            "beta2": cfgGen["beta2TauyHerschelAndBulkley"],
+        }
+    elif cfgGen["frictModel"].lower() == "bingham":
+        reportST["Friction model"] = {
+            "type": "columns",
+            "model": "Bingham",
+            "alpha1": cfgGen["alpha1EtaBingham"],
+            "beta1": cfgGen["beta1EtaBingham"],
+            "alpha2": cfgGen["alpha2TauyBingham"],
+            "beta2": cfgGen["beta2TauyBingham"],
+        }
 
     # check if secondary release area
     if secRelAreaFlag == "Yes":
@@ -1855,6 +1886,9 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         "voellmyminshear",
         "coulombminshear",
         "spatialvoellmy",
+        "obrienandjulien",
+        "herschelandbulkley",
+        "bingham"
     ]
     frictModel = cfgGen["frictModel"].lower()
     frictType = frictModelsList.index(frictModel) + 1

avaframe/com1DFA/com1DFACfg.ini

@@ -284,6 +284,7 @@ centeredPosition = 1
 # or use an implicit method (0)
 explicitFriction = 0
 # friction type. Available options are:
+### Granular friction models
 # samosAT - standard setup
 # samosATSmall - setup for release volumes < 25k m³
 # samosATMedium - setup for release volumes < 60k m³
@@ -294,6 +295,10 @@ explicitFriction = 0
 # spatialVoellmy
 # CoulombMinShear
 # wetsnow
+### Rheological models
+# OBrienAndJulien
+# HerschelAndBulkley
+# Bingham
 # Please note that each type has their own/separate parameters:
 # https://docs.avaframe.org/en/latest/theoryCom1DFA.html#friction-model
 frictModel = samosATAuto
@@ -352,6 +357,56 @@ TIni = -10
 entTempRef = -10.0
 # J/kgK (ice)
 cpIce = 2050.0
+#+++++++++DEBRIS FLOW properties
+# density of sediment [kg/m³]
+rhoSediment = 2650
+# mean sediment size [m]
+## TODO  adjusting
+sizeSediment = 0.002
+# maximum concentration of sediment particles [m³/m³]
+cvMaxSediment = 0.615
+# volumetric sediment concentration [m³/m³]
+## TODO adjusting
+cvSediment = 0.5
+#++++++++++++++Parameter Rheological models
+#++++++++++++O´Brien and Julien
+# viscosity parameter
+## TODO adjusting, default value O´Brien and Julien (1985)
+alpha1EtaObrienAndJulien = 0.650
+## TODO adjusting, default value O´Brien and Julien (1985)
+beta1EtaObrienAndJulien = 16.81
+# yield shear stress parameter
+## TODO adjusting, default value O´Brien and Julien (1985)
+alpha2TauyObrienAndJulien = 0.00886
+ ## TODO adjusting, default value O´Brien and Julien (1985)
+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
+# flow exponent
+## TODO adjusting
+nHerschelAndBulkley = 1.5
+#++++++++++++Bingham
+# viscosity parameter
+## TODO adjusting, default value O´Brien and Julien (1985)
+alpha1EtaBingham = 0.650
+ ## TODO adjusting, default value O´Brien and Julien (1985)
+beta1EtaBingham = 16.81
+# yield shear stress parameter
+ ## TODO adjusting, default value O´Brien and Julien (1985)
+alpha2TauyBingham = 0.00886
+ ## TODO adjusting, default value O´Brien and Julien (1985)
+beta2TauyBingham = 13.11
 
 #++++++++++++ Resistance model
 # default setup: