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

В научной литературе широко исследованы гибридные методы, сочетающие МКЭ и МГЭ для решения задач акустики, взаимодействия жидкости и конструкций, а также рассеяния волн. В работе [3] представлен итерационный подход для анализа взаимодействия насыщенных пористых и герметичных упругих сред комбинацией МКЭ-МГЭ. Метод позволяет учитывать нелинейности в подобластях МКЭ. В работе [4] разработана явная процедура прямой связи МКЭ-МГЭ для нелинейной динамики. Метод позволяет анализировать подобласти независимо и использовать разные временные шаги, что повышает гибкость и эффективность связанных анализов. В работе [5] представлен метод связи МКЭ-МГЭ для анализа взаимодействия жидкости и конструкции во временной области. Метод использует оператор Пуанкаре-Стеклова и позволяет получать априорные оценки погрешности для модельных задач. В работе [6] проанализирована перекрывающаяся непрерывная модель для задачи распространения волн Гельмгольца в неограниченных областях. Метод позволяет точно сохранять условие излучения и демонстрирует эффективность для гетерогенных сред. В работе [7] разработана стратегия взаимодействия МКЭ-МГЭ для акустического рассеяния на гармониках времени в средах с переменной скоростью звука. Метод использует растворную переменную на границе связи и демонстрирует квази-оптимальную сходимость. В работе [8] предложен гибридный численный метод для решения задач изменения фазы с движущимися границами. Метод сочетает радиальные базисные функции и сеточный подход, что позволяет решать задачи растворения и затвердевания. В работе [9] исследованы апостериорные оценки ошибок для метода виртуальных элементов (МВЭ) и МКЭ на основе МГЭ. Метод позволяет вычислять расстояние между градиентом конечных элементов и его постобработанной версией, а также обнаруживать точки необычайной точности. В работе [10] предложен регуляризатор на основе условий излучения на поверхности (OSRC) для стабилизации связи МКЭ-МГЭ. Метод улучшает сходимость итеративных решателей, особенно на высоких частотах, и снижает вычислительные затраты. В работе [11] разработана структура глубокого обучения для прогнозирования передаточных функций в структурно-акустических моделях. Метод позволяет значительно сократить время вычислений по сравнению с традиционными МКЭ-вычислениями. В работе [12] предложена структура акустического анализа и оптимизации формы для трехмерных двукратно периодических многослойных структур. Метод использует изогеометрический метод граничных элементов (ИГ-МКЭ) и позволяет оптимизировать форму акустических структур. В работе [13] представлен полуаналитический метод конечных элементов для моделирования 2D направленных волн в композитных ламинатах. Метод позволяет анализировать рассеяние волн на расслоениях и получать коэффициенты отражения и пропускания. В работе [14] проведено сравнение высокопорядковых и плосковолновых базисных функций для решения задач Гельмгольца. Показано, что базис NURBS превосходит базис Лагранжа, а также проанализированы ошибки и численная дисперсия для различных базисов. В работе [15] предложен гибридный метод МКЭ-МГЭ для решения акустико-структурных задач. Метод использует дискретный сдвиговый зазор для построения МКЭ без блокировок и коллокацию МГЭ на основе формулировки Бертона-Миллера. В работе [16] разработана симметричная схема связи МКЭ-МГЭ для рассеяния переходных акустических волн на неоднородных анизотропных препятствиях. Метод позволяет дискретизировать прошедшую и рассеянную волны и демонстрирует устойчивость в гильбертовом пространстве.

Значительное число научных работ посвящено разработке численных методов для анализа слоистых, анизотропных и функционально-градиентных сред. В работе [17] предложен модифицированный метод конечных элементов с масштабируемой границей (МКЭМГ) для анализа динамического сопротивления анизотропных слоистых сред. Метод позволяет решать задачи для кросс-анизотропных и общих анизотропных сред, а также упрощает расчеты для кросс-анизотропных слоистых сред. В работе [18] МКЭМГ применен для решения динамической реакции трехмерного бесконечного грунта. Алгоритм непрерывного деления позволяет получать матрицу динамической жесткости слоистого грунта в частотной области, а также решать задачи для заглубленных фундаментов. В работе [19] разработан 3D подход, сочетающий МКЭ и МГЭ в частотной области для анализа динамического отклика конструкций на слоистых трансверсально-изотропных полупространствах. Метод использует функции Грина и позволяет учитывать сингулярное поведение напряжений на границах раздела. В работе [20] предложена новая методология для анализа продольно-инвариантных задач взаимодействия грунт-конструкция в эластодинамике. Метод использует 2.5D связь МКЭ-МГЭ и позволяет повысить вычислительную эффективность за счет использования метода фундаментальных решений для моделирования распространения волн. В работе [21] разработан метод для вычисления 2.5D эластодинамических функций Грина для слоистых полупространств в декартовых координатах. Метод позволяет получать точные решения для однородных и слоистых сред. В работе [22] предложена 2.5D методология МКЭ-МГЭ для анализа динамической реакции сегментированных тоннелей в насыщенном грунте. Метод позволяет учитывать влияние соединений тоннеля на динамические напряжения и давление поровой воды, что недооценивается в однородных моделях. В работе [23] представлен метод конечных элементов на основе границ (МКЭГ) для анализа напряжений в анизотропных функционально-градиентных материалах. Метод использует кусочно-однородные слои и позволяет учитывать эффекты анизотропии и градиентности без построения сетки области.

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

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

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

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

Поведение всех составляющих конструкции и N-слойного основания считается линейно-упругим. При этом свойства среды в ее j-й составляющей описываются плотностью ρj и коэффициентами Ламе λj, μj (j = 0, 1, …, N – 1). Контакт всех сред предполагается жестким, определяющим равенство нормальных и касательных элементов векторов напряжений и перемещений на границах сопряжения.

Движение среды определяется решением системы уравнений Ламе относительно амплитудных функций перемещений точек среды [26]:

(1)

Здесь система координат (x, y, z) связана с поверхностью многослойного полупространства, как показано на рис. 1.

Не ограничивая общности постановки задачи, считаем, что колебания возбуждаются распределенными в ограниченной области Ω усилиями на плоской границе полупространства при x = 0:

При этом линейный параметр α в формуле (1) соответствует характерному размеру этой области. Замыкает постановку задачи требование удовлетворения решений принципу предельного поглощения [27].

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

Вид матриц фундаментальных решений опирается на решение вспомогательной задачи о колебаниях однородного полупространства (x > 0) с упругими параметрами j-го слоя при заданном векторе напряжений X1(y, z) на его поверхности x = 0 в локальной системе координат.

Для этого выразим поле перемещений в преобразованном по Фурье матричном виде:

(2)

Аналогично может быть представлен вектор напряжений:

Вектор C(α, β) определяется из граничного условия на поверхности x = 0 в форме:

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

Для построения матриц фундаментальных решений считаем, что в точке r0 = (x0, y0, z0) j-го упругого слоя 0 < x < hj действует сосредоточенное усилие:

Тогда смещения точки с координатами (x, y, z), удовлетворяющие уравнению (1) с добавлением в левой части слагаемого для объемной силы qδ(r – r0) представляют собой интеграл вида:

(3)

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

(4)

Для построения матриц фундаментальных решений в преобразованном по Фурье виде используется принцип суперпозиции решений трех задач:

(5)

где , соответствующее формуле (2), представляет собой поле излучаемых и отраженных волн от границы слоя x = 0; — поле излучаемых и отраженных волн от границы слоя x = hj; — поле источника колебаний в безграничном пространстве; легко получить из путем перехода к новой системе координат: x → h – x, изменении знаков у ортов ex, ez и введении нового вектора напряжений вместо .

Компоненты матрицы можно найти применением формул переразложения сферических волн по набору плоских в форме двойного интеграла Фурье [28], что приводит к результату:

,

,

,

,

,

,

,

.

Аналогичное представление имеют и компоненты напряжений . В случае, если источник колебаний отсутствует в выбранном j-м слое, третье слагаемое в представлении (5) исключается.

Для окончательного определения фундаментальных решений необходимо исключить вектор-функции для каждой компоненты многослойного полупространства. Это реализуется применением условий непрерывности векторов смещений и напряжений u(j+1) = u(j), t(j+1) = t(j) в форме преобразований Фурье при переходе через границы x = xj, j = 0, 1 ,…, N – 2 раздела j-й и (j + 1)-й компонент слоистой среды, а также условий свободной границы при x = 0. Полученная система линейных алгебраических уравнений и ее определитель определяют дисперсионные свойства многослойного полупространства.

С использованием динамической теоремы взаимности далее получим граничные интегральные уравнения на поверхности S для узловых значений векторов перемещений и напряжений, имеющие следующий вид [29]:

(6)

Здесь up, τk — компоненты векторов напряжений и перемещений соответственно; — компоненты построенных матриц фундаментальных перемещений и напряжений в точке r от действия источника в точке r0; Rk — компоненты вектора внешних усилий на границе полупространства; nk — компоненты вектора внешней нормали к S.

Путем дискретизации системы (6) методом линейных граничных элементов получена СЛАУ относительно неизвестных узловых перемещений uS и сил FS:

(7)

Рассматриваемая задача сводится к совместному решению СЛАУ МГЭ и уравнения движения по МКЭ [30]:

(8)

Здесь M — матрица масс; C — матрица сопротивлений; K — матрица жесткости конструкции; u — амплитудные функции вектора узловых перемещений; F — вектор узловых нагрузок; ω — частота установившихся колебаний. Отметим, что векторы u и F содержат переменные для узлов, лежащих на общей границе S, и являющиеся неизвестными величинами.

После совместного решения СЛАУ (7, 8), выполняется подстановка полученного решения в известные представления волновых полей (2, 3, 4), что позволяет рассчитать динамические характеристики напряженно-деформированного состояния многослойного основания и рассмотреть его взаимодействие со связанным полузаглубленным объектом.

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

Рассматривается несколько случаев постановки задачи — при наличии и отсутствии строительной конструкции, а также при различных конфигурациях жесткости сред основания. Физические свойства слоев основания и элементов конструкции отражены в таблицах 1 и 2 соответственно, где ρ — плотность; E — модуль Юнга; ν — коэффициент Пуассона.

Нагрузка q = 1 кПа приложена на поверхности в точке (x0; y0; z0) = (0; 4,7; 0), один из нижних углов куба расположен в точке (x1; y1; z1) = (0; –1; –1), координаты точки наблюдения (x; y; z) = (0; 1; 0).

Для решения поставленной задачи используется комбинация разработанного авторами программного продукта, реализующего рассмотренный механизм определения характеристик напряженно-деформированного состояния многослойного основания методом граничных элементов и комплекса ANSYS Mechanical для анализа конструкции и окружающего грунта методом конечных элементов. Численный эксперимент проводился на ПК с процессором Intel Core i5-13600k и 32 Гб оперативной памяти с частотой 6000 МТ/c.

Результаты исследования. Получены данные об амплитудно-частотных характеристиках точки наблюдения (рис. 3).

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

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

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

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

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

Рассмотрим реализацию различных методик распараллеливания вычислений в программе на языке Fortran. Используем три технологии: Fortran Coarray, Message-Passing Interface (MPI) и OpenMP.

Coarray — это технология, в рамках которой программа на Fortran реплицируется на нескольких процессорах (образах). Данная технология внедрена в Fortran 2008 для поддержки параллельного программирования. Она предоставляет простой и интуитивно понятный синтаксис для выполнения параллельных вычислений без явной передачи сообщений. В Fortran Coarray — это переменная с одной или несколькими объявленными дополнительными соразмерностями. Каждый образ имеет собственную локальную копию данных Coarray [31].

MPI, или интерфейс передачи сообщений, представляет собой стандартизированную и портативную систему передачи сообщений, предназначенную для облегчения связи и координации между параллельными процессами в среде параллельных вычислений с распределенной памятью. MPI обычно используется в приложениях высокопроизводительных вычислений (HPC) для обеспечения связи между несколькими процессами, работающими на разных узлах кластера. В Fortran MPI предоставляет набор вызовов подпрограмм, которые обеспечивают прямую связь между параллельными процессами [32].

OpenMP (открытая многопроцессорная обработка) — это API, который поддерживает кроссплатформенное параллельное программирование с общей памятью на C, C++ и Fortran. OpenMP упрощает разработку параллельных программ, предоставляя набор директив компилятора и библиотечных процедур для выражения параллелизма в коде [33]. Эта технология особенно хорошо подходит для распараллеливания на уровне циклов и областей. OpenMP контролирует безопасность потоков, включая обработку частных и общих переменных. К недостаткам OpenMP можно отнести необходимость учета потенциальных условий конкуренции за данные и ручной синхронизации.

В исходное программное решение внесены модификации для обеспечения параллельных вычислений. В ходе эксперимента выполнено распределение нагрузки на 2, 4, 6 и 8 потоков (образов), результаты представлены в таблицах 3, 4 и на рис. 5.

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

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

В связи с тем, что решение данной задачи является ресурсоемким процессом и требует значительных временных затрат, применение технологий параллельных вычислений для оптимизации расчетов является естественным этапом развития программного продукта. Анализ затрат времени на вычисление фундаментальных решений с использованием технологий Coarray, MPI и OpenMP показал высокую степень сходства прироста производительности для всех трех методов, что свидетельствует о зрелости данных технологий в рамках инструментария языка Fortran. При переходе от последовательной реализации к двухпоточной наблюдается сокращение времени выполнения расчетов на ~30 %. Увеличение числа потоков на два демонстрирует схожий относительный прирост производительности. Данная тенденция прослеживается и далее. Таким образом, дальнейшее увеличение числа потоков и, соответственно, задействованных ядер процессора, приведет к снижению затрат времени на расчеты, но в абсолютных значениях разница во времени для большого числа потоков будет невелика. Об этом свидетельствует и рост затраченного процессорного времени с увеличением числа потоков для технологий MPI и OpenMP, который связан с увеличением затрат времени на синхронизацию параллельных процессов. Отличия в поведении реализации программного модуля на базе Coarray могут быть обусловлены повышенными затратами времени на инициализацию параллельных образов или, потенциально, связаны с особенностями выполнения программного средства на данном оборудовании. Более ранние исследования показывают, что дальнейшее увеличение числа образов Coarray приведет к схожим затратам процессорного времени относительно других технологий.

Заключение. В соответствии с полученными результатами можно сделать вывод, что использованная технология матричного построения фундаментальных решений для многослойных сред с поверхностными и заглубленными объектами является эффективным инструментом, позволяющим точно учесть дисперсионные свойства многослойной среды, условия излучения энергии на бесконечность при взаимодействии объекта и основания, а также естественным образом распараллелить вычислительные процессы при совместном применении методов граничных интегральных уравнений и конечных элементов. Применение технологий Coarray, MPI и OpenMP для параллельных вычислений позволяет кратно снизить время работы программного модуля, не приводя к значительному увеличению затрат ресурсов. Применение различных технологий дает схожий прирост производительности программного модуля, а различия в затратах процессорного времени связаны с особенностями реализации технологии распараллеливания. Выбор оптимальной технологии зависит от конкретной реализации приложения и ограничений выбранной платформы. Так, в случае компиляции разработанного модуля в качестве DLL-библиотеки и интеграции в комплексное решение необходимо учитывать, что для работы технологии Coarray базовое приложение также должно быть реализовано на Fortran, а инициализация MPI-приложений происходит через утилиту mpiexec. В качестве альтернативы рассмотренным технологиям предлагается использовать OpenCL, OpenACC и CUDA, позволяющих проводить эффективные вычисления с использованием графических ускорителей. Исследование методик применения данных технологий для анализа сложных систем строительных конструкций и многослойных оснований является одним из приоритетных направлений дальнейших исследований.