calcUstar.py

import numpy as np
from .uniqueFinder import uniqueFinder
from .Rxy import Rxy
from .outlierRemove import outlierRemove
from .velocityProfile import velocityProfile


def calcUstar(data, Hs=None, bottomWall=None, topWall=None, vonKarman=None):
    '''
    This function is to calculate the friction velocity using output of Rxy(data).
    Inspired by Nezu's book on Turbulence. This method should not be used for finding the shear stress at the carpet. An alternative way is to find the slope using log-wall by keeping kappa = 0.412.
    Input : measurement data

    Output : U* - friction velocity
             hs - water depth at which the shear stress is zero (R. N. Parthasarathy and M. Muste (1995))
             topWall - If it is True, Method of (R. N. Parthasarathy and M. Muste (1995)) is used here to calcuate the friction velocity at the carpet.
    '''

    if vonKarman==None:
        data = Rxy(data)

        dataWithoutOutliers = outlierRemove(data,'Rxy')

        # Fitting line
        slope, intercept = np.polyfit(dataWithoutOutliers['depth[cm]'], dataWithoutOutliers['uv_prime'], 1)

        Ustar = np.sqrt(intercept) ## Ref. paper of R. N. Parthasarathy (1995)
        hs  = -1*intercept/slope  ## hs = water depth at which the shear stress is zero.
        #y=np.linspace(0,data[-2])
        #Ustar = np.sqrt(np.abs(slope)*data[-2])
        print('U* = ', Ustar)

    if topWall:  ## Method of (R. N. Parthasarathy and M. Muste (1995)) is used here to calcuate the friction velocity at the carpet.
        depthUnderCarpet = data[1]-0.02 # average water depth - particle's diameter
        ## See eq.3 of paper of R. N. Parthasarathy (1995)
        Ustar_topWall = np.sqrt(Ustar**2 * depthUnderCarpet/hs) ## friction velocity at the top wall
        print("U* at wall", Ustar_topWall)
        return Ustar_topWall

    if bottomWall == True and vonKarman == True:
        U = velocityProfile(data) # Use velocityProfile function to calcuate the velocity profile in the vertical direction
        cutoffPoint = np.argmax(U[0].reset_index()['U']) ## To find the depth of maximum velocity
        logZ = np.log(U[0].reset_index()['depth[cm]'][:cutoffPoint]/100/1e-6) #  log(y/nu)
        slope, intercept = np.polyfit(logZ, U[0].reset_index()['U'][:cutoffPoint], 1) ## Use polyfit of numpy to find the best fitting line
        kappa = 0.412 # von Karman's constant
        Ustar = kappa*slope # m/sec
        A = (intercept-Ustar/kappa*np.log(Ustar))/Ustar

        return Ustar, A


    if vonKarman: ## To calcuate shear stress at the carpet by using log-wall law by using von Karman's constant = 0.412. We used fitting line method to find out the shear stress.
        U = velocityProfile(data) # Use velocityProfile function to calcuate the velocity profile in the vertical direction
        ks = 0.01 # half of the diameter 0.02/2 [meter]
        cutoffPoint = np.argmax(U[0].reset_index()['U'])-1  ## To find the depth of maximum velocity, we first find the index of maximum velocity. It is subtracted with -1 for later purpose.
        depthUnderCarpet = U[1]-1 ## Water depth under the carpet = average water depth - radius of the particle (1 cm). My vertical coordinate is pointing downwards from the half of the carpet thickness.
        newDepth = depthUnderCarpet-U[0].reset_index()['depth[cm]'][:cutoffPoint:-1] ## Since we changed vertical refernce line, we calculated the new water levels up to the point of max; vel.
        logZ = np.log(newDepth)
        slope, intercept = np.polyfit(logZ, U[0].reset_index()['U'][:cutoffPoint:-1], 1) ## Use polyfit of numpy to find the best fitting line
        kappa = 0.412 # von Karman's constant
        Ustar = kappa*slope # shear velocity
        y0 = np.exp(-intercept*kappa/Ustar) # roughness characteristic height
        B = 1/kappa*np.log(ks/(y0/100)) # 0.02 is diameter of the particle. It is 'k' roughness height.

        return Ustar, B

    if Hs:
        return Ustar, hs

    return Ustar