Акустический журнал, 2019, T. 65, № 4, стр. 448-459

ВВЕДЕНИЕ

Дифракция волн на периодических структурах представляет большой практический интерес [1–5]. В литературе имеется также много работ, посвященных дифракции на неровных поверхностях (в том числе периодических) [6–13]. В указанных работах рассматриваются асимптотические и численные методы решения задачи дифракции на неровных поверхностях, например, метод касательной плоскости, метод малых наклонов, метод поверхностных интегральных уравнений и другие. В большинстве работ рассматриваются двумерные задачи дифракции. Однако более адекватной моделью является трехмерная двоякопериодическая поверхность. В настоящей работе рассматривается задача дифракции плоской волны на акустически мягкой или жесткой двоякопериодической волнистой поверхности. Задача решается двумя методами: модифицированным методом дискретных источников (ММДИ) и методом диаграммных уравнений (МДУ). По-существу, данная работа является обобщением работ [11] и [13], посвященных двумерной задаче дифракции на периодической поверхности. Отметим, что МДУ применим к так называемым слабоневыпуклым поверхностям, в то время как ММДИ может быть применен к существенно менее “пологим” поверхностям. Также отметим, что одним из преимуществ МДУ является, во-первых, простота алгебраизации задачи. Во-вторых, в некоторых частных случаях матричные элементы соответствующей алгебраической системы выражаются в явном виде, то есть не требуют большого объема вычислений.

Как известно (см., например, [13]), в случае применения ММДИ задача дифракции сводится к решению интегрального уравнения относительно некоторой неизвестной функции, распределенной на центральном периоде вспомогательной поверхности, расположенной под исходной (волна падает сверху). При этом ядро интегрального уравнения выражается через периодическую функцию Грина (ФГ). Известно, что двойной ряд периодической ФГ сходится чрезвычайно медленно при условии, что точка наблюдения и точка источника близки друг к другу. С использованием формулы суммирования Пуассона ряд ФГ преобразуется в ряд по гармоникам Флоке. Однако сходимость этого ряда не лучше, чем у исходного, при указанном выше условии [5]. С целью преодоления этой трудности в работе применялся алгоритм вычисления периодической ФГ, предложенный в работе [5], в которой рассмотрена задача дифракции на двоякопериодической решетке, состоящей из тел вращения.

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

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим математическую постановку задачи дифракции на двоякопериодической поверхности. Введем декартову систему координат, в которой уравнение поверхности есть $z = f(x,y),$ причем $f(x + {{d}_{x}},y)$ = $f(x,y + {{d}_{y}}) = f(x,y).$ Начало координат выберем в середине какого-либо периода поверхности (назовем его центральным периодом поверхности и обозначим через ${{S}_{{00}}}$). Считаем, что структура облучается плоской волной вида:

где U ⁰ – давление падающей волны, $k$ – волновое число, $(r,{\theta },{\varphi })$ – сферические координаты, ${{{\theta }}_{0}},{{{\varphi }}_{0}}$ – углы падения. Предполагаем, что область над поверхностью (при $z > f(x,y)$) заполнена однородной средой так, что рассеянное поле ${{U}^{1}}$ в этой области удовлетворяет однородному уравнению Гельмгольца:

На поверхности выполнено условие Дирихле:

либо Неймана где ${\partial \mathord{\left/ {\vphantom {\partial {\partial n}}} \right. \kern-0em} {\partial n}}$ – производная по нормали, составляющей острый угол с осью z, $U = {{U}^{0}} + {{U}^{1}}$ – полное поле над поверхностью.

Вторичное поле в области над поверхностью удовлетворяет условиям периодичности Флоке:

где ${{h}_{x}} = k{{d}_{x}}\sin {{{\theta }}_{0}}\cos {{{\varphi }}_{0}},$ ${{h}_{y}} = k{{d}_{y}}\sin {{{\theta }}_{0}}\sin {{{\varphi }}_{0}}$ – параметры Флоке. На бесконечности полное поле удовлетворяет условию излучения: где ${{\kappa }_{{pq}}} = {{u}_{p}}{{i}_{x}} + {{v}_{q}}{{i}_{y}} + {{\Gamma }_{{pq}}}{{i}_{z}},$ ${{u}_{p}} = \frac{{{{h}_{x}} + 2{\pi }p}}{{{{d}_{x}}}},$ ${{v}_{q}} = \frac{{{{h}_{y}} + 2{\pi }q}}{{{{d}_{y}}}},$ ${{\Gamma }_{{pq}}} = \sqrt {{{k}^{2}} - u_{p}^{2} - v_{q}^{2}} ,$ причем знак квадратного корня выбирается из условия неположительности его мнимой части. В формуле (7) ${{z}_{{\max }}}$ – максимальное значение координаты $z$ точки на поверхности, $A_{{pq}}^{{}}$ – неизвестные коэффициенты.

2. РЕШЕНИЕ ЗАДАЧИ ПРИ ПОМОЩИ ММДИ

В соответствии с ММДИ представим волновое поле над периодической поверхностью в виде [13–15]:

В этой формуле $\Sigma $ – центральный период вспомогательной поверхности, расположенной под исходной поверхностью (выбор вспомогательной поверхности описан ниже), а $J$ – неизвестная функция, заданная на поверхности $\Sigma $. Функция $G$ представляет собой периодическую ФГ, которая имеет вид [5]:

где

Таким образом, рассеянное поле, записанное в виде (8), удовлетворяет условиям периодичности Флоке. Далее подставим формулу (8) в краевое условие (3) или (4). В результате задача сводится к решению интегрального уравнения

либо соответственно в случае условия Дирихле или Неймана на периодической поверхности. Ядра интегральных уравнений (11) и (12) выражаются через периодическую ФГ $G(r,r{\text{'}}),$ которая вычисляется по следующей формуле [5]: где причем $R,\Theta ,\Phi $ – сферические координаты вектора $r - r{\text{',}}$ $P_{s}^{l}(x)$ – присоединенные функции Лежандра, ${{j}_{s}}(x)$ – сферические функции Бесселя. Величины $W_{{sl}}^{{}}$ в (13) выражаются по формулам, приведенным в работе [5].

Сделаем замечание по поводу разрешимости интегрального уравнения (11) или (12). Как показано в работе [16] (см. также [15, 17]), имеет место следующая теорема существования (рассмотрим ее на примере решения задачи дифракции на компактном теле с граничным условием Дирихле на поверхности тела).

Пусть простая замкнутая поверхность S такова, что $k$ не является собственным значением внутренней однородной задачи Дирихле для области внутри S. Тогда интегральное уравнение вида (11) $\left( {{\text{г д е }}\,\,G \equiv {{G}_{0}}(R),\,\,R = \left| {{\mathbf{\vec {r}}} - {\mathbf{\vec {r}}}'} \right|} \right)$ разрешимо в том и только том случае, если $\Sigma $ охватывает все особенности решения ${{U}^{1}}(\vec {r})$ краевой задачи с граничным условием Дирихле на поверхности рассеивателя. В этом случае уравнение (11) имеет единственное решение.

Аналогичная теорема будет справедлива и для задачи дифракции на периодической поверхности. В последнем случае вспомогательная поверхность $\Sigma $ должна располагаться выше особенностей продолжения волнового поля под поверхность (плоская волна падает сверху). Отметим, что в приведенной теореме не указывается способ построения вспомогательной поверхности, нужно лишь, чтобы эта поверхность охватывала все особенности продолжения волнового поля внутрь границы тела или под периодическую поверхность. В работах [13–15] было показано, что для получения наиболее эффективных и устойчивых численных алгоритмов необходимо выбирать вспомогательную поверхность при помощи аналитической деформации границы рассеивателя.

Рассмотрим выбор вспомогательной поверхности на примере решения двумерной задачи дифракции плоской волны на неровной периодической поверхности, заданной уравнением $z = f(x).$ В соответствии с техникой ММДИ вводим комплексную переменную [13–15]:

где $\delta > 0.$ В случае если величина $\delta $ равна нулю, то переменная $\zeta \in C,$ где $C$ контур на комплексной плоскости $\zeta $, соответствующий контуру периодической поверхности (назовем его $S$) и конгруэнтный ему. Если начать увеличивать δ, то $C$ будет деформироваться и лежать в области, расположенной под контуром волнистой поверхности $S$. При этом мы получим новый контур $\Sigma $, который можно выбрать в качестве вспомогательного контура. Заметим, что такого рода деформация контура $S$ возможна до тех пор, пока отображение $\zeta (x) = x + if(x)$ остается взаимно однозначным. Если рассматривать параметры $x$ и δ в формуле (15) как некоторые криволинейные координаты, то соответствующие координатные линии будут ортогональны. Это означает, что при увеличении параметра δ вспомогательные источники движутся по траекториям, ортогональным к контуру $S$. Данный факт обеспечивает устойчивость методики, основанной на ММДИ (то есть возможность решать уравнение первого рода с гладким ядром).

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

где ${{f}_{l}}(x)$ и ${{g}_{t}}(y)$ – периодические аналитические функции. Например, уравнение поверхности $z = f(x,y)$ может быть представлено отрезком двойного ряда Фурье. Для построения вспомогательной поверхности рассмотрим следующую величину: где $i,j,k$ – мнимые кватернионные единицы, $\delta $ – положительный параметр. Таким образом, величина $\zeta $ представляет собой функцию от двух кватернионных переменных [18, 19]. При этом, если записать кватернион $\zeta $ в виде то в качестве координат точки на вспомогательной поверхности следует положить

Например, для поверхности f(x, y) = = $A\sin {{\tau }_{1}}x\sin {{\tau }_{2}}y,$ где ${{\tau }_{1}} = \frac{{2\pi }}{{{{d}_{x}}}},$ ${{\tau }_{2}} = \frac{{2\pi }}{{{{d}_{y}}}},$ используя формулы

а также известные формулы для произведения мнимых кватернионных единиц [18, 19], будем иметь следующие параметрические уравнения вспомогательной поверхности

Заметим, что при дифференцировании уравнений вспомогательной поверхности по переменной δ мы получим, что в первом приближении по δ точки вспомогательной поверхности будут двигаться по ортогональным к поверхности $\Sigma $ траекториям. Этот факт обеспечивает устойчивость алгоритма, как и в случае применения стандартного ММДИ.

Необходимо также заметить, что предлагаемый подход переходит в двумерный вариант ММДИ при условии

Действительно, в этом случае

Если обозначить $f(x + j\delta )$ = $u(x,\delta ) + jv(x,\delta ),$ то в соответствии с (19) получим

Этот же результат получается из формулы (15), так как действительная и мнимая части комплексной переменной $\zeta (x)$ совпадают с выражениями для ${{x}_{\Sigma }}$ и ${{z}_{\Sigma }}$ в формуле (24) (при этом, очевидно, $f(x - i\delta ) = u(x,\delta ) - iv(x,\delta )$).

Для численного решения интегрального уравнения (11) или (12) применим метод коллокации. С этой целью поверхностные интегралы в (11) и (12) заменяем на двойные интегралы по прямоугольнику $[{{ - {{d}_{x}}} \mathord{\left/ {\vphantom {{ - {{d}_{x}}} 2}} \right. \kern-0em} 2},{{{{d}_{x}}} \mathord{\left/ {\vphantom {{{{d}_{x}}} 2}} \right. \kern-0em} 2}] \times [{{ - {{d}_{y}}} \mathord{\left/ {\vphantom {{ - {{d}_{y}}} 2}} \right. \kern-0em} 2},{{{{d}_{y}}} \mathord{\left/ {\vphantom {{{{d}_{y}}} 2}} \right. \kern-0em} 2}].$ При этом сделаем замену неизвестной функции $J$ по формуле [5]:

В формуле (25) $r{\text{'}}$ – радиус-вектор точки на вспомогательной поверхности, точка означает дифференцирование по параметру x или y. Далее выберем в области интегрирования прямоугольную сетку:

Заменим двойные интегралы суммами Римана и приравняем полученные выражения в точках коллокации, которые выберем следующим образом:

В результате перейдем от интегрального уравнения к следующей алгебраической системе относительно неизвестных значений функции $I$ в точках $r_{{nm}}^{{}}$ на вспомогательной поверхности центрального элемента поверхности:

Здесь

в случае граничного условия (3) и в случае условия (4). В формуле (28) ${{c}_{{nm}}} = I({{r}_{{nm}}})\sigma ,$ где $\sigma $ – площадь элементарных ячеек, на которые разбивается область интегрирования в уравнении (11) или (12). Точки $r_{{nm}}^{{}}$ в (29) и (30) выбираются с использованием формул (17)–(19).

В качестве выходных характеристик задачи рассмотрим амплитуды пространственных гармоник, обусловленных дифракцией на неровной поверхности. Как описано в работе [5], амплитуды плоских волн, распространяющихся в верхнем полупространстве, будут иметь вид:

3. РЕШЕНИЕ ЗАДАЧИ ПРИ ПОМОЩИ МДУ

Рассмотрим решение задачи дифракции плоской волны на периодически неровной поверхности при помощи МДУ [11, 15, 20]. Основная идея МДУ состоит в том, что исходная задача дифракции сводится к алгебраической системе не относительно значений поля или его нормальной производной на неровной поверхности (или на вспомогательной поверхности), а относительно амплитуд отраженных плоских волн. С целью решения задачи дифракции при помощи МДУ используем интегральное представление рассеяного поля:

в случае условия Дирихле, либо в случае условия Неймана. Соотношения (32) и (33) получаются из формулы Грина. При этом запишем периодическую ФГ, используя формулу суммирования Пуассона, в виде

Считая вначале, что поверхность рэлеевская (полагаем $\left| {z - z{\text{'}}} \right| = z - z{\text{'}},$ то есть $z > z{\text{'}}$), подставим выражение (34) в формулу (32) или (33). В результате перемены местами порядка интегрирования и суммирования, получим следующее представление рассеянного поля [11, 15]:

где ${{g}_{{pq}}} = g({{u}_{p}},{{v}_{q}}),$ причем в случае абсолютно мягкой и в случае абсолютно жесткой поверхности. В формулах (36) и (38) $g(u,v)$ – диаграмма центрального периода поверхности, $\Gamma = \sqrt {{{k}^{2}} - {{u}^{2}} - {{v}^{2}}} .$ Таким образом, мы получаем представление (35) для волнового поля в области над периодической поверхностью в виде суммы уходящих плоских волн. Подставим далее полное поле $U = {{U}^{0}} + {{U}^{1}},$ где рассеянное поле имеет вид (35) в формулу (36) или (38). В результате будем иметь где в случае условия Дирихле, в случае условия Неймана. В формулах (41) и (43) учтено, что падающее поле имеет вид ${{U}^{0}} = \exp ( - i{{u}_{0}}x$ – $i{{v}_{0}}y + i{{\Gamma }_{{00}}}z).$ Положив в формуле (39) $u = {{u}_{m}},$ $v = {{v}_{n}},$ $m,n = 0, \pm 1, \pm 2,...$, получим алгебраическую систему в которой $A_{{mn}}^{{pq}} = {{F}_{{pq}}}({{u}_{m}},{{v}_{n}}),$ $g_{{mn}}^{0} = {{g}^{0}}({{u}_{m}},{{v}_{n}}).$ Как видно из формулы (35) и из системы (44), исходная задача дифракции сводится к нахождению амплитуд отраженных плоских волн, то есть величин ${{g}_{{mn}}}.$ Этот факт очень удобен при вычислениях, так как именно эти величины представляют интерес в приложениях. Отметим, что при численных расчетах система (44) решалась методом редукции, причем $m,p = \overline { - {{M}_{1}},{{M}_{1}}} ,$ $n,q = \overline { - {{N}_{1}},{{N}_{1}}} .$

Сделаем замечание относительно применимости МДУ к исследуемой задаче дифракции. Рассмотрим случай абсолютно мягкой поверхности. Можно показать, что имеют место следующие оценки матричных элементов $A_{{mn}}^{{pq}},$ получаемые применением метода перевала к соответствующим интегралам (см., например, [11, 15]):

где $\lambda = \frac{{2\pi }}{k}$ – длина волны. При этом пары комплексных величин $({{x}_{{01}}},{{y}_{{01}}})$ и $({{x}_{{02}}},{{y}_{{02}}})$ являются корнями следующих систем уравнений

В формулах (49) и (50) максимумы и минимумы берутся по всем корням систем уравнений (51), соответствующих как верхним, так и нижним знакам.

Сделаем замену неизвестных в системе (44) по формуле

где

Подставим формулу (52) в алгебраическую систему (44). В результате получим

где

Из последней формулы и из оценок (45)–(48) видно, что с ростом индексов новые матричные элементы экспоненциально убывают при $p \to \infty $ или $q \to \infty $ при условии

В случае $m \to \infty $ или $n \to \infty $ убывание имеет степенной характер. В случае знака равенства в формуле (56) матричные элементы системы (54) имеют порядок $O\left( {{1 \mathord{\left/ {\vphantom {1 {{{{\left| p \right|}}^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}}}} \right. \kern-0em} {{{{\left| p \right|}}^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}}}} \right)$ или $O\left( {{1 \mathord{\left/ {\vphantom {1 {{{{\left| q \right|}}^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}}}} \right. \kern-0em} {{{{\left| q \right|}}^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}}}} \right)$ при $p \to \infty $ или $q \to \infty $ соответственно. Таким образом, система (44) будет разрешима методом редукции при ограничениях (56). Заметим, что область применимости МДУ, определяемая условием (56), шире, чем, например, область применимости методики, основанной на предположении о справедливости гипотезы Рэлея для неровной поверхности [11, 15].

Для численной оценки матричных элементов $A_{{mn}}^{{pq}}$ и правых частей $g_{{mn}}^{0},$ определяемых по формулам (40)–(43), необходимо использовать метод интегрирования по частям и замену $\xi = \frac{{2\pi }}{{{{d}_{x}}}}x,$ $\eta = \frac{{2\pi }}{{{{d}_{y}}}}y.$ В результате соответствующие интегралы сводятся к следующим:

при ${{\Gamma }_{{mn}}} \ne {{\Gamma }_{{pq}}},$ либо при ${{\Gamma }_{{mn}}} = {{\Gamma }_{{pq}}},$

В формулах (57)–(59) обозначено

в случае абсолютно мягкой поверхности и

в случае абсолютно жесткой.

Как указано во Введении, в некоторых частных случаях интегралы (57) и (59) могут быть вычислены аналитически. Например, если поверхность задана уравнением:

то двойные интегралы в формулах (57) и (59) примут вид [21]: где через $\Omega $ обозначено ${{\Gamma }_{{mn}}} - {{\Gamma }_{{pq}}},$ либо ${{\Gamma }_{{mn}}} + {{\Gamma }_{{00}}},$ ${{J}_{s}}(x)$ – функции Бесселя.

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ

В качестве примера неровной двоякопериодической поверхности была рассмотрена синусоидальная поверхность вида (62). Ниже всюду исследуется случай абсолютно мягкой поверхности. Рассмотрим вначале результаты, касающиеся проверки корректности предлагаемых методов. На рис. 1 представлено распределение невязки краевого условия (3) на центральном периоде поверхности (62), полученной при помощи ММДИ. Невязка определяется по формуле

Здесь $r_{{\nu \mu }}^{{{\text{**}}}}$ – радиус-векторы точек на поверхности центрального периода, расположенные посередине между точками коллокации. Параметры поверхности были следующими: ${{d}_{x}} = 2{{d}_{y}},$ ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.5,$ ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda } = 1.2,$ ${{\theta }_{0}} = 45^\circ ,$ ${{\varphi }_{0}} = 0,$ где $d = \min ({{d}_{x}},{{d}_{y}}).$ Число дискретных источников (точек коллокации) было выбрано равным $N = 55,$ $M = 35.$ Параметр δ выбирался так, чтобы вспомогательная поверхность не имела точек самопересечения, и был равен $\delta = 0.2.$ Как видно из рисунка, максимальный уровень невязки не превосходит 0.01, то есть достаточно мал. В качестве еще одной проверки результатов, получаемых при помощи ММДИ, была найдена точность выполнения закона сохранения энергии, который в рассматриваемом случае имеет вид:

где ${{P}_{0}} = k\cos {{\theta }_{0}}$ – мощность падающей плоской волны (предполагается, что амплитуда падающей волны равна 1 Па). В формуле (65) суммирование проводится по всем распространяющимся гармоникам Флоке. Результаты показывают, что точность выполнения закона сохранения при выбранных параметрах модели составила $1.39 \times {{10}^{{ - 4}}}.$ Таким образом, имеет место высокая точность расчета амплитуд отраженных гармоник при наличии достаточно больших значений отношения ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d}$ неровной поверхности.

В таблице 1 представлены результаты расчета амплитуд трех первых отраженных гармоник в зависимости от угла падения ${{\theta }_{0}}$ плоской волны, полученные при помощи МДУ и ММДИ. Параметры задачи имели значения ${{d}_{x}} = {{d}_{y}},$ ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09,$ ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda } = 1.2,$ ${{\varphi }_{0}} = 0.$ Как было указано выше, МДУ применим для достаточно “пологих” поверхностей, поэтому в данном случае отношение ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d}$ невелико. Отметим, что для достижения данной точности при использовании МДУ требуется относительно небольшое число базисных функций: ${{M}_{1}} = {{N}_{1}} = 15.$ Из таблицы видно, что результаты, полученные обоими методами, отличаются не более чем на $1.4 \times {{10}^{{ - 4}}}.$ Данный факт подтверждает корректность обоих методов. В таблице 2 приведены значения точности выполнения закона сохранения энергии, полученные при помощи ММДИ и МДУ для двух значений отношения ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d}$. Остальные параметры задачи были такие же, как и для табл. 1. Как видно из таблицы, в случае ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09$ относительная разность левой и правой частей в формуле (65) не превосходит 5 × 10^–5 при использовании обоих методов. В то же время при увеличении отношения ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d}$ до 0.11 эта величина имеет максимальное значение, равное примерно 0.1 при использовании МДУ, а в случае применения ММДИ точность выполнения закона сохранения энергии остается примерно такой же, как и для отношения ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09$. Расчеты показывают, что предельное значение отношения ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d}$, при котором МДУ дает приемлемую точность в широком диапазоне длин волн, составляет примерно 0.1. Этот факт согласуется с условиями применимости МДУ (56).

В табл. 3 и 4 приведены результаты, касающиеся внутренней сходимости методов на основе ММДИ и МДУ соответственно. Входные данные были следующие: ${{d}_{x}} = 2{{d}_{y}},$ ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda } = 1.2,$ ${{\theta }_{0}} = 45^\circ ,$ ${{\varphi }_{0}} = 0.$ Отношение ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.5$ для таблицы 3 и ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09$ для табл. 4. Как видно из приведенных данных, результаты расчетов стабилизируются в 4-м десятичном знаке (см. результаты расчета величины $\left| {{{A}_{{ - 2,0}}}} \right|$). Из таблицы также следует, что с увеличением числа базисных функций увеличивается точность выполнения закона сохранения энергии.

На рис. 2 и 3 изображены зависимости модулей амплитуд трех первых гармоник Флоке от параметра ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda }.$ Параметры задачи были следующие: ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09,$ углы падения волны ${{\theta }_{0}} = 45^\circ ,$ ${{\varphi }_{0}} = 0.$ Кривые 1–3 соответствуют амплитудам ${{A}_{{00}}},{{A}_{{ - 1,1}}},{{A}_{{ - 2,0}}}.$ Кривые на рисунках получены при помощи МДУ, который позволяет рассматривать случай достаточно высоких частот. Как видно из приведенных рисунков, имеются изломы кривых, соответствующие критическим значениям волнового параметра ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda },$ на которых возникают новые распространяющиеся (незатухающие) пространственные гармоники. Из рисунков следует, что модуль коэффициента отражения плоской волны, то есть величина $\left| {{{A}_{{00}}}} \right|,$ монотонно убывает, а величина $\left| {{{A}_{{ - 2,0}}}} \right|$ монотонно возрастает с увеличением параметра ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda },$ то есть происходит перераспределение энергии в высшие гармоники Флоке.

Рис. 4–7 иллюстрируют зависимости модулей отраженных гармоник указанных выше номеров от угла падения плоской волны ${{\theta }_{0}}$. Рисунки 4 и 5 соответствуют случаю “пологой” поверхности (при ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09$), а рисунки 6 и 7 – случаю, когда отношение ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.5.$ Параметры задачи были следующие: ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda } = 1.2,$ ${{\varphi }_{0}} = 0.$ Кривые 1–3 на рисунках вновь соответствуют амплитудам ${{A}_{{00}}},{{A}_{{ - 1,1}}},{{A}_{{ - 2,0}}}.$ Как следует из рисунков, в случае “пологой” поверхности зависимости модулей амплитуд ${{A}_{{ - 1,1}}},{{A}_{{ - 2,0}}}$ монотонны, при этом $\left| {{{A}_{{00}}}} \right|$ меняется незначительно. В случае, когда ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.5,$ имеется “провал” зависимости величины $\left| {{{A}_{{00}}}} \right|,$ то есть модуля коэффициента отражения плоской волны, соответствующий углу падения, примерно равному 42° (при ${{d}_{x}} = {{d}_{y}}$) и 45° (при ${{d}_{x}} = 2{{d}_{y}}$). При этом величина $\left| {{{A}_{{ - 1,1}}}} \right|$ имеет максимальное значение для этого угла падения волны.

ЗАКЛЮЧЕНИЕ

На основе ММДИ и МДУ разработаны два подхода для решения трехмерной задачи дифракции плоской волны на неровной двоякопериодической поверхности. Выведены соотношения, позволяющие рассчитывать амплитуды гармоник Флоке для абсолютно мягкой и абсолютно жесткой поверхности двумя методами. Численные результаты приведены для абсолютно мягкой синусоидальной поверхности. Для тестирования метода на основе ММДИ построено распределение невязки краевого условия на центральном периоде поверхности. Показано, что максимальный уровень невязки не превосходит ${{10}^{{ - 2}}}$ для значения относительной высоты поверхности, равной ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.5.$ Проведена оценка точности выполнения закона сохранения энергии при помощи обеих методик. Показано, что для значения параметра ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda } = 1.2$ точность выполнения закона сохранения составляет 4 × 10^–5 (при ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09$) либо 1.4 × 10^–4 (при ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.5$) при использовании ММДИ и 4 × 10^–5 при использовании МДУ (при ${a \mathord{\left/ {\vphantom {a d}} \right. \kern-0em} d} = 0.09,$ то есть в том случае, когда этот метод применим). Проведено сравнение результатов, полученных двумя методами, для случая “пологой” поверхности. Показано, что относительная разность результатов, полученных с использованием ММДИ и МДУ, не превосходит 1.4 × 10^–4. Приведены результаты, касающиеся внутренней сходимости обоих методов.

В работе построены зависимости модуля амплитуд трех первых распространяющихся гармоник Флоке от волнового параметра ${d \mathord{\left/ {\vphantom {d \lambda }} \right. \kern-0em} \lambda }$ и угла падения плоской волны для синусоидальной поверхности. Показано, что в случае “пологой” поверхности модуль коэффициента отражения плоской волны монотонно убывает с ростом волнового параметра, таким образом происходит перераспределение энергии в высшие гармоники. Показано также, что в случае достаточно большого значения амплитуды неровной поверхности имеется “провал” зависимости модуля коэффициента отражения плоской волны от угла падения.

Работа выполнена при поддержке Российского фонда фундаментальных исследований, проекты № 18-02-00961 и № 19-02-00654.