CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

| Download

Computational notebook for computing ks from Cf using the method in Monty et al. (2016) to be used alongside the manuscript "Effects of fetch length on turbulent boundary layer recovery past a step-change in surface roughness"

Views: 18
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Tweet
Share
Share
Kernel: Python 3 (system-wide)

This code was made in support of the article "Effects of fetch length on turbulent boundary layer recovery past a step-change in surface roughness" by Martina Formichetti1, Dea D. Wangsawijaya, Sean Symon, and Bharathram Ganapathisubramani from the Department of Aeronautical and Astronautical Engineering, University of Southampton, University Road, Southampton, SO17 1BJ, UK.

It provides a tool to calculate the equivalent sandgrain roughness height, ksk_s, given a known friction coefficient, CfC_f, in lign with the method adopted by the authors of Monty, J. P. et al. (2016) "An assessment of the ship drag penalty arising from light calcareous tubeworm fouling" Biofouling 32(4), 451–464.

# load necessary modules
import numpy as np
from scipy.optimize import minimize
from IPython.display import Latex, display

First we load the experimental data and define the necessary constants:

# Smooth wall wss (measured in BLWT)
nu = 1.5259e-05
B = 4.3474
PI = 0.3386
kappa = 0.39

# Roughness wss (P24 grit sandpaper)
# Roughness height k = 6σ, per Gul & Ganapathisubramani (2022)
# k = 1.626 mm, ks = 2.3312 mm (from the highest Re)
# From Alicona scanner, k = 1.695 mm for the current P24 batch
ks0 = 2.3312E-3  # in m
BFR = 8.5
epsilon = 1e-9 #stability factor avoids division by 0
x = 6.7 #effective fetch of the longest fetch case i.e. fully rough homogeneous roughness

# load your data or sub it below where u0=mean streamwise velocity,rho=density,
# nui=kinematic viscosity,q=dynamic pressure,utau=friction velocity,Cf=friction coefficient
u0 = [10.3105047395910, 15.3773793135073, 20.5083531748649, 25.6487506373946, 30.7611612595695, 35.9244836186404, 41.1042657303576]
rho = [1.19094491559810, 1.19372394996242, 1.19400529473343, 1.19441581902224, 1.19415794261944, 1.19408836992398, 1.19315583861218]
nui = [1.52975528487435e-05, 1.52332457531582e-05, 1.52268440926290e-05, 1.52169861538793e-05, 1.52216476130015e-05, 1.52231857218278e-05, 1.52437306108320e-05]
q = [63.3024792529917, 141.136513783458, 251.095123115035, 392.878233846369, 564.985348813404, 770.526395903804, 1007.95455643058]
utau = [0.558750450525668, 0.840222031932186, 1.12506801135256, 1.40885978604175, 1.68977813566214, 1.97230432501546, 2.25588902949122]
Cf = [0.00587363934676032, 0.00597106015366755, 0.00601902874984967, 0.00603438632339124, 0.00603509236308117, 0.00602832676915149, 0.00602409278258317]
wss = []
# Append the mean dictionaries to the wss list
wss.append({
    'u0': np.array(u0),
    'rho': np.array(rho),
    'nu': np.array(nui),
    'q': np.array(q),
    'utau': np.array(utau),
    'Cf': np.array(Cf),
})
wss = wss[0]
L = 1 #in this example I only used one fetch length = 1 delta

Create arrays for the lines of constant length, ks/xk_s/x, and the lines of constant Reynolds, U∞/νU_\infty/\nu 

# For plotting Cf at constant Uinf/nu and constant ks/x
Re_const = {
    'ks0': ks0,
    'nu': nu,
    'u0': np.array([10, 130]),  # the range of measurements is between 10-40 mps
}
Re_const['Re_unit'] = Re_const['u0'] / Re_const['nu']

# Constant L (or ks(x))
npoints = 2 #keeping the number of points minimal to avoid memory problems, do increase these for better accuracy
# L_const.L = [1.4 8.8]; #length in m
L_const = {
    'L': np.logspace(np.log10(1.4), np.log10(8.8), npoints),
}
L_const['ks_on_x'] = ks0 / L_const['L']

Define a function to model a TBL mean velocity profile in inner units, this one is based on the wake function of Jones, et al (2001)

def wake_jones(kappa, B, PI, delta_plus):

    yplus = np.arange(delta_plus + 1)
    eta = yplus / delta_plus

    uplus = np.zeros_like(yplus, dtype=float)
    mask = yplus > 0
    uplus[mask] = 1 / kappa * np.log(yplus[mask]) + B - 1 / (3 * kappa) * eta[mask] ** 3 + 2 * PI / kappa * eta[mask] ** 2 * (3 - 2 * eta[mask])
    uplus[0] = 0  # force wall to zero

    S = 1 / kappa * np.log(delta_plus) + B - 1 / (3 * kappa) + 2 * PI / kappa

    return yplus, uplus, S

Define a function for the iterative process to obtain as a function of and described in Monty, et al, 2016. 

def solve_Cf_rough(Cf_guess, kappa, delta_plus, ks0, Re_unit, BFR, PI):
    rhs = 1 / kappa * np.log(delta_plus) - 1 / kappa * np.log(ks0 * Re_unit) + BFR - 1 / (3 * kappa) + 2 * PI / kappa
    lhs = np.sqrt(2 / Cf_guess) + 1 / kappa * np.log(np.sqrt(Cf_guess / 2))
    err = np.abs(lhs - rhs)
    return err

Define the bounds for CfC_f, RexRe_x, and ΔU+\Delta U^+

# delta_plus (Re_tau)
lb = 1E2
ub = 1E8
delta_plus = np.logspace(np.log10(lb), np.log10(ub), npoints)
Rex_lim = [5E5, 5E7]
Cf_lim = [2.5E-3, 7E-3]

Use the wake function to obtain the mean velocity profile in inner units of a rough and a smooth wall and the wall normal viscous coordinate. Use these to obtain ReθRe_\theta, RexRe_x and CfC_f for a smooth wall.

# find C_f vs Re_theta
Re_theta = np.zeros_like(delta_plus)
Cf = np.zeros_like(delta_plus)
for j in range(len(delta_plus)):
    yplus, uplus, S = wake_jones(kappa, B, PI, delta_plus[j])
    Re_theta[j] = np.trapz(uplus - uplus ** 2 / S, yplus)
    Cf[j] = 2 / S ** 2

# find C_f vs Re_x
Rex = np.zeros_like(delta_plus)
for i in range(1, len(delta_plus)):
    Rex[i] = np.trapz(2 / Cf[:i + 1], Re_theta[:i + 1])

Re_const['Cf'] = np.zeros((len(delta_plus), len(Re_const['Re_unit'])))
Re_const['Rex'] = np.zeros((len(delta_plus), len(Re_const['Re_unit'])))

Use the function adopting the method in Monty et al. 2016 to obtain CfC_f for constant RexRe_x for the fully rough flow. A small epsilon was added to avoid division by zero and log of zero issues.

for ii in range(len(Re_const['Re_unit'])):
    # find C_f vs Re_theta
    for j in range(len(delta_plus)):
        # solve for Cf
        optfun = lambda Cf_guess: solve_Cf_rough(Cf_guess, kappa, delta_plus[j], Re_const['ks0'],
                                                 Re_const['Re_unit'][ii], BFR, PI)
        #keeping iterations low to avoid memory ptoblems, increase for better accuracy
        res = minimize(optfun, 0.001, bounds=[(epsilon, 0.01)], method='SLSQP', options={'ftol': 1e-4, 'maxiter': 1})
        Re_const['Cf'][j, ii] = res.x.item()

        # solve for theta
        ks_plus0 = Re_const['ks0'] * Re_const['Re_unit'][ii] * np.sqrt(Re_const['Cf'][j, ii] / 2)
        deltau_plus = 1 / kappa * np.log(ks_plus0) + B - BFR
        yplus, uplus, S = wake_jones(kappa, B, PI, delta_plus[j])
        uplus -= deltau_plus
        S -= deltau_plus
        Re_theta[j] = np.trapz(uplus - uplus ** 2 / S, yplus)

    # find C_f vs Re_x
    for i in range(1, len(delta_plus)):
        Re_const['Rex'][i, ii] = np.trapz(2 / Re_const['Cf'][:i + 1, ii], Re_theta[:i + 1])

L_const['Rex'] = np.zeros((len(Re_const['Re_unit']), len(L_const['L'])))
L_const['Cf'] = np.zeros((len(Re_const['Re_unit']), len(L_const['L'])))

Compute CfC_f for the fully rough flow for constant L=ks/xL=k_s/x

for j in range(len(L_const['L'])):
    # lookup table
    for ii in range(len(Re_const['Re_unit'])):
        L_const['Rex'][ii, j] = L_const['L'][j] * Re_const['Re_unit'][ii]
        L_const['Cf'][ii, j] = np.interp(L_const['Rex'][ii, j], Re_const['Rex'][:, ii], Re_const['Cf'][:, ii])
--------------------------------------------------------------------------- NameError Traceback (most recent call last) /tmp/ipykernel_1003/1308400285.py in <cell line: 1>() ----> 1 for j in range(len(L_const['L'])): 2 # lookup table 3 for ii in range(len(Re_const['Re_unit'])): 4 L_const['Rex'][ii, j] = L_const['L'][j] * Re_const['Re_unit'][ii] 5 L_const['Cf'][ii, j] = np.interp(L_const['Rex'][ii, j], Re_const['Rex'][:, ii], Re_const['Cf'][:, ii]) NameError: name 'L_const' is not defined

Compute ksk_s,ks+k_s^+, and ΔU+\Delta U^+

delta99 = 149.5
fetch = (1 * 300 - 150) / delta99

u0_lim = [20, 40]  # only use the Cf between 20-40 mps

ks = []
ks_plus = []
deltau_plus = []

u0 = wss['u0']
Cf = wss['Cf']
utau = wss['utau']
nu = wss['nu']

# Find indices based on the limits
id1 = np.argmin(np.abs(u0 - u0_lim[0]))
id2 = np.argmin(np.abs(u0 - u0_lim[1]))

# Compute the mean of Cf in the range [id1, id2]
Cf_mean = np.mean(Cf[id1:id2])
idopt = np.argmin(np.abs(np.mean(L_const['Cf'], axis=0) - Cf_mean))
ks = L_const['ks_on_x'][idopt] * x * 1E3
ks_plus = ks * 1E-3 * utau / nu
ks_plus = ks_plus[0]
deltau_plus = 1 / kappa * np.log(ks_plus) + B - BFR
display(Latex(r'Fetch L = 1$\delta_2$:'))
display(Latex(f'1. $k_s$ [mm]: {ks}'))
display(Latex(f'2. $k_s^+$ [-]: {ks_plus}'))
display(Latex(f'3. $\Delta U^+$ [-]: {deltau_plus}'))
display('-' * 50)