# This script computes 2D potential flow around an airfoil represented by a distribution of point sources.
# It expects three files in the working directory:
# - NACA0012_x.txt : x-coordinates of source points along the airfoil surface (1D)
# - NACA0012_y.txt : y-coordinates of source points along the airfoil surface (1D)
# - NACA0012_sigma.txt : source strengths for each point (same length as x/y)
#
# It builds a 51x51 grid over a rectangular domain, superposes the velocity
# induced by all sources with a uniform free stream in +x, then computes
# streamfunction, potential, and Cp. It reports the max Cp and its indices,
# and renders three figures as requested.
#
# If the files are not present, the script will print a clear message telling you
# what to upload and where. No plots will be created in that case.
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
# ---------------------------
# User-adjustable parameters
# ---------------------------
U_inf = 1.0 # Free-stream speed (non-dimensionalized)
# Domain for visualization (choose a sensible window around a unit chord airfoil)
# You can change these if your data requires a different view.
x_min, x_max = -1.0, 2.0
y_min, y_max = -1.0, 1.0
# Grid resolution
nx = ny = 51
# Filenames
fx = Path("NACA0012_x.txt")
fy = Path("NACA0012_y.txt")
fs = Path("NACA0012_sigma.txt")
def load_column_vector(path: Path) -> np.ndarray:
# Loads a text file as a 1D float array (whitespace- or newline-separated)
return np.loadtxt(path, dtype=float).reshape(-1)
# ---------------------------
# Load inputs (or exit early)
# ---------------------------
if not (fx.exists() and fy.exists() and fs.exists()):
missing = [p.name for p in [fx, fy, fs] if not p.exists()]
print("Missing required input file(s):", ", ".join(missing))
print("Please upload these files with exactly these names into this session and re-run this cell:")
print(" - NACA0012_x.txt")
print(" - NACA0012_y.txt")
print(" - NACA0012_sigma.txt")
else:
x_src = load_column_vector(fx)
y_src = load_column_vector(fy)
sigma = load_column_vector(fs)
if not (len(x_src) == len(y_src) == len(sigma)):
raise ValueError("Input files must have the same length: len(x) == len(y) == len(sigma).")
n_sources = len(x_src)
# ---------------------------
# Build computational grid
# ---------------------------
x = np.linspace(x_min, x_max, nx)
y = np.linspace(y_min, y_max, ny)
X, Y = np.meshgrid(x, y, indexing="ij") # X,Y shapes: (nx, ny)
# Initialize velocity, potential, and streamfunction with uniform flow
u = np.full_like(X, U_inf, dtype=float) # u_x component
v = np.zeros_like(X) # u_y component
phi = U_inf * X.copy() # velocity potential
psi = U_inf * Y.copy() # streamfunction
# ---------------------------
# Superpose contributions from point sources
# ---------------------------
# 2D point source at (x_i, y_i) with strength sigma_i (per unit depth):
# u = (sigma_i / (2π)) * (x - x_i) / r^2
# v = (sigma_i / (2π)) * (y - y_i) / r^2
# φ = (sigma_i / (2π)) * ln(r)
# ψ = (sigma_i / (2π)) * arctan2(y - y_i, x - x_i)
#
# We add a small epsilon to r^2 to avoid singularities exactly at source points.
twopi = 2.0 * np.pi
eps2 = 1e-12
for i in range(n_sources):
dx = X - x_src[i]
dy = Y - y_src[i]
r2 = dx*dx + dy*dy + eps2
inv_r2 = 1.0 / r2
contrib = sigma[i] / twopi
u += contrib * dx * inv_r2
v += contrib * dy * inv_r2
# Potential and streamfunction
phi += contrib * 0.5 * np.log(r2) # 0.5*ln(r^2) == ln(r)
psi += contrib * np.arctan2(dy, dx)
# ---------------------------
# Pressure coefficient Cp
# ---------------------------
V2 = u*u + v*v
Cp = 1.0 - V2 / (U_inf**2)
# Find maximum Cp and its indices (in the Cp array)
flat_idx = np.argmax(Cp)
i_max, j_max = np.unravel_index(flat_idx, Cp.shape) # i: x-index, j: y-index
Cp_max = Cp[i_max, j_max]
print(f"Max Cp: {Cp_max:.6f}")
print(f"Indices of max Cp (i, j) with indexing='ij': ({i_max}, {j_max})")
print(f"Coordinates of max Cp (x, y): ({X[i_max, j_max]:.6f}, {Y[i_max, j_max]:.6f})")
# ---------------------------
# Plot 1: Streamlines + airfoil profile
# ---------------------------
fig1 = plt.figure(figsize=(7, 5), dpi=120)
plt.streamplot(x, y, u.T, v.T, density=1.2, linewidth=0.8) # transpose fields to (ny, nx) for streamplot
plt.plot(x_src, y_src, linewidth=2.0) # airfoil surface points
plt.xlabel("x")
plt.ylabel("y")
plt.title("Streamlines and NACA0012 Profile")
plt.axis("equal")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.tight_layout()
plt.show()
# ---------------------------
# Plot 2: Velocity potential φ + airfoil profile
# ---------------------------
fig2 = plt.figure(figsize=(7, 5), dpi=120)
cf = plt.contourf(X, Y, phi, levels=40)
plt.plot(x_src, y_src, linewidth=2.0)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Velocity Potential φ and NACA0012 Profile")
plt.axis("equal")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.colorbar(cf, shrink=0.85, label="φ")
plt.tight_layout()
plt.show()
# ---------------------------
# Plot 3: Pressure coefficient Cp with marker at max
# ---------------------------
fig3 = plt.figure(figsize=(7, 5), dpi=120)
cf2 = plt.contourf(X, Y, Cp, levels=40)
plt.plot(x_src, y_src, linewidth=2.0)
plt.scatter([X[i_max, j_max]], [Y[i_max, j_max]], s=50, marker="o") # location of max Cp
plt.xlabel("x")
plt.ylabel("y")
plt.title("Pressure Coefficient Cp with Max Location")
plt.axis("equal")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.colorbar(cf2, shrink=0.85, label="Cp")
plt.tight_layout()
plt.show()
