НЕЛИНЕЙНАЯ ^^^^^^^^^^^^^^ АКУСТИКА
534.2
МОДЕЛИРОВАНИЕ ТРЕХМЕРНЫХ НЕЛИНЕЙНЫХ ПОЛЕЙ УЛЬТРАЗВУКОВЫХ ТЕРАПЕВТИЧЕСКИХ РЕШЕТОК
© 2011 г. П. В. Юлдашев1, В. А. Хохлова12
1 Московский государственный университет им. М.В. Ломоносова 119991 Москва, Ленинские горы 2Центр промышленного и медицинского ультразвука, Лаборатория прикладной физики университета шт. Вашингтон, Сиэтл, шт. Вашингтон, 98105, США E-mail: {petr,vera}@acs366.phys.msu.ru Поступила в редакцию 20.12.10 г.
Развита новая численная модель для описания трехмерных нелинейных полей, создаваемых современными ультразвуковыми терапевтическими решетками. Модель основана на решении уравнения Вестервельта; разработанный алгоритм позволяет рассчитывать нелинейные поля периодических волн при наличии ударных фронтов, локализованных вблизи фокуса. В численном эксперименте исследована роль нелинейных эффектов при фокусировке в воде для двумерной решетки, состоящей из 256 элементов, в диапазоне интенсивностей на элементах до 10 Вт/см2. Полученные решения показали, что для характерных режимов эксплуатации современных решеток нелинейные эффекты играют важную роль и в профиле волны в фокусе происходит образование ударных фронтов.
Ключевые слова: многоэлементные ультразвуковые решетки, нелинейность, дифракция, фокусированные пучки.
УДК
I. ВВЕДЕНИЕ
В современных устройствах неинвазивной хирургии с использованием мощного фокусированного ультразвука (ШБи) все больше начинают применяться двумерные фазированные решетки, состоящие из многих элементов, расположенных на сегменте сферической поверхности [1—4]. Каждый элемент такой решетки управляется независимо, что позволяет электронным образом перемещать фокус в пространстве, создавать сложную конфигурацию поля в виде нескольких фокусов, минимизировать нагрев акустических препятствий (например, ребер) при сохранении эффективности воздействия в фокусе [3—7]. При помощи решеток можно улучшать качество фокусировки в неоднородной ткани с использованием методов обращения времени, а также отслеживать область воздействия, которая смещается за счет дыхания [1, 8].
При описании полей, создаваемых НШи излучателями, разработке протоколов облучения и предсказании соответствующих биоэффектов в ткани важным инструментом исследования является численный эксперимент [9, 10]. Уровни интенсивности в фокусе НШи-излучателей достигают нескольких десятков тысяч Вт/см2, при этом за счет эффектов акустической нелинейности в профиле ультразвуковой волны образуются ударные фронты, что принципиальным образом меняет эффективность теплового воздействия уль-
тразвука на ткань и может приводить к новым биологическим эффектам нетеплового характера [11, 12].
Нелинейные эффекты в фокусированных ультразвуковых полях с учетом образования разрывов подробно исследовались для аксиально-симметричных пучков [9, 10, 12, 13]. Трехмерные задачи решались для пучков, создаваемых излучателями диагностического ультразвука, т.е. при слабой фокусировке и слабом проявлении нелинейных эффектов [14, 15]. Исследовать нелинейные эффекты в трехмерных полях, создаваемых многоэлементными ШБи-решетками, при учете формирования разрывов до сих пор не удавалось. В этой задаче соединяются сразу несколько трудностей: сложная дифракционная структура ближнего поля и большие углы фокусировки, что требует использования точных дифракционных моделей и мелкой пространственной сетки численной схемы, а также сильные нелинейные эффекты, требующие учета большого количества гармоник либо мелкой временной сетки. Для моделирования в трехмерной геометрии с учетом образования ударных фронтов в общем случае требуются существенные затраты памяти и машинного времени, превышающие возможности современных 8МР-компьютеров, т.е. компьютеров с общей памятью.
В данной работе впервые удалось получить численное решение для нелинейного поля тера-
(а) P/Po
-50 0 50 х, мм
Рис. 1. Пространственные распределения безразмерной амплитуды давления р/ро, рассчитанные на плоскости вблизи края излучателя при г = 20 мм (а), и в плоскости % = 0 (б). Здесь р00 = р00с00и00, где и00 — амплитуда колебательной скорости на поверхности элемента.
певтической решетки в режиме развитых разрывов. Был разработан экономный по объему используемой памяти алгоритм, позволяющий проводить такие расчеты на 8МР-компьютерах. В качестве примера была рассмотрена решетка с частотой 1.2 МГц, состоящая из 256 элементов радиусом 3.3 мм, расположенных на сферической чашке радиусом 68 мм, с отверстием для диагностического датчика радиусом 9 мм, фокусным расстоянием 120 мм и интенсивностью акустического поля вблизи излучающего элемента до 10 Вт/см2 [3, 16].
II. ТЕОРЕТИЧЕСКИЙ МЕТОД
1. Уравнение Вестервельта
Моделирование поля решетки проводилось на основе уравнения Вестервельта [17], которое в со-
провождающей системе координат можно записать в виде:
= С0 Ар + + дЦ (1)
дтдг 2 2р 0с0 дт 2с0 дт
Здесь р — акустическое давление, г — выделенное направление вдоль оси пучка,т = I - г/с0, I —
время, Ар = д2р/дг2 + д2р/ду2 + д2р/дх2, х и у — поперечные к г пространственные координаты; р0, с0, в и 8 — плотность, скорость звука, коэффициент нелинейности и коэффициент поглощения в среде соответственно. Расчеты проводились в воде, соответствующие физические параметры в уравнении (1) были равны: р0 = 1000 кг/м3,
с0 = 1500 м/с, р = 3.5, 5 = 4.33 х 10-6 м2/с. Начало координат соответствовало центру сферического сегмента, на котором были расположены одиночные излучатели, так что точка х = 0, у = 0, г = Р соответствовала геометрическому фокусу решетки. Уравнение (1), описывающее распространение нелинейных волн в вязкой среде в положительном направлении оси г, широко использовалось ранее для моделирования слабонелинейных и слабофокусированных полей, создаваемых диагностическими датчиками [18, 19].
2. Граничные условия для численного алгоритма
Для решения уравнения Вестервельта (1), записанного в эволюционной форме по координате г, требуется задание граничных условий на некоторой начальной плоскости (х,у, г = г0). Поскольку элементы решетки расположены на поверхности сферической чашки, то вначале рассчитывалось поле на плоскости г0 = 2 см от центра решетки вблизи ее края, расположенного на расстоянии г = 1.85 см от центра, с использованием интеграла Рэлея:
р(г) = -Р0с0 ± ^Щ'Щ^Ьз', (2)
2я л \г - г\
5 1 1
где к = ю/с0 — волновое число, ю = 2п/, / — частота ультразвука, и(г') — комплексная амплитуда колебательной скорости на поверхности излучателя 5.
Сложная дифракционная структура поля давления, создаваемого 256 элементами на указанной плоскости (х,у,г = г0), показана на рис. 1а. Чтобы учесть влияние нелинейных эффектов в ближнем поле решетки, далее граничное условие переносилось линейным образом на плоскость (х, у, г = 0), касательную к центру сферической поверхности решетки, рассчитывая интеграл Рэлея описанным ниже методом углового спектра. Полученное поле "виртуальных" источников на плоскости (х, у, г = 0) показано на рис. 1б. Числен-
ные решения рассчитывались, используя оба граничных условия, наблюдаемые различия в результатах моделирования обсуждаются ниже.
3. Численный алгоритм
Численное решение уравнения (1) строилось последовательно, переходя от плоскости (х, у, г = г1) к плоскости (х, у, г = г1 + Аг) с шагом Дг, следуя методу расщепления по физическим факторам [9, 10, 13, 14, 18, 19]. Отметим здесь лишь основные особенности применения этого метода для данной задачи. Уравнение (1) разбивалось на более простые уравнения для дифракции:
Хт
д р _е0
= - Ар,
нелинейности:
и поглощения:
дтдг 2
др _ в др2 дг 2р0с03 дт '
др _ 5 д2р дг 2с0 дт2
(3)
(4)
(5)
(6)
р(т,х,у,г) = ^ р„(х,у,г)ехр(-/лют).
(7)
П=-Х т
Дифракционный оператор (3) рассчитывался для амплитуд каждой из гармоник методом углового спектра [18—21]. Согласно этому методу, амплитуда давления р„(х,у,г) п-ой гармоники в плоскости (х,у) с помощью быстрого преобразования Фурье (БПФ) разлагается в двумерный спектр р„(кх, ку, г) по пространственным частотам (кх, ку). Компоненты углового спектра на следующем шаге рп(кх, ку, г + Аг) получаются умножением на соответствующий фазовый множитель
рп(кх, ку, г + Аг) = рп(кх, ку, г) х х ехр[/Аг^кП - к\ - к] - /Агкп],
(8)
При известном распределении давления в плоскости (х, у, г = г1), комбинируя определенным образом решения этих уравнений последовательно друг за другом, можно получить решение уравнения (1) в плоскости (х, у, г = г1 + Аг) со вторым порядком точности, т.е. ошибкой аппроксимации на каждом шаге 0(Аг ), а на всем интервале расчета по координате г — 0(Дг ) [13, 18, 19]. В данной работе метод расщепления применялся в том варианте, когда каждый шаг по координате г начинается и завершается оператором дифракции, рассчитываемым на половинном шаге сетки. Если обозначить действие оператора дифракции на шаге Аг как Г 0 Дг, а совместное действие операторов нелинейности и поглощения как Г х+Дг, то схема применения метода расщепления выглядит следующим образом:
р(т, х, у, г + Аг) =
= ГД Аг/2Г N+Л,АгГДАг/2р(т,х,у,г) ■
По мере изменения расстояния г и увеличения амплитуды высокочастотных компонент спектра волны величина шага Аг изменялась с целью обеспечения необходимой точности решения.
Моделирование проводилось с использованием преимуществ как спектрального, так и временного представлений акустического поля. Для перехода к спектру, решение уравнения (1) представлялось в виде разложения решения в конечный ряд Фурье [10]:
где кп = лю/с0 — волновое число п-ой гармоники.
Для подавления отраженных волн, появление которых связано с использованием конечной пространственной области для поля рп(х, у, г) и периодическими граничными условиями по координатам х и у, использовался метод обнуления компонент спектра на пространственных частотах вне круга радиуса ктах, т.е. в области
кI + к2у > к2ах [18]. Радиус круга для каждой из гармоник определялся соотношением ктях =
= кп/д/грГор Ак2 /п2 + 1, где А к = 412п/Ь — шаг по пространственной частоте в случае квадратной области с размером Ь, гргор — пройденное расстояние, различающееся для каждой гар
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.