CC BY
6
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
2022, Т. 164, кн. 1 С. 68-84
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
ОРИГИНАЛЬНАЯ СТАТЬЯ
УДК 519.63
doi: 10.26907/2541-7746.2022.1.68-84
ЭФФЕКТИВНЫЙ ЧИСЛЕННЫЙ МЕТОД ДЛЯ ОПРЕДЕЛЕНИЯ ЗАХВАЧЕННЫХ МОД АКУСТИЧЕСКИХ ВОЛНОВОДОВ
Р.З. Даутов
Казанский (Приволжский) федеральный университет, г. Казань, 420008, Россия
Предложен эффективный приближенный метод определения всех захваченных мод уравнения Гельмгольца на основе метода конечных элементов и точных нелокальных краевых условий. Рассмотрен бесконечный двумерный канал с параллельными стенками на бесконечности, который может содержать препятствия произвольной формы, при этом предполагается, что частоты захваченных мод лежат ниже некоторого порогового значения. Предлагаемая дискретная задача представляет собой алгебраическую задачу на собственные значения для симметричных положительно определенных разреженных матриц, одна из которых зависит нелинейно от спектрального параметра. Разработан быстрый итерационный метод решения подобных задач. Приведены результаты численных расчетов.
Ключевые слова: акустический волновод, захваченная мода, дискретный и непрерывный спектр, метод конечных элементов, нелинейная спектральная задача
В настоящей работе предлагается и численно исследуется приближенный метод на основе метода конечных элементов (МКЭ) для решения задачи о распространении волн, описываемых уравнением Гельмгольца, в неограниченном в одном направлении двумерном акустическом волноводе.
В рамках классической волновой теории захваченные моды (trapped modes) представляют собой гармонические во времени колебания на некоторой четко определенной частоте, локализованные вблизи границы волновода или ловушеч-ной структуры в неограниченной области. Затухание по мере удаления от объекта захвата обычно является экспоненциальным. В зависимости от контекста захваченные моды известны как акустические резонансы, волны Рэлея-Блоха, краевые волны, поверхностные волны, квантово-механические связанные состояния и т. д. В настоящей работе мы рассматриваем захваченные моды в акустическом контексте, но можно провести прямые аналогии с соответствующими проблемами в волнах на воде и в квантовых волноводах [1—6].
Мы будем рассматривать волны в двумерном неограниченном канале (волноводе) П, который имеет параллельные стенки на бесконечности и может включать необязательно связное ограниченное препятствие ш произвольной формы. Рассмотрим гармонические во времени колебания с угловой частотой ш = Хе, где c - скорость волны. Амплитуда колебания поля давления в волноводе, и = u(x, y), удовлетворяет уравнению Гельмгольца Au + Хи = 0 в П\ш при соответствующих
Аннотация
Введение
однородных граничных условиях. Захваченная мода (или акустический резонанс) -это решение задачи с конечной энергией, то есть такое, что
J (и2(х,у) + \Уи(х,у)\2) йхйу < ж. (1)
П\ш
Положительное значение волнового числа А, для которого существует ненулевое решение уравнения Гельмгольца, называется частотой захваченной моды. Таким образом, математически частоты захваченных мод соответствуют собственным числам оператора Лапласа в неограниченной области, и их существование в решающей степени зависит от природы границ области и типа граничных условий.
Неограниченность волновода представляет основную трудность как при исследовании разрешимости задачи и определении числа решений, так и при разработке и теоретическом исследовании численных методов ее решения. Охарактеризируем кратко эти трудности.
В силу неограниченности области спектр а оператора Лапласа (—Д), рассматриваемого как самосопряженный, полуограниченный снизу оператор в пространстве ^(П^), не может быть чисто дискретным и содержит непрерывную компоненту ас = [Л1, то), где пороговое значение Л1 зависит от ширины канала и краевых условий на бесконечности. Если выполнено условие
\Уи\2 йх йу
А\ = ^- < Л1, (2)
«ещо} / 2
П\ш
то также не пуст дискретный спектр а^. В этом случае существует конечное число собственных чисел, включая А1 , определяющих частоты захваченных мод ниже порогового значения Л1. Остальные собственные числа, если они есть, включены в непрерывный спектр, что существенно усложняет задачу их поиска. В условии (2) V есть множество функций с конечной энергией (пространство Соболева Н 1(П\^)), удовлетворяющих краевым условиям Дирихле (если они имеются).
Если условие (2) не выполнено (то есть А1 = Л1), то дискретный спектр либо пуст, либо вложен в непрерывный спектр. В последнем случае нередко подобные задачи могут быть переформулированы так, что оператор новой задачи имеет конечное число собственных значений (являющихся также собственными числами исходной задачи) ниже нового порогового значения Л2. Например, в [7] был предложен прием постановки искусственных краевых условий Дирихле на плоскости геометрической симметрии канала для создания такого искусственного порогового значения Л2 (см. также [8]). Этот прием нашел применение в многочисленных исследованиях и привел к построению примеров захваченных волн в разнообразных задачах математической физики (см., например, [9-11].
В настоящей работе мы предполагаем, что исходная задача уже сформулирована таким образом, что условие (2) выполнено, и остановимся на самом численном методе определения всех захваченных мод и соответствующих им частот ниже порогового значения.
Для исследования захваченных мод применяются различные численные схемы. Имеется ряд методов, которые используют особую форму препятствий и приводят к эффективным методам вычислений. Отметим, например, статью [1], в которой предложен метод, использующий линейную комбинацию соответствующим образом модифицированных мультипольных потенциалов, а также [12].
Трудности, с которыми сталкиваются при применении методов общего назначения, связаны с необходимостью обрезать бесконечную область, применяя соответствующее граничное условие на некоторой конечной фиктивной границе. Однако простые схемы могут порождать паразитные резонансы, устранение которых может представлять сложную задачу при автоматическом и эффективном исследовании широкого пространства параметров. Среди методов общего назначения отметим метод граничных элементов, основанный либо на вещественной, либо комплексной функции Грина [13, 14]. В этих работах авторы методом граничных элементов исходную задачу первоначально сводят к алгебраической задаче на собственные значения вида Л(Л)х = 0, матрица которой нелинейно зависит от спектрального параметра Л, а затем предлагают решать уравнение (¿е1;(Л(А)) = = 0 . Отметим также метод разложений по собственным функциям во внутренней области [15].
Предлагаемый нами метод имеет две составляющих:
1) первоначально известным способом исходная спектральная задача в неограниченной области эквивалентно сводится к задаче в ограниченной области с использованием нелокальных граничных П1М операторов и стандартно аппроксимируется методом конечных элементов;
2) для решения полученной на первом шаге алгебраической задачи на собственные значения, которая имеет вид К и + Б (Л)и = ЛМи ,1 предлагается итерационный метод, позволяющий определить все ее решения ниже порогового значения. Разреженные симметричные положительно определенные матрицы жесткости К и масс М имеют большую размерность, «граничная» неотрицательно определенная матрица Б на порядок более разрежена, чем матрица К, непрерывно и монотонно зависит от спектрального параметра Л.
В настоящей работе мы ограничиваемся численным исследованием точности предложенного метода конечных элементов и эффективности предложенного итерационного метода. Результаты проведенных вычислительных экспериментов позволяют утверждать, что метод прост, эффективен и пригоден для определения всех захваченных мод и соответствующих им частот ниже порогового значения для волноводов общего вида.
1. Формулировка задачи
С учетом вышесказанного относительно существования захваченных волн рассмотрим задачу поиска антисимметричных мод, считая, что геометрия задачи в декартовых координатах х, у имеет ось симметрии у = 0, а исходная задача сразу формулируется в половине области [7].
1.1. Исходная задача в неограниченной области. Пусть П = {(х,у) : 0 < у < 7(х), х € Я} - верхняя половина волновода, Г = {у = 7(х), х € Я} -его твердая граница, Го = {(х, 0) : х € Я} - линия симметрии (7 есть постоянная функция в случае волновода с прямолинейными стенками). Если в волноводе имеется препятствие (необязательно связное), то пусть ш обозначает его часть, расположенную в П, а дш - границу ш. Будем предполагать, что ш С {(х,у) : 0 < у < 7(х), \х\ < а}, а также, что х) = ! при \х\ > а. Функция 7(х) и граница препятствия считаются кусочно-гладкими (липшицевыми). Таким образом, = П\ш - область, занятая акустической средой.
1 Подобные задачи будем называть для краткости нелинейными задачами на собственные значения, а задачи вида Ах = ХЫх — линейными.
Определим два участка границы дО« области О« (несвязные): Гд = Го\д^ -точки дО«, лежащие на оси симметрии, и Гк - множество остальных точек дО« (лежащих на твердой стенке Г и границе препятствия). Определим также пороговое значение Л = (п/2й)2 .
Исходная задача формулируется следующим образом: найти ненулевую функцию и и число А € (0, Л) такие, что
Ди + Аи = 0, (х,у) € О«, (3)
ди
— =0, (х,у) € Гк, (4)
и = 0, (х,у) € Гд, (5)
и ^ 0, х ^ \у\ < й. (6)
Здесь д/дп - производная по внешней нормали п = п(х) к границе О«.
Из-за возможного сингулярного поведения градиента решения вблизи особенностей границы препятствия задачу (3)-(6) следует понимать в обобщенном смысле. В связи с этим введем дополнительные обозначения. Пусть Ь2(С) есть стандартное пространство Лебега со скалярным произведением и нормой
(и, у)п = ! и(х) • V (х) йх, ||и||ь2(д) = (и, и)Д2, д
Н 1(В) есть пространство Соболева с нормой || • ||#1 (д) (см. (1)):
1М1я1(д) = J(и2 + \уи\2) йхйу. д
Через V« обозначим множество функций из Н 1(Ото), удовлетворяющих условию Дирихле (5). Отметим, что функции из V« также удовлетворяют условию (6).
Под решением задачи (3)-(6) будем понимать пару (и, А) € ^«\{0} х (0, Л) такую, что
(V«) (Уи, Vv)^ = А(и^)Пте VV € V«,.
1.2. Задача в ограниченной области. Сведем эквивалентно задачу (V«) к задаче в ограниченной области, используя операторы точных нелокальных краевых условий (см., например, [16, 17]). При получении новой задачи будем следовать [18, 19] (см. также [20]). Для этого фиксируем Ь > а и введем фиктивные границы
Г± = {(±Ь,у):0 < у < й},
и представим О« в виде объединения непересекающихся областей О_, О и О+ так, что средняя область О содержит препятствие,
О = {(х,у) € О« : \х\ < Ь, 0 <у < й},
а области О± имеют прямолинейные границы и полуограничены. Положим
Vо = {V € V« : v(x, у) = 0, (х, у) € О}.
Пусть фиксирована пара (и, А) € V« х (0, Л). Функцию и\ назовем метагармо-ническим продолжением функции и\^ в область О_ У О+, если
Дих + Аих = 0, (х, у) € О_и О+,
ди\(х,й) . .
их(х, 0)=0, ^ , ' =0, \х\ > Ь, ду
их = и, (х,у) е г_иГ+,
и\ ^ 0, х ^ \у\ < (.
Легко видеть, что и\ удовлетворяет интегральному тождеству
(Уих, и = Х(их,у)п_и п+ У V е Уо,
и просто находится с помощью метода разделения переменных:
¿2 с-(и) ехр(Vх^ - X (х + Ь)) фп(у), (х, у) е О-
и\(х,у) = <
п=1
(7)
^с+(и) ехр(^/ Vп - X (Ь - х)) фп(у), (х, у) е О+,
-п=1
где (фп(у),^п) - решения задачи Штурма-Лиувилля —ф" = иф, ф(0) = ф'(() = 0, ип = (22(1^' фп(у) = *ш("пу\ с±(и) = 2 ! и(х,у)фп(у) (у. (8)
г±
Отметим, что величины -у/^ — X вещественны, поскольку X е (0, Л), Л = VI. Определим основное гильбертово пространство функций
У = {V е Н1(О): v\y=o = 0}, (и^)у = (Уи, ^)ып). (9)
Отметим, что на У нормы и ||и||у = (и, и)у эквивалентны в силу нера-
венства Фридрихса ||и||£2(о) < С ||и|V .
Определим также функционал в : К х У х У ^ К:
( /у2_ X
в(Х, = ^ -2- (с-п(и) с~п^) + с+п(и) с+ п^))
п=1 2
и введем обозначение
а^,и^) = (Уи, Vv)^ + (11)
Рассмотрим искомую задачу: найти пары (и, X) е У\{0} х (0, Л) такие, что (V) а^,и^) = X(u,v)q У V е У.
Она эквивалентна исходной задаче в следующем смысле.
Теорема 1. Пусть (и*(х),\) -решение исходной задачи (Рж), и(х) - сужение и*(х) на область О. Тогда (и(х)^) - решение ззадачи (V). Обратно: пусть (и(х)^) - решение задачи (V), функция и*(х) в области О совпадает с и(х), а в области О— У О+ равен его метагармоническому продолжению и\(х). Тогда (и*(х)^) - решение исходной задачи (Рж).
Доказательство. Пусть (и*(х)^) - решение (Рж), и = и*\п. Выбирая в (Рж) функцию v равной нулю в О, получаем, что и* = и\ в О— У О+. Поэтому тождество (Рж) можно переписать в следующем виде:
(Уи, Vv)n + в^,и^)= X(u,v)Q У v е Уто, (12)
где
s(X,u,v) = (Vux, Vv)n_u п+ - X(u\,v)n_u п+. (13)
После интегрирования по частям в первом слагаемом s будем иметь
s(X,u,v) = -(Aux + Xux,v)Q_u q+ +
f dux j f dux . f dux j f dux , + vdy - vdy = —— vdy - —— v dy.
J dx J dx J dx J dx
Г_ Г + Г_ Г+
Используя здесь формулу (7), придем к (10). Учитывая (10) в (12), получим, что (u,X) есть решение задачи (V), поскольку тождество (12) не изменится, если пространство Vœ заменить на V.
Обратно: пусть (u(x),X) - решение задачи (V). По u(x) определим ее метагар-моническое продолжение ux(x) и определим функцию u*(x), которая в области О совпадает с u(x), а в области О- У О+ равен ux(x). Ясно, что u* G Vœ. Согласно (13) справедливо равенство
s(X,u,v) = (Vu*, Vv)Q_u q+ - X(u*,v)Q_u q+.
С учетом этой формулы в (V) после объединения интегралов получим, что (u*(x),X) - решение исходной задачи (Vœ). □
2. Существование и число решений
Отметим ряд свойств функционала a, определенного формулой (11).
Лемма 1. При каждом X G [0, Л] билинейный функционал a является:
a) симметричным, то есть a(X,u,v) = a(X,v,u) для любых u, v G V;
b) положительно определенным и ограниченным. Точнее, найдется положительная постоянная C такая, что \\u\\V < a(X,u,u) < C\\u\\V для любого u G V.
Кроме того, при фиксированном u G V отношение Рэлея X ^ R(X, u) = = a(X,u,u)/(u,v)a является непрерывной и невозрастающей функцией на [0, Л], R(0,u) > c> 0 с постоянной c, не зависящей от u.
Доказательство. Свойство a) очевидным образом следует из определения a, как и оценка снизу в b), поскольку s(X,u,u) > 0. Докажем оценку сверху в b). Обозначим через H 1/2(Г -) пространство следов функций из H 1(О) на Г- . Хорошо известно (см., например, [21, p. 135]), что
œ
M2hi/2(г_) п (c-n(u))2 < C\\u\\l.
n=1
Аналогично, \\u\\hi/2(r+) < C\\u\\V . Искомая оценка теперь следует из оценки
a(X,u,u) < a(0, u, u) < — ^ n ((c n(u))2 + (c+n(u))2).
Свойства отношения Рэлея вытекают из свойств функций Л ^ \]— X и неравенства Фридрихса. □
Запишем задачу (V) в операторной форме в виде
А(Л)и = ЛМи, Л е (0, Л), и е V\{0}, (14)
где, согласно лемме 1, А(А) : V ^ V - самосопряженный ограниченный оператор, M : V ^ V - компактный оператор (в силу компактности вложения V С L2(Q)).
Непрерывность и монотонность оператор-функции А ^ А(А) позволяет достаточно просто как исследовать разрешимость задачи (14), так и определить число решений.
C этой целью определим обычную задачу на собственные значения, линейно зависящую от спектрального параметра. При фиксированном А е [0, Л] рассмотрим задачу определения пар (¡,u), ¡л = ¡(А), u = и(А) таких, что
A(A)u = ¡Mu, ¡> 0, u е V\{0}. (15)
Эта задача имеет счетное множество положительных собственных значений конечной кратности, которые занумеруем в невозрастающем порядке:
0 < ц(А) < ¡2(А) < ■■■ < ¡п(А) < ■■■ , lim ¡¡п(А) = ж,
n—
указывая каждое собственное значение с учетом его кратности. Собственные функции uk (А), к = 1, 2,..., соответствующие этим собственным значениям, образуют полную систему в пространстве V. Согласно теореме Рисса-Фишера
¡к(А) = min max И(А, v), к = 1, 2,..., (16)
VkCV veVk
где минимум берется по всем к -мерным подпространствам пространства V .
Ясно, что (А^(А)) есть решение задачи (14) тогда и только тогда, когда А является решением уравнения
¡(А) = А, А е [0,Л]. (17)
Теорема 2. Пусть волновод имеет прямолинейные стенки, ф\(у) определяется формулой (8) и выполнено условие
Jdi(y) dxdy> 0, Di(y) = \ф\(у)\2 — Л2\ф-1(у)\2. (18)
oj
Тогда .задача (V) имеет решение, причем число ее собственных чисел ( c учетом кратности) равно числу собственных чисел задачи A^)u = ¡Mu (c учетом кратности) меньших, чем Л.
Доказательство. Рассмотрим уравнение
¡k(А) = А, А е [0,Л], к > 1, (19)
где ¡k (А) определяется формулой (16). С учетом свойств отношения Рэлея, указанных в лемме 1, несложно проверить, что функция А ^ ¡k(А) является непрерывной и невозрастающей функцией на [0, Л] (см., например, [22]), причем R(0,u) > const > 0. Поэтому уравнение (19) имеет единственное решение, если ¡k (Л) < Л. При к = 1 это условие будет выполнено, если, например, R(A, ф\) < Л, иначе, если
выполнено условие Ii = (ф\, ф{)а — Л а(Л, ф1, ф1) > 0.
d
Так как система {фkортогональна, j Di(y) dy = 0, то в(Л, ф1, ф{) = 0 и
о
Ii = (ф1,фг)п — Л(Уфъ Уф1 = —J Di(y) dxdy =
Q
= —J Di(y) dxdy + J Di (y) dxdy = j Di(y) dx dy > 0,
По ш
где По = {(x,y) : \x\ < L, 0 < y < d}. Следовательно, пл > 1 при условии (18). □
Отметим, что критерий разрешимости (18), как следует из доказательства, не зависит от краевых условий функции фх. Он был получен другим методом в [11]. В рассматриваемом нами случае, согласно (8),
= А2 J сов(2^1у) ё,хё,у = А2 J сов ^.хё,у > 0.
Легко видеть, что это условие выполнено, например, если ш имеет ненулевую толщину и ш = {(х, у) : х € [—а, а], 0 < ^1(х) < у < Ч2(х) < ¿}.
3. Метод решения задачи
3.1. Дискретизация задачи. Пусть к > 0 есть малый параметр. Для дискретизации задачи (V) используем метод конечных элементов на основе линейных треугольных элементов. Пусть Тн = {ех, е2,. .., еп} - семейство треугольников (конечных элементов), образующих, как это принято в МКЭ, правильную триангуляцию области О (то есть два треугольника либо не имеют общих точек, либо имеют общую вершину или общую сторону),
Он = е, к = тах &ат(е).
ееТн
Множество вершин треугольников из Тн будем называть сеткой узлов на О. Все узлы пронумеруем в некотором порядке, через а^ обозначим координату г-го узла, через шн обозначим множество узлов сетки, не лежащих на оси абцисс (линии симметрии волновода). Считаем, что любые угловые точки О являются узлами сетки, прямолинейные участки границы О являются объединением сторон элементов.
Пусть Р\ = {ох + в2х + сзу, е^ € Я} есть множество линейных полиномов. Определим пространство конечных элементов Ун, аппроксимацию пространства V:
Ун = {ун € С (Он) : ук\е € Рь е € Тн, Ун\у=о = 0}.
Базис Лагранжа {^¿}1=1 в нем определяется так, что ^¿(а^) = , 1 < г, у < N, где N = ^ - число элементов шн. Функция из ин € Ун имеет представление
ин(х,у) = ^2 иг¥г(х,у), и = ик(а^). (20)
N ¿=1
Вектор-столбец ин = (их,щ,... ,UN)Т принято называть вектором узловых параметров (значений) ин. Усекая бесконечную сумму, определим функционал
т ^ х/туП — А
(А, и, у) = -2-(с-п(и) с-п(у) + с+п(и) с+ п(у)). (21)
п=1 2
Дискретный аналог задачи (V) формулируется следующим образом: найти пары (uh,АNm) € Ун\{0} х (0, Л) такие, что
(Рн) (Уин, Уу)о_к + 8т^т,ин,у) = АNm(uh,v)nh У у € Ун.
Теорема 2 применима и к дискретной задаче (Рн), поскольку границы Г± , на которых определяются коэффициенты Фурье с±, аппроксимируются точно.
Сформулируем задачу (Рн) на матричном языке. Положим у = ^ в (Рн) и используем разложение (20) вместо ин . Придем к матричной задаче: найти пары (ин, АNm) € Я1^\{0} х (0, Л) такие, что
AN ин + SNm(АNm)Uн = АNmMNUн. (22)
Здесь симметричные разреженные матрицы определяются соотношениями
An = {(VVj, ^ViW}Nj=v MN = {(vj,ViW }Nj=v
SNm(^) =ч ^ dVVn Л ^c n(^j ) С n(Vi)+ C+n (Vj) C+n(Vi))
N
/ иП
- {c-n(Vj ) C-n(Vi ) + C+n (Vj ) C+n\
U=1 ) i,j=i
Матрицы An и Mn - матрицы жесткости и масс, являются стандартными в МКЭ, точно и эффективно вычисляются. Матрицу SNm(X) удобно вычислять следующим образом. Один раз точно вычисляются матрицы
S± = {sik = c±Vi), 1 < i < N, 1 < к < M},
и при заданном Л определяется диагональная матрица D(X) размера M х M с элементами dyVf-X/2. Тогда SNm(Л) = S—D^)S—T + S+D^)S+T.
3.2. Метод PINVIT для линейной задачи. Итерационный метод NPINVIT, который будет определен далее для решения нелинейной задачи (22), является обобщением метода, известного как предобусловленный метод обратных итераций или предобусловленный градиентный метод. Метод PINVIT предназначен для определения группы m > 1 младших собственных чисел и соответствующих им векторов линейной задачи Ax = ЛMx с симметричными и положительно определенными матрицами A и M размера N. В нем собственные пары определяются последовательно друг за другом.
Будем считать, что предобуславливатель B есть симметричная положительно определенная матрица такая, что 7i(Bx,x) < (Ax,x) < 72(Bx,x) для любого x € Rn , итерационный параметр т = 2/ (71 + 72).
Опишем вычисления собственной пары с номером s при условии, что s — 1 первая пара (Лi, xi) уже найдена [24, р. 485]. Используем обозначения: (Л, x) - текущая итерация, x ортогонален найденным векторам, ||x||m = 1; (Л,з:) и/или (Л, sc) - промежуточная итерация; (Л',x') - следующая итерация собственной пары (Лs,xs),
112
R(x) = (Ax,x)/(Mx,x) - отношение Рэлея. Здесь ||x||m = (x,x)M , (x,y)M = = (Mx,y).
I. Начальная итерация. Выберем M-ортогональный найденным векторам xi, i = 1,..., s — 1, вектор x, ||x||m = 1. Положим Л = (Ax, x);
II. Итерационный цикл. До выполнения условия сходимости находим:
1. r = Ax — ЛMx, Л = B— 1 r;
s — 1
2. v = Л aixi, ai = (C,xi)M;
i=i
3. Л = x — tv, x' = Л/||Л||м;
4. Л' = (Ax',x').
Здесь на шаге 2 (при s > 2) осуществляется исчерпывание из вектора v найденных ранее собственных векторов. Этот шаг позволяет поддерживать итерации ортогональными найденным собственным векторам. В [24, р. 485] доказано, что этот метод имеет геометрическую скорость сходимости, зависящую от распределения собственных чисел и отношения 71/72. Кроме того, итерационные приближения к собственным числам образуют монотонно убывающую последовательность. Все это делает метод привлекательным для сеточных задач большой размерности, когда удается определить B так, чтобы отношение 71/72 не зависело от шага сетки.
3.3. Итерационный метод для нелинейной задачи. Рассмотрим обобщение итерационного метода, рассмотренного выше, для определения всех собственных пар (х, А) € И1^ х Я+ следующей нелинейной задачи на собственные значения:
А(А)х = АМх, А € (0, Л), \\х\\м = 1. (23)
Его можно рассматривать и как обобщение метода, предложенного в статье [22], для определения первого собственного числа задачи вида (23).
Будем считать, что М и А(А) (при каждом А € [0, Л]) являются симметричными и положительно определенными матрицами большой размерности, отношение Рэлея А ^ Я(А,х) = (А(А)х,х)/(Мх,х) непрерывно и не возрастает на [0, Л] при каждом х € RN, \\х\\м = 1. Кроме того, предположим, что имеются т собственных чисел (с учетом кратности) следующей задачи:
А(Л)х = /Мх, / € (0, Л), \\х\\м = 1. (24)
Обозначим через х+ соответствующие собственные векторы.
Рассмотренная ранее задача (22) является частным случаем (23). Как и при доказательстве теоремы 2, убеждаемся, что при сделанных предположениях задача (23) имеет ровно т решений.
Будем считать, что выбран симметричный положительно определенный предо-буславливатель В, а также, что известны функции 71 и 72 такие, что
71(А)(Вх,х) < (А(А)х, х) < 72(А)(Вх,х) Ух € RN, А € [0, Л]. (25)
Как и в линейном случае, выберем итерационный параметр т(А) = 2/(^х(А) +
+ 12(А)).
Предлагаемый метод NPINVIT состоит из двух шагов:
1. Предварительный шаг. Находим все решения (p+¿,x+¿), г = 1, 2,...т, задачи (24);
2. Основной шаг. Последовательно, друг за другом определяем собственные пары А, х¿), г =, 2,..., задачи (23).
Опишем шаг 2, считая, что в — 1 первая собственная пара уже найдена.
I. Начальная итерация. Положим х = х+ и найдем А, решая уравнение
(и1): А = (А(А)х, х), А € [0, Л].
II. Итерационный цикл (х, А) ^ (х', А'). Находим
1. г = А(А)х — А Мх, у = В-1г, х = х — т(А)у;
8-1
2. ж = X — а3х3 , где величины а3 определяются из равенств = 0,
3 = 1
г =1, 2,..., в — 1, г = А(А)х — АМх;
3. х' = х/\\х\\м;
4. А' найдем как решение уравнения (^2) : А' = (А(А')х',х'). Отметим два отличия этого метода и метода РШУ1Т.
1) Для определения нового приближения к собственному числу требуется решать скалярные уравнения (и{), (^2). При А(А) = К + Б (А) имеем (А(А)х,х) = = (Кх, х) + (Б(А)х, х). Вектор Кх вычисляется ранее при определении невязки. На итерациях решения скалярного уравнения необходимо вычислять лишь (Б(А)х, х) .
И) Поскольку собственные векторы нелинейной задачи не ортогональны, то исчерпывание ранее найденных векторов нельзя провести, как в методе РШУ1Т. Пусть Т(А) = А(А) — АЫ. Тогда на шаге 2 получаем
з-1
(г,х£) = (Т (Х)х,х£) = (Т (Х)х£,х) = (т (Х)х£,х — а.3 х^ =0. (26)
3 = 1
Таким образом, для определения коэффициентов а приходим к системе уравнений 8-1
х з , х £ J аз
(Т (А)хз, х£) аз = (Т (Х)х£, х), г = 1, 2, . . . , в — 1.
3 = 1
Если А(А) не зависят от А, то из (26) получим, что
(А£ — Х)(Ых£ ,х) = 0, г = 1, 2,..., в — 1,
то есть новое приближение х' является Ы-ортогональным всем ранее найденным векторам и методы РШУ1Т и МРШУ1Т совпадают, если А(А) не зависят от А.
4. Результаты вычислений
Работоспособность предложенного метода изучим на примерах решения ряда задач, варьируя основные параметры метода. Во всех вычислениях предобуслав-ливатель В в (25) выбирался одинаково как В = АN + Б^т(Ао), Ао = Л/2 (см. (22)). Сравнивая функции ^— А и ^— Ао, получаем
71(а) = * л т™ № — А)/("1 — ао),
1<и<т
12(х) = — а)/(^1 — а0).
\/ 1<и<т
Итерации ОТШУ1Т заканчивались при условии |Ак+1 — Ак\/Ак+1 < 10-6, где Ак есть к-я итерация А.
4.1. Одно цилиндрическое тело в прямолинейном канале. Рассмотрим прямолинейный канал с полушириной с! = 1, в котором расположено круговое препятствие радиуса а = 0.1. В этом случае имеется одна захваченная мода, причем при малых а для частоты справедлива формула [1]
А - (п/2(1 — п4а4/8))2. (27)
Расчетную область О заключим в прямоугольник [—Ь,Ь] х [0,1], Ь = 0.175. На подробной сетке с числом узлов N = 5529 было получено А = 2.461906, при этом абсолютная погрешность формулы (27) оказалась равной 1.3 • 10-4. При изменении N — 50 + 5500 число итераций возрастало от 5 до 6. Грубая сетка узлов, график собственной функции на твердой стенке представлены на верхних рисунках рис. 1. Из графика х ^ и1(х, !) видно, что собственная мода медленно убывает на бесконечности (решение из малой области О было метагармонически продолжено во внешность расчетной области).
Предложенный нами метод зависит от двух основных расчетных параметров: числа узлов сетки N и числа удерживаемых Фурье гармоник т. Для определения зависимости погрешности собственных чисел от числа узлов сетки N были проведены вычисления на последовательности сеток. Результаты представлены в табл. 1.
-0.6 -0.4 -0.2 0 0.2 0.4 0.6
1од N
N=5529 А=2.4619
-5 0 5 10 15 X
N=1492
1 2 3 4 5 6 7 8 т
0
10
10
10
10
10
10
10
10
2
10
Рис. 1. Задача 4.1. Грубая триангуляция области П и график захваченной моды на твердой стенке (верхние рисунки). Зависимость погрешности N ^ \Х — (т = 10) и т ^ \Х — А^т| (N = 1492) (нижние рисунки)
Табл.1
Зависимость первого собственного числа от N при т = 10
N 50 117 256 433 605 1492 5529
Хмт 2.46414 2.46240 2.462023 2.461954 2.461935 2.461913 2.461906
Соответствующий график погрешности N ^ А — А^т в двойных логарифмических координатах представлен на левом нижнем рис. 1. Он соответствует ожидаемой оценке |А — А^т\ < C/N. Зависимость погрешности от т приведена на правом нижнем рис. 1. Из графика следует, что погрешность экспоненциально убывает с ростом числа Фурье гармоник в определении граничного оператора.
4.2. Два цилиндрических тела в прямолинейном канале. Рассмотрим прямолинейный канал с полушириной с! = 1, в котором расположены два круговых препятствия радиуса а = 0.5, расстояния между которыми равно 2. В этом случае имеются две захваченные моды, симметричная и антисимметричная относительно оси симметрии х = 0. Графики собственных функций на твердой стенке и собственные числа представлены на рис. 2. Расчетная область П была заключена в прямоугольник [—Ь,Ь] х [0,1], Ь = 1.6. При изменении N ~ 120 + 10000 число итераций возрастало от 4 до 5. При всех N достаточно было выбрать т = 3.
4.3. Два цилиндрических тела в каверне. Рассмотрим расчетную область, изображенную на первом рис. 3. Радиус кругов равен 0.6, Ь = 1.8. В этом случае имеется 3 решения: первая и вторая собственные функции являются четными, а третья - нечетной относительно оси симметрии х = 0. На всех сетках и для
-1 — -20
1.5
1
0.5 0
-0.5 -1
Рис. 2. Задача 4.2. Графики собственных функций на твердой стенке
3.5 3 2.5 2 1.5 1
0.5 0
N=197 N1=304
3.5 3 2.5 2 1.5 1
0.5 0
Рис. 3. Задача 4.3. Грубая триангуляция области П, собственные числа и функции
N=1120 А1 = 1.6873
N=1120 А2=2.3157
0
10
20
20
N=2552 А1=0.18737
5
N=2552 Л3=2.3088
N=2552 Л3=2.3088
0
5
всех мод достаточно было выбрать число Фурье гармоник т = 3. Для определения А1 понадобилось 4 итерации, А2 - 7 итераций, А3 - 14 итераций.
Приведенные результаты вычислений позволяют с уверенностью утверждать, что предложенный в работе вариант метода конечных элементов и итерационный метод его реализации являются подходящими для вычисления захваченных мод акустических волноводов.
Благодарности. Работа выполнена за счет средств Программы стратегического академического лидерства Казанского (Приволжского) федерального университета («ПРИОРИТЕТ-2030»).
Литература
1. Callan M.A., Linton C.M., Evans D.V. Trapped modes in two-dimensional waveguides // J. Fluid Mech. - 1991. - V. 229. - P. 51-64. - doi: 10.1017/S0022112091002938.
2. Evans D.V., Linton C.M., Ursell F. Trapped mode frequencies embedded in the continuous spectrum // Q. J. Mech. Appl. Math. - 1993. - V. 46, No 2. - P. 253-274. - doi: 10.1093/qjmam/46.2.253.
3. Evans D.V., Porter R. Trapping and near-trapping by arrays of cylinders in waves // J. Eng. Math. - 1999. - V. 35. - P. 149-179. - doi: 10.1023/A:1004358725444.
4. Exner P., Seba P. Bound states in curved quantum waveguides //J. Math. Phys. - 1989. -V. 30, No 11. - P. 2574-2580. - doi: 10.1063/1.528538.
5. Postnova J., Craster R.V. Trapped modes in elastic plates, ocean and quantum waveguides // Wave Motion. - 2008. - V. 45, No 4. - P. 565-579. - doi: 10.1016/j.wavemoti.2007.11.002.
6. Caspers F., Scholz T. Measurement of trapped modes in perforated waveguides // Part. Accel. - 1989. - V. 51. - P. 251-262.
7. Evans D.V., Levitin M., Vassiliev D. Existence theorems for trapped modes //J. Fluid Mech. - 1994. - V. 261. - P. 21-31. - doi: 10.1017/S0022112094000236.
8. Linton C.M., McIver M., McIver P., Ratcliffe K., Zhang J. Trapped modes for off-centre structures in guides // Wave Motion. - 2002. - V. 36, No 1. - P. 67-85. - doi: 10.1016/S0165-2125(02)00006-9.
9. Linton C.M., McIver P. Embedded trapped modes in water waves and acoustics // Wave Motion. - 2007. - V. 45, No 1-2. - P. 16-29. - doi: 10.1016/j.wavemoti.2007.04.009.
10. Nazarov S.A. Properties of spectra of boundary value problems in cylindrical and quasicylindrical domains // Maz'ya V. (Ed.) Sobolev Spaces in Mathematics II. International Mathematical Series, V. 9. - N. Y.: Springer, 2009. - P. 261-309. - doi: 10.1007/978-0-387-85650-6_12.
11. Nazarov S.A. Variational and asymptotic methods for finding eigenvalues below the continuous spectrum threshold // Sib. Math. J. - 2010. - V. 51, No 5. - P. 866-878. -doi: 10.1007/s11202-010-0087-3.
12. Evans D.V., Linton C.M. Trapped modes in open channels // J. Fluid Mech. - 1991. -V. 225. - P. 153-175. - doi: 10.1017/S0022112091002008.
13. McIver M., Linton C.M., McIver P., Zhang J., Porter R. Embedded trapped modes for obstacles in two-dimensional waveguides // Q. J. Mech. Appl. Math. - 2001. - V. 54, No 2. - P. 273-293. - doi: 10.1093/qjmam/54.2.273.
14. Sargent C.V., Mestel A.J. Trapped modes of the Helmholtz equation in infinite waveguides with wall indentations and circular obstacles // IMA J. Appl. Math. - 2019. - V. 84, No 2. - P. 312-344. - doi: 10.1093/imamat/hxy060.
15. Levitin M., Marletta M. A simple method of calculating eigenvalues and resonances in domains with infinite regular ends // Proc. R. Soc. Edinburgh, Sect. A: Math. - 2008. -V. 138, No 5. - P. 1043-1065. - doi: 10.1017/S0308210506001144.
16. Keller J.B., Givoli D. Exact non-reflecting boundary conditions //J. Comput. Phys. -1989. - V. 82, No 1. - P. 172-192. - doi: 10.1016/0021-9991(89)90041-7.
17. Givoli D. Non-reflecting boundary conditions //J. Comput. Phys. - 1991. - V. 94, No 1. -P. 1-29. - doi: 10.1016/0021-9991(91)90135-8.
18. Dautov R.Z., Karchevskii E.M. On a spectral problem of the theory of dielectric waveguides // Comput. Math. Math. Phys. - 1999. - V. 39, No 8. - P. 1293-1299.
19. Dautov R.Z., Karchevskii E.M. Existence and properties of solutions to the spectral problem of the dielectric waveguide theory // Comput. Math. Math. Phys. - 2000. -V. 40, No 8. - P. 1200-1213.
20. Dautov R.Z., Karchevskii E.M., Kornilov G.P. A numerical method for finding dispersion curves and guided waves of optical waveguides // Comput. Math. Math. Phys. - 2005. -V. 45, No 12. - P. 2119-2134.
21. Kress R. Linear Integral Equations. - N. Y.: Springer, 1999. - XIV, 367 p.
22. Solov'ëv S.I. Preconditioned iterative methods for a class of nonlinear eigenvalue problems // Linear Algebra Its Appl. - 2006. - V. 415, No 1. - P. 210-229. - doi: 10.1016/j.laa.2005.03.034.
23. Dautov R.Z., Lyashko A.D., Solov'ev S.I. The bisection method for symmetric eigenvalue problems with a parameter entering nonlinearly // Russ. J. Numer. Anal. Math. Modelling. - 1994. - V. 9, No 5. - P. 417-427. - doi: 10.1515/rnam.1994.9.5.417.
24. D'yakonov E.G. Optimization in Solving Elliptic Problems. - CRC Press, 1996. - 590 p. -doi: 10.1201/9781351075213.
Поступила в редакцию 01.02.2022
Даутов Рафаил Замилович, доктор физико-математических наук, ведущий научный сотрудник НИЛ «Интеллектуальные биомиметические и природосообразные системы» Института математики и механики им. Н.И. Лобачевского Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: rafail.dautov@gmail.com
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online) UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)
2022, vol. 164, no. 1, pp. 68-84
ORIGINAL ARTICLE
doi: 10.26907/2541-7746.2022.1.68-84
An Efficient Numerical Method for Determining Trapped Modes in Acoustic Waveguides
R.Z. Dautov
Kazan Federal University, Kazan, 420008 Russia E-mail: rafail.dautov@gmail.com
Received February 1, 2022 Abstract
An efficient numerical method for determining all trapped modes of the Helmholtz equation based on the finite element method and exact nonlocal boundary conditions is proposed. An infinite two-dimensional channel with parallel walls at infinity, which may contain obstacles
of arbitrary shape, is considered. It is assumed that the frequencies of the trapped modes lie below a certain threshold value. The discrete problem is an algebraic eigenvalue problem for symmetric positive definite sparse matrices, one of which depends nonlinearly on the spectral parameter. A fast iterative method for solving such problems is introduced. The results of the numerical calculations are presented.
Keywords: acoustic waveguide, trapped mode, discrete and continuous spectrum, finite element method, nonlinear spectral problem
Acknowledgments. This study was supported by the Kazan Federal University Strategic Academic Leadership Program ("PRI0RITY-2030").
Figure Captions
Fig. 1. Problem 4.1. Rough triangulation of the O region and the graph displaying the trapped mode on a solid wall (upper figures). Dependence of the error N ^ \X — XNm\ (m = 10) and m ^ \X — XNm\ (N = 1492) (lower figures).
Fig. 2. Problem 4.2. The plots of eigenfunctions on a solid wall.
Fig. 3. Problem 4.3. Rough triangulation of the O region, eigenvalues and functions.
References
1. Callan M.A., Linton C.M., Evans D.V. Trapped modes in two-dimensional waveguides. J. Fluid Mech., 1991, vol. 229, pp. 51-64. doi: 10.1017/S0022112091002938.
2. Evans D.V., Linton C.M., Ursell F. Trapped mode frequencies embedded in the continuous spectrum. Q. J. Mech. Appl. Math., 1993, vol. 46, no. 2, pp. 253-274. doi: 10.1093/qjmam/46.2.253.
3. Evans D.V., Porter R. Trapping and near-trapping by arrays of cylinders in waves. J. Eng. Math., 1999, vol. 35, pp. 149-179. doi: 10.1023/A:1004358725444.
4. Exner P., Seba P. Bound states in curved quantum waveguides. J. Math. Phys., 1989, vol. 30, no. 11, pp. 2574-2580. doi: 10.1063/1.528538.
5. Postnova J., Craster R.V. Trapped modes in elastic plates, ocean and quantum waveguides. Wave Motion, 2008, vol. 45, no. 4, pp. 565-579. doi: 10.1016/j.wavemoti.2007.11.002.
6. Caspers F., Scholz T. Measurement of trapped modes in perforated waveguides. Part. Accel., 1989, vol. 51, pp. 251-262.
7. Evans D.V., Levitin M., Vassiliev D. Existence theorems for trapped modes. J. Fluid Mech., 1994, vol. 261, pp. 21-31. doi: 10.1017/S0022112094000236.
8. Linton C.M., McIver M., McIver P., Ratcliffe K., Zhang J. Trapped modes for off-centre structures in guides. Wave Motion, 2002, vol. 36, no. 1, pp. 67-85. doi: 10.1016/S0165-2125(02)00006-9.
9. Linton C.M., McIver P. Embedded trapped modes in water waves and acoustics. Wave Motion, 2007, vol. 45, nos. 1-2, pp. 16-29. doi: 10.1016/j.wavemoti.2007.04.009.
10. Nazarov S.A. Properties of spectra of boundary value problems in cylindrical and quasi-cylindrical domains. In: Maz'ya V. (Ed.) Sobolev Spaces in Mathematics II. International Mathematical Series. Vol. 9. New York, Springer, 2009, pp. 261-309. doi: 10.1007/978-0-387-85650-6_12.
11. Nazarov S.A. Variational and asymptotic methods for finding eigenvalues below the continuous spectrum threshold. Sib. Math. J., 2010, vol. 51, no. 5, pp. 866-878. doi: 10.1007/s11202-010-0087-3.
12. Evans D.V., Linton C.M. Trapped modes in open channels. J. Fluid Mech., 1991, vol. 225, pp. 153-175. doi: 10.1017/S0022112091002008.
13. McIver M., Linton C.M., McIver P., Zhang J., Porter R. Embedded trapped modes for obstacles in two-dimensional waveguides. Q. J. Mech. Appl. Math., 2001, vol. 54, no. 2, pp. 273-293. doi: 10.1093/qjmam/54.2.273.
14. Sargent C.V., Mestel A.J. Trapped modes of the Helmholtz equation in infinite waveguides with wall indentations and circular obstacles. IMA J. Appl. Math., 2019, vol. 84, no. 2, pp. 312-344. doi: 10.1093/imamat/hxy060.
15. Levitin M., Marletta M. A simple method of calculating eigenvalues and resonances in domains with infinite regular ends. Proc. R. Soc. Edinburgh, Sect. A: Math., 2008, vol. 138, no. 5, pp. 1043-1065. doi: 10.1017/S0308210506001144.
16. Keller J.B., Givoli D. Exact non-reflecting boundary conditions. J. Comput. Phys., 1989, vol. 82, no. 1, pp. 172-192. doi: 10.1016/0021-9991(89)90041-7.
17. Givoli D. Non-reflecting boundary conditions. J. Comput. Phys., 1991, vol. 94, no. 1, pp. 1-29. doi: 10.1016/0021-9991(91)90135-8.
18. Dautov R.Z., Karchevskii E.M. On a spectral problem of the theory of dielectric waveguides. Comput. Math. Math. Phys., 1999, vol. 39, no. 8, pp. 1293-1299.
19. Dautov R.Z., Karchevskii E.M. Existence and properties of solutions to the spectral problem of the dielectric waveguide theory. Comput. Math. Math. Phys., 2000, vol. 40, no. 8, pp. 1200-1213.
20. Dautov R.Z., Karchevskii E.M., Kornilov G.P. A numerical method for finding dispersion curves and guided waves of optical waveguides. Comput. Math. Math. Phys., 2005, vol. 45, no. 12, pp. 2119-2134.
21. Kress R. Linear Integral Equations. New York, Springer, 1999. XIV, 367 p.
22. Solov'ev S.I. Preconditioned iterative methods for a class of nonlinear eigenvalue problems. Linear Algebra Its Appl., 2006, vol. 415, no. 1, pp. 210-229. doi: 10.1016/j.laa.2005.03.034.
23. Dautov R.Z., Lyashko A.D., Solov'ev S.I. The bisection method for symmetric eigenvalue problems with a parameter entering nonlinearly. Russ. J. Numer. Anal. Math. Modell., 1994, vol. 9, no. 5, pp. 417-427. doi: 10.1515/rnam.1994.9.5.417.
24. Dyakonov E.G. Optimization in Solving Elliptic Problems. CRC Press, 1996. 590 p. doi: 10.1201/9781351075213.
I Для цитирования: Даутов Р.З. Эффективный численный метод для определения I захваченных мод акустических волноводов // Учен. зап. Казан. ун-та. Сер. Физ.-ма-\ тем. науки. - 2022. - Т. 164, кн. 1. - С. 68-84. - doi: 10.26907/2541-7746.2022.1.68-84.
For citation: Dautov R.Z. An efficient numerical method for determining trapped modes / in acoustic waveguides. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matema-\ ticheskie Nauki, 2022, vol. 164, no. 1, pp. 68-84. doi: 10.26907/2541-7746.2022.1.68-84.
(In Russian)
