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 [5] should be noted separately. Paper [6] described the combination of photo- and piezoelectric effects to create efficient compact energy sources. New materials designed for the application under specific conditions are being studied in science and industry. In [7], the creation of a lead-free piezo-active composition suitable for operation at various temperatures was considered.

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 [8]. When solving a generalized eigenvalue problem, methods based on matrix inversion are widely used. However, they are not applicable to large-dimensional matrices. In the presented research, this limitation was overcome as follows:

Materials and Methods. Principally, the proposed approach was designed to solve static problems of electroelasticity when implementing the averaging method [9], which was used to calculate the effective properties of piezo composites. In this regard, only stiffness matrices were presented at the stage of constructing global FEM matrices. In this study, we additionally implemented the logic of constructing mass matrices and developed software interfaces (application programming interface, API) for exchanging data on eigenvalue tasks with pre- and postprocessing modules. The developed routines were evaluated in terms of performance and applicability for large-scale tasks. Numerical experiments were carried out to validate the algorithms created for such small-dimensional problems that provide obtaining a solution by general methods in the MATLAB computing package. Next, testing was performed on tasks with a large number of unknowns and taking into account the parallelization of individual operations.

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

,,(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 [10]. The author of this modification is T. S. Martynova. The description of the development is not given in this article. Of the operations used in the method, the most expensive from the point of view of computational resources was the solution to a system of linear algebraic equations (SLAE), needed for performing a spectral transformation.

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 [11]; however, their applicability to the stiffness matrix obtained when solving a three-dimensional problem using FEM requires a separate study.

To solve the SLAE (system of linear algebraic equations), an iterative symmetric LQ method (SYMLQ [12]) adapted to the storage formats listed above was used.

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.

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.

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.