1. На главную
  2. Электронные версии
  3. Известия РАН. Серия физическая
  4. T. 86, № 2, 2022

Известия РАН. Серия физическая, 2022, T. 86, № 2, стр. 257-262

Расчет акустической ловушки для упругого сферического рассеивателя большого волнового размера

А. А. Крохмаль 1, *, Н. Е. Крохмаль 1, О. А. Сапожников 1

1 Федеральное государственное бюджетное образовательное учреждение высшего образования “Московский государственный университет имени М.В. Ломоносова”
Москва, Россия

* E-mail: aa.krokhmal@physics.msu.ru

Поступила в редакцию 01.10.2021
После доработки 11.10.2021
Принята к публикации 22.10.2021

Полный текст (PDF)

Аннотация

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

ВВЕДЕНИЕ

В большом количестве прикладных задач возникает необходимость создания акустической ловушки [1, 2] – области пространства, в которой с помощью ультразвукового пучка может быть удержана твердотельная частица. Чтобы убедиться, что данная конфигурация акустического поля является ловушкой для некоторой частицы, необходимо рассчитать действующую на эту частицу акустическую радиационную силу. Если размер частицы существенно превышает длину волны или сравним с ней, то расчет акустической радиационной силы становится непростой вычислительной задачей, требующей решения задачи рассеяния акустической волны, так как существующие упрощения для расчета радиационной силы (например, приближение Горькова [3]) подходят только для рассеивателей малых волновых размеров.

В работе [4] был предложен подход по расчету радиационной силы, действующей со стороны произвольного акустического пучка на упругую сферу в жидкости. Падающий на рассеиватель пучок представляется в виде суммы плоских волн разных направлений (метод углового спектра). Для каждой такой волны применяется классическое решение для рассеяния на упругой сфере и затем результирующее рассеянное поле представляется в виде суперпозиции рассеянных волн от всех компонент углового спектра падающего пучка. Сумма падающей и рассеянной волн позволяет рассчитать компоненты тензора радиационных напряжений на замкнутой поверхности, окружающей рассеиватель, и в итоге вычислить акустическую радиационную силу.

Хотя предложенный в [4] метод позволяет с хорошей точностью рассчитать радиационную силу для произвольного акустического поля и для любого материала и размера сферического рассеивателя, то есть применим в большом диапазоне практических задач, численная реализация этого метода может вызвать трудности ввиду громоздкости формул и необходимости одновременного учёта большого количества параметров. Чтобы облегчить вычислительную процедуру и анализ результатов, был реализован сервис с графическим интерфейсом на языке Python, позволяющий проводить расчет акустической радиационной силы, задавая посредством дружественного интерфейса параметры акустического поля и рассеивателя.

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

МЕТОД РАСЧЕТА АКУСТИЧЕСКОЙ РАДИАЦИОННОЙ СИЛЫ

Для расчета акустической радиационной силы необходимо знать акустическое поле, в которое помещен рассеиватель. Рассмотрим акустический источник, излучающий волновой пучок вдоль оси $z$. Предполагается, что источник является гармоническим по времени. Акустическое поле может быть охарактеризовано голограммой – двумерным распределением давления ${\text{на}}$ некоторой плоскости $xy$ перед источником. Это распределение может быть рассчитано на основе модели вибрации поверхности источника с использованием интеграла Рэлея или может быть измерено экспериментально. Выберем в качестве начала декартовой системы координат пересечение плоскости голограммы с осью z. Таким образом, нам известно распределение комплексной амплитуды давления падающей волны в плоскости $xy$ при $z = 0{\text{:}}$ ${{p}_{i}}\left( {x,y,z = 0} \right),$ и, следовательно, угловой спектр в этой же плоскости $S\left( {{{k}_{x}},~{{k}_{y}}} \right),$ где ${{k}_{x}},~{{k}_{y}}$ – компоненты волнового вектора $k = \{ {{k}_{x}},~{{k}_{y}},~{{k}_{z}}$ = $\sqrt {{{k}^{2}} - k_{x}^{2} - k_{y}^{2}} \} .$ Угловой спектр связан с падающим полем давления ${{p}_{i}}\left( {x,y,0} \right)$ двумерным пространственным преобразованием Фурье по координатам $x,y.$ Комплексная амплитуда акустического давления в произвольной точке пространства находится из углового спектра следующим образом [5, 6]:

(1)
$\begin{gathered} {{p}_{i}}\left( {x,y,z} \right) = \\ = \,\,\frac{1}{{4{{\pi }^{2}}}}\iint\limits_{k_{x}^{2} + k_{y}^{2} < {{k}^{2}}} {d{{k}_{x}}d{{k}_{y}}S\left( {{{k}_{x}},~{{k}_{y}}} \right){{e}^{{i{{k}_{x}}x + i{{k}_{y}}y + i\sqrt {{{k}^{2}} - k_{x}^{2} - k_{y}^{2}} ~z}}}}.{\text{\;}} \\ \end{gathered} $

В работе [2] показано, что если акустический пучок падает на упругий сферический рассеиватель, то комплексную амплитуду полного поля (суммы падающего и рассеянного полей) в случае, если рассеиватель находится в точке $\vec {r} = \left\{ {x,y,z} \right\},$ можно представить в виде следующего разложения:

(2)
$p = \frac{1}{\pi }\sum\limits_{n = 0}^\infty {{{i}^{n}}\left\{ {{{j}_{n}}\left( {kr} \right) + {{c}_{n}}h_{n}^{{\left( 1 \right)}}\left( {kr} \right)} \right\}} \sum\limits_{m = - n}^n {{{H}_{{nm}}}{{Y}_{{nm}}}\left( {\theta ,~\varphi } \right)} ,$
где jn, $h_{n}^{{\left( 1 \right)}}$ – сферические функции Бесселя и Ханкеля, ${{Y}_{{nm}}}$ – сферическая функция, ${{c}_{n}}$ – множители, зависящие от параметров рассеивателя и окружающей среды, ${{H}_{{nm}}}$ – матричные коэффициенты, характеризующие падающее поле, $r = \left| {\vec {r}} \right|.$ Коэффициенты ${{H}_{{nm}}}$ можно найти на основе углового спектра падающей волны $S\left( {{{k}_{x}},~{{k}_{y}}} \right)$ и комплексно-сопряженных сферических гармоник ${{Y}_{{nm}}}{\text{:}}$

(3)
${{H}_{{nm}}} = \int\limits_{k_{x}^{2} + k_{y}^{2} \leqslant {{k}^{2}}} {d{{k}_{x}}d{{k}_{y}}S\left( {{{k}_{x}},~{{k}_{y}}} \right)Y_{{nm}}^{*}\left( {{{\theta }_{k}},~{{\varphi }_{k}}} \right)} .$

Заметим, что ${{H}_{{nm}}}$ рассчитывается для такого углового спектра, которому соответствует положение рассеивателя в координатах $\vec {r} = \left\{ {0,0,0} \right\}.$ Если рассеиватель помещен в другую точку пространства $\vec {r} = \left\{ {x,y,z} \right\},$ то соответствующий такому положению угловой спектр находится умножением исходного углового спектра на соответствующий пропагатор:

(4)
$\tilde {S}\left( {{{k}_{x}},~{{k}_{y}}} \right) = S\left( {{{k}_{x}},~{{k}_{y}}} \right){{e}^{{i{{k}_{x}}x + i{{k}_{y}}y + i\sqrt {{{k}^{2}} - k_{x}^{2} - k_{y}^{2}} ~z}}}.$

После нахождения коэффициентов ${{H}_{{nm}}}$ они используются для того, чтобы рассчитать декартовы компоненты вектора акустической радиационной силы, действующей на рассеиватель в данной точке пространства:

(5)
${{F}_{x}} = \frac{1}{{8{{\pi }^{2}}\rho {{c}^{2}}{{k}^{2}}}}{\text{Re}}\left\{ {\sum\limits_{n = 0}^\infty {{{\psi }_{n}}} \sum\limits_{m = - n~}^n {{{A}_{{nm}}}} \left( {{{H}_{{nm}}}H_{{n + 1,~m + 1}}^{*} - {{H}_{{n, - m}}}H_{{n + 1,~ - m - 1}}^{*}} \right)} \right\},$
(6)
${{F}_{y}} = \frac{1}{{8{{\pi }^{2}}\rho {{c}^{2}}{{k}^{2}}}}{\text{Im}}\left\{ {\sum\limits_{n = 0}^\infty {{{\psi }_{n}}} \sum\limits_{m = - n~}^n {{{A}_{{nm}}}} \left( {{{H}_{{nm}}}H_{{n + 1,~m + 1}}^{*} + {{H}_{{n, - m}}}H_{{n + 1,~ - m - 1}}^{*}} \right)} \right\},$
(7)
${{F}_{z}} = - \frac{1}{{4{{\pi }^{2}}\rho {{c}^{2}}{{k}^{2}}}}{\text{Re}}\left\{ {\sum\limits_{n = 0}^\infty {{{\psi }_{n}}} \sum\limits_{m = - n~}^n {{{B}_{{nm}}}} {{H}_{{nm}}}H_{{n + 1,~m}}^{*}} \right\},$
где ${{A}_{{nm}}},~{{B}_{{nm}}}~$ и ${{\psi }_{n}}$ представляются следующим образом:

(8)
${{\psi }_{n}} = \left( {1 + 2{{c}_{n}}} \right)\left( {1 + 2c_{{n + 1}}^{*}} \right) - 1,$
(9)
${{A}_{{nm}}} = \sqrt {\frac{{\left( {n + m + 1} \right)\left( {n + m + 2} \right)}}{{\left( {2n + 1} \right)\left( {2n + 3} \right)}}} ,$
(10)
${{B}_{{nm}}} = \sqrt {\frac{{\left( {n + m + 1} \right)\left( {n - m + 1} \right)}}{{\left( {2n + 1} \right)\left( {2n + 3} \right)}}} .$

Как видно из приведенных соотношений (2), (4)–(7), везде фигурирует суммирование по $n$ от нуля до бесконечности. В действительности величина х-компонент радиационной силы довольно быстро убывает, и степень убывания зависит в том числе от множителя ${{\psi }_{n}},$ который зависит от параметра $ka$ – волнового размера рассеивателя, где $a$ – радиус рассеивателя. Чем больше величина $ka,$ тем больше членов суммирования нужно учитывать для точного расчета акустической радиационной силы: если $ka < 1,$ то требуется всего два члена ряда, и выражения сводятся к приближению Горькова для рассеивателей малого волнового размера [3].

ОЦЕНКА ТОЧНОСТИ РАСЧЕТА

Время расчета квадратично зависит от количества итерационных членов и линейно зависит от количества пространственных точек, поэтому для оптимального времени работы алгоритма нужно определить минимальное количество удерживаемых членов ряда, которое обеспечит приемлемую точность расчета. В реализованном алгоритме количество членов ряда было ограниченно величиной ${{N}_{{max}}} = ka + 5,$ где k – волновое число в среде, a – радиус рассеивателя. Было показано, что в таком случае численная ошибка расчета не превышает 10–5.

Можно представить ${{F}_{z}}$ из выражения (7) как функцию от величины верхнего индекса суммирования $N$ следующим образом:

(11)
${{F}_{z}}\left( N \right) = \sum\limits_{n = 0}^N {{{f}_{n}}} ~,$
где ${{f}_{n}}$ – n–е компоненты из выражения (7). Тогда будем считать количество членов $N$ ряда достаточным, если Fz(N) достигает асимптоты, а величина ${{\left| {{{F}_{z}}\left( {N - 1} \right) - {{F}_{z}}\left( N \right)} \right|} \mathord{\left/ {\vphantom {{\left| {{{F}_{z}}\left( {N - 1} \right) - {{F}_{z}}\left( N \right)} \right|} {{{F}_{z}}\left( N \right)}}} \right. \kern-0em} {{{F}_{z}}\left( N \right)}}$ не превышает 10–5.

Для иллюстрации были рассчитаны величины z-компоненты радиационной силы, действующей на помещенные в фокус сферические рассеиватели разных волновых размеров. Для оценки точности в качестве источника поля моделировался вогнутый фокусирующий излучатель в виде сферической чаши диаметром 10 см и радиусом кривизны 7 см, который излучал ультразвуковой пучок в воду на частоте 1 МГц. Параметры рассеивателя были заданы следующими: плотность ${{\rho }_{{sc}}}$ = 2040 кг/м3, скорости продольных и поперечных волн ${{c}_{l}}$ = 4540 и ${{c}_{t}}$ = 2130 м/с, соответственно, радиус рассеивателя от 1 до 6 мм, что соответствовало волновому размеру $ka$ от 5 до 26 (т.е. рассеиватели имели большой волновой размер). Параметры рассеивателя соответствовали материалу моногидрат оксалата кальция, который используется в клинических испытаниях по выталкиванию ультразвуком почечных камней.

На рис. 1а показана зависимость ${{F}_{z}}\left( N \right)$ для каждого размера рассеивателя. Как видно из графика, чем меньше волновой размер рассеивателя, тем быстрее величина ${{F}_{z}}$ выходит на асимптоту. Для того, чтобы оценить необходимое количество удерживаемых членов ряда N для достижения необходимой точности расчета g, были построены зависимости вида $g\left( N \right) = 1 - \left| {\frac{{{{F}_{z}}\left( {N - 1} \right)}}{{{{F}_{z}}\left( N \right)}}} \right|$ (рис. 1б). Видно, что для всех размеров рассеивателя при достижении ${{N}_{{max}}} = ka + 5$ (отмечено точкой на каждом графике) точность счета достигает значения 10–5. Аналогичные расчеты для поперечных компонент акустической радиационной силы также показали достаточность приведенного критерия для ${{N}_{{max}}}.$

Рис. 1.

Зависимость ${{F}_{z}}\left( N \right)$ (а) и $g\left( N \right)$ (б) для каждого размера рассеивателя. Рисунки демонстрируют, что для рассеивателей разных размеров количества удерживаемых членов ряда ${{N}_{{max}}} = ka + 5$ достаточно для высокой точности численных расчетов акустической радиационной силы.

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

Для того, чтобы удерживать в некоторой области пространства довольно крупный рассеиватель, предлагается создать структуру звукового поля в виде “чаши”, в центре которой амплитуда акустического давления имеет минимум, а по краям – максимум. Характерная частота ультразвукового пучка предполагается равной 2 МГц, тогда как характерные размеры рассеивателя могут достигать нескольких миллиметров в радиусе, поэтому зона акустической ловушки должна быть достаточно широкой. В такой постановке задачи возникла идея сформировать стены “чаши” акустической ловушки с помощью каустики, подобрав правильную фазу сигнала для каждого элемента излучателя.

Каустику в виде полусферы можно создать, если задать излучающую поверхность с кривизной, соответствующей эвольвенте окружности. В таком случае лучи, исходящие перпендикулярно от излучающей поверхности, образуют каустику в виде полусферы радиуса r, а параметрическое уравнение эвольвенты будет описываться уравнениями:

(12)
$x = r\left( {\cos \phi + \phi \sin \phi } \right),$
(13)
$y = r\left( {\sin \phi - \phi \cos \phi } \right),$
где $r$ – радиус окружности каустики, $\phi $ – полярный угол точки касания касательной и окружности.

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

Для задания акустической ловушки шириной 1 см был произведен расчет фазы для элементов существующего кольцевого многоэлементного излучателя. Он состоит из 12 колец с минимальным радиусом 2 см и максимальным радиусом 5 см и с радиусом кривизны 7 см. На рис. 2а показано распределение фазы излучаемого сигнала на каждом элементе преобразователя, а на рис. 2б – структура результирующего поля, рассчитанная с помощью интеграла Рэлея для излучения ультразвукового пучка частоты 2 МГц в воду. Как видно из последнего рисунка, с помощью фазирования элементов преобразователя удалось достигнуть создания каустики в виде “чаши”, в которой потенциально может удерживаться рассеивателя довольно большого волнового размера.

Рис. 2.

Распределение фазы на элементах преобразователя (а) и соответствующее поле акустического давления, формирующее каустику в виде “чаши” – акустическую ловушку для рассеивателя большого волнового размера (б).

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

РАСЧЕТ АКУСТИЧЕСКОЙ РАДИАЦИОННОЙ СИЛЫ

Указанный выше сервис создан на языке программирования Python 3, при этом графический интерфейс программы реализован с помощью фреймворка веб-разработки Flask и доступен в любом браузере. В качестве базы данных, в которой хранятся результаты расчетов и загруженные модели акустического пучка, используется PostgreSQL. Кроме того, имеется Docker образ данного сервиса, который может быть запущен на любой операционной системе, а также на облачном сервере.

В качестве источника поля была использована теоретически рассчитанная с помощью интеграла Рэлея голограмма ультразвукового пучка [5] – распределение комплексной амплитуды давления в плоскости, перпендикулярной направлению распространения пучка, представляющее собой таблицу N × N элементов в формате .mat (геометрия излучателя и плоскости голограммы показана на рис. 3). При практическом использовании в качестве источника информации об ультразвуковом пучке можно также использовать экспериментально измеренную голограмму этого же излучателя, что делает расчет акустической радиационной силы более точным.

Рис. 3.

Схематическое изображение плоскости расчета голограммы фокусирующего излучателя.

Для иллюстрации работоспособности теоретической модели в качестве рассеивателя был рассмотрен стальной шарик радиуса 1 мм с продольной скоростью звука 5740 м/с, сдвиговой скоростью звука 3092 м/с и плотностью 7800 кг/м3.

РЕЗУЛЬТАТЫ РАСЧЕТА

В разработанном сервисе после расчета акустической радиационной силы результаты выводятся в виде интерактивных графиков (рис. 4), компоненты радиационной силы нормируются на мощность звукового поля $W$ и скорость звука в среде c. Нормировка радиационной силы необходима для тех случаев, когда доподлинно неизвестна амплитуда давления на поверхности излучателя, но было проведено измерение мощности источника, например, методом акустического взвешивания и известна его мощность. Нормирование производится по следующей формуле: ${{F}_{{norm}}} = {{Fc} \mathord{\left/ {\vphantom {{Fc} W}} \right. \kern-0em} W}.$

Рис. 4.

Нормированные x- (а) и z-компоненты (б) радиационной силы на плоскости xz, действующие на стальной сферический рассеиватель радиуса 1 мм. Пунктиром обозначена зона акустической ловушки, стрелками показано направление x- и z-компоненты радиационной силы.

Как показывает распределение x-компоненты радиационной силы в плоскости xz (рис. 4), поперечная компонента радиационной силы будет притягивать рассеиватель к оси пучка в зоне акустической ловушки, продольная же z-компонента всегда положительна, действует на рассеиватель по направлению распространения пучка и может быть уравновешена силой тяжести, действующей на рассеиватель. Таким образом, предложенная конфигурация поля может оказаться эффективной акустической ловушкой для рассеивателей большого волнового размера.

ЗАКЛЮЧЕНИЕ

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

Работа выполнена при поддержке Российского фонда фундаментальных исследований (проекты №№ 20-32-90093, 20-02-00139) и Фонда развития теоретической физики “Базис”.

Список литературы

  1. Крохмаль А.А., Сапожников О.А., Цысарь С.А. и др. // Уч. зап. физ. фак. МГУ. 2020. Т. 1. С. 2010902.

  2. Крохмаль А.А., Сапожников О.А., Кудан Е.В. и др. // Изв. РАН. Сер. физ. 2021. Т. 85. № 6. С. 883; Krokhmal A.A., Sapozhnikov O.A., Kudan E.V. et al. // Bull. Russ. Acad. Sci. Phys. 2021. V. 85. No. 6. P. 681.

  3. Горьков Л.П. // Докл. АН СССР. 1961. Т. 140. № 1. С. 88.

  4. Sapozhnikov O.A., Bailey M.R. // J. Acoust. Soc. Amer. 2013. V. 133. No. 2. P. 661.

  5. Sapozhnikov O.A., Tsysar S.A., Khokhlova V.A., Krei-der W. // J. Acoust. Soc. Amer. 2015. V. 138. No. 3. P. 1515.

  6. Schafer M.E., Lewin P.A. // J. Acoust. Soc. Amer. 1989. V. 85. No. 5. P. 2202.

  7. https://github.com/nkrokhmal/rfc_calculation.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Известия РАН. Серия физическая

Архивы выпусков Информация о журнале Отправить рукопись в журнал