«РАТАУШКО ЯН ЮРЬЕВИЧ ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ УПРУГИХ И ПОРОУПРУГИХ ТРЁХМЕРНЫХ ТЕЛ НА ОСНОВЕ СОВМЕСТНОГО ПРИМЕНЕНИЯ МЕТОДОВ ГРАНИЧНЫХ ЭЛЕМЕНТОВ И РУНГЕ-КУТТЫ Специальность 01.02.04 – механика деформируемого ...»

В пороупругой динамике находит применение также непрямой подход метода ГИУ [179, 180]. В [196-198] с его помощью рассматривается задача о жёстком включении в пороупругое полупространство. В [152] на базе непрямой формулировки МГЭ для пороупругости проводится исследование рассеяния сдвиговых волн.

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

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

Методика исследований основана на граничных интегральных уравнениях прямого подхода трёхмерных изотропных линейных теорий упругости и пороупругости в преобразованиях по Лапласу; на описании пороупругой среды моделью Био с четырьмя базовыми функциями – три компоненты перемещений упругого скелета и поровое давление; на получении решений во времени на основе шагового метода численного обращения преобразования Лапласа на узлах методов Рунге-Кутты; на компьютерном моделировании искомых решений методом граничного элементав сочетании с методом коллокации; локальной поэлементной аппроксимацией на основе согласованной модели интерполирования Гольдштейна [28].

Научную новизну работы составляют: гранично-элементное моделирование краевых задач смешанного типа динамики трёхмерных упругих и пороупругих тел в сочетании с шаговым методом численного обращения преобразования Лапласа на узлах семейства схем Рунге-Кутты; применение в компьютерном моделировании трёхмерной динамики согласованной модели аппроксимации Гольдштейна на обобщённых четырёхугольных элементах; решение на основе применения метода гранично-временных элементов совместно с методами Рунге-Кутты волновых задач о действии силы в виде функции Хевисайда по времени на упругое и пороупругое призматические тела, полупространство, полупространство, ослабленное полостью, и полупространство с выемкой; исследование возбуждения медленной волны в пороупругой среде с помощью шаговой схемы МГЭ на узлах методов Рунге-Кутты.

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

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

Основные положения, выносимые на защиту трёхмерных упругих и пороупругих тел на основе совместного применения метода гранично-временных элементов с шаговой схемой на узлах методов Рунге-Кутты.

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

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

Моделирование решений следующих волновых задач на основе шаговой схемы на узлах методов Рунге-Кутты:

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

о действии вертикальной силы в виде функции Хевисайда по времени на поверхность однородного и слоистого пороупругих полупространств;

о действии вертикальной силы в виде функции Хевисайда по времени на поверхность пороупругого полупространства, ослабленного полостью, и пороупругого полупространства с выемкой.

Результаты диссертационной работы докладывались на X Всероссийском съезде по фундаментальным проблемам теоретической и прикладной механики (Н.Новгород, 2011), XXV Международной конференции «Математическое моделирование в механике деформируемых сред и конструкций. Методы граничных и конечных элементов» (СанктПетербург, 2013), VII Всероссийской (с международным участием) конференции по механике деформируемого твердого тела (Ростов-на-Дону, 2013), XVIII Нижегородской сессии молодых ученых – естественные, математические науки (Н. Новгород, 2013), форуме молодых ученых (Н. Новгород, 2013), XX Международном симпозиуме «Динамические и технологические проблемы механики конструкций и сплошных сред»

им. А.Г. Горшкова (Вятичи, 2014), V Международной инженерной конференции Университета Кхон Кэна KKU-IENC2014 «Engineering and Technological Responses to Global Challenges» (Кхон Кэн, Таиланд, 2014), III Международном симпозиуме «Physics and Mechanics of New Materials and Underwater Applications (PHENMA 2014)» (Кхон Кэн, Таиланд, 2014), VIII Всероссийской молодежной научно-инновационной школе (Саров, 2014).

Опубликовано 16 работ [37-45, 57-60, 134-136], из них 16 по теме диссертации. В изданиях, рекомендуемых ВАК РФ для опубликования результатов кандидатских диссертаций, опубликовано 5 работ в соавторстве [38, 40, 42, 43, 57]. Результаты работ принадлежат Ратаушко Я.Ю., кроме постановок задач и постпроцессорных представлений.

Структура и объем работы Диссертационная работа состоит из введения, трёх глав, заключения, списка литературы из 218 наименований. Общий объем диссертации составляет 179 страниц машинописного текста, включая 235 рисунков и 2 таблицы.

На различных этапах работа поддерживалась грантами Президента РФ для государственной поддержки ведущих научных школ (№ НШ-2843.2012.8 2012-2013гг.; № НШ-593.2014.8 2014-2015гг.); средствами ФЦП "Научные и научно-педагогические кадры инновационной России" на 2012 - 2014 годы (ГК №14.1337.21.1249 от 14 сентября 2012г.);

грантами РФФИ (проекты № 12-01-00698-а, № 13-08-00742-а, № 14-08-00811-а, № 14-08мол_а).

Введение содержит описание основного направления исследований в диссертационной работе; краткое состояние вопроса по модели пороупругости Био, подходам ГИУ и МГЭ и их использованию для решения динамических задач, в частности пороупругости; обоснование актуальности темы работы; формулировки целей работы и основных положений, выносимых на защиту. Также приведена информация об объёме и структуре работы, о конференциях и публикациях, где были представлены результаты диссертационной работы.

В главе I представлены система уравнений Ламе для упругой среды и система уравнений модели Био пороупругой среды, постановки краевых задач в преобразованиях по Лапласу для обеих моделей. Изложены основные принципы метода квадратур свёрток, а также его модификации с использованием методов Рунге-Кутты. Приведены методика численного обращения интегрального преобразования Лапласа с помощью как традиционного шагового метода, так и шагового метода на узлах семейства схем РунгеКутты, модификация методики с переменным шагом интегрирования при расчёте коэффициентов квадратурной формулы, а также модификация схемы «с ключом» на основе формул интегрирования сильно осциллирующих функций. Для исследований выбрана шаговая схема численного обращения преобразования Лапласа. Проведено сравнение подходов к численному обращению преобразования Лапласа на основе шаговых схем на узлах метода Эйлера и других методов семейства Рунге-Кутты на примере численно-аналитических решений задач о действии осевой нагрузки на упругий и пороупругий стержни. На примере задачи о действии осевой нагрузки на пороупругий стержень проведено численно-аналитическое моделирование эффекта возбуждения медленной продольной волны в пороупругой среде с использованием шаговой схемы на узлах методов Рунге-Кутты для получения оригиналов динамических откликов порового давления и потока; приведено сравнение решений с результатами других авторов и сделаны выводы о превосходстве полученных результатов над ранее опубликованными.

В главе II представлена система ГИУ прямого подхода в изображениях по Лапласу и гранично-элементная методика решения ГИУ, предполагающая использование шаговой схемы численного обращения преобразования Лапласа. Запись системы ГИУ в изображениях по Лапласу позволяет для получения оригиналов искомых функций построить шаговую схему, опирающуюся как на теорему о свёртках оригиналов, так и на теорему об интегрировании оригинала. В работе выбран подход, опирающийся на теорему об интегрировании оригинала, и на её основе построены шаговые схемы на узлах методов Рунге-Кутты для решения дискретных аналогов краевых задач теории упругости и пороупругости. Для гранично-элементной дискретизации использованы четырёхугольные восьмиузловые биквадратичные элементы, применяется метод коллокации.

Аппроксимация обобщённых граничных функций построена по согласованной модели.

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

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

полупространстве. Решены задачи о действии вертикальной силы на поверхность однородного пороупругого полупространства и пороупругого полупространства с фиктивной границей (постановка для составного тела). Задача использована для выбора рабочей гранично-элементной сетки: приведено исследование сходимости решений на разных пространственных сетках и проведено сравнение с решением другого автора.

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

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

диссертации.

Глава I обращения преобразования Лапласа 1.1. Математические модели 1.1.1. Упругая среда Начально-краевую задачу теории упругости составляет система уравнений Ламе:

где L0 G ( K )grad div – эллиптический оператор изотропной упругой статики, дополненная условиями на границе:

и начальными условиями:

Применение преобразования Лапласа сводит краевую задачу теории упругости к задаче эллиптического типа:

где u – граница Дирихле, а – граница Неймана.

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

где V f – общий объём пор (наполнителя), V – общий объём рассматриваемого тела.

Наполнитель может свободно перемещаться внутри пор; закрытые поры считаются частью упругого скелета. Насыщение наполнителем означает совпадение суммарного объёма пор и объёма наполнителя и гарантирует равенство V V f V s, где V s – объём упругого скелета без учёта пор.

Идея заключается в том, чтобы получить описание такой двухфазной среды через характеристики каждого компонента ( K s, K f – соответственно упругие объёмные модули скелета и наполнителя) и усреднённые характеристики пористой среды ( K, G – упругие константы пористой среды). В работе [85] Био вводятся суммарное напряжение и коэффициент активного напряжения :

На основе этих формул получено определяющее соотношение в перемещениях упругого скелета и поровом давлении:

Более поздняя работа [86] рассматривает также альтернативное определяющее соотношение, содержащее перемещения упругого скелета и относительные смещения наполнителя. Таким образом, существует несколько разных формулировок теории Био с четырьмя или шестью неизвестными [87, 89, 217]. Однако в работе [89] Bonnet G. показал, что для описания среды достаточно набора компонент перемещения упругого скелета и порового давления.

Для линейной модели пороупругости справедливо Взаимодействие между упругим скелетом и наполнителем описывается набором постоянных, R, где Вводится также уравнение неразрывности для изменения объёма наполнителя:

где – изменение объёма наполнителя на единицу объёма, – поровый поток.

В дополнение к уравнению баланса движения наполнителя вводится общее уравнение баланса среды:

где общая плотность вычисляется как:

а также используется обобщённый закон Дарси:

Здесь – коэффициент поровой проницаемости, а – плотность присоединённой массы – мнимая массовая плотность, введённая Био в [87] для описания динамического воздействия между текучим наполнителем и упругим скелетом. Коэффициент C зависит от геометрического расположения пор и частоты возбуждения. Bonnet G. and Auriault J.-L. [90] было установлено C 0.66 на низких частотах возбуждения для совокупности сферических пор, расположенных вплотную друг к другу. На высоких частотах предлагается функциональная зависимость коэффициента C от частоты возбуждения [88, 90].

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

где Исключая vi, получим окончательную систему дифференциальных уравнений полной модели Био в преобразованиях по Лапласу:

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

Полагая несжимаемую модель, т.е. = 1 и R, можем получить следующую систему уравнений движения для несжимаемого материала:

Поровое давление для такой модели определяется только деформациями скелета и более не образует степень свободы. Также находит применение упрощённая модель Био [216], основанная на упрощённом уравнении баланса:

и упрощённом законе Дарси:

что позволяет сформулировать дифференциальные уравнения движения во времени:

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

1.2. Метод квадратур свёрток и шаговая схема численного обращения преобразования Лапласа Прямое и обратное интегральные преобразования Лапласа соответственно определяются формулами:

где s i – комплексная переменная, введённая в полуплоскости 0.

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

1.2.1. Традиционный метод квадратур свёрток [156, 157] Метод квадратур свёрток рассматривает интеграл свёртки:

Интеграл свёртки (1) можно преобразовать следующим образом [82, 188]:

дифференциального уравнения первого порядка После дискретизации времени t на N равных временных шагов t интеграл свёртки при каждом дискретном значении времени t nt определяется следующим образом:

Решение дифференциального уравнения (3) может быть аппроксимировано посредством линейного многошагового метода где x(nt, s) xn есть приближенное решение x(t, s) при заданном дискретном времени nt. Однако данное представление многошагового метода непригодно для получения дискретных значений xn с последующей подстановкой их в уравнение (4). Поэтому уравнение (5) умножается на z n ( z C ) и суммируется по n от 0 до :

Исходя из условия, что заданы однородные начальные значения, сделаем следующие преобразования:

Так как g (t 0) 0, то аналогично сделаем преобразования в (6) (подчёркнутая сумма). Таким образом, следующие суммы могут быть вынесены за скобку, и выражение (6) преобразуется следующим образом:

Введём обозначение Функция (z ) определяется многошаговым методом (5) и является его преобразуя (7), можем получить представление функции x(t, s) в виде формальных степенных рядов [188]:

Прежде, чем полученное представление x(nt, s) может быть подставлено в (4), умножим его на z n и просуммируем по n. Окончательно, в интеграле свёртки:

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

Для корректного использования данного приёма требуется, чтобы функция f удовлетворяла следующим условиям [82, 188]:

то есть функция f должна быть ограничена в бесконечности, в связи с чем интеграл по полуокружности обращается в 0 при R. Таким образом, дальнейшее интегрирование будет осуществляться с использованием теоремы о вычетах. По определению обратного преобразования Лапласа все особенности функции f должны располагаться в левой полуплоскости относительно прямой, проходящей через c. В таком случае единственная интегрирования, это:

Следуя теореме о вычетах [156], комплексный интеграл (8) представим следующим образом:

Для получения представления y(nt ) правую часть (9) запишем в виде ряда с f ( ( z ) t ) в степенной ряд:

Коэффициенты n (t ) из (10) определяются двумя способами:

сравнением коэффициентов, когда функция f перегруппируема в виде ряда;

использованием для f интегральной формулы Коши:

где R – радиус окружности области аналитичности функции f ( ( z ) t ).

Подстановка степенного ряда (10) и соответствующие преобразования по формуле произведения рядов Коши приводят к следующей двойной сумме:

Подставляя это выражение в (9), имеем:

Сравнение соответствующих коэффициентов перед одинаковыми степенями приводит к следующему результату [188, 215]:

Уравнение (12) представляет аппроксимацию интеграла свёртки (1), построенную на основе изображения по Лапласу f (s) и на основе дискретных значений другой свёртываемой функции g (kt ).

Единственная аппроксимация, которая использовалась при выводе формулы (12), заключалась в использовании линейного многошагового метода (5) для нахождения непосредственно. Многошаговый метод должен быть порядка точности p 1, кроме того, полуплоскости относительно прямой c i, c i, то есть:

Однако если функция f (s) аналитична и ограничена в области где, критерий устойчивости может быть ослаблен до A( ) устойчивости.

Выражаясь в терминах характеристической функции (z ), имеем [82, 188]:

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

Если функция f ( ( z ) t ) не может быть представлена в виде аналитического степенного ряда, то коэффициенты степенного ряда должны быть посчитаны с помощью интегральной формулы Коши (11). Вычисления интеграла в формуле (11) для определения весов интегрирования n представимы численно. После приведения (11) к полярным координатам z Rei получим:

Этот интеграл аппроксимируется с L равными шагами 2 L :

Математическое доказательство сходимости и устойчивости метода квадратур свёрток могут быть найдены в работах [156, 157]. При условии того, что функция f (s) погрешность вычисления n порядка ( ).

1.2.2. Метод квадратур свёрток на основе методов Рунге-Кутты Для решения уравнения (3) рассмотрим метод Рунге-Кутты, представленный с помощью таблицы Бутчера функцию устойчивости Для корректной формулировки метода квадратур свёрток должны быть выполнены следующие условия [82]:

1. Метод Рунге-Кутты должен быть A-устойчивым;

Будем считать, что bT A1 (0,...,0,1), тогда метод L-устойчив, и значит, Aустойчив. В методах семейства Рунге-Кутты вычисления производятся не только в точках t n nt, но и в промежуточных t n сl t, l 1,2,..., m. Применяя метод Рунге-Кутты для решения уравнения (3), получим [82]:

Значения g (t n ci t ) записаны в вектор g n. Для вывода метода квадратур полиномов. Для метода Рунге-Кутты можем записать:

Умножив на z C и просуммировав, получим:

Учитывая A и L – устойчивости, окончательно получим, что характеристическая функция имеет следующий вид [82]:

Подставляя (14) в (13), получим окончательную формулу для решения дифференциального уравнения (3):

Теперь (2) можно переписать в виде:

Теперь запишем разложение функции f ( ( z ) t ) в степенной ряд:

Коэффициенты n (t ) определяются из интегральной формулы Коши:

Интеграл для определения весов интегрирования n представим численно. После преобразования (17) к полярным координатам z Rei получим:

или, аппроксимируя с L равными шагами 2 L :

Подставляем (16) в (15) и окончательно получаем [80]:

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

Заменим интеграл (18) квадратурной суммой, весовые множители которой определяются с помощью изображения по Лапласу f и линейного многошагового метода (с учётом результатов, полученных в [156, 157, 188]):

где N – количество шагов по времени, L – количество расчётных узлов для численного интегрирования по аргументу.

Аппроксимация, используемая при выводе формул (19), (20), основана на применении линейного многошагового метода (с характеристической функцией (z ) ) для решения возникающей в процессе преобразования интеграла (18) задачи Коши для дифференциального уравнения первого порядка:

На выбор многошагового метода накладываются следующие условия: метод устойчивым. Функция должна быть ограничена в правой полуплоскости относительно прямой c i, c i, то есть:

Если функция f (s) аналитична и ограничена в области arg( s c), где, критерий устойчивости может быть ослаблен до A( ) устойчивости.

К соответствующим примерам многошаговых методов относятся методы дифференцирования назад порядка p 6 : для A-устойчивого метода дифференцирования назад второго порядка ( 90 ) можем записать:

При условии того, что функция f (s) в уравнении (20) вычисляется с некоторой порядка ( ).

1.2.4. Модификация шаговой схемы с переменным шагом интегрирования по аргументу Модификация формулы (20) для вычисления n с переменным шагом и линейной аппроксимацией функции выглядит следующим образом:

Для случаев, когда f s / ein – сильно осциллирующая функция, в сочетании с (22) целесообразно использовать комбинированную формулу, учитывающую специфику интегрирования таких функций [10]: с линейной аппроксимацией функции f s s :

или с квадратичной аппроксимацией функции f s s :

Формулы (22) – (24) позволяют перераспределять расчётные узлы по промежутку изменения для получения большей точности результата при сохранении затрат. Общая схема на основе таких формул называется схемой «с ключом», где w2 или w3 – ключ.

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

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

Метод Рунге-Кутты должен быть A-устойчивым;

устойчивости, [1] (1,...,1)T ;

Если принять bT A1 (0,...,0,1), то метод автоматически L-устойчив.

Применяя метод Рунге-Кутты вместо линейного многошагового метода для решения задачи Коши (21), получим с учётом [82]:

Формулы (22) – (24) сохраняют свой формальный вид и в случае схемы обращения преобразования Лапласа на основе метода Рунге-Кутты, но используют соответствующую матричную характеристическую функцию.

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

сформулированным условиям, выберем схемы Радо (Radau IIA) и Лобатто (LobattoIIIC) [80].

1.3. Численно-аналитические результаты 1.3.1. Задача о действии осевой силы на упругий стержень Рассматривается решение задачи о действии осевой силы на упругий стержень (рис. 2).

Изображение и оригинал функции перемещения представляются следующими аналитическими формулами [188]:

Отклики перемещений, вызванные силой F= 1 Н / м 2, наблюдаются в точке y 3м стержня длиной l 3м. Численные результаты получены для материала с параметрами: Е = 2.11 1011 H / м2 ; = 0 ; = 7850кг / м3 ( с = 5184.5 м / с ). Выберем отрезок времени 0.01 с, что составляет чуть более 4 периодов функции перемещения по времени.

Для обращения преобразования Лапласа использован коэффициент R=0.997. Общее число узлов по аргументу (учитывая двукратное использование каждого i для двухэтапной схемы Рунге-Кутты) обозначим через, число узлов по времени – через N.

На рис. 3 а, б приведён вид соответственно действительной ( ) и мнимой ( ) частей спектра перемещений при N=1200.

Результаты обращения при N=/2=500 и N=/2=1000 представлены на рис. 4, соответственно. Решения получены на основе формул шагового метода обращения преобразования Лапласа на узлах методов Эйлера, Лобатто и Радо. Для сравнения приведено решение, рассчитанное по аналитической формуле. На рис. 6 более крупным планом представлен вид двух последних пиков.

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

На рис. 7, 8 представлено сравнение модификаций шагового метода на основе схем семейства Рунге-Кутты с модификацией, использующей переменный шаг интегрирования по аргументу ; для последней рассматривается кусочно-равномерная сетка на промежутках 0, / 2, / 2, 3 / 2, 3 / 2, 2. Число узлов на промежутках выбирается исходя из их процентного распределения 47.5, 5, 47.5, и сохранения общего количества узлов при расчётах по разным схемам. Рис. 7 соответствует N=/2=500, рис.

8 – N=/2=1000. Рис. 8 изображает вид двух последних пиков решения крупным планом.

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

Перераспределение узлов i по расчётному промежутку вышеописанным способом для традиционного подхода не позволяет намного превзойти результаты, полученные с помощью подхода на основе схемы Радо, лишь ускоряя сходимость с измельчением сеток. На малом числе узлов подход на основе методов семейства РунгеКутты отстаёт в точности аппроксимации от подхода с переменным шагом; с дальнейшим измельчением сеток схемы на основе методов Рунге-Кутты предоставляют лучшие результаты. Также при N=/2=2000 график, полученный с помощью схемы с переменным шагом, проявляет характерные осцилляции [10] (рис. 10а), поэтому необходимо использовать модификацию метода на основе формул интегрирования сильно осциллирующих функций (рис. 10б), в то время как остальные схемы использованы без изменений.

На рис. 11, 12 представлено сравнение результатов, полученных с помощью подходов на основе методов семейства Рунге-Кутты и на основе метода Эйлера при использовании переменного шага интегрирования по аргументу. Рассматривается кусочно-равномерная сетка, построенная с теми же условиями, что и ранее. Рис. соответствует N=/2=500, рис. 12 – N=/2=1000. Рис. 13 изображает вид двух последних пиков решения крупным планом.

схема Эйлера, переменный шаг схема Лобатто, переменный шаг схема Эйлера, переменный шаг схема Лобатто, переменный шаг схема Эйлера, переменный шаг схема Лобатто, переменный шаг схема Радо, переменный шаг аналитическое решение Сочетание схем Рунге-Кутты с переменным шагом уже при N=/2=500 допускает некоторые осцилляции, поэтому для расчётов использована модификация на основе формул интегрирования сильно осциллирующих функций. С учётом этого можно утверждать, что шаговый метод на основе схемы Радо на любой из рассмотренных сеток представляет лучший результат при приблизительно равных вычислительных затратах.

Традиционный шаговый метод с переменным шагом обнаруживает появление осцилляций только при измельчении сетки до N=/2=2000. На рис. 14 показаны результаты, полученные на этой сетке с его помощью, а также методом на основе схемы Радо с переменным шагом без учёта и с учётом формул интегрирования сильно осциллирующих функций.

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

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

показано, что с помощью формулы переменного шага интегрирования по аргументу «с ключом» можно снизить вычислительные затраты за счёт уменьшения количества расчётных узлов L. Схема на основе метода Радо позволяет сверх того добиться сходной точности приближения решения на меньшем числе шагов по времени N. На рис. показаны отклики перемещений, полученные традиционным шаговым методом при N=L=2000 и комбинированной формулой с кусочно-равномерным шагом на интервалах 0; / 3, / 3; 5 / 3, 5 / 3; 2 при N=2000 и L 1 =333, L 2 =41, L 3 =333 (решение из [10]) – кривые совпадают. Также приведено решение шаговым методом на основе схемы Радо при N=300, 1 =330, 2 =40, 3 =330.

решение из [10], N=2000, L=707 схема Радо, N=300, = При 1 =330, 2 =40, 3 =330, N=800 подход на основе схемы Радо уже выгодно отличается от подхода на основе метода Эйлера меньшим запаздыванием решения относительно аналитического: см. рис. 16.

1.3.2. Задача о действии осевой силы на пороупругий стержень Рассмотрим решение задачи о действии осевой силы на пороупругий стержень.

аналитическими формулами [188]:

где и, следуя [112], Оригиналы могут быть вычислены следующим образом:

где Отклики перемещения и порового давления, вызванные силой F= 1 Н / м 2, наблюдаются соответственно в точках y 3м и y 0 стержня длиной l 3м. Численные результаты получены для материала с параметрами: K = 4.8 109 Н / м2, G = 7.2 109 H / м2, = 1.9 1010 м4 /( Н с). Рассмотрим отрезок времени 0.02 с, что составляет чуть более периодов функций по времени. Для обращения преобразования Лапласа использован коэффициент R=0.997.

На рис. 17, 18 приведён вид действительной и мнимой частей спектра соответственно перемещений и поровых давлений при N=1200.

Рис. 17а Рис. 17б Рис. 18а Рис. 18б Рассмотрим кусочно-равномерную сетку на промежутках 0, / 2, / 2, 3 / 2, 3 / 2, 2 для двух модификаций шагового метода с использованием переменного шага по аргументу (построенных соответственно на основе метода Эйлера и схемы Радо), и равномерную сетку с сохранением общего числа узлов для модификации на узлах схемы Радо. Для схем с переменным шагом распределение узлов по промежуткам задаётся процентным соотношением 140/3, 20/3, 140/3. Результаты обращения для отклика перемещений при N=/2=300, N=/2=600, N=/2=1200 представлены на рис. 19, 20, соответственно. Для сравнения приведено аналитическое решение. При расчетах с применением схемы Радо и переменного шага интегрирования использована модификация метода на основе формул интегрирования сильно осциллирующих функций.

схема Эйлера, переменный шаг схема Радо, постоянный шаг схема Эйлера, переменный шаг схема Радо, постоянный шаг схема Эйлера, переменный шаг схема Радо, постоянный шаг На рис. 22-24 крупным планом показан вид двух первых и двух последних пиков графика для N=/2=300, N=/2=600, N=/2=1200 соответственно.

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

На рис. 25, 26 представлены результаты, полученные с помощью схемы Радо для отклика перемещения соответственно с постоянным и переменным шагами интегрирования, при увеличении числа узлов.

Результаты обращения для отклика давления при N=/2=600, N=/2=1200, N=/2=2400 представлены на рис. 27-29 соответственно (применены шаговые методы на основе метода Эйлера и схемы Радо с постоянным и переменным шагом интегрирования).

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

На рис. 30, 31 представлены результаты, полученные для отклика давления с помощью схемы Радо соответственно с постоянным и переменным шагом интегрирования, при увеличении числа расчётных узлов: N=/2=300, N=/2=600, N=/2=1200.

Нужно отметить, что качественный характер расхождения результата, предоставленного шаговым методом, с аналитическим решением может варьироваться в зависимости от параметров расчётного материала. На рис. 32, 33 показано соответственно сравнение откликов перемещения и давления, полученных методом на узлах схемы Радо с постоянным шагом интегрирования при N=/2=600, для = 2458 кг / м3 и = 2181 кг / м3.

Чёрными линиями приведены соответствующие аналитические решения.

Здесь в случае разных расчётных плотностей имеем на одной и той же сетке преобладающее занижение (голубая кривая) и превышение (синяя кривая) отклика давления, и занижение амплитуды (голубая кривая), отсутствующее в другом случае (синяя кривая), для отклика перемещения. Этот эффект не является характерной чертой шагового метода на узлах схем Рунге-Кутты, а заметен также и для традиционного шагового метода: см. рис. 34, 35. Численные результаты получены при N=L=600 для Несмотря на различия, во всех случаях расхождение численных результатов обращения с аналитическим решением нивелируется с измельчением сеток.

1.3.3. Моделирование медленной продольной волны в одномерном случае Рассмотрим эффект появления медленной продольной волны в пороупругой среде на примере одномерного пороупругого стержня (рис. 36) длиной l = 9 м с параметрами f =1000 кг / м3, K f = 3.3 109 Н / м2 под действием осевой силы F 1H / м 2. Для расчётов используется численное обращение преобразования Лапласа шаговым по времени методом на узлах схемы Радо семейства Рунге-Кутты. Применяется шаговая схема с ключом, порождаемая квадратурными формулами для сильно осциллирующих функций; интегрирование по аргументу производится с переменным шагом: узлов распределены по промежуткам 0, / 2, / 2, в отношении 95% / 5%.

На рис. 37-44 представлены численно-аналитические решения для откликов порового давления и потока в точке A ( l1 = 1.5 м ). На рис. 37, 38 приведены результаты исследования сходимости шаговой схемы на узлах метода Радо при расчёте соответственно порового давления и потока в точке A для коэффициента проницаемости = 1.9 10 7 м4 /( Н с). На рис. 39-41 показано сравнение откликов давления, а на рис. 42откликов потока, полученных с помощью рассматриваемой шаговой схемы, с результатами [56] на основе шаговой схемы на узлах метода Эйлера при значениях коэффициента проницаемости материала = 1.9 10 6 м4 /( Н с), = 1.9 10 7 м4 /( Н с) и = 1.9 10 10 м4 /( Н с) соответственно.

решение из [56], =N=2000 аналитическое решение =N=2000 решение из [56], =N= =N=2000 решение из [56], =N= решение из [56], =N=2000 аналитическое решение продемонстрирован как на откликах порового давления, так и на откликах порового потока. В работах [46, 47] аналогичные исследования были проведены на основе метода Дурбина. В работе [56] представлены более точные по сравнению с [46, 47] результаты благодаря использованию шагового метода численного обращения преобразования Лапласа, но тем не менее не позволяющие обнаружить детализированный ступенчатый характер изменения отклика порового потока (см. рис. 38, 42, 43, отрезок времени 0В качестве основных результатов главы I можно отметить следующее. Глава посвящена описанию математических моделей. Записаны системы уравнений теории упругости и пороупругости, постановки краевых задач в преобразованиях по Лапласу. Для описания пороупругой среды выбрана полная модель Био. Представлена методика численного обращения преобразования Лапласа на основе шагового метода, как традиционного, так и на узлах семейства схем Рунге-Кутты, приведены модификации с переменным шагом интегрирования и на основе применения квадратур сильно осциллирующих функций. Рассмотрены численно-аналитические решения задач о действии осевой силы на упругий и пороупругий стержни; на примере задачи о пороупругом стержне продемонстрирован эффект возбуждения медленной продольной волны в пористой среде. Проведено сравнение результатов шаговых схем на узлах методов Эйлера и Рунге-Кутты (Лобатто, Радо), в т.ч. с результатами других авторов. В качестве далее используемого метода семейства Рунге-Кутты выбран метод Радо.

Глава II ГИУ и гранично-элементная методика 2.1. Граничные интегральные уравнения Решение краевых задач динамики трёхмерных упругих / пороупругих тел сводится к решению граничных интегральных уравнений. ГИУ для задач динамической следующий вид [11]:

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

2.2. Гранично-элементная дискретизация Граничные интегральные уравнения решаются численно на базе метода граничных элементов. Для дискретизации ГИУ используется метод коллокации.

Методика численного решения ГИУ строится на основе регуляризованных граничных интегральных уравнений в преобразованиях по Лапласу вида [11]:

особенностями, x – точка коллокации.

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

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

где (k, l ) – глобальный номер узла, l – локальный номер узла в элементе k, N l ( ) – функции формы, в качестве которых выбраны следующие интерполяционные полиномы:

В рассматриваемых условиях естественный базис (a1, a2 ), метрический тензор g и единичная нормаль n на элементе записываются следующим образом:

Узлы интерполяции неизвестных граничных функций являются подмножеством геометрических узлов наложенной гранично-элементной сетки. Локальная аппроксимация строится по согласованной интерполяционной модели Р.В. Гольдштейна [28]. Согласно модели Гольдштейна обобщённые граничные перемещения аппроксимируются по билинейным граничным элементам, в то время как обобщённые усилия представляются на элементах постоянными. Для метода коллокации выберем множество узлов, совпадающее с множеством узлов аппроксимации исходных граничных функций. В итоге для каждой контактирующей подобласти получим систему линейных алгебраических уравнений вида:

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

где s – параметр преобразования Лапласа.

Ключевым аспектом расчёта коэффициентов A, B, D является разная модель поведения фундаментальных и сингулярных решений по координатам в рассматриваемом элементе: возможны особенности типа O(1/ r ) и O(1/ r 2 ), в этом случае соответственно используется алгоритм устранения или понижения порядка особенности [11]. После квадратурной формулы Гаусса:

где оценка погрешности, n – порядок квадратурной формулы – определяется исходя из неравенства учитывающего выражение для оценки погрешности квадратуры [4] и в допущении малых изменений якобиана J ( ) в пределах граничного элемента, где * – наперёд заданная погрешность. Если неравенство удовлетворяется на элементе, выбирается минимальный допустимый порядок интегрирования; в противном случае (возможен, т.к. на порядок рассматриваемый элемент рекурсивно переразбивается на четыре новых элемента до тех пор, пока на каждом из элементов неравенство не будет выполнено (см. рис. 46). Такой процесс разбиения конечен, т.к. производная степени 2n от ( x y ) 2 при фиксированных x и граничном элементе ограничена, а размеры подэлементов при каждом новом разбиении стремятся к нулю.

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

компьютерном моделировании строятся непосредственно матрица [K] (полностью заполненная, несимметричная) и вектор правой части {Y} разрешающей системы, минуя поэтапное получение дискретных аналогов ГИУ.

2.3. Программная реализация На основе представленной схемы метода гранично-временных элементов на узлах методов Рунге-Кутты был создан программный комплекс для расчёта задач трёхмерной динамической теории упругости и пороупругости. В него входят два основных исполняемых приложения LBEM и CQM_for_BEM, написанных на языке Фортран.

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

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

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

На рис. 47 представлена схема организации взаимодействия между основными модулями программы. Стрелками обозначены направления потоков данных модулей.

Модуль SYST обеспечивает доступ к банку данных и конфигурационному файлу программы. Модуль FREQ позволяет осуществлять обмен данными с управляющим пакетом CQM_for_BEM (или другим) и передавать необходимые конфигурационные данные в модуль INTG. Обработка параметров материалов: считывание параметров из конфигурационного файла и их приведение к безразмерным величинам производится в модуле GMAT. Обезразмеренные константы материалов передаются в модуль INTG.

Параметры обезразмеривания считываются модулем SCAL, который состоит из подпрограмм, используемых для обезразмеривания величин другими модулями. Данные и подпрограммы модуля SCAL используются также модулями BCON и GMAT. Модуль BCON отвечает за обработку граничных условий задачи и передает информацию о неизвестных величинах в модуль NUMB для вычисления количества неизвестных и их нумерации. Если решаемая задача имеет физическую симметрию, то в конфигурационном файле могут задаваться плоскости симметрии задачи, которые обрабатываются модулем REFL. Данные обработки используются модулем GEON для построения полной геометрии задачи. Данные REFL используются также в модуле INTG. Этот модуль замыкает цепочку, ведущую к получению решения задачи. Его основной задачей является вычисление коэффициентов дискретного аналога и решение соответствующей системы линейных алгебраических уравнений. Результатом вычислений становится набор значений изображений искомых функций по Лапласу для каждого узла/элемента гранично-элементной сетки.

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

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

Программа CQM_for_BEM. Программа CQM_for_BEM является исполняемым файлом и для корректной работы требует файл конфигурации, передаваемый в качестве значения одного из параметров командной строки (–config). Конфигурационный файл содержит информацию, необходимую для управления пакетом LBEM, такую, как данные об узлах интегрирования по аргументу параметра преобразования Лапласа и о временной дискретизации, а также параметры для получения оригинала искомых функций как завершающего этапа вычислений. Корректный проход вычислений с помощью программного комплекса подразумевает двукратный запуск пакета CQM_for_BEM:

первый раз в режиме управления пакетом LBEM посредством заполнения необходимых данных об аргументах изображения искомых функций по Лапласу в зарезервированный файл (это достигается путём передачи в качестве аргумента командной строки параметра –fill); второй раз в режиме прохода шаговой по времени схемы численного обращения преобразования Лапласа, результатом которого является получение оригиналов искомых функций (режим работы по умолчанию). Предусмотрены два варианта формирования итоговых результатов: в виде отдельных файлов для каждой компоненты искомых функций и каждого узла/элемента гранично-элементной сетки (используется для построения графиков динамических откликов), либо в виде одного файла, содержащего все итоговые данные (используется программой визуализации расчётного объекта). Также предусмотрена возможность вывода отдельных данных о спектрах искомых функций, что позволяет, не заканчивая решение задачи, скорректировать параметры схемы или исходные данные для дальнейшего получения точных результатов. Управление режимом вывода спектров функций достигается с помощью параметра командной строки –spectre, который может принимать значения no, function, method или all. Аналогично возможно отключить расчёт оригиналов функций, если в данный момент в качестве результата расчётов интересен только спектр: это достигается параметром –results, который может принимать значения no или yes.

Эксплуатация программного комплекса. Алгоритм решения задачи с помощью программного комплекса на основе метода гранично-временных элементов сводится к трём шагам:

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

Запуск пакета LBEM. Предварительно необходимо заполнить соответствующий конфигурационный файл.

Запуск пакета CQM_for_BEM в режиме шаговой схемы с уже имеющимся конфигурационным файлом.

Предусмотрено несколько наиболее часто используемых сценариев запуска общего комплекса расчётов, выполненных в виде batch-скриптов: один общий проход схемы на основе двух конфигурационных файлов; только непосредственный расчёт оригиналов искомых функций на основе заранее решённой в изображениях по Лапласу задачи и конфигурационного файла для пакета CQM_for_BEM; последовательный итерационный проход схемы на основе массива конфигурационных файлов для LBEM и массива соответствующих конфигурационных файлов для CQM_for_BEM с копированием результатов, полученных на каждой итерации, в назначенные папки.

Также имеются специализированные средства обработки итоговых данных, такие, как несколько вариантов графопостроителя и программа визуализации view-design, позволяющая на основе входных данных об узлах сетки, обходах элементов и откликах обобщённых перемещений/напряжений во времени построить как статическую (в некоторый момент времени), так и динамическую визуализацию расчётного объекта.

Предусмотрены различные режимы визуализации: трёхмерный, полутоновой и комбинированный.

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

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

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

debug.bem //имя отладочного лог-файла 13 //число модулей в конфигурационном файле BGN SYST //начало блока SYST, системные входные параметры 28886080 //используемый стэк динамической памяти в блоках по 4 байта 90000000 //общий размер памяти для хранения данных в блоках по 4 байта END //конец блока BGN MISC //вспомогательные параметры 1 //число задаваемых параметров MEPA R 1.000000000e-07 //порядок формулы интегрирования определяется наперёд заданной погрешностью END BGN KERN //параметры нагрузки 1 //флаг обращения нагрузки 1 //закон нагружения 1.0000000E+00 //коэффициенты нагружения по соответствующему закону 2.0000000E+00 //коэффициенты нагружения по соответствующему закону END BGN SCAL //параметры обезразмеривания 001 2.110000000e+11 3.000000000e+00 1.000000000e-03 //напряжения, расстояния, время END BGN GMAT //характеристики материалов NMAT 002 //число подобластей Material 1 001 1 2.110000000e+11 0.000000000e+00 7850.000000e+ 0.000000000e+00 1.000000000e+00 //характеристики материала первой подобласти Material 2 001 1 2.110000000e+11 0.000000000e+00 7850.000000e+ 0.000000000e+00 1.000000000e+00 //характеристики материала второй подобласти END BGN REFL //плоскости физической симметрии NPLSYM = 0 //число плоскостей симметрии //первая матрица симметрии -1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+ //вторая матрица симметрии 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 -1.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+ END BGN GEON //геометрия гранично-элементной сетки NKRG 2 //количество подобластей isrg 1 //номер подобласти alfad 1.000000000e+00 //признак конечности nel 6 //начальное количество граничных элементов //обходы граничных элементов //координаты граничных узлов 1 +0.50000000e+00 +0.15000000e+01 -0.50000000e+ 2 -0.50000000e+00 +0.15000000e+01 -0.50000000e+ 3 -0.50000000e+00 +0.15000000e+01 +0.50000000e+ 4 +0.50000000e+00 +0.15000000e+01 +0.50000000e+ 5 +0.50000000e+00 +0.00000000e+00 -0.50000000e+ 6 -0.50000000e+00 +0.00000000e+00 -0.50000000e+ 7 -0.50000000e+00 +0.00000000e+00 +0.50000000e+ 8 +0.50000000e+00 +0.00000000e+00 +0.50000000e+ 1 //признак дополнительного разбиения элементов ISETKA 1 //матрица дополнительных разбиений по рёбрам элементов isrg 2 //вторая подобласть alfad 1.000000000e+ 1 +0.50000000e+00 +0.30000000e+01 -0.50000000e+ 2 -0.50000000e+00 +0.30000000e+01 -0.50000000e+ 3 -0.50000000e+00 +0.30000000e+01 +0.50000000e+ 4 +0.50000000e+00 +0.30000000e+01 +0.50000000e+ 5 +0.50000000e+00 +0.15000000e+01 -0.50000000e+ 6 -0.50000000e+00 +0.15000000e+01 -0.50000000e+ 7 -0.50000000e+00 +0.15000000e+01 +0.50000000e+ 8 +0.50000000e+00 +0.15000000e+01 +0.50000000e+ END BGN GAMR // CAMR // END BGN BCON //граничные условия ICNT 0 //признак контакта твёрдого тела с жидкостью LAMR //список кодов контакта по элементам IBCP //список признаков контакта в виде граничных усилий по элементам IBCU //список признаков контакта в виде граничных перемещений по узлам P //список значений граничных усилий 1 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 288 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ U //список значений граничных перемещений 1 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 290 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ LAMR IBCP IBCU 1 0.00000000e+00 0.00000000e+00 -1.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 37 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 288 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 1 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 290 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+ END BGN FREQ //информация о параметрах преобразования Лапласа CQM //признак использования метода квадратур свёрток либо шагового метода численного обращения преобразования Лапласа 4 //признак передачи управления стороннему модулю END BGN NUMB //обработка неизвестных значений граничных функций END BGN ANAL //аналитическое решение 0 //признак использования модуля; отключён 9.0000000e+00 //параметры расчётного тела 7.5000000e+ END BGN INTG //численное интегрирование и решение системы алгебраических уравнений 0 //вариант наименования выходного файла; дата и время отключены u_1 //название выходного файла с изображениями обобщённых перемещений p_1 //название выходного файла с изображениями обобщённых напряжений END Конфигурационный файл программы CQM_for_BEM содержит набор параметров используемой шаговой схемы для управления программой LBEM. В таблице представлен пример входного файла для задачи о действии силы на торец составного упругого призматического тела. Символами // отделены комментарии о значении параметров.

02_u_1 //имя файла для извлечения изображений обобщённых перемещений/напряжений //параметры временной сетки 1.e-2 500 //расчётное время; число шагов по времени //параметры сетки по аргументу параметра преобразования Лапласа 3 //число промежутков равномерного разбиения 0. 0.5 5.783185307179586 6.283185307179586 //границы промежутков 238 24 238 //количество разбиений на промежутках dontmergeresults //признак вывода результатов для различных компонент/узлов в разные файлы; флаг mergeresults включает режим объединения всех результатов в одном файле //параметры пространственной сетки 290 1 //количество расчётных узлов; количество узлов для прогона шаговой схемы 74 //список номеров узлов, для которых нужно получить оригиналы граничных функций 3 1 //количество компонент обобщённых граничных функций; количество компонент для прогона шаговой схемы 2 //список номеров компонент граничной функции, для которых нужно получить оригинал 1.e-3 3. 3. 3.//коэффициенты обезразмеривания по времени и всем компонентам обобщённых граничных функций radau //выбор метода, на основе которого строится шаговая схема; тж euler, lobatto и т.д.

3. //значение ключа в квадратуре сильно осциллирующих функций;

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

тж quadratic и т.д.

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

При включённом режиме симметрии по аргументу параметра преобразования Лапласа расчёты ведутся только на промежутке [0,], поэтому достаточно задать разбиение, покрывающее только этот промежуток. Программа автоматически выберет значения, достаточные для итогового вычисления спектра изображений на всей комплексной окружности, и впоследствии вычислит его путём симметричного отображения относительно. Режим симметрии позволяет экономить приблизительно половину расчётного времени без потери точности результата.

2.4. Задача о действии торцевой силы на упругое призматическое тело Рассматривается задача о действии торцевой силы на однородное упругое призматическое тело с жёстко закреплённым концом (рис. 49). Задача имеет наблюдаются в точке A консоли длиной l 3м. Для расчётов взяты следующие параметры материала: Е = 2.11 1011 H / м2 ; = 0 ; = 7850кг / м3 ( с = 5184.5 м / с ).

Граничные условия (боковая поверхность):

известные:

неизвестные: u1, u2 u3.

Результаты вычислений при N=/2=500, N=/2=1000 представлены на рис. 50, соответственно, где N – количество шагов по времени, – общее количество узлов по аргументу (с учётом двукратного использования узлов). Использованы модификации шагового метода обращения преобразования Лапласа на основе метода Эйлера с переменным шагом интегрирования по ; на основе схемы Радо с постоянным шагом и распределением узлов 47.5, 5, 47.5. Здесь и далее переменный шаг используется совместно с модификацией схемы на основе квадратур сильно осциллирующих функций. Граничноэлементная сетка содержит 224 квадратных элемента. Для сравнения приведён отклик перемещений, рассчитанный по аналитической формуле.

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

На сетке N=/2=500 отчётливо виден выигрыш от использования переменного шага интегрирования. Однако на сетке N=/2=1000 решения, полученные с помощью схемы Радо, начиная примерно с середины рассмотренного промежутка времени, искажены нарастающими осцилляциями. Это объясняется численной ошибкой, внесённой при расчётах значений обращаемой функции и обусловленной достаточно большим шагом пространственной сетки и достаточно малым – временной. Подтвердим этот факт.

На рис. 52 представлено сравнение решения, полученного на основе метода Эйлера, с аналогичными результатами (/2=1000) при меньшем числе шагов по времени для схемы на основе метода Радо: N=500 и N=800.

схема Радо, переменный шаг, N=800 схема Радо, переменный шаг, N=500, схема Эйлера, переменный шаг, N=1000 аналитическое решение Видно, что схема на основе метода Радо позволяет добиться гладкости графика за счёт увеличения шага по времени, при этом не проигрывая в точности приближения решения традиционному шаговому методу.

Численное обращение преобразования Лапласа шаговым методом на узлах схемы Рунге-Кутты можно рассматривать как обращение преобразования традиционным шаговым методом, применённое к некой линейной комбинации спектров f s s, порождённых целочисленными и дробными узлами схемы Рунге-Кутты. Так как эта комбинация определяется выбранным методом Рунге-Кутты, на основе которого строится схема обращения, для краткости будем называть её спектром, порождённым соответствующим методом. На рис. 53, 54 для сетки N=/2=1000 приведены соответственно действительная и мнимая части спектра, порождённого схемой Радо.

Видно, что они значительно искажены на участках 1 / 2, / 2 по сравнению со спектром, построенным с использованием аналитической формулы для f s на той же сетке (см. рис. 55, 56).

Следует заметить, что в спектре, полученном на основе метода Эйлера в рамках гранично-элементной схемы, появляются искажения, не влияющие на результат обращения преобразования Лапласа: см. рис. 57, 58.

Сравним спектры функции f s, порождённые целочисленными и дробными узлами схемы Радо в рамках гранично-элементной схемы и совместно с расчётами f s по аналитической формуле: рис. 59-66. На рис. 59-60 представлены действительные части, на рис. 61-62 мнимые части спектров на целочисленных узлах, полученные в рамках ГЭсхемы и аналитически, соответственно. Аналогично на рис. 63-64 представлены действительные и на рис. 65-66 мнимые части спектров на дробных узлах, полученные в рамках ГЭ-схемы и аналитически, соответственно.

Рис. Рис. Рис. Рис. Рис. 63а Рис. 63б Рис. Рис. Таким образом, можно утверждать, что возникновение осцилляций на рис. обусловлено численной ошибкой при расчёте f s в наборе аргументов s, порождённом дробными узлами схемы Радо. Численная ошибка расчёта образа функции в наборе аргументов, порождённом целочисленными узлами, слабо влияет на характер спектра. Это совпадает результатами, полученными при помощи схемы на основе метода Эйлера (т.е. с рассмотрением только целочисленных узлов).

Обезразмеривание исходных данных (в качестве характерных величин взяты Е* = 2.11 1011 H / м2, l* = 3 м, t* = 6 103 c ) позволяет ослабить влияние численной ошибки на вид отклика перемещений. На рис. 67 приведены решения без учёта и с учётом обезразмеривания при вычислениях.

На рис. 68 представлены графики отклика перемещений, полученные с помощью схемы Радо при N=/2=1000 на кусочно-равномерной сетке на промежутках 0, / 2, одинаковым процентным распределением узлов по промежуткам 47.5, 5, 47.5. Синяя кривая соответствует равномерному распределению /2=2000 узлов по общему промежутку 0, 2.

схема Радо, переменный шаг, без учёта обезразмеривания схема Радо, переменный шаг, обезразмеривание схема Радо, переменный шаг, N=/2=1000, 0, / 2, / 2, 3 / 2, 3 / 2, схема Радо, переменный шаг, N=/2=1000, 0, 1 / 2, 1 / 2, 2 1 / 2, схема Радо, постоянный шаг, N=/4=1000 аналитическое решение Равномерное распределение, достигнутое без уменьшения количества узлов на улучшить результат и избежать паразитных осцилляций при наличии некой погрешности расчёта f s, в отличие от случая численного обращения преобразования Лапласа с применением аналитической формулы f s. Однако другое перераспределение узлов i по расчётному промежутку способно снизить влияние численной ошибки и получить лучший результат.

На рис. 69, 70 приведены соответственно действительная и мнимая части спектра, порождённого схемой Радо при N=/2=1000 на кусочно-равномерной сетке на На рис. 71 представлен отклик перемещения, полученный схемой на основе метода Эйлера при N=1000 и количестве узлов i соответственно 950, 100, 950 на промежутках 0, / 2, / 2, 3 / 2, 3 / 2, 2, и схемой на основе метода Радо при N=700 и том же распределении i и N=1000 и количестве i 996, 8, 996 на промежутках Как видно из рис. 68, 71, перераспределения узлов оказывается достаточно для подавления нежелательных осцилляций решения.

На рис. 72 представлен отклик перемещения, полученный при N=2000. На временной сетке схема на основе метода Эйлера испытывает ощутимое влияние численной ошибки (красная кривая. Приведены результаты для схемы на основе метода Эйлера при количестве узлов i соответственно 1900, 200, 1900 на промежутках 0, / 2, / 2, 3 / 2, 3 / 2, 2, количестве узлов i 1996, 8, 1996 на 0, 0.15, 0.15, 2 0.15, 2 0.15, 2, а также схемы на основе метода Радо при количестве i 1996, 8, 1996 на 0, 0.15, 0.15, 2 0.15, 2 0.15, 2.

Дальнейшее измельчение сетки по времени не даёт ощутимого улучшения результата (см. рис. 73). Разница между рассматриваемыми схемами решения (рис. 72, 73) за счёт использования дробных узлов практически не приносит рационального вклада в результат, а подавляется ошибкой вычислений. На рис. 73 приведены решения, полученные основе метода Радо при N=1000 и количестве i 996, 8, 996 на 0, 1 / 2, 1/ 2, 2 1/ 2, 2 1/ 2, 2, и на основе метода Эйлера и метода Радо при N=2000 и практически неразличимы.

На рис. 74, 75 показаны решения, полученные на основе метода Радо при N=/2=1000 и процентном распределении узлов i 47.5, 5, 47.5 по промежуткам 0, / 2, / 2, 3 / 2, 3 / 2, 2 на разных пространственных сетках: из 224 (сетка 1), 504 (сетка 2) и 896 (сетка 3) элементов (см. рис. 76). Рис. 75 более крупным планом изображает вид первого и третьего пиков решения.

схема Эйлера, N=2000, i : 1996 0, 0.15, 8 0.15, 2 0.15, 1996 2 0.15, схема Радо, N=2000, i : 1996 0, 0.15, 8 0.15, 2 0.15, 1996 2 0.15, сетка 1 сетка 2 сетка 3 аналитическое решение Измельчение гранично-элементной сетки обеспечивает сходимость схемы. На сетке из 896 элементов шаговая схема на основе метода Радо (см. рис. 77, зелёная кривая) сходится быстрее схемы на основе метода Эйлера (красная кривая).

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