figure12.ipynb

Kernel: Python 3 (system-wide)

import numpy as np
from scipy.integrate import solve_ivp, trapezoid
from scipy.optimize import brentq
import matplotlib.pyplot as plt
import scienceplots
import pandas as pd
import io
plt.style.use(['science', 'no-latex'])
from cycler import cycler
import matplotlib as mpl

Generalised Actuator Disk (AD) model: Thrust and Power Coefficients

Note that variables are dimensionless: distances are normalised by disk diameter $D$, and velocities are normalised by the incoming velocity $U_0$.

Equation labels in the code refer to the model formulation in the accompanying manuscript.

#--------------------------------
# MAIN FUNCTION
#--------------------------------
def new_AD(E1, E2, I, *,
           a=None, ct=None,
           L_up=5, L_down=30,
           tol=5e-3, r=0.05, max_iter=500,
           n_eval_up=1000, n_eval_down=10000):
    """
Exactly one of (a, ct) must be provided. The code will compute the other missing value (a or ct), as well as upstream and downstream variations of axial velocity U, pressure P, and cross-sectional width of CV, sigma.

    Parameters
    ----------
    E1, E2 : float
        Entrainment coefficients in Eqs. (2.5)–(2.6).
    I : float
        Ambient turbulence intensity in Eq. (2.6).
    a : float, optional (keyword-only)
        Induction factor at the disk; if supplied, Ct is obtained iteratively.
    ct : float, optional (keyword-only)
        Thrust coefficient; if supplied, a is obtained iteratively.
    L_up : float
        Upstream extent (in diameters).
    L_down : float
        Downstream extent (in diameters).
    tol : float
        Convergence tolerance for iteration (|new − old| ≤ tol).
    r : float
        Relaxation factor in the iteration update (0 < r ≤ 1).
    max_iter : int
        Maximum number of iterations permitted for each case.
    n_eval_up : int
        Number of points for the upstream region.
    n_eval_down : int
        Number of points for the downstream region.

    Returns
    -------
    x : ndarray
        Axial coordinates (dimensionless), concatenated as [upstream, downstream].
    sigma : ndarray
        Cross-sectional width of the control voluem σ/D.
    U : ndarray
        Axial velocity U/U0.
    Ue : ndarray
        Entrainment velocity U_e/U0 from Eq. (2.7).
    result_scalar : float
        The computed scalar: `Ct` if `a` was provided, or `a` if `ct` was provided.

    Notes
    -----
    • Arguments after the single '*' are keyword-only (Python syntax).
    """

    # Exactly one of a, ct must be set
    if (a is None) == (ct is None):
        raise ValueError("Provide exactly one of 'a' or 'ct', not both or neither.")

    R_hat = 0.5  # disk radius in units of D

    # ----------------------------------------------------------------------
    # Entrainment model (dimensionless) — Eqs. (2.5)–(2.7)
    # ----------------------------------------------------------------------
    def entrainment_velocity(U, x):
        Ue_w = E1 * (1.0 - U)  # Eq. (2.5)
        Ue_b = E2 * I * x / np.sqrt(x**2 + R_hat**2)   # Eq. (2.6)
        n = 4
        return (Ue_w**n + Ue_b**n)**(1/n) # Eq. (2.7)

    # ----------------------------------------------------------------------
    # Downstream: integrate sigma(x), U(x) given (a, Ct) via Eqs. (2.17a–b)
    # ----------------------------------------------------------------------
    def downstream_sigma_U_given_ct(current_a, ct_val, L=L_down):
        def rhs(x, y):
            sigma, U = y
            Ue = entrainment_velocity(U, x)
                        
            # Eq. (2.17b)
            dUdx = (1.0 / (2.0 * U)) * (
                (R_hat**2) / (x**2 + R_hat**2)**(1.5) * ((2.0 * current_a - current_a**2) - ct_val)
                + (8.0 / sigma) * Ue * (1.0 - U))
            
            # Eq. (2.17a)
            dsigmadx = (1.0 / (2.0 * U)) * (4.0 * Ue - sigma * dUdx)
      
            return [dsigmadx, dUdx]

        y0 = [1.0, 1.0 - current_a]           # sigma(0)=1, U(0)=1−a
        x_span = (0.0, L)                     # downstream domain in diameters
        x_eval = np.linspace(x_span[0], x_span[1], n_eval_down)

        # Runge–Kutta (RK45) integration of Eqs. (2.17a–b)
        sol = solve_ivp(rhs, x_span, y0, t_eval=x_eval, method='RK45')

        x = sol.t
        sigma = sol.y[0]
        U = sol.y[1]
        Ue = entrainment_velocity(U, x)
        return x, sigma, U, Ue

    # ----------------------------------------------------------------------
    # Upstream: closed-form sigma(x), U(x) from Eqs. (2.18a–b)
    # ----------------------------------------------------------------------
    def upstream_sigma_U(current_a, L=L_up):
        x = np.linspace(-L, 0.0, n_eval_up)
        # Eq. (2.18a)
        U = np.sqrt((current_a**2 - 2.0 * current_a)
                    * (x / np.sqrt(x**2 + R_hat**2) + 1.0) + 1.0)
        # Eq. (2.18b)
        sigma = np.sqrt((1.0 - current_a) / U)
        return x, sigma, U

    # ----------------------------------------------------------------------
    # Iteration method to find the Ct-a relation
    # Case A (known a): Eq. (2.23) for Ct (dimensionless).
    # Case B (known Ct): inverted Eq. (2.23) for a (dimensionless).
    # The iteration counter is printed for diagnostics.
    # ----------------------------------------------------------------------
    tol_now = np.inf
    count = 0
    convergence=True
    
    if a is not None:
        ct = 4.0 * a * (3.0 - a) / (3.0 * (1.0 + a)) # Initial guess for Ct (Steiros & Hultmark 2018)
        
        while tol_now > tol:
            x_pos, sigma_pos, U_pos, Ue_pos = downstream_sigma_U_given_ct(a, ct)
            dx = x_pos[1] - x_pos[0]
            
            # Mask to calculate integral Y up to 3D only
            mask = x_pos < 3.0
            integrand = (1.0 - x_pos / np.sqrt(x_pos**2 + R_hat**2)) * np.gradient(sigma_pos**2, dx)
            
            # Eq. (2.24)
            Y = trapezoid(integrand[mask], x_pos[mask])

            # Eq. (2.23)
            ct_new = 2.0 * a + (1.0 / Y - 1.0) * a**2

            # Convergence and relaxed update
            tol_now = abs(ct - ct_new)
            ct = r * ct_new + (1.0 - r) * ct

            count += 1
            if count > max_iter:
                ct = np.nan
                convergence=False
                break
        if convergence:
            print(f"a = {a:.2f} and E1 = {E1:.2f} --> CT={ct:.2f} - converged after {count:.0f} iterations. " )
        else:
            print(f"a = {a:.2f} and E1 = {E1:.2f} --> NOT converged after {count:.0f} iterations." )
        
        result_scalar = ct

    else:
        
        a = 0.5 * ct # Initial guess for a

        while tol_now > tol:
            x_pos, sigma_pos, U_pos, Ue_pos = downstream_sigma_U_given_ct(a, ct)
            dx = x_pos[1] - x_pos[0]
            
            # Mask to calculate integral Y up to 3D only
            mask = x_pos < 3.0
            integrand = (1.0 - x_pos / np.sqrt(x_pos**2 + R_hat**2)) * np.gradient(sigma_pos**2, dx)
            
            # Eq. (2.24)
            Y = trapezoid(integrand[mask], x_pos[mask])

            # Eq. (2.23) for a
            a_new = Y * (-1.0 + np.sqrt((Y + (1.0 - Y) * ct) / Y)) / (1.0 - Y)

            # Convergence and relaxed update
            tol_now = abs(a - a_new)
            a = r * a_new + (1.0 - r) * a

            count += 1
            if count > max_iter:
                a = np.nan
                convergence=False
                break
                
        if convergence:
            print(f"CT = {ct:.2f} and E1 = {E1:.2f} --> a={a:.2f} - converged after {count:.0f} iterations." )
        else:
            print(f"CT = {ct:.2f} and E1 = {E1:.2f} --> NOT converged after {count:.0f} iterations." )                

        result_scalar = a

    # ----------------------------------------------------------------------
    # Assemble upstream + downstream
    # ----------------------------------------------------------------------
    x_neg, sigma_neg, U_neg = upstream_sigma_U(a)
    x = np.concatenate([x_neg, x_pos])
    sigma = np.concatenate([sigma_neg, sigma_pos])
    U = np.concatenate([U_neg, U_pos])
    Ue = np.concatenate([np.zeros_like(x_neg), Ue_pos])

    return x, sigma, U, Ue, result_scalar

########### DEFINE INPUT PARAMETERS & ANALYTICAL MODELS ######################
I = 0.05       # incoming turbulence intensity

A = np.arange(0.01, 0.85, 0.05)  # induction factor values

# CT from analytical models
CT_froude = 4 * A * (1 - A)
CT_Steiros = 4 * A * (3 - A) / (3 * (1 + A))
CT_Glauert = CT_froude.copy()
CT_Glauert[A > 0.4] = 0.889 - (0.0203 - (A[A > 0.4] - 0.143)**2) / 0.6427

########### LOAD EXTERNAL DATA & PROCESS SHADED REGIONS ######################
try:
    mit_ct = pd.read_csv('MIT_ct.csv', header=None, names=['a', 'CT'])
    mit_cp = pd.read_csv('MIT_cp.csv', header=None, names=['a', 'CP'])
    Burton_ct = pd.read_csv('Burton_ct.csv', header=None, names=['a', 'CT'])
    nrel_ct = pd.read_csv('NREL_ct.csv', header=None, names=['a', 'CT'])
    nrel_cp = pd.read_csv('NREL_cp.csv', header=None, names=['a', 'CP'])
except FileNotFoundError:
    print("Warning: One or more reference CSV files not found. They will be skipped in plotting.")
    mit_ct = mit_cp = Burton_ct = nrel_ct = nrel_cp = pd.DataFrame(columns=['a', 'CT', 'CP'])

# --- Process Dehtyriov 2023 Data for Shaded Backgrounds ---
dehtyriov_csv = """0.949859596068689,1.01654846335697,0.999768993531818,2.35933806146572,0.999363982191501,2.99763593380614
0.881331677286963,1.01654846335697,0.917371686407219,2.21749408983451,0.953753705103742,2.87943262411347
0.822956042769197,1.01654846335697,0.838784485965607,2.0709219858156,0.901810250687019,2.74231678486997
0.763311372718436,1.01654846335697,0.724709291860172,1.85342789598108,0.84227058357634,2.57683215130023
0.704935738200669,1.01654846335697,0.648660162484549,1.70685579196217,0.758655242346785,2.35460992907801
0.646560103682903,1.01654846335697,0.573883068725924,1.55555555555555,0.694051433440136,2.17021276595744
0.585646398099146,1.01654846335697,0.487702655674358,1.37588652482269,0.638312872760437,2.01418439716312
0.505700159604468,1.01182033096926,0.398999171976815,1.1725768321513,0.577507170200765,1.84397163120567
0.427037957062797,0.983451536643026,0.317903901309236,0.978723404255319,0.49010272287624,1.59338061465721
0.347139719912157,0.903073286052009,0.244419843755625,0.789598108747045,0.372304424523886,1.2434988179669
0.285010980307448,0.817966903073286,0.160813502778077,0.553191489361702,0.259582268303512,0.893617021276595
0.208970851183833,0.657210401891253,0.0531674886896833,0.203309692671394,0.170929786034008,0.609929078014184
0.145621077390166,0.496453900709219,0,0,0.0873564459804874,0.321513002364066
0.0911575524114673,0.330969267139479,,,0,0
0.0417761697327523,0.156028368794326,,,,
0,0,,,,"""

deh_df = pd.read_csv(io.StringIO(dehtyriov_csv), header=None)

# Extract, drop NaNs, and sort each curve by induction factor (a)
lower_curve = deh_df.iloc[:, [0, 1]].dropna().sort_values(by=0)
far_curve = deh_df.iloc[:, [2, 3]].dropna().sort_values(by=2)
near_curve = deh_df.iloc[:, [4, 5]].dropna().sort_values(by=4)

# Create a dense grid of 'a' for smooth shading up to the plot limit (0.75)
a_grid = np.linspace(0, 0.75, 200)

# Interpolate the CT values onto the common a_grid
Ct_lower = np.interp(a_grid, lower_curve[0], lower_curve[1])
Ct_far = np.interp(a_grid, far_curve[2], far_curve[3])
Ct_near = np.interp(a_grid, near_curve[4], near_curve[5])

# Compute corresponding CP boundaries (CP = CT * (1 - a))
Cp_lower = Ct_lower * (1 - a_grid)
Cp_far = Ct_far * (1 - a_grid)
Cp_near = Ct_near * (1 - a_grid)

############### PLOTTING #######################
fig, axs = plt.subplots(1, 2, figsize=(15, 7), sharex=False)
f1 = 24  

# Subplot 1: Shaded regions for CT
axs[0].fill_between(a_grid, Ct_lower, Ct_far, color='blue', alpha=0.1)
axs[0].fill_between(a_grid, Ct_far, Ct_near, color='red', alpha=0.1)

# Subplot 1: CT vs a
axs[0].plot(A, CT_froude, label="Froude $C_T=4a(1-a)$", color='black', linewidth=3, linestyle='-')
axs[0].plot(A, CT_Steiros, label="Steiros & Hullmark (2018)", color='black', linewidth=1.5, linestyle='-', marker='.', markersize=8)
if not mit_ct.empty:
    axs[0].plot(mit_ct['a'], mit_ct['CT'], label='Liew et al. (2024)', color='black', linewidth=1.5, linestyle='-', marker='x', markersize=8)
if not nrel_ct.empty:
    axs[0].scatter(nrel_ct['a'], nrel_ct['CT'], label='Martinez-Tossas et al. (2020)', color='black', marker='+', s=100, linewidths=1.5)

axs[0].set_xlabel("Induction factor $a$", fontsize=f1)
axs[0].set_ylabel("$C_T$", fontsize=f1)
axs[0].grid(False)

# Subplot 2: Shaded regions for CP
axs[1].fill_between(a_grid, Cp_lower, Cp_far, color='blue', alpha=0.1)
axs[1].fill_between(a_grid, Cp_far, Cp_near, color='red', alpha=0.1)

# Subplot 2: CP vs a
axs[1].plot(A, CT_froude * (1 - A), label="Froude", color='black', linewidth=3, linestyle='-')
axs[1].plot(A, CT_Steiros * (1 - A), label="Steiros & Hullmark (2018)", color='black', linewidth=1.5, linestyle='-', marker='.', markersize=8)
if not mit_cp.empty:
    axs[1].plot(mit_cp['a'], mit_cp['CP'], label='Liew et al. (2024)', color='black', linewidth=1.5, linestyle='-', marker='x', markersize=8)
if not nrel_cp.empty:
    axs[1].scatter(nrel_cp['a'], nrel_cp['CP'], label='Martinez-Tossas et al. (2020)', color='black', marker='+', s=100, linewidths=1.5)

axs[1].set_xlabel("Induction factor $a$", fontsize=f1+2)
axs[1].set_ylabel("$C_P=C_T(1-a)$", fontsize=f1+2)
axs[1].grid(False)

# Compute CT and CP from new_AD model
EE = [(0.02, 0), (0.05, 0), (0.1, 0), (0.1, 0.6)]
CT = np.zeros_like(A)

for E1, E2 in EE:
    for i, a in enumerate(A):
        # Result scalar is located at index 4
        CT[i] = new_AD(E1, E2, I, a=a)[4]
    
    axs[0].plot(A, CT, label=f'$E_1 = ${E1:.2f}, $E_2$ = {E2:.2f}', linewidth=3, linestyle='--')
    CP = CT * (1 - A)
    axs[1].plot(A, CP, label=f'$E_1 = ${E1:.2f}, $E_2$ = {E2:.2f}', linewidth=3, linestyle='--')
    
axs[0].text(0.05, 0.9, "(a)", transform=axs[0].transAxes, fontsize=20, weight="bold")
axs[1].text(0.05, 0.9, "(b)", transform=axs[1].transAxes, fontsize=20, weight="bold")
 
for ax in axs:
    ax.tick_params(axis='both', labelsize=f1)
    ax.set_xlim([0, 0.75])
  
# Clean up legends to avoid duplication
handles, labels = axs[0].get_legend_handles_labels()

# Place the legend at the top, centered, with a visible border
fig.legend(handles, labels, 
           fontsize=f1, 
           loc='upper center',          
           bbox_to_anchor=(0.5, 1.08),  
           ncol=3, 
           frameon=True,                
           edgecolor='black',           
           framealpha=1.0)              

axs[0].set_ylim([-0.05, 1.5])
axs[1].set_ylim([-0.05, 0.75])

# Adjust tight_layout to close the gap
plt.tight_layout(rect=[0, 0, 1, 0.82]) 
plt.show()

a = 0.01 and E1 = 0.02 --> CT=0.04 - converged after 1 iterations. 
a = 0.06 and E1 = 0.02 --> CT=0.22 - converged after 1 iterations. 
a = 0.11 and E1 = 0.02 --> CT=0.38 - converged after 2 iterations. 
a = 0.16 and E1 = 0.02 --> CT=0.53 - converged after 12 iterations. 
a = 0.21 and E1 = 0.02 --> CT=0.66 - converged after 15 iterations. 
a = 0.26 and E1 = 0.02 --> CT=0.77 - converged after 16 iterations. 
a = 0.31 and E1 = 0.02 --> CT=0.86 - converged after 13 iterations. 
a = 0.36 and E1 = 0.02 --> CT=0.93 - converged after 3 iterations. 
a = 0.41 and E1 = 0.02 --> CT=0.99 - converged after 13 iterations. 
a = 0.46 and E1 = 0.02 --> CT=1.04 - converged after 16 iterations. 
a = 0.51 and E1 = 0.02 --> CT=1.06 - converged after 12 iterations. 
a = 0.56 and E1 = 0.02 --> CT=1.08 - converged after 10 iterations. 
a = 0.61 and E1 = 0.02 --> CT=1.08 - converged after 12 iterations. 
a = 0.66 and E1 = 0.02 --> NOT converged after 501 iterations.
a = 0.71 and E1 = 0.02 --> NOT converged after 501 iterations.
a = 0.76 and E1 = 0.02 --> NOT converged after 501 iterations.
a = 0.81 and E1 = 0.02 --> NOT converged after 501 iterations.
a = 0.01 and E1 = 0.05 --> CT=0.04 - converged after 1 iterations. 
a = 0.06 and E1 = 0.05 --> CT=0.22 - converged after 13 iterations. 
a = 0.11 and E1 = 0.05 --> CT=0.38 - converged after 13 iterations. 
a = 0.16 and E1 = 0.05 --> CT=0.52 - converged after 7 iterations. 
a = 0.21 and E1 = 0.05 --> CT=0.65 - converged after 1 iterations. 
a = 0.26 and E1 = 0.05 --> CT=0.76 - converged after 9 iterations. 
a = 0.31 and E1 = 0.05 --> CT=0.86 - converged after 12 iterations. 
a = 0.36 and E1 = 0.05 --> CT=0.94 - converged after 12 iterations. 
a = 0.41 and E1 = 0.05 --> CT=1.01 - converged after 8 iterations. 
a = 0.46 and E1 = 0.05 --> CT=1.06 - converged after 6 iterations. 
a = 0.51 and E1 = 0.05 --> CT=1.11 - converged after 13 iterations. 
a = 0.56 and E1 = 0.05 --> CT=1.13 - converged after 11 iterations. 
a = 0.61 and E1 = 0.05 --> CT=1.15 - converged after 9 iterations. 
a = 0.66 and E1 = 0.05 --> CT=1.16 - converged after 8 iterations. 
a = 0.71 and E1 = 0.05 --> NOT converged after 501 iterations.
a = 0.76 and E1 = 0.05 --> NOT converged after 501 iterations.
a = 0.81 and E1 = 0.05 --> NOT converged after 501 iterations.
a = 0.01 and E1 = 0.10 --> CT=0.04 - converged after 6 iterations. 
a = 0.06 and E1 = 0.10 --> CT=0.21 - converged after 22 iterations. 
a = 0.11 and E1 = 0.10 --> CT=0.36 - converged after 24 iterations. 
a = 0.16 and E1 = 0.10 --> CT=0.50 - converged after 23 iterations. 
a = 0.21 and E1 = 0.10 --> CT=0.63 - converged after 20 iterations. 
a = 0.26 and E1 = 0.10 --> CT=0.75 - converged after 14 iterations. 
a = 0.31 and E1 = 0.10 --> CT=0.85 - converged after 1 iterations. 
a = 0.36 and E1 = 0.10 --> CT=0.94 - converged after 12 iterations. 
a = 0.41 and E1 = 0.10 --> CT=1.02 - converged after 16 iterations. 
a = 0.46 and E1 = 0.10 --> CT=1.09 - converged after 17 iterations. 
a = 0.51 and E1 = 0.10 --> CT=1.14 - converged after 16 iterations. 
a = 0.56 and E1 = 0.10 --> CT=1.19 - converged after 14 iterations. 
a = 0.61 and E1 = 0.10 --> CT=1.22 - converged after 11 iterations. 
a = 0.66 and E1 = 0.10 --> CT=1.25 - converged after 7 iterations. 
a = 0.71 and E1 = 0.10 --> CT=1.26 - converged after 4 iterations. 
a = 0.76 and E1 = 0.10 --> CT=1.27 - converged after 4 iterations. 
a = 0.81 and E1 = 0.10 --> CT=1.28 - converged after 310 iterations. 
a = 0.01 and E1 = 0.10 --> CT=0.03 - converged after 25 iterations. 
a = 0.06 and E1 = 0.10 --> CT=0.19 - converged after 34 iterations. 
a = 0.11 and E1 = 0.10 --> CT=0.35 - converged after 30 iterations. 
a = 0.16 and E1 = 0.10 --> CT=0.50 - converged after 26 iterations. 
a = 0.21 and E1 = 0.10 --> CT=0.63 - converged after 21 iterations. 
a = 0.26 and E1 = 0.10 --> CT=0.75 - converged after 15 iterations. 
a = 0.31 and E1 = 0.10 --> CT=0.85 - converged after 1 iterations. 
a = 0.36 and E1 = 0.10 --> CT=0.94 - converged after 12 iterations. 
a = 0.41 and E1 = 0.10 --> CT=1.02 - converged after 16 iterations. 
a = 0.46 and E1 = 0.10 --> CT=1.09 - converged after 17 iterations. 
a = 0.51 and E1 = 0.10 --> CT=1.14 - converged after 17 iterations. 
a = 0.56 and E1 = 0.10 --> CT=1.19 - converged after 14 iterations. 
a = 0.61 and E1 = 0.10 --> CT=1.22 - converged after 11 iterations. 
a = 0.66 and E1 = 0.10 --> CT=1.25 - converged after 7 iterations. 
a = 0.71 and E1 = 0.10 --> CT=1.26 - converged after 4 iterations. 
a = 0.76 and E1 = 0.10 --> CT=1.27 - converged after 4 iterations. 
a = 0.81 and E1 = 0.10 --> NOT converged after 501 iterations.