Introduction. The widespread use of piezoelectric materials in various industries stimulates the study of their physical characteristics and determines the urgency of such research. In this case, modal analysis makes it possible to determine the operating frequency and the coefficient of electromechanical coupling of piezoelectric elements of various devices. These indicators are of serious theoretical and applied interest. The study was aimed at the development of numerical methods for solving the problem of determining resonance frequencies in a system of elastic bodies. To achieve this goal, we needed new approaches to the discretization of the problem based on the finite element method and the execution of the software implementation of the selected method in C# on the .net platform. Current solutions were created in the context of the ACELAN-COMPOS class library. The known methods of solving the generalized eigenvalue problem based on matrix inversion are not applicable to large-dimensional matrices. To overcome this limitation, the presented scientific work implemented the logic of constructing mass matrices and created software interfaces for exchanging data on eigenvalue problems with pre- and postprocessing modules.

Materials and Methods. A platform was used to implement numerical methods .net and the C# programming language. Validation of the research results was carried out through comparing the values found with solutions obtained in well-known SAE packages (computer-aided engineering). The created routines were evaluated in terms of performance and applicability for large-scale tasks. Numerical experiments were carried out to validate new algorithms in small-dimensional problems that were solved by known methods in MATLAB. Next, the approach was tested on tasks with a large number of unknowns and taking into account the parallelization of individual operations. To avoid finding the inverse matrix, a modified Lanczos method was programmatically implemented. We examined the formats for storing matrices in RAM: triplets, CSR, СSC, Skyline. To solve a system of linear algebraic equations (SLAE), an iterative symmetric LQ method adapted to these storage formats was used.

Results. New calculation modules integrated into the class library of the ACELAN-COMPOS complex were developed. Calculations were carried out to determine the applicability of various formats for storing sparse matrices in RAM and various methods for implementing operations with sparse matrices. The structure of stiffness matrices constructed for the same task, but with different renumbering of nodes of a finite element grid, was graphically visualized. In relation to the problem of the theory of electroelasticity, data on the time required to perform basic operations with stiffness matrices in various storage formats were summarized and presented in the form of a table. It has been established that the renumbering of grid nodes gives a significant increase in performance even without changing the internal structure of the matrix in memory. Taking into account the objectives of the study, the advantages and weaknesses of known matrix storage formats were named. Thus, CSR was optimal when multiplying a matrix by a vector, SKS was optimal when inverting a matrix. In problems with the number of unknowns of the order of 103, iterative methods for solving a generalized eigenvalue problem won in speed. The performance of the software implementation of the Lanczos method was evaluated. The contribution of all operations to the total solution time was measured. It has been found that the operation of solving SLAE takes up to 95% of the total time of the algorithm. When solving the SLAE by symmetric LQ method, the greatest computational costs were needed to multiply the matrix by a vector. To increase the performance of the algorithm, parallelization with shared memory was resorted to. When using eight threads, the performance gain increased by 40–50%.

Discussion and Conclusion. The software modules obtained as part of the scientific work were implemented in the ACELAN-COMPOS package. Their performance for model problems with quasi-regular finite element grids was estimated. Taking into account the features of the structures of the stiffness and mass matrices obtained through solving the generalized eigenvalue problem for an electroelastic body, the preferred methods for their processing were determined.

Введение. Широкое использование пьезоматериалов в различных отраслях стимулирует изучение их физических характеристик и обусловливает актуальность таких изысканий. В рассматриваемом случае модальный анализ позволяет определить рабочую частоту и коэффициент электромеханической связи пьезоэлементов различных устройств. Эти индикаторы представляют серьезный теоретический и прикладной интерес. Цель исследования — разработка численных методов для решения задачи определения частот резонанса в системе упругих тел. Для достижения цели нужны новые подходы к дискретизации задачи на основе метода конечных элементов и выполнение программной реализации выбранного метода на языке С# на платформе .net. Актуальные решения созданы в контексте библиотеки классов комплекса ACELAN-COMPOS. Основанные на обращении матриц известные методы решения обобщенной задачи на собственные значения неприменимы к матрицам большой размерности. Для преодоления этого ограничения в представленной научной работе реализована логика построения матриц масс и созданы программные интерфейсы для обмена данными о задачах на собственные значения с модулями пре- и постпроцессинга.

Материалы и методы. Для реализации численных методов задействовали платформу .net и язык программирования C#. Валидация результатов исследования проводилась путем сравнения найденных значений с решениями, полученными в известных CAЕ-пакетах (англ. computer-aided engineering — компьютеризированная инженерия). Созданные подпрограммы оценивались с точки зрения производительности и применимости для задач большой размерности. Проводились численные эксперименты с целью валидации новых алгоритмов в задачах малой размерности, которые решаются известными методами в MATLAB. Далее подход тестировали на задачах с большим числом неизвестных и с учетом распараллеливания отдельных операций. Чтобы избежать нахождения обратной матрицы, программно реализовали модифицированный метод Ланцоша. Рассмотрели форматы хранения матриц в оперативной памяти: триплеты, CSR, СSC, SKyline. Для решения системы линейных алгебраических уравнений (СЛАУ) задействовали итерационный симметричный метод LQ, адаптированный к этим форматам хранения.

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

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

Introduction. Devices made of piezoelectric materials have been widely used, actively studied and improved for a long time. Medical ultrasound devices (diagnostic equipment, ultrasonic scalpels) [1–4] and mobile energy generators [

In the study on piezoelectric elements, a key role is played by the modal analysis stage, which enables to establish the resonance and antiresonance frequencies of the device. These data:

The study was aimed at the creation of numerical methods for solving the problem of determining resonance frequencies for a system of elastic bodies. Achieving the stated goal required solving two tasks. The first was to develop methods of discretization of the problem based on the finite element method (FEM). The second was to carry out a software implementation of the selected method in C# on the .net platform. All known programs take into account the context of the ACELAN-COMPOS class library [

Materials and Methods. Principally, the proposed approach was designed to solve static problems of electroelasticity when implementing the averaging method [

The mathematical model of the problem being solved consists of the defining relations [

,,(1)

, , (2)

, (3)

Here, σ— stress tensor; ρj — body density; ε — strain tensor; u — displacement vector; D — electric displacement vector; E — electric-field vector; — body force vector; — electric potential; , , — damping coefficients;, , — tensors of elastic constants, piezoelectric modules and dielectric permittivity; index j — body number in the model.

The discretization is performed by replacing:

Here, — shape function matrix for the displacement field; Nϕ — shape function vector for electric potential; — global vectors of the corresponding nodal degrees of freedom.

In this case, the original problem (1–3) takes the form:

. (4)

Here, matrices M and K are global matrices of mass and stiffness, respectively, and the vector is a general vector of unknowns:

a = [U, Ф].

In the problem of the theory of electroelasticity:

(5)

Matrices , and — symmetric. In the case of harmonic oscillations at natural frequency , it is possible to write:

, denoting the corresponding eigenvector by .

Consider free oscillations if . In this case, task (4) is represented as:

. (6)

Thus, the original problem is reduced to a generalized eigenvalue problem (6). For nonzero , inequality (6) is solved by finding the matrix inverse to . However, at the same time, the sparse matrix becomes full, i.e., the method is unsuitable for large matrices. Therefore, it is needed to use other methods that do not require finding the inverse matrix. To solve this problem, a modified Lanczos method was programmatically implemented in this paper [

Matrices M and K — sparse, with a small number of nonzero elements. Several formats are used to store such matrices in RAM:

The coordinate format involves storing triples (triplets) of values (i, j, k), representing coordinates (i, j) and values (k) of nonzero elements. CSR is sometimes referred to as CRS or Yale format. It involves storing a sparse matrix in the form of three arrays. Consider matrix of N size with NZ nonzero elements. We describe the possible organization of its storage. All nonzero elements must be placed in one array of NZ size. The positions of these elements in the columns should be placed in another array of NZ size, and the third array of N size should be used to store the indices of the first elements of the rows. Similarly, the storage in CSV format is implemented.

The SKS format assumes the storage of a variable-width matrix band that includes all nonzero elements. In this case, zeros are allowed. The efficiency of this format depends on the renumbering of the matrix rows. Methods for reducing the size of the tape are described in [

To solve the SLAE (system of linear algebraic equations), an iterative symmetric LQ method (SYMLQ [

Research Results. At the beginning of the study, we chose the optimal storage format for sparse matrices. The coordinate format enabled to quickly add and change an element of the matrix. These operations were needed at the stage of assembling the global matrix and taking into account the boundary conditions. In addition, for ill-conditioned matrices, to which refers, a preliminary transformation is often used for normalization. It is also convenient to perform it in the coordinate format. However, this format is ineffective when it comes to algebraic operations.

CSR is ill-adapted for changing the structure of the matrix: by adding a nonzero element, you need to insert into two arrays. In this case, the matrix is multiplied by a vector quite easily and efficiently.

SKS has similar problems with the addition of nonzero elements and is highly dependent on the renumbering of unknowns in the problem. We focus on the example of a quasi-regular grid, which is used in the ACELAN-COMPOS package to work with representative volumes of composites. The width of the band containing all nonzero elements can be predetermined and depends on the number of nodes and the type of final element. In the general case of an arbitrary finite element grid, it is difficult to estimate the size of the band in advance.

Four methods of numbering unknowns were used in numerical experiments. Figure 1 shows the structure of stiffness matrices constructed for the same task, but with different renumbering of nodes of a finite element grid.

Fig. 1. Stiffness matrix structure with various methods of node numbering:a — unknowns are ordered by nodes; b — first, displacement nodes are ordered, then potentials; c — nodes are sorted by layers of the FE grid, and unknowns — by nodes;b, d — nodes are sorted by layers of the FE grid, and unknowns — as in example

Thus, the grid was a cube with regular partitioning by eight-node finite elements. A model matrix of 500 lines was used for the illustrations. The matrices shown in Figures 1 a and 1 b, were not subjected to additional renumbering of nodes and differed only in the numbering of degrees of freedom. In 1 a:

.

In 1 b, the unknowns responsible for the potential distribution were collected at the end of the vector:

.

The unknowns in matrices 1 c and 1 d were numbered similarly, but the nodes of the finite element grid were pre-numbered according to their coordinates through alternately sorting all nodes by each of the coordinates. This technique is widely used to build more efficient SLAE solution modules, as it enables to work with the matrix in a suitable band format, convenient for parallelization. Similar external modules were implemented for the ACELAN-COMPOS complex [13–15]; however, in this work, only formats for storing sparse matrices of a general form were used.

Table 1 summarizes the data on the time required to perform basic operations with matrices in various formats.

Table 1

Time to perform basic operations with the stiffness matrix in the problem of the theory of electroelasticity.

19,652 rows

Storage format | Operation | Elapsed time, ms | |||

1 a | 1 b | 1 c | 1 d | ||

CSR | Conversion from coordinate format | 123 | 132 | 97 | 117 |

CSR | Multiplication by vector, 100 operations | 260 | 260 | 260 | 260 |

SKS | Conversion from coordinate format | 690 | 703 | 124 | 268 |

SKS | Multiplication by vector, 100 operations | 60,558 | 61,450 | 7,616 | 22,113 |

The experimental results showed that the conversion operation from a coordinate storage format to a compact one took little time. At the same time, the renumbering of grid nodes to form a block-tape matrix made it possible to get a noticeable increase in performance even without changing the internal structure of the matrix in memory. CSR format turned out to be optimal in terms of the efficiency of the matrix-vector multiplication operation. When the matrix was inverted, SKS format was more efficient, but for problems with the number of unknowns of the order of 103, iterative methods for solving a generalized eigenvalue problem worked noticeably faster.

Further, the performance of the software implementation of the Lanczos method was experimentally evaluated. The contribution of all operations to the total solution time was measured. As a result, it was found that the operation of solving SLAE took up to 95% of the total time of the algorithm. In the course of the algorithm, a Krylov subspace was constructed, and depending on its dimension, the number of SLAE that needed to be solved changed. Note that the dimension of the Krylov subspace was chosen based on heuristics with respect to the number of desired eigenvalues. Here, the SLAE differed only in the right part, so that the requirements for allocated memory remained low. Among the basic operations used in the course of solving SLAE by the SYMLQ method, the greatest computational costs were needed for multiplying the matrix by the vector.

To increase the performance of the algorithm, the simplest parallelization with shared memory was implemented. Blocks of rows were allocated for the CRS format. They were transmitted to separate threads that calculated the corresponding components of the resulting vector. The performance gain was 40–50% when using 8 threads. At the same time, for matrices of the order of 103 elements, the increase was about 40%, and for matrices of the order of 104 — about 50%.

Discussion and Conclusions. Within the framework of this study, a method for solving a generalized eigenvalue problem for matrices obtained by modeling electroelastic bodies was implemented. Software modules were created in C# for constructing mass matrices by the finite element method and performing auxiliary operations within the framework of the Lanczos method (working with Krylov subspace vectors, reorthogonalization, finding eigenvectors). The computational complexity was mainly due to the operations of multiplying sparse matrices by a vector. In this regard, numerical experiments were carried out to determine the optimal formats for storing matrices, the optimal structure of the matrix obtained as a result of renumbering the nodes of the FE grid and degrees of freedom in the nodes. A version of the SYMMLQ iterative algorithm using parallel computing was developed. The final scheme of work included three points. First, global matrices were constructed in coordinate format with a renumbering algorithm (Fig. 1 c). Secondly, the data was converted to CRS format. Thirdly, the data was processed by the Lanczos method, which included the SYMLQ method for solving SLAE. The results of the work were included in the ACELAN-COMPOS software package.

The authors declare that there are no conflicts of interest present.