diff --git a/pyrpod/plume/RarefiedPlumeGasKinetics.py b/pyrpod/plume/RarefiedPlumeGasKinetics.py index e571e27..47be556 100644 --- a/pyrpod/plume/RarefiedPlumeGasKinetics.py +++ b/pyrpod/plume/RarefiedPlumeGasKinetics.py @@ -34,15 +34,193 @@ ' = simplified analytical results """ +import warnings + import numpy as np -import sympy as sp from scipy import integrate -import matplotlib.pyplot as plt +from scipy.special import erf #define constants AVOGADROS_NUMBER = 6.0221e23 GAS_CONSTANT = 8.314 + +def get_K_factor(Q, S_0): + ''' + Scaled special factor exp(-S_0^2) * K [Cai & Wang 2012, Eq. 10]. + + The exp(-S_0^2) prefactor of the field solutions (Eqs. 5-8, 14) + is combined analytically with Eq. 10's exp(Q * S_0^2) into + exp(-S_0^2 * (1 - Q)). Since Q <= 1 both exponentials are <= 1, + so this never overflows for large speed ratios, while the plain + exp(Q * S_0^2) would. The [1 + erf(...)] factor is bounded by 2 + and needs no scaling. + + Parameters + ---------- + Q : float + special factor Q (full or simplified), 0 < Q <= 1 + S_0 : float + molecular speed ratio at the nozzle exit + + Returns + ------- + float + exp(-S_0^2) * K(Q, S_0) + ''' + erf_term = (1 + erf(S_0 * np.sqrt(Q))) * np.exp(-S_0 ** 2 * (1 - Q)) + term1 = Q * S_0 * np.exp(-S_0 ** 2) + term2 = (0.5 + Q * S_0 ** 2) * np.sqrt(np.pi * Q) + return Q * (term1 + term2 * erf_term) + + +def get_M_factor(Q, S_0): + ''' + Scaled special factor exp(-S_0^2) * M [Cai & Wang 2012, Eq. 11]. + + Overflow-safe exponential combination as in get_K_factor. + + Parameters + ---------- + Q : float + special factor Q (full or simplified), 0 < Q <= 1 + S_0 : float + molecular speed ratio at the nozzle exit + + Returns + ------- + float + exp(-S_0^2) * M(Q, S_0) + ''' + erf_term = (1 + erf(S_0 * np.sqrt(Q))) * np.exp(-S_0 ** 2 * (1 - Q)) + term1 = (1 + Q * S_0 ** 2) * np.exp(-S_0 ** 2) + term2 = S_0 * (1.5 + Q * S_0 ** 2) * np.sqrt(np.pi * Q) + return Q ** 2 * (term1 + term2 * erf_term) + + +def get_N_factor(Q, S_0): + ''' + Scaled special factor exp(-S_0^2) * N [Cai & Wang 2012, Eq. 12]. + + Overflow-safe exponential combination as in get_K_factor. + + Parameters + ---------- + Q : float + special factor Q (full or simplified), 0 < Q <= 1 + S_0 : float + molecular speed ratio at the nozzle exit + + Returns + ------- + float + exp(-S_0^2) * N(Q, S_0) + ''' + erf_term = (1 + erf(S_0 * np.sqrt(Q))) * np.exp(-S_0 ** 2 * (1 - Q)) + term1 = S_0 * Q ** 2 * (1.25 + Q * S_0 ** 2 / 2) * np.exp(-S_0 ** 2) + term2 = 0.5 * np.sqrt(np.pi * Q ** 3) + term3 = 0.75 + 3 * Q * S_0 ** 2 + Q ** 2 * S_0 ** 4 + return term1 + term2 * term3 * erf_term + + +def get_Q_full(r, epsilon, X, Z): + ''' + Full special factor Q [Cai & Wang 2012, Eq. 9] in closed form. + + Eq. 9 defines Q = cos^2(psi) * [sum_n P_n(sin(psi) sin(epsilon)) + * (r/sqrt(X^2+Z^2))^n]^2 with P_n the Legendre polynomials. The + series is the Legendre generating function sum_n P_n(x) t^n = + 1/sqrt(1 - 2*x*t + t^2) evaluated at x = sin(psi)sin(epsilon), + t = r/sqrt(X^2+Z^2), which collapses (with cos^2(psi) = + X^2/(X^2+Z^2)) to + + Q = X^2 / (X^2 + Z^2 - 2*Z*r*sin(epsilon) + r^2) + + so no series truncation or convergence loop is needed. The + denominator equals X^2 + (Z - r*sin(epsilon))^2 + + (r*cos(epsilon))^2 >= X^2 > 0, hence 0 < Q <= 1 always, which + also guarantees the overflow-safe exponential combination in + get_K_factor. On the centerline (Z = 0) this reduces to + Q(r) = X^2/(X^2 + r^2). + + Parameters + ---------- + r : float or ndarray + radial integration variable over the exit disk, 0 <= r <= R_0 (m) + epsilon : float or ndarray + angular integration variable, -pi/2 <= epsilon <= pi/2 (rad) + X : float + axial coordinate of the field point, X > 0 (m) + Z : float + transverse coordinate of the field point (m) + + Returns + ------- + float or ndarray + special factor Q at (r, epsilon) for field point (X, 0, Z) + ''' + return X ** 2 / (X ** 2 + Z ** 2 - 2 * Z * r * np.sin(epsilon) + r ** 2) + + +def get_far_field_velocity_normalized(S_0): + ''' + Far-field centerline velocity asymptote, lim U_1 * sqrt(beta_0) + as X -> infinity [Cai & Wang 2012, Eqs. 22 and 24]. A function of + the exit speed ratio only: far enough from the exit the problem + degenerates to a small-hole effusion flow, so no nozzle geometry + factors remain. + + Evaluated by dividing numerator and denominator through by + [1 + erf(S_0)] * exp(S_0^2), which would itself overflow for + S_0 >~ 26; the residual factor exp(-S_0^2) / (1 + erf(S_0)) is + bounded for S_0 > 0. + + Parameters + ---------- + S_0 : float + molecular speed ratio at the nozzle exit, S_0 > 0 + + Returns + ------- + float + lim (X -> infinity) U_1(X, 0, 0) * sqrt(beta_0) + ''' + inv_E = np.exp(-S_0 ** 2) / (1 + erf(S_0)) + sqpi = np.sqrt(np.pi) + return S_0 + (inv_E + sqpi * S_0) / (S_0 * inv_E + (0.5 + S_0 ** 2) * sqpi) + + +def get_far_field_temp_ratio(S_0): + ''' + Far-field centerline temperature asymptote, lim T_1/T_0 as + X -> infinity [Cai & Wang 2012, Eqs. 23-24]. A function of the + exit speed ratio only (see get_far_field_velocity_normalized). + + Implemented as -2/3 * G^2 + 4*N(Q=1)/(3*K(Q=1)) with G of + Eq. 24, i.e. the Q -> 1 limit of Eq. 17; this reads Eq. 23's + printed denominator "3 S_0 + (1/2 + S_0^2) sqrt(pi) [...]" with + the 3 distributing over the whole denominator, 3*K(Q=1). The + equivalence was verified independently to 40 digits. Overflow + safety as in get_far_field_velocity_normalized. + + Parameters + ---------- + S_0 : float + molecular speed ratio at the nozzle exit, S_0 > 0 + + Returns + ------- + float + lim (X -> infinity) T_1(X, 0, 0) / T_0 + ''' + inv_E = np.exp(-S_0 ** 2) / (1 + erf(S_0)) + sqpi = np.sqrt(np.pi) + G = get_far_field_velocity_normalized(S_0) + num = S_0 * (5 + 2 * S_0 ** 2) * inv_E + 2 * sqpi * (0.75 + 3 * S_0 ** 2 + S_0 ** 4) + den = 3 * (S_0 * inv_E + (0.5 + S_0 ** 2) * sqpi) + return -(2 / 3) * G ** 2 + num / den + + def get_maxwellian_pressure(rho_inf, U, S, sigma, theta, T, T_w): ''' Rarefied Gas Dynamics - Shen - eq. 4.19 @@ -79,7 +257,7 @@ def get_maxwellian_pressure(rho_inf, U, S, sigma, theta, T, T_w): p1 *= np.exp(- (S * np.cos(theta)) ** 2) p2 = (2 - sigma) * ((S * np.cos(theta)) ** 2 + 0.5) p2 += (S * np.cos(theta) * (sigma / 2) * np.sqrt(np.pi * T_w / T)) - p2 *= 1 + sp.erf(S * np.cos(theta)) + p2 *= 1 + erf(S * np.cos(theta)) p = p1 + p2 p *= (rho_inf * U ** 2) / (2 * S ** 2) return p @@ -111,7 +289,7 @@ def get_maxwellian_shear_pressure(rho_inf, U, S, sigma, theta): tau1 = np.exp(- (S * np.cos(theta)) ** 2) tau2 = np.sqrt(np.pi) * S * np.cos(theta) - tau2 *= (1 + sp.erf(S * np.cos(theta))) + tau2 *= (1 + erf(S * np.cos(theta))) tau = tau1 + tau2 tau *= -(sigma * rho_inf * np.sin(theta) * U ** 2) / (2 * np.sqrt(np.pi) * S) return tau @@ -150,7 +328,7 @@ def get_maxwellian_heat_transfer(rho_inf, S, sigma, theta, T, T_r, R, gamma): the pressure exerted on the surface element (N / m^2) ''' q = (S ** 2) + (gamma / (gamma - 1)) - (((gamma + 1) * T_r) / (2 * (gamma - 1) * T)) - q *= np.exp(- (S * np.cos(theta)) ** 2) + (np.sqrt(np.pi) * (S * np.cos(theta)) * (1 + sp.erf(S * np.cos(theta)))) + q *= np.exp(- (S * np.cos(theta)) ** 2) + (np.sqrt(np.pi) * (S * np.cos(theta)) * (1 + erf(S * np.cos(theta)))) q -= 0.5 * np.exp(- (S * np.cos(theta)) ** 2) q *= sigma * rho_inf * R * T * np.sqrt(R * T / (2 * np.pi)) return q @@ -206,7 +384,7 @@ class Simons: Number density from continuity equation with constant mass flux across different spherical surfaces. Return the ratio of number density at an analyzed point outside of the exit vs at the exit. ''' - def __init__(self, gamma, R, T_c, P_c, R_0, r): + def __init__(self, gamma, R, T_c, P_c, R_0, r, kappa=None): ''' Simple constructor, saves parameters to self. @@ -223,8 +401,15 @@ def __init__(self, gamma, R, T_c, P_c, R_0, r): R_0 : float nozzle exit radius (m) r : float - distance from the evaluated point + distance from the evaluated point to the nozzle exit center (m) + kappa : float, optional + plume beaming exponent for the cosine-law decay + function [Cai & Wang 2012, Eq. 26]. Defaults to + Boyton's kappa = 2/(gamma - 1) [Cai & Wang 2012, + Sec. II.C, refs. 22-23]. Paper-cited alternatives: + kappa = 2 (Ashkenas & Sherman, ref. 20) and + kappa = 1/(gamma - 1) (Albini, ref. 21). Returns ------- @@ -237,12 +422,11 @@ def __init__(self, gamma, R, T_c, P_c, R_0, r): self.P_c = P_c self.r = r self.R_0 = R_0 + if kappa is None: + kappa = 2 / (gamma - 1) #Boyton 1967/68 from Cai2012 [22][23] + self.kappa = kappa self.set_normalization_constant() - #self.rho_throat = self.get_nozzle_throat_density() - #self.theta_max = self.get_limiting_turn_angle() - #self.A = self.get_normalization_constant() - #self.U_t = self.get_limiting_velocity() def get_nozzle_throat_density(self): ''' @@ -259,8 +443,10 @@ def get_nozzle_throat_density(self): ''' #ideal gas law P_throat = rho_throat * R * T_throat #therefore: rho_throat = P_throat / (R * T_throat) - P_throat = self.P_c * 0.5283 #from Isentropic Flow Tables @ M = 1 - T_throat = self.T_c * 0.8333 #from Isentropic Flow Tables @ M = 1 + #isentropic ratios at M = 1: T*/T_c = 2/(gamma+1), + #P*/P_c = (2/(gamma+1))^(gamma/(gamma-1)) + T_throat = self.T_c * (2 / (self.gamma + 1)) + P_throat = self.P_c * (2 / (self.gamma + 1)) ** (self.gamma / (self.gamma - 1)) rho_throat = P_throat / (self.R * T_throat) return rho_throat @@ -281,28 +467,39 @@ def get_limiting_turn_angle(self): def get_plume_angular_density_decay_function(self, theta): ''' - From Cai 2012. Solve for the density decay function at a given off-centerline angle. - Expression for kappa is chosen from Boyton 1967/68. + Solve for the density decay function at a given off-centerline + angle [Cai & Wang 2012, Eq. 26], f(theta) = + cos^kappa(pi*theta/(2*theta_max)) with the beaming exponent + kappa chosen at construction (default Boyton 2/(gamma-1)). + + For theta >= theta_max the plume model region is empty (the + gas cannot turn past the limiting angle), so the decay + function is 0. Without this guard the cosine goes negative + and fractional kappa would produce complex numbers. Parameters ---------- theta : float plume centerline off-angle of current position (rad) - + Returns ------- float evaluation of plume angular density decay function at theta ''' - kappa = 2 / (self.gamma - 1) #Boyton 1967/68 from Cai2012 [22][23] theta_max = self.get_limiting_turn_angle() - f = (np.cos((np.pi / 2) * (theta / theta_max))) ** kappa + if theta >= theta_max: + return 0.0 + f = (np.cos((np.pi / 2) * (theta / theta_max))) ** self.kappa return f def set_normalization_constant(self): ''' From Lumpkin 1999. Setter for normalization constant. - Integrates theta from 0 to the limiting turn angle. + Numerically integrates sin(theta) * cos^kappa(pi*theta/(2*theta_max)) + from 0 to the limiting turn angle. The integrand is real and + non-negative on [0, theta_max] for any kappa, so no complex + arithmetic can arise. Parameters ---------- @@ -313,17 +510,13 @@ def set_normalization_constant(self): None. ''' theta_max = self.get_limiting_turn_angle() - theta = sp.symbols('theta') - f = (sp.cos((sp.pi / 2) * (theta / theta_max))) ** (2 / (self.gamma - 1)) - integrand = sp.sin(theta) * f - integral = sp.integrate(integrand, (theta, 0, theta_max)) #integrate from 0 to max turning angle - A = 0.5 * np.sqrt((self.gamma - 1) / (self.gamma + 1)) / (integral.evalf()) - - # if A is complex, take the real part - if type(A) == sp.Add: - A = sp.re(A) + kappa = self.kappa + + def integrand(theta): + return np.sin(theta) * np.cos((np.pi / 2) * (theta / theta_max)) ** kappa - self.A = A + integral, _ = integrate.quad(integrand, 0, theta_max) + self.A = 0.5 * np.sqrt((self.gamma - 1) / (self.gamma + 1)) / integral return def get_sonic_velocity(self): @@ -340,7 +533,7 @@ def get_sonic_velocity(self): float sonic velocity (m/s) ''' - T_throat = self.T_c * 0.8333 #from Isentropic Flow Tables @ Mach = 1 + T_throat = self.T_c * (2 / (self.gamma + 1)) #isentropic ratio at M = 1 sonic_velocity = np.sqrt(self.gamma * self.R * T_throat) return sonic_velocity @@ -374,25 +567,68 @@ def get_static_pressure(self, rho_ratio): def get_num_density_ratio(self, theta): ''' - Number density from continuity equation with constant mass flux across different spherical surfaces. - Return the ratio of number density at an analyzed point outside of the exit vs at the exit. + Number density from continuity equation with constant mass flux + across different spherical surfaces [Cai & Wang 2012, Eq. 25]. + + NOTE: this ratio is THROAT-referenced, n/n_s: the returned + value normalizes by the number density at the nozzle throat + (rho_s in the paper), not at the nozzle exit. For a ratio + comparable to the gas-kinetic classes (which normalize by the + exit density n_0), use get_num_density_ratio_exit. + + Returns 0.0 for theta >= theta_max (empty plume region). Parameters ---------- theta : float Angle off centerline of the current point being analyzed (rad). - + Returns ------- float - the ratio of number density at an analyzed point outside of the exit vs at the exit + number density at the analyzed point normalized by the + nozzle THROAT number density, n/n_s ''' + theta_max = self.get_limiting_turn_angle() + if theta >= theta_max: + return 0.0 f = self.get_plume_angular_density_decay_function(theta) - #??? make own function for rho_ratio??? rho_ratio = self.A * ((self.R_0/self.r) ** 2) * f #rho_ratio = density / nozzle throat denisty aka rho / rho_s n_ratio = rho_ratio return n_ratio + def get_num_density_ratio_exit(self, theta, exit_mach): + ''' + Exit-referenced number density ratio n/n_0. + + Rescales the throat-referenced Eq. 25 result by the isentropic + density ratio between the throat (M = 1) and the nozzle exit + (M = exit_mach): + + n_s / n_0 = [(1 + (gamma-1)/2 * M_e^2) / ((gamma+1)/2)]^(1/(gamma-1)) + + This makes the Simons model directly comparable to the + gas-kinetic classes (SimplifiedGasKinetics, + CollisionlessGasKinetics), which normalize by the exit + number density n_0. + + Parameters + ---------- + theta : float + Angle off centerline of the current point being analyzed (rad). + exit_mach : float + Mach number at the nozzle exit plane. + + Returns + ------- + float + number density at the analyzed point normalized by the + nozzle EXIT number density, n/n_0 + ''' + throat_to_exit = ((1 + (self.gamma - 1) / 2 * exit_mach ** 2) + / ((self.gamma + 1) / 2)) ** (1 / (self.gamma - 1)) + return self.get_num_density_ratio(theta) * throat_to_exit + # TODO save plume constants into self class SimplifiedGasKinetics: ''' @@ -600,8 +836,11 @@ def set_Q_simple(self): def set_K_simple(self): ''' - Setter for simplified special factor K. This simplification is just the - substitution of Q for Q'. + Setter for simplified special factor K [Cai & Wang 2012, Eq. 10] + with Q substituted by Q'. Stored scaled by exp(-S_0^2) for + overflow safety (see get_K_factor); ratio methods (U, W, T) + are unaffected since the scaling cancels, and the density + method uses the scaled factor directly. Parameters ---------- @@ -611,21 +850,15 @@ def set_K_simple(self): ------- None. ''' - #K_simple = Q_simple * ((Q_simple * S_0) + ((0.5 + (Q_simple * S_0 ** 2)) * np.sqrt(np.pi * Q_simple) * - #(1 + sp.erf(S_0 * np.sqrt(Q_simple))) ** (Q_simple * S_0 ** 2))) - term1 = self.Q_simple * self.S_0 - term2 = 0.5 + self.Q_simple * self.S_0 ** 2 - term3 = np.sqrt(np.pi * self.Q_simple) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(self.Q_simple))) * np.exp(self.Q_simple * self.S_0 ** 2) - K_simple = self.Q_simple * (term1 + term2 * term3 * term4) - self.K_simple = K_simple + self.K_simple = get_K_factor(self.Q_simple, self.S_0) return - + def set_M_simple(self): ''' - Setter for simplified special factor M. This simplification is just the - substitution of Q for Q'. + Setter for simplified special factor M [Cai & Wang 2012, Eq. 11] + with Q substituted by Q'. Stored scaled by exp(-S_0^2) for + overflow safety (see get_K_factor). Paramters --------- @@ -635,21 +868,15 @@ def set_M_simple(self): ------- None. ''' - #M_simple = (Q_simple ** 2) * ((Q_simple * S_0 ** 2) + 1 + (S_0 * (1.5 + (Q_simple * S_0 ** 2)) * - #np.sqrt(np.pi * Q_simple)) * (1 + sp.erf(S_0 * np.sqrt(Q_simple))) ** (Q_simple *S_0 ** 2)) - term1 = 1 + self.Q_simple * self.S_0 ** 2 - term2 = self.S_0 * (1.5 + self.Q_simple * self.S_0 ** 2) - term3 = np.sqrt(np.pi * self.Q_simple) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(self.Q_simple))) * np.exp(self.Q_simple * self.S_0 ** 2) - M_simple = self.Q_simple ** 2 * (term1 + term2 * term3 * term4) - self.M_simple = M_simple + self.M_simple = get_M_factor(self.Q_simple, self.S_0) return - + def set_N_simple(self): ''' - Setter for simplified special factor N. This simplification is just the - substitution of Q for Q'. + Setter for simplified special factor N [Cai & Wang 2012, Eq. 12] + with Q substituted by Q'. Stored scaled by exp(-S_0^2) for + overflow safety (see get_K_factor). Parameters ---------- @@ -659,14 +886,7 @@ def set_N_simple(self): ------- None. ''' - #N_simple = S_0 * (Q_simple ** 2) * (1.25 + (Q_simple * S_0 ** 2) / 2) - #N_simple += (0.5 * np.sqrt(np.pi * Q_simple ** 3)) * (0.75 + 3 * Q_simple * S_0 **2 + Q_simple ** 2 * S_0 ** 4) * (1 + sp.erf(S_0 * np.sqrt(Q_simple))) ** (Q_simple * S_0 ** 2) - term1 = self.S_0 * self.Q_simple ** 2 * (1.25 + self.Q_simple * self.S_0 ** 2 / 2) - term2 = 0.5 * np.sqrt(np.pi * self.Q_simple ** 3) - term3 = 0.75 + 3 * self.Q_simple * self.S_0 ** 2 + self.Q_simple ** 2 * self.S_0 ** 4 - term4 = (1 + sp.erf(self.S_0 * np.sqrt(self.Q_simple))) * np.exp(self.Q_simple * self.S_0 ** 2) - N_simple = term1 + term2 * term3 * term4 - self.N_simple = N_simple + self.N_simple = get_N_factor(self.Q_simple, self.S_0) return @@ -684,10 +904,9 @@ def get_num_density_ratio(self): float number density at a point (X, 0, Z) vs number density at the nozzle exit ''' - # num_density_ratio = n_1s(X, 0, Z) / n_0 + # num_density_ratio = n_1s(X, 0, Z) / n_0 [Cai & Wang 2012, Eq. 14] + # K_simple already carries the exp(-S_0^2) prefactor (see set_K_simple) num_density_ratio = (self.K_simple / (2 * np.sqrt(np.pi)) * (self.R_0 / self.X) ** 2) - num_density_ratio *= np.exp(-(self.S_0 ** 2)) - num_density_ratio = float(num_density_ratio) return num_density_ratio def get_U_normalized(self): @@ -705,9 +924,8 @@ def get_U_normalized(self): float returns U normalized ''' - # U_normalized = U_1s (X, 0, Z) * sqrt(beta) + # U_normalized = U_1s (X, 0, Z) * sqrt(beta) [Cai & Wang 2012, Eq. 15] U_normalized = self.M_simple / self.K_simple - U_normalized = float(U_normalized) return U_normalized def get_W_normalized(self): @@ -725,9 +943,8 @@ def get_W_normalized(self): float returns W normalized ''' - # W_normalized = W_1s (X, 0, Z) * sqrt(beta) + # W_normalized = W_1s (X, 0, Z) * sqrt(beta) [Cai & Wang 2012, Eq. 16] W_normalized = (self.M_simple / self.K_simple) * (self.Z / self.X) - W_normalized = float(W_normalized) return W_normalized def get_temp_ratio(self): @@ -744,11 +961,9 @@ def get_temp_ratio(self): float ratio of temperature at a point (X, 0, Z) to the temperature at the nozzle exit ''' - # T_ratio = T_1s / T_0 + # T_ratio = T_1s / T_0 [Cai & Wang 2012, Eq. 17] T_ratio = ((-2 * self.M_simple ** 2) / (3 * self.Q_simple * self.K_simple ** 2)) T_ratio += (4 * self.N_simple / (3 * self.K_simple)) - T_ratio = float(T_ratio) - #print(f'S0 = {S_0}, Qs = {Q_simple}, Ks = {K_simple}, Ms = {M_simple}, Ns = {N_simple}, T_ratio = {T_ratio}') return T_ratio def get_num_density_centerline(self): @@ -766,10 +981,10 @@ def get_num_density_centerline(self): number density at a point on the centerline vs the number density at the nozzle exit ''' - p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) + # n_1(X, 0, 0) / n_0 [Cai & Wang 2012, Eq. 18] + p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) p2 = self.R_0 / np.sqrt(self.X ** 2 + self.R_0 ** 2) - n_ratio = 0.5 + 0.5 * sp.erf(self.S_0) - (p1 * np.exp(-self.S_0 ** 2 * p2 ** 2) / 2) * (1 + sp.erf(p1 * self.S_0)) - n_ratio = float(n_ratio) + n_ratio = 0.5 + 0.5 * erf(self.S_0) - (p1 * np.exp(-self.S_0 ** 2 * p2 ** 2) / 2) * (1 + erf(p1 * self.S_0)) return n_ratio @@ -791,11 +1006,11 @@ def get_velocity_centerline(self): velocity at a point on the centerline (X, 0, 0) normalized with the parameter beta at the exit ''' - p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) + # U_1(X, 0, 0) * sqrt(beta_0) [Cai & Wang 2012, Eq. 19] + p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) p2 = self.R_0 / np.sqrt(self.X ** 2 + self.R_0 ** 2) n_ratio = self.get_num_density_centerline() - U_ratio = 1 / (2 * n_ratio) * ((p2 ** 2 * np.exp(- self.S_0 ** 2) / np.sqrt(np.pi)) + (self.S_0 * (1 + sp.erf(self.S_0))) - (np.exp(- p2 ** 2 * self.S_0 ** 2) * p1 ** 3 * self.S_0 * (1 + sp.erf(p1 * self.S_0)))) - U_ratio = float(U_ratio) + U_ratio = 1 / (2 * n_ratio) * ((p2 ** 2 * np.exp(- self.S_0 ** 2) / np.sqrt(np.pi)) + (self.S_0 * (1 + erf(self.S_0))) - (np.exp(- p2 ** 2 * self.S_0 ** 2) * p1 ** 3 * self.S_0 * (1 + erf(p1 * self.S_0)))) return U_ratio def get_temp_centerline(self): @@ -803,6 +1018,13 @@ def get_temp_centerline(self): Method to calculate the temperature at a point (X, 0, 0) outside of the nozzle. This temperature is normalized over the temperature at the nozzle exit. + Evaluates Eq. 21's integral of N(Q(r)) * r over the exit + radius by numerical quadrature with the exact centerline + Q(r) = X^2/(X^2 + r^2) [Cai & Wang 2012, Eqs. 9, 21]. (An + earlier version treated N as constant across the exit, + integral = 0.5 * N(Q') * R_0^2, which is only valid in the + far field X >> R_0 where Q(r) ~ Q' = 1.) + Parameters ---------- None. @@ -810,20 +1032,20 @@ def get_temp_centerline(self): Returns ------- float - temperature on the point on the centerline (X, 0, 0) + temperature on the point on the centerline (X, 0, 0) vs temperature at the nozzle exit ''' + # T_1(X, 0, 0) / T_0 [Cai & Wang 2012, Eq. 21] + # get_N_factor carries the exp(-S_0^2) prefactor of Eq. 21 n_ratio = self.get_num_density_centerline() U1 = self.get_velocity_centerline() - ''' - r = sp.symbols("r") - f = N * r - integral = sp.integrate(f, (r, 0, R_0)) - temp_ratio = (4 * np.exp(- S_0 ** 2)) / (3 * n_ratio * np.sqrt(np.pi) * X ** 2) * integral.evalf() - (U1 ** 2 / (3/2)) - ''' - integral = 0.5 * self.N_simple * self.R_0 ** 2 - temp_ratio = (4 * np.exp(- self.S_0 ** 2)) / (3 * n_ratio * np.sqrt(np.pi) * self.X ** 2) * integral - (U1 ** 2 / (3/2)) - temp_ratio = float(temp_ratio) + + def integrand(r): + Q = get_Q_full(r, 0.0, self.X, 0.0) + return get_N_factor(Q, self.S_0) * r + + integral, _ = integrate.quad(integrand, 0, self.R_0) + temp_ratio = 4 / (3 * n_ratio * np.sqrt(np.pi) * self.X ** 2) * integral - (U1 ** 2 / (3/2)) return temp_ratio def get_pressure(self): @@ -946,191 +1168,307 @@ def get_heat_flux(self): heat_flux = get_maxwellian_heat_transfer(rho_inf, S, self.sigma, self.theta, T, self.T_w, self.R, self.gamma) return heat_flux -''' -import math -from scipy.special import legendre -class CollisionlessPlume: - def __init__(self, U_0, R, T_0, n_iters, conv_tol): - self.n_iters = n_iters - self.conv_tol = conv_tol - self.U_0 = U_0 - self.R = R - self.T_0 = T_0 - self.S_0 = self.get_speed_ratio() +class CollisionlessGasKinetics(SimplifiedGasKinetics): + ''' + Full collisionless analytical plume model [Cai & Wang 2012, + Sec. II.A, Eqs. 5-12]: a free jet expanding from a round exit + into vacuum, evaluated at a point (X, 0, Z) in front of the + nozzle (X > 0). Unlike the parent SimplifiedGasKinetics (which + substitutes the far-field Q' of Eq. 13), this class integrates + the exact special factor Q of Eq. 9 -- via its closed form, see + get_Q_full -- over the finite exit disk, so it remains valid in + the near field. + + The field solutions (Eqs. 5-8) are double integrals over + r in [0, R_0] and epsilon in [-pi/2, pi/2]. They are evaluated + with tensor-product Gauss-Legendre quadrature: the integrands + are analytic on the compact rectangle (the denominator of Q is + bounded below by X^2 > 0), so Gauss-Legendre converges + geometrically. The order is doubled (40 -> 80 -> 160) until two + successive orders agree to QUAD_RTOL; the finer result is kept. + In practice order 40 already reaches machine precision except + very near the nozzle lip (X -> 0, Z ~ R_0), where the density + field has a singularity [Cai & Wang 2012, Sec. III]. + + Centerline closed forms (Eqs. 18-21) and the Maxwell gas-surface + interface are inherited from SimplifiedGasKinetics; the exact + analytical solutions reduce to Eqs. 18-21 on the centerline, so + the inherited methods are exact there. The inherited surface + methods (get_pressure, get_shear_pressure, get_heat_flux) + dispatch to this class's overridden field getters, so they are + fed by the full-model n, U, W, T. + + Attributes + ---------- + (all of SimplifiedGasKinetics, plus) + + I_K : float + integral of r * exp(-S_0^2) * K over the exit disk [Eq. 5] + + I_M : float + integral of r * exp(-S_0^2) * M over the exit disk [Eq. 6] + + I_W : float + integral of (Z - r sin(epsilon)) * r * exp(-S_0^2) * M over + the exit disk [Eq. 7] + + I_N : float + integral of r * exp(-S_0^2) * N over the exit disk [Eq. 8] + + Methods + ------- + get_num_density_ratio() + + get_U_normalized() + + get_W_normalized() + + get_Vr_normalized() + + get_temp_ratio() + + get_pressure_ratio() + + (get_pressure / get_shear_pressure / get_heat_flux and the + centerline closed forms are inherited from SimplifiedGasKinetics) + ''' + + QUAD_ORDERS = (40, 80, 160) + QUAD_RTOL = 1e-9 + + def __init__(self, distance, theta, thruster_characteristics, T_w, sigma): + ''' + Mirrors SimplifiedGasKinetics(distance, theta, + thruster_characteristics, T_w, sigma) exactly; X = d*cos(theta) + and Z = d*sin(theta) are derived internally. Additionally + precomputes the four field integrals of Eqs. 5-8. + ''' + super().__init__(distance, theta, thruster_characteristics, T_w, sigma) + self.set_field_integrals() + return - def get_speed_ratio(self): - S_0 = self.U_0 / np.sqrt(2 * self.R * self.T_0) - return S_0 + def _compute_field_integrals(self, order): + ''' + Evaluate the four exit-disk integrals of Eqs. 5-8 with a + tensor-product Gauss-Legendre rule of the given order per + axis, over r in [0, R_0] and epsilon in [-pi/2, pi/2]. - #solves a summation series of legengre polynomials of the first kind - #from degree 0 to degree n or until the convergance tolerance, tol is met - #n: int, x: float (value to evaluate function over), tol: float - def get_Q(self, X, Z, r, epsilon): - - # TEST TODO TEST TODO TEST + The K, M, N factors carry the exp(-S_0^2) prefactor of the + field solutions (see get_K_factor), keeping every term + bounded for large speed ratios. - psi = np.arctan(Z/X) - leg_x = np.sin(psi) * np.sin(epsilon) + Parameters + ---------- + order : int + number of Gauss-Legendre nodes per axis - sum = 0 - - P_n_minus_2 = 1 #P_0 = 1 - P_n_minus_1 = leg_x #P_1 = x - sum += P_n_minus_2 + (P_n_minus_1 * r / np.sqrt(X ** 2 + Z ** 2)) + Returns + ------- + tuple of float + (I_K, I_M, I_W, I_N) + ''' + nodes, weights = np.polynomial.legendre.leggauss(order) + # map [-1, 1] to [0, R_0] (radial) and [-pi/2, pi/2] (angular) + r = 0.5 * self.R_0 * (nodes + 1) + w_r = 0.5 * self.R_0 * weights + eps = 0.5 * np.pi * nodes + w_eps = 0.5 * np.pi * weights + + R, E = np.meshgrid(r, eps, indexing='ij') + W2D = np.outer(w_r, w_eps) + + Q = get_Q_full(R, E, self.X, self.Z) + K = get_K_factor(Q, self.S_0) + M = get_M_factor(Q, self.S_0) + N = get_N_factor(Q, self.S_0) + + I_K = np.sum(W2D * R * K) + I_M = np.sum(W2D * R * M) + I_W = np.sum(W2D * (self.Z - R * np.sin(E)) * R * M) + I_N = np.sum(W2D * R * N) + return I_K, I_M, I_W, I_N + + def set_field_integrals(self): + ''' + Setter for the field integrals I_K, I_M, I_W, I_N of + Eqs. 5-8, with quadrature-order doubling until two + successive orders agree to QUAD_RTOL (see class docstring + for the accuracy rationale). Warns if the finest order is + reached without convergence (only possible extremely close + to the nozzle-lip singularity). - - for degree in range (2, self.n_iters + 1, 1): - - P_n = legendre(degree)(leg_x) - #print(f'degree: {degree}; x: {leg_x}; P_n: {P_n}') - sum += P_n - if abs(P_n) < self.conv_tol and degree % 2 == 0: + Parameters + ---------- + None. + + Returns + ------- + None. + ''' + previous = None + for order in self.QUAD_ORDERS: + current = self._compute_field_integrals(order) + if previous is not None and self._integrals_converged(previous, current): break - - #my recurrence legendre solver - #solve for polynomial of degree of current iter - - P_n = (2 * (degree - 1) + 1) * leg_x * (P_n_minus_1) - P_n -= (degree - 1) * P_n_minus_2 - P_n /= degree + previous = current + else: + warnings.warn( + 'CollisionlessGasKinetics quadrature did not converge to ' + 'rtol={} at order {} for point (X={}, Z={}); using finest ' + 'result.'.format(self.QUAD_RTOL, self.QUAD_ORDERS[-1], + self.X, self.Z), + RuntimeWarning) - sum += P_n * (r / np.sqrt(X ** 2 + Z ** 2)) ** degree + self.I_K, self.I_M, self.I_W, self.I_N = current - P_n_minus_2 = P_n_minus_1 - P_n_minus_1 = P_n + return - if abs(P_n) < self.conv_tol: - break - - sum = 0.712 - Q = (np.cos(psi) ** 2) * (sum ** 2) + def _integrals_converged(self, previous, current): + ''' + Convergence test between two quadrature orders. I_K, I_M and + I_N are strictly positive, so a relative test applies; I_W + can legitimately vanish (centerline, Eq. 20), so its + difference is measured against the natural magnitude of its + integrand, (|Z| + R_0) * I_M. - return Q - - def get_K(self, X, Z, r, epsilon): - #K = Q * ((Q * S_0) + ((0.5 + (Q * S_0 ** 2)) * sqrt(pi * Q) * - #(1 + erf(S_0 * sqrt(Q))) ** (Q * S_0 ** 2))) - Q = self.get_Q(X, Z, r, epsilon) - term1 = Q * self.S_0 - term2 = 0.5 + Q * self.S_0 ** 2 - term3 = np.sqrt(np.pi * Q) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(Q))) * np.exp(Q * self.S_0 ** 2) - K = Q * (term1 + term2 * term3 * term4) - return K - - def get_M(self, X, Z, r, epsilon): - #M = (Q ** 2) * ((Q * S_0 ** 2) + 1 + (S_0 * (1.5 + (Q * S_0 ** 2)) * - #sqrt(pi * Q_simple)) * (1 + erf(S_0 * sqrt(Q))) ** (Q *S_0 ** 2)) - Q = self.get_Q(X, Z, r, epsilon) - term1 = Q * self.S_0 ** 2 - term2 = 1 + self.S_0 * (1.5 + Q * self.S_0 ** 2) - term3 = np.sqrt(np.pi * Q) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(Q))) * np.exp(Q * self.S_0 ** 2) - M = Q ** 2 * (term1 + term2 * term3 * term4) - return M - - def get_N(self, X, Z, r, epsilon): - #N = S_0 * (Q ** 2) * (1.25 + (Q * S_0 ** 2) / 2) + - #(0.5 * sqrt(pi * Q ** 3)) * (0.75 + 3 * Q * S_0 **2 + Q ** 2 * S_0 ** 4) * - #(1 + erf(S_0 * sqrt(Q))) ** (Q * S_0 ** 2) - Q = self.get_Q(X, Z, r, epsilon) - term1 = self.S_0 * Q ** 2 * (1.25 + Q * self.S_0 ** 2 / 2) - term2 = 0.5 * np.sqrt(np.pi * Q ** 3) - term3 = 0.75 + 3 * Q * self.S_0 ** 2 + Q ** 2 * self.S_0 ** 4 - term4 = (1 + sp.erf(self.S_0 * np.sqrt(Q))) * np.exp(Q * self.S_0 ** 2) - N = term1 + term2 * term3 * term4 - return N - - def get_num_density_ratio(self, X, Z, R_0): - n_ratio = np.exp(-self.S_0 ** 2) / (X ** 2 * np.pi ** (3/2)) - n = 15 - e_a, e_b = -np.pi, np.pi - e_h = (e_b - e_a) / n - - e_vals = np.arange(e_a, e_b, e_h) - e_vals = np.append(e_vals, e_b) - dbl_integral = 0 - - epsilon = e_a - for idx, epsilon in enumerate(e_vals): - if epsilon == e_a or idx == len(e_vals) - 1: - integral = self.get_integral(X, Z, R_0, epsilon) - dbl_integral += integral - elif ((epsilon / e_h) % 3 == 0): - integral = self.get_integral(X, Z, R_0, epsilon) - dbl_integral += 2 * integral - else: - integral = self.get_integral(X, Z, R_0, epsilon) - dbl_integral += 3 * integral - dbl_integral *= 3 * e_h / 8 - n_ratio *= dbl_integral - print(f'S0: {self.S_0}; X: {X}; Z: {Z}; {n_ratio}') - return n_ratio + Parameters + ---------- + previous : tuple of float + integrals from the coarser rule + current : tuple of float + integrals from the finer rule - def get_integral(self, X, Z, R_0, epsilon): - - r_a, r_b = 0, R_0 - n = 15 - r_h = (r_b - r_a) / n - - r_vals = np.arange(r_a, r_b, r_h) - r_vals = np.append(r_vals, r_b) - integral = 0 - r = r_a - for idx, r in enumerate(r_vals): - if r == r_a or idx == len(r_vals) - 1: - integral += r * self.get_K(X, Z, r, epsilon) - elif ((r / r_h) % 3 == 0): - integral += 2 * r * self.get_K(X, Z, r, epsilon) - else: - integral += 3 * r * self.get_K(X, Z, r, epsilon) - integral *= (3 * r_h / 8) - - return integral + Returns + ------- + bool + True when every integral has converged to QUAD_RTOL + ''' + rtol = self.QUAD_RTOL + I_K0, I_M0, I_W0, I_N0 = previous + I_K1, I_M1, I_W1, I_N1 = current + scale_W = (abs(self.Z) + self.R_0) * abs(I_M1) + return (abs(I_K1 - I_K0) <= rtol * abs(I_K1) + and abs(I_M1 - I_M0) <= rtol * abs(I_M1) + and abs(I_N1 - I_N0) <= rtol * abs(I_N1) + and abs(I_W1 - I_W0) <= rtol * scale_W) -import numpy as np -import matplotlib.pyplot as plt - -# Assuming you have an instance of CollisionlessPlume named plume - -# Define the range of X / (2 * R_0) values -x_values = np.linspace(0.3, 10, 100) # Adjust the range as needed - -# Values of S_0 to plot -u0_values = [2] -R_0 = 5 -# Create subplots -plt.figure(figsize=(10, 6)) -plt.title('Normalized analytical number density along centerline') -plt.xlabel('X / D') -plt.ylabel('num_density_ratio') - -# Plot for each S_0 value -for u0 in u0_values: - num_density_ratios = [] - plume = CollisionlessPlume(u0 * 15.15255, 0.287, 400, 999, 0.00001) - # Calculate num_density_ratio for each X / (2 * R_0) - for x_ratio in x_values: - X = x_ratio * 2 * R_0 # Calculate X from the ratio - Z = 0 # Fixed Z value - num_density_ratio = plume.get_num_density_ratio(X, Z, R_0) - num_density_ratios.append(num_density_ratio) - - # Plot the results - plt.plot(x_values, num_density_ratios, label=f'S_0 = {u0}') - -# Show legend -plt.legend() - -# Show the plot -plt.show() - -#test = CollisionlessPlume.get_num_density_ratio() - #def get_U_normalized(): - #def get_W_normalized(): - #def get_temp_ratio(): - -''' \ No newline at end of file + def get_num_density_ratio(self): + ''' + Number density at (X, 0, Z) normalized by the exit number + density, n_1/n_0 [Cai & Wang 2012, Eq. 5]. The exp(-S_0^2) + prefactor is carried inside I_K (see get_K_factor). + + Parameters + ---------- + None. + + Returns + ------- + float + n_1(X, 0, Z) / n_0 + ''' + return self.I_K / (np.pi ** 1.5 * self.X ** 2) + + def get_U_normalized(self): + ''' + Macroscopic x-velocity at (X, 0, Z) normalized by + sqrt(beta_0), i.e. U_1 * sqrt(beta_0) [Cai & Wang 2012, + Eq. 6]. The shared prefactor and n_0/n_1 reduce the ratio to + I_M / I_K. + + Parameters + ---------- + None. + + Returns + ------- + float + U_1(X, 0, Z) * sqrt(beta_0) + ''' + return self.I_M / self.I_K + + def get_W_normalized(self): + ''' + Macroscopic z-velocity at (X, 0, Z) normalized by + sqrt(beta_0), i.e. W_1 * sqrt(beta_0) [Cai & Wang 2012, + Eq. 7]. Eq. 7's integrand factor is read as + (Z - r sin(epsilon)); the printed "(Z - r sin(theta))" is a + typo -- sin(epsilon) is the convention consistent with Q + (Eq. 9) and it makes W vanish on the centerline as Eq. 20 + requires. Eq. 7's extra 1/X (prefactor 1/X^3 vs 1/X^2) + reduces the ratio to I_W / (X * I_K). + + Parameters + ---------- + None. + + Returns + ------- + float + W_1(X, 0, Z) * sqrt(beta_0) + ''' + return self.I_W / (self.X * self.I_K) + + def get_Vr_normalized(self): + ''' + Radial (spherical, from the nozzle exit center) velocity + component in the Y = 0 plane, normalized by sqrt(beta_0): + V_r = (X*U + Z*W)/sqrt(X^2 + Z^2) [Cai & Wang 2012, + Figs. 16-18]. + + Parameters + ---------- + None. + + Returns + ------- + float + V_r(X, 0, Z) * sqrt(beta_0) + ''' + U = self.get_U_normalized() + W = self.get_W_normalized() + return (self.X * U + self.Z * W) / np.sqrt(self.X ** 2 + self.Z ** 2) + + def get_temp_ratio(self): + ''' + Temperature at (X, 0, Z) normalized by the exit temperature, + T_1/T_0 [Cai & Wang 2012, Eq. 8]. With beta_0 = 1/(2*R*T_0), + the kinetic term -(U_1^2 + W_1^2)/(3*R*T_0) equals + -(2/3) * [(U_1 sqrt(beta_0))^2 + (W_1 sqrt(beta_0))^2]. + + Parameters + ---------- + None. + + Returns + ------- + float + T_1(X, 0, Z) / T_0 + ''' + U = self.get_U_normalized() + W = self.get_W_normalized() + return -(2 / 3) * (U ** 2 + W ** 2) + (4 / 3) * self.I_N / self.I_K + + def get_pressure_ratio(self): + ''' + Flowfield static pressure at (X, 0, Z) normalized by the exit + static pressure: p_1/p_0 = (n_1/n_0) * (T_1/T_0) from the + ideal gas law with the LOCAL temperature. Cai & Wang 2012 + (p. 64, Fig. 13 discussion) explicitly warn that computing + the local pressure as n(X, 0, Z) * k * T_0 is theoretically + invalid, because the flowfield temperature is always lower + than the exit temperature T_0. + + Parameters + ---------- + None. + + Returns + ------- + float + p_1(X, 0, Z) / p_0 + ''' + return self.get_num_density_ratio() * self.get_temp_ratio() diff --git a/tests/plume/data/digitized/README.md b/tests/plume/data/digitized/README.md new file mode 100644 index 0000000..9f43dd1 --- /dev/null +++ b/tests/plume/data/digitized/README.md @@ -0,0 +1,29 @@ +# Digitized reference data for the plume verification figures + +Drop digitized curves from Cai & Wang 2012 (JSR 49(1), DOI +10.2514/1.A32046) here and the manual-run figure scripts +(`tests/plume/plume_verification_test_04` ... `_27`) will overlay them +automatically on their next run. No code changes are needed. + +## File convention + +- Name: `_.csv`, where `` is the paper figure + (e.g. `fig19`) and `` names the dataset. Examples: + - `fig03_dsmc.csv` — DSMC centerline density (S0 = 2), Fig. 3 + - `fig19_dsmc.csv` — DSMC centerline density, Kn = 100, Fig. 19 + - `fig25_dsmc_kn100.csv`, `fig25_dsmc_kn0p1.csv`, + `fig25_dsmc_kn0p01.csv` — the three DSMC mass-flux curves, Fig. 25 + - contour figures: one CSV per digitized contour polyline, e.g. + `fig06_dsmc_0p001.csv` (the DSMC n/n0 = 0.001 line of Fig. 6) +- Content: two comma-separated columns `x,y` with one header row. + x/y are in the figure's plotted units (X/D, Z/D, theta in degrees, + normalized quantity values). +- Series suffixes are prettified for legends: `dsmc` -> `DSMC`, + `kn0p1` -> `Kn=0.1`. +- Lower-half-plane series in the split contour figures (Figs. 6-8, + 11-18) should be digitized with positive Z/D values taken from the + paper's lower half; enter them as printed (negative y) — overlays are + drawn exactly as given. + +Each figure script's docstring lists the exact filenames it looks for. +Files that are absent are silently skipped. diff --git a/tests/plume/output/.gitignore b/tests/plume/output/.gitignore new file mode 100644 index 0000000..d6b7ef3 --- /dev/null +++ b/tests/plume/output/.gitignore @@ -0,0 +1,2 @@ +* +!.gitignore diff --git a/tests/plume/plume_figure_utils.py b/tests/plume/plume_figure_utils.py new file mode 100644 index 0000000..194772e --- /dev/null +++ b/tests/plume/plume_figure_utils.py @@ -0,0 +1,559 @@ +# ======================== +# PyRPOD: tests/plume/plume_figure_utils.py +# ======================== +# Shared helpers for the manual-run verification figure scripts +# plume_verification_test_04 ... _27, which reproduce Figs. 2-25 of +# Cai, C. and Wang, L., "Numerical Validations for a Set of Collisionless +# Rocket Plume Solutions," JSR 49(1), 2012, DOI 10.2514/1.A32046. +# +# Design notes +# ------------ +# * All figures use the paper's validation conditions: argon, +# D = 0.2 m, gamma = 5/3, R = 208.13 J/(kg K). Every plotted quantity +# is normalized (n/n0, U*sqrt(beta0), T/T0, ...), so the absolute +# T_0 and n_0 values are inert; the ones below match the existing +# verification tests. ve is set from the exit speed ratio, +# ve = S_0 * sqrt(2*R*T_0). +# * Single-source-of-truth policy: the simplified and Simons overlays +# delegate point-by-point to SimplifiedGasKinetics / Simons (they are +# closed-form and cheap), so any model change is reflected in the +# figures automatically. Only the FULL model keeps a vectorized +# re-implementation of its Eq. 5-8 quadrature here (per-point class +# instantiation would cost hours per contour sweep instead of +# seconds); as a guard, every call cross-checks sample points against +# CollisionlessGasKinetics and raises if they disagree beyond 1e-6, +# so that path cannot silently diverge from the physics module +# either -- a model change makes the figure scripts fail loudly +# until this file is updated to match. +# * The models require X > 0; grids start at X_MIN_OVER_D and angular +# sweeps stop at THETA_MAX_DEG < 90 deg. +# * Digitized-data convention (see tests/plume/data/digitized/README.md): +# any CSV named '_.csv' (two columns x,y, one header +# row) is overlaid automatically by overlay_digitized(); absent files +# are silently skipped, so the DSMC series appear once the paper's +# curves are digitized. +# * Simons overlays are EXIT-referenced (n/n_0) via the isentropic +# throat-to-exit conversion with Me = S_0*sqrt(2/gamma) (from +# S_0 = U0/sqrt(2RT0) and M = U0/sqrt(gamma*R*T0)). Following the +# paper's remark that the cosine-law curves with different kappa +# coincide on the centerline (p. 65: "the cosine-law model with +# different kappa values actually provides the same value"), the +# normalization constant A is always taken from the Boyton default +# kappa while the plotted decay exponent may vary (Figs. 22-24). + +import sys +import time +from pathlib import Path + +import numpy as np +import matplotlib +matplotlib.use('Agg') # save-only figures; must precede pyplot import +import matplotlib.pyplot as plt # noqa: E402 (backend must be set first) + +_THIS_DIR = Path(__file__).resolve().parent +_REPO_ROOT = _THIS_DIR.parents[1] +if str(_REPO_ROOT) not in sys.path: # direct `python tests/plume/...py` runs + sys.path.insert(0, str(_REPO_ROOT)) + +from pyrpod.plume.RarefiedPlumeGasKinetics import ( # noqa: E402 + CollisionlessGasKinetics, + SimplifiedGasKinetics, + Simons, + get_K_factor, + get_M_factor, + get_N_factor, + get_Q_full, +) + +# Paper validation conditions (Sec. III: argon, D = 0.2 m) +R_SPECIFIC = 208.13 # J / (kg K) +T_0 = 500.0 # K (inert for normalized results) +GAMMA = 5.0 / 3.0 +N_0 = 1.0e20 # m^-3 (inert for normalized results) +D_NOZZLE = 0.2 # m +R_0 = D_NOZZLE / 2 # m +T_W = 300.0 # K (unused by plotted quantities) +SIGMA = 1.0 +BOLTZMANN = 1.380649e-23 # J / K + +OUTPUT_DIR = _THIS_DIR / 'output' +DIGITIZED_DIR = _THIS_DIR / 'data' / 'digitized' + +X_MIN_OVER_D = 0.05 # models require X > 0 +THETA_MAX_DEG = 89.5 # angular sweeps stop short of 90 deg (X -> 0) +QUAD_ORDER = 80 # per-axis Gauss-Legendre order for grid fields +GRID_N = 161 # nodes along the longer contour-grid axis + +_GL_NODES, _GL_WEIGHTS = np.polynomial.legendre.leggauss(QUAD_ORDER) + + +# --------------------------------------------------------------------------- +# model factories +# --------------------------------------------------------------------------- + +def thruster_characteristics(S_0): + """Thruster dict for the gas-kinetic classes at exit speed ratio S_0.""" + ve = S_0 * np.sqrt(2 * R_SPECIFIC * T_0) + return {'d': D_NOZZLE, 've': ve, 'R': R_SPECIFIC, + 'gamma': GAMMA, 'Te': T_0, 'n': N_0} + + +def make_plume(cls, distance, theta, S_0): + """Instantiate SimplifiedGasKinetics or CollisionlessGasKinetics.""" + return cls(distance, theta, thruster_characteristics(S_0), T_W, SIGMA) + + +def exit_mach(S_0): + """Me from S_0: S_0 = U0/sqrt(2RT0), M = U0/sqrt(gamma R T0).""" + return S_0 * np.sqrt(2.0 / GAMMA) + + +def make_simons(distance, S_0, kappa=None): + """Simons instance with chamber conditions isentropically consistent + with the exit state at Me(S_0). T_c/P_c only rescale absolute + density; they are inert for the plotted ratios.""" + Me = exit_mach(S_0) + stag = 1 + (GAMMA - 1) / 2 * Me ** 2 + T_c = T_0 * stag + P_c = (N_0 * BOLTZMANN * T_0) * stag ** (GAMMA / (GAMMA - 1)) + return Simons(GAMMA, R_SPECIFIC, T_c, P_c, R_0, distance, kappa=kappa) + + +def simons_density_field(X, Z, S_0, kappa_f=None): + """Exit-referenced Simons density n/n_0 at points (X, 0, Z) [m], + delegated point-by-point to Simons.get_num_density_ratio_exit so + that any change to the Simons model is reflected here. + + Following the paper (p. 65: the cosine-law curves coincide at + theta = 0 for every kappa), the normalization constant A is always + the Boyton-kappa value; kappa_f optionally varies only the plotted + decay exponent f = cos^kappa_f (Figs. 22-24), implemented by + overriding A on a kappa_f-built instance with the Boyton value. + A single instance is reused with its evaluation radius updated per + point (reconstructing would re-run the A quadrature 10^4 times). + """ + X = np.asarray(X, dtype=float) + Z = np.asarray(Z, dtype=float) + shape = np.broadcast(X, Z).shape + Xf = np.broadcast_to(X, shape).ravel() + Zf = np.broadcast_to(Z, shape).ravel() + + simons = make_simons(1.0, S_0, kappa=kappa_f) + if kappa_f is not None: + simons.A = make_simons(1.0, S_0).A # paper's shared Boyton A + Me = exit_mach(S_0) + + out = np.empty(Xf.size) + for i in range(Xf.size): + simons.r = float(np.hypot(Xf[i], Zf[i])) + theta = float(np.arctan2(np.abs(Zf[i]), Xf[i])) + out[i] = simons.get_num_density_ratio_exit(theta, Me) + return out.reshape(shape) + + +# --------------------------------------------------------------------------- +# centerline profiles (exact closed forms, Eqs. 18/19/21) +# --------------------------------------------------------------------------- + +def centerline_profiles(x_over_D, S_0, quantities=('n',)): + """Exact analytical centerline values at X/D array; quantities from + {'n', 'U', 'T'} -> Eqs. 18 / 19 / 21.""" + getters = {'n': 'get_num_density_centerline', + 'U': 'get_velocity_centerline', + 'T': 'get_temp_centerline'} + out = {q: np.empty(len(x_over_D)) for q in quantities} + for i, xd in enumerate(x_over_D): + plume = make_plume(SimplifiedGasKinetics, xd * D_NOZZLE, 0.0, S_0) + for q in quantities: + out[q][i] = getattr(plume, getters[q])() + return out + + +# --------------------------------------------------------------------------- +# simplified-model field evaluation (delegated to SimplifiedGasKinetics) +# --------------------------------------------------------------------------- + +def simplified_field_values(quantity, X, Z, S_0): + """Simplified-model quantity at 1-D point arrays X, Z [m], + delegated point-by-point to SimplifiedGasKinetics (Eqs. 13-17) so + that any change to the simplified model is reflected here. The + class is closed-form, so the per-point loop costs well under a + second even on full contour grids. + + quantity: 'n', 'U', 'W', 'T', plus the derived compositions + 'Vr' = (X*U + Z*W)/sqrt(X^2+Z^2) and 'p' = n*T (the class exposes + no methods for those). + """ + getters = {'n': 'get_num_density_ratio', 'U': 'get_U_normalized', + 'W': 'get_W_normalized', 'T': 'get_temp_ratio'} + if quantity not in getters and quantity not in ('Vr', 'p'): + raise KeyError(quantity) + X = np.asarray(X, dtype=float).ravel() + Z = np.asarray(Z, dtype=float).ravel() + out = np.empty(X.size) + for i in range(X.size): + d = float(np.hypot(X[i], Z[i])) + theta = float(np.arctan2(Z[i], X[i])) + plume = make_plume(SimplifiedGasKinetics, d, theta, S_0) + if quantity in getters: + out[i] = getattr(plume, getters[quantity])() + elif quantity == 'Vr': + out[i] = (X[i] * plume.get_U_normalized() + + Z[i] * plume.get_W_normalized()) / d + else: # 'p' + out[i] = plume.get_num_density_ratio() * plume.get_temp_ratio() + return out + + +# --------------------------------------------------------------------------- +# vectorized full-model field evaluation (Eqs. 5-8) +# --------------------------------------------------------------------------- + +def full_field_values(quantity, X, Z, S_0, _chunk=256): + """Full-model field quantity at 1-D point arrays X, Z [m]. + + quantity: 'n', 'U', 'W', 'Vr', 'T', or 'p'. Same integrals as + CollisionlessGasKinetics at its converged order (80); each call + verifies sample points against the class to 1e-6 (see module + docstring). + """ + X = np.asarray(X, dtype=float).ravel() + Z = np.asarray(Z, dtype=float).ravel() + + r = 0.5 * R_0 * (_GL_NODES + 1) + w_r = 0.5 * R_0 * _GL_WEIGHTS + eps = 0.5 * np.pi * _GL_NODES + w_eps = 0.5 * np.pi * _GL_WEIGHTS + Rn, En = np.meshgrid(r, eps, indexing='ij') + W2D = np.outer(w_r, w_eps) + sinE = np.sin(En) + + n_pts = X.size + I_K = np.empty(n_pts) + I_M = np.empty(n_pts) + I_W = np.empty(n_pts) + I_N = np.empty(n_pts) + for start in range(0, n_pts, _chunk): + sl = slice(start, min(start + _chunk, n_pts)) + Xc = X[sl][:, None, None] + Zc = Z[sl][:, None, None] + Q = get_Q_full(Rn, En, Xc, Zc) + K = get_K_factor(Q, S_0) + M = get_M_factor(Q, S_0) + N = get_N_factor(Q, S_0) + I_K[sl] = np.sum(W2D * Rn * K, axis=(1, 2)) + I_M[sl] = np.sum(W2D * Rn * M, axis=(1, 2)) + I_W[sl] = np.sum(W2D * (Zc - Rn * sinE) * Rn * M, axis=(1, 2)) + I_N[sl] = np.sum(W2D * Rn * N, axis=(1, 2)) + + n = I_K / (np.pi ** 1.5 * X ** 2) + U = I_M / I_K + W = I_W / (X * I_K) + T = -(2 / 3) * (U ** 2 + W ** 2) + (4 / 3) * I_N / I_K + values = {'n': n, 'U': U, 'W': W, 'T': T, + 'Vr': (X * U + Z * W) / np.sqrt(X ** 2 + Z ** 2), + 'p': n * T} + _verify_against_class(values, X, Z, S_0) + return values[quantity] + + +def _verify_against_class(values, X, Z, S_0, rtol=1e-6): + """Cross-check sample points against CollisionlessGasKinetics. + + W legitimately vanishes on the centerline, so its deviation is + measured against the natural velocity scale max(|W|, |U|) instead + of a purely relative test. + """ + for idx in {0, X.size // 2, X.size - 1}: + d = float(np.hypot(X[idx], Z[idx])) + th = float(np.arctan2(Z[idx], X[idx])) + ref = make_plume(CollisionlessGasKinetics, d, th, S_0) + U_ref = ref.get_U_normalized() + W_ref = ref.get_W_normalized() + checks = {'n': (ref.get_num_density_ratio(), None), + 'U': (U_ref, None), + 'W': (W_ref, max(abs(W_ref), abs(U_ref))), + 'T': (ref.get_temp_ratio(), None)} + for q, (expected, scale) in checks.items(): + scale = abs(expected) if scale is None else scale + got = values[q][idx] + if abs(got - expected) > rtol * scale: + raise AssertionError( + f'vectorized field {q} diverged from ' + f'CollisionlessGasKinetics at (X={X[idx]}, Z={Z[idx]}): ' + f'{got} vs {expected}') + + +def full_field_grid(quantity, x_over_D, z_over_D, S_0): + """Full-model quantity on the tensor grid; returns shape + (len(z_over_D), len(x_over_D)) for use with plt.contour.""" + Xg, Zg = np.meshgrid(np.asarray(x_over_D) * D_NOZZLE, + np.asarray(z_over_D) * D_NOZZLE) + vals = full_field_values(quantity, Xg.ravel(), Zg.ravel(), S_0) + return vals.reshape(Xg.shape) + + +def simplified_field_grid(quantity, x_over_D, z_over_D, S_0): + """Simplified-model quantity on the tensor grid (delegates to + simplified_field_values); returns shape + (len(z_over_D), len(x_over_D)) for use with plt.contour.""" + Xg, Zg = np.meshgrid(np.asarray(x_over_D) * D_NOZZLE, + np.asarray(z_over_D) * D_NOZZLE) + vals = simplified_field_values(quantity, Xg.ravel(), Zg.ravel(), S_0) + return vals.reshape(Xg.shape) + + +def default_contour_axes(x_max=10.0, z_max=10.0): + """Paper-style contour grid: X/D in [X_MIN, x_max], Z/D in [0, z_max] + (fields are symmetric in Z; mirror for the lower half-plane).""" + x = np.linspace(X_MIN_OVER_D, x_max, GRID_N) + z = np.linspace(0.0, z_max, (GRID_N + 1) // 2) + return x, z + + +# --------------------------------------------------------------------------- +# digitized-data overlay (see data/digitized/README.md) +# --------------------------------------------------------------------------- + +def load_digitized(fig_stem): + """dict {series_suffix: (x, y)} from data/digitized/_*.csv.""" + datasets = {} + for path in sorted(DIGITIZED_DIR.glob(f'{fig_stem}_*.csv')): + data = np.genfromtxt(path, delimiter=',', skip_header=1) + data = np.atleast_2d(data) + datasets[path.stem[len(fig_stem) + 1:]] = (data[:, 0], data[:, 1]) + return datasets + + +def _pretty_series_label(suffix): + tokens = [] + for token in suffix.split('_'): + if token.lower() == 'dsmc': + tokens.append('DSMC') + elif token.lower().startswith('kn'): + tokens.append('Kn=' + token[2:].replace('p', '.')) + else: + tokens.append(token) + return ' '.join(tokens) + + +def overlay_digitized(ax, fig_stem, style='markers', mirror_z=False, **kw): + """Overlay every digitized series for fig_stem; returns #series. + + style 'markers' suits profile figures; 'line' suits digitized + contour polylines. mirror_z flips y (for lower-half-plane DSMC + slots in the split contour figures). + """ + count = 0 + for suffix, (x, y) in load_digitized(fig_stem).items(): + y = -y if mirror_z else y + label = _pretty_series_label(suffix) + if style == 'line': + ax.plot(x, y, color='k', lw=1.0, label=label, **kw) + else: + ax.plot(x, y, linestyle='none', marker='d', mfc='none', + color='k', ms=4, label=label, **kw) + count += 1 + return count + + +# --------------------------------------------------------------------------- +# figure assembly helpers +# --------------------------------------------------------------------------- + +def max_rel_diff(candidate, reference): + """max |candidate - reference| / |reference| over finite entries.""" + candidate = np.asarray(candidate, dtype=float) + reference = np.asarray(reference, dtype=float) + mask = np.isfinite(candidate) & np.isfinite(reference) & (reference != 0) + return float(np.max(np.abs(candidate[mask] - reference[mask]) + / np.abs(reference[mask]))) + + +def annotate_error(ax, text, loc='lower left'): + """Small model-difference annotation box on the axes.""" + positions = {'lower left': (0.02, 0.02, 'left', 'bottom'), + 'lower right': (0.98, 0.02, 'right', 'bottom'), + 'upper left': (0.02, 0.98, 'left', 'top')} + x, y, ha, va = positions[loc] + ax.text(x, y, text, transform=ax.transAxes, fontsize=7, ha=ha, va=va, + bbox=dict(boxstyle='round', fc='white', ec='0.6', alpha=0.85)) + + +def save_figure(fig, name): + """Save PNG to tests/plume/output and return the path.""" + OUTPUT_DIR.mkdir(parents=True, exist_ok=True) + path = OUTPUT_DIR / f'{name}.png' + fig.savefig(path, dpi=200, bbox_inches='tight') + plt.close(fig) + return path + + +def run_script(generate_figure): + """Manual-run entry point shared by all figure scripts.""" + t0 = time.perf_counter() + path = generate_figure() + print(f'saved {path} ({time.perf_counter() - t0:.1f} s)') + + +# --------------------------------------------------------------------------- +# shared builders for the Kn-triplet figures (identical analytic content; +# separate figures so per-Kn DSMC overlays land in the right file) +# --------------------------------------------------------------------------- + +def density_contour_figure(kn_label, fig_stem): + """Figs. 6-8: number-density contours, S_0 = 2. Upper half-plane: + analytical (solid) + simplified (dashed); lower half-plane: Simons + (dashed) + DSMC digitized slot (solid lines when data present).""" + S_0 = 2.0 + levels = [0.001, 0.01, 0.1] + x, z = default_contour_axes() + n_full = full_field_grid('n', x, z, S_0) + n_simp = simplified_field_grid('n', x, z, S_0) + Xg, Zg = np.meshgrid(x * D_NOZZLE, z * D_NOZZLE) + n_simons = simons_density_field(Xg, Zg, S_0) + + fig, ax = plt.subplots(figsize=(6, 7)) + cs = ax.contour(x, z, n_full, levels=levels, colors='k', + linewidths=1.2) + ax.clabel(cs, fmt='%g', fontsize=7) + cs = ax.contour(x, z, n_simp, levels=levels, colors='k', + linewidths=0.9, linestyles='dashed') + ax.clabel(cs, fmt='%g', fontsize=7) + cs = ax.contour(x, -z, n_simons, levels=levels, colors='k', + linewidths=0.9, linestyles='dashed') + ax.clabel(cs, fmt='%g', fontsize=7) + n_dig = overlay_digitized(ax, fig_stem, style='line', mirror_z=False) + + ax.axhline(0.0, color='k', lw=0.6) + ax.set_xlim(0, 10) + ax.set_ylim(-10, 10) + ax.set_xlabel('X/D') + ax.set_ylabel('Z/D') + ax.set_title(f'Normalized number density, Kn = {kn_label}, ' + f'$S_0$ = 2.0') + handles = [plt.Line2D([], [], color='k', lw=1.2, label='Analytical (top)'), + plt.Line2D([], [], color='k', lw=0.9, ls='--', + label='Simplified (top)'), + plt.Line2D([], [], color='k', lw=0.9, ls='--', + label='Simons (bottom)')] + if n_dig: + handles.append(plt.Line2D([], [], color='k', lw=1.0, + label='DSMC (bottom, digitized)')) + ax.legend(handles=handles, fontsize=7, loc='upper right') + far = x >= 1.0 # near the exit the simplified (R_0/X)^2 form diverges + annotate_error( + ax, + f'max |n_S/n_A - 1| = {max_rel_diff(n_simp[:, far], n_full[:, far]):.2g}\n' + f'max |n_Sim/n_A - 1| = {max_rel_diff(n_simons[:, far], n_full[:, far]):.2g}\n' + '(over grid with X/D $\\geq$ 1; largest in the plume tails)') + return save_figure(fig, fig_stem) + + +def top_bottom_contour_figure(quantity, levels, kn_label, fig_stem, + title_quantity, x_max=10.0, z_max=10.0, + z_axis_label='Z/D', fmt='%g'): + """Figs. 11-18 pattern: analytical contours in the upper half-plane, + lower half-plane reserved for the DSMC digitized overlay.""" + S_0 = 2.0 + x, z = default_contour_axes(x_max, z_max) + field = full_field_grid(quantity, x, z, S_0) + + fig, ax = plt.subplots(figsize=(6, 7)) + cs = ax.contour(x, z, field, levels=levels, colors='k', linewidths=1.1) + ax.clabel(cs, fmt=fmt, fontsize=7) + n_dig = overlay_digitized(ax, fig_stem, style='line') + + ax.axhline(0.0, color='k', lw=0.6) + ax.set_xlim(0, x_max) + ax.set_ylim(-z_max, z_max) + ax.set_xlabel('X/D') + ax.set_ylabel(z_axis_label) + ax.set_title(f'{title_quantity}, Kn = {kn_label}, $S_0$ = 2.0') + ax.text(0.7 * x_max, 0.55 * z_max, 'Analytical', fontsize=9) + ax.text(0.7 * x_max, -0.55 * z_max, + 'DSMC' if n_dig else 'DSMC (pending digitized data)', + fontsize=9) + return save_figure(fig, fig_stem) + + +def centerline_density_profile_figure(kn_label, fig_stem): + """Figs. 19-21: centerline density -- analytical (Eq. 5 integral, + circles), simplified (the Eq. 18 closed form from the paper's + simplified-solutions section II.B, triangles), Simons (solid line, + exit-referenced), DSMC digitized slot. As in the paper, the + analytical and simplified series coincide on the centerline (the + Eq. 5 integral reduces exactly to Eq. 18) and converge to 1 at + X = 0.""" + S_0 = 2.0 + x_over_D = np.linspace(X_MIN_OVER_D, 10.0, 200) + r = x_over_D * D_NOZZLE + analytical = full_field_values('n', r, np.zeros_like(r), S_0) + simplified = centerline_profiles(x_over_D, S_0, ('n',))['n'] + simons = simons_density_field(r, np.zeros_like(r), S_0) + + fig, ax = plt.subplots(figsize=(6, 4.5)) + ax.plot(x_over_D, simons, 'k-', lw=1.2, label='Simons') + ax.plot(x_over_D, analytical, linestyle='none', marker='o', mfc='none', + color='k', ms=4, markevery=4, label='Analytical') + ax.plot(x_over_D, simplified, linestyle='none', marker='>', color='k', + ms=4, markevery=(2, 4), label='Simplified') + overlay_digitized(ax, fig_stem) + + ax.set_xlim(0, 10) + ax.set_ylim(0, 1.2) + ax.set_xlabel('X/D') + ax.set_ylabel('$n/n_0$') + ax.set_title(f'Centerline density profiles, Kn = {kn_label}, ' + f'$S_0$ = 2.0') + ax.legend(fontsize=8) + tail = x_over_D >= 1.0 + annotate_error( + ax, + f'Eq. 5 vs Eq. 18: max rel diff = ' + f'{max_rel_diff(simplified, analytical):.1g}\n' + f'max |n_Sim/n_A - 1| = {max_rel_diff(simons[tail], analytical[tail]):.2g} ' + f'(X/D $\\geq$ 1)', loc='lower right') + return save_figure(fig, fig_stem) + + +def angular_density_profile_figure(kn_label, fig_stem): + """Figs. 22-24: density along r/D = 10 -- analytical, simplified, + Simons kappa = 1.5/2/3 (single Boyton-A normalization, see module + docstring), DSMC digitized slot.""" + S_0 = 2.0 + r = 10.0 * D_NOZZLE + theta = np.deg2rad(np.linspace(0.0, THETA_MAX_DEG, 90)) + X = r * np.cos(theta) + Z = r * np.sin(theta) + analytical = full_field_values('n', X, Z, S_0) + simplified = simplified_field_values('n', X, Z, S_0) + theta_deg = np.rad2deg(theta) + + fig, ax = plt.subplots(figsize=(6, 4.5)) + for kappa, ls in [(2.0, '-'), (1.5, '--'), (3.0, '-.')]: + simons = simons_density_field(X, Z, S_0, kappa_f=kappa) + ax.plot(theta_deg, simons, 'k', ls=ls, lw=1.1, + label=f'Simons $\\kappa$={kappa:g}') + ax.plot(theta_deg, analytical, linestyle='none', marker='o', mfc='none', + color='k', ms=4, markevery=2, label='Analytical') + ax.plot(theta_deg, simplified, linestyle='none', marker='>', color='k', + ms=4, markevery=(1, 2), label='Simplified') + overlay_digitized(ax, fig_stem) + + ax.set_xlim(0, 90) + ax.set_ylim(0, 0.015) + ax.set_xticks([0, 30, 60, 90]) + ax.set_xlabel(r'$\theta$') + ax.set_ylabel('$n/n_0$') + ax.set_title(f'Density profiles along r/D = 10, Kn = {kn_label}, ' + f'$S_0$ = 2.0') + ax.legend(fontsize=8) + simons_boyton = simons_density_field(X, Z, S_0) + annotate_error( + ax, + f'max |n_S/n_A - 1| = {max_rel_diff(simplified, analytical):.2g}\n' + f'max |n_Sim($\\kappa$=3)/n_A - 1| = ' + f'{max_rel_diff(simons_boyton, analytical):.2g}', + loc='lower left') + return save_figure(fig, fig_stem) diff --git a/tests/plume/plume_unit_test_02.py b/tests/plume/plume_unit_test_02.py new file mode 100644 index 0000000..8022e3a --- /dev/null +++ b/tests/plume/plume_unit_test_02.py @@ -0,0 +1,74 @@ +# ======================== +# PyRPOD: tests/plume/plume_unit_test_02.py +# ======================== +# Verifies the closed form of the special factor Q (Cai & Wang 2012, +# Eq. 9) against a direct truncation of the printed Legendre series. +# +# Eq. 9: Q = cos^2(psi) * [sum_n P_n(sin(psi)sin(eps)) * t^n]^2 with +# t = r/sqrt(X^2+Z^2). The series is the Legendre generating function, +# so Q collapses to X^2/(X^2 + Z^2 - 2*Z*r*sin(eps) + r^2). The series +# converges like t^n, so with the test points restricted to t <= 0.6 a +# 50-term truncation is accurate to ~t^51 ~ 5e-12, and rel=1e-9 both +# passes robustly and would catch any error in the closed form. + +import numpy as np +import pytest +from scipy.special import eval_legendre + +from pyrpod.plume.RarefiedPlumeGasKinetics import get_Q_full + +pytestmark = [pytest.mark.plume, pytest.mark.unit] + +N_TERMS = 50 + + +def q_series(r, epsilon, X, Z, n_terms=N_TERMS): + """Truncated Eq. 9 exactly as printed in the paper.""" + psi = np.arctan2(Z, X) + x = np.sin(psi) * np.sin(epsilon) + t = r / np.sqrt(X ** 2 + Z ** 2) + series = sum(eval_legendre(n, x) * t ** n for n in range(n_terms)) + return np.cos(psi) ** 2 * series ** 2 + + +# (r, epsilon, X, Z) with t = r/sqrt(X^2+Z^2) <= 0.6 +SERIES_POINTS = [ + (0.05, 0.0, 1.0, 0.5), + (0.10, 0.7, 1.0, 0.5), + (0.10, -1.2, 1.0, 0.5), + (0.30, 1.5, 0.8, 0.6), + (0.60, -0.4, 1.0, 0.2), + (0.05, 0.3, 0.1, 0.05), + (0.50, 1.0, 2.0, 1.5), +] + + +@pytest.mark.parametrize("r,epsilon,X,Z", SERIES_POINTS) +def test_closed_form_matches_legendre_series(r, epsilon, X, Z): + """Closed form == 50-term truncation of the printed Eq. 9 series.""" + assert get_Q_full(r, epsilon, X, Z) == pytest.approx( + q_series(r, epsilon, X, Z), rel=1e-9) + + +def test_centerline_reduction(): + """On the centerline (Z = 0), Q reduces to X^2/(X^2 + r^2) for any + epsilon (needed by the Eq. 21 centerline temperature quadrature).""" + X = 0.7 + for r in [0.0, 0.05, 0.1]: + for epsilon in [-1.0, 0.0, 1.3]: + assert get_Q_full(r, epsilon, X, 0.0) == pytest.approx( + X ** 2 / (X ** 2 + r ** 2), rel=1e-14) + + +def test_q_bounded_between_zero_and_one(): + """0 < Q <= 1 everywhere in the integration domain: the denominator + is X^2 + (Z - r sin(eps))^2 + (r cos(eps))^2 >= X^2. This bound is + what makes the exp(-S_0^2*(1-Q)) combination overflow-safe, so it + is asserted over a grid that includes the worst case Z ~ r*sin(eps).""" + r = np.linspace(0.0, 0.1, 41) + epsilon = np.linspace(-np.pi / 2, np.pi / 2, 41) + R, E = np.meshgrid(r, epsilon) + for X, Z in [(1.0, 0.5), (0.05, 0.08), (0.01, 0.0), (0.3, 3.0)]: + Q = get_Q_full(R, E, X, Z) + assert np.all(Q > 0.0) + assert np.all(Q <= 1.0) diff --git a/tests/plume/plume_unit_test_03.py b/tests/plume/plume_unit_test_03.py new file mode 100644 index 0000000..a9d48de --- /dev/null +++ b/tests/plume/plume_unit_test_03.py @@ -0,0 +1,136 @@ +# ======================== +# PyRPOD: tests/plume/plume_unit_test_03.py +# ======================== +# Tests for the Simons cosine-law plume model (Cai & Wang 2012, Sec. II.C, +# Eqs. 25-26) after the gamma-generalization fixes: +# - isentropic throat ratios computed from gamma instead of the +# hardcoded gamma=1.4 table values 0.5283 / 0.8333, +# - theta >= theta_max returns an empty-region density of 0.0, +# - parameterizable beaming exponent kappa. + +import numpy as np +import pytest + +from pyrpod.plume.RarefiedPlumeGasKinetics import Simons + +pytestmark = [pytest.mark.plume, pytest.mark.unit] + +R_SPECIFIC = 287.0 # J / (kg K) +T_C = 500.0 # K +P_C = 1.0e5 # N / m^2 +R_0 = 0.1 # m +R_FIELD = 2.0 # m + + +def make_simons(gamma, kappa=None): + return Simons(gamma, R_SPECIFIC, T_C, P_C, R_0, R_FIELD, kappa=kappa) + + +def test_throat_ratios_reproduce_gamma_1_4_table_values(): + """The closed forms T*/T_c = 2/(gamma+1), P*/P_c = (2/(gamma+1))^ + (gamma/(gamma-1)) must reproduce the isentropic-table values + 0.8333 / 0.5283 for gamma = 1.4 to 4 decimal places.""" + gamma = 1.4 + t_ratio = 2 / (gamma + 1) + p_ratio = t_ratio ** (gamma / (gamma - 1)) + assert t_ratio == pytest.approx(0.8333, abs=5e-5) + assert p_ratio == pytest.approx(0.5283, abs=5e-5) + + # And the class must use them: density/sonic velocity at gamma=1.4 + # match the old hardcoded-table implementation to table precision. + simons = make_simons(gamma) + rho_old_table = (P_C * 0.5283) / (R_SPECIFIC * T_C * 0.8333) + assert simons.get_nozzle_throat_density() == pytest.approx(rho_old_table, + rel=1e-4) + a_old_table = np.sqrt(gamma * R_SPECIFIC * T_C * 0.8333) + assert simons.get_sonic_velocity() == pytest.approx(a_old_table, rel=1e-4) + + +def test_gamma_2_yields_real_finite_normalization(): + """Regression for the old 'broken when gamma = 2' failure: the + numerically integrated normalization constant is a plain real, + positive, finite float for fractional/even kappa alike.""" + simons = make_simons(2.0) + assert isinstance(simons.A, float) + assert np.isfinite(simons.A) + assert simons.A > 0.0 + + +@pytest.mark.parametrize("gamma", [1.4, 5.0 / 3.0, 2.0]) +def test_density_zero_beyond_limiting_angle(gamma): + """theta >= theta_max is a physically empty region: both the decay + function and the density ratio must return exactly 0.0 (previously + the cosine went negative and fractional kappa produced complex + numbers).""" + simons = make_simons(gamma) + theta_max = simons.get_limiting_turn_angle() + for theta in [theta_max, 1.01 * theta_max, 2.0 * theta_max]: + f = simons.get_plume_angular_density_decay_function(theta) + n = simons.get_num_density_ratio(theta) + assert f == 0.0 + assert n == 0.0 + assert not isinstance(f, complex) + assert not isinstance(n, complex) + + +def test_density_positive_inside_limiting_angle(): + simons = make_simons(2.0) + theta_max = simons.get_limiting_turn_angle() + n = simons.get_num_density_ratio(0.99 * theta_max) + assert np.isreal(n) + assert 0.0 < n < simons.get_num_density_ratio(0.0) + + +def test_default_kappa_is_boyton(): + """Default beaming exponent is Boyton's kappa = 2/(gamma-1) + [Cai & Wang 2012, Sec. II.C, refs. 22-23].""" + gamma = 1.4 + simons = make_simons(gamma) + assert simons.kappa == pytest.approx(2 / (gamma - 1)) + + +@pytest.mark.parametrize("kappa", [2.0, 2.5, 1 / 0.4]) +def test_kappa_parameterization(kappa): + """Custom kappa (e.g. Ashkenas & Sherman's 2, Albini's 1/(gamma-1)) + is used by the decay function: f = cos^kappa at the half angle, + and f(0) = 1 for every kappa (paper Sec. II.C).""" + simons = make_simons(1.4, kappa=kappa) + assert simons.kappa == kappa + assert simons.get_plume_angular_density_decay_function(0.0) == 1.0 + theta_max = simons.get_limiting_turn_angle() + theta = 0.5 * theta_max + expected = np.cos((np.pi / 2) * (theta / theta_max)) ** kappa + f = simons.get_plume_angular_density_decay_function(theta) + assert f == pytest.approx(expected, rel=1e-12) + + +def test_kappa_changes_normalization_constant(): + """A comes from mass-flow continuity over the chosen decay function, + so it must respond to kappa: a larger kappa concentrates the plume, + requiring a larger A.""" + a_kappa_2 = make_simons(1.4, kappa=2.0).A + a_kappa_5 = make_simons(1.4, kappa=5.0).A + assert a_kappa_5 > a_kappa_2 > 0.0 + + +def test_exit_referenced_density_at_mach_1_equals_throat_referenced(): + """With exit Mach 1 the exit IS the throat, so n/n_0 == n/n_s.""" + simons = make_simons(1.4) + theta = 0.3 + assert simons.get_num_density_ratio_exit(theta, 1.0) == pytest.approx( + simons.get_num_density_ratio(theta), rel=1e-12) + + +def test_exit_referenced_density_scaling(): + """n/n_0 = (n/n_s) * (n_s/n_0) with the isentropic throat-to-exit + density ratio [(1+(gamma-1)/2*Me^2)/((gamma+1)/2)]^(1/(gamma-1)). + The throat is denser than a supersonic exit, so the factor is > 1.""" + gamma = 1.4 + exit_mach = 5.0 + simons = make_simons(gamma) + theta = 0.3 + factor = ((1 + (gamma - 1) / 2 * exit_mach ** 2) + / ((gamma + 1) / 2)) ** (1 / (gamma - 1)) + assert factor > 1.0 + assert simons.get_num_density_ratio_exit(theta, exit_mach) == pytest.approx( + simons.get_num_density_ratio(theta) * factor, rel=1e-12) diff --git a/tests/plume/plume_verification_error_summary.py b/tests/plume/plume_verification_error_summary.py new file mode 100644 index 0000000..b99b6d7 --- /dev/null +++ b/tests/plume/plume_verification_error_summary.py @@ -0,0 +1,111 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_error_summary.py +# ======================== +# Model-vs-model analog of Cai & Wang 2012 Table 1: maximum relative +# differences in density and mass flux between the analytical, +# simplified and Simons (exit-referenced, Boyton kappa) plume models +# along the centerline and along the angular curve r/D = 10, S0 = 2.0. +# +# Table 1's reference columns are DSMC results, which PyRPOD does not +# have yet; those columns are emitted as "pending digitized data" and +# will be filled from tests/plume/data/digitized/ in a follow-up. As +# the paper notes, "the largest density relative error actually +# happens within the curve tails with very small values" -- each row +# therefore also reports the restricted-range maximum. +# +# The filename intentionally avoids pytest's collection patterns +# (test_*.py / *_test_*.py): this is a manual-run generator -- +# python tests/plume/plume_verification_error_summary.py +# It writes tests/plume/output/model_error_summary.csv and .md. + +import numpy as np + +import plume_figure_utils as u + +S_0 = 2.0 + + +def _row(curve, quantity, comparison, values, reference, mask, mask_note): + full = u.max_rel_diff(values, reference) + restricted = u.max_rel_diff(values[mask], reference[mask]) + return {'curve': curve, 'quantity': quantity, 'comparison': comparison, + 'max_rel_diff': full, + 'max_rel_diff_restricted': restricted, + 'restriction': mask_note, + 'vs_DSMC': 'pending digitized data'} + + +def build_rows(): + rows = [] + + # centerline, X/D in [X_MIN, 10]. As in the paper's Figs. 19-21, + # "simplified" on the centerline is the Eq. 18 closed form (paper + # Sec. II.B) and "analytical" is the Eq. 5 integral; they coincide + # to quadrature accuracy. + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + r = x_over_D * u.D_NOZZLE + n_analytical = u.full_field_values('n', r, np.zeros_like(r), S_0) + n_simplified = u.centerline_profiles(x_over_D, S_0, ('n',))['n'] + n_simons = u.simons_density_field(r, np.zeros_like(r), S_0) + core = x_over_D >= 1.0 + rows.append(_row('centerline', 'density', + 'simplified (Eq. 18) vs analytical (Eq. 5)', + n_simplified, n_analytical, core, 'X/D >= 1')) + rows.append(_row('centerline', 'density', 'Simons vs analytical', + n_simons, n_analytical, core, 'X/D >= 1')) + + # angular curve r/D = 10 (simplified series delegates to + # SimplifiedGasKinetics via plume_figure_utils) + theta = np.deg2rad(np.linspace(0.0, u.THETA_MAX_DEG, 90)) + X = 10.0 * u.D_NOZZLE * np.cos(theta) + Z = 10.0 * u.D_NOZZLE * np.sin(theta) + nA = u.full_field_values('n', X, Z, S_0) + nS = u.simplified_field_values('n', X, Z, S_0) + nSim = u.simons_density_field(X, Z, S_0) + # mass flux with the same Fig. 25 normalization, n * Vr*sqrt(beta0) + fluxA = nA * u.full_field_values('Vr', X, Z, S_0) + fluxS = nS * u.simplified_field_values('Vr', X, Z, S_0) + core = theta <= np.deg2rad(60.0) + rows.append(_row('angular r/D=10', 'density', 'simplified vs analytical', + nS, nA, core, 'theta <= 60 deg')) + rows.append(_row('angular r/D=10', 'density', 'Simons vs analytical', + nSim, nA, core, 'theta <= 60 deg')) + rows.append(_row('angular r/D=10', 'mass flux', 'simplified vs analytical', + fluxS, fluxA, core, 'theta <= 60 deg')) + return rows + + +def write_summary(): + rows = build_rows() + u.OUTPUT_DIR.mkdir(parents=True, exist_ok=True) + columns = ['curve', 'quantity', 'comparison', 'max_rel_diff', + 'max_rel_diff_restricted', 'restriction', 'vs_DSMC'] + + csv_path = u.OUTPUT_DIR / 'model_error_summary.csv' + with open(csv_path, 'w', encoding='utf-8', newline='') as fh: + fh.write(','.join(columns) + '\n') + for row in rows: + fh.write(','.join( + f'{row[c]:.4g}' if isinstance(row[c], float) else str(row[c]) + for c in columns) + '\n') + + md_path = u.OUTPUT_DIR / 'model_error_summary.md' + with open(md_path, 'w', encoding='utf-8', newline='') as fh: + fh.write('# Model-vs-model error summary (Table 1 analog)\n\n') + fh.write(f'Cai & Wang 2012 conditions, S0 = {S_0}. DSMC reference ' + 'columns are pending digitized data.\n\n') + fh.write('| ' + ' | '.join(columns) + ' |\n') + fh.write('|' + '---|' * len(columns) + '\n') + for row in rows: + fh.write('| ' + ' | '.join( + f'{row[c]:.4g}' if isinstance(row[c], float) else str(row[c]) + for c in columns) + ' |\n') + fh.write('\nNote: unrestricted maxima occur in the plume tails ' + 'where the reference density is very small (the paper ' + 'makes the same observation about Table 1).\n') + return csv_path, md_path + + +if __name__ == '__main__': + for path in write_summary(): + print(f'wrote {path}') diff --git a/tests/plume/plume_verification_test_02.py b/tests/plume/plume_verification_test_02.py new file mode 100644 index 0000000..54a7d61 --- /dev/null +++ b/tests/plume/plume_verification_test_02.py @@ -0,0 +1,214 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_02.py +# ======================== +# Pinning tests for the SimplifiedGasKinetics class (Cai & Wang 2012, +# "Numerical Validations for a Set of Collisionless Rocket Plume Solutions", +# JSR 49(1), DOI 10.2514/1.A32046). +# +# Reference values were computed independently from the paper's equations +# with mpmath at 40 significant digits (Eqs. 13-19 and 22-24 implemented +# directly from the printed formulas, not from pyrpod code). They pin the +# known-correct behavior of the simplified far-field model so that later +# refactors (sympy removal, overflow-safe exponentials) can be verified +# to preserve behavior. + +import numpy as np +import pytest + +from pyrpod.plume.RarefiedPlumeGasKinetics import ( + SimplifiedGasKinetics, + get_far_field_temp_ratio, + get_far_field_velocity_normalized, +) + +pytestmark = [pytest.mark.plume, pytest.mark.verification] + +# Fixed thruster gas properties used for all cases (argon-like). +R_SPECIFIC = 208.13 # J / (kg K) +T_0 = 500.0 # K +GAMMA = 5.0 / 3.0 +N_0 = 1.0e20 # m^-3 +D_NOZZLE = 0.2 # m (R_0 = 0.1 m) +T_W = 300.0 # K +SIGMA = 1.0 + +# Relative tolerance for pinned closed-form values: the implementation +# evaluates the same closed-form expressions in float64 (machine epsilon +# ~1e-16, tens of flops => accumulated rounding well below 1e-12), while +# the references carry 17 correct digits. 1e-9 leaves ample headroom for +# rounding-level refactors (e.g. sympy.erf -> scipy.special.erf, combined +# exponentials) yet fails loudly on any actual formula change. +RTOL_PINNED = 1e-9 + + +def make_plume(distance, theta, S_0): + """Build a SimplifiedGasKinetics instance with exit speed ratio S_0.""" + ve = S_0 * np.sqrt(2 * R_SPECIFIC * T_0) + thruster_characteristics = { + 'd': D_NOZZLE, 've': ve, 'R': R_SPECIFIC, + 'gamma': GAMMA, 'Te': T_0, 'n': N_0, + } + return SimplifiedGasKinetics(distance, theta, thruster_characteristics, + T_W, SIGMA) + + +# (distance, theta, S_0) -> (n/n0, U*sqrt(beta0), W*sqrt(beta0), T/T0) +# per Eqs. 14-17 with Q' = X^2/(X^2+Z^2) (Eq. 13). +FIELD_CASES = [ + ((1.0, 0.2, 2.0), (0.036338635060152093, 2.3636427215997389, + 0.47913410002529834, 0.27393936134733628)), + ((2.0, 0.7, 2.0), (0.0010317934247266036, 1.5832262838600719, + 1.3335331025390796, 0.2543577108315116)), + ((5.0, 1.2, 2.0), (4.4418364832646119e-6, 0.54854957955069592, + 1.4109526908580327, 0.20360334802441553)), + ((1.5, 0.4, 1.0), (0.0046772248683472122, 1.5085528987977848, + 0.63780593571949983, 0.21734975912503903)), + ((3.0, 0.9, 3.0), (1.0977476844642729e-5, 1.4507999313815155, + 1.8282374555518849, 0.27008477291007084)), +] + +# (X, S_0) -> (n/n0, U*sqrt(beta0)) per Eqs. 18-19 (centerline). +CENTERLINE_CASES = [ + ((0.5, 2.0), (0.15923868831612018, 2.4062627987321269)), + ((1.0, 2.0), (0.043598168447392599, 2.434502010083779)), + ((5.0, 2.0), (0.0017976259018792059, 2.4441625633761652)), + ((1.0, 1.0), (0.014624592798873514, 1.6841200399767633)), + ((1.0, 3.0), (0.089793762869630585, 3.3010882334220839)), +] + +# S_0 -> (U_inf*sqrt(beta0), T_inf/T0) per Eqs. 22-24, evaluated with +# mpmath. (Eq. 23's denominator is 3*[S0 + (1/2+S0^2)*sqrt(pi)* +# (1+erf(S0))*exp(S0^2)], i.e. 3*K(Q=1), consistent with +# -2/3*G^2 + 4N(1)/(3K(1)).) +ASYMPTOTE_CASES = { + 1.0: (1.689948557881606, 0.22268161973032418), + 2.0: (2.444572023854802, 0.2754744452640276), + 3.0: (3.3157896665554836, 0.30193859122021722), +} + + +@pytest.mark.parametrize("point,expected", FIELD_CASES) +def test_field_solution_pinned(point, expected): + """Pin Eqs. 14-17 (off-centerline simplified field solutions).""" + distance, theta, S_0 = point + n_ref, U_ref, W_ref, T_ref = expected + plume = make_plume(distance, theta, S_0) + assert plume.get_num_density_ratio() == pytest.approx(n_ref, rel=RTOL_PINNED) + assert plume.get_U_normalized() == pytest.approx(U_ref, rel=RTOL_PINNED) + assert plume.get_W_normalized() == pytest.approx(W_ref, rel=RTOL_PINNED) + assert plume.get_temp_ratio() == pytest.approx(T_ref, rel=RTOL_PINNED) + + +@pytest.mark.parametrize("point,expected", CENTERLINE_CASES) +def test_centerline_solution_pinned(point, expected): + """Pin Eqs. 18-19 (centerline density and velocity closed forms).""" + X, S_0 = point + n_ref, U_ref = expected + plume = make_plume(X, 0.0, S_0) + assert plume.get_num_density_centerline() == pytest.approx(n_ref, + rel=RTOL_PINNED) + assert plume.get_velocity_centerline() == pytest.approx(U_ref, + rel=RTOL_PINNED) + + +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_far_field_asymptote_functions_pinned(S_0): + """Pin the Eq. 22-24 module-level asymptote functions against the + independent mpmath evaluation (tolerance rationale as RTOL_PINNED).""" + U_ref, T_ref = ASYMPTOTE_CASES[S_0] + assert get_far_field_velocity_normalized(S_0) == pytest.approx( + U_ref, rel=RTOL_PINNED) + assert get_far_field_temp_ratio(S_0) == pytest.approx( + T_ref, rel=RTOL_PINNED) + + +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_centerline_velocity_approaches_asymptote(S_0): + """Centerline U converges to the Eq. 22/24 far-field asymptote. + + Convergence is O((R_0/X)^2): at X/R_0 = 1000 the residual is ~4e-7 + relative, so rel=1e-5 passes with an order-of-magnitude margin while + still requiring genuine convergence. + """ + U_inf = get_far_field_velocity_normalized(S_0) + X = 1000 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, S_0) + assert plume.get_velocity_centerline() == pytest.approx(U_inf, rel=1e-5) + + +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_centerline_temperature_approaches_asymptote(S_0): + """Centerline T converges to the Eq. 23 far-field asymptote. + + T = -2/3*U^2 + 4/3*/ subtracts two O(S_0^2) quantities, so the + O((R_0/X)^2 * S_0^2) residual of each term is amplified by ~U^2/T + (a factor ~25 at S_0 = 3) in the relative error of T. At + X/R_0 = 5000 the amplified residual is <~1e-5 for S_0 <= 3, so + rel=1e-4 passes with margin for both the constant-N approximation + and the exact Eq. 21 quadrature, while still requiring convergence. + """ + T_inf = get_far_field_temp_ratio(S_0) + X = 5000 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, S_0) + assert plume.get_temp_centerline() == pytest.approx(T_inf, rel=1e-4) + + +def old_temp_centerline_approximation(plume): + """The pre-fix Eq. 21 evaluation: N treated as constant across the + exit radius, integral ~ 0.5 * N(Q'=1) * R_0^2 (N_simple carries the + exp(-S_0^2) prefactor).""" + integral = 0.5 * plume.N_simple * plume.R_0 ** 2 + return (4 / (3 * plume.get_num_density_centerline() * np.sqrt(np.pi) + * plume.X ** 2) * integral + - plume.get_velocity_centerline() ** 2 / 1.5) + + +def test_temp_centerline_quadrature_matches_old_in_far_field(): + """(a) Far field: with X >> R_0, Q(r) = X^2/(X^2+r^2) ~ 1 across the + whole exit, so the exact quadrature and the constant-N approximation + must agree. Measured relative difference at X/R_0 = 1000 is ~5e-5 + (it scales as (R_0/X)^2 amplified by the -2/3 U^2 cancellation, see + the asymptote test above); rel=1e-3 gives ~20x margin.""" + X = 1000 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, 2.0) + assert plume.get_temp_centerline() == pytest.approx( + old_temp_centerline_approximation(plume), rel=1e-3) + + +def test_temp_centerline_quadrature_diverges_from_old_in_near_field(): + """(b) Near field: at X = 2*R_0 the exact Q(r) varies strongly over + the exit (Q(R_0) = 0.8), and the old approximation errs by ~87%, + yielding an unphysical T > T_0 (collisionless expansion can only + cool the gas below the exit temperature -- paper Fig. 13 argument). + Require a divergence of at least 50% to prove the fix changes the + near field meaningfully.""" + X = 2 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, 2.0) + T_new = plume.get_temp_centerline() + T_old = old_temp_centerline_approximation(plume) + assert abs(T_new / T_old - 1) > 0.5 + assert T_new < 1.0 + assert T_old > 1.0 + + +@pytest.mark.parametrize("point,expected", [ + ((1.0, 0.2, 2.0), (0.36855461276327334, -0.056815859301478394, + 143.3923778236429)), + ((2.0, 0.7, 2.0), (0.005299929030826453, -0.003007472717182616, + 1.9141355894449432)), +]) +def test_surface_interaction_pinned(point, expected): + """Pin the Maxwell gas-surface interface outputs (off-centerline). + + Unlike the field-solution pins above, these reference values were + generated from the current implementation (not independently), so + they pin implementation behavior for refactor safety rather than + correctness against the paper. Tolerance rationale as RTOL_PINNED. + """ + distance, theta, S_0 = point + p_ref, tau_ref, q_ref = expected + plume = make_plume(distance, theta, S_0) + assert float(plume.get_pressure()) == pytest.approx(p_ref, rel=RTOL_PINNED) + assert float(plume.get_shear_pressure()) == pytest.approx(tau_ref, + rel=RTOL_PINNED) + assert float(plume.get_heat_flux()) == pytest.approx(q_ref, rel=RTOL_PINNED) diff --git a/tests/plume/plume_verification_test_03.py b/tests/plume/plume_verification_test_03.py new file mode 100644 index 0000000..fa3f5ba --- /dev/null +++ b/tests/plume/plume_verification_test_03.py @@ -0,0 +1,166 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_03.py +# ======================== +# Verification tests for the full collisionless analytical plume model +# CollisionlessGasKinetics (Cai & Wang 2012, Eqs. 5-12), checked against +# the paper's own internal anchors: +# (a, b) the exact field integrals must reduce to the Eq. 18/19 +# closed forms on the centerline, +# (c) W must vanish on the centerline (Eq. 20), +# (d) the full and simplified models must agree in the far field +# within a bound derived from Eq. 30, +# (e) qualitative physics of Figs. 3, 5, 13 (monotonic centerline +# density decay, cooling below T_0), +# (f) convergence to the Eq. 22-24 far-field asymptotes. + +import numpy as np +import pytest + +from pyrpod.plume.RarefiedPlumeGasKinetics import ( + CollisionlessGasKinetics, + SimplifiedGasKinetics, + get_far_field_temp_ratio, + get_far_field_velocity_normalized, +) + +pytestmark = [pytest.mark.plume, pytest.mark.verification] + +R_SPECIFIC = 208.13 # J / (kg K), argon-like +T_0 = 500.0 # K +GAMMA = 5.0 / 3.0 +N_0 = 1.0e20 # m^-3 +D_NOZZLE = 0.2 # m (R_0 = 0.1 m) +R_0 = D_NOZZLE / 2 +T_W = 300.0 # K +SIGMA = 1.0 + +# The class doubles the quadrature order until successive orders agree +# to QUAD_RTOL = 1e-9; comparisons against exact closed forms therefore +# use a slightly looser 1e-8. +RTOL_QUAD = 1e-8 + + +def make_plume(cls, distance, theta, S_0): + ve = S_0 * np.sqrt(2 * R_SPECIFIC * T_0) + thruster_characteristics = { + 'd': D_NOZZLE, 've': ve, 'R': R_SPECIFIC, + 'gamma': GAMMA, 'Te': T_0, 'n': N_0, + } + return cls(distance, theta, thruster_characteristics, T_W, SIGMA) + + +CENTERLINE_POINTS = [(0.2, 1.0), (0.5, 2.0), (1.0, 2.0), (5.0, 2.0), + (1.0, 3.0)] + + +@pytest.mark.parametrize("X,S_0", CENTERLINE_POINTS) +def test_centerline_density_matches_eq18(X, S_0): + """(a) The Eq. 5 double integral evaluated at (X, 0, 0) must equal + the Eq. 18 closed form (they were derived from the same solution).""" + plume = make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + assert plume.get_num_density_ratio() == pytest.approx( + plume.get_num_density_centerline(), rel=RTOL_QUAD) + + +@pytest.mark.parametrize("X,S_0", CENTERLINE_POINTS) +def test_centerline_velocity_matches_eq19(X, S_0): + """(b) The Eq. 6 integral at (X, 0, 0) must equal the Eq. 19 closed + form.""" + plume = make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + assert plume.get_U_normalized() == pytest.approx( + plume.get_velocity_centerline(), rel=RTOL_QUAD) + + +@pytest.mark.parametrize("X,S_0", CENTERLINE_POINTS) +def test_centerline_W_vanishes(X, S_0): + """(c) Eq. 20: W(X, 0, 0) = 0. The quadrature can only leave + rounding noise, so |W| is required to be negligible against U. On + the centerline V_r must also reduce to U.""" + plume = make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + U = plume.get_U_normalized() + assert abs(plume.get_W_normalized()) < 1e-12 * U + assert plume.get_Vr_normalized() == pytest.approx(U, rel=1e-12) + + +# Far-field points (distance, theta) with S_0 = 2. +FAR_FIELD_POINTS = [(20.0, 0.4), (50.0, 0.8), (100.0, 0.2)] + + +@pytest.mark.parametrize("distance,theta", FAR_FIELD_POINTS) +def test_far_field_agrees_with_simplified_within_eq30_bound(distance, theta): + """(d) Eq. 30 bounds the pointwise error of the simplified Q': + |Q' - Q|/Q <= max(R_0^2, 2*Z*R_0)/(X^2 + Z^2) over the exit disk. + The K, M, N factors depend on Q through powers up to Q^(5/2), the + factor exp(S_0^2 * Q), and bounded erf terms, so their logarithmic + sensitivity |d ln K / d ln Q| is at most ~(S_0^2 + 4). The relative + difference between the full and simplified field quantities is + therefore bounded by (S_0^2 + 4) * max(R_0^2, 2*Z*R_0)/(X^2 + Z^2).""" + S_0 = 2.0 + full = make_plume(CollisionlessGasKinetics, distance, theta, S_0) + simple = make_plume(SimplifiedGasKinetics, distance, theta, S_0) + + bound = ((S_0 ** 2 + 4) + * max(R_0 ** 2, 2 * full.Z * R_0) / (full.X ** 2 + full.Z ** 2)) + + assert full.get_num_density_ratio() == pytest.approx( + simple.get_num_density_ratio(), rel=bound) + assert full.get_U_normalized() == pytest.approx( + simple.get_U_normalized(), rel=bound) + assert full.get_W_normalized() == pytest.approx( + simple.get_W_normalized(), rel=bound) + assert full.get_temp_ratio() == pytest.approx( + simple.get_temp_ratio(), rel=bound) + + +def test_centerline_density_decays_monotonically(): + """(e) Paper Fig. 3: for fixed S_0 the centerline density decreases + monotonically with X.""" + S_0 = 2.0 + X_values = [0.15, 0.3, 0.6, 1.2, 2.5, 5.0, 10.0, 20.0] + densities = [make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + .get_num_density_ratio() for X in X_values] + assert all(a > b for a, b in zip(densities, densities[1:])) + assert all(n > 0 for n in densities) + + +@pytest.mark.parametrize("distance,theta", [(0.2, 0.0), (1.0, 0.0), + (1.0, 0.5), (5.0, 1.0)]) +def test_temperature_below_exit_temperature_downstream(distance, theta): + """(e) Paper Figs. 5, 13: the collisionless expansion only cools the + gas, so T_1 < T_0 everywhere downstream; consequently the static + pressure ratio p_1/p_0 = (n/n_0)(T/T_0) lies below the density + ratio (the paper's p. 64 argument against using n*k*T_0).""" + plume = make_plume(CollisionlessGasKinetics, distance, theta, 2.0) + T_ratio = plume.get_temp_ratio() + assert 0.0 < T_ratio < 1.0 + assert 0.0 < plume.get_pressure_ratio() < plume.get_num_density_ratio() + + +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_centerline_converges_to_far_field_asymptotes(S_0): + """(f) Full-model centerline U and T converge to the Eq. 22-24 + asymptotes at large X. Tolerances follow the O((R_0/X)^2) + convergence-rate arguments in plume_verification_test_02.py + (rel=1e-5 for U at X/R_0 = 1000; rel=1e-4 for T at X/R_0 = 5000 + due to the -2/3 U^2 cancellation amplification).""" + plume_U = make_plume(CollisionlessGasKinetics, 1000 * R_0, 0.0, S_0) + assert plume_U.get_U_normalized() == pytest.approx( + get_far_field_velocity_normalized(S_0), rel=1e-5) + + plume_T = make_plume(CollisionlessGasKinetics, 5000 * R_0, 0.0, S_0) + assert plume_T.get_temp_ratio() == pytest.approx( + get_far_field_temp_ratio(S_0), rel=1e-4) + + +def test_surface_interaction_uses_full_model_and_is_finite(): + """The inherited Maxwell gas-surface interface dispatches to the + overridden full-model field getters; off the centerline it must + return finite, physically-signed values (p > 0, q finite, shear + opposing the flow direction).""" + plume = make_plume(CollisionlessGasKinetics, 1.0, 0.4, 2.0) + p = float(plume.get_pressure()) + tau = float(plume.get_shear_pressure()) + q = float(plume.get_heat_flux()) + assert np.isfinite(p) and p > 0.0 + assert np.isfinite(tau) and tau < 0.0 + assert np.isfinite(q) diff --git a/tests/plume/plume_verification_test_04.py b/tests/plume/plume_verification_test_04.py new file mode 100644 index 0000000..5978a31 --- /dev/null +++ b/tests/plume/plume_verification_test_04.py @@ -0,0 +1,44 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_04.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 2: plume boundaries (n/n0 = 0.001 +# contour of the full analytical model) for exit speed ratios +# S0 = 1, 2, 3; X/D in [0, 30], Z/D in [0, 15]. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_04.py +# The figure is saved to tests/plume/output/. +# Digitized overlays looked for: fig02_*.csv (none published -- Fig. 2 +# has no DSMC data in the paper). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x = np.linspace(u.X_MIN_OVER_D, 30.0, 181) + z = np.linspace(0.0, 15.0, 91) + + fig, ax = u.plt.subplots(figsize=(7, 4.5)) + for S_0, ls in [(1.0, '-'), (2.0, '--'), (3.0, '-.')]: + n = u.full_field_grid('n', x, z, S_0) + cs = ax.contour(x, z, n, levels=[0.001], colors='k', + linestyles=ls, linewidths=1.2) + ax.clabel(cs, fmt='0.001', fontsize=7) + ax.plot([], [], 'k', ls=ls, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig02') + + ax.set_xlim(0, 30) + ax.set_ylim(0, 15) + ax.set_xlabel('X/D') + ax.set_ylabel('Z/D') + ax.set_title('Plume boundaries at different exit speed ratios ' + '(Fig. 2)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig02_plume_boundaries') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_05.py b/tests/plume/plume_verification_test_05.py new file mode 100644 index 0000000..430f282 --- /dev/null +++ b/tests/plume/plume_verification_test_05.py @@ -0,0 +1,38 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_05.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 3: normalized analytical number +# density along the centerline (Eq. 18) vs X/D for S0 = 1, 2, 3. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_05.py +# Digitized overlay looked for: fig03_dsmc.csv (the paper's S0 = 2 +# DSMC circles). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + for S_0, ls in [(1.0, '-.'), (2.0, '--'), (3.0, '-')]: + n = u.centerline_profiles(x_over_D, S_0, ('n',))['n'] + ax.plot(x_over_D, n, 'k', ls=ls, lw=1.2, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig03') + + ax.set_xlim(0, 10) + ax.set_ylim(0, 1.2) + ax.set_xlabel('X/D') + ax.set_ylabel('$n_1/n_0$') + ax.set_title('Normalized analytical number density along centerline ' + '(Fig. 3)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig03_centerline_density') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_06.py b/tests/plume/plume_verification_test_06.py new file mode 100644 index 0000000..829d5dd --- /dev/null +++ b/tests/plume/plume_verification_test_06.py @@ -0,0 +1,37 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_06.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 4: normalized analytical U-velocity +# along the centerline (Eq. 19), U1*sqrt(beta0) vs X/D, S0 = 1, 2, 3. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_06.py +# Digitized overlay looked for: fig04_dsmc.csv (S0 = 2 DSMC circles). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + for S_0, ls in [(1.0, '-.'), (2.0, '--'), (3.0, '-')]: + U = u.centerline_profiles(x_over_D, S_0, ('U',))['U'] + ax.plot(x_over_D, U, 'k', ls=ls, lw=1.2, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig04') + + ax.set_xlim(0, 10) + ax.set_ylim(1, 4) + ax.set_xlabel('X/D') + ax.set_ylabel(r'$U_1\sqrt{\beta_0}$') + ax.set_title('Normalized analytical U-velocity along centerline ' + '(Fig. 4)') + ax.legend(fontsize=8, loc='lower right') + return u.save_figure(fig, 'fig04_centerline_velocity') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_07.py b/tests/plume/plume_verification_test_07.py new file mode 100644 index 0000000..3ca4a96 --- /dev/null +++ b/tests/plume/plume_verification_test_07.py @@ -0,0 +1,38 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_07.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 5: normalized analytical temperature +# along the centerline (Eq. 21, exact quadrature), T1/T0 vs X/D for +# S0 = 1, 2, 3. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_07.py +# Digitized overlay looked for: fig05_dsmc.csv (S0 = 2 DSMC circles). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + for S_0, ls in [(1.0, '-.'), (2.0, '--'), (3.0, '-')]: + T = u.centerline_profiles(x_over_D, S_0, ('T',))['T'] + ax.plot(x_over_D, T, 'k', ls=ls, lw=1.2, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig05') + + ax.set_xlim(0, 10) + ax.set_ylim(0.2, 1.2) + ax.set_xlabel('X/D') + ax.set_ylabel('$T_1/T_0$') + ax.set_title('Normalized analytical temperature along centerline ' + '(Fig. 5)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig05_centerline_temperature') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_08.py b/tests/plume/plume_verification_test_08.py new file mode 100644 index 0000000..bf5baec --- /dev/null +++ b/tests/plume/plume_verification_test_08.py @@ -0,0 +1,23 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_08.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 6: normalized number-density contours, +# Kn = 100, S0 = 2.0. Upper half-plane: analytical (solid) + simplified +# (dashed); lower half-plane: Simons (dashed) + DSMC slot. The analytic +# curves are Kn-independent (collisionless); Figs. 6-8 differ only in +# their DSMC data, so each gets its own script/overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_08.py +# Digitized overlays looked for: fig06_*.csv (e.g. fig06_dsmc_0p001.csv, +# one CSV per digitized DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.density_contour_figure('100', 'fig06') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_09.py b/tests/plume/plume_verification_test_09.py new file mode 100644 index 0000000..5dae292 --- /dev/null +++ b/tests/plume/plume_verification_test_09.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_09.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 7: normalized number-density contours, +# Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 6 (collisionless +# models are Kn-independent); the Kn distinction lives in the DSMC +# overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_09.py +# Digitized overlays looked for: fig07_*.csv (one CSV per digitized +# DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.density_contour_figure('0.1', 'fig07') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_10.py b/tests/plume/plume_verification_test_10.py new file mode 100644 index 0000000..f5ab6b3 --- /dev/null +++ b/tests/plume/plume_verification_test_10.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_10.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 8: normalized number-density contours, +# Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 6 (collisionless +# models are Kn-independent); the Kn distinction lives in the DSMC +# overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_10.py +# Digitized overlays looked for: fig08_*.csv (one CSV per digitized +# DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.density_contour_figure('0.01', 'fig08') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_11.py b/tests/plume/plume_verification_test_11.py new file mode 100644 index 0000000..deff98a --- /dev/null +++ b/tests/plume/plume_verification_test_11.py @@ -0,0 +1,50 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_11.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 9: relative density error between the +# analytical and Simons (cosine-law) solutions, |n_A/n_Simons - 1| +# [Eq. 28], Kn = 100 (analytic content is Kn-independent), S0 = 2.0, +# X/D and Y/D in [0, 10]. The Simons field is exit-referenced with the +# Boyton kappa (see plume_figure_utils docstring); the paper's A +# normalization additionally depends on exit Mach number, so absolute +# error magnitudes may differ somewhat from the printed figure. +# +# The paper notes the upper-nozzle-lip region is a density singularity +# whose error pattern "shall be neglected". +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_11.py +# Digitized overlays looked for: fig09_*.csv. + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + S_0 = 2.0 + x, y = u.default_contour_axes() + n_full = u.full_field_grid('n', x, y, S_0) + Xg, Yg = np.meshgrid(x * u.D_NOZZLE, y * u.D_NOZZLE) + n_simons = u.simons_density_field(Xg, Yg, S_0) + error = np.abs(n_full / n_simons - 1) + + fig, ax = u.plt.subplots(figsize=(6, 5)) + levels = [0.1, 0.4, 1, 2, 4, 5] + cs = ax.contour(x, y, error, levels=levels, colors='k', linewidths=1.0) + ax.clabel(cs, fmt='%g', fontsize=7) + u.overlay_digitized(ax, 'fig09', style='line') + + ax.set_xlim(0, 10) + ax.set_ylim(0, 10) + ax.set_xlabel('X/D') + ax.set_ylabel('Y/D') + ax.set_title('Relative density error: analytical vs Simons (Fig. 9)') + u.annotate_error( + ax, f'max |n_A/n_Sim - 1| = {np.nanmax(error):.2g} ' + '(peak at the nozzle-lip singularity)', loc='lower right') + return u.save_figure(fig, 'fig09_error_analytical_vs_simons') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_12.py b/tests/plume/plume_verification_test_12.py new file mode 100644 index 0000000..1e2f62f --- /dev/null +++ b/tests/plume/plume_verification_test_12.py @@ -0,0 +1,47 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_12.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 10: relative density error between the +# analytical and simplified analytical solutions, |n_A/n_As - 1| +# [Eq. 29], Kn = 100 (analytic content is Kn-independent), S0 = 2.0, +# X/D and Y/D in [0, 10]. +# +# The paper notes the upper-nozzle-lip region is a density singularity +# whose error pattern "shall be neglected". +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_12.py +# Digitized overlays looked for: fig10_*.csv. + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + S_0 = 2.0 + x, y = u.default_contour_axes() + n_full = u.full_field_grid('n', x, y, S_0) + n_simp = u.simplified_field_grid('n', x, y, S_0) + error = np.abs(n_full / n_simp - 1) + + fig, ax = u.plt.subplots(figsize=(6, 5)) + levels = [0.005, 0.03, 0.05, 0.1, 0.2, 0.5, 1, 2] + cs = ax.contour(x, y, error, levels=levels, colors='k', linewidths=1.0) + ax.clabel(cs, fmt='%g', fontsize=7) + u.overlay_digitized(ax, 'fig10', style='line') + + ax.set_xlim(0, 10) + ax.set_ylim(0, 10) + ax.set_xlabel('X/D') + ax.set_ylabel('Y/D') + ax.set_title('Relative density error: analytical vs simplified ' + '(Fig. 10)') + u.annotate_error( + ax, f'max |n_A/n_As - 1| = {np.nanmax(error):.2g} ' + '(peak at the nozzle-lip singularity)', loc='lower right') + return u.save_figure(fig, 'fig10_error_analytical_vs_simplified') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_13.py b/tests/plume/plume_verification_test_13.py new file mode 100644 index 0000000..09c2385 --- /dev/null +++ b/tests/plume/plume_verification_test_13.py @@ -0,0 +1,25 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_13.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 11: normalized pressure contours +# p1/p0 = (n1/n0)(T1/T0), Kn = 100, S0 = 2.0. Analytical field in the +# upper half-plane; lower half-plane reserved for the DSMC digitized +# overlay. Analytic content is Kn-independent; Figs. 11-12 differ only +# in their DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_13.py +# Digitized overlays looked for: fig11_*.csv (one CSV per digitized +# DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'p', [1e-5, 1e-4, 1e-3, 1e-2], '100', 'fig11', + 'Normalized pressure contours') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_14.py b/tests/plume/plume_verification_test_14.py new file mode 100644 index 0000000..a2a2bfe --- /dev/null +++ b/tests/plume/plume_verification_test_14.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_14.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 12: normalized pressure contours +# p1/p0, Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 11; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_14.py +# Digitized overlays looked for: fig12_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'p', [1e-5, 1e-4, 1e-3, 1e-2], '0.1', 'fig12', + 'Normalized pressure contours') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_15.py b/tests/plume/plume_verification_test_15.py new file mode 100644 index 0000000..493cf1b --- /dev/null +++ b/tests/plume/plume_verification_test_15.py @@ -0,0 +1,25 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_15.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 13: normalized temperature contours +# T1/T0, Kn = 100, S0 = 2.0; X/D in [0, 2.5], Y/D in [-2.5, 2.5]. +# Analytical field in the upper half-plane; lower half-plane reserved +# for the DSMC digitized overlay. This is the figure the paper uses to +# argue p = n*k*T0 is invalid (T1 < T0 everywhere downstream). +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_15.py +# Digitized overlays looked for: fig13_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'T', [0.3, 0.4, 0.6, 0.9], '100', 'fig13', + 'Normalized temperature distribution', + x_max=2.5, z_max=2.5, z_axis_label='Y/D') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_16.py b/tests/plume/plume_verification_test_16.py new file mode 100644 index 0000000..f65f68f --- /dev/null +++ b/tests/plume/plume_verification_test_16.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_16.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 14: normalized U-velocity contours +# U1*sqrt(beta0), Kn = 100, S0 = 2.0. Analytical field in the upper +# half-plane; lower half-plane reserved for the DSMC digitized overlay. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_16.py +# Digitized overlays looked for: fig14_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'U', [0.6, 1.0, 1.8, 2.4], '100', 'fig14', + 'Normalized U-velocity distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_17.py b/tests/plume/plume_verification_test_17.py new file mode 100644 index 0000000..2c9a9e6 --- /dev/null +++ b/tests/plume/plume_verification_test_17.py @@ -0,0 +1,25 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_17.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 15: normalized "V-velocity" contours, +# Kn = 100, S0 = 2.0. In the plotted XOZ plane (Y = 0) the y-component +# V is identically zero by axisymmetry; the transverse component shown +# in the paper's Fig. 15 corresponds to the model's W (Eq. 7), which is +# what is plotted here. Analytical field in the upper half-plane; lower +# half-plane reserved for the DSMC digitized overlay. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_17.py +# Digitized overlays looked for: fig15_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'W', [0.1, 0.5, 1.0], '100', 'fig15', + 'Normalized V-velocity distribution (transverse W at Y = 0)') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_18.py b/tests/plume/plume_verification_test_18.py new file mode 100644 index 0000000..d42e5ee --- /dev/null +++ b/tests/plume/plume_verification_test_18.py @@ -0,0 +1,26 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_18.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 16: normalized radial velocity Vr +# contours, Kn = 100, S0 = 2.0. Analytical field in the upper +# half-plane; lower half-plane reserved for the DSMC digitized overlay. +# Analytic content is Kn-independent; Figs. 16-18 differ only in their +# DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_18.py +# Digitized overlays looked for: fig16_*.csv. + +import plume_figure_utils as u + +VR_LEVELS = [1.0, 1.8, 2.0, 2.2, 2.4, 2.8] + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'Vr', VR_LEVELS, '100', 'fig16', + 'Normalized velocity $V_r$ distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_19.py b/tests/plume/plume_verification_test_19.py new file mode 100644 index 0000000..150635d --- /dev/null +++ b/tests/plume/plume_verification_test_19.py @@ -0,0 +1,24 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_19.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 17: normalized radial velocity Vr +# contours, Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 16; the +# Kn distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_19.py +# Digitized overlays looked for: fig17_*.csv. + +import plume_figure_utils as u + +from plume_verification_test_18 import VR_LEVELS + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'Vr', VR_LEVELS, '0.1', 'fig17', + 'Normalized velocity $V_r$ distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_20.py b/tests/plume/plume_verification_test_20.py new file mode 100644 index 0000000..5e8b6e7 --- /dev/null +++ b/tests/plume/plume_verification_test_20.py @@ -0,0 +1,24 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_20.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 18: normalized radial velocity Vr +# contours, Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 16; the +# Kn distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_20.py +# Digitized overlays looked for: fig18_*.csv. + +import plume_figure_utils as u + +from plume_verification_test_18 import VR_LEVELS + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'Vr', VR_LEVELS, '0.01', 'fig18', + 'Normalized velocity $V_r$ distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_21.py b/tests/plume/plume_verification_test_21.py new file mode 100644 index 0000000..ec85e17 --- /dev/null +++ b/tests/plume/plume_verification_test_21.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_21.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 19: centerline density profiles, +# Kn = 100, S0 = 2.0 -- analytical (Eq. 18), simplified (Eq. 14), +# Simons (exit-referenced cosine law) and the DSMC digitized slot. +# Analytic content is Kn-independent; Figs. 19-21 differ only in their +# DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_21.py +# Digitized overlay looked for: fig19_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.centerline_density_profile_figure('100', 'fig19') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_22.py b/tests/plume/plume_verification_test_22.py new file mode 100644 index 0000000..3d428c7 --- /dev/null +++ b/tests/plume/plume_verification_test_22.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_22.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 20: centerline density profiles, +# Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 19; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_22.py +# Digitized overlay looked for: fig20_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.centerline_density_profile_figure('0.1', 'fig20') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_23.py b/tests/plume/plume_verification_test_23.py new file mode 100644 index 0000000..2f9ff55 --- /dev/null +++ b/tests/plume/plume_verification_test_23.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_23.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 21: centerline density profiles, +# Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 19; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_23.py +# Digitized overlay looked for: fig21_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.centerline_density_profile_figure('0.01', 'fig21') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_24.py b/tests/plume/plume_verification_test_24.py new file mode 100644 index 0000000..cc3858a --- /dev/null +++ b/tests/plume/plume_verification_test_24.py @@ -0,0 +1,24 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_24.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 22: density profiles along r/D = 10, +# Kn = 100, S0 = 2.0 -- analytical, simplified, Simons kappa = 1.5/2/3 +# and the DSMC digitized slot. Following the paper (p. 65), the Simons +# curves share one normalization constant A (Boyton kappa) so they +# coincide at theta = 0; only the plotted decay exponent varies. +# Analytic content is Kn-independent; Figs. 22-24 differ only in their +# DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_24.py +# Digitized overlay looked for: fig22_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.angular_density_profile_figure('100', 'fig22') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_25.py b/tests/plume/plume_verification_test_25.py new file mode 100644 index 0000000..3932df9 --- /dev/null +++ b/tests/plume/plume_verification_test_25.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_25.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 23: density profiles along r/D = 10, +# Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 22; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_25.py +# Digitized overlay looked for: fig23_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.angular_density_profile_figure('0.1', 'fig23') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_26.py b/tests/plume/plume_verification_test_26.py new file mode 100644 index 0000000..93a2a1d --- /dev/null +++ b/tests/plume/plume_verification_test_26.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_26.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 24: density profiles along r/D = 10, +# Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 22; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_26.py +# Digitized overlay looked for: fig24_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.angular_density_profile_figure('0.01', 'fig24') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_27.py b/tests/plume/plume_verification_test_27.py new file mode 100644 index 0000000..cf7ee1a --- /dev/null +++ b/tests/plume/plume_verification_test_27.py @@ -0,0 +1,59 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_27.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 25: normalized mass flux along +# r/D = 10 vs theta, S0 = 2.0 -- the analytical curve plus three DSMC +# digitized slots (Kn = 100, 0.1, 0.01). +# +# Normalization note: the paper states the flux is "normalized by the +# inflow mass flux at the exit plane" (p. 66), which would suggest +# rho*Vr / (rho0*U0) = (n/n0)*(Vr*sqrt(beta0))/S0. That convention +# peaks at ~0.0137 at theta = 0, half the paper's ~0.027. The plotted +# magnitudes of Fig. 25 are instead reproduced exactly by +# (n/n0)*(Vr*sqrt(beta0)), i.e. rho*Vr normalized by +# rho0*sqrt(2*R*T0) = rho0*U0/S0. That convention is adopted here +# (verified: theta = 0 value 0.0273 vs the paper's ~0.027); it should +# be re-calibrated against the digitized DSMC data when available. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_27.py +# Digitized overlays looked for: fig25_dsmc_kn100.csv, +# fig25_dsmc_kn0p1.csv, fig25_dsmc_kn0p01.csv. + +import numpy as np + +import plume_figure_utils as u + + +def mass_flux_curve(theta_deg, S_0=2.0, r_over_D=10.0): + """(n/n0)*(Vr*sqrt(beta0)) along the angular curve (see header).""" + theta = np.deg2rad(theta_deg) + r = r_over_D * u.D_NOZZLE + X = r * np.cos(theta) + Z = r * np.sin(theta) + n = u.full_field_values('n', X, Z, S_0) + Vr = u.full_field_values('Vr', X, Z, S_0) + return n * Vr + + +def generate_figure(): + theta_deg = np.linspace(0.0, u.THETA_MAX_DEG, 90) + flux = mass_flux_curve(theta_deg) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + ax.plot(theta_deg, flux, linestyle='none', marker='o', mfc='none', + color='k', ms=4, markevery=2, label='Analytical') + u.overlay_digitized(ax, 'fig25', style='line') + + ax.set_xlim(0, 90) + ax.set_ylim(0, 0.03) + ax.set_xticks([0, 30, 60, 90]) + ax.set_xlabel(r'$\theta$') + ax.set_ylabel('Mass flux') + ax.set_title('Normalized mass flux along r/D = 10 (Fig. 25)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig25_mass_flux') + + +if __name__ == '__main__': + u.run_script(generate_figure)