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

Моделирование поведения сложных сред на основе комбинированного дискретно-континуального подхода

С.Г. Псахье, А.Ю. Смолин, Ю.П. Стефанов, П.В. Макаров, Е.В. Шилько, М.А. Чертов, Е.П. Евтушенко

Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия

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

1. Введение

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

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

Стремительное развитие вычислительной техники позволило в течение последних десятилетий добиться определенных успехов в применении дискретных численных методов для описания гранулированных и сыпучих сред [1-6]. При этом в большинстве работ используются уравнения движения, характерные для метода частиц [6], а межчастичные силы взаимодействия вычисляются в рамках моделей жестких или мягких сфер (rigid body dynamics). Несмотря на достигнутый прогресс необходимо отметить, что дискретные подходы, основанные на этой концепции, по своей природе ориентированы на моделирование сыпучих материалов. Формализм таких методов, как правило, не позволяет достаточно корректно описывать поведение таких «слабо связанных» гранулированных систем как грунты, для которых необходимо учитывать как «сплош-

© Псахье С.Г., Смолин А.Ю., Стефанов Ю.П., Макаров П.В., Шилько Е.В., Чертов М.А., Евтушенко Е.П., 2003

ность» среды, так и возможность ее дробления с последующим перемешиванием частиц и разуплотнением (дилатансией). Свободным от данного недостатка является численный метод, активно развиваемый в последние годы — метод подвижных клеточных автоматов [7-13] (МСА-метод—method of movable cellular automata). Хотя данный метод и основан на дискретном подходе, однако в его основе лежат уравнения движения, отличные от классических, используемых в методе частиц. Кроме того, в [10] показано, что при стремлении характерного размера автомата к нулю формализм метода МСА позволяет перейти к классическим соотношениям механики сплошной среды. Основным достоинством данного метода является возможность явным образом моделировать как формирование несплошностей различного типа (от генерации отдельных повреждений до распространения магистральных трещин), так и эффекты перемешивания масс. Это делает его уникальным инструментом для изучения отклика геологических сред, где при динамических воздействиях подобные процессы протекают крайне интенсивно.

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

2. Методика расчета

Поскольку в данном подходе используются два прин-

ципиально различных метода, то вначале приведем

краткое изложение каждого из них. Затем подробно рассмотрим вопрос совмещения используемых методов.

2.1. Континуальный метод

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

, j = puui, (1)

уравнение неразрывности р + Риц = уравнение энергии E = ор Ej

(3)

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

— (и,-2 г

1

Р а = 2(ui, j+u л )>

(4)

®ij =- (ui,j - uj,i )>

где — компоненты тензора скорости деформаций Коши; иг- — компоненты вектора перемещения; Юу — компоненты тензора скоростей вращения.

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

Р.. = Р • + Р"

ij ij ij.

(5)

Упругое поведение среды описывается гипоупругим законом:

CTj = - р8 ij + sij,

-¡Dj = 2^fPe -3Рkk8j I,

Dsj

Dt

'■ sij - sik ® jk - sjk ® ik >

P=. V

(6)

(7)

(8) (9)

Здесь <5^ — компоненты тензора напряжений; sij — компоненты девиатора тензора напряжений; Р — среднее давление; К и ц — модули сжатия и сдвига соответственно; р — плотность.

Неупругая деформация определяется в соответствии с заданными законом течения и поверхностью нагружения:

/ (5, е, X) = 0,

ЭФ (10)

е =-------

j

где f — поверхность текучести; X — параметр состояния и Ф — пластический потенциал.

Соотношения, используемые для интегрирования уравнений в континуальной области, эквиваленты конечно-разностным уравнениям, приведенным в [14]. Для того чтобы описать особенности реализации условий сопряжения континуальной и дискретной областей, нам потребуются лишь выражения для численного интегрирования уравнения движения. Согласно [14] уравнение движения в конечно-разностном виде записывается как

Л

Рис. 1. Шаблон построения разностного аналога уравнения движения. Показана нумерация узлов и ячеек

* О _ *о ■

Дьп

2фо

[(о^)П (уП - УП) + (о„Уа(у1 - Уп) +

+()Пн(УП - УП) + (^х*)Пу(УП - Уп) -

- (Оху )П (хП - хП) - (а*У )Л(хП - хП) -

- (Оху )Пи(хП - хП) - (оху )2у(хП - хП)], (11)

у О _ УО + ~-[(оуу )П (х2 - х3 ) + (0уу )П1(хЪ - хП ) +

2Фо

+ ( о уу )Пп( хП - хП) + (Оуу )Пу(хП - хП) -

- (оху )П (у2 - уП ) - (оху )£(узп - уП ) -

- (оху )Пп(Уп - уП) - (Оху )Пу(уП - уП )], хп+п _ хп + х Дьп+1/2

п+1 п . • л^п+1/2

Уо _ Уо + У оД •

одес. х _ хп+п/2. х' _ х2-12 • у _ уп+1/2. у' _ у2-1/2*

ОдесЬ .Хо _ хо ’ хо _ хО ’ у о уо ’ у О_ уО ’

2ф — масса, соответствующая узлу О, которая вычисляется как среднее по окружающим ячейкам:

фо _ 4 х (р а) к

4 к _ 1,1У

где А — площадь ячейки.

Шаблон построения разностного аналога показан на рис. 1.

2.2. Метод подвижных клеточных автоматов

В методе подвижных клеточных автоматов моделируемая среда представляется в виде ансамбля взаимодействующих частиц или клеточных автоматов, имеющих конечный размер. Метод наследует все преимущества классического метода клеточных автоматов, при этом добавляется способность автоматов к пространственному движению и вводятся новые параметры состояния, отнесенные к паре автоматов. Таковыми, в частности, являются перекрытие и сила в паре [7-9]. Для определения критерия переключения состояний вводится параметр межавтоматного перекрытия:

где г — расстояние между центрами соседних элементов; гЦ определяется как г/ _ (dl + dj)/2, dj — размер автомата (рис. 2). Существуют два типа состояний (отношений) пар:

связанные — Н < Н^,

несвязанные — Н > Нтах.

В простейшем случае связанными называют пары автоматов, между которыми имеется химическая связь, а несвязанными — пары, между которыми химическая связь отсутствует. Изменение отношения пары определяется относительным движением автоматов, и среда, образованная такими парами, может рассматриваться как бистабильная. Следуя модели Винера-Розенблюта, распределенная бистабильная активная среда может быть описана уравнением:

дну

Д

_ / (Н) + х С (у, 1к) 1( Нк) +х С ( у, II) I (На),

к Ф у IФI

где С(у, ¡ку!)) — коэффициенты, связанные с переносом параметра перекрытия h от одной пары автоматов к другой; I(Нк(1)) — явная функция Н1 к(у), которая определяет перераспределение Н1 к(1) между парами у, ¡к иjl. Функция состояния пары /(Н) имеет смысл относительной скорости VП автомата у.

В линейном приближении функция I(Н1 к(1)) может быть записана как

1 (Н‘к(у]) _ у(ау, к(у,))

Vn ' ’

Рис. 2. Параметры перекрытия пар: у, ¡к и]1

где ^(ау Л(у,)) определяется взаимным расположением (ориентацией) пар автоматов у, ¡к и ]1; а у к ,) — некоторый параметр взаимной ориентации.

Эволюция среды подвижных клеточных автоматов будет описываться следующими уравнениями движения для трансляционной составляющей:

dt2

_|А+ЛI р +

т

т

+ Х С (І1, 1к)¥(а у к)

к ф у

1 1к .

—р + т

+

X С (l1, у,Ма у у )-17 ру , (12)

,Ф1 т

где р1} — центральная сила парного взаимодействия.

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

ё2 0 у ( Л -у1 \

(!+у

Ji Ji

г +

+ X ^ (1у 1к) Тк +

к Ф у J

+ Х ^ (l1,у1) ду т 11 ■

,Ф1

(13)

Здесь 0у — угол относительного разворота (он является параметром переключения, также как и Н'} для трансляционной составляющей); д1 (у^ — расстояние от центра автомата ¡(у) до точки контакта с автоматом у(г') (плечо момента); ту — тангенциальная сила парного взаимодействия; S(ij, ¡к(у1)) — некоторые коэффициенты, связанные с переносом параметра 0 от одной пары к другой (они аналогичны коэффициентам С(у, ¡к(у1)) в уравнениях для трансляционной составляющей движения). Следует заметить, что уравнения (12) и (13) полностью аналогичны уравнениям движения для многочастичного подхода и в общем случае могут быть записаны как:

т

J¡

ёь2

ё20' ёь2

_ FД +

= Х К

X FІJ

(14)

где Rl(у^ — координаты автоматов; 01 (у^ — углы их поворота; FQ — эффективная сила, действующая на автомат ¡ и обусловленная взаимодействием ее соседей с другими автоматами, F1 _ ру + Т1; ру ная составляющая парной силы взаимодействия; т’ тангенциальная составляющая этой силы; К _

-централь-

у

_ д1 (пу XТ1); единичный вектор пу определяется как пу _ ^у - Ri)/(ду + ду )■

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

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

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

2.3. Совмещение континуалъного и дискретного методов

Из описания представленных методов видно, что в обоих методах уравнения движения могут быть записаны для некоторых масс через действующие на них напряжения или силы. Эти особенности дают основания полагать, что совмещение данных методов может быть осуществлено достаточно корректно. Для этого в настоящей работе моделируемая среда считается состоящей из континуальной и дискретной областей. Между ними определяется некоторая граница сопряжения и полагается, что она принадлежит обеим областям. При этом каждому узлу расчетной сетки, лежащему на границе, ставится в соответствие некоторый автомат (элемент МСА-метода).

Рассмотрим границу сопряжения областей, описываемых различными методами. В общем случае, между двумя соседними узлами сетки, лежащими на границе, может располагаться несколько клеточных автоматов (рис. 4, а). Тогда размер автоматов должен быть кратным шагу сетки и не превышать его: D _ Д/п, где D — размер автомата; Д — шаг сетки; п — целое число. В

Рис. 3. Основные типы функции отклика автоматов

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

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

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

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

Область подвижных клеточных автоматов

Узлы с сопряженными автоматами на границе сопряжения двух методов

Область конечноразностной сетки

Рис. 4. Варианты сопряжения континуальной и дискретной расчетных областей: размер подвижных автоматов меньше шага сетки (а); размер подвижных автоматов равен шагу сетки (б)

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

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

(2Ф X) о = -(^24 + ^3! + а ™ у42 + Ст^у™) +

І ґ~І „І . II . _ІІІ III . -IV 13 ч

+ (ахух24 + ахух31 + ахух42 + ахух13 Л

_ ґ~І І І «.II II . «.III III . «.IV IV ч

(2ф у) о = (ауух24 + ауух31 + ауух42 + ауух13 )

^І _і_ «-II і _пі III . ~І3,Д3ч

- (ахуУ24 + ахуУ31 + ахуУ42 + ахуУ13 )’

(15)

где О^х, ОУу , О^у — компоненты напряжений в соответствующих ячейках; хктп _ хкт ~ х2, У^п _ У^т ~ у2 '; к= I, ..., IV — номер ячейки; п, ш = 1, ..., 4 — номер узла (рис. 1).

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

Независимо от трактовки и используемого похода — конечно-разностного или конечно-элементного — при построении разностного аналога уравнений их форма и смысл слагаемых остаются одинаковыми [14-18]. В левой части стоит произведение массы на компоненту ускорения узла расчетной сетки, а справа—компонента равнодействующей силы, приложенной к узлу:

X рх > т у = X РУ ’

(16)

=1,^

=1,^

F1х _ (-°хх ДУ + 0ху^ ),

Fiy _ (о1ууДх1 - о ^Ду1 )■

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

X

Т7І- /«.I . ^ІІ „II . ^ІІІ,,ІІІ .

- - (ахху24 + ахху31 + ахху42 + ахху13 ) +

і=І,І3

І ґ~І „Д . ~ІІ „ДІ . ^.ІІІ ІІ^ , «.ІЗ^ч + (ахух24 + ахух31 + ахух42 + ахух13 ),

ЕТ7Ї ґ~І ^І . ~ІІ ^ІІ . ^ІІІ^ІІІ . ГСч

РУ = (ауух24 + ауух31 + ауух42 + ауух13 ) -

і=І,І3

^І ,Д _і_ «-II і _підп . ~І3,Д3ч

- (ахуУ24 + ахуУ31 + ахуУ42 + ахуУ13 )•

(17)

ху7 24 ху-У 31 ху7 42 ху

Для всех граничных узлов слагаемые от недостающих ячеек заменяются силами, действующими на соответствующие автоматы со стороны МСА-области. Расчет движения таких узлов осуществляется исходя из уравнения (16), куда в зависимости от геометрии по-

верхности раздела будут входить силы, полученные в МСА-области. Далее, после перемещения автомата, неразрывно связанного с узлом сетки, соответствующее воздействие из области сетки передается во внутреннюю область МСА.

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

Ет-,І ( „III III . „III III \ . рх - (-ахху42 + ахух42 ) +

X

+ (-а х3у13 + а х3х13) + ^.

ру - (а уу1 х42 -а “ у4І)+

+ (аУ3х^ -ах3у13) + ¥'у,

(18)

где ¥'х, ¥’у — компоненты вектора силы, действующего со стороны МСА-области на автомат, лежащий на границе раздела.

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

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

Таким образом, схематически алгоритм совместного расчета выглядит следующим образом.

1. В континуальной области решается система уравнений (1)—(10), осуществляется расчет напряженно-деформированного состояния, рассчитываются скорости и координаты.

2. Перед расчетом скоростей узлов, лежащих на границе раздела, осуществляется вызов подпрограммы, реализующей метод МСА. В нее передаются координаты и скорости совмещенных узлов автоматов с предыдущего шага, а также шаг интегрирования.

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

4. Новые данные о граничных автоматах (узлах) возвращаются в сеточный метод. После этого в нем осу-

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

5. Осуществляется согласование величины нового шага интегрирования по времени.

3. Результаты расчетов

Использованные методы и программы их численной реализации ранее неоднократно применялись авторами для описания динамических процессов, в том числе и распространения упругих волн [11-13, 19-21 ]. Для проверки возможности объединения описанных выше подходов в рамках единой методики и тестирования разработанных алгоритмов были рассмотрены задачи о распространении упругих волн в среде со свободной поверхностью. Одна часть моделируемой среды рассматривалась как сплошная, а другая как дискретная. При этом их механические характеристики полагались одинаковыми. Поскольку таким образом задавалась фактически однородная среда, то граница сопряжения не должна проявляться как граница раздела различных сред.

Был рассмотрен случай, когда имеется только одна линейная граница сопряжения двух методов (рис. 6, а). Размеры автоматов совпадали с размерами ячеек расчетной сетки. Совмещение движения узлов и автоматов по границе раздела осуществлялось на этапе расчета скорости с учетом сил, действующих на совмещенные с этими узлами автоматы. Были рассмотрены варианты, отличающиеся видами упаковки автоматов в дискретной области. В первом случае использовалась плотная упаковка автоматов (рис. 4), во втором — квадратная (рис. 5).

А -►

-►

-► GRID МСА

-►

-►

0

X

Рис. 6. Схема расчетной области с условиями нагружения (а), а также изменение скорости в последовательных точках образца при распространении плоской волны (б) в образце, состоящем из двух частей

Рис. 7. Картины волнового поля при симметричном расположении источника. На верхнем рисунке представлено поле скоростей в виде векторов, а на нижнем — то же поле после компьютерной обработки в виде теневой картины. Различные типы волн обозначены буквами: Р — продольная, S — поперечная, С — коническая и R — Рэлея

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

На рис. 6, б видно, что прохождение волны через границу не сопровождалось формированием отраженной волны, также не заметно искажения формы импульса после его прохождения в дискретную область. Аналогичное поведение наблюдалось и после отражения волны от жесткой тыльной поверхности дискретной (МСА) области и обратного прохождения ею границы со стороны дискретной области в континуальную.

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

На втором этапе была рассмотрена более сложная задача, а именно: задача о генерации и распространении в среде со свободной поверхностью волн всех типов. Для этого участок поверхности упругого полупространства подвергался кратковременному действию локальной вертикальной нагрузки. Были рассмотрены два случая: 1) когда источник располагался симметрично (находился точно на линии сопряжения GRID- и MCA-

Рис. 8. Картины волнового поля при несимметричном расположении источника

областей); 2) источник располагался несимметрично (был смещен относительно линии сопряжения GRID-и МСА-областей). Во всех случаях расположения источника относительно границы сопряжения анализировались детали распространения продольной, поперечной и рэлеевской волн, а также симметрия поля скоростей смещений. Тестирование проводилось как для квадратной, так и для плотной (гексагональной) упаковки автоматов в дискретной области.

В результате воздействия в среде на некотором расстоянии от источника формируются продольная (растяжения-сжатия) Р и поперечная (сдвига) S волны, распространяющиеся с различными скоростями (рис. 7, 8). Скорость распространения продольной и поперечной волн можно рассчитать как

и VS2 = М

Р

соответственно.

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

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

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

Во всех рассмотренных случаях результаты расчетов не показали существенного искажения волновых фронтов на границе сопряжения (рис. 7, 8). Имеет место лишь незначительное отличие в форме импульсов, а заметных вторичных (отраженных, преломленных и конических) волн не возникало. Следует отметить, что влияние упаковки автоматов в МСА-области проявлялось лишь в слабом изменении формы импульсов.

4. Заключение

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

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

Работа выполнена при финансовой поддержке РФФИ (гранты №№ 01-05-64482, 02-05-65346-а), Гранта Президента РФ поддержки ведущих научных школ РФ № НШ-2324.2003.1, Программы фундаментальных исследований Отделения энергетики, машиностроения, механики и процессов управления РАН №3.12 (проект № 5), Федеральной целевой научно-технической программы «Исследования и разработки по приоритетным направлениям развития науки и техники» на 2002-2006 гг. (гос. контракт №2 41.002.11.2424 от 31 января 2002 г.), проекта 13.12 Программы фундаментальных исследований Президиума РАН на 2003 год № 13, проекта 6.5.2 Программы специализированных отделений РАН.

Авторы выражают благодарность А.В. Димаки и Р.А. Бакееву за разработку графического интерфейса и средств визуализации и обработки результатов расчетов.

Литература

1. Cundall P.A., Strack O.D.L. A discrete numerical model for granular assemblies // Geotechnique. - 1979. - V. 29. - No. 1. - 1979. - P. 4765.

2. Herrmann H.J. Simulatig granular media on the computer // 3rd Granada Lectures in Computational Physics. - Ed. by P.L. Garrido, J. Marro. - Heidelberg: Springer, 1995. - P. 67-114.

3. Luding S. Granular materials under vibration: Simulation of rotating spheres // Phys. Rev. E. - 1995. - V. 52. - No. 4. - P. 4442-4457.

4. Poschel T. Granular material flowing down an inclined chute: A mole-

cular dynamics simulation // J. Phys. VII. - 1993. - No. 3. - P. 27^0.

5. Walton O.R. Numerical simulation of inclined chute flows of monodisperse, inelastic, friction spheres // Mechanics of Materials. - 1993. -No. 16.- P. 239-247.

6. Herrmann H.J., Luding S. Modeling granular media on the computer // Continuum Mech. Thermodyn. - 1998. - No. 10. - P. 189-231.

7. Псахъе С.Г., Хори Я., Коростелев С.Ю., Смолин А.Ю., Дмитриев А.И., Шилъко Е.В., Алексеев С.В. Метод подвижных клеточных автоматов как инструмент для моделирования в рамках физической мезомеханики // Изв. вузов. Физика. - 1995. - Вып. 38. -№ 11. - С. 58-69.

8. Псахъе С.Г., Коростелев С.Ю., Смолин А.Ю., Дмитриев А.И., Шилъко Е.В., Моисеенко Д.Д., Татаринцев Е.М., Алексеев С.В. Метод подвижных клеточных автоматов как инструмент физической мезомеханики материалов // Физ. мезомех. - 1998. - Т. 1. - №2 1. -С. 95-108.

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

10. Псахъе С.Г., Чертов М.А., Шилъко Е.В. Интерпретация параметров метода подвижных клеточных автоматов на основе перехода к континуальному описанию // Физ. мезомех. - 2000. - Т. 3. -№ 3. - С. 93-96.

11. Psakhie S.G., Horie Y, Ostermeyer G.P., Korostelev S. Yu., Smo-linA.Yu., Shilko E.V., Dmitriev A.I., Blatnik S., Spegel M., Zavsek S. Movable cellular automata method for simulating materials with meso-structure // Theor. and Appl. Frac. Mech. - 2001. - V. 37. - Nos. 1-

3.- P. 311-334.

12. Псахъе С.Г., Ружич В.В., Смекалин О.П., Шилъко Е.В. Режимы отклика геологических сред при динамических воздействиях // Физ. мезомех. - 2001. - Т. 4. - № 1. - С. 67-71.

13. Голъдин С.В., Псахъе С.Г, Дмитриев А.И., Юшин В.И. Переупаковка структуры и возникновение подъемной силы при динамическом нагружении сыпучих грунтов // Физ. мезомех. - 2001. -Т.4. - № 3. - С. 97-103.

14. Уилкинс М.Л. Расчет упругопластических течений / Вычислительные методы в гидродинамике / Под ред. Б. Олдера, С. Фернба-ха, М. Ротенберга. - М.: Мир, 1967. - С. 212-263.

15. Нох В. Ф. СЭЛ — совместный эйлерово-лагранжев метод для расчета нестационарных двумерных задач // Вычислительные методы в гидродинамике / Под ред. Б. Олдера, С. Фернбаха, М. Ротенберга. - М.: Мир, 1967. - С. 128-184.

16. Johnson G.R. Dynamic response of axisymmetric solids subjected to impact and spin // AIAA Journal. - 1979. - V. 17. - No. 9. - P. 975979.

17. Корнеев А.Н., Николаев А.П., Шиповский И.Е. Приложение метода конечных элементов к задачам соударения твердых деформируемых тел // Численные методы решения задач теории упругости и пластичности: Матер. VII Всес. конф. - Новосибирск: ИТПМ АН СССР, 1982. - С. 122-129.

18. Gulidov A.I., Fomin V.M., Shabalin I.I. Mathematical simulation of fracture in impact problems with formation of fragments // International Journal of Fracture. - 1999. - V. 100(2). - P. 121-131.

19. Немирович-Данченко М.М., Стефанов Ю.П. Применение конечно-разностного метода в переменных Лагранжа для расчета волновых полей в сложнопостроенных средах // Геология и геофизика. -1995.- Т. 36. - № 11. - С. 96-105.

20. Стефанов Ю.П. Численное исследование поведения упругоидеальнопластических тел, содержащих неподвижную и распространяющуюся трещины, под действием квазистатических и динамических растягивающих нагрузок // Физ. мезомех. - 1998. -Т. 1. - №2. - С. 81-93.

21. Stefanov Yu.P. Wave dynamics of cracks and multiple contact surface interaction // Theor. and Appl. Frac. Mech. - 2000. - V. 34/2. - P. 101108.

Simulation of the behavior of complex media based on the combined discrete-continuous approach

S.G. Psakhie, A.Yu. Smolin, Yu.P. Stefanov, P.V. Makarov, E.V Shilko, M.A. Chertov, and E.P. Evtushenko

Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia

The paper presents an approach to simulating the behavior of heterogeneous media based on a combination of the discrete and continuous descriptions. The possibilities of this approach are investigated on the example of a problem about the elastic wave propagation in a region with the free boundary. The considered region is assumed to consist of two parts. In one part the simulation is carried out by the finite-difference method (continuous approach) and in the other by the method of movable cellular automata (discrete description). The calculation results show that the approach proposed allows an adequate description of elastic wave propagation in complex media. Both methods give identical results and on the boundary of their conjugation there are no energy loss or generation and any non-physical elastic wave distortions. The results obtained allow us to conclude that the combined discrete-continuous approach holds much promise for solving the problems related with the numerical investigation of the behavior of complex media with marked differences in physical and mechanical properties.