Численное моделирование движения гранулированной среды в подвижных сосудах Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 16, № 2, 2011

Численное моделирование движения гранулированной среды в подвижных сосудах

И.О. Вогульский1, Н. А. Вогульская2 1 Институт, вычислительного моделирования СО РАН, Красноярск, Россия 2 Сибирский федеральный университет, Красноярск, Россия e-mail: bogul@icm.krasn.ru

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

Ключевые слова: моделирование, численные методы, гранулированная среда.

Настоящая работа является результатом сотрудничества авторов с инженерами-конструкторами, занимающимися созданием универсальных высевающих аппаратов вибрационного типа. Полученные инженерные решения поддержаны рядом публикаций [1, 2] и подтверждены патентами на изобретение [3, 4]. Главный критерий оптимальности работы конструкции — равномерность высева семян из отверстий лотка высевающего устройства. Модель, предложенная в работе, описывает движение зерен в мельницах и другие процессы, связанные с движением гранулированных сред.

Первые исследования по равномерности истечения гранул были связаны с попытками описания сыпучей среды в рамках модели потенциального течения идеальной жидкости [5, 6]. Это позволило получить ряд достаточно важных результатов, однако возникла необходимость имитационного, дискретного моделирования изучаемых процессов.

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

В основе предлагаемой модели движения гранулированной среды заложено упругое взаимодействие частиц следующей структуры. Будем считать, что каждая частица представляет собой абсолютно твердое тело (круг в плоском случае, шар в пространстве) массы га, радиуса г, окруженное достаточно тонкой упругой оболочкой. Коэффициент жесткости этой оболочки неизвестен, но его связь со сжатием е определяется из условия: если частица неподвижно лежит на поверхности (рис. 1), то сила упругости Еу полностью компенсирует вес Р:

се = тд. (1)

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

У

х

Рис. 1. Схема сил, действующих на частицу Рис. 2. Схема взаимодействия упругих частиц с упругой оболочкой

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

Модель заведомо нелинейная. Если расстояние между центрами частиц меньше диаметра — существует линейное упругое взаимодействие, если больше — взаимодействие отсутствует.

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

1. Алгоритм решения задачи

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

Сначала рассмотрим простейший вариант, не учитывающий силы трения. Пусть сосуд содержит N элементов, координаты центров которых х^ у^ % = 1, 2, 3,..., N известны.

Расчеты сил начнем с определения элементов, "ближайших" к г-му. Можно сразу вычислить расстояние между центрами элементов % и ^ = г.

и отобрать в качестве "ближайших" те, для которых Rij < д. Однако с точки зрения экономии вычислений более выгодной является следующая процедура: найдем модуль разности координат 1х-1 — xj• 3 = 1, 2, 3,..., N 3 = Если эта величина больше (I, элемент с номером ] больше не рассматривается. В противном случае вычисляются величины ^ — у?-1, 3 = 1, 2, 3,..., N 3 = Величину Rij определяем только в случае У — yj | < и в случае Rij < д, считаем элемент с номером ] лежащим в "ближайшем окружении". Таких элементов не может быть много. Запоминаем их номера и далее работаем только с ними.

Вычислим силу, действующую на /-П элемент со стороны у-го. Так как расстояние

меньше диаметра возникает упругая сила fij (см. рис. 2), равная

fij = (<& - Rij)с = (<1 — Rij)

тд ~2е

(2)

и действующая в направлении точки О по линии 00 . Единичный вектор в направлении 00 равен

xj xi Уj Уi

в1

Ri j Ri j

Единичный вектор е2, перпендикулярный е1, соответственно равен

yj yi xj xi

е2

Ri j Ri j

(3)

(4)

Следовательно, со стороны у-го элемента действует сила с компонентами

тд (xj — Хг)

2е Щ j

& = —— Rij)

тд (уз - уг) 2е ^ j

(5)

Полные компоненты силы, действующей на /-н элемент, получаются суммированием величин (5) по всему окружению, в направлении у дополнительно учитывается — тд

Достаточно естественно записывается и взаимодействие элемента со стенками и дном лотка. Предположим, закон движения стенок лотка известен: левая граница движется по закону 01 = 01(1), правая — по за кону 0р = 01 + Ь, где Ь — постоянная длина лотка. Если расстояние (х^ — 01) меньше г (рис. 3, о), то возникает горизонтальная сила

Гх = (хг — 01 — г)

тд

В случае "проникновения" элемента в правую стенку (рис. 3, в)

т? ( п х \тУ

гх = -(-Г - Ор + Хг)-.

£

Если координата у,1 становится меньше г (рис 3, б), возникает вертикальная компонента

=(г — уО

тд

х = С1(()

у = О

х = вр

Рис. 3. Схема взаимодействия частиц со стенками и дном лотка

£

При необходимости можно учесть и вертикальное движение лотка по известному закону.

После описанных действий (для каждого элемента окружения) суммарные силы, действующие на элемент с номером /. вычислены.

Далее необходимо численно решать задачу Коши для уравнений Ньютона, Вопрос о выборе численного алгоритма для решения данной задачи не тривиален, что требует пояснения. Дело в том, что даже линейная задача об упругом взаимодействии тяжелых точек (в простейшем случае шарик на пружине) не является асимптотически устойчивой. Спектр оператора перехода в лучшем случае чисто мнимый, а в данной задаче вполне может содержать положительные вещественные части. Кроме того, система уравнений является жесткой. Поэтому применение достаточно простых численных методов типа метода Эйлера первого порядка точности [7] исключается (область устойчивости этого метода вообще не содержит точек мнимой оси), С другой стороны, вычисление сил требует достаточно большого (в плоском случае порядка М2) количества операций, в силу чего использование многостадийных явных методов [8] также не желательно.

Поясним последнее. Хорошо известны эффективные многостадийные методы с автоматическим контролем точности и устойчивости (см., например, [8]), выигрывающие по сравнению с классическими за счет большего шага по времени. Однако в нашей задаче можно провести простую оценку максимального шага, при котором имеет смысл рассчитывать процесс. Частота вибрации лотка в конструкции ^ 10 Гц, ЯМ ТТЛ итудя, ~ 1 За один шаг можно "разрешить продвинуться" стенке лотка на треть-четверть радиуса. Элементарные оценки дают в этом случае шаг по времени не больше 0,002 с,

В качестве компромисса используем схему Рунге — Кутты второго порядка точности [7], по которой можно (с полученным шагом по времени) вычислить устойчивое решение.

Запишем закон Ньютона для /-го элемента в виде системы обыкновенных дифференциальных уравнений первого порядка

ОО 1 V11,

■ (6) VI = Р -

Уг г, (7)

ь2г = Ру -

Здесь и Ру — действующие па г-й элемент суммарные силы в направлениях х и у, ь1г, у2г — компоненты скорости центра масс в направлениях х и у соответственно, а диссипативные (вязкие) члены и где ^ > 0 — коэффициент вязкости, введены искусственным образом для повышения устойчивости решения. Непосредственно при численном эксперименте значение вязкости ^ выбирается чисто эмпирически так, чтобы обеспечивалась устойчивость, но значение ^ было, по возможности, минимальным.

На каждом шаге по времени ДЬ на первой стадии вычисляются промежуточные величины па шаге ДЬ/2. По методу Эйлера

= % + ^Г1 = ^ + "

уГ1 = Уг + ^2» ^Г1 = ь2пг + - ^2?). (8)

На второй стадии совершается переход на следующий шаг по времени Ар х?+1 = х? + АгйГ', г>1Г1 = VI? + А* - ,

уГ1 = У? + д, = < + Д* - ) • (9)

Здесь Рх+2 и Ру + 2 — силы, вычисленные для и

Начальные условия можно выбрать достаточно произвольно. Проще всего взять ровные вертикальные "столбцы" из элементов, где верхний ряд "вдавлен" во второй на величину £, второй в третий — на 2£ и т. д. Это дает возможность легко вычислить вертикальные координаты гранул в начальный момент времени. Если лоток неподвижен, такая система теоретически должна находиться в равновесии все время, что является одним из самых важных тестов устойчивости численного решения динамической задачи, Когда лоток приходит в движение, элементы среды за несколько шагов по времени принимают "плотную упаковку", обеспечивающую минимум потенциальной энергии системы.

Учет сил трения в рассматриваемой задаче — наиболее сложный вопрос в моделировании процесса движения гранулированной среды. Ясно, что нужно различать "сухое" трение, когда движение соседних частиц и стенок сосуда происходит без проскальзывания, и трение скольжения (максимально возможное), возникающее при "проскальзывании" соседних элементов.

Будем считать, что между элементами в процессе движения возникает сила трения, пропорциональная силе упругого сдавливания и направленная против относительной скорости движения точки контакта.

Рассмотрим два "соседних" элемента с номерами г и у (рис, 4), Мы уже выяснили, что сила давления элемента у на элемент г

и направлена против вектора е1 (3), Возникающая при этом сила трения Д] равна

где V — коэффициент трения скольжения, и в зависимости от проекций скоростей точек А и А (см. рис, 4) па линию (1), идущую вдоль вектора е2, направлена либо в направлении вектора е2, либо в противоположном. Для точки А , принадлежащей элементу у,

Vа' = V + ^ х ОА = ^и, 1*2^ + х ОА.

Рис. 4. Схема взаимодействия двух элементов с учетом сил трения

Для точки А, принадлежащей элементу г,

^Л = ^ + Ш х О А = ) + Ш3 х О А.

Обозначим направляющие косинусы вектора в1 как

ех „ , еу .

кг з кг з

Тогда -в2 = {ву, -вх).

Проекция скоростей точек А и А есть их скалярное произведение на вектор — в2:

РЯ\Л |(1) = {У1г,у2г){ву, —вх) + ШгК = у1гву — у2^вх + ШгГ,

РК\А' |(2) = {ь1з,ь2]){ву, —вх) + Ш3К = VIзву — ь2звх + Ш3г.

Здесь Шг и Шз — угловые скорости вращения г-го и ^'-го элементов. Относительная скорость вдоль линии (1)

W = РК\А !(1) — РК\А' !(!) = — VIз)ву — — v2з)вх + {Шг — Шз)г. (10)

С точки зрения программирования удобнее всего ввести величину {гп)гз, которая принимает значения {гп)гз = 1, если W < 0, и {гп)гз = — 1, если W > 0, В этом случае к силам (5) необходимо добавить соответствующие компоненты сил трения

Рхз !х3 + Ргз вх{гп)гз,

Ру3 = 1-у3 — ву{гп)гз,

Мг = —М0 • {т)гз. (11)

Взаимодействие с дном и боковыми стенками записывается подобным образом: если 0 < уг < г, то {гп)г = 1 при {у^ — из — Шз г) < 0 и {гп)г = — 1 при {vN — из — Шз г) > 0, Здесь vN = С¡{Ь) — скорость движения дна лотка. Далее силы вычисляются следующим образом:

Рх3 !х3 + Ртр {гп)г,

Р г3 = / г3

у

Мг = —М0 • {гп)г. (12)

Учет трения обязательно приводит к учету вращения элемента, поэтому в уравнениях (11) и (12) необходимо вычислять также и момент силы.

После того как суммарные (с учетом трения) силы подсчитаны, численное решение уравнений проводится тем же методом (8), (9), Однако размерность системы возрастает, кроме системы четырех уравнений (6) и (7) необходимо интегрировать уравнение, описывающее вращение элемента под номером г вокруг оси, проходящей через центр тяжести:

Ршъ = Мг. (13)

тг2

Здесь ,] = —---момент инерции круга относительно оси, проходящей через центр.

Решение уравнений (12) производится численным алгоритмом (8), (9), Записывать полное уравнение вращения не обязательно, так как силы и моменты сил, действующие на рассматриваемый элемент, зависят только от координат центра круга и угловой скорости.

2. Вычислительные эксперименты

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

На рис, 5 приведено численное решение задачи для удаленного от начального момента времени, когда система уже достаточно плотно упакована. Здесь рассмотрено 2400 элементов в подвижном лотке и порядка 100 в бункере в начальный момент времени, В дне лотка имеются три отверстия размером в три диаметра элемента, через которые происходит высыпание частиц. Соответственно опустошается и бункер, В численном эксперименте длина лотка принята 50 см, диаметр частиц — 0,8 см. Величина е составляет 0,01 радиуса частицы. Частота горизонтальных гармонических колебаний лотка 9 Гц, амплитуда 1 см, шаг по времени 0,002 с. Около отверстий расположены счетчики, показывающие количество выпавших частиц.

На рис, 6 светлые гранулы образуют зоны уплотнения, темные — зоны разрыхления среды.

Рис. 5. Результат работы программы

гш* а

Рис. 6. Зоны уплотнения

Предложенный подход применялся для моделирования поведения гранулированной среды в конструкциях роторного типа. Была проведена оптимизация конструкций

и вибрационных режимов их работы. Получены результаты для трехмерного варианта

задачи с гранулами в виде шаров.

Список литературы

[1] Вишняков A.A. Универсальное почвообрабатывающе-посевные машины. Красноярск: Изд-во КрасГАУ, 2004. 202 с.

[2] Вогульская H.A., Вогульский И.О., Вишняков A.A. Оптимизация конструкционных параметров высевающего аппарата вибрационного типа // Вестник КрасГАУ. 2009. № 1(28). С. 110-112.

[3] Пат. № 2192111 Россия, МКИ А 01 С7/16. Вибрационный высевающий аппарат овощной сеялки / A.A. Вишняков, A.C. Вишняков, A.A. Вишняков, П.Г. Ляшок, И.О. Вогульский. Опубл. в ВН., 10.11.2002, № 31.

[4] Пат. № 2310311 Россия, МПК А 01 С7/16. Высевающий аппарат сеялки / A.A. Вишняков, A.C. Вишняков, И.О. Вогульский, H.A. Вогульская, Д.А. Каркавин, В.А. Козлов. Опубл. в ВН., 20.11.2007, № 32.

[5] Вогульская H.A., Вогульский И.О., Вишняков A.A. Вертикальная вибрация жидкости в сосуде // Ресурсосберегающие технологии: Сб. науч. тр. / Прил. к "Вестнику КрасГАУ". 2005'№ 3. С. 115-120.

[6] Никифоров А.Л., Вишняков A.A., Вогульский О.И. Движение сухих семян в вибрирующем лотке сеялки // Международная научно-практ. конф. "Земледельческая механика в растениеводстве". \!.. 2001. С. 106-110.

[7] Калиткин H.H. Численные методы. М.: Наука, 1978. 512 с.

[8] Новиков Е.А. Явные методы для жестких систем. Новосибирск: Наука, 1997. 195 с.

Поступила в редакцию 15 августа 2010 г., с доработки — 22 октября 2010 г.