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

ПРИБОРОСТРОЕНИЕ. ИНФОРМАТИКА

УДК 621.937

ПРИМЕНЕНИЕ МЕТОДА ПЕРЕДАТОЧНЫХ МАТРИЦ ДЛЯ ИССЛЕДОВАНИЯ ГИБКИХ УЛЬТРАЗВУКОВЫХ ВОЛНОВОДОВ

Канд. техн. наук СТЕПАНЕНКО Д. А., канд. техн. наук, доц. МИНЧЕНЯ В. Т.

Белорусский национальный технический университет

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

Теоретические основы методики расчета.

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

икТ{х)а'У + - а) + = 0;

{ С^бЗДт!')' - С^бЗДа)' + ю2р£(х)т1 = 0,

где Е - модуль упругости материала волновода; J(x) = кс?(х)/64 - осевой момент инерции поперечного сечения волновода; d(x) - диаметр волновода; т](х) - амплитуда поперечных смещений; р - плотность материала волновода; ш = = 2л/- циклическая частота колебаний; /- частота колебаний; S(x) = 0,25лй?2(х) - площадь поперечного сечения волновода; а(х) - амплитуда угловых смещений поперечного сечения волновода; Ks - коэффициент формы поперечного сечения волновода; G - модуль сдвиговой упругости материала волновода.

Для того чтобы система уравнений (1) имела однозначно определенное решение, необходимо задать четыре граничных условия, например значения а(х), т](х), а'(х), т]' (х) во входном сечении волновода х = 0 либо четырех линейно независимых линейных комбинаций этих величин. Если обозначить эти комбинации через у\(х), у2(х), Уз(х), у4(х), то согласно методу передаточных матриц выполняется соотношение

у(х) = Т(х)у(0), (2)

где у(х) = (у,(х) у2(х) у3(х) Ул(х))Т] Т(х) -передаточная матрица.

Докажем соотношение (2) для: у\(х) =

= а(х); у2(х) = т1(х); уъ(х) = а'(х); у4(х) = ц'(х). В этом случае система уравнений (1) может быть представлена в виде эквивалентной системы дифференциальных уравнений первого порядка

У'(*) =

О 0 10

О 0 0 1

aifi(x)-a2 0 -fi(x) -ûi/i(x)

0,5 f2(x) -a2/ai 1 -0,5 f2(x))

у(х) = А(х)у(х),

(3)

где Мх) = 16/d*(x); f2(x) = 4d'(x)/d(x) ; 0l =Kt

Рассмотрим вектор-функции вида =

fj\ t t:\ . - . = yv-'(.ic»/vv'mi». где i = 1, ..., 4, a вектор-

функции уЮ(х) = (й®(х) ^«f

представляют собой частные решения системы уравнений (3) при граничных условиях у(0) =

= (5,i 5,2 5,з 5,4 ) , где 5,у - символ Кроне-

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

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

т(х) = det

Л*)

У(г\х) У$\х) У$\х)

у?\х) У?\х) У?\х) У?\х)

У?\х)

У4\х)

У?\х) У?\х)

у?Чх)

У?\х)

АЧх)

Для вронскиана выполняется тождество Якоби - Лиувилля [8]

W(x) = W(Q) exp || Sp( А(х))й6с | :

= W(0) exp |-1,5 j f2 (x)dx | = W(G)\

m

d(x)

\6

где Бр обозначает след матрицы (сумму элементов главной диагонали).

Из тождества Якоби - Лиувилля следует, что вронскиан либо тождественно равен нулю, либо всюду отличен от нуля. В первом случае система вектор-функций будет линейно зависимой, а во втором — линейно независимой. В силу граничных условий, принятых при определении вектор-функций у('Нх). вронскиан

в точке х = 0 представляет собой определитель единичной матрицы и выполняется равенство

£?/£; а2 = ш2р/£ .

Щ0) = 1. Таким образом, вронскиан будет всюду отличен от нуля, а система вектор-функций у(')(х) - линейно независимой.

Так как вектор-функции у^<Х> пропорциональны вектор-функциям у&{х). они также линейно независимы, и частное решение у(х) системы уравнений (3) при произвольных граничных условиях у, (0), у2(0), у3(0), у4(0) можно представить в виде разложения

у (д;) = (0)у« (х) + (0)у (2) (л:) + +>"з (0)У(3) (*) + У л (0)у(4) (х), которое можно представить в матричной форме у(х) = [у0)(х) у(2)(х) уЩх) У(4) (*)]у(0) = = Т(х)у(0). (4)

Уравнение (4) по своей форме аналогично (2), что доказывает справедливость последнего. Из (4) также следует, что столбцы передаточной матрицы представляют собой в рассматриваемом случае определенные выше вектор-функции у(0 (х).

Описание методики расчета. Рассмотрим далее использование метода передаточных матриц для проектирования комбинированного гибкого волновода, конструкция которого приведена на рис. 1.

Рис. 1

Волновод состоит из трех участков: ступеней 1 и 2 цилиндрической формы, имеющих постоянную площадь поперечного сечения, и переходного участка 3 с переменной площадью поперечного сечения, служащего для плав-

ного сопряжения ступеней, снижающего концентрацию напряжений. Общая передаточная матрица Т(Х) комбинированного волновода

может быть определена как произведение передаточных матриц отдельных участков

Т(Х) = Т2(^)Тг(Л1)ВД), (5)

где и Х2 - длины ступеней волновода; АЬ -длина переходного участка, Ь = Ь1+Ь2 + АЬ -общая длина волновода; Т! и Т2 - передаточные матрицы первой и второй ступеней; Тг -передаточная матрица переходного участка.

В дальнейшем будет рассматриваться задача определения длин ступеней волновода, обеспечивающих его резонанс на заданной частоте / (задача синтеза волновода с заданной резонансной частотой), в связи с чем в (5) можно опустить аргументы, являющиеся постоянными величинами, и представить его в виде

Т(Х1,Х2) = Т2(Х2)ТД1(Х1). (6)

Граничные условия на концах волновода имеют вид:

а(0) = 0; г|(0) = 0; а'(0) = 0;

т!'(1)-а(Х) = 0, (7)

в связи с чем представляется рациональным задать элементы вектор-функции у(х) в (2) в виде:

Зл(лг) = сс(х); у2(х) = ц(х); уз(х) = а'(х);

У4(х) = ц'(х)-а(х).

Для определения передаточной матрицы волновода с постоянной площадью поперечного сечения (передаточных матриц ступеней комбинированного волновода) запишем аналитические решения (1):

Т](х) = Ху 8Ш(К]Х) + Х2 С08(К]Х) +

+Х38Ь(к2х) + ХАс\\{к2х); а(х) = Ух 8т(К]Х) + У2 ««(к^х) +

а(0) Л(0) а'(0) Л'(0)-а(0)

К1

р<в2 Кгв К! 0

0

рш2 к,0 К!

0 1

+}з8Ь(к2х) + У4сЬ(к2х), (8)

где Х1 и У1 (г= 1, ..., 4) - коэффициенты, зависящие от граничных условий; К[ и к2 - волновые числа, определяемые соотношением

К12 =

1

+

1 + -

со

2с2

- + „

\2

1 —

ю4 16ю2 4с4 сЮ2

(9)

где -О - диаметр волновода (диаметр ступени); с = у[Е/р - скорость продольной стержневой

ультразвуковой волны в материале волновода.

Рассмотрим значения элементов вектор-функции у(х) для входного сечения волновода х = 0:

^ 1 \ рш2

а(0) = У2 + У4 =

+

к2 +■

К1-

рш2

К,в К!

х1 +

К,в к2 т](0) = Хг+ Хг,

X,

а'(0) = к1У1+к2У3 =

-кf +

рш2

Х? +

+

щ +

рш2

Хл

л'(0) - а(0) = кЛ + к2Х3 -У2-У4 = рш2 рш2

-Х3.

(10)

К,в К! К,вк2

Здесь учтено, что между коэффициентами X,- и У{ существует линейная зависимость [7]. Действительно, так как (1) имеет второй порядок по производным функций а(х) и г](х), то решения (8) содержат четыре линейно независимых коэффициента, а остальные четыре коэффициента являются их линейными комбинациями. В матричной форме соотношения (10) могут быть записаны в виде

к2 +

-к( +

рш2

ад

рш2 0

К,вк2

0 1

0 к2 + р(°2

ад

рш2 0

'Х^

Х2

хг

Выполнив обращение матрицы, можно найти связь между коэффициентами Х^ и значениями элементов вектор-функции у(х) при х = 0:

х2 хг

\ХА у

К! 0

0 к?+-рЮ

К,в

к2 0

0 к?-

2 Р®

к«0

0 К! -1

0 к2

1

рш2

о

' 11Л ч Р®2

О

а(0) Л(0) а'(0) Л'(0)-а(0),

к^ + к2

-Т«

а(0) Л(0) а'(0) П'(0)-а(0)

(П)

Рассмотрим далее значения элементов вектор-функции у(х) для выходного сечения волноводах = Ь\ а(Ь) = 8т(к,1) + У2 С08(к,1) + Уз8Ь(к21) + У4с11(к2Х) =

(

к, -

рю

2 Л

(

Хх С08(к1Х) +

-К! +

рю

2 Л

Х2 8т(к1Х) +

+

рю

2 Л

к, +-

Х3сЬ(к2Ь)+

рю

2 Л

к, +•

Х48Ь(к2Х) ;

г|(Х) = Хх вц^к^) + Х2 совСк^) + Х38Ь(к2Х) + Х4сЬ(к2Ь);

а'(Ь) = со8(к1Х)-к1Г2 8т(к!Х) + к2Г3сЬ(к2Х) + к2Г48Ь(к2Х) =

Г 1 \

ч

К.О у

Г 2 \

ч

К.О у

Х2 сов( к^) +

' 2 ^

Х38Ь(к2Х) +

' 2 ^

Х4сЬ(к2Х);

Т|'(Х) — сх.(Х) = (к^ - Г2) сов( к^) - (щХ2 + ) вт( к^) +

+ (к2Х3 - У4)сЬ(к2/.) + (к2Х4 - У3)8Ь(к2/.) = Х] совСк^) -

рю

2 2

-X, втГ к,/,) - -Х,сЬ(к- -ХдвЬГк,/,) .

2 1 К,в к2 3 2 К,в к2 4 2

В матричной форме эти соотношения могут быть записаны в виде

а(1)

а \Ь) ч'Щ-аЩ

= Т(2)(1)

'Х^

х2

Х3 хч

или с учетом уравнения (11)

а(1) а \Ь)

Матрица

Т(£)

к^ + к2

Т(2)(1)Т0)

а(0) Л(0) а'(0) л'(0)-а(0).

= , 1 Т(2)(1)Т(1) (12)

2 2 щ+щ

представляет собой передаточную матрицу волновода с постоянным поперечным сечением (ступени комбинированного волновода).

Для определения передаточной матрицы rTt переходного участка необходимо в передаточной матрице, определенной выражением (4), к первому столбцу добавить четвертый столбец и в полученной матрице из четвертой строки вычесть первую строку. Методика вычисления передаточной матрицы, входящей в (4), была описана ранее [3, 4] и в связи с этим не рассматривается. Она основана на численном решении последовательности задач Коши для системы уравнений (3) с помощью метода Рун-ге-Кутта. Функция (1(х), описывающая профиль переходного участка, задается в виде многочлена.

С учетом выражения (12) общая передаточная матрица комбинированного волновода будет определяться соотношением

Т(Х1,Х2) = -

1

(к?+к2)(к2+к I)

где волновые числа К1 и к2 соответствуют первой ступени волновода и определяются выражением (9), если положить в нем £> = а волновые числа к3 и К4 соответствуют второй ступени волновода и определяются выражением (9), в котором Б = /)2; и Б2 - диаметры первой и второй ступеней волновода.

С учетом граничных условий (7) можно записать следующую систему уравнений относительно неизвестных величин а(Ь), т](£), а'(0), т]'(0):

' 0 ^

Л(^) 0 = Т(Ц,Ь2) 0 а'(0)

, 0 ,

или

-1 0 ТЬС^ь^г) Ти (1ц, 1а)

0 -1 ТгзС^ь^а) Т24{Ь1, Ьг)

0 0 а'(0)

0 0 ^иС^Ь^г), [чХО))

Га(Х)]

а'(0) = 0.

(13)

Разложив определитель матрицы Т'(Ь[,Ь2) по элементам первого столбца, можно свести его к определителю третьего порядка, разложив который по элементам второго столбца, получим равенство

Г34 (Ь\, Ь2 ) ^43 (ЬиЪ) Т44 Щ

Для того чтобы система (13) имела нетривиальное решение, ее определитель должен обращаться в нуль, что дает условие резонанса

(14)

В ранее предложенной авторами модели [3, 4] условие резонанса представлялось в виде

определителя шестого порядка. Неизвестными переменными являлись коэффициенты Хц, Х2, Х5, ..., Х8, первые два из которых являются коэффициентами аналитического решения вида (8) для амплитуды поперечных смещений первой ступени волновода, а коэффициентыХ$, ...,Х8 входят в аналитическое решение для амплитуды поперечных смещений второй ступени волновода. Коэффициенты Х3 и как следует из граничных условий (7), пропорциональны коэффициентам Хц и Х2. Применение метода передаточных матриц позволило представить условие резонанса в виде определителя второго порядка, что дает ряд преимуществ:

1) упрощается вычисление определителя и снижаются возникающие при этом погрешности;

2) для системы уравнений (13) заранее известно, какие уравнения системы являются линейно зависимыми (это уравнения относительно неизвестных г|'(0) и а'(0)). Это позволяет задать одну из переменных, например Т1'(0), произвольным образом, выразить через нее переменную а'(0), а остальные переменные (т](£) и а(Ь)) выразить в виде линейных комбинаций переменных г|'(0) и а'(0). В ранее предложенной модели [3, 4] неизвестен ранг матрицы Аг, входящей в условие резонанса. Вследствие этого неизвестно, какие из коэффициентов Хи Х2, Х5, ..., Х8 могут быть заданы произвольным образом. Численное определение ранга матрицы Аг затрудняется наличием в ней быстрорастущих экспоненциальных членов, в результате чего определитель матрицы и определители подматриц более низкого порядка сильно отклоняются от нуля даже в точках весьма близких к резонансу. Это затрудняет проверку условия обращения в нуль определителей всех подматриц с порядком выше чем гаг^А,.. Отсутствие информации о ранге матрицы Аг и наличии линейной зависимости между уравнениями системы АД = О (X = {Хх Х2 Х5 Хь Х7 Х8)Т) делает невозможным решение этой системы относительно коэффициентов Хъ Х2, Х5, ..., Х8 и определение собственных форм колебаний, которые связаны с этими коэффициентами уравнениями (8).

Определение недостающих граничных условий г|'(0), а'(0) для входного сечения волно-

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

Численные результаты и их обсуждение. Численная реализация описанной выше методики с помощью программы МаИгСас! позволила построить резонансные кривые изгибных колебаний для частоты /=25 кГц. Эти кривые представляют собой совокупность точек на плоскости параметров (Ь\, Ь2), для которых выполняется условие резонанса (14). Полученные резонансные кривые с высокой точностью совпадают с кривыми, построенными на основании ранее разработанной авторами модели [3-6], и в связи с этим не приводятся. В качестве примера результатов расчета собственных форм колебаний на рис. 2 приведен график зависимости амплитуды поперечных смещений от координаты х для волновода с длинами ступеней Ь\ = 0,0895 ми12 = 0,0209 м, который имеет на частоте /=25 кГц резонанс на первой моде продольных колебаний и 11-й моде изгибных колебаний.

1,00 0,82 0,64 6 0,46 * 0,28 I 0,10 | -0,08 1 -0,26 -0,44 -0,62

'О 0,012 0,036 0,060 0,084 0,110

дс, м

Рис. 2

Сплошная линия соответствует результатам расчета с использованием метода передаточных матриц, а точки - результатам конечно-элементного моделирования с помощью программы АКБУБ. Вертикальные линии х = 0,09 м и х = 0,096 м соответствуют границам переходного участка волновода. Методика моделирования с использованием программы АКБУБ была описана [3, 4, 9] и в связи с этим не рассматривается. Как видно из приведенного графика, результаты расчета с использованием

обоих методов совпадают с высокой точностью (среднеквадратическое отклонение около 0,009), что подтверждает корректность и эффективность предложенной методики. Как отмечалось выше, ранее разработанная авторами модель не позволяла производить расчет собственных форм колебаний.

выводы

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

2. На конкретном численном примере путем сравнения с результатами конечно-элементного моделирования с использованием программы ANSYS показана корректность и эффективность предложенной методики расчета.

ЛИТЕРАТУРА

1. Zhou, G. The performance and design of ultrasonic vibration system for flexural mode / G. Zhou // Ultrasonics. -2000. - Vol. 38. - P. 979-984.

2. Квашнин, С. E. Ультразвуковые электроакустические преобразователи и волноводы-инструменты для медицины / С. Е. Квашнин. -М.: Изд-во МГТУ, 1995. -43 с.

3. Stepanenko, D. A. Modeling of flexible waveguides for ultrasonic vibrations transmission: Longitudinal and flexural vibrations of non-deformed waveguide / D. A. Stepa-panenko, V. T. Minchenya // Ultrasonics. - 2010. - Vol. 50. -P. 424—430.

4. Минченя, В. Т. Линейные колебания двухступенчатого волновода-концентратора для ультразвукового тромболизиса / В. Т. Минченя, Д. А. Степаненко // Доклады НАН Беларуси. - 2009. - Т. 53, № 6. - С. 114-119.

5. Стенаненко, Д. А. Собственные колебания ультразвуковых волноводов для минимально-инвазивной хирургии / Д. А. Степаненко, В. Т. Минченя, А. В. Чигарев // Теоретическая и прикладная механика. — 2010. — Вып. 25. — С. 276-281.

6. Минченя, В. Т. Анализ резонансных явлений в гибких волноводных системах с применением теорий Эйлера-Бернулли и Тимошенко / В. Т. Минченя, Д. А. Степаненко, Е. Н. Юрчик // Приборостроение-2009: материалы МНТК. - Минск, 2009. - С. 247-248.

7. Han, S. M. Dynamics of transversely vibrating beams using four engineering theories / S. M. Han, H. Benaroya, T. Wei // Journal of Sound and Vibration. - 1999. - Vol. 225. -P. 935-988.

8. Ланкастер, П. Теория матриц / П. Ланкастер. - М.: Наука, 1982.-272 с.

9. Bubulis, A. Semi-automatic modal analysis of flexible ultrasonic waveguides in ANSYS / A. Bubulis, V. T. Minchenya, D. A. Stepanenko // Приборостроение-2009: материалы МНТК. - Минск, 2009. - С. 145-146.

Поступила 30.03.2010