diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 0877190d7..7e8b8698d 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/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 diff --git a/avaframe/com1DFA/checkCfg.py b/avaframe/com1DFA/checkCfg.py index 205abb4cb..ef8289a43 100644 --- a/avaframe/com1DFA/checkCfg.py +++ b/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"]: diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index f7f6aa3fd..ccac63e7d 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/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 diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 58391766b..d7085e8f5 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/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: diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index ede4dc84e..488ccf91b 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -9,6 +9,7 @@ # Local imports import avaframe.com1DFA.DFAfunctionsCython as DFAfunC import avaframe.com1DFA.DFAtools as DFAtls +import avaframe.com1DFA.DFAToolsCython as DFAtlsC def test_getNeighborsC(capfd): @@ -783,6 +784,226 @@ def test_computeCohesionForceC(): assert force["forceSPHY"][2] < 0 assert force["forceSPHY"][1] == -force["forceSPHY"][2] +def test_computeForceC(): + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "alpha1EtaObrienAndJulien": "0.650", + "beta1EtaObrienAndJulien": "16.81", + "alpha2TauyObrienAndJulien": "0.00886", + "beta2TauyObrienAndJulien": "13.11", + "alphaObrienAndJulien": "0.01", + "rhoSediment": "2650.0", + "sizeSediment": "0.002", + "cvMaxSediment": "0.615", + "cvSediment": "0.5", + "rho": "2000.0", + "rhoEnt": "100.0", + "gravAcc": "9.81", + "velMagMin": "1.0e-6", + "depMin": "1.0e-6", + "interpOption": "2", + "explicitFriction": "0", + "reprojMethodForce": "2", + "reprojectionIterations": "5", + "thresholdProjection": "0.001", + "subgridMixingFactor": "100.00", + "viscOption": "0", # 0: no viscosity -> ux as defined below + "enthRef": "20000", + "musamosat": "0.155", + "tau0samosat": "0", + "Rs0samosat": "0.222", + "kappasamosat": "0.43", + "Bsamosat": "4.13", + "Rsamosat": "0.05", + "musamosatsmall": "0.22", + "tau0samosatsmall": "0", + "Rs0samosatsmall": "0.222", + "kappasamosatsmall": "0.43", + "Bsamosatsmall": "4.13", + "Rsamosatsmall": "0.05", + "musamosatmedium": "0.17", + "tau0samosatmedium": "0", + "Rs0samosatmedium": "0.222", + "kappasamosatmedium": "0.43", + "Bsamosatmedium": "4.13", + "Rsamosatmedium": "0.05", + "entEroEnergy": "5000", + "entShearResistance": "0", + "entDefResistance": "0", + "xsivoellmy": "4000.", + "muvoellmy": "0.155", + "xsivoellmyminshear": "4000.", + "muvoellmyminshear": "0.155", + "tau0voellmyminshear": "70", + "mucoulomb": "0.155", + "mucoulombminshear": "0.155", + "tau0coulombminshear": "70", + "alpha1EtaHerschelAndBulkley": "0.650", + "beta1EtaHerschelAndBulkley": "16.81", + "alpha2TauyHerschelAndBulkley": "0.00886", + "beta2TauyHerschelAndBulkley": "13.11", + "nHerschelAndBulkley": "1.5", + "alpha1EtaBingham": "0.650", + "beta1EtaBingham": "16.81", + "alpha2TauyBingham": "0.00886", + "beta2TauyBingham": "13.11", + "curvAccInFriction": "1", + "curvAccInTangent": "0", + "curvAccInGradient": "0", + "mu0wetsnow": "0.155", + "xsiwetsnow": "4000.", + } + + # set particles-parameters; dictionary at t - see 'com1DFA.DFAfunctionsCython.pyx.computeForceC()' + # Minimal particles setup + particles = { + "nPart": 1, + "x": np.array([0.0]), + "y": np.array([0.0]), + "z": np.array([0.0]), + "ux": np.array([1.0]), + "uy": np.array([0.0]), + "uz": np.array([0.0]), + "m": np.array([1.0]), + "h": np.array([1.0]), + "dt": 0.1, + "ID": np.array([1]), + "totalEnthalpy": np.array([0.0]), + "indXDEM": np.array([0], dtype=np.int32), + "indYDEM": np.array([0], dtype=np.int32), + } + + # set dem-parameters; dictionary - see 'com1DFA.DFAfunctionsCython.pyx.computeForceC()' + # Minimal dem setup + dem = { + "header": { + "cellsize": 5, + "nrows": 2, + "ncols": 2, + }, + "rasterData": np.ones((2, 2)), + # normal vectors of surface + "Nx": np.zeros((2, 2)), + "Ny": np.zeros((2, 2)), + "Nz": np.ones((2, 2)), + "areaRaster": np.ones((2, 2)), # Area of grid cells + "outOfDEM": np.zeros((4,), dtype=bool), # no data mask (used to find out of dem particles) + } + #dem["headerNeighbourGrid"] = dem["header"] + + # set fields parameters; dictionary - see 'com1DFA.DFAfunctionsCython.pyx.computeForceC()' + # Minimal fields setup + nrows = dem["header"]["nrows"] + ncols = dem["header"]["ncols"] + + fields = { + "Vx": np.zeros((nrows, ncols)), + "Vy": np.zeros((nrows, ncols)), + "Vz": np.zeros((nrows, ncols)), + "entrMassRaster": np.zeros((nrows, ncols)), + "entrEnthRaster": np.zeros((nrows, ncols)), + "muField": np.zeros((nrows, ncols)), + "xiField": np.zeros((nrows, ncols)), + "detRaster": np.zeros((nrows, ncols)), + "cResRaster": np.zeros((nrows, ncols)), + } + + + # initialize test array friction force + test_tauArray = np.zeros(3) # saving tau of the three models + test_forceFrict = np.zeros(particles['nPart']) + + # resistanceType + resistanceType = 1 # default; since cResRaster = 0 there is no resistance + + # apply function + # looping over models + forceFrictArray = np.zeros(3) + # i = index of result in array, j = frictType + for i, j in enumerate(range(10,13)): + particles, force,*_ = DFAfunC.computeForceC(cfg["GENERAL"], particles, fields, dem, j, resistanceType) + forceFrictArray[i] = force['forceFrict'][0] + + ## Check output types and some expected values + assert "forceX" in force + assert "forceY" in force + assert "forceZ" in force + assert "forceFrict" in force + assert isinstance(force["forceX"], np.ndarray) + assert isinstance(force["forceFrict"], np.ndarray) + + # Mass should not change (no entrainment/detrainment) + assert np.isclose(particles["m"][0], 1.0) + # Output velocities should be unchanged (no artificial viscosity, no mass change) + assert np.isclose(particles["ux"][0], 1.0) + assert np.isclose(particles["uy"][0], 0.0) + assert np.isclose(particles["uz"][0], 0.0) + + ## Testing rheologcial models for given input parameters + # substitution of shear rate gamma + test_uMag = DFAtlsC.norm(particles['ux'][0], particles['uy'][0], particles['uz'][0]) + test_shearRate = 3 * test_uMag / particles['h'] + + ## O'Brien & Julien + # read input parameters + test_alpha1EtaObrienAndJulien = float(cfg['GENERAL']['alpha1EtaObrienAndJulien']) + test_beta1EtaObrienAndJulien = float(cfg['GENERAL']['beta1EtaObrienAndJulien']) + test_alpha2TauyObrienAndJulien = float(cfg['GENERAL']['alpha2TauyObrienAndJulien']) + test_beta2TauyObrienAndJulien = float(cfg['GENERAL']['beta2TauyObrienAndJulien']) + test_cvSediment = float(cfg['GENERAL']['cvSediment']) + test_cvMaxSediment = float(cfg['GENERAL']['cvMaxSediment']) + test_rho = float(cfg['GENERAL']['rho']) + test_alphaObrienAndJulien = float(cfg['GENERAL']['alphaObrienAndJulien']) + test_rhoSediment = float(cfg['GENERAL']['rhoSediment']) + test_sizeSediment = float(cfg['GENERAL']['sizeSediment']) + # compute tau + test_etaObrienAndJulien = test_alpha1EtaObrienAndJulien * math.exp(test_beta1EtaObrienAndJulien * test_cvSediment) + test_tauyObrienAndJulien = test_alpha2TauyObrienAndJulien * math.exp(test_beta2TauyObrienAndJulien * test_cvSediment) + test_lmObrienAndJulien = 0.4 * particles['h'] + test_lambdaBagnold = 1.0 / (math.pow(test_cvMaxSediment/test_cvSediment, 1.0 / 3.0) - 1.0) + test_cObrienAndJulien = test_rho * test_lmObrienAndJulien * test_lmObrienAndJulien + test_alphaObrienAndJulien * test_rhoSediment * test_lambdaBagnold * test_lambdaBagnold * test_sizeSediment * test_sizeSediment + test_tauObrienAndJulien = test_tauyObrienAndJulien + test_etaObrienAndJulien * test_shearRate + test_cObrienAndJulien * (test_shearRate * test_shearRate) + test_tauArray[0] = test_tauObrienAndJulien[0] + + ## 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_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_shearRate[0], test_nHerschelAndBulkley) + test_tauArray[1] = test_tauHerschelAndBulkley + + ## Bingham + # read input parameters + test_alpha1EtaBingham = float(cfg['GENERAL']['alpha1EtaBingham']) + test_beta1EtaBingham = float(cfg['GENERAL']['beta1EtaBingham']) + test_alpha2TauyBingham = float(cfg['GENERAL']['alpha2TauyBingham']) + test_beta2TauyBingham = float(cfg['GENERAL']['beta2TauyBingham']) + # compute tau + test_etaBingham = test_alpha1EtaBingham * math.exp(test_beta1EtaBingham * test_cvSediment) + test_tauyBingham = test_alpha2TauyBingham * math.exp(test_beta2TauyBingham * test_cvSediment) + test_tauBingham = test_tauyBingham + test_etaBingham * test_shearRate + test_tauArray[2] = test_tauBingham[0] + + ## compute friction force + # deduce area - see 'com1DFA.DFAfunctionsCython.pyx.computeForceC()' + test_areaPart = particles['m'] / (test_rho * particles['h']) + # adding bottom shear resistance contribution - see 'com1DFA.DFAfunctionsCython.pyx.computeForceC()' + test_forceBotTang = - test_areaPart * test_tauArray + # explicit friction = 0 + test_velMagMin = float(cfg['GENERAL']['velMagMin']) + if test_uMag < test_velMagMin: + test_uMagRes = test_velMagMin + else: + test_uMagRes = test_uMag + test_forceFrictArray = test_forceFrict - test_forceBotTang/test_uMagRes + + assert (forceFrictArray == test_forceFrictArray).all() """ TODO: When calling pytest, the following function raises an error ("Fatal Python error: Aborted") diff --git a/docs/_static/Overview_rheological_models.png b/docs/_static/Overview_rheological_models.png new file mode 100644 index 000000000..c39bda107 Binary files /dev/null and b/docs/_static/Overview_rheological_models.png differ diff --git a/docs/references_all.bib b/docs/references_all.bib index f42f8fc7b..62dca87b2 100644 --- a/docs/references_all.bib +++ b/docs/references_all.bib @@ -1825,6 +1825,15 @@ @article{Bi1998 year = {1998} } +@article{Bi1919, + title = {Paint, a plastic material and not a viscous liquid; the measurement of its mobility and yield value}, + number = {19}, + journal = {Proc. Am. Soc. Test. Mater.}, + author = {Bingham, E.C. and Green, H.}, + year = {1919}, + pages = {640--64}, +} + @article{BiLiTu2001, abstract = {Vertical profiles of snowdrift density and concurrent meteorological data were measured over Antarctic snow and blue ice during a 1 month @@ -6239,6 +6248,22 @@ @article{Ge1982 year = {1982} } +@article{GiFlSaHe2020, + title = {Comparing single‐phase, non‐{Newtonian} approaches with experimental results: {Validating} flume‐scale mud and debris flow in {HEC}‐{RAS}}, + volume = {46}, + issn = {0197-9337, 1096-9837}, + shorttitle = {Comparing single‐phase, non‐{Newtonian} approaches with experimental results}, + doi = {10.1002/esp.5044}, + language = {en}, + number = {3}, + urldate = {2025-07-28}, + journal = {Earth Surface Processes and Landforms}, + author = {Gibson, S. and Floyd, I. and Sánchez, A. and Heath, R.}, + month = dec, + year = {2020}, + pages = {540--553}, +} + @article{GiThLaBoSl2012, author = {A. Giacomini and K. Thoeni and C. Lambert and S. Booth and S.W. Sloan}, doi = {10.1016/j.ijrmms.2012.07.030}, @@ -7616,6 +7641,22 @@ @article{HeBiSiGu2011 year = {2011} } +@article{HeBu1926, + title = {Konsistenzmessungen von {Gummi}-{Benzollösungen}}, + volume = {39}, + copyright = {http://www.springer.com/tdm}, + issn = {0303-402X, 1435-1536}, + doi = {10.1007/BF01432034}, + language = {de}, + number = {4}, + urldate = {2025-07-28}, + journal = {Kolloid-Zeitschrift}, + author = {Herschel, Winslow H. and Bulkley, Ronald}, + month = aug, + year = {1926}, + pages = {291--300}, +} + @article{HeGuZa2008, author = {J. Heierli and P. Gumbsch and M. Zaiser}, file = {HeGuZa2008_anticrack_science.pdf}, @@ -9730,6 +9771,19 @@ @article{JoSiAgRo2001 year = {2001} } +@book{Ju2010, + edition = {2}, + title = {Erosion and {Sedimentation}}, + copyright = {https://www.cambridge.org/core/terms}, + isbn = {978-0-521-83038-6 978-0-521-53737-7 978-0-511-80604-9}, + urldate = {2025-09-02}, + publisher = {Cambridge University Press}, + author = {Julien, Pierre Y.}, + month = jun, + year = {2010}, + doi = {10.1017/CBO9780511806049}, +} + @article{JuHu2007, author = {H. Juliussen and O. Humlum}, file = {JuHu2007_PPP.pdf}, @@ -9855,6 +9909,21 @@ @article{KaEnTa1998 year = {1998} } +@article{KaZhHaHe2025, + title = {Classic, modern, and physics-based rheological laws for geophysical granular flows in a landslide hazard chain}, + volume = {269}, + issn = {00128252}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0012825225001655}, + doi = {10.1016/j.earscirev.2025.105204}, + language = {en}, + urldate = {2025-09-02}, + journal = {Earth-Science Reviews}, + author = {Kamali Zarch, Mohsen and Zhang, Limin and Haeri, S. Mohsen and He, Jian}, + month = oct, + year = {2025}, + pages = {105204}, +} + @article{KePeBeAm2006, author = {G.M. Keevil and J. Peakall and J.L. Best and K.J. Amos}, file = {KePeBeAm2006_marinegeology.pdf}, @@ -14206,6 +14275,28 @@ @book{Ob2011b year = {2011} } +@article{ObJu1993, + title = {Two‐{Dimensional} {Water} {Flood} and {Mudflow} {Simulation}}, + volume = {119}, + issn = {0733-9429, 1943-7900}, + doi = {10.1061/(ASCE)0733-9429(1993)119:2(244)}, + language = {en}, + number = {2}, + urldate = {2025-07-28}, + journal = {Journal of Hydraulic Engineering}, + author = {O'Brien, J. S. and Julien, P. Y. and Fullerton, W. T.}, + month = feb, + year = {1993}, + pages = {244--261}, +} + +@article{ObJu1985, + title = {Physical properties and mechanics of hyperconcentrated sediment flows}, + journal = {Proc. ASCE HD Delineation of landslides, flash flood and debris flow Hazards}, + author = {O’brien, JS and Julien, PY}, + year = {1985}, +} + @article{ObKiSc2009, author = {M. Oberguggenberger and J. King and B. Schmelzer}, file = {ObKiSc2009.pdf}, @@ -18172,6 +18263,21 @@ @book{Ta2007 year = {2007} } +@article{Ta1978, + title = {Mechanical {Characteristics} of {Debris} {Flow}}, + volume = {104}, + issn = {0044-796X, 2690-2524}, + doi = {10.1061/JYCEAJ.0005046}, + language = {en}, + number = {8}, + urldate = {2025-08-14}, + journal = {Journal of the Hydraulics Division}, + author = {Takahashi, Tamotsu}, + month = aug, + year = {1978}, + pages = {1153--1169}, +} + @article{TaMa2004, abstract = {Mechanical properties of snow were investigated by means of a vibration response technique in a frequency range from 10 Hz to 1 MHz and a diff --git a/docs/theoryCom1DFA.rst b/docs/theoryCom1DFA.rst index c25129259..da3fce612 100644 --- a/docs/theoryCom1DFA.rst +++ b/docs/theoryCom1DFA.rst @@ -459,9 +459,10 @@ has to be introduced in order to completely solve the equations. Friction Model ~~~~~~~~~~~~~~~~~ -The problem can be solved by introducing a constitutive equation which -describes the basal shear stress tensor :math:`\tau^{(b)}` as a function -of the flow state of the avalanche. +The problem can be solved by introducing constitutive equations which +describe the behaviour of a material due to stress. +For granular materials like avalanches the basal shear stress tensor :math:`\tau^{(b)}` +is expressed as a function of the flow state of the avalanche (:ref:`theoryCom1DFA:Granular friction models`). .. math:: \tau^{(b)}_i = f(\sigma^{(b)},\overline{u},\overline{h},\rho_0,t,\mathbf{x}) @@ -478,18 +479,39 @@ With &t \qquad &\text{time}\\ &\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`). + +.. math:: + \tau^{(b)}_i = f(\eta_m,\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}\\ + &\dot\gamma \qquad &\text{shear rate}\\ + &\overline{h} \qquad &\text{average flow thickness}\\ + &\rho_m \qquad &\text{density of the solid-fluid mixture}\\ + &t \qquad &\text{time}\\ + &\mathbf{x} \qquad &\text{position vector}\end{aligned} + Several friction models already implemented in the simulation tool are -described here. +described in the following sections. +Granular friction models +"""""""""""""""""""""""""" Mohr-Coulomb friction model -""""""""""""""""""""""""""""""" +++++++++++++++++++++++++++++ The Mohr-Coulomb friction model describes the friction interaction between twos solids. The bottom shear stress simply reads: .. math:: \tau^{(b)} = \tan{\delta}\,\sigma^{(b)} + :label: Mohr-Coulomb friction model :math:`\tan{\delta}=\mu` is the friction coefficient (and :math:`\delta` the friction angle). The bottom shear stress linearly increases with the normal stress component :math:`\sigma^{(b)}` (:cite:`Zw2000,BaSaGr1999,WaHuPu2004,Sa2007`). @@ -503,12 +525,13 @@ granular flow. Because of its relative simplicity, this friction model is also very convenient to derive analytic solutions and validate the numerical implementation. Chezy friction model -"""""""""""""""""""""""" +++++++++++++++++++++++ The Chezy friction model describes viscous friction interaction. The bottom shear stress then reads: .. math:: \tau^{(b)} = c_{\text{dyn}}\,\rho_0\,\bar{u}^2 + :label: Chezy friction model :math:`c_{\text{dyn}}` is the viscous friction coefficient. The bottom shear stress is a quadratic function of the velocity. (:cite:`Zw2000,BaSaGr1999,WaHuPu2004,Sa2007`). @@ -517,14 +540,17 @@ This model enables to reach more realistic velocities for avalanche simulations. The draw back is that the avalanche doesn't stop flowing before the slope inclination approaches zero. This implies that the avalanche flows to the lowest local point. +.. _voellmyfrict: + Voellmy friction model -"""""""""""""""""""""" ++++++++++++++++++++++++ Anton Voellmy was a Swiss engineer interested in avalanche dynamics :cite:`Vo1955`. He first had the idea to combine both the Mohr-Coulomb and the Chezy model by summing them up in order to take advantage of both. This leads to the following friction law: .. math:: \tau^{(b)} = \tan{\delta}\,\sigma^{(b)} + \frac{g}{\xi}\,\rho_0\,\bar{u}^2 + :label: Voellmy friction model where :math:`\xi` is the turbulent friction term. This model is described as Voellmy-Fluid :cite:`Sa2004,Sa2007`. @@ -535,25 +561,27 @@ For this option, raster files with values for :math:`\mu` and :math:`\xi` need t VoellmyMinShear friction model -""""""""""""""""""""""""""""""" ++++++++++++++++++++++++++++++++ In order to increase the friction force and make the avalanche flow stop on steeper slopes than with the Voellmy friction relation, a minimum shear stress can be added to the Voellmy friction relation. This minimum value defines a shear stress under which the snowpack doesn’t move, and induces a strong flow deceleration. This expression of the basal layer friction model also resembles the one used in the swiss RAMMS model, where the Voellmy model is modified by adding a yield stress supposed to account for the snow cohesion (https://ramms.slf.ch/en/modules/debrisflow/theory/friction-parameters.html). .. math:: \tau^{(b)} = \tau_0 + \tan{\delta}\,\sigma^{(b)} + \frac{g}{\xi}\,\rho_0\,\bar{u}^2 + :label: VoellmyMinShear friction model .. _samosatfrict: SamosAT friction model -"""""""""""""""""""""""" ++++++++++++++++++++++++ SamosAT friction model is a modification of some more classical models -such as Voellmy model :ref:`theoryCom1DFA:Voellmy friction model`. The basal shear stress tensor :math:`\tau^{(b)}` +such as Voellmy model :ref:`voellmyfrict`. The basal shear stress tensor :math:`\tau^{(b)}` is expressed as (:cite:`Sa2007`): .. math:: \tau^{(b)} = \tau_0 + \tan{\delta}\,\left(1+\frac{R_s^0}{R_s^0+R_s}\right)\,\sigma^{(b)} + \frac{\rho_0\,\overline{u}^2}{\left(\frac{1}{\kappa}\,\ln\frac{\overline{h}}{R} + B\right)^2} + :label: SamosAT friction model With @@ -585,7 +613,7 @@ The default configuration also provides two additional calibrations for small- avalanches. A further constraint is the altitude of runout below 1600m msl for both. Wet snow friction type -"""""""""""""""""""""""" ++++++++++++++++++++++++ .. Note:: @@ -597,6 +625,7 @@ This approach is based on the Voellmy friction model but with an enthalpy depend .. math:: \tau^{(b)} = \mu\,\sigma^{(b)} + c_\text{dyn}\,\rho_0\,\bar{u}^2 + :label: Wet snow friction model where, @@ -604,6 +633,7 @@ where, .. math:: \mu = \mu_0\,\exp(-enthalpy/enthRef) + :label: mu enthalpy The total specific enthalpy of the particles is initialized based on their initial temperature, specific heat capacity, @@ -612,7 +642,147 @@ Throughout the computation, the particles specific enthalpy is then computed fol .. math:: enthalpy = totalEnthalpy - g\,z - 0.5\,\bar{u}^2 + :label: enthalpy + +Rheological Models +""""""""""""""""""" + +.. Note:: This documentation about rheological models is currently under development and therefore not complete yet! + Also the parameters are not calibrated yet! + +General +++++++++ +Debris flows are gravity-driven masses of poorly sorted and water saturated sediments whose dynamics are strongly influenced by +solid and fluid forces :cite:`Iv1997`. Unlike the granular friction models presented in the previous section, which consider the solid +fraction of a material, the rheological models implemented in :py:mod:`com1DFA` are intended to simulate the visco-turbulent behaviour of the solid-fluid mixtures. +A rheological constitutive law describes the shear stress applied at a given shear rate (:numref:`Overview-rheological-models-fig`). +In general, these laws can be subdivided into: Newtonian rheologies, in which shear stress increases linearly with shear rate (e.g. clear water), +and non-Newtonian fluids, in which a yield shear stress has to be exceeded or shear stress increases non-linearly with shear rate (e.g. debris flows, mudflows). + +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 + :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). +The flow index :math:`n` describes the rheological behaviour of the mixture as :cite:`KaZhHaHe2025`: + +- :math:`n = 0` a rate-independent solid-like behaviour, +- :math:`n = 1` a Newtonian fluid-like behaviour, +- :math:`0 < n < 1` a shear-thinning non-Newtonian fluid-like behaviour, +- :math:`1 < n < 2` a shear-thickening non-Newtonian fluid-like behaviour, +- :math:`n = 2` an intertial mixture (fluid or solid). + +: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: + +.. _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` | + +-------------------------------------+-----------------------------+------------------------+------------------------+ + + +:numref:`Overview-rheological-models-fig` illustrates the behaviour of the graphs of the rheological models due to shear rate. + +.. _Overview-rheological-models-fig: + +.. figure:: _static/Overview_rheological_models.png + :width: 90% + + Overview rheological models + + +It is well known that the bulk viscosity :math:`\eta_m` and the yield shear stress :math:`\tau_y` are both a function of the +volumetric sediment concentration :math:`C_v`. The dependencies are expressed by following equations :cite:`ObJu1985, ObJu1993`: + +.. math:: + \eta_m = \alpha_1 \cdot e^{\beta_1 \cdot C_v} + :label: eta-cv + +.. math:: + \tau_y = \alpha_2 \cdot e^{\beta_2 \cdot C_v} + :label: tauy-cv + +where :math:`\alpha_1` and :math:`\beta_1`, :math:`\alpha_2` and :math:`\beta_2`, respectively, are empirical coefficients which +are determined in lab experiments. + +Since the governing equations to be solved are depth-averaged, the shear rate :math:`\dot\gamma = \frac{du}{dz}` +cannot be computed directly. Assuming a parabolic vertical flow velocity distribution the integration over the flow thickness results +in following substitution :cite:`Iv1997, DeIv2001, GiFlSaHe2020`: + +.. math:: + \dot{\gamma} = \frac{3 \cdot \overline{u}}{\overline{h}} + :label: shearRate_substitution + +where :math:`\overline{u}` is the depth-averaged flow velocity and :math:`\overline{h}` is the flow thickness. + +O´Brien and Julien ++++++++++++++++++++ +The quadratic rheological model proposed by O´Brien and Julien :cite:`ObJu1985` accounts for various shear stress components, including +cohesive yield stress, viscous stress, turbulent stress, and dispersive stress, which arise from sediment particle +collisions under high deformation rates. +The resulting shear stress reads: + +.. math:: + \tau^{(b)} = \tau_y + \eta_m \cdot \dot{\gamma} + C \cdot \dot{\gamma}^2 + :label: oBrienAndJulien + +The turbulence and dispersive shear stresses, which are both functions of the second power of the shear rate, +are taken into account by the factor :math:`C` which incorporates : + +.. math:: + C = \rho_m \cdot l_m^2 + \alpha_i \cdot \rho_s \cdot \lambda^2 \cdot d_s^2 + :label: dispersiveShearStress + +.. math:: + \frac{1}{\lambda} = \left(\frac{C_m}{C_v}\right)^{1/3} - 1 + :label: sedimentConcentration + +where :math:`\rho_m` is the mass density of the solid-fluid mixture :math:`[kg \cdot m^{-3}]`, :math:`l_m` is the Prandtl mixing length :math:`[m]` +(approximate :math:`0.4 \cdot h`), :math:`\alpha_i` is a coefficient (= 0.01, :cite:`Ta1978`), :math:`\rho_s` is the mass density +of sediment :math:`[kg \cdot m^{-3}]`, :math:`d_s` is the sediment size :math:`[m]`, :math:`\lambda` is the linear sediment concentration :cite:`Ba1954`, +:math:`C_m` is the maximum concentration of sediment particles (= 0.615, :cite:`Ba1954`) and :math:`C_v` is the volumetric sediment concentration. + +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. + +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 + :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. + +Bingham +++++++++++++++++ +After exceeding a threshold shear stress the rheological model proposed by Bingham :cite:`Bi1919` describes a linear relationship +between shear stress and shear rate. The equation reads: + +.. math:: + \tau^{(b)} = \tau_y + \eta_m \cdot \dot{\gamma} + :label: bingham +The Bingham model is well-suited to homogeneous suspensions of fine particles at low shear rates (e.g. mudflows, :cite:`Ju2010`). Dam ~~~