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

Расчет выталкивающей силы является фундаментальным методом, необходимым для динамического моделирования объектов, плавающих или погруженных в жидкостях. Взаимодействие между твердыми телами и жидкостями называется проблемой взаимодействия твердого тела и жидкости, подробный обзор этой темы можно найти в [BPD20]. Соединение твердого тела с жидкостью - исключительно сложная задача, включающая два основных аспекта: определение движения тела, погруженного в жидкость, и, наоборот, оценка воздействия движущегося тела на жидкость.

В этом кратком исследовании мы сосредоточимся исключительно на первой проблеме. Существующие решения часто работают на сетках [BR13] [Cla13] или непосредственно на сетке, представляющей объект [KKKT06] [Roy23]. Эти методы решают более сложные задачи, чем наши, но требуют больших вычислительных затрат. Напротив, мы нашли множество реализаций, которые могли быть вычислены в режиме реального времени на центральном процессоре, но основывались на крайне неточных приближениях.

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

Рассмотрим твердое трехмерное тело, (частично) погруженное в объем жидкости. Предположим, что наше тело жесткое и имеет однородную плотность. Мы также предполагаем, что жидкость несжимаема, ее поверхность находится в состоянии покоя и достаточно велика. Для дальнейшего упрощения задачи мы проигнорируем все силы сопротивления как в воздухе, так и в жидкости.

В этом случае движение объекта определяют две силы: сила тяжести и выталкивающая сила, как показано на рисунке 1. Силу тяжести можно рассчитать как:

где m - масса объекта, g - гравитационная постоянная, а j = (0,1,0) - единичный вектор, направленный вверх. Согласно принципу Архимеда, выталкивающая сила рассчитывается как:

где ρ - плотность жидкости, а V_B - объем погруженной части объекта, который равен объему вытесненной жидкости.

Рисунок 1: Плавающее движение. Сила тяжести направлена вниз от центра тяжести, который является центроидом всего тела. Выталкивающая сила направлена вверх от центра плавучести, который является центроидом погруженной части тела.

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

В среде динамического моделирования масса и центр масс объекта обычно известны заранее. Их следует указать или предварительно рассчитать, поскольку они остаются постоянными во время движения. Напротив, V_B и C_B изменяются в каждый момент времени. Для моделирования движения объекта важно понимать, как вычислять V_B и C_B для данного жесткого преобразования.

Рассмотрим множество Ω ⊂ R³, которое определяет точки твердого тела в пространстве. Объем и центроид тела определяются как интегралы по объему. Например, объем Ω и x-координата центроида равны:

соответственно. Согласно теореме Гаусса о дивергенции (при соответствующих условиях), эти объемные интегралы могут быть заменены поверхностными интегралами, вычисленными по границе [MT12], например:

Если поверхность тела представлена треугольной сеткой, интегралы могут быть вычислены треугольник за треугольником. Таким образом, и объем, и центроид могут быть определены взвешенным суммированием элементов объема со знаком V_i и элементов центроида C_i для каждого треугольника Δ_i:

Элементы объема и центроида могут быть рассчитаны по формулам:

как мы можем прочитать, например, в [Cat06]. Здесь a_i, b_i, c_i обозначают три вершины треугольника Δ_i по порядку. Мы предположили, что треугольники имеют согласованную ориентацию, а нормаль, обращенная наружу, задается как (b_i-a_i)×(c_i−a_i). Обратите внимание, что элементы объема V_i представляют собой подписанные объемы тетраэдра, определяемые a_i, b_i, c_i и произвольно выбранной точкой отсчета o. Аналогично, элемент центроида C_i соответствует центру масс того же тетраэдра.

Предположим, что текущее положение и ориентация движущегося тела могут быть описаны жестким преобразованием T. Рассмотрим произвольную вершину x сетки в состоянии покоя. Для каждого временного шага в моделировании существует поворот R и перемещение t, которые описывают текущее местоположение вершины x: x' = T(x)= Rx + t. В каждый момент моделирования мы будем выполнять следующие шаги.

1. 1 Обновите текущее жесткое преобразование, Т.
2. 2 Классифицируйте треугольники преобразованного тела (непогруженные, частично погруженные, полностью погруженные).
3. 3 Вычислите вклад каждого треугольника в выталкивающую силу.
4. 4 Передайте силу и крутящий момент физическому симулятору.

Наши алгоритмы были реализованы с использованием платформы разработки в реальном времени Unity. Сама Unity включает физический движок, поэтому мы полагались на него при моделировании динамики твердого тела. В частности, мы использовали встроенный компонент Rigidbody для присвоения массы объекту. Кроме того, в нашем моделировании этот компонент рассчитывает влияние силы тяжести и сопротивления.

Шаг 1, определяющий текущее жесткое преобразование, может быть выполнен с помощью встроенного преобразования Unity.Функция TransformPoint. Затем шаг 2 можно выполнить простым способом. Шаг 4 может быть выполнен просто с помощью встроенного Rigidbody в Unity.Добавьте функции Force и Rigidbody.Добавьте функции Torque. Таким образом, большая часть динамического моделирования выполняется Unity, в то время как наш специально разработанный компонент изменяет движение объектов, добавляя силы и крутящие моменты.

Наиболее важной задачей, требующей решения, является Шаг 3: вычисление выталкивающей силы. Далее мы представляем 5 различных алгоритмов, направленных на решение этой проблемы. Первый обеспечивает быструю, но грубую аппроксимацию, второй предлагает более точную аппроксимацию, а третий вычисляет точные значения, последние два предлагают дальнейшие возможности развития. Разница между этими алгоритмами заключается в том, какой вклад треугольников учитывается при их вычислениях, как показано на рисунке 2.

Рисунок 2: Быстрые алгоритмы учитывают только полностью погруженные треугольники (синий). Улучшенный алгоритм также включает в себя вклад треугольников на поверхности жидкости, которые окружают затопленную область сверху (зеленый). Точные алгоритмы включают, в дополнение к вышесказанному, вклад частично погруженных треугольников (красный).

1. Быстрый приближенный алгоритм

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

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

2. Улучшенный приближенный алгоритм

Быстрый алгоритм может привести к неточным результатам. В таких случаях вклад пропущенных треугольников становится значительным. Важно отметить, что это упущение касается не только частично погруженных треугольников. Формулы для объема и центроида основаны на интегралах по замкнутым поверхностям, что означает, что вклад многоугольников, образованных на пересечении поверхности жидкости (плоскости) и сетки, также должен быть включен в наши расчеты.

Простой подход включал бы вычисление этого пересечения, триангуляцию результирующих полигонов (с потенциальными дырами), а затем учет вклада этих фрагментов поверхности. Однако эта процедура была бы дорогостоящей с точки зрения вычислений [Roy23].

Ключевая идея для решения этой проблемы заключается в следующем. Напомним, что при расчетах как для объемных, так и для центроидных элементов выбор начала координат o является произвольным. Если в каждый момент времени мы корректируем начало координат так, чтобы оно лежало на поверхности жидкости, то объемный вклад всех элементов поверхности на плоскости становится равным нулю, поскольку четыре вершины любых тетраэдров копланарны. Следовательно, пусть o - проекция текущего центра тяжести объекта на плоскость y = 0 и вычислите объем и элементы центра тяжести относительно этой контрольной точки.

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

3. Точный алгоритм

Точный алгоритм представляет собой небольшое усовершенствование улучшенного алгоритма, в котором также учитывается вклад треугольников, пересекающих поверхность жидкости. Исключая краевые случаи (когда одна или несколько вершин треугольника лежат точно на плоскости y = 0), возможны два сценария: частично погруженный треугольник имеет либо одну, либо две вершины ниже поверхности жидкости. В первом случае мы рассматриваем вклад подтреугольника под поверхностью жидкости. Во втором случае четырехугольник, образованный под поверхностью жидкости, разделяется на два треугольника, и их вклады добавляются к существующим, см. Рисунок 3.

Рисунок 3: Непогруженные (белые), полностью погруженные (синие) и частично погруженные (красные) треугольники. Справа: два основных случая для частично погруженных треугольников: если только вершина b находится ниже поверхности жидкости, учитывается вклад треугольника pbq; если c и a находятся ниже поверхности, также учитываются вклады треугольников pqc и apc.

4. Алгоритмы, основанные на поверхности

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

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

Следуя предыдущему подходу, вычисления можно выполнять треугольник за треугольником, суммируя вклад сил и крутящих моментов для каждого треугольника. После простых вычислений мы получаем следующую формулу для элементов силы: F_i = ρg(C_i)_yA_i, где (C_i)_y обозначает y-координату центра тяжести треугольника, A_i - элемент ориентированной поверхности.

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

Точный алгоритм учитывает вклад каждого частично или полностью погруженного треугольника (включая те, которые лежат на поверхности жидкости), обеспечивая точные результаты, не считая ошибок в числовых расчетах. Мы проверили эту точность на двух простых примерах, а также оценили производительность двух других алгоритмов с точки зрения точности.

Для тестов мы исследовали погружение единичной сферы и единичного куба. Центроиды объектов были концептуально размещены в точке (0,h,0), и мы определили, как изменяется величина выталкивающей силы в зависимости от h. Для простоты мы предположили, что ρg = 1. Основные геометрические соображения привели к прямолинейному предварительному выводу. В случае куба для h ∈ (-1/2, 1/2) получаем:

Для сферы задача была лишь немного сложнее, требуя интегрирования полиномиальных выражений. В этом случае для h ∈ (-1, 1) имеем:

Во время наших тестов мы использовали кубические сетки двух разных уровней детализации. Куб из 10³ состоял из 1200 треугольников, в то время как куб из 100² состоял из 120 000 треугольников. Дополнительно мы протестировали две сетки, приближенные к сфере: сфера, аппроксимированная 50² сегментами, содержала 5000 треугольников, тогда как сфера с 200² сегментами содержала 80 000 треугольников.

Рисунок 4: Тесты точности наших алгоритмов в случае погружения единичного куба и единичной сферы.

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

Мы провели тесты производительности с использованием нескольких моделей, три из которых представлены здесь. Модель torusknot содержит 1760 треугольников, модель Stanford Bunny - 5,002, а модель Armadillo - 30,000 треугольников. Во время измерений мы зафиксировали время срабатывания 5 алгоритмов, разработанных для вычисления выталкивающих сил.

В начале теста объект находился почти полностью над поверхностью жидкости, а к концу он был почти полностью погружен в воду. После выполнения 100 измерений мы рассчитали среднюю стоимость вычисления выталкивающей силы и крутящего момента. Измерения проводились на настольном компьютере со следующей конфигурацией: Intel Core i5-9400, 8 ГБ оперативной памяти, Intel UHD Graphics 630, работающем под управлением Windows 10 Pro. Результаты приведены в таблице 1.

Даже для детализированной модели Armadillo измеренное время отклика в наихудшем случае составило примерно 10 мс. Это обеспечивает скорость моделирования 100 кадров в секунду на тестовом компьютере, обеспечивая производительность сложных моделей в режиме реального времени.

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

Из предыдущих методов [BPD20] заслуживает особого внимания. Авторы рассматривали более сложную проблему, чем наша; однако их расчет выталкивающей силы сопоставим с нашим. Для достижения этой цели они аппроксимируют погруженный объем аналогично нашим методам аппроксимации на основе объема. Аналогичное измерение было проведено для проверки точности, подтвердив, что их метод лишь приближается к теоретической кривой.

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

Компонент Unity, реализующий продемонстрированные алгоритмы, можно найти в репозитории [Fab25] на GitHub.

Поддерживается Стипендиальной программой EKÖP-24 для повышения квалификации в университетах Министерства культуры и инноваций за счет средств Национального фонда исследований, разработок и инноваций.