In [ ]:

# script to evaluate the raw data provided in zenodo (https://zenodo.org/record/5222697)
# unpack the zip archive to a folder and set the following path accordingly for the respective case:
dataFolder = Path('rheology_simulation_data/mono')

# then execute all cells of this notebook, which will
# - read the vertical velocity profiles (-> velocity profile over z)
# - read the particle data and extract relevant and horizontally averaged vertical profiles (-> solid volume fraction profile over z)
# - apply averaging in time over the given time span. This time span should be adjusted for each case to just include the quasi-static regime.
# - apply a vertical running average to coarse-grain the profiles
# - evaluate derived profiles like macroscopic friction coefficient, particle pressure, viscous number
# - plot all these profiles

# further processing on the data can then be done to assess the rheology
# more infos on that can be found in 
# Rettinger, Eibl, Rüde, Vowinckel - "Rheology of mobile sediment beds in laminar shear flow: effects of creep and polydispersity"
# see either the preprint on arXiv (https://arxiv.org/abs/2103.04103) 
# or the published version in the Journal of Fluid Mechanics
# note that in the original work, data of z < 5 * averageDiameter is discarded due to particle layering effects

# Author of this script: Christoph Rettinger (christoph.rettinger@fau.de), 2021

In [ ]:

# general setup parameters
Lx = 1024.
Ly = 512.
Lz = 480.

densityRatio = 1.5
uTop = 0.03
Rep = 4.
Sh = 0.5
hbRef=340

runningAvgFilterWidth = 11

In [ ]:

# read fluid profile data
fluidProfileFolder = dataFolder / 'fluid'
fluidProfileFiles = list(fluidProfileFolder.glob('*0.txt'))

fluidHeights = np.loadtxt(fluidProfileFolder / 'height.txt')

def readFluidData(files):
    timesteps = []
    data = []
    for file in files:
        fileName = file.name
        timestep = int(fileName[0:fileName.find(".txt")])
        timesteps.append(timestep)

        profileData = np.transpose(np.loadtxt(file, usecols=[0])) # ux
        data.append(profileData)
        
    return np.array(timesteps), np.array(data)

fluidTimesteps, fluidData = readFluidData(fluidProfileFiles)

In [ ]:

# read and evaluate particle data
particleFolder = dataFolder / 'particle'
particleFiles = list(particleFolder.glob('*0.txt'))

def readAndEvaluateParticleData(files, heightArray, horizontalArea):
    
    def evaluateProfiles(data, heightArray, horizontalArea):
        
        gridSpacing = heightArray[1] - heightArray[0] # dz
        
        fractionProfile = np.zeros(len(heightArray))
        velocityProfile = np.zeros(len(heightArray))

        def sphereShellVolume(r, R, h):
            r1 = max(min(r-h/2,R),-R)
            r2 = max(min(r+h/2,R),-R)
            a1 = np.sqrt(R**2 - r1**2)
            a2 = np.sqrt(R**2 - r2**2)
            hNew = r2-r1

            #see wikipedia "Kugelschicht"
            V = np.pi * hNew / 6 * (3*a1**2 + 3*a2**2 + hNew**2)
            return V

        for j,p in enumerate(data):
            zp = p[0]
            radius = p[1]
            velX = p[2]
            idxBegin = max(int((zp-radius) / gridSpacing - 1),0)
            idxEnd = max(int((zp + radius) / gridSpacing + 2),0)

            for i in range(idxBegin, idxEnd):
                z = heightArray[i]
                weight = sphereShellVolume(z-zp, radius, gridSpacing)
                fractionProfile[i] += weight
                velocityProfile[i] += weight * velX

        # the resulting warning can safely be ignored here
        velocityProfile = np.where(fractionProfile>0, velocityProfile / fractionProfile, np.zeros(len(heightArray)))
        fractionProfile /= (gridSpacing * horizontalArea)
    
        return (fractionProfile,velocityProfile)
        
    
    timesteps = []
    fractionData = []
    velocityData = []
    
    for i,file in enumerate(files):
        fileName = file.name
        timestep = int(fileName[0:fileName.find(".txt")])
        timesteps.append(timestep)

        particleData = np.loadtxt(file, usecols=[2,3,4]) # z, radius, ux 
        if i == 0:
            meanDiameter = np.mean(2. * particleData[:,1])
            
        fractionProfile,velocityProfile = evaluateProfiles(particleData, heightArray, horizontalArea)
        
        fractionData.append(fractionProfile)
        velocityData.append(velocityProfile)
        
    return meanDiameter, np.array(timesteps), np.array(fractionData), np.array(velocityData)

# this step may take a substantial amount of time, in the order of hours
dAvg, particleTimesteps, particleFractionData, particleVelocityData = readAndEvaluateParticleData(particleFiles, fluidHeights, Lx * Ly)

In [ ]:

# definition of further reference quantities
viscosity = uTop * dAvg / Rep
wallShearStressRef = viscosity * uTop / (Lz - hbRef)
gravitationalAcceleration = wallShearStressRef / (Sh * (densityRatio - 1)*dAvg)
normalizeShearStress = (densityRatio - 1) * dAvg * gravitationalAcceleration
tref = dAvg / uTop

In [ ]:

# apply time averaging

# this value has to be a valid timestep (= file name)
beginTemporalAveraging = 1500000 # mono
#beginTemporalAveraging = 5000000 # poly10
#beginTemporalAveraging = 5000000 # poly50
#beginTemporalAveraging = 5000000 # poly100

print("Applying temporal averaging for t in [{},{}] = [{:.1f},{:.1f}] * tref".format(beginTemporalAveraging, particleTimesteps[-1], beginTemporalAveraging/tref, particleTimesteps[-1]/tref ))

beginIdx = np.where(particleTimesteps==beginTemporalAveraging)[0][0]

phiTempAvg = np.mean(particleFractionData[beginIdx:-1], axis=0)
velTempAvg = np.mean(fluidData[beginIdx:-1], axis=0)

In [ ]:

# apply running average for coarse graining

def applyRunningAverage(data, averagingWindowSize):
    # https://stackoverflow.com/questions/13728392/moving-average-or-running-mean/43200476#43200476

    if averagingWindowSize > 1:
        dataAvg = uniform_filter1d(data, size=averagingWindowSize, mode='nearest')
    else:
        dataAvg = data

    # override border values since filter is not able to extrapolate
    dataAvg[0:averagingWindowSize//2] = data[0:averagingWindowSize//2]
    dataAvg[-1-averagingWindowSize//2:] = data[-1-averagingWindowSize//2:]

    return dataAvg

phi = applyRunningAverage(phiTempAvg, runningAvgFilterWidth)
vel = applyRunningAverage(velTempAvg, runningAvgFilterWidth)

In [ ]:

# evaluate the remaining profiles and plot all profiles as a function of normalized height (z)

def plotAllProfiles(height, phiProfile, velocityProfile, 
                    densityRatio, gravitationalAcceleration, viscosity, dAvg):
    
    def getTotalPressure(phi):
        return np.sum(phi) * (densityRatio - 1)*gravitationalAcceleration
    
    def getConfinementPressureProfile(phi, gridSpacing):
        Pp = np.cumsum(phi[::-1])[::-1]*(densityRatio - 1)*gravitationalAcceleration*gridSpacing
        return Pp
        
    
    def getPlotProfilesFromVelocityProfile(height, phi, flowProfile):
        
        gridSpacing = height[1]-height[0] # = dz
        
        shearRate = np.diff(flowProfile,prepend=0) / gridSpacing # gamma        
        tau = viscosity * np.max(shearRate) # = visc * du/dz at top wall, use max here to avoid problems with previous averaging
        ShActual = tau / ( (densityRatio - 1) * gravitationalAcceleration * dAvg )
        print("Actual Shields number = {:.3f}".format(ShActual))
        
        Pp = getConfinementPressureProfile(phi, gridSpacing)
                
        mu = np.array([ tau / x if x > 1e-15 else np.nan for x in Pp]) # friction coefficient
        
        J = np.array([viscosity * abs(x[0]) / x[1] if x[1] > 1e-15 else np.nan for x in zip(shearRate, Pp)]) # viscous number
        
        return (shearRate, Pp, mu, J, tau)
    
    figWidth = 9
    figHeight= 6
    heightNormalizationFactor = 1. / dAvg
    linestyle='.'
    markersize=3
    
    shearRate, Pp, mu, J, tau = getPlotProfilesFromVelocityProfile(height, phiProfile, velocityProfile)

    # solid volume fraction profile
    plt.figure(figsize=(figWidth,figHeight))
    plt.plot(phi, height*heightNormalizationFactor,linestyle,markersize=markersize)
    plt.xlabel(r"$\langle \phi \rangle_t$")
    plt.ylabel(r"$z/\bar{d}_p$")
    plt.xscale("linear")
            
    # (dimensionless) velocity profile
    plt.figure(figsize=(figWidth,figHeight))
    plt.plot(velocityProfile/uTop, height*heightNormalizationFactor,linestyle,markersize=markersize)
    plt.xlabel(r"$\langle u_f \rangle_{V,t}/U_w$")
    plt.ylabel(r"$z/\bar{d}_p$")
    plt.xscale("log")
    
    # (dimensionless) shear rate profile -> as Stokes number
    plt.figure(figsize=(figWidth,figHeight))
    plt.plot(densityRatio * np.abs(shearRate) * dAvg**2 / viscosity, height*heightNormalizationFactor,linestyle,markersize=markersize)
    plt.xlabel(r"$St = \rho_p\bar{d}_p^2 |\dot{\gamma}| / \eta_f$")
    plt.ylabel(r"$z/\bar{d}_p$")
    plt.xscale("log")
    
    # (dimensionless) pressure profile
    plt.figure(figsize=(figWidth,figHeight))
    Ptot = getTotalPressure(phi)
    plt.plot(Pp/Ptot, height*heightNormalizationFactor,linestyle,markersize=markersize)
    plt.xlabel("$p_p/P_\\mathit{tot}$")
    plt.ylabel(r"$z/\bar{d}_p$")
    plt.xscale("linear")
    plt.gca().ticklabel_format(axis='x',style='sci',scilimits=(0,0))
    
    # mu (friction coefficient) profile
    plt.figure(figsize=(figWidth,figHeight))
    plt.plot(mu, height*heightNormalizationFactor,linestyle,markersize=markersize)
    plt.xlabel("$\mu$")
    plt.ylabel(r"$z/\bar{d}_p$")
    plt.xscale("log")
    
    # J (viscous number) profile
    plt.figure(figsize=(figWidth,figHeight))
    plt.plot(J, height*heightNormalizationFactor,linestyle,markersize=markersize)
    plt.xlabel("$J$")
    plt.ylabel(r"$z/\bar{d}_p$")
    plt.xscale("log")


    plt.show()
    
    return mu, J, tau
    
muProfile, JProfile, tau = plotAllProfiles(fluidHeights, phi, vel, densityRatio, gravitationalAcceleration, viscosity, dAvg)