CC BY
248
Моделирование случайной упаковки шаров
И.Г. Дик, Е.Н. Дьяченко1, Л.Л. Миньков1
Университет Эрланген-Нюрнберг, Эрланген, 91054, Германия 1 Томский государственный университет, Томск, 634050, Россия
Цель данной работы — проиллюстрировать возможности численного моделирования случайных упаковок. Генерация случайной упаковки заключается в последовательном запуске частиц сферической формы со случайно выбранной координатой на верхней грани куба с ребром единичной длины по направлении его дна. Предусмотрена возможность заполнять куб полидисперс-ной совокупностью частиц, задаваемой в виде любой монотонно возрастающей аналитической функции распределения частиц по размерам. Учитывается, что при столкновении сферы с одной из уже осевших частиц между ними происходит адгезионное взаимодействие, в результате которого рассматриваемая сфера может либо прилипнуть к уже осевшей, либо изменить вектор своей скорости и продолжить движение в зависимости от сил, действующих вдоль поверхности контакта.
Подтверждено известное из эксперимента наличие колебания локальной пористости у стенок контейнера. Проанализирована структура пор для случаев наличия и отсутствия адгезионной силы между шариками упаковки. Обнаружена анизотропия структуры упаковки, связанная с действием адгезии. Подтверждено известное из экспериментов повышение пористости упаковки при наличии адгезии. Показана роль «размытости» распределения частиц по размерам на величину образуемой ими упаковки. Проанализировано и дано сравнение с экспериментом влияния состава разноразмерных шаров на пористость упаковки. Показана возможность моделирования работы фильтра очистки воды, в частности, его блокирования.
Simulation of a random sphere packing
J.G. Dueck, E.N. Dyachenko1, and L.L. Minkov1
Friedrich-Alexander-University of Erlangen-Nuremberg, Erlangen, 91054, Germany 1 Tomsk State University, Tomsk, 634050, Russia
The aim of this paper is to illustrate the possibilities of numerical simulation of random packings. The generation of a random packing consists in the sequential starting of spherical particles with a randomly selected coordinate on the upper cube face having unit edge length down to its bottom. We provide for the possibility to fill the cube with a polydisperse mixture of particles which is given by any monotonically increasing analytic particle size distribution function. It is taken into account that after a sphere collides with one of the already settled particles there occurs adhesive interaction between them. As a result, the considered sphere can either stick to the settled particle or change the velocity vector and continue its motion depending on the forces acting along the contact surface.
The experimentally known fluctuation of local porosity near the container walls is confirmed. The structure of pores is analysed for the cases with and without the adhesive force between packing spheres. We have revealed anisotropy of the packing structure related with the adhesive forces. The experimentally known increase of packing porosity in the presence of adhesion is confirmed. The role of particle size distribution spreading by the value of the formed packing is shown. The effect of the composition of different size spheres on porosity is analysed and compared with the experiment. We demonstrate the possibility to simulate water filter operation and in particular filter blocking.
1. Введение
Характеристики пористой структуры играют существенную роль в некоторых природных явлениях и важны во многих отраслях техники (композитные материалы, гетерогенный катализ, фильтрация суспензий, транспорт и складирование зернистых материалов, седименты). Как правило, речь идет о нерегулярных, статистических упаковках. Такие упаковки в виде случайных ансамблей частиц естественным образом возникают в ходе смешивания, осаждения, напыления, спекания и т.д. Как пра-
вило, для их характеристики ограничиваются двумя макропараметрами: пористостью и удельной поверхностью, через которые должны быть выражены остальные интересные для практики величины (например, размер пор).
Часто именно пористость (плотность) и средний размер пор являются ответственными за прочностные характеристики композитных материалов [1]. Характеристики упаковок играют определяющую роль при доработке материалов после проведения основных про-
© Дик И.Г., Дьяченко Е.Н., Миньков Л.Л., 2006
цессов порошковой металлургии, например, при пропитке их легкоплавким металлом [2, 3]. Существенным для управления такими пост-процессами является наличие детальной информации об упаковке. Основным источником нахождения более детальных характеристик, включающих статистику размеров пор, их форму, влияние внешних сил на плотность упаковки, гранулометрического состава смеси и т.д., является эксперимент. Теоретическое изучение случайных упаковок на основе моделирования силового взаимодействия частиц, ведущего к возможности их самоорганизации и образования специфических мезомасштабных структур [4], естественно делать, привлекая статистические методы, проводя компьютерный эксперимент [5-12].
Цель данной работы заключается в иллюстрации возможности численного моделирования случайных упаковок из шаров различных размеров.
2. Численная модель
Генерация случайной упаковки заключается в последовательном запуске частиц сферической формы со случайно выбранной координатой на верхней грани куба с ребром единичной длины (рис. 1) в направлении его дна. Начальная скорость оседания для всех частиц задается одинаковая.
В рассматриваемой модели предусмотрена возможность заполнять куб совокупностями частиц трех видов: монодисперсной, бидисперсной и полидисперсной. Последняя может задаваться в принципе любой монотонно возрастающей функцией распределения частиц по размерам, в виде аналитической функции Qi (г), где / означает род распределения: / = 0 — числовая доля частиц радиуса меньше г, / = 2 — доля площади поверхности частиц радиуса меньше г, / = 3 — объемная доля этих частиц. В частности, в качестве функции распределения удобно применять так называемую функцию ИИ^В [13]:
Qi (r) =1 - exP
r
rm
Случай монодисперсной совокупности частиц можно рассматривать как предел при а ^ ж При этом гт будет отвечать радиусу частиц монодисперсной совокупности. Аналогично, предел композиции функций (1) в виде:
Qi(r) = 1 -ф exp
- (1 -ф)ехр
(2)
при ах и а2 ^ ж дает бидисперсную совокуп-
ность частиц с радиусами а (доля частиц этой фракции Ф) и V
Кроме начальной координаты каждой запускаемой частицы варьируется и ее радиус. Как правило, на практике задается функция объемного распределения Qз (г). Для расчетов необходим переход к штучному распределению, что осуществляется переходом:
/
s - 3 Qs> *
Qo(r) = -
ds
\
s - 3 dQäis) ds
ds
(3)
Выбирая случайным образом Q0(r) между 0 и 1, находим из (3) радиус шара г.
В дальнейшем производится слежение за движением запущенного шара. Чтобы исключить расчет движения частицы от верхней грани куба до первого ее контакта с уже сложившейся упаковкой, запуск шара реально производится с некоторой точки внутри куба над упаковкой, в которой движущийся шар достаточно близко находится от уровня упакованных частиц. Координата 2 подконтрольного шара принимает значение Z = z + + г + R, где z — координата центра самого высокосидя-щего из уже осевших шаров; г — радиус движущегося шара и R — радиус частицы из упаковки с центром в z. Координаты X и У определяются генератором случайных чисел, который с равной вероятностью принимает значения от 0.5 - (X - г) до 0.5 + (X - г), где X — наперед задаваемое число между г и 0.5. Установка числа X позволяет задавать размер источника запускаемых частиц. При X = г запускаемые частицы находятся на линии X = 0.5 и У = 0.5.
При столкновении шара с одной из уже осевших частиц между ними происходит адгезионное взаимодействие, в результате которого рассматриваемый шар может либо прилипнуть к уже осевшей частице, либо изменить вектор своей скорости и продолжить движение. Условие остановки (прилипания) частицы определяется равновесием сил, действующих вдоль поверхности контакта (рис. 2).
Таким образом может моделироваться взаимовлияние различных иерархических уровней [4]: микроуровня (здесь — характер и силы взаимодействия отдельных частиц), мезоуровня (строение пор, определяющее конфигурацию силовых векторов) и макроуровня (пористость, плотность материала, зависящяя от количества и величины пор).
Предполагается, что на частицу, находящуюся в точке (х, у, z), могут действовать следующие силы:
1. Сила адгезии, которую можно записать в виде [14]: 4гЯ
ad
С
r + R
где CH — константа Г амакера; r и R — радиусы контактирующих шаров.
2. Сила гидродинамического увлечения частицы потоком (в случае моделирования образования осадка при фильтрации), например, записанная в виде [15]:
Fh =6пцги,
где п — вязкость жидкости; v — скорость потока.
3. Сила тяжести частицы. Ее проекция на плоскость касания есть mg cos а, где а — угол между прямой, проходящей через центры взаимодействующих сфер, и горизонтальной плоскостью. Учет тяжести свода упаковки производится прибавлением к весу частицы веса вышележащих слоев m + m(1 - z)(е), где (е) — средняя пористость слоя между z и 1. В случае зависимости е от z (имеет место при наличии адгезионных сил) для нахождении применяются итерации.
В общем случае условие равновесия можно записать в виде критического значения для угла: cos а < A, где A — коэффициент адгезии, представляющий собой отношение сил, удерживающих частицу на поверхности и пытающихся сдвинуть ее с места. В программе условие предела покоя можно формулировать достаточно широким образом в виде некоторой аналитической функции.
Рис. 2. Схема сил, действующих на частицу, контактирующую с частицей из упаковки
Если частица не прилипла, то она продолжает свое движение по касательной к поверхности осевшего шара со скоростью, определяемой по формуле: V’= = V - п^п), где V’ — скорость частицы после столкновения; V — скорость частицы до столкновения; п — единичный вектор, выходящий из центра осевшего шара к центру движущегося шара.
Условием остановки движения частицы, кроме прилипания к уже осевшей частице, является касание нижней грани куба либо достижение устойчивого положения с тремя (или более) точками опоры при касании других частиц или граней куба. После остановки сфера не может изменять свое положение.
Упаковка частиц продолжается до полного заполнения куба.
3. Результаты исследования упаковок
Представление результатов проведем для различных случаев. Радиусы шаров в приведенных ниже примерах даны в долях длины ребра заполняемого куба.
3.1. Монодисперсные упаковки
На рис. 3 показан пример упаковки шаров одинакового радиуса 0.05, сгенерированной при нулевом уровне адгезионных сил (А = 0). На рисунке показан вид сверху на шары, лежащие ниже уровня z = 0.6. Шары изображены тем темнее, чем глубже от поверхности z = 0.6 они расположены.
На рис. 4 дан вид при z = 0.6. Черные круги — сечения шаров. Незакрашенная площадь дает пористость в сечении (просветность), которая статистически равна пористости объема. Исследуя поверхностную пористость различных сечений можно убедиться в том, что локальная пористость вблизи стенок является периодической функцией расстояния до стенки. Иллюстрация этому дана на рис. 5. Там же нанесены экспериментально измеренные значения локальности пористости [15]. У самой стенки частицы имеют только точки касания со стенкой и просветность равна 1. Следующий максимум уже гораздо ниже единицы находится при х = 2г. При дальнейшем удалении от стенки пристеночный порядок теряется и значение пористости выходит на постоянное значение 8 ~ 0.42.
Чтобы исключить влияние стенки, в дальнейшем при вычислении объемной пористости берется куб, образованный из исходного единичного удалением поверхностных слоев с каждой стенки толщиной примерно 4-5 диаметров шаров.
Для каждой конкретной упаковки результат является статистическим. Для получения устойчивой осреднен-ной характеристики упаковки исследуется набор упаковок с одинаковыми параметрами. Обычно достаточно рассчитать порядка 10 упаковок.
В трехмерной упаковке (в отличие от двухмерной [9-12]) дать определение границ поры затруднительно.
Рис. 3. Упаковка монодисперсных шаров
0.9
0.7
0.5
0.3
\ 7
|—1 1-І п □ п П і—і і-і
□ 1-І 1-І □ □ 1-І—
х/2 г
Рис. 5. Пористость монодисперсной упаковки в зависимости от расстояния до стенки. Квадратные символы — экспериментальные данные [14]; линии-расчеты: сплошная — в направленииX, пунктирная — в направлении 2
мерная пористость) равны между собой. Аналогично и объемная пористость равна одномерной и двухмерной [5].
Пример вычисления функции распределения DL по длинам пор L дан на рис. 6. Видно, что распределения длин прерывистости вдоль вертикальной и горизонтальной осей практически одинаковы (изотропия). Максимальная длина пор соответствует примерно радиусу шариков, образующих упаковку. Доля DL таких пор достигает 65 %. Пор с длиной, большей 4 радиусов, уже крайне мало.
Зная доли длин пор DL i для отдельных классов Li, можно вычислить и среднее значение для длины поры в сечении упаковки (ь) = X DL iLi. Для приведенного примера оно равно 0.075, несколько больше чем длина наиболее часто встречающейся поры 0.051.
На структуру упаковки может существенно влиять адгезия между частицами.
Рис. 4. Сечение в упаковке монодисперсных шаров при z = 0.6
Поэтому для нахождения статистических характеристик пор используется лучевой метод. При этом вдоль плоскости некоторого выбранного сечения упаковки в произвольно выбранном направлении исследуется луч, проходящий от одной стенки куба до другой. Часть пути луч проходит внутри шаров (или в сечении — внутри кругов), а часть — вне. Тот участок пути, который луч проходит вне шара (прерывистость), считаем порой, а длину этого участка — длиной поры. Проанализировав множество таких лучей, можно построить функцию распределения по длинам пор L.
Нетрудно показать, что поскольку площадь пор в произвольном сечении, например z = С, может быть вычислена обычным интегрированием, например, вдоль
лучей в направлении
у, т.е. |Ь(х, С)ёх, то
прерывис-
тость (одномерная пористость) и просветность (двух-
0.05064 0.10128 0.15192 0.20256 0.2532 1_
Рис. 6. Доли длин прерывистости для лучей, пущенных из точки (0, 0.5, 0.5) в направлении оси X и точки (0.5, 0.5, 0) в направлении оси 2
Рис. 7. Упаковка монодисперсных шаров. Вид на глубине z = 0.6
Рис. 9. Распределение длин прерывистости для лучей, пущенных из точки (0, 0.5, 0.5) в направлении оси X и точки (0.5, 0.5, 0) в направлении оси 2. Уровень адгезии 0.5
Рис. 8. Сечение в упаковке монодисперсных шаров при z = 0.6
из них конактирует с 5-6 соседними, то при уровне адгезии А = 0.5 в основном с 2-3.
Заметим, что с ростом адгезии пористость растет — факт известный из экспериментов по образованию фильтрационных и седиментационных слоев [16-18]. На основании наших расчетов пористость может быть описана линейной функцией 8 = 8 0 + 0.4 А.
При возможно наиболее высоком А = 1 пористость может достичь значения порядка 0.8, что хорошо совпадает с экстраполяцией пористости при стремлении размера частиц в упаковке к нулю [14, 15].
3.2. Полидисперсные упаковки
Из рис. 5 видно, что вычисленные значения пористости упаковки несколько выше измеренных в приведенном численном эксперименте. Обычно для пористости шаров приводят значения между 0.38 и 0.43 [13-
На рис. 7 приведен вид упаковки, а на рис. 8 сечение упаковки при z = 0.6 для случая значения коэффициента адгезии А=0.5 (критический угол покоя а = п/ 3).
Сравнивая рис. 7, 8 с рис. 3, 4, видим, что при наличии эффекта слипания частиц, поры приобретают несколько вытянутую в вертикальном направлении форму (анизотропия). Это же подтверждает рис. 9, из которого очевидна анизотропия упаковки. Особенно следует подчеркнуть наличие длин прерывистости в направлении оси 2, многократно превышающих радиус слоеобразующих шаров.
Имеется возможность получить функцию распределения Dk координационного числа к — количества соприкосновений данной сферы с другими. Соответствующий график дан на рис. 10. Очевидно, что это распределение зависит от уровня адгезии. Если для упаковки из адгезионно невзаимодействующих шаров каждый
Рис. 10. Распределение координационного числа в упаковках с различным уровнем адгезионной силы
1.00
0.75
О
0.50 -
0.25
0.00 0.01
1 1 Y772 п НИ— а = 1.5, в = 0.35 \ / а - ц с-пяа LY
а = ос, 8 = 0.42 г /*
/1 П 9 у 1
г~г * * -Сг /
0.1
Г/Ггг
10
Рис. 11. Функции распределения частиц по размерам при различных т и соответствующие значения пористости упаковки
15]. Причины этого разброса могут быть различны. Рассмотрим одну из них.
На рис. 11 показаны кривые функций распределения частиц по размерам при различных значениях а. Там же даны и соответствующие значения пористости. Видно, что «размытость» функции распределения проявляется как уменьшение пористости.
Этот эффект уменьшения пористости смеси крупных и мелких шариков можно проиллюстрировать исследованием бидисперсных упаковок. На рис. 12 показана зависимость пористости бидисперсной упаковки от объемной доли мелкой фракции при различных соотношениях радиусов шариков. Точки — экспериментальные значения, взятые из [15, 19].
Видно, что качественно рассчитанные и измеренные зависимости совпадают. Систематическое завышение рассчитанных значений пористости можно, в частности, объяснить тем, что в реальных экспериментах, вероятно, трудно выдержать строго размеры частиц.
Интересное приложение моделирования связано с расчетом засорения фильтров, применяемых, например,
0.9
0.8
0.7
0.6
0.5
\ ф ш 9 ——
W 1 °
о\ /о
1 ■ ■О г1 /г2 = 0 ■О Г|/Г2 = 0 О eXD П/Г; .2 .5 , = 0.2
Ф ехр Г|/Г2 = 0.5
0.00
0.25
0.50
Ф
0.75
1.00
Рис. 12. Зависимость пористости бидисперсной упаковки от объемной доли мелкой фракции при различных соотношениях радиусов шариков. Точки — экспериментальные значения
для очистки воды. Фильтр представляет собой пористый слой бумаги, ткани или относительно крупнодисперсной засыпки. Пропуская через такой слой суспензию, можно задержать твердую фазу. Мелкие твердые частицы, удерживаясь в матрице фильтра, постепенно забивают поры, что приводит к повышению энергии, необходимой для прокачки суспензии, так что через определенное время фильтр приходится менять или чистить.
Мы ограничимся иллюстрацией расчета забивания фильтра, представляя его упаковкой шаров, генерированной описанным выше методом. Такого рода естественные фильтры применяются, например, при очистке шахтных вод [20, 21]. Представленный ниже пример может быть хорошей моделью и для фильтров с регулярной структурой. Моделируется прохождение достаточно мелких частичек через зернистый слой.
Рис. 13. Блокирование пор насыпи из частиц с радиусом, равным 0.05, мелкими частицами радиуса 0.015
Рис. 14. Блокирование пор насыпи из частиц с радиусом, равным 0.05, мелкими частицами радиуса 0.01
На рис. 13 и 14 показаны конечные состояния фильтра, после того как процесс блокирования пор закончился. Мелкие частички, примерно втрое меньшие размера пор, заблокировали фильтр, тогда как 2273 мелких частиц остались в фильтре при общем числе запущенных в него частиц 2 654 (т.е. лишь 15 % частиц смогли пройти через фильтр). При этом большинство частиц скопилось в верхнем слое фильтра. Таким образом, большая часть объема фильтра осталась незадействованной.
Для более мелких частиц характерна большая глубина проникания частиц в фильтр. В приведенном примере из 6 777 запущенных частиц фильтром на момент нарушения проходимости было задержано 2 777 (доля прошедших частиц равна 60 %).
В показанных примерах предполагалось, что прилипания мелких частиц к матрице фильтра не происходит. Учет адгезионной силы, естественно, существенно сказывается на работе фильтра, повышая, с одной стороны, его селектирующие характеристики, укорачивая, с другой, объем суспензии, переработанный до полного блокирования.
4. Выводы
Разработаны метод и программа моделирования образования насыпного слоя, седиментационного осадка, фильтрационного слоя или слоя напыления в виде случайной упаковки шаров.
Подтверждено известное из экспериментов наличие колебания локальной пористости у стенок контейнера. Проанализирована структура пор для случаев наличия и отсутствия адгезионной силы между шариками упаковки. Обнаружена анизотропия структуры упаковки, связанная с действием адгезии. Подтверждено известное из экспериментов повышение пористости упаковки при наличии адгезии.
Показана роль «размытости» распределения частиц по размерам на величину образуемой ими упаковки. Проанализировано и дано сравнение с экспериментом влияния состава разноразмерных шаров на пористость упаковки.
Показана возможность моделирования работы фильтра, в частности, его блокирования.
Другим приложением результатов является объяснение влияния состава дисперсного материала и физикохимических свойств поверхности частиц на плотность упаковок и соответственно прогноз характеристик твердого продукта, например осадков в фильтр-прессе.
Литература
1. Арзамасов Б.Н., Крашенинников А.И., Пастухова Ж.П., Рах-штадт А.Г. Научные основы материаловедения - М.: Изд-во МГТУ им. Н.Э. Баумана, 1994. - 336 с.
2. Dueck J., Niebling F, Neesse Th., Otto A. Infiltration as post-processing
of laser sintered metal parts // Powder Technology. - 2004. - No. 145. -P. 62-68.
3. Dueck J., Niebling F., Neesse Th., Otto A. Smoothing the path for rapid manufacturing // Metal Powder Report. - 2004. - No. 10. -P. 28-32.
4. Макаров П.В. Подход физической мезомеханики к моделированию процессов деформации и разрушения // Физ. мезомех. -1998. - T. 1.- № 1. - P. 61-81.
5. Аюкаев Р.И., Воробьев В.А., Кивран В.К. и др. Применение ЭВМ в
исследовании физико-структурных свойств пористых материалов. - Куйбышев: Куйб. инж.-строит. ин-т им. И. Микояна, 1976. -155 с.
6. Берлин А.А., Балабаев Н.К. Имитация свойств твердых тел и жидкостей методами компьютерного моделирования // Соросовс-кий образовательный журнал. - 1997. - № 11. - С. 85-92.
7. Винокуров Г.Г Компьютерное моделирование статистических характеристик макроструктуры газотермических покрытий // Физика и химия обработки материалов. - 2002. - № 3. - С. 29-32.
8. Хазан Г.Л., Злыгостев С.Н. Моделирование полидисперсной струк-
туры // Расплавы. - 1998. - № 4. - С. 49-51.
9. Дик ИГ, Дьяченко Е.Н., Миньков Л.Л. Моделирование образования
двухмерного насыпного слоя // Фундаментальные и прикладные проблемы современной механики: Доклады IV Всероссийской научной конференции, Томск, 5-7 октября 2004 г. - Томск: Изд.-во ТГУ, 2004. - С. 295-297.
10. Дик И.Г, Югов В.И. О моделировании структуры насыпного слоя // Инж. физ. журн. - 2005. - Т. 78. - № 2. - С. 36-43.
11. Дьяченко Е.Н., МиньковЛ.Л., Дик И.Г. Компьютерное моделирование случайных упаковок дисков // Изв. вузов. Физика. - 2005. -Вып. 11. - С. 83-91.
12. Okubo T., Odagaki T. Random packing of binary hard discs // J. Phys.: Condense Matter. - 2004. - V. 16. - No. 37. - P. 6651-6659.
13. Stiess M. Mechanische Verfahrenstechnik 1. - Berlin-Hidelbeig-New York: Springer-Verlag, 1995. - 369 p.
14. Israelachvili J. Intermolecular and Surface Forces. - San Diego: Academic Press Inc., 1995. - 450 p.
15. BrauerH. Grundlagen der Einphasen- und Mehrphasenstromungen. -Aarau: Sauerlander, 1971. - 955 p.
16. DueckJ.G., PurevjavD., KilimnikD.Yu. A contribution to the theory of porosity of fine-grained sediments // J. Eng. Phys. Thermophys. -2004. - V. 77. - No. 1. - P. 93-102.
17. Dueck J.G., Zvetanov E., Neesse Th. Porosity prediction for finegrained filter cake // Chem. Eng. & Technol. - 2000. - V. 23. -No. 1. - P. 18-22.
18. Dueck J., PurevjavD., Neesse Th. Surface force effects on the filter cake characteristics // Trans. Filt. Soc. - 2002. - V. 2. - No. 4. -P. 94-98.
19. JescharR. Druckverlust in Mehrkornschuttungen aus Kugeln // Archiv fur das Eisenhuttenwesen. - 1964. - No. 2. - P. 1-18.
20. ВасенинИ.М., ДьяченкоН.Н., Дьяченко Л.И. Осветление шахтных вод на слое песка // Изв. вузов. Горный журнал. - 2004. - № 6. -P. 51-55.
21. Бэр Я., Заславски Д., Ирмей С. Физико-математические основы фильтрации воды. - М.: Мир, 1971. - 280 с.
