Моделирование упругих столкновений

Абсолютно упругий удар — это модель соударения, при которой полная кинетическая энергия системы сохраняется. В классической механике при этом пренебрегают деформациями тел. Соответственно, считается, что энергия на деформации не теряется, а взаимодействие распространяется по всему телу мгновенно. Хорошим приближением к модели абсолютно упругого удара является столкновение бильярдных шаров или упругих мячиков.

Википедия

В статье мы будем рассматривать моделирование упругих столкновений на примере атомов одноатомного газа в двумерном пространстве.

Реализация

Я сначала расчет столкновений написал сам, потом спросил Claude, он мне показал более элегантную версию. Привожу его версию с небольшими корректировками.
Код на Javascript.

/**
 *  @param {Atom} atom1
 *  @param {Atom} atom2
 */
processCollision(atom1, atom2) {
  // vector from atom1 to atom2
  const dx = atom2.pos.x - atom1.pos.x
  const dy = atom2.pos.y - atom1.pos.y
  const distance = Math.sqrt(dx*dx + dy*dy)

  const minNoIntersectionDistance = atom1.radius + atom2.radius
  if (distance >= minNoIntersectionDistance) return

  // cos and sin of vector angle
  const cosA = dx / distance
  const sinA = dy / distance

  // relative velocity
  const dvx = atom1.velocity.x - atom2.velocity.x
  const dvy = atom1.velocity.y - atom2.velocity.y

  // relative velocity to normal
  const dvn = dvx * cosA + dvy * sinA

  // if atoms go to opposite directions, there is no collision
  // NOTE: for negative time the equality sign should be opposite
  if (dvn < 0) return

  // impulse of elastic collision
  const impulseCoef = (2 * dvn) / (atom1.mass + atom2.mass)

  // update velocities
  const dvx1 = -impulseCoef * atom2.mass * cosA
  const dvy1 = -impulseCoef * atom2.mass * sinA
  const dvx2 = impulseCoef * atom1.mass * cosA
  const dvy2 = impulseCoef * atom1.mass * sinA

  atom1.velocity.add(dvx1, dvy1)
  atom2.velocity.add(dvx2, dvy2)
}

Классы

class Point {
  constructor(x, y) {
    this.x = x
    this.y = y
  }

  add(x, y) {
    this.x += x
    this.y += y
  }
}

class Vector {
  constructor(x, y) {
    this.x = 0
    this.y = 0
    this.add(x, y)
  }

  add(x, y) {
    this.x += x
    this.y += y
    this.module = Math.sqrt(this.x * this.x + this.y * this.y)
  }

  cos() {
    return this.module < 0.0000001 ? 1 : this.x / this.module
  }

  sin() {
    return this.module < 0.0000001 ? 0 : this.y / this.module
  }
}

class Atom {
  /**
   * @param {Point} position
   * @param {number} radius
   * @param {number} mass
   * @param {Vector} velocity
   */
  constructor(position, radius, mass, velocity) {
    this.pos = position
    this.radius = radius
    this.mass = mass
    this.velocity = velocity
  }
}

Вызывается примерно так:

calcStep(dt) {
  for (let atom of this.atomList) {
    this.moveAtom(atom, dt)
  }

  for (let atom of this.atomList) {
    this.processCollisions(atom)
  }
}

/**
 * @param {Atom} atom
 * @param {number} dt
 */
moveAtom(atom, dt) {
  ...

  atom.pos.add(atom.velocity.x * dt, atom.velocity.y * dt)

  ...
}

processCollisions(atom) {
  const possibleCollisions = this.getPossibleCollisions(atom)

  for (let neighborAtom of possibleCollisions) {
    this.processCollision(atom, neighborAtom)
  }
}

Что здесь происходит:
Получаем вектор направления из первого объекта на второй, это ось столкновения.
Получаем его длину, проверяем, есть ли пересечение объектов.
Находим вектор скорости первого объекта относительно второго.
Проецируем его на ось столкновения.
Если его длина отрицательная, значит атомы разлетаются, столкновение рассчитывать не надо.
Мы нашли проекцию скорости вдоль линии столкновения, делаем расчет по формуле упругого столкновения для одной оси с учетом масс.

Здесь все компоненты обрабатываются в декартовых координатах, то есть раскладываются на компоненты в виде проекций на оси координат. Можно считать в полярных координатах “длина, направление”, но проще вычисления не будут. Там все равно нужно учитывать только проекцию скорости на ось столкновения, а проекцию в перпендикулярном направлении оставлять как есть, и потом по ним рассчитывать финальную длину и направление векторов.

Схема:

Примечания

Ось столкновения называется нормаль, потому что она идет перпендикулярно линии касания поверхностей.

Нет особых причин, почему относительная скорость вычисляется именно как скорость первого относительно второго, наоборот тоже можно, поменяются только знаки в нескольких местах. Но я бы сказал, на рисунке так получается нагляднее, когда относительная скорость указана для первого атома.

Если просто напрямую считать финальные скорости по формуле из Википедии, то может быть так, что атомы будут каждый раз менять направление скорости на противоположное, потому что они не успевают разлететься и все еще пересекаются, хотя уже разлетаются. Поэтому нужно специально проверять разлет.

Two-dimensional collision
Dot product

${\begin{aligned}\mathbf {v} '_{1}&=\mathbf {v} _{1}-{\frac {2m_{2}}{m_{1}+m_{2}}}\ {\frac {\langle \mathbf {v} _{1}-\mathbf {v} _{2},\,\mathbf {x} _{1}-\mathbf {x} _{2}\rangle }{\|\mathbf {x} _{1}-\mathbf {x} _{2}\|^{2}}}\ (\mathbf {x} _{1}-\mathbf {x} _{2}),\\\mathbf {v} '_{2}&=\mathbf {v} _{2}-{\frac {2m_{1}}{m_{1}+m_{2}}}\ {\frac {\langle \mathbf {v} _{2}-\mathbf {v} _{1},\,\mathbf {x} _{2}-\mathbf {x} _{1}\rangle }{\|\mathbf {x} _{2}-\mathbf {x} _{1}\|^{2}}}\ (\mathbf {x} _{2}-\mathbf {x} _{1})\end{aligned}}$ $\langle \mathbf {a}, \mathbf {b} \rangle =\left\|\mathbf {a} \right\|\left\|\mathbf {b} \right\|\cos \theta$

В коде можно использовать отрицательное время, все объекты будут правильно двигаться в обратном направлении, но в проверке на разлет if (dvn < 0) return нужно поменять знак сравнения на обратный (“больше нуля”), иначе группы атомов будут собираться вокруг некоторого общего центра.

В этом коде столкновения рассчитываются попарно сразу для обоих атомов. Это значит, что если атом сталкивается сразу с 2 другими, при расчете второго столкновения будет использована уже новая скорость после первого.
Чтобы это исправить, можно добавить в атомы специальное свойство dv, которое будет хранить изменение скорости, и переносить его в вектор скорости только после обработки всех столкновений.
Но это все равно не даст правильную обработку. Например, когда один атом влетает сразу в 3 других атома, каждое столкновение будет считаться независимо, и импульс от атома 1 учтется полностью 3 раза, а должен распределяться на 3 атома. Claude и Google AI мне подсказали, что там нужно решать сложную систему уравнений.

Здесь можно почитать подробнее про столкновения нескольких объектов.

Интересно, как столкновения работают во Вселенной на низком уровне? Думаю, попарно, так как это значительно проще, просто на макроскопическом уровне различия незаметны.

Если квант времени будет слишком большой, то объекты могут пролетать сквозь друг друга. Поэтому квант времени нужно делать поменьше.

Теоретически можно вместо этого проверять столкновения не только в финальной точке, но и во всех точках пространства до нее. Это означает введение кванта пространства.
У этого есть интересное следствие. Чем больше скорость объекта, тем больше квантов пространства он проходит за квант времени. Чем больше квантов пространства он проходит, тем больше нам нужно проверить столкновений, и тем меньше у нас остается части от этого кванта времени на обработку других взаимодействий. То есть все процессы объекта начинают идти медленнее, собственное время объекта замедляется. Как в теории относительности.

Результат

Посмотреть результат и полный код можно тут.

GIF 8Mb

Гравитация

Можно добавить гравитацию. Просто добавляем нужную величину к вертикальной компоненте скорости всех атомов - ускорение, умноженное на время кадра.

calcStep(dt) {
  // гравитация
  for (let atom of this.atomList) {
    atom.velocity.add(0, -this.gravity * dt)
  }

  for (let atom of this.atomList) {
    this.moveAtom(atom, dt)
  }

  for (let atom of this.atomList) {
    this.processCollisions(atom)
  }
}

GIF 1Mb

推荐订阅源

Все публикации подряд на Хабре