Численное решение трехмерной стационарной задачи дифракции акустических волн Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 15, № 1, 2010

Численное решение трехмерной стационарной задачи дифракции акустических волн*

Н.Е. Ершов, Л. В. Илларионова, С. И. Смагин Вычислительный центр ДВО РАН, Хабаровск, Россия e-mail: illarionova_l@list.ru

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

Ключевые слова: уравнение Гельмгольца, трехмерная задача дифракции, метод граничных интегральных уравнений.

Введение

Задачи распространения стационарных акустических колебаний в средах, содержащих трехмерные включения, имеют большое значение для многих областей науки и практики. Они встречаются, например, в дефектоскопии, в акустике океана и атмосферы, в геофизике. Для решения этих задач применялись конечно-разностные и проекционно-сеточные методы численного решения [1-5], а также численные алгоритмы, основанные на преобразовании исходных дифференциальных задач к эквивалентным им интегральным уравнениям при помощи прямого [6, 7] и непрямого [8-10] вариантов метода интегральных уравнений.

Данная работа является продолжением публикаций [11 -14], посвященных развитию непрямого метода интегральных уравнений для численного решения задач дифракции. С помощью представления решения в виде потенциалов простого слоя исходная задача сводится к системе двух граничных интегральных уравнений Фредгольма первого и второго рода со слабыми особенностями в ядрах на поверхности включения. Впервые такой подход был предложен В.А. Цецохо [8] и использовался для численного решения осесимметричных и плоских задач дифракции электромагнитных и упругих колебаний

*Работа выполнена при финансовой поддержке ДВО РАН (грант № 09-П-СО-005), РФФИ (гранты № 08-01-00947 и 09-07-98509) и Программы Президиума РАН № 2. © ИВТ СО РАН, 2010.

в работах [8, 9, 15]. Применение данного подхода позволяет снизить размерность задачи и локализовать поиск неизвестных функций на границе включения, а затем рассчитать искомые потенциалы волновых полей в любых точках пространства как внутри, так и вне включения. При этом ядра полученной системы интегральных уравнений имеют более простой вид, чем в интегральных уравнениях прямого варианта. Корректная разрешимость полученной системы, а также ее эквивалентность исходной задаче дифракции в классической и обобщенной постановках исследованы в работах [10, 12].

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

В настоящее время при численном решении граничных интегральных уравнений наибольшее распространение получил метод граничных элементов, к преимуществу которого следует отнести его теоретическую прозрачность [7, 16]. Коэффициенты получающихся при этом алгебраических систем выражаются в виде многомерных интегралов от сложных выражений с особенностями при совпадении аргументов, вычисление которых требует значительных ресурсов ЭВМ. Применение других, более экономичных методов, хорошо зарекомендовавших себя при решении одномерных уравнений, к которым прежде всего можно отнести методы коллокации и квадратур, сдерживается серьезными трудностями их теоретического обоснования, особенно для уравнений первого рода [8, 9]. В данной работе применяется численный метод, сочетающий в себе простоту реализации метода коллокации с возможностью полного теоретического обоснования. Этот метод был предложен и впервые использовался в упрощенной форме для решения граничных уравнений первого рода в работах [14, 17].

Аппроксимация интегральных уравнений системой линейных алгебраических уравнений осуществляется при помощи разбиения единицы на границе включения, связанного с системой узловых точек, а также согласованного с порядком дискретизации метода осреднения слабо сингулярных ядер интегральных операторов. В его основе лежит возможность приближения поверхностных потенциалов выражениями для объемных потенциалов, плотности которых локализованы вблизи граничной поверхности. Возникающие при дискретизации многократные интегралы вычисляются аналитически. Это позволяет получать в явном виде формулы для аппроксимации граничных интегральных операторов с особенностями в ядрах и использовать их для вычисления коэффициентов систем линейных алгебраических уравнений. При этом не требуется предварительная триангуляция поверхности, а достижение необходимой точности, как показали результаты численных экспериментов, происходит при сравнительно небольшом количестве точек дискретизации [14-18].

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

1. Постановка задачи

Рассмотрим стационарную задачу дифракции акустических волн в однородной безграничной среде с включением (см, [7, 19])

ДФг(е) + k(e)$i(e) =0 В 0i(e), (1)

Л дФг дФе дФо Pi§i - реФе = реФ0, - = -7^- на S, (2)

j — гкеФе = о(|ж|-1) при |ж| —> сю. (3)

Здесь Hi — ограниченная область из R3 со связной границей S G C(в > 0e = К3\Г2г, п = п{х) — нормаль к S, направленная вне Qj, в точке х = (х1, х2, х3). Константы Yi(ey, Pi(e), Ci(e) характеризуют коэффициент поглощения, плотность и скорость звука сред, заполняющих 0^), Yi(e) > 0 Pi(e),ci(e) > 0 ш ~ круговая частота колебаний, Ф0 — заданная комплексная амплитуда потенциала смещений акустических волн в 0e, Фi, Фе — комплексные амплитуды потенциалов смещений проходящих и дифрагированных акустических волн в 0e, ki(e) — волновые чпела 0^), k2(e) = + ¿Yi(e)) /c2(e) •

Классическим решением задачи дифракции (1)-(3) называют функции Фi G C2(0i)n C1(Qi), Фе G C2(Qe) П Cl(Qe U S), удовлетворяющие уравнениям Гельмгольца (1), условиям сопряжения (2) и условию Зоммерфельда (3).

Решение задачи (1)-(3) будем искать в виде потенциалов простого слоя

Ф^е)(ж) = (Ai(e))(x) = j Gi(e)(x - y)qi(e)(y) dSy, X G 0i(e). (4)

S

Здесь qi(e) — неизвестные плотности вспомогательныx источников, Gi(e) — фундаментальные решения уравнения (1),

exp(ifci(e)r) _

^i(e){X) - ^ , г - \х\.

Функции, определенные формулами (4), удовлетворяют уравнениям (1) и условию излучения на бесконечности (3). Подставив их в условия сопряжения (2) и воспользовавшись формулами для скачка нормальных производных потенциала простого слоя [6], получаем систему интегральных уравнений Фредгольма первого и второго рода со слабыми особенностями в ядрах:

PiAiqi - peAeqe = реФо,

Уг + Яе R R _ <9Ф0 ня q ^

—---Ь пiQi — Beqe — —— на Ь, (о)

2 дп

где

д

(Bi{e)qi(e))(x) = / — Gi(e)(a; - y)qiiye){y)dSy, х G S.

Sx

Теорема 1 [7, 19]. Задача дифракции (1)-(3) имеет не более одного классического решения.

Теорема 2 [20]. Задача дифракции (1)-(3) имеет не более одного обобщенного решения.

Теорема 3 [10, 21], Пусть Б € С1'13 (I + в > 1); к2 не является собственной частотой задачи

ДФ + к2еФ = 0вП,, Ф = 0 на Б.

Тогда система интегральных уравнений (5) корректно разрешима при любой правой части, Ф0 € С1'13(Б), дФ0/дп € С1-1'в(Б) в пространстве пар плотностей qi,qe € С1-1,13(Б), при этом, формулы (4) дают решение задачи, дифракции (1)-(3).

Вопрос о существовании и единственности обобщенного решения интегральных уравнений задачи дифракции (1)-(3) рассмотрен в [12],

2. Аппроксимация интегральных уравнений

Методику аппроксимации задачи (5) системой линейных алгебраических уравнений рассмотрим на примере следующего уравнения:

a(x)q(x) + J К(x,y)q(y)dБy = /(х), х € Б, (6)

где а, К, / — заданные функции; ^ ^ ^^^^^^^^^^^ плотность; ядро К можно представить в виде суммы К = В+Т; Т — гладкая на Б функция; В — функция с особенностью, которая может иметь вид

1*-»Г'. Р= 1,2,3.

Пусть {Б,}^ — покрытие Б гада Б, = {х € Б : |х — х,| < к,}, где Xi — точки на Б, а {у,}^ — множество функций, образующих разбиение единицы на Б, подчиненное этому покрытию, т, е,

N

вирр у, С Б,, 0 < у,(х) < 1, = 1 V х € Б.

i= 1

В качестве у, будем использовать удобные для вычисления функции

'Ш = , у*(х) = ( I1 " " ж ^' (7)

Е I0'

Приближенное решение уравнения (6) будем искать на сетке {х,}^1;

Хг = = [ Х(Рг(х)<ЛБх, <Рг= [ У^х)<ЛБх, 1= 1,2,. . . (8)

Я Я

узлами которой являются "центры тяжести" функции у,. При этом будем предполагать, что для всех г = 1, 2,..., N выполняются неравенства

0 <к* < |х, — х, | , г = з, к* < к, < к, к/к, < ро < то, з = 1, 2,...,^ (9) где к*, к — положительные числа, зависящие от N р0 не зависит от N.

Вместо неизвестной функции ^^ ^^^адной на Б, будем искать обобщенную функцию q5s, действующую по правилу

, п)«з = п)я, V п € На (К3) , а > 1/2,

где (•, — отношение двойственности на Н-а+1/2(Б) х На-1/2(Б), обобщающее ека-

Н0(Б)

странств На(Б) см., например, в [22]),

= У ^Б, Vq,п € Н0(Б).

Я

Обобщенную функцию q5s будем приближать выражением

N

И (10)

,=1

где qi — подлежащие определению компоненты неизвестного вектора ф = ..., qN)

ф,(х) = (па2)-3/2 ехр (— (х — х)2/о2) х € К3,

а, — радиусы осреднения, которые выбираются из условий хорошей аппроксимации соответствующих интегральных уравнений [17],

/2>г [ <Рг(у)

—г 3

ИЛИ

Я

Аппроксимация обобщенной функции q$s, заданной па всем прострапстве К3, в виде (10) представляет собой ее замену линейной комбинацией бесконечно дифференцируемых в К3 функций 'ф, (х), близких к 0 вне достаточно малой окрестности поверхноети Б. Таким образом, если

= {я, гПг= Г]{х)ф1{х)(1х,

ТО

И

J п(х)у,(х^Б ^ п,

Я

N \ N „ N

^ Г] \ / ф)ф^х)йх = ^(д, Уг)яЩ -

л=1 / К3 ,= 1 К3 ,= 1

q, ПгУА ~(q,n)s Vq € Н-а+1/2(Б),п € На (К3)

,= 1 / Я

В частности, можно показать, что для достаточно гладких функций п

= У г](х)ф^х)с1х = г](х^ + О (Л2) = ^ J ц(х)(р^х)с18 + О (Л2) .

К3 ' 5

Следовательно, для любых таких функций п и д € Н 1(Б) справедливы равенства

N

N

N

N

^ n)s = (ni ^ vdsqi — ni),^i)s((q— %)(n— n^^^s--

i=l

S i=1 N

i=l

i=1

= +О (И2) / 7]{х)фг{х)йх + О (И2) = ( + 0{к2).

г=1 г=1 ¿3 \ г=1 / Д3

Аппроксимация плотности потенциала простого слоя объемной плотностью вида (10) позволяет получить весьма простые формулы для вычисления главных частей интегралов с особенностями в ядрах и теоретически обосновать применяемый численный метод с "саморегуляризацией" для решения слабосингулярных интегральных уравнений первого рода [17].

Уравнение (6) с учетом вышеизложенного будем аппроксимировать системой линейных алгебраических уравнений

N

«¿Ф + X! + = /¿, г = 1, 2,..., N.

3 = 1

Здесь

а = а^йБ, = / / В(х,у)фг(х)фз(у)йуйх, Тз = Т(х^,хз), ^ = / ¡^¿Б.

Интегралы Bij вычисляются аналитически с использованием формул

д дхР

lim

ht^ о

|x - У1 f ф%{х)ф,{у)

I \х-У\

г.

ij

dydx =---—г^-erf (7у)(2 - /л^),

(rij)

(Н)

/ir>i (x)

—zz—g(x) dS

S

0 V i = j, p = 1, 2, 3,

полученных в [13] для вычисления пределов при rij — 0, Здесь erf (x) — интеграл ошибок [23], g(x) — любая непрерывная в окрестности узловой точки xi и интегрируемая в R3 функция,

2 , 2ч 1/2 _ т _ ехр(-72) _ _ _

<Jij — \CTi ~Г CTj ) , ^ij — J- 7T1/2erf (7 ••) ' ^ ~~ r4/<7ij> r4 ~ \x% Xi\-

Интегралы (11) для i = j находятся предельным переходом при xj — xi в выражениях для правых частей соответствующих равенств

lim

Xj >xt г.

erf (Yij) 2

lim

x

— xPp

ij

n1/2Oii Xj ^Xi (rij)

-erf (7y)(2 — ¡iij) = 0, p = 1,2,3. (12)

3

3. Численное решение задачи

При аппроксимации системы интегральных уравнений ядра интегральных операторов представляются в виде суммы двух слагаемых, одно их которых содержит особенность при совпадении аргументов, а другое гладкое. Интегралы, содержащие слагаемые с особенностями в ядрах, вычисляются аналитически по формулам (11), (12), а интегралы без особенностей — по квадратурным формулам прямоугольников, В результате получаем систему линейных алгебраических уравнений для отыскания приближенного решения д^, д^ системы (5) следующего вида:

У^ (fijiPiAikjgij - peAekjgej) = реФ0(хк), j=i

N ^ , s

(gik + дек)/2 + Ihkj'lij - Bekjgej) = —!k= 1,2,...

(13)

Qn , .. N.

j=i

Здесь

a - Air. ^ _ erf Ы) , exp(ikrkj) - 1

Акз АЫ + 4vr rk3 '

Bkj = = [erf + exp(ikrkj) • (1 - ikrkj)], 1 3

Ukj = - (Xk ~ rf) Пк> Пк = П(Хк) ПРИ k i"

rkj p=l

Akk = °kk (1 + ik/okk )/4n, Bkk = 0.

Для повышения точности вычисления коэффициентов системы (13) при j = k, которые по модулю больше остальных, и нахождения приближенного решения системы (5) с меньшей погрешностью применяются тождества, полученные с помощью теоремы Гаусса [6]:

— ± i ^-P-dSy = ^ ^ ^ ± [ (Я(У)#" + dSy. (14)

2 J dnx 47гг y 2 J \4yyJdnx чу J dny J 4тгг y K '

S S

Улучшения здесь связаны с тем, что тождества позволяют выделить главную часть соответствующего интеграла и избавиться от особенности в подынтегральном выражении при совпадении аргументов. При этом последние интегралы в правой части равенства (14) аппроксимируются по следующим правилам:

[a(v)J-±dS~ v f--dS n.-lr rj k-TN

s j=l,j=k S

Приближенное решение задачи (l)-(3) в любых точках областей П^ и Пе находим, заменяя интегралы (4) выражениями

Фг(е}(Х) = ^ j=l

N Г 1 — erf (rj)

Cj(e) (x, Xj) —

qi(e)j<Pj, X e Qj(e).

4. Аппроксимация интегральных уравнений на "звездной" поверхности

Существует достаточно широкий класс замкнутых поверхностей Б, допускающих взаимно однозначное гладкое отображение на трехосный эллипсоид Г, Это позволяет свести задачу решения системы (5) на такой поверхности к решению системы на поверхно-Г

этом могут выбираться исходя из возможности использования на нем почти равномерных сеток,

ГГ

эллипсоид, который в декартовой системе координат Ох1х2х3 задается уравнением

3 /.п\ 2

т. íУ-

i \ap p=i ч '

Зададим сетку Xj £ Г и локальные координаты г/j = (у], у2, у3) (j = 1, N) так, чтобы в

окрестности узловой точки Xj £ Г уравнение поверхности Г имело в ид y3 = fj (yj,y2).

Сетку на эллипсоиде Г выберем, равномерно располагая узловые точки на эллип-

X2

ГГ

вой точки Xj = (x],x2,x3) осуществим проектированием ее на касательную плоскость, проходящую через Xj. Пусть nj = (n],n2,n3) — нормаль к этой плоскости. Выберем направляющий вектор Ip (p £ {1, 2, 3}) системы Ox1x2x3 из условия

|n7-IJ = min |n7-Im|.

J p me[l,3] j

Орты ej1; ej2j ej3 локальной системы координат Oyiy2y3 определим выражениями

Ip Tljijljlp) I Ip - Hjirijlp) I

eji = 77-,, г M, ei2 = rij x eji, ej3 = щ.

Рассмотрим теперь вопрос об аппроксимации интегральных уравнений на "звездной" поверхности S £ C 1'в (в > 0), Такая поверхность состоит из точек, которые в декартовых координатах определяются формулами

(xl,x2,x3) = R(ip,e)(a\ cos y sin в, a2 sin y sin в, a3 cos в), 0 < y < 2п, 0 < в < п,

(15)

где R £ C([0, 2п] x [0,п]) — некоторая функция, удовлетворяющая условиям:

а) R(y,9) > 0 для всех y £ [0, 2п], в £ [0,п],

б) = lim И. И.„ = 0, р. 0,1. *=1,2.

дтр Т1^2п дтр Г2^0,п дтр

Обозначим через Г эллипсоид с центром в точке x = 0 и полуосями ai, a2, a3, т.е. в декартовых координатах точки Г определяются формулами

(xl,x2,x3) = (ai cos y sin в, a2 sin y sin 9, a3 cos 9) , 0 < y < 2п, 0 < в < п.

Заданные на поверхности S интегральные операторы системы (5), используя соотношение (13), преобразуем следующим образом [13, 14]:

У К(x,y)q(y) dSy = j К(x,y)q(y)J(yy) dTy = ^ J К(x,yJ)<j{y)q{y)J(yy) dTy =

S г j=1 г

N Г -

= £ / К'(Х>,УЩ(X)X(у')|gj(X)| dy)dy2, (16)

j=1 D

где К(x, y) — ядра интегралов системы (5); q — неизвестные функции, определенные в точках на поверхности S; К К', у X, j <j — выражения для К, q и < в координатах поверхности Г и локальных координатах yj = (у],у2, y|); J — модуль якобиана преобразования координат точек поверхности S к координатам точек поверхности Г (у),у2) _ локальные координаты поверхности Г в окрестности j y p = y p(yj) (p = 1, 2, 3), yj £ Dj Dj — прообраз rj в системе yj; |gj(yj)| — дискриминант метрической формы поверхности в системе yj.

Sx

определим из соотношений

np(x) = Pp^2(Pi)2)-1/2, p = 1, 2, 3,

p i=i

где

dx' dxm dx' dR(d,y)^t . dX

p, = E w'щ, = + ^ W

t,m=1

dx'

T— = atvkt} k = 1,2, i= 1,2,3, дтк

Ti = д, T2 = y, Vi,i = ^2^4, Vi,2 = ^2^3, Vi,3 = -^i, V2,1 = V2,2 = V2,3 = 0,

= sin d, w2 = cos д ш3 = sin y, ш4 = cos y, здесь eptm — символ Леви — Чивита,

5. Вычисление поверхностных интегралов

Интегралы Тр^, у., входящие в коэффициенты системы (13), учитывая соотношения (15), (16), запишем в виде функций точек ¡г "стандартной" поверхности Г:

Ъ = / пШШГу, <р. = I 3 = 1, 2,..., N. (17)

г г

Их значения могут быть найдены численным интегрированием.

При вычислении интегралов у перейдем сначала к параметрическому представлению Г в системе локальных координат ), а затем введем полярные координаты и приведем Тр^ к виду

ъ = /Iф;Щ1 - \%/Ьз\2)3(Юу,. =

О; В,-

г 1 п2п

= (Ъ)2 / (1 - Р2)3Р / ф(р,ф)йфйр. (18)

ио ио

Здесь — выражепия у в локальных координатах, связанных с узловой точкой ф3 _ непрерывная функция, гладкость которой определяется гладкостью поверхности Б и функций, образующих разбиение единицы. Так как интегралы (18) от любого одночлена, содержащего нечетную степень какой-либо координаты, равняются нулю [24], то после замены переменной г = р2 получаем

Уз

2п ^ Мф

о

Оптимальные веса Ск и узлы гк для (19) при различпых М имеются в [25], Положим Мф = 4Ыг. Тогда алгебраические порядки точности квадратурных формул по г и ф будут одинаковыми и равными 4МХ — 1,

Интегралы <р. содержат особенности в ядрах при у —> х^. Их можно преобразовать к виду

* = / = / ~ =

О, 3 3 В,-

1 2п

= ъА (1 — г)3 Ф"(гк3,ф) йфйг. (20)

~з

о

Приближенные значения интегралов (20) получим заменой этих интегралов по г и ф квадратурными формулами, имеющими вид

д. (">

В (21) ик и гк выберем такие, которые обеспечивают максимальный алгебраический порядок точности при интегрировании с весом (1 — г)3. Тогда при Мф = 4М порядок точности формулы (20) будет равным 4МХ — 1 (см, [24, 25]),

6. Результаты численных экспериментов

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

Применение предварительного сведения системы граничных интегральных уравнений к системе уравнений по "стандартной" поверхности в виде трехосного эллипсоида существенно упрощает процедуру построения сеток и делает алгоритм более универсальным, Ранее такой подход использовался для построения сеток при численной реализации алгоритма решения методом коллокации интегральных уравнений пространственной задачи дифракции акустических волн на упругом теле [13, 14],

Для численного решения системы (13) с не очень большим числом узлов дискретизации N (М < 500) применялся метод Гаусса с выбором ведущего элемента в строке, а при N > 500 — обобщенный метод минимальных невязок (GMR.ES) [26],

Введем обозначения: — приближенное решение задачи, полученное в ре-

^ г ЦФГ - Ф*1к2(П;)

з\льтате численных расчетов, Ф¿, Фе — точное решение, дг = ---г- х 100,

||Фг|и2(П;)

ЦФ^ - Фв|1ь2(Пе)

ое =- --—1 х 100 — относительные погрешности решения, || ■ \\Ь2 — сеточная

|ФеУь2(П е)

£2-норма, Пе — шар, содержащий Пг, с радиусом ~ 10Я, где Я — характерный размер области Пг,

Для проверки правильности работы алгоритмов и созданного программного комплекса решались следующие тестовые задачи (примеры 1-3),

Пример 1. Первая внутренняя (внешняя) краевая задача. Найти функцию Фг (Фе), удовлетворяющую в Пг (в Пе) уравнению (1) и краевому условию

Фг(х) = ехр(гкгх3) I Фе(х) = ехр(гкег)/(4пг) ) , х € Б,

Пг

с центром в точке (0, 0, 0) и полуосями (0.75,1, 0.5),

V = 1, 7г(е) = 0, Рг = 1, Ре = 2, Сг(е) = 1. Здесь и далее г = |х|. Точное решение имеет вид

Фг(х) = ехр(гкгх3), х € Пг ^Фе(х) = ехр(гкег)/(4пг), х € П^ . (22)

На рис, 1, о приведены графики зависимости относительных погрешностей 5гъ 5е внутренней и внешней задачи соответственно от числа N точек дискретизации для первой краевой задачи.

а б

Пример 2. Вторая внутренняя (внешняя) краевая задача. Она заключается в отыскании функции Ф^ (Фе), удовлетворяющей в Qj (в Qe) уравнению (1) и краевому условию

дФ (х) I дФ (х) ik r_1 3 Л \

—^-= %кггг?{х) exp(«fc¿x3) I —^-= —т—^— exp(¿fcer) } прхр 1 , х Е S,

dn у dn 4пг3 у

а также условию излучения (3) в случае внешней задачи; Qi — шар с центром в начале координат единичного радиуса,

Ш =1, Yi(e) =0, Pi =1, Pe = 2, Cj(e) = 0, 05.

Точное решение определяется формулами (22), На рис, 1, б" приведены графики зависимости относительных погрешностей ^ и Se внутренней и внешней задачи соответственно от числа N точек дискретизации для второй краевой задачи.

Пример 3. Задача дифракции (1)^(3). Пусть QQ — внутренность эллипсоида S с центром (0, 0, 0) и полуосями (0.75,1, 0.5), Параметры сред ш = 1, рi = 1, pe = 2, Ъ(е) = 0, Положим

Ргф*-РеФ: 5£о = 5Ф:_5Ф: pe on on on

где

Ф*=ехР(^3), Ф: = (о,о,0.3).

4n|x _ y|

Тогда решением задачи (1)-(3) являются функции Фi = Ф*, Фe = Ф*. На рис. 2 показаны графики зависимости относительной погрешности Si (Se) от числа N точек дискретизации для вариантов ci = 0.8, ce = 1; ci = 0.6, ce = 0.5.

Пример 4. Рассматривается задача дифракции (1)^(3). Поверхность вклю-S=S1

где Si состоит из точек x = (x1,x2,x3), определяемых формулами (x1, x2, x3) = Е(в) (ai cos у sin в, a2 sin у sin в, a3 cos в) , в G [0, п], у G [0, 2п),

ct = 0.3, се= 1 Ci — 0.6, Се= 0.5

0.125 -

0.650 -

0.550 -I

Я(6) = 2 - (в/1.1)1Л ^(п - в)/1.8^ , «1 = 0.75, а2 = 1, аз = 0.5, (23) показана па рис, 3, б. Источник акустических волн — плоская волна вида

Фо(х) = ехр(гкех3).

Параметры среды и включения следующие:

Сг = 0.1, Се = 0.18, Рг = 1, Ре = 3, V = 1, ^¿{е) = 0.

Количество точек дискретизации N = 970,

На рис, 4 изображены линии уровня и проективная поверхность функции

П(Ф) = + Фо|-|Фо|

на квадрате |х1 '21 < 4.5, х3 = 0, Рисунки 4, а и б" позволяют судить о местоположении, размерах и форме тела, а также о физических параметрах вмещающей среды и включения. Видно, что область наибольших значений функции п(Ф) находится над телом

Рис. 3. Поверхности 51 (а) и Б2 (б)

и ее размеры вдоль оси x2 больше, чем вдоль оси x1, Точка максимума расположена над центром тела. Плоскости x1 = 0 и x2 = 0 являются плоскостями симметрии задачи.

Пример 5. Отличия от примера 4: S = S2, поверхпость S2 определяется формулами (4), где

Я(в) = -0.5 (eos 29 + (1.1 - sin2 2в)1/2)1/2 + 1.25, ai = 0.75, a2 = 1, аз = 0.5; параметры среды

Ci = 0.4, Ce = 0.27, Pi =1, Pe =1, Yi(e) = 0.

S2

Как и для примера 4, рис, 5, о и б" позволяют определить местоположение, размер и форму тела. Изменение формы тела и значений параметров включения приводит к тому, что здесь линии уровня имеют несколько иной вид.

-4 -2 0 2 4

б

Рис. 5. Линии уровня п(Ф) (а) и проективная поверхность п(Ф) Для примера 5 (б)

Пример 6. Отличается от примера 5 тем, что включение имеет другую ориентацию в пространстве. Углы наклона осей включения, отложенные от положительного направления осей Ох2 и Ох3 против часовой стрелки, равны п/4. Результаты расчетов представлены на рис, 6, Здесь видна симметричность линий уровня относительно оси

Ох2. Область наибольших значений функции п(Ф) в отличие от рис, 4, о и ¿¡"смещена в

Ох2

Пример 7. В отличие от примера 5 падающее поле создается точечным источником

вида

Фо(х)

ехр(гке|х - у|)

|х - у1 :

расположенным в точке у = (0, -2, 2); параметры сред следующие:

Сг = 0.4, Се = 0.2, Рг = 3, Ре = 2, ^г(е) = 0.

Рис. 7. Линии уровня п(Ф) (я) и проективная поверхность п(Ф) Для примера 7 (б)

На рис, 7 изображены линии уровня и проективная поверхность функции п(Ф) на квадрате |xx'21 < 4.5, x3 = 0, В этом случае линии уровня симметричны относительно оси Ox2 и указывают на расположение источника акустических волн. Исчезает локальный максимум в центре тела включения.

Пример 8. Отличается от примера 7 только положением источника, который расположен в точке y = (-2, 2, -2), Результаты расчетов представлены на рис, 8, Видно, что изменение координат источника приводит к изменению симметрии линии уровня (в данном случае они симметричны относительно плоскости Ox1 = -Ox2),

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

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

[1] Самарский A.A. Теория разностных схем. М.: Наука, 1977. 656 с.

[2] Марчук Г.И., Агошков В.И. Введение в проекционно-сеточные методы. М.: Наука, 1981. 416 с.

[3] Завадский В.Ю. Моделирование волновых процессов. М.: Наука, 1991. 248 с.

[4] ihlenburg F. Finite Element Analysis of Acoustic Scattering. Springer, 1998. 239 p.

[5] Braess D. Finite Elemente. Berlin; Heidelberg: Springer-Verlag, 2003. 359 p.

[6] Купрадзе В.Д., Гегелия Т.Г., Башелейшвили М.О., Бурчуладзе Т.В. Трехмерные задачи математической теории упругости и термоупругости. М.: Наука, 1976. 664 с.

[7] Колтон Д., Кресс Р. Методы интегральных уравнений в теории рассеяния. М.: Мир, 1987. 311 с.

[8] Цецохо В.А. Задача об излучении электромагнитных волн в слоистой среде с осевой симметрией // Вычислительные системы. Новосибирск: Изд-во ВЦ СО АН СССР, 1964. Вып. 12. С. 52-78.

[9] Цецохо В.А., Воронин В.В., Смагин С.И. О решении задач дифракции потенциалами простого слоя // Докл. АН СССР. 1988. Т. 302, № 2. С. 323-327.

[10] Смагин С.И. Об одной системе интегральных уравнений теории дифракции // Дифферент уравнения. 1990. Т. 26, № 8. С. 1432-1437.

[11] Смагин С.И., Цецохо В.А. Численное решение задачи дифракции электромагнитных волн на трехмерном включении // Журн. вычисл. математики и мат. физики. 1991. Т. 31, № 5. С. 718-734.

[12] Каширин A.A., Смагин С.И. Обобщенные решения интегральных уравнений скалярной задачи дифракции // Дифференц. уравнения. 2006. Т. 42, № 1. С. 79-90.

[13] Ершов Н.Е., Смагин С.И. Приближенное решение пространственных задач акустики и упругости методом потенциалов // Математические модели, методы и приложения. Хабаровск: Изд-во ХГПУ, 2002. С. 45-115.

[14] Ершов Н.Е., Смагин С.И. Решение пространственных задач акустики и упругости методом потенциалов // Дифференц. уравнения. 1993. Т. 29, № 9. С. 1517-1525.

[15] Дробница В.В., Цецохо В.А. Метод расчета плоского электромагнитного поля в средах со слоем переменной толщины // Математические проблемы геофизики. Новосибирск: Изд-во ВЦ СО АН СССР. 1971. Вып. 2, С. 251-284.

[16] Atkinson К.Е. The Numerical Solution of Integral Equations of the Second Kind. Cambridge, UK: Cambridge Univ. Press, 1997.

[17] Смагин С.И. Численное решение интегрального уравнения I рода со слабой особенностью на замкнутой поверхности // Докл. АН СССР. 1988. Т. 303, № 5. С. 1048-1051.

[18] Белоносов A.C., Цецохо В.А. Вычислительный алгоритм и процедуры сглаживания функций, заданных приближенно в узлах нерегулярной сетки на плоскости // Некорректные задачи математической физики и проблемы интерпретации геофизических наблюдений. Новосибирск: Изд-во ВЦ СО АН СССР, 1976. С. 96-104.

[19] Купрадзе В.Д. Граничные задачи теории колебаний и интегральные уравнения. М.: Л.: 11. III. 1950.

[20] Илларионова Л.В. Задача оптимального управления для стационарных уравнений дифракции акустических волн //Журн. вычисл. математики и мат. физики. 2008. Т. 48, № 2. С. 106-117.

[21] Смагин С.И. Интегральные уравнения задач дифракции. Владивосток: Дальнаука, 1995. 203 с.

[22] Овэн Ж.-П. Приближенное решение эллиптических краевых задач. М.: Мир, 1977. 384 с.

[23] Градштейн И.С., Рыжик U.M. Таблицы интегралов, сумм, рядов и произведений. М.: Гос. изд-во физ.-мат. литературы, 1963. 1100 с.

[24] Крылов В.И. Приближенное вычисление интегралов. М.: Наука, 1967. 324 с.

[25] Крылов В.И., Шульгина Л.Т. Справочная книга по численному интегрированию. М.: Наука, 1966. 395 с.

[26] Saad Y., Schultz М. (¡MRKS: A generalized minimal residual algorithm for solving nonsvmmetric linear systems // SIAM J. Sei. Statist. Comput. 1986. Vol. 7. P. 856-869.

Поступила в редакцию 28 марта 2008 г., в переработанном виде — 28 сентября 2009 г.