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

УДК 621.937

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

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

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

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

Объектом исследования является волновод, состоящий из двух цилиндрических участков (ступеней) с длинами Ь и ¿2, имеющих постоянную площадь поперечного сечения и связанных переходным участком длиной АЬ с плавно изменяющейся площадью поперечного сечения (рис. 1).

Форма переходного участка установлена опытно-экспериментальным путем [5] и описывается полиномиальной функцией [1, 2]. Длины ступеней Ь1 и Ь2 выбираются таким образом, чтобы обеспечить резонанс на продольной моде колебаний первого порядка при частоте / = 25 кГц, для чего могут быть использованы

[ СКшОБ (х) + Р( х, t)) дХ ] + р£ (х)

Б

расчетные резонансные кривые [1, 2]. Длины Ь\ и Ь2 связаны однозначной зависимостью, что позволяет при построении диаграммы устойчивости задавать только одну их них. Для каждой из резонансных геометрических конфигураций волновода, характеризуемых длиной первой ступени Ьь необходимо определить критические амплитуды сг (¿1) продольных

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

Рис. 1

Сведение задачи к уравнению Матьё. Как

показано ранее [1, 2], наиболее точным является моделирование изгибных колебаний гибких волноводов на основе теории Тимошенко, согласно которой изгибные колебания волновода переменного сечения под действием осевой вынуждающей силы описываются уравнением

Б (п( х, t) а (х, t) )Т = 0,

(1)

где п(х, t) - поперечное смещение; а(х, t) -

угловое смещение поперечного сечения волно-

еет

-К^ + £ (х)|-

ах дх

\

К&Б (х)|

дх

д_ дх

(

(

3 (х)

Е -

Р( х, t) ^

£(х) ) дх

\

-КшОБ (х) + р (х) ^

дt2

где Ко - коэффициент формы поперечного сечения; О - модуль сдвиговой упругости, Б(х) - площадь поперечного сечения; Р(х, X) = = Р0 (х) с о - осевая вынуждающая сила; Р0(х) - амплитуда вынуждающей силы; О -циклическая частота вынуждающей силы; р -плотность материала волновода; /(х) - осевой момент инерции поперечного сечения; Е - модуль продольной упругости.

Коэффициент формы для волновода с круглым поперечным сечением определим следующим образом:

= 6(1 + У) 5 7 + 6у ,

где V - коэффициент Пуассона материала волновода.

В соответствии с методом Бубнова-Галёр-кина будем искать приближенное решение (1) в виде разложения по собственным формам свободных изгибных колебаний волновода

(П (х, X) а(х, X))Г (X) Уш (х), (2)

т

где базисные вектор-функции

V т ( х ) = ( Ут\( х) Ут х) )Т

представляют собой собственные векторы обобщенной задачи на собственные значения

(Бо -ю2Б(х)) х) = 0, (3)

9т (X) - подлежащие определению функции времени.

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

Оператор Б0 и матрица Б(х) в уравнении (3) могут быть получены из оператора Б при условии Р(х, X) = 0 (свободные колебания без вынуждающей силы) и гармоническом законе изменения решений (д2/дх2 ^-ю2 ) путем выделения членов, содержащих собственные частоты ю:

Б0

1

^ ( КоОБ{х)^ 1 -Ко О ( ^ + Б(х) ^ ^

йх

йх

йх йх

КО (х)— —\Е/(х)— |-КОБ (х) ах ах V ах

Б( х) =

( Б (х) 0 0 /(х)

Л

(4)

Оператор Б0 является эрмитовым

| иТ (х)Б0 х)йх =| VТ (х)Б0и( х)йх, (5) 0 0 где и(х), х) - любые вектор-функции, удовлетворяющие граничным условиям, описываю-

щим закрепление концов волновода; Ь - общая длина волновода.

Так как оператор Б0 является эрмитовым, а элементы матрицы Б(х) неотрицательны, то собственные векторы задачи (3) удовлетворяют обобщенному условию ортогональности [7, раздел 15.4.6]

Ь

| уШ (х)Б(х) V„ (х)йх = 8т„,

0

где 8т„ - символ Кронекера.

С учетом выражений (4) запишем это условие в виде

Ь

I (Б (х)у

ш1 (х)У„1 (х) + /(х)Уш2 (х)У„2 (х))йх = 0.

0

В соответствии с методом Бубнова-Галёр-кина функции 9ш (X) должны удовлетворять условию

Ь

I УТ (х)Б (г| (х, X) а (х, X) )т йх = 0. (6)

0

Ошибка аппроксимации, возникающая при замене точного решения (1) приближенным решением (2) и входящая в условие (6), определяется вектором

1 (бт +Ю2твт )Б(х)Уш:(х) + ^

Б

П (х, X) а (х, X)

= 1

+ бш соо Qx I р (х) от |

р йх V йх )

(б т ) / (х)Ут2 (х) -

0ш cos 0.x й

йх

Р)(х)/(х) йУт2

Б(х) йх

)

Подставляя ошибку аппроксимации в (6), получим с учетом условия ортогональности уравнение вида

Р

ё„ + < 0„ + cos Qt ^ c„mem = 0, (7)

m

коэффициенты которого имеют вид

|а( х) ij (х)

dVm2 dv„2 _s(х) dm dvn_ dx dx dx dx

dx

Cnm

pj (S (x)v2( x) + J (x)v22( x))dx

,(8)

где c(x) = P0(x)/S(x) - амплитуда осевых механических напряжений, возникающих при воздействии вынуждающей силы.

В случае одномодовых колебаний, соответствующем рассмотрению одного слагаемого в разложении (2), уравнение (7) принимает вид

Уп + (an _ 2q„ cos 2z)yn = 0, (9)

где введены новые переменные z = Qt/2; ап = = 4ю2/Q2; qn = _2Cnn/Q2; yn(z) = 0n(2z/Q); n - порядок изгибной моды колебаний, возникающей при потере устойчивости.

Уравнение (9) известно как уравнение Ма-тьё [3, 4, 8]. Устойчивость его решений определяется параметрами a и q и может быть представлена в виде областей устойчивости на плоскости параметров (a, q), графическое изображение которых известно как диаграмма Айнса-Стретта (Ince-Strutt) [9, с. 223] (рис. 2).

-35 -30 -25

Ф1

F

Fo

Fi

F2

Рис. 2

Границы областей устойчивости, выделенных штриховкой, соответствуют 2п-периоди-ческим решениям (9) и определяются характеристическим уравнением, записываемым с помощью цепных дробей. Для кривых ^2п+1 и Ф2п+1, п = 0,1, 2,... (верхние и нижние границы

нечетных областей неустойчивости) характеристическое уравнение имеет вид [8, с. 468]

±д +1 - а_... _

| 9 - а | 25 - а

д 2|

| (2r +1)2 _ a

_... = 0,

где использована сокращенная форма записи цепных дробей

ai

ai a2

bi +

- = b0 +—1 + —1 +... a2 | Ь | b2

b2 + ...

При д ^ 0 уравнение (9) принимает вид уравнения гармонических колебаний, которое имеет 2п-периодическое решение при условии а = п2, где п - натуральное число. Поэтому при д ^ 0 границы областей устойчивости определяются равенством О«2юп/п, т. е. при п = 1 (главная область неустойчивости) циклическая частота вынуждающей силы должна быть равна удвоенной резонансной циклической частоте изгибных колебаний. Более точное уравнение границ главной области неустойчивости при д ^ 0 можно получить, воспользовавшись характеристическим уравнением

О « 2юп^/1 ±д.

Определение параметров уравнения Ма-тьё с помощью метода конечных элементов.

Для определения параметров а и д в (9) необходимо для каждой из геометрических конфигураций волновода рассчитать резонансные частоты и формы изгибных колебаний, а также распределение амплитуды о(х) осевых механических напряжений по длине волновода. Как показано ранее [1, 2], резонансные частоты, рассчитанные на основе теории Тимошенко, с погрешностью не более 0,2 % совпадают с резонансными частотами, рассчитанными методом конечных элементов, в связи с чем для определения резонансных частот и форм колебаний использовалась программа А№У8. Для решения последовательности задач, соответствующих множеству геометрических кон-фигу-

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

создана программа на языке APDL (ANSYS Parametric Design Language). Расчет каждой геометрической конфигурации описывается с помощью цикла, включающего в себя следующие действия: построение геометрической модели и наложение граничных условий, разбиение модели на конечные элементы, модальный анализ (определение резонансных частот и форм колебаний), постобработку результатов с вычислением параметров a и q. Блок-схема программы приведена на рис. 3. Программа разработана на основе ранее созданной авторами программы для автоматизированного моделирования гибких ультразвуковых волноводов [10].

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

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

Рис. 3. Физико-механические параметры материала (Е, р, V), диаметры ступеней, геометрические параметры переходных участков; N - число исследуемых геометрических конфигураций (число точек дискретизации параметров Ь\ и Ь2); * - комментарий ко 2-му блоку схемы

Собственные частоты fn(L1), п = 1, 2, ..., 9 изгибных колебаний, найденные в результате модального анализа, использовали для расчета параметра a в соответствии с формулой an (L1) = 4 /n2(Lj)/ f2. При постобработке использовали интерполяцию результатов анализа на линии (пути), в качестве которых выбирались продольная ось волновода (путь P1, рис. 1) и его образующая (путь P2), в результате чего формировались путевые переменные (path variables), над которыми могут выполняться математические операции. При постобработке результатов анализа продольных колебаний использовался путь P1, а в качестве путевой переменной ux1 - осевая составляющая ux1( x) смещения точек этого пути, соответствующая амплитуде 2(x) продольных колебаний. Так как величины смещений определяются при модальном анализе в произвольном масштабе, необходимо задать масштабный множитель, например потребовав, чтобы амплитуда 2(0) продольных колебаний входного сечения была одинаковой для всех геометрических конфигураций (например, 2(0) = 20 = 1 мкм). Дифференцирование масштабированной переменной ux\ использовалось для определения амплитуды механических напряжений o(x) (путевой переменной stress) по формуле с(x) = —^L.

ux1 (0) dx

В отличие от an (L1) параметр qn помимо длины L1 будет также зависеть от амплитуды 2(0) колебаний входного сечения, т. е. qn = qn( L1,2,(0)). Так как параметр qn связан с амплитудой 2(0) линейной зависимостью, достаточно вычислить его значение qn0(L1) = qn (L1,20)для амплитуды 2(0) = 20. Также вычислялись значения коэффициента усиления продольных колебаний K(L1) = |ux1(L)/ux1(0)| (путевой переменной K).

При постобработке результатов анализа изгиб-ных колебаний использовались пути P1 и P2. В качестве первой путевой переменной uy1 использовалась перпендикулярная оси волновода составляющая uy1(x) смещения точек пути P1, соответствующая амплитуде n(x) изгибных ко-

лебаний (собственной функции vn1( x)), а в качестве второй переменной ux2 - составляющая ux 2( x) смещения точек пути P2, параллельная оси волновода, которая использовалась для определения амплитуды a(x) угловых смещений (собственной функции vn2( x), путевой переменной angle) в соответствии с формулой a(x) = arctg(2ux2(x)/d(x)) « 2ux2(x)/d(x), где d(x) - диаметр волновода. Четвертьволновая продольная мода исключалась из рассмотрения путем проверки условия | ux1(L)/uy1(L) | < 0,1,

которое выполняется лишь для изгибных колебаний. Вычисление параметра q производили путем численного интегрирования. Значения параметров an (L1) и qn0(L1) для всех девяти изгибных мод, а также коэффициента усиления K (L1) записывали для каждой геометрической конфигурации в текстовые файлы в формате ASCII для последующего считывания программой MathCad при построении диаграммы устойчивости.

Построение диаграммы устойчивости. При построении диаграммы устойчивости для каждого значения длины L1 и каждой моды изгибных колебаний из файла считывался параметр an (L1). Затем из характеристического уравнения находилось соответствующее ему критическое значение qn cr(L1) параметра q, соответствующее критической амплитуде 2« cr(L1) колебаний входного сечения, т. е. qn cr(L1) = = q(L1,2n cr). Так как параметр q связан с амплитудой колебаний входного сечения линейной зависимостью, q(A, 2 n cr)/q(A, 20) = 2n cr/ 20, откуда следует 2ncr(L1) = 20qncr(A)/qM(L). Параметр qn0(L1) считывался из файла. Допустимую амплитуду продольных колебаний выходного сечения вычисляли путем умножения критической амплитуды колебаний входного сечения 2n cr(L1) на коэффициент усиления K (L1), который считывался из файла.

На рис. 4 приведен расчетный график нижней границы области неустойчивости продольных колебаний, соответствующей потере

устойчивости с возникновением изгибной моды колебаний 8-го порядка.

Из рис. 4 следует, что при значениях длины Ь[, равных 0,0585 и 0,0878 м, продольные колебания волновода являются неустойчивыми даже при близких к нулю значениях амплитуды колебаний входного сечения, что соответствует точке (0, 1) на диаграмме Айнса-Стретта. В реальных условиях из-за наличия демпфирования потеря устойчивости будет происходить при более высоких значениях амплитуды [4].

20

15

, 10

-

•

'.г

0,02 0,04 0,06 0,08 0,10 ¿1, м

Рис. 4

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

^ог( ¿1) = тт %п сг( ¿1).

п

Результирующая диаграмма устойчивости приведена на рис. 5.

0,02 0,04 0,06 0,08 0,10 Li, м

Рис. 5

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

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

В Ы В О Д Ы

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

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

Л И Т Е Р А Т У Р А

1. Stepanenko, D. A. Modeling of flexible waveguides for ultrasonic vibrations transmission: Longitudinal and flex-ural vibrations of non-deformed waveguide / D. A. Stepanenko, V. T. Minchenya // Ultrasonics. - 2010. - Vol. 50. - P. 424^30.

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

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

4. Болотин, В. В. Динамическая устойчивость упругих систем / В. В. Болотин. - М.: Гос. изд-во техн.-теорет. лит., 1956. - 600 с.

5. Евразийский патент EA 005704 B1, МПК A61B 17/22, 17/32, C25F 3/16. Волновод для внутрисосудистой тромбоэктомии тромбов и тромбоэмболов и способ его изготовления / А. Г. Мрочек [и др.]. - № 200300259; заявл. 11.02.2003; опубл. 28.04.2005; приоритет 24.01.2003 BY 20030052.

6. Dynamic stiffness for piecewise non-uniform Timo-shenko column by power series - part I: Conservative axial force / A. Y. T. Leung [et al.] // International Journal for Numerical Methods in Engineering. - 2001. - Vol. 51. -

P. 505-529.

7. Корн, Г. Справочник по математике для научных работников и инженеров. Определения, теоремы, формулы / Г. Корн, Т. Корн. - М.: Наука, 1970. - 720 с.

8. Анго, А. Математика для электро- и радиоинженеров / А. Анго. - М.: Наука, 1965. - 780 с.

9. Светлицкий, В.А. Механика стержней: в 2-х ч. Ч. II: Динамика / В. А. Светлицкий. - М.: Высш. шк., 1987. - 304 с.

5

0

8

6

4

2

0