Содержание к диссертации
Введение
1 Математическое моделирование динамических процессов в средах с микроструктурой на многопроцессорных вычислительных системах 12
1.1 Обзор исследований по моделированию сред с микроструктурой 12
1.2 Обзор методов численного решения динамических задач . 15
1.3 Обзор параллельных вычислительных технологий 21
2 Модель Коссера 25
2.1 Уравнения моментной теории упругости 25
2.2 Корректность задачи 31
2.3 Одномерные движения 33
2.3.1 Плоские продольные волны 33
2.3.2 Поперечные волны с вращением частиц 34
2.3.3 Волны кручения 37
2.3.4 Задача о простом сдвиге 39
2.4 Резонансные спектры поперечных возмущений 40
2.5 Плоское деформированное состояние 47
2.6 Редуцированная среда Коссера 50
3 Численное моделирование 53
3.1 Вычислительный алгоритм 53
3.1.1 Метод расщепления 53
3.1.2 Одномерная ENO-схема 55
3.1.3 Предельная реконструкция 57
3.1.4 Алгоритм реализации граничных условий 59
3.2 Программный комплекс 61
3.3 Алгоритм сжатия файлов 67
4 Результаты численных расчетов 73
4.1 Плоские задачи 73
4.2 Поверхностные волны Рэлея 79
4.3 Резонансные воздействия 82
4.4 Пространственные задачи Лэмба 83
Заключение 88
Литература 89
Приложения 100
- Обзор методов численного решения динамических задач
- Поперечные волны с вращением частиц
- Алгоритм реализации граничных условий
- Поверхностные волны Рэлея
Введение к работе
Актуальность работы. В 2009 году исполнилось 100 лет со дня опубликования работы братьев Коссера, в которой была предложена новая математическая модель сплошной среды. В отличие от классической теории упругости, в этой модели каждая материальная точка наделяется свойствами твердого тела - для нее учитываются вращательные степени свободы. Математическая модель Коссера служит для описания напряженно-деформированного состояния структурно неоднородных материалов: композитов, гранулированных, порошкообразных, сыпучих, микроразрушенных и микрополярных сред. Структура - один из важнейших показателей качества материалов, непосредственно влияющий на их прочностные характеристики. В зависимости от типа материала и масштаба исследований в практических задачах требуется учитывать структуру нано-, микро- или мезоуровня. Особую актуальность математические модели материалов со структурой получили в последнее время в связи с развитием нанотехнологий.
При численном решении задач деформирования в рамках теории Коссера необходимо согласовывать размер ячеек используемых сеток с характерным размером неоднородности, представляющим собой малую величину. В результате дискретизации получаются задачи большой размерности, для реализации которых недостаточно вычислительных ресурсов персонального компьютера или рабочей станции с последовательной архитектурой. Методы моделирования с использованием высокопроизводительных распределенных вычислений оказываются едва ли не единственным способом получения информации об исследуемых процессах.
Целью исследования является разработка и реализация вычислительного алгоритма для решения динамических задач моментной теории упругости, описывающей процессы распространения волн напряжений и деформаций в средах с микроструктурой, на многопроцессорных вычислительных системах.
Для достижения поставленной цели решаются следующие задачи:
Приведение полной системы уравнений моментной теории упругости к симметрической t - гиперболической форме, позволяющей применить к решению задач эффективные вычислительные алгоритмы.
Разработка параллельной версии алгоритмов, ориентированных на использование многопроцессорных вычислительных систем.
Создание комплекса прикладных программ для исследования процессов распространения упругих волн в средах с микроструктурой на кластерных системах с распределенной памятью.
В качестве метода исследования используется вычислительный эксперимент, включающий в себя следующие этапы: математическая формулировка задачи, построение численного алгоритма решения, программная реализация алгоритма, проведение расчетов, анализ полученных результатов.
Новые научные результаты, выносимые на защиту:
Разработан параллельный вычислительный алгоритм для решения динамических задач моментной теории упругости, основанный на расщеплении пространственной задачи на серию одномерных задач.
Алгоритм реализован на языке программирования Fortran-90 с использованием библиотеки передачи сообщений MPI в виде комплекса программ для многопроцессорных вычислительных систем.
На основании серии расчетов показано, что в моментной упругой среде существует собственная резонансная частота, зависящая только от инерционных свойств частиц микроструктуры и от параметров упругости материала.
Личный вклад. Результаты, составляющие основное содержание диссертации, получены автором самостоятельно. В совместных работах соавторам принадлежит постановка задачи, автором диссертации проведены необходимые численные расчеты и обработка полученных результатов.
Достоверность и обоснованность результатов обеспечивается применением строгих математических методов исследования; использованием при
компьютерном моделировании тестовых задач, допускающих точное аналитическое решение.
Научная новизна работы заключается в том, что в ней впервые выполнена численная реализация модели моментной среды Коссера на многопроцессорных вычислительных системах в пространственной постановке и показано, что в такой среде существует собственная резонансная частота.
Практическая ценность работы состоит в создании комплекса прикладных программ, который может быть использован для численного исследования волновых процессов в средах с микроструктурой в задачах сейсмики и акустики, а также в учебном процессе при подготовке специалистов по математической обработке геофизической информации.
Апробация работы. Основные результаты исследований были представлены на научных конференциях: XLIII международная научная студенческая конференция «Студент и научно-технический прогресс» (Новосибирск, 2005); VI школа-семинар «Распределенные и кластерные вычисления» (Красноярск, 2006); VII, VIII, IX Всероссийские конференции молодых ученых по математическому моделированию и информационным технологиям (Красноярск, 2006, Новосибирск, 2007, Кемерово, 2008); XV, XVI зимние школы по механике сплошных сред (Пермь, 2007, 2009); Конкурс-конференция молодых ученых Красноярского научного центра (Красноярск, 2007); II международная научная конференция «Параллельные вычислительные технологии» (Санкт-Петербург, 2008); V Российско-германская школа по параллельным вычислениям (Новосибирск, 2008); Всероссийская молодежная школа по параллельному программированию (Новосибирск, 2009); Международная конференция «Математические и информационные технологии» (Будва, Черногория, 2009); Первая Всероссийская конференция «Проблемы механики и акустики сред с микро- и наноструктурой» (Нижний Новгород, 2009).
Публикации. По теме диссертации опубликовано 8 печатных работ, из них 2 статьи в изданиях по списку ВАК.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 08-01-00148), Комплексной программы фундаментальных исследований Президиума РАН № 17 «Параллельные вычисления на многопроцессорных вычислительных системах», № 14 «Фундаментальные проблемы информатики и информационных технологий» и № 2 «Интеллектуальные информационные технологии, математическое моделирование, системный анализ и автоматизация», Междисциплинарного интеграционного проекта СО РАН № 40 и Красноярского краевого фонда науки (грант 17G029).
Объем и структура работы. Диссертация состоит из введения, четырех глав, заключения, списка литературы и приложений. Общий объем диссертации составляет 105 страниц, включая 42 рисунка и 3 таблицы. Список используемой литературы содержит 96 наименований.
Обзор методов численного решения динамических задач
Многообразие численных алгоритмов решения задач динамики возникает из выбора схем расчета и сетки. Выделяются следующие основные группы алгоритмов: разностные методы, в том числе методы, основанные на точном и приближенном решении задачи Римана, методы с выделениями разрывов, методы сквозного счета; методы конечных элементов; методы частиц.
Численный метод должен сохранять свойство монотонности в тех областях, где искомые решения имеют большие перепады значений, и должен обладать высоким порядком точности в тех областях, где решение является гладким. Доказанная для линейного уравнения переноса теорема С. К. Годунова гласит, что среди линейных разностных схем не существует монотонных с порядком выше первого.
Для решения многомерных задач математической физики применяют процедуры расщепления, заключающиеся в сведении решения исходной задачи, описывающей сложный физический процесс, к решению последовательности задач, описывающих процессы более простой физической структуры [26,27].
Схема Годунова [28] первого порядка точности с расщеплением по пространственным переменным использует дивергентную запись уравнений. В ней все функции, определенные на сетке точек, предполагаются кусочно-постоянными. Важным свойством этой схемы является ее монотонность. Повышение точности схемы Годунова развивается по нескольким направлениям: выделение сильных разрывов, которые могут двигаться по разностной сетке, построение схем переменного порядка точности (гибридных и TVD-схем), схем с адаптивными к решению разностными сетками и другие [29].
Гибридные разностные схемы являются одним из промежуточных этапов развития численных методов. Под термином «гибридность» понимается возможность численного метода менять свои свойства, в частности, порядок аппроксимации. Гибридные схемы, обеспечивающие повышение точности, были предложены Федоренко [30]. В этих схемах сильные разрывы рассчитываются по схеме первого порядка, а гладкие решения - по схеме второго порядка точности. Данный подход получил развитие в работах [31, 32], где предложен метод коррекции потоков или метод FCT (Flux Correct Transport). На первом шаге вычисления решения используется монотонная схема первого порядка точности. На втором шаге полученное численное решение должно быть модифицировано с тем, чтобы повысить его порядок до второго по времени и пространству. На этом шаге в численном решении не должны возникать новые локальные экстремумы, а также возрастать (уменьшаться) значения локальных максимумов (минимумов), которые имели место в начале этого шага. Такие условия эквивалентны условию неувеличения полной вариации численного решения, или условию TVD (total variation diminishing), более слабому, чем требование монотонности схемы [33]. На основе TVD-подхода развита специальная процедура кусочно-линейной (кусочно-полиномиальной) реконструкции сеточных функций внутри ячейки, позволяющая повысить точность в разностных схемах Годунова.
Модификации схемы Годунова, повышающие порядок ее аппроксимации до второго на гладких решениях (Схема Лакса-Вендрофа [34], Русанова [35] и гибридные схемы), уже не обладают свойством монотонности или имеют смешанный порядок.
Для поиска точного решения задачи Римана, необходимого для реализации многих схем, обладающих монотонностью, сначала следует построить основные элементарные решения: движущийся разрыв и волны Римана. Затем точное решение общего вида строится на основе их комбинации.
К методам, основанным на приближенном решении задачи Римана, относятся метод Куранта-Изаксона-Риса (КИР), Роу и Ошера [36-39]. Методы КИР и Роу основаны на решении задачи Римана для локально линеаризованной системы уравнений. В этом случае решение состоит только из движущихся разрывов, которые разделяют области с постоянными значениями величин. В методе Ошера приближенное решение задачи Римана строится для квазилинейной системы уравнений. При. этом решение является комбинацией только волн Римана.
Поперечные волны с вращением частиц
В работе [45] предложен и обоснован класс экономичных разностных схем для динамических задач теории упругости. В основу положены эквивалентные сопряженно-факторизованные постановки задачи в перемещениях и в напряжениях. Доказано, что на решениях задачи в напряжениях тензор несовместности равен нулю. Это позволило в достаточно общем случае (криволинейная система координат, основные типы краевых условий) построить самосопряженный положительно-определенный сеточный оператор задачи в напряжениях и получить основные теоремы сходимости.
В основе метода конечных элементов лежит идея расчета сложной конструкции путем расчленения ее на отдельные простейшие части (конечные элементы), напряженно-деформированное состояние которых сравнительно легко описывается, и затем объединения их вновь в единую конструкцию с выполнением условий равновесия и непрерывности поля перемещений [46-49]. Основное преимущество этого метода по сравнению с методом конечных разностей состоит в возможности сгущения расчетных сеток в локальных подобластях. Метод конечных элементов сочетает идеи, лежащие в основе разностных и вариационных методов. С одной стороны, подобно разностным схемам, на рассматриваемую область наносится сетка разбиения. С другой стороны, при теоретическом обосновании этой схемы применяется подход, использующий вариационную интерпретацию метода.
В настоящее время метод конечных элементов широко применяется для исследования самых разнообразных физических процессов, в задачах статики и динамики [47,50]. Подтверждением развития метода служат многочисленные пакеты прикладных программ.
Метод частиц состоит в представлении тела совокупностью взаимодействующих частиц (материальных точек или твердых тел), описываемых законами классической механики. Одним из наиболее хорошо разработанных вариантов этого метода является метод молекулярной динамики [51], интенсивно использующийся для исследования физико-химических свойств материалов. В классической молекулярной динамике в качестве частиц выступают атомы и молекулы, составляющие материал. Взаимодействие частиц описывается посредством потенциалов взаимодействия, основным свойством которых является отталкивание при сближении и притяжение при удалении. Перед началом моделирования задается некоторое начальное распределение частиц в пространстве (исходная структура материала) и начальное распределение скоростей частиц (механическое и тепловое движение системы в исходном состоянии). Далее задача сводится к решению задачи Коши для системы обыкновенных дифференциальных уравнений [52-56]. В силу ограниченности радиуса взаимодействия между частицами, этот метод допускает почти полное распараллеливание процессов, происходящих в смежных областях пространства, что позволяет эффективно применять данный метод при моделировании процессов сильного деформирования и разрушения на многопроцессорных вычислительных системах.
При решении задач механики неограниченная область задания функции заменяется расчетной областью конечных размеров. В связи с этим, возникает проблема постановки условий на искусственных (фиктивных) границах этой области. Формулировка таких условий представляет собой сложную проблему, которую можно считать в какой-то степени решенной только для уравнений линейной теории упругости и вязкоупругости.
В общем случае нужно либо отображать неограниченную область на конечную путем преобразования системы координат, либо располагать границы расчетной области достаточно далеко от приложенных воздействий [57]. Другой путь состоит в том, чтобы на близких искусственных границах ставить такие условия, которые не искажали бы решение внутри области, т.е. обеспечивали бы полное или частичное поглощение приходящих возмущений. Граничные условия, которые подавляют или уничтожают приходящие волны, называются неотражающими, поглоща- ющими или излучательными. В работе [58] изучаются такие условия для линейных гиперболических систем. Неотражающие граничные условия должны быть простыми по реализации. В противном случае не получится выигрыша во времени счета по сравнению с простым расширением области.
Разнообразие существующих подходов к конструированию неотражающих искусственных граничных условий принято разделять на три группы: локальные условия, нелокальные условия и поглощающие слои.
В методе «идеально-согласованных поглощающих слоев» PML (Perfectly Matched Layer) расчетная область окружается слоем некоторой специальной модельной среды, имеющим конечную толщину [59,60]. Параметры слоя выбираются таким образом, чтобы при прохождении через него волн, распространяющихся РІЗ области решений, не возникало отраженной волны. На внешней границе такого слоя ставится корректное краевое условие, удобное для аппроксимации в выбранной разностной схеме.
Алгоритм реализации граничных условий
Рассмотрим алгоритм вычисления вектора U на правой границе области решения задачи (на левой границе этот вектор рассчитывается аналогично). Граничные условия представимы в виде где m+ - количество строго положительных значений Q, соответствующих «уходящим» характеристикам; Dj и qi - заданные векторы и коэффициенты. Недостающие т — т+ уравнений для отрицательных и равных нулю значений Ch (h = т+ + 1,..., т) имеют вид где Г]г - выражения в правой части равенства (3.5). Очевидно, что к (3.7) приводятся линейные граничные условия общего вида, в частности, условия в скоростях или в напряжениях. Граничные условия (3.7) корректны, если т векторов D/ и Y A образуют линейно независимую систему. В этом случае приведенные уравнения позволяют найти все инварианты її для I = 1, ...,т, после чего U определяется из разложения по базису Y/.
Кроме рассмотренных выше типов граничных условий на практике часто используются условия на искусственных границах, возникающих из-за необходимости урезания расчетной области. Граничные условия, моделирующие беспрепятственное прохождение волн через границу, называются неотражающими, поглощающими или излучательными [39]. Такие условия обеспечивают затухание волн, приходящих в расчетную область. При численном исследовании распространения волн в неограниченном пространстве использовался простейший вариант неотражающих условий, состоящий в следующем.
Пусть в момент времени t в точку х на искусственной границе приходит волна, движущаяся в направлении вектора п. В конкретных задачах этот вектор можно определить, например, приближая волну плоской сферической или цилиндрической, распространяющейся от заданного источника. Вектор п остается постоянным вплоть до прихода в точку новой волны и по мере прохождения волны он изменяется. В первом приближении можно считать, что в окрестности точки х волна падает на границу под прямым углом. В этом случае вектор п совпадает с нормалью к границе, а искусственные граничные условия записываются в форме уравнений для инвариантов, которые переносят возмущения в направлении, противоположном п: Такие уравнения использовались при реализации граничных условий для одномерных задач метода расщепления совместно с уравнениями, задающими значения инвариантов Д на приходящих характеристиках. Получаемая система линейных алгебраических уравнений представляет собой частный случай системы (3.7), ее решение может быть осуществлено с помощью прямых методов, не требующих итераций.
Изложенный вычислительный алгоритм реализован в виде комплекса прикладных программ для численного решения плоских и пространственных задач динамики моментных сред на многопроцессорных вычислительных системах. Программирование выполнено по технологии SPMD (Single Program - Multiple Data) на языке Fortran—90 с применением библиотеки передачи сообщений MPI (Message Passing Interface) [78,79].
Программный комплекс позволяет проводить расчеты распространения упругих волн, вызванных внешними воздействиями, в массиве среды, составленном из разнородных блоков с криволинейными границами. Он состоит из программы-препроцессора, основной программы и программы-постпроцессора.
Информация, необходимая для расчетов, представляется в виде текстовых файлов. Один из таких файлов содержит феноменологические параметры материалов в определенной системе координат (например, СИ), в другом хранятся сведения о блочной структуре массива, идентификатор материала в блоке и пространственные размерности сеток. Третий файл содержит информацию о внешних механических воздействиях: зона приложения нагрузки, характер (распределенная, сосредоточенная), вид (периодическая, импульсная и т.п), длительность нагружения.
Поверхностные волны Рэлея
Такую нагрузку можно создать, например, через жесткую пластинку, приклеенную к образцу, с помощью периодически повторяющегося поперечного удара. Остальные грани являются неотражающими.
В результате воздействия в среду распространяется серия волн на-гружения и разгрузки, представляющих собой области плавного изменения решения с четкими передними и задними фронтами. Это видно из рис. 4.1, в котором для двух моментов времени t = 167 мс (слева) и t — 372 мс (справа) изображены линии уровня касательного напряжения т\2 (оттенкам серого соответствует диапазон изменения напряжения от —0.6 до 0.3). В первом случае (рис. 4.1, слева) наблюдается одна волна, инициируемая первым импульсом нагрузки на границе, во втором (рис. 4.1, справа) - три волны от действия трех импульсов. В отличие от одномерного решения, в котором касательное напряжение всюду положительно (рис. 2.3), в плоской задаче из-за боковой разгрузки появляются зоны отрицательного значения.
На рисунках, соответствующих угловой скорости частиц UJ$ (рис. 4.2) и моментному напряжению шіз (рис. 4.3) для тех же моментов времени, наблюдаются осцилляции, характерный масштаб которых может быть
На рисунках 4.4 и 4.5 представлены аналогичные результаты для нормального напряжения Гц и линейной скорости V2- Линии уровня напряжения лежат в диапазоне от -0.3 до 0.35, а линии уровня скорости -от —5 Ю-7 до 8 Ю-7. Соответствующие волны генерируются в точках, ограничивающих зону приложения нагрузки на левой стороне массива.
Аналогичные результаты о действии нормальной нагрузки приведены на рис. 4.6 и 4.7 в момент времени 167 мс.
Расчеты проводились для квадрата со стороной 0.1 м. Характерный масштаб микроструктуры материала равен 0.15 мм. Используемая равномерная сетка составляла 1000 х 1000 ячеек с шагом 0.1 мм. Шаг по времени составил 372 мкс. Время действия импульса равно 1 мс, время между соседними импульсами также равно 1 мс. Время расчета задачи составило 6 ч 45 м на 10 процессорах.
На рис. 4.8 - 4.9 приведены результаты численного решения задачи Лэмба о мгновенном действии сосредоточенной силы на поверхности полупространства при плоской деформации. Размерность сетки составила 1000 х 500 ячеек. Шаг по времени, допустимый по условию Куранта-Фридрихса-Леви, равен 372 мкс. В первом случае импульсная нагрузка приложена в центре верхней границы расчетной области в направлении нормали (рис. 4.8), а во втором - по касательной к границе (рис. 4.9).напряжения Гц, касательного напряжения ті2, скорости vi и угловой скорости и3
На рисунках слева направо изображены линии уровня нормального и касательного напряжения, линейной и угловой скорости, соответственно. На линиях уровня четко прослеживаются все волны, характерные для решения задачи Лэмба в рамках классической теории упругости [89,90]: падающие продольная и поперечная волны с круговыми фронтами, две головные поперечные волны в виде симметричных отрезков прямых, касающихся полуокружности меньшего радиуса, а также быстро затухающие с глубиной поверхностные волны Рэлея - яркие точки на границе, движущиеся вслед за падающей поперечной волной. Отличие от классической теории упругости состоит в том, что в модели моментнои упругой среды решение имеет ярко выраженный колебательный характер.
Рассмотрим полупространство, поверхность которого свободна от нагрузок в случае отсутствия массовых сил и моментов. Оси декартовых координат Х\ и х2 направлены по поверхности, а ось х3 _ вглубь полупространства. Решение системы (2.1), описывающее поверхностную волну Рэлея [91,92], представимо в виде где г = у—1; к - волновое число; / - круговая частота; щ(х3), щ(х3), 2(#з) _ зависящие от глубины функции затухания, а физический смысл имеют только вещественные части данных комплекснозначных функций.
