CC BY
19
ПРОМЫШЛЕННОЕ РЫБОЛОВСТВО, АКУСТИКА
УДК 539.3
С.М. Балабаев, Н.Ф. Ивина
Дальневосточный государственный технический рыбохозяйственный университет,
690087, г. Владивосток, ул. Луговая, 52б
МОДИФИЦИРОВАННЫЙ МЕТОД КОНЕЧНЫХ И ГРАНИЧНЫХ ЭЛЕМЕНТОВ ДЛЯ АНАЛИЗА ГИДРОАКУСТИЧЕСКИХ ПЬЕЗОПРЕОБРАЗОВАТЕЛЕЙ
Рассмотрен модифицированный метод конечных и граничных элементов для компьютерного моделирования гидроакустических пьезопреобразователей произвольной формы и размеров. Введено понятие матрицы импеданса излучения. Выполнен анализ основных параметров цилиндрического водозаполненного пьезопреобразователя с радиальной поляризацией. Определены оптимальные геометрические размеры пьезопреобразователя. Рассмотрено влияние паразитной изгибной моды.
Ключевые слова: гидроакустический пьезопреобразователь, метод конечных элементов, метод граничных элементов.
S.M. Balabaev, N.F. Ivina MODIFIED FINITE ELEMENT-BOUNDARY ELEMENT METHOD FOR ANALYSIS OF HYDROACOUSTIC TRANSDUCERS
A modified finite element-boundary element method for computer simulation of hydroacoustic piezoelectric transducers of arbitrary shape and size is considered. The concept of the radiation impedance matrix is introduced. Analysis of basic parameters of free-flooded piezoelectric cylinder transducer with radial polarization is executed. The optimal geometrical dimensions of piezoelectric transducer are determined. The influence of the parasitic bending mode is considered.
Key words: hydroacoustic piezoelectric transducer, finite element method, boundary element method.
Введение
Пьезокерамические преобразователи (пьезопреобразователи) являются обязательной и важной частью различных акустических устройств, применяемых в гидроакустике, рыбопоисковой технике, радиотехнике, автоматике, медицине. Классические задачи анализа пьезопреобразователей решены на основе их одномерных моделей. Основной недостаток одномерных моделей состоит в том, что они дают удовлетворительные результаты только для основной моды колебаний пьезопреобразователей простейших форм при определенных ограничениях на соотношение их резонансных и нерезонансных размеров. Для преобразователей с произвольным соотношением размеров задача становится принципиально трехмерной и описывается не обыкновенными дифференциальными уравнениями, а уравнениями электроупругости в частных производных со сложными граничными условиями. К настоящему времени аналитические решения этих задач в точной постановке не получены. Ограниченные возможности аналитических методов приводят к тому, что многие вопросы разработки пьезопреобразователей можно решить только опытным путем, что связано с большими затратами времени и материалов.
Компьютерное моделирование на основе современных численных методов: метода конечных элементов (МКЭ) и метода граничных элементов (МГЭ) - позволяет анализиро-
вать реальную конструкцию при варьировании ее параметров и таким образом значительно ускорить и удешевить разработку высокоэффективных пьезопреобразователей. Высокая производительность современных персональных компьютеров изменила соотношение между численными и аналитическими методами в пользу первых. Современные численные методы следует считать наиболее эффективным и перспективным математическим аппаратом для анализа реальных пьезопреобразователей. Важными преимуществами МКЭ и МГЭ являются эффективность, относительная простота, универсальность: многие подпрограммы не зависят от геометрической формы преобразователя. Таким образом, авторы солидарны с мнением известного сторонника численных методов в акустике В.Ю. Завадского [1]: «Высокая производительность компьютеров изменила соотношение между вычислениями при использовании теоретических решений в виде формул для волновых задач и численными алгоритмами, схемами, разностными методами, применяемыми непосредственно к исходным уравнениям задачи. Конечно, формулы сохранили свое значение как более наглядные средства представления решения простых идеализированных задач. В то же время стало ясно, что в более сложных задачах необходимо использовать всю мощность компьютера и применять универсальные алгоритмы».
Работы по применению современных численных методов, давшие мощный толчок компьютерному моделированию пьезопреобразователей, были начаты в 70-х гг. прошлого столетия практически одновременно и независимо в США и Японии. В статье [2] изложен один из методов учета жидкой акустической нагрузки, основанный на применении граничного интегрального уравнения Гельмгольца, фактически МГЭ с кусочно-постоянной аппроксимацией искомых функций, хотя сам этот термин в то время еще не сформировался, а был введен позже.
Объекты и методы исследований
Авторами в статье [3] разработан комбинированный метод конечных и граничных элементов для анализа пьезопреобразователей произвольной формы и размеров. В комбинированном методе преобразователи моделируются пьезоэлектрическими конечными элементами второго порядка. Для учета акустического излучения используется граничное интегральное уравнение Гельмгольца, которое интегрируется численно с использованием граничных элементов с квадратичной аппроксимацией. Использование таких элементов значительно увеличивает точность расчетов и уменьшает размеры глобальных матриц, а следовательно, и время вычислений. Применение комбинированного метода позволяет моделировать пьезопреобразователи в точной постановке и учесть все эффекты, которые принципиально не могут быть учтены при использовании аналитических методов. Как показано в статье [3], комбинированный метод конечных и граничных элементов приводит к системе двух матричных уравнений большой размерности относительно узловых смещений и давления. В данной работе представлена модификация комбинированного метода с введением импедансной матрицы, что позволит записать окончательное выражение в виде одного матричного уравнения большой размерности.
Результаты и их обсуждение
Как показано в монографии авторов [4], колебания пьезопреобразователя произвольной геометрической формы и размеров описываются матричной системой линейных алгебраических уравнений большой размерности
{[нии ]- с2 [м ]+ с[с Ми. }= [ит}+ [^} (1)
Обозначения в (1) соответствуют принятым в [4]: [Иии ] - глобальная матрица жесткости с учетом электрических граничных условий; с - круговая частота; [м] - глобальная
матрица массы; г - мнимая единица; [С] - демпфирующая матрица для учета механических потерь; {иг} - вектор узловых смещений; \Ит} - глобальная матрица-столбец пьезоэлектрической «жесткости»; {р} - вектор узловых сил.
Для решения задач излучения нужно определить вектор узловых сил, учитывающий реакцию жидкости на колебания пьезопреобразователя, и описать излучаемое акустическое поле. Для этих целей используется граничное интегральное уравнение Гельмгольца, которое определяет давление в произвольной точке пространства, если известны распределения давления и нормальной компоненты колебательной скорости на замкнутой поверхности £ преобразователя [5]:
С (X )Р (X )= ДО (Р (0 )дО (X, б)/дп + гюрО (X, б ) ()), (2)
где P(X), Р(0) - давления в точке X акустического поля и в точке 0 на поверхности 5; р -плотность жидкости; Уп(0) - колебательная скорость, нормальная к поверхности 5; п -внешняя по отношению к излучающей поверхности единичная нормаль; G(X,Q) - функция Грина, которая в общем случае трехмерного пространства равна
О (X, 0) = О (Я) = ехр( -гкЯ)/ Я.,
где Я - расстояние между точками Xи 0; к - волновое число, к = ю/с, с - скорость звука в жидкости. Коэффициент С(^) в уравнении (2) определяется в зависимости от расположения точки X относительно излучающей поверхности 5. С(^=4 п, если точка X находится в пространстве вне поверхности 5; С(А')=0, если точка X находится внутри 5.
Для интегрирования уравнения (2) излучающая поверхность 5 разбивается на конечное количество граничных элементов, число узловых точек которых зависит от порядка аппроксимирующих функций. Порядок аппроксимации в пределах граничного элемента учитывается выбором соответствующих функций формы N (х1, х2), которые определены в
локальных координатах х}, х2. Тогда нормальную скорость и давление в пределах граничного элемента с номером к можно выразить следующим образом:
I+1 I +1
V* = Е , Рк = Е ир, (3)
1 =1 1 =1
где I - порядок аппроксимации функций; ¥}к, Р1к - узловые значения функций. Интегралы в уравнении (2) можно представить как сумму интегралов по граничным элементам.
При выполнении интегрирования в пределах каждого граничного элемента удобно использовать локальные координаты; тогда компоненты элементных матриц можно ввести следующим образом:
ак]т = Л N 1дО / дп ёег [/ ]х1 дх 2; Ьк]т = ДО N }О ёе! [/ ]дх1 дх 2;
5к 5 к
скт = ДО дЯ-1/дп ёе! [/ ]х 1 дх 2, (4)
5 к
где ёе/ - якобиан.
Если точку т располагать последовательно во всех узловых точках, то получим систему линейных уравнений относительно граничных узловых значений давления и нормальной скорости
Цатр -(4п +1ст)Рт ЬтУЧ, (5)
к ] т т ]
где т=1+М; М- количество узловых точек и количество уравнений системы (5).
Определим векторы давления {Р} и нормальной скорости {¥„} как совокупность соответствующих узловых значений и перепишем систему уравнений (5) в виде
[л\р}=-Ю)р[В ]{Гп}, (6)
где [Л] и [В] - глобальные матрицы МГЭ, полученные специальным суммированием соответствующих элементных матриц (4). Уравнение (6) можно использовать для решения задач Неймана и Дирихле, определяя решение во внутренней или внешней области.
Для заданной точки X можно вычислить давление в акустическом поле по известным граничным значениям давления и нормальной скорости в соответствии с уравнением (2)
Р = (4п)-1 ({Л } {Р}+ шр{Вх }Т V }), (7)
где {Л1}т, {В1}Т - матрицы-строки для заданной точки акустического поля. Выражение (7) применяется для определения поля изучения в ближней и дальней зоне, а также границы ближней зоны.
Введем понятие матрицы импеданса излучения, ее можно получить из системы уравнений (6) по определению как отношение силы к скорости узловых точек:
[7 ]=-1ар[Б ][Л ]-1 [В ], (8)
где [5] - диагональная матрица площади.
Импеданс излучения зависит от геометрической формы пьезопреобразователя и частоты. При использовании классических одномерных моделей импеданс - скалярная величина, которая характеризует реакцию акустической среды и равна отношению силы, действующей на излучающую поверхность, к ее колебательной скорости. Если излучающая поверхность представлена граничными элементами с узловыми точками, то можно определить собственный и взаимный импеданс каждой точки (точнее небольшой области вблизи этой точки) и сформировать импедансную матрицу.
Реакцию акустической среды выразим через граничное давление
}=-№ ][5 ]Р }, (9)
где [Ж] - матрица направляющих косинусов нормали к поверхности излучения. Давление на границе определим из уравнения (6):
{р}=-1юр[л]~1 [В ]Уп}. (10)
Выразив нормальные скорости на поверхности излучателя через узловые смещения
[V„ } = ia[W ]T [и,}, (11)
из (10) получим
[P}=a2pW [B ][W ]T [u,}. (12)
Перепишем (9) с учетом (8) и (12):
[F }=-ia[W ][Z ][W ]T [u,}. (13)
Подставив выражение (13) в систему (1), получим окончательную систему уравнений модифицированного комбинированного метода конечных и граничных элементов с учетом акустического излучения в виде одного матричного уравнения большой размерности
([Huu ]-a2 [M ]+ ia[C ]+ ia>[W ][Z ][W ]t }=[Huv }. (14)
Из системы (14) можно определить узловые смещения пьезопреобразователя, затем по выражению (11) - нормальную скорость узловых точек, выходящих на поверхность излучения, и давление в этих точках (12). Далее определяется давление в дальней и ближней зоне по выражению (7).
На основе изложенного алгоритма разработан пакет программ для анализа осесиммет-ричных гидроакустических преобразователей. В осесимметричных (фактически двумерных) задачах применяются кольцевые конечные и граничные элементы, поэтому разбивать на элементы необходимо только осевое сечение преобразователя. Для расчетов используются конечные и граничные элементы второго порядка (с квадратичной аппроксимацией). Для демонстрации возможностей разработанного алгоритма и пакета программ рассмотрим цилиндрический радиально поляризованный пьезопреобразователь с открытой внутренней полостью (в отечественной терминологии - водозаполненный, в американской -free-flooded - «свободно-затопленный»).
Водозаполненный пьезопреобразователь был исследован экспериментально в пионерской работе [6]. Было показано, что ниже основного «рабочего» радиального резонанса пьезоцилиндра наблюдается еще один резонанс, в дальнейшем названный объемным. Этот резонанс обусловлен излучением внутренней цилиндрической поверхности и резонансом внутреннего объема жидкости и части окружающей жидкости.
Эта работа стимулировала большой интерес гидроакустиков к разработке математических моделей водозаполненного пьезоцилиндра, являющегося широкополосным двухрезо-нансным гидроакустическим преобразователем, анализу и оптимизации его параметров. Как отмечено в статье [7], эта задача является «интересной и содержательной задачей акустики, имеющей богатую историю». Заметим также, что для первого автора знакомство с этой задачей произошло в январе 1972 г. при подготовке дипломного проекта в одном ленинградском «почтовом ящике» (т.е. это - первая и, наверное, последняя «акустическая любовь»).
Приближенные аналитические модели водозаполненного пьезоцилиндра были разработаны рядом зарубежных и отечественных акустиков [5, 7, 8]. Однако они не учитывают многие принципиально важные эффекты.
В статье [8] разработана численно-аналитическая модель водозаполненного пьезоци-линдра с радиальной поляризацией. Для моделирования пьезоцилиндра использована теория оболочек, а для численного интегрирования граничного интегрального уравнения Гельмгольца - известная программа SHIP, ранее разработанная под руководством автора статьи [8]. Как показано авторами данной статьи в работе [9], теория оболочек, а следовательно, и предложенная комбинированная модель в статье [8] не учитывают изгибные моды, поэтому полученная на ее основе информация будет неполной.
Рассмотрим водозаполненный пьезоцилиндр с радиальной поляризацией, выполненный из пьезокерамики ЦТБС-3, параметры которой соответствуют справочным данным [10], на внешней и внутренней цилиндрических поверхностях которого нанесены сплошные серебряные электроды. Размеры цилиндра: a - внешний радиус, l - высота, h - толщина стенки. В дальнейшем размеры определяются безразмерными отношениями: l / a и h / a, для краткости - просто l и h .
Для анализируемого пьезопреобразователя рассчитывались частотные характеристики чувствительности в режиме излучения в точке дальнего поля в радиальном направлении (на плоскости симметрии цилиндра, в дБ относительно уровня 1 мкПа/В при 1 м), нормированные диаграммы направленности в вертикальной плоскости, частотные зависимости безразмерной активной составляющей акустической мощности, а также распределения модулей нормальной (радиальной) компоненты колебательной скорости и давления на внешней и внутренней цилиндрических поверхностях пьезоцилиндра.
Ранее [3] для водозаполненного тангенциально поляризованного пьезоцилиндра было введено понятие оптимальных геометрических размеров. Под ними понимаются размеры, при которых чувствительность излучения на объемном и радиальном резонансах примерно одинакова, а неравномерность частотной характеристики не превышает 3 дБ.
Частотные характеристики чувствительности в режиме излучения: а - h = 0,2, 1 - l = 0,7, 2 - l = 0,8, 3 - l = 0,9; б - h = 0,15, 1 - l = 1,0, 2 - l = 1,1, 3 - l = 1,2 The frequency sensitivity characteristics in the radiation mode:
а - h = 0,2, 1 - l = 0,7, 2 - l = 0,8, 3 - l = 0,9; б - h = 0,15, 1 - l = 1,0, 2 - l = 1,1, 3 - l = 1,2
На рисунке, а приведены частотные характеристики чувствительности излучения v радиально поляризованного пьезоцилиндра с толщиной стенки h = 0,2 при высотах, близких к оптимальным. По оси абсцисс отложена безразмерная частота kta, kt - волновое число поперечной волны в пьезокерамике. Из этого рисунка следует, что оптимальная вы-
сота равна примерно 0,85. При высоте пьезоцилиндра меньше оптимальной чувствительность на объемном резонансе больше, чем на радиальном, а при высоте цилиндра больше оптимальной - наоборот. Максимальная чувствительность в режиме излучения для ради-ально поляризованного цилиндра равна примерно 137 дБ, т.е. меньше, чем у тангенциально поляризованного пьезоцилиндра [3]. Это объясняется тем, что секционированный пье-зоцилиндр работает на продольном пьезоэффекте, а радиально поляризованный - на поперечном. Поэтому в первом случае коэффициент электромеханической связи больше, чем во втором. Это хорошо известно из приближенной одномерной теории и из практики.
Отметим, что в известной статье [6] измерена чувствительность пьезоцилиндра из керамики цирконата титаната свинца (конкретный тип не указан) с размерами: наружный радиус 1,59 см, толщина стенки 0,32 см, высота цилиндра 1,21 см, что соответствует принятым безразмерным отношениям: I = 0,76, к = 0,2. Из соответствующей экспериментальной кривой следует, что чувствительность на радиальном резонансе меньше, чем на объемном, так как высота цилиндра меньше оптимальной. Следовательно, чувствительность можно повысить, увеличив высоту цилиндра до оптимальной.
Для пьезоцилиндра с меньшей толщиной стенки почти равномерная частотная характеристика чувствительности излучения получается при большей высоте (рисунок, б). На рисунке, б видно, что плавный ход кривой правее радиального резонанса искажают небольшие дополнительные экстремумы. Они связаны с проявлением изгибной по высоте цилиндра моды. При увеличении высоты цилиндра частота изгибной моды приближается к частоте радиального резонанса. У изгибной моды наблюдаются две узловые окружности колебательной скорости по высоте пьезоцилиндра. При этом на концах пьезоцилиндра колебательная скорость достигает большой величины. Это приведет к большим механическим напряжениям в пьезоцилиндре в местах узловых окружностей, а следовательно, может привести к разрушению пьезокерамического цилиндра при подаче на него большого электрического напряжения. По этой причине нужно избегать попадания паразитной из-гибной моды в рабочую полосу частот пьезопреобразователя. Этого можно достичь увеличением стенки пьезоцилиндра.
Выводы
Разработан модифицированный комбинированный метод конечных и граничных элементов для компьютерного моделирования и анализа гидроакустических пьезопреобразо-вателей произвольной формы и размеров. Для учета влияния жидкой акустической нагрузки введена матрица импеданса излучения. Введение импедансной матрицы позволило записать окончательную систему уравнений модифицированного комбинированного метода конечных и граничных элементов с учетом акустического излучения в виде одного матричного уравнения большой размерности. Выполнен анализ основных параметров цилиндрического водозаполненного пьезопреобразователя с радиальной поляризацией; определены его оптимальные геометрические размеры. Показано, что характеристики такого пье-зопреобразователя может искажать паразитная изгибная мода, которая может привести к разрушению пьезокерамического цилиндра при подаче на него большого электрического напряжения.
В заключение отметим, что пьезоцилиндр с внутренним твердым заполнением обладает характеристиками, подобными характеристикам водозаполненного цилиндра. При заполнении пьезоцилиндра легким материалом с небольшим акустическим сопротивлением, близким к акустическому сопротивлению воды (эбонит, сферопластик и т. п.), он также обладает широкополосной частотной характеристикой с двумя максимумами. Это показано авторами в работах [11, 12]. Таким образом, пьезоцилиндры с внутренним твердым заполнением, выполненным из легких материалов, также можно использовать в широкополосных гидроакустических антеннах.
Список литературы
1. Завадский В.Ю. Моделирование волновых процессов. М.: Наука, 1991. 248 с.
2. Smith R.R., Hunt J.T., Barach D. Finite element analysis of acoustically radiating structures with applications to sonar transducers // Journ. Acoust. Soc. Amer. 1973. Vol. 54, № 5. P.1277-1288.
3. Балабаев С.М., Ивина Н.Ф. Анализ пьезопреобразователей комбинированным методом конечных и граничных элементов // Акуст. журн. 1996. Т. 42, № 2. С. 172-178.
4. Балабаев С.М., Ивина Н.Ф. Компьютерное моделирование и анализ собственных колебаний пьезопреобразователей методом конечных элементов. Владивосток: Дальрыб-втуз, 2007. 242 с.
5. Шендеров Е.Л. Излучение и рассеяние звука. Л.: Судостроение, 1989. 304 с.
6. McMahon G.W. Performance of open ferroelectric ceramic cylinders in underwater transducers // Journ. Acoust. Soc. Amer. 1964. Vol. 36, № 3. P.528-533.
7. Вовк И.В. Излучение звука конечной решеткой, состоящей из открытых цилиндрических пьезокерамических оболочек // Акуст. журн. 1992. Т. 38, № 3. С.427-434.
8. Rogers P.H. Mathematical model for a free-flooded piezoelectric cylinder transducer // Journ. Acoust. Soc. Amer. 1986. Vol. 80, № 1. P. 13-18.
9. Балабаев С.М., Ивина Н.Ф. Анализ собственных колебаний пьезокерамических цилиндров произвольных размеров // Прикл. механика. 1989. Т. 25, № 10. С. 37-41.
10. Пьезокерамические преобразователи / под ред. С.И. Пугачева. Л.: Судостроение, 1984. 256 с.
11. Балабаев С.М., Ивина Н.Ф. Акустическое излучение цилиндрического пьезопреоб-разователя с внутренним твердым заполнением // Акуст. журн. 1999. Т. 45, № 4. С. 450-453.
12. Balabaev S.M., Ivina N.F. Computer modeling and analysis of a cylindrical piezoce-ramic transducer with interior solid filling // Journ. Acoust. Soc. Amer. 1999. Vol. 106, № 4, Pt. 2. P. 2196.
Сведения об авторах: Балабаев Сергей Михайлович, доктор физико-математических наук, профессор, e-mail: [email protected];
Ивина Наталья Федоровна, доктор технических наук, доцент, e-mail: ivinanata@ yandex.ru.
