同一极坐标下不同粒径粒子轨迹仿真异常：路径呈圆形原因排查

问题描述 #

我尝试用Matplotlib通过以下Python代码在同一极坐标图中绘制不同粒径粒子的运动轨迹，预期轨迹和参考图一致，但实际运行后出现以下问题：

• 仅能看到部分粒径的轨迹，前四个粒径的轨迹因重叠不可见
• 所有轨迹都变成了圆形
• 单独运行不同粒径的仿真时，轨迹各不相同
• 调整粒径参数后，只有特定粒径的轨迹可见

仿真代码 #

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi

# Constants
μ = 379312077e8
R = 60268e3
g_10 = 21141e-9
Ω = 1.7e-4
ρ = 1e3
b_values = [0.95e-10, 0.94e-10, 0.26e-10, 0.25e-10, 0.23e-10, 0.1e-10]
V = 10
ε = 8.85e-12
j_2 = 1.629071e-2

# Function to define the ODEs
def odes(t, p,b):
    r, x, θ, y, ϕ, z = p

# Constants for the given b
    m = (4/3) * pi * (b**3) * ρ
    q = 4 * pi * ε * b * V
    β = q / m

# Other constants
    μ_0 = 4 * pi * 1e-7
    B_θ = (μ_0 / (4 * pi)) * (R / r)**3 * g_10 * sin(θ)
    B_r = (μ_0 / (4 * pi)) * 2 * (R / r)**3 * g_10 * cos(θ)

# Defining the ODEs
    drdt = x
    dxdt = r * (y**2 + (z + Ω)**2 * sin(θ)**2 - β * z * sin(θ) * B_θ) - (μ / r**2) * (1 - (3/2) * j_2 * (R / r)**2 * (3 * cos(θ)**2 - 1))
    dθdt = y
    dydt = (-2 * x * y + r * (z + Ω)**2 * sin(θ) * cos(θ) + (β * r * z * sin(θ) * B_r)) / r + (3 * μ / r**2) * j_2 * (R / r)**2 * sin(θ) * cos(θ)
    dϕdt = z
    dzdt = (-2 * (z + Ω) * (x * sin(θ) + r * y * cos(θ)) + (β * (x * B_θ - r * y * B_r)) / (r * sin(θ)))

    return np.array([drdt, dxdt, dθdt, dydt, dϕdt, dzdt])
 # Define the event function
def event_func(t, p,b):
    return p[0] - R

event_func.terminal = True
event_func.direction = -1

# Time window
t_span = (0, 1000)
t = np.linspace(t_span[0], t_span[1], 1001)

# Initial conditions
r0 = 1.12 * R
p0 = np.array([1.12 * R, 0.0, 90.0 * (pi / 180), 0.0, 0.0, 0.0203 * (pi / 180)])

# Create a figure with a polar projection
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})

# Define colors for plotting
colors = ['r', 'g', 'b', 'm', 'y', 'k']

# Iterate through different values of b
for b, color in zip(b_values, colors):
    
# Solve IVP for the current value of b
    sol = solve_ivp(odes, t_span, p0, args=(b,), t_eval=t,events=(event_func,), method="DOP853",
                dense_output=True, vectorized=False, atol=1e-4, rtol=1e-5)

    # Extract r and ϕ values
    r = sol.y[0, :]
    ϕ = sol.y[4, :]

    # Plot r against ϕ for the current value of b
    ax.plot(ϕ, r, label=f'b = {b:.2e}', color=color)

# Set radial limits
ax.set_rmax(1.5 * R)

# Set radial ticks
ax.set_rticks([R, 1.12 * R, 1.18 * R, 1.2 * R, 1.3 * R])

# Add a circle at the center
circle = plt.Circle((0, 0), R, transform=ax.transData._b, facecolor='yellow',
                fill=True, edgecolor='red', linewidth=0.5, zorder=2, alpha=0.5)
ax.add_artist(circle)

# Set plot title and legend
ax.set_title("Polar Plot of r vs Azimuthal Angle (ϕ)")
ax.legend()
plt.show()

问题根源与修复方案 #

1. 事件函数参数传递错误 #

你的event_func需要接收b参数，但在solve_ivp的events参数中仅传入了函数名，未绑定当前循环的b值，导致事件触发逻辑完全失效。粒子不会在到达r=R时停止积分，而是一直运行到t=1000，最终轨迹趋近于圆形稳态。

修复方式：
在循环中用lambda绑定b参数，确保事件函数能正确获取当前粒径：

sol = solve_ivp(odes, t_span, p0, args=(b,), t_eval=t,
                events=lambda t, p: event_func(t, p, b),
                method="DOP853", dense_output=True, vectorized=False,
                atol=1e-8, rtol=1e-10)

2. 积分精度不足 #

部分粒径对应的运动方程刚性较强，默认的atol=1e-4和rtol=1e-5精度无法捕捉轨迹的细微差异，导致不同粒径的轨迹被“抹平”成相似的圆形。

修复方式：
调小精度阈值，比如设置atol=1e-8、rtol=1e-10，确保积分结果能保留轨迹的独特特征。

3. 额外优化建议 #

• 给不同轨迹添加差异化线型（如linestyle=['-', '--', ':', '-.', '-', '--']），即使颜色相近也能区分
• 打印sol.t_events确认每个粒径的粒子都触发了终端事件，而非运行到最大时间
• 单独输出每个粒径的r和ϕ数组，验证积分结果是否符合预期

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