StarDEM: эффективный метод дискретных элементов для частиц звездообразной формы

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

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

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

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

В нескольких работах [LCH20; Wan*21; LZZH22] моделируются твердые частицы звездообразной формы. Однако их оценка перекрытия между двумя частицами, хотя и эффективна для частиц, близких к сферическим (обычно представляющих собой природные зерна или горные породы), не очень подходит для частиц с сильной вогнутостью или удлиненной формы. Наш метод StarDEM фокусируется на стабильном и эффективном моделировании звездообразных частиц и может работать с большим разнообразием несферических частиц, в том числе с сильно невыпуклой геометрией.

Наш вклад заключается в следующем:

• Мы предлагаем использовать аппроксимацию расстояния первого порядка для эффективной оценки расстояния между точкой и границей звездообразной формы.
• Мы реализуем эффективное 2D-моделирование на CPU и GPU, предназначенное для этого изготовленного гранулированного материала, в частности материалов, состоящих из сильно невыпуклых частиц.

Мы предоставляем открытый исходный код для нашей реализации здесь https://github.com/schreckc/StarDEM.

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

Метод дискретных элементов (DEM) является эффективным способом моделирования большого количества частиц, каждая из которых представлена в виде отдельного объекта с определенным положением и ориентацией. В оригинальной статье Кандалла и соавт. [CS79] жесткие сферы моделируются с использованием мягких столкновений: частицам разрешается слегка накладываться друг на друга. Учитывая две сферы S₁ и S₂ с центрами c₁, c₂ и радиусами r₁, r₂, нормальное усилие между ними зависит от глубины перекрытия

и контактной нормали

Также могут быть вычислены силы трения, а также тангенциальные силы, которые возникают из-за разницы тангенциальных скоростей между двумя частицами в точке контакта xc (центр перекрытия). Мы отсылаем читателя к [CS79] за дополнительной информацией. Затем в точке xc прикладывается сила F12 к S₁, и результирующие сила и импульс, действующие на центр тяжести S₁, равны

В cилы приложены аналогичным образом на секциях S₂ с

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

За прошедшие годы многие статьи расширили оригинальную DEM, и, в частности, некоторые распространили метод на несферические гладкие частицы; мы отсылаем читателя к [Fen23] за подробным недавним обзором. Наиболее известный метод представляет форму частицы в виде скопления сфер [BYM05]. Этот метод использует преимущества эффективности обработки столкновений сфер. Однако для точного представления гладких, сложных форм требуется большое количество сфер на частицу, что делает метод дорогостоящим в этих случаях.

Также были исследованы другие методы, использующие неявные функции, такие как суперэллипсоид и поля расстояний со знаком (SDF). Среди них есть несколько работ, посвященных звездообразным частицам, радиальная функция которых представлена сферическими гармониками (или рядами Фурье в 2D). Звездообразные формы представляют собой широкий диапазон нетривиальной невыпуклой геометрии (см., Например, рисунок 1), которая может привести к интересному физическому поведению.

Форма звезды определяется как фигура, содержащая любой сегмент между центром фигуры и любой граничной точкой. Затем можно определить двумерную звездообразную форму с гладкой границей, используя центр c и дифференцируемую 2π-периодическую радиальную функцию r(α), связывающую радиус с углом α вокруг центра.

где α(p) − угол, образованный (p - c), и направление, соответствующее α = 0. Таким образом, точка p находится внутри фигуры, если f(p) < 0, и снаружи, если f(p) > 0. Граница может быть выбрана путем равномерной выборки угла α:

Эти свойства позволяют эффективным алгоритмам обнаруживать контакты между двумя звездообразными формами [CH21; GB13] путем нахождения отобранных узлов одной частицы внутри другой.

Только в нескольких статьях рассматривается моделирование DEM звездообразных частиц. Ближе к нам Ван и др. [Wan*21] рассмотрим схему обнаружения контактов для трехмерных звездообразных форм. Они предлагают вычислять глубину перекрытия, используя радиальное расстояние, которое является приличным приближением евклидова расстояния между одной точкой и границей формы звезды, только если форма достаточно близка к сферической. Работа Лая и др. [LZZH22] моделирует более общие формы с SDF, но в случае сферических гармоник они страдают от той же проблемы, поскольку они приближают SDF к радиальному расстоянию. Наконец, Лай и др. [LCH20] предлагают оптимизацию замкнутой формы для нахождения точек перекрытия и пересечения границ двух звездчатых форм в 2D. Однако в их работе не рассматриваются случаи, когда между двумя частицами имеется несколько пересечений, что часто бывает с невыпуклыми формами.

Аппроксимация расстояния первого порядка

Чтобы вычислить перекрытие, мы должны оценить расстояние d между точкой и границей формы звезды S, определяемой центром c и радиальной функцией r. В предыдущей работе использовалось радиальное расстояние вдоль луча, проходящего через p и центр: d_r = f(p). Это разумное приближение до тех пор, пока нормаль к границе близка к направлению луча из центра. К сожалению, это не относится к вытянутой или сильно невыпуклой частице.

Мы предлагаем вместо этого использовать приближение расстояния первого порядка, используемое Таубином [Tau94] для растрирования неявной кривой. В этом приближении используется разложение неявной функции по Тейлору первого порядка. В нашем случае граница определяется неявной функцией f (Уравнение 1). Расстояние между точкой p и границей может быть аппроксимировано с помощью:

где ∇f является аналитическим градиентом f. Это расстояние асимптотически эквивалентно евклидову расстоянию: ошибка при евклидовом расстоянии приближается к нулю, когда p приближается к границе. Таким образом, ошибка невелика при приближении к границе (рис. 3), где вычисляются мягкие столкновения, учитывая, что предполагается, что они допускают только небольшое перекрытие.

Это расстояние приближения первого порядка между точкой и границей S таким образом, может быть вычислено с использованием только одной оценки f и его градиента. Нормаль поверхности вблизи точки p ∇f затем может быть аппроксимирована нормализованным градиентом ||∇f||. Расстояния первого порядка и радиальные расстояния оцениваются за постоянное время. Оценка первого лишь немного дороже, чем второго.

Действительно, это радиальное расстояние, деленное на норму градиента f, которая также используется для вычисления нормали. Несмотря на то, что возможны аппроксимации d более высокого порядка [Tau94], мы нашли наше моделирование стабильным, используя только аппроксимации первого порядка.

Обнаружение столкновений

Мы используем схему обнаружения столкновений, аналогичную Wang et al. [Wan*21]. Ограничивающий круг центрирован в геометрическом центре частицы. Его радиус равен максимальному радиусу частицы. Если существует пересечение между ограничивающими окружностями двух частиц S₁ и S₂, область возможного контакта ограничивается областью пересечения двух дисков. Обозначим точки пересечения двух окружностей I и J (см. рисунок 2).

Как Ван и др. [Wan*21], затем мы используем схему "узлы-поверхность" для вычисления возможных контактов между S₁ и S₂. Мы находим выбранные узлы из S₁, которые находятся внутри границы S₂, используя знак f. Луч, начинающийся из центра формы звезды, пересекает ее границу только один раз. Итак, нам нужно только проверить граничные узлы x из S₁, соответствующие углу α(x), который принадлежит интервалу [α(I), α(J)], ограниченному точками пересечения I и J (фиолетовая область на рисунке 2 (a)). Каждый набор последовательных узлов внутри S₂ соответствует отдельному контакту между S₁ и S₂.

Насколько нам известно, предыдущие работы ([LZZH22; Wan*21]) с использованием методов столкновения узлов с поверхностью используют только узлы из S₁ (красные узлы на рисунке 2) для разрешения столкновения. Однако это приводит к асимметричному столкновению, при котором силы изменяются, если S₁ и S₂ меняются местами. Следовательно, мы предлагаем также рассмотреть граничные узлы из S₂. Следовательно, мы предлагаем также рассмотреть граничные узлы из S₂. Возможно, это связано с использованием таких геометрических функций, как: узлы первого и последующего узлов 1 отдельного контакта (красные узлы на рисунке 2(а) представляют два контакта) определяют угловую область в S₂ В (синяя область в (в)). Мы выбираем узлы из S₂ внутри этой области (синие узлы). Хотя это выполняется быстро для центрального процессора, выполнение этого становится дорогостоящим для графического процессора. Итак, мы рассматриваем только красные узлы из S₁ для реализации графического процессора.

Разрешение столкновений

Для каждого контакта между S₁ и S₂ мы выбираем в качестве точки контакта x_c самый глубокий узел (желтые узлы на рисунке 2 (c)) от одной из частиц внутри другой. Глубина d перекрытия для этого контакта тогда равна расстоянию между x_c и границей другой частицы. d вычисляется с использованием расстояния первого порядка аппроксимации, описанного выше. Нормаль n контакта оценивается как градиент f в точке x_c. Затем используются значения x_c, d и n , вычисленные для каждого контакта, как описано в разделе 2.

Обратите внимание, что, хотя в большинстве методов столкновения используется одна точка соприкосновения для одного отдельного контакта, в предыдущих работах с использованием методов "узлы-поверхность", как правило, рассматривались все узлы S₁ внутри S₂ в качестве точек соприкосновения и добавлялись силы для каждой из них. Мы решили использовать самую глубокую точку, поскольку это заставляет силу зависеть только от глубины перекрытия и не зависит от количества выборочных узлов внутри. Однако относительная физическая точность двух методов требует дальнейшего изучения. Кроме того, силу и импульс необходимо вычислять только один раз для каждого контакта методом самой глубокой точки.

Эффективная реализация

Мы специально оптимизируем реализацию нашего метода для моделирования гранулированного материала, состоящего из идентичных частиц.

Учитывая, что все частицы имеют одинаковый размер, мы используем пространственную сетку с ячейками размером 2r_max. Каждой частице присваивается ячейка в зависимости от ее положения. Для любой данной частицы нам нужно только изучить потенциальное столкновение с частицами в пределах квадрата ячейки 3x3, окружающего положение частицы.

Расстояние первого порядка и, следовательно, радиальная функция r и ее производная r' часто вычисляются и могут быть дорогостоящими, например, при использовании многих коэффициентов Фурье. Чтобы обойти эту проблему, мы предварительно вычисляем r и r' в таблицах подстановки. Эта замена вычисления функции доступом к таблице упрощает процесс и обеспечивает постоянный доступ по времени. Важно отметить, что количество выборок для справочной таблицы отличается от количества, используемого в методе nodesto-surface. Учитывая, что мы выполняем выборку одномерной кривой и что нам нужно только один раз вычислить и сохранить таблицы, мы можем позволить себе выборку с высоким разрешением. В наших примерах мы используем 20000 образцов, гарантируя точные результаты.

Реализация на графическом процессоре аналогична. Мы строим сетку с использованием атомарных операций, создавая связанный список частиц в каждой записи сетки. Таблица 2D-поиска -- "спрайт" частицы -- используется для проверки того, находятся ли точки S₁ в S₂ и одновременного получения расстояния столкновения и нормали.

Результаты и тайминги нашего процессора были рассчитаны на 12-ядерном процессоре AMD Ryzen 9 3900X. Мы распараллелили моделирование с использованием OpenMP. В наших примерах используются частицы радиусом около 5 мм и весом 0,2 г. Нормальная жесткость составляет 1e3 Нм⁻¹, а касательная жесткость 0.5e3 Нм⁻¹. Граница частицы выбирается из 100 узлов. В нашем методе может использоваться любая дифференцируемая функция r. В нашем примере мы используем сумму Фурье для представления r:

Это можно обобщить в 3D с использованием сферических гармоник. Пожалуйста, обратитесь к прилагаемому видео, чтобы ознакомиться с результатами моделирования, описанными в примерах, приведенных в этом разделе.

Аппроксимация первого порядка в зависимости от радиального расстояния

Мы сравниваем аппроксимацию первого порядка d_fo евклидова расстояния с численно вычисленной основной истиной d_gt. Ошибка e = |d_fo - d_gt|/d_gt остается низкой вблизи границы (рисунок 3, внизу). Например, в случае фигуры на рисунке 3 с радиусом r_max = 1 ошибка на расстоянии < 0.01 границы остается < 0.17.

Радиальное расстояние d_r = f(p) = ||p - c|| - r(α(p)) только близко к евклидову расстоянию, когда нормаль к границе близка к направлению луча из центра. По сравнению с расстоянием истинности на земле ошибка может стать очень большой для вытянутых областей даже вблизи границы (рисунок 3, вверху). Ошибка на расстоянии < 0.01 от границы может достигать 2.7. Как показано на видео, это может вызвать нестабильность при моделировании.

Сравнение с скоплением сфер

Мы реализовали метод Кэндалла и Стрэка для сфер [CS79] и метод Белла для скопления сфер [BYM05]. Реализации те же, за исключением обработки коллизий. Мы сравниваем методы с использованием упаковки: в коробку попадает 3154 частицы (рисунок 4). Каждая частица имеет площадь 30 мм².

Мы измеряем среднее время вычислений в одном потоке для обработки коллизий (обнаружение и вычисление перекрытий). В случае сфер (a) использование StarDEM с r(α) = 1 выполняется примерно в девять раз медленнее, чем использование прямых сфер (9,6 мс против 1,1 мс за один шаг столкновения). Однако при использовании скопления сфер для представления более сложной формы стоимость использования сфер квадратично возрастает с увеличением количества сфер, необходимых для каждой частицы.

Для формы креста, использующей форму звезды с r(α)= 1+cos(4α) (b), время вычисления составляет 26 мс, в то время как при использовании скопления сфер (c) увеличивается до 12 мс для 9 сфер и 41 мс для 17 сфер. Более сложные формы (например, (d), (e)) потребовали бы большого количества сфер и увеличения вычислительных затрат для группы сфер, в то время как временные интервалы для нашего моделирования остаются аналогичными ((d) 30 мс, (e) 29 мс).

Сравнение с [LCH20]

Мы сравнили 2D-тайминги с методом Лай и др. [LCH20]. Они сообщают, что время вычисления составляет около 5 секунд для 1e4 возможных столкновений между парами частиц, что приводит к фактическому столкновению примерно в половине случаев. Для сравнения, для последних 100 кадров нашего примера упаковки на рисунке 4 (e) мы обнаруживаем 84.9e4 возможных столкновения (пересекающиеся ограничивающие круги), в результате чего 51e4 пары фактически сталкиваются и 84.7e4 отдельных контакта. Это вычисляется за 3 секунды для одного потока и 0,6 секунды для 8 потоков.

Другие результаты

Мы продемонстрировали наш метод на различных примерах с разными формами. Мы используем временной шаг ∆t = 0.1 мс. Сообщаемые значения времени для одного шага моделирования усредняются по всему моделированию и вычисляются в 8 потоках. Частицы, используемые для этих значений времени, соответствуют тем, которые показаны на рисунках, однако для разных форм время обычно довольно схоже.

Мы показываем пример неплотной упаковки (рисунок 1 (слева)) с частицами, падающими в коробку (3154 частицы, один этап моделирования занимает в среднем 8 мс). На рисунке 1 (посередине) показаны частицы внутри вращающегося барабана (900 частиц, 2 мс). Частицы, выходящие через узкий зазор, имеют тенденцию образовывать перемычки, как показано на рисунке 1 (справа) (279 частиц, 1 мс) -- на видео также показана силовая сеть, созданная в этом примере. Последний пример в видео показывает, как колонна сворачивается и образует кучу (429 частиц, 2 мс).

Графический процессор

На NVidia GTX 4080 полный этап моделирования с соприкасанием 1 миллиона вогнутых частиц занимает всего 7,5 мс (см. Видео).

Хотя в этой статье мы проводили оптимизацию для идентичных частиц, суть нашего метода также может быть использована для гетерогенных частиц. В случае большого количества различных форм затраты памяти на использование справочных таблиц могут стать непомерно высокими (тем более для 2D справочных таблиц, которые мы использовали для графического процессора или которые использовались бы в аналогичном методе с использованием растрированного SDF). Эта проблема станет более острой при переходе к трехмерным звездообразным частицам.

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

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

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

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