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

56 Вестник СамГУ — Естественнонаучная серия. 2007. №4(54).

УДК 539.3

СИММЕТРИЧНЫЕ СХЕМЫ КОНЕЧНО-ЭЛЕМЕНТНОГО АНАЛИЗА ПЬЕЗОЭЛЕКТРИЧЕСКИХ УСТРОЙСТВ С УЧЕТОМ ВНЕШНИХ ЭЛЕКТРИЧЕСКИХ ЦЕПЕЙ И АКУСТИЧЕСКИХ НАГРУЗОК

© 2007 А.В. Белоконь, А.В. Наседкин, А.С. Даниленко1

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

Введение

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

1 Белоконь Александр Владимирович (niimpm@math.rsu.ru), Наседкин Андрей Викторович (nasedkin@math.rsu.ru), Даниленко Александр Сергеевич (alex_serf@mail.ru), НИИ механики и прикладной математики им. Воровича И.И., кафедра математического моделирования Южного федерального университета, 344090, Россия, г. Ростов-на-Дону, пр. Стачки, 200/1.

элементов ансамблируются с используемыми в ANSYS локальными матрицами диэлектрических проницаемостей по методологии, описанной в [3, 4].

Разработанный на кафедре математического моделирования ЮФУ (РГУ) пакет ACELAN [5, 6] до недавнего времени не имел возможностей расчета пьезоустройств, включенных во внешнюю цепь, и в настоящей работе предлагаются подходы, позволяющие присоединить к связанным конечно-элементным уравнениям пьезоэлектричества и акустики ACELAN элементы внешних электрических цепей. Как описывалось в [7], пакет ACELAN поддерживает технологию формирования всех конечно-элементных систем, решателей систем линейных алгебраических уравнений (СЛАУ), обобщенных задач на собственные значения, процедур интегрирования по времени для нестационарных задач с симметричными матрицами седловой структуры, что отличает принятую технологию конечно-элементного пьезоэлектрического анализа в ACELAN от используемой в ANSYS. Данная структура симметричных матриц и алгоритмов поддерживается и в настоящей работе при добавлении базовых элементов электрических цепей.

1. Континуальные постановки задач акустоэлектроупругости

Рассмотрим некоторый пьезопреобразователь П, представленный набором областей Пj = Прк; к = 1,2,...,Ир; у = к со свойствами пьезоэлектрических матриалов и набором областей Пу = Пет; т = 1,2,..., Ие; j = Ир + т со свойствами упругих материалов. Будем считать, что физико-механические процессы, происходящие в средах Прк и Пет, можно адекватно описать в рамках теорий пьезоэлектричества (электроупругости) и упругости.

Для пьезоэлектрических сред Пу = Прк предположим, что выполняются следующие полевые уравнения и определяющие соотношения:

где аау, вау, — неотрицательные коэффициенты демпфирования, а

остальные обозначения стандартны для теории электроупругости [8], за исключением дополнительного индекса ”у ”, указывающего на принадлежность к среде Пу с номером у.

Для сред Пу = Пет с чисто упругими свойствами будем учитывать только механические поля, для которых примем аналогичные (1)—(4) полевые уравнения и определяющие соотношения в пренебрежении электрическими полями и эффектами пьезоэлектрической связности

руи + аауруи -V- о = у V- Б = 0,

(1)

(2)

(3)

(4)

о = сЕ ••(£ + вау*) - с] ■ Е,

Б + Б = еу ■,(£ + £) + еу ■ Е>

£ = (Vu + Vuг )/2; Е = ^ф,

руи + аууи -V- о = fу,

(5)

о = Су • <е + |3^у£); £ = (Уи + Уит)/2. (6)

Наконец, пьезоэлектрическое устройство может быть нагружено на рабочие акустические среды Пу = Па/; I = 1,2,...,Иа] у = I + Ир + Ие. Для этих областей будем использовать уравнения акустики с учетом линейных диссипативных эффектов

——¿р + V • V; V = У\|/, (7)

р Ус2

р у\ = У-о; о = - р I + ЬуУу, (8)

где ру — равновесное значение плотности; су — скорость звука; Ьу—диссипативный коэффициент для среды Пу = Па1; р — звуковое давление; V — вектор скорости; у — потенциал скоростей; о — тензор напряжений; I — единичный тензор.

Отметим, что принятые здесь модели (1)—(8) отличаются специфическими способами учета демпфирования. Для упругих сред в (5), (6) имеем хорошо известную модель учета демпфирования по Релею. Эта модель в

(1)—(4) обобщена на пьезоэлектрические среды. Подробное обсуждение модели (1)—(4) содержится в [6, 9, 10], а обсуждение модели (7), (8) для акустической среды с диссипацией — в [11].

К (1)—(8) необходимо добавить соответствующие главные и естественные граничные условия, условия согласования полей по границам контакта различных сред, импедансные условия для ”усечения” акустических областей, а также начальные условия для нестационарных задач [11].

2. Аппроксимации МКЭ

Для решения динамических задач акустоэлектроупругости будем использовать МКЭ в классической лагранжевой формулировке.

Выберем согласованную конечно-элементную сетку, задаваемую в областях Пу аппроксимирующих области Пу. На этой сетке неизвестные полевые функции и, ф и у аппроксимируем в форме

и * ^(х) ■ т ф * N(х) ■ ф(о, у * ^(х) ■ т (9)

где ^ — матрица функций формы для поля перемещений и; ^^ — вектор-строки функций формы для полей электрического потенциала ф и потенциала скоростей в акустической среде у; и(г), Ф(0, ¥(0 — глобальные векторы соответствующих узловых степеней свободы.

Аппроксимация МКЭ (9) обобщенных постановок динамических задач

(1)—(8), включающих основные главные и естественные граничные условия, приводит к следующей системе дифференциальных уравнений:

М ■ а + С ■ а + К ■ а = Р, (10)

( Мии 0 п N ну Снн 0 1? \ ну

М= 0 0 0 , С = ^К1ф 0 0

Й т \ ^ну 0 —Муу , Йт ну 0 - С уу

К ии Киф 0

К= Кт ^ иф —Кфф 0

0 0 —Куу ,

Е = №и, Еф + КйЕ ф 0]т, а = [и,

(11)

относительно вектора неизвестных а. Здесь Смм = £у(айуМииу + |3^/Кмму), где Мииу и Кииу — структурные конечно-элементные матрицы масс и жесткости соответственно, а остальные входящие в (10), (11) подматрицы описаны в

[11]. Заметим, что вектор Рф является здесь вектором узловых электрических зарядов, взятых с обратным знаком.

В случае нестационарных задач к системе (10), (11) следует присовокупить начальные условия а(0) = ао; ¿1 (0) = ао, получающиеся из (9) и соответствующих континуальных начальных условий.

т

3. Моделирование элементов внешних электрических цепей

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

Резистивный элемент с номером т, узлами г, у и с сопротивлением Ят описывается законом Ома

Фг - Фу = ЯтЬ], (12)

где Фг, Фу — значения электрического потенциала в узлах г, у; 1гу — ток, текущий от узла г к узлу у по соответствующей ветви. Поскольку

1гу = й гу, (13)

где йгу — заряд в узле г, перетекающий в узел у, то можно записать уравнение для резистивного элемента с номером т в виде

КЯт . фет = (дет, (14)

чг=£(-1 "!)= •» ! ! «■={§!}• <•»

причем йуг = -йгу.

Для емкостного элемента с номером т с электрической емкостью Ст элементное уравнение можно представить в форме

кст • Фет = дст, (16)

а для индуктивного элемента с номером т и с индуктивностью Ьт — в

виде

V ^ - VI

где в (16), (17) элементные матрицы даются формулами

К1т • Фет = (3!т, (17)

к? = с.| -М: кг = ^(_| I- (И)

Конечно-элементные уравнения электрической цепи удобно формировать с использованием топологических понятий теории цепей и первого закона Кирхгофа для токов [12]. Согласно этому закону, в электрической цепи для каждого его узла г алгебраическая сумма токов всех ветвей, присоединенных к узлу, равна нулю: £кI-к = 0. Из закона Кирхгофа следует, что имеют место уравнения

^ О-к = О- = сош^-, (19)

к

где О- — заданный заряд в узле с номером г, определяемый из начальных условий.

Уравнение (19) может служить основой для ансамблирования элементных уравнений (14), (16), (17) и для их объединения с конечно-элементными уравнениями пьезоустройств (10), (11). Основная проблема состоит в сборке глобальных уравнений по единым по физическому смыслу правым частям, в частности, по сборке уравнений по электрическим зарядам. Как видно из (10), (11), (14), (16), (17), ситуация здесь осложняется тем, что в правые части КЭ уравнений входят как значения электрических зарядов, так и их производные различного порядка.

4. Симметричные формы разрешающих уравнений

Для задач об установившихся колебаниях, когда Ри = Риехр (гшг); Рф = Рфехр(гшг); а = аехр(гшг), из (10), (11) легко получить СЛАУ относи-

- ф

тельно вектора амплитуд а

К • а = Рс; Рс = |_Ри, Рф, 0]Т

с симметричной матрицей

К *^иис Киф Киус

= Кг ^иф К — ффс 0

КТ У ^иус 0 —Куус ,

(20)

(21)

где

Ки¥с = — ш2Йиш + гшК,

п = и, у, 1

иу

'■иу>

Кффе “ (1 + ^)к-

Отметим, что конечно-элементные уравнения для элементов внешней электрической цепи (14), (16), (17) можно записать в виде

(22)

Кет

Я '

фк"‘ = — (2Ят;

1

Кет

с

фк"‘ = —(3 ст;

Кет

ь '

ф*,п = —(3ьт,

-шгет _ тт^ет. к* - •

Кет _ т^ет.

с = Кс ;

(23)

Поскольку в (20), (22) в Рф, -ОЬ™ (Ь = С, Я, L) фигурируют узловые электрические заряды с обратным знаком, то проблем со сборкой глобальных конечно-элементных уравнений пьезоустройств со внешними цепями на основе закона (19) с случае установившихся колебаний не возникает.

В случае нестационарных задач непосредственно собрать КЭ уравнения

(10), (14), (16), (17) не удается. Заметим, что для численного интегрирования по времени уравнений движения в МКЭ обычно используются некоторые разностные схемы с весами. В ANSYS и в ACELAN применяются схемы Ньюмарка, однако, в существенно различных формулировках.

Для нестационарных задач акустоэлектроупругости, как было принято в ACELAN и для задач электроупругости [6, 7], используем метод Ньюмарка в альтернативной формулировке [13]. Для этого на V j-м временном слое tj = jx (т = At — шаг по времени) к (10), (11) применим осредняющие операторы Yj, согласно формулам

2

Yjb = ^ p*bj+i-k, b = a, FM, F9, (24)

k=0

Yjb = (yaj+i - (2y - 1)bj - (1 - Y)bj-i)/T, Yjb = (bj+i - 2bj + bj-i)/T2, (25)

где

во = в, Pi = Yi - 2в, P2 = Y2 + в, Yi = i/2 + Y, Y2 = i/2 - Y- (26)

Здесь в, Y — параметры методы Ньюмарка, безусловно устойчивого при в ^ (i/2 + y)2/4; y ^ i/2, и не обладающего аппроксимизационной вязкостью при в ^ i/4; y = i/2.

Соотношения (25) можно также представить в форме

Y i

Ур = рУр - УМ1) + УМ', Yjb = ^(Yjb - Y-У), (27)

Yjbp = Yibj + Y2bj-i ; Yjbp = (bj - bj-0/т, (28)

позволяющей для всех осредняющих операторов выделить единое слагаемое

Yja со значениями на j-м временном слое.

При Zd = 0 в (3) возможна достаточно простая процедура ансамблиро-вания осредненных по времени КЭ уравнений пьезоустройств и элементов внешней цепи. Именно, продифференцировав нужное число раз по времени, уравнение (10) для Ф и уравнения (14), (16) вместе с (17) можно записать в виде

К:ф • tj - Кфф • Ф = F?, (29)

-K^m • фет = -QCm, -KRm • Фem = -QRm, -К^т • Фет = -Qm (30)

Теперь, применяя осредняющие операторы (24), (25) к (10), (11) с заменой второго векторного уравнения на (29) и к (30), можно объединить (ансамблировать) полученные уравнения по однотипным правым частям, и при этом внутренние узловые электрические заряды в правых частях сократятся, а останутся только активные узловые источники электрических зарядов.

Полученную систему уравнений можно записать в форме, подходящей

для ACELAN, с симметричной матрицей блочной седловой структуры

Klf ■ Uj+i + Кнф • Ф/+1 + К{/ ■ ^.,+| = Kff aj, aj_!J, (31)

е//(т2

(32)

(33)

(34)

(35)

(36)

где

Подчеркнем, что при ансамблировании уравнений (31)—(36) внутренние узловые заряды в т27уЁф, т2У^Ьп (Ь = С, Я, Ь) сокращаются, и в правых частях остаются только известные величины.

Недостатком данного подхода является ограничение ^ = 0, а также дифференцирование по времени конечно-элементных уравнений (правда, только для применения осредняющих операторов), что может однако повлиять на точность аппроксимаций.

Альтернативный подход состоит в непосредственном применении осредняющих по времени операторов в форме (24), (27), (28) к уравнениям (10),

(11), (14), (16), (17) и в выделении в правых частях слагаемых Рф, Qbш по (27). При ансамблировании полученных уравнений будут сокращаться внутренние заряды в форме — YjРф, Qb'n, однако внутренние электрические заряды из -Еф^ —Рфу_1), ОЬп, QЬnj-l), получаемые из (28), останутся. Для их определения можно привлечь дополнительные аппроксимации уравнений (10), (11), (14), (16), (17) на временных слоях j и j — 1.

Фактически именно близкий по идеологии подход, описанный в [4] при использовании метода Ньюмарка в совсем другой формулировке, принят в

квадратного корня [7], что позволяет на каждом временном слое решать только СЛАУ с нижними и верхними треугольными матрицами.

Отметим, что такие важнейшие процедуры, как повороты узловых степеней свободы и учет главных граничных условий, необходимые для формирования систем (20), (31)—(36), также могут быть реализованы в симметричных формах без нарушения структур матриц МКЭ [7].

Реализация описанных выше симметричных алгоритмов решения задач акустоэлектроупругости осуществлена в последних версиях АСЕЬА^ Были проведены тестовые расчеты ряда гидроакустических пьезопреобразователей и пьезоустройств, включенных во внешние цепи, со сравнением результатов с аналогичными, полученными по конечно-элементной программе

ANSYS [3].

Эффективная матрица жесткости Ке// ансамблированной системы (31)—(36) имеет блочную седловую структуру, т.к. КЦ/ > 0; Кфф > 0; > 0; Ке// = Ке//Т. Эта матрица может быть факторизована по методу

ANSYS. Например, конечно-элементный расчет реального пьезоэлектрического устройства в виде многослойного пьезотрансформатора Розена, включенного в развитую электрическую цепь, был осуществлен в [14] и показал эффективность предлагаемой формы конечно-элементного пьезоэлектрического анализа с элементами внешней цепи.

Работа выполнена при частичной поддержке РФФИ (грант 05-01-00752).

Литература

[1] Данов, Г.А. Пьезоэлектрические трансформаторы / Г.А. Данов. - М., 2003. - 320 с.

[2] Клигман, Е.П. Динамические характеристики тонкостенных электро-упругих систем / Е.П. Клигман, В.П. Матвеенко, Н.А. Юрлова // Изв. РАН. МТТ. - 2005. - №2. - C. 179-187.

[3] ANSYS 10.0, Theory Reference. ANSYS Inc., 2005 © SAS IP, Inc.

[4] Wang, J.S. A finite element-electric circuit coupled simulation method for piezoelectric transducer / J.S. Wang, D.F. Ostergaard // Proc. IEEE Ultrasonics Symp. - 1999. - P. 1105-1108.

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

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

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

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

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

[10] Наседкин, А.В. Особенности учета демпфирования в конечно-элементном пьезоэлектрическом анализе / А.В. Наседкин // Материалы Межд.

научн.-практич. конф. ’’Фундамент. проблемы пьезоэлектрич. приборостроения” (”Пьезотехника-2000”), Москва, 27 ноября-1 декабря 2000 г. -М., 2000. - С. 154-158.

[11] Наседкин, А.В. К расчету по МКЭ пьезопреобразователей, нагруженных на акустическую среду / А.В. Наседкин // Известия ВУЗов. Северо-Кавказский регион. Естеств. науки. - 1999. - №1. - С. 48-51.

[12] Ушаков, В.Н. Электротехника и электроника / В.Н. Ушаков. - М.: Радио и связь, 1997. - 328 с.

[13] Зенкевич, О. Конечные элементы и аппроксимация / О. Зенкевич, К. Морган. - М.: Мир, 1986. - 318 с.

[14] Наседкин А.В. Конечно-элементный расчет многослойного пьезотрансформатора Розена с двумя выходами / А.В. Наседкин // Мат. моделирование, биомеханика и информационные технологии в современном университете: тез. докл. Всерос. шк.-семинара, Ростов-на-Дону: 22-26 мая 2006 г.. - Ростов-н/Д: Тера Принт, 2006. - С. 36-37.

Поступила в редакцию 15/У/2007; в окончательном варианте — 15/У/2007.

SYMMETRICAL SCHEMES OF THE FINITE-ELEMENT ANALYSIS FOR PIEZOELECTRIC DEVICES WITH EXTERNAL ELECTRIC CIRCUITS AND ACOUSTIC

LOADS

© 2007 A.V. Belokon, A.V. Nasedkin, A.S. Danilenko2

The finite element technologies for piezoelectric devices modelling with electric circuits and acoustic load are considered. For finite element package ACELAN the basic finite elements of electric circuits and symmetric forms of finite element matrices and algorithms are constructed. The harmonic and transient operating regimes of piezoelectric devices are considered.

Paper received 15/V/2007. Paper accepted 15/V/2007.

2Belokon Alexandr Vladimirovich (niimpm@math.rsu.ru), Nasedkin Andrey Victorovich (nasedkin@math.rsu.ru), Danilenko Aleksey Sergeevich (alex_serf@mail.ru), Dept. of Mathematical Modelling, Research Institute of Mechanics and Applied Mathematics, Southern Federal University, Rostov-on-Don, 344090, Russia.