CC BY
93
№ 4 (32), 2014
Физико-математические науки. Математика
УДК 517.3
Е. Д. Деревянчук, Е. Ю. Смолькин, А. А. Цупак
МЕТОД ГАЛЕРКИНА РЕШЕНИЯ СКАЛЯРНОЙ ЗАДАЧИ РАССЕЯНИЯ ПРЕПЯТСТВИЕМ СЛОЖНОЙ ФОРМЫ1
Аннотация.
Актуальность и цели. Целью работы является численное исследование скалярной задачи рассеяния плоской волны препятствием сложной формы, состоящим из объемных тел и бесконечно тонких акустически мягких экранов.
Материалы и методы. Задача рассматривается в квазиклассической постановке. Исходная краевая задача для уравнения Гельмгольца в неограниченном пространстве сводится к системе интегральных уравнений по ограниченным многообразиям размерности 2 и 3. Для нахождения численного решения задачи применяется метод Галеркина с использованием финитных кусочно-постоянных базисных функций.
Результаты. Разработан и программно реализован численный метод решения системы интегральных уравнений скалярной задачи дифракции, проведен ряд вычислительных экспериментов.
Выводы. Предложенный численный метод является эффективным способом приближенного решения задач дифракции на препятствиях сложной геометрической формы; он может применяться и для решения более широкого круга задач.
Ключевые слова: скалярная задача дифракции, интегральные уравнения, метод Галеркина, базисные функции, условие апроксимации
E. D. Derevyanchuk, E. Yu. Smol'kin, A. A. Tsupak
THE GALERKIN METHOD FOR SOLVING THE SCALAR PROBLEM OF SCATTERING BY AN OBSTACLE OF COMPLEX SHAPE
Abstract.
Background. The aim of this work is to numerically study the scalar problem of scattering by an obstacle of complex shape consisting of solid bodies and infinitely thin acoustically soft screens.
Material and methods. The problem is considered in the quasiclassical statement; the original boundary value problem for the Helmholtz equation in unbounded space is reduced to a system of integral equations over bounded manifolds of dimension 2 and 3. To find the numerical solution to the problem of diffraction the Galerkin method is applied using finite piecewise constant basis functions.
Results. The projection method for solving the system of integral equations of the scalar diffraction problem was developed and implemented; several computational experiments were performed.
Conclusions. The proposed numerical method is an effective way to find approximate solutions to the scalar problems of diffraction on obstacles of complex shape; it can also be used to solve problems of a wider range.
1 Работа выполнена при поддержке гранта Российского научного фонда, проект № 14-11-00344.
Physical and mathematical sciences. Mathematics
57
Известия высших учебных заведений. Поволжский регион
Key words: scalar diffraction problem, integral equations, Galerkin method, basis functions, approximation condition
Введение
Настоящая работа посвящена методу интегральных уравнений для численного исследования скалярной задачи дифракции плоской волны на системе, состоящей из двух- и трехмерных рассеивателях сложной формы.
Впервые скалярная задача рассеяния плоской гармонически зависящей от времени волны была рассмотрена в [1,2] - в качестве рассеивающей структуры рассматривалась, как и в данной статье, система тел и «мягких» экранов. Краевая задача рассмотрена в квазиклассической. В [2-5] рассмотрена задача, в которой рассеивателями являются дополнительно акустически «жесткие» экраны. В [5] представлено наиболее полное доказательство теоремы единственности решения краевой задачи. Далее краевая задача сводится к системе слабосингулярных интегральных уравнений. В [1-4] доказаны теоремы о гладкости решения этой системы, об эквивалентности интегральных уравнений задаче в исходной дифференциальной постановке. Наконец, установлена обратимость матричного интегрального оператора в пространствах Соболева.
Для нахождения приближенных решений рассматриваемой задачи применяется метод Бубнова - Галеркина с выбором кусочно-постоянных базисных функций как на трехмерном, так и на двухмерном многообразии. Для таких функций доказана выполнимость условия аппроксимации. Кроме того, простота базисных функций позволяет без труда обосновать и программно реализовать численный метод для решения задач рассеяния на объемных телах и тонких экранах сложной формы. В работе приведено описание базисных функций, в том числе на неплоских экранах. Заданных параметрически. Приведены расчетные формулы матричных элементов согласно методу Галеркина; приведены результаты расчетов для различных рассеивающих структур.
1. Постановка краевой задачи и система интегральных уравнений
Пусть Q - связная ориентируемая незамкнутая и ограниченная поверхность класса Cв Ж3 . Край dQ := Q \ Q поверхности есть гладкая
кривая класса Cбез точек самопересечения; трубчатые окрестности края экранов имеют вид
дП;,8 := У BS (x), (1)
xedQ
3
где Bg (x) - открытый шар в Ж .
Пусть Q - ограниченная область, границы которой dQ = Q \ Q -кусочно-гладкая замкнутая ориентируемая поверхность, состоящая из конечного числа поверхностей класса C1. Рассматриваемое тело может быть
неоднородным и описывается общем случае функцией k(x) е C(Q);
3
свободное пространство однородно с волновым числом ke; всюду в Ж выполняются условия Rek(х) > 0 и Imk(x) ^ 0.
58
University proceedings. Volga region
№ 4 (32), 2014
Физико-математические науки. Математика
Предполагаем, что экран Q является акустически мягкий - полное полеудовлетворяет однородному краевому условияю Дирихле на поверхности экрана Q (за исключением точек края экранов).
Требуется определить полное поле и, представимое суммой Uq + us падающей плоской волны, гармонически зависящей от времени, и рассеянного системой препятствий поля.
Задача дифракции сводится к решению системы слабосингулярных интегральных уравнений
u(x)- j(((y)-kke)(x,y)u(y)dy- jG(x,y)q(y)dsy = Uqq(x), xe Q,
Q Q
- j (((У) - ke )( x, y)u(y)dy - jG(x, y)q>(y)dsy = Uqq (x), x ей (2) Q Q
, ik |x— y\
1 e e' 7
Здесь G(x, y) =----------, Uqq и Uq q - сужения падающего поля на
4п | x - y |
Q и Q соответственно; и - неизвестное полное поле в Q , а ф - неизвестная плотность поверхностьного потенциала на Q. Решение системы (2) - пара функций (и,ф).
Краткая операторная форма системы (2) имеет вид
L(u, ф) =
Матричный оператор Le(Q) хЙ-1/2(й) в Le(Q) уй но обратим.
(I - A K Л
K2
S
лл (
V ф У
u0 \q Uq \q
(3)
1/2
L рассматривается как отображение из (Q). В работе [1] доказано, что L непрерыв-
2. Метод Галеркина
Приближенные решения fN =(un, фт) рассматриваемой задачи будем искать в виде
n т
un (x) = Yjc'uvi(x), x е Q; Фт (x) = Vi(x) x eQ,
i=1 i=1
(4)
где Си, сф - неизвестные коэффициенты, vi (x), Vi (x) - базисные функции.
Определим финитные кусочно-постоянные функции vi на объемном теле, считая, что Q - параллелепипед: Q = {x:at <xt <a\,i = 1,2,3}. Разобьем Q на элементарные параллелепипеды
nii3 = {x : x1i < x1 < x1i +b x2,i2 < x2 < x2,i2+2, x3,i3 < x3 < x3,i3+1} ,
a} - a1 a2 - a2 a3 - a3
x1i = a1 +----------1 x2,i2 = a2 +-----------i2, x3,i3 = a3 +--------i3,
1 nj n2 Пз
Physical and mathematical sciences. Mathematics
59
Известия высших учебных заведений. Поволжский регион
где ik =0,..,щ -1, n := «««3; тогда
v; (x) = vi1i2i3(x)
^, x G п / ?2 /'3,
0, x й Пг'/з.
Важно отметить, что при измельчении сетки и на Q необходимо одновременное увеличение значений всех величин. Это требование является необходимым и достаточным для выполнения условия аппроксимации в L2.
Определим базисные функции V; = V/j,;, (ik =0,...,mfc -1) с носителями Ю; на плоском экране Q=P = {x: bj< xj< b,, i = 1,2}. Разобьем P на элементарные прямоугольники:
п/1/2/3 = {x : x1,i1 < x1 < x1,i1 +Ъ x2,i2 < x2 < x2,i2 + 2},
, bj- bb b2- b2.
x1,i1 = b1 +~-1 /1, x2,i2 = a2 +—-2 ^
1 «1 2 «2
где ik =0,...,mk -1, m := m1m2; тогда
Vi (X) = V/j/, (x)
1, X G Я// , 0, x й п/1/2.
Если же экран Q не плоский, то в качестве носителей можно выбирать, например, элементы покрытий {а);-}, i = 1,.,m (m = mjm,) экрана, пересекающихся по кусочно-гладким кривым и таких, что ^)ю; = Q и
i
lim max Mq (ю; ) = 0, (5)
m^» ;
где Mq - мера на экране.
В работе реализован важный случай, когда экран Q может быть задан параметрически: x = x(tj,t,), (tj,12)g P с Ш , где P - компакт (например, прямоугольник). Тогда V; = Vi;i2 - функции-«ступеньки», носители которых
суть образы в Q прямоугольников П; при заданной параметризации. Потребуем, чтобы lim maxM_D (П/j ;2) = 0; тогда из условия гладкости экрана
m^» ij ,/'2
получим lim maxMQ (supp(Vi ,;2)) = 0.
mij ,/'2
В [3] доказаны следующие утверждения.
Утверждение 1. Введенные функции V; и V; удовлетворяют условию
аппроксимации в L,(Q) и H (Q) соответственно.
Утверждение 2. Метод Галеркина с выбранными базисными функциями сходится.
60
University proceedings. Volga region
№ 4 (32), 2014
Физико-математические науки. Математика
Неизвестные коэффициенты clu, сф в разложениях (8) будем искать согласно методу Галеркина из системы линейных алгебраических уравнений
[L]u = uq ,
где u = (Си,..,спи ,сф’--’Сер ) - вектор-столбец неизвестных коэффициентов;
[L] - основная матрица СЛАУ; uq - столбец правой части.
Расширенную матрицу системы запишем в блочном виде
La = '[ A] [ Ki] [u0,Q]
[ K2] [ S ] [u0,Q ]_
Элементы матрицы вычисляются согласно методу Галеркина по формулам (для удобства в определении матричных элементов принята сквозная (одинарная) нумерация базисных функций и их носителей вместо тройной и двойной нумерации для тела и экрана соответственно):
[А]. = ус1(Пг.)8у - j j (k2(y)-k0)G(x,y)dydx, п.п.
. J
[K ]j =- j j G(x, y)dsydx, [K2 ]j =- j j (k2(y) -k0)G(x, y)dydsx,
[S]ij =- j j G(x, y)dsydsx, [u0 Q ]. = j u0,Q (x)dx, [u0 Q ]. = ju0 Q (x)dsx
aiWj
“i
Здесь уо1(П. ) - объем носителя i -й базисной функции на теле Q.
Если экран Q является плоским, то интеграл по поверхности представляет собой двойной интеграл Римана по плоской области. Если экран не плоский и задан параметрически, то поверхностный интеграл выражается через двойной стандартным образом.
Метод Галеркина является достаточно трудоемким, прежде всего из-за необходимости вычисления большого количества многомерных интегралов. В реализации алгоритма, предложенной авторами, была учтена симметрия задачи в случае, когда тело является однородным (или описывается кусочнопостоянной функцией k(x)), а экран - плоским. Так как ядра интегральных операторов зависят от разности аргументов, то основная матрица СЛАУ является блочно-теплицевой, это позволяет сократить объем вычислений матричных элементов (число различных элементов N' ~ N, N - порядок матрицы) и провести эксперименты с достаточно большим числом базисных функций. В табл. 1 проведено сравнение времени (в секундах), затрачиваемого на заполнение основной матрицы СЛАУ без учета и с учетом симметрии.
Для приближенного вычисления многомерных интегралов используется кубатурная формула центральных прямоугольников.
Отметим, что все интегралы в блоках [KJ, [K2], [uoi], [uo2]
являются собственными интегралами от гладких функций, ведь G е C, так
Physical and mathematical sciences. Mathematics
61
Известия высших учебных заведений. Поволжский регион
как экран и тело не имеют по предположению общих точек, а k(x),
Uq(x)e C(столь жесткие условия используются в [1], [3] при теоретическом исследовании задачи, а в расчетах могут быть существенно ослаблены).
Таблица 1
n, m t, с tsymm, C
n = 8, m = 4 0,171 0,016
n = 64, m = 16 19,39 0,436
n = 216, m = 36 177,9 2,496
n = 512, m = 64 518 11,12
n = 1000, m = 100 2003 79,63
Если носители базисных функций на теле (на экране) имеют непустое пересечение, то интегралы в блоках [A] ([S']) являются слабосингулярными и могут быть приближенно вычислены по формуле прямоугольников с отбрасыванием тех узлов кубатурной формулы, которые являются особыми точками функции G( x, y).
В табл. 2 приведены результаты приближенного вычисления интеграла
J J G(x, y)dydx ПП j
при различном числе узлов формулы центральных прямоугольников; в проведенном эксперименте в качестве Q выбран единичный куб, n = 4,
i = j = 1; N - число элементов разбиения П по каждому измерению (общее число узлов кубатурной формулы для шестикратного интеграла без учета особых точек составляет N6 ).
Таблица 2
N Значение интеграла Время, с
4 (0,0004078; 0,00005708) 0,003
8 (0,0004275; 0,00005787) 0,213
12 (0,0004311; 0,00005794) 1,833
16 (0,0004324; 0,00005796) 10,34
20 (0,0004329; 0,00005797) 39,64
Из табл. 2 видно, что вычисление интегралов с приемлемой точностью является достаточно трудоемким. Однако количество несобственных интегралов в матрице в сравнении с общим числом матричных элементов невелико. Учет структуры матрицы позволяет вычислять фиксированное количество различных слабосингулярных интегралов независимо от порядка матрицы. Так, для тела это случаи совпадения носителей, а также пересечения их по одной из шести граней, одному из двенадцати ребер или одной из восьми вершин; для экрана - случаи совпадения носителей, пересечения по одной из четырех сторон или вершин прямоугольного носителя.
62
University proceedings. Volga region
№ 4 (32), 2014
Физико-математические науки. Математика
3. Описание вычислительных экспериментов
Результаты расчетов
На рис. 1-3 представлены результаты приближенного решения системы (3) задачи дифракции на двух рассеивателях: кубе и ограниченном экране. На рис. 1(а)-3(а) изображаются значения модуля плотности (| ф(х) |) на экране и модуля решения в объемном теле (\u(x) |) в различных горизонтальных сечениях куба - при Х3 =1, Х3 = 0,66, Х3 = 0,33 и Х3 = 0. В проведенных вычислительных экспериментах 1-3 в методе Галеркина использованы описанные выше кусочно-постоянные базисные функции, количество базисных функций выбиралось одинаково, а именно положено n = 16, mk = 40. В подписях к рисункам указаны диапазоны изменения абсолютных величин полей на экране и теле. Для рис. 1,а, 2,а, 3,а дополнительно показаны две шкалы значений модуля решений на экране и теле.
а)
в)
Рис. 1. Общий вид модуля решения на системе рассеивателей (а), значения | и | на верхней грани (Х3 =1) куба Q (б) и значения | ф | на прямоугольном экране ^ (в);
диапазоны изменения - (0,87;1,58) и (2,26;35,59) соответственно; U0 = e™3
Physical and mathematical sciences. Mathematics
63
Известия высших учебных заведений. Поволжский регион
1.6
1.3
0.7
0.4
а)
в)
Рис. 2. Общий вид модуля решения на системе рассеивателей (а), значения | и | на верхней грани ( x3 =1) куба Q (б) и значения | ф | на прямоугольном экране Q (в);
диапазоны изменения - (0,44; 1,51) и (2,28; 39,38) соответственно; и0 = еХз
Фиксируем значения волнового числа - ke = 1. Будем рассматривать безразмерные координаты (нормированные на ke ).
На рис. 1 и 2 представлены результаты решения задачи, в которой рассеивающая структура образована однородным единичным кубом Q :={x:Xj е [0;1]} (k(x) = const = 3) и расположенного над ним плоского квадратного экрана.
В первом эксперименте экран Q := {x : Х1, Х2 е [0;1], Х3 =2} расположен строго над телом; падающая волна распространяется вдоль оси 0x3 . На рис. 1 видна симметрия решений относительно оси x1 = x2=0,5, ожидаемая для рассмотренных условий экспериментов.
При сдвиге экрана относительно тела (рис. 2) наблюдается нарушение симметрии решения как на экране (представлен вид на него сверху), так и на теле - в точках тела, расположенных «ближе» к экрану, поле ослабевает сильнее, что также соответствует физической сути задачи.
Из результатов, представленных на рисунках, видно, что значения | ф(x) | возрастают при стремлении x к границе экрана - это соответствует условию принадлежности решения ф(x) уравнения на экране с краем пространству Соболева H (Q).
64
University proceedings. Volga region
№ 4 (32), 2014
Физико-математические науки. Математика
в)
Рис. 3. Общий вид модуля решения на системе рассеивателей (а), значения | и | на нижней грани (Х3 = -0,5 ) куба Q (б) и значения | ф | на цилиндре Q (в);
диапазоны изменения - (0,42; 0,45) и (0,76; 10,17) соответственно; uo = e*3
На рис. 3 приведен результат расчета для системы рассеивателей «куб-цилиндр» (цилиндр изображен не полностью из соображений наглядности). Здесь рассеиватели заданы следующим образом: однородный куб
2 2
Q = {x: Xj е [—0,5; 0,5]}, цилиндр Q = {x: *i + *2 =1*3 е[—1,2]}, падающая
ix3
волна распространяется вдоль оси 0X3 : U0 = e 3. Возрастание решения на экране проявляется уже в различной степени на нижней и верхней границах поверхности цилиндра: со стороны набегания волны U0 (на нижней границе) решение растет быстрее (см. рис 3,в), однако значения | ф| у обоих краев существенно меньше, чем в предыдущих экспериментах - сказывается отсутствие углов у краев цилиндрического экрана. Симметрия поля в Q сохраняется, однако качественное поведение поля вблизи нижней границы экрана иное - в горизонтальных оно растет при стремлении точки к краям сечения (см. рис. 3,б).
В последнем вычислительном эксперименте рассеивающая структура состоит из неполного тора. Численное решение такой задачи дифракции
Physical and mathematical sciences. Mathematics
65
Известия высших учебных заведений. Поволжский регион
приводит к решению СЛАУ с матрицей, отвечающей только «поверхностной» части матричного оператора. Расширенная матрица имеет вид
[[s] ад],
матричные элементы определяются описанными ранее соотношениями.
На рис. 4,а,б изображаются значения модуля плотности поверхностного потенциала (| ф(х) |) на тороидальном экране при различном числе базисных
функций - расчеты проведены при т-у = m2 = 30, и Ш\ = m2 = 70.
а)
б)
Рис. 4. Значения | ф | на торе О при mt = 30 (а), mt = 70 (б); диапазоны изменения значений |ф| - (0,74;6,7) и (0,78; 10,9) соответственно; u0 = e*3
Заключение
В статье рассмотрена скалярная задача рассеяния плоской волны рассеивателями сложной структуры и метод Галеркина для ее приближенного решения. Описаны базисные функции на объемном теле и параметрически заданном ограниченном экране, приведены расчетные формулы для элемен-
66
University proceedings. Volga region
№ 4 (32), 2014
Физико-математические науки. Математика
тов матрицы системы линейных уравнений. Приведены в графическом виде результаты расчетов: рассмотрены случаи однородного и неоднородного тел, а также плоского и параметрически заданного неплоского экрана. Проведенные эксперименты позволили выявить симметрию решения задачи, а также влияние экрана, неоднородности тела и направления падающего поля на величину значений решения и его качественное поведение.
Список литературы
1. Медведик, М. Ю. Скалярная задача дифракции плоской волны на системе непересекающихся экранов и неоднородных тел / М. Ю. Медведик, Ю. Г. Смирнов, А. А. Цупак // Журнал вычислительной математики и математической физики. - 2014. - Т. 54, № 8. - С. 1319-1331.
2. Medvedik, M. Yu. Scalar Problem of Plane Wave Diffraction by a System of Nonintersecting Screens and Inhomogeneous Bodies / M. Yu. Medvedik, Yu. G. Smirnov, and A. A. Tsupak // Computational Mathematics and Mathematical Physics. - 2014. -Vol. 54, № 8. - P. 1280-1292.
3. Смирнов, Ю. Г. Метод интегральных уравнений в скалярной задаче дифракции, состоящей из «мягкого» и «жесткого» экранов и неоднородного тела / Ю. Г. Смирнов, А. А. Цупак // Дифференциальные уравнения. - 2014. - Т. 50, № 9. -С. 1164-1174.
4. Smirnov, Yu. G. Method of Integral Equations in the Scalar Problem of Diffraction on a System Consisting of a «Soft» and a «Hard» Screen and an Inhomogeneous Body / Yu. G. Smirnov, A. A. Tsupak // Differential Equations. - 2014. - Vol. 50, № 9. -P. 1150-1160.
5. Цуп ак, А. А. О единственности решения задачи дифракции акустической волны на системе непересекающихся экранов и неоднородных тел / А. А. Цупак // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2014. - № 1 (29). - С. 30-38.
References
1. Medvedik M. Yu., Smirnov Yu. G., Tsupak A. A. Zhurnal vychislitel'noy matematiki i matematicheskoyfiziki [The journal of calculus mathematics and mathematical physics]. 2014, vol. 54, no. 8, pp. 1319-1331.
2. Medvedik M. Yu., Smirnov Yu. G. and Tsupak A. A. Computational Mathematics and Mathematical Physics. 2014, vol. 54, no. 8, pp. 1280-1292.
3. Smirnov Yu. G., Tsupak A. A. Differentsial'nye uravneniya [Differential equations]. 2014, vol. 50, no. 9, pp. 1164-1174.
4. Smirnov Yu. G., Tsupak A. A. Differential Equations. 2014, vol. 50, no. 9, pp. 11501160.
5. Tsupak A. A. Izvestiya vysshikh uchebnykh zavedeniy. Povolzhskiy region. Fiziko-matematicheskie nauki [University proceedings. Volga region. Physics and mathematics sciences]. 2014, no. 1 (29), pp. 30-38.
Деревянчук Екатерина Дмитриевна
аспирант, Пензенский государственный университет (Россия, г. Пенза, ул. Красная, 40)
Derevyanch.uk Ekaterina Dmitrievna
Postgraduate student, Penza State University (40 Krasnaya street,
Penza, Russia)
E-mail: [email protected]
Physical and mathematical sciences. Mathematics
67
Известия высших учебных заведений. Поволжский регион
Смолькин Евгений Юрьевич аспирант, Пензенский государственный университет (Россия, г. Пенза, ул. Красная, 40)
E-mail: [email protected]
Цупак Алексей Александрович кандидат физико-математических наук, доцент, кафедра математики и суперкомпьютерного моделирования, Пензенский государственный университет (Россия, г. Пенза, ул. Красная, 40)
E-mail: [email protected]
Smol'kin Evgeniy Yur'evich Postgraduate student, Penza State University (40 Krasnaya street, Penza, Russia)
Tsupak Aleksey Aleksandrovich Candidate of physical and mathematical sciences, associate professor, sub-department of mathematics and supercomputer modeling, Penza State University (40 Krasnaya street, Penza, Russia)
УДК 517.3 Деревянчук, Е. Д.
Метод Галеркина решения скалярной задачи рассеяния препятствием сложной формы / Е. Д. Деревянчук, Е. Ю. Смолькин, А. А. Цупак // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2014. - № 4 (32). - С. 57-68.
68
University proceedings. Volga region
