О вращательном движении в методе подвижных клеточных автоматов Текст научной статьи по специальности «Физика»

УДК 531.133, 531.16, 539.37

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

А.Ю. Смолин1,2, Н.В. Роман2, С.А. Добрынин1, С.Г. Псахье1,2

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

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

Ключевые слова: моделирование, метод частиц, механическое движение, вращение

On rotation in the movable cellular automaton method

A.Yu. Smolin1,2, N.V. Roman2, S.A. Dobrynin1 and S.G. Psakhie1,2

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

2 Tomsk State University, Tomsk, 634050, Russia

The paper deals with the problem of rotation consideration in the framework of the particle method. The rotation taken either as an independent degree of freedom or as rotation of the surrounding allows the movable cellular automaton method to correctly describe a classical continuum. For more complex media, e.g., micropolar ones, both ways of rotation consideration should be realized with the independent degree of freedom taken as an independent rotation.

Keywords: simulation, particle method, mechanical motion, rotation

1. Введение

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

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

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

© Смолин А.Ю., Роман Н.В., Добрынин С.А., Псахье С.Г., 2009

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

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

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

2. Учет вращения в методе подвижных клеточных автоматов как самостоятельной степени свободы

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

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

нием вокруг оси пары (Юу • г у = 0), то для Юу получим следующее выражение:

у

ю у =-

Г у х (у.) - V I) = п у х (V у -

у у

где Пу — единичный вектор направления от центра автомата i к центру автомата j.

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

ю у =■

qц

qц

= (юу - ю)х пу

Угловая величина сдвиговой деформации каждого из автоматов в момент времени Т определится интегралом

Уу =1юУЬ(0& = 0у -0„

0

(1)

где 0у — вектор поворота оси пары автоматов у относительно начального положения.

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

Ц = -(°7 у + &уу у,) х пу > (2)

где О, — модуль сдвига материала автомата i. В выражении (2) представлена парная тангенциальная сила, обусловленная сдвиговой деформацией каждого автомата пары.

Из выражений (1) и (2) видно, что при вращении элементов с одинаковой скоростью, но в противоположных направлениях (0, = -0у), силы сопротивления сдвигу не будут препятствовать такому движению. Однако среда, описываемая автоматами, в данном случае будет испытывать деформацию и должна ей сопротивляться. В случае линейно-упругого поведения момент сил, препятствующих относительному вращению, должен быть прямо пропорционален разности векторов поворота элементов в паре 0, - 0у. Возникает вопрос: Чему должен быть равен коэффициент этой пропорциональности, какому упругому модулю? Если рассматривать пару элементов как брусок, то такая деформация

где VА и V0 — линейные скорости соответственно точек А и О в произвольной инерциальной системе координат; г0А — вектор, проведенный из точки О в точку А. В случае вращения пары автоматов как целого (без деформации, в том числе сдвиговой) это соотношение можно переписать следующим образом (рис. 1):

vу - vi = ю у х Гу,

где Гу = Гу - г, ; Юу — мгновенная скорость вращения пары как целого. Если умножить векторно слева обе части последнего уравнения на Гу и пренебречь враще-

Рис. 1. Вращение пары и вращение каждого автомата

v Л = О + Ю х Г0А.

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

3. Учет вращения как согласованное движение элементов среды

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

= 2 юу/к

где Юу — скорость вращения оси пары автоматов у'; N— количество взаимодействующих соседей.

В трехмерном случае повороты описываются ортогональными тензорами или векторами, имеющими три независимые компоненты. Задачу об определении среднего вращения можно сформулировать так. Имеются несколько точек у = 1, ..., N определяемых векторами Гу, проведенными из неподвижной точки I. Каждая точка у движется относительно точки I со скоростью Vij. Такому движению можно сопоставить вращение некоторого твердого тела, образуемого у-ми точками, вокруг точки I. Реальное движение каждой точки у будет складываться из среднего движения точки как твердого тела

и деформации этого тела. Нужно определить скорость вращения такого твердого тела Шу.

Как известно [12], кинетический момент твердого тела при вращении вокруг оси, заданной единичным вектором п Ю, определяется как

= 2 ту 'ОуР у = 2 туЮр2 = Ую,

К = 2 гу х ту V у

у у у

где ру = |п юх гу — расстояние от оси вращения до точки у; 3— момент инерции тела относительно оси п ю. Тогда требование, чтобы движение системы точек у = 1, ..., N в среднем соответствовало вращению некоторого твердого тела, можно записать как равенство соответствующих кинетических моментов:

ю 2шур2 = 2 Гу х ту Vу.

Учитывая, что массы всех автоматов одинаковы и Гу хVу = ГуЮу, искомая средняя скорость вращения окружения автомата в трехмерном случае запишется как

Ш,

= 2

у

'2 |п

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

4. Преимущества учета вращения как самостоятельной степени свободы

Несмотря на сходство получаемых результатов следует отметить, что явный учет вращений автоматов позволяет моделировать более сложные среды с неявным учетом их структуры. Например, в работе [14] показано, что тензор напряжений, в общем случае, для таких моделей является несимметричным, как и должно быть в случае среды Коссера со стесненным вращением [7]. Кроме того, на основе расчетов методом дискретных элементов (с явным учетом вращений) можно построить определяющие соотношения для микрополярной среды, описывающей движение гранулированных материалов в рамках континуума [15].

Учет поворота одновременно двумя рассмотренными способами дает возможность реализации среды Кос-сера в методе подвижных клеточных автоматов, если при этом собственное вращение автоматов рассматривать как независимое. Естественно, что такой подход приведет к необходимости введения дополнительных упругих модулей, которые рассматриваются в теории микрополярных сред и, вообще говоря, могут быть определены экспериментально [7]. Как известно, в этом случае в среде возможно распространение новых типов упругих волн, в том числе так называемых волн продольного вращения. На рис. 2 показано распростране-

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

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

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

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

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

ранения упругих волн в полупространстве от точечного источника на его поверхности (задача Лэмба), полученные с учетом вращения элементов и без него (рис. 3). В силу симметрии задачи рассматривалась расчетная область, представляющая собой 1/4 куба, ее размеры составляли 0.25х0.25х0.25 м. Размер автомата — 0.0025 м, материал — сталь (р = 7 800 кг/м3, ур = = 5.95 км/с, vS = 3.19 км/с). На плоскостяхX = 0 и У = = 0 задавались условия симметрии, остальные грани куба считались свободными. Импульс прикладывался в точке, совпадающей с началом координат (показано на рис. 3, а стрелкой), в виде «купола» синусоиды длительностью 5 мкс. Все результаты на рис. 3, в-е представлены для 40-й мкс в плоскости У = 0.

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

Рис. 3. Поле скоростей подвижных клеточных автоматов при решении задачи Лэмба для различных вариантов модели и упаковок автоматов: а, е — кубическая упаковка с учетом вращения элементов, б — плотная упаковка на плоскости без учета вращения, в — ГЦК-упаковка без учета вращения, г — кубическая упаковка без учета вращения, д — ГЦК-упаковка с учетом вращения элементов

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

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

Вторым примером, показывающим важность учета вращений в методе подвижных клеточных автоматов, являются результаты теста на простой сдвиг, показанные на рис. 4. Здесь представлены данные только для ГЦК-упаковки автоматов, потому что соответствующие

Рис. 4. Схема деформирования модельного образца и диаграммы нагружения при простом сдвиге. Обозначения кривых соответствуют ГЦК-упаковке с учетом вращения элементов (1) и без учета вращения (2)

диаграммы для кубической упаковки совпадают с ними. Видно, что без учета вращений угол наклона рассчитываемых диаграмм в два раза превышает значение модуля сдвига G = 79 ГПа. Это объясняется тем, что в этом случае заданное смещение Ах целиком соответствует сдвиговой деформации г у = у = Ах^, в то время как согласно формуле для компонент тензора малых деформаций г у = 1/2(диу/ дх + дих/ ду) при простом сдвиге должно быть г у = 0.5у, т.е. при простом сдвиге образец испытывает сдвиговую деформацию равную 0.5у и одновременно поворачивается на такой же угол как целое. Действительно, в случае учета вращений все автоматы образца при его простом сдвиге оказываются повернутыми на угол равный 0.5у и сдвиговая деформация во всех парах имеет то же самое значение.

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

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

Работа выполнена при финансовой поддержке междисциплинарного интеграционного проекта СО РАН № 32, РФФИ (грант № 07-08-00192), а также проекта 13.13.3 Программы фундаментальных исследований ОЭМПУ РАН № 13.13.

Литература

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

2. Панин B.E. Основы физической мезомеханики // Физ. мезомех. -

1998. - Т. 1. - № 1. - С. 5-22.

3. Пcaxьe C.r., Дмumрueв A.И., Шилько E.B., Cмoлuн A.Ю., Коро^е-лев C.Ю. Метод подвижных клеточных автоматов как новое направление дискретной вычислительной механики. I. Теоретическое описание // Физ. мезомех. - 2000. - Т. 3. - № 2. - С. 5-15.

4. Пcaxьe C.r., Хори Я., Коро^елев C.Ю., Cмoлuн И.Ю., Дмum-риев A.И., Шилько E.B., Aлeкceeв C.B. Метод подвижных клеточных автоматов как инструмент для моделирования в рамках физической мезомеханики // Изв. вузов. Физика. - 1995. - Т. 38. -№11.- С. 58-69.

5. Пcaxьe C.r., Чeрmoв M.A., Шилько E.B. Интерпретация параметров

метода подвижных клеточных автоматов на основе перехода к континуальному описанию // Физ. мезомех. - 2000. - Т. 3. - № 3. -С. 93-96.

6. Cosserat E., Cosserat F. Theorie des Corps Deformables. - Paris: A. Hermann et fils, 1909. - 226 p.

7. Eрoфeeв B.И. Волновые процессы в твердых телах с микрострукту-

рой. - М.: Изд-во МГУ, 1999. - 328 с.

8. Пcaxьe C.r., Cмoлuн A.Ю., Cmeфaнoв Ю.П., Maкaрoв ПМ., Шилько E.B., Чeрmoв M.A., Eвmyшeнкo E.П. Моделирование поведения сложных сред на основе комбинированного дискретно-континуального подхода // Физ. мезомех. - 2003. - Т. 6. - № 6. - С. 1121.

9. Уилкинс MM. Расчет упругопластических течений // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 212-263.

10. Heгрескул C.И., Пcaxьe C.r., Коро^елев C.Ю., Панин B.E. Моделирование зернистых сред методом элементной динамики. -Томск, 1989. - 27 с. / Препринт ТНЦ СО АН СССР № 39.

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

12. Ландау Л.Д., Лифшиц E.M. Теоретическая физика: В 10 т. Т.1. Механика. - М.: Наука, 1988. - 216 c.

13. Aлeкcaндрoв A.B., Пomanoв B.Д., Державин Б.П. Сопротивление материалов: Учебн. для вузов. - М.: Высш. школа, 2001. - 560 с.

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

15. Walsh S.D.C., Tordesillas A. A thermomechanical approach to the development of micropolar constitutive models of granular media // Acta Mech. - 2004. - V. 167. - No. 3-4. - P. 145-169.

16. Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических явлений. - М.: Наука, 1966. -686 с.

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

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

Смолин Алексей Юрьевич, к.ф.-м.н., доцент, старший научный сотрудник ИФПМ СО РАН, доцент ТГУ, asmolin@ispms.tsc.ru

Роман Никита Витальевич, магистрант ТГУ, nroman@kvadro.net

Добрынин Сергей Александрович, аспирант ИФПМ СО РАН, laesod@usgroups.com

Псахье Сергей Григорьевич, д.ф.-м.н., профессор, директор ИФПМ СО РАН, профессор ТГУ, sp@ispms.tsc.ru