CC BY
54
12. Leawitt R.P. Pressure broadening and shifting in microwave and infrared spectra of molecules of arbitrary symmetry: An iireducible tensor approache. // J. Chem. Phys. - 1980. - V. 73. -P. 5432—5450.
13. Черкасов М.Р. Формализм квантовомеханического оператора Лиувилля в расчетах релаксационных параметров // Томск: ИОА СО АН СССР, 1975. - 47 c. (Препринт № 26).
14. Ma Q., Tipping R.H., Boulet C. Irredusible correlation functions of the S matrix in the coordinate representation: Application in calculating Lorentzian half-widths and shifts // J. Chem. Phys. - 2006. -V. 124. - 014109 (14 pages).
15. Варшалович Д.А., Москалев А.Н., Херсонский В.К. Квантовая теория углового момента. - Л.: Наука, 1975. - 436 c.
16. Bykov A.D., Stroinova V.N., Smekalina E.L. Application ofvariatio-nal technique to relaxation parameters calculation of highly vibratio-nally excited CO molecule // Atomic and molecular pulsed lasers: Abstracts of papers of 8th Intern. Conf. - Tomsk, 2007. - P. 108.
17. Bykov A.D., Stroinova V.N., Smekalina E.L. Application of variational technique to relaxation parameters calculation of highly vibratio-nally excited CO molecule // Atmospheric and ocean optics: Proc.s of XIV Intern. Symp. - Moscow, 2007. - V. 6936. - 693604 (5 pages).
18. Булдаков М.А., Черепанов В.Н. Полуэмпирические функции дипольного момента молекул СО и NO // Оптика атмосферы и океана. - 2004. - Т. 17. - № 1. - С. 42-46.
19. Maroulis G. Accurate higher electric multipole moments for carbon monoxide // J. Chem. Phys. Lett. - 2001. - V. 334. - P. 214-219.
20. Maroulis G. Electric polarizability and hyperpolarizability of carbon monoxide // J. Phys. Chem. - 1996. - V. 100. - P. 13466-13473.
21. Gordov E.P., Lykosov V.N., Fazliev A.Z. Web portal on environmental sciences «ATMOS» // Advances in Geophysical Sciences. -2006. - V. 8. - № 1. - P. 33-38.
Поступила 19.06.2008 г.
УДК 531.355
РАЗВИТИЕ ФОРМАЛИЗМА МЕТОДА ЧАСТИЦ ДЛЯ РАСЧЕТА УСЛОВИЙ НА ГРАНИЦЕ РАЗДЕЛА ТВЕРДОГО ТЕЛА С ЖИДКОЙ СРЕДОЙ
А.Ю. Панченко, Е.В. Шилько, С.В. Астафуров, С.Ю. Коростелев, С.Г. Псахье
Институт физики прочности и материаловедения СО РАН, г. Томск E-mail: lion@usgroups.com
Для численного моделирования методом частиц поведения твердых тел в жидкой среде предложен базовый формализм граничных условий, позволяющих эффективным образом осуществлять учет «механического» влияния жидкости. Несмотря на ряд ограничений, развитый алгоритм реализации условий на границе твердое тело - жидкость продемонстрировал применимость для решения широкого класса задач, связанных, в частности, с изучением поведения плитных сред различной природы (ледовых покровов водоемов, фрагментов литосферы и т. д.), покоящихся на жидкоподобном основании. При этом простота и эффективность предложенного подхода делает возможным реализацию граничных условий данного типа в рамках различных методов частиц, относящихся к классу методов дискретных элементов.
Ключевые слова:
Методы частиц, метод подвижных клеточных автоматов, граничные условия, граница «твердое тело - жидкость».
1. Введение
Интенсивное развитие вычислительной техники, наблюдаемое на протяжении последних трех десятилетий, стимулировало создание новых, а также развитие формализма существующих методов компьютерного моделирования процессов деформации и разрушения материалов на различных масштабных уровнях. В целом совокупность таких методов может быть условно разбита на два класса: методы механики сплошной среды (представлен -ные, главным образом, методом конечных разностей [1] и методом конечных элементов [2]) и методы частиц. Последние также развиваются в рамках двух основных направлений: методов динамики частиц, в рамках которого дискретный элемент определяется на основе использования функции распределения плотности [3], и методов дискретных элементов с явным учетом размеров и формы элементов (методы «кундалловского типа») [4].
Адаптация методов частиц для решения задач механики сложных гетерогенных материалов и сред в значительной степени связана с разработкой эффективных алгоритмов расчета условий на контактных границах твердого тела, моделируемого ансамблем элементов конечного размера. При этом особенности задания граничных условий определяются не только природой взаимодействующих тел, но и характеристиками (в частностями, скоростями) деформационных процессов.
Важным классом задач механики структурнонеоднородных сред является изучение поведения твердого тела (или ансамбля тел), частично или полностью погруженного в среду с нулевым или близким к нулю сопротивлением сдвигу. К таким системам относятся, например, фрагменты литосферы [5] или ледовые покровы водоемов [6], характеризующиеся наличием стесненных граничных условий в плоскости «залегания» (подобные среды
называют квазидвумерными блочными или плитными).
Следует отметить, что при решении динамических задач формализм граничных условий для твердых тел произвольной геометрии, контактирующих с жидкостью, получил свое развитие, главным образом, в рамках методов механики сплошной среды [2, 7]. В то же время, актуальным является перенос данной методологии на дискретные методы механики, в частности, методы «кундалловского типа» (в случае методов динамики частиц могут быть использованы подходы, изложенные в [3, 8]).
В настоящей работе предложен достаточно простой и универсальный алгоритм расчета условий на поверхности контакта твердого тела, моделируемого ансамблем частиц, с жидкой средой. Данный алгоритм учитывает силу Архимеда и диссипативную силу сопротивления жидкоподобной среды движению тела. Предложенные граничные условия реализованы в рамках формализма метода подвижных клеточных автоматов [9, 10], который можно рассматривать как представителя методов частиц «кундалловского типа».
2. Формализм граничных условий
для поверхностей погруженного твердого тела,
контактирующих с жидкой средой
2.1. Основные приближения
В рамках метода подвижных клеточных автоматов деформируемое твердое тело моделируется ансамблем дискретных элементов (клеточных автоматов). Подвижный клеточный автомат, являясь дискретным объектом «кундалловского типа», имеет конечные пространственные характеристики (размер, форму, массу, момент инерции и т. д.). В настоящей работе изложена реализация алгоритма граничных условий в двумерной версии метода. Рассмотрено приближение, в рамках которого исходное (недеформированное) тело описывается регулярной плотной упаковкой автоматов (рис. 1, а) одинакового исходного размера й на плоскости (реализация алгоритма для другого типа упаковки, а также распространение его на трехмерный случай производится по аналогии).
В исходном (недеформированном) состоянии автоматы имеют форму призмы высотой й и объема Уса=^3й 3/2 (рис. 1, б). В рамках данного рассмотрения ограничимся случаем хрупких твердофазных материалов, изменение формы и размеров которых в результате деформационных процессов является незначительным. Для описания упругопластического отклика подвижного клеточного автомата используется модель, аналогичная модели малых упругопластических деформаций [11].
Реализация предлагаемого алгоритма расчета условий на внешней поверхности погруженной в жидкость части твердого тела связана с использованием следующих «внешних» сил, действующих
на дискретные элементы (подвижные клеточные автоматы) в дополнение к потенциальным и вязким сил взаимодействия автоматов [10]:
1. Сила тяжести, ответственная за погружение тела. Данная сила действует на все элементы (подвижные клеточные автоматы), составляющие твердое тело: Fgrav=m.g, где т - масса автомата, g - ускорение свободного падения.
2. Гидростатическое давление жидкости, действующее только на автоматы, составляющие часть поверхности твердого тела, контактирующей с жидкостью.
3. Сила сопротивления, ответственная за диссипацию кинетической энергии тела, движущегося в жидкости. Диссипативная сила также действует только на поверхностные автоматы твердого тела, контактирующие с жидкостью. В первом приближении в качестве такой силы может использоваться только сила лобового сопротивления жидкости.
Рис. 1. Определение: а) размеров и б) объема подвижных клеточных автоматов в случае исходной плотной упаковки на плоскости
Отметим, что две последние силы действуют также на поверхностные автоматы твердого тела, контактирующие с воздушной средой, однако здесь не учитываются ввиду относительной малости.
2.2. Расчет сил гидростатического давления жидкости
на поверхностные элементы
Как объект конечного размера, клеточный автомат, моделирующий часть поверхности твердого тела, может быть погружен в жидкость полностью или частично. При этом на каждый малый участок поверхности автомата й£, взаимодействующей с жидкостью, действует гидростатическое давление жидкости Рн=рж-&Н (рж - плотность жидкости, Н -глубина), направленное по нормали к этому участку. Суммарная сила определяется как интеграл соответствующих сил й¥Н=рН;й8 по поверхности S взаимодействия автомата с жидкостью. В случае одиночного полностью погруженного подвижного клеточного автомата она равна ?Л°=рл;-4х3Я3.п (где К=й/2 - эффективный радиус автомата, Я - единичный вектор, направленный противоположно силе тяжести). В то же время, у автомата, моделирующего некоторую часть поверхности твердого тела, есть взаимодействующие с ним соседи (рассматриваемые пары автоматов могут находиться в связанном состоянии, либо несвязанном, но контактирующем [9]). На соответствующие «грани» (пло-
щадки контакта с «соседями») гидростатическое давление не действует. Поэтому результирующее значение силы гидростатического давления Р1оШ\ действующей на полностью погруженный^поверх-ностный автомат /, будет равно разности РА° и «несуществующих» вкладов АР/ (рис. 2):
Ff
= Ff -^aF' =
j=1
■Рж8
• 4л/зЛ3 ч
■Y^p^gSj
j=1
(У, - y,)
q + ff
m, (1)
где у1 и у - координаты центров масс взаимодействующих автоматов / и у; Н>0 - глубина погружения центра масс автомата /; N - количество взаимодействующих с ним соседей автомата /; - рас-
стояние от центра масс автомата до поверхности контакта с автоматом у; - площадь поверхности
контакта в паре /-у; г=д+у - расстояние между центрами масс автоматов / и у; ту1 - единичный вектор, направленный от центра масс автомата у к центру масс автомата /. Вычисленное значение силы 6ша1 используется как дополнительный вклад в правой части уравнения трансляционного движения подвижного клеточного автомата / [9,10].
Рис. 3. Различные конфигурации частично погруженного автомата. Горизонтальной линией отмечен уровень жидкости
В случае, когда одиночный (изолированный) подвижный клеточный автомат частично погружен в жидкость, выталкивающая (Архимедова) сила равна Р0=рж^Уп (Уп - объём погруженной части автомата). Тогда ввиду призматической шестиугольной формы автомата существуют две конфигурации (рис. 3), для которых Уп легко описать кусочнозаданной функцией Н, таблица.
Так как подвижный клеточный автомат может вращаться [9], в общем случае он будет располагаться под некоторым углом в к положению 1 (рис. 3), и его объём Уп будет являться суперпозицией «крайних» значений 1 и 2 при одинаковой величине Н: Уп=аУ1+рУ2. В первом приближении можно положить в=(1-а)=2в/п (0°<в<30°).
В случае неизолированного (поверхностного) автомата из выталкивающей силы РА°=рлУпп вычитаются вклады АР/, соответствующие поверхностям контакта с «соседями», аналогично тому, как показано в (1). При этом, если площадки контакта погружены не полностью, аРау рассчитывается по формуле (рис. 4):
~ Рж8к • 1 (— -) % + |X - X |>г
AFf
I R
где '> = VJ+i
(У, - У j)
- погружен-
ная часть «ребра контакта» автоматов.
Рис. 2. Направления действия сил гидростатического давления (указаны стрелками) на свободные поверхности автоматов ¡(х,у) и\(х,у)
Рис. 4. Определение погруженной части «ребра контакта» автоматов i и j. Горизонтальной линией отмечен уровень жидкости. Стрелки указывают направления действия гидростатических сил
2.3. Расчет силы лобового сопротивления жидкости
В качестве диссипативной силы, ответственной за торможение движущихся в жидкости тел, в расчетах использовалась сила лобового сопротивления. Как известно, при различных значениях числа
Таблица. Функциональные зависимости объема погруженной части автомата
Объём V для конфигурации 1 Объём V2 для конфигурации 2
V =< 2WS •(2 W 2R „ „2 —¡^ + 4R л/3 10 R + R i л/з 1 R Л2 , ( 2R R = + h , h = —¡=, —т= 3 ) V л/3) ( R А , Г R R Л 3 + h ’ h = /Г’ /г 1 V V3 ) L V3 л/3) h-R 1(6R-2л/э • h) h = |~R^1 , л/3 JV 7 La/3 л/3 J Г 2 R -^(3R + h)(R + h)’ h = (-R’0] V2 =< ^3 2tf-J3 + 2Rh ^(4 R h )’ h = (0, R) L v3
Рейнольдса (Яе) используется линейная, либо квадратичная зависимость силы лобового сопротивления: Ра/=-кш о или Ра=-кщг и о, где 5 - вектор скорости тела (в данном случае - подвижного клеточного автомата). Первое приближение справедливо при Яе<40 (для рассматриваемой задачи это соответствует и и<10 3 м/с), второе - при 104>Яе>40 [12]. При моделировании деформационных процессов в хрупком твердом теле, сопровождающихся возникновением повреждений и трещин, значения скоростей автоматов могут изменяться в широких пределах. Поэтому в зависимости от локального значения числа Рейнольдса для рассматриваемого автомата должно использоваться то или другое приближение.
Общая формула для расчета «квадратичной» силы сопротивления жидкости движению погруженного тела имеет вид Рй=-0,5Схрж8 и и, где £ - площадь лобовой поверхности тела в направлении вектора скорости движения и, Сх - коэффициент, учитывающий форму тела. В зависимости от формы объекта Сх изменяется от 0,45 до 1,45. В первом приближении для клеточного автомата величина коэффициента Сх может быть положена постоянной и равной 1,528 (значение для цилиндра).
В случае изолированного полностью погруженного автомата площадь лобовой поверхности варьируется между двумя «крайними» значениями, соответствующими конфигурациям 1 и 2, описанным в пункте 2.2 (рис. 3):
^ = Яф + Б2(1 -в), (2)
где соответствует положению 1, а £2=2£0 -
положению 2 (£0 - площадь поверхности контакта двух недеформированных автоматов).
В рамках используемого приближения Сх=1,528 при расчете силы лобового сопротивления движению частично погруженного неизолированного (поверхностного) автомата его призматическая форма заменяется цилиндром с радиусом основания Ле#=£/(4Л). При этом площадь поперечного сечения полностью погруженного цилиндра (сечения в направлении, перпендикулярном к плоскости расположения автоматов) совпадает с 8йщ, вычисляемой согласно (2). Для расчета площади лобового сопротивления частично погруженного автомата осуществляется следующая последовательность операций:
• строятся две прямые, перпендикулярные вектору скорости (рис. 5); одна из них касающется основания цилиндра в точке, к которой направлен вектор скорости (линия ЬоШ на рис. 5), а другая проходит через центр масс автомата (линия Ьк);
• находятся координаты точек пересечения основания цилиндра (круга радиуса Ле#) с прямой (точки 1 и 2 на рис. 5), а также с горизонтальной линией поверхности жидкости (точки 5 и 6 на рис. 5).
• вычисляются длины перпендикуляров из точек
1, 2, 5 и 6 к линии Ьш (длины соответствующих отрезков обозначены на рис. 5 как гх, г2, ¿2);
• находятся минимальные отрезки в парах (ад) и (ад), которые и определяют длину 1йщ отрезка линии ЬоШ, соответствующего ребру лобовой поверхности автомата в плоскости (Х,У); на рис. 5 этот отрезок ограничен точками 7 и 4.
• площадь лобовой поверхности определяется
как 1йщ 2 Я.
В случае, когда какой-либо взаимодействующий сосед у автомата I закрывает часть его лобовой поверхности, площадь «закрытой» части вычитается из £га£:
N
^ ^ Ш 2К - £ 2Я • 1 сои(щ ).
j=1
Здесь ^ и 1?МК - площадь и длина ребра лобовой поверхности клеточного автомата I; индекс у нумерует взаимодействующих соседей автомата , находящихся перед его лобовой частью; - число таких соседей; определение угла <ру приведено на рис. 6.
Уровень жидкости ,7 \6
Рис. 5. Определение лобовой поверхности частично погруженного в жидкость цилиндра радиусом Reff. Горизонтальной линией отмечен уровень жидкости
i
Рис. 6. Определение угла щмежду вектором скорости автомата / и линией, соединяющей центры масс в паре \-\
Для малых чисел Рейнольдса (Яе<102) используем формулу, описывающую силу лобового сопротивления для шара радиуса г (закон Стокса):
Р^=-к1пи=-6пц/. и, где ¡лж - вязкость жидкости. Поскольку подвижный клеточный автомат имеет другую форму, при нахождении км используем следующее приближение: кип=Лцжк4Г. Из равенства «линейной» и «квадратичной» сил сопротивления при Яе=40 коэффициент Л для клеточного автомата г определяется как А,=15,28 £?тгЩ/ и в общем случае зависит от времени.
Таким образом, в рамках «цилиндрического» приближения выражение для силы лобового сопротивления движению в жидкости автомата г со скоростью и! имеет вид:
Р1 а/ = -0,5Схц. ц., при Яе > 40 Р. = -15,28цж^ и,. 1^., при Яе < 40 '
Вычисленное значение Р/, как и Ра°ш, используется в качестве дополнительного вклада в правой части уравнения трансляционного движения подвижного автомата г.
3. Тестовые расчеты и апробирование граничных условий
Корректность предложенного алгоритма расчета граничных условий для методов частиц «кундал-ловского типа», реализованного в рамках метода подвижных клеточных автоматов, проверялась в ряде тестовых численных расчетов. В частности, рассматривались следующие задачи:
• всплытие полностью или частично погруженных в воду тел (моделируемых ансамблем клеточных автоматов) симметричной и несимметричной формы с различной плотностью, ориентированных различным образом;
• торможение полностью или частично погруженных в воду тел, обладающих начальным импульсом.
Анализ результатов моделирования показал их качественное и количественное согласие с аналитическими оценками, а также результатами наблюдений и натурных экспериментов, проводимых на блоках льда озера Байкал. В частности, в случае ледовых блоков толщины 0,75 м равновесная глубина погружения в воду составляла около 90 % от их толщины, что соответствует аналитическим расчетам, результатам наблюдений и известным литературным данным [13]. Погруженные тела, ось симметрии которых была ориентирована перпендикулярно линии уровня жидкости, всплывали вертикально; в других случаях, по мере всплытия, они начинали поворачиваться как целые вплоть до переворачивания в положение, соответствующее наиболее низкому положению их центра тяжести. Характеристики диссипации кинетической энергии тел в воде также соответствовали аналитическим оценкам. Полученные результаты свидетельствуют о корректности предложенных граничных условий и возможности их использования для различных задач, связанных с моделированием методами частиц поведения твердых тел в жидкости.
Отметим, что развитые граничные условия были использованы авторами при численном моделировании методом подвижных клеточных автоматов деформационных процессов на конвергентных границах в модельной плитной среде (а именно, ледовом покрове крупного водоема) и образования деформационных структур субдукционного типа [14]. Значительный интерес к экспериментальному и теоретическому изучению фрагментов ледового покрова Байкала, проявляемый в последние годы, связан с тем, что с точки зрения реологии, блочного строения и динамики поведения он представляет собой среду, качественно подобную литосфере [6], хотя и более простую. Это дает возможность апробировать на ледовом покрове различные модели, впоследствии применяемые для теоретического изучения крупных фрагментов земной коры.
Моделирование с применением описанных граничных условий позволило получить некоторые важные результаты, касающиеся возможных механизмов образования и режимов функционирования деформационных структур субдукционного типа [14]. В частности, показано, что в результате образования таких структур величина силы сопротивления сжатию снижается на порядки. Это оказывает значительное влияние на кинематику плитной системы в целом. Результаты моделирования и данные экспедиционных исследований 2006-2008 гг. на ледовом покрове озера Байкал [6, 14] подтвердили возможность реализации механизма инициирования субдукции на «подготовленных» конвергентных границах в литосфере, связанного с процессами фрагментации земной коры.
Заключение
В работе предложна достаточно простая и эффективная реализация алгоритма расчета граничных условий на поверхности контакта твердого тела, моделируемого методом частиц, с жидкой средой. Предложенные граничные условия, учитывающие влияние гидростатического давления жидкости и диссипативную силу сопротивления движению твердого тела, реализованы в рамках одного из методов частиц «кундалловского типа», а именно метода подвижных клеточных автоматов. Результаты тестовых расчетов подтвердили корректность принятых граничных условий, что позволило использовать их для компьютерного моделирования деформационных процессов в одном из представителей класса специфических квазидвумерных блочных геологических сред, покоящихся на жидкоподобном основании (в ледовом покрове водоемов).
Изложенный алгоритм должен рассматриваться только как первое (базовое) приближение, область применения которого ограничивается изучением относительно медленных деформационных процессов и малых деформаций твердого тела. Его дальнейшее развитие может быть связано, в частности, с реализацией в рамках других математических моделей упруго-пластического деформи-
рования твердого тела, а также с учетом динамических процессов в жидкой среде (распространением в жидкости колебаний, что можно осуществить, например, по аналогии с [2], и наличием потоков). Предложенный в работе подход к реализации граничных условий в рамках методов частиц может быть модифицирован для случая контакта твердого тела со средой, обладающей низким, однако от-
СПИСОК ЛИТЕРАТУРЫ
1. Уилкинс М.Л. Расчет упруго-пластических течений // Вычислительные методы в гидродинамике / Под ред. Б. Олдера, С. Фернбаха, М. Ротенберга. - М.: Мир, 1967. - С. 212-263.
2. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 271 с.
3. Monaghan J.J. Smoothed particle hydrodynamics // Reports on progress in physics. - 2005. - V. 68. - № 8. - P. 1703-1759.
4. Cundall PA., Strack O.D.L. A discrete numerical model for granular assemblies // Geotechnique. - 1979. - V. 29. - № 1. - P. 47-65.
5. Добрецов Н.Л., Кирдяшкин А.Г., Кирдяшкин А.А. Глубинная геодинамика. - Новосибирск: Изд-во СО РАН, 2001. - 409 с.
6. Добрецов Н.Л., Псахье С.Г., Ружич В.В. и др. Ледовый покров озера Байкал как модельная среда для изучения тектонических процессов в земной коре // Доклады РАН. - 2007. - Т. 412. -№ 5. - С. 656-660.
7. Аксенов А., Коньшин В. Анализ задач взаимодействия «жидкость - конструкция» с использованием программных комплексов ABAQUS и FlowVision // САПР и графика. - 2006. -№ 9. - С. 20-24.
8. Liu M.B., Liu G.R. Restoring particle consistency in smoothed particle hydrodynamics // Applied numerical mathematics. - V. 56. -№ 1. - P. 19-36.
9. Псахье С.Г., Дмитриев А.И., Шилько Е.В. и др. Метод подвижных клеточных автоматов как новое направление дискретной
личным от нуля модулем сдвига. Такие условия являются более корректными, например, при моделировании деформационных процессов в литосфере.
Работа выполнена в рамках интеграционных проектов СО РАН № 27и 87, проекта № 3 программы Президиума РАН № 16, а также при поддержке грантов РФФИ (№ 06-05-64792-а) и Фонда содействия отечественной науке.
вычислительной механики. I. Теоретическое описание // Физическая мезомеханика. - 2000. - Т. 3. - № 2. - С. 5-15.
10. Синтез и свойства нанокристаллических и субструктурных материалов / Под ред. А.Д. Коротаева. - Томск: Изд-во Том. унта, 2007. - 368 с.
11. Работнов Ю.И. Механика деформируемого твердого тела. -М.: Наука, 1988. - 712 с.
12. Стрелец М.Х., Травин А.К., Шур М.Л., Спаларт Ф.Р Метод моделирования отсоединенных вихрей для расчета отрывных турбулентных течений: предпосылки, основная идея и примеры применения // Научно-технические ведомости. - 2004. -№ 2. - С. 22-33.
13. Сокольников В.М. Вертикальные и горизонтальные смещения и деформации сплошного ледяного покрова Байкала // Труды Байкальской лимнологической станции АН СССР. - 1960. -Т. 18. - С. 291-350.
14. Псахье С.Г, Шилько Е.В., Астафуров С.В. и др. Модельные исследования процессов возникновения и развития деформационных структур субдукционного типа в ледовом покрове озера Байкал // Физическая мезомеханика. - 2008. - Т 11. -№1. - С. 55-65.
Поступила 28.04.2008г.
