CC BY
69
ЧЕБЫШЕВСКИЙ СБОРНИК Том 18 Выпуск 3
УДК 539.3:539.62:004.94 Б01 10.22405/2226-8383-2017-18-3-439-460
МЕТОД ПОДВИЖНЫХ КЛЕТОЧНЫХ АВТОМАТОВ КАК НАПРАВЛЕНИЕ ДИСКРЕТНОЙ ВЫЧИСЛИТЕЛЬНОЙ МЕХАНИКИ1
С. Г. Псахье2, А. Ю. Смолин3, А. И. Дмитриев4, Е. В. Шилько5, С. Ю. Коростелев6 (г. Томск)
Аннотация
В статье рассмотрены основные положения метода подвижных клеточных автоматов, который предназначен для моделирования процессов деформирования и разрушения материалов и сред в рамках метода частиц на различных масштабах. Изначально метод частиц в механике материалов применялся только для моделирования поведения материалов на микроуровне в виде метода молекулярной динамики. Дальнейшее его развитие привело к целой группе методов под общим названием метода дискретных элементов, которые в основном применяются для моделирования сыпучих и гранулированных материалов на макроуровне. Рассматриваемый в работе метод подвижных клеточных автоматов разработан для моделирования процессов деформирования и разрушения материалов на различных масштабах: на мезоскопиче-ском масштабе с явным учётом структуры материала и на макроскопическом масштабе в рамках среды с эффективными свойствами. В работе изложены важные отличия и преимущества данного подхода по сравнению с другими методами современной дискретной вычислительной механики. Эти преимущества обусловлены прежде всего тем,
1 Исследование выполнено за счёт проекта 111.23.2.4 Программы фундаментальных исследований СО РАН на 2017-2020 годы.
2 Псахье Сергей Григорьевич, член-корреспондент Российской академии наук, директор Института физики прочности и материаловедения Сибирского отделения РАН, [email protected]
3 Смолин Алексей Юрьевич, ведущий научный сотрудник Института физики прочности и материаловедения Сибирского отделения РАН, [email protected]
4Дмитриев Андрей Иванович, профессор кафедры физики металлов Томского государственного университета, ведущий научный сотрудник Института физики прочности и материаловедения Сибирского отделения РАН
5 Шилько Евгений Викторович, профессор кафедры физики металлов, ведущий научный сотрудник Института физики прочности и материаловедения Сибирского отделения РАН
6 Коростелев Сергей Юрьевич, старший научный сотрудник Института физики прочности и материаловедения Сибирского отделения РАН
что представленный здесь подход основывается на двух базовых методах дискретного моделирования: методе частиц и методе клеточных автоматов. Использование формализма клеточных автоматов иозволя-ет явно описывать как процессы зарождения и развития повреждений (разрушения), так и залечивания трещин и микросварки. Кроме того, в рамках этого же формализма возможно описание процессов теплопроводности, химических и фазовых превращений. Вторым важным преимуществом метода подвижных клеточных автоматов является многочастичный характер взаимодействия его элементов. В результате использования многочастичного взаимодействия удаётся избавиться от искусственного влияния упаковки частиц и локальности их взаимодействия в точках контакта на поведение моделируемого материала, что наиболее важно для моделирования его упруго-пластического течения. В плане дальнейшего развития рассмотренного подхода в работе предложены способы описания в рамках метода частиц контактного взаимодействия поверхностей различных тел на микро- и мезоскопическом масштабах.
Ключевые слова: вычислительная механика, метод частиц, клеточные автоматы, деформация и разрушение.
Библиография: 31 название.
MOVABLE CELLULAR AUTOMATON METHOD AS A TREND IN DISCRETE COMPUTATIONAL MECHANICS
S. G. Psakhie, A. Yu. Smolin, A. I. Dmitriev, E. V. Shilko, S. Yu. Korostelev (Tomsk)
Abstract
The paper presents the basics of movable cellular automaton method aimed for simulating deformation and fracture of materials and media at different scales. Initially, the particle method has been employed in mechanics of materials only at microscale as molecular dynamics. Its further development has been led to a group of methods which are usually called as discrete element method and used for simulation of loose and granular materials at the macroscale. The presented method of movable cellular automata was developed for simulating deformation and fracture of materials at different scales: at mesoscale with an explicit account for material structure, and at macroscale within the framework of a media with effective properties. The main advantages and differences of the approach compared with the other methods of discrete computational mechanics are considered. These advantages, first of all, are determined by the fact that the considered approach is based on two basic methods of discrete simulation: particle method and cellular automaton method. Employing
the formalism of cellular automata allows explicit description of both processes of damage generation and evolution as well as of crack healing and microwelding. More of that, it is possible to describe heat transfer, chemical reactions and phase transitions as well. The second important advantage of the movable cellular automaton method is the many-body type of interaction among its elements. The use of many-body interaction allows us to avoid artificial effect of the particle packing and locality of their interaction on the resulting behavior of the modeled material that is extremely important for modeling elastic-plastic matereials. As a further development of the considered approach, two techniques are discussed which enable to describe contact interaction of solid bodies surfaces at the micro-and mesoscopic scales within the framework of the particle method.
Keywords: computational mechanics, particle method, cellular automata, deformation and fracture.
Bibliography: 31 titles.
1. Введение
Механика — это наука о движении в самом общем смысле этого понятия. Если говорить о механике материалов и сред, то объектом её описания может быть как некий континуум, так и множество взаимодействующих частиц (выделенных элементов среды). Эти два различных подхода были предложены соответственно Огюетеном Луи Коши и Анри Навье ещё в начале 1820-х годов [1]. Преимущественное развитие континуального подхода на протяжении последующих более чем ста пятидесяти лет было обусловлено более широкими возможностями аналитического описания.
Важный этап развития механики, — разработка и создание численных методов, — был связан с появлением вычислительной техники. Вычислительную механику можно смело рассматривать как самостоятельный раздел механики, дополняющий теоретические и экспериментальные исследования. Для неё характерны как теоретический, так и экспериментальный аспекты. По этой причине компьютерное моделирование часто называют компьютерным экспериментом.
Специфика численного решения уравнений механики сплошных сред требует дискретизации изучаемого объекта, что приводит к ряду дополнительных проблем, связанных в первую очередь со сходимостью и устойчивостью получаемых решений. Большинство этих проблем, а также проблема описания процессов возникновения и развития нееплошноетей в рамках континуума были успешно решены [2, 3]. Успехи современной континуальной вычислительной механики трудно переоценить. Её инженерные приложения совершили своеобразную революцию в технических науках и проектировании [4].
Дискретная вычислительная механика вначале развивалась главным образом на атомном уровне как молекулярная динамика. Необходимо отме-
тить, что хотя молекулярная динамика и получила исключительно широкое развитие во многих областях современной науки начиная от физики и кончая биологией, она, безусловно, не исчерпывает всех возможностей дискретного подхода, 40-45 лет назад появились работы, в которых метод частиц использовался для описания различных сред, в том числе гранулированных и сыпучих. Это в первую очередь работы таких зарубежных учёных, как P. A. Cundall, Н, J, Herrmann, О, Е, Walton, S, Luding, Т. Pöschel, D, Greenspan, G, P. Ostermeyer и др. [5, 6, 7, 8, 9, 10,11, 12, 13]. Фактически все эти работы были основаны на технике и уравнениях движения молекулярной динамики, с той лишь разницей, что добавлялись уравнения для вращения частиц. Т. е. использовались классические уравнения движения Ньютона-Эйлера. В рамках этого подхода информация о взаимодействии частиц (в том числе атомов) передается в радиусе ближайших соседей за один шаг интегрирования уравнений движения (т. е. практически мгновенно) [16, 17]. Если для описания поведения систем на атомном уровне это обосновано в рамках адиабатического приближения, то для частиц на мезо- и макроуровне это является более грубым приближением и требует использования искусственных приёмов при интегрировании уравнений движения.
В настоящее время наиболее известным представителем дискретной вычислительной механики является так называемый метод дискретных элементов (discrete element method — DEM). Точнее говоря, под этим термином понимается целая группа методов, в которых моделируемый объект представляется набором взаимодействующих частиц различной формы (элементов). Главное отличие различных представителей этих методов заключается в аппроксимации деформации этих элементов. И здесь можно выделить два подхода: отдельных элементов (distinct elements) и однородно деформируемых элементов (simply deformable elements) [14, 15]. В первом случае дискретные элементы рассматриваются как абсолютно жёсткие тела, связанные между собой пружинами с заданными жёсткостями и вяз костями. Второй метод является более строгим и учитывает возможность деформирования образца в линейном приближении распределения смещений в объёме элемента. Фактически данное приближение подразумевает, что напряжения и деформации однородно распределены в объёме элемента.
Таким образом, подходы, основанные на классических уравнениях движения Ньютона-Эйлера, хотя и включают в себя уравнения для вращения, фактически описывают поведение точечных масс в парном приближении. Данный подход хорошо подходит для описания упруго-хрупких материалов и сред, однако при описании больших упругих, а также пластических деформаций сказывается искусственное влияние упаковки частиц и локальность их взаимодействия в точках контакта.
Новый вид уравнений движения, позволяющий преодолеть данное противоречие был получен в рамках метода подвижных клеточных автоматов [18, 19, 20, 21, 22].
2. Описание моделируемой системы
В рамках метода подвижных клеточных автоматов (movable cellular automata — MCA) моделируемая система представляет собой ансамбль взаимодействующих автоматов (элементов), имеющих конечный размер. Концепция метода МСА основана на введении нового типа состояний в подходе классических клеточных автоматов — состояния пары автоматов. Это позволило сделать принципиально важный шаг — перейти к использованию пространственной переменной как параметра переключения, В качестве такого параметра было выбрано перекрытие пары автоматов (рис, 1):
где г4 — расстояние между автоматами г и ] (верхние индексы в данной работе будут обозначать номера автоматов); нижний индекс 0 обозначает начальный момент времени (недеформируемая конфигурация, расстояние в которой вычисляется как г0 = (¿г + в?)/2, где параметр с1г характеризует размер автомата г).
Рис, 1: Пространственные параметры пары подвижных клеточных автоматов г и j.
В простейшем случае существует два состояния пар автоматов: связанное (linked) и несвязанное (unlinked). Связанное состояние означает наличие химических связей между элементами (они образуют единый фрагмент консолидированного тела), а несвязанное — отсутствие каких-либо химических связей (элементы принадлежат раздельным фрагментам разрушенного тела или различным телам). Следует подчеркнуть, что данный параметр имеет пространственную размерность. Это приводит к качественно новому свойству клеточных автоматов — способности изменять своё пространственное
(1)
hIJ
r 'J
положение и как результат, — изменять пространственную конфигурацию всей системы. Таким образом, изменение состояния пар автоматов определяется (управляется) относительными смещениями (в общем случае также поворотами) автоматов, образующих пару. Поскольку названные состояния являются устойчивыми, то среда, образуемая парами таких автоматов, может рассматриваться как бистабильная.
Рассмотрим пару автоматов г]. Состояние данной пары однозначно определяется двумя типами полей (¿). Первый тип связан с состоянием самой пары ii — № (£), второй — с состояниями соседних пар, а именно: пар г к — Нгк (¿) и jl — Ь?1 (¿) Следует отметить, что эти поля по существу определяют локальные деформации (деформации элементов среды) и, следовательно, определяют распределение и потоки упругой энергии в моделируемой среде. Производная по времени от параметра перекрытия может быть определена подобно тому как это делается в модели Винера,-Розенблюта [23]. В этом случае рассматриваемая система может быть описана на основе соотношения:
= №) + шты11) + ттп1), (2)
к=з 1=г
где функция ) имеет смысл нормальной составляющей относительной скорости автоматов г и ] (У^); С (г], гк) — коэффициент, связанный с переносом параметра Ы от пары г к к паре г]] 1(Нгк) — явная функция от Ык, которая определяет перераспределение параметра Нгк между п арами: ц и г к (рис. 2).
Рис. 2: Параметры перекрытия соседних пар элементов 13, г к и jl. В рамках линейного приближения функция 1(Нгк) может быть представ-
лена как:
1(кк) = ф(агз'гк )У1к, (3)
В общем случае множитель ф(аг^'гк) определяется взаимным расположением автоматов иг к и является комбинацией тригонометрических функций, а аг^'гк — набор значений соответствующих углов, определяемых относительной упаковкой.
Таким образом, изменение параметра перекрытия автоматов к^ определяется следующими факторами:
1, относительной скоростью между автоматами г и
2, изменением перекрытий Дк4' г-го автомата с его соседями;
3, изменением перекрытий Дк1 ]-го автомата с его соседями.
3. Уравнения движения подвижных клеточных автоматов
Как следует из соотношения (2) система, состоящая из подвижных клеточных автоматов, будет описываться следующими уравнениями движения для трансляционного поведения:
А2у.г] , 11 ч 1 1
^т = (—++Е С * )—р* +Е С шшс?»1 )—у.
at2 \тг т? / ^ тг т?
k=j 1=г
(4)
ayv 1
В уравнении (4) вклад —был представлен как —:ргк, где ргк — сила
at т\
нормального взаимодействия между автоматами i и к, тг — масса автомата
г. Коэффициент С(ij,ik) как и в соотношении (2), связан с переносом параметра h от пары г к к паре ij. Член ф(аг:>'гк) определяется коэффициентом Пуассона и соответственно связан со взаимным расположением пар элементов ij и г к.
Уравнения движения для поворотов могут быть записаны в виде:
dw = (т + Т У' + Е5 (»■*) рк + Е5 Ш о Р1- и
dt jjJ кфз Jг i= J
Здесь скорость относительного поворота элементов (сам относительный
h
j
Коэффициент S (ij, г к) связан с переносом относительного вращения от па-отражена в работе [24],
Хорошо видно, что если положить С(1],1к(]1)) = 1 и Б(1],1к(]1)) = 1, то уравнения (4) и (5) полностью эквивалентны уравнениям Ньютона-Эйлера для многочастичного взаимодействия.
Как отмечалось выше, члены с множителями С(%], гк(]1)) определяют перенос параметра к от пар г к ми к паре Таким образом, уравнения движения (4) и (5) явно учитывают некоторое «запаздывание» 5Ьг^'гк и б^'Э1
размерами автоматов, продольной и поперечной скоростями звука, а также взаимным расположением пары и пар г к ]1. Очевидно, что эффекты, связанные с генерацией и распространением ударных волн, требуют специального рассмотрения.
4. Явный учёт отклика, обусловленного изменением объёма
Уравнения движения (4) и (5) могут быть записаны в следующем виде
[22]:
m
J1
:d2Rг
dûг ~ddt =
/ , FpLir + Fh
^Mij
(6)
где Кг — вектор положения, шг — скорость вращения, тг мент инерции автомата г, соответственно, ^
pair
pij + Tij
масса, J% — мо— общая парная
сила взаимодействия автоматов г и j, Р^ — объёмно-зависящая сила, дей-
автоматов, M%j
дгз (пгз х + Хсл _ вращательный момент, действующий - расстояние от цнентра автомата г до точки его взаимодействия с автоматом пг^ = ( К — Кг) /ггг — единичный вектор, направленный
ij
ij
расстояние между центрами
автоматов, К^ — вращательный момент, вызванный относительным вращением автоматов пары [21, 22],
В рамках подхода модели «погруженной частицы», аналогичной модели погруженного атома в молекулярной динамике [25], объёмно-зависящая сила может быть представлена как [21, 22]
FX
2 Gj
) Р jAlj п
(7)
где индекс ](] = 1... Nг) нумерует автоматы, взаимодействующие с автоматом ц Рз — давление j-тo соседнего автомата; А4 — площадь контакта г-го автомата с С^ — модуль сдвига, К — модуль сжатия материала автома-
(линейном) случае эта зависимость имеет вид:
О-7' — О-7
р1 = к 0 . °о /8)
где 030 — начальный (равновесный) объём автомата; О-7 — его текущий объём,
Д
но определить исходя из соответствующих изменений расстояний от центра
Д к
ДО3 = ^ ДдзкАзк. (9)
5. Функция отклика подвижного клеточного автомата
Как известно из механики сплошной среды формоизменение при деформировании описывается девиатором тензора деформаций или вторым ин-
напряжения сдвига характеризуются вторым инвариантом тензора напряжений (интенсивностью напряжений а). Как показано в работах [21, 22], в приближении однородно деформируемого элемента можно вычислить средние значения этих величин в подвижном клеточном автомате через силы в точках его взаимодействия с соседними автоматами, В случае линейно-упругого тела значения интенсивности напряжений пропорциональны значениям интенсивности деформаций с коэффициентом пропорциональности, равном утроенному модулю сдвига, В случае пластического течения часто постулируется наличие единой кривой упрочнения, как а = f(e) для произвольного типа нагружения. Экспериментально это подтверждено для простого нагружения,
а = ( )
то называть функцией отклика. Она является аналогом определяющих соотношений механики деформируемого твёрдого тела и характеризует силовой отклик материала автомата при деформировании.
Можно выделить четыре основных типа функции отклика автоматов, используемых при описании различных материалов и сред (рис, 3), В простейшем случае межавтоматное взаимодействие полагается упругим и линейным. Вид такой функции показан на рис, 3, а. Следует отметить, что в этом случае нагрузка и разгрузка идут по одной и той же линии. Переключение пары автоматов из состояния linked в состояние unlinked происходит когда их перекрытие достигает значения, соответствующего нагрузке ас.
Для того чтобы учесть генерацию повреждений на масштабном уровне, меньшем чем размер автомата, линейная функция оклика должна быть модифицирована, Простейшая такая модификация показана на рис, 3, б, где приведён вид функции отклика типичной для бетонов. Хорошо известно,
а
а
а
б
г
г
а
н
а
0
8
8
Рис. 3: Примеры функций отклика: а^ - предел упругого поведения (точка
что для подобных материалов начиная с некоторой нагрузки (выше точки деградации эффективный модуль упругости понижается в результате генерации микроповреждений. Это означает, что линейный отклик имеет место при нагрузке в диапазоне [0, а^] в диапазоне же [а^, ас] генерируются повреждения и функция отклика ведёт себя нелинейно. Кроме того, по этой же причине разгрузка определяется новым модулем упругости (см. пунктирные линии на рис. 3, б). Из рис. 3, б хорошо видно, что каждой точке такой функции отклика (например, точкам 1 и 2) из диапазона [а^, ] соответствует свой модуль упругости, который и будет определять поведение материала при разгрузке.
Функции отклика показанные на рис. 3, а и б могут быть использованы для моделирования разрушения хрупких материалов типа керамики, бетона, и т. п.
Примеры функций отклика для необратимого поведения материала приведены на рис. 3, в и г. Случай, показанный на рис. 3, в, соответствует пластической деформации, а функция отклика, приведённая на рис. 3, г — комбинации пластического течения и процессов деградации материала (генерации микроповреждений). Следует подчеркнуть, что примеры, показанные на рис. 3, в и г носят иллюстративный характер. Определение функции отклика для описания необратимого поведения, особенно пластической деформации, является сложной задачей и подразумевает достаточно глубокие знания о закономерностях инициирования и развития пластического течения [27, 26].
Таким образом, для моделирования поведения образцов, структур, кон-
деградации); ас — прочность.
струкций (или их частей) в рамках метода подвижных клеточных автоматов необходимо знание конкретных функций отклика, примеры которых показаны па рис, 3, В случае композитных систем необходимо знание функций отклика для всех составляющих компонент. Подобная информация может быть получена как экспериментальным путем, так и па основе теоретического рассмотрения,
6. Метод отрезков и описание взаимодействия контактирующих поверхностей
При моделировании контактирующих поверхностей в рамках метода частиц возникает серьёзная проблема, связанная с корректным заданием их структуры и взаимодействия. Очевидно, что при рассмотрении тела как па-бора частиц, поверхности являются искусственно шероховатыми. Характеристики шероховатости, как можно видеть из рис, 4, связаны с размерами элементов.
Рис, 4: Искусственная шероховатость поверхностей двух соприкасающихся тол, рассматриваемых как набор частиц,
В общем случае при описании реальной поверхности можно выделить три масштабных уровня: макро- (масштаб моделируемого тона как целого), микро- (масштаб отдельных «бугорков», «впадин», повреждений и небольших трещин) и мезо- (промежуточный уровень, комбинирующий особенности макро- и микромасштабных уровней) 128, 29|,
Следовательно, существует два способа дня моделирования поверхности контакта:
1. Задать размеры автоматов настолько малыми, чтобы непосредственно (топологически) описать реальный рельеф поверхности. Это позволяет исследовать процессы на всех трёх перечисленных масштабных уровнях, В этом случае поверхностная шероховатость, неоднородность материала, другие структурные дефекты тина повреждений, трещин и прочих задаются явным образом через топологию поверхностного слоя.
2, Задать размеры автоматов достаточно большими дня того, чтобы рассматривать контактирующие поверхности как плоские с косвенным заданием шероховатости через шероховатость меньшого масштабного
уровня. Такое приближение позволяет охватить макро- и мезомасштаб-ные уровни. В этом случае, взаимодействие двух контактирующих поверхностей можно описать путем введения сил, зависящих от микропараметров поверхностей.
Даже с учётом мощностей современных компьютеров первый способ позволяет исследовать только малые области контактирующих поверхностей (см. рис. 5, а), поскольку для описания реальных материалов необходимо огромное число элементов. Следовательно, такое приближение может быть использовано для изучения механизмов очень сложных процессов в поверхностном слое материала, например возникновение и накопление повреждений, распространение трещин, изменение поверхностного профиля, структуры и состава материала возле поверхности и другие. Результаты подобных исследований могут быть очень полезными для понимания причин и тенденций явления на мезомасштабном уровне поскольку характерные размеры моделируемой системы соответствуют данному масштабу. Но этого недостаточно для анализа и описания взаимодействия между двумя поверхностями для реального (макро-) масштабного уровня.
\б\ мезомасштаб
VkM
Рис. 5: Два способа описания контактирующих поверхностей: на микромасштабном уровне шероховатость поверхности задается явным образом (а); опосредованное (мезомасштабное) описание шероховатости на основе приближения отрезков в методе МСА (б).
В рамках второго способа описания контактирующие поверхности рассматриваются как совокупность отрезков в двумерном случае и площадок в трёхмерном. В двумерном случае для описания взаимодействия между подвижными автоматами контактирующих поверхностей был разработан новый подход названный методом отрезков. Схема данного подхода приведена на рис. 5,6. Здесь автоматы in j соответствуют нижней контактирующей
поверхности, автомат к соответствует верхней поверхности. Пара г] рассматривается как отрезок, взаимодействующий с автоматом к. Вводятся два типа взаимодействия между автоматом и отрезком: нормальное (действующее по нормали к отрезку) потенциальное взаимодействие и тангенциальная сила трения.
Нормальное взаимодействие определяется высотой к между автоматом к и отрезком г — ]. Нормальная сила прикладывается к каждому автомату в треугольнике к — г — ] в следующем соотношении. Сила, действующая на автомат к со стороны отрезка г —
/ 1 к(^) т^к(ц) \
р на (Н^)) = Еуп ) . (10)
Сила, действующая со стороны автомата к на автомат г как часть отрезка г — Г-
Р]3 / ь^Л _
Р(н^)=щн- ). (11)
Сила, действующая со стороны автомата к на автомат ] как часть отрезка г — Г-
Е1г — кЦ(гЛ'
р(Л)к (кк(гз)\ = ^ I Нп — )
V п ) ^ )'
т К, I I ^ )
Но
Равновесное расстояние (высота) к0 для взаимодействия «отрезок - автомат» определяется как И0 = \/3Н где К — радиус автомата. Можно видеть, что Рк(гз) = р Ш)к + р (3г
Тангенциальное взаимодействие двух контактирующих поверхностей описывается диссипативпой силой трения, В рамках предложенного подхода ме-зофрагменты соприкасающихся поверхностей рассматриваются как плоские. Тогда тангенциальное взаимодействие «отрезок - автомат» можно записать следующим образом (по аналогии с (10)—(12)):
^(Уг ) = Я,, (13)
^)к(Уг) = ^, С")
тр 11
= ^(15)
где ^ — локальная сила трения.
Сопротивление относительному тангенциальному смещению определяется характеристиками шероховатости поверхностей на микромасштабном уровне, В данной работе общее выражение для локальной силы трения на мезомасштабном уровне имеет вид:
^ = ^оУт + / ),
(16)
где первое слагаемое в правой части описывает вязкое трение, а второе — трение, как результат взаимодействия контактирующих областей поверхностей, Здесь Ут — относительная тангенциальная скорость; — коэффициент вязкого трения; $ — некоторая функция; — коэффициент трения при контакте сухих поверхностей; — локальное значение силы нормального сжатия; а0 — предел упру гости; К^ — специальная переменная, эффективно отражающая локальное сопротивление тангенциальному смещению контактирующих поверхностей благодаря их зацеплению на микроуровне,
В общем случае параметры Ки являются сложными функци-
ями различных факторов микроуровня, таких как локальный профиль контакта поверхностей ж, локальное состояние поверхностей и история их взаимодействия Ь, относительная тангенциальная скорость Ут локальное значение нормального сжимающего напряжения а^са1 локальная температура Т1оса1 и т. д. Переменная К^ также является функцией средней площади контакта Б.
Для учёта зависимости перечисленных параметров уравнения (16) от микро- и мезоуровневых переменных могут быть использованы различные теоретические подходы и результаты экспериментов. Использование различных форм зависимости позволит описать и исследовать в деталях различные эффекты, возникающие в процессе трения, на мезо- и макромасштабном уровнях [29, 30],
Таким образом, предложенная модель даёт возможность непосредственно изучать взаимодействие контактирующих поверхностей на двух последовательных масштабных уровнях: мезо- и макро-, а также непосредственно исследовать влияние параметров микроуровня.
Использование мезоекопичееких моделей трения и износа является очень важным для получения практически важной обобщенной информации об изменениях поверхности (например, сглаживании рельефа) и характеристиках трения.
На практике часто используются макрохарактеристики трения, такие как коэффициент трения двух поверхностей
(17)
где — сила сопротивления контактирующих тел относительному тангенциальному движению; ^ — приложенная нормальная сила, В зависимости от мезомодели трения можно получить различные формы зависимости от времени и частоты трения, пройденного пути, начальной температуры и состояния поверхностей, а также других используемых на практике характеристик.
7. Эффективный учёт микрошероховатости поверхности
Для переменной К^, описывающей микрошероховатость, очень сложная реальная зависимость К^(х) (см. рис. 6, а) аппроксимируется суперпозицией двух кусочно-периодических функций (см. рис. 6,6). Периодические фрагменты аппроксимирующей функции описываются фрагментами синусоид
А™
Кп — — бЫ -К
( . ( х^-12\
\*—-г) + Т 81Ч ^-г--)
лп ^2
Ж?- 12
^с — Х> х^ — X/
п-1 1
п-1
+ кп Аг, + Кп Аг.
(18)
Здесь п — номер временного шага; Ацу и 11(2) — амплитуда н полупериод синусоиды, соответственно; х — текущее значение координаты на оси, направленной вдоль линии контакта поверхностей (тангенциальное направление), индексы 1 и 2 отмечают номер синусоиды; А — размер представительного мезообъема.
Рис. 6: Действительная (а) и аппроксимирующая (б, в) зависимости мезоко-эффициента Км от характеристик зацепления поверхностей на микроуровне (то есть, рельефа контакта).
Первая из синусоид (номер 1) в выражении (18) рассматривается как первое (грубое) приближение реальной зависимости Км (х). Вторая синусоида позволяет получать более точное (тонкое) приближение к реальному профилю (см. рис. 6,б).
Каждая синусоидальная компонента сложной кривой Км(х) имеет свое собственное значение амплитуды А, которое коррелирует с заданными средними значениями амплитуд первой и второй синусоид (см. рис. 6, в).
2
Для оценки характерной длины регулярного (периодического) участков кривой К^(х) были введены параметры регулярности d\ ъ d2 (см, рис, 6, в). Каждый периодический (синусоидальный) участок кривой имеет свою собственную длину, которая коррелирует со средним заданным значением dl(2).
Очевидно, что значения параметров dl(2) как и АЦ^ и £1(2) могут быть получены из экспериментов или теоретических оценок на основе других моделей [28, 30, 31],
В общем случае кривая К^(х) колеблется около некоторой средней величины большей нуля. Эта величина не зависит от FN и определяется только профилем контакта самих поверхностей. Вообще говоря, условие К^ > 0 определяется тем, что на мезоуровне на характерной длине R укладывается множество элементарных единиц поверхности («выпуклостей», «впадин» и т, д.). Так на микроуровне при рассмотрении взаимодействия двух сцепленных «выпуклостей» неизбежно реализуется механизм stick-slip (прерывистого скольжения) и, следовательно, здесь К^ = 0, При взаимодействии двух клеточных автоматов на мезоуровне вне зависимости от взаимодействия элементарных единиц всегда будет существовать некоторое зацепление по-
R
элементарной единицы поверхности), и, следовательно, Ктт > 0, Также за-
R
Км (по отношению к К^, будет уменьшаться до нуля.
Таким образом, при увеличении отношения размера клеточного автомата R к характерному размеру элементарных единиц поверхности величина К^
R
макромаештабный уровень. Относительное значение амплитуды колебаний при этом должно спадать до нуля.
8. Заключение
В работе изложены основы метода подвижных клеточных автоматов и его отличия от других методов частиц, предназначенных для моделирования деформации материалов и сред. Главным существенным отличием МСА является запись уравнений движения с учётом многочастичности взаимодействия элементов. Подробно изложены способы описания поверхностей контактирующих тел в рамках метода МСА, Следует отметить, что данный метод не является закрытым и позволяет применять различные подходы и модели для описания моделируемых сред. Так, в зависимости от функции отклика клеточного автомата, существует возможность описания процессов в твёрдом теле как на мезомасштабном уровне (с явным учётом анизотропии представительных объёмов), так и на макроуровне (в приближении изотропности автоматов),
Важным преимуществом метода МСА в сравнении с методами механики сплошной среды является возможность прямого моделирования процессов разрушения. Такая возможность напрямую вытекает из постулатов метода,
поэтому не нуждается в искусственных построениях,
К сожалению, в рамках одной статьи невозможно охватить все составляющие метода МСА, В частности, не было рассмотрено использование этого метода для многоуровневого моделирования структурированных материалов и сред. Следует отметить, что в настоящее время созданы и апробированы два основных подхода. Первый основан на совместном использовании метода МСА (мезоуровень) с численными методами континуальной механики (макроуровень). Второй подход предполагает поэтапное моделирование методом МСА представительных элементов всех структурных уровней. Для этого проводятся серии расчётов начиная с самого нижнего уровня, где все структурные элементы учитываются явно. На основе анализа результатов этих расчётов определяются функции отклика подвижных клеточных автоматов следующего структурного уровня и функции случайного распределения параметров функций отклика. Полученная таким образом информация позволяет проводить серии расчётов на представительных образцах следующего структурного уровня, и так далее. Подобный подход был успешно применён для моделирования наноетруктурной керамики и металлокерами-ческих композитов.
Необходимо отметить, что перспективы развития метода МСА являются достаточно широкими. Так, ведётся работа по развитию на основе МСА так называемого метода гибридных клеточных автоматов, предназначенного для решения двумерных и трёхмерных задач деформирования флюидо- и газонасыщенных сред.
Таким образом, благодаря своим преимуществам метод МСА открывает новые возможности для непосредственного моделирования таких сложных процессов, как перемешивание масс, возникновение и накопление повреждений, распространение трещин, фазовые переходы и химические реакции и т, д.
СПИСОК ЦИТИРОВАННОЙ ЛИТЕРАТУРЫ
1, Darrigol О, Between hydrodynamics and elasticity theory: the first five births of the Navier-Stokes equation // Arch, Hist, Exact Sei, 2002, Vol, 56, P. 95150.
2, Левин В, А, Нелинейная вычислительная механика прочности, Т. 1, Модели и методы. Образование и развитие дефектов. Под ред. В, А, Левина, М,: Физматлит, 2015, 454 с,
3, Левин В, А,, Вершинин A.B. Нелинейная вычислительная механика прочности, Т. 3, Численные методы. Реализация на высокопроизводительных вычислительных системах. Под ред. В, А, Левина, М,: ФИЗМАТЛИТ, 2015. 543 с.
4, Морозов Е, М,, Левин В, А,, Вершинин А. В, Прочностной анализ, Фц ionic в руках инженера, М,: UEES, 2015, 400 с.
5, Cundall P. A,, Straek O.D.L. A discrete numerical model for granular assemblies // Geoteehnique, 1979, Vol, 29, No, 1, P. 47-65,
6, Cundall P. A, A computer simulations of dense sphere assembles // Micromechanics of granular materials / Eds, by M, Satake and J.T. Jenkins, Amsterdam: Elsever Sci. Publ., 1988. P. 113-123.
7, Herrmann H, J, Simulating granular media on the computer // 3rd Granada lectures in computational physics / Eds, by P.L. Garrido and J, Marro, Heidelberg: Springer, 1995, P. 67-114,
8, Hemmingsson J,, Herrmann H, J,, Roux S, On stress networks in granular media //J. Phvs. I. 1997. Vol. 7. P. 291-302.
9, Walton О. E. Numerical simulation of inclined chute flows of monodisperse, inelastic, frictional spheres // Mechanics of Materials. 1993. Vol. 16. P. 239247.
10. Luding S. Granular materials under vibration: Simulations of rotating spheres // Phvs. Rev. E. 1995. Vol. 52. No. 4. P. 4442.
11. Poschel T. Granular material flowing down an inclined chute: A molecular dynamic simulation //J. Phvs. II. 1993. V. 3. P. 27.
12. Greenspan D. Particle modeling in science and technology // Coll. Math. Societatis Janos Bolvai, 1988. №. 50. P. 51.
13. Ostermaver G. P. Friction models with discrete layers // Z. Angew. Math. Mech. 2000. Vol. 80. P. 61-64.
14. Mustoe G. G. W. A generalized formulation of the discrete element method // Engineering computations. 1992. Vol. 9. P. 181-190.
15. Bicanic N. Discrete element methods. In: E. Stein, R. Borst, T.J.R. Hughes (Eds.), Encyclopedia of computational mechanics. Vol. 1: Fundamentals, Wiley, Chichester, 2004. P. 311-337.
16. Хокни P., i In нуд Дж, Численное моделирование методом частиц. М,: Мир, 1987. 640 с.
17. Поттер Д. Вычислительные методы в физике. М,: Мир, 1975. 218 с.
18. Псахье С. Г., Хори Я., Коростелев С. К).. Смолин А. К).. Дмитриев А. И., Шилько Е.В., Алексеев С. В. Метод подвижных клеточных автоматов, как инструмент для моделирования в рамках физической мезомеханики // Изв. вузов. Физика. 1995. JV2 11. С. 58-69.
19. Псахье С. Г., Коростелев С.Ю., Смолин А. К).. Дмитриев А. И., Шилько Е.В., Моисеенко Д. Д., Татаринцев Е.М., Алексеев С. В. Метод подвижных клеточных автоматов как инструмент физической мезомеханики материалов // Физическая мезомеханика. 1998. Т. 1. № 1. С. 95-108.
20. Дмитриев А, И,, Короетелев С.Ю., Оетермайер Г, П., Пеахье С, Г., Смолин А. К).. Шилько Е, В, Метод подвижных клеточных автоматов, как инструмент для моделирования на мезоуровне // Известия РАН, Механика твердого тела, 1999, JV2 6, С, 87-94,
21. Psakhie S,, Shilko Е,, Smolin A,, Astafurov S,, Oveharenko V, Development of a formalism of movable cellular automaton method for numerical modeling of fracture of heterogeneous elastic-plastic materials // Frattura ed Intégrité Strutturale, Vol. 24. 2013. P. 26-59.
22. Shilko E.V., Psakhie S.G., Schmauder S,, Popov V. L., Astafurov S.V., Smolin A, Yu, Overcoming the limitations of distinct element method for multiscale modeling of materials with multimodal internal structure // Computational Materials Science, 2015, Vol, 102, P. 267-285,
23. Mikhailov A. S. Foundations of Synergetics I. Distributed Active Systems. Berlin: Springer, 1994. 112 p.
24. Смолин А. Ю., Роман H. В., Добрынин С. А., Пеахье С. Г. О вращательном движении в методе подвижных клеточных автоматов // Физическая мезомеханика. 2009. Т. 12, № 2. С. 17-22.
25. Foiles S. M., Baskes M. I., Daw M.S. Embeded-atom-method functions for the f.c.c. metals Cu, Ag, Au, Ni, Pd, Pt and their alloys // Physical Review B. 1986. Vol. 33. No. 12. P. 7983-7991.
26. Psakhie S.G., Zolnikov K.P., Krvzhevich D.S., Protodefect as a basis of multilevel nanoscale plasticity of crystal materials // American Institute of Physics. 2008. Vol. 999 P. 20-31.
27. Корчуганов А.В., Зольников К.П., Крыжевич Д.С., Чернов В.М., Пеахье С.Г. Моделирование зарождения пластической деформации в механичееки-нагруженных кристаллитах при радиационном воздействии // ВАНТ. Сер. Термоядерный синтез. 2015. Т. 38, Вып. 1. С. 42-48.
28. Крагельекий И. В. Трение и износ. М,: Машгиз,, 1962. 384 с.
29. Tworzydlo W. W,, Cecot W,, Oden J, T., Yew С, H, Computational miero-and macroscopic models of contact and friction: formulation, approach and applications // Wear. 1999. Vol. 220. P. 113-140.
30. Eozman M.G., Urbakh M,, Klafter J. Stick-slip dynamics of interfacial friction // Phvsica A. 1998. Vol. 249. P. 184-189.
31. Eaharijaona F,, Eoizard X., Stebut J. Usage of 3D roughness parameters adapted to the experimental simulation of sheet-tool contact during a drawing operation // Tribologv International, 1999, Vol, 32, P. 59-67,
REFERENCES
1. Darrigol, O, 2002, "Between hydrodynamics and elasticity theory: the first five births of the Navier-Stokes equation", Arch. Hist. Exact Sci., vol, 56, pp. 95-150.
2. Levin, V. A. 2015, Nonlinear computational 'mechanics of strength. V. 1. Models and methods. Nucleation and development of defects. Eds. V. A. Levin, Phvsmatlit, Moscow. (inRussian)
3. Levin, V. A. & Vershinin, A. V. 2015, Nonlinear computational 'mechanics of strength. V. 3. Numerical methods. Implementation for high-performance computers. Eds. V. A. Levin, Phvsmatlit, Moscow. (inRussian)
4. Morozov, E. M., Levin, V. A. & Vershinin, A. V. 2015, Strength analysis. Fidesis for engineers, URRS, Moscow. (inRussian)
5. Cundall, P. A. & Strack, O. D. L. 1979, "A discrete numerical model for granular assemblies", Geotechnique, vol. 29. No. 1. pp. 47-65.
6. Cundall, P. A. 1988, "A computer simulations of dense sphere assembles", Micromechanics of granular materialsEds. Eds. by M. Satake and J. T. Jenkins. Amsterdam: Els ever Sci. Publ, pp. 113-123.
7. Herrmann, H. J. 1995, "Simulating granular media on the computer", 3rd Granada lectures in computational physics. Eds. by P.L. Garrido and J. Marro. Heidelberg: Springer, pp. 67-114.
8. Hemmingsson, J., Herrmann, H.J. & Roux, S. 1997, "On stress networks in granular media", J. Phys. J, vol. 7. pp. 291-302.
9. Walton, O. R. 1993, "Numerical simulation of inclined chute flows of monodisperse, inelastic, frictional spheres", Mechanics of Materials, vol. 16. pp. 239-247.
10. Luding, S. 1995, "Granular materials under vibration: Simulations of rotating spheres", Phys. Rev. E, vol. 52. No. 4. pp. 4442.
11. Posehel, T. 1993, "Granular material flowing down an inclined chute: A molecular dynamic simulation", J. Phys. II, vol. 3. pp. 27.
12. Greenspan, D. 1988, "Particle modeling in science and technology", Coll. Math. Societatis Janos Bolyai., №. 50. pp. 51.
13. Ostermaver, G. P. 2000, "Friction models with discrete layers", Z. Angew. Math. Mech,., vol. 80. pp. 61-64.
14. Mustoe, G. G. W. 1992, "A generalized formulation of the discrete element method", Engineering computations, vol. 9. pp. 181-190.
15, Bieanie, N. 2004, "Discrete element methods", In: E. Stein, R. Borst, T.J.R. Hughes (Eds.), Encyclopedia of computational 'mechanics. Vol. 1: Fundamentals, Wiley, Chichester, pp. 311-337,
16, Hockney, E, W, & Eastwood, J.W, 1981, Numerical Simulation by the Particle Method, McGraw-Hill International Book Company,
17, Potter, D, 1973, Computational Physics, John Wiley & Sons Ltd., London,
18, Psakhie, S.G., Horie, Y,, Korostelev, S.Yu,, Smolin, A.Yu,, Dmitriev, A,I,, Shilko, E.V, & Alekseev, S.V., 1995, "Method of movable cellular automata as a tool for simulation within the framework of mesomeehanies", Russian Physics Journal, vol, 38, №11, pp. 1157-1168,
19, Psakhie, S.G., Korostelev, S.Yu,, Smolin, A.Yu,, Dmitriev, A, I,, Shilko, E.V., Moisevenko, D,D,, Tatarintsev, E, M, & Alexeev, S.V, 1998, "Movable cellular automata method as a tool for physical mesomeehanies of materials", Physical Mesomeehanies, vol, 1, No, 1 pp. 89-102,
20, Dmitriev, A, I,, Korostelev, S.Yu,, Ostermever, G.P., Psakhie, S.G., Smolin, A, Yu. & Shilko, E.V. 1999, "Movable cellular automata method as a tool for simulation at the mesolevel", Proc. of RAS. Mechanics of solids, vol. 6. pp. 87-94.
21, Psakhie, S,, Shilko, E,, Smolin, A., Astafurov, S. & Oveharenko, V. 2013, "Development of a formalism of movable cellular automaton method for numerical modeling of fracture of heterogeneous elastic-plastic materials", Frattura ed Integrita Strutturale, vol, 24, pp. 26-59,
22, Shilko, E. V., Psakhie, S.G., Schmauder, S,, Popov, V.L., Astafurov, S.V. & Smolin, A.Yu, 2015, "Overcoming the limitations of distinct element method for multiscale modeling of materials with multimodal internal structure", Computational Materials Science, vol. 102. pp. 267-285.
23, Mikhailov, A. S. 1994, Foundations of Synergetics I. Distributed Active Systems. Springer, Berlin.
24, Smolin, A.Y., Eoman, N.V., Dobrvnin, S.A. & Psakhie, S.G. 2009, "On rotation in the movable cellular automaton method", Physical Mesomeehanies vol. 12. № 2. pp. 17-22.
25, Foiles, S.M., Baskes, M.I. & Daw, M.S. 1986, Embeded-atom-method functions for the f.c.c. metals Cu, Ag, Au, Ni, Pd, Pt and their alloys", Physical Review B, vol. 33. No. 12. pp. 7983-7991.
26, Psakhie, S.G., Zolnikov, K.P. & Krvzhevieh, D.S. 2008, "Protodefeet as a basis of multilevel nanoscale plasticity of crystal materials", American Institute of Physics, vol. 999. pp. 20-31.
27. Korchuganov, А. V,, Zolnikov, К, P., Kryzhevich, D, S, Chernov, V, М. & Psakhie, S, G, 2016, "MI) simulation of plastic deformation nucleation in stressed crystallites under irradiation", Physics of Atomic Nuclei, vol, 79, iss. 7. pp. 1193-1198.
28. Kragelskii, I. V. 1962, Friction and Wear, Mashgiz, Moscow.
29. Tworzvdlo, W. W., Cecot, W., Oden, J. T. & Yew, С. H. 1999, "Computational micro- and macroscopic models of contact and friction: formulation, approach and applications", Wear, vol. 220. P. 113-140.
30. Eozman, M. G,, Urbakh, M. & Klafter, J. 1998, "Stick-slip dynamics of interfacial friction", Physica A, vol. 249. pp. 184-189.
31. Eaharijaona, F,, Eoizard, X. & Stebut, J. 1999. "Usage of 3D roughness parameters adapted to the experimental simulation of sheet-tool contact during a drawing operation", Tribology International, vol. 32, pp. 59-67.
Получено 22.05.2017
принято в печать 14.09.2017
