Моя модель следующая:
Планер запускается с начальной скоростью и углом атаки (θ в код), действующие силы здесь — это его вес, сопротивление и подъемная сила, отсутствие тяги вообще.
Я уравнял модель и вот что я нашел:
Код: Выделить всё
dX/dt = -B(X^2+Y^2)(Cl*sin(θ)+Cd*cos(θ) dY/dt = - g + B(X^2+Y^2)(Cl*cos(θ)-Cd*sin(θ) where X = Vx and Y = Vy, B = 0.5*µ*S/m and Cl, Cd the lift and drag coefficientsСначала я импортирую данные Cl/Cd из файла CSV и интерполирую их.
Я определяю свои функции dX/dt, dY/dt
Я определяю свою функцию Эйлера для решения
Код: Выделить всё
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as interpolate
# Data
rho = 1.2
m = 0.06
S = 0.05
g = 9.81
B = 0.5*rho*S*(1/m)
# Data on angle of attack, lift coefficient, and drag coefficient
angles = []
CL_raw = []
CD_raw = []
file = open("C:/Cours Prépa/TIPE/xf-naca0015-il-50000.csv")
L = file.readlines()
for i in range(11, 100):
angles.append(float(L[i].split(",")[0]))
CL_raw.append(float(L[i].split(",")[1]))
CD_raw.append(float(L[i].split(",")[2]))
# Interpolation of the data
f_CD = interpolate.interp1d(angles, CD_raw)
f_CL = interpolate.interp1d(angles, CL_raw)
# Creating arrays with higher resolution values
angles_interp = np.linspace(min(angles), max(angles), 1000)
CD_interp = f_CD(angles_interp) # Using interpolation function for CD
CL_interp = f_CL(angles_interp) # Using interpolation function for CL
# Functions to calculate derivatives dX/dt and dY/dt based on the inclination
def dX_dt(X, Y, θ):
CD = f_CD(θ)
CL = f_CL(θ)
return -B*(X**2 + Y**2)*(CL*np.sin(np.radians(θ)) + CD*np.cos(np.radians(θ)))
def dY_dt(X, Y, θ):
CD = f_CD(θ)
CL = f_CL(θ)
return -g + B*(X**2 + Y**2)*(CL*np.cos(np.radians(θ)) - CD*np.sin(np.radians(θ)))
# Main function for Euler's method
def euler(dt):
# Initial conditions
X = [15] # Initial horizontal velocity in m/s
Y = [0] # Initial vertical velocity in m/s
x = [0] # Initial horizontal position in m
y = [10] # Initial vertical position in m
θ = [np.rad2deg(np.arctan(Y[0]/X[0]))] # Initial angle of attack
v = [np.sqrt(X[0]**2 + Y[0]**2)]
elapsed_time = 0
i = 0
# Euler's method
while y[i] >= 0:
X.append(X[i] + dX_dt(X[i], Y[i], θ[i])*dt)
Y.append(Y[i] + dY_dt(X[i], Y[i], θ[i])*dt)
x.append(x[i] + X[i]*dt)
y.append(y[i] + Y[i]*dt)
v.append(np.sqrt(X[i]**2 + Y[i]**2))
elapsed_time += dt
θ.append(np.rad2deg(np.arctan(Y[i]/X[i])))
if θ[i] > angle_limit or θ[i] < -angle_limit:
print("Stall")
return x, y, v, elapsed_time, θ
i += 1
return x, y, v, elapsed_time, θ
# Simulation parameters
dt = 0.001 # Time step in seconds
angle_limit = 10 # Stall angle limit in degrees
x, y, v, elapsed_time, θ = euler(dt)
print(elapsed_time)
time_array = np.arange(0, elapsed_time, dt)
# Plotting the graph
fig = plt.figure(figsize=(15, 10))
plt.subplot(3, 2, 1)
plt.plot(angles, CD_raw, 'o', label='Raw data')
plt.plot(angles_interp, CD_interp, '-', label='Interpolation')
plt.xlabel('Angle of Attack (degrees)')
plt.ylabel('Drag Coefficient (CD)')
plt.title('Variation of CD with Angle of Attack')
plt.legend()
plt.grid(True)
plt.subplot(3, 2, 2)
plt.plot(angles, CL_raw, 'o', label='Raw data')
plt.plot(angles_interp, CL_interp, '-', label='Interpolation')
plt.xlabel('Angle of Attack (degrees)')
plt.ylabel('Lift Coefficient (CL)')
plt.title('Variation of CL with Angle of Attack')
plt.legend()
plt.grid(True)
plt.subplot(3, 2, 3)
plt.plot(x, y)
plt.ylim([0, 11])
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.title('Aircraft Trajectory')
plt.grid(True)
# Plotting the graph of the aircraft's inclination relative to the ground
plt.subplot(3, 2, 4)
plt.plot(time_array, θ)
plt.xlabel('Time')
plt.ylabel('Inclination (degrees)')
plt.title('Aircraft Inclination Relative to Ground')
plt.grid(True)
plt.subplot(3, 2, 5)
plt.plot(time_array, v)
plt.xlabel('Time')
plt.ylabel('Speed magnitude v')
plt.title('Speed Magnitude Over Time')
plt.grid(True)
plt.tight_layout()
plt.show()
Я перепробовал множество начальных условий, временных шагов и т. д., но неважно что я делаю, угол атаки всегда быстро меняется, чего не происходит при реальном тестировании.
Дайте мне знать, если у вас есть идеи
Спасибо
Когда я запускаю код, угол атаки внезапно падает, что бы я ни пытался.
Я перепробовал множество начальных условий, временных шагов и т. д., но независимо от того, что я делаю, угол атаки всегда быстро расходится, чего не происходит при реальном тестировании.
Дайте мне знать, если у вас есть идеи
Спасибо
Подробнее здесь: https://stackoverflow.com/questions/783 ... ers-method
