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

где р = const, Р >0, значение Р определяется с помощью испытаний; t- время работы конструкции (до первого отказа).

Используя (11), (12), находим время tkq эксплуатации конструкции Vq (имеющей коэффициент запаса прочности пкг) tkr = Р пкг ИЛИ tkr = Рпо (1 - pkq )1(\-ро).

Литература

1. Афанасьев Н.Н. Статистическая теория усталостной прочности металлов. - М.: Изд-во Акад. наук УкрССР, 1953.

2. Волков С.Д. Статистическая теория прочности. - М.: Машгиз, 1960.

3. Вейбулл В. Усталостные испытания и анализ их результатов. - М.: Машгиз, 1964.

4. Болотин В.В. Статистические методы в строительной механике. - М.: Госстроиздат, 1965.

5. Когаев В.П. Расчеты на прочность при напряжениях переменных во времени. - М.: Машиностроение, 1993.

6. Ломакин В.А. Статистические задачи механики твердых деформируемых тел. - М.: Наука, 1970.

7. Гмурман В.Е. Теория вероятностей и математическая статистика. - М.: Высш. шк., 2003.

8. Болотин В.В. Ресурс машин и конструкций. - М.: Машиностроение, 1990.

9. Москвичев В.В. Основы конструкционной прочности технических систем и инженерных сооружений. -Новосибирск: Наука, 2002.

10. Немировский Ю.В., Янковский А.П. Рациональное проектирование армированных конструкций. - Новосибирск: Наука, 2002.

11. Матвеев А.Д. Определение коэффициентов запаса прочности конструкций с учетом распределения напряжений // Вестн. КрасГАУ. - 2007. - №3. - С. 159-168.

УДК 631.331.62-66 Н.А. Богульская, И.О. Богульский

МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ НЕОДНОРОДНОЙ ГРАНУЛИРОВАННОЙ СРЕДЫ С ГРАНУЛАМИ НЕПРАВИЛЬНОЙ ФОРМЫ

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

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

N.A. Bogutskaya, I.O. Bogulskyi MOTION MODELING INHOMOGENEOUS GRANULAR MEDIUM WITH GRANULES OF INCORRECT FORMS

The mathematical model of granular medium motion is presented in the article. Granules are represented as ensembles of regular absolutely rigid bodies with elastic membranes, making it possible to describe their interactions and movement at a sufficiently large time interval.

Key words: mathematical modeling, granular medium, optimization.

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

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

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

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

cs = mg. (1)

Здесь с - жесткость упругой оболочки.

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

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

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

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

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

Сначала рассмотрим простейший вариант, не учитывающий силы трения.

Пусть сосуд содержит N элементов, координаты центров которых х., У, 1=1, 2, 3, ..., N известны.

Вычисление сил начнем с определения элементов, «ближайших» к /-му. Можно сразу вычислить расстояние между центрами элементов / и ] Ф /, = ^(х. — х. )2 + (у. — у. )2 и отобрать в качестве

«ближайших» те, для которых Я < с1. Однако более выгодной, с точки зрения экономии вычислений, яв-

ляется следующая процедура: вычислим модуль разности координат

,]=1, 2, 3, А/, у Ф г. Если

эта величина больше ^ элемент с номером ] больше не рассматривается. В противном случае вычисляются

х.-х.

г 3

величины

у,-У,

у,-У,

< (I, вычисляем величину К и в

, ¡=1, 2, 3, А/, j Фi. Только в случае

случае Я < считаем элемент с номером\ лежащим в «ближайшем окружении».

* ]

Таких элементов не может быть много. Мы запоминаем их номера и далее работаем только с ними. Вычислим силу, действующую на /-й элемент со стороны )-го.

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

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

Г \

Х.—Х. у.-у.

Єї =

(3)

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

Є2

г \

у - у х.—х.

^ ^ I ] I

Я ’Я

V ч ч У

(4)

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

г] =-(с!-Я ХА (5)

V ' 2е Я

и

/г7 =-(с1-Я )тё^у1~у \. у 4 2е Я

ч

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

Достаточно естественно записывается и взаимодействие элемента со стенками и дном лотка. Предположим, закон движения стенок лотка известен: левая граница движется по закону 01=01(1), а правая 6р=6/+/_, где /. - постоянная длина лотка.

Если расстояние (х. — 01) становится меньше г, возникает горизонтальная сила

^ =(Х1-а-г)^.

Х ' £

В случае «проникновения» элемента в правую стенку - сила

77 / Г1 і Л т&

є

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

77 / лmg

Ру=(г-У)----------■

є

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

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

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

Поэтому применение достаточно простых численных методов, типа метода Эйлера первого порядка точности [3], исключается (область устойчивости этого метода вообще не содержит точек мнимой оси). С

другой стороны, вычисление сил требует достаточно большого (порядка N2 в плоском случае) количества операций, поэтому использование многостадийных явных методов [4] тоже не очень желательно.

Это тоже требует пояснения. Хорошо известны эффективные многостадийные методы с автоматическим контролем точности и устойчивости (к примеру, [4]), которые выигрывают по сравнением с классическими за счет большего шага по времени. Однако в нашей задаче можно провести простую оценку максимального шага, при котором имеет смысл рассчитывать процесс. Частота вибрации лотка в конструкции » 10Гц, амплитуда «1мм. За один шаг можно «разрешить продвинуться» стенке лотка на треть-четверть радиуса. Элементарные оценки дают в этом случае шаг по времени не больше, чем 0.002 сек.

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

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

г х = и.,

: ' (6) и = р — ц и .

г х I

У = V.,

■к (7)

V = Ь - ц V .

г у ' г

Здесь Г и Г - действующие на /-й элемент суммарные силы в направлениях х и у, и., V. - компоненты скорости центра масс в направлениях х и у соответственно, а диссипативные (вязкие) члены [Л и. и //V , где ¡л>0 - коэффициент вязкости, введены искусственным образом для повышения устойчивости решения. Непосредственно при численном эксперименте значение вязкости /л выбирается чисто эмпирически, и таким образом, чтобы обеспечивалась устойчивость, но значение ¡л было бы, по возможности, минимальным.

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

Хп+1 =хп+— и\

г г 2

йг=й;+(8)

уГ=у:+^у%

уй+1 = Vй + —Vй),.

г г 2 1

На второй стадии совершается переход на следующий шаг по времени А?:

1

хп+1=хп+А ги+~\

I I I 5

и"+1 = и" + - ¡и ), (9)

1

уГ=у:+

1 1

' ПА— ПА—

у;+ =у;+А№ 2-//уг 2).

11 11

~ П-\— ~ ПА— "ПЛ— ^пл—

Здесь ^ 2 и Р 2 -силы, вычисленные для значений X. 2 и у. 2.

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

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

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

Рассмотрим два «соседних» элемента с номерами «/» и «у». Мы уже выяснили, что сила давления элемента «у» на элемент «/» равна

Р ,=(с1-Я г; ^ 1}) 28

и направлена против вектора в\, (3). Возникающая при этом сила трения Ртр равна

* ]

Ртр = ур ,

ч Ч '

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

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

у а = у, + со, ■ОА = (м.,у.) + со1 •ОА. Для точки А, принадлежащей элементу «/»

у а = V, + О] • О А = (и у ) + со] - О А.

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

х -х У --У-

і і У і У і

е = —-----, е =—----------

Я у Я

ч Ч

Тогда -е2=(еу,-ех).

Проекция скоростей точек А и А есть их скалярное произведение на вектор — вг: Пр Уа пл=(и ..у )(е -е ) + (Х>Я = ие -у е +0)Г,

± (\) \ г? гу? х у г 1 у гх

ПрУл =(и ,у )(е-е) + 6)Я = и е -у е +СО Г.

г (1) ^ J, J■'к■ У ’ ] ] У ] X ]

Здесь ш. и й)] - угловые скорости вращения /-го и 7-го элементов. Относительная скорость вдоль

линии (1)

w = Пру А т-Пр Va т = 0иг -Ui)ev - (v. - У.)ех + (щ ~ СО .)г

(1) V Z j / у V i у/ X V i j.

(10)

С точки зрения программирования, удобнее всего ввести величину (т)и , которая принимает значения — 1 , если СО < 0 И = — 1 если со > 0.

В этом случае к силам (5) необходимо добавить соответствующие компоненты сил трения

FiJ = fi¡ + Fmpe (in) ,

x Jx IJ y\ Jij'

= fu -Fmpe (in) ,

y Jy 4 yy

M =-M0 -(in). ..

l l v Jij

(11)

Взаимодействие с дном и боковыми стенками также записывается подобным образом: если (О < у. < 2), то (ш). = 1 если (Уд, - и} - сот) < 0 и (ш)г = -1, если (у^ - и} - сот) > 0.

Здесь Уд, - скорость движения дна лотка.

' N

После этого

Fij - fij + F (in).,

x J x mp v si*

pij _ fij y J X 1

M = -M° ■ (in)i.

(12)

Учет трения обязательно приводит к учету вращения элемента, поэтому в уравнениях (11) и (12) необходимо вычислять еще и момент. Здесь М{). = • г, если элемент «проникает» в левую границу

(ш)г =

1, если уг - (Ог < О, -1, если У -cor > 0.

1 1

Тогда

F'1 = flJ +F8,

х J х 5

¥ч = fü +Fmp(in)

У J У ¡ J v y¡ ’

М = -М° • (ш)г

Для случая правой границы:

(ш)г =

1, если У — СОТ < О,

1 1

-1, если у -cor > 0.

1 1

FIJ = fij +F8,

X J X 5

/•' =/ + Fmp (in).,

y J y ' J v y¡ ’

М = -М° • (ш)г

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

J сог = М.. (13)

ш г2

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

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

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

Следующий вычислительный эксперимент проводился для того, чтобы выбрать оптимальный вариант расположения отверстий в дне вибрирующего сосуда, обеспечивающий равномерность истечения гранул через отверстия. Будем считать часть сыпучей среды уплотненной, если частицы касаются всех «соседей» и разрыхленной в противном случае. Для равномерности выпадения необходимо, чтобы отверстия находились в зоне разрыхления.

Применяем созданную имитационную модель для определения зон разрыхления и уплотнения. При вибрации стенки начинают колебаться и создают зоны уплотнения. Свободное истечение частиц в зонах уплотнения нарушается.

На рисунке 1 - расположение зон разрыхления и уплотнения. Уплотнение материала в данном случае вызвано: у левой и правой стенок лотка - колебаниями, у горловины бункера - давлением столба материала, находящегося в бункере. Отверстия в дне для данной конструкции расположены в зоне разрыхления. В этом случае свободное истечение частиц через отверстия не нарушается.

Рис. 1. Расположение зон разрыхления; светлыми изображены частицы, расположенные в зоне уплотнения материала, темными - в зоне разрыхления 3. Представление гранул в виде ансамблей шаров

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

лительные затраты. В работе предлагается представлять гранулы неправильной формы в виде ансамбля шаров (рис.2).

Рис. 2. Модели гранул

Был выполнен вычислительный эксперимент, в котором рассматривался упрощенный вариант удлиненного тела - «диполь»: каждая гранула представлялась в виде двух шаров, жестко соединенных стержнем (рис. 3). При умеренной протяженности пары взаимодействие таких элементов вычисляется не сложнее, чем взаимодействие отдельных шаров и, в то же время, учитывается ориентация элементов.

Рис. 3. Модель гранулы в виде «диполя»

Предположим, координаты центра гранулы Оі = {хі,уі). Тогда координаты шаров, ее составляю-

щих

хи1=х,+1со^щ, уи1=у, + Нтщ, хи2=х,-1со^щ, у2=у-Ытд?г

Пусть т - масса, И - радиус, момент инерции ,1 = т(И2 + 2/2 ). Рассмотрим взаимодействие гранул. Обозначим

-\/(А /.а хі ./) +(.У^,к У и) .^ — 1^2 (У],к~Уи)

"'і, ] ,к, I

Силы взаимодействия могут быть вычислены по формулам

/л = -

/'.і'./.А'

2

моменты сил

C(ci-Riilk)

К =—ор С

/. /. I.k

c((/-Riilk)

Fv=-------(yj.k -Уи)’

2Ru.,.k

M

(X./ -2.V, + xJM)(г/А. -yu)-(yu -2у + yJk)(xJk -,v(/)

Уравнения движения с учетом вязкости Рэлея

2mxi=Fx-2ju\x, 2myi=Fy-2Juvy, Jqjj = М -2/л!2ф,

где Ff = -2juvx, FK = -2juv Мк = -2ц11а

?R

X ’ у

k

Рис. 4. Результат работы программы Вычислительный эксперимент подтверждает работоспособность представленной модели.

Литература

1. Богульская Н.А., Богульский И.О., Вишняков А.А. Имитационный подход к моделированию движения гранулированных сред // Вестн. КрасГАУ. - Красноярск, 2005. - № 9. - С. 214-218.

2. Богульский И.О., Богульская Н.А. Численное моделирование движения гранулированной среды в

подвижных сосудах // Вычислительные технологии. - 2011. - №2.

3. Калиткин Н.Н. Численные методы. - М.: Наука, 1978. - 512 с.

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