Об одном алгоритме синтеза упаковок частиц в рамках аристотелевской механики

УДК 51.72, 539.3

Об одном алгоритме синтеза упаковок частиц в рамках аристотелевской механики

А.Ф. Ревуженко, С.В. Клишин, О.А. Микенина

Институт горного дела им. Н.А. Чинакала СО РАН, Новосибирск, 630091, Россия

Предлагается динамический способ формирования трехмерных упаковок отдельных твердых сферических частиц. Частицы осаждаются внутри заданной области с фиксированной жесткой границей. Скорость движения каждой частицы пропорциональна ее весу и силам, возникающим от контакта с границей и соседними частицами. Проведен размерный анализ уравнений движения частиц. Для равновесной упаковки вычисляются ее средняя плотность и распределение координационного числа. Исследуется зависимость данных характеристик от коэффициента вязкости, гранулометрического состава и способа задания начальных условий (численный аналог способа засыпки материала в емкость).

Ключевые слова: механика Аристотеля, гранулированный материал, упаковка, метод молекулярной динамики, координационное число

Algorithm for particle packing based on Aristotelian mechanics

A.F. Revuzhenko, S.V. Klishin, and O.A. Mikenina

Chinakal Institute of Mining SB RAS, Novosibirsk, 630091, Russia

A dynamic algorithm is proposed for three-dimensional packing of spherical solid particles. The particles are deposited inside a specified region with a fixed rigid boundary. The velocity of each particle is proportional to its weight and forces due to contact of the particle with the boundary and neighbor particles. Dimensional analysis of the equations of particle motion is performed. The average density and coordination number distribution for equilibrium packing are calculated. The dependence of these characteristics on viscosity, granulometric composition, and representation of initial conditions (numerical analogue of material supply into a specified volume) is studied.

Keywords: Aristotelian mechanics, granular material, packing, molecular dynamics, coordination number

1. Введение

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

тивных лазерно-порошковых технологиях объемного формообразования (технологиях объемного синтеза), производстве керамики, цементной промышленности, порошковой металлургии, при создании композиционных материалов, состоящих из жидких и твердых фаз, и т.д.

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

Исследования в данной области имеют богатую историю и восходят к классическим работам по упаковке

© Ревуженко А.Ф., Клишин С.В., Микенина О.А., 2014

шаров (гипотеза Кеплера (1611 г.), проблема Гильберта (1900 г.), задача Фон-Коссена и др.

Теоретическая постановка задачи о размещении непересекающихся одинаковых шаров в пространстве звучит следующим образом: найти такой способ расположения шаров, при котором покрыта наибольшая доля этого пространства. В случае монодисперсных упаковок на плоскости наилучшим заполнением является размещение центров кругов в вершинах «паркета», образованного правильными шестиугольниками, в котором каждый круг окружен шестью другими [1]. Координационное число данной упаковки хп = 6, а относительная плотность £ = 0.9069. Здесь и ниже £ — отношение суммарного объема частиц к общему объему. Для монодисперсных упаковок в трехмерном пространстве максимальная плотность упаковки £ = 0.7405, а хп = 12.

Экспериментальные исследования плотности упаковок, состоящих из случайно засыпанных под действием силы тяжести сферических частиц одинаковых размеров, изготовленных из различных материалов, представлены в работах [1-3]. Было показано, что в случайной упаковке существует предельная плотность, значение которой составляет £ = 0.64, а среднее координационное число хп изменяется в пределах от 6.92 до 9.14. Видно, что значения плотности значительно меньше, чем известные значения для идеальных упаковок.

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

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

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

Представленные динамические способы создания упаковок предполагают учет всех соударений частиц во время процесса уплотнения, до тех пор пока не будет достигнуто равновесное распределение контактных сил по всему объему. Поскольку движение каждой частицы должно быть описано в течение всего процесса уплотнения, то эти методы требуют большого количества вычислений, иногда сравнимого с проведением последующего численного эксперимента. Применению динамических методов при различных граничных условиях посвящена обширная литература [4-10].

Геометрический подход к созданию плотных исходных полидисперсных упаковок частиц является альтернативой динамическим алгоритмам. Основная особенность этих схем состоит в том, что упаковки создаются на основе теоретических геометрических расчетов, без имитации динамики движения и соударения частиц. Как правило, процесс упаковывания происходит последовательно для каждой частицы. Положение новой частицы рассчитывается на основе известных положений уже упакованных частиц таким образом, что она касается определенных ранее уложенных частиц в достаточном числе контактов. Эти алгоритмы разработаны главным образом для двух измерений, но в настоящее время разработаны алгоритмы и для случая трехмерного пространства [11-14]. Преимуществом такого подхода является вычислительная эффективность, особенно при укладывании большого числа частиц. С другой стороны, процесс упаковывания часто приводит к большому количеству дислокаций между частицами или вблизи границ, причем значения среднего координационного числа таких упаковок получаются довольно низкими. Обзоры геометрических методов построения исходных упаковок даны в работах [15, 16]. Численные исследования плотности и координационных чисел упаковок, созданных на основе указанных подходов, показали удовлетворительное совпадение с экспериментальными данными.

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

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

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

2. Постановка задачи

Пусть иг- — перемещение центра частицы с номером i, — сила, действующая на данную частицу со стороны частицы с номером j, t — время. Как отмечалось, в рамках механики Аристотеля скорости пропорциональны силе:

Эи, ^ 4 3

=? F- + 3п rYe'

где Ц, — эффективная вязкость, которую будем считать постоянной, т.е. Ц, = const; у, — удельный вес i-й частицы; e — единичный вектор в направлении силы тяжести. Величина Ц, зависит от вязкости жидкости, в которой находится частица, а также от ее размера и формы. Например, для частицы в форме шара формула Стокса дает

du,

6п r, v—- = F,,

г dt '

где v — динамическая вязкость жидкости; F, — сила, действующая на частицу с номером i (сила трения, также называемая силой Стокса). В формуле (2) предполагается, что F, = const и движение частицы является установившимся. Следовательно, эффективная вязкость для i-й частицы равна Ц, = 6щ v.

Предположим далее, что величина перекрытия частиц с номерами i и j равна 8-. В первом (линейном) варианте будем считать, что

F = k- 8-n- пРи 8- ^

F- = 0 при 8--< 0, где k- = const; n- — единичный вектор, направленный вдоль линии, соединяющей центры частиц i и j. Второе условие (3) означает, что частицы между собой контактировать перестали. Параметр k- зависит от размеров и свойств материала частиц. Задача сводится к решению

уравнений (1), (3) при заданных начальных условиях (заданы начальное положение всех частиц и их радиусы). Нагружение сводится к заданию направления действия внешней силы е (у >0 в уравнении (1)).

Вначале сделаем размерный анализ данной задачи. Параметры среды, фигурирующие в уравнениях (1) и (3), дают естественные масштабы для переменных размерности длины, времени и силы (отмеченные индексом 0):

г0, г0 = ц/к, Р0 = кг0

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

Х = -

kr 0

= Е8 -n - +^e> Fa =8n -.

4/3 пу(г0)3'

Если использовать указанные масштабы и соответствующие безразмерные переменные, то придем к следующей системе:

ди

дг

Здесь безразмерные переменные обозначены так же, как и размерные. Механический смысл числа X очевиден: это отношение характерной силы отталкивания, которая возникает вследствие перекрытия частиц, к силе тяжести. Величина г0 — это время релаксации силы, возникающей при перекрытии частиц. Релаксация происходит вследствие удаления частиц друг от друга по закону Аристотеля (без учета силы тяжести).

Интересно отметить, что размерные параметры задачи дают еще одну характерную величину, которая явно в уравнениях не фигурирует. Это величина М0:

М0 =ц2Д,

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

Ее смысл можно интерпретировать следующим образом. В равенстве (5) не фигурирует удельный вес частиц у. Поэтому рассмотрим размерные уравнения без учета веса у и в одномерном варианте: уш' = Р, Р = = k8. Штрихом обозначена производная по времени. Отсюда

ц = Р/и, к = Р/8.

Подстановка (6) в (5) дает

М У2 = Р8 2 2'

Справа в равенстве (7) стоит энергия, которая высвободится при уменьшении перекрытия частиц от величины 8 до 0, а слева — кинетическая энергия тела, обладающего скоростью и'. Величина М 0 — это масса, обеспечивающая закон сохранения энергии в указанной ситуации.

3. Реализация алгоритма

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

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

Был проведен ряд численных экспериментов для различного гранулированного состава частиц на основе представленных двух способов засыпки. Эксперименты показали отличие кинематической картины процесса засыпки, предшествующей равновесному состоянию в рамках представленной модели (1), (3) и классической постановки метода дискретных элементов на основе ньютоновской механики [17]. В качестве примера на

рис. 1 приведены этапы засыпки полидисперсной системы, состоящей из 80 000 частиц. Здесь цветом показаны вертикальные компоненты векторов скоростей отдельных частиц. Видно, что в материале образуется фронт, движущийся в противоположном направлении от внешнего воздействия, ниже которого частицы практически покоятся. Этапы численного моделирования в рамках ньютоновской механики представлены на рис. 2. Для создания равновесной упаковки и устранения колебательных процессов в данном случае вводились либо небольшая по величине контактная вязкость, либо малое сухое трение между частицами. В отличие от представленной в настоящей работе аристотелевской схемы видно, что вследствие ударов частиц о дно емкости или о частицы, которые находятся на свободной поверхности, происходит сжатие материала, при котором частицы меняют направление на противоположное и процесс усадки приобретает периодический характер, затухающий со временем. При этом в материале отсутствует выраженная граница, за которой частицы находятся в состоянии равновесия.

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

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

vz, м/с 0.00

-0.25 -0.50 -0.75 -1.00

Рис. 1. Вертикальная компонента вектора скорости в фиксированные моменты времени при засыпке полидисперсной системы сферических частиц в цилиндрическую емкость в рамках аристотелевской механики

Рис. 2. Этапы процесса засыпки полидисперсного набора частиц в цилиндрическую емкость в рамках ньютоновской механики (а-в); конечное равновесное состояние (г). Стрелками показано направление движения материала

цами и свободной поверхностью. Поскольку свободная поверхность может приобретать случайную форму, для ее определения использовался алгоритм построения поверхностных сеток, представленный в работе [18]. Результат применения данного алгоритма показан на рис. 3. Здесь можно отметить, что при указанных способах засыпки и заданных выше отношениях размеров емкости к размерам частиц (более 50) в проведенных численных экспериментах конфигурация свободной поверхности была близкой к плоской. Поэтому объем, занятый частицами, можно вычислять как площадь основания контейнера, умноженную на высоту засыпки. Сравнение результатов показало, что полученный таким способом объем отличается не более чем на 5 % от объема, полученного с помощью алгоритма построения сеток.

Для проведенных численных экспериментов относительная плотность упаковки составила 0.62-0.64 для монодисперсных наборов частиц, что соответствует максимальным значениям для случайных засыпок. В случае полидисперсных упаковок плотность составила 0.65-0.71, увеличиваясь с увеличением отношения максимального радиуса частиц к минимальному (т.е. с расширением спектра размеров частиц).

Вторая характеристика — это число контактов, которые имеет каждая фиксированная частица с соседними частицами либо с границей (координационное число). Упаковка, состоящая из N частиц, дает N координационных чисел х, * = 1, 2, ..., N. Среднее координационное

число определяем, как

* 1 N

х = Т7 ^ х'

N¿=1

Определим плотность распределения координационных чисел / (х) как отношение количества частиц с фиксированным числом контактов к общему количеству частиц. На рис. 4 приведены характерные диаграммы плотности распределения fкоординационных чисел монодисперсных и полидисперсных наборов частиц. Здесь сплошной линией показан график плотности,

Рис. 3. Определение объема, занятого дискретными элементами на основе алгоритмов построения поверхностных сеток: цилиндрическая емкость

Рис. 4. Плотность распределения координационных чисел. Радиусы частиц равны между собой

полученной на основе механики Аристотеля, а штриховой — на основе механики Ньютона.

4. Выводы

Алгоритмы создания упаковок частиц в рамках механики Аристотеля являются достаточно эффективными и имеют ряд преимуществ по сравнению с алгоритмами, реализуемыми на основе ньютоновской механики. В случае монодисперсных систем частиц относительная плотность приближается к максимально возможной для случайных упаковок (0.62-0.64). Расширение спектра размеров частиц приводит к увеличению относительной плотности до значений 0.65-0.71.

Работа выполнена при финансовой поддержке РФФИ (проект № 13-05-00432).

Литература

1. Дересевич Г. Механика зернистой среды // Проблемы механики: сб. статей / Под ред. X. Драйдена, Т. Кармана. - М.: Изд-во иностранной литературы, 1961. - Вып. 3. - С. 92-152. Deresiewicz H. Mechanics of Granular Matter, V. 3, Advances in Applied Mechanics / Ed. by H.L. Dryden, Th. Karman. - New York: Academic Press, 1958. - P. 233-303.

2. Scott G.D., Kilgour D.M. The density of random close packing of spheres // British J. Appl. Phys. (J. Phys. D). - 1969. - V. 2. - No. 6. -P. 863-866.

3. Bernal J.D. The Bakerian Lecture, 1962. The Structure of Liquids // Proc. Royal Soc. London. A. Math. Phys. Sci. - 1962. - V. 280. -No. 1382. - P. 299-322.

4. Guises R., Xiang J., Latham J.-P., Munjiza A. Granular packing: numerical simulation and the characterization of the effect of particle shape // Granul. Matter. - 2009. - V. 11. - No. 5. - P. 281-292.

5. Makse H.A., Johnson D.L., Schwartz L.M. Packing of compressible granular materials // Phys. Rev. Lett. - 2000. - V. 84. - P. 4160-4163.

6. SilbertL.E., Erta§ D., Grest G.S., Halsey T.C., Levine D. Geometry of frictionless and factional sphere packings // Phys. Rev. E. - 2002. -V. 65. - No. 3.

7. Potyondy D.O., Cundall P.A. A bonded-particle model for rock // Int. J. Rock Mech. Min. Sci. - 2004. - V. 41. - No. 8. - P. 1329-1364.

8. Salot C., Gotteland Ph., Villard P. Influence of relative density on granular materials behavior: DEM simulations of triaxial tests // Granul. Matter. - 2009. - V. 11. - No. 4. - P. 221-236.

9. Ergenzinger C., Seifried R., Eberhard P. A discrete element model to describe failure of strong rock in uniaxial compression // Granul. Matter. - 2011. - V. 13. - No. 4. - P. 341-364.

10. Fukumoto Y., Sakaguchi H., Murakami A. The role of rolling friction in granular packing // Granul. Matter. - 2013. - V. 15. - No. 2. -P. 175-182.

11. Потураев В.Н., Стоян Ю.Г., Шуляк И.А., Пономаренко Л.Д., Сансиев В.Г. О моделировании зернистой среды вычислительными методами // ФТПРПИ. - 1989. - № 2. - С. 3-8. Potuarev V.N., Stoyan Yu.G., Shulyak I.A., Ponomarenko L.D., San-siev V.G. Modeling of a granular medium by computer methods // Soviet Mining. - 1989. - V. 25. - No. 2. - P. 101-106.

12. Рыжков Ю.А., Гоголин В.А., Карпенко Н.В. Моделирование структуры массивов из кусковых и зернистых материалов (плоская задача) // ФТПРПИ. - 1992. - № 1. - С. 8-14.

Ryzhkov Yu.A., Gogolin, V.A., Karpenko N. V. Modelling the structure of solid masses of lump and granular materials (plane problem) // J. Min. Sci. - 1992. - V. 28. - No. 1. - P. 6-12.

13. Collins C.R., Stephenson K. A circle packing algorithm // Comp. Geom. - 2003. - V. 25. - No. 3. - P. 233-256.

14. Bagi K. An algorithm to generate random dense arrangements for discrete element simulations of granular assemblies // Granul. Matter. - 2007. - V. 7. - No. 1. - P. 31-43.

15. Stoyan D. Simulation and characterization of random systems of hard particles // Image Analysis Stereol. - 2002. - P. 41-48.

16. Stoyan D. Random systems of hard particles: models and statistics // Chinese J. Stereol. Image Analysis. - 2002. - V. 7. - No. 1. - P. 113.

17. Клишин С.В., Ревуженко А. Ф. Исследование задачи Янсена методом дискретных элементов в трехмерной постановке // ФТПРПИ. -2014. - № 3. - С. 10-16.

Klishin S. V., Mikenina O.A., Revuzhenko A.F. Deformation of granular material around a rigid inclusion // J. Min. Sci. - 2014. - V. 50. -No. 2. - P. 229-234.

18. Stukowski A. Computational analysis methods in atomistic modeling of crystals // JOM. - 2014. - V. 66. - No. 3. - P. 399-407.

Поступила в редакцию

__19.06.2014 г.

Сведения об авторах

Ревуженко Александр Филиппович, д.ф.-м.н., проф., зав. отд., зав. лаб. ИГД им. H.A. Чинакала СО РАН, revuzhenko@yandex.ru Клишин Сергей Владимирович, к.т.н., снс ИГД им. H.A. Чинакала СО РАН, sv.klishin@gmail.com Микенина Ольга Александровна, к.ф.-м.н., мнс ИГД им. H.A. Чинакала СО РАН, olgarev@yandex.ru