Влияние поглощающего слоя на численное решение параболического уравнения Текст научной статьи по специальности «Физика»

УДК: 621.371: 519.6

Ю.П. Акулиничев, П.В. Абрамов, И.Н. Ваулин

Влияние поглощающего слоя

на численное решение параболического уравнения

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

Последнее десятилетие ознаменовалось разработкой ряда методов непосредственного численного решения волновых уравнений, описывающих распространение радиоволн в неоднородной тропосфере над неровной поверхностью Земли [1]. Исходными данными для расчета ожидаемой напряженности поля в прямоугольной области * дальность — высота над поверхностью земли», расположенной в заданном направлении от источника, являются характеристики источника излучения (мощность, форма диаграммы направленности антенны, её высота над поверхностью земли и т.п.), высотные профили индекса преломления воздуха во всех точках трассы, рельеф и электрические свойства подстилающей поверхности. Одним из таких уравнений является двумерное параболическое уравнение (ПУ) для комплексной огибающей U(X,Z) напряженности E(X,Z) = U(X,Z)exp(ikX) электрического или магнитного поля [1]

ш dU^Z) + d2U(XZ) + fe2 г 10_6iV(Xj 2) + 2/ал ЩХг Z) = 0t (1)

ал dz L J

где X — координата на оси, направленной вдоль трассы по «выпрямленной» поверхности Земли (горизонтальная дальность); Z — высота точки наблюдения над поверхностью Земли; N(X,Z) — поле индекса преломления тропосферы; аэ = 8,5-106 м — эквивалентный радиус Земли для стандартной радиоатмосферы.

На практике для численного решения (1) в области О <Х < D, 0 < Z < Н вводится расчетная сетка с ячейками AX-AZ, где АХ = D/M и ДZ = H/N — величины шагов квантования непрерывных функций U(X,Z) и N(X,Z) по дальности и высоте соответственно. Задав во всех узлах сетки на линии X = 0, что соответствует, например, апертуре передающей антенны, значения функции ¡У„(0) = U(0,nAZ) и, используя один и тот же алгоритм, например метод преобразования Фурье (ПФ) [1], последовательно, шаг за шагом, находят значения этой функции во всех узлах сетки Un( 1) = U(AX,nAZ), Un(2) = U(2AX,nAZ), ... при X = АХ, X = 2АХ и так далее до X = D.

При численном решении ПУ методом ПФ на каждом шаге по дальности производятся вычисления вида

U(X + АХ, Z) = w(Z) exp [ikАХ ■ (N(X, Z) - 1<Г6 + Z / 2аэ)] F'1 {ВД^ {*У(Х, Z)}} ,

где F{U(X,Z)} = S(ß) = k J U(X,Z)exv(-ik&Z)dZ - Фурье-образ; #(ß) = exp(-0,5 ■ ¿MXß2) —

частотная характеристика, описывающая распространение плоской волны на расстояние АХ под углом ß к оси ОХ; w(Z) — характеристика пропускания искусственного поглощающего слоя, который вводится на каждом шаге по дальности для имитации нелокального граничного условия на верхней границе области расчета. Эта верхняя граница должна представлять на конечной дистанции условия излучения Зоммерфельда [6] — она должна быть абсолютно прозрачной, позволяющей всей энергии, приходящей снизу к границе, уходить в бесконечность.

Поглощающий слой задаем в виде модифицированного окна Хеннинга (рис. 1):

ш(п) = 1 -Ас+Ас- сое

71(п-АГс)

для ЛГС < п < N , (2)

где Ас — амплитуда функции поглощения (для типичного окна Хеннинга Ас = 0,5). Область 0 < п < Ыс, в которой шп = 1, считается рабочей.

Конкурирующей методу ПФ является схема Кранка-Николсон (К.-Н.) [1, 2], суть которой заключается в следующем. В (1) производные заменяются соответствующими конечными разностями, в результате чего составляется система уравнений вида

g*(Un+1(m +1) + ип_х{т +1)) + (1 - 2ё*)ип{т +1) = й(ип+1(т) + и^т)) + (1 - 2ё)ип(т), (3) где g = АД;е/(2гяА^2) — весовой коэффициент схемы К.-Н.; * — знак комплексного сопряжения.

Зная значения поля на т.-м шаге по дальности и решая систему уравнений (3) методом прогонки [2], находят все значения поля на следующем, то есть (т+1)-м, шаге.

Для метода Фурье поглощающий слой пока является единственной возможностью задания верхнего граничного условия, а для схемы К.-Н. существуют другие методы [3-5].

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

Влияние Земли учитываться не будет, поэтому граничное условие вводится одинаково сверху и снизу. Ширину поглощающего слоя будем считать фиксированной и равной четвертой части рабочей области (см. рис. 1).

Полоса расчета

Рис. 1. Распределение коэффициентов поглощения в полосе расчета при амплитудах функции поглощения: 1 — Ас = 0,5; 2 — Д, = 0,3; 3 — Д. = 0,1

Поглощающий слой для двух методов вводится одинаково: умножаем поле в узлах сетки на коэффициенты поглощения (2). Погрешности оценим, сравнивая результаты расчета поля Со^ в безграничной среде {И —» идеальные условия) с полем ву, рассчитанным при конечных значениях N с использованием поглощающего слоя. Для этого будем вычислять интегральную квадратическую погрешность (ИКП) (детерминированный аналог среднеквадрати-ческой ошибки) ст (4), скалярное произведение г (детерминированный аналог коэффициента корреляции) (5), также удобно пользоваться величиной ог, полученной из г(6). Чем меньше о и точнее проведен расчет, тем лучше поглощающий слой:

а =

N-1

1=0_

ЛГ-1

2 N

7=0

(4)

г = -

N-1

X Ofio] /=о_

Г*"*, .И

х ЬГ х N

V i=o > U=o

(5)

(6)

Сначала приведем результаты, полученные для метода ПФ. В расчете использовались следующие величины параметров: Х= 0,1 м, AZ = 1 м. Для тестового расчета (N —» взяли 131 072 = 217 точек по вертикали. Из теоремы отсчетов для спектра получили максимальную дальность, на которой начинает проявляться наличие границ. Она составляет 1 310 720 м. Чтобы надежно исключить влияние граничных условий, расчетная дальность была взята много меньше — 150 км. Результаты расчетов показаны на рис. 2.

1.5

Нориирсв энная i амплитуда поля

0.5

О

\ 1 1 J 1

1 1 1 1

О 2000 4000 6000 8000 1 10 Ширина полосы

Рис. 2. Тестовый расчет для Фурье-метода (дальность 150 км)

Ширина четной функции Грина (поле, возбужденное элементарным источником) составляет примерно 15000 точек, что соответствует в нашем случае 15 км. Теперь уменьшим полосу расчета до N = 4096 и введем поглощающий слой снизу и сверху, шириной в четверть полосы расчета (см. рис. 1). О степени совпадения графиков будем судить по величине ст.

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

1.5

« «

« ч

к о **

со £

к £

О s

К I

0.5 "

1 ........1 1 " "Г 2 ...................................."I.................................... '

m \ 3 1

\ 1 ...... 1 L V

500

1000 1500 Полоса расчета

2000

2500 3000

N

Рис. 3. Тестовый расчет функции Грина и расчет с поглощающим слоем на дальности 150 км при амплитуде функции поглощения 0,04: 1 — расчет с поглощающим слоем; 2 — тестовый расчет; 3 — характеристика поглощающего слоя

На рис. 3 показаны результаты расчета модуля напряженности поля с использованием поглощающего слоя при N = 4096 и АХ = 200 м, тестового расчета и функция затухания в слое. При этом |г| = 0,9998, а сгг = 0,02. Это связано и с тем, что фазы функций в полосе расчета практически не отличаются. На рис. 4 показана разность фаз для идеальной и расчетной функций.

ч -1-I-1-1-1-

0 500 1000 1500 2000 2500 3000

Рис. 4. Разность фаз идеальной и расчетной функций: 1 — в рабочей области; 2 — за пределами рабочей области

Практически аналогичные оценки были сделаны для схемы К.-Н. Количество уравнений в системе равно количеству точек по вертикали без учета двух крайних точек, задающих верхнее и нижнее граничные условия.

В расчете использовались следующие параметры: А. = 0,1 м; АЯ =2м; АХ = 8 м; количество шагов 7000. Для тестового расчета (идеальный случай) взяли 1000 точек по вертикали. В этом случае можно считать, что граничные условия не повлияют на решение (рис. 5), так как распределение поля еще не достигает границ.

амплитуда

поля 0.1

J

400 600 800 1000 ширина полосы

Рис. 5. Тестовый расчет для схемы К.-Н. (дальность 56 км)

Затем количество точек уменьшили до 200 и ввели поглощающий слой. Результаты вычислений показаны на рис. 6.

Для этого случая |г| в полосе расчета составляет 0,995 (стг= 0,1). Причем, как видно из рисунка, отличие проявляется только по краям, в центре полосы графики хорошо совпадают.

При амплитуде поглощающего слоя А,. = 0,05 получаем |г| = 0,999996 (стг = 0,003). Если амплитуду уменьшать далее, то корреляция уменьшается. Фаза в этом случае очень хорошо совпадает во всей полосе.

Отметим, что при значительном уменьшении полосы расчета и соответственно объема вычислений результаты очень незначительно отличаются от тестовых расчетов (рис. 7, 8). То есть слой весьма точно имитирует отсутствующее в расчетах бесконечное пространство, куда уходят радиоволны.

Рис. 6. Тестовый расчет при ширине 1000 точек и расчет с поглощающим слоем при ширине полосы расчета 200 точек (Д. = 0,5): 1 — расчет с поглощающим слоем; 2 — тестовый расчет; 3 — характеристика поглощающего слоя, доходит до единицы (см. рис. 1)

250

Рис. 7. Модуль поля для тестового расчета и расчета с поглощающим слоем (амплитуда функции поглощения 0,05)

Разность фаз 0

Рис. 8. Разность фаз полей точечных источников для тестового расчета и расчета с поглощающим слоем (амплитуда функции поглощения 0,05)

Следует также сказать об амплитуде функции поглощения. Вначале предполагалось, что ее величина на каждом шаге будет 0,5; но, как показали результаты, чем больше шагов выполняется, тем меньше должна быть амплитуда поглощающего слоя на каждом шаге. Это

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

В табл. 1-3 приведены значения ИКП, допускаемой при расчете напряженности поля методом ПФ с использованием поглощающего слоя (25 % от полосы расчета) при различных значениях параметров вычислительной схемы.

Интересна особенность табл. 1: когда уменьшается полоса расчета, то увеличивается |г| (хотя предположительно должно быть наоборот). Причем, как показал дополнительный расчет, это не связано с неоптимальностью амплитуды функции поглощения. Это «эффект большой дальности», который был отмечен в [8]. Это влияние отсутствия выбросов на краях функции Грина, то есть, уменьшая полосу расчета, прячем нежелательные выбросы в поглощающем слое.

Уменьшение точности расчета при уменьшении шага дальности (табл. 3) — результат влияния поглощающего слоя; видно, что для АХ — 400 м значение Д, = 0,04 уже не является оптимальным. Как показал дополнительный расчет, оптимально подобранная амплитуда для каждого значения АХ приводит к тому, что результаты будут отличаться незначительно. Все эти данные показывают, что ошибка, вызванная неидеальностью поглощающего слоя, если он занимает четверть полосы расчета по высоте, искажает результаты расчетов напряженности поля не более чем на 0,1-0,3 дБ на расстояниях до 300 км.

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

Таблица 1

Зависимость |г| и аг от ширины полосы расчета N при фиксированной дальности

N 1024 2048 4096

0,04 0,04 0,04/0,05/0,02

И 0,999929 0,999828 0,999704/0,999657/0,999663

стг, дБ 6,2-Ю"4 1,5-10~3 2,6-10~3/3-10~3/3-10_3

Таблица 2

Зависимость |г| и стг от дальности при фиксированной полосе расчета N = 2048, шаге по дальности 200 м и амплитуде функции поглощения 0,04 / 0,05 / 0,02

D, м 50 150 300

Л 0,04 / 0,05 / 0,02 0,04 0,04

И 0,999765 / 0,999727 / 0,999538 0,999828 0,999937

» td 2-Ю"3 / 2,4-Ю"3/ 4-10"3 1,5-Ю-3 5,5-Ю"4

Даже при достаточно экономном задании значений параметров возможно обеспечить суммарную ИКП вычислений в пределах 0,5...1 дБ, что меньше величин ошибок, обусловленных другими причинами, главной из которых является недостаточно детальное знание структуры высотной зависимости индекса преломления на всех участках трассы.

Таблица 3

Зависимость |г| и ог от шага по дальности при фиксированной дальности 150 км, полосе расчета N — 2048 и амплитуде функции поглощения 0,04 / 0,08

АХ, м 50 100 200 400

К 0,04 0,04 0,04 0,04 / 0,08

И 0,999275 0,999598 0,999828 0,996336/0,999825

ог, дБ 6,ЗЮ~3 3,5-10~8 1,5'10~3 3,2-10~2 / 1,4-10~3

Заключение

До настоящего времени не было количественных данных о влиянии граничного условия в виде поглощающего слоя на точность численного решения параболического уравнения. В данной работе подобраны значения параметров функции поглощения, приведена методика их нахождения. Проведена оценка точности расчетов при использовании поглощающего слоя толщиною 1/4 от ширины области вычислений. Было установлено, что уменьшение вычислительных затрат при использовании поглощающего слоя с определенными параметрами может проводиться практически без потерь точности. При этом можно использовать обычное окно Хеннинга, хотя при достаточно большом количестве шагов по дальности (> 1500) и достаточно малом количестве узлов сетки по высоте (< 512) следует уменьшить амплитуду поглощающего слоя Ас примерно в 10 раз. При этом интегральная квадратическая погрешность, обусловленная несовершенством поглощающего слоя, не должна превышать 0,1 дБ.

Литература

1. Levy M. Parabolic equation methods for electromagnetic wave propagation / M. Levy. -London : The IEE, 2000. - 336 p.

2. Самарский A.A. Введение в численные методы / А.А. Самарский. - СПб. : Лань, 2005. -288 с.

3. Arnold A. Discrete Transparent Boundary Conditions for Wide Angle Parabolic Equations in Underwater Acoustics / A. Arnold, M. Ehrhardt // Journal of Computational Physics. - 1998. -Vol. 145, № 2. - P. 611-638.

4. Ehrhardt. M. Discrete transparent boundary conditions for the Schrôdinger equation: fast calculations, approximation and stability / M. Ehrhardt, A. Arnold, I. Sofronov // Communs Math. Sci. - 2003. - Vol. 1. - P. 501-556.

5. Ehrhardt M. Solutions to the Discrete Airy Equation: Application to Parabolic Equation Calculations / M. Ehrhardt, R.E. Mickens // J. Comput. Appl. Math. - 2004. - Vol. 172, Issue 1. -P. 183-206.

6. Нефёдов Е.И. Дифракция электромагнитных волн на диэлектрических структурах / Е.И. Нефёдов. - М. : Наука, 1979. - 272 с.

7.Tirkas Р.А. Higher Order Absorbing Boundary Conditions for the Finite-Difference Time-Domain Method / P.A. Tirkas, C.A. Balanis, R.A. Renaut // IEEE Transactions on Antennas and Propagation. - 1992, Oct. - Vol. 40, Ns 10. - P. 1215-1222.

8. Акулиничев Ю.П. Предельная точность схемы Кранка - Никол сон при численном решении параболического волнового уравнения / Ю.П. Акулиничев, М.Е. Ровкин, И.Н. Ваулин // Распространение радиоволн: сб. докл. XXI всерос. науч. конф. (Йошкар-Ола, 25-27 мая 2005 г.). В 2 т. - Йошкар-Ола : МарГТУ, 2005. - Т. 2. - С. 272-276.

Акулиничев Юрий Павлович

Д-р. техн. наук, профессор, профессор каф. радиотехнических систем ТУСУРа

Тел.: (3822) 41 34 67

Эл. почта: ayp@rts.tomsk.ru

Абрамов Павел Викторович

Бывший младший научный сотрудник НИИ радиотехнических систем при ТУСУРе Ваулин Иван Николаевич

Научный сотрудник НИИ радиотехнических систем при ТУСУРе, Тел.: (3822) 41 38 92 Эл. почта: iyme@mail.ru

Yu.P. Akulinichev, P.V. Abramov, I.N. Vaulin

The absorbing layer influence on the parabolic equation solutions

For many radiowave propagation problems the bottom boundary of the domain is the physical ground/ air interface, often represented by a surface impedance boundary condition, while the top boundary is a computational artifact needed to limit the integration domain in height. The top boundary must represent at a finite distance the Sommerfeld radiation condition - it must be perfectly transparent, letting all the energy coming from below the boundary escape to infinity. Though the absorbing layers are the popular truncating technique there are no investigation, as far as it's known, of its influence on the numerical parabolic equation solution. This paper is devoted for these questions._