Kernel: Python 3 (system-wide)
In [0]:
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 plt.style.use(['science', 'no-latex']) from cycler import cycler import matplotlib as mpl
Generalised Actuator Disk (AD) model
Note that variables are dimensionless: distances are normalised by disk diameter D, and velocities are normalised by the incoming velocity U0.
Equation labels in the code refer to the model formulation in the accompanying manuscript.
In [0]:
#-------------------------------- # 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): W_hat = E1 * (1.0 - U) # Eq. (2.5) A_hat = E2 * I # Eq. (2.6) n = 4 return (W_hat**n + A_hat**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) # 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) 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] # Eq. (2.24) Y = trapezoid( (1.0 - x_pos / np.sqrt(x_pos**2 + R_hat**2)) * np.gradient(sigma_pos**2, dx), x_pos) # 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] # Eq. (2.24) Y = trapezoid( (1.0 - x_pos / np.sqrt(x_pos**2 + R_hat**2)) * np.gradient(sigma_pos**2, dx), x_pos) # 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
In [0]:
########### DEFINE INPUT PARAMETERS ###################### I = 0.05 # incoming turbulence intensity E1 = 0.1 # wake-shear entrainment coefficient E2 = 0.3 # Background-turbulence entrainment coefficient A=[.2,.3,.4,.5,.6] # Induction factor values CT = np.zeros_like(A)
In [0]:
############### PLOTTING ####################### f_base = 14 # font size matplotlib.rcParams.update({ "font.size": f_base, "axes.labelsize": f_base, "axes.titlesize": f_base, "xtick.labelsize": f_base-2, "ytick.labelsize": f_base-2, "legend.fontsize": f_base-1, "lines.linewidth": 1.5, "grid.linestyle": ":", "grid.alpha": 0.6, "axes.grid": True, "pdf.fonttype": 42, "ps.fonttype": 42 }) linestyles = [ "-", # solid "--", # dashed "-.", # dash-dot ":", # dotted (0, (3, 1, 1, 1)) # custom dash: short dash, short gap, dot, short gap ] markers = ["o", "s", "D", "^", "v", "P", "X", "*"] colors = [ "#4C72B0", # muted blue "#55A868", # muted green "#C44E52", # muted red "#8172B2", # muted purple "#CCB974" # muted yellow-brown ] # generate figure fig, axs = plt.subplots(3, 1, figsize=(6, 5), sharex=True) line_handles, line_labels = [], [] for i, a in enumerate(A): x, sigma, U, Ue, ct = new_AD(E1, E2, I,a=a) # Pressure relation (Eq. 2.14) normalized by 0.5 rho U0^2 P = (1-(1-a)**2)*(1+x/np.sqrt(x**2+(1/2)**2))*np.heaviside(-x, .5) + \ (1-(1-a)**2-ct)*(1-x/np.sqrt(x**2+(1/2)**2))*np.heaviside(x, .5) # Plot line, = axs[1].plot(x, sigma, markevery=100, label=fr"$a={a:.2f},C_T={ct:.2f}$",linestyle=linestyles[i],color=colors[i]) axs[0].plot(x, U,linestyle=linestyles[i],color=colors[i]) axs[2].plot(x, P,linestyle=linestyles[i],color=colors[i]) line_handles.append(line) line_labels.append(fr"$a={a:.2f},C_T={ct:.2f}$") # Labels axs[1].set_ylabel(r"$\sigma/D$") axs[0].set_ylabel(r"$U/U_0$") axs[2].set_ylabel(r"$P/(0.5\,\rho\,U_0^2)$") axs[2].set_xlabel(r"$x/D$") # Panel letters for ax, label in zip(axs, ["(a)", "(b)", "(c)"]): ax.text(0.01, 0.85, label, transform=ax.transAxes, va="top", ha="left", fontsize=f_base, fontweight="bold") # Axes limits, ticks, and grids for ax in axs: ax.set_xlim([-2, 10]) ax.minorticks_on() ax.grid(True, which="major") ax.grid(True, which="minor", alpha=0.25) axs[0].set_ylim([0, 1.1]) axs[1].set_ylim([0, 2.0]) axs[2].set_ylim([-.9, .9]) # Legend leg = fig.legend( line_handles, line_labels, loc="center left", bbox_to_anchor=(1.0, 0.5), borderaxespad=0, frameon=False, ) plt.setp(leg.get_title(), fontsize=f_base) # Layout & save plt.tight_layout(rect=[0, 0, 0.98, 1]) # leave room for the legend plt.show()
In [0]:
In [0]:
In [0]:
In [0]: