автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование деформационных динамических процессов в средах со сложной структурой

На правах рукописи

Челноков Федор Борисович

Численное моделирование

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

Специальность 05.13.18 Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Москва - 2005

Работа выполнена на кафедре информатики Московского физико-технического института (государственного университета)

Научный руководитель:

доктор физико-математических наук, профессор Петров Игорь Борисович

Официальные оппоненты:

доктор физико-математических наук, профессор Кондауров Владимир Игнатьевич,

Ведущая организация:

Институт математического моделирования РАН

заседании диссертационного совета К 212.156.02 при Московском физико-техническом институте (государственном университете) по адресу: 141700, Московская обл., г. Долгопрудный, Институтский пер., д. 9.

С диссертацией можно ознакомиться в библиотеке МФТИ.

кандидат физико-математических наук Антоненко Максим Николаевич

Защита состоится «

2005 г. в 3

часов на

Автореферат разослан с 20 » СЙкЯ _2005 г.

Ученый секретарь диссертационного совет кандидат физико-математических наук

С .рактеристика работы

Актуальность темы

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

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

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

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

Цель работы

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

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

Научная новизна

1. Реализован комплекс программ для исследования динамических волновых задач в неоднородных телах, в том числе, в многослойных, перфорированных, кавернозных и трещиноватых средах.

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

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

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

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

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

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

Практическая ценность

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

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

Защищаемые положения

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

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

3. Комплекс программ применен для исследования природы распространения волн в трехмерных перфорированных конструкциях.

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

5. Спектральное исследование матриц коэффициентов уравнений позволило добиться ускорения счета и упрощения программного кода.

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

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

Публикации

Научные результаты диссертации опубликованы в работах [1-15]. Апробация

Результаты работы были доложены, обсуждены и получили одобрение специалистов на следующих научных конференциях:

• международная конференция «Математика. Компьютер. Образование» (Пущино, 2005) [11];

• электронная конференция МЭИ «Топливо и энергетика» (Москва, 2004) [12];

• научные конференции МФТИ «Современные проблемы фундаментальных и прикладных наук» (Долгопрудный, 2002 - 2004) [13-15].

Структура и объем диссертации

Диссертация состоит из введения, 8-ми глав и заключения. Общий объем диссертации составляет 251 страницу. Список литературы содержит ссылки на 75 публикаций.

Содержание работы

Введение

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

Глава 1. Численное решение уравнений теории линейной упругости

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

рЬ = У-Т,

(1)

т = А(У + +

Здесь р и А, ц — константы упругого материала: плотность и параметры Л яме соответственно. Переменными являются V — скорость движения среды в данной точке, и Т — симметричный тензор напряжений Коши. V — оператор градиента, I — единичный тензор, (V • Т)* = £ ®

— оператор тензорного произведения: (а®6)и = а1Ь}. Две скорости звука выражаются через константы материала следующим образом:

/а + 2/х П1

В вычислительной математике при анализе и построении разностных схем для таких систем уравнений как (1) принято приводить ее запись к канонической форме

где А; — квадратные матрицы с постоянными коэффициентами, а вектор й включает все переменные системы (ЗБ):

т т

и = {у, Т} = Х1уь ЬХХ1¿т ^уу> ¿гг} )

— векторы выбранного базиса на плоскости, не обязательно ортогонального или нормированного; -Л- = • V — производные вдоль направлений,

задаваемых базисными векторами. Для системы (1) показывается, что

_ -у

действие А, на произвольный вектор {¿7, Т} имеет вид

А,

¿(Й.-Т)

А(п, • и)1 + /х(п, 0 V + г7 ® п,)

(3)

где {Ьпх, и^г) — биортогональный базис по отношению к , £2), а векторы п, имеют единичную длину. В последующем индексы г опускаются, и для каждого вектора п, в отдельности рассматривается пара векторов (ЗБ), образующих с ним совместно ортонормированный базис и также обозначаемых (для каждого г) п^о), П\,П1- Кроме того, вводятся симметричные матрицы Иу = N¿4 с равной пространству размерностью как

N>.7 = ^® Щ + ® (*>= 0,1,2)

Для каждой матрицы А существует разложение:

ПА = АП,

(4)

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

А = ¿сИа§{с1, -сьс2, -с2,с2, -с2,0,0,0},

(5)

т, и'2

■Ш7 №8 ■Шд

= п

п-(«Т ^Тп)

"1 ■ (" Т фТЙ) ^•(^¿Тп)

N12 : т (N11-N22)^ (I - ТгУ^оо) : Т

Столбцы матрицы Г2 1 имеют следующий вид:

(6)

а*'1-2 = -

ТР[(С1 - с3)^о + с31]

&'7 =

Г,3,4 =

о

4N12

щ

¿Г'8 = г

¿.,5,6

п2

Т2 С2^02

N11 - №

22

¿5*.® = ±

О

1-1Ч

оо

Полученное разложение используется для специализации, в частно-с т. схемы Куранта - Изаксона - Риса, в одну из форм записи которой входит произведение П-1|Л|Г&:

£

С2# + (сх - сг)(п • у)п = I (Т : N00)1(01 - 2с2 - с3)^о + с31]+ . +с2[(Т • п) ® п + п ® (Т • п)]

Для граничных узлов, на которых поставлено условие вида

В гГ+1 = 6

(где В — прямоугольная матрица, Ь — вектор правых частей), выводится формула их двухэтапного расчета:

'Здесь '«*" определятся в результате первого этапа из значений » местах пе-]>ессчения уходя1Цнх внутрь тела характеристик с предшествующим временным слоем; прямоугольная матрица О*'""1 составляется из столбцов !П~г, «ютвстствующих характеристикам выходящим наружу.

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

Глава 2. Построение нерегулярной треугольной сетки

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

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

= 1?" - П'<ои1 (ВГГ'ои*у1 {Ы?п - Ь).

(8)

сетки, согласованной с заданными контурами. Эффективность алгоритма достигается благодаря линейной сложности триангулирования монотонных многоугольников, на которые можно условно разделить любую исходную фигуру. По определению, монотонным является такой многоугольник, пересечение которого любой вертикальной прямой либо пусто, либо является непрерывным. Доказывается оценка трудоемкости 0(п log' ri), где п — суммарное количество вершин во всех входных контурах.

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

Предлагается новый алгоритм построения квазиравномерной треугольной сетки, который исходит из заданной триангуляции тела, внося в нее минимальные изменения. В качестве исходной берется сетка с предшествующего шага по времени, вершины которой смещены на произведение скорости тела в этих точках и величины шага интегрирования. Мелкость сетки и допустимая степень ее неравномерности оиределякпч я двумя числовыми параметрами алгоритма: Z*llu > 0 и а > 1 соответственно.

Минимальность вносимых изменений в сетку (введение и удаление вершил) обеспечивается за счет увеличения параметра cv (в расчетах использовалось а = 1.1). Идея заключается в том, что каждое ребро в триангуляции может иметь длину от некоторого минимального размера 1„,¡„ до размера большего в 2а раз. При этом допустимые пределы на размеры объектов сетки увеличиваются за счет се равномерности.

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

Шаги алгоритма поддержания заданной плотности сетки:

1. Построение триангуляции Делоне из данной изначально и поддержание этого свойства в дальнейшем.

2. Сокращение граничных вершин до тех пор, пока имеются граничные ребра короче lbmin = /(£,,„,а).

3. Введение дополнительных вершип в центрах граничных ребер длии-нее 'max —

4. Сокращение внутренних вершин, одно из ребер которых короче [¡ши.

5. Введение новых внутренних вершин в центрах описанных окружностей треугольников сетки, если их радиус больше R'nax = nl'niu. а центр попадает внутрь тела.

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

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

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

Рис. 1. Слева направо: начальная триангуляция неодносвязной области, триангуляция Делоне, квазиравпомерпая треугольная сетка.

Глава завершается несколькими примерами сеток, полученных в результате выполнения описанных алгоритмов (рис. 1).

Глава 3. Контакт тел в динамических задачах

В главе рассматриваются особенности численного моделирования в задачах с несколькими деформируемыми телами. За время моделирования тела могут неоднократно касаться и отрываться друг от друга. Также возможен контакт между разными частями одного и того же тела.

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

подход, предлагаемый и данной диссертации, основывается па построении триангуляции Делоне пространства между телами (рис. 2). Отмечается ряд аспектов, в силу которых имеет смысл предпочесть именно его:

• Возможность использования одних методов, как для построения сетки внутри тел, так и для определения контакта между ними.

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

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

Рис. 2. Удар по упругой решетчатой конструкции. Слева: взаимодействующие тела, справа: триангуляция пространства между заключенными в ограниченный объем телами, используемая для установления контактом.

В главе исследуется вопрос, как выбрать допустимый шаг интегрирования, чтобы за его период одно быстролетящее тело не проникло вглубь другого или не перескочило бы через него, не «заметив» факт контакта.

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

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

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

Глава 4. Интерполяция в треугольнике

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

Для того чтобы иметь возможность однозначно определять коэффициенты полинома степени N в треугольнике ABC проведем прямые, параллельные его сторонам, которые делят каждую из его сторон на N частей (рис. 3). В результате исходный треугольник подразделяется на JV2 меньших треугольников, равных друг другу и подобных всему треугольнику. KojHrrecTBO точек внутри треугольника, в которых прямые пересекаются между собой и со сторонами треугольника равно + 2)(N +1), что совпадает с количеством коэффициентов у искомого полинома. Именно эти точки пересечения выбираются в качестве опорных.

Рис. 3. Равномерное распределение опорных точек в треугольнике.

Обозначим координаты вершин треугольника символами Га, fg, re, а ссылаться на опорные точки будем при помощи ra¡,c. Индекс состоит из грех частей, каждая из которых описывает место пересечения прямой, проходящей через опорную точку и параллельной определенной стороне, с другой стороной треугольника (рис. 3).

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

111

SA = -[г<: -fB,f-fв]. SB = ^[гА ~гс, Г-re], Se - ^[гв-ГА,Г-ГЛ],

где S, — площадь треугольника, не содержащего вершины г ABC. Квадратными скобками обозначено векторное произведение, результат кото-

poro ортогонален плоскости и считается скаляром. Используем также относительные площади тех же треугольников:

Sa SB Se с Ir- , - -, , . 1

SA — Sb ~ Sc~~s'' = 2^Гв~Га, с~Га^' sA+sb+sC = 1-

Значение полинома v(r) необходимо выразить через значения vai>c, которые он принимает в опорных точках:

v{r) = wabc(r)vabc, а,Ь,с

где wabc(r) — вес опорной точки ГаЬс, также являющийся полиномом степени N. Поскольку в опорных точках полином должен обращаться в заг данные значения, то каждый вес равняется единице в одной опорной точке и нулю во всех остальных опорных точках. Следующая функция удовлетворяет поставленным условиям:

wMf) = Н^51;^ ~ ^,, №е{А,в,С},о<т<Щ ПГ=1(^(гаЬс) -

при условии правильного подбора величин {Г,} и {n¿}, описанного в диссертации. Приводятся примеры весов для полиномов первых четырех порядков, а также формулы градиента интерполяционных полиномов и их интегралы по треугольнику.

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

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

1. Определить пробную квадратичную функцию по вышеописанным формулам из известных значений в 6-ти опорных точках (рис. 4).

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

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

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

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

Глава 5. Численный метод для бесструктурных треугольных сеток

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

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

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

системы уравнений (2), для решения которой используется расщепление но направлениям: последовательно рассматриваются две системы вида й + = 0. Если сделать замену переменных v — flu. где О определяется по формуле (4), то оказывается, что необходимо решить ряд независимых уравнений переноса Vj + Xj ^ = 0. Эффективное вычисление п и о-1, необходимых для перехода к уравнениям переноса и обратно, реализуется благодаря результатам главы 1. Решение уравнений переноса может быть выписано в следующей форме:

Vj(fk, tn+1) = Vj(ffc - AjrCu *»),

где в правой части равенства стоит восстановленное значение с временного слоя tn в точке пересечения его характеристикой уравнения, проходящей также через (г*, £n+1), fk — положение узла, в котором тциггн ' решение. Если сетка была ростроена методом, предложенным в главе 2, то можно легко выбрать такой максимальный шаг т, при котором все подобные точки будут принадлежать смежным с даппым узлом треугольникам (рис. 5). Для граничных узлов не все уравнения переноса могут быть" разрешены, поэтому выясняется тип граничных условий, возможно, с использованием алгоритмов главы 3, после чего выполняется соответствующая корректировка. Восстановление решения в треугольнике сводится к применению одного из способов интерполяции, изложенных в главе 4. выбор интерполяции определяет порядок и свойства численного метода.

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

Главк 6. Распространение упругих волн в неоднородных массивных породах

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

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

Рис. 6. Подробный план сетки рядом с круговыми кавернами и трещинами.

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

Рис. 7. Вертикальная составляющая скорости точек среды. Слева: плоская волна, распространяющаяся в сторону кавернозной зоны. Справа: отраженные и рассеянные волны, полученные в результате математического моделирования.

На рис. 7 представлен результат моделирования прохождения сигнала Берлаге частотой 300 Гц сквозь скопление из 177 каверн диаметром 1 м.

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

Глава 7. Распространение упругих волн в перфорированных средах

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

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

а. Однородный материал. Ь. Перфорированная среда.

Рис. 8. Распространение волн в трехмерных телах кубической формы. Показано центральное сечение и изоповерхности модуля скорости.

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

Глава 8.

Высокоскоростной удар по многослойной преграде

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

Рис. 9. Слева: вихревые структуры в поле скоростей в верхнем слое преграды.

Справа: образованные трещины и зоны раздробленного стекла после удара.

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

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

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

Заключение

Заключение содержит в развернутом виде формулировку основных результатов и выводов диссертации.

Основные результаты и выводы диссертации

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

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

3. Численными методами исследовано распространение упругих волн в перфорированных средах, включая трехмерную постановку.

4. Проведено математическое моделирование зависимости характера отраженной волны от кластеризации неоднородностей в породе.

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

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

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

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

Список работ по теме диссертации

1. Челноков Ф. В., Мацисвский Н. С. Реконструкция кусочно-полино-мналыюго непрерывного поля па плоскости по значениям в опорных точках треугольной сетки // Процессы и методы обработки информации: Сб. ст. / Моск. физ.-тех. ин-т. - М., 2005. — С. 201 - 208.

2. Истрой И. В., Челноков Ф. В. Численное исследование волновых процессов и процессов разрушения в многослойных преградах // Ж. вы-числ. матем. и матем. физ. — 2003. — Т. 43, № 10. — С. 1562 - 1579.

3. Петров И. В.. Челноков Ф. В., Чнбриков В. В. Численное исследование волновых процессов в перфорированных деформируемых средах // Матем. моделирование. - 2003. — Т. 15, № 10. — С. 89 - 94.

4. Петров И. В., Челноков Ф. В. Численная проверка прочности железобетонной наружной оболочки иод действием динамической наг грузки // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ии-т. — М., 2003. -- С. 4 - 13.

5. Пет,ров И. В., Ртвелиашвили Д. П., Челноков Ф. Б. Численный расчет разрушения бетонных конструкций с учетом влияния гравитации // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ин-т. • - М., 2003. - С. 14 - ]8.

(». Агапов П. И., Челноков Ф. Б. Сравнительный анализ разностных схем для численного решения двумерных задач механики деформируемого твердого тела // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ин-т. М., 2003. — С. 19 -- 27.

7. Петров И. В., Челноков Ф. В., Чибриков В. В. Расчет волновых процессов и процессов разрушения в пористых средах // Обработка ин-

формации и моделирование: Сб. ст. / Моск. фи ¡.-тех. ин-т. — М., 2002. - С. 137 - 147.

8. Агапов П. И., Петров И. Б., Челноков Ф. Б. Численное исгледовахшс задач механики деформируемого твердого тела в неоднородных областях интегрирования // Обработка информации и моделирование: Сб. ст. / Моск. ф из .-тех. ин-т. - М., 2002. — С. 148 - 157.

9. Агапов П. И., Обухов А. С., Петров И. Б., Челноков Ф. Б. Численное решение динамических задач биомеханики агючио-характеристическим методом // Компьютерные модели и прогресс медицины: Сб. ст. / РАН. - М.: Наука, 2001. - С. 275 - 300.

10. Агапов П. И., Обухов А. С., Петров И. Б., Челноков Ф. Б. Компьютерное моделирование биомехшгачееких процессов еоточпо-характеристическим методом // Управление и обработка информации: модели процессов: Сб. ст. / Моск. физ.-тсх. ин-т. - М..'2001. С. 95 - 114.

11. Челноков Ф. Б. Расчет прохождения импульса в перфорированных трехмерных телах // Труды 12-ой международной конференции «Математика. Компьютер. Образование». — Пущине, 2005. — С. 1Г>7.

12. Петров И. Б., Челноков Ф. Б. Численное исследование прочности железобетонной наружной оболочки под действием динамической нагрузки // Электронная конференция «Топливо и энергетика». — М.: МЭИ, 2004. - С. 56.

13. Челноков Ф. Б. Исследование характера распространения упругих возмущений в перфорированных трехмерных телах // Труды ХЕУГ1 научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук». Часть VII «Управление и прикладная математика». — М.-Долгопрудный:МФТИ. 2004. — С. 32.

14. Агапов П. И., Челноков Ф. Б. Сопоставление нескольких р&шостных схем для численного ]х!шепия двумерных гиперболических уравнений в частных производных // Труды ХЬУ1 научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук». М.-Долгопрудный: МФТИ, 2003. - С. 48.

15. Петров И. Б., Челноков Ф. Б. Моделирование высокоскоростного удара по многослойной преграде с учетом формирования областей разрушения // Труды XIУ научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук». Часть VII «Прикладная математика и экономика». — М.-Долгопрудпый: МФТИ, 2002. — С. 43.

Челноков Федор Борисович

Численное моделирование деформационных

динамических процессов в средах со сложной структурой

Подписано в печать 05.09.2005. Формат 60x84 1/16. Печать офсетная. Усл. печ. л. 1,0. Уч.-изд. л. 1,0. Тираж 70 экз. Заказ N ф-

Государственное образовательное учреждение высшего профессионального образования Московский физико-технический институт (государственный университет) Отдел автоматизированных издательских систем «ФИЗТЕХ-ПОЛИГРАФ» 141700, Московская обл., г. Долгопрудный, Институтский пер., д. 9

»16728

РНБ Русский фонд

2006-4 11267

;

Введение

0.1. Способы дискретизации уравнений механики.

0.2. Способы построения сетки в области интегрирования

0.2.1. Квадратная регулярная сетка.

0.2.2. Регулярная сетка из четырехугольников, сопряженная с границами области интегрирования.

0.2.3. Нерегулярная треугольная сетка.

0.2.4. Изменение сетки при деформировании тел.

Глава 1. Численное решение уравнений упругости.

1.1. Математическая модель.

1.2. Выбор системы координат.

1.3. Обобщение записи дифференциальных уравнений

1.4. Спектральное исследование системы.

1.4.1. Прямая задача

1.4.2. Сопряженная задача.26 *

1.4.3. Нормировка собственных векторов.

1.4.4. Нулевые собственные значения.

1.4.5. Матрицы Л, П, Г2-1 .-.30 ->

1.5. Покоординатное расщепление.

1.6. Разностные схемы.

1.7. Сеточно-характеристические схемы.

1.8. Расчет на границе области интегрирования.

1.8.1. Заданная внешняя сила.

1.8.2. Заданная скорость границы.

1.8.3. Смешанные условия

1.8.4. Условия поглощения и симметрии

1.8.5. Решение на границе при наличии правой части.

1.9. Контакт между двумя телами.54 •

1.9.1. Полное слипание.

1.9.2. Свободное скольжение.

1.10. Интегрирование уравнений акустики.

1.11. Двумерные уравнения упругости.

1.12. Эйлерова сетка и границы из маркеров

Глава 2. Построение нерегулярной треугольной сетки

2.1. Представление триангуляции в программе.

2.1.1. Наиболее компактный формат

2.1.2. Расширенные структуры данных для ускорения поиска

2.2. Триангуляция невыпуклого многоугольника с полостями

2.3. Оптимальная триангуляция Делоне.

2.4. Поддержание заданной плотности сетки

2.4.1. Сокращение граничных вершин.

2.5. Обоснование корректности алгоритма.

2.6. Размеры внутренних треугольников сетки.

2.7. Допустимые размеры длины граничного ребра.

2.7.1. Минимальный угол границы тела.

2.7.2. Обработка треугольников с двумя граничными ребрами

2.8. Трудоемкость поиска треугольника.

2.9. Момент вырождения сетки при движении вершин с постоянной скоростью.

2.10. Примеры работы алгоритмов.

Глава 3. Контакт между телами в динамических задачах

3.1. Поиск сегментов контактирующих границ.

3.1.1. Структуры многомерного поиска.

3.1.2. Триангуляция пространства между телами.

3.2. Проверка сблизившихся узлов.

3.3. Несовпадение узлов в контактирующих телах

3.4. Примеры расчетов с контактом нескольких тел.

Глава 4. Интерполяция в треугольнике.

4.1. Реконструкция полинома заданного порядка.

4.2. Кусочно-линейная интерполяция.

4.3. Градиент интерполяционного полинома.

4.4. Вычисление интеграла полинома в треугольнике.

4.5. Монотонная квадратичная реконструкция.

4.6. Борьба с ростом вариации при интерполяции.

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

4.7.1. Выполнение законов сохранения импульса и энергии

Глава 5. Численный метод для бесструктурных сеток

5.1. Уравнение переноса.

5.2. Гиперболическая система уравнений.

5.3. Сравнение одномерных схем на решении уравнения переноса

Глава 6. Распространение упругих волн в массивных породах

6.1. Введение.

6.2. Начальное состояние среды.

6.3. Граничные условия.

6.3.1. Поверхности трещин.

6.4. Примененная неравномерная треугольная сетка.

6.5. Исследование энергии в области интегрирования.

6.6. Равномерность распределения полостей.

6.7. Оценка вариации плотности тела со случайным распределением круговых полостей.

6.7.1. Полости одного размера.

6.7.2. Случайное распределение размеров полостей.

6.8. Детали численных экспериментов.

6.9. Анализ результатов расчетов.

Глава 7. Распространение волн в перфорированных средах

7.1. Двумерная постановка задачи.

7.2. Трехмерная постановка задачи

Глава 8. Высокоскоростной удар по многослойной преграде

8.1. Постановка задачи.

8.2. Изменение скорости и положения шарика со временем

8.3. Подвижная расчетная сетка.

8.3.1. Перестройка сетки как задача оптимизации.

8.4. Учет разрушения материалов.

8.4.1. Результаты расчетов.

8.4.2. Увеличение рассчитываемого периода соударения за счет фрагментации.

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

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

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

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

Краткие обзоры-классификации используемых численных методов и видов сеток можно найти в работах [5,6].

0.1. Способы дискретизации уравнений механики

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

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

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

Метод конечных разностей (finite-difference method) основывается на . замене производных в поставленных дифференциальных уравнениях на конечные, недифференциальные комбинации (в простейшем случае — разности) значений, определенных в точках сетки. Приближенное решение состоит из значений U^ (в двумерном случае), представляющих поточечное приближение к решению в узлах сетки (Xi, yj) во время tn:

Ujj ~ и(х{, yjt tn).

Метод конечных объемов (finite-volume method), с другой стороны, основывается на приближении средних значений и в ячейках сетки

ГУ] +1/2 r*i+l/2 u(x,y,tn)dxdy.

Jyj-l/2 Jxi-l/2

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

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

Рис. 1. Используемые значения в методах конечных разностей (слева) и конечных объемов (справа), стрелками показаны потоки между ячейками.

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

Метод конечных элементов (finite element method), примененный, например, в [7-9], относится к вариационным методам [10]. Решение в нем ищется среди конечного множества функций, определенных согласованно с текущей сеткой. Метод конечных элементов наиболее часто используется в случае нерегулярности сетки из-за простоты его обобщения на этот случай.

Методы конечных объемов и конечных элементов высоких порядков аппроксимации с упрощенной, аналитической реконструкцией для неструктурированных сеток получили названия методов спектральных объемов (spectral volume method) [11] и спектральных элементов (spectral element method) [12] соответственно. Некоторые соображения, лежащие в их основе, были заимствованы в данной работе для построения сеточно-характе-ристического метода второго порядка на неструктурированных сетках.

Методы частиц [13] и гладких частиц (smoothed particle hydrodynamics) [14-16], или метод свободных точек [17] применяются при моделировании значительных деформаций тел и разрушений в них, приводящих к отколу множества обломков.

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

0.2. Способы построения сетки в области интегриро

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

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

• способностью меняться по мере деформирования моделируемых тел и

• наличию структуры (регулярности) в сетке.

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

0.2.1. Квадратная регулярная сетка вания

Рис. 2. Квадратная регулярная сетка в окрестности исследуемого объекта.

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

Рис. 3. Декартова решетка вблизи границы области интегрирования. Слева отмечен несимметричный шаблон конечно-разностных схем, справа — ячейка неправильной формы конечнообъемных методов,

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

0.2.2. Регулярная сетка из четырехугольников, сопряженная с

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

Рис. 4. Возможные способы построения регулярной сетки в треугольной области.

Существуют алгоритмы [35, 36], позволяющие так разместить точки регулярной сетки (body-fitted grid), чтобы они не выходили за пределы тела, а образованные ячейки между линиями сетки были максимально близки по размерам друг к другу и похожи по форме на квадрат. При этом сетка и шаблоны численных методов не отличаются принципиально внутри и на границе области интегрирования. На рис. 4 представлены возможные виды регулярной сетки для области треугольной формы. Предлагаемые алгоритмы обычно решают задачу оптимизации итерационными методами, последовательно приближаясь к оптимальному значению. Как следствие, полученная сетка существенно зависят от начального приближения, а при плохом выборе не всегда обеспечивается сходимость алгоритма границами области интегрирования и отыскание глобального экстремума. Сложность построения регулярной сетки в телах сложной геометрии заставляет некоторых исследователей [6] возвращаться к декартовой решетке с ее упомянутыми недостатками.

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

• отсутствие острых углов между пересекающимися сеточными линия

• малая кривизна сеточных линий,

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

На рис. 5 представлена еще одна сетка, через каждую точку которой проходит ровно две сеточные линии, однако взаимнооднозначное отображение на решетку отсутствует, и, следовательно, сетка не может считаться регулярной. Такого рода сетки возникают при разбиении исходной области на ряд четырехугольных частей, в каждой из которых сетки строятся независимо, после чего они «сшиваются». Если конечно-разностному численному методу для расчета узла важно лишь наличие четырех расположенных симметрично соседних узлов, то такого рода сетка может быть использована вместо регулярной. Однако, например, для методов, в которых уравнения «расщепляются» вдоль сеточных линий, эта сетка не подойдет. ми,

Рис. 5. Сетка, части которой являются регулярными.

0.2.3. Нерегулярная треугольная сетка

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

Чаще всего используются нерегулярные сетки (рис. 6) с простей- ' шей формой ячейки — треугольником, хотя иногда применяют четырехугольные и смешанные (треугольно-четырехугольные) нерегулярные сетки. Ячейки более сложной формы используются намного реже.

0.2.4. Изменение сетки при деформировании тел

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

В механике сплошной среды принято записывать динамические уравнения в одном из двух видов переменных [5]: эйлеровых — фиксированных в пространстве и лагранжевых — движущихся вместе с точками тела. Аналогично выделяют и два основных типа сеток, но кроме них используют и другие [19]. Рассмотрим далее каждый из них отдельно.

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

• Точки лагранжевой сетки смещаются вместе с точками тела, поэтому границы тела всегда совпадают с сеточными линиями. Однако при наличии сдвиговых деформаций в теле, ячейки лагранжевой сетки могут постепенно вырождаться и пересекать друг друга. Численные методы при вырождении сетки обнаруживают неустойчивость и счет приходится прекращать. Поэтому неизбежен подход, называемый лагран-жева сетка с перестройкой, в котором, как только детектируется приближение сетки к вырожденному состоянию, производится построение новой сетки, а значения в новых узлах интерполируются из прежних значений. Процедура интерполяции сама по себе приводит к потере точности решения, поэтому желательно перестраивать сетку как можно реже. В противном случае имеет смысл обратиться к подвижной сетке. Задачи о соударении в лагранжевых переменных рассматривались в работах [37,38].

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

• По-координатная комбинация разных типов движения сетки рассматривается в статье [19]. Например, если правая и левая границы области интегрирования неподвижны, то дг-координата узлов может быть эйлеровой и не храниться в памяти программы, а у-координата — лагранжевой.

• Совместная эйлерово-лагранжева сетка представлена в работе [39]. Подобная сетка имеет неподвижные точки внутри тела и движущиеся вместе с поверхностью тела точки на ее границе. Основная сложность работы с такими сетками — обеспечить сопряжение двух типов узлов в приграничной полосе. В другой работе [40] исследование задач о соударении в эйлеровых переменных приводило к необходимости использовать метода маркеров.

• Иной вариант совмещения эйлерового и лагранжевого подхода предлагается в методе «частиц в ячейке» Ф. Харлоу [13,17]. Метод использует две сетки частиц: лагранжеву сетку для описания частиц, представляющих элементы жидкости, движущихся по неподвижной эйлеровой сетке. Другой пример совмещения неподвижной сетки с подвижными частицами дает работа [41].

• Неструктурированная лагранжева сетка с перестройкой (изменением «соседства») используется в методике «Медуза», также изложенной в [17]. Область интегрирования разбивается на построенные вокруг заданного набора точек многоугольные ячейки («глобулы»), во многом схожие с ячейками диаграммы Вороного [42]. Заметную трудность о в этой методике представляет расчет граничных ячеек, поскольку их стороны не совпадают с контуром области интегрирования.

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

Основные результаты и выводы диссертации

Аналитическим образом произведено спектральное исследование матриц коэффициентов системы уравнений теории упругости, выписанной в произвольной прямолинейной системе координат. В компактной форме получены выражения для всех собственных значений этих матриц Л (1.22), их левых собственные векторов Q (1.23) и векторов взаимного к ним базиса fi-1 (1.24). В предшествующих работах такие выражения были известны только для декартовой системы координат [55,75], а в прочих случаях опре- . делялись численно [19,22,72].

В работе предлагается использовать явное представление сеточно-ха-рактеристических схем, основываясь на произведенном спектральном исследовании, поскольку в записи таких схем входят громоздкие выражения f относительно A, Q, Q-1. В полученной упрощенной записи не требуется решения системы линейных уравнений, обращения и даже перемножения матриц. Полученные выражения приведены к виду, инвариантному относительно размерности пространства (справедливы в 2D и 3D), тогда как запись A, ft, ft'1 отлична для двумерного и трехмерного пространств [75]. В результате чего вместе с упрощением программы достижима более высокая скорость ее работы, а также исключаются численные ошибки, связанные с решением возможно обусловленных систем линейных уравнений.

Для граничных узлов помимо явной записи было предложено использовать двухэтапный метод, причем первый этап не зависит от граничных условий, а второй — от порядка аппроксимации. Разделение на этапы чрезвычайно удобно в программной реализации, поскольку позволяет отдельно отлаживать компактные модули, отвечающие только за перенос значений вдоль характеристик с тем или иным порядком точности либо только за корректировку для конкретного граничного условия. Метод требует решения лишь системы из m-линейных уравнений, где га-число выходящих из области характеристик, тогда как классический подход [19,22] требует решения полной системы из n-уравнений, где п-число переменных в задаче.

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

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

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

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

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

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

Заключение

1. Седов Л. И. Механика сплошной среды. — М.: Наука, 1970.

2. Новацкий В. К. Теория упругости. — М.: Мир, 1975.

3. Новацкий В. К. Волновые задачи теории пластичности. — М.: Мир, 1978.

4. Партон В. 3Перлин 77. И. Методы математической теории упругости. — М.: Наука, 1981.

5. Кондауров В. И., Фортов В. Е. Основы термомеханики конденсированной среды. — М.: МФТИ, 2002.

6. Бураго 77. Т., Кукуджанов В. 77. Решение упругопластических задач методом конечных элементов, пакет прикладных программ «Астра» // Препринт ИПМ АН СССР. 1988. - № 280.

7. O'Brien J. F., Hodgins J. К. Graphical modeling and animation of brittle fracture // Proceedings of ACM SIGGRAPH. 1999. — Pp. 137 - 146.

8. O'Brien J. F., Hodgins J. K. Animating fracture // Communications of the ACM. 2000. - Vol. 43, no. 7. - Pp. 69 - 75.

9. Рябенький В. С. Введение в вычислительную математику. — М.: Физматлит, 2000.

10. Wang Z. J. Spectral (finite) volume method for conservation laws on unstructured grids // Journal of Computational Physics. — 2002. — Vol. 178. Pp. 210 - 251.

11. Penrose D., ed. Sourcebook of Parallel Computing. — Elsevier Science (USA), 2003.

12. Харлоу Ф. X. Численный метод частиц в ячейках для задач гидродинамики // Вычисл. методы в гидродинамике. — М. :Мир, 1967. — С. 316 -342.

13. Блажевич Ю. В., Иванов В. Д., Петров И. Б., Петвиашвили И. В. Моделирование высокоскоростного соударения методом гладких чи-стиц // Матем. моделирование. — 1999. — Т. 11, № 1. — С. 88 -100.

14. Parshikov A. N., Medin S. A. Smoothed particle hydrodynamics using interparticle contact algorithms // Journal of Computational Physics. — 2002. no. 180. - Pp. 358 - 382.

15. Блажевич Ю. В., Петров И. Б., Сабельник А. Е. Моделирование динамических процессов разрушения пористых конструкций в проблеме безопастности жилищных сооружений. — 2002.http://cs. mipt. ru/docs/whitepapers / petrovl0052002.pdf.

16. Бабенко К. ред. Теоретические основы и конструирование численных алгоритмов задач математической физики. — М.: Наука, 1979.

17. Магомедов К. М., Холодов А. С. О построении разностных схем для уравнений гиперболического типа на основе характеристический соотношений // Ж. вычисл. матем. и матем. физ. — 1969. — Т. 9, № 2. — С. 373 386.i

18. Петров И. Б., Холодов А. С. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом // Ж. вычисл. матем. и матем. физ. — 1984. Т. 24, № 5. - С. 722 - 739.

19. Петров И. Б., Холодов А. С. О регуляризации разрывных численных решений уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. 1984. - Т. 24, № 8. - С. 1172 - 1188.

20. Петров И. Б. Волновые и откольные явления в слоистых оболочках конечной толщины // Изв. АН СССР. Механ. твердого тела. — 1986.- №4.-С. 118-124.I

21. Магомедов К. М., Холодов А. С. Сеточно-характеристические численные методы. — М.: Наука, 1988.

22. Петров И. Б., Тормасов А. Г., Холодов А. С. О численном изучении нестационарных процессов в деформируемых средах многослойной структуры // Изв. АН СССР. Механ. твердого тела. — 1989. — 4. — С. 89 95.

23. Коротин П. И., Петров И. Б., Холодов А. С. Численное решение некоторых задач о воздействии тепловых нагрузок на металлы // Изв. АН СССР. Механ. твердого тела. — 1989. — № 5. — С. 63 69.

24. Коротин П. Н., Острик А. В., Петров И. Б. Численное исследование волновых процессов при объемном энергопоглощении в мишенях конечной толщины // Докл. АН СССР. — 1989. — Т. 308, JV« 5. — С. 1065- 1070.

25. Коротин 77. П., Петров И. Б., Холодов А. С. Численное моделирование поведения упругих и упругопластических тел под воздействием мощных энергетических потоков // Матем. моделирование. — 1989. — Т. 1, № 7. — С. 1 12.

26. Иванов В. Д., Кондауров В. ИПетров И. В., Холодов А. С. Расчет динамического деформирования и разрушения упругопластических тел сеточно-характеристическими методами // Матем. моделирование. — 1990. Т. 2, № 11. - С. 10 - 29.

27. Петров И. БТормасов А. Г. О численном исследовании трехмерных задач обтекания волнами сжатия препятствия или полости в упруго-пластическом пространстве // Докл. АН СССР. — 1990. — Т. 314, № 4. С. 817 - 820.

28. Жуков Д. С., Петров И. В., Тормасов А. Г. Численное и экспериментальное изучение разрушения твердых тел в жидкости // Изв. АН СССР. Механ. твердого тела. — 1991. — № 3. — С. 183 190.

29. Петров И. В., Тормасов А. Г. Численное исследование косого соудаiрения жесткого шарика с двухслойной упругопластической плитой // Матем. моделирование. — 1992. — Т. 4, № 3. — С. 20 27.

30. Иванов В. Д., Петров И. Б. Моделирование деформаций и разрушений в мишенях под действием лазерного излучения // Труды института общей физики. — 1992. — Т. 36. — С. 247 266.

31. Иванов В. Д., Петров И. БТормасов А. Г., Холодов А. С., Пашу-тин Р. А. Сеточно-характеристический метод расчета динамического деформирования на нерегулярных сетках // Матем. моделирование. — 1999. Т. 11, № 7. - С. 118 - 127.

32. Kholodov Y. A Monotone High-Order Accuracy Scheme for Hyperbolic CFD Problems // APS Meeting Abstracts. — 2000. — Pp. B4+.

33. Иваненко С. А., Чарахчьян А. А. Криволинейные сетки из выпуклых, четырехугольников // Ж. вычисл. матем. и матем. физ. — 1988. — Т. 28, № 4. С. 503 - 514.

34. Петров И. Б., Челноков Ф. Б. Численное исследование волновых процессов и процессов разрушения в многослойных преградах // Ж. вычисл. матем. и матем. физ. — 2003. — Т. 43, № 10. — С. 1562 1579.

35. Меньшиков Г. П., Одинцов В. А., Чудов Л. А. Внедрение цилиндрического ударника в конечную плиту // Изв. АН СССР. Механ. твердого тела. 1976. — № 1. — С. 125 - 130.

36. Нох В. Ф. СЭЛ совместный эйлерово-лагранжев метод для расчета нестационарных двумерных задач // Вычисл. методы в гидродинамике.- М. :Мир, 1967. С. 128 - 184.

37. Гриднева В. А., Корнеев А. И., Трушков В. Г. Численный метод расчета напряженного состояния и разрушения плиты конечной толщины при ударе бойками различной формы // Изв. АН СССР, Механ. твердого тела. — 1977. — № 1. — С. 146 157.

38. Кукуджанов В. Н. Метод расщепления упругопластических уравнений // Известия РАН. Механика твердого тела. — 2004. — № 1. — С. 98 108.

39. Lax P. D. Weak solutions of nonlinear hyperbolic equations and their numerical computations // Comm. Pure and Appl. Math. — 1954. — Vol. 7, no. 1. — Pp. 159 193.

40. Lax P. D., Wendroff B. System of conservation laws // Comm. Pure and Appl. Math. 1960. - Vol. 13, no. 2. - Pp. 217 - 237.

41. Courant R., Isaacson E., Rees M. On the solutions of nonlinear hyperbolic differential equations by finite differences // Comm. Pure and Appl. Math.- 1952. — Vol. 5, no. 5. — Pp. 243 254.

42. Агапов 77. И., Челноков Ф. Б. Сравнительный анализ разностных схем для численного решения двумерных задач механики деформируемого твердого тела // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ин-т. — М., 2003. — С. 19 27.

43. Петров И. В., Челноков Ф. Б. Численная проверка прочности железобетонной наружной оболочки под действием динамической нагрузки // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ин-т.- М., 2003. С. 4-13.

44. Агапов П. И., Обухов А. СПетров И. Б., Челноков Ф. Б. Численное решение динамических задач биомеханики сеточно-характеристическим методом // Компьютерные модели и прогресс медицины: Сб. ст. / РАН. — М.: Наука, 2001. С. 275 - 300.

45. Попов И. В., Поляков С. В. Построение адаптивных нерегулярных треугольных сеток для двумерных многосвязных невыпуклых областей // Матем. моделирование. — 2002. — Т. 14, № 6. — С. 25 35.

46. Shewchuk J. R. Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates // Discrete & Computational Geometry. — 1997. — Vol. 18, no. 3. — Pp. 305-363.

47. Кормен Т., Лейзерсон, Pueecm P. Алгоритмы. Построение и анализ. — М.: МЦНМО, 2000.

48. Shewchuk J. R. Delaunay refinement algorithms for triangular mesh generation // Computational Geometry: Theory and Applications. — 2002. — Vol. 22(1-3), no. 5. — Pp. 21 74. http://www.cs.berkcley.edu/~jrs/papers/2dj.ps.

49. Агапов П. И., Петров И. Б., Челноков Ф. Б. Численное исследование задач механики деформируемого твердого тела в неоднородных областях интегрирования // Обработка информации и моделирование: Сб. ст. / Моск. физ.-тех. ин-т. — М., 2002. — С. 148 157.

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

51. Тер-Крикоров А. М., Шабунин М. И. Курс математического анализа. Издание второе, переработанное. — М.: МФТИ, 1997.

52. Warming R. F., Beam R. М. Upwind second-order difference schemes and applications in unsteady aerodynamic low // AIAA 2nd CFD Conf. — Hartford, Connecticut, 1974. — P. 17.

53. Федоренко P. П. Введение в вычислительную физику. — М.: МФТИ, 1994.

54. Русанов В. В. Разностные схемы третьего порядка точности для сквозного счета разрывных решений // Докл. АН СССР. — 1968. — Т. 180, № 6. С. 1303 - 1305.

55. Файзуллин И. С., Куценко Н. В. О возможности применения рассеянных волн для изучения трещиноватости геосреды по данным численного моделирования // Геофизика. — 2004. — № 5. — С. 5 9.

56. Saenger E. H., Kruger O. S., Shapiro S. A. Effective elastic properties of randomly fractured soils: 3D numerical experiments // Geophysical Prospecting. 2004. — Vol. 52. - Pp. 183 - 195.

57. Левянт В. В., Антоненко М. Н., Антонова И. Ю. Исследование методами численного моделирования сейсмического поля, обусловленного рассеиванием на зонах диффузной кавернозности и трещиноватости / / Геофизика. 2004. — № 2. — С. 8 - 20.

58. Бочаров П. П., Печинкин А. В. Теория вероятностей. Математическая статистика. — М.: Гардарика, 1998.

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

60. Майчен Дж.} Сак С. Метод расчета «Тензор» // Вычисл. методы в гидродинамике. — М. :Мир, 1967. — С. 185 211.

61. Кондауров В. И. Континуальное разрушение нелинейно-упругих тел // Матем. моделирование. — 1988. — Т. 52, № 2. — С. 302 310.

62. Белоцерковский О. М. Численное моделирование в механике сплошных сред. — М.: Физматлит, 1994.

63. Белоцерковский О. М. Вычислительная механика. Современные проблемы и результаты. — М.: Наука, 1991.

64. Петров И. Б., Тормасов А. Г. О численном решении пространственных задач соударения // Матем. моделирование. — 1990. — Т. 2, № 2. — С. 58 72.

65. Yngve G. D., O'Brien J. F., Hodgins J. K. Animating explosions // Proceedings of ACM SIGGRAPH. 2000. - Pp. 29 - 36.

66. Андрущенко В. А., Головешкин В. А., Холин H. Я. Вихревые движения твердых сред в динамических задачах теории упругости // Инж.-физ. журнал. 1999. - Т. 72, № 4. - С. 803 - 810.

67. Антоненко М. Н. Численное моделирование распространения упругих волн в неоднородной среде: Диссертация кандидата физ.-мат. наук. — М.: ИАП, 2004.