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

УДК 51.72, 539.3, 539.4, 533.15

Развитие подхода к моделированию деформирования и разрушения иерархически организованных гетерогенных, в том числе контрастных, сред

С.Г. Псахье1’2’3, Е.В. Шилько1,2, А.Ю. Смолин1,2, А.В. Димаки1,

А.И. Дмитриев12, Иг.С. Коноваленко1, С.В. Астафуров1, С. Завшек4

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

2 Национальный исследовательский Томский государственный университет, Томск, 634050, Россия

3 Национальный исследовательский Томский политехнический университет, Томск, 634050, Россия

4 Веленская угольная шахта, Веленье, 3320, Словения

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

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

Approach to simulation of deformation and fracture of hierarchically organized heterogeneous media, including contrast media

S.G. Psakhie1,2,3, E.V. Shilko1,2, A.Yu. Smolin1,2, A.V. Dimaki1, A.I. Dmitriev1,2,

Ig.S. Konovalenko1, S.V. Astafurov1 and S. Zavshek4

1 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia 2 National Research Tomsk State University, Tomsk, 634050, Russia 3 National Research Tomsk Polytechnic University, Tomsk, 634050, Russia 4 Velenje Coal Mine, Velenje, 3320, Slovenia

The paper concerns the development of a formalism of the movable cellular automaton method for simulation of consolidated heterogeneous elastoplastic media at different scales. Using the developed formalism as the base, an approach was formulated for construction of structural models that describe mesoscopic response (including fracture) of heterogeneous media with regard to hierarchical organization of their internal structure. In the approach, the effect of structural scale levels higher than the level under consideration is taken into account by a technique combining the particle method and conventional methods of continuum mechanics. The effect of lower structural scale levels is taken into account by determining integral characteristics of the response of representative volumes of lower scale and by specifying appropriate values of the interaction parameters of particles. The proposed formalism was advanced for description of contrast heterogeneous media whose components can assume different aggregate states. The capabilities of the particle method for description of hierarchically organized media are illustrated by studying the mechanisms of the response and fracture of materials with a developed porous structure, including contrast media.

Keywords: heterogeneous medium, hierarchy of structural scale levels, numerical simulation, movable cellular automaton method, elasticity, plasticity, fracture

1. Введение тел [1, 2] деформируемое твердое тело является иерар-

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

концепции структурных уровней де формации твердых схематически представлена ансамблем последовательно

D Псахье С.Г., Шилько Е.В., Смолин А.Ю., Димаки А.В., Дмитриев А.И., Коноваленко Иг.С., Астафуров С.В., Завшек С., 2011

сменяющих друг друга и взаимосвязанных масштабов структуры и деформационных механизмов. На каждом масштабном уровне гетерогенная среда характеризуется наличием структурных элементов нескольких основных типов. При этом каждый из этих элементов при более детальном рассмотрении также имеет сложную внутреннюю композицию, состоящую из элементов более низкого ранга [3]. Иерархическая организация структуры твердых тел приводит к тому, что их поведение под нагрузкой представляет собой сложный нелинейный процесс, в ходе которого происходит последовательное и самосогласованное вовлечение «доступных» деформационных механизмов с различной энергией вовлечения [4, 5]. Конечным этапом этого процесса является вовлечение механизмов наибольшего ранга, приводящих к разрушению системы. Поэтому не только физическое, но и математическое описание отклика гетерогенных материалов должно осуществляться в рамках концепции их многомасштабности.

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

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

Перспективным классом численных методов механики деформируемого твердого тела, «генетически» приспособленным для моделирования разрушения, являются методы частиц. В настоящее время эти методы широко применяются для изучения механического или термомеханического отклика твердых тел и жидкостей на различных масштабных уровнях вплоть до атомарного. Следует отметить, что в современном понимании термин «методы частиц» является собирательным и включает в себя весьма разноплановые численные методы, как относящиеся к классическим представителям дискретного подхода в механике, так и беcсеточные алгоритмы численного решения уравнений механики сплошных сред (например метод частиц в ячейках [8], метод сглаженных частиц [9], SPAM [10] и т.д.). Более того, в настоящее время к методам частиц зачастую относят и некоторые современные реализации традиционных численных методов, такие как точечный метод конечных элементов (particle finite element method) [11].

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

d2R, dt2

Ji

d20i

= E (if + I),

j=i

Ni

= E Mj, j=i

(1)

&

где ^ и 0, — радиус-вектор и угол поворота частицы г; т, и — ее масса и момент инерции; Щ и Щ7 — силы центрального (т.е. в направлении, соединяющем центры масс) и тангенциального взаимодействия частиц г и j; Му — момент силы; Ni — количество взаимодействующих «соседей» частицы. Интегральные реологические свойства ансамбля элементов определяются видом и параметрами потенциала (потенциальных сил) межчастичного взаимодействия.

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

сы) и, соответственно, учитывается взаимодействие только с ближайшим окружением.

Наиболее известным представителем данной группы методов является так называемый метод дискретных элементов [13, 14]. В настоящее время он широко используется для изучения деформационных процессов в гранулированных (сыпучих) и слабосвязанных средах, в частности реологических особенностей таких систем, процессов разрушения и перемешивания [15, 16]. В то же время до последнего времени возможности применения метода дискретных элементов для изучения механических явлений в консолидированных средах ограничивались главным образом хрупкими пористыми материалами [14, 17], что связано с недостаточным развитием математических моделей взаимодействия дискретных элементов. В частности, подавляющее число моделей в рамках метода дискретных элементов основано на использовании парных (двухчастичных) потенциалов/сил взаимодействия элементов. Данное упрощение может сопровождаться возникновением ряда искусственных проявлений отклика ансамбля частиц, не присущих моделируемой среде. Среди них важно отметить затруднительность корректного моделирования накопления необратимых деформаций в материалах, пластичность которых обеспечивается механизмами масштаба кристаллической решетки.

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

В рамках формализма методов частиц авторами развит подход к построению структурных моделей для мезоскопического описания отклика (включая разруше-

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

2. Формализм метода подвижных клеточных автоматов для описания механических явлений в гетерогенных упругопластических средах

2.1. Общий подход к описанию сред с различной реологией на основе многочастичного взаимодействия

В рамках развитого подхода к построению многочастичного взаимодействия выражение для силы, действующей на клеточный автомат (дискретный элемент) i со стороны окружения, записывается в следующей форме:

т ^=р = 1 у+. (2)

& 7=1

Она представляется как суперпозиция парных составляющих зависящих от пространственного положения/перемещения автомата i относительно ближайшего соседаj, и объемнозависящей составляющей связанной с совместным влиянием окружения.

При моделировании локально изотропных сред с различной реологией объемнозависящий вклад ¥га можно выразить через давление р в объеме автомата (дискретного элемента) i в следующей форме [19, 20]:

N

Ра =-А 1 пу, (3)

7=1

где Sij — площадь поверхности взаимодействия автоматов i иу; пу — единичный вектор, ориентированный вдоль линии, соединяющей центры масс автоматов; А — материальный параметр.

В данной постановке правую часть выражения (2) можно свести к сумме сил взаимодействия пар автоматов и однозначно разделить их нормальную р и тангенциальную р7 составляющие:

И,

Р = 1 (У - АР^Пу) =

7=1

N

= Е (7 гТ7), (4)

7=1

где и -Рр7іг Т — парные силы центрального и тан-

генциального взаимодействия, зависящие соответственно от величины перекрытия автоматов ^ и относительного сдвигового смещения Та-, вычисляемого с учетом вращения обоих автоматов [17, 22]; t7• — единичный вектор, ориентированный перпендикулярно линии, соединяющей центры масс автоматов. Отметим, что

хотя правая часть (4) формально соответствует записи взаимодействия элементов в классических моделях метода дискретных элементов (1) [13-17], ее принципиальным отличием является многочастичность центрального взаимодействия автоматов.

Величина давления р (или, что то же, среднего напряжения) в объеме автомата определяется на основе вычисления компонент тензора усредненных напряжений [17]. В терминах центральной F.Щ и тангенциальной FTj сил взаимодействия выражение для компонент тензора усредненных напряжений в объеме автомата имеет вид:

°«р= V х

i

N г - - п

хЕ 9ij I Fn cos 0у a cos 0j в ± FT cos % ,а sm 0- в I, (5) j=i

где а, в — осиX, Y, Zлабораторной системы координат; Vi — текущее значение объема автомата i; qij — расстояние от центра масс элемента i до центральной точки поверхности взаимодействия с соседом j; 0j a — угол между линией, соединяющей центры масс взаимодействующих автоматов i и j, и осью а лабораторной системы координат (рис. 1).

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

P _-сг„

axx + ayy +

(6)

Вычисление компонент ст'ар позволяет рассчитывать и другие инварианты тензора напряжений в объеме автомата, в частности интенсивность напряжений:

4t ^^12{(aXx - )2 + (&yy - )2 + (aZz - aXx)2 +

+ 6[(aXy )2 + (ayz )2 + (aXz )2]}^2.

yz

(7)

Рис. 1. Пример определения угла 0уа между линией, соединяющей центры масс взаимодействующих автоматов г и у, и осью а лабораторной системы координат (а = X). Здесь у, к, I, т, п — индексы автоматов, взаимодействующих с рассматриваемым автоматом г. Центр системы координат транслирован в центр автомата г

Как следует из (1), (4), (5), конкретный вид выражений для рЩ и р/ (и в частности для рД^,„ и рД^,т) определяется реологическими характеристиками моделируемой среды.

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

2.2.1. Удельные параметры взаимодействия

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

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

Ac J _

Ar,

(d + dj )/2 dj 2 djj 2

Aq, Aq,

J + J _ Ac + Ac

+ • _ ACi(J) +ACJ(i),

Al\

Ay _ shear

ZshearAt _Av

i(

"7 г г “Г>-( 7) +АУ 7(')’

Чу ',7

где символ А здесь и далее обозначает приращение соответствующего параметра за шаг по времени А? численного интегрирования уравнений движения (1); '7 — расстояние между центрами масс автоматов i иу (рис. 2); di — размер автомата г; У,1еаг — скорость относительного сдвигового перемещения автоматов пары г— (рассчитываемая с учетом их вращения [17]). Пространственные переменные е, (7) и у, (7) будем называть нормальной и сдвиговой деформациями автомата г в паре —'.

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

Силы центрального рЩ и тангенциального Рт7 взаимодействия автоматов г и у рассматриваются в удельных единицах:

\РП = СТ;Д;

FiJ _Т S

iJSiJ.

(9)

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

Рис. 2. Схема определения параметров пространственного отношения пары подвижных клеточных автоматов і и у: расстояние между центрами масс г и расстояние между центром масс автомата и центральной точкой поверхности взаимодействия и q 7)

а7 = а7

г(е„) - АР,,

(10)

(11)

7

|т- =Трак(У..).

" 7 7 '*у’

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

2.2.2. Упругопластическое взаимодействие автоматов

Напряженно-деформированное состояние изотропной линейно-упругой среды описывается обобщенным законом Гука:

кш= 2^еаа + (1 - Ж/К К

[ТаР = °Уав’

где а, в =Х, Y, Z; ааа и еаа — диагональные компоненты тензоров напряжений и деформаций; тар и уар — недиагональные компоненты; К—модуль всестороннего сжатия; ст^еал =(ахх + ауу + а22 V3 — среднее напряжение; G — модуль сдвига.

Как видно из (10) и (11), форма и содержание выражений, описывающих центральное и тангенциальное взаимодействие подвижных клеточных автоматов, аналогичны соотношениям закона Гука для диагональных и недиагональных компонент тензора напряжений. Это приводит к очевидной идее записать конкретные выражения для силового отклика автомата г на воздействие со стороны соседа у путем простой переформулировки соотношений закона Гука:

= 2^е,( 7) + (1 - 2Gi /К, К

К- = щ У,( ]),

(12)

где О и К, — модули сдвига и всестороннего сжатия материала автомата г; е, 7 и у, 7 — нормальная и сдвиговая деформации автомата г в паре г'—'; среднее напряжение атеап вычисляется в соответствии с (6).

Предложенные соотношения (12) для силы реакции автомата г на воздействие со стороны соседау не являются произвольными. Так, подстановкой соотношений (12) в (5) легко показать, что предложенные выражения для сил межавтоматного взаимодействия автоматически обеспечивают выполнение закона Гука для средних значений компонент тензоров напряжений а‘ар и деформаций еа в в объеме элемента г'. Отметим, что еа в определяются по аналогии с (5) с использованием компонент е,7 и у7 вектора относительного перемещения элементов в парах г'—'.

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

ных перемещений в паре автоматов выражения для удельных сил взаимодействия записываются в следующем виде:

Аа7 = Аа 7, = 20, Ае,( 7) +

1 -

Щ

К,

Аа,,

= 207Ае 7(,) +

1-20

К7

Аа 7

теап >

(13)

Ае

:(7) dil2 + Ае 7 (,) <17/ 2 Ач7

\Ат7 = Ат7, = ЩАУц 7) = 207АУ j(1), 1АУ,(7) ^/2 + АУ7(,) dj|2 = 7Аt•

(14)

(15)

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

Приведенные выше соотношения являются трехмерными. В то же время в ряде задач оправданным является использование двумерной постановки. При этом обычно используются приближения плосконапряженного или плоскодеформированного состояний. В рамках рассматриваемой модели это означает, что компоненты 0х2 и 0 тензора усредненных напряжений равны нулю (предполагается, что движение объектов происходит в плоскости ХУ), а компонента аг22 определяется следующим образом:

1 -2о7к,- ч

Аа122 =-----г±т-^ (Аахх + Аа' )

22 2 + 2С,/К, уу

для плоскодеформированного состояния.

Аа,22 = 0

для плосконапряженного состояния.

Соотношения (13)—( 15) дают возможность моделировать фрагменты локально изотропной линейно-упругой среды ансамблем подвижных клеточных автоматов.

Для описания отклика упругопластических сред на основе изложенной модели взаимодействия реализована теория пластического течения с критерием Мизеса. Для этого проведена адаптация алгоритма Уилкинса [7, 23] к концепции дискретных элементов. Алгоритм Уилкинса, как правило, формулируется в терминах девиа-торов напряжений Г)а (рис. 3):

4= 4м, (16)

где М — коэффициент сброса (приведения) напряжений.

В терминах напряжений для рассматриваемых компонент тензора усредненных напряжений в объеме клеточного автомата г алгоритм Уилкинса записывается в следующей форме:

|(а'аа ) = (а'аа-атеап)М, + <

[Кв ),=аавм, >

Рис. 3. Схематическое представление приведения напряжений согласно алгоритму Уилкинса. Здесь Ое] — интенсивность напряжений, вычисляемая по окончании решения упругой задачи на текущем временном шаге

где а, в =Х, У, Zи а ф в; (а'аа) и (а'ар) — приведенные компоненты тензора усредненных напряжений;

и а,

ав

исходные значения компонент, полу-

ченные в результате решения упругой задачи (13)—( 15) на текущем временном шаге; М, =а|з]/ст|п4 — текущее значение коэффициента Мдля автомата г (рис. 3); а^ — текущий радиус предельной поверхности для автомата г; а|п( вычисляется с помощью (7) по окончании решения упругой задачи на текущем временном шаге.

По аналогии с упругой задачей для корректировки удельных нормальных и тангенциальных сил взаимодействия автоматов применяются соотношения, полученные непосредственной переформулировкой алгоритма Уилкинса для усредненных напряжений (17):

К = (а,7 -атеап)М, +а;

[Т/ = 'ТуМ,,

(18)

где а7 и Tj• — приведенные значения удельных сил. Легко показать, что при такой записи подстановка соотношений (18) в выражения (5) автоматически обеспечивает приведение компонент тензора усредненных напряжений в объеме подвижного клеточного автомата г к кругу текучести.

Необходимо отметить, что приведенные значения удельных сил а'ij и 17 для автомата г в общем случае отличаются от тех же а7 и для автомата у. Ввиду

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

[4 = (а,7 ^еапЖк +°

к

теап,

,7 =т ,7Мк,

(19)

где Мк = min{MiMj}, а индекс к в (19) принимает соответствующее значение из пары {г, у}.

При реализации алгоритма Уилкинса в двумерной постановке учтены его особенности для плоскодеформи-рованного и плосконапряженного приближений [23, 24]:

(а22 ) = (а22 -атеап)М, +атеап

для плоскодеформированного состояния,

ЬМ

Ла>22) =а2.

Ь,М« + (1 - М,) для плосконапряженного состояния, где ^’122 — Z-компонента тензора девиатора усредненных напряжений, Ь, = (1 + 4 01^1), М,11 — коэффициент приведения напряжений в плосконапряженной постановке, который вычисляется на основе итерационного метода Ньютона с М, в качестве исходного приближения [24].

Реологические свойства материала подвижного клеточного автомата г определяются заданием единой кривой упрочнения = Ф(ё^) (данная зависимость так-

же носит название функции отклика автомата).

Численное интегрирование уравнений движения (1) подвижных клеточных автоматов осуществлялось с использованием скоростной схемы Верле [12], модифицированной введением предикторного цикла для оценки величины а‘ав на текущем шаге.

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

2.2.3. Моделирование разрушения

Принципиальным преимуществом методов дискретного подхода в механике является способность непосредственного моделирования явлений разрушения материала (в том числе множественного) и связывания (сцепления) фрагментов посредством изменения состояния пары частиц («связанная» пара о «несвязанная» пара, рис. 4, а). Критерием переключения состояния пары является, как правило, предельное значение силы взаимодействия или относительного перемещения [17, 22]. Возможности развитого подхода для описания взаимодействия подвижных клеточных автоматов в многочастичном приближении позволяют применять различные многопараметрические «силовые» критерии разрушения (Губера-Мизеса, Друкера-Прагера и т.д.) в качестве критериев разрыва межэлементных связей.

В рамках развитой модели вычисление этих критериев для пар связанных клеточных автоматов г и у основано на определении локальных компонент тензора напряжений на поверхности взаимодействия пары (далее будем обозначать этот тензор как а^ф') в локальной системе координатX'У' (рис 4, б). В соответствии

аа

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

с (5) в данной системе координат значения компонент

У У

и а

X у

аУу - аУ' У У пі см у' = = У

аХу -аХу - а, па' а х = ҐІ = X

для автоматов і и у тождественно равны друг другу и численно равны удельным силам центрального и тангенциального взаимодействия автоматов а,,

у

и т,,:

(21)

где а' = х', у' па' — направляющие косинусы; /а — проекция вектора напряжений (удельных сил взаимодействия автоматов і и у) на ось а'. Такими же являются значения этих компонент и на поверхности взаимодействия автоматов (а,-,- = аУУ и т, = аХУ). В то же

У уу У ху

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

а:- =

Ял +а ххЧи

а +а г’г’Чу

(22)

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

(23)

1СТт1 >СТо

[°.5(А + 1)стШ, + 1-5(а-Г)<еап >стс, где стс — соответствующее пороговое значение для рассматриваемой пары (величина, характеризующая прочность химической связи); а — отношение прочностей материала на сжатие и растяжение; стЩ и сттеап вычисляются по аналогии с (6), (7).

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

только при Сту < 0) и наличие ограничения на величину силы тангенциального взаимодействия, определяемое применяемой моделью трения поверхностей элементов.

3. Развитие подхода к учету взаимосвязи различных масштабных уровней в иерархически организованных средах

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

1. Учет «вышележащих» структурных уровней твердого тела осуществляется через «погружение» расчетной области, моделируемой ансамблем подвижных клеточных автоматов, в объем, характерные масштабы которого могут быть на несколько порядков величины выше линейных размеров ансамбля и достигать разме-

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

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

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

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

3.1. Совмещение различных численных методов для эффективного учета влияния более высоких структурных уровней

Развитый способ учета «вышележащих» структурных уровней твердого тела основан на «погружении»

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

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

СТ/ у = РЦ, (24)

уравнение неразрывности

р + РЦ= 0 (25)

уравнение энергии

(26)

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

Е/= У 2 (ц. / + и у.), (27)

где Еу — компоненты тензора деформаций Коши; и. — компоненты вектора перемещения. Полные деформации состоят из упругих и пластических:

Еу =Е/ +Е/. (28)

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

(29)

Е = а у4 г, >

Ву

Ві

'-Яу - Я гкЮ ук - ЯукЮгк >

(30)

(31)

(32)

где а у — компоненты тензора напряжений; я, — компоненты девиатора тензора напряжений; ю, — компоненты тензора скоростей вращения:

Юу = !/2(мг,І - иу,г);

Р—среднее давление; К и ц — модули сжатия и сдвига соответственно.

f (СТ е, Х) = 0, Еу =-

(33)

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

ЭФ дсту

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

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

„ ЫП

хо = хо-т—

2Фо

+(СТх )ш ( у'П - уп)+(СТх )ПЛуП - уП) -

-(СТу )П (ХП - ХП ) - (СТху )1(£П - Х4П ) -

-(СТху )Ш(Х4 - Х1 ) - (СТ ху )1У(Х1 - Х2 ) ],

'\(СТхх)1 (у2 у3 ) + (СТхх)11(у3 у4 ) +

АГ

У0 = У +

0 о 2ф,

\(СТУУ )1 (х2 Х3) + (стуу )11( Х3 Х4) +

о

(34)

+(стуу)ш(хП -хП) + (стуу&(хп - хП) --(Стху )П( УП - Уп) - (СТху )п( У3 - Уп) --(СТху )Ш(у4 - у1 ) - (СТХУ )1У(у1 - у2) ],

х0+1 = х0 + ХоАе+У2, Уп01 = у0 + УоА1п+уп, где Х0 =

= Хп+У2

= ХП -1/2

= • п+1/2

= •п-1/2.

„ , Хо = хо~'"> Уо = Уо~'Уо = Уо"2Фо

масса, соответствующая узлу О, которая вычисляется как среднее по окружающим ячейкам: фо = -У Е (рДк-

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

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

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

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

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

(2фх)о = -(стХху24+стХхУэУ +

.-Ш III , -IV ]^\ ./-I I ,

+СТхху42 + СТхху13 ) + (СТхуХ24 +

(35)

, „.II II , „.III III , „.IV IV Ч +СТхуХ31 + СТхуХ42 + ху х13 ),

(2фУ)о = (стУух24 + СТ°Х3П1 +

, „.III III , „.IV IV ч /„.I I ,

+СТУУХ42 + СТууХ13 ) (СТхуу24 +

, _П II . _Ш III . _:у IV ч +СТхуу31 + СТхуу42 + СТху у13 ),

где стХХ, о*, стХ,, —компоненты напряжений в соот-хх УУ ХУ ^ к к к к к к ветствующих ячейках; хтп = хт - ХП , Утп = Ут - У п и т.д. Верхний индекс означает номер ячейки, а нижний индекс для координат — номер узла (рис. 5).

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

(36)

гХ = Е К ’ ту = Е К

I =Ш

I =1, N

КХ =-СтХх V + СтХуАХ.

К =&ууАХ -стХуV.

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

Ет-І Ґ-І I , _П II ,

Рх = _(а;схУ24 + аххУ31 +

і

, „.III III , „.IV , Ґ0Л I ,

+ ахху42 +а хху13 ) + (ахуХ24 +

, II II , III III , IX 13 ч

+ ахух31 + ахух42 + ахух13 ),

_Ш III

„ х4.

ху .IV IV ч

+ а уух31 +

(37)

Е ру = (а^24-

і

, „.III III , „.IV IV ч /„.I I ,

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

, _П II . _Ш III . „IV IV ч + ахуу31 + ахуу42 + ахуу13 ).

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

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

V Ы -( „III,,III , „Ш ПК ,

Е рх = ( ахху42 + ахух42 ) +

+ (^ув +а5; х-) + Р.,

Е ру = (а уу х42 -а ху у42)+

(38)

+ (а:

■IV х IV уух13

„IV р'

— а ху у13 ) + Ру,

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

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

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

1. В континуальной области решается система уравнений (24)-(33), осуществляется расчет напряженно-

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

деформированного состояния, рассчитываются скорости и координаты.

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

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

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

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

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

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

1) когда источник располагался симметрично, т.е. точно на линии сопряжения сеточной и клеточно-автоматной областей;

2) источник располагался несимметрично (был смещен относительно линии сопряжения сеточной и клеточно-автоматной областей).

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

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

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

GRID [R] I MCA

Рис. 7. Теневая картина волнового поля при симметричном расположении источника. Левая область (GRID) описывалась конечно-разностным методом, правая (MCA) — методом подвижных клеточных автоматов. Типы волн обозначены буквами: P — продольная, S — поперечная, C — коническая и R — Рэлея

Рис. 8. Теневая картина поля скоростей при расположении источника, смещенного относительно границы сопряжения вправо. Обозначения аналогичны рис. 7

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

Заметим, что в тестовых расчетах были использованы как плотная, так и квадратная упаковки автоматов в дискретной области. При этом вид упаковки автоматов практически не сказывается на параметрах импульсов, что особенно важно, так как позволяет сделать вывод об отсутствии искусственных эффектов, наведенных симметрией расчетной модели. Размер автоматов варьировался от АН до Ай/8, где АН — шаг сетки. При этом также не наблюдалось искажений волновых фронтов.

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

3.2. Формализм описания гетерогенных сред с учетом более низких структурных уровней

Предлагаемый формализм учета «нижележащих» структурных уровней иерархически организованного твердого тела основан на определении соответствующих функций отклика автоматов. Сущность данного формализма можно эффективно продемонстрировать на примере пористых керамических материалов с иерархической поровой структурой, которые являются практически важным классом гетерогенных сред [29]. В зависимости от технологии изготовления, в них могут присутствовать поры нескольких масштабных уровней в различных соотношениях и вариантах их пространственного распределения. Например, стенки макропор (ячеек) могут включать поры меньшего размера [30, 31]. Экспериментальные исследования зачастую дают информацию об испытываемом материале только на макроуровне, в то время как крайне важным является достоверное знание о свойствах материала и на иных структурных уровнях, в частности для оценки их вклада в процессы деформации и разрушения. Обычно подобные проблемы исследуются на основе численного моделирования. Однако в рамках одноуровневого подхода явный учет особенностей структуры и поведения материала на каждом из масштабных уровней представляется невозможным. Это связано как со значительными вычислительными затратами при явном задании поро-

вой структуры материала в модели, так и с иерархическим строением порового пространства [30, 31]. Последнее обстоятельство само по себе предполагает использование многоуровнего подхода при описании подобных систем. Расчеты в данной работе проводились для модельного материала со свойствами спеченной керамики ZrO2 со средним размером пор, превышающим размер зерна, и двумя максимумами на гистограмме распределения пор по размерам [30, 31].

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

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

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

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

На микромасштабном уровне предлагаемой модели представительный объем определяется на основе анализа сходимости упругих и прочностных характеристик модельных пористых образцов по мере увеличения их размеров. Для решения этой задачи в работе моделировалось механическое поведение шести групп плоских пористых керамических образцов в условиях как сдвигового нагружения, так и одноосного сжатия. Внутри каждой группы образцы характеризовались одинаковыми размерами, но различным пространственным расположением пор. Каждая группа содержала шесть образцов. Рассматривались квадратные образцы, у которых сторона Н составляла 20, 60, 100, 150, 200 и 250 мкм. Принято допущение, что все поры рассматриваемой керамики и, соответственно, модельного материала одинаковы и имеют форму сферы. Их размер, в соответствии с максимумом на гистограмме распределения реальной керамики, составлял 3 мкм [30, 31]. Размер клеточного автомата выбран в соответствии со средним размером зерна и составлял 1 мкм. Поровая структура образцов задавалась путем удаления в случайном порядке одиночных автоматов, а также шести их ближайших соседей. Величина пористости образцов составляла 7 % [30, 31]. Исходная структура одного из модельных образцов представлена на рис. 9.

В случае сдвигового нагружения механическая нагрузка прикладывалась путем задания одинаковой скорости в горизонтальном направлении верхнему слою автоматов при жестком закреплении автоматов нижнего слоя образца (рис. 9, а). На начальном этапе скорость движения автоматов верхнего слоя нарастала по синусоидальному закону от 0 до 1 м/с, а затем полагалась постоянной (рис. 9, в). Такая схема использовалась для

V* = УХ(Ц, V, = О м/с

V* = Уу(і)

V» = 0 м/с

Рис. 9. Начальная структура модельного образца со стороной h = 60 мкм и схема приложения механической нагрузки при сдвиговом нагружении (а) и одноосном сжатии (б); закон изменения скорости движения автоматов верхнего слоя образца (в)

устранения искусственных динамических эффектов и обеспечения плавного и быстрого выхода процесса деформирования образца на квазистационарный режим. По горизонтальной оси использовались периодические граничные условия. При одноосном сжатии скорость автоматов верхнего слоя в вертикальном направлении нарастала так же, как и при сдвиге, по синусоидальному закону (до 1 м/с), а нижнего — была задана равной нулю. Автоматам верхнего и нижнего слоев образца были также разрешены горизонтальные смещения, а боковые поверхности образца были свободны (рис. 9, б). Задачи решались в условиях плоской деформации. Функция отклика автоматов имела линейный вид и соответствовала диаграмме нагружения моделируемой керамики с пористостью 2 % [30, 31]. Модуль сдвига для клеточного автомата G составлял 59.2 ГПа, коэффициент Пуассона V = 0.3. В качестве критерия разрыва межавтоматных связей использовался критерий разрушения по интенсивности касательных напряжений [32].

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

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

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

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

Результаты расчетов, отражающие относительное изменение эффективных упругих свойств монолитных (сплошных) образцов в зависимости от их размера, представлены на рис. 10. Видно, что для ЕеВ существует нелинейная сходимость к Е0, а при размерах образца 60 мкм и более величина (Е0 - Ее{г)/Е0 не превосходит 2 %. Указанная точность более чем достаточна для целей данной работы, в связи с чем монолитный образец со стороной Н = 60 мкм можно считать представительным объемом для рассматриваемой сплошной модельной среды.

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

тие стс и на сдвиг тс) в зависимости от их размеров при различных видах механического нагружения, представлены на рис. 11. Для этих характеристик также получена нелинейная сходимость в рассматриваемом диапазоне изменения размеров образцов. У образцов со стороной Н = 150 мкм величины относительного отклонения ЕеВ, Geff, стс и тс от соответствующих средних по группе величин (ЕеВ), ^ей-), (стс), (тс) составляют

0.51, 0.38, 6 и 7.6 %, что не превосходит требуемого предела (3 и 15 %) и достаточно для решения поставленной задачи.

Таким образом, показано, что пористые образцы со стороной Н = 150 мкм и более могут быть рассмотрены в качестве представительных объемов рассматриваемой модельной среды на микромасштабном уровне. На основе проведенных расчетов стало возможным продолжение исследований на более высоком, мезомасштаб-

Рис. 10. Отклонение эффективного упругого модуля модельного образца Ее£у от упругого модуля клеточного автомата Е0 в зависимости от размеров образца

- (ЕеА»/<Ее«) о.аз-

0.02- о

0.01- 0.00- -0.01- о А А

8 о і А

-0.02-

-0,03-

50

100

і.......I

150

200

50

100

150

200

в, мкм

ш

И, мкм

- <тс»/{те) . 0.1 - Гв"| . і ; ; — 0.030 0.02Т * « 0.010 * ♦ т « . . 0і* т * 8 : • , . .

0.0 - ъ —; >— —’ і- °-°°- 0 т * * ^ -001- г -- 1 3 і ♦ А О

-0.1 -

о -0.02-

-0.03-

0.2 - > 1 т

50

100

150

200

Ь, мкм

Рис. 11. Отклонение упругого модуля Ее1Еу, модуля сдвига Ge^^, прочности на сжатие ас и на сдвиг тс модельных образцов со стороной к от соответствующих средних по группе величин (Ее&), {Ое((), (^с) и (Тс) в случае одноосного сжатия (а, б) и сдвигового нагружения (в, г)

ном уровне, при этом средние по группе значения (Ее1Г) и (стс) были приняты в качестве параметров функции отклика автоматов на этом масштабном уровне. Небольшие отклонения этих величин (9 и 30 % соответственно) от экспериментально определяемых значений (связанные с двумерной постановкой задачи, предположением о виде напряженного состояния и неполным соответствием морфологии пор в модели и в керамике) были откорректированы введением поправочного коэффициента. Переход от упругого модуля, определяемого по рассчитанной диаграмме нагружения в условиях плоской деформации ЕццС, к модулю Юнга осуществлялся на основе соотношения Е = ЕццС (1 -V2) [33].

3.2.2. Проведение расчетов на мезоуровне с явным учетом поровой структуры

На мезоуровне расчеты выполнены для девяти пористых плоских квадратных образцов со стороной к = 22.5 мм. Диаметр автомата, в соответствии с размером определенного представительного объема, составлял 150 мкм. Неявный перенос информации о поро-вой структуре и определяемых ею эффективных прочностных и упругих свойствах материала с микромасштабного на мезомасштабный уровень осуществлялся посредством использования для автоматов мезоуровня функций отклика, определенных на основе расчетов для представительных объемов материала на микроуровне.

Функция отклика автоматов макроуровня имела линейный вид, а ее параметры (максимальное значение удельной силы сопротивления нагружению и упругий параметр, соответствующий модулю Юнга) составляли 846 МПа и 112 ГПа. Учет поровой структуры образцов осуществлялся явным образом, так же как и на первом (микромасштабном) уровне. В соответствии с гистограммами распределения пор по размерам пористость образцов составляла 28 %, размер пор — 450 мкм [30, 31]. Схема и характер приложения механической нагрузки аналогичны используемым в задаче об опреде-

Рис. 12. Начальная структура образца со стороной h = 22.5 мм, размером пор 450 мкм и величиной пористости 28 %

лении представительного объема. Начальная структура образца представлена на рис. 12.

3.2.3. Проверка адекватности модели на мезоуровне

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

1) линейный вид диаграммы нагружения модельных образцов с наличием горизонтального плато, соответствующего квазивязкому поведению пористых хрупких материалов с пористостью выше 20 %;

2) качественное взаимное соответствие картин разрушения модельных образцов и реальной керамики;

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

Типичные для всех модельных образцов диаграммы нагружения при различных способах приложения к ним механической нагрузки представлены на рис. 13. На диаграммах можно выделить несколько характерных участков: первый, линейный участок, соответствующий упругому деформированию образца (характерен для хрупких материалов при любом значении пористости); далее незначительный по протяженности восходящий участок со срывами кривой и почти горизонтальное плато с множественными колебаниями напряжения (только при одноосном сжатии, рис. 13, а). Этим участкам соответствуют повторяющиеся процессы генерации повреждений по всему образцу, локального растрескивания и затем—упругого деформирования материала и т.д. Последний участок диаграммы — ниспадающий, соответствует развитию макротрещины или системы макротрещин, а также генерации отдельных множественных повреждений. При сдвиговом нагружении после перечисленных участков на диаграммах наблюдается восходящий участок со «срывами» кривой и затем еще один ниспадающий участок (рис. 13, б). В силу стесненных условий деформирования развитие

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

Заметим, что горизонтальное плато на диаграмме сжатия хрупких образцов (рис. 13, а) говорит об их квазивязком разрушении, которое реализуется только при величинах пористости образца более 20 % [34]. Степень проявления этих свойств пропорциональна протяженности данного участка (плато) и различна для каждого из модельных образцов. Уменьшение длины данного участка диаграммы говорит о приближении разрушения к хрупкому типу. Естественно предположить, что как возникновение квазивязкого характера разрушения образцов, так и степень его выраженности определяются некоторым критическим значением локальной пористости. Последняя связана, в частности, с величиной общей пористости образца, а также с формой и размером пор. Стоит заметить, что в данном случае квазивязкий характер разрушения полностью обусловлен структурным фактором, поскольку в модели не учитываются ни фазовые переходы, ни перестройки решетки материала.

Сравнение рассматриваемых диаграмм с соответствующими диаграммами хрупких пористых тел в условиях сдвига и одноосного сжатия [30, 34, 35] показало их хорошее качественное соответствие. Таким образом, первый критерий адекватности модели на макроскопическом уровне, касающийся общего вида диаграммы нагружения материала, выполняется.

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

Рис. 13. Диаграммы нагружения модельных образцов со стороной к = 22.5 мм при одноосном сжатии (а) и сдвиговом нагружении (б)

Рис. 14. Картины разрушения образцов с пористостью 28 % в момент сдвиговом нагружении (б)

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

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

Для проверки выполнения третьего критерия были определены средние по группе модельных образцов эффективные значения их упругого модуля (Ее1Г) и максимальной удельной силы сопротивления нагружению (стс) при одноосном сжатии и простом сдвиге. Далее они сравнивались с соответствующими величинами, найденными из натурных экспериментов. Было показано, что отклонение (стс) и (Ее1Г) модельных образцов от экспериментальных данных не превосходит 30 и 12 % соответственно, что является достаточно хорошей точностью при моделировании высокопористых сред в плоском приближении. Это говорит о хорошем количественном соответствии расчетов и эксперимента и о выполнении третьего критерия адекватности модели.

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

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

образования в них системы макротрещин при одноосном сжатии (а) и

3.2.4. Влияние перколяционного перехода в пористой структуре керамики на ее механические свойства

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

В данной работе было проведено трехмерное компьютерное моделирование механического поведения хрупких пористых образцов при одноосном сжатии с использованием метода подвижных клеточных автоматов. Функция отклика используемых автоматов соответствовала керамике ZrO2(Y2O3) со средним размером пор, соизмеримым с размером зерна [31]. Размер автоматов, в соответствии с диаграммой распределения пор по размерам, составлял 1 мкм. Рассматривались образцы в форме параллелепипеда с квадратным основанием со стороной 10, 15, 20, 30 и 40 мкм и высотой в два раза больше стороны основания. Нагрузка прикладывалась аналогично двумерным задачам. Пористость задавалась удалением одиночных автоматов из исходной структуры автоматов, расположенных в ГЦК-упаковке. Генерировались образцы со значениями пористости от 0 до 50 % с шагом 5 %. Заметим, что максимальная

Рис. 15. Периодические структуры при удалении одиночных шаров (черные кружки) из их плотной упаковки на плоскости

пористость, которую можно получить удалением одиночных автоматов из двумерной плотной упаковки при условии изолированности пор, составляет 1/3 ~ 33.3 %. При этом поры располагаются регулярно и соответствуют структуре, представленной на рис. 15, а. Соответствующая максимальная пористость для ГЦК-упа-ковки в трехмерном случае равна 1/4 = 25 % и соответствует регулярной структуре, представленной на рис. 15, б (показана плоскость 111). При случайном удалении одиночных автоматов из плотной упаковки максимальное значение пористости, которое можно получить без слияния пор в кластеры, в двумерном случае составляет приблизительно 22.4 %. В 3D-случае слияние наступает при пористости приблизительно 17.6%. Заметим, что предел перколяции в трехмерной ГЦК-структуре для задачи узлов составляет 19.8%, а в плотной упаковке на плоскости — 69.62 % [37].

Таким образом, большинство рассматриваемых образцов (с пористостью 20 % и более) всегда содержали кластеры сообщающихся пор, причем при пористости более 25 % пористость является проницаемой. Согласно [38, 39] структура пористого пространства может определять не только прочностные, но также и упругие характеристики образцов. Важно отметить, что наличие кластеров сообщающихся пор следует рассматривать как самостоятельный элемент структуры, дополнительно к одиночным порам. Следовательно, с увеличением пористости можно ожидать изменения характера зависимости упругих и прочностных свойств модельного материала от пористости при переходе предела перколяции.

На первом этапе был определен представительный объем рассматриваемого модельного материала. Для этого анализировалась сходимость упругих и прочностных характеристик сплошных образцов (нулевая пористость) с увеличением их размеров при постоянном размере автоматов. Расчеты показали, что представительными (отклонение упругого модуля не превышало 3 %) можно считать образцы с размером основания 10 мкм. Очевидно, что представительный размер пористых образцов будет различным для различных значений пористости. Так, проведенные расчеты показали, что для пористости менее 15 % представительными (разброс в эффективном модуле упругости образцов с различными

Рис. 16. Диаграммы нагружения керамических образцов с различной пористостью. Цифры п около кривых соответствуют значениям пористости в десятках процентов С = п • 10 %

вариантами размещения пор — не более 3 %) являются образцы с основанием 20 мкм, для пористости в интервале от 15 до 35 % — образцы с размером основания 30 мкм, от 35 до 50 % — образцы с размером основания 40 мкм (максимальные отклонения модуля 6 %).

Диаграммы нагружения моделируемых трехмерных образцов (рис. 16) качественно хорошо согласуются с двумерными расчетами [38]. Количественные различия объясняются разницей в напряженно-деформируемом состоянии и степени влияния пористости на эффективные характеристики пористых сред в случае двумерных и трехмерных задач. Следует подчеркнуть, что в трехмерном случае образцы со стохастически распределенной по пространству пористостью, также как и в двумерных расчетах, способны демонстрировать ква-зивязкий характер разрушения. Из рис. 16 видно, что в трехмерном случае величина пористости таких образцов должна превышать 40 %.

Рассмотрим зависимость эффективного упругого модуля исследуемого материала от его общей пористости. На рис. 17, а квадратами нанесены средние величины по пяти представительным образцам с различными вариантами размещения пор. На этой зависимости явно можно выделить два характерных участка, связанных со структурой пористого пространства: первый соответствует одиночным порам (5-20 %), второй — кластерам сообщающихся пор (20-50 %). На рис. 17, б для сравнения представлены экспериментальные данные, взятые из [30]. Видно, что рассчитанные данные хорошо соответствуют экспериментальным. Наклоны аппроксимирующих прямых, показывающие степень влияния пористости на упругие характеристики материала на характерных участках, несколько отличаются от экспериментальных. С учетом результатов [38, 39] это можно объяснить тем, что стохастический выбор элементов при генерации пор в образце, а следовательно, и ориентации кластеров сообщающихся пор нивелирует влияние пористой структуры на механические свойства среды. Кроме того, по данным [30, 31] перко-ляционные переходы в керамике ZrO2 вызывают изме-

1пЕ-

4.0-

3.5-

3.0-

| а 1пЕ •ч б N. •

\- 4.0- •4и

\ 3.5- •м

з.о- і і і і і і і і і

-3.0

-2.5

-2.0

-1.5

-1.0 ІпС

-3.0

-2.5

-2.0 -1.5

-1.0

ІпС

Рис. 17. Зависимости в логарифмических координатах упругого модуля керамики от пористости: данные расчетов (а), экспериментальные данные (б)

нение микроструктуры. В частности, при непрерывной пористой структуре напряжения, инициируемые в керамике, ограничивают рост кристаллитов.

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

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

3.3. Развитие формализма гибридных клеточных автоматов для моделирования отклика контрастных сред

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

1пэ

7.0-

6.5-

6.0

5.5-

1—\ \± ІпБ, I . \б_

N. • .

у 7.0-

IV

\ • в4'

\1

\ 6.5 -

\ •

\ 6.0-

•\«

■ \

5.5- \

1 1 1 1 1 1 1 1 1 • \ ■ 1 I 1 I 1 I 1 I

-3.0

-2.5

-2.0

-1.5

-1.0

ІпС

-3.0 -2.5 -2.0 -1.5 -1.0 ІпС

Рис. 18. Зависимости в логарифмических координатах предела прочности керамики от пористости: данные расчетов (а), экспериментальные данные (б)

различного масштаба. В настоящее время общепризнанным является представление о важной роли газового компонента как фактора изменения локального напряженно-деформированного состояния в газоносных средах [40-43]. Особенно значимым влияние газа, содержащегося в объеме материала, становится в толще горных пород, где, начиная с определенных глубин, динамика упругого восстановления деформированной среды при ослаблении степени ее стеснения определяется, при прочих равных условиях, «не только механической упругостью пород, но и в значительной мере давлением свободного газа, формирующимся в соответствии с законами фильтрации и газовой кинетики, инициируемыми, в свою очередь, снижением механических напряжений в деформирующейся среде» [41]. Необходимо также отметить, что в условиях значительных литостатического и бокового давлений газовый компонент может оказывать существенное влияние и на характер разрушения горных пород.

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

3.3.1. Общий формализм гибридного метода клеточных автоматов

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

ставлен ансамблем подвижных клеточных автоматов, а другой — сеткой классических автоматов (рис. 19, а).

В рамках данного метода расчетный шаг состоит из двух этапов. Первый из них называется «механическим» и представляет собой шаг, выполняемый на слое подвижных клеточных автоматов. На этом этапе решаются уравнения движения подвижных автоматов, моделируются процессы разрушения и массопереноса под воздействием механических нагрузок. Второй этап связан с описанием эволюции газовой/жидкой фаз в трещинно-поровом пространстве подвижных автоматов (процессов переноса газа в порах и каналах твердого каркаса) и окружающей среде. Здесь моделирование осуществляется на сетке классических клеточных автоматов. На основе рассчитанных значений концентраций газа в порах вычисляется давление, действующее на твердый каркас со стороны газовой фазы. Геометрия макроскопических пор (имеющих размер, сопоставимый с размером подвижного автомата) проецируется на сетку классических клеточных автоматов со слоя подвижных клеточных автоматов, как показано на рис. 19, б. Указанный подход позволяет объединить решения механической и газодинамической задач и описать поведение многофазной гетерогенной среды. Необходимо отметить, что точность метода зависит от соотношения между размерами классического и подвижного автоматов. Чем меньше данное соотношение, тем больше точность как описания формы и размеров макропор, так и проецирования свойств со слоя подвижных клеточных автоматов на слой сеточных автоматов и обратно.

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

Конкретные соотношения, определяющие эволюцию твердофазного каркаса и газовой/жидкой фаз, а также характер их взаимного влияния определяются

Классические клеточные автоматы

Рис. 19. Слои подвижных и классических клеточных автоматов (а); проекция пары подвижных автоматов на неподвижную сетку (б); определение сорта классического автомата (в)

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

3.3.2. Модель переноса газов в твердофазном каркасе

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

1. Открытая пористость (или фильтрационный объем), образованная сообщающейся системой пор, каналов и трещин, имеющих выход на внешнюю поверхность. Характеризуется величиной общей пористости у, средним размером пор и их распределением по размеру, просветностью фг, характерным диаметром каналов просачивания d и др. [43, 44].

2. Закрытая пористость, образованная ансамблем изолированных микропор, не имеющих связи с фильтрационным объемом материала. Характеризуется величиной общей пористости у0, средним размером микро-пор и их распределением по размеру.

3. Непосредственно твердофазный каркас, вмещающий фильтрационный объем и ансамбль закрытых пор. Важными характеристиками каркаса являются средний размер монолитного фрагмента R, а также растворимость V, отражающая его абсорбционную способность и связывающая объемную концентрацию абсорбированного газа с концентрацией газа на поверхности рассматриваемого монолитного блока [43].

Математический формализм, применяемый для описания переноса газа в такой контрастной среде, зависит от наличия открытой/закрытой пористости и состояния заключенного в объеме материала газа. В настоящей модели для каждого компонента смеси газов используются приближения идеального газа и изотермических условий (T = const). В связи со сказанным, математические соотношения, описывающие перенос компонента газовой смеси в такой контрастной среде, включают уравнение фильтрации

др

Y— = dt дх

где р — плотность рассматриваемого газа; к — коэффициент проницаемости среды; ц — динамическая вязкость газа; Р — парциальное давление газа; уравнение диффузии газа из закрытых пор (и твердофазного каркаса) в фильтрационный объем (или обратно):

У о д- = Оев (4°)

дt дх

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

Связь между концентрацией газа в твердофазном каркасе (в том числе в закрытых микропорах) и его плотностью в фильтрационном объеме устанавливается законом Генри:

у оР^Р, V << 1. (41)

В качестве уравнения идеального газа используется уравнение Менделеева-Клапейрона:

РУ = NRgT, (42)

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

Проницаемость пористой среды удовлетворительно описывается формулой [43]

k = yd2.

(43)

3.3.2.2. Численная реализация уравнений переноса. Уравнения фильтрации (39) и диффузии (4°) решаются на сетке классических клеточных автоматов с использованием неявной разностной схемы.

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

1) элемент окружающей среды (автомат не принадлежит ни одному подвижному клеточному автомату и представляет собой область свободного пространства, заполненную газом или газовой смесью);

2) автомат принадлежит одному из клеточных автоматов и представляет собой фрагмент контрастной сре-

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

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

В связи с этим на сетке классических клеточных автоматов выделяются граничные условия двух типов:

1) контрастная среда - окружающее пространство,

2) «монолитная» среда - контрастная среда.

Для каждого компонента газовой смеси (далее нумеруется индексом {), заключенного в фильтрационном объеме твердофазного каркаса, распределение его парциального давления р на сетке классических клеточных автоматов рассчитывается на основе неявной численной схемы с разделением по пространственным переменным:

> ™

Y к-

т}П+1/2 т)П

Pi,l — Pi,l

Y к-

At

к1+1 + к1_ Ц1+1 Ц1 ?

к1 -1 + h_ Ц1 -1 Ц1

Pn+1 — р

Pi, m Pi, m

4Ax

(Pn + pn )(pn+V2 — pn+V2) — (p, l+1 + Pi, l) (p, l+1 Pi, l )

(P%1 + P?M)( P^'1 — P£f)

J

n+1/2

(44)

At

km+1 +

Цт+1 Mm

—1 + km

4Ax

(pn+yi + pn+\/2

(pi, m+1 r pi,

)(pn+1 — pn+1) — )(Pi, m+1 pi, m )

'm—1

Mm

(pn+V2 + pn+V2)(pn+1 — pn+1 )

(Pi, m i, m—1 )(Pi, m pi, m—1)

(45)

Здесь индексы I и т нумеруют ячейки квадратной сетки в горизонтальном X и вертикальном У направлениях (здесь и далее рассматривается двумерная модель); п — номер шага численной схемы; Ах — размер классического клеточного автомата. Уравнения (44) и (45) решаются методом прогонки с учетом граничных условий. Результатом их решения являются новые значения давления газа i в классическом клеточном автомате (I, т). Количество газа i на шаге (п + 1) вычисляется из уравнения состояния (42):

лп+1 = P+V

1 RgT '

(46)

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

pn+1

Ptot

= Е P

n+1

(47)

В случае присутствия в моделируемой системе смеси нескольких газов уравнения (39) и (4°) решаются независимо для каждого компонента.

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

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

3.3.2.3. Учет взаимного влияния слоев классических и подвижных клеточных автоматов. Одним из аспектов влияния слоя классических клеточных автоматов на ансамбль подвижных автоматов является учет вклада давления газов в напряженное состояние твердофазного каркаса. Данный учет может осуществляться различными способами (один из них предложен в [45]). В рамках настоящей модели используется более простой способ. Влияние газов, содержащихся в фильтрационном объеме и закрытых порах подвижного клеточного автомата, на его напряженное состояние учитывается через дополнительный вклад в величину среднего напряжения в его объеме:

pg = pci Yo + Po Y, (48)

где Po и Pcl — полные давления газовых смесей в фильтрационном объеме и закрытых порах подвижного автомата. Отметим, что Po и Pcl определяются усреднением соответствующих значений по классическим клеточным автоматам, относящимся к данному подвижному автомату. Таким образом, полная величина среднего напряжения в объеме подвижного клеточного автомата j записывается в следующем виде:

rrJ — rrJ

(49)

mean, mech g ,

где S,mean mech связано с влиянием окружающих клеточных автоматов и вычисляется по формуле (6). Отметим, что выражение (48) подразумевает однородное распределение пор, трещин и фильтрационных каналов в объеме подвижного клеточного автомата. Кроме того, развитая модель не учитывает влияния молекул газа, абсорбированных в кристаллической решетке, на величину упругой энергии твердофазного каркаса. Тем не менее, как показано в [46], такое приближение является вполне допустимым.

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

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

N

= Е (Fj + Ff) + Fgg,

(50)

к=1

где — суммарная сила со стороны внешних газов на поверхностный автомат:

N Ы

р = dJ А х££Рь Пьа. (51)

а=1Ь =1

Здесь индекс а нумерует классические клеточные автоматы, принадлежащие подвижному автомату j и граничащие с классическими автоматами свободного внешнего пространства (Щ — число таких автоматов); индекс Ь нумерует классические автоматы свободного внешнего пространства, граничащие с автоматами а (— число таких автоматов, в случае квадратной упаковки на плоскости Щ < 3); Рь — суммарное давление газов в объеме классического автомата Ь; п Ьа — единичный вектор, направленный перпендикулярно площадке, разделяющей автоматы Ь и а, в направлении автомата а (т.е. внутрь подвижного клеточного автомата).

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

'jk

C0S 0jk,а C0S 0jk,р±

±c0s 0jkSm 0jk>p)-

N Ng +qja Axdj EEc0s(

a=1b =1

Ьа ,aPb C0S(

nba >P

(52)

где qia — расстояние между центром масс подвижного автомата j и центром принадлежащего ему классического автомата a; C0S 0ja a — косинус угла между линией, соединяющей центр масс подвижного автомата j с центром принадлежащего ему классического автомата a, и осью а лабораторной системы координат;

C0S 0Пь р -косинус угла между вектором nba и осью

Р лабораторной системы координат.

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

f Y = Yin (1 + 3 Emean ) = Yin (1 + 3 °mean /KX

1.Y0 = Y0 (1 + 3Emean ) = Y0n (1 + 3 °mean /KX

(53)

где у ;п и у о — значения открытой и закрытой пористостей в объеме недеформированного подвижного клеточного автомата. Отметим, что выражение (53) не учитывает возможности коллапса пор/каналов в контрастной среде, но является вполне обоснованным приближением в случае достаточно высоких значений Р0 и Рс1.

3.3.3. Оценка возможностей гибридного метода клеточных автоматов для моделирования связанных задач механики контрастных сред

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

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

3.3.3.1. Верификация модели переноса газов в пористых твердофазных средах. Тестирование развитой модели переноса газообразных компонентов контрастной среды в твердофазном каркасе показало ее корректность и применимость для моделирования отклика таких систем в сложных условиях нагружения. Это может быть продемонстрировано на примере моделирования установления стационарного режима фильтрации газа в плоскопараллельном пористом пласте (рис. 2°). В случае идеального газа установившееся распределение давления газа по длине пласта описывается следующим аналитическим выражением:

0.1 м

Рис. 20. Структура моделируемого пласта (а) и его проекция на сетку классических клеточных автоматов (б)

р(х) р2 - р1_р1 х, (54)

где р1 — постоянное давление на левой границе пласта (в толще массива); р2 — постоянное давление на правой границе (в галерее); L — длина моделируемого участка пласта; х — координата в направлении простирания пласта. Отметим, что данное решение получено в приближении фильтрации только в направлении простирания пласта (что эквивалентно непроницаемым верхней и нижней границе приведенного на рис. 20 слоя).

При численном моделировании использовались следующие параметры: протяженность L = 0.1 м, р1 = = 5.6 МПа (характерное значение давления в пласте лигнита на одном из нижних уровней в Беленском угольном бассейне), р2 = 0.1 МПа (имитация атмосферного давления на галерее). Рассматривалась фильтрация углекислого газа С02 (молярная масса — 0.044 кг/моль, коэффициент динамической вязкости |Л = 1 • 10-5 Па-с). Б настоящем расчете диффузия газа в твердом каркасе и эволюция ансамбля подвижных клеточных автоматов не учитывались, что обусловлено необходимостью сопоставления численного решения с аналитическим выражением (54). Соотношение размеров подвижного клеточного автомата и классического автомата в проводимых расчетах составляло 8 : 1.

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

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

тете (Любляна, Словения) под руководством профессора И. Пездича.

Б данных экспериментах лигнит литотипа дробного детрита, предварительно измельченный и дегазированный, помещался в сосуд обьемом 105 мл и утрамбовывался. Масса детрита составляла 60 г. После утрамбовывания лигнита сосуд закрывался и в течение 10-15 с заполнялся газом (С02, СН4, N или их смесью) до достижения рабочего давления, равного нескольким десяткам атмосфер. Примерно через 24 ч после начала эксперимента в сосуде достигаются равновесные условия (адсорбция газа лигнитом прекращается). После этого давление в сосуде сбрасывается до атмосферного, и горловина снова закрывается для изучения динамики десорбции газа из предварительно насыщенного образца. Постоянный мониторинг изменения давления в сосуде в ходе эксперимента позволяет оценить сорбционную способность угля по отношению к находящемуся в сосуде газу.

На рис. 22, а приведена структура модельного двумерного образца с линейными размерами 5.1 х 5.1 см, обьемное содержание детрита в котором соответствует экспериментам. Следует отметить, что в настоящих расчетах учитывались процессы фильтрации и диффузии С02 в обьеме лигнита. Применялось приближение равнораспределения давления газа в свободном (т.е. не занятом лигнитом) обьеме сосуда. Использовались экспериментально определенные значения открытой у = 2 % и закрытой у о = 40 % пористости дробного детрита,

0.0 0.02 0.04 0.06 0.08 0.10

Расстояние, м

Рис. 21. Распределение давления газа по простиранию пласта

коэффициента растворимости 0.05, коэффициента диффузии Deff = 2-10-7 см2/с и характерного размера монолитных блоков R ~ 10-4 м. Б то же время однозначное экспериментальное определение величины характерного диаметра фильтрационных каналов d является трудноосуществимым. Б связи с этим данный параметр являлся модельным и подбирался из условия наилучшего согласия результатов расчетов с экспериментальными данными на начальной стадии (в течение первых минут) адсорбции/десорбции.

Результаты моделирования показали хорошее качественное и количественное согласие с экспериментальными кривыми изменения давления газа в свободном обьеме сосуда. Б качестве примера на рис. 22, б приведены данные лабораторного и численного экспериментов по адсорбции С02 при начальном давлении газа в сосуде 5.63 МПа. Полученные зависимости характеризуются наличием двух участков: 1) начальный участок быстрого падения давления в сосуде, определяемого фильтрацией С02 в открытом трещинно-поровом пространстве детрита (первые 100-150 с) и 2) основной по протяженности участок, на котором скорость адсорбции С02 контролируется диффузией в закрытое поровое пространство. Такое поведение модельной системы находится в хорошем соответствии с современными теоретическими представлениями о физической природе процессов сорбции газов пористыми материалами. Полученные результаты свидетельствуют о корректности формулировки задачи переноса газового компонента в контрастных средах.

Наилучшее согласие расчетных и экспериментальных данных получено при величине характерного диаметра фильтрационных каналов d ~ 2.5 • 10-9 нм. Следует отметить, что такое характерное значение d обеспечивает хорошее согласие результатов моделирования и с данными более ранних лабораторных экспериментов, проведенных на образцах различных литотипов веленского лигнита [47]. Полученный результат демон-

стрирует, что несмотря на широкий спектр размеров фильтрационных каналов в образцах лигнита (10-910-5 м) определяющую роль в процессе фильтрации газа играют наиболее мелкие из них, имеющие характерный диаметр нанометрового диапазона. Физическая интерпретация данного результата может быть следующей. Крупные фильтрационные каналы (в том числе трещины) в лигните являются разрозненными, т.е. не образующими связанную сеть во всем обьеме (что не исключает наличия малых областей, где они обьедине-ны в единую систему). Их соединение в фильтрационную систему, пронизывающую весь обьем лигнита, обеспечивается мостиками из мелких каналов с характерным диаметром порядка нескольких нанометров. Б этом случае эффективная скорость фильтрации в лигните контролируется пропускной способностью таких каналов, даже если обьемная доля более крупных каналов является преобладающей (поскольку только пренебрежимо малая часть фильтрационных каналов имеет выход на свободную поверхность).

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

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

3.3.3.2. Берификация гибридного метода клеточных автоматов для решения связанных задач деформирования и разрушения контрастных сред. Быше отмечалось, что газовая/жидкая фаза, содержащаяся в трещин-

Рис. 22. Структура моделируемой системы (а) и зависимость давления газа в сосуде от времени (б)

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

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

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

Б настоящей работе проведено изучение влияния скорости изменения напряженного состояния образцов лигнита на их способность сохранять целостность и несущую способность. Исследовались образцы дробного детрита, наиболее совершенного и высокопрочного из литотипов лигнита, обладающего выраженными упругохрупкими свойствами [48]. Структура моделируемой системы приведена на рис. 23, а. Б качестве газа использовался С02.

На первом этапе нагружения образца (при фиксированном положении поршня и опорной плиты, соответствующем ненагруженному образцу), давление газа СО2 медленно повышалось до Р^ = 4.2 МПа. При этом после окончания релаксационных процессов давление образца на поршень составляло 3.13 МПа (рис. 23, б). Отметим, что на этой стадии материал оставался практически неповрежденным. На втором этапе образец претерпевал одноосное сжатие по стандартной процедуре, описанной выше, до достижения деформации £ тх = 1.15 %, после чего положение поршня фиксировалось (при этом максимальное значение осевой нагрузки сттх = 10.2 МПа соответствует повышенным, т.е. так называемым динамическим, вертикальным напряжениям в пластах веленского лигнита). Отметим, что на второй стадии в образце возникло значительное количество мезоповреждений, моделируемых разорванными межавтоматными связями (около 15 % от их общего

а, МПа 8 - Ч \б_

6 -

4 -

2 - 1 11 III

Время

Рис. 23. Схематическое изображение структуры моделируемой системы (а) и изменения удельного значения силы реакции образца дробного детрита на поршень на различных стадиях численного эксперимента (б): I — стадия линейного возрастания ру до 4.2 МПа; II — стадия сжатия образца до £т^ = 1.15 %; III — стадия бокового разгружения

Рис. 24. Структуры межавтоматных связей образца дробного детрита по окончании заключительной стадии бокового разгружения с различными скоростями: <т1 = 0.25 (а), 1.25 (б), 10 мс (в). Точками показаны центры масс подвижных клеточных автоматов. Линиями соединены центры масс связанных автоматов

числа). На третьей стадии осуществлялось боковое раз-гружение образца путем снижения давления С02 по линейному закону в течение определенного интервала времени £ш1.

Результаты расчетов показали, что при ^ип1 > 10-3 с (1 мс) моделируемый образец дробного детрита сохраняет целостность (рис. 24) и несущую способность (рис. 25). При этом снижение прочности образцов вблизи порогового значения <ип1 —1 -10-3 с не превышало 35 % (рис. 25, кривая 2). Дальнейшее увеличение £ш1 сопровождается постепенным увеличением остаточной прочности до предельного значения, составляющего порядка 0.75сттх. Отметим, что основные процессы формирования мезоповреждений и трещин происходят динамически в течение интервала разгрузки ^ип1. Позднее имеют место лишь остаточные релаксационные процессы, не сопровождающиеся заметным увеличением поврежденности образца, что находит отражение на диаграммах нагружения, стабилизирующихся вскоре после достижения ^ип1 (рис. 25). Схожие закономерности были получены и в других сериях расчетов при различных р и етх.

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

ствием этого является более интенсивное развитие системы повреждений и трещин в материале.

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

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

Рис. 25. Зависимость удельной силы реакции образца дробного детрита на поршень от времени на стадии бокового разгружения: гш1 = = 0.25(1), 1.25(2), 10 мс (3)

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

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

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

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

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

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

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

Работа выполнена в рамках проектов Ш.20.2.2, Ш.20.2.3 и Vn.64.L8 Программы фундаментальных исследований СО РАН, а также при поддержке проекта программы Президиума РАН 23.2, гранта РФФИ № 09-05-00968-а, грантов Президента РФ №°№ МК-130.2010.5 и МК-5260.2010.8 и государственного контракта № 16.523.12.3002.

Литература

1. Структурные уровни пластической деформации и разрушения / Под ред. В.Е. Панина. - Новосибирск: Наука, 1990. - 255 с.

2. Панин В.Е., Лихачев В.А., Гриняев Ю.В. Структурные уровни деформации твердых тел. - Новосибирск: Наука, 1985. - 229 с.

3. Садовский М.А. Естественная кусковатость горной породы // ДАН СССР. - 1979. - Т. 247. - № 4. - С. 829-831.

4. Панин В.Е. Методология физической мезомеханики как основа построения моделей в компьютерном конструировании материалов // Изв. вузов. Физика. - 1995. - Т. 38. - № 11.- С. 6-25.

5. Панин В.Е., Гриняев Ю.В., Псахье С.Г. Физическая мезомеханика: достижения за два десятилетия развития, проблемы и перспективы // Физ. мезомех. - 2004. - Т. 7. - Спец. вып. - Ч. 1. - С. I-25-I-40.

6. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975.- 271 с.

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

8. E-vans M. W., Harlow F.H. The Partide-in-СеП Method for Hydrodynamic Calculations / Los Alamos Scientific Laboratory Report No. LA-2139. - Los Alamos, 1957.

9. Monaghan JJ. Smoothed particle hydrodynamics // Rep. Progr. Phys. -

2005. - V. 68. - No. 8. - P. 1703-1759.

10. Hoover W.G., Hoover C.G. Computational physics with particles // Amer. J. Phys. - 2008. - V. 76. - No. 4-5. - P. 481-492.

11. Onate E., Idelson S.R., Del Pin F., Aubry R. The particle finite element method. An overview // Int. J. Comput. Meth. - 2004. - V. 1. -No. 2. - P. 267-307.

12. Rieth M. Nano-engineering in Science and Technology. - Singapore: World Scientific, 2003. - 151 p.

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

14. Jing L., Stephansson O. Fundamentals of Discrete Element Method for Rock Engineering: Theory and Applications. - Oxford: Elsevier,

2007. - 450 p.

15. Sibille L., Nicot F., Donze F.V., Darve F. Material instability in granular assemblies from fundamentally different models // Int. J. Numer. Anal. Meth. Geomech. - 2007. - V. 31. - No. 3. - P. 457-481.

16. Martin C.L., BouvardD. Study of the cold compaction of composite powders by the discrete element method // Acta Mater. - 2003. -V. 51.- No. 2. - P. 373-386.

17. Potyondy D.O., Cundall P.A. A bonded-particle model for rock // Int. J. Rock Mech. Min. Sci. - 2004. - V. 41. - No. 8. - P. 1329-1364.

18. Daw M.S., Foiles S.M., Baskes M.I. The embedded-atom method: A review of theory and applications // Mat. Sci. Rep. - 1993. - V. 9. -No. 7-8. - P. 251-310.

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

20. Psakhie S.G., Horie Y., Ostermeyer G.-P. et al. Movable cellular automata method for simulating materials with mesostructure // Theor. Appl. Fract. Mech. - 2001. - V. 37. - No. 1-3. - P. 311-334.

21. Psakhie S.G., Shilko E.V., Smolin A. Yu., Dmitriev A.I., Korostelev S.Yu. Computer Aided Study of Reaction-Assisted Powder Mixture Shock Compaction at Meso-scale. New Computational Technique // Proc. US-Russian Workshop «Shock Induced Chemical Processing». - St.-Petersburg, 1996. - P. 21-34.

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

23. Wilkins M.L. Computer Simulation of Dynamic Phenomena. - Berlin: Springer-Verlag, 1999. - 246 p.

24. Wilkins M.L., Gutman M.W. Plane Stress Calculations with a Two Dimensional Elastic-Plastic Computer Program. - University of California, Lawrence Livermore Laboratory, 1976. - 15 p. / Preprint UCRL-77251.

25. Псахье С.Г., Смолин А.Ю., Стефанов Ю.П., Макаров П.В., Шилько Е.В., Чертов М.А., Евтушенко Е.П. Моделирование поведения сложных сред на основе комбинированного дискретно-континуального подхода // Физ. мезомех. - 2003. - Т. 6. - № 6. - С. 1121.

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

27. Дмитриев А.И., Смолин А.Ю., Попов В.Л., Псахье С.Г. Многоуровневое моделирование процессов трения и износа на основе численных методов дискретной механики и феноменологической теории // Физ. мезомех. - 2008. - Т. 11. - № 4. - С. 15-24.

28. Смолин А.Ю., Добрынин С.А., Псахье С.Г. Факторы, определяющие генерацию упругих волн при трении. Моделирование на основе дискретно-континуального подхода // Изв. ТПУ. - 2009. -Т. 314. - № 2. - С. 76-82.

29. Global Roadmap for Ceramics: Proc. 2nd Int. Congress on Ceramics (ICC2) / Ed. by A. Belosi, G.N. Babini. - Verona (Italy): Institute of Science and Technology for Ceramics, National Research Council,

2008. - 833 p.

30. Буякова С.П. Свойства, структура, фазовый состав и закономерности формирования пористых наносистем на основе ZrO2 / Дис. ... докт. техн. наук. - Томск: ИФПМ СО РАН, 2008. - 309 с.

31. Кульков С.Н., Буякова С.П., Масловский В.И. Структура, фазовый состав и механические свойства керамик на основе диоксида циркония // Вестник ТГУ. - 2003. - № 13. - С. 34-57.

32. Качанов Л.М. Основы механики разрушения. - М.: Наука, 1974. -312 с.

33. Седов Л.И. Механика сплошной среды. Т. 2. - М.: Наука, 1970. -568 с.

34. Коноваленко И.С. Теоретическое исследование деформации и разрушения пористых материалов медицинского назначения и биомеханических конструкций / Дис. ... канд. физ.-мат. наук. -Томск: ИФПМ СО РАН, 2007. - 174 с.

35. Кульков С.Н., Масловский В.И., Буякова С.П., Никитин Д.С. Негуковское поведение пористого диоксида циркония при активной деформации сжатием // ЖТФ. - 2002. - Т. 72. - № 3. - С. 3842.

36. Kovacik J. Correlation between shear modulus and porosity in porous materials // J. Mater. Sci. Lett. - 2001. - V. 20. - P. 1953-1955.

37. Клеман М., Лаврентович О.Д. Основы физики частично упорядоченных сред. - М.: Физматлит, 2007. - 680 c.

38. Смолин А.Ю., Коноваленко Иг.С., Кульков С.Н., Псахье С.Г. О возможности квазивязкого разрушения хрупких сред со стохастическим распределением пор // Письма в ЖТФ. - 2006. - Т. 32. -№ 17. - С. 7-14.

39. Коноваленко Иг.С., Смолин А.Ю., Коростелев С.Ю., Псахье С.Г. О зависимости макроскопических упругих свойств пористых сред от параметров стохастического пространственного распределения пор // ЖТФ. - 2009. - Т. 79. - № 5. - С. 155-158.

40. Христианович С.А., Коваленко Ю.Ф. Об измерении давления газа в угольных пластах // ФТПРПИ. - 1988. - № 3. - С. 3-24.

41. Полевщиков Г.Я. Динамические газопроявления при проведении подготовительных и вскрывающих выработок в угольных шахтах. - Кемерово: ИУУ СО РАН, 2003. - 326 c.

42. Полевщиков Г.Я., Козырева Е.Н. Газокинетический паттерн разрабатываемого массива горных пород // Горный информационноаналитический бюллетень. - 2002. - № 11. - С. 117-120.

43. Алексеев А.Д., Василенко Т.А., Гуменник К.В. и др. Диффузионнофильтрационная модель выхода метана из угольного пласта // ЖТФ. - 2007. - Т. 77. - № 4. - С. 65-74.

44. Христианович С.А. Об основах теории фильтрации // ФТПРПИ. -1989.- № 5. - С. 3-18.

45. Borisenko A.A. Effect of gas pressure on stresses in coal strata // J. Min. Sci. - 1985. - V. 21. - No. 1. - P. 88-92.

46. Алексеев А.Д., Синолицкий В.В. Теория абсорбции газа пористыми веществами // Физика и техника высоких давлений. - 1983. -№12.- С. 103-106.

47. PezdicJ., MarkicM., LeticM. et al. Laboratory simulation of ad-sorption-desorption processes on different lignite lithotypes from the Velenje lignite mine // RMZ - Materials and Geoenvironment. -1999. - V. 46/3. - P. 555-568.

48. Psakhie S.G., Zavshek S., Jezershek J. et al. Computer-aided examination and forecast of strength properties of heterogeneous coal-beds // Comput. Mater. Sci. - 2000. - V. 19. - No. 1-4. - P. 69-76.

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

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

Псахье Сергей Григорьевич, д.ф.-м.н., проф., дир. ИФПМ СО РАН, проф. ТГУ, зав. каф. ТПУ, sp@ispms.tsc.ru Шилько Евгений Викторович, д.ф.-м.н., внс ИФПМ СО РАН, проф. ТГУ, shilko@ispms.tsc.ru Смолин Алексей Юрьевич, д.ф.-м.н., доц., снс ИФПМ СО РАН, доц. ТГУ, asmolin@ispms.tsc.ru Димаки Андрей Викторович, к.т.н., нс ИФПМ СО РАН, dav@ispms.tsc.ru

Дмитриев Андрей Иванович, д.ф.-м.н., доц., внс ИФПМ СО РАН, проф. ТГУ, dmitr@ispms.tsc.ru Коноваленко Игорь Сергеевич, к.ф.-м.н., мнс ИФПМ СО РАН, igkon@ispms.tsc.ru Астафуров Сергей Владимирович, к.ф.-м.н., нс ИФПМ СО РАН, astaf@ispms.tsc.ru Завшек Симон, Ph.D., нач. отд. Веленской угольной шахты, Simon.zavsek@rlv.si