球状星团的NBody模拟

Question

我正试图用蛙跳方案编写一个n体模拟程序来模拟一个球状星团，但是我遇到了这样一个问题：粒子从（我认为）接近其他粒子的系统中被甩出，这导致了势能的大幅增加。我试着写一些代码来解释碰撞，但是如果我也使用了软化因子，我不确定应该使用什么样的碰撞半径

初始位置在球体内随机生成，所有初始速度均为0

如果将软化因子设置为生成球体的0.1*半径，并将碰撞半径设置为太阳半径（不会导致碰撞），则会得到以下结果：

从能量图中，我们可以看到，当物体最初彼此靠近时，势能会出现一个巨大的峰值。随着时间的推移，我的能量逐渐下降到0，这仅仅是由于数值计算

有没有办法阻止这些星星被抛出去。我试图研究初始团簇形成平衡态需要多长时间

球体生成：

sun_mass = 1.989e30


N = 100
mass = sun_mass
M = np.full([N],mass)
R = 1e13
epsilon = 0.1 * R
collision_radius = 7e8
V = np.zeros([N,3])
M = np.full([N],mass)
P = np.zeros([N,3])


t = 50 * 365 * 24 * 60 * 60
dt = 30 * 24 * 60 * 60



print(t/dt)

np.random.seed(54321)

for i in range(N):
    
    phi = np.random.uniform(0,2*np.pi)
    costheta = np.random.uniform(-1,1)
    u = np.random.uniform(0,1)
    
    theta = np.arccos( costheta )
    r = R * (u) **(1/3)
    
    
    x = r * np.sin( theta) * np.cos( phi )
    y = r * np.sin( theta) * np.sin( phi )
    z = r * np.cos( theta )

    P[i] = (x,y,z)

节目：

def programe(position, mass, velocity, softening, time, dt, R, collision_radius):
    
    no_of_time_steps = (round(time/dt))

    all_positions = []
    all_velocities = []
    #print(all_positions)
    #print(len(all_positions[0]))
    kinetic_energy = []
    potential_energy = []
    total_energy = []
    
        
    for i in range(no_of_time_steps):
        
        position, mass, velocity = detect_collisions(position, mass, velocity, collision_radius)
        
        all_positions.append(position)
        all_velocities.append(velocity)
    
        'graph'
        plots = np.round(np.linspace(0,no_of_time_steps,num=500))
        for k in range(len(plots)):
            if i == plots[k]:
                print("test")
                #print(i)
                graph(position, R, k)
        
        'energies'
        kinetic_energy.append(calc_kinetic_energy(velocity, mass))
        potential_energy.append(calc_potential_energy(position, mass))
        total_energy.append(calc_total_energy(position, velocity, mass))
        
        'leap frog'
        velocity = calc_next_v_half(position, mass, velocity, softening, dt)    
        position = calc_next_position(position, mass, velocity, dt)    
        velocity = calc_next_v_half(position, mass, velocity, softening, dt)
        

    
    all_positions = np.array(all_positions)
    graphing(all_positions, time, dt, kinetic_energy, potential_energy, total_energy, no_of_time_steps, R)
    #print(len(mass))

    return(all_positions, all_velocities, kinetic_energy, potential_energy, total_energy)

碰撞检测：

def indexOf(Array, item):
    for i in range(len(Array)):
        if (Array[i] == item).all():
            return i
        
def detect_collisions(position, mass, velocity, collision_radius):
    i = 0
    
    newP = position
    newM = mass
    newV = velocity
    #print(len(position), len(newM))
    while i < len(newM):
        if newM[i] == 0:
            i += 1
            continue
        j = i + 1
        while j < len(newP):
           #print(j)
            if newM[j] == 0:
                j += 1
                continue
            p1 = position[i]
            p2 = position[j]
            
            if calc_seperation(p1,p2) < collision_radius:
                index1 = indexOf(position, p1)
                index2 = indexOf(position, p2)
                print('collision', index1, index2)
                newM, newV, newP = handle_collision(newM, newV, newP, [index1, index2])
            j += 1
        i += 1
    
    return(newP, newM, newV)

def handle_collision(M, V, P, indexes):
    if M[indexes[0]] > M[indexes[1]]:
        primary = indexes[0]
        secondary = indexes[1]
    else:
        primary = indexes[1]
        secondary = indexes[0]

    primaryMomentum = M[primary] * V[primary]
    secondaryMomentum = M[secondary] * V[secondary]
    newMomentum = primaryMomentum + secondaryMomentum
    
    newMass = M[primary] + M[secondary]
    newVelocity = newMomentum / newMass
    M[primary] = newMass
    V[primary] = newVelocity
    M[secondary] = 0
    V[secondary] = 0
    P[secondary] = 0

    newM = np.delete(M, secondary, axis=0)
    newV = np.delete(V, secondary, axis=0)
    newP = np.delete(P, secondary, axis=0)

    return (newM, newV, newP)