ubuntu2404
=>PYTHONTEX#py#default#default#0#code#####68#
import numpy as np
import matplotlib.pyplot as plt
def naca_4digit(m, p, t, c=1.0, n=100):
"""
Generate NACA 4-digit airfoil coordinates
m: maximum camber (as fraction of chord)
p: position of maximum camber (as fraction of chord)
t: thickness (as fraction of chord)
c: chord length
n: number of points
"""
# Cosine spacing for better resolution near leading edge
beta = np.linspace(0, np.pi, n)
x = c * (1 - np.cos(beta)) / 2
# Thickness distribution
yt = 5 * t * c * (0.2969 * np.sqrt(x/c) -
0.1260 * (x/c) -
0.3516 * (x/c)**2 +
0.2843 * (x/c)**3 -
0.1015 * (x/c)**4)
# Camber line
yc = np.zeros_like(x)
dyc_dx = np.zeros_like(x)
if p > 0: # Cambered airfoil
# Forward of maximum camber
mask1 = x <= p * c
yc[mask1] = (m * c * x[mask1] / (p * c)**2 *
(2 * p * c - x[mask1]) / c**2)
dyc_dx[mask1] = 2 * m / p**2 * (p - x[mask1]/c)
# Aft of maximum camber
mask2 = x > p * c
yc[mask2] = (m * c * (c - x[mask2]) / (1 - p)**2 / c**2 *
(1 + x[mask2]/c - 2*p))
dyc_dx[mask2] = 2 * m / (1 - p)**2 * (p - x[mask2]/c)
# Surface coordinates
theta = np.arctan(dyc_dx)
xu = x - yt * np.sin(theta)
yu = yc + yt * np.cos(theta)
xl = x + yt * np.sin(theta)
yl = yc - yt * np.cos(theta)
return xu, yu, xl, yl, x, yc
# Generate NACA 2412 airfoil
xu, yu, xl, yl, x_camber, y_camber = naca_4digit(0.02, 0.4, 0.12)
# Plot airfoil
plt.figure(figsize=(10, 6))
plt.plot(xu, yu, 'b-', linewidth=2, label='Upper surface')
plt.plot(xl, yl, 'r-', linewidth=2, label='Lower surface')
plt.plot(x_camber, y_camber, 'k--', linewidth=1, label='Camber line')
plt.xlabel('x/c')
plt.ylabel('y/c')
plt.title('NACA 2412 Airfoil Geometry')
plt.grid(True, alpha=0.3)
plt.legend()
plt.axis('equal')
plt.tight_layout()
plt.savefig('naca2412_geometry.pdf', dpi=300, bbox_inches='tight')
plt.show()
print(f"NACA 2412 Airfoil Parameters:")
print(f"Maximum camber: 2% at 40% chord")
print(f"Maximum thickness: 12% chord")
print(f"Leading edge radius: {5 * 0.12**2:.4f}c")
=>PYTHONTEX#py#default#default#1#code#####151#
# Simplified thin airfoil theory calculations
def thin_airfoil_theory(alpha_deg, camber_max, camber_pos):
"""
Calculate lift coefficient using thin airfoil theory
alpha_deg: angle of attack in degrees
camber_max: maximum camber as fraction of chord
camber_pos: position of maximum camber as fraction of chord
"""
alpha_rad = np.radians(alpha_deg)
# Lift curve slope (per radian)
a0 = 2 * np.pi
# Zero-lift angle of attack for cambered airfoil
alpha_L0 = -2 * camber_max * (1 - 2 * camber_pos) / np.pi * np.pi/180
# Lift coefficient
cl = a0 * (alpha_rad - alpha_L0)
return cl, alpha_L0 * 180/np.pi
# Performance analysis for NACA 2412
alpha_range = np.linspace(-5, 15, 21)
cl_values = []
alpha_L0_theory = None
for alpha in alpha_range:
cl, alpha_L0 = thin_airfoil_theory(alpha, 0.02, 0.4)
cl_values.append(cl)
if alpha_L0_theory is None:
alpha_L0_theory = alpha_L0
cl_values = np.array(cl_values)
# Plot lift curve
plt.figure(figsize=(10, 6))
plt.plot(alpha_range, cl_values, 'bo-', linewidth=2, markersize=6,
label='Thin Airfoil Theory')
plt.axhline(y=0, color='k', linestyle='--', alpha=0.5)
plt.axvline(x=alpha_L0_theory, color='r', linestyle='--', alpha=0.7,
label=f'Zero-lift AoA = {alpha_L0_theory:.2f}°')
plt.xlabel('Angle of Attack (degrees)')
plt.ylabel('Lift Coefficient, $C_L$')
plt.title('NACA 2412 Lift Curve - Thin Airfoil Theory')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('naca2412_lift_curve.pdf', dpi=300, bbox_inches='tight')
plt.show()
print(f"Aerodynamic Analysis Results:")
print(f"Zero-lift angle of attack: {alpha_L0_theory:.2f} degrees")
print(f"Lift curve slope: {2*np.pi:.3f} per radian "
f"({2*np.pi*180/np.pi:.1f} per degree)")
print(f"At alpha = 5 degrees: CL = {cl_values[np.argmin(np.abs(alpha_range - 5))]:.3f}")
=>PYTHONTEX#py#default#default#2#code#####234#
import pandas as pd
# Create table of important dimensionless parameters
parameters = {
'Parameter': ['Reynolds Number', 'Mach Number', 'Prandtl Number',
'Lift Coefficient', 'Drag Coefficient', 'Pressure Coefficient'],
'Symbol': ['$Re$', '$Ma$', '$Pr$', '$C_L$', '$C_D$', '$C_p$'],
'Definition': ['$\\frac{\\rho U L}{\\mu}$', '$\\frac{U}{a}$', '$\\frac{\\mu c_p}{k}$',
'$\\frac{L}{\\frac{1}{2}\\rho U^2 S}$', '$\\frac{D}{\\frac{1}{2}\\rho U^2 S}$',
'$\\frac{p - p_\\infty}{\\frac{1}{2}\\rho U^2}$'],
'Physical Meaning': ['Inertial/Viscous forces', 'Flow/Sound speed', 'Momentum/Thermal diffusivity',
'Lift/Dynamic pressure', 'Drag/Dynamic pressure', 'Local/Dynamic pressure']
}
df = pd.DataFrame(parameters)
print("Key Dimensionless Parameters in Aerospace CFD:")
print("=" * 40)
for i, row in df.iterrows():
print(f"{row['Parameter']:<20} {row['Symbol']:<8} {row['Physical Meaning']}")
=>PYTHONTEX#py#default#default#3#code#####260#
def turbojet_performance(Ma0, gamma=1.4, cp=1005, T0=288.15, pi_c=8, pi_b=0.95,
T04=1400, eta_c=0.85, eta_t=0.88, eta_n=0.98):
"""
Simplified turbojet engine performance analysis
Ma0: flight Mach number
gamma: specific heat ratio
cp: specific heat at constant pressure (J/kg/K)
T0: ambient temperature (K)
pi_c: compressor pressure ratio
pi_b: burner pressure ratio
T04: turbine inlet temperature (K)
eta_c: compressor efficiency
eta_t: turbine efficiency
eta_n: nozzle efficiency
"""
R = cp * (gamma - 1) / gamma
# Station 0: Ambient conditions
T00 = T0 * (1 + (gamma - 1) / 2 * Ma0**2)
p00 = 101325 * (T00 / T0)**(gamma / (gamma - 1))
# Station 2: Compressor exit (ideal)
T02_ideal = T00 * pi_c**((gamma - 1) / gamma)
T02 = T00 + (T02_ideal - T00) / eta_c
p02 = p00 * pi_c
# Station 3: Burner exit
T03 = T04 # Assuming complete combustion
p03 = p02 * pi_b
# Station 4: Turbine exit
# Work balance: turbine work = compressor work
T04_ideal = T03 - (T02 - T00) / eta_t
T04_actual = T03 - (T02 - T00) / eta_t
pi_t = (T04_actual / T03)**(gamma / (gamma - 1))
p04 = p03 * pi_t
# Station 9: Nozzle exit
p9 = 101325 # Assume perfectly expanded
T09_ideal = T04_actual * (p9 / p04)**((gamma - 1) / gamma)
T9 = T04_actual - eta_n * (T04_actual - T09_ideal)
# Exhaust velocity
V9 = np.sqrt(2 * cp * (T04_actual - T9))
# Flight velocity
V0 = Ma0 * np.sqrt(gamma * R * T0)
# Specific thrust and specific fuel consumption (simplified)
F_specific = V9 - V0 # N per kg/s of air
return {
'V0': V0,
'V9': V9,
'F_specific': F_specific,
'T02': T02,
'T04': T04_actual,
'pi_c': pi_c,
'pi_t': pi_t
}
# Performance analysis for different flight conditions
mach_numbers = np.linspace(0, 2.0, 21)
results = []
for Ma in mach_numbers:
perf = turbojet_performance(Ma)
results.append([Ma, perf['F_specific'], perf['V9'], perf['V0']])
results = np.array(results)
# Plot performance
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.plot(results[:, 0], results[:, 1], 'b-o', linewidth=2, markersize=4)
plt.xlabel('Flight Mach Number')
plt.ylabel('Specific Thrust (N·s/kg)')
plt.title('Turbojet Specific Thrust vs Flight Mach')
plt.grid(True, alpha=0.3)
plt.subplot(2, 2, 2)
plt.plot(results[:, 0], results[:, 2], 'r-', linewidth=2, label='Exhaust Velocity')
plt.plot(results[:, 0], results[:, 3], 'g--', linewidth=2, label='Flight Velocity')
plt.xlabel('Flight Mach Number')
plt.ylabel('Velocity (m/s)')
plt.title('Velocities vs Flight Mach')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(2, 2, 3)
efficiency = results[:, 1] * results[:, 3] / (0.5 * (results[:, 2]**2 - results[:, 3]**2))
plt.plot(results[:, 0], efficiency, 'purple', linewidth=2)
plt.xlabel('Flight Mach Number')
plt.ylabel('Propulsive Efficiency')
plt.title('Propulsive Efficiency vs Flight Mach')
plt.grid(True, alpha=0.3)
plt.subplot(2, 2, 4)
thrust_ratio = results[:, 1] / results[0, 1] # Normalized to Ma=0
plt.plot(results[:, 0], thrust_ratio, 'orange', linewidth=2)
plt.xlabel('Flight Mach Number')
plt.ylabel('Normalized Specific Thrust')
plt.title('Thrust Variation with Altitude')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('turbojet_performance.pdf', dpi=300, bbox_inches='tight')
plt.show()
print("Turbojet Performance Summary:")
print(f"Sea level static specific thrust: {results[0, 1]:.1f} N·s/kg")
print(f"Cruise (Ma=0.8) specific thrust: {results[16, 1]:.1f} N·s/kg")
print(f"Maximum analyzed Mach: {results[-1, 0]:.1f}")
=>PYTHONTEX#py#default#default#4#code#####388#
def static_stability_analysis(cg_position, ac_position, cl_alpha, cm_alpha_wing,
cm_alpha_tail, tail_volume):
"""
Analyze longitudinal static stability
cg_position: center of gravity position (fraction of MAC)
ac_position: aerodynamic center position (fraction of MAC)
cl_alpha: lift curve slope (per radian)
cm_alpha_wing: wing pitching moment curve slope
cm_alpha_tail: tail pitching moment curve slope
tail_volume: tail volume coefficient
"""
# Static margin
static_margin = ac_position - cg_position
# Total pitching moment curve slope
cm_alpha_total = cm_alpha_wing + cm_alpha_tail * tail_volume
# Neutral point
neutral_point = ac_position - cm_alpha_total / cl_alpha
# Stability criterion
is_stable = static_margin > 0
return {
'static_margin': static_margin,
'neutral_point': neutral_point,
'cm_alpha_total': cm_alpha_total,
'is_stable': is_stable
}
# Example aircraft stability analysis
aircraft_data = {
'cg_range': np.linspace(0.20, 0.35, 100), # CG positions from 20% to 35% MAC
'ac_position': 0.25, # Aerodynamic center at 25% MAC
'cl_alpha': 2 * np.pi, # Lift curve slope
'cm_alpha_wing': -0.1, # Wing contribution
'cm_alpha_tail': -0.8, # Tail contribution
'tail_volume': 0.5 # Tail volume coefficient
}
stability_results = []
for cg in aircraft_data['cg_range']:
result = static_stability_analysis(
cg, aircraft_data['ac_position'], aircraft_data['cl_alpha'],
aircraft_data['cm_alpha_wing'], aircraft_data['cm_alpha_tail'],
aircraft_data['tail_volume']
)
stability_results.append([cg, result['static_margin'], result['is_stable']])
stability_results = np.array(stability_results)
# Plot stability analysis
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(stability_results[:, 0] * 100, stability_results[:, 1] * 100,
'b-', linewidth=2)
plt.axhline(y=0, color='r', linestyle='--', linewidth=2, label='Neutral Stability')
stable_mask = stability_results[:, 1] > 0
plt.fill_between(stability_results[stable_mask, 0] * 100,
stability_results[stable_mask, 1] * 100,
alpha=0.3, color='green', label='Stable Region')
plt.xlabel('CG Position (% MAC)')
plt.ylabel('Static Margin (% MAC)')
plt.title('Longitudinal Static Stability')
plt.grid(True, alpha=0.3)
plt.legend()
plt.subplot(1, 2, 2)
# Control authority analysis (simplified)
elevator_deflection = np.linspace(-20, 20, 41)
cm_delta_e = -0.02 # Elevator effectiveness (per degree)
trim_moments = []
for delta_e in elevator_deflection:
cm_trim = cm_delta_e * delta_e
trim_moments.append(cm_trim)
plt.plot(elevator_deflection, trim_moments, 'g-', linewidth=2,
label='Elevator Moment')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.5)
plt.xlabel('Elevator Deflection (degrees)')
plt.ylabel('Pitching Moment Coefficient')
plt.title('Control Authority')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('flight_dynamics_stability.pdf', dpi=300, bbox_inches='tight')
plt.show()
print("Flight Dynamics Analysis:")
print(f"Aerodynamic center: {aircraft_data['ac_position']*100:.1f}% MAC")
print(f"Stable CG range: {stability_results[stable_mask, 0].min()*100:.1f}% to {stability_results[stable_mask, 0].max()*100:.1f}% MAC")
print(f"Maximum static margin: {stability_results[:, 1].max()*100:.1f}% MAC")
=>PYTHONTEX#py#default#default#5#code#####504#
def heat_equation_1d(nx=50, nt=100, dx=0.02, dt=0.001, alpha=0.01):
"""
Solve 1D heat equation using finite differences
Demonstrates numerical methods used in CFD
"""
x = np.linspace(0, 1, nx)
T = np.zeros((nt, nx))
# Initial condition: Gaussian temperature distribution
T[0, :] = np.exp(-50 * (x - 0.5)**2)
# Boundary conditions (Dirichlet)
T[:, 0] = 0
T[:, -1] = 0
# Time stepping using explicit finite differences
r = alpha * dt / dx**2 # Stability parameter
for n in range(1, nt):
for i in range(1, nx-1):
T[n, i] = T[n-1, i] + r * (T[n-1, i+1] - 2*T[n-1, i] + T[n-1, i-1])
return x, T
# Solve heat equation
x, T_solution = heat_equation_1d()
# Plot temperature evolution
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
time_steps = [0, 20, 50, 99]
for t in time_steps:
plt.plot(x, T_solution[t, :], label=f't = {t*0.001:.3f}s')
plt.xlabel('Position (x)')
plt.ylabel('Temperature')
plt.title('1D Heat Equation Solution')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
X, T_time = np.meshgrid(x, np.linspace(0, 0.099, 100))
plt.contourf(X, T_time, T_solution, levels=20, cmap='hot')
plt.colorbar(label='Temperature')
plt.xlabel('Position (x)')
plt.ylabel('Time (s)')
plt.title('Temperature Evolution')
plt.tight_layout()
plt.savefig('heat_equation_solution.pdf', dpi=300, bbox_inches='tight')
plt.show()
print("Numerical Method Demonstration:")
print(f"Grid points: {len(x)}")
print(f"Time steps: {T_solution.shape[0]}")
print(f"Stability parameter r = {0.01 * 0.001 / 0.02**2:.3f} (should be <= 0.5)")
print(f"Maximum temperature at t=0: {T_solution[0, :].max():.3f}")
print(f"Maximum temperature at final time: {T_solution[-1, :].max():.3f}")
=>PYTHONTEX:SETTINGS#
version=0.18
outputdir=pythontex-files-main
workingdir=.
workingdirset=false
gobble=none
rerun=default
hashdependencies=default
makestderr=false
stderrfilename=full
keeptemps=none
pyfuture=default
pyconfuture=none
pygments=true
pygglobal=:GLOBAL||
fvextfile=-1
pyconbanner=none
pyconfilename=stdin
depythontex=false
pygfamily=py|python3|
pygfamily=pycon|pycon|
pygfamily=sympy|python3|
pygfamily=sympycon|pycon|
pygfamily=pylab|python3|
pygfamily=pylabcon|pycon|