Физика «на пальцах»: солверы физических движков

Автор: Александр Санников

Данная статья объясняет принципы работы солверов современных физических движков.

В этом очерке я постараюсь объяснить, как работают солверы современных на момент написания статьи движков игровой физики: PhysX, ODE, Bullet, Box2D. За строгой математикой любопытный читатель всегда может обратиться к многочисленным существующим пейперам, ссылки на которые можно найти в нашем сообществе и лекции, прочитанные в нём же. Я же постараюсь абстрагироваться от строгой математики и попытаюсь объяснить, как это работает, что же означает каждое действие солвера. В любом случае, понимая, как он работает, проще вести и диагностику и расширение функциональности, например, разработать новый, специфический тип джойнта. Также рекомендую к прочтению более обзорную статью о разработке своего импульсного физического движка.

Формулировка задачи
Описание метода Sequential Impulses
Warmstarting (горячий старт)
Случай упругого соударения
Усовершенствования алгоритма
Pseudovelocities (псевдоскорости)
Другие типы соединений (джойнтов)
Трение
Заключение
Код

Формулировка задачи

Прежде всего сформулируем в чём же состоит задача солвера. На входе у него есть особым образом заданные ограничения на движения тел, его задача состоит в том, чтобы воздействовать на тела таким образом, чтобы все эти ограничения не нарушались или нарушались как можно меньше.

К сожалению, солвер «понимает» далеко не любые типы ограничений, а только заданные особым образом — ограничивающие ровно одну степень свободы системы.

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

Так какие же степени свободы мы хотим ограничить, например, при столкновении двух тел? Рассмотрим простой случай, когда тела касаются в единственной точке point, тела взаимодействуют по нормали norm и «отскок» нулевой. Контакт будет считаться удовлетворённым, если проекция относительной скорости двух тел на нормаль контакта будет нулевой. Интуитивное объяснение - при столкновении тел их скорость меняется таким образом, что контактирующие точки перестают сближаться.

Рассмотрим на простом примере двух столкнувшихся шаров без учёта вращения:

скорость первого шара: body1->velocity,
скорость второго шара: body2->velocity,
относительная скорость в точке contact->point можно рассмотреть, например, в системе отсчёта второго шара, она будет равна:

  //все операции гораздо удобнее проводить над векторами,
  //а не над их отдельными компонентами
  Vector3 relativeVelocity = body1->velocity - body2->velocity;

а теперь просто выразим проекцию относительной скорости на нормаль контакта:

  //* - скалярное произведение
  float velocityProjection = relativeVelocity * contact->norm;

Теперь, если мы потребуем от солвера, чтобы velocityProjection было равно нулю, мы тем самым ограничим одну степень свободы системы — движение тел вдоль нормали контакта. Важно понять, что, удовлетворив этому требованию, мы лишь запретим контактирующим телам пролезать всё глубже и глубже друг в друга, никак не разрешая коллизию в случае, если они пересеклись. О том, как разрешать именно взаимопроникновения, мы рассмотрим в разделе "Pseudovelocities", ниже.

Как же учесть угловую составляющую скорости? Из курса механики известно, что если тело вращается вокруг центра масс body->pos с угловой скоростью body->angularVelocity, то скорость его точки point можно выразить как

  //^ - cross product, векторное произведение
  Vector3 pointVelocity = 
      body->velocity + (point - body->pos) ^ body->angularVelocity;

Таким образом распишем выражение для нашей ограниченной степени свободы с учётом угловой скорости:

  Vector3 relativeVelocity = 
      body1->velocity + (contact->point - body1->pos) ^ body1->angularVelocity
    - body2->velocity - (contact->point - body2->pos) ^ body2->angularVelocity
  float velocityProjection = relativeVelocity * contact->norm;
  //полагаем velocityProjection = 0

Преобразуем, объединив в одно выражение:

  float velocityProjection = 
      body1->velocity * contact->norm 
    + ((contact->point - body1->pos) ^ body1->angularVelocity) * contact->norm
    - body2->velocity * contact->norm 
    - ((contact->point - body2->pos) ^ body2->angularVelocity) * contact->norm;
  //полагаем velocityProjection = 0

Далее применяем нехитрое векторное преобразование (a ^ b) * c == (c ^ a) * b и ещё раз преобразуем формулу:

  float velocityProjection = 
      contact->norm * body1->velocity 
    + (contact->norm ^ (contact->point - body1->pos)) * body1->angularVelocity
    - contact->norm * body2->velocity
    - (contact->norm ^ (contact->point - body2->pos)) * body2->angularVelocity;
    //полагаем velocityProjection = 0

и замечаем, что эту формулу можно переписать в виде:

  Vector3 n1 = contact->norm;
  Vector3 w1 =  (contact->norm ^ (contact->point - body1->pos));
  Vector3 n2 = -contact->norm;
  Vector3 w2 = -(contact->norm ^ (contact->point - body2->pos));
  float velocityProjection =
      n1 * body1->velocity + w1 * body1->angularVelocity 
    + n2 * body2->velocity + w2 * body2->angularVelocity;

То есть получается, что ограничение, вносимое в систему из двух тел контактом, можно описать четырьмя векторами: n1, w1, n2, w2. справедливости ради заметим, что записанные подряд координаты этих самых векторов образуют строку Той Самой Матрицы Якоби, или Якобиана

j = (n1, w1, n2, w2),

а ограничение можно описать как j * v = 0, где v - это вектор обобщённых скоростей:

v = (body1->velocity, body1->angularVelocity, body2->velocity, body2->angularVelocity)

Описание метода Sequential Impulses

Осталось всего-то ничего, понять, что нужно сделать с телами, чтобы этому самому velocityProjection = 0 удовлетворить. Другими словами - что должно произойти с контактирующими телами, чтобы они перестали друг в друга проникать? Продолжаем рассуждать из физических соображений:

Взаимодействие столкнувшихся тел осуществляется передачей импульса в направлении contact->norm некоторой амплитуды lambda к точке contact->point к первому телу и такого же импульса противоположенного направления второму телу. Другими словами, тела при столкновении обмениваются имульсом contact->norm * lambda. Работу солвера в случае одного контакта можно свести к нахождению этой самой лямбды.

Ещё немного книжной механики. Распишем, как изменится линейная и угловая скорость тела, если приложить к его точке point импульс в направлении norm и с модулем lambda:

линейная составляющая высчитывается просто: делением вектора импульса на массу тела. Масса каждого тела для удобства хранится инвертированная, body->invMass = 1.0f / body->mass;

  Vector3 velocityAfterCollision =
      body->velocity + norm * lambda * body->invMass;

С угловой составляющей чуть сложнее. Для простоты положим, что эллипсоид инерции симметричен, выражен в шар, то есть задаётся единственным числом — моментом инерции вокруг произвольной оси, проходящей через центр масс тела. Нехилое допущение, но желающие без труда могут преобразовать формулу для произвольного тензора инерции.
Говоря простым языком, момент инерции в нашем случае - число, полностью аналогичное массе, только масса связывает силу с ускорением, а момент инерции связывает момент силы с угловым ускорением. Для тестовых расчётов момент инерции вполне можно брать порядка массы объекта, помноженной на квадрат его характерного размера(т.е. длины или, например, высоты). Если ошибиться в расчёте момента инерции раз в 10, то физика работать, просто несколько менее естественно, а неподготовленный глаз может и не заметить отличий. Момент инерции, как и в случае с массой, храним инвертированным: body->invInertia = 1.0f / body->inertia;

  Vector3 angularVelocityAfterCollision = body->angularVelocity 
        + ((norm * lambda) ^ (point - body->pos)) * body->invInertia;

Почему мы храним инвертированную массу и инвертированный момент инерции? На самом деле просто для удобства задания бесконечных масс, если body->mass устремить в бесконечность, то body->invMass будет равна нулю. Поэтому вместо обработки бесконечных величин мы просто храним обратные, которые в этом случае равны нулю.

Что же замечаем? Изменения линейной и угловой скоростей можно выразить через те самые n1, w1, n2, w2:

  float velocityProjection =
      n1 * (body1->velocity + n1 * body1->invMass * lambda) 
    + w1 * (body1->angularVelocity + w1 * body1->InvInertia * lambda)
    + n2 * (body2->velocity + n2 * body2->invMass * lambda) 
    + w2 * (body2->angularVelocity + w2 * body2->InvInertia * lambda);

Преобразуем, получаем:

float velocityProjection = n1 * body1->velocity + n2 * body2->velocity 
        + w1 * body1->angularVelocity + w2 * body2->angularVelocity
        + (n1 * n1 * body1->invMass  + w1 * w1 * body1->InvInertia 
        + n2 * n2 * body2->invMass  + w2 * w2 * body2->InvInertia) * lambda;

Напомню, что контакт считается решённым, если velocityProjection обращается в ноль. Получаем легко решающееся уравнение:

  float a = n1 * body1->velocity + n2 * body2->velocity
    + w1 * body1->angularVelocity
    + w2 * body2->angularVelocity;

  float b = (n1 * n1 * body1->invMass 
    + w1 * w1 * body1->InvInertia 
    + n2 * n2 * body2->invMass  
    + w2 * w2 * body2->InvInertia);

  a + b * lambda = 0;

Уравнение такое можно решить, очевидно: lambda = –a / b; таким образом мы нашли имульс, которым обмениваются тела при введении ограничения степени свободы (n1, w1, n2, w2). Важно отметить, что сила реакции опоры может быть направлена только от поверхности контакта, то есть никак не вовнутрь. Таким образом, если lambda получается отрицательной, это значит, что тела не расталкиваются, а слипаются. Этого допустить никак нельзя, и если в каком-то контактном джойнте(и вообще любом, где телам разрешено передавать импульс только в одном направлении), импульс(lambda) получается отрицательным, мы его просто обнуляем. Далее, применив этот импульс к обоим телам, находим корректированные скорости после контакта.

Если перед решением джойнта мы замечаем, что проекция относительной скорости тел в точке контакта на нормаль этого контакта уже положительная(физический смысл - тела разлетаются), то контакт заведомо можно пропустить, так как он уже считается решённым.

Вот, собственно, и все чудеса, такая нехитрая последовательность действий позволяет солверу удовлетворить одному ограничению степени свободы. В чём же тогда проблема? Не слишком-то это было трудно. А проблема, как нетрудно догадаться, в том, что ограничений таких мноого (по меньшей мере по одному на каждый контакт) и удовлетворить им нужно всем одновременно. Чтобы сделать это, нужно каким-то образом решить систему уравнений, состоящую из количества уравнений, равного количеству ограниченных степеней свободы. Обычно эта величина порядка нескольких тысяч. Большинство методов решения системы линейных уравнений работает за сложность порядка O(n³), что просто убило бы производительность, прямой подход «в лоб» применяется редко. Зато существуют итеративные методы, постепенно делающие решение всё более и более точным. Так получилось, что если решать полученную систему уравнений одним из таких методов, называемых Projected Gauss-Seidel Solver, PGS, то последовательность производимых при этом действий математически эквивалентна последовательному решению джойнтов, буквально как описывалось выше. Этот подход широко распространён и называется Sequential Impulses.

То есть, если мы будем просто решать ограниечение за ограничением, не обращая внимание, что оно, вообще говоря, может быть и не одно и будем делать это достаточно долго, то полученный результат будет сходиться к результату, как если бы их решали в системе. За строгими доказательствами — в другие источники, но, уж на слово поверьте, там действительно ничего слишком сложного.

Warmstarting (горячий старт)

Существует множество эвристик, существенно улучшающих сходимость алгоритма Sequential Impulses. Одна из самых важных называется warmstarting, заключается в том, что мы ищем решение системы уравнений каждую итерацию процесса моделирования не с нуля, а используя результаты с предыдущего шага. Эвристика использует свойство так называемой временной когерентности (time coherence) — предполагается, что напряжения, действующие в ограничениях, меняются от кадра к кадру не слишком значительно.

На практике реализуется так:

• сначала применяем аккумулированные на прошлой итерации импульсы во всех ограничениях
• продолжаем решать ограничения, как ни в чём не бывало находя импульсы lambda, и просто прибавляем их к аккумулированному значению.

при этом также немного меняется стратегия, по которой мы накладываем ограничения на величину передаваемого импульса. если в случае без warmstarting'а ограничение на то, что тела могут только разлетаться, не слипаясь, схематично выглядело так:

lambda = ComputeLambda();
if(lambda < 0.0f) return;
ApplyImpulse(lamda);

то с добавлением аккумулирования импульса, стратегия немного меняется:

lambda = ComputeLambda();
accumulatedImpulse += lambda.
if(accumulatedImpulse < 0.0f)
{
  lambda += (0.0f - accumulatedImpulse);
  accumulatedImpulse = 0.0f;
}
ApplyImpulse(lamda);

Физический смысл произошедшего в том, что мы теперь ограничиваем только накопленный импульс, а lambda считаем его изменением. и если мы видим, что какое-то изменение lambda сделало накопленный импульс отрицательным, мы урезаем это изменение таким образом, чтобы оно обращало накопленный импульс ровно в ноль.

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

Случай упругого соударения

Далее вспомним, чему будет равняться скорость, если тела обменяются через данный контакт импульсом с модулем lambda:

  float velocityProjection =
      n1 * (body1->velocity + n1 * body1->invMass * lambda) 
    + w1 * (body1->angularVelocity + w1 * body1->InvInertia * lambda)
    + n2 * (body2->velocity + n2 * body2->invMass * lambda) 
    + w2 * (body2->angularVelocity + w2 * body2->InvInertia * lambda);

Если в случае с абсолютно неупругим соударением мы ставили условие velocityProjection >= 0. Теперь поставим условие следующим образом: velocityProjection = -bounce * initialVelocityProjection. Здесь initialVelocityProjection - величина, равная:

  float initialVelocityProjection =
      n1 * body1->velocity
    + w1 * body1->angularVelocity
    + n2 * body2->velocity
    + w2 * body2->angularVelocity;

Чтобы понять, что это значит, вернёмся к аналогии с двумя сталкивающимися шариками. То есть мы хотим, чтобы при bounce = 1 тела после решения контакта начали бы отдаляться с такой же скоростью, с какой они до его решения сближались, а при bounce = 0, чтобы они относительно друг друга остановились. В остальном алгоритм решения каждого контакта остаётся практически без изменений:
1) Считаем величину b, которая является коэффициентом при lambda и зависит от нормали, точки контакта, массы и момента инерции каждого из тел. Формула для неё остаётся той же самой.
2) Считаем велиичну a, которая зависит от скоростей контактирующих тел и всё тех же параметров контакта. Теперь, когда мы ввели коэффициент bounce, принимает такой вид:

  float a = (n1 * body1->velocity + n2 * body2->velocity
    + w1 * body1->angularVelocity
    + w2 * body2->angularVelocity) + bounce * initialVelocityProjection;

Выражение для initialVelocityProjecton - выше.
3) находим lambda = -a / b;
4) применяем lambda к обоим телам.

Важно отметить, что если используется несколько итераций солвера, то математически корректно будет сначала посчитать для каждого контакта initialVelocityProjection, а потом прогонять несколько итераций солвера, не изменяя это значение.

Усовершенствования алгоритма

Приведённая последовательность действий лишь частично предотвратит взаимопроникновение тел. К сожалению, из-за ошибок интегрирования и итеративного не совсем точного решения системы, тела всё-таки будут частично проникать друг в друга, это может быть особенно заметно, если на тела продолжительно действует сближающая их сила. Простейший случай: бокс, стоящий на столе, при таком подходе будет постепенно в него погружаться. Чтобы предотвратить это явление, существует ряд методик. Одна из самых простых называется velocity-based, основана на том, что мы считаем контакт решённым не при условии velocityProjection >= c а при модифицированном: velocityProjection >= c + contact->depth * ERP.

ERP (error reduction parameter) играет роль чего-то вроде архимедовой силы — чем глубже тела проникают друг в друга, тем быстрее мы хотим, чтобы они искусственно расталкивались. физическая аналогия: что-то в духе архимедовой силы, выталкивающей пенопласт из воды. Грубоватая, правда, аналогия, потому что, во-первых, архимедова сила зависит от объёма погружённой в жидкость части тела (мы учитываем глубину), и архимедова сила влияет на ускорение, а мы воздействуем напрямую на скорость.

К сожалению, это простой, но не слишком эффективный подход, из-за него глубоко проникшее тело «выплёвывается» из того, в которое оно проникло и может неестественно далеко улететь. Куда более практичный подход — использование так называемых псевдоимпульсов и псевдоскоростей:

Pseudovelocities (псевдоскорости)

Псевдоскорость — искусственно введённый в систему параметр, служащий для расталкивания тел. В случае контакта, вместо использования одного уравнения мы поступаем следующим образом:

1) составляем как описывалось выше систему контактов без velocity-based добавки:
velocityProjection >= c;
2) решаем, как ни в чём не бывало.
3) обнуляем псевдоскорости всех тел.
4) составляем ещё одну систему, на этот раз не для скоростей, а для псевдоскоростей:
pseudoVelocityProjection >= contact->depth * ERP;
5) решаем в точности таким же подходом систему для псевдоскоростей, предполагая, что это настоящие скорости тел.
6) интегрируем движение каждого тела для обычных скоростей по закону
velocity += acceleration * timeStep; pos += velocity * timeStep;
7) интегрируем движение каждого тела для псевдоскоростей по закону
pos += pseudoVelocity * timeStep;

Физический смысл подхода не слишком прозрачен, но состоит примерно в том, что мы наделяем тела дополнительными мгновенными скоростями, которые обнуляются, как только тела расталкиваются друг из друга. В результате чего они раздвигаются, но не разлетаются. Вообще говоря, они обнуляются каждый кадр, а не только в момент, когда одно тело вытолкнется из другого, но сути это не меняет. Хочется отметить, что попытки ввести в систему для псевдоскоростей warmstarting добром редко увенчиваются, потому что система и без него сходится достаточно быстро, а стабильность warmstarting понижает значительно. Это обозначает, что warmastarting имеет смысл применять для обычной системы импульсов, а систему для псевдоскоростей решать каждый кадр заново, без warmstarting'а.

#gauss-seidel, #pgs, #основы, #солверы

5 января 2010 (Обновление: 9 авг 2020)

Комментарии [644]