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

Вычислительные технологии

Том 14, № 6, 2009

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

О. В. Садовская, В.М. Слдовский Учреждение Российской академии наук Институт вычислительного моделирования СО РАН, Красноярск, Россия

e-mail: [email protected]

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

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

Введение

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

*Работа выполнена при поддержке РФФИ (грант № 08-01-00148), Междисциплинарного интеграционного проекта СО РАН № 40 и Комплексной программы фундаментальных исследований Президиума РАН № 2 "Интеллектуальные информационные технологии, математическое моделирование, системный анализ и автоматизация".

© ИВТ СО РАН, 2009.

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

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

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

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

1. Основные уравнения моделей

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

pV = V•о + рд, 2 а : о = Vv + Vv*. (1)

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

Принципиально важно то, что система (1) относится к симметрическим ¿-гиперболическим системам [6]. В декартовых координатах она записывается в матричной форме

„ди АЖ

вк

ди

джк

+ С,

(2)

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

и = (^1,^2,^3,^11,^22,^33,^23,^13,^12), матрицы А и Вк имеют блочный вид

0 о ) 0 вк

0 А, Акт 1 , вк = I В, 0

0 Ат/ \вк 0

А

Здесь А^ = р5 (5 — единичная матрица размерности 3 х 3),

0 0

а1111 а 1122 а1133 Ак = I а2211 а2222 а2233 ,а3311 а3322 а3333 /

а2323 а2313 а2312 Ат = 2 I а1323 а1313 а1312 , а1223 а1213 а12121

А,

1123 а1113 а 1112 \ а2223 а2213 а2212 I а3323 а3313 а3312)

5 523 0 0 0 5 512 513

0 5 513 0 ) , вт = 1512 0 523

0 0 512 513 5 523 0

5,

ч

вк

1, если к = г и к = 0, если к = г или к =

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

д дЬ0(и) _ д дЬ(и)

+ С(и).

д* ди джк ди

Система (3) обладает дополнительным законом сохранения

(3)

1(и дЬ0(и) дП ди

- ь0(и)

А

дж

и ■ ^^^ - Ь(и)) + и ■ С(и).

Для нее при условии строгой выпуклости Ь0(и) корректно поставлены задача Коши и краевая задача с диссипативными граничными условиями [7]. Эта система описывает обобщенные решения с сильными разрывами скоростей и напряжений (ударными волнами) и, кроме того, к ее решению могут быть применены эффективные численные методы сквозного счета, адаптированные к расчету разрывов [8, 9].

Для модели линейной теории упругости производящие потенциалы — квадратичные формы: 2 Ь0 = и ■ Аи и 2 Ьк = и ■ Вк и, первая из которых положительна. В более

к

к

сложной модели упругопластической среды, учитывающей необратимые деформации, тензор напряжений подчиняется ограничению о £ Г (Г — выпуклое множество допустимых напряжений), а тензор скоростей деформации е = (Vv + Vv*)/2 разлагается в сумму упругой составляющей, вычисляемой через скорость изменения тензора напряжений по закону Гука, и пластической составляющей ер = е — а : о. Для пластической составляющей выполняется неравенство Мизеса

(О — о): ер < 0, 0, о £ Г. (4)

В канонической форме система приводится к вариационному неравенству

ф — и).(°д-Ш — — а(иЛ > 0, и, и £ г. (5)

1 ; Уд^ ди дхк ди 1 ') " ' v у

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

Для описания волновых процессов в грунтах применяется математическая модель идеально сыпучей среды, которая в зависимости от уровня напряжений проявляет упругие или пластические свойства при сжатии и не сопротивляется растяжению [11]. В рамках этой модели тензор малых деформаций разлагается в сумму упруго-сыпучей и пластической составляющих: е = ес + ер. Определяющие соотношения для первого слагаемого приводятся к неравенству Хаара и Кармана:

(о — о): (а : о — ес) > 0, о, о £ К, (6)

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

Можно показать, что единственным решением неравенства (6) является тензор вп, где в — тензор условных напряжений, вычисляемый по линейному закону Гука а : в = ес, а индекс п означает проекцию на конус К по евклидовой норме |о| = у/о : а : о. Зависимость о от ес не является взаимно-однозначной, что свидетельствует о механической некорректности модели — по заданным напряжениям невозможно, вообще говоря, определить деформированное состояние среды. Регуляризацией (6) служит уравнение

о = Лв + (1 — Л) вп, 0 < Л < 1,

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

— — — аоо) > 0, ^ £ Г, (7)

V = Ли +(1 — Л) ип, и = ! V — Л . (8)

ЛЛ

Здесь V и и — вектор-функции, первая из которых составлена из отличных от нуля компонент вектора скорости частиц среды v и тензора действительных напряжений о, а во вторую вместо о входит тензор условных напряжений в.

Для описания упругого деформирования материалов с микроструктурой служит математическая модель моментного континуума Коссера [12, 13]. Основное отличие этой модели от классической линейной теории упругости состоит в том, что наряду с поступательными степенями свободы в ней учитывается вращение частиц. В рамках геометрически линейного приближения уравнения поступательного и вращательного движения моментной среды

р г; = V •о + рд, ] и = V •т — 2 оа + (9)

записываются относительно векторов линейной скорости v и угловой скорости и, несимметричных тензоров напряжений о и моментных напряжений т. Здесь ] — суммарный момент инерции частиц в единице объема, д — вектор массовых моментов, оа — вектор антисимметричной составляющей тензора напряжений (о — о*)/2. Как обычно, в декартовой координатной системе антисимметричный тензор

10 — из и2

и= 1 из 0 —и1

— и2 и1 0

отождествляется с вектором и, координатами которого являются числа и>1, и2, из.

В модели Коссера для изотропной моментной упругой среды постулируются тензор-но-линейные определяющие уравнения

о = Л (6 :Л*) 6 + 2 ^Л' + 2 а Ла, т = в (6 : М*) 6 + 2£М* + 2 пМа (10)

с шестью константами — параметрами материала: Л, а, в, С и П

Полная система уравнений моментной теории упругости (9), (10) также приводится к симметрической системе

Аж = вк ^ + «и + а

д£ дхк

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

и = (Vl, V2, Vз, оИ, о22, озз, о23, оз2, оз1, о1з, о12, о21,

иь и2, из, тп, т22, тзз, т2з, тз2, тз1, т1з, т12, т21).

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

- К)(Л» — ВкЦ! - «К — ^ > 0, (Л, К £ Г, (11)

где вектор-функции V и и связаны между собой уравнениями (8).

2. Вычислительная технология

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

Рассмотрим вычислительную технологию на примере решения вариационного неравенства наиболее общего вида (11), являющегося частным случаем неравенства (7). Явный по времени алгоритм численной реализации неравенства строится с помощью метода расщепления по физическим процессам следующим образом: сначала на каждом временном слое решается задача деформирования разномодульной упругой среды, а затем полученное решение корректируется для учета пластических свойств. В начальный момент времени задаются некоторые начальные данные и(0,ж) = и0(ж). Граничные условия формулируются в скоростях или в напряжениях.

Пусть и — вектор-функция, которая определяется через искомое решение краевой задачи из системы

= в'^ + «V + С

д* джк

и, кроме того, в фиксированный момент времени * подчиняется условию [/(¿,ж) = и(*,ж). Тогда из исходного неравенства после аппроксимации производной по времени конечной разностью на интервале + Д*) в каждой ячейке пространственной сетки возникает задача

(V" - V) А (и - /) > 0, V, V €

которая с учетом (8) приводится к виду

(V" - V) А (V - (1 - Л) V* - Л /) > 0.

Отсюда, по определению проекции,

V = ((1 - Л) V* + Л/)П (12)

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

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

правлениях х2 и х1. Формулы расщепления выглядят так:

ди1 д(1

А—^ = В1 д— + а1, и1 (/, х) = и (/,х), д/ дх1

ди2 д(2

^ = В2 д— + а2, и2(/,х) = и1 (/ + Д//2,х), д/ дх2 диз д(

а-^ = Взд— + аз, из(/,х) = и2(/ + Д//2,х),

д/ дх.

з

ди 4

А^- = «V4, и4(/,х) = из (/ + Д//2, х), ди/5 дV5

^ = Вз д— + аз, и5(/ + Д//2, х) = и4(/ + Д/,х),

ди6 дV6

^ = В2 д— + а2, и 6(/ + Д//2, х) = и 5(/ + Д/,х),

/ х2

ди7 дУ7 А-^ = В1 д— + а1, и7(/ + Д//2, х) = и6(/ + Д/,х).

Искомое значение и(/ + Д/, х) равно и7(/ + Д/,х). Векторы ак выбираются исходя из принципа физической адекватности одномерных моделей, а1 + а2 + аз = а. Например, при учете силы тяжести, действующей в направлении оси х1, вектор массовых сил а целиком переносится на первый и последний этапы. При расчете двумерной (плоской или осесимметричной) задачи в методе расщепления отсутствуют третий и пятый этапы, относящиеся к направлению хз.

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

„ ит — и „ ит + и

А-гт-= «

Д/ 2

(г — номер шага по времени). В общем случае для нелинейной системы применяется следующая схема:

и^+1 — и® и^+1 + и® / \

А и Д/ и = « Г+1/2, Г+1/2 = Л и 2+ и + (1 — Л) (а (ип)т + (1 — а)(ип)*),

где а — параметр, который в каждой ячейке сеточной области определяется из соображений консервативности как решение нелинейного уравнения

(а(и- Г« + (1 — а)(и- )')Аи'+' — и' = ^ + ' >' А (ЦП^+Д—^.

Д/ 2 Д/

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

^А — иг+1 = ^А + ^ + (1 — Л) Д/« (а (ип)г+1 + (1 — а)(ип)*),

решение которой строится по методу последовательных приближений [4].

Оставшиеся шесть одномерных систем уравнений на этапах расщепления решаются с помощью явной монотонной ЕМО-схемы типа "предиктор-корректор" с кусочно-линейными распределениями скоростей и напряжений по ячейкам [9]. В случае постоянных матриц-коэффициентов на шаге "корректор" используются соотношения

и+1/2 = иг+1/2 + Д А-1 (в' + Ск

2 V Дхк

Здесь индекс г + 1/2 относится к центру ячейки пространственной разностной сетки, верхний индекс соответствует текущему временному слою, нижний — предыдущему слою. Вектор Vг+1/2 вычисляется через иг+1/2 по формуле (8). Если матрицы переменные, то в качестве разностной производной по X' берутся соответствующие слагаемые консервативной аппроксимации. Для замыкания схемы необходимо доопределить ее соотношениями шага "предиктор". Пусть У, и с, — полная система левых собственных векторов и собственных чисел матрицы в'А-1: У, в' = с, У, А, У, А У, = Умножая уравнения системы слева на вектор У,, перейдем к системе дифференциальных уравнений, которые для модели упругой среды представляют собой уравнения на характеристиках:

УА^ = с, У + У С' •

д* дхк

После аппроксимации получим

(1?+1/2)± = 4*1/2 ± а,т/2 Д2Хк + (с, в, + У Ск)г+1/2

где а, ¿+1/2 и вг»+1/2 — производные от коэффициентов разложения и и V по базису У,: /, = (У, А)^+1/2 и и = (У, А)^+1/2 V, полученные с помощью итерационной процедуры предельной реконструкции. Индексами " -" и "+" отмечены значения этих коэффициентов на левой и правой границах одной и той же ячейки. Процедура предельной реконструкции позволяет повысить точность численного решения и состоит в построении монотонных кусочно-линейных сплайнов, приближающих /, и с минимальными разрывами на границах соседних ячеек сетки. Во внутренних узлах расчетной области на шаге "предиктор" величины и находятся по формуле осреднения и = (и+ + и—)/2 через значения и , относящиеся к разным сторонам границы между ячейками и удовлетворяющие системе нелинейных алгебраических уравнений:

(У,А)г+1/2 и+ = 4—+1/2 для С > 0,

(У,А)г-1/2 и— = /,+—1/2 для С, < 0,

А V+ = .

В этой системе уравнения с матрицами А, представляют собой условия непрерывности вектора скорости и вектора напряжений при переходе через границу, а число уравнений с учетом данных условий равно числу неизвестных величин. Условия непрерывности реализуются численно методом последовательных приближений. Значения V пересчи-тываются через и в соответствии с (8).

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

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

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

Каждый из узлов кластера при выполнении программы-препроцессора упаковывает исходные данные в двоичный файл прямого доступа — файл вещественных чисел, в который поблочно записываются параметры материала, часть сетки, приходящаяся на этот узел, и начальные значения решения, и файл целых чисел, содержащий соответствующие им адреса (указатели) — порядковые номера первых элементов. Вещественные файлы такой же структуры в дальнейшем создаются основной программой для организации контрольных точек при счете и для последующего анализа полученных результатов. Суммарный размер этих файлов может значительно превышать объем оперативной памяти отдельного процессора. Каждый процессор при старте основной программы считывает целочисленный файл и соответствующий ему вещественный файл. Далее целочисленный массив редуцируется — параметры сетки и указатели принимают в нем индивидуальные значения для данного узла. Балансировка вычислительной нагрузки достигается за счет равномерного распределения сеточной области между узлами кластера. Если размерность сетки какого-либо блока больше средней размерности в расчете на один узел кластера, то этот блок обслуживается несколькими узлами, и наоборот, один и тот же узел обслуживает несколько блоков, если их суммарная размерность не превышает средней. Из соображений минимизации количества пересылок используются 1D-, 2D- или 3Д-разбиения расчетной области.

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

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

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

3. Алгоритм сжатия файлов

Идея сжатия файлов числовых данных для компактного хранения и передачи решения заимствована в методе разделения переменных, в соответствии с которым заданная скалярная функция и, зависящая от вектора х = (х^ х2,..., хп), разлагается в бесконечный ряд. Элементы этого ряда представляют собой произведения одномерных функций:

и(х) = 5] 71X 1'1(х1)Х1'2 (х2)...Хг'п(х„). 1=1

Дискретный аналог ряда с разделенными переменными выглядит так:

х

Хда.-ХГ. (13)

¿n

1=1

Базисные функции Xг'2, ..., X 1'п определяются из соображений наилучшего приближения и. При вычислении первого слагаемого ряда решается задача условной минимизации квадратичной функции [14]

М1 М2 мп

ЕЕ-ЕК * - тх! Х2 ...X")2 (14)

¿1 = 1 ¿2 = 1 ¿п = 1

относительно переменных 7, X,¿11, X,2 , ..., X™. Предполагается, что выполнены следующие условия нормировки:

М1 М2 Мп

Е№ )2 = E(X2 )2 =... = E(х: )2 = 1.

Точка минимума (14) удовлетворяет системе уравнений, которая может быть получена с учетом условий нормировки методом множителей Лагранжа:

Ml M2 Mn

-1 V"2 vn

¿n '

¿1 = 1 ¿2 = 1 ¿n = 1

Y = E E - E Uii¿2-inX1X2... Xin;, (15)

"¿1 ¿2 ...¿n

Mi

Mí-1 M¡+1

X,

E-E

11 = 1 гг—1 = 1 »¡+1 = 1

Mn

^ ] Í2 ••• »n Xi1

in = 1

X i—1 X i+1 X™

" гг—1 гг+1"' »n

\

M¡ / M1

ЕЕ-

гг = 1 у г1 = 1

(16)

M¡—1 M¡+1

E E

»г—1=1 »¡+1=1

Mn

E

ь- ■ ■ X1

ь»1 »2 ••• tn Xt1 '

TAÍ— 1 тЛ + 1 ' Xí¡—1 '

»n

= 1 »n = 1

Эта система используется при численной реализации алгоритма. Параметр 7 определяется по формуле (15) в завершение расчета, когда базисные функции уже найдены. Для нахождения базисных функций применяется один из вариантов метода последовательных приближений. Учитывая условие нормировки, начальное приближение для Х^ принимается равным 1 /\[М\. Все последующие приближения рассчитываются по формуле (16), в которой значения Х?+, ..., X™ берутся с предыдущего шага, а в качестве X,*, ..., X1— участвуют вновь вычисленные значения. Таким образом, реализуется полунеявная процедура вычислений. Второе и все оставшиеся слагаемые ряда (13) находятся аналогично. При определении второго слагаемого дискретная функция ,2 заменяется погрешностью ,2 — X,2 ...X™ , полученной после вычисления на предыдущем этапе, и т.д. Условие окончания счета формулируется в терминах относительной среднеквадратичной погрешности приближения.

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

b(x1,x2,x3) = e h(r ro) /r2, r

2 , 2 1 2

X1 + X2 + X3 ,

описывающей локализацию решения на сфере радиуса г0. Безразмерный параметр локализации Н = 20, радиус сферы равен наибольшей диагонали параллелепипеда. В данном случае были вычислены 12 членов ряда. Относительная погрешность, характеризующая потЁрянзнеуденннфорэма^Ерасчравсна 2&азалж,эффищшрщитный фрйцас^яулйется экономичным в том смысле, что требуемое для заданной точности число итераций метода последовательных приближений практически не зависит от размерности дискретной задачи.

2

2

Рис. 1. Локализация решения на сфере и результат сжатия

4. Результаты расчетов на кластерах

Разработанный комплекс программ применялся к решению методических задач, ориентированных на приложения в геофизике (сейсмике), с целью демонстрации возможностей распределенных вычислений в этих задачах. Численные расчеты проводились на кластерах МВС-1000М Института вычислительного моделирования СО РАН (г. Красноярск) и МВС-100К Межведомственного суперкомпьютерного центра РАН (г. Москва).

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

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

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

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

Сосредоточенная нормальная нагрузка стц = — р*1 $(ж) ^(¿) действует мгновенно в одной ячейке верхней границы расчетной области. На фронтальной границе области решения ставятся условия симметрии, на остальных частях боковой поверхности и нижней границе — неотражающие условия, верхняя граница представляет собой свободную поверхность. Область распределяется между 68 процессорами по принципу равномерной загрузки, каждым из процессоров решается часть задачи на подсетке размерностью 50 х 50 х 50 ячеек. На рис. 3, а изображены поверхности уровня напряжения ац, первый из них относится к моменту, когда отраженная волна еще не достигла свободной поверхности, второй — к моменту, когда преломленная волна обгоняет падающую продольную волну. На рис. 3, б приведена сейсмограмма, описывающая поведение вертикальной компоненты перемещения и в зависимости от времени. Приемники расположены на верхней границе расчетной области вдоль линии, проходящей через точку приложения нагрузки параллельно горизонтальной оси. Видны падающие продольные и поперечные волны, отраженные и преломленные волны. Годографы падающих продольных и преломленных волн представляют собой прямолинейные отрезки, годографы отраженных волн — кусочно-параболические сплайны, меняющие свою кривизну в точках перехода от отражения от поверхности раздела к отражению от угла.

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

Рис. 3. Задача Лэмба для среды с жестким включением: а — поверхности уровня нормального напряжения в различные моменты времени, б — сейсмограмма отраженных волн

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

На рис. 5 приведены результаты расчетов для задачи о периодическом воздействии сосредоточенной нормальной нагрузки ац = — p* #(x)sin(2 nvt) с частотой v. Изображены поверхности уровня напряжения при частоте внешнего воздействия, равной частоте собственных колебаний вращательного движения частиц v* = 1/T, где T = n^j/a — период колебаний частиц среды, а также при v = 0.5 v* и v = 1.5 v*. В отличие от предыдущих задач здесь нижняя граница расчетной области закреплена, поэтому от нее внутрь области распространяются отраженные волны. Расчеты проводились для пенистого полиуретана на 64 процессорах, 200 шагов по времени были посчитаны примерно за 10 ч.

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

Рис. 4. Задача Лэмба для моментной упругой среды: схемы нагружения и сейсмограммы падающих волн

Рис. 5. Задача о периодическом воздействии сосредоточенной нагрузки для моментной среды слева направо: v = 0.5 v*, V = V* и v = 1.5 v*

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

Список литературы

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

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

[3] Lakes R. Experimental methods for study of Cosserat elastic solids and other generalized elastic continua // Continuum Models for Materials with Micro-Structure / Eds. H. Muhlhaus, J. Wiley. New-York, 1995. P. 1-22.

[4] Садовская О.В., САдовский В.М. Математическое моделирование в задачах механики сыпучих сред. М.: Физматлит, 2008. 368 с.

[5] РАвотнов Ю.Н. Механика деформируемого твердого тела. М.: Наука, 1979. 744 с.

[6] Годунов С.К., РомЕнский Е.И. Элементы механики сплошных сред и законы сохранения. Новосибирск: Научная книга, 1998. 280 с.

[7] Годунов С.К. Уравнения математической физики. М.: Наука, 1979. 391 с.

[8] Численное решение многомерных задач газовой динамики / Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. М.: Наука, 1976. 400 с.

[9] Куликовский А.Г., ПогоРЕлов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 607 с.

[10] САдовский В.М. Разрывные решения в задачах динамики упругопластических сред. М.: Физматлит, 1997. 208 с.

[11] Садовская О.В. Метод сквозного счета для исследования упругопластических волн в сыпучей среде // Журн. вычисл. математики и мат. физики. 2004. Т. 44, № 10. С. 1909-1920.

[12] Cosserat E., Cosserat F. Theorie des corps deformables // Chwolson's Traite Physique: 2nd ed. Paris, 1909. P. 953-1173.

[13] ПАЛЬмов В.А. Основные уравнения теории несимметричной упругости // Прикл. математика и механика. 1964. Т. 28, вып. 3. С. 401-408.

[14] Варыгина М.П., Киреев И.В., Садовская О.В., САдовский В.М. Программное обеспечение для анализа волновых движений в моментных средах на многопроцессорных вычислительных системах // Вестник Сибирского государственного аэрокосмического унта им. акад. М.Ф. Решетнева. 2009. Вып. 2 (23). С. 104-108.

[15] САдовский В.М. К теории распространения упругопластических волн в сыпучих средах // Докл. РАН. 2002. Т. 386, № 4. С. 487-489.

Поступила в редакцию 20 октября 2009 г.