欧拉法求解滑翔机耦合微分方程：攻角发散问题求助

问题背景 #

我3D打印并测试了一款滑翔机，尝试用欧拉法模拟其无推力飞行过程，模型考虑重力、升力与阻力作用，推导得到耦合微分方程：

dX/dt = -B(X²+Y²)(Clsin(θ)+Cdcos(θ))
dY/dt = -g + B(X²+Y²)(Clcos(θ)-Cdsin(θ))
其中 X=Vx，Y=Vy，B=0.5ρS/m，Cl、Cd为升力系数和阻力系数

我从CSV文件导入NACA0015的攻角、Cl、Cd数据并完成插值，编写Python代码实现欧拉法求解，但无论调整初始条件、时间步长等参数，攻角都会快速发散，与实际测试结果不符。

完整代码 #

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()

请求解决思路 #

恳请各位提供解决思路，谢谢！

内容的提问来源于stack exchange，提问作者DarrOw