Computationally efficient solution for finding the current density on the illuminated and shadow sides of an infinitely thin circular disk

全文:

Введение

Одной из канонических задач теории дифракции в векторном (электромагнитном) случае является задача дифракции на бесконечно тонком идеально проводящем круглом диске [1]. Ее исследованию в отношении внешней и внутренней задач электродинамики для стороннего плоского монохроматического поля произвольной поляризации и произвольного направления падения посвящено большое число работ [2–10] и др. Их основу составляют следующие решения: 1) аналитические, формируемые при применении метода разделения переменных и представления волновых уравнений в вырожденной эллиптической системе координат с получением решения в виде ряда по сфероидальным функциям [2; 3]; 2) асимптотические, реализуемые в приближениях физической оптики [4], физической теории дифракции [5], геометрической теории дифракции и ее модификаций [6; 7]; 3) численные, основанные на применении метода моментов [8] или его модификаций [9; 10].

Указанные решения в исследовании дифракционной задачи на бесконечно тонком идеально проводящем круглом диске при применении метода конформных отображений [11; 12] потенциально позволяют обобщить получаемый результат на более сложные геометрические структуры [13–15]. Вместе с тем перечисленные методы в существующих реализациях не позволяют вычислительно эффективно выделить поверхностную плотность тока на освещенной и теневой сторонах круглого диска.

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

1. Постановка и решение задачи дифракции

Пусть $\Omega \subset {ℝ}^2=(x_3=0)\subset {ℝ}^3$ – бесконечно тонкий идеально проводящий диск радиусом $R$ с центром в начале координат и границей $\partial \Omega =\overline{\Omega }\backslash \Omega ,$ а ${\mathbf{E}}^0, {\mathbf{H}}^0$ – падающее стороннее плоское монохроматическое поле (рис. 1).

Рис. 1. Геометрическое представление задачи дифракции плоской волны на бесконечно тонком идеально проводящем круглом диске

Fig. 1. Geometric representation of the problem of diffraction of a plane wave on an infinitely thin perfectly conducting circular disk

Задачу дифракции ${\mathbf{E}}^0, {\mathbf{H}}^0$ на $\Omega$ сведем к определению рассеянного электромагнитного поля (ЭМП) [16]:

$\mathbf{E},\mathbf{H}\in C^2({ℝ}^3\backslash \overline{\Omega })\times \underset{\delta >0}{\cap }C\left(\overline{{ℝ}_{+}^3}\backslash \partial {\Omega }_{\delta }\right)\underset{\delta >0}{\cap }C\left(\overline{{ℝ}_{-}^3}\backslash \partial {\Omega }_{\delta }\right),$ (1)

удовлетворяющего условиям [1; 13]:

$\nabla \times \mathbf{H}=-i\beta \mathbf{E},\nabla \times \mathbf{E}=i\beta \mathbf{H},\mathbf{x}\in {ℝ}^3\backslash \overline{\Omega };$ (2)

${\left.{\mathbf{E}}_{\tau }\right|}_{\Omega }={\left.-{\mathbf{E}}_{\tau }^0\right|}_{\Omega };{\left.{\mathbf{E}}_{\tau }^0\right|}_{\Omega }\in C^{\infty }\left(\overline{\Omega }\right);\mathbf{E},\mathbf{H}\in L_{loc}^2\left({ℝ}^3\right);$

$\begin{array}{l}\mathbf{E}\mathbf{,}\mathbf{H}=o\left(r^{-1}\right),r:=\left|\mathbf{x}\right|\to \infty ,\mathrm{Im}\beta >0;\\ \mathbf{H}\times {\mathbf{e}}_r-\mathbf{E}=o\left(r^{-1}\right);\mathbf{E}\times {\mathbf{e}}_r+\mathbf{H}=o\left(r^{-1}\right);\\ \mathbf{E}\mathbf{,}\mathbf{H}=O\left(r^{-1}\right),r\to \infty при\mathrm{Im}\beta =\mathrm{0,}\end{array}$

где

$\mathbf{x}=\left\{x_1,x_2,x_3\right\};$ ${\mathbf{e}}_r=\mathbf{x}/\left|\mathbf{x}\right|;$

${\beta }^2={\omega }^2\mu (\epsilon +i\sigma {\omega }^{-1});$

$\omega >0$ – угловая частота; $\epsilon >0$ и $\mu >0$ – абсолютные диэлектрическая и магнитная проницаемости среды; индекс $\tau$ обозначает тангенциальную составляющую поля на $\Omega$,

$\partial {\Omega }_{\delta }=\left\{\mathbf{x}:\left|\mathbf{x}-\mathbf{y}\right|<\delta ,\mathbf{y}\in \partial \Omega \right\}.$

Доказательство существования и единственности решения задачи (1)–(2) при $\mathrm{Im}\beta \ge 0$ известно из [16, с. 45]. Используя векторные потенциалы и граничное условие для ${\mathbf{E}}_{\tau }$ из (2), представим (1)–(2) в виде

$\mathcal{L}\mathbf{J}=-\mathbf{f},$ (3)

где

$\mathcal{L}\mathbf{J}=\nabla \mathcal{A}(\nabla \cdot \mathbf{J})+{\beta }^2\mathcal{A}\mathbf{J};$ $\mathbf{f}=i\omega (\epsilon +i\sigma {\omega }^{-1}){\left.{\mathbf{E}}_{\tau }^0\right|}_{\overline{\Omega }};$

$\mathcal{A}\mathbf{J}\left(x\right)=\underset{\Omega }{\int }G(\mathbf{x},\mathbf{y})\mathbf{J}\left(\mathbf{y}\right)d\mathbf{y};$ $\mathbf{x}=\left\{x_1,x_2,x_3\right\}\in \overline{\Omega };$

$\mathbf{J}\left(\mathbf{y}\right)\cdot \mathbf{n}=0;$ $G(\mathbf{x},\mathbf{y})=\frac{1}{4\pi }\frac{e^{i\beta \left|\mathbf{x}-\mathbf{y}\right|}}{\left|\mathbf{x}-\mathbf{y}\right|}$

– функция Грина; $\mathbf{J}\left(\mathbf{y}\right)={\mathbf{J}}_{}^e\left(\mathbf{y}\right)+{\mathbf{J}}_{}^h\left(\mathbf{y}\right)$ – поверхностная плотность тока на $\Omega$ представленная суммой безвихревых ${\mathbf{J}}_{}^e\left(\mathbf{y}\right)$ и вихревых ${\mathbf{J}}_{}^h\left(\mathbf{y}\right)$ токов [17]; $\mathbf{n}=\left(\mathrm{0,0,1}\right)$ – орт нормали к $\Omega$.

Решение задачи (3) выполним в проекционной постановке метода Галеркина [18] при разложении:

$\tilde{\mathbf{J}}\left(\mathbf{x}\right)=\sum _{j=1}^{\mathcal{M}}{c}_j{\psi }_j\left(\mathbf{x}\right),$ (4)

искомой функции J по базису ${\mathbf{\psi }}_j\left(\mathbf{x}\right).$

С учетом первой формулы Грина [1] и свойств дифференциальных операторов [19] при удовлетворении ${\psi }_j$ граничным условиям ${\left.{\mathbf{\psi }}_j\cdot \mathbf{\nu }\right|}_{\partial \Omega }=\mathrm{0,}$ ($\mathbf{v}$ – нормаль к $\partial \Omega$) выражение (3) при требовании ортогональности

$\underset{\Omega }{\int }\mathcal{N}\left(x\right){\psi }_i\left(x\right)dx=0$

невязки

$\mathcal{N}\left(\mathbf{x}\right)=\mathcal{L}\left[\sum _{j=1}^{\mathcal{M}}{c}_j{\psi }_j\left(\mathbf{x}\right)\right]-f\left(\mathbf{x}\right)$

к ${\psi }_i\left(x\right)$ примет вид

$\sum _{j=1}^{\mathcal{M}}{c}_j\cdot \left[\underset{\Omega }{\int }\nabla {\psi }_i\left(\mathbf{x}\right)\cdot \right.\underset{\Omega }{\int }G(\mathbf{x},\mathbf{y})\nabla {\psi }_j\left(\mathbf{y}\right)d\mathbf{y}d\mathbf{x}-$ (5)

$\left.-{\beta }^2\underset{\Omega }{\int }{\psi }_i\left(\mathbf{x}\right)\cdot \underset{\Omega }{\int }G\left(\mathbf{x},\mathbf{y}\right){\psi }_j\left(\mathbf{y}\right)d\mathbf{y}d\mathbf{x}\right]=$

$=\underset{\Omega }{\int }f\left(\mathbf{x}\right)\cdot {\psi }_i\left(\mathbf{x}\right)d\mathbf{x}.$

Эффективность решения (5) существенным образом зависит от выбора ${\mathbf{\psi }}_j$ для $\Omega$. При этом исследование задачи предлагается формировать из функций ${\mathbf{\psi }}_j$, способных в последующем численном решении (5) и аппроксимации (4) обеспечить возможность представить $\tilde{\mathbf{J}}\left(\mathbf{x}\right)$ в виде суммы безвихревых ${\tilde{\mathbf{J}}}^e\left(\mathbf{x}\right)$ и вихревых ${\tilde{\mathbf{J}}}^h\left(\mathbf{x}\right)$ токов.

2. Решение дифракционной задачи с учетом разделения поверхностной плотности тока на вихревую и безвихревую составляющие

Особенность предложенного решения состоит в формировании гладкой аппроксимации $\mathbf{J}\left(\mathbf{x}\right):$

$\tilde{\mathbf{J}}\left(\mathbf{x}\right)={\tilde{\mathbf{J}}}_e^{}\left(\mathbf{x}\right)+{\tilde{\mathbf{J}}}_h^{}\left(\mathbf{x}\right),$

заданной конечными суммами:

$\tilde{\mathbf{J}}\left(\mathbf{x}\right)=\sum _{j=1}^{\stackrel{⌢}{N}}{c}_j^e\cdot {\mathbf{\psi }}_j^e\left(\mathbf{x}\right)+\sum _{\overline{j}=1}^{\overline{N}}{c}_{\overline{j}}^h\cdot {\mathbf{\psi }}_{\overline{j}}^h\left(\mathbf{x}\right),$

где $\stackrel{⌢}{N},$ $\overline{N}$ – количество базисных функций каждого типа.

Согласно [20], для $x\in \Omega$ имеет место двумерный аналог декомпозиции Гельмгольца, который для поля $H$ имеет вид

$\mathbf{H}={\nabla }_{\perp }\cdot {\phi }_{}^D+{{\nabla }^{\text{'}}}_{\perp }\cdot {\phi }_{}^N,$

где

${\nabla }_{\perp }=\left(\frac{\partial }{\partial x_1},\frac{\partial }{\partial x_2}\right),$

${{\nabla }^{\text{'}}}_{\perp }=\left(-\frac{\partial }{\partial x_2},\frac{\partial }{\partial x_1}\right),$

${\phi }_{}^N,$ ${\phi }_{}^D$ – неизвестные скалярные функции, удовлетворяющие условиям Неймана:

${\left.{\nabla }_{\perp }{\phi }_{}^N\cdot \nu \right|}_{\partial \Omega }=0$ (6)

и Дирихле:

${\left.{\phi }_{}^D\right|}_{\partial \Omega }=0$ (7)

соответственно.

Тогда с учетом граничных условий на $\partial \Omega$ решение задачи (5) на $\Omega$ определим:

$\tilde{\mathbf{J}}\left(\mathbf{x}\right)=\sum _{j=1}^{\stackrel{⌢}{N}}{c}_j^e\cdot {\nabla }_{\perp }{\phi }_j^N\left(\mathbf{x}\right)+\sum _{\overline{j}=1}^{\overline{N}}{c}_{\overline{j}}^h\cdot {{\nabla }^{\text{'}}}_{\perp }{\phi }_{\overline{j}}^D\left(\mathbf{x}\right).$ (8)

В качестве ${\phi }_j^N$ и ${\phi }_{\overline{j}}^D$ предлагается применять двумерные модифицированные функции Бесселя первого рода действительного переменного ${\stackrel{⌢}{J}}_j^{}\left(\mathbf{x}\right)$ и многочлены Цернике ${\stackrel{⌢}{Z}}_j^{}\left(x\right),$ удовлетворяющие условиям (6) и (7) соответственно.

Удовлетворяющие граничному условию Неймана (6) функции ${\stackrel{⌢}{J}}_j^{}\left(\mathbf{x}\right)$ зададим в виде

${\stackrel{⌢}{J}}_j^{}\left(x\right)={\stackrel{⌢}{J}}_j^{}(\rho ,\phi )=$

$=\left\{\begin{array}{l}{J}_{}^{\alpha }\left({\lambda }_{(j+1)/2}^{\alpha }\rho \right)\cos \left(k\phi \right),j\mathrm{mod}2=\mathrm{1,}\\ J_{}^{\alpha }\left({\lambda }_{j/2}^{\alpha }\rho \right)\sin \left(k\phi \right),j\mathrm{mod}2=0;\end{array}\right.$ (9)

где

$\rho =\rho \left(x\right)=\sqrt{x_1^2+x_2^2};$

$\phi =\phi \left(x\right)=\arg (x_1+i{x}_2);$

$J_{}^{\alpha }\left(\rho \right)$ – функции Бесселя перового рода порядка $\alpha \ge 1;$ $\lambda =\left\{{\lambda }_i^{\alpha }\right(k\left)\right\}$ – упорядоченное по возрастанию ${\lambda }_i^{}<{\lambda }_{i+1}^{}$ (рис. 2) множество всех нулей

${\left.\frac{d{J}_{}^{\alpha }\left(\rho \right)}{d\rho }\right|}_{\rho ={\lambda }_i^{\alpha }}=0;$

$k$ – порядковый номер нуля функции

${J^{\text{'}}}^{\alpha }\left(\rho \right)=\frac{d{J}_{}^{\alpha }\left(\rho \right)}{d\rho };$

$i\in \left[\mathrm{1,}\stackrel{⌢}{N}/2\right]$ – порядковый номер элемента множества $\lambda$.

Рис. 2. Примеры графиков ${J^{\text{'}}}^{\alpha }\left(\rho \right)$ и расположение первых 12 элементов множества нулей $\lambda =\left\{{\lambda }_i^{\alpha }\right(k\left)\right\}$

Fig. 2. Examples of graphs of ${J^{\text{'}}}^{\alpha }\left(\rho \right)$ and location of the first 12 elements of the set of zeros $\lambda =\left\{{\lambda }_i^{\alpha }\right(k\left)\right\}$

Модификацию ${\stackrel{⌢}{Z}}_{\overline{j}}^{}\left(\mathbf{x}\right)$ функций Цернике $Z_{\overline{j}}^{}\left(\mathbf{x}\right)$ при удовлетворении граничному условию (7) определим как

${\stackrel{⌢}{Z}}_{\overline{j}}^{}\left(\mathbf{x}\right)={\stackrel{⌢}{Z}}_n^m(\rho ,\phi )=\left\{\begin{array}{c}{R}_n^m\left(\rho \right)\cos \left(m\phi \right),m\ge 0;\\ R_n^m\left(\rho \right)\sin \left(-m\phi \right),m<\mathrm{0,}\end{array}\right.$ (10)

$R_n^m\left(\rho \right)=\rho (1-{\rho }^2)P_{\frac{n-m}{2}}^{\left(\mathrm{2,}m\right)}(2{\rho }^2-1),$

где

$P_t^{(\alpha ,\beta )}\left(x\right)=\frac{\Gamma (\alpha +t+1)}{t!\Gamma (\alpha +\beta +t+1)}\times$

$\times \sum _{s=0}^t\left(\begin{array}{c}t\\ s\end{array}\right)\left(\frac{\Gamma (\alpha +\beta +t+s+1)}{\Gamma (\alpha +s+1)}\right){\left(\frac{x-1}{2}\right)}^2$

– многочлены Якоби порядка t; $n\in \mathbb{N};$ $m\in \mathbb{Z}/\left\{0\right\};$ $n\ge \left|m\right|;$ $\frac{n-m}{2}\in \mathbb{Z};$ $\overline{j}=\mathrm{0,5}(n^2-1-n\mathrm{mod}2+\mathrm{sgn}(m\left)\right)+\left|m\right|.$

Примеры графического представления предложенной модификации функции Бесселя первого рода действительного переменного ${\stackrel{⌢}{J}}_j^{}\left(x\right)$ и многочленов Цернике ${\stackrel{⌢}{Z}}_j^{}\left(x\right),$ удовлетворяющие условиям (6) и (7) соответственно, приведены на рис. 3.

Рис. 3. Пример визуализации многочленов ${\stackrel{⌢}{J}}_j^{}\left(x\right)$ (а) и ${\stackrel{⌢}{Z}}_{\overline{j}}^{}\left(x\right)$ (б)

Fig. 3. Example of polynomials visualization ${\stackrel{⌢}{J}}_j^{}\left(x\right)$ (a) and ${\stackrel{⌢}{Z}}_{\overline{j}}^{}\left(x\right)$ (b)

С учетом поставленной задачи дифракции и предложении аппроксимации поверхностной плотности тока в виде (8) зададим (5) виде системы уравнений:

$\sum _{{\overline{j}}^{\text{'}}=1}^{\overline{N}}{c}_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\cdot \left\{\begin{array}{l}\\ \end{array}\right.\underset{\Omega }{\int }\nabla \cdot {\psi }_{\overline{j}}^h\left(\mathbf{x}\right)\underset{\Omega }{\int }G(\mathbf{x}\mathbf{,}\mathbf{y})\nabla \cdot {\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)d\mathbf{y}dx-$ (11)

$-{\beta }^2\underset{\Omega }{\int }{\psi }_{\overline{j}}^h\left(\mathbf{x}\right)\underset{\Omega }{\int }G(\mathbf{x},\mathbf{y}){\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)d\mathbf{y}d\mathbf{x}\left.\begin{array}{l}\\ \end{array}\right\}+$

$\begin{array}{l}+\sum _{j^{\text{'}}=1}^{\stackrel{⌢}{N}}{c}_{j^{\text{'}}}^{e^{\text{'}}}\cdot \left\{\begin{array}{l}\\ \end{array}\right.\underset{\Omega }{\int }\nabla \cdot {\psi }_{\overline{j}}^h\left(x\right)\underset{\Omega }{\int }G(x,y)\nabla \cdot {\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx-\\ -{\beta }^2\underset{\Omega }{\int }{\psi }_{\overline{j}}^h\left(x\right)\underset{\Omega }{\int }G(x,y){\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx\left.\begin{array}{l}\\ \end{array}\right\}=\\ =\underset{\Omega }{\int }f\left(x\right){\psi }_{\overline{j}}^h\left(x\right)dx,\end{array}$

$\sum _{{\overline{j}}^{\text{'}}=1}^{\overline{N}}{c}_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\cdot \left\{\begin{array}{l}\\ \end{array}\right.\underset{\Omega }{\int }\nabla \cdot {\psi }_j^e\left(x\right)\underset{\Omega }{\int }G(x,y)\nabla \cdot {\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)dydx-$ (12)

$-{\beta }^2\underset{\Omega }{\int }{\psi }_j^e\left(x\right)\underset{\Omega }{\int }G(x,y){\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)dydx\left.\begin{array}{l}\\ \end{array}\right\}+$

$\begin{array}{l}+\sum _{j^{\text{'}}=1}^{\stackrel{⌢}{N}}{c}_{j^{\text{'}}}^{e^{\text{'}}}\cdot \left\{\begin{array}{l}\\ \end{array}\right.\underset{\Omega }{\int }\nabla \cdot {\psi }_j^e\left(x\right)\underset{\Omega }{\int }G(x,y)\nabla \cdot {\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx-\\ -{\beta }^2\underset{\Omega }{\int }{\psi }_j^e\left(x\right)\underset{\Omega }{\int }G(x,y){\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx\left.\begin{array}{l}\\ \end{array}\right\}=\\ =\underset{\Omega }{\int }f\left(x\right){\psi }_j^e\left(x\right)dx.\end{array}$

Решение системы уравнений (11), (12) позволит определить поверхностную плотность тока $\stackrel{\mathbf{~}}{\mathbf{J}}\mathbf{\left(}\mathbf{x}\mathbf{\right)}$ с требуемой точностью при обеспечении возможности выделения тока на освещенной ${\stackrel{\mathbf{~}}{\mathbf{J}}}^{осв}\left(\mathbf{x}\right)$ и теневой ${\tilde{\mathbf{J}}}^{тен}\left(\mathbf{x}\mathbf{\right)}$ сторонах диска.

3. Расчет поверхностной плотности тока на освещенной и теневой сторонах Ω

В [16] доказано, что с каждой стороны бесконечно тонкого идеально проводящего плоского экрана известны тангенциальная составляющая магнитного поля $H_{\tau }$ и нормальная составляющая электрического поля $E_n$:

$\underset{x_3\to \pm 0}{\lim }{\mathbf{E}}_n\left(\mathbf{x}\right)=\mp \frac{1}{2}\frac{i}{\omega \epsilon }{\nabla }_{\perp }\mathbf{J}\left(\mathbf{x}\right);$ (13)

$\underset{x_3\to \pm 0}{\lim }{\mathbf{H}}_{\tau }\left(\mathbf{x}\right)=\pm \frac{1}{2}\mathbf{J}\left(\mathbf{x}\right)\times \mathbf{n}.$

Согласно принципу физического эквивалента [23], поля ${\mathbf{H}}^0$ и $\mathbf{H}$ на поверхности идеально проводящего тела можно заменить эквивалентным поверхностным током ${\mathbf{J}}^{eq}$:

${\mathbf{J}}_{}^{eq}=\tilde{\mathbf{n}}\times ({\mathbf{H}}^0+\mathbf{H}),$ (14)

где $\tilde{\mathbf{n}}$ – нормаль к поверхности тела.

С учетом того что поле ${\mathbf{H}}^0$ известно, а $\mathbf{H}={\left.{\mathbf{H}}_{\tau }\left(\mathbf{x}\right)\right|}_{x_3\to \pm 0}$ поддается вычислению для каждой из сторон $\Omega$, подставив (13) в (14) при непрерывности падающего поля ${\mathbf{H}}^0$ в пределе получим:

${\mathbf{J}}_{}^{осв}\left(\mathbf{x}\right)=\underset{x_3\to +0}{\lim }{\mathbf{J}}^{eq}\left(\mathbf{x}\right)=$

$\begin{array}{l}=\underset{x_3\to +0}{\lim }\left[n\times H^0\left(x\right)+n\times H_{\tau }\left(x\right)\right]=\\ =\mathbf{n}\times {\mathbf{H}}^0\left(\mathbf{x}\right)+\mathbf{n}\times \left(\frac{1}{2}\mathbf{J}\left(\mathbf{x}\right)\times \mathbf{n}\right)=\\ =\mathbf{n}\times {\mathbf{H}}^0\left(\mathbf{x}\right)+\frac{1}{2}\mathbf{J}\left(\mathbf{x}\right)\approx \\ \approx \mathbf{n}\times {\mathbf{H}}^0\left(\mathbf{x}\right)+\frac{1}{2}\stackrel{\mathbf{~}}{\mathbf{J}}\left(\mathbf{x}\right)={\tilde{\mathbf{J}}}_{}^{осв}\left(\mathbf{x}\right);\end{array}$

${\mathbf{J}}_{}^{тен}\left(\mathbf{x}\right)=\underset{x_3\to -0}{\lim }{\mathbf{J}}^{eq}\left(\mathbf{x}\right)=$

$\begin{array}{l}=\underset{x_3\to -0}{\lim }\left[-\mathbf{n}\times {\mathbf{H}}^0\left(\mathbf{x}\right)-\mathbf{n}\times {\mathbf{H}}_{\tau }\left(\mathbf{x}\right)\right]=\\ =-\mathbf{n}\times {\mathbf{H}}^0\left(\mathbf{x}\right)-\mathbf{n}\times \left(-\frac{1}{2}\mathbf{J}\left(\mathbf{x}\right)\times \mathbf{n}\right)=\\ =-\mathbf{n}\times {\mathbf{H}}^0\left(\mathbf{x}\right)+\frac{1}{2}\mathbf{J}\left(\mathbf{x}\right)\approx \\ \approx -\mathbf{n}\times {\mathbf{H}}^0\left(\mathbf{x}\right)+\frac{1}{2}\stackrel{\mathbf{~}}{\mathbf{J}}\left(\mathbf{x}\right)={\tilde{\mathbf{J}}}_{}^{тен}\left(x\right).\end{array}$

Таким образом, предложенный способ позволяет преодолеть ограничения (5) при асимптотическом выделении тока с каждой из сторон $\Omega$. Для сформированных представлений составим численную схему решения задачи (11), (12) при определении $J_{}^{осв}\left(x\right),$ $J_{}^{тен}\left(x\right)$ и выделении алгоритмических особенностей.

4. Особенности алгоритмической реализации

Точность решения системы (11), (12) во многом зависит от численного вычисления одиночных и двойных интегралов. С этой целью представим область интегрирования $\Omega$ многоугольником $S\in {ℝ}^2,$ $({х}_3=0),$ состоящим из

$S=\underset{s=1}{\overset{M}{\cup }}{S}_s$

треугольных областей таких, что

${\overline{S}}_m\cap {\overline{S}}_{m^{\text{'}}}=\varnothing$ $(m,m^{\text{'}}\in \{\overline{\mathrm{1,}M}\},m\ne m^{\text{'}}).$

Разбиение $\Omega$ на треугольные элементы произведем триангуляцией Делоне (рис. 4) [6].

Рис. 4. Вариант представления диска $\Omega$ многоугольником S. Точками обозначены узлы численного интегрирования

Fig. 4. Variant of the disk representation $\Omega$ by a polygon S. The dots indicate the nodes of numerical integration

Окончательная система линейных алгебраических уравнений (СЛАУ) относительно неизвестных коэффициентов ${с}_{}^{e^{\text{'}}},$ ${с}_{}^{h^{\text{'}}}$ примет вид

$\sum _{s=1}^M\left[\sum _{s^{\text{'}}=1}^M\left[\sum _{{\overline{j}}^{\text{'}}=1}^{\overline{N}}{c}_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\cdot \right.\right.\left\{\begin{array}{l}\\ \end{array}\right.\underset{S_s}{\int }\nabla \cdot {\psi }_{\overline{j}}^h\left(\right)\underset{S_{s^{\text{'}}}}{\int }G(,)\nabla \cdot$

$\begin{array}{l}{\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)dydx-{\beta }^2\underset{S_s}{\int }{\psi }_{\overline{j}}^h\left(x\right)\underset{S_{s^{\text{'}}}}{\int }G(x,y){\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)dydx\left.\begin{array}{l}\\ \end{array}\right\}+\\ +\sum _{j^{\text{'}}=1}^N{c}_{j^{\text{'}}}^{e^{\text{'}}}\cdot \left\{\begin{array}{l}\\ \end{array}\right.\underset{S_s}{\int }\nabla \cdot {\psi }_{\overline{j}}^h\left(x\right)\underset{S_{s^{\text{'}}}}{\int }G(x,y)\nabla \cdot \\ \cdot {\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx-{\beta }^2\underset{S_s}{\int }{\psi }_{\overline{j}}^h\left(x\right)\underset{S_{s^{\text{'}}}}{\int }G(x,y){\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx\left.\left.\left.\begin{array}{l}\\ \end{array}\right\}\right]\right]=\\ =\sum _{s=1}^M\underset{S_s}{\int }f\left(x\right){\psi }_{\overline{j}}^h\left(x\right)dx,\end{array}$

$\sum _{s=1}^M\left[\sum _{s^{\text{'}}=1}^M\left[\sum _{{\overline{j}}^{\text{'}}=1}^{\overline{N}}{c}_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\cdot \right.\right.\left\{\begin{array}{l}\\ \end{array}\right.\underset{S_s}{\int }\nabla \cdot {\psi }_j^e\left(x\right)\underset{S_{s^{\text{'}}}}{\int }G(x,y)\nabla \cdot$

$\begin{array}{l}\cdot {\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)dydx-{\beta }^2\underset{S_s}{\int }{\psi }_j^e\left(x\right)\underset{S_{s^{\text{'}}}}{\int }G(x,y){\psi }_{{\overline{j}}^{\text{'}}}^{h^{\text{'}}}\left(y\right)dydx\left.\begin{array}{l}\\ \end{array}\right\}+\\ +\sum _{j^{\text{'}}=1}^N{c}_{j^{\text{'}}}^{e^{\text{'}}}\cdot \left\{\begin{array}{l}\\ \end{array}\right.\underset{S_s}{\int }\nabla \cdot {\psi }_j^e\left(x\right)\underset{S_{s^{\text{'}}}}{\int }G(x,y)\nabla \cdot \\ \cdot {\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx-{\beta }^2\underset{S_s}{\int }{\psi }_j^e\left(x\right)\underset{S_{s^{\text{'}}}}{\int }G(x,y){\psi }_{j^{\text{'}}}^{e^{\text{'}}}\left(y\right)dydx\left.\left.\left.\begin{array}{l}\\ \end{array}\right\}\right]\right]=\\ =\sum _{s=1}^M\underset{S_s}{\int }f\left(x\right){\psi }_j^e\left(x\right)dx.\end{array}$

Ввиду того что $\Omega \subset {ℝ}^2=$$(x_3=0)\subset {ℝ}^3,$ $\nabla \cdot {\psi }_j^h\left(\right)=0$ и $\nabla \cdot {\psi }_{j^{\text{'}}}^{h^{\text{'}}}\left(x\right)=0.$

В качестве узлов численного интегрирования внутри каждого $S_s$ используются корни многочленов Дубинера – Курнвиндера [21]:

$Q_{n,k}(x_1,x_2)=P_{n-k}^{(2k+\mathrm{1,0})}\left(x_1\right){(1-x_1)}^k{P}_k^{\left(\mathrm{0,0}\right)}\left(\frac{x_2}{1-x_1}\right),$ (19)

где n – максимальная степень многочлена Якоби $P_n^{(\alpha ,\beta )}\left(x\right);$ $0\le k\le n.$

Координаты узлов интегрирования определим по правилу (рис. 5):

${\mathbf{x}}_{a,b}=\left\{\left(x_1,x_2(1-x_1)\right):x_1={\left(X_1^{a,b}\right)}_b,\right.$

$\begin{array}{l}{x}_2={\left(X_2^{a,b}\right)}_{\left(n-b-a\right)};a\in \overline{\left[\mathrm{0,}n-1\right]};\\ \left.b\in \overline{\left[\mathrm{0,}a\right]};a,b,n\in {\mathbb{N}}_0\right\};\\ X_1^{a,b}=\left\{x_s:J_{n-(a-b)}^{\left(2(a-b)+\mathrm{1,0}\right)}\left(x_s\right)=0;\right.\\ s\in \overline{\left(\mathrm{0,}n-(a-b)\right)};\\ x_0<x_1<\mathrm{...}<x_{n-(a-b)};x_s\in ℝ\left.\begin{array}{l}\\ \end{array}\right\};\\ X_2^{a,b}=\left\{x_p:J_{n-b}^{\left(\mathrm{0,0}\right)}\left(x_p\right)=0;p\in \overline{(\mathrm{0,}n-b)};\right.\\ x_0<x_1<\mathrm{...}<x_{n-b};x_p\in ℝ\left.\begin{array}{l}\\ \end{array}\right\}.\end{array}$

Рис. 5. Корни многочлена Дубинера – Курнвиндера при $n=4.$ Нумерация корней $x_{a,b}$ отображена в виде $(row,col),$ где row – номер строки, col – номер столбца: a – на единичном треугольнике; б – на произвольном треугольнике

Fig. 5. Roots of the Dubiner–Koornwinder polynomial at $n=4.$ The numbering of the roots $x_{a,b}$ is shown in the form $(row,col),$ where row is the row number, col is the column number: a – on a unit triangle; b – on an arbitrary triangle

В ходе вычисления (11) и (12) возникает ситуация, когда $s=s\text{'}$ и $x=y$ при $\left|x-y\right|\to 0$ и $G(x,y)\to \infty .$ Для исключения сингулярности в $G(x,y)$ при $s=s\text{'}$ интегрирование производится в полярной системе координат [22]. Обозначим $S_s=S_{s\text{'}}=T_s$ тогда в интегралах вида

$G_{j_1,j_2}^1=\underset{T_s}{\int }\nabla \cdot {\psi }_{j_1}^{}\left(x\right)\underset{T_s}{\int }\frac{e^{i\beta \left|x-y\right|}}{\left|x-y\right|}\nabla \cdot {\psi }_{j_2}^{}\left(y\right)dydx;$ (20)

$G_{j_1,j_2}^2=\underset{T_s}{\int }{\psi }_{j_1}^{}\left(x\right)\underset{T_s}{\int }\frac{e^{i\beta \left|x-y\right|}}{\left|x-y\right|}{\psi }_{j_2}^{}\left(y\right)dydx,$

(где $j_1, j_2$ – порядковые номера базисных функций) преобразуем по $dy$ в локальную относительно $T_s$ полярную систему координат $\text{}y(\phi ,\rho )$ с центром в $x$. В таком представлении (20) преобразуется:

$G_{j_1,j_2}^1=\underset{T_s}{\int }\nabla \cdot {\psi }_{j_1}^{}\left(x\right)\cdot \sum _{q=0}^2\underset{{\phi }_q}{\overset{{\phi }_{q1}}{\int }}\underset{0}{\overset{r_q\left(\phi \right)}{\int }}{e}^{i\beta \rho }\nabla \cdot$

$\begin{array}{l}\cdot {\psi }_{j_2}^{}\left(y(\phi ,\rho )\right)d\rho d\phi dx;\\ G_{j_1,j_2}^2=\underset{T_s}{\int }{\psi }_{j_1}^{}\left(x\right)\cdot \sum _{q=0}^2\underset{{\phi }_q}{\overset{{\phi }_{q1}}{\int }}\underset{0}{\overset{r_q\left(\phi \right)}{\int }}{e}^{i\beta \rho }{\psi }_{j_2}^{}\left(y(\phi ,\rho )\right)d\rho d\phi dx,\end{array}$

где $q$ – номер вершины $T_s$;

${\phi }_m=\text{arctg}\frac{x_2^q-x_2}{x_1^q-x_1};$

$y(\rho ,\phi )=(\rho \cos \phi +x_1,\rho \sin \phi +x_2);$

$r_q\left(\phi \right)=\frac{{R^{\text{'}}}_q}{\cos (\phi -{\varphi }_q)};$ ${\varphi }_q=arctg\frac{{x^{\text{'}}}_2^q-x_2}{{x^{\text{'}}}_1^q-x_1};$

$x{\text{'}}^q=\left(\frac{b_q{d}_q-a_q{c}_q}{{\left(a_q\right)}^2+{\left(b_q\right)}^2},\frac{-a_q{d}_q-b_q{c}_q}{{\left(a_q\right)}^2+{\left(b_q\right)}^2}\right)$

– координаты пересечения $△x^q{x}^{q_1}x$ со стороной $\left(x^q{x}^{q_1}\right);$

$a_q=x_2^{q_1}-x_2^q;$ $b_q=x_1^{q_1}-x_1^q;$

$c_q=x_1^{q_1}{x}_2^q-x_2^{q_1}{x}_1^q;$ $d_q=b_q{x}_1-a_q{x}_2;$

$q_1=(q+1)\mathrm{mod}3;$ ${R^{\text{'}}}_q=\left|{x^{\text{'}}}^l-x\right|;$ $x^q=(x_1^q,x_2^q)$

– координаты вершин $T_s$

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

5. Результаты верификации сформированных решений

Для наглядной демонстрации предпочтительности сформированного решения выполним серию вычислительных экспериментов, которые предполагают получение плотности тока на освещенной и теневой сторонах диска для различных порядков аппроксимации полиномов и размеров диска, заданного соотношением $D/\lambda$, $\lambda$ – длина волны.

Эталонная модель, используемая для верификации полученных результатов, разработана в САПР ANSYS HFSS. Она имеет форму идеально проводящего цилиндра высотой 0,0001 м и радиусом D. В качестве источника ЭМП задана плоская монохроматическая волна с частотой 1,5 ГГц и фронтом, параллельным диску. Установки программы: режим моделирования: HFSS IE Solver; Maximum Number of Passes: 20; Maximum Residual Error: 0,0002. Сравнение полученных результатов для различных соотношений $D/\lambda =\left\{\mathrm{1,2,3,4}\right\}$ $(\lambda =\text{const})$ приведены на рис. 6, 7.

Рис. 6. Сравнение $\left|{\tilde{J}}_{}^{осв}\left(x\right)\right|$ и $\left|{\tilde{J}}_{}^{тен}\left(x\right)\right|$ в САПР Ansys HFSS (а) и разработанного решениея (б) при различных соотношениях $D/\lambda$

Fig. 6. Comparison of $\left|{\tilde{J}}_{}^{осв}\left(x\right)\right|$ and $\left|{\tilde{J}}_{}^{тен}\left(x\right)\right|$ in Ansys HFSS CAD (a) and the developed solution (b) at different $D/\lambda$

Рис. 7. Совмещение графиков $\left|{\tilde{J}}_{}^{осв}\left(x\right)\right|$ и $\left|{\tilde{J}}_{}^{тен}\left(x\right)\right|$ для различных соотношений $D/\lambda$. Оранжевый – САПР Ansoft HFSS, цветной – разработанное решение

Fig. 7. Matching plots $\left|{\tilde{J}}_{}^{осв}\left(x\right)\right|$ and $\left|{\tilde{J}}_{}^{тен}\left(x\right)\right|$ for different $D/\lambda .$ Orange – Ansoft HFSS CAD, colored – developed solution

Оценка апостериорной сходимости оценивалась относительно модуля

$\left|\tilde{\mathbf{J}}\left(\mathbf{x}\right)\right|=\sqrt{{\tilde{\mathbf{J}}}_1\left(\mathbf{x}\right)\cdot \overline{{\tilde{J}}_1\left(\mathbf{x}\right)}+{\tilde{J}}_2\left(\mathbf{x}\right)\cdot \overline{{\tilde{J}}_2\left(\mathbf{x}\right)}+{\tilde{J}}_3\left(\mathbf{x}\right)\cdot \overline{{\tilde{J}}_3\left(\mathbf{x}\right)}}$

по нормам:

${||\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|||}_{L_2}=\sqrt{\underset{\Omega \backslash \delta \Omega }{\int }{\left(\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|\right)}^{2}dx},$

${||\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|||}_C=\underset{x\in \Omega \backslash \delta \Omega }{\max }\left|\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|\right|,$

где $\delta \Omega :=\left\{x:\left|x-y\right|<\delta ,y\in \partial \Omega \right\}$ при $\delta >0.$

В качестве эталона $\overline{\mathbf{J}}\left(\mathbf{x}\right)$ принято решение, сформированное в САПР Ansoft HFSS.

Из представленных на рис. 8 графических зависимостей следует, что не все базисные функции вносят одинаковый вклад в решение тестовой задачи. В этой связи на графиках наблюдаются «скачки» (рис. 8). Отдельные увеличения ошибки (особенно для малых $D/\lambda$ по норме $C$ при росте числа базисных функций связано с ошибочным нахождением в САПР Ansoft HFSS поверхностной плотности тока $\overline{\mathbf{J}}\left(\mathbf{x}\right)$ вблизи $\partial \Omega$ (рис. 7). Также для улучшения сходимости при малых соотношениях $D/\text{λ}$ требуется дополнительная регуляризация СЛАУ. В реализованном алгоритме ее решение выполнено прямым методом.

Рис. 8. Зависимость ${||\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|||}_{L_2}$ (а) и ${||\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|||}_C$ (б) от числа базисных функций

Fig. 8. Dependence of ${||\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|||}_{L_2}$ (a) and ${||\left|\tilde{J}\left(x\right)\right|-\left|\overline{J}\left(x\right)\right|||}_C$ (b) on the number of basis functions

В целом применение составленного алгоритмического решения при исследовании задачи дифракции на $\Omega$ и выделении $J_{}^{осв}\left(x\right),$ $J_{}^{тен}\left(x\right)$ обеспечивает экспоненциальную сходимость по норме $L_2$ и полиномиальной по норме в $C$.

Заключение

Полученные результаты позволяют сделать вывод о предпочтительном применении предложенной модификации функций Бесселя и Цернике при решении задачи дифракции на $\Omega$. Основное достоинство предлагаемой схемы состоит в разделении суммарной поверхностной плотности тока $\tilde{\mathbf{J}}\left(\mathbf{x}\right)$ на две составляющие: на освещенной ${\tilde{\mathbf{J}}}^{осв}\left(\mathbf{x}\right)$ и теневой ${\tilde{\mathbf{J}}}^{тен}\left(\mathbf{x}\right)$ сторонах $\Omega$. Применение векторных базисных функций ${\psi }_j^e\left(x\right)$ и ${\psi }_{\overline{j}}^h\left(x\right),$ основанных на многочленах (9) и (10), удовлетворяющих граничным условиям (6) и (7), позволяет учесть влияние ${\mathbf{E}}^{\mathbf{0}}$ на ${\tilde{\mathbf{J}}}_e^{}\left(\mathbf{x}\right)$ и ${\tilde{\mathbf{J}}}_h^{}\left(\mathbf{x}\right)$. В свою очередь, поле ${\mathbf{H}}^0$ формирует дополнительные поверхностные токи равной амплитуды, но разного направления с каждой из сторон $\Omega$, что в совокупности с принципом эквивалентности позволяет разделить $\tilde{\mathbf{J}}\left(\mathbf{x}\right)$ на ${\tilde{\mathbf{J}}}_{}^{осв}\left(\mathbf{x}\right)$ и ${\tilde{\mathbf{J}}}_{}^{тен}\left(\mathbf{x}\right)$

Следует уточнить, что применение стандартных норм $L_2$ и $C$ при $x\in \Omega$ для оценки сходимости в рассматриваемой задаче является некорректным, что обуславливается граничными условиями $\mathbf{J}\left(x\right)$ вблизи $\partial \Omega$ при возникающей сингулярности [16]. В этой связи в нормах $L_2$ и $C$ для апостериорной оценки сходимости результатов численного решения выбрана модификация при $x\in \Omega \backslash \delta \Omega .$ В последующих исследованиях для априорной оценки сходимости предполагается выбирать пространства Соболева [24].

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

作者简介

Dmitry Ketukh

编辑信件的主要联系方式.
Email: kdk_hcc@mail.ru

member of the Academy of the Federal Guard Service of the Russian Federation. Author of 10 scientific publications.

Research interests: radio engineering, antenna systems, computational electrodynamics, digital signal processing

参考

补充文件

附件文件

动作

1. JATS XML

下载

2. Fig. 1. Geometric representation of the problem of diffraction of a plane wave on an infinitely thin perfectly conducting circular disk

下载 (52KB)

索引源数据

3. Fig. 2. Examples of graphs of and location of the first 12 elements of the set of zeros

下载 (506KB)

索引源数据

4. Fig. 3. Example of polynomials visualization (a) and (b)

下载 (2MB)

索引源数据

5. Fig. 4. Variant of the disk representation by a polygon S. The dots indicate the nodes of numerical integration

下载 (789KB)

索引源数据

6. Fig. 5. Roots of the Dubiner–Koornwinder polynomial at The numbering of the roots is shown in the form where row is the row number, col is the column number: a – on a unit triangle; b – on an arbitrary triangle

下载 (205KB)

索引源数据

7. Fig. 6. Comparison of and in Ansys HFSS CAD (a) and the developed solution (b) at different

下载 (1MB)

索引源数据

8. Fig. 7. Matching plots and for different Orange – Ansoft HFSS CAD, colored – developed solution

下载 (3MB)

索引源数据

9. Fig. 8. Dependence of (a) and (b) on the number of basis functions

下载 (356KB)

索引源数据