CC BY
26
Вычислительные технологии
Том 13, Специальный выпуск 4, 2008
Численное решение интегральных уравнений трехмерных задач акустики*
А. А. Каширин Вычислительный центр ДВО РАН, Хабаровск, Россия e-mail: elomer@mail.ru
Three-dimensional acoustic problems are considered. The problems are formulated in terms of weakly singular integral equations of the first and second kinds. A discretization of these equations is proposed with the help of the special smoothing method for the integral kernels. Approximate solutions of these equations are found using the variational iterative methods. The results of numerical experiments are presented.
Введение
Математическое моделирование процессов распространения стационарных акустических волн в средах с трехмерными включениями играет важную роль в различных областях науки и техники и приводит к постановке задач математической физики, которые принято называть задачами дифракции или задачами рассеяния. Точные аналитические решения таких задач могут быть построены только в исключительных случаях, поэтому основным направлением их исследования является прямое компьютерное моделирование.
Использование компьютера предполагает предварительное построение дискретного аналога рассматриваемой задачи, которое может осуществляться различными способами. При выборе способа дискретизации необходимо учитывать следующие свойства задач дифракции: решения задач ищутся в неограниченных областях, зависят от трех пространственных переменных, медленно убывают на бесконечности, могут быть быстро осциллирующими функциями и должны удовлетворять условиям излучения на бесконечности. Перечисленные особенности приводят к тому, что алгоритмы численного решения задач дифракции, основанные на дифференциальных формулировках, приводящих к конечно-разностным или проекционно-ееточным схемам, являются малоэффективными. Этот факт стимулировал развитие интегральных методов решения задач дифракции. Исследованиями в этой области занимались U.M. Гюнтер, В.Д. Купрадзе, С.Г. Михлин, А.Г. Свешников, В.И. Дмитриев, С.М. Белоносов, В.А. Цецохо, В.В. Воронин, С. Muller, D. Colton, Е. Kress, J.-С. Nedelec, W, McLean и др.
Применяемые в этих работах методы позволяют сводить задачи дифракции к различным системам из двух интегральных уравнений с двумя неизвестными функциями
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-96024) и Дальневосточного отделения РАН (гранты № 06-1-П14-053, № 06-II-CO-01-001, № 06-III-B-01-018).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
(плотностями). Более перспективным с вычислительной точки зрения представляется развивающийся с конца 1980-х — начала 1990-х годов подход, когда такие задачи формулировались в виде одного слабо сингулярного интегрального уравнения Фредгольма первого или второго рода с одной неизвестной плотностью (С,И, Смагин, H.K. Kleinman, P.A. Martin, Р, Ola), При этом в рамках указанного подхода каждая задача допускает различные эквивалентные формулировки,
В данной работе для численного решения скалярных задач дифракции использовались интегральные уравнения первого рода. Для численного решения полученных уравнений применялся согласованный с шагом сетки специальный метод осреднения интегральных операторов со слабыми особенностями в ядрах, который был разработан и использовался ранее для решения интегральных уравнений краевых задач акустики. Решения диекретизованных уравнений находились итерационными методами вариационного типа. Проверка правильности и точности численного метода проводилась с применением трехмерных задач акустики, имеющих известные аналитические решения.
1. Интегральные уравнения скалярной задачи дифракции
Задача 1 (обобщенная постановка скалярной задачи дифракции). Найти в ограниченной области ^ трехмерного евклидова, пространства Я3 и в неограниченной области Пе = Я3\0г, разделенных замкнутой поверхностью Г класса, Гёльдера, Сг+Ц ((г + в) > 1), комплекснозначные функции ^¿(е) € Н ^^¿(е)), удовлетворяющие интегральным тождествам,
J ^¿(е)У^(е) <1х - А^) J Рг(е)Фг(е) ¿Х = 0 ^¿(е) € Н (П*(е)) , (1)
^г(е) ^г(е)
условиям сопряжения на границе раздела, сред из Пг м Пе
- = <^>г ^ € Н-1/2(Г), (2)
{'Ч/РгМ^г - Ре^е)г = <^,Ре^1>Г У'П € Н 1/2(Г), а, также условию излучения на, бесконечности, для, ре
д^е/д |х| — гкеРе = О (|х|-1) , |х| ^ то, (3)
если на Г заданы функции р0 € Н1/2 (Г) м € Н-1/2(Г).
Здесь и далее <•, ^)г — отношение двойственности на Н 1/2(Г) х Н-1/2(Г), обобщающее скалярное произведение в Н0(Г):
</,д)г = / V/, д € Н0(Г),
г
д — комплексно-сопряженная к д функция; N = д/дп, п(х) — вектор единичной внеш-
ГХ
Pi(e) = (pi(e)U (W + Íji(e))) 1 , k\e) = W + ^(e)) /ci(e), Im(ki(e)) > 0,
ш — круговая частота колебаний, сг(е), рг(е), 7г(е) _ скорости звука, плотности и коэффициенты поглощения сред, заполняющих области Пг(е), сг(е) > 0 рг(е) > 0 7г(е) > О, В работе [1] доказана
Теорема 1. Задача 1 имеет не более одного решения. Введем обозначения
(Аг(е)?) (х) = (^(е)^ О^)г , (Вг(е)?) (х) = (Я^е)^ О^)г ,
(Вг'(е)^) (х) = (^^О^г^Сх 0)г , &г(е) ^ = ехр (¿^(е) |х - /(4п |х - у|)
и рассмотрим потенциалы
<£е(х) = (Ае?) (х), X е Пе, (4)
<£г(х) = Рег (А» (Ж^е + Ы + ) (х) - (Вг' (<£е + Ы + ) (х), X € П»,
где ? е Н-1/2(Г) — неизвестная плотность; рег = ре/рг; знаком "+" отмечаются предельные значения соответствующих выражений на Г когда х — Г го области Пе,
Ядрами потенциалов (4) являются фундаментальные решения уравнений Гельм-гольца и их нормальные производные, поэтому они удовлетворяют тождествам (1) и условию излучения на бесконечности (3), Кроме того, выполнение для них первого из условий сопряжения (2) автоматически влечет за собой выполнение второго условия сопряжения. Подставляя потенциалы (4) в первое условие сопряжения, получим эквивалентное задаче дифракции слабо сингулярное интегральное уравнение Фредгольма первого рода:
<(С?) ,^>г = <^>г ^ е Н-1/2(Г), (5)
где
(С?) (X) = ((0.5 (Ае + РегАг) + (Вг'Ае - Рег^Ве)) (х); ^о(х) = -0.5^о(х) + Рег (А^) (х) - (В»^о) (х).
Справедлива [1]
Теорема 2. Пусть е Н 1/2(Г), е Н-1/2(Г), 7е > 0 или ш не является собственной частотой задачи
Д^ + = 0, х е Пг, = О, х е Г. (6)
Тогда уравнение (5) корректно разрешимо в классе плотностей ? е Н-1/2(Г) м формулы, (4) дают решение задачи 1.
Задача 1 допускает еще одну эквивалентную формулировку в виде интегрального уравнения Фредгольма первого рода со слабой особенностью в ядре. Будем искать ее решение в виде
<£г(х) = (Аг?) (х), х е Пг, (7)
^е(х) = (Ае(^1 - Р^М^)-) (х) - ^(^о - <£г)-) (х), х е Пе,
где рге = рг/ре, знаком "-" отмечаются предельные значения соответствующих выражений на Г, когд а х — Г из области Пг, В этом случае задача 1 сводится к уравнению
<(Аг) ,^>г = <^о,^>г ^ е Н-1/2(Г), (8)
(Д?) (х) = ((0.5 (Аг + РгеАе) + (ргеАеВг - Ве^)) (х)
и имеет место [1]
Теорема 3. Пусть € Н 1/2(Г), ^ € Н-1/2(Г), 7е > 0 или ш не является собственной частотой задачи (6). Тогда уравнение (8) корректно разрешимо в классе плотностей д € Н-1/2(Г) и формулы (7) дают решение задачи 1.
В тех случаях, когда необходимо рассчитать волновое поле в области Qe, задачу 1 предпочтительно сводить к уравнению (5), которое допускает расчет отраженного поля по более простой формуле. По аналогичной причине в случаях, когда необходимо рассчитать проходящее волновое поле в области предпочтительно использовать уравнение (8),
2. Метод численного решения
Излагаемый метод решения интегральных уравнений со слабыми особенностями в ядрах впервые предложен в работе [2]. Его идея заключается в том, что неизвестная плотность находится в виде линейной комбинации гладких финитных функций, образующих разбиение единицы на границе включения. Такой подход не требует предварительной триангуляции поверхности и одинаково просто реализуется как на регулярных, так и на нерегулярных сетках. При дискретизации уравнений поверхностные интегралы приближенно заменяются выражениями, содержащими интегралы по Я3, которые затем вычисляются аналитически, что позволяет достаточно просто находить коэффициенты систем линейных алгебраических уравнений (СЛАУ), аппроксимирующих интегральные уравнения.
Следуя работе [3], кратко опишем общую схему реализации метода. Построим покрытие поверхности Г системой (Гт)м=1 окрестностей узловых точек xm £ Г, лежащих внутри сфер радиуса hm с центрам и в xm, и обозначим че рез {^m} подчиненное ему разбиение единицы, В качестве будем использовать функции
^m(x) = ^m(x) (Y, ж(x)l , ^m(x) = {О1 - rm/hm)3' rm — hm'
У k=i / L 0 ' m — hm;
где x £ Г rm = |x — xm|, £ С1(Г) при Г £ r + ß > 1,
Будем предполагать, что для всех m = 1, 2,..., M выполняются неравенства
О < h' < |xm — xn| , m = n, n = 1, 2,..., M, h' < (2n)1/2am < hm < h, h/h' < q0 < oo,
где h, h' — положительные числа, зависящие от M, q0 не зависит от M,
2\ 1/2 / f Vm(V)
Вместо заданной на Г неизвестной функции д будем искать обобщенную функцию д$г, действующую по правилу (д$г,^)д3 = (д, ^)г. Приближать эту функцию будем выражением
м
q (x) ir (x) ~ ^ gn^n(x)
nn
n=1
где qn — неизвестные коэффициенты; фп(x) = (паП) 3/2 exp (—(x — xn)2/^n).
Умножим обе части интегрального уравнения на фт и проинтегрируем по Г, Тогда при достаточно больших М интегральный оператор первого рода аппроксимируем по формуле [4]
» м
/ (Аг(еП)^т(1Г и т = 1, 2,...,М, (9)
г
ут
п=1
а интегральный оператор второго рода — по формуле [3]
м
(ад + (Вг(е)д))^Г/ и фт X] В^, т = 1, 2,...,М. (10)
п=1
г
Здесь
Атп = о.бвгспехр нто ы - ™ ы)/г™, т=п,
вго = фп (4п)-1 (1 - 0"%))2 + 0.5 (^5(2))^-1 , = 0.5кг( ,
' тп |хт хп1 , ^"тп (^т, + ^п) , Ттп гmn/^'mn,
*1 = ^тп - ®Ттп, ^2 = ^тп + ¿7тп, ^ (г) = ехр (-г2) 1 + 2п 1/2г / ехр (¿2) <И
о
Атт = вгТеТ [^г(е)™ + 2 (оттП1/2)-1 + 2п1/2атт (1 - 2 ) V3) ^т^
3
Вг(е) — (4пГтп) Птп ехр [}кг(е)гтп) (^кг(е)гтп ^ фп птп — ^ ^ п1т (х1т ^п^
1=1
Вгтет = -Св(хт) при а = 0.5, Вгтет = -1 - Ов(жт) при а = -0.5, Г , ^ _ у^ п*тпфп _ 1 Г д 1 _1
Интегральные операторы в левых частях уравнений (5) и (8) являются композицией интегральных операторов (9) и (10), аппроксимируем их по формулам
м
(Сд) фт^Г и ^ (А5пВгтп - РегАГпВетп) дп, т =1, 2,..., М, (11)
п=1
м
(Яд) фт^Г (РгеА5пВгтп - АГпВетп) дп, т =1, 2,..., М, (12)
п=1
а правые части уравнений (5) и (8) — по формулам
м
^о фт^Г И фт ( Св(Жт)фо(Хт) + ^ фп [РегА™ф1(Хп) - фо(Хп)Вгпт] ) , т =1, 2, ..., М,
п=1
г
£
J ipоpmdr и ipm^a(xm), m = 1, 2,..., M. г
Решая соответствующие СЛАУ, найдем приближенные значения плотностей в точках дискретизации для каждого из интегральных уравнений. После этого искомые решения задач акустики могут быть одинаково просто и точно вычислены как в дальней, так и в ближней зоне,
3. Результаты численных экспериментов
Правильность и точность численного метода проверялись с использованием тестовых задач, имеющих точные аналитические решения, В качестве последних выбирались сформулированные в виде интегральных уравнений внутренние и внешние краевые задачи Дирихле и Неймана для уравнения Гельмгольца [5] в областях, ограниченных единичной сферой и эллипсоидом с полуосями (0.4,1, 0.2), и задача о рассеянии плоской акустической волны на единичном шаре [6]. Интегральные операторы в этих задачах имеют вид (9)—(12), волновые числа ki(e) выбирались из диапазона от 0 до 21, число точек дискретизации M = 500 ... 16 000,
Для рассматриваемых диекретизованных уравнений исследовалась возможность их приближенного решения итерационными методами вариационного типа. Тестировались метод минимальных невязок (MINRES), модификации метода сопряженных градиентов для невырожденных СЛАУ (CGNR, CGNE), метод биеопряженных градиентов (BiCG), обобщенный метод минимальных невязок (GMRES) и стабилизированный метод биеопряженных градиентов (BiCGStab) [7].
Численные эксперименты показали, что величина погрешности в рассматриваемых задачах при больших M имеет порядок не хуже h2 ~ M-1, При этом относительная погрешность вычислений в норме сеточного пространства H^(Qi(e)) для всех задач принадлежит интервалу от 10-5 до 10-3 в зависимости от исходных данных. Все итерационные методы, кроме MINRES, позволяют достаточно быстро находить искомые решения СЛАУ, причем число итераций, необходимых для отыскания решений с за-
M
GMRES,
Список литературы
[1] Каширин A.A., Смагин С.И. Обобщенные решения интегральных уравнений скалярной задачи дифракции // Дифференц. уравнения. 2006. Т. 42, № 1. С. 79-90.
[2] Смагин С.И. Численное решение интегрального уравнения I рода со слабой особенностью на замкнутой поверхности // Докл. АН СССР. 1988. Т. 303, № 5. С. 1048-1051.
[3] Ершов Н.Е., Смагин С.И. Численное решение трехмерной стационарной задачи дифракции акустических волн на трехмерном упругом включении. Владивосток, 1989 (Препр. АН СССР, ДВО, ВЦ).
[4] Каширин A.A. Исследование и численное решение интегральных уравнений трехмерных стационарных задач дифракции акустических волн: Дис. ... канд. физ.-мат. наук. Хабаровск, 2006. 118 с.
[5] Nedelec J.-С. Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems. Berlin: Spinger-Verlag, 2001. 316 p.
[6] Селезов И.Т., Кривонос Ю.Г., Яковлев В.В. Рассеяние волн локальными неодно-родностями в сплошных средах. Киев: Наук, думка, 1985. 136 с.
[71 Saad Y. Iterative Methods for Sparse Linear Systems. Boston: PWS Publ. Co., 2000. 460 p.
Поступила в редакцию 20 февраля 2008 г.
