CC BY
0
Вычислительные технологии, 2022, том 27, № 2, с. 19-36. © ФИЦ ИВТ, 2022 Computational Technologies, 2022, vol. 27, no. 2, pp. 19-36. © FRC ICT, 2022
ISSN 1560-7534 elSSN 2313-691X
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
D01:10.25743/ICT.2022.27.2.003
Моделирование акустического рассеяния
от множества звукопроницаемых сфер в трехмерном
пространстве
Э.Ш. Насибуллаева
Институт механики им. P.P. Мавлютова УФИЦ РАН, 450054, Уфа, Россия Контактный автор: Насибуллаева Эльвира Шамилевна, e-mail: elviraSanrb. ru
Поступила 26 октября 2021 г., доработана 4 февраля 2022 г., принята в печать 11 февраля 2022 г.
Для решения задачи акустического рассеяния от множества звукопроницаемых сфер различных радиусов, произвольно расположенных в трехмерном пространстве, при внешнем воздействии проведена оптимизация численной методики с целью минимизации используемого объема оперативной памяти и машинного времени. Верификация компьютерной модели показала хорошее соответствие численным данным других исследователей. Проведенная серия вычислительных экспериментов позволяет получить количественную оценку числа усечения рядов при разложении, которая может быть использована при исследовании ресурсоемких задач. С помощью представленной оптимизированной методики построены диаграммы распределения давления вне и внутри системы сфер, которые позволяют определить зоны повышения и понижения давления.
Ключевые слова: звукопроницаемая сфера, акустическое рассеяние, метод разложения по мультиполям, вычислительный эксперимент, монопольный источник излучения, диаграмма распределения давления.
Цитирование: Насибуллаева Э.Ш. Моделирование акустического рассеяния от множества звукопроницаемых сфер в трехмерном пространстве. Вычислительные технологии. 2022; 27(2):19-36. D01:10.25743/ICT.2022.27.2.003.
Введение
Изучение рассеяния звуковых волн на неоднородностях малых размеров, отличающихся плотностью и/или сжимаемостью от основной среды, — одна из важнейших задач акустики. Прежде всего это связано с тем, что на данном явлении основываются многие практические применения акустических волн, которые позволяют, например:
• в гидролокации [1, 2] определять местоположение косяков рыб, отдельных крупных рыб или других подводных объектов и невидимых подводных препятствий;
• в приборах неразрушающего контроля [3] измерять как геометрические параметры при одностороннем доступе к изделию, так и физико-механические свойства металлов и металлоизделий без их разрушения;
• в медицинских сканерах [4, 5] визуализировать внутренние органы и процессы, протекающие в тканях, для диагностических целей;
• при зондировании атмосферы и океана [1, 2, 6] получать информацию о поверхности Земли и объектах на ней, об атмосфере, океане и верхнем слое земной коры;
• при создании позиционируемого ЗБ-звука [7] придавать звуковой модели реализм и усиливать ощущения при восприятии звука слушателем. Таким образом, задача акустического рассеяния от препятствий малых размеров при внешних воздействиях является актуальной.
Случай одиночной сферы (сферическая твердая частица, пузырек или капля) достаточно хорошо исследован (см., например, работы [8-12]), Можно выделить несколько статей [13-15], посвященных рассеянию звуковых волн от пары сфер, где рассеянные поля от сфер хотя и взаимодействуют друг с другом, но данное взаимодействие является достаточно простым для того, чтобы его можно было подробно изучить, В последние два десятилетия основное внимание сосредоточено на рассеянии от множества сфер при внешних воздействиях [16-20], Так, в [16] представлена основанная на быстром методе мультиполей численная техника, которая разработана для решения задачи акустического рассеяния от множества звуконепроницаемых (волна не проходит через поверхность) сфер с произвольным акустическим импедансом. Вычислительные эксперименты, проведенные с помощью данной техники и метода граничных элементов (в пакете программ COMET), показали, что предложенная в работе [16] методика расчета позволяет достигать высокой точности получаемых результатов при минимальных затратах машинного времени,
В работе [17] представлены теоретические исследования рассеяния от множества резонансных воздушных пузырьков под действием падающей плоской волны (т, е, рассмотрен случай звукопроницаемых сфер, когда волна проходит через поверхность раздела) , а также численные расчеты обратного сечения рассеяния для простых кластеров из двух или трех пузырьков в различных конфигурациях и для разных углов падения плоской волны,
В работе [18] обобщена численная техника разложения потенциала поля по муль-типолям на случай акустического рассеяния от множества звукопроницаемых сфер с центрами, расположенными на одной оси, при прохождении сферической волны от монопольного источника излучения, а в работе [19] — при падении плоской волны под произвольным углом к оси, соединяющей центры сфер. Диаграммы, построенные в данных работах, позволяют наглядно продемонстрировать полную картину распределения давления вне и внутри системы сфер,
В работе [20] предложен альтернативный подход к теории многократного рассеяния, в которой рассеянное отдельной частицей поле имеет монопольный, дипольный и ротационный тип рассеяния звука. Однако данный подход из-за осреднения звукового поля не позволяет определить давление в конкретной точке пространства, следовательно, такие эффекты, как повышение или понижение давления в некоторых зонах, не могут быть обнаружены,
В настоящей работе рассматривается задача акустического рассеяния от множества звукопроницаемых сфер различных радиусов, произвольно расположенных в трехмерном пространстве, при внешнем воздействии с целью получения распределения давления вне и внутри сфер в широком диапазоне параметров системы. Для численного решения данной задачи методика расчета, представленная в работе [16], адаптирована на случай звукопроницаемых сфер. При большом числе сфер в конфигурации системы так же, как и при увеличении их взаимодействия (вызванного, например, более сильным воздействием внешнего поля или более плотным расположением сфер в конфигурации), решение задачи требует значительных вычислительных ресурсов и машинного времени. Поэтому в работе проведена оптимизация алгоритма построения матрицы ко-
эффициенгов при повторном разложении с цепью минимизации используемого объема оперативной памяти. Выполнен анализ затрат вычислительных ресурсов в зависимости от числа усечения рядов при разложении.
1. Постановка задачи и основные уравнения
Рассматриваются N сфер различных радиусов а 1; а2, ..., а^ с центрами, которые имеют декартовы координаты г^ = (х^ , ), V = 1,...,М, и расположены в бесконечном трехмерном пространстве, заполненном однородной средой, которая характеризуется плотностью р0 и скоростью звука с0. На рис, 1 представлены схема задачи и основные обозначения. Моделирование осуществляется при следующих предположениях:
• сферы звукопроницаемы, т.о. через их границу проходит волна, которая распро-
страняется в другой среде, имеющеи плотность pv и скорость звука cv, v
1,
N:
• система сфер подвергается внешнему воздействию сферической волны от монопольного источника излучения, расположенного в произвольной точке МБ(х3) трехмерного пространства;
• центры сфер неподвижны и радиальное движение сферической поверхности отсутствует.
Задача рассеяния звука от системы звукопроницаемых сфер сводится к решению уравнений Гельмгольца для комплексных потенциалов вне сфер ф(г) и внутри у-й сфе-
ры [10, 18]
V2^ + = 0, V24S + kl ф^ = 0
с граничными условиями вида
(^(r) - 4S) = 0, (
1 дф Ро дг
0,
1
Pv дг
)
0, v = 1,
Ж
Здесь к0 и к.и — волновые числа для внешней и внутренних сред.
MS (Xs,Vs,Zs)
z O
MC
Рис. 1. Схема и основные обозначения в различных системах отсчета: MS — монопольный источник излучения; МС — расчетная точка
Fig. 1. Scheme and basic notations in various reference systems: MS is monopole radiation source; MC is calculation point
s
■u
V
В соответствии с методом повторного разложения по мультпполям [16, 17] потенциал внешнего поля представляется в форме
ф(г) = фы (г) + фюл(г),
где фуп (г) — потенции падающего поля; ^са^г) — потенциал поля рассеяния, удовлетворяющий условию излучения Зоммерфельда [21]:
Иш г
(
дф.
scat
дг
- гкоф.
Q^scat^
0.
Далее проводятся разложение всех потенциалов и повторное разложение потенциала поля рассеяния по мультпполям, их вид подробно описан, например, в работе [18], где рассмотрен случай системы звукопроницаемых сфер, центры которых расположены на одной оси, В настоящей работе задача не является осесимметрнчной, поэтому при повторном разложении мультиполей ¿^(гг,) необходимо использовать общую формулу [16]
S?(rv) = ЕЕ (s\RYin«w)Rt(rw), V,w = 1,...,N, v = w.
l=Q s=-l
Здесь (в(г^) — коэффициенты перехода при повторном разложении; Д™(гад) — регулярные фундаментальные решения уравнения Гельмгольца в сферических координатах, связанные с и>-й сферой [22],
Для упрощения записи дальнейших соотношений введем следующие обозначения:
и = (ь - 1Ж + I)2 + (1 + I)2 - (I - з), й=(1 + I)2 - (I - з), г = (<ш - 1)(щг + I)2 + (п + I)2 - (п - т), 1=(п + I)2 - (п - т), Ъ = уМ - ь(ь + 1)/2 - (М - 1»),
где по умолчанию индексы принимают значения
(1)
I = 0,1 ,...,щг, з = -1,...,1, п = 0,1 ,...,щг, т = -п,...,п, =1,...,Ж
(2)
Задача в общем случае сводится к решению системы линейных уравнений отноеи-
тельно неизвестных коэффициентов А( ) , записанной в матричном виде
ЬЛ = Б
с матрицей и векторами, скомпонованными следующим образом: Л = {Л} = (Л(1),... , Л(м ))т = (|Д((1)5}Т, ■
(3)
|ДГ)s }T)T
D = {Dt} = (D(1), ••• , D(N))T = ({^1)m}T, ••• , [DlN )m}T)
L = {Lut} =
/Li11) L(12) L(21) L(22)
где элементы определяются как
l(1n )\
L(2N) L(NN)
( /r(11)
¿(21)
ut,
¿(1-2)
utt ¿(2-2) ut tt
4N2)}
L
L
(2 N) ut tt
( N N) ut tt
T
д = ^ = ^(w)m = — jn (ко aw )jn (kw aw ) Kw jn (ко 0>w )jn (kw Q>w ) ^,„1 (r' )
n hn(koaw )jn(k h'n(koaw)jn(kwaw) m'n w '
Lut = L^"rv) = <( hi(koaw)j[(kwaw) - kwh[(koaw)ji(kwaw)
ji(koaw)j'i(kwaw) - kwj[(koaw)ji(kwaw^ , (4) -(ь\Щы (r'vw) дая v = w,
ut
8й1 дая v = w.
Здесь знак "Т" операцию транспонирования; jn(z) и Кп(г) — сферические
функции Бесселя и Хапкеля первого типа соответственно [22]; к,ш = (к0рт)/(ктр0) — введенный параметр; С^11(г1ш) — коэффициенты разложения падающего поля около г = гШ; — символ Кронекера, Отметим, что при численной реализации система уравнений (3) должна быть конечной, поэтому при разложении в ряды проводится их усечение по / и п при некотором фиксированном значении щг. Вопрос о выборе числа щг является важным и будет подробно рассмотрен в подразд, 2,2,
2. Программная реализация техники разложения по мультиполям
При численной реализации техники разложения по мультиполям разработан программный код на языке программирования Fortran 90 (GCC) в среде MSYS2 (MinGW-w64) с подключением библиотеки LAPACK [23] для решения системы (3) с помощью функции zcgesv (для комплексных матриц общего вида). Для вычисления специальных сферических функций hi(z), ji(z) и присоединенных функций Лежандра Pf (z) и их производных были адаптированы программные коды [24] для специальных цилиндрических функций и полиномов Лежандра, написанные на языке Fortran 77 ,
Поскольку матрица L и векторы А и D в матричном уравнении (3) имеют размеры N(ntr + 1)2 х N(ntr + 1)2 и N(ntr + 1)2 соответственно, при большом числе сфер N и числе усечения рядов щг при разложении данная матрица и векторы занимают значительный объем памяти, С целью минимизации используемого объема оперативной памяти проведена оптимизация алгоритма [16] для определения коэффициентов перехода при повторном разложении (S\R)f™(r'vw) (см, подразд, 2,1), а для минимизации машинного времени выполнен анализ подходов к определению числа ntr (см, подразд, 2,2),
2.1. Оптимизация алгоритма определения коэффициентов перехода при повторном разложении (S\R)^(r'vw)
Благодаря соотношению симметрии
(s\R)in«w) = (-1)l+n(S\R)t:^(rlwv) (5)
коэффициенты повторного разложения достаточно вычислить для w = v+1,..., N. Тогда вместо двухмерной матрицы (S\R) размер а N (щг + 1)2 х N (щг + 1)2-, представленной в работе [16], с учетом введенных обозначений (1) и (2) можно ввести трехмерную матрицу (S\R) = {(S\R)tyt}, в которой элементы определяются следующим образом:
(5\R),ui = (S\R)in(r'vw), W = V + 1,..., N.
Матрица имеет размер N(N — 1)/2 х (щг + 1)2 х (щг + 1)2-, что уменьшает объем выде-
50(N — 1)
ляемой памяти под данную матрицу на -—- процентов,
С учетом равенства
(s\R)n<w) = (-l)l+n(S\R)-r-s «w)
получим следующий алгоритм заполнения матрицы (S|R) с помощью рекуррентных соотношений для w = v + l,..., N:
• рекуррентный процесс начинается с определения коэффициентов (S\R)S°(r'vw) и (S\R)on(r'vW) 110 формулам
(5 = V^(—l)1 s-s (r'vw)
(s\R)
vlu
V^SS(r'vw) (при I =0);
• вычисляются вспомогательные коэффициенты для п,т = —ntr,... ,ntr
bm
(n — m — 1)(п — m) (n — m — 1)(п — m) ^0, \m\ > n;
0 < m < n,
—n<m< 0,
• для m = 0,...,ntr — 1, I = 1,..., ntr коэффициенты вида (S\R)Syn\ и (S \Д) вычисляются по формулам
щ = I2 — (I — s), u2 = (l + 2)2 — (I — s), t0 = m2 + 1, ti = (m +1)2 + 1, t2 = (m + 2)2,
sm |s|n
(s\R)
bj S(S\ ^)V-Vq(VQ —1) — tf+1(S\ V('V2—2)(Vi — 1)
V'Vv2 = 7- mm— 1
"m+1
_ tf(S\R)v(uQ + 2)tQ — fy + 1 1 (S\R)vU2tQ
(S\R)
vuti
h-m—1 um+1
а для m = 1,... ,ntr — 1, I = 1,..., ntr7 s = —I + 1,... ,1 — 1 — по формулам u3 = (l + 1)2 — (I + s), to = m2 + 1, ¿3 = (m + 1)2, (s\R)vt3u3 = (—l)l+m(s\R)vutQ, \R)vtQ'U3 = (—l)l+m(s\R)vut3; • вычисляются вспомогательные коэффициенты для n,m = —ntr,... ,ntr
/
(n + l + \m\)(n + l — \m\)
, n > \m\,
(2n + l)(2n + 3)
0, \m\ > n;
• вычисляются остальные коэффициенты матрицы по формулам
щ = I2 — (I — l — s), u2 = (l + 2)2 — (I + l — s),
U3
(I + l)2 — (I + s),
t0 = n2 — (n — l — m), t2 = (n + 2)2 — (n +l — m), t3 = (n + 2)2 — (n + l + m),
(S\R)vut2 =
(S\R)v t3u3
ari-1 (S\R)vuto aS (S\R)vu2t + aS-1(S\R)vuot
(—l)n+w (S \R)vut2 nP ml = (n + l)
для п = 0,... ,щг — 1 I = п + 1, . . . ,ЩГ, 8 = —I + 1,... ,1 — 1. В результате с учетом соотношения (5) для определения коэффициентов Ь при V < т получим выражения
(wv)
ut
(4)
в
Т (ши) _ 31(ко аш) ][ ( к
) - КшЗ[ (ко аш )31 (кш^ш) ,с
ш Ц(коаш)31 (кшОьш) - КшЦ(коаш)]г(к ш ш )
-(иш) = 1+п31 (коау) Л( кцац) - КуЦ(коау)]1(куау) Ц(коау)Л(куЧу) - КуЦ(коау(куЧу)
Таким образом, все коэффициенты основной матрицы Ь системы (3) определены, 2.2. Выбор числа членов усечения рядов при разложении
Вопрос правильного выбора числа щг членов при усечении рядов важен, поскольку при малом числе щг точность расчетов будет низкой, а при большом — возрастет не только точность, но и время расчета.
Анализ литературы по данной проблеме показал, что можно выделить два подхода к определению числа щг.
Подход 1. Усечение рядов основывается на сравнении двух последовательных значений суммы искомого ряда Е (при п = щг и п = щг + 1) — как только их относительная погрешность 6 = | (Еп. - Еп.+\) /Еп. | • 100 % становится меныне 8^Х, дальнейший расчет суммы ряда прекращается и принимается значение щг = п^ [9],
Подход 2. Выполняется усечение всех рядов в каждом разложении при фиксированном числе щг, определяемом с помощью одной из эвристических формул:
щг
к^ + (2/3 (коау)1/з + 1
(6)
= ( [е ко г'уш] для малых ког'уш, , ,
Пгг \ [еког'уш/2] для больших ког'уш.
Здесь [г] — целая часть числа г; е — точность, заданная при разложении в ряд Фу-
формула (7) приведена в работе [16].
Заметим, что вывод формулы (7) основывается на том, что асимптотическое разложение сферических функций Ханкеля кп (г) при больших п и фиксированных г может привести к росту данных функций, начиная с ^ ~ /2
функции Ханкеля Цп(ког'уш) входят в матрицу коэффициентов перехода при повторном разложении (8|И.), то начиная с щг, определенного формулой (7) при наименьшем значении г'уш (= 1,... V = и>), возможен экспоненциальный рост решения вместо его стабилизации. Следовательно, данная формула позволяет получить только верхнюю оценку числа щг, а не значение, при котором происходит стабилизация решения.
В работе [27] представлено сравнение эвристических формул (6) и (7) на примере двух звуконепроницаемых сфер, однако данное сравнение ограничено только одной точкой. В работе [28] проведено сравнение двух подходов для случая трех звукопроницаемых сфер при сильном их взаимодействии. Показано, что в рассмотренном диапазоне параметров (для коа > 5) оптимально применение комбинированного подхода: с помощью подхода 2 по формулам (6) и (7) вычислить минимальное значение щг', затем, начиная с этой величины, применять подход 1 до достижения необходимой точности 8/гХ. Для анализа данных подходов к определению числа щг для конфигураций, состоящих из множества звукопроницаемых сфер, при значениях волнового радиуса 1 < коа < 5 необходимо провести серию вычислительных экспериментов (см. подразд. 3.3).
3. Вычислительный эксперимент
Проведен ряд вычислительных экспериментов дня верификации численной методики расчета, профилирования программного кода и анализа подходов к определению числа усечения рядов при их разложении. На конкретных примерах с помощью оптимизированной методики расчета продемонстрированы возможности получения результатов дня различных значений параметров системы.
3.1. Верификация оптимизированной численной методики
Большинство работ по изучению акустического рассеяния от множества сфер численно исследуют системы, содержащие от одной до трех сфер. Рассмотрим конфигурацию, представленную в работе |16|, поскольку она — одна из немногих, где проводились расчеты для систем с N > 3. В данной работе рассмотрен рассеивающий слой, который состоял из N = 11 х 11 = 121 звуконепроницаемых жестких сфер одного радиуса а = 1, равномерно распределенных в плоскости Оуг с наименьшим расстоянием между центрами соседних сфер 61 (конфигурация представлена на рис. 2). Монопольный источник излучения МБ находится на оси Ох на расстоянии й = 10а от плоскости, содержащей центры сфер, а точка расчета МС расположена симметрично относительно данной плоскости.
В настоящей работе исследуются звукопроницаемые сферы, дня которых имеет место согласование со случаем малых жестких сфер |10|, когда отношения упру гостей вещества сферы [х-и = Р'ос-) и внешней среды (%0 = РосО), и плотностей внутренней и внешней сред стремятся к бесконечности, т.е. /\0 ^ ж и р-и/р0 ^ то, Тогда для возможности верификации с результатами работы |16| выбирается внешняя среда с физическими свойствами воздуха (р0 = 1.205 кг/м3, с0 = 343.1 м/с), а внутренняя — со свойствами воды (ри = 998 кг/м3, си = 1484 м/с), поскольку в этом случае для всех сфер имеем
Рис. 2. Рассеивающий слой из 121 равномерно распределенных в плоскости Oyz звукопроницаемых сфер одного радиуса а: 51 — наименьшее расстояние между центрами сфер; MS = М(—10а, 0, 0) — монопольный источник излучения; МС = М(10а, 0, 0) — расчетная точка
Fig. 2. A scattering layer of 121 sound-permeable spheres of the same radius a which uniformly distributed in the plane Oyz: 51 is the least distance between the centers of the spheres; MS = M(—10a, 0, 0) is monopole radiation source; МС = M(10a, 0, 0) is calculation point
Х'о/хо ^ 1й Ру /Ро ^ 1. Наименьшее расстояние между центрами соседних сфер 51 = 5а. Выберем систему координат так, чтобы плоскость Оуг содержала центры сфер, ось Ох проходила через центр рассеивающего слоя, монопольный источник излучения ложа;: на оси Ох на расстоянии й = 10а от начала координат (в точке МБ = М(-10а, 0, 0)), тогда точка расчета МС = М(10а, 0, 0) (рис, 2),
Дня монопольного источника излучения падающее поло определяется но формуле
е г&о|г-г3|
ф^уесг) = -Уо~,-Г,
4^|г - г,|
а коэффициенты разложения падающего ноля имеют следующий вид |18|: С™п(г'у) = -Уо1коК( Ыг3 - гу ^^(вг^ ,<Рг.-к),
где Уо — производительность монопольного источника. Здесь и далее временной множитель е-ш* (ш = 2тг/ — угловая частота; f — частота внешнего поля) опущен. На рис, 3, а представлены результаты расчета функции
STF = 20 lg
\мс
Фш \
щ с
(8)
при щг = 5 для волновых радиусов коа = 1, 3 и 5, полученные в работе [16] (большие голубые закрашенные круги) и в настоящей работе (красная линия с закрашенными кругами). Относительные погрешности, отнесенные к расчетным данным работы |16|, равны 0,19% для коа = 1, 0.92% для коа = 3 и 4.64% для коа = 5. Увеличение погрешности с ростом волнового радиуса связано с тем, что чем больше значение коа, тем больше число щг, при котором достигается устойчивое решение (см. подразд, 3.4), следовательно, тем больше отклонение конечного результата от случая малой жесткой сферы.
12345 12345
ко a ко a
Рис. 3. Зависимость функции STF (8) в точке МС(10а, 0, 0) от волнового числа коа для 51 = Ьа{а) и 51 = 3а{б): 1 — данные для случая жестких сфер при щг = 5, полученные в работе [16]; 2 — данные настоящей работы для случая капель воды в воздухе щг = 5; 3 — капли воды в воздухе (Сф I); 4 воздушные пузырьки в воде (Сф II); 5 капли дихлорэтана в воде (Сф III)
Fig. 3. The dependence of the STF function (8) on the wave number k0a at the point MC(10a, 0, 0) for 51 = 5a {a) and 51 = 3a (б): 1 is data for the case of rigid spheres at щг = 5 obtained in [16]; 2 is data of present work for the case of water drops in air at щг = 5; 3 are water drops in air (Sph I); 4 are air bubbles in water (Sph II); 5 are dichloroethane drops in water (Sph III)
Таким образом, результаты настоящей работы, полученные для капель воды в воздухе, в пределах небольшой погрешности хорошо согласуются с результатами, представленными в работе [16] для малых жестких сфер,
3.2. Профилирование программного кода
В табл. 1 в зависимости от числа усечения рядов ntr представлены значения таких параметров, как:
• объем памяти VL, занимаемый самой большой расчетной матрицей — основной матрицей системы (3);
• используемый объем Vram оперативной памяти;
• общее время tcpu выполнения программного кода;
• доля времени 5eq = teq/tcpu-l00 %, приходящаяся на решение системы уравнений (3);
• относительные погрешности и ^з для функции 1р/рт |> вычисленные в точках М1(10а, 0, M2(0,a, 0) и М3(0, 0.5а, 0) соответственно.
Рассматривалась конфигурация, представленная на рис, 2, в случае капель дихлорэтана (pv = 1252.6 кг/м3, cv = 1034 м/с) в воде при к0а = 5 и 51 = 5а. Расчеты проводились в ОС Windows 10 на процессоре Intel Core i7-4702MQ (4 ядра и 8 потоков по 2,2 ГГц, кэш 6 Мб) с 32 Гб ОЗУ, Затраченное на выполнение программы общее время tcpv определялось консольной командой time; профилирование кода (погрешность 5eq) выполнялось с помощью команды cpu_time; используемый объем ОЗУ Vram определялся с помощью диспетчера задач Windows,
Анализ данных, приведенных в табл. 1, показывает, что основное время при выполнении программы занимает решение системы уравнений (3), Увеличение числа ntr даже на единицу приводит к значительному увеличению как объема оперативной памяти (в 3,11 раза для ntr = 1 и в 1.34 раза для ntr = 12), так и машинного времени (в 4 раза для ntr = 1 и в 1.64 раза для ntr = 12).
Таблица 1. Использование вычислительных ресурсов и относительные погрешности в зависимости от параметра ntr
Table 1. Analysis of the usage of computing resources and the relative errors depending on the parameter ntr
ntr VL, Мб Vrami Мб tcpui С 5eg, % Si, % ¿2, % ¿3, %
1 3.57 12.7 0.2 76.9 - - -
2 18.10 55.3 0.8 83.7 29.628 10.318 467.333
3 57.19 172.5 3.7 92.0 89.038 41.470 14.163
4 139.63 419.6 13.5 94.1 49.562 4.397 29.127
5 289.53 868.8 45.2 95.2 16.852 67.555 25.143
6 536.39 1608.5 113.0 96.3 28.429 326.222 14.079
7 915.06 2743.2 249.3 96.9 29.711 51.019 0.021
8 1465.75 4393.1 507.3 97.1 27.795 82.504 5.489
9 2234.04 6695.0 945.6 97.5 4.823 4.029 2.587
10 3270.86 9801.3 1750.1 97.9 9.358 9.353 0.420
11 4632.50 13880.9 2835.8 98.2 3.299 1.405 0.226
12 6380.64 19118.4 4642.6 98.4 0.498 0.188 0.009
13 8582.29 25727.8 7636.1 98.7 0.002 0.042 0.005
3.3. Применение подхода 1 к определению числа усечения рядов п4г при разложении
В табл. 1 также представлены значения относительных погрешностей Sj (] = 1, 2, 3) для последовательных значений щг на трех примерах расчета модуля нормированного давления |р/р-ш| в следующих точках пространства: точка вне сфер МС = М\{\0а, 0, 0), граничная точка М2(0, а, 0) и точка внутри сферы М3(0, 0.5а, 0), Поскольку связь между акустическим давлением р и потенциалом ф определяется по формуле
дф . Р = Р^ = шР
значения р и р[п выражаются через соответствующие потенциалы с точностью до постоянной величины %шр.
С помощью подхода 1 к определению числа щг и, например, выполнения условия < = 1% 0 = 1, 2, 3) го табл. 1 получаем, что щг = 11 для точек и М2, щг = 6 для точки М3, однако стабилизация решения для М3 устанавливается при щг = 9, На рис, 4 приведены зависимости функции |р/р;п | от числа щг в данных трех точках.
3.4. Параметрический анализ
На рис, 3 представлены значения функции STF (8) в точке MC = М1(10а, 0, 0) в зависимости от волнового радиуса k0a = 1, 2, 3, 4, 5 при 61 = 5а (см. рис. 3, а) и 61 = 3а (см. рис. 3, б) дня различных сред вне и внутри сфер: капни воды в воздухе (Сф I), воздушные пузырьки в воде (Сф II), капни дихлорэтана в воде (Сф III). Соответствующие значения чисел усечения рядов приведены в табл. 2. Данные значения получены с помощью подхода 1 с использованием двойного контроля точности расчетов, а именно, относительная погрешность рассчитывалась дня двух функций — STF и модуля нормированного давления |р/р;п Необходимая точность считалась достигнутой при 5fiX = 1%, Применялся двойной контроль, поскольку были обнаружены параметры, при которых для одной из функций относительная погрешность могла принимать значения меньше 6fiX до достижения стабилизации искомой функции.
ntr
Рис. 4. Зависимость модуля нормированного давления [р/рт | от числа щг для конфигурации, представленной на рис. 2, в точках: 1 — Mi(10a, 0, 0); 2 — М2(0,а, 0); 3 — Мэ(0, 0.5a, 0). Система сфер: капли дихлорэтана в воде; основные параметры: к0а = 5, SI = 5а Fig. 4. The dependence of the normalized pressure modulus |p/p-m | on the truncation number щг for the configuration shown on Fig. 2 at the points: 1 — M1(10a, 0,0); 2 — M2(0,a, 0); 3 — М3(0, 0.5a, 0). Sphere system: dichloroethane drops in water; main parameters: k0a = 5, SI = 5a
Таблица 2. Число усечения рядов ntr, полученное с помощью подхода 1 при 5fix = 1 %, в зависимости от волнового радиуса к0а для рассеивающего слоя из 121 сферы (см. рис. 2) в точке МС для двух значений параметра 51: Сф I — капли воды в воздухе; Сф II — воздушные пузырьки в воде; Сф III — капли дихлорэтана в воде
Table 2. The series truncation number щг obtained using the approach 1 at Sfix = 1 % as function of the wave radius k0a fa a scattering layer of 121 spheres (see Fig. 2) at the point MC for two values of the parameter 51: Sph I — water drops in air; Sph II — air bubbles in water; Sph III — dichloroethane drops in water
к0а Ccj t> 1 Сф II Сф III Формула (6) Формула (7)
51 = 3а 51 = 5а 51 = 3а 51 = 5а 51 = 3а 51 = 3а 51 = 3а 51 = 5а
1 2 3 4 3 2 1 6 4 6
2 5 5 8 5 4 4 9 8 13
3 9 8 7 7 7 7 11 12 20
4 10 8 9 8 10 8 13 16 27
5 8 12 12 10 11 12 14 20 33
Из рис. 3 видно влияние основных геометрических и физических параметров системы, таких как:
• расстояние между центрами сфер 51, поскольку чем ближе сферы друг к другу, тем больше их взаимодействие;
• волновой радиус к0а = (ш/с0)а, так как при увеличении параметра к0а и фиксированном радиусе сферы а усиливается и воздействие внешнего поля на систему сфер, а при фиксированной частоте внешнего поля f увеличивается радиус самих сфер, а следовательно, и их отклик на внешнее воздействие;
• физические параметры сред вне и внутри сфер (плотность р и скорость звука с), поскольку отношение данных величин существенно влияет на характер рассеянного поля. Рассеяние капель воды в воздухе, как уже отмечалось, согласуется со случаем жестких сфер {\v/Хо ^ ж и pv/р0 ^ то), рассеяние от воздушных пузырьков в воде вне области резонанса {к0а ^ 0.013) — со случаем мягких сфер [10] (x0/xv ^ ж и p0/pv ^ то), а случай капель одной жидкой среды в другой является промежуточным (в рассмотренном примере \v/х0 ~ 0.6 и pv/р0 & 1.3).
Проведем анализ данных, представленных в табл. 2. Сравнение двух подходов к определению числа ntr для параметров системы 51, к0а, рис показало следующее:
• значения ntr для рассмотренных случаев Сф I-Сф III и разных 51 меньше числа ntr, полученного с помощью формулы (6), которая не зависит ни от физических свойств сред, ни от конфигурации системы. Таким образом, для волновых радиусов к0а < 5 формула (6) завышает значения ntr, что приводит к существенному увеличению времени расчета;
• все значения ntr для случаев Сф I-Сф III не превышают соответствующих значений, полученных по формуле (7), которая задает верхнюю границу числа ntr',
• логично предположить, что чем больше расстояние между сферами, тем меньше их взаимодействие, следовательно, тем меньше будет число ntr. Однако имеются параметры, для которых данное предположение не выполняется, например для Сф I при к0а = 1 и 5, а также для Сф III при к0а = 5, что затрудняет формулировку закономерности для числа ntr в зависимости от 51, Значит, в каждом случае необходимо проводить дополнительное исследование для определения числа усечения рядов при разложении;
• аналогично, чем больше волновой радиус к0а, тем больше число щг. Но и здесь имеются исключения — для Сф I значение при к0а = 5 и 51 = 3а меньше, чем при к0а = 3 и 4, а при 51 = 5а данные значения при к0а = 3 и 4 равны, что также затрудняет формулировку зависимости щг от к0а и требует проведения дополнительных исследований. Таким образом, в случае 1 < к0а < 5 необходимо использовать подход 1 для получения наименьшего значения щг-, при котором достигается заданная точность, поскольку эвристические формулы для данных волновых радиусов завышают значение числа щг, что приводит к существенному увеличению времени расчета. Для расчетов областей, когда требуется привлечение значительных вычислительных ресурсов, рекомендуется предварительно вычислить значение щг в нескольких характерных точках, например внешних, внутренних и граничных. Выбирается наибольшее значение числа щг из полученных дня данных характерных точек при условии, что оно не превосходит значения, вычисленного но формуле (7).
На рис. 5 и 6 представлены диаграммы распределения нормированного давления 1р/рш | для конфигурации капель воды в воздухе (см. рис. 2) в плоскостях, параллельных рассеивающему слою, при значениях параметров 51 = 3а и к0а = 1. Поскольку система симметрична относительно осей Ох и Оу, расчеты проведены только в первой четверти. Число расчетных точек плоскости 91 х 91. Для диаграмм на рис. 5 требуемая точность (5/гх = 1 %) достигается при щг = 3 (для всех расчетных точек плоскости), а на рис. 6 — при щг = 6 (предварительные расчеты в характерных точках), так как для получения стабильного решения внутренних точек требуется больше членов в разложении рядов.
Плоскость х = -10а (рис. 5, а) содержит монопольный источник излучения МБ. На распределение давления влияют волны, рассеянные системой сфер в обратном к падающей сферической волне направлении. В плоскости х = 10а (см. рис. 5, б) на распределение давления влияют волны, прошедшие сквозь рассеивающий слой. Распределение давления на самом слое (в плоскости, содержащей центры сфер х = 0) представлено на рис. 6. Капли воды в воздухе мало пропускают сквозь себя падающую волну, что б.лиз-
у/а у/а
Рис. 5. Диаграммы распределения нормированного давления [p/pm | вне рассеивающего слоя (см. рис. 2) при 51 = 3а и коа = 1 в первой четверти плоскостей х = — 10а (а) их = -10а (б) Fig. 5. Diagrams of the normalized pressure distribution lp/p-m | outside the scattering layer (see Fig. 2) which consists of water drops in the air at 51 = 3a and k0a = 1 in the first quarter of planes x = —10а (a) and x = —10а (6)
0 3 6 9 12 15
у/а
Рис. 6. Диаграмма распределения нормированного давления |р/ртI в плоскости х = 0 при 51 = 3а и коа = 1 в первой четверти
Fig. 6. Diagrams of the normalized pressure distribution Ip/p-m | in the plane x = 0 at 51 = 3a and k0a = 1 in the first quarter
ко к случаю звуконепроницаемых жестких сфер. На диаграммах хорошо видны зоны повышения и понижения давления как результат наложения падающей и рассеянных от системы сфер волн. Таким образом, рассеивающий слой в рассмотренном примере играет роль сложной двумерной дифракционной решетки.
Заключение
Дня задачи акустического рассеяния от множества звукопроницаемых сфер разных радиусов, произвольно расположенных в трехмерном пространстве, при воздействии акустическим нолем проведена оптимизация алгоритма построения матрицы коэффициентов повторного разложения по мультиполям с целью минимизации используемого объема оперативной памяти. Верификация оптимизированного алгоритма показана хорошее соответствие численным данным, полученным другими исследователями.
Анализ использования вычислительных ресурсов и результаты профилирования программного кода показали, что для волновых радиусов 1 < к0а < 5 наименьшее значение числа усечения рядов при разложении щг, для которого достигается заданная точность, дает подход 1. Подход 2 завышает значение щг, что приводит к существенным затратам машинного времени. Проведенная серия вычислительных экспериментов позволила получить количественную оценку ntr, которая может быть использована при исследовании ресурсоемких задач по акустическому рассеянию от звукопроницаемых сфер.
Диаграммы, построенные с помощью рассмотренного подхода к изучению рассеяния от множества сфер, в отличие от подходов, основанных па осредпеппых уравнениях, позволяют наглядно продемонстрировать картину распределения давления вне и внутри системы сфер, в том числе дня определения зон повышения и понижения давления.
В дальнейшем планируется проведение исследований акустического рассеяния дня различных конфигураций системы в более широком диапазоне изменения физических и геометрических параметров.
Благодарности. Работа выполнена при поддержке госзадания № 0246-2019-0052.
Список литературы
[1] Исимару А. Распространение и рассеяние волн в случайно-неоднородных средах. М.: Мир; 1981: 280.
[2] Сташкевич А.П. Акустика моря. Ленинград: Судостроение; 1966: 356.
[3] Алешин Н.П., Щербинский В.Г. Радиационная, ультразвуковая и магнитная дефектоскопии металлоизделий. М.: Высшая школа; 1991: 271.
[4] Применение ультразвука в медицине: Физические основы (под ред. К. Хилла). М.: Мир; 1989: 589.
[5] Демин И.Ю., Прончатов-Рубцов Н.В. Современные акустические методы исследований в биологии и медицине / Учеб. метод, пособие. Нижний Новгород: Издательство ННГУ; 2007: 121.
[6] Каллистратова М.А. Радиоакустическое зондирование атмосферы. М.: Наука; 1985: 197.
[7] Технология создания позиционируемого 3D звука. Адрес доступа: https: //www. ixbt. com/ multimedia/3dsound-tech.html (дата обращения 21.10.2021).
[8] Medwin H., Clay C.S. Fundamentals of acoustical oceanography. San Diego: Academic Press; 1998: 739.
[9] Duda R.O., Martens W.L. Range dependence of the response of a spherical head model. Journal of the Acoustical Society of America. 1998; 104(5):3048-3058. D01:10.1121/1.423886.
[10] Гринченко В.Т., Вовк И.В., Мацыпура В.Т. Основы акустики. Киев: Наукова думка; 2009: 867.
[11] Насибуллаева Э.Ш. Исследование рассеяния от звуконепроницаемой одиночной сферы при внешнем воздействии. Труды Института механики им. P.P. Мавлютова Уфимского научного центра РАН. 2017; 12(1):73-82. D01:10.21662/uim2017.1.011.
[12] Насибуллаева Э.Ш. Исследование акустического рассеяния от одиночной звукопроницаемой сферы. Многофазные системы. 2018; 13(4):79—91. D01:10.21662/mfs2018.4.012.
[13] Kapodistrias G., Dahl P.H. Effects of interaction between two bubble scatterers. Journal of the Acoustical Society of America. 2000; (107):3006-3017. DOhlO.1121/1.429330.
[14] Румелиотис Д.А., Котсис А.Д. Рассеяние звуковых волн на двух сферических телах, одно из которых имеет малый радиус. Акустический журнал. 2007; 53(1):38-49.
[15] Насибуллаева Э.Н1. Исследование акустического рассеяния от пары звуконепроницаемых сфер при внешнем воздействии. Многофазные системы. 2019; 1 I(1): 11 51. DOI: 10.21662/mfs2019.1.006.
[16] Gumerov N.A., Duraiswami R. Computation of scattering from N spheres using multipole reexpansion. Journal of the Acoustical Society of America. 2002; 112(6):2688—2701. DOI:W. 1121/1.1517253.
[17] Skaropoulos N.C., Yagridou H.D., Chrissoulidis D.P. Interactive resonant scattering by a cluster of air bubbles in water. Journal of the Acoustical Society of America. 2003; 113(6):3001—3011. DOLIO.1121/1.1572141.
[18] Насибуллаева Э.Н1. Численный анализ акустического рассеяния от звукопроницаемых сфер при внешнем воздействии. Вестник УГАТУ. 2021; 25(2):93—101. D01:10.54708/19926502_2021_2529293.
[19] Насибуллаева Э.Ш. Моделирование акустического рассеяния от коаксиальных сфер при внешнем воздействии. Многофазные системы. 2019; 14(2): 115^124. D01:10.21662/mfs2019.2.016.
[20] Кобелев Ю.А. К теории многократного рассеяния звуковых волн на сферических частицах в жидких и упругих средах. Акустический журнал. 2011; 57( I): 113 I 19.
[21] Владимиров B.C. Уравнения математической физики. М.: Наука; 1981: 512.
[22] Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука; 1974: 832.
[23] LAPACK — Linear Algebra PACKage. Available at: https://www.netlib.org/lapack
(accessed at October 21, 2021).
[24] Zhang S., Jin J. Computation of special functions. N.Y.: Wiley-Interscience; 1996: 717.
[25] Antoine X., Chniti C., Ramdani K. On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders. Journal of Computational Physics. 2008; 227(3):175 I 1771.
[26] Chew W.C., Jin J.M., Michielssen E., Song J. Fast and efficient algorithms in computational electromagnetics. Norwood: Artech House Antennas and Propagation Library; 2001: 932.
[27] Amamou M.L. A theoretical and numerical resolution of an acoustic multiple scattering problem in three-dimensional case. Acoustical Physics. 2016; 62(3):280—291. DOI: 10.1134/S1063771016030015.
[28] Насибуллаева Э.Ш. Определение числа членов при усечении рядов для численного решения задачи акустического рассеяния от множества звукопроницаемых сфер. Многофазные системы. 2020; 15(3 1):170 1S2. D01:10.21662/mfs2020.3.128.
Вычислительные технологии, 2022, том 27, № 2, с. 19-36. © ФИЦ ИВТ, 2022 ISSN 1560-7534
Computational Technologies, 2022, vol. 27, no. 2, pp. 19-36. © FRC ICT, 2022 elSSN 2313-691X
MATHEMATICAL MODELLING
DOI: 10.25743/ICT.2022.27.2.003 Simulation of acoustic scattering from a set of sound-permeable spheres in 3D space
Nasibullaeva Elvira Sh.
Mavlyutov Institute of Mechanics UFRC RAS, 450054, Ufa, Russia Corresponding author: Nasibullaeva Elvira Sh., e-mail: elvira@anrb.ru Received October 26, 2021, revised February 4, 2022, accepted February 11, 2022.
Abstract
The problem of acoustic scattering from a sound-permeable spheres set, arbitrarily located in 3D space, under the acoustic field action is considered. In the case of a spheres large number in the system configuration, as well as with an increase in their interaction, the problem solution requires significant computing resources and computer time. Therefore the algorithm optimization for constructing the reexpansion coefficient matrix is carried out in order to minimize the amount of RAM used, and the analysis of the costs of computing resources is carried out depending on the truncated series number. The optimized algorithm verification showed a good correspondence with the numerical data obtained by other researchers.
The purpose is to determine the pressure distribution outside and inside the spheres in a wide range of system parameters.
Methodology. When solving the Helmholtz equations, a numerical technique based on the fast multipole method was used. For the numerical solution, the method developed for the case of soundproof spheres has been adapted for the sound-permeable one.
As a finding of the conducted series of computational experiments, a quantitative estimate of the truncated series number was obtained and can be applied to solution of resource-intensive problems. The optimized calculation method for various values of system parameters is demonstrated for particular examples.
Originality. Diagrams constructed using the considered approach, in contrast to approaches based on averaged equations, have represented a complete pressure distribution picture outside and inside the spheres system.
Keywords: sound-permeable sphere, acoustic scattering, multipole expansion method, computational experiment, monopole radiation source, pressure distribution diagram.
Citation: Nasibullaeva E.Sh. Simulation of acoustic scattering from a set of sound-permeable spheres in three-dimensional space. Computational Technologies. 2022; 27(2):19—36. D01:10.25743/ICT.2022.27.2.003. (In Russ.)
Acknowledgements. The work was supported by the state assignment No. 0246-2019-0052.
References
1. Ishimaru A. Wave propagation and scattering in random media. 1st edition. N.Y.: Academic Press; 1978: 339.
2. Stashkevich A.P. Akustika morya [Acoustics of the sea]. Leningrad: Sudostroenie; 1966: 356. (In Russ.)
3. Aleshin N.P., Shcherbinskiy V.G. Radiatsionnaya, ul'trazvukovaya i magnitnaya defektoskopii metalloizdeliy [Radiation, ultrasonic and magnetic flaw detection of metal products]. Moscow: Vysshaya Shkola; 1991: 271. (In Russ.)
4. Physical principles of medical ultrasonics. Eds.: C.R. Hill. Chichester: Ellis Horwood; 1986: 495.
5. Demin I.Yu., Pronchatov-Rubtsov N.V. Sovremennye akusticheskie metody issledovaniy v biologii i meditsine [Modern acoustic research methods in biology and medicine]. Nizhni Novgorod: Izdadelstvo NNGU; 2007: 121. (In Russ.)
6. Kallistratova M.A. Radioakusticheskoe zondirovanie atmosfery [Radio-acoustic sounding of the atmosphere]. Moscow: Nauka; 1985: 197. (In Russ.)
7. Technology for creating positioned 3D sound. Available at: https://www.ixbt.com/multimedia/ 3dsound-tech.html (accessed at October 21, 2021). (In Russ.)
8. Medwin H., Clay C.S. Fundamentals of acoustical oceanography. San Diego: Academic Press; 1998: 739.
9. Duda R.O., Martens W.L. Range dependence of the response of a spherical head model. Journal of the Acoustical Society of America. 1998; 104(5) :3048-3058. D01:10.1121/1.423886.
10. Grinchenko V.T., Vovk I.V., Matsypura V.T. Osnovy akustiki [Basics of acoustics]. Kiev: Naukova Dumka; 2009: 867. (In Russ.)
11. Nasibullaeva E.Sh. Investigation of scattering from soundproof single sphere under external influence. Proceedings of the Mavlyutov Institute of Mechanics. 2017; 12(l):73-82. DOI: 10.21662/uim2017.1.011. (In Russ.)
12. Nasibullaeva E.Sh. The study of acoustic scattering from a single sound-permeable sphere. Multiphase Systems. 2018; 13(4):79-91. D01:10.21662/mfs2018.4.012. (In Russ.)
13. Kapodistrias G., Dahl P.H. Effects of interaction between two bubble scatterers. Journal of the Acoustical Society of America. 2000; (107):3006-3017. D01:10.1121/1.429330.
14. Roumeliotis J.A., Kotsis A.D. Acoustic scattering from two spheres, one with a small radius. Acoustical Physics. 2007; 53(l):33-43. D0110.1134/S1063771007010046.
15. Nasibullaeva E.Sh. Investigation of acoustic scattering from a pair soundproof spheres under external influence. Multiphase Systems. 2019; I I( I): 11 51. D01:10.21662/mfs2019.1.006. (In Russ.)
16. Gumerov N.A., Duraiswami R. Computation of scattering from N spheres using multipole reexpansion. Journal of the Acoustical Society of America. 2002; 112(6):2688-2701. D01:10.1121/l.1517253.
17. Skaroponlos N.C., Yagridou H.D., Chrissoulidis D.P. Interactive resonant scattering by a cluster of air bubbles in water. Journal of the Acoustical Society of America. 2003; 113(6):3001-3011. D01:10.1121/l.1572141.
18. Nasibullaeva E.Sh. Numerical analysis of acoustic scattering from soundpermeable spheres under external influence. Vestnik UGATU. 2021; 25(2):93-101. (In Russ.)
19. Nasibullaeva E.Sh. Numerical simulation of acoustic scattering from coaxial soundpenetrable spheres. Multiphase Systems. 2019; 14(2):115-124. D01:10.21662/mfs2019.2.016. (In Russ.)
20. Kobelev Yu.A. On the theory of multiple scattering of sound waves by spherical particles in liquid and elastic media. Akusticheskij Zhurnal. 2011; 57(4):443-449. (In Russ.)
21. Vladimirov V.S. Equations of mathematical physics. N.Y.: Marcell Dekker Incoporated; 1971: 427.
22. Korn G.A., Korn Th.M. Mathematical handbook for scientists and engineers. McGraw Hill Book Company; 1968: 943.
23. LAPACK — Linear Algebra PACKage. Available at: https://www.netlib.org/lapack (accessed at October 21, 2021).
24. Zhang S., Jin J. Computation of special functions. N.Y.: Wiley-Interscience; 1996: 717.
25. Antoine X., Chniti C., Ramdani K. On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders. Journal of Computational Physics. 2008; 227(3): 1754-1771.
26. Chew W.C., Jin J.M., Michielssen E., Song J. Fast and efficient algorithms in computational electromagnetics. Norwood: Artech House Antennas and Propagation Library; 2001: 932.
27. Amamou M.L. A theoretical and numerical resolution of an acoustic multiple scattering problem in three-dimensional case. Acoustical Physics. 2016; 62(3):280-291. D01:10.1134/S1063771016030015.
28. Nasibullaeva E.Sh. Terms number determination at the series truncation for the numerical solution of the problem of acoustic scattering from a sound-permeable spheres set. Multiphase Systems. 2020; 15(3 I): 176 1X2. DOI: 10.21662/mfs2020.3.128. (In Russ.)
