CC BY
26
Шубин С.Н., Фрейдин А.Б. Алгоритм построения пространственных стохастических микроструктур дисперсно-упрочненных композитов с включениями эллипсоидальной формы // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2016. - № 4. - С. 317-337. DOI: 10.15593/perm.mech/2016.4.19
Shubin S.N., Freidin A.B. A new algorithm for generating a random packing of ellipsoidal inclusions to construct composite microstructure. PNRPU Mechanics Bulletin. 2016. No. 4. Рр. 317-337. DOI: 10.15593/perm.mech/2016.4.19
ВЕСТНИК ПНИПУ. МЕХАНИКА
№ 4,2016 PNRPU MECHANICS BULLETIN
http://vestnik.pstu.ru/mechanics/about/inf/
001 10.15593/регш.шесЬ/2016.4.19 УДК 620.22
АЛГОРИТМ ПОСТРОЕНИЯ ПРОСТРАНСТВЕННЫХ СТОХАСТИЧЕСКИХ МИКРОСТРУКТУР ДИСПЕРСНО-УПРОЧНЕННЫХ КОМПОЗИТОВ С ВКЛЮЧЕНИЯМИ ЭЛЛИПСОИДАЛЬНОЙ ФОРМЫ
С.Н. Шубин1, А.Б. Фрейдин1'2
1 Институт проблем машиноведения Российской академии наук, Санкт-Петербург, Россия 2Санкт-Петербургский политехнический университет Петра Великого, Санкт-Петербург, Россия
О СТАТЬЕ АННОТАЦИЯ
Статья посвящена построению пространственных микроструктур дисперсно-упрочненных композитов. Особенностью таких материалов является наличие связной компоненты (матрицы) и равномерно распределенных по объему частиц наполнителя (включений), которые могут иметь разную форму и размер. В работе рассматриваются включения эллипсоидальной формы, при этом их распределение по объему композита и ориентация считаются стохастическими. Излагается новый алгоритм построения таких микроструктур, основанный на методе, идея которого заключается в итерационном процессе, минимизирующем области пересечения частиц. Вначале случайным образом задаются координаты центров включений и их ориентация, т.е. на этом шаге частицы могут пересекаться. В процессе работы алгоритма пересечения частиц также допускается, при этом на каждом шаге области пересечения минимизируются с помощью задания нового положения включений, которое определяется следующим образом: если два включения пересекаются, то выбирается реперная точка внутри области пересечения и затем новое положение частиц задается таким образом, что реперная точка становится точкой касания. Для того чтобы определить взаимное расположение двух частиц (ни одной общей точки, одна общая точка, множество общих точек) и в случае пересечения выбрать репер-ную точку, используется метод, основанный на решении обобщенной проблемы поиска собственных значений, сформулированной для двух эллипсоидов. Разработанный алгоритм позволяет сформировать микроструктуры с объемными долями включений, близкими к максимально достижимым, для эллипсоидов произвольной формы, в том числе для дисков и вытянутых волокон с большим (более 10) отношением полуосей. Анализ матрицы ориентации показывает, что распределение ориен-таций включений в сформированных микроструктурах изотропно.
© ПНИПУ
Получена: 01 ноября 2016 г. Принята: 12 декабря 2016 г. Опубликована: 30 декабря 2016 г.
Ключевые слова:
композит, стохастическая
микроструктура,
представительный
элемент объема,
эллипсоидальное
включение, матрица
ориентаций.
© Шубин Сергей Николаевич - аспирант, e-mail: [email protected]
Фрейдин Александр Борисович - доктор физико-математических наук, зав. отделом, e-mail: [email protected] Sergey N. Shubin - Postgraduate Student, e-mail: [email protected]
Alexander B. Freidin - Doctor of Physical and Mathematical Sciences, Head of Department, e-mail: [email protected]
A NEW ALGORITHM FOR GENERATING A RANDOM PACKING OF ELLIPSOIDAL INCLUSIONS TO CONSTRUCT COMPOSITE MICROSTRUCTURE
S.N. Shubin1, A.B. Freidin1,2
11nstitute for Problems in Mechanical Engineering of Russian Academy of Sciences, Saint-Petersburg, Russian Federation
2Peter the Great St. Petersburg Polytechnic University, Saint-Petersburg, Russian Federation
ABSTRACT
The subject of the work is a microstructure of a composite which consists of a continuum matrix and a set of isolated particles homogeneously distributed inside the matrix. It is assumed that the reinforcing particles have ellipsoidal shapes, while distribution and orientation are random. The main point of the work is a new computationally-efficient algorithm to generate microstructure of such a composite. In the algorithm the existing "concurrent" method based on an overlap elimination is extended to ellipsoidal shapes of the particles. It begins with randomly distributed and randomly oriented ellipsoidal particles which can overlap each other. During the performance of the algorithm intersections between particles are allowed and at each step the volumes of intersections are minimized by moving the particles. The movement is defined for each pair of particles based on the volume of the intersection: if two particles are overlapped, then the reference point inside the intersection is chosen and then two particles are moved in such a way that the reference point becomes the tangent point for both particles. To define the relative configuration of two particles (separate, tangent or overlapping) and to choose reference point inside the intersection volume the technique based on formulating the problem in four dimensions and then analyzing the roots of the characteristic equation are applied. The algorithm is able to generate close packed microstructures containing arbitrary ellipsoids including prolate and oblate ellipsoids with high aspect ratios (more than 10). The generated packings have a uniform distribution of orientations.
©PNRPU
Введение
Синтез представительного элемента объема (ПЭО) композиционного материала играет ключевую роль при использовании численных методов моделирования в задачах определения эффективных свойств материала, физических полей на микроуровне, несущей способности адгезионного слоя, условий зарождения микротрещин и многих других. Очевидно, что для таких задач принципиально важно смоделировать микроструктуру в наибольшей степени приближенной к реальному внутреннему строению материала.
Данная работа посвящена алгоритму построения пространственных микроструктур дисперсно-упрочненных композитов. Такие материалы характеризуются наличием двух компонентов: связующей однородной матрицы и равномерно распределенных в ней частиц наполнителя (включений). При этом распределение включений в объеме матрицы, а также их форма и размер считаются стохастическими. Формирование внутренней структуры таких композитов связано со случайным расположением непересекающихся включений различной геометрической формы в некотором объеме. В качестве модели реальных геометрических форм включений в данной работе используются эллипсоиды, в предельных случаях переходящих в форму шаров (три полуоси эллипсоида равны друг другу), дисков (одна из полуосей много меньше двух других) и вытянутых волокон (одна из полуосей много больше двух других).
В последние десятилетия активно развиваются методы компьютерного построения микроструктуры, которые условно можно разделить на «последовательные» и «парал-
ARTICLE INFO
Received: 01 November 2016 Accepted: 12 December 2016 Published: 30 December 2016
Keywords:
composite, stochastic microstructure, representative volume element, ellipsoidal inclusion, matrix of orientations.
лельные». Работа последовательных методов основывается на поштучном добавлении включений в структуру ПЭО, в то время как в параллельных алгоритмах реализуется некоторая процедура распределения заданного количества включений в случайную конфигурацию. Критериями оценки работы алгоритмов могут выступать максимально достижимая объемная доля включений и, в случае нешаровой формы включений, равномерность распределения ориентаций включений.
Максимально достижимая объемная доля зависит от формы и размеров включений. Имеется ряд экспериментальных работ по определению этой величины. Наибольшее количество работ посвящено включениям шаровой формы. Максимально плотная упаковка шаров равного радиуса составляет -\/18 « 0,7405 и достигается при образовании регулярной гексагональной укладки (гипотеза Кеплера, сформулированная в 1611 году [1] и доказанная в 1998 [2]). Для нерегулярной укладки шаров максимальная доля точно не определена, как и само это понятие (см. обсуждение этого вопроса в работе [3]). Обычно эту величину оценивают в 0,64 по серии экспериментов, проведенных в работах [4, 5]. В случае эллипсоидальной формы включений максимальная объемная доля при нерегулярной укладке быстро убывает с увеличением отношения полуосей [6-13]. Так, для вытянутых волокон существует экспериментальная оценка максимально достижимой объемной доли [9, 12]:
ра = 5,3,
где р - объемная доля включений; а - отношение длины волокон к диаметру. Подобная зависимость в логарифмической форме представлена в работе [11]. Для эллипсоидов дискообразной формы максимально достижимая объемная доля также убывает с увеличением отношения полуосей [13-17]. Отметим, что для общего случая эллипсоидов произвольной формы и произвольного распределения по размерам предельная объемная доля на сегодняшний день не определена. Для оценки этой величины могут быть использованы методы компьютерного построения микроструктуры.
Среди методов компьютерного построения микроструктуры в литературе наиболее часто встречается метод последовательного типа, впервые представленный в работе [18]. Его суть заключается в следующем. Положение первого включения задается случайным образом. Случайным же образом задаются положения следующих включений, вносимых в структуру последовательно, однако, при этом проверяется пересечения каждого добавляемого включения с ранее размещенными, и если пересечения отсутствуют, то включение добавляется в структуру ПЭО. В противном случае выбирается новое положение до тех пор, пока не найдется такое, в котором пересечения отсутствуют. Процесс добавления включений в структуру ПЭО продолжается до достижения необходимой объемной доли. Данный алгоритм для включений шаровой формы равного радиуса позволяет достичь 0,385 объемной доли [19]. Для вытянутых волокон эллипсоидальной формы качественно зависимость объемной доли от отношения полуосей повторяет экспериментальные оценки, но достигаемая максимальная объемная доля значительно ниже [19].
Принципиально другая идея построения ПЭО предложена в работе [20] для включений круглой формы (плоская задача). Алгоритм основан на методе молекулярной динамики и начинается со стохастического задания векторов, определяющих положения и скорости заданного количества частиц нулевого объема (точек). Новые положения частиц определяются путем интегрирования уравнений движения, при этом размеры частиц увеличиваются с течением времени и отслеживаются два типа событий: столкновение
частиц друг с другом и столкновение частиц с границами ПЭО. Процесс интегрирования по времени итерационный, а шаги интегрирования определяются временем до ближайшего события. Процесс движения частиц в ПЭО продолжается до достижения необходимой объемной доли. Для включений шаровой формы определить время столкновения двух движущихся включений можно из простых аналитических формул, однако ситуация принципиально меняется, когда включения имеют эллипсоидальную форму. В этом случае для определения времени столкновения приходится применять численные методы, что существенно замедляет работу алгоритма. Так, в работе [21, 22] время столкновения определяется из решения оптимизационной задачи по поиску экстремума нелинейной функции, а в работе [23] - из решения нелинейного уравнения с переменными коэффициентами. Таким образом, реализация данного алгоритма, относящегося к параллельному типу, требует больших вычислительных ресурсов, однако позволяет сформировать ПЭО с объемной долей, близкой к предельным значениям. В работе [24] демонстрируется формирование ПЭО для шаровой формы включений равного радиуса с объемной долей 0,637, а в [23] приводятся примеры ПЭО для эллипсоидальной формы включений со значениями объёмных долей, близких к экспериментальным оценкам.
К параллельному же типу относится алгоритм, представленный в работе [25] для включений шаровой формы. Как и в предыдущем алгоритме, сначала случайным образом задаются векторы положения центров включений, при этом включения имеют заданный размер, определяющий объемную долю. Случайное расположение приводит к тому, что в начале работы алгоритма включения пересекаются (если пересечений нет, то ПЭО сформирован и работа алгоритма завершена). Затем запускается итерационный процесс, в котором для пересекающихся включений задаются векторы перемещений, минимизирующих области пересечения. Процесс продолжается до тех пор, пока не останется пересечений. Для включений шаровой формы удается достичь объемной доли в 0,64 [25]. В работе [26] данный алгоритм расширен для случая включений цилиндрической формы, и достигаемые объемные доли близки к экспериментальным оценкам. Для включений эллипсоидальной формы идея алгоритма не реализована.
Выше перечислены алгоритмы, наиболее часто встречающиеся в литературе. Не вдаваясь в детали, приведем ссылки на еще два оригинальных, но сложных в исполнении и требующих значительных вычислительных ресурсов алгоритма, представленных в работах [27] и [28]. Отметим, что на сегодняшний день только один алгоритм, основанный на методе молекулярной динамики, позволяет сформировать ПЭО для эллипсоидальных включений с предельными значениями объемных долей, однако его реализация является вычислительно затратной.
Основная цель данной работы - расширить алгоритм, разработанный для включений шаровой и цилиндрической форм, минимизирующий области пересечения [25, 26], на случай включений эллипсоидальной формы. Идея, заложенная в предлагаемом методе, заключается в следующем. Сначала случайным образом задаются положения и ориентации включений заданного размера. Затем осуществляется поиск пар пересекающихся включений. Для каждой такой пары находится точка внутри области пересечения, и затем каждое из включений перемещается таким образом, что эта точка становится точкой касания. В разработанном алгоритме учтена возможность использования сформированных микроструктуры в задачах численной гомогенизации, для которых предпочтительнее использовать периодические ПЭО [29]. В этом случае распределение включений внутри ПЭО явля-
ется стохастическим, однако для каждого включения, пересекающего границу ПЭО, есть периодический образ на противоположной грани ПЭО. Данная опция является вспомогательной и может быть с легкостью исключена из алгоритма в случае необходимости.
Статья построена так, чтобы на ее основе можно было реализовать предлагаемый алгоритм, поэтому все шаги описаны достаточно детально. В следующем за введением разделе описывается математический аппарат, необходимый для геометрического моделирования включений в ПЭО. Третий раздел посвящен непосредственно реализации алгоритма. Разбираются основные шаги алгоритма, и делаются замечания по их численной реализации. В четвертом разделе представлены результаты работы алгоритма: приводятся примеры ПЭО с различными формами включений, оценивается равномерность распределения включений по ориентациям, и проводится сравнение с другими алгоритмами.
1. Включения в представительном элементе объема
Данный раздел является вспомогательным и служит для сохранения целостности изложения. Сведения, приведенные здесь, общеизвестны, однако геометрическое моделирование вносит определенную специфику в описание местоположения и взаимного расположения включений в ПЭО, поэтому для понимания основной части работы мы достаточно детально останавливаемся на данном вопросе.
Всюду в статье прописными латинскими буквами обозначены скаляры, прописными жирными буквами - векторы, заглавными жирными буквами - матрицы.
1.1. Положение включения в представительном элементе объема
В статье рассматривается ПЭО в форме куба, стороны Ь внутри которого случайным образом расположены включения эллипсоидальной формы (рис. 1). Введем ортогональную систему координат Ое1е2ез так, чтобы точка О совпадала с одной из вершин куба, а орты направим вдоль сторон, образующих эту вершину. Рассмотрим эллипсоид, расположенный произвольным образом относительно Ое1е2е3. Обозначим через О' центр эллипсоида и введем локальную систему координат О'е'1е'2е'3, оси которой направлены вдоль главных осей эллипсоида. Тогда уравнение эллипсоида в локальной системе координат примет канонический вид
х'г Ах' = 0,
0 0 0 ^ 1/Ь2 0 0 ; 0 1/с2 0 ; 0 0 -1 ,
а, Ь, с - полуоси эллипсоида. Здесь для описания местоположения точек в пространстве использованы однородные координаты, т. е. вектор х определяется заданием четырех чисел х = (хх, х2, х3, м?)г. Напомним, что х = (хх/ м>, х2/ м>, х3/ м?)г есть вектор соответствующей точки в декартовой системе координат.
где
А' =
^1/ а2 0 0 0
х'\
WO'
/ r
Введение однородных координат продиктовано эффективностью работы для выражения преобразований точек при численном моделировании, так как в этом случае вращение и смещение может быть описано с помощью матричного умножения. Так, переход от локальных координат к глобальным осуществляется матрицей М:
х = Мх',
где
O
Рис. 1. Эллипсоидальное включение в представительном элементе объема Fig. 1. An ellipsoidal inclusion in the representative volume element
M
Г R 0T
Я - матрица поворота; г - радиус-вектор точки 0'. Тогда в глобальной системе координат уравнение эллипсоида примет вид
где
xT Ax = 0,
A = M A' M"
здесь М т обозначает обратную матрицу от Мт. В силу того, что Я ортогональная матрица, обратная матрица М-1 может быть вычислена по формуле
М 1
Г RT
RT r ^
0T
1
Матрица Я определяет поворот на угол 9 вокруг единичного вектора т. Запишем эту матрицу через кватернион, который определяется с помощью задания скаляра а и вектора V по следующему правилу:
( ) Г 9 . 9Л
q = (а, V) = 1 соб-^,т вт-^ I.
Заметим, что так как |т| = 1, то = ^а2 + = 1, т. е. кватернион нормирован. Тогда матрицу поворота можно представить в следующем виде:
Я = (2а2 -1)Е + 2vvт + 2а8, где Е - единичная матрица, а матрица 8 выражается через компоненты вектора V:
Л
Г 0 "V3 V2
S = V3 0
v"V2 V1 0
где - 1-й компонент вектора V. Напомним также правило умножения двух кватернионов, т.е. правило сложения поворотов. Пусть ориентация эллипсоида задана кватернионом
X
2
X
q1 = (а1, у1), и пусть необходимо совершить поворот этого эллипсоида кватернионом q2 = (а2, у2), тогда кватернион, определяющий новую ориентацию, вычисляется по формуле
ql2 = (аШ У12) = (а1а2 - У1Гv2, а1у2 + а2У1 - V Х У2)-1.2. Взаимное расположение двух включений
Пусть имеются два эллипсоида хт А'х = 0 и хтА1 х = 0 . Составим характеристическое уравнение
ёе1;(ЯА' - А1) = р^4 + р2Х3 + р3Х2 + р4Х + р5 = 0. (1)
Это полином четвертой степени, и поэтому всегда имеется четыре корня. В работе [30] показывается, что определить взаимное расположение двух эллипсоидов можно исходя из анализа корней характеристического уравнения (1). Сначала отметим, что данное уравнение в случае эллипсоидов всегда имеет два положительных корня, поэтому анализ основывается на оставшихся двух корнях. Возможны следующие варианты:
а) два некратных отрицательных корня: А' и А1 не имеют общих точек;
б) два кратных отрицательных корня: А' и А1 касаются друг друга с внешней стороны обоих эллипсоидов (т.е. имеют только одну общую точку);
в) два комплексно сопряженных корня: А' и А1 имеют общую область;
г) два кратных положительных корня: А' и А1 имеют общую область, и их границы имеют одну или две точки касания;
д) два некратных положительных корня: А' и А1 пересекаются или один из эллипсоидов лежит внутри другого.
Других возможных комбинаций корней нет.
Вычислить коэффициенты характеристического полинома можно по формулам, приведенным в работе [31]:
Р1 = -ЗА^; Р2 = 52§3СП + 5153С22 + 5152С33 -515253С44;
Р3 = 5152 (С33С44 - С34С43 ) ^ 5253 (С11С44 - С14С41 ) ^ 5153 (С22С44 - С24С42 ) ^ +51(С23С32 -С22С33) ^52(С13С31 -С11С33) ^53(С12С21 -С11С22);
р4 = 51 (-С22С33С44 ^ С22С34С43 ^ С33С42С24 ^ С44С32С23 - С32С24С43 - С42С23С34) +§2(-СПС'33С44 + С11С34С43 + С33С14С41 + С44С13С31 - С31С14С43 - С41С13С34 ) +53(-С11С22С44 ^ С11С24С42 ^ С22С14С41 ^ С44С12С21 - С21С14С42 - С41С12С24 ) +С„С22С33 -С11С23С32 -С22С13С31 -С33С12С21 + С21С13С32 + С31С12С23;
р5 = ёе1(А1),
где
51 = 1/а'2, 52 = 1/V2, 53 = 1/ С2, С = М'тА3М'.
Таким образом, взаимное расположение двух эллипсоидов сводится к определению корней полинома четвертой степени, что численно легко реализуется и достигается выполнением счетных алгебраических действий.
Отметим еще одно важное свойство задачи (1) на обобщенные собственные значения. В работе [19] показывается, что: в случае б) и г) соответствующий(ие) собствен-ный(ые) вектор(ы) указывает(ют) на точку(и) касания; в случае в) вещественная часть соответствующих собственных векторов (собственные векторы комплексно сопряженные) указывает на точку внутри области пересечения; в случае д) один из соответствующих собственных векторов указывает на точку внутри области пересечения. Данное свойство будет использовано ниже.
1.3. Взаимное расположение включения и границ представительного элемента объема
Пусть в глобальной системе координат уравнение эллипсоида имеет вид
хт Ах = 0.
Матрица А симметричная и может быть представле в виде
А =
А 1 Л в -а 2
1 т
V 2
где В - матрица 3 х 3 ; а - вектор 3 х1; F - скаляр. Тогда в декартовой системе координат уравнение эллипсоида запишется в виде
хт Вх + ат х + F = 0, (2)
здесь уже - декартовы координаты. Отметим, что в симметричная положительно определенная матрица, так как последнее уравнение есть уравнение эллипсоида. Границы ПЭО определяются уравнениями
х, = р, к = 1,2,3, р = 0, Ь.
То есть, например, для первой переменной имеем две плоскости х1 = 0 и х1 = Ь . Для того чтобы установить, пересекает ли эллипсоид плоскость хк = Р, введем в рассмотрение вектор на этой плоскости
х = (х,, хт )т, I < т, IФ к Ф т.
Тогда точка на плоскости может быть представлена в виде
х = Рх + р, (3)
где Р - матрица 3 х 2, составленная из двух ортов осей I и т,
Р = (е,, ет),
а р - вектор, натянутый на орт ек,
Р = рек.
Так, например, для плоскости х1 = Ь уравнение (3) примет вид
V х3
( 0 0 Л ( х1 л ( Ь Л
1 0 1 + 0
V 0 1, V х2 ) V 0 ,
Подставив (3) в (2), получим
хт в х+ат х+р = 0, (4)
где
В = Рт ВР,
а = 2Рт вр+рт а,
Р = рт вр+ат р+р.
Заметим, что матрица В получается из В путем вычеркивания к строки и к столбца. А так как матрица В положительно определенная, то и В тоже положительно определенная.
Сделаем в (3) замену переменных х = х * -х0. Имеем
х *т Вх *+(2хтВ - ат )х *+Р+х0Вх0 - атх0 = 0.
Выберем х 0 так, чтобы выражение в скобках равнялась нулю:
х„ =
1В -а.
(5)
Тогда уравнение (5) примет вид
х *т Вх * +§ = 0,
где скаляр £ определен выражением
1
(6)
£ = р—ат в-1а.
4
Итак, вид пересечения эллипсоида с плоскостью хк = Р определяется уравнением (6),
а так как В положительно определенная матрица, то возможны три случая:
1) если £ < 0, эллипсоид пересекает плоскость и пересечение есть эллипс;
2) если £ = 0, эллипсоид пересекает плоскость в одной точке;
3) если £ > 0, эллипсоид не пересекает плоскость.
2. Алгоритм формирования представительного элемента объема
Опишем основные шаги алгоритма. Сначала случайным образом задается местоположение центров эллипсоидов в кубе стороной Ь. Ориентация эллипсоидов также задается случайным образом. Размеры эллипсоидов выбираются исходя из объёмной доли, которую должен иметь ПЭО. На каждом шаге алгоритма проверяется два типа пересечения: пересечение
эллипсоидов с границами куба и взаимные пересечения эллипсоидов. Для первого типа пересечения создаются периодические образы эллипсоидов с противоположной стороны, чтобы сохранялась периодичность ПЭО. Во втором случае эллипсоиды раздвигаются так, чтобы минимизировать область пересечения. Работа алгоритма завершается при отсутствии взаимного пересечения эллипсоидов. Рассмотрим более детально каждый из типов пересечения.
2.1. Пересечение включением границ представительного элемента объема
Пусть эллипсоид пересекает границу ПЭО. В этом случае необходимо создать его периодический образ с противоположной стороны. Расстояние между центрами исходного эллипсоида и его периодического образа задается вектором
I Р = 0,
(7)
-I, р = 1.
Отметим, что периодический образ имеет такую же пространственную ориентацию, что и исходный эллипсоид. Ясно, что количество периодических образов эллипсоида зависит от того, сколько границ области он пересекает. Если эллипсоид пересекает две плоскости хк = Рк и х1 = Рг, то помимо двух периодических образов, заданных выражением (7), необходимо создать третий, центр которого определяется вектором
^ К . (8)
Наконец, если эллипсоид пересекает три плоскости хк = Рк, х1 = Рг и хт = Рт, то периодических образов будет семь. Местоположение центров шести из них определены векторами (7) и (8), а седьмой - вектором
Щк К )е, +Щт К. 2.2. Взаимное пересечение включений
Пусть два включения А' и А1 пересекаются. Идея предлагаемого метода заключается в поиске нового положения включений, в котором бы они имели только одну общую точку. Поиск такого положения предлагается вести итерационным способом, постепенно уменьшая область пересечения. Для этого на каждом шаге внутри области пересечения выбирается реперная точка. Далее задается вектор перемещения точек эллипсоидов, совпадающих с реперной точкой таким образом, чтобы реперная точка оказалась на гранях каждого из эллипсоидов и одновременно приближенно (степень приближенности будет определена ниже) их точкой касания. Новое положение эллипсоидов определяется с помощью поворота и перемещения эллипсоидов так, чтобы точки эллипсоидов, совпадающих с реперной точкой, переместились на заданные векторы.
2.2.1. Выбор реперной точки внутри области пересечения
Отметим, что выбор реперной точки не определен однозначно. Ясно, что удачный выбор может ускорить последующую работу алгоритма. С другой стороны, не стоит забывать о времени, потраченном на поиск «оптимальной» точки. Например, в качестве та-
кой точки можно выбрать центр масс области пересечения. Однако поиск центра масс может занять достаточно много машинного времени, что сильно замедлит работу алгоритма. Поэтому предлагается в качестве реперной выбрать точку, определяемую собственными векторами для задачи на обобщенные собственные значения приведенной в п. 2.3, так как на этом шаге алгоритма собственные значения уже известны и вычисление соответствующих собственных векторов не требует больших вычислительных затрат.
Как уже было сказано выше, для пересекающихся эллипсоидов возможны три случая корней характеристического полинома. Если корни комплексно сопряженные, то и собственные векторы будут также комплексно сопряженными. В этом случае вещественная часть собственных векторов указывает на точку внутри области пересечения. Если имеется два некратных положительных корня, то один из собственных векторов указывает на точку внутри области пересечения. Наконец, для двух кратных положительных корней собственный(е) вектор(ы) будет(ут) указывать на точку(и) касания. В этом случае выбор точки внутри области пересечения предлагается провести следующим образом. Введем в рассмотрение прямую, проходящую через точку касания хс (в случае двух точек выбирается произвольным образом одна из них) и параллельную вектору нормали п к поверхностям эллипсоидов в этой точке:
х = хс + пг, г е Я.
Для определения вектора нормали п продифференцируем уравнение одного из эллипсоидов:
п = -V(хтвх + атх + р)| = -2Вхс - а.
1х=хс
Мы выбрали знак минус, так как нас будет интересовать нормаль внутрь области. Найдем точки пересечения прямой с каждым из эллипсоидов. Для этого подставим уравнение прямой в уравнения для каждого из эллипсоидов (2):
(хс + пг)т вк (хс + пг)+акт (хс + пг)+Рк = 0, к = I, ].
Раскрывая скобки и приводя подобные, получим два квадратных уравнения
а кг2 + ркг + ук = 0,
где
ак = ет Вк е, рк = 2ет Вке + акт е, ук = хстВкхс + акт хс + Рк.
Каждое из двух уравнений имеет нулевой и положительный корни. Обозначим положительные корни через г1 и г1 соответственно для эллипсоидов I и 1 . Положение репер-ной точки определим выражением
хт = хс + п ш1и(гг, г1 )/2.
2.2.2. Определение вектора перемещения
Построение вектора перемещения идентично для обоих эллипсоидов, поэтому рассмотрим только один из них. По идее метода точка эллипсоида, совпадающая с реперной точкой, должна переместиться так, чтобы реперная точка оказалась на поверхности эл-
липсоида, при этом векторы нормали обоих эллипсоидов в реперной точке должны быть коллинеарными. Таким образом, необходимо определить вектор нормали в новом, пока неизвестном, положении эллипсоида. Для этого введем в рассмотрение функцию
/ (х) = -хт Бх-ат х- ^.
Эта функция имеет максимум в середине эллипсоида и в точке максимума равна 1. Она равна нулю на границе эллипсоида и меньше нуля вне его. Найдем градиент этой функции:
V/ (х) = -2Бх- а. (9)
Определим приближенно вектор нормали для нового положения эллипсоида в репер-ной точке хт вдоль градиента функции /(х) :
е= V/ (хт) V/ (хт )|"
На рис. 2 изображены поверхности уровня функции /(х). Ясно, что чем меньше область пересечения, тем ближе поверхность
о. т
уровня, на которой лежит точка х , к поверхности эллипсоида, и тем меньше угол Рис. 2. Линии уровня функции / (х)
между вектором е и нормалью к поверхно-
Рщ. 2. 18о1те8 ог шпсйоп/ (х)
сти эллипсоида в точке ближайшей к хт .
Искомый вектор перемещения сонаправим с вектором е. Величину перемещения определим как ближайшее расстояние от реперной точки до поверхности эллипсоида вдоль прямой, определяемой вектором е . Для этого запишем уравнение прямой в параметрическом виде
х = хт + et, г е Я.
Найдем точки пересечения прямой с эллипсоидом. Для этого подставим уравнение прямой в уравнение эллипсоида (2):
(хт + ег )Т Б(хт + ег) + Ат (хт + ег) + ^ = 0.
Раскрывая скобки и приводя подобные, получим квадратное уравнение
аг2 +рг + у = 0,
где
а = еТ Бе, р = 2ет Бе + Ат е, у = хтТ Бхт + Ат хт + ^.
Из этого уравнения находим два вещественных корня г1 и г2 (будем считать, что г1 < г2). Двум значениям параметра г соответствуют две точки:
х1 = хт + ег1, х2 = хт + ег2.
Очевидно, что так как хт находится внутри эллипсоида, то корни г1 и г2 разных знаков. Так как перемещение должно осуществляется вдоль вектора е, то модуль перемещения будет определяться отрицательным корнем. Величину перемещения точки хт определим как
II // m 1\Т / m 1\ /, T , L I
u = y¡ (x - x ) (x - x ) =yjt1e t1e = í^ . Итак, вектор перемещения точки xm
u = |ii|e. (10)
Такому определению вектора перемещения соответствует простая механическая интерпретация, в которой функция f(x) представляет собой потенциал некоторой силы, «выталкивающей» эллипсоид из области пересечения.
2.2.3. Определение перемещения и поворота эллипсоида
После того как вектор перемещения точки эллипсоида, совпадающей с реперной, найден, необходимо определить новое положение эллипсоида. В предельном случае, когда эллипсоид представляет собой шар, новое положение определяется простым переносом центра шара на вектор u . В общем случае определение нового положения неоднозначно. Например, можно перенести центр эллипсоида или осуществить поворот эллипсоида вокруг некоторой оси. В данной статье предлагается комбинация этих двух возможностей. Для этого представим вектор u в виде суммы двух векторов:
u = ucec + (u - ucec X (11)
где единичный вектор ec указывает направление от точки xm к центру эллипсоида,
г - x m
ec = i , (12)
// m \T / m \
V(r - x ) (г - x )
где г - радиус-вектор, указывающий на центр эллипсоида; uc определяется скалярным произведением векторов u и ec :
uc = eT u.
Новое положение центра эллипсоида определим выражением
= г + ucec (13)
и зададим поворот эллипсоида так, чтобы смещение точки, совпадающей с xm , определилось вектором (u - ucec). Для этого определим вектор, относительно которого совершается поворот векторным произведением:
m = (u - ucec ) Х ec
V(u - ucec )Т (u - ucec )
Отметим, что |m| = 1. Угол поворота определим приближенно:
V(u - ucec)Т (u - ucec)
V(r - xm )Т (г - xm )
Таким образом, поворот определяется кватернионом
q
е . е
cos—, m sin — 2 2
2
На рис. 3 проиллюстрирован итерационный процесс, переводящий пересекающиеся эллипсоиды в касающиеся. На первом шаге выбирается реперная точка внутри области пересечения; вычисляется градиент функции /(х) по формуле (9); с помощью выражения (10) определяются векторы перемещения точек эллипсоидов, совпадающих с реперной
точкой; задается новое положение и пространственная ориентация эллипсоидов, определяемые формулами (13) и (14). На втором шаге реперная точка оказывается расположенной на гранях эллипсоидов, но она не является точкой касания, так как векторы нормалей к двум эллипсоидам в этой точке были определены приближенно, более того, они не были коллинеарны, и как следствие эллипсоиды все еще имеют общую область, однако она меньше, чем на шаге 1. Далее итерационный процесс продолжается с постепенным уменьшением области пересечения, при этом определяемые в реперной точке векторы нормалей постепенно становятся коллинеарными.
3. Результаты
3.1. Примеры упаковок
На рис. 4 и рис. 5 приведены примеры ПЭО для различных форм включений с максимальной объемной долей, полученных с помощью разработанного алгоритма, а в табл. 1 проведено сравнение с результатами двух алгоритмов последовательного [19] и параллельного [23] типов. Из таблицы видно, что, как уже отмечалось во введении, максимальная объемная доля последовательного алгоритма значительно ниже значений, получаемых параллельными алгоритмами. Сравнение с результатами работы [23] показывает, что максимальная объемная доля разработанного алгоритма не меньше, а для некоторых соотношений полуосей больше.
На практике нередко встречаются композитные материалы, в которых размеры включений и соотношение полуосей являются стохастическими величинами. Для оценки форм включений таких материалов используется анализ, основанный на двумерном изображении включений. Выходными данными анализа являются распределения включений
Точка касания
Рис. 3. Итерационный процесс по уменьшению области пересечения Fig. 3. Stages in the elimination of overlaps
а б в
Рис. 4. Представительные элементы объема для включений эллипсоидальной формы с соотношением полуосей меньшим единицы: (а) a = b = с/2, объемная доля 0,66; (б) a = b = с/10, объемная доля 0,33; (в) a = b = с/30, объемная доля 0,10 Fig. 4. Fig. 4. Representative volume elements contain inclusions of prolate ellipsoids: (а) a = b = c/2, volume fraction of 0,66; (b) a = b = c/10, volume fraction of 0,33; (с) a = b = c/30, volume fraction of 0,10
а б в
Рис. 5. Представительные элементы объема для включений эллипсоидальной формы с соотношением полуосей большим единицы: (а) a = b = 2с, объемная доля 0,65; (б) a = b = 10с, объемная доля 0,33; (в) a = b = 30с, объемная доля 0,10 Fig. 5. Representative volume elements contain inclusions of oblate ellipsoids: (а) a = b = 2с, volume fraction of 0,65; (б) a = b = 10с, volume fraction of 0,33; (в) a = b = 30с, volume fraction of 0,10
Таблица 1
Сравнение максимальных объемных долей полученных с помощью разных алгоритмов
Table 1
Comparison of maximum volume fraction obtained with various algorithms
Отношение полуосей эллипсоида a / с Результаты, приведенные в работе[19] Результаты, приведенные в работе [23] Результаты данной работы
1 / 2 0,40 0,60 0,66
1 / 10 0,25 0,30 0,33
1 / 30 - 0,10 0,10
2 0,41 0,60 0,65
10 0,26 0,30 0,33
30 - 0,10 0,10
35
30
25
20
15
и 10
35
x= 30 0х
£ 25 и
<D
j 20
в 15
H X
и
10
0 7,2 14,4 21,6 28,8 Диаметр/Сторона ПЭО, x 10-2
36
0,25 0,4 0,55 0,7 0,85 Отношение полуосей
а б в
Рис. 6. Представительный элемент объема с включениями, имеющими заданное распределение по размерам и отношениям полуосей: а - представительный элемент объема; б - распределение по размерам; в - распределение по отношению полуосей. Объемная доля включений 36 %,
количество включений 300 Fig. 6. Fig. 6. Representative volume element of polydisperse ellipsoids with the particular size and aspect ratio distributions: а - representative volume element; b - size distribution; c - respect of semiaxes distribution. Volume fraction 36 %, number of ellipsoids 300
по осредненному диаметру (диаметр круга с площадью, равной площади проекции включения) и отношению полуосей. На рис. 6, б и в приведен пример таких распределений, а на рис. 6, а изображен соответствующий ПЭО. Таким образом, разработанный алгоритм позволяет сформировать микроструктуру с эллипсоидами случайных размеров и формы.
3.2. Распределение эллипсоидов по ориентации
Одна из важных характеристик ПЭО - это распределение по ориентациям эллипсоидов. Ясно, что для случайного расположения включений распределение должно быть близко к равномерному (изотропному). Для оценки близости функции вероятности распределения ориентации включений к постоянной величине введем в рассмотрение матрицу ориентации, предложенную в работе [32],
ТШпП =(^п) , т П = 1,2,3
где w1, w2 и - три компоненты единичного вектора направленного вдоль одной из осей эллипсоида; означает осреднение по всем эллипсоидам в ПЭО. Отметим, что
матрица ориентации симметрична и след матрицы равен единице. В случае изотропного распределения тензор ориентации равен (1/3) Е, где Е единичная матрица. Для оценки
близости матрицы ориентации к (1/3) Е введем в рассмотрение отношение р = ^тах /^тт, где ^тах и ^тт обозначают максимальное и минимальное собственные значения матрицы ориентации. По близости числа р к единице можно оценить близость распределения ори-ентаций к изотропному. Зависимость числа р от количества включений для двух типов эллипсоидов показана на рис. 7. Для каждого количества включений проведена серия из 10 реализаций алгоритма. На графиках отмечены средние значения числа р с 95 % доверительным интервалом. В обоих случаях объемная доля эллипсоидов равна 0,1.
5
5
0
0
1
Из графиков видно, что р^ 1 при N , при этом сужается и доверительный интервал. Это говорит о том, что разработанный алгоритм формирует действительно стохастический ПЭО.
О 100 200 300 400 0 100 200 300 400
Количество включений Количество включений
а б
Рис. 7. Отношение собственных чисел матрицы ориентации от количества включений.
Каждая точка обозначает среднее значение по 10 реализациям алгоритма: (a) a = b = 0,1с; (б) a = b = 10с. На графиках также отображен 95 % доверительный интервал Fig. 7. Ratio р as a function of the ellipsoids number. Each point represents the mean value obtained for 10 random realizations: (a) a = b = 0,1с; (b) a = b = 10с. There is also 95% confidence interval shown
Приведем еще одну оценку близости распределения ориентаций включений в сформированных ПЭО к изотропному. Напомним, что работа алгоритма начинается со случайного задания месторасположения и ориентации включений. Для того чтобы убедиться, что алгоритм не выстраивает включения в определенном направлении, можно сравнить число р в начале работы алгоритма и в конце. При этом так как вначале ориентация произвольная, то соответствующее значение р можно считать в некотором смысле «эталонным» для данного количества частиц. В табл. 2 приведены значения числа р в начале и в конце работы алгоритма для ПЭО, приведенных на рис. 4 и 5. Из таблицы видно, что выстраивание включений вдоль определенного направления не происходит и числа р входят в коридор, показанный на рис. 7.
Таблица 2
Значение р в начале и в конце работы алгоритма для ПЭО, приведенных на рис. 4 и 5
Table 2
Ratio p at the beginning and end of the algorithm performance for the RVE, shown in fig. 4 and 5
Отношение Количество Значение р в начале Значение р в конце
полуосей a/с эллипсоидов работы алгоритма работы алгоритма
1 / 2 30 2,17 1,53
1 / 10 90 1,44 1,47
1 / 30 175 1,25 1,19
2 30 1,70 2,03
10 60 1,36 1,99
30 90 1,96 1,53
Заключение
В статье изложен новый алгоритм формирования пространственных микроструктур дисперсно-упрочненных композиционных материалов с упрочняющими частицами в форме эллипсоидов. Разработанный алгоритм позволяет сформировать ПЭО:
- с равномерным стохастическим распределением эллипсоидов по ориентациям и по внутреннему объему;
- с заданным распределением эллипсоидов по размерам и отношениям полуосей, т.е. эллипсоиды внутри ПЭО могут иметь неодинаковые размеры;
- с объемной долей не меньших, а в некоторых случаях превосходящих значений, достигаемых известными на сегодняшний день алгоритмами.
Разработанный алгоритм достаточно прост в реализации и не требует больших вычислительных ресурсов. Определение собственных чисел и собственных векторов матриц 4x4 - это наиболее ресурсоемкая операция алгоритма. Однако в вычислительном смысле решение этой задачи представляет собой счетную комбинацию простых алгебраических действий. Эти действия значительно менее ресурсоемки по сравнению с решением нелинейного уравнения с переменными коэффициентами, необходимыми для реализации алгоритма, основанного на методе молекулярной динамики [23] (напомним, что на сегодняшний день из всех известных алгоритмов только этот позволяет достичь объемных долей, близких к получаемым с помощью приведенного в данной работе). Отметим, что предложенный алгоритм допускает касание частиц. Дополнительное искусственное увеличение размеров превращает области касания в переходные слои заданной толщины.
В завершение отметим, что сформированные микроструктуры дисперсно-упрочненных композитов могут быть использованы в широком спектре задач физики и механики композиционных материалов. Так, например, можно исследовать влияния распределений формы (отношения полуосей) и размеров включений на эффективные упругие, электрические и другие свойства композита, на зарождение микродефектов и многое другое.
Библиографический список
1. Kepler J. Strena seu de nive sexangula. - Francofurti ad Moenum: apud Tampach, 1611. - 21 p.
2. Hales T.C. A proof of the Kepler conjecture // Annals of mathematics. - 2005. - Vol. 162. -P.1065-1185.
3. Torquato S., Truskett T.M., Debenedett P.G. Is random close packing of spheres well defined? // Physical Review Letters. - 2000. - Vol. 84 - P. 2064-2067.
4. Scott G.D. Packing of spheres // Nature. - 1960. - Vol. 188 - P. 908-909.
5. Scott G.D., Kilgour D.M. The density of random close packing of spheres // J. Phys. D: Appl. Phys. - 1969. - Vol. 2 - P. 863-866.
6. Milewski J.V. A study of the packing of milled fibreglass and glass beads // Composites. - 1973. -Vol. 4 - P. 258-265.
7. Milewski J.V. The combined packing of rods and spheres in reinforcing plastics // Ind. Engng. Chem. Prod. Res. Dev. - 1978. - Vol. 17. -P. 363-366.
8. Nardin M., Papirer E. Contribution a l'etude des empilements au hasard de fibres et/ou de particules spheriques // Powder Technology - 1985. - Vol. 44. - P. 131-140.
9. Evans K.E., Gibson A.G. Prediction of the maximum packing fraction achievable in randomly oriented short-fibre composites // Composites Science and Technology - 1986. - Vol. 25. - P. 149-162.
10. Evans K.E., Ferrar M.D. The packing of thick fibers // J. Phys. D: Appl. Phys. - 1989. -Vol. 22. - P. 354-360.
11. Parkhouse J.G., Kelly A. The Random Packing of Fibres in Three Dimensions // Proc. R. Soc. Lond. A. - 1995. - Vol. 451. - P. 737-746.
12. Philipse A.P. The random contact equation and its implications for (colloidal) rods in packings, suspensions, and anisotropic powders // Langmuir. - 1996. - Vol. 12. - P. 1127-1133.
13. Effects of grain shape on packing and dilatancy of sheared granular materials / Wegner S. [et al.] // Soft Matter. - 2014. - Vol. 10. - P. 5157-5167.
14. Improving the density of jammed disordered packings using ellipsoids / Donev A. [et al.] // Science. - 2004. - Vol. 303. - P. 990-993.
15. Experiments on random packings of ellipsoids / Man W. [et al.] // Physical Review Letters. -2005. - Vol. 94 (19). - P. 198001-1 - 198001-4.
16. X-ray tomography study of the random packing structure of ellipsoids / Xia C. [et al.] // Soft Matter. - 2014. - Vol. 10. - P. 990-996.
17. Local Origin of Global Contact Numbers in Frictional Ellipsoid Packings / Schaller F.M. [et al.] // Physical Review Letters. - 2015. - Vol. 114(15). - P. 158001-1 - 158001-5.
18. Bennett C. H. Serially depositied amorphous aggregates of hard spheres // J. Appl. Phys. -1972. - Vol 32(6). - P. 2727-2734.
19. Sherwood J. D. Packing of spheroids in three-dimensional space by random sequential addition // J. Phys. A: Math. Gen. - 1997. - Vol. 30. - P. 839-843.
20. Lubaehevsky B.D., Stillinger F.H. Geometric properties of random disk packings // Journal of Statistical Physics. - 1990. - Vol. 60. - P. 561-583.
21. Donev A., Torquato S., Stillinger F.H. Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. I. Algorithmic details // Journal of computational physics. -2005. - Vol. 202 (2) - P. 737-764.
22. Donev A., Torquato S., Stillinger F.H. Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. II Application to ellipses and ellipsoids // Journal of computational physics. - 2005. - Vol. 202 (2). - P. 765-793.
23. Ghossein E., L'evesque M. Random generation of periodic hard ellipsoids based on molecular dynamics: a computationally-efficient algorithm // Journal of Computational Physics. - 2013. - Vol. 253. -P.471-490.
24. Lubachevsky B.D., Stillinger F.H., Pinson E.N. Disks vs. spheres: contrasting properties of random packings // Journal of Statistical Physics. - 1991. - Vol. 64. - P. 501-524.
25. Jodrey W.S. Computer simulation of close random packing of equal spheres // Physical Review A. - 1985. - Vol. 32. - No. 4. - P. 2347-2351.
26. Williams S.R., Philipse A.P. Random packings of spheres and spherocylinders simulated by mechanical contraction // Physical Review E. - 2003. - Vol. 67. - P. 051301-1 - 051301-9.
27. Stillinger F.H., Weber T.A. Inherent structure of liquids in the hard-sphere limit // J. Chem. Phys. - 1985. - Vol. 83(9). - P. 4767-4775.
28. Zinchenko A. Algorithm for random close packing of spheres with periodic boundary conditions // J. Comp. Phys. - 1994. - Vol. 114. - P. 298-307.
29. Determination of the size of the representative volume element for random composites: statistical and numerical approach / Kanit T. [et al.] // International Journal of Solids and Structures. -2003. - Vol. 40. - P. 3647-3679.
30. Chan K. A simple mathematical approach for determining intersection of quadratic surfaces // Multiscale optimization methods and applications. - New York, 2006. - P. 271-298.
31. Continuous collision detection for ellipsoids / Choi Y. [et al.] // Visualization and computer Graphics, IEEE Transactions. - 2009. - Vol. 15. - P. 311-325.
32. Advani S.G., Tucker C.L. The use of tensors to describe and predict fiber orientation in short fiber composites // Journal of rheology. - 1987. - Vol. 31. - P. 751-784.
References
1. Kepler J. Strena seu de nive sexangula. Francofurti ad Moenum, apud Tampach, 1611, 21 p.
2. Hales T.C. A proof ofthe Kepler conjecture, Annals of mathematics, 2005, vol. 162, pp. 1065-1185.
3. Torquato S., Truskett T.M., Debenedett P.G. Is random close packing of spheres well defined? Physical Review Letters, 2000, vol. 84, pp. 2064-2067.
4. Scott G.D. Packing of spheres. Nature, 1960, vol. 188, pp. 908-909.
5. Scott G.D., Kilgour D.M. The density of random close packing of spheres. J. Phys. D: Appl. Phys, 1969, vol. 2, pp. 863-866.
6. Milewski J.V. A study of the packing of milled fibreglass and glass beads. Composites, 1973, vol. 4, pp. 258-265.
7. Milewski J.V. The combined packing of rods and spheres in reinforcing plastics. Ind. Engng. Chem. Prod. Res. Dev., 1978, vol. 17, pp. 363-366.
8. Nardin M., Papirer E. Contribution a l'etude des empilements au hasard de fibres et/ou de particules spheriques. Powder Technology, 1985, vol. 44, pp. 131-140.
9. Evans K.E., Gibson A.G. Prediction of the maximum packing fraction achievable in randomly oriented short-fibre composites. Composites Science and Technology, 1986, vol. 25, pp. 149-162.
10. Evans K.E., Ferrar M.D. The packing of thick fibers. J. Phys. D: Appl. Phys., 1989, vol. 22, pp. 354-360.
11. Parkhouse J.G., Kelly A. The Random Packing of Fibres in Three Dimensions. Proc. R. Soc. Lond. A., 1995, vol. 451, pp. 737-746.
12. Philipse A.P. The random contact equation and its implications for (colloidal) rods in packings, suspensions, and anisotropic powders. Langmuir, 1996, vol. 12, pp. 1127-1133.
13. Wegner S., Stannarius R., Boese A., Rose G., Szabo B., Somfai E., Börzönyi T. Effects of grain shape on packing and dilatancy of sheared granular materials. Soft Matter, 2014, vol. 10, pp. 5157-5167.
14. Donev A., Cisse I., Sachs D., Variano E.A., Stillinger F.H., Connelly R., Torquato S., Chaikin P.M. Improving the density of jammed disordered packings using ellipsoids. Science, 2004, vol. 303, pp.990-993.
15. Man W., Donev A., Stillinger F.H., Sullivan M.T., Russel W.B., Heeger D., Inati S., Torquato S., Chaikin P.M. Experiments on random packings of ellipsoids. Physical Review Letters, 2005, vol. 94(19), pp. 198001-1 - 198001-4.
16. Xia C., Zhu K., Cao Y., Sun H., Kou B., Wang Y. X-ray tomography study of the random packing structure of ellipsoids. Soft Matter, 2014, vol. 10, pp. 990-996.
17. Schaller F.M. Neudecker M., Saadatfar M., Delaney G.W., Schröder-Turk G.E., Schröder M. Local Origin of Global Contact Numbers in Frictional Ellipsoid Packings. Physical Review Letters, 2015, vol. 114(15), pp. 158001-1 - 158001-5.
18. Bennett C. H. Serially depositied amorphous aggregates of hard spheres. J. Appl. Phys., 1972, vol 32(6), pp. 2727-2734.
19. Sherwood J.D. Packing of spheroids in three-dimensional space by random sequential addition. J. Phys. A: Math. Gen., 1997, vol. 30, pp. 839-843.
20. Lubaehevsky B.D., Stillinger F.H. Geometric properties of random disk packings. Journal of Statistical Physics, 1990, vol. 60, pp. 561-583.
21. Donev A., Torquato S., Stillinger F.H. Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. I. Algorithmic details. Journal of computational physics, 2005, vol. 202 (2), pp. 737-764.
22. Donev A., Torquato S., Stillinger F.H. Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. II Application to ellipses and ellipsoids. Journal of computational physics, 2005, vol. 202 (2), pp. 765-793.
23. Ghossein E., L'evesque M. Random generation of periodic hard ellipsoids based on molecular dynamics: a computationally-efficient algorithm, Journal of Computational Physics, 2013, vol. 253, pp. 471-490.
24. Lubachevsky B.D., Stillinger F.H., Pinson E.N. Disks vs. spheres: contrasting properties of random packings. Journal of Statistical Physics, 1991, vol. 64, pp. 501-524.
25. Jodrey W.S. Computer simulation of close random packing of equal spheres, Physical Review A, 1985, vol. 32, no. 4, pp. 2347-2351.
26. Williams S.R., Philipse A.P. Random packings of spheres and spherocylinders simulated by mechanical contraction. Physical Review E, 2003, vol. 67, pp. 051301-1 - 051301-9.
27. Stillinger F.H., Weber T.A. Inherent structure of liquids in the hard-sphere limit. J. Chem. Phys, 1985, vol. 83(9), pp. 4767-4775.
28. Zinchenko A. Algorithm for random close packing of spheres with periodic boundary conditions, J. Comp. Phys, 1994, vol. 114, pp. 298-307.
29. Kanit T., Forest S., Galliet I., Mounoury V., Jeulin D. Determination of the size of the representative volume element for random composites: statistical and numerical approach. International Journal of Solids and Structures, 2003, vol. 40, pp. 3647-3679.
30. Chan K. A simple mathematical approach for determining intersection of quadratic surfaces. Multiscale optimization methods and applications. New York, 2006, pp. 271-298.
31. Choi Y., Chang J., Wang W., Kim M., Elber G. Continuous collision detection for ellipsoids. Visualization and computer Graphics, IEEE Transactions, 2009, vol. 15, pp. 311-325.
32. Advani S.G., Tucker C.L. The use of tensors to describe and predict fiber orientation in short fiber composites. Journal of rheology, 1987, vol. 31, pp. 751-784.
