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

            Метод дискретного элемента (DEM), изначально разработанный Канделом и Стрэком [5], является наиболее мощным инструментом для расчета динамики большого количества частиц размера микрон и более. DEM основан на том, что каждая твердая частица рассматривается как отдельный объект и её движение и взаимодействия рассчитываются уравнениями Ньютона и Эйлера (Ур. 1-2).

		(1)
		(2)

            Сила F и момент M действующие на частицу i массы m_i которая имеет позицию с радиус вектором r_i а также угол поворота в пространстве φ_i. Момент инерции J_i представлен:

(3)

где L — это длина измерения взятая от центра массы (радиус сферической частицы), c — константа инерции (c = 2/5, если объект является сферой и c = 1/2, если объект является цилиндром или диском вращающимся вокруг его центра). Условие контакта частиц следующее:

(4)

где s_ij — пересечение частиц, R_i и R_j — радиусы частиц (рис. 1). Если условие контакта выполняется, то сила взаимодействия рассчитывается. Контакт частиц описывается моделями коллизий.
            В последнее время были разработаны различные коммерческие и открытые бесплатные системы для DEM моделирования. Их сравнение приведено в табл. 1. Последующий анализ показал, что ни одна из существующих систем с открытым кодом не может быть эффективно применена для моделирования процессов сыпучих материалов. Главной целью представленного проекта является разработка новой параллельной кросс-платформенной системы для DEM моделирования процессов твердых тел с хорошим соотношением с экспериментальными данными.

Таблица 1

Сравнение существующих сред DEM моделирования

	EDEM	PFC 3D	Chute maven	ELFEN	Passage	Ascalaph	Yade	LAMMPS	Pasimodo	Esys-Particle
Коммерческая	+	+	+	+	+	+	-	-	-	-
Распараллеливание	+	+	-	+	-	+	+	+	+	+
Совместимость с CFD	+	-	-	-	-	-	+	+	-	-
Совместимость с FEM	+	-	-	+	-	-	+	+	-	-
Импорт CAD	+	-	+	+	-	-	-	+/-	-	-
Различные контактные модели	+	+	-	+	+	+	-	+	+	+
Установка форм частиц	+	+	-	+	+	+	+	+	+	-
3-х мерность	+	+	+	+	+	+	+	+	+	-
Информация о частицах	+	+	+	+	+	+	+	+	+	-

Рисунок 1 — Схема коллизии модели упругих сфер

            Существуют две основные стратегии столкновения частиц: упругая и твердая модель сфер. Модель твердых сфер наиболее часто используется для моделирования разреженных энергетических потоков частиц. Для выполнения более реалистичного моделирования используется модель упругих сфер (рис.1) [8].
            Силы, которые действующие на частицу, могут быть описаны в зависимости от поведения материала частицы во время столкновения. Сила взаимодействия делится на нормальную (ур. 5, 6, 8) и тангенциальную составляющую
(ур. 10, 12). Нормальная сила делится на эластичную и рассеивающую часть.

(5)

где k и — эластичная и рассеивающая константы. s — пересечение в направлении нормали. Эластичная часть нормальной силы, представленная Герцем [1] и рассеивающая часть , представленная Бриллиантовым [6] приведены в ур. 6:

(6)

где E* — эффективный модуль Юнга частиц (индекс i и j), R* — эффективный радиус, n — коэффициент Пуассона, Ad — функция вязкости материала:

(7)

h — рассеивающий коэффициент сталкивающихся частиц. E — модуль Юнга частицы.
В ур. (6-7) E_i = E_j = E*, n_i = n_j.
Нелинейная модель нормальной силы была представлена Кувабара и Коно [3]:

(8)

Линейная модель нормальной силы была представлена Уолтоном и Брауном [7]:

(9)

где l, u — индексы загруженного и разгруженного состояния частиц, k_l, k_u — жесткость пружины во время загрузки и разгрузки. s₀ — максимальная результирующая деформация взаимодействующих частиц.
Модель тангенциальной силы была представлена Хаффом и Вернером [1]:

(10)

где — относительная скорость центров сфер (ур. 11). — коэффициенты рассеивания для нормальной и тангенциальной силы, m — коэффициент трения.

(11)

Линейная модель тангенциальной силы была также представлена Ди Майо и Ди Ренцо [2]:

(12)

где k^t — тангенциальная жесткость. s^t — смещение в тангенциальном направлении:

(13)

            Наилучшее соотношение с экспериментальными данными показала расширенная линейная модель нормальной силы представленная Уолтоном и Брауном [3]. Более того, использование линейной модели ускоряет моделирование по сравнению с нелинейными моделями. Линейная тангенциальная модель представленная Ди Майо и Ди Ренцо [2] является одной из наиболее соответствующих экспериментальным данным.

            Интегрирование уравнений движения Ньютона может быть выполнено с помощью схемы интегрирования Gear. Алгоритм Gear (рис. 2) подходит для решения данной проблемы, в основном в силу его численной стабильности [1].

Рисунок 2 — схема интегрирования Gear

Первым шагом (в цикле) в алгоритме предсказываются позиции частиц, скорости, и производные ко времени t+Δt используя разложение в ряд Тейлора (ур. 14).

		(14)

После чего предсказанные координаты и производные используются для расчета и . Из сил и моментов вычисляются линейные и угловые ускорения и (ур. 15). Чтобы ускорить вычисления был разработан алгоритм фиксирования коллизий основанный на Verlet списках.

(15)

            В этом алгоритме силы рассчитываются между частицами, которые являются близкими соседями, т.е. дистанция между ними меньше verlet дистанции (ур. 16)

< verlet дистанция

(16)

            Verlet список составляется для каждой частицы и содержит индексы соседних частиц. Двойное включение в список Verlet исключается, используя ограничение, что частица с индексом i имеет только соседей с индексом
j < i. Список необходимо составлять заново в случае, если частица переместилась на расстояние, которое равно половине verlet дистанции.
             Пространство моделирование делится сеткой на ячейки. Размер каждой ячейки не должен превышать размера наибольшей частицы. Условие (ур. 16) проверяется только для частиц в соседних ячейках.
            В случае если частицы движутся по направлению друг к другу с максимальной скоростью, они могут столкнуться не присутствуя в списке Verlet. Таким образом необходимо сохранять коордиднаты частиц и списки сосдей когда составляется новый список Verlet. В худшем случае моделирование должно быть возвращено к шагу, где был составлен последний удачный список и продолжить с увеличенной Verlet дистанцией.
            В дополнение к алгоритмам оптимизации может быть использовано распараллеливание для ускорения вычислений.

            Разработка стратегии распараллеливания играет важную роль, т.к. необходимо моделировать большие количества частиц. Для этого объем моделирования разбивается на сеть пересекающихся ячеек. После чего расчеты распараллеливаются с учетом разбиения объема моделирования. Все частицы в ячейке или в группе ячеек рассчитываются независимо на отдельном процессоре.
            Эффективность стратегии распараллеливания зависит от минимизации двух основных критериев: дисбаланса загрузки и объема передаваемых данных. Минимизируя дисбаланс загрузки, расчеты должны быть равномерно распределены между процессорами. Это может быть достигнуто с помощью не постоянного размера сетки, где размер ячейки пропорционален количеству частиц в ячейке (рис. 4). Минимизирование объема передаваемых данных между процессорами может быть достигнуто используя сетку с движущимися границами. Миграция частиц между процессорами может быть минимизирована, если сетка движется в том же направлении, что и поток частиц.
            Основная часть компьютерного времени в моделировании сыпучих материалов расходуется на векторные операции. Использование SSE (Streaming SIMD Extensions) инструкций может сделать расчеты с векторами и матрицами до 4х раз быстрее.

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

Рисунок 4 – DEM-моделирование заполнения и освобождения бункера в 2D и 3D
Параметры анимации: количество кадров: 6; количество циклов: 7; Размер: 39 Кб; Сделано в AVI->GIF Converter & GIF-Animator

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

            1. T. Poeschel, T. Schwager. Computational granular dynamics. Models and algorithms. Springer, 2005.
            2. H. Kruggel-Emden, S. Wirtz, V. Scherer. A study on tangentional force laws applicable to the discrete element method for materials with viscoelastic or plastic behaviour, Chemical Engineering Science 63, 2008, 1523-1541.
            3. H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, V. Scherer. Review and extension of normal force models for DEM, Powder Technology 130, 2007, 157-173.
            4. S. Antonyuk, M. Khanal, J. Tomas, S. Heinrich and L. Moerl. Impact breakage of spherical granules: experimental study and DEM simulation, Chemical Engineering and Processing 45, 2006, 838-856.
            5. P.A. Cundall, O.D.L. Strack. A discrete numerical model for granular assemblies. Geotechnique 29, 1979, 47-65.
            6. N.V. Brilliantov, F. Spahn, J.-M. Hertzsch, T. Poeschel. Model for collision in granular gases, Physical review E, Vol. 53, 1996.
            7. Walton, O.R., Braun, R.L. Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. Journal of Rheology, 30, 1986, 949- 980.
            8. Antonyuk, S., Palis, S., Heinrich, S. Breakage behaviour of agglomerates and crystals by static loading and impact, Powder Technology, Special Issue: Agglomeration, in Press, February 2010.
            9. Poschel, T., Saluena, C., Schwager, T., 2001. Scaling properties of granular materials. Physical Review E 64 (1), (Art. No. 011308 Part 1).
            10. Schaefer, J., Dippel, S., Wolf, D.E., 1996. Force schemes in simulations of granular materials. Journal De Physique I 6 (1), 5–20.

---

            При написании данного автореферата магистерская работа еще не завершена. Окончательное завершение: декабрь 2010 г. Полный текст работы и материалы по теме могут быть получены у автора или его руководителя после указанной даты.