К расчету амплитудно-частотных характеристик задач об установившихся колебаниях на основе кластерных технологий в ACELAN Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 10, № 1, 2005

К РАСЧЕТУ АМПЛИТУДНО-ЧАСТОТНЫХ

ХАРАКТЕРИСТИК ЗАДАЧ ОБ УСТАНОВИВШИХСЯ КОЛЕБАНИЯХ НА ОСНОВЕ КЛАСТЕРНЫХ ТЕХНОЛОГИЙ

В ACELAN*

К.Е. ВАСИЛЪЧЕНКО НИИ механики и прикладной математики Ростовского государственного университета, Ростов-на-Дону, Россия

e-mail: [email protected]

А. В. НАСЕДКИН

Ростовский государственный университет, Ростов-на-Дону, Россия

e-mail: [email protected]

А. Н. СОЛОВЬЕВ Донской государственный технический университет, Ростов-на-Дону, Россия e-mail: [email protected]

Algorithms of parallel calculations are implemented into ACELAN, a software package based on finite elements method. Specific problems include direct problems of structural vibrations of permanent form and inverse boundary problems in application to elasticity and electroelasticity theories. We discuss an implementation of these algorithms using MPI standard. Numerical examples are presented.

Введение

Современные конечно-элементные программные комплексы, такие, например, как ANSYS [1], предоставляющие возможности расчета связанных физических полей в элементах конструкций, предлагают пользователю широкий арсенал средств и приемов. К последним относятся, в частности, интерактивный режим, предусматривающий непосредственное описание задачи во время сеанса работы с использованием интерфейса пакета, а также пакетный режим, в котором пользователь должен предварительно изучить некоторый специализированный язык команд и затем описать исследуемую задачу с его помощью.

* Работа выполнена при частичной поддержке Российского фонда фундаментальных исследований (гранты № 03-07-90411, № 02-01-01124).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

Возможности использования компьютерных кластеров, которые также предусматривают современные программные комплексы, позволяют решать задачи с большим объемом вычислений. К таким задачам относятся расчеты пьезоэлектрических преобразователей [2], которые находят свое применение во многих технических устройствах и представляют собой существенно трехмерные составные конструкции, нагруженные на газообразные или жидкие среды. Динамические режимы функционирования преобразователей обычно являются установившимися или, что несколько реже, нестационарными. Конечно-элементное моделирование таких процессов [3] приводит к большому объему вычислений, которые целесообразно проводить на компьютерных кластерах.

В последние годы интенсивно развиваются методы решения обратных граничных и геометрических задач механики деформируемого твердого тела, которые, в частности, могут быть сведены к неклассическим краевым задачам [4]. Некоторые из этих методов [5, 6] предполагают решение последовательности прямых задач, конечно-элементные модели которых не изменяются от шага к шагу, что допускает простые схемы распараллеливания таких вычислений.

Специализированный конечно-элементный комплекс ACELAN [2, 3, 7, 8], который разрабатывается с 1996 года в Ростовском государственном университете, предназначен для исследования напряженно-деформированного состояния составных конструкций, состоящих из упругих, электроупругих и акустических тел. "Дружелюбный" русскоязычный интерфейс комплекса ACELAN позволяет применять его широкому кругу пользователей с различным уровнем подготовки как в механике деформируемого твердого тела, так и в вычислительной технике. Вместе с этим последние версии комплекса (ACELAN 0.8.5.5) имеют встроенную консоль команд и пакетный режим работы с файлами, созданными в текстовых редакторах и описывающими исследуемую задачу на специально разработанном языке AceLang [9].

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

При построении в ACELAN амплитудно-частотных характеристик (АЧХ) и при реализации нестационарного анализа пьезоэлектрических преобразователей основной ресурс времени тратится на пошаговое преобразование и решение систем линейных алгебраический уравнений (СЛАУ). В случае установившихся колебаний и построения АЧХ распараллеливание можно осуществить путем разделения частотного интервала и использования на каждом компьютере "нераспараллеленных решателей" прежних версий. При нестационарном анализе распараллеливаются преобразование и решение СЛАУ на каждом временном слое. Алгоритмы параллельных вычислений решения СЛАУ можно найти в [10]. При статическом анализе время, затраченное на решение СЛАУ, сопоставимо со временем построения конечно-элементных объектов, включающего следующие основные операции: формирование матрицы смежности графа геометрической модели и графа узловых неизвестных; перенумерацию узлов; сборку матриц жесткости, масс, демпфирования и учет главных и естественных граничных условий. Алгоритмы описанных подпрограмм, основанные на поэлементном построении КЭ-объектов, допускают простой способ рас-

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

1. Континуальные и конечно-элементные модели

При расчетах пьезоэлектрических устройств из-за наличия пьезоэффекта обязателен учет связанности механических и электрических полей. В случаях взаимодействия пьезоэлек-триков с акустической средой требуется также и учет акустических полей. В комплексе ACELAN для описания движения упругих и пьезоэлектрических тел принята модель линейной теории упругости и электроупругости [11]. Затухание в упругой среде описывается в рамках общепринятой методики КЭ-пакетов [1], потери в электроупругой среде учитываются по оригинальной методике, предложенной в [12]. Движение жидких и газообразных сред описывается в рамках акустического приближения [13]. К уравнениям движения описанных сред добавляются соответствующие главные и естественные граничные условия, условия согласования полей по границам контакта различных сред, импедансные граничные условия для "усечения" акустических областей, а также начальные условия для нестационарных задач.

В общем случае при частичной (по пространственным координатам) конечно-элементной аппроксимации линейных задач акустоэлектроупругости можно получить векторные уравнения, в которых учтены естественные граничные условия, вида [3]

M ■ a + C ■ a + K ■ a = F, (1)

где

M

Muu 0 0 0 »T

иф 0

И«ф 0 — Мфф

C

K

Kuu K

KT -K

uu

T —

up AV PP 00

0 0

—Кфф

F

F

C ии

CdKUp

T

Ииф

Fu

p + ZdF 0

0 И-иф 00 0 —Сфф

p

a

U Ф Ф

(2)

(3)

и для нестационарных задач к (1) следует добавить начальные условия.

Здесь неизвестный вектор узловых степеней свободы а = а(£) может содержать узловые механические перемещения и, электрические потенциалы Ф и потенциалы скоростей акустической среды Ф; — коэффициент демпфирования для пьезоэлектрического материала [3, 12]. Блоки матриц М, С, К и вектора Е формируются при учете взаимодействия полей, которые указаны в нижних индексах. Обозначения в нижних индексах соответствуют следующим полевым функциям: и — поле механических перемещений, (р — поле электрического потенциала, ф — потенциал поля скоростей акустической среды.

Входящие в (2), (3) блоки Мии, Мфф и Сфф симметричны и положительно определены,

а блоки Kuu, Kpp и K

фф симметричны и неотрицательно определены.

Следующий этап формирования разрешающей системы уравнений метода конечных элементов (МКЭ) связан с учетом главных граничных условий и других связей. Для удобства дальнейших операций представим (1) в компактной форме [14]

p dk

dtfc a™ = Fra' (4)

k=0

Здесь p = 2, а при n = 0 A00 = K, A01 = C, A02 = M, F0 = F и a0 = a.

Первый шаг n =1 этого этапа связан с поворотом систем координат узловых степеней свободы, например, для удобства наложения ограничений на нормальные un или касательные uT компоненты векторов перемещений на границах, не параллельных координатным осям глобальной системы. Следует отметить [14], что эти преобразования не меняют структуру и свойства соответствующих объектов.

Теперь для учета связей (при n = 1) разобьем все множество степеней свободы a™ на множество главных степеней свободы, из которых составим вектор ac, и множество подчиненных степеней свободы, из которых составим вектор as. В соответствии с этим разбиением представим КЭ-объекты из (4) в новой блочной форме

A™k=( Aik tk) • f™={ £} • an={:}' (5)

Рассмотрим уравнения связей вида

as = TsC ■ ac + s, (6)

где rank Tsc = dim ac.

Уравнения (6) включают следующие важнейшие типы связей:

— задание отдельных степеней свободы;

— связи групп степеней свободы в одну, что в задачах расчета пьезоустройств необходимо для определения эквипотенциальных участков (электродов);

— произвольные линейно независимые уравнения связей.

Учитывая уравнения (5), (6) и гипотезу идеальности связей, получаем первую форму уравнений МКЭ:

:dtk'

Е

A-ra a = FC. (7)

k=0

VT Л T , тлТ

Ak = Acck = Acck + Acsk Tsc + TscAcsk + TscAssk Tsc (8)

p dk c * С I J-scA's — (TscAssk + Acsk) dfk '

FC = Fc + TSCFs -J](TSCAssk + Acsk) dtkS' (9)

&=0

Переход от (4) к (7)-(9) можно рассматривать как процедуру конденсации степеней свободы при неоднородной зависимости (6). Система (7) имеет меньшую размерность по сравнению с размерностью матриц Лга& в (4), что может потребовать перестройки форм хранения КЭ-объектов и некоторых дополнительных технических процедур.

Для сохранения размерности КЭ-матриц добавим к системе (7) уравнения связей (6). После процедур симметризации полученной системы и выравнивания порядков коэффициентов получим вторую форму уравнений МКЭ:

p ,ik

i ac = Fc- L (io)

k=0

dtk as F

s

д'' _ ( Асс0 + Асс0 Ас«0 \ д'' — Асст 0 \ ™ _ 1 о.

Ао — 1 асто А;«, ,)' Ат ч 0 о), т — 1,2; (11)

ЕС' - 1_ЕС + ЕС; Е«]т; Е« - А'^в, ЕС - -Т^Е«; (12)

А' = -Тт А' А' = -Т А'т А' = А (13)

лс«0 — 1 БС^ввО, лсс0 — 1 8С^-СВ 0, — л ««0, (13)

где А^0 — диагональная матрица, составленная из диагональных компонент блока А««0 матрицы Ап0.

Уравнения МКЭ в форме (10) имеют размерность исходных матриц Ага& и представлены в симметричном виде, удобном для применения численных процедур. Отметим, что в (7), (10) структура уравнений МКЭ задач акустоэлектроупругости представлена в общем виде. Реально форматы хранения матриц и векторов МКЭ зависят от алгоритмов решения задач (7) или (10) и стратегий формирования конечно-элементных объектов и работы с ними.

В АСЕЬАМ в процедурах статического, гармонического и нестационарного анализа используются уравнения в форме (10) в сочетании с полностью "перемешанной" их формой, когда в векторе а последовательно для каждого узла записываются компоненты различных степеней свободы.

2. Параллельные алгоритмы расчета АЧХ задач об установившихся колебаниях

Конечно-элементные модели установившихся колебаний с частотой f составных упругих, электроупругих и акустических тел получаются из (1) или (10) с помощью замены

а — а(х) ехр^и), Е — Е(х) ехр^и1), и — 2пf. (14)

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

-и2 М ■ а + ■ а + К ■ а — Е, (15)

которая в силу структуры (2), (3) приводится к симметричной форме

Aeff (и) ■ а — Ее// (16)

после умножения второго блока уравнений на множитель (1 + ^и^)-1.

Рассмотрим далее задачу построения в АСЕЬАМ АЧХ на интервале частот и £ [и1, и2]. Тогда кроме этого частотного интервала пользователю необходимо задать количество N точек, равномерно распределенных в этом интервале. Отметим, что для получения адекватных экспериментальным данным АЧХ N обычно должно иметь порядок нескольких десятков и более в зависимости от размера частотного интервала и его положения относительно резонансных частот. Основной ресурс времени в этом случае связан с преобразованием матрицы системы (16) (шаги п — 1, 2 описанного выше алгоритма) и ее решением на каждом частотном шаге. Пусть в распоряжении пользователя имеется кластер, состоящий из АрГОС одинаковых процессоров со скоростью передачи информации ус. Процесс подготовки задания, триангуляция области производится пользователем на центральном

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

Применяемый в ACELAN алгоритм распараллеливания вычислений в "решателе" может быть представлен следующей схемой.

1. Центральный процессор считывает данные о задаче и рассылает их на процессоры кластера за время Tin.

2. Каждый процессор за время Tprep производит подготовку конечно-элементных объектов — матриц и векторов в уравнении (15), которые в дальнейшем хранятся в памяти каждого процесса.

3. В цикле по частотному интервалу каждый процесс обрабатывает свой набор частот (шаги n = 1, 2 разд. 2, решение СЛАУ (16)) и после каждого решения посылает найденный вектор a центральному процессу, который осуществляет вывод этих данных в соответствующий файл. Обозначим время обработки одного частотного варианта через T1, время вывода одного варианта — через Tout = Vr /vout и время передачи-приема вектора a — через Tr = Vr/vc, где vout — скорость вывода информации в файл, Vr — объем информации в векторе а, который пропорционален некоторому параметру S, характеризующему "сложность" задачи (количество узлов конечно-элементной сетки или степеней свободы).

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

T = Tprep + N (Ti + Tout). (17)

Будем считать для простоты, что N кратно Nproc. Тогда время Tc выполнения задания на рассматриваемом кластере определится соотношением

Tc = Tin + Tprep + — (N - N^-) + Ti + N Tout. (18)

vc Nproc Nproc

Ускорение E алгоритма можно определить посредством отношения

E = (19)

Tc

При большом N ускорение алгоритма, кроме числа процессоров Nproc, определяется тремя скоростями vout, vc и скоростью обработки одного частотного варианта vs = V /T1 следующим приближенным соотношением:

E ^- -1 + Vsout-, (20)

Nproc (1 - Vsc) + Vsc + Vsout

где Vsout = vs/vout, Vsc = vs/vc.

Очевидно, что при большом числе процессоров Nproc имеет место асимптотическое равенство

E ^ 1 + Vsout , (21)

Vsc + vsout

несложный анализ которого указывает пути совершенствования предложенного алгоритма. При известных скоростях решения, вывода и обмена соотношение (20) может служить основой для проектирования кластера заданного ускорения. Отметим, что если vc ^ vs и vout ^ vs, то ускорение кластера в (20) приближается к числу процессоров Nproc. Такая ситуция имеет место для задач с "большой сложностью", так как время вывода и обмена растет пропорционально объему V, а время решения — значительно быстрее.

3. Численный пример реализации кластерного алгоритма

На базе лаборатории открытых систем НИИ механики и прикладной математики РГУ разработан вычислительный гетерогенный кластер, состоящий фактически из трех подкла-стеров. Каждый подкластер состоит из узлов примерно одинаковой производительности, при этом производительность узлов из разных подкластеров различается. 1-й подкластер (основной) состоит из 22 вычислительных узлов на базе Intel Celeron 2500 МГц с включением в него четырех узлов меньшей производительности на базе процессоров Pentium-II с тактовой частотой 300-400 МГц, что, однако, не влияет на общую его производительность, поскольку в нем доминируют узлы высокой производительности. При этом узлы, входящие в этот подкластер, укомплектованы сетевыми адаптерами, поддерживающими технологию Gigabit Ethernet. Это позволяет не только ускорить собственно расчетную часть, но и обеспечить высокую скорость обмена данными между узлами кластера. Второй под-кластер состоит из восьми узлов, тактовая частота процессора каждого узла находится в пределах 600-800 МГц. Третий подкластер состоит из 11 узлов на базе процессоров Pentium-II 450 MHz. В качестве операционной системы для кластеров используются Linux Mandrake 9.0/ALT Linux Junior 2.0 и Windows NT/2000/XP. На каждом из узлов обеспечивается двойная загрузка — Linux и Windows. Таким образом, фактически имеются два кластера, которые, однако, могут совместно работать. Разработка программ ведется на языках C/C++ и Фортран с включением библиотек MPI. Каждый из кластеров управляется своим отдельным сервером, причем серверы не предполагают двойной загрузки. В Linux-кластере в качестве среды параллельного программирования используется библиотека LAM, реализующая интерфейс MPI. В Windows-кластере используется библиотека WMPI 1.3, реализующая тот же интерфейс MPI. Она же позволяет организовать расчет задачи в гетерогенной среде одновременно на Linux и Windows.

В численном эксперименте рассматривался (рис. 1) цилиндрический многослойный

Рис. 1. Половина осевого разреза контейнера.

контейнер и изучалось влияние наличия дефектов сплошности конструкции и их размеров на АЧХ электрического потенциала свободных электродов пьезокерамических элементов, расположенных на боковой поверхности и нижнем торце контейнера. Пьезоэлементы принимают упругие волны, возбуждаемые пьезокерамическим излучателем, находящимся на верхней крышке контейнера. В качестве дефекта сплошности рассматривалось отверстие в центре нижнего торца контейнера, частично заполненное менее жестким материалом, чем внешняя оболочка контейнера, выполненная из стали. В качестве параметра, контролирующего сплошность конструкции, была выбрана разность электрического напряжения. (Эта величина, наряду с электрическим импедансом, наиболее просто может быть измерена в натурном эксперименте.) Конечно-элементная сетка для конструкции без дефектов содержала 1346 узлов и 587 конечных треугольных квадратичных шестиузловых элементов.

Первым шагом численного эксперимента было нахождение диапазона частот, в котором АЧХ электрического потенциала для конструкции с дефектом и без него различались бы наиболее существенно (так, чтобы это различие можно было обнаружить не только в численном, но и в натурном эксперименте). Поэтому был выбран достаточно широкий частотный интервал f е [0; 15 000] (f в Гц), а число точек, в которых проводились вычисления, равнялось N = 1501, т. е. шаг на частотном интервале составлял 10 Гц.

На рис. 2 представлены АЧХ (f е [10 000; 14 500]) электрического потенциала ^ свободного электрода приемника, расположенного на нижнем торце контейнера, цифрой 1 отмечены результаты для конструкции без дефекта, а цифрой 2 — с дефектом (отверстие радиуса 0.01 м во внутреннем контейнере и заполнителе в центре нижнего торца, причем в районе дефекта материал внешней оболочки заменен на материал заполнителя). Как следует из этого рисунка, искомым частотным интервалом является интервал [11 890; 12 200], для которого характерны колебания с наибольшими амплитудами вертикальных смещений в окрестности дефекта.

На рис. 3 изображена (черно-белая копия окна постпроцессора ACELAN) форма колебаний конструкции без дефекта при f = 12 130 Гц (тонами серого цвета задано распределение вертикальных смещений). На рис. 4 представлены аналогичные результаты для конструкции с дефектом при f = 11890 Гц, когда радиус отверстия равен 0.01 м.

В табл. 1 представлены результаты, отражающие зависимость собственной частоты колебаний от поперечного размера (радиуса) дефекта. Здесь следует отметить, что частоты, соответствующие рис. 3 и 4, близки к резонансным частотам, расположенным в крайних столбцах табл. 1.

(р

■240-------

■320;-------

10000 11000 12000 13000 14000 /

Рис. 2. Амплитудно-частотная характеристика электрического потенциала.

Рис. 3. Форма колебания контейнера без дефекта. Рис. 4. Форма колебания — с отверстием.

В рассмотренном примере отношение средних скоростей решения и вывода составило ^о^ = 'Уз/'Уо^ = 0.465, при этом относительный разброс не превысил 10 %. Время Тс передачи-приема одного и того же объема информации существенным образом зависело от конфигурации кластера и в рассмотренном примере изменялось от Тс = 5.0 • 10-4 с до Тс = 1.6 с. Таким образом, отношение скоростей решения и передачи-приема изменялось от узс = г^/г^ = 1.2 • 10-3 до узс = 3.7. В табл. 2 представлены результаты, отражающие ускорение кластерного алгоритма для АрГОС = 6, в зависимости от числа точек N по частоте, в которых проводились вычисления. При этом ускорение, вычисленное по формуле (20), равнялось Е = 2.32 для = 1.2 • 10-3.

В табл. 3 при = 100 представлены результаты, отражающие зависимость ускорения алгоритма от числа процессоров ^ргос.

Результаты численных экспериментов и анализ соотношений (20) и (21) открывают пути совершенствования предлагаемого алгоритма, так как, например, при уменьшении отношения на порядок при ^ргос = 6 ускорение алгоритма возрастает с 2.32 до 4.89. (Последнего можно добиться, изменив схему вывода результатов, применяя "бесформатный" вывод и запись в единый файл.)

Таблица!

Радиус, м 0 0.001 0.05 0.01

Частота, Гц 12135 12 050 12 015 11905

Таб лица

100 500 800

Е — численный эксперимент 1.94 2.14 2.21

Т абл ица3

^ргос 2 3 4 6

Е — численный эксперимент 1.45 1.73 1.82 1.94

Е — формула (20), рзс = 1, 2 • Ю-3 1.52 1.83 2.05 2.32

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

Список литературы

[1] ANSYS. Theory Reference. Rel.7.0 / Ed. P. Kothnke. ANSYS Inc. Houston, 2003.

[2] Белоконь А.В., Наседкин А.В., НикитАЕВ А.В. и др. Новая версия пакета ACELAN для для проведения расчетов пьезоизлучателей и пьезоприемников акустических волн // Пьезотехника-2002. Междунар. научно-практ. конф. "Фундаментальные проблемы пьезоэлектрического приборостроения". Тверь, 17-21 сент. 2002: Сб. докл. Тверь: ТвГУ, 2002. С. 171-179.

[3] Белоконь А.В., Наседкин А.В., Соловьев А.Н. Новые схемы конечно-элементного динамического анализа пьезоэлектрических устройств // Прикл. математика и механика. 2002. Т. 66, вып. 3. С. 491-501.

[4] Ватульян А.О., Ворович И.И., Соловьев А.Н. Об одном классе граничных задач в динамической теории упругости // Там же. 2000. Т. 64, вып. 3. C. 373-380.

[5] Козлов В.А., Мазья В.Г., Фомин А.В. Итерационный метод решения задачи Коши для эллиптических уравнений // Журн. вычисл. математики и мат. физики. 1991. Т. 31. С. 45-52.

[6] Weikl W., Andra H., Schnack E. An Alternating iterative algorithm for the reconstruction of internal cracks in a three-dimensional solid body // Inverse Problem. 2001. Vol. 17. P. 1957-1975.

[7] Акопов О.Н., Белоконь А.В., НАдолин К.А. и др. Статический анализ пьезоэлектрических устройств в ACELAN. I. Структура и возможности. II. Численные эксперименты // Современные проблемы механики сплошной среды: Тр. IV Междунар. конф. Ростов-на-Дону, 27-28 окт. 1998. Т. 1. Ростов-на-Дону: Изд-во СКНЦ ВШ, 1999. С. 14-22.

[8] Белоконь А.В., Петушков А.Л., Соловьев А.Н. О реализации кластерных технологий при конечно-элементном моделировании пьезоэлектрических устройств в ACELAN // Пьезотехника-2003. Междунар. научно-практ. конф. "Фундаментальные проблемы пьезоэлектрического приборостроения". Москва, 26-29 ноября. 2003: Сб. докл. М.: МИРЭА, 2003. С. 239-240.

[9] Петушков А.Л., Соловьев А.Н. О реализации пакетного режима работы конечно-элементного комплекса ACELAN // Современные проблемы механики сплошной среды: Тр. VIII Междунар. конф. Ростов-на-Дону, 14-18 окт. 2002. Т. 2. Новая книга, 2003. С. 149-153.

[10] Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. М.: Мир, 1991. 275 с.

[11] Партон В.З., Кудрявцев Б.А. Электромагнитоупругость пьезоэлектрических и электропроводных тел. М.: Наука, 1988. 472 с.

[12] Наседкин А.В. Новая модель учета демпфирования для конечно-элементного пьезоэлектрического анализа // Современные проблемы механики и прикладной математики: Мат. шк.-семинара. Воронеж, 25-30 сент. Ч. 2. Воронеж: ВГУ, 2000. С. 319-323.

[13] КрАСильников В.А., Крылов В.В. Введение в физическую акустику. М.: Наука, 1984. 400 с.

[14] Акопов О.Н., Белоконь А.В., Нлдолин К.А. и др. Симметричные седловые алгоритмы конечно-элементного анализа составных пьезоэлектрических устройств // Мат. моделирование. 2001. Т. 13, № 2. С. 51-60.

Поступила в редакцию 11 марта 2004 г.