В этой статье из цикла «Физика на пальцах» речь пойдёт об одном из подходов, с успехом применяемом в самых разных областях математического моделирования: начиная от динамики жидкости и заканчивая физикой взаимодействия твёрдых тел.
Рассмотрим обычный подход, ещё называемый velocity-based, к моделированию, например, движения материальной точки. Точку в таком подходе принято описывать её мгновенным положением, скоростью и ускорением:
И если мы хотим предсказать, как изменится положение точки, когда пройдёт интервал времени dt, мы считаем, что на протяжении этого интервала времени частица движется с постоянным ускорением acceleration. Так как dt мал, мы можем посчитать приближённое изменение скорости и положения, численно проинтегрировав по Эйлеру. Это нетрудно себе представить: скорость изменяется на малую величину acceleration * dt, а положение — на velocity * dt.
+ Для интересующихся мат.частью
− Скрыть
Пусть точка движется по некоторой траектории x(t). Разложим положение точки в момент времени t + dt в ряд Тейлора:
x(t + dt) = x(t) + x'(t) * dt + o(dt2)
x'(t + dt) = x'(t) + x''(t) * dt + o(dt2)
То есть ошибка одного шага интегрирования - o(dt2). Общее количество шагов интегрирования - порядка 1/dt, таким образом погрешность всего счёта - o(dt2) / dt = o(dt), то есть метод Эйлера обладает первым порядком точности.
К первому равенству можно прибавить x''(t) * dt2 / 2, тогда координата будет считаться со вторым порядком точности, но скорость по-прежнему с первым.
Отличие же Position-Based подхода от Velocity-Based состоит в том, что мы не храним скорость и ускорение в явном виде. Вместо этого мы храним положение точки на предыдущем кадре и на текущем. Например, так:
А передвигаем точку из интуитивно понятного принципа: если с предыдущего кадра до текущего прошёл малый интервал времени, и за это время частица сместилась на малое расстояние deltaPos, то, скорее всего, за следующий такой же интервал времени частица сместится на такое же расстояние. За объяснением, откуда взялся член с ускорением, придётся открыть спойлер:
+ Строгое обоснование
− Скрыть
Разложим положение точки x(t + dt) и x(t - dt) в текущий и предыдущий моменты времени в ряд:
x(t + dt) = x(t) + x'(t) * dt + x''(t) * dt2 / 2 + x'''(t) * dt3 / 6 + o(t4)
x(t - dt) = x(t) - x'(t) * dt + x''(t) * dt2 / 2 - x'''(t) * dt3 / 6 + o(t4)
Сложим два равенства:
x(t + dt) + x(t - dt) = 2 * x(t) + x''(t) * dt2 + o(t4)
Переносим в правую часть, получаем:
x(t + dt) = 2 * x(t) - x(t - dt) + x''(t) * dt2 + o(t4)
Видим, что метод имеет третий порядок точности, что на два порядка выше базовой версии интегратора Эйлера и на один порядок выше, чем если учитывать в интеграторе Эйлера ускорение.
Обратим внимание, что для успешной работы метода dt всегда должен быть постоянным, так как мы считаем, что от предыдущего шага интегрирования до текущего прошло столько же времени, сколько от текущего до следующего.
В чём же профит?
Что же нам даёт этот метод? Рассмотрим, например, взаимодействие большого количества круглых частиц друг с другом:
Сперва расширим класс нашей position-based точки, добавив ей радиус, теперь это будет полноценная частица:
В классическом честном velocity-based подходе нам бы пришлось считать время до ближайшего соударения двух шаров, смещать всю систему на это малое расстояние по времени, высчитывать их углы отскока, заниматься неинтересными и вычислительно дорогими вещами. В нашем же подходе мы делаем следующее: просто расталкиваем на равное расстояние шары, которые на текущей итерации пролезли немного друг в друга:
for(int i = 0; i < particles.size() - 1; i++)
{
for(int j = i + 1; j < particles.size(); j++)
{
Vector3 p1 = particles[i].position;
Vector3 p2 = particles[j].position;
if((p1 - p2).SquareLength()< sqr(particles[i].radius + particles[j].radius))// расстояние между центрами шаров меньше// суммы их радиусов -> они пересекаются
{
// шары удобнее всего расталкивать вдоль прямой, проходящей через их центры
Vector3 penetrationDirection = (p2 - p1).Normalize();
// глубина проникновения - это сумма радиусов двух шаров// минус расстояние между их центрамиfloat penetrationDepth = particles[i].radius + particles[j].radius
- (p2 - p1).Length();
// расталкиваем шары в противоположенных направлениях
particles[i].position += -penetrationDirection * penetrationDepth * 0.5f;
particles[j].position += penetrationDirection * penetrationDepth * 0.5f;
}
}
}
Данный цикл лучше прогнать несколько раз, так как при выталкивании шара A из шара B, он может ещё глубже залезть в шар C, но растолкнув их попарно несколько раз, в пределе все коллизии разрешатся.
Далее, если мы вызовем для всех шаров функцию ::Move, столкнувшиеся шары «магическим» образом разлетятся в нужном направлении с нужными скоростями!
for(int i = 0; i < particles.size(); i++)
{
paraticles[i].Move(dt);
}
+ Почему это работает
− Скрыть
Метод сохраняет порядок точности даже при "жёстком" взаимодействии, например, при столкновении шаров. Смещение позиции можно интерпретировать как некоторое ускорение x''(t), такое, что x''(t) * dt2 нейтрализует коллизию. То есть нам не важно, какое конкретно ускорение получают шарики при столкновении, нам важно, что это ускорение устраняет их коллизию за один шаг интегрирования.
Несколько мутный аспект такого расчёта - скорость шаров при отскоке меняется непредсказуемо. Энергия гарантированно не повысится, то есть шары разлетятся со скоростью меньшей, чем слетались, но конкретная величиная отскока не определена. Энтузиасты могут попытаться избавиться от этого ограничения, модифицировав prevPos в соответствии с законом сохранения импульса, но доказательство корректности такого подхода выходит за рамки статьи.
С виду простой метод возволяет моделировать вполне внушительное количество взаимодействующих объектов:
На этом просчитывается взаимодействие 5000 шариков, расчёт всей физики занимает чуть больше 2мс на одном ядре процессора Core2Duo E6600 2.4GHz.
Но и это ещё не всё. Предположим, мы хотим соединить какие-то пары шариков жёсткими связями. Заведём тип данных, описывающий эту самую связь:
После каждого шага по времени мы варварски передвигаем частицы на расстояние, соответствующее их смещению с предыдушего кадра к текущему. Разумеется, из-за этого связи могут нарушаться, растягиваясь или сжимаясь. Придерживаясь нашего position-based подхода, всё, что от нас требуется — просто сжать или растянуть связь до её собственной длины. Реализуем это в методе Solve():
Обратим внимание на параметр stiffness. Если он равен нулю, то после «решения» связи частицы не сдвинутся с места. Если он равен единице, частицы примут своё положение, как если бы связь мгновенно приняла «удобную» для себя длину. Этот параметр означает что-то вроде жёсткости связи. Его максимальное значение — единица, соостветствует максимальной жёсткости, допустимой численным методом.
Грубо говоря, всё, что от нас требуется для реализации какого-то физического явления — это каждый кадр ставить все частицы в положение, удовлетворяющее этому явлению и более-менее близкое к текущему. Если это затруднительно сделать, как, например, в случае с большим количеством связей не тривиально найти положение шариков, когда не нарушена ни одна связь, мы просто выбираем «менее плохую» конфигурацию, постепенно решая, каждую связь отдельно от других. Всё остальное получается как бы само собой — шарики расталкиваются, конструкция из связей получает жёсткость.
Further Improvements
Одной из неприятных особенностей метода является то, что энергия системы при каждом расталкивании пары шаров гарантированно не увеличивается. Более того, нельзя точно сказать, уменьшается ли она. В описанной реализации отскочут ли шары друг от друга, как мячи, или слипнутся, как пластилиновые, сказать заранее нельзя. Благо, в хаосе из тысяч частиц это не заметно, но если мы реализуем, например, бильярд, на самотёк это дело пускать нельзя.
Нужно найти способ, как явно задавать скорость частиц после соударения. Один из подходов — менять prevPosition, но я считаю его антилогичным. Вместо этого я предлагаю поступить следующим образом:
Вместо того, чтобы хранить пару prevPosition + position, можно хранить пару position + deltaPosition, где deltaPosition — это разница между position и prevPosition. Соотвественно, каждый раз изменяя значение position, нам нужно в точности на такую же величину изменить и deltaPosition, в этом нетрудно убедиться взяв листок бумаги, нарисовав предыдущую и текущую позицию точки и посмотрев, как изменяется вектор, их соединяющий, если двигать текущую позицию.
Теперь, чтобы корректно обработать отскок пары шаров, мы не только расталкиваем их (обратим внимание, теперь уже не простым присвоением текущей позиции, а методом Push()), но и пересчитываем новые скорости из закона сохранения импульса:
for(int i = 0; i < particles.size() - 1; i++)
{
for(int j = i + 1; j < particles.size() - 1; j++)
{
Vector3 p1 = particles[i].position;
Vector3 p2 = particles[j].position;
if((p1 - p2).SquareLength()< sqr(particles[i].radius + particles[j].radius))// расстояние между центрами шаров// меньше суммы их радиусов -> они пересекаются
{
// шары удобнее всего расталкивать вдоль прямой,// проходящей через их центры
Vector3 penetrationDirection = (p2 - p1).Normalize();
// глубина проникновения - это сумма радиусов двух шаров// минус расстояние между их центрамиfloat penetrationDepth = particles[i].radius + particles[j].radius
- (p2 - p1).Length();
//расталкиваем шары в противоположенных направлениях
particles[i].Push(-penetrationDirection * penetrationDepth * 0.5f);
particles[j].Push( penetrationDirection * penetrationDepth * 0.5f);
float bounce = 0.0f; //коэффициент отскока, [0..1]//относительная скорость шариков
Vector3 relativeVelocity = particles[i].deltaPosition - particles[j].deltaPosition;
//тела обмениваются импульсом вдоль нормали контактаfloat exchangeVelocity = (1.0f + bounce) * (relativeVelocity * penetrationDirection);
if(exchangeVelocity >0)
{
particles[i].deltaPosition += penetrationDirection * exchangeVelocity * 0.5f;
particles[j].deltaPosition -= penetrationDirection * exchangeVelocity * 0.5f;
}
}
}
}
Для интересующихся, откуда был взят код пересчёта скоростей после отскока рекомендую обратиться к статье «Физика «на пальцах»: солверы физических движков» по разрешению контактов в импульсных солверах.
Ничто принципиально не мешает добавить в метод угловые скорости, трение и тела более сложной формы. Но, к сожалению, после таких улучшений подход выглядит уже далеко не так лаконично и красиво.
Немного модифицировав этот алгоритм, чтобы связи образовывались и рвались на определённом расстоянии между частицами, можно добиться интересных эффектов, например, моделировать что-то вроде вязкой жидкости:
Если добавить третье измерение, построить сетку из частиц и выправлять слишком «сложившиеся» её ячейки, получим что-то вроде ткани:
Если мы выделим группу частиц и скажем, что объём, который они ограничивают, должен сохраняться (этакий бурдюк с несжимаемой жидкостью), и будем каждый кадр просто увеличивать его (растягивая вдоль координатных осей) во столько раз, чтобы сохранить объём, получим так называемый baloon:
Наконец, если свяжем частицы не отрезками, а тетраэдрами, мы получим объёмные связи. «Выправляя» их от кадра к кадру можно получить нетривиальные явления, например, таким подходом удобно моделировать деформируемые тела:
или распространение ударного волнового фронта:
Приятная находка
Я ещё давно заметил, что если интегрировать движение тела таким образом:
Объяснение просто - "волшебная" формула при постоянном шаге является не интегратором по Эйлеру первого порядка, а переходит в position-based третьего порядка, в формуле достаточно сделать замену velocity = deltaPosition / dt и она полностью перейдёт в ту, что описана в Further Improvements.
Ограничения метода
Разумеется, чудес не бывает и не бывает одновременно простого, быстрого и общего метода моделирования физики. За простоту реализации и широкий спектр моделируемых явлений position-based подход заплатил своё:
Любые попытки изменять шаг по времени до добра не доведут. Если оригинальная схема имеет третий порядок точности по времени, то использование нефиксированного шага сбрасывает порядок точности до первого. Если в игре не фиксированный fps, лучше вызывать физику не каждый кадр а, например, тогда, когда с момента предыдущего обновления физики прошло какое-то время, например, 1/50 секунды.
Метод вовсе не дружит с continuous collision detection. Если неаккуратно выбрать размеры или скорости тел, они могут безнаказанно пролетать друг через друга. Если избавиться от быстро летящих частиц нельзя, то можно попытаться обрабатывать их другим, не position-based подходом.