CC BY
60
Известия Института математики и информатики УдГУ
2015. Вып. 2 (46)
УДК 519.63
© А. С. Караваев, С. П. Копысов, А. С. Сармакеева
ПОСТРОЕНИЕ АГЛОМЕРАТОВ СФЕР В МЕТОДЕ ДИСКРЕТНЫХ ЭЛЕМЕНТОВ
В работе предложен алгоритм построения агломерата сфер при аппроксимации геометрии тел произвольной формы в методе дискретных элементов. Для оценки получаемой поверхности произвольных тел введена мера, характеризующая микрорельеф поверхности сформированного агломерата.
Ключевые слова: метод дискретных элементов, агломерат сфер, оценка погрешности.
Введение
Метод дискретных элементов (МДЭ), широко используемый для исследования деформационных процессов в гранулированных (сыпучих), слабосвязанных средах, впервые был предложен в работе [1], где среда представлялась через ансамбль взаимодействующих сфер. Методом дискретных элементов также моделируют динамику твердых тел сложной формы [2], и в этом случае каждое тело состоит из множества жесткосвязанных сфер и представляется как муль-тисферный агломерат.
Существует несколько алгоритмов построения агломератов сфер. Построение агломерата из поверхностной треугольной сетки исследовано в [3]. При таком построении каждая сфера определяется четырьмя узлами поверхностной сетки путем решения, соответствующего уравнения второго порядка. В этом случае неизбежно появление сфер, находящихся за пределами модели, и необходимо применение специальной двухшаговой процедуры для их удаления.
В [4] исходной информацией для построения является множество точек, распределенных на поверхности моделируемого тела. Объем заполняется сферами разного радиуса, а построение начинается от случайно выбранной точки на поверхности, далее диаметр сферы растет вдоль направления внутренней нормали точки до достижения им максимального размера, когда сфера пересекает другую точку на поверхности модели. Процесс выбора новых точек регулируется расстоянием до уже построенных предыдущих сфер агломерата. К недостаткам этого алгоритма можно отнести сложность поиска максимального диаметра вписанной сферы и, как следствие, непростую программную реализацию метода. Отметим, что качество модели возрастает с увеличением плотности точек на поверхности. Развитием данного подхода можно считать [5], где также используется построение сфер из множества точек, описывающих поверхность моделируемого тела, которая предварительно покрывается сеткой с кубических ячейками. Добавление сфер в агломерат идет последовательно, начиная со сфер, охватывающих максимальное число кубических ячеек внутри тела. Однако построение сетки, содержащей кубические ячейки, для сложных тел является нетривиальной задачей. В этом случае в качестве исходных данных может быть использована геометрическая сетка, полученная путем компьютерной томографии реального физического тела [6]. Сетка характеризуется наличием внутренних и граничных узлов. Центры и радиусы сфер агломерата вычисляются путем анализа минимального расстояния от внутренних узлов до узлов на поверхности. Производительность предлагаемого метода существенно снижается для тел с развитой криволинейной поверхностью, когда для описания модели потребуется значительное число граничных узлов.
В работе [7] рассматривается метод построения агломерата из объемных данных в формате «знакоопределенных полей расстояний», когда каждый воксель сохраняет наименьшее расстояние до поверхности объекта с учетом знака. Построение агломерата осуществляется последовательно, начиная со сфер, имеющих максимальный радиус, аппроксимация модели контролируется количеством сфер, а также плотностью их наложения друг на друга.
Также, основываясь на объемных данных в формате «знакоопределенных полей расстояний», в работе [8] предложен метод построения шестигранных конечноэлементных сеток и алгоритм конвертации начальной поверхностной треугольной сетки, описывающей геометрическую
модель, в формате «знакоопределенных нолей расстояний». Метод характеризуется простой программной реализацией, высокой производительностью, а также возможностью применения параллельных вычислений для больших моделей.
В данной работе исследуются возможности предложенных в [8] алгоритмов построения агломератов сфер из объемных данных. Исходные объемные данные получались как из вексельного представления, так и преобразованием поверхностной сетки, описывающей геометрическую модель. Приводятся результаты работы алгоритмов для различного числа сфер и плотности их распределения. Вводятся меры качества для оценки полученных агломератов сфер.
§ 1. Построение агломерата сфер
Для представления модели в виде агломерата сфер V применяется декартова сетка векселей V в формате «знакоопределенных полей расстояний», представляющая собой трехмерный массив векселей, содержащих наименьшее расстояние до границы модели с учетом знака. Любая сфера р с центром в вокселе V и радиусом г, равным значению вокселя ¡и, гарантированно находится внутри моделируемого объекта. Обозначим через Л список вокселей внутри модели, при этом будем полагать, что элементы в нем отсортированы по значениям ¡и от большего к меньшему.
Построение агломерата V происходит последовательно из отсортированного списка вок-Л
р € V, построенная из вокселя V € Л, определяется центром, расположенным в вокселе, и радиусом Тр, равным значению вокселя ¡и.
Воксель V € Л пропускается алгоритмом, если полученная из него сфера р' будет целиком содержатся в одной из уже построенных сфер агломерата V, то есть Зр € V: р' С р. Кроме того, сфера не строится, если расстояние между V' € Л и центром рр одной из сфер р € V меньше определенной доли ее радиуса Зр € V: 11V — рр||/ гр < Я, где Я — значение, характеризующее расстояние между центрами сфер агломерата V. В случае Я ^ 1 центр рр любой сферы не может располагаться внутри других сфер, при этом для Я < 1 возможно расположение точки рр внутри другой сферы па расстоянии не ближе Я относительно ее радиуса. Уменьшение параметра Я позволяет создать более плотное наложение сфер друг на друга, тем самым обеспечив сглаженную поверхность модели. Предлагаемый алгоритм включает следующие основные шаги (см. алгоритм 1).
а) б) в)
Рис. 1. Модель шестигранного тела, представленная агломератом сфер: а) Я = 0.85, Р = 200; б) Я = 0.85, Р = 1000; в) Я = 0.5, Р = 15000
В реализации алгоритма предусмотрена возможность определения максимального числа элементов Ртах, содержащихся в V, так как с уменьшением параметра Я число сфер Р в агломерате возрастает.
На рисунках 1 и 2 представлены аппроксимации моделей куба и расчетной шестигранной сетки модели позвонка человека соответственно. Характеристики сферических агломератов приведены в таблице 1, где и V т — объемы агломерата V и шестигранной сетки Т,
Алгоритм 1 Построение агломерата сфер P(dim, R) трехмерного тела
Require: Сетка вокселей V размерности
dim V = dimz х dimy х dimx на расчетной области
1: Инициализировать список Л = 0
2: Инициализировать агломерат P = 0
3: for z = 1; z < dimz; z++ do
4: for y =1; y < dimy; y++ do
5: for x = 1; x < dimx; x++ do
6: if fV( ) > 0 then
7: Добавить V(x y z) в список Л
8: end if
9: end for
10: end for
11: end for
Л
13: for all v e Л do
14: if P > Pmax then
15: break
16: end if
17: Инициализировать сферу < : pv' = v,rv' = fv
18: if 3ip G P, <f' С <f V 3ip G P, < 7г then
19: continue
20: end if
21: Добавить <' в P
22: end for
построенных для некоторой рассчитываемой области. Параметр Qv находится как соотношение Vр/Vт и оценивает точность аппроксимации объема модели, вычисленного для агломерата сфер и соответствующей шестигранной сетки. Значение Qr — гладкость получаемой поверхности. Отметим, что относительная разница объемов агломерата и модели не превышает 10%, что дает достаточно точную аппроксимацию области, при этом значение Qv уменьшается с увеличением количества сфер в агломерате.
Таблица 1. Характеристики мультисферных агломератов.
Характеристика P = 200 Шестигранник P = 1000 P = 15000 P = 200 Позвонок P = 1000 P = 15000
R 0.85 0.85 0.5 0.85 0.85 0.5
V P 58.572 61.1457 63.018 31601.4 35338.8 37948.4
V T 64 64 64 39068.6 39068.6 39068.6
Qv 0.91 0.95 0.98 0.81 0.90 0.97
Qr 0.92 0.82 0.12 0.82 0.56 0.03
§ 2. Масса и центр инерции агломерата
Объем V и, соответственно, масса агломератов сфер m = pV вычисляются на сетке равномерно распределенных точек методом Монте-Карло. Центр инерции pc агломерата определяется как центр инерции системы точек, попавших внутрь агломерата при вычислении его объема.
Для увеличения точности вычисляемого объема шаг распределения точек h = k • h, k e (0,1] выбирался не превышающим шаг h декартовой тетки вокселей V, из которой был построен агломерат. В рассмотренных далее примерах принималось k = 0.5. При известных координатах p0) Pdtm первого и последнего вокселя V и сетки размерности dim V = dimz х dimy х dimx
P
алгоритме 2.
Рис. 2. Модели отдельного позвонка позвоночного столба: а) конечно-элементная сетка шестигранных ячеек; агломераты сфер для МДЭ: б) Я = 0.85, Р = 200; в) Я = 0.85, Р = 1000; г) Я = 0.5, Р = 15000
Представленный алгоритм позволяет с заданной точностью определять инерциальные свойства агломерата сфер. Отметим, что наряду с простотой метода он обладает значительными вычислительными затратами. При практической реализации алгоритма необходимо использовать технологии параллельных вычислений.
§ 3. Оценка качества построения поверхности агломерата сфер
Другими значимыми параметрами агломератов при рассмотрении динамики соударения твердых тел являются параметры формы тел такие, как [9] сферичность, определяемая как отношение радиуса наибольшей вписанной в тело сферы к минимальному радиусу сферы, описанной вокруг тела; эллиптичность, вычисляемая как отношение оередненного радиуса кривизны поверхности к радиусу максимальной сферы, вписанной в тело; гладкость/шероховатость поверхности, которая определяется совокупностью неровностей на базовой длине.
При построении агломерата сфер микрорельеф поверхности формируется за счет множества непересекающихся частей сфер. По аналогии с вычислением объема для оценки гладкости/шероховатости поверхности агломерата предлагается следующий алгоритм.
Определим декартову сетку V с шагом распределения Н = к • Н, к € (0,1] (в рассмотренных
P
VP
h « k • h
V = 0, Pc = 0, ni = 0, П2 = 0 for z = 1; z < dimz; z++ do for y =1; y < dimy; y++ do for x = 1; x < dimx; x++ do
Вычислить p = (hx + p0[1], hy + po[2], hz + _Po[3]) out = true for all < e P do
if ||pv - p|| < rv then out = false break end if end for out
П2 = П2 + 1
else
ni = ni + 1
pPc = pPc + pP
end if end for end for end for
P = Pdim - P0 1/ = У1ЩЗ]
Pc = Pc/ni
далее в примерах принималось к = 0.5). На основе V построим поверхностную четырехугольную сетку Т для линии уровня 1а согласно алгоритму [8]. Гладкость поверхности агломерата Qr определим отношением числа узлов N поверхностной сетки, которые не по пали внутрь Р к общему числу узлов N содержащихся в множестве вершин N сетки Т.
На рисунке 3, а) серым изображена линия уровня 1а/ для а = 0, совпадающая с границей аппроксимируемой модели куба. В таком случае 1а/ целиком включает в себя агломерат Р, обозначенный более темным цветом, и пересекается с ним только в отдельных точках. С увеличением а до а = 1.5 линия уровня 1ап продвигается «внутрь» объекта, что приводит к появлению областей пересечений, как показано на рисунке 3, б). Построенная четырехугольная сетка Т ДЛЯ 1а'' изображена на рисунке 3, в). Вершины сетки, которые не попали в области пересечений Р и Т (обозначены темным), составляют множество N1. Представим алгоритм оценки
а) линия уровня Та/ б) линия уровня /а» в) сетка T да я /а«
Рис. 3. Мера гладкости поверхности агломерата, качества поверхности агломерата сфер так, как показано в алгоритме 3:
Алгоритм 3 Вычисление гладкости поверхности агломерата P
Require: Сетка векселей V, агломерат сфер P
1: h = k • h,N = 0
2: Построить сетку вокселей V с шагом h
3: На основе V построить се тку Т для лини и ур овня 1а, соответствующей значению а [8]
4: for all p eNci do
5: out = true
6: for all ^ eP do
7: if - p| < rv then
8: out = false
9: break
10: end if
11: end for out
13: N = IN +1
14: end if
15: end for
16: Qr = N/N
Оценки гладкости Qr для агломератов куба и позвонка представлены в таблице 1. Как видно из таблицы, для фиксированного значения а параметр Qr уменьшается с возрастанием количества сфер в агломерате и уменьшением R В случае увеличения а значение гладкости улучшается и с некоторого момента становится равным идеальному значению Qr = 0, однако при уменьшении а наблюдается увеличение Qr. Значение Qr ~ 1 достигается при зна-а=0
увеличения точности работы алгоритма возможно увеличение множества вершин N С Т, что достигается перестроением четырехугольной сетки, а также добавлением произвольных точек в ее ячейках.
Полученные агломераты использовались в вычислительных экспериментах дискретно-элементного моделирования динамики соударения твердых тел.
Список литературы
1. Cundall P., Starck О. A discrete numerical model for granular assemblies // Geotechnique. 1979. Vol. 29. .V» 1. P. 47-65. DOI: 10.1680/geot.l979.29.1.47
2. Luding S. Introduction to discrete element methods: basic of contact force models and how to perform the micro-macro transition to continuum theory // European Journal of Environmental and Civil Engineering. 2008. Vol. 12. №7-8. P. 785-826. DOI: 10.1080/19648189.2008.9693050
3. Price M., Murariu V., Morrison G. Sphere clump generation and trajectory comparison for real particles // Proceedings of 5th International Conference on Discrete Element Methods. LCPC. Paris, 2007. P. 239-246.
4. Ferellec J.-F., McDowell G. R. A method to model realistic particle shape and inertia in DEM // Granular Matter. 2010. Vol. 12. № 5. P. 459-467. DOI: 10.1007/sl0035-010-0205-8
5. Garcia X., Latham J .-P., Xiang J., Harrison J.P. A clustered overlapping sphere algorithm to represent real particles in discrete element modelling // Geotechnique. 2009. Vol. 59. № 9. P. 779-784.
DOI: 10.1680/ geot.8.T.037
6. Gao R., Du X., Zeng Y., Li Y., Yan J. A new method to simulate irregular particles by discrete element method // Journal of Rock Mechanics and Geotechnical Engineering. 2012. Vol. 4. № 3. P. 276-281.
7. Amberger S., Friedl M., Goniva C., Pirker S., Kloss C. Approximation of objects by spheres for multisphere simulations in DEM // European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012). Vienna, 2012. P. 5446-5458.
8. Караваев А.С., Копысов С.П. Метод построения неструктурированных шестигранных сеток из объемных данных // Компьютерные исследования и моделирование. 2013. Вып. 5. С. 11-24.
9. Cho G.-C., Dodds J., Santamarina J.С. Particle Shape Effects on Packing Density, Stiffness, and Strength: Natural and Crushed Sands // Journal of Geotechnical and Geoenvironmental Engineering. 2006. Vol. 132. № 5. P. 591-602.
Поступила в редакцию 30.09.2015
Караваев Александр Сергеевич, аспирант, Институт механики УрО РАН, 426067, Россия, г. Ижевск,
ул. Т. Барамзиной, 34.
E-mail: [email protected]
Копысов Сергей Петрович, д. ф.-м. н., заведующий лабораторией, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: [email protected]
Сармакеева Анастасия Семеновна, младший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: [email protected]
A. S. Karavaev, S.P. Kopysov, A. S. Sarmakeeva Construction of clumps for discrete element method
Keywords: discrete element method, clump, error estimation.
MSC: 65M55, 65M75
In this article an algorithm of clumps construction under approximation of the arbitrary shaped bodies geometry in the discrete element method is reviewed. To estimate the surface of arbitrary bodies a measure characterizing the surface microrelief of the formed clump is introduced.
REFERENCES
1. Cundall P., Strack O. A discrete numerical model for granular assemblies, Geotechnique, 1979, vol. 29, no. 1, pp. 47-65. DOI: 10.1680/geot.1979.29.1.47
2. Luding S. Introduction to discrete element methods: basic of contact force models and how to perform the micro-macro transition to continuum theory, European Journal of Environmental and Civil Engineering, 2008, vol. 12, no. 7-8, pp. 785-826. DOI: 10.1080/19648189.2008.9693050
3. Price M., Murariu V., Morrison G. Sphere clump generation and trajectory comparison for real particles, Proceedings of 5th International Conference on Discrete Element Methods, LCPC, Paris, 2007, pp. 239-246.
4. Ferellec J.-F., McDowell G.R. A method to model realistic particle shape and inertia in DEM, Granular-Matter, 2010, vol. 12, no. 5, pp. 459-467. DOI: 10.1007/s10035-010-0205-8
5. Garcia X., Latham J.-P., Xiang J., Harrison J.P. A clustered overlapping sphere algorithm to represent real particles in discrete element modelling, Geotechnique, 2009, vol. 59, no. 9, pp. 779-784.
DOI: 10.1680/geot.8.T.037
6. Gao R., Du X., Zeng Y., Li Y., Yan J. A new method to simulate irregular particles by discrete element method, Journal of Rock Mechanics and Geotechnical Engineering, 2012, vol. 4, no. 3, pp. 276-281.
7. Amberger S., Friedl M., Goniva C., Pirker S., Kloss C. Approximation of objects by spheres for mul-tisphere simulations in DEM, European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012), Vienna, 2012, pp. 5446-5458.
8. Karavaev A.S., Kopysov S.P. The method of constructing unstructured hex mesh for volumetric data, Komp'yuternye issledovaniya i modelirovanie, 2013, vol. 5, pp. 11-24 (in Russian).
9. Cho G.-C., Dodds J., Santamarina J.C. Particle shape effects on packing density, stiffness, and strength: natural and crushed sands, Journal of Geotechnical and Geoenvironmental Engineering, 2006, vol. 132, no. 5, pp. 591-602.
Received 30.09.2015
Karavaev Aleksandr Sergeevich, Post-Graduate Student, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: [email protected]
Kopysov Sergei Petrovich, Doctor of Physics and Mathematics, Head of Laboratory, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: [email protected]
Sarmakeeva Anastasiya Semenovna, Junior Researcher, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: [email protected]
