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

УДК 51.72, 539.3

Численный метод построения континуальной модели деформирования твердого тела, эквивалентной заданной модели дискретных элементов

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

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

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

Ключевые слова: континуальная модель, сплошная среда, гранулированный материал, аффинное преобразование, численный анализ, метод молекулярной динамики, метод дискретных элементов

Numerical method for constructing a continual deformation model equivalent to a specified discrete element model

A.F. Revuzhenko and S.V. Klishin

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

A class of boundary value problems of solid and granulated solid mechanics is identified in which uniform strain and strain rate distributions are possible. Boundary conditions appropriate to these problems are applied for statements of numerical experiments. An affine transformation method is proposed for constructing a continual model equivalent to a specified discrete element model. An example of using the method to construct a three-dimensional model of granular media is given.

Keywords: continual model, continuous medium, granulated material, affine transformation, numerical analysis, molecular dynamics method, discrete element method

1. Введение

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

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

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

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

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

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

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

С другой стороны, в настоящее время альтернативой континуальным моделям являются модели молекулярной динамики [1, 2]. В последние десятилетия в связи с развитием вычислительной техники эти методы приобрели большую популярность [3-6]. В рамках этого подхода среда сплошной не предполагается, а задаются форма частиц, составляющих среду, их начальная конфигурация (упаковка), начальные скорости и законы взаимодействия частиц друг с другом, а также с внешними границами. Затем решается система уравнений, описывающая либо квазистатическое деформирование ансамбля частиц, либо их движение с учетом сил инерции.

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

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

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

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

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

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

3. Метод аффинных преобразований

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

Известно, что идеальными в качестве базовых являются эксперименты, в которых реализуется однородное распределение деформаций. Рассмотрим данный вопрос подробнее, следуя [7, 8]. Предположим, что среда является сплошной и процесс ее деформирования раз-

вивается в пространстве Ох1 х2х3 и времени /. Пусть и,, V, — компоненты смещений и их скоростей; е^ и в^ — компоненты тензоров деформаций и скоростей деформаций:

du:

Vi =-

dt

1

= — 2

"ч

du dUj

dx.

eij т

dv г

dxi dv■

1 duk duk

2 dx dx,-

(1)

dx,■ dxг

Здесь и ниже все индексы принимают значения 1, 2, 3, по немым индексам проводится суммирование. Уравнения движения и определяющие уравнения запишем в следующем виде:

доу' ^ dvг

—- + рГ = р—,

дх,- dt

= Ф</ [e

ijV^kb akl'

(2) (3)

где о у — компоненты тензора напряжений; р — плотность; Г, — компоненты массовой силы; фу — заданные функции.

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

и, = и, (х1, х2, х3, /), vi = VI (х1, х2, х3, /), (4)

=°ii (^ x2 ' x3' t).

(5)

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

Действительно, рассмотрим распределение скоростей, линейное по координатам: dx■

V, = — = аЛ (Г) х + а 12 (Г)х2 + аг3 (Г)х3. (6)

dt

Коэффициенты ау заданы как функции времени. Пусть в исходном состоянии тело ограничено замкнутой поверхностью S0. Зададим на скорости (6). Тогда при выполнении некоторых ограничений можно утверждать, что и внутри области распределение скоростей будет иметь вид (6).

Вначале рассмотрим необходимые условия. Предположим, что решение для скоростей имеет вид (6). Из (1) видно, что распределение деформаций и их скоростей будет однородным, т.е. компоненты е^ и ву зависят от времени, но не зависят от координат хк.

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

Обратимся к уравнениям движения. Так как напряжения от координат xk не зависят, то уравнения (2) примут следующий вид:

77 dv

PF =Р-Г7. dt

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

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

Теорема об аффинном деформировании. Если деформируемое тело: 1) является сплошным и однородным; 2) нагружается при отсутствии массовых и инерционных сил (т.е. квазистатически) и, кроме того, 3) является устойчивым (имеет место единственность решений), то задание аффинного преобразования границы тела гарантирует, что и все его внутренние элементы будут испытывать аффинные деформации независимо от реологии тела.

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

Основная идея настоящей работы состоит в том, чтобы базовые физические эксперименты заменить на численные, сохранив при этом всю общепринятую методологию построения континуальных моделей. Сделаем сопоставление экспериментальной (индекс exp) и численной (индекс cal) методик построения континуальных моделей.

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

Шаг IIexp: определяем класс нагружений, который должен быть реализован. Данный класс зависит от возможностей оборудования и поставленных задач. Шаг IIcal: определяем параметры нагружения и конкретный

тип аффинной деформации, который должен быть реализован, т.е. выбираем функции а^ (/) в преобразовании (6).

Шаг Шехр: изготавливаем образцы материала необходимой формы. Шаг Шса1: задаем исходную конфигурацию расчетной области деформирования.

Шаг 1^хр: проводим серии испытаний. Шаг ГУса1: проводим серии численных расчетов. Проверяем выполнение критерия существования континуальной модели (рассмотрен ниже). Если критерий не выполняется, то признается, что эквивалентной континуальной модели не существует. В случае выполнения критерия переходим к шагу ^а1.

Шаг ^хр: обрабатываем экспериментальные результаты с целью формулировки математической модели среды и определения значений параметров, которые фигурируют в данной модели. Шаг ^а1: применяем процедуру осреднения и определяем связи между макрохарактеристиками среды. На данной основе строим математическую модель и вычисляем макропараметры среды, которые фигурируют в модели.

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

4. Критерий существования эквивалентной континуальной модели

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

Нетрудно ввести количественные характеристики для допустимых и недопустимых отклонений ломаной от прямой. Однако и без этого можно оценить ситуацию в целом визуально (пример приведен на рис. 2 работы [9]).

Если предположить, что для заданной дискретной модели эквивалентная ей континуальная модель существует, то указанный выше критерий отсюда последует. Это означает, что выполнение критерия — это только необходимое условие эквивалентности. Достаточность критерия будем проверять на конкретных примерах.

Доказательство достаточности в общем случае, по-видимому, возможно, но здесь не рассматривается.

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

5. Пример использования метода аффинных преобразований

Остановимся на модели деформирования сыпучей среды на основе метода дискретных элементов в трехмерной постановке. Будем действовать последовательно согласно изложенному выше алгоритму. Суть метода дискретных элементов применительно к сформулированной задаче заключается в дискретизации исследуемой конечной области на отдельные элементы (частицы), заданием их механических свойств и геометрических параметров. На основе этих условий численно интегрируется система уравнений движения частиц. В данной работе рассматривается сухая несвязная среда, причем силы взаимодействия между отдельными частицами возникают только при их контакте. Частицы при деформировании среды не меняют своей формы и не разрушаются. Форма частиц ограничена сферами с заданным распределением их радиусов.

5.1. Шаг 1са1. Начальная упаковка

Создание начальной плотной упаковки частиц (сфер) с заданным распределением размеров является одной из основных задач на начальном этапе моделирования методом дискретных элементов. Обычно предполагается, что частицы или близки по размеру (типа монодисперсных сфер) или, наоборот, сильно различаются (мелкие сферы вложены в полости укладки крупных). В настоящей работе выбраны близкие по размеру частицы, а именно: радиусы сфер выбираются случайным образом на основе равномерного распределения и отличаются друг от друга в 1-1.25 раза. Это позволяет применить алгоритм, при котором на первом этапе задается начальное распределение частиц по координатам и размерам, а на втором имитируется рост их диаметров с учетом контактного взаимодействия (алгоритм последовательного роста частиц). Таким образом, начальная упаковка «разбухает», заполняя собой исследуемую область, причем конечная плотность упаковки (суммарная доля

объемов частиц по отношению к объему всей области) составляет порядка 0.65-0.75 [10, 11].

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

Силовое взаимодействие в точке контакта между сферическими частицами с номерами - и ] можно представить в векторной форме следующим образом:

р.. = Е"п- + ЕхХ-

где Е/ и Е/ — нормальная и касательная компоненты контактной силы; П/ — единичный вектор, определяющий плоскость контакта между двумя сферами, а именно вектор нормали к плоскости пересечения сфер; X/ — единичный вектор, принадлежащий плоскости контакта.

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

Рассмотрим две составляющие перекрытия — по нормали 5П и по касательной 5/. Первая составляющая определяется по формуле

Щ = г + г/ -1/ > 0

где г и Г/ — радиусы - -й и/-й частиц соответственно; I/ — расстояние между их центрами. В приращениях по времени Дt выражение для Д5П формулируется следующим образом:

Д5П =-Д/.

Приращение касательной составляющей определяется скалярным произведением по формуле

Д5/ = (Ди,-Ди /) • Ц,

где вектор X/ направлен вдоль вектора относительного

смещения частиц в плоскости контакта; Ди,- и Ди, —

1 ]

приращения смещений частиц.

Выберем упругое взаимодействие частиц на основе закона Герца. Тогда компоненты контактной силы будут определяться следующим образом:

еп=ф (5/ ^

(7)

(8)

где

1 -V? 1 -V

Е,-

2Щ

(2 + )(1 -V2 + (? + V;)(1 ^;)

Е

Е,

— = - + — Г г; '

Здесь Е, и Е; — модули упругости частиц , и;; V, и V; — коэффициенты Пуассона --й и ;-й частиц соответственно.

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

На первом шаге нагружения для определенного контакта определяются значения Е/ и Е/ в соответствии с (7), (8) и проверяется выполнение неравенства | Е/1 > tg ф | Е/1, где ф — заранее заданный угол внешнего трения скольжения между частицами. Если это неравенство не выполняется, то значение силы Е/ остается определенным в соответствии с (8). Если неравенство выполнено, то частицы начинают проскальзывать друг по другу и на текущем шаге нагружения касательная составляющая силы отталкивания вычисляется по закону:

Е/ = tg Ф Е/. (9)

После определения значений нормальной и касательной компонент силы, действующих на каждом контакте, происходит переход к следующему шагу нагру-жения и данная схема повторяется для существующих и вновь образовавшихся контактов. Если в какой-то момент нагружения для нормальной составляющей перекрытия выполняется условие 5П < 0, то такой контакт считается исчерпанным и в дальнейшем не рассматривается. Шаг по времени Дt выбирается достаточно малым, чтобы число новых контактов, перешедших в предельное состояние на очередном шаге нагружения, было относительно небольшим.

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

1

5.2. Шаг 11са1. Выбор типа аффинной деформации

Для систематического исследования необходимо рассмотреть серии трехосных нагружений без поворота осей тензора деформаций. Затем перейти к сложным нагружениям с изломами траекторий, а также к на-гружениям с непрерывным поворотом осей тензора деформаций. В случае необходимости сюда включаются и различные знакопеременные пути нагружений. Выбор пути нагружения сводится к выбору матрицы аффинного преобразования (6). Таким образом, речь идет о достаточно обширной программе исследований. Ее реализация в полном объеме выходит за рамки одной статьи.

Ниже ограничимся только одним наиболее интересным классом нагружений, а именно сложным на-гружением с непрерывным поворотом осей тензора деформаций [7, 8]. Указанный тип деформирования представляет собой суперпозицию течений Куэтта между параллельными пластинами. Это течение обладает замечательным свойством: если на границах реализовы-вать постоянные скорости, то внутри может возникать только линейное распределение скоростей независимо от реологии среды. Такой деформации отвечает простой сдвиг со скоростью у. Предположим, что сдвиг реализуется в течение времени At. Затем тело поворачивается на угол юДt, снова реализуется простой сдвиг и т.д. Таким образом, мы приходим к сложному нагружению, которое складывается из простых сдвигов, направления которых меняются с постоянной скоростью ю. В течении Куэтта, если оно осуществляется неограниченно долго, траектории точек (вне оси симметрии) уходят на

бесконечность. В рассматриваемой задаче за счет поворота (если его скорость больше скорости сдвига w > g) траектории материальных точек остаются в ограниченной области пространства. Последнее обстоятельство открывает путь для экспериментальной реализации подобных течений. В представленной работе ограничимся численным экспериментом. Рассмотрим вначале кинематику деформирования для плоского случая.

Суммирование деформаций сдвига показывает, что движение материальной точки в плоскости Ox1x2 описывается следующей системой уравнений:

dx dx2

V = —1 = -(w-g)X2, v2 = wxi, (10)

dt dt

где v1, v2 — компоненты вектора скорости; w и у —

константы. При w > у решение (10) имеет вид:

x1 (t) = x° cos (— t) -—x° sin (— t), w

w (11)

x2(t) =—xfsinG— t) + x° cos(—t), Ц_

где — = yjw (w - g) и начальные координаты

x1(0) = x0, x2(0) = x0. (12)

Материальная точка (12) движется по эллиптической траектории

x2 (t) + ^ x2(t) = (x0)2 (x20)2. (13)

w w

На рис. 1, а показана траектория точки с начальными координатами

x1(0) = a, x2(0) = 0. (14)

Рис. 1. Схема численного эксперимента (а); расположение частиц в начальный момент времени (б) и после 1.75 полных оборотов частиц на границе (в)

5.3. Шаг IIIcaj. Задание исходной конфигурации области деформирования.

Выбираем в качестве границы исследуемого тела траекторию (13). Полуоси эллипса равны a и ^ю/(ю-g) a. Закон движения точки по эллиптической «орбите» обладает замечательной особенностью: за равное время радиус-вектор r точки, движущейся со скоростью v, «заметает» одинаковые площади:

| vxr | = wa2 + (w-g)b2 = const. (15)

Таким образом, в данном случае реализуется закон Кеплера. В частности, линейные скорости в точках A и B (рис. 1, а) равны соответственно

Vb Ю-g

vA = ю а, =- (ю-g) , — =---.

vA ю

5.4. Шаг IVcat. Численные расчеты методом дискретных элементов

Пусть в пространстве OXj Х2 Х3 задана область — цилиндр, ось которого ориентирована вдоль оси Ox3 и сечением которого в горизонтальной плоскости Oxj x2 является эллипс с соотношением полуосей 0.8 (рис. 1, а).

Здесь необходимо отметить следующее. Очевидно, что в классической постановке метода дискретных элементов существует большая зависимость кинематики движения среды как от начальной упаковки частиц, так и от упругих постоянных. Задавая движение точек границы по закону (11), весьма высока вероятность получения либо разрыхленных участков (где отсутствует взаимодействие частиц), либо участков с предельными силовыми значениями Fn и В рассматриваемой задаче вектор граничной скорости направлен вдоль границы. Для вязкого трения этого условия было бы достаточно для того, чтобы среда и внутри области вовлекалась в процесс деформирования. В случае же сухого трения это не всегда так. При недостаточных по величине нормальных напряжениях на границе, сухого трения будет недостаточно для реализации движения материала. Проще говоря, на контактах с границей возможно проскальзывание и численные эксперименты это подтверждают. Практически в лабораторных экспериментах [7, 8] данная проблема решалась за счет веса слоя материала, расположенного выше. В представленном численном эксперименте роль веса выполняет постоянная нагрузка, приложенная к верхнему и нижнему торцам образца и действующая вдоль оси Ox3.

Рассмотрим данный вопрос подробнее. Формулы (10)-(15) относятся к плоской деформации, в то время как в настоящей работе рассматривается случай трехмерного деформирования: образец представляет собой прямой эллиптический цилиндр, который заполнен сферическими частицами, образующими определенную упаковку. Поэтому возникает вопрос об условиях на торцах и условиях для третьей компоненты смещения

на боковой поверхности цилиндра. Примем, что на торцах задано следующее распределение скоростей: V и v2 заданы формулами (10), а вертикальная компонента v3 зависит от времени и равна vz (t) = v* на верхнем торце и vz (t) = 0 на нижнем. На боковой поверхности третья (вертикальная) компонента скорости задается линейной по координате z и согласованной с вертикальной компонентой на торцах образца. Такие краевые условия в точности являются аффинными и содержатся в классе (6).

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

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

5.5. Проверка критерия существования континуальной модели

Один из главных вопросов численного моделирования на основе метода дискретных элементов связан с выбором размеров (радиусов) частиц. На рис. 1, в приведена картина деформирования образца после 1.75 оборотов частиц на границе. Размер частиц выбирался таким образом, что отношение длины большой полуоси эллипса к среднему радиусу частицы превышало 100 (всего в модели рассматриваются 180 000 частиц). Видно, что первоначальная прямая AC (рис. 1, а), перейдя в B'D , сохранила свою прямолинейную форму, т.е. при данном наборе частиц критерий существования континуальной модели выполняется. Это означает, что любое материальное волокно конечной длины, которое было прямолинейным в начальный момент времени, осталось прямолинейным и в последующие моменты времени (деформации соответствует аффинное преобразование). Это свойство использовано как простой и наглядный тест для решения вопроса о том, сколько частиц необходимо взять в деформируемой области, чтобы численное решение дискретной задачи соответствовало континуальному решению.

5.6. Шаг ¥са1. Определение компонент деформаций, скоростей деформаций и напряжений на макроуровне

Согласно (1), компоненты тензора скоростей деформаций и поворота равны

дю2

dvi

eii = ^ = 0, e22

д x1

д Хо

— 0,

е12

dv1 dv2 1 ■ + - 2

dxo д х

i J

Q — 1

2

dv2 д х1

dv1 д x

2 J

Y

— — --. 2

Оси тензора скоростей деформаций направлены под углом 45° к осям координат, главные значения скоростей деформаций равны

e =

^ — -

Рассмотрим вопрос о величине самих деформаций. Закон движения точки (12) имеет вид:

х1 (t) = a cos (pt), х2 (t) = — a sin (pt).

p

Материальная точка стартует из положения A (см. рис. 1, а) и достигает положения B за время t *:

cos (p t *) = 0, sin (p t *) = 1, p t * = П,

t * = П 1

2 ^ю (ю-у)

За это время короткая ось эллипса растягивается до длинной оси, а длинная ось сжимается до короткой. Таким образом, за время t* деформации области (если их оценивать по изменению длин осей эллипса) достигают максимального значения.

Определим компоненты нелинейного тензора деформаций Коши-Грина:

<2 , л2"

еп =-

дщ 1

дх0

+

е 22 = '

ди1

ди

дх

V дХ1° J ди1

дщ

дх

V дХ2 J

дх

V^! J ди

V дХ2 J

дх

(16)

е19 =

12 2

ди1 \дх2

ди2 Л

a? J

1

+— 2

ди1 ди1

дщ дщ

дх? дх? дх? дх

Из (11) заключаем, что

и1 = х? (cos(pt) -1) -—х? sin(pt), —

и2 — — х? sin(pt) + х^соэ^t) -1). p

Подставляя (17) в (16), получим

(17)

1

е11 =_

11 2

/

\

V

—Т -1

p2

sin (pt),

У

/

е 22 —

22 2

\

4 -1

V— J

е12 = ~

/

— VP

\

cos (pt),

sin(pt) cos(pt).

При t — t * имеем

еп — -

1

/

—

p

-1

е 22 = '

--1

—

12

= ?.

Таким образом, главные оси тензора деформаций совпали с осями хх и х2 (рис. 1, а). Так как

2 ю ,

2 = - >

ц2 ю-у

то в момент t = t *

2

—

= = 1 е1 "е11 = 2

—2 -1 p2

> Q,

j

= =-1

е 2 — е 22 —

22 2

2

1 2

V — J

< ?,

где и 8 2 — главные деформации.

Таким образом, выше описана кинематика деформирования среды в континуальном приближении. Как отмечалось, особенность краевых условий (11) состоит в том, что, какие бы определяющие уравнения среды мы не брали, кинематика деформирования от этих уравнений не зависит.

Перейдем теперь к процедуре осреднения напряжений и определению компонент тензора напряжений на макроуровне.

Для наглядности рассмотрим плоский случай. Пусть внутри области деформирования задана произвольная площадка длиной Ь и нормалью п. Эта площадка разбивает область на две подобласти Ь+ и Ь~ (рис. 2). На каждую -ю частицу (выделены серым цветом), которую пересекает площадка Ь, со стороны Ь+ действует сила Р, равная сумме векторов ^ и рп. Примем, что напряжения равны

Рис. 2. Определение касательных и нормальных напряжений на площадках

Рис. 3. Изменение высоты образца в зависимости от значений угла внешнего трения между частицами ф = 10° (1), 20° (2), 30° (3), 40° (4) (а); дилатансионные кривые, полученные в лабораторном эксперименте при сложном нагружении (б)

т п = |1 Г ,, о „ =1X ,

L i L ■

(18)

где тп, оп — касательная и нормальная компоненты напряжения соответственно; суммирование происходит по всем частицам, попавшим на площадку L. Эти же рассуждения аналогично переносятся на случай трехмерного пространства. Разница состоит в том, что длина L в выражении (18) заменяется значением площади площадки.

Воспользовавшись (18), найдем напряжения оп, о22 и т12 на фиксированных площадках ЛЛ'С'С и ВВББ (рис. 1, а), ориентированных вдоль осей Ох1 и Ох2 в процессе нагружения. Как показали численные эксперименты, значения напряжений стабилизируются уже после 0.5 оборота частиц на границе. Направление площадок главных напряжений определялось по следующей формуле:

2о12

^(2 к) = -

о - о

(19)

22

Сыпучий материал, как известно, обладает двумя основными свойствами: дилатансией и внутренним трением. Переход от двумерного моделирования к трехмерному позволяет исследовать не только внутреннее трение, но и дилатансию, которая оказывает существенное влияние на напряженно-деформированное состояние сыпучей среды. Для этой цели была проведена серия численных экспериментов по нагружению образца на основе закона движения (11). В каждом из численных экспериментов высота образца являлась одинаковой, изменялся коэффициент трения между частицами. На рис. 3, а представлено изменение высоты образца Н в процессе деформирования для различных значений коэффициента трения (углов внешнего трения ф) между частицами. Видно, что при совершении четверти оборота (а > п/4) дилатансия себя практически исчерпывает.

На рис. 3, б представлены дилатансионные кривые, полученные в лабораторных экспериментах по сложному нагружению сыпучих материалов [12]. Здесь нас

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

При нагружении образца до полутора оборотов (а = 3 п/ 2) процесс останавливался и нагрузка прикладывалась в противоположном направлении. Видно, что вначале высота образца стабилизируется (рис. 4, участок 1 кривой), затем, при изменении направления нагружения, уменьшается до некоторого предела (участок 2) (материал уплотняется), а затем опять увеличивается, достигая своего устойчивого состояния (участок 3).

Рассмотрим теперь, согласно [8], соотношение, описывающее угол внутреннего трения Ф материала:

. лДОц -о??)2 + 4о Ф = агсзт^^-1-^-=

о11 + о22

Здесь в числителе приведен инвариант т = (о2 -0)72 (максимальное касательное напряжение), в знаменателе — о = (о2 + а^/2 (среднее сжимающее напряжение); о1, о2 — главные напряжения. Графики, по-

Рис. 4. Изменение высоты образца при изменении направления вращения точек границы

Рис. 5. Изменение угла внутреннего трения Ф в зависимости от значений угла внешнего трения между частицами ф = 10° (7), 20° (2), 30° (3), 40° (4)

6. Выводы

Показано, что при определенных ограничениях аффинное преобразование границы деформируемого тела гарантирует аффинное преобразование и любых внутренних объемов тела независимо от его реологии.

На этой основе разработан метод аффинных преобразований, который позволяет по заданной модели дискретных элементов построить эквивалентную континуальную модель тела.

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

Работа выполнена при финансовой поддержке РФФИ (грант № 10-05-91 -002-АНФ).

казывающие значения угла Ф, полученные в численном эксперименте для различных углов внешнего трения между частицами, приведены на рис. 5.

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

Из результатов численных экспериментов можно сказать, что через 0.25-0.5 оборота границы напряжения, действующие внутри области, переходят в стационарное состояние. В то же состояние переходит и упаковка частиц. В этом случае деформации сдвига не сопровождаются изменением объема упаковки, т.е. дилатансия отсутствует, поэтому здесь имеет место условие несжимаемости. Тот факт, что напряжения ап, а22, а12 становятся постоянными, означает, что среда находится в предельном состоянии. Таким образом, можно сказать, что рассмотренная на основе метода дискретных элементов модель соответствует следующей континуальной модели:

5ап_+d0ii=0 don+d022=0 dui= 0

дх1 дх2

дх1 дх2

дх1 дх2

VÖ

Oil - а22 )2 + 4aÎ2 = - sin Ф(ап + а22 ),

дщ/дх2 + дщ Iдх1 2а12 Л „ arctg———2-Ц—^ = arctg-± 0.1.

дщ/ дх1 - ди2 / дх2

а11 - а

22

Литература

1. Haile J.M. Molecular Dynamics Simulation: Elementary Methods. -Chichester: Wiley, 1997. - 512 p.

2. RapaportD.C. The Art of Molecular Dynamics Simulation. - Cambridge:

Cambridge University Press, 2004. - 564 p.

3. Механика — от дискретного к сплошному / Под ред. В.М. Фомина. -

Новосибирск: Изд-во СО РАН, 2008. - 344 c.

4. Псахье С.Г., Смолин А.Ю., Стефанов Ю.П., Макаров П.В., Чертов М.А. Моделирование поведения сложных сред на основе совместного использования дискретного и континуального подходов // ПЖТФ. - 2004. - Т. 30. - № 17. - C. 7-13.

5. Псахье С.Г., Остермайер Г.П., Дмитриев А.И., Шилько Е.В., Смолин А.Ю., Коростелев С.Ю. Метод подвижных клеточных автоматов как новое направление дискретной вычислительной механики. I. Теоретическое описание // Физ. мезомех. - 2000. -Т. 3. - № 2. - С. 5-13.

6. Кривцов А.М., Кривцова Н.В. Метод частиц и его использование в

механике деформируемого твердого тела // Дальневосточный математический журнал ДВО РАН. - 2002. - Т. 3. - № 2. - С. 254-276.

7. Ревуженко А.Ф. Механика сыпучей среды. - Новосибирск: ЗАО ИПП «Оффсет», 2003. - 373 с.

8. Revuzhenko A.F. Mechanics of Granular Media. - Berlin-Heidelberg: Springer-Verlag, 2006. - 380 p.

9. Ревуженко А.Ф., Клишин С.В. Об одном методе сопоставления моделей молекулярной динамики с континуальными моделями сплошных сред // Проблемы и достижения прикладной математики и механики. - Новосибирск: Параллель, 2010. - С. 577-585.

10. Kolonko M., Raschdorf S., Wäsch D. A hierarchical approach to simulate the packing density of particle mixtures on a computer // Granular Matter. - 2010. - V. 12. - No. 6. - P. 629-643.

11. Labra C., Onate E. High-density sphere packing for discrete element method simulations // Commun. Numer. Meth. En. - 2009. - V. 25. -No. 7. - P. 837-849.

12. Бобряков А.П., Ревуженко А.Ф. Об одном методе испытания неупругих материалов // Изв. АН СССР. МТТ. - 1990. - № 4. -С.178-182.

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

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

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