CC BY
104
УДК 534.231
МОДЕЛИРОВАНИЕ АКУСТИЧЕСКИХ ПОЛЕЙ, СОЗДАВАЕМЫХ ПРИ УЛЬТРАЗВУКОВОЙ АНГИОПЛАСТИКЕ, МЕТОДОМ ГРАНИЧНЫХ ЭЛЕМЕНТОВ
Кандидаты техн. наук, доценты СТЕПАНЕНКО Д. А., МИНЧЕНЯВ. Т.
Белорусский национальный технический университет
Эффективность разрушения материала тромба при ультразвуковой ангиопластике зависит от степени развития кавитации. Интенсивность кавитации, в свою очередь, определяется значением амплитуды акустического давления, создаваемого головкой волновода в окружающей жидкой среде (крови). Для возникновения кавитации величина интенсивности ультразвука, соответствующая амплитудному значению акустического давления, должна превышать порог кавитации. Теоретические и экспериментальные исследования показывают, что амплитудный показатель акустического давления, создаваемого волноводом в жидкости, при постоянном значении амплитуды колебательных смещений зависит от формы рабочего окончания (головки) волновода [1, 2]. В частности, теоретические расчеты с помощью метода конечных элементов (МКЭ) показывают, что при диаметре головки 1 мм, частоте ультразвука 23,5 кГц и амплитуде колебательных смещений 60 мкм сферическая головка создает максимальную амплитуду акустического давления 300 кПа, плоская головка - 550 кПа, а вогнутая сферическая головка - 1100 кПа [1]. Увеличение диаметра головки приводит к росту давления. Влияние формы головки на амплитуду акустического давления также косвенно подтверждается результатами экспериментальных исследований скорости разрушения тромбов in vitro [2]. В частности, диаграмма зависимости убыли массы тромба от интенсивности ультразвука для различных форм головки при постоянном времени ультразвукового воздействия приведена на рис. 1.
Представленные результаты позволяют расположить волноводы с различными формами головки в порядке убывания их эффективности следующим образом: волновод с плоской головкой с отверстием (WFHO), волновод с плос-
Наука итехника, № 6, 2013
кой головкой (^РН), волновод со сферической головкой (^БН), волновод со сферической головкой с отверстием (WSHO).
100 80 6040 -20 0
2 4
4,2 16,2
Интенсивность, Вт/см
46,2 2
Рис. 1. Зависимость убыли массы тромба
от интенсивности ультразвука для различных форм головки волновода: 1 - ВСФ; 2 - ВПФ; 3 - ВСФО; 4 - ВПФО
Теоретические расчеты параметров акустического поля, создаваемого волноводом в жидкости, обычно проводятся с помощью МКЭ [1, 3, 4]. Однако необходимо отметить, что решение подобных задач (внешних задач акустики) с помощью МКЭ не является оптимальным с точки зрения производительности вычислений. Дифференциальное уравнение Гельмголь-ца, пространственная дискретизация которого производится при решении внешних задач акустики с помощью МКЭ, требует задания граничных условий на поверхности излучателя и на бесконечности (условия излучения Зом-мерфельда). Необходимость учета граничных условий на бесконечности требует конечно-элементной дискретизации большой (по сравнению с размерами излучателя) области пространства вокруг излучателя с целью имитации реального бесконечно-протяженного пространства. Это приводит к значительному увеличению числа конечных элементов и продолжительности расчета. Более эффективным способом решения внешних задач акустики является применение метода граничных элементов
4
(МГЭ) [5, 6]. Данная статья посвящена изучению возможности применения МГЭ для теоретического анализа акустических полей, создаваемых при ультразвуковой ангиопластике.
Методика моделирования. Моделирование акустического поля волновода с помощью МГЭ производили с использованием программы с открытым кодом OpenBEM*, реализованной в виде библиотеки функций для MatLab [7]. По сравнению с МКЭ, МГЭ имеет следующие преимущества:
1) гранично-элементная дискретизация необходима только на поверхностях излучателей и отражателей акустических волн. Это позволяет сократить общее число дискретных элементов и вычислительную сложность задачи;
2) условия излучения Зоммерфельда автоматически удовлетворяются, так как они учтены при переходе от дифференциального уравнения Гельмгольца к интегральному уравнению Кирхгофа - Гельмгольца, пространственная дискретизация которого используется при решении акустических задач с помощью МГЭ.
Для оценки эффективности применения МГЭ для моделирования ультразвуковой ангиопластики была рассмотрена упрощенная задача расчета акустического поля, создаваемого волноводом в полубесконечной жидкой среде. В случае необходимости рассматриваемая задача может быть приближена к реальным условиям экспериментов in vivo и in vitro путем введения дополнительных границ с импеданс-ными граничными условиями, имитирующих стенки кровеносного сосуда или пробирки. Схематическое изображение рассматриваемой модели приведено на рис. 2а.
В случае продольных колебаний волновода 1 граничные условия (распределение нормальной составляющей колебательной скорости) на границе 2 его контакта с жидким полупространством 3 будут осесимметричными. В совокупности с геометрической осевой симметрией волновода это дает возможность использования двумерной осесимметричной модели. Решатель осесимметричных задач реализован в OpenBEM в виде модуля AxiBEM, основанного на осесимметричной формулировке уравнения Кирхгофа - Гельмгольца [8].
* Авторы выражают благодарность профессору Vicente Cutanda Henriquez (Syddansk Universitet, Дания), любезно предоставившему программу OpenBEM.
Практическая применимость представленной на рис. 2а модели ограничена присутствием бесконечной границы 4 между жидким полупространством 3 и полубесконечной газообразной средой 5. Вследствие своей бесконечности такая граница не может быть реализована в геометрической модели задачи. Необходимость рассмотрения границы 4 может быть устранена в случае использования метода зеркального источника, применяемого для анализа излучения акустических волн в полупространство, ограниченное импедансной плоскостью с импедансом X [9].
а б
р2С2
I-/- 1
— " р2С2 —
р1С1
— р2С2
1'
Рис. 2. Модели излучения волновода: а - в полубесконечной жидкой среде; б - с зеркальным излучателем в бесконечной жидкой среде
Основная идея метода состоит во введении дополнительного фиктивного излучателя 1', являющего зеркальным изображением первичного излучателя 1 относительно границы полупространства (рис. 2б). В случае X = 0 (акустически мягкая плоскость) излучатели должны создавать противофазные акустические давления, а в случае Z = (акустически жесткая плоскость) - синфазные давления. Поскольку в рассматриваемой задаче р,^ « р2с2, где р, и С - плотность и скорость звука для газообразной среды; р2 и с2 - плотность и скорость звука для жидкой среды, граница 4 может рассматриваться как акустически жесткая (идеально отражающая). Так как для гармонических волн
,/2л/р
где Уо - вектор амплитуды колебательной скорости; р0 - амплитуда акустического давления;
Наука итехника, № 6, 2013
f - частота колебаний; р - плотность среды; ] -мнимая единица, в случае синфазных акустических давлений ^-составляющие векторов колебательной скорости первичного и зеркального источников будут иметь противоположные знаки. Результирующее акустическое поле, создаваемое первичным и зеркальным источниками в однородной бесконечной среде без им-педансной плоскости, будет идентично акустическому полю, создаваемому первичным источником в полупространстве, ограниченном импедансной плоскостью.
Подобный метод зеркального источника используется в электростатике для расчета электрического поля точечного заряда, находящегося вблизи проводящей плоскости [10, с. 56]. Длину Ь погруженной в жидкость части волновода определяли из условия четвертьволнового резонанса на частоте f = 25 кГц. Граничные условия задавали в виде распределения нормальной составляющей \п вектора амплитуды
колебательной скорости на границе волновода с жидкостью. В случае продольных колебаний для цилиндрической части волновода О < у < , где Ь\ - длина цилиндрической части, можно принять граничное условие л>п (_)') = 0. Это условие основано на допущении, что величина радиальной деформации волновода значительно меньше величины продольных деформаций. Для участка 1л<у<Ь (головка волновода) принимали граничное условие
уп(У) =1 уо II "у (У) I 88Р(иДЯ)»
где | у0 | - модуль вектора амплитуды колебательной скорости; п (у) - ^-составляющая
внешней единичной нормали к границе волновода.
В этом условии пренебрегают изменением амплитуды колебательной скорости на протяжении рассматриваемого участка, так как Ь—ЬуКкЬ. Модуль вектора амплитуды колебательной скорости рассчитывали по формуле
1^0|=2тгД0,
где с,п = 10 мкм - амплитуда колебательных смещений.
Наука итехника, № 6, 2013
Такое же граничное условие с противоположным знаком было принято для головки зеркального излучателя.
При решении трехмерных задач с помощью OpenBEM возможно использование генерации граничных элементов с помощью внешних программ (например, применяя программное обеспечение для моделирования с помощью МКЭ, в частности конечные элементы типа SHELL в ANSYS) с последующим импортом. Для этих целей предусмотрены функции readnodes ('nlist.nod') и readelements('elist.ele'), где nlist. nod - имя файла, содержащего номера и координаты узлов, elist.ele - имя файла, содержащего номера элементов и номера принадлежащих им узлов. Эти функции также могут применяться с двумерными моделями, однако в решении двумерных задач OpenBEM использует трехузловые конечные элементы, а конечные элементы типа BEAM, применяемые для построения сетки на границах в ANSYS, являются двухузловыми.
В связи с этим мы использовали генерацию граничных элементов с помощью внутренних инструментов OpenBEM. Геометрия границы осесимметричного излучателя описывается в OpenBEM в виде последовательности прямолинейных и дуговых сегментов с помощью матрицы segments, содержащей информацию о координатах граничных точек и радиусах кривизны сегментов, а также сведения о количестве граничных элементов, используемых для разбиения каждого из сегментов. Генерация граничных элементов осуществляется с помощью функции nodegen (segments). При решении рассматриваемой задачи сегменты генерировали с помощью ANSYS путем разбиения границы излучателя на конечные элементы типа BEAM. Координаты граничных точек сегментов импортировали в OpenBEM с помощью функций readnodes и readelements.
Для подтверждения достоверности результатов, полученных с помощью МГЭ, задачу также решали с помощью МКЭ с использованием модуля Acoustics Module программы COMSOL Multiphysics. Задачу рассматривали как осесимметричную проблему о вынужденных колебаниях волновода, взаимодействующего с жидкостью (тип анализа "Acoustic-Structure Interaction —> Frequency response analy-
sis"). Геометрическая модель состояла из двух областей: области волновода и области жидкости, ограниченной акустически жесткой плоскостью (поверхность раздела жидкости и газообразной среды) и полусферической поверхностью с излучательными граничными условиями для сферических волн. Излучательные граничные условия позволяют исходящим волнам покидать область моделирования без отражений [11; 12, с. 79]. Вектор амплитуды колебательного ускорения частиц жидкости вблизи поверхности волновода был принят равным вектору амплитуды нормального ускорения поверхности (граничное условие типа Normal acceleration). Колебания волновода рассматривали с учетом давления, создаваемого жидкостью на его поверхность (нагрузка типа Fluid load). Вынуждающая нагрузка была приложена в виде заданного колебательного смещения рабочей
поверхности головки волновода (граничное условие типа Prescribed displacement). Амплитуда смещения была принята равной 10 мкм, а частота - 25 кГц. Перемещения узловой плоскости волновода, совпадающей с поверхностью раздела жидкой и газообразной сред, были ограничены по всем степеням свободы.
Результаты моделирования и их обсуждение. Результаты расчета амплитуды акустического давления на поверхности волновода со сферической головкой с помощью МГЭ приведены на рис. 3 а. Нумерация узлов граничных элементов изменена таким образом, чтобы они располагались в порядке увеличения значений дуговой координаты, измеряемой вдоль границ зеркального и первичного излучателей в направлении от точки (0, —L) к точке (0, L). Левая часть графика соответствует зеркальному излучателю, а правая - первичному.
я 1 10
С
я 1 10
я 1 • 10:
я s с
S <
100
10
100
200 300 Номер узла
400
Расчетное распределение давления
Поверхность
раздела «жидкость -твердое тело»
Колебательное смещение
Элементы FLUID29
0,066 0,064 0,062 0,060 0,058 0,056 0,054 0,052 0,050
0,8
0
0,004 0,008 0,012 0,016
24,6
11,9
-0 ,8
-13,5
-26,2
-38,9
-51,6
-64,3
-77,0
-89,7 -96,0
Рис. 3. Результаты расчета амплитуды акустического давления: а - с помощью МГЭ, поверхностное распределение; б - с помощью МКЭ, пространственное распределение согласно [1]; в - с помощью МКЭ, пространственное распределение
82
Наука итехника, № 6, 2013
б
а
в
Графики амплитуды акустического давления и его фазы (рис. 4) являются симметричными, то есть первичный и зеркальный излучатели создают синфазные акустические давления. Это согласуется с принятыми граничными условиями. На рис. 3б и 3в для сравнения представлены результаты расчета акустического поля волновода с помощью МКЭ согласно [1] и настоящему исследованию. Из сравнения приведенных графиков следует, что во всех случаях имеется основной максимум амплитуды акустического давления (основной лепесток диаграммы направленности), соответствующий точке (0, Ь), и вторичный максимум амплитуды акустического давления (боковой лепесток диаграммы направленности), соответствующий точке на задней части поверхности головки (пу (V) < 0). Значение амплитуды акустического
давления, рассчитанное с помощью МГЭ, согласуется со значением, рассчитанным с помощью МКЭ (300 кПа - при амплитуде колебательных смещений Е, = 60 м согласно [1]; 82,6 кПа - при амплитуде 10 мкм, расчет с помощью МГЭ; 96,0 кПа - при амплитуде 10 мкм, расчет с помощью COMSOL).
Моделирование с помощью МКЭ показывает, что акустические давления в основном и боковом лепестках диаграммы направленности являются противофазными (т. е. сжатие среды у передней части поверхности головки соответствует расширению у задней части). Это согласуется с результатами расчета фазы давления с помощью МГЭ, приведенными на рис. 4.
1
1
k 2 / h
100 200 300 400 Номер узла
Рис. 4. Графики амплитуды акустического давления и его фазы, рассчитанные с помощью МГЭ
Как следует из приведенного графика, между основным и вторичным максимумами амплитуды акустического давления имеется фазовое смещение, равное п. С целью упрощения положений максимумов амплитуды график 1 фазы скомбинирован с графиком 2 (рис. 4) амплитуды в произвольном масштабе.
Наука итехника, № 6, 2013
В Ы В О Д
На основе анализа простой тестовой задачи показано, что МГЭ может быть использован в качестве эффективного инструмента для моделирования акустических полей, генерируемых в ультразвуковых хирургических процедурах, в частности ультразвуковой ангиопластике. МГЭ имеет ряд преимуществ по сравнению с МКЭ и может быть использован в качестве альтернативы традиционно используемому моделированию с помощью МКЭ либо дополняющего метода. Результаты моделирования с помощью МГЭ и МКЭ хорошо согласуются друг с другом как качественно, так и количественно.
Л И Т Е Р А Т У Р А
1. Wylie, M. P. A linear finite element acoustic fluid-structure model of ultrasonic angioplasty in vivo / M. P. Wylie, G. B. McGuinness, G. P. Gavin // International Journal for Numerical Methods in Biomedical Engineering. - 2010. - Vol. 26. - P. 828-842.
2. Тун, Ц. Эффективность восстановления проходимости пораженных атеросклерозом артерий ультразвуковыми волноводами различных модификаций in vitro: дис. ... канд. мед. наук / Ц. Тун. -Минск, 2006. - 110 с.
3. Development of piezoelectric ultrasonic thrombolysis device for blood clot emulsification / T. Li [et al.] // ISRN Materials Science. - 2012. - Article ID 106484.
4. Horn-type piezoelectric ultrasonic transducer: modelling and applications / T. Li [et al.] // Advances in Piezoelectric Transducers / Ed. by F. Ebrahimi. - N.Y.: InTech, 2011. - P. 3-26.
5. Juhl, P. M. The boundary element method for sound field calculations: PhD thesis / P. M. Juhl. - Copenhagen, 1993. - 195 p.
6. Kirkup, S. The boundary element method in acoustics / S. Kirkup. - S.l.: Integrated Sound Software, 1998. - 148 p.
7. OpenBEM: Open source Matlab codes for the Boundary Element Method [Электронный ресурс]. -Режим доступа: http://www.openbem.dk/.
8. A boundary integral approach for acoustic radiation of axisymmetric bodies with arbitrary boundary conditions valid for all wave numbers / W. Wang [et al.] // Journal of the Acoustical Society of America. - 1997. - Vol. 101. - P. 1468-1478.
9. Brick, H. A half-space BEM for the simulation of sound propagation above an impedance plane / H. Brick, M. Ochmann // Proc. of Acoustics'08. - Paris, 2008 [Электронный ресурс]. - Режим доступа: http:// webistem. com / acoustics2008 / acoustics2008/cd 1/data/ articles/001156.pdf.
10. Тамм, И. Е. Основы теории электричества / И. Е. Тамм. - М.: Наука, 1989. - 504 с.
11. Bayliss, A. Boundary conditions for the numerical solution of elliptic equations in exterior regions / A. Bayliss [et al.] // SIAM Journal of Applied Mathematics. - 1982. - Vol. 42. - P. 430-451.
12. COMSOL Multiphysics. Version 3.5. Acoustics Module. User's Guide. - 2008.
Поступила 23.10.2013
