Figure 9.ipynb

Kernel: Python 3 (system-wide)

Model for predicting skin friction from effects of pressure gradient histories on skin-friction and mean flow of high Reynolds number turbulent boundary layers over smooth and rough walls

This code presents the work given in the paper to calulate the skin friction assuming the ZPG conditions along with the pressure gradient history. The model is based on the $\Delta \beta$ parameter which is defined as the following -

$\Delta \beta = \left(\frac{\delta^*}{\tau_w}\right)_{DS} \left[\frac{1}{x_{DS}-x_{US}}\int_{x_{US}}^{x_{DS}}\left(\frac{dP}{dx}\right) w(x) dx\right]~~{\rm where}~w(x) = \frac{x-x_{US}}{x_{DS}-x_{US}}$

It is based on the weighted average of the pressure gradient and the local conditions at the measurement station. Furthermore, we can predict skin friction based on the following equation -

$\sqrt{\frac{2}{C_f^{PG}}} - \sqrt{\frac{2}{C_f^{ZPG}}} = \frac{1}{\kappa} ln \left(\frac{y_0^{ZPG}}{y_0^{PG}}\right) + \frac{2}{\kappa} (\Pi^{PG} - \Pi^{ZPG}) + \frac{1}{\kappa} ln\left(\frac{\delta_{PG}^+}{\delta_{ZPG}^+}\right)$

We can assume the last two terms to be negligible at matched $Re_\tau$ since this work found $y_0$ to be independent of the pressure gradient history. This model uses a calibration based on the smooth wall data to predict the rough wall data, this function is based on the relationship between $\Pi_{PG} - \Pi_{ZPG}$ and $\Delta \beta$. An iterative function is then used to solve for the $C_f$, $\Pi$ and $\delta^*$.

This first cell imports the required functions for the model

import numpy as np  # Numerical computing
import matplotlib.pyplot as plt  # Plotting
from scipy.optimize import curve_fit, fmin  # Curve fitting and optimization
import pickle  # Object serialization
from scipy.special import erf  # Error function

This next cell imports the date required for the model as well as defines the velocity and angle of attack arrays along with the symbols used for plotting.

# Import the data required for the model
with open('Cf_predictive_model_data.pkl', 'rb') as f:
    Model_Data = pickle.load(f)  # Load model data from a pickle file

# Define arrays of variables
AOA_array = [-8, -4, 0, 4, 8]  # Angles of Attack 500mm (AOA)
AOA_array_400mm_RW = [-10, -8, -4]  # AOA for 400mm Right Wing (RW)
AOA_array_400mm_SW = [-10, -8]  # AOA for 400mm Side Wing (SW)

Vel_array_400mm_RW = [20, 25, 30]  # Velocity array for 400mm RW
Vel_array_400mm_SW = [10, 20, 30]  # Velocity array for 400mm SW
Vel_array_RW = [10, 15, 20, 25, 30]  # Velocity array for RW
Vel_array_zpg = [15, 20, 25, 30, 35, 40]  # Velocity array for zero pressure gradient (zpg)
Vel_array_SW_all = [10, 15, 20, 25, 30]  # Velocity array for all SW
Vel_array_SW = [10, 20, 30]  # Velocity array for SW

Vel_array_upstream_SW = [10, 20, 30]  # Velocity array for upstream SW

# Symbols and colors for plotting
symbol = ['s', 'v', 'D', '^', 'H']  # Symbols for 500mm cases
symbol_400mm = ['>', 'd', '<']  # Symbols for 400mm cases
symbol_400mm_SW = ['>', 'd']  # Symbols for 400mm SW cases

colors_SW = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']  # Colors for SW

colors_SW_400mm = ['tab:brown', 'tab:cyan']  # Colors for 400mm SW
colors_RW_400mm = ['tab:brown', 'tab:cyan', 'tab:pink']  # Colors for 400mm RW

kappa = 0.39  # Constant value for kappa

This cell defines some common functions required for the code. The first is the functional fit assumed between the difference in PI and the $\Delta \beta$. The next function defines the fit used for the wake fit based on Lewkowicz, A. 1982. The next two define the musker function for which we use the average value of a found from the fitting of the 500mm pressure gradient cases. For a smooth wall, the velocity profile is defined as using fit_composite_yplus, which combines Musker with the normal outer region func. Finally, a weight function is defined to be linear; however, it can be varied to use the error function if required.

def deltabeta_func(deltabeta, A):
    return A * deltabeta  # Linear relationship between deltabeta and A used for curve fitting

def W_func(eta):
    term1 = 2 * eta**2 * (3 - 2 * eta)
    term2 = (1 / np.pi) * eta**2 * (1 - eta) * (1 - 2 * eta)
    return term1 - term2  # Computes the weighting based on eta which is y/delta

def musker_func(yplus, a):
    alpha = (-1 / kappa - a) / 2
    beta = np.sqrt(-2 * a * alpha - alpha**2)
    R = np.sqrt(alpha**2 + beta**2)
    term1 = (1 / kappa) * np.log((yplus - a) / (-a))
    term2a = R**2 / (a * (4 * alpha - a))
    term2b = (4 * alpha + a) * np.log(-(a / R) * (np.sqrt((yplus - alpha)**2 + beta**2) / (yplus - a)))
    term2c = (alpha / beta) * (4 * alpha + 5 * a) * (np.arctan((yplus - alpha) / beta) + np.arctan(alpha / beta))
    term2 = term2a * (term2b + term2c)
    Uplusinner = term1 + term2
    return Uplusinner  # Computes inner layer (Musker profile) for turbulent flow over a smooth wall

def fit_composite_yplus(yplus, a, PI, retau):
    musk = musker_func(yplus, a) + np.exp(-np.log(yplus / 30)**2) / 2.85
    eta = yplus / retau
    W = 2 * eta**2 * (3 - 2 * eta) - (1 / np.pi) * eta**2 * (1 - eta) * (1 - 2 * eta)
    outer = PI / kappa * W
    return musk + outer  # Combines Musker and outer-layer profile

downstream_hw = 9.03 # defines where pressure gradient history for the integral ends
upstream_reference = 9.03 - 1.25 * 3  # defines where pressure gradient history starts from 

def weight_func(x):
    # Linear weight
    weight = (x - upstream_reference) / (downstream_hw - upstream_reference)

    # Error function-based weighting (commented out alternative)
    # mean, sigma = 0.5, 0.15
    # xnorm = (x - x[0]) / (x[-1] - x[0])
    # weight = 0.5 * (1 + erf((xnorm - mean) / (sigma * np.sqrt(2))))

    return weight  # Computes weight for the intergral of the history, it can be changed to use ERF function to show effect of weighting

The following cell produces a calibration function using the smooth wall data between the $\Pi_{PG} - \Pi_{ZPG}$ and $\Delta \beta$ and then plots the result. We calibrate the function by calculating $\delta^*$ from the composite velocity profile.

downstream_hw = 9.03  # Downstream reference limit where local conditions are to be evalualted
upstream_reference = 9.03 - 1.25 * 3  # Upstream reference limit

Vel = str(20)  # Velocity key as string

# Create figure and axes for plotting using a subplot mosaic
fig, axs = plt.subplot_mosaic([['a)']])

def beta_func(beta, A):
    return A * beta  # Defines a linear relationship for beta

# Initialize arrays for storing results
beta_all = []
pi_all = []

# Loop through AOA_array to process 500mm smooth wall (SW) data
for j in range(len(AOA_array)):
    AOA = str(AOA_array[j])
    Key_AOA = f"{AOA}_500mm_SW"  # Construct key for 500mm SW data
    Key_ZPG = 'ZPG_SW'  # Key for Zero-Pressure-Gradient (ZPG) SW data

    Tap_positons_m = Model_Data[Vel][Key_AOA]['Pressure Taps (m)']
    indexes_fit = np.where((Tap_positons_m < downstream_hw) & (Tap_positons_m > upstream_reference))  # Select fitting range

    weight = weight_func(Tap_positons_m[indexes_fit])  # Compute weight
    grad = Model_Data[Vel][Key_AOA]['ddeltaP/dx'][indexes_fit]  # Pressure gradient

    # Calculate difference in PI between AOA and ZPG
    delta_PI = Model_Data[Vel][Key_AOA]['PI'] - Model_Data[Vel][Key_ZPG]['PI']
    tau0 = Model_Data[Vel][Key_AOA]['rho'] * Model_Data[Vel][Key_AOA]['utau']**2  # Compute tau0

    # Compute composite profile for U+ and boundary layer integral parameters
    zplus = np.logspace(np.log10(5), np.log10(Model_Data[Vel][Key_AOA]['Re_tau']), 100)
    Uplus_profile = fit_composite_yplus(zplus, -10.36590192, Model_Data[Vel][Key_AOA]['PI'], Model_Data[Vel][Key_AOA]['Re_tau'])
    Cf = Model_Data[Vel][Key_AOA]['Cf']
    zfit = zplus * Model_Data[Vel][Key_AOA]['nu'] / Model_Data[Vel][Key_AOA]['utau']
    deltastar = np.trapz((1 - (Uplus_profile / np.sqrt(2 / Cf))), zfit) + 0.5 * (1 - (Uplus_profile[0] / np.sqrt(2 / Cf))) * zfit[0]

    # Calculate Δβ integral using weighted gradient
    Delta_beta_int = (deltastar / tau0) * (1 / (Tap_positons_m[indexes_fit][-1] - Tap_positons_m[indexes_fit][0])) * \
        np.trapz(grad * weight, Tap_positons_m[indexes_fit])
    
    # Plot results
    axs['a)'].plot(Delta_beta_int, delta_PI, marker=symbol[j], color=colors_SW[j], markerfacecolor='none', linestyle='none')

    # Store results in arrays
    beta_all.append(Delta_beta_int)
    pi_all.append(delta_PI)

# Plot reference point at (0, 0) for ZPG results
axs['a)'].plot(0, 0, marker='o', color='k', markerfacecolor='none', linestyle='none')

# Loop through AOA_array_400mm_SW to process 400mm smooth wall (SW) data
for j in range(len(AOA_array_400mm_SW)):
    AOA = str(AOA_array_400mm_SW[j])
    Key_AOA = f"{AOA}_400mm_SW"  # Construct key for 400mm SW data
    Key_ZPG = 'ZPG_SW'  # Key for Zero-Pressure-Gradient (ZPG) SW data

    Tap_positons_m = Model_Data[Vel][Key_AOA]['Pressure Taps (m)']
    indexes_fit = np.where((Tap_positons_m < downstream_hw) & (Tap_positons_m > upstream_reference))  # Select fitting range

    weight = weight_func(Tap_positons_m[indexes_fit])  # Compute weight
    grad = Model_Data[Vel][Key_AOA]['ddeltaP/dx'][indexes_fit]  # Pressure gradient

    delta_PI = Model_Data[Vel][Key_AOA]['PI'] - Model_Data[Vel][Key_ZPG]['PI']
    tau0 = Model_Data[Vel][Key_AOA]['rho'] * Model_Data[Vel][Key_AOA]['utau']**2  # Compute tau0

    # Compute composite profile for U+ and boundary layer integral parameters
    zplus = np.logspace(np.log10(5), np.log10(Model_Data[Vel][Key_AOA]['Re_tau']), 100)
    Uplus_profile = fit_composite_yplus(zplus, -10.36590192, Model_Data[Vel][Key_AOA]['PI'], Model_Data[Vel][Key_AOA]['Re_tau'])
    Cf = Model_Data[Vel][Key_AOA]['Cf']
    zfit = zplus * Model_Data[Vel][Key_AOA]['nu'] / Model_Data[Vel][Key_AOA]['utau']
    deltastar = np.trapz((1 - (Uplus_profile / np.sqrt(2 / Cf))), zfit) + 0.5 * (1 - (Uplus_profile[0] / np.sqrt(2 / Cf))) * zfit[0]

    # Calculate Δβ integral using weighted gradient
    Delta_beta_int = (deltastar / tau0) * (1 / (Tap_positons_m[indexes_fit][-1] - Tap_positons_m[indexes_fit][0])) * \
        np.trapz(grad * weight, Tap_positons_m[indexes_fit])

    # Plot results
    axs['a)'].plot(Delta_beta_int, delta_PI, marker=symbol_400mm[j], color=colors_SW_400mm[j], markerfacecolor='none', linestyle='none')

    # Store results in arrays
    beta_all.append(Delta_beta_int)
    pi_all.append(delta_PI)

# Fit Δβ func to data Δβ and ΔPI using curve fitting
pi_all = np.array(pi_all)
beta_all = np.array(beta_all)
fit_SW_only, _ = curve_fit(beta_func, beta_all, pi_all)

def beta_to_pi_diff(beta):
    return beta_func(beta, *fit_SW_only)  # Function to predict ΔPI from Δβ

# Generate fitted curve
beta_fit = np.linspace(np.min(beta_all) - 0.2, np.max(beta_all) + 0.2, 20)
axs['a)'].plot(beta_fit, beta_to_pi_diff(beta_fit), 'k--')  # Plot fitted curve

# Configure plot appearance
axs['a)'].grid()
axs['a)'].set_xlabel(r'$\Delta \beta$')
axs['a)'].set_ylabel(r'$\Pi_{PG} - \Pi_{ZPG}$')
axs['a)'].set_aspect(1.0 / axs['a)'].get_data_ratio(), adjustable='box')
axs['a)'].set_title(r'Smooth Wall Data Only')
plt.show()

# Print fit parameters for SW data only
print(fit_SW_only)

[0.94422564]