Введение. Устройства из пьезоматериалов широко используются, давно активно изучаются и совершенствуются. Отдельно следует отметить медицинские ультразвуковые приборы (оборудование для диагностики, ультразвуковые скальпели) [1–4] и мобильные генераторы энергии [5]. Исследование [6] описывает комбинирование фото- и пьезоэлектрического эффектов для создания действенных компактных источников энергии. В науке и промышленности изучаются новые материалы, рассчитанные на эксплуатацию в специфических условиях. В [7] рассматривается создание бессвинцового пьезоактивного состава, подходящего для эксплуатации при различных температурах.

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

Цель исследования — создание численных методов решения задачи определения частот резонанса для системы упругих тел. Достижение заявленной цели требует решения двух задач. Первая: разработать методы дискретизации задачи на основе метода конечных элементов (МКЭ). Вторая: провести программную реализацию выбранного метода на языке С# на платформе .net. Все известные программы учитывают контекст библиотеки классов комплекса ACELAN-COMPOS [8]. При решении обобщенной задачи на собственные значения широко применяются методы, основанные на обращении матриц. Однако они неприменимы к матрицам большой размерности. В представленной научной работе это ограничение преодолевается следующим образом:

Материалы и методы. В первую очередь предлагаемый подход призван решать статические задачи электроупругости при реализации метода осреднения [9], который задействуют для расчета эффективных свойств пьезокомпозитов. В связи с этим на этапе построения глобальных матриц МКЭ представлены только матрицы жесткости. В данном исследовании дополнительно реализовали логику построения матриц масс и разработали программные интерфейсы (application programming interface, API, англ. — интерфейс прикладного программирования) для обмена данными о задачах на собственные значения с модулями пре- и постпроцессинга. Разработанные подпрограммы оценивались с точки зрения производительности и применимости для задач большой размерности. Проводились численные эксперименты с целью валидации созданных алгоритмов для таких задач малой размерности, которые позволяют получить решение общими методами в вычислительном пакете MATLAB. Далее выполнялось тестирование на задачах с большим числом неизвестных и с учетом распараллеливания отдельных операций.

Математическая модель решаемой задачи состоит из определяющих соотношений [9]:

,,(1)

, , (2)

, (3)

Здесь σ — тензор напряжений; ρj — плотность тела; ε — тензор деформаций; u — вектор перемещений; D — вектор электрической индукции; E — вектор напряженности электрического поля; — вектор массовых сил; — электрический потенциал; , , — коэффициенты демпфирования; , , — тензоры упругих констант, пьезомодулей и диэлектрических проницаемостей; индекс j — номер тела в модели.

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

Здесь — матрица функций формы для поля перемещений;

Nϕ — вектор функций формы для электрического потенциала;

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

В таком случае исходная задача (1–3) приобретает вид:

. (4)

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

a = [U, Ф].

В задаче теории электроупругости:

(5)

Матрицы , и — симметричные. В случае гармонических колебаний на собственной частоте можно записать:

, обозначив через соответствующий собственный вектор.

Рассмотрим свободные колебания если . В этом случае задача (4) представляется в виде:

. (6)

Таким образом, исходная задача сводится к обобщенной задаче на собственные значения (6). Для ненулевого vi неравенство (6) решается нахождением матрицы, обратной K. Однако при этом разреженная матрица становится заполненной, то есть метод непригоден для матриц больших размеров. Поэтому нужно использовать другие методы, не требующие нахождения обратной матрицы. Для решения этой задачи в данной работе программно реализован модифицированный метод Ланцоша [10]. Автор этой модификации — Т.С. Мартынова. Описание разработки в данной статье не приводится. Из используемых в методе операций наиболее затратной с точки зрения вычислительных ресурсов оказалось решение системы линейных алгебраических уравнений (СЛАУ), необходимое для выполнения спектральной трансформации.

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

Координатный формат предполагает хранение троек (триплетов) значений (i, j, k), представляющих собой координаты (i, j) и значения (k) ненулевых элементов. CSR иногда называют CRS или Йельским форматом. Он предполагает хранение разреженной матрицы в виде трех массивов. Рассмотрим матрицу размера с ненулевыми элементами. Опишем возможную организацию ее хранения. Все ненулевые элементы нужно разместить в одном массиве размера. Позиции этих элементов в столбцах разместить в другом массиве размера , а третий массив размера использовать для хранения индексов первых элементов строк. Аналогично реализуется хранение в формате CSC.

Формат SKS предполагает хранение ленты матрицы переменной ширины, включающей в себя все ненулевые элементы. В этом случае допускается хранение нулей. Эффективность этого формата зависит от перенумерации строк матрицы. Методы сокращения размера ленты описаны в [11], однако требует отдельного исследования их применимость к матрице жесткости, получаемой при решении трехмерной задачи с использованием МКЭ.

Для решения СЛАУ задействовали итерационный симметричный метод LQ (Symmetric LQ Method, SYMMLQ [12]), адаптированный к перечисленным выше форматам хранения.

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

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

SKS имеет аналогичные проблемы с добавлением ненулевых элементов и сильно зависит от перенумерации неизвестных в задаче. Остановимся на примере квазирегулярной сетки, которая используется в пакете ACELAN-COMPOS для работы с представительными объемами композитов. Ширина ленты, содержащей все ненулевые элементы, может быть определена заранее и зависит от числа узлов и типа конечного элемента. В общем случае произвольной конечноэлементной сетки сложно заранее оценить размер ленты.

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

Итак, сетка представляла собой куб с регулярным разбиением восьмиузловыми конечными элементами. Для иллюстраций использовалась модельная матрица из 500 строк. Матрицы, представленные на рис. 1 а и 1 б, не подвергались дополнительной перенумерации узлов и отличаются только нумерацией степеней свободы. В 1 а:

.

В 1 б неизвестные, отвечающие за распределение потенциала, собраны в конце вектора:

.

Неизвестные в матрицах 1 в и 1 г занумерованы аналогично, но узлы конечноэлементной сетки предварительно перенумерованы согласно их координатам путем поочередной сортировки всех узлов по каждой из координат. Такой прием широко используется для построения более эффективных модулей решения СЛАУ, так как позволяет работать с матрицей в специальном ленточном формате, удобном для распараллеливания. Для комплекса ACELAN-COMPOS реализованы подобные внешние модули [13–15], однако в данной работе использовались только форматы хранения разреженных матриц общего вида.

В таблице 1 сведены данные о времени, необходимом на выполнение базовых операций с матрицами в различных форматах.

Результаты экспериментов показали, что операция преобразования из координатного формата хранения в компактные занимает мало времени. При этом перенумерация узлов сетки для формирования блочно-ленточной матрицы позволяет получить заметный прирост производительности даже без изменения внутренней структуры матрицы в памяти. Формат CSR оказался оптимальным с точки зрения эффективности операции умножения матрицы на вектор. При обращении матрицы формат SKS более эффективен, однако для задач с числом неизвестных порядка 103 заметно быстрее работают итерационные методы решения обобщенной задачи на собственные значения. Далее экспериментально оценили производительность программной реализации метода Ланцоша. Измерили вклад всех операций в общее время решения. В результате установили, что операция решения СЛАУ занимает до 95 % от общего времени работы алгоритма. В ходе работы алгоритма строится подпространство Крылова, и в зависимости от его размерности меняется число СЛАУ, которые необходимо решать. Отметим, что размерность подпространства Крылова выбиралась на основе эвристик относительно числа искомых собственных значений. При этом СЛАУ отличаются только правой частью, благодаря чему требования к выделяемой памяти остаются невысокими. Среди базовых операций, применяемых в ходе решения СЛАУ методом SYMMLQ, наибольшие вычислительные затраты нужны для умножения матрицы на вектор. Для увеличения производительности алгоритма реализовано простейшее распараллеливание с общей памятью. Для формата CRS выделялись блоки строк. Они передавались в отдельные потоки, которые вычисляли соответствующие компоненты результирующего вектора. Прирост производительности составил 40–50 % при использовании 8 потоков. При этом для матриц порядка 103 элементов прирост был около 40 %, а для матриц порядка 104 — около 50 %.

Обсуждение и заключение. В рамках данной работы реализован метод решения обобщенной задачи на собственные значения для матриц, получаемых при моделировании электроупругих тел. Созданы программные модули на языке C# для построения матриц масс методом конечных элементов и выполнения вспомогательных операций в рамках метода Ланцоша (работа с векторами подпространства Крылова, переортогонализация, нахождение собственных векторов). Вычислительная сложность обусловлена в основном операциями умножения разреженных матриц на вектор. В связи с этим проводились численные эксперименты по определению оптимальных форматов хранения матриц, оптимальной структуры матрицы, получаемой в результате перенумерации узлов КЭ-сетки и степеней свободы в узлах. Разработана версия итерационного алгоритма SYMMLQ с использованием параллельных вычислений. Итоговая схема работы включает три пункта. Во-первых, строятся глобальные матрицы в координатном формате с алгоритмом перенумерации (рис. 1 в). Во-вторых, данные преобразуются в формат CRS. В-третьих, данные обрабатываются методом Ланцоша, который включает метод SYMMLQ для решения СЛАУ. Результаты работы включили в программный пакет ACELAN-COMPOS.