Александр Зеленин.
Граничные условия для FDTD.
В счетном объеме каждый вектор Е или Н вычисляется через 4 соседних вектора. Так происходит по всему объему. Но на границах самые последние векторы Е имеют: на гранях параллелепипеда счетного объема только три соседних вектора Н из четырех необходимых; на ребрах - два. Поэтому точно вычислить поле Е на границах невозможно.
Проблема вычисления граничных полей решается различными способами.
1. Условия PEC - идеальный проводник.
Условия PEC таковы, что граничные вектора Е никогда не вычисляются, а, значит, всегда равны нулю. Как известно, поле Е всегда равно нулю в идеальном проводнике, поэтому такие границы ведут себя как идеальный проводник: электромагнитные волны 100 % отражаются обратно в счетный объем.
2. Условия симметрии.
В некоторых случаях поле Е или поле Н может быть симметричным относительно некоторой плоскости. Тогда в этой плоскости можно задать условие симметрии и тем самым вдвое уменьшить счетный объем. При этом рядом с данной плоскостью симметрии будет проходить граница счетного объема с условиями симметрии.
Симметрия может быть четной или нечетной.
При нечетной симметрии плоскость симметрии проходит внутри счетного объема параллельно грани на расстоянии пол-ячейки от грани. Условия нечетной симметрии для симметрии по Е получаются простым переносом значений ближайших к границе векторов Е на саму границу, а для симметрии по Н получаются так же, но при этом вектор Е меняет свой знак. Ниже приведены примеры нечетной симметрии для Е и Н.
Пример: нечетная симметрия по Е:
// ПЛОСКОСТЬ Y=1/2 симметрична по E
For I:=1 to NX-1 do
For K:=1 to NZ do EX[I,1,K]:=EX[I,2,K];
For I:=1 to NX do
For K:=1 to NZ-1 do EZ[I,1,K]:=EZ[I,2,K];
Пример: нечетная симметрия по Н:
// ПЛОСКОСТЬ Y=1/2 симметрична по H
For I:=1 to NX-1 do
For K:=1 to NZ do EX[I,1,K]:=-EX[I,2,K];
For I:=1 to NX do
For K:=1 to NZ-1 do EZ[I,1,K]:=-EZ[I,2,K];
Условия четной симметрии несколько сложнее. Плоскость симметрии проходит на расстоянии целой ячейки от границы, поэтому кроме поля Е на границе необходимо помнить о прилегающих к границе векторах Н. При этом симметрии по Е и Н отличаются друг от друга только знаками переносимых значений.
Пример: четная симметрия по Н:
// ПЛОСКОСТЬ Y=1 симметрична по H, вычисление поля Е
begin
For K:=1 to NZ do
For I:=1 to NX-1 do EX[I,1,K]:=-EX[I,3,K];
For K:=1 to NZ-1 do
For I:=1 to NX do EZ[I,1,K]:=-EZ[I,3,K];
end;
// ПЛОСКОСТЬ Y=1 симметрична по H, вычисление поля Н
begin
For K:=1 to NZ-1 do
For I:=1 to NX do HX[I,1,K]:=HX[I,2,K];
For K:=1 to NZ do
For I:=1 to NX-1 do HZ[I,1,K]:=HZ[I,2,K];
end;
Пример: четная симметрия по Е:
// ПЛОСКОСТЬ Y=1 симметрична по Е, вычисление поля Е
begin
For K:=1 to NZ do
For I:=1 to NX-1 do EX[I,1,K]:=EX[I,3,K];
For K:=1 to NZ-1 do
For I:=1 to NX do EZ[I,1,K]:=EZ[I,3,K];
end;
// ПЛОСКОСТЬ Y=1 симметрична по Е, вычисление поля Н
begin
For K:=1 to NZ-1 do
For I:=1 to NX do HX[I,1,K]:=-HX[I,2,K];
For K:=1 to NZ do
For I:=1 to NX-1 do HZ[I,1,K]:=-HZ[I,2,K];
end;
При четной симметрии граничные поля Е и Н устанавливаются не одновременно, а на соответствующих полушагах по времени.
3. Простые условия поглощения (АВС)
Для условий поглощения значения векторов электрического поля на границе вычисляются на основании известных полей в приграничных слоях. Причем берутся приграничные поля не только на текущем шаге по времени, но и на предыдущих шагах. Все эти условия пытаются спрогнозировать поле на границе.
В литературе встречается описание целого ряда простых условий поглощения разных авторов. Всего их наберется с десяток. Практически чаще всего применяют условия Мура (Mur) и Лиао (Liao). Остальные не применяют или используют очень редко из-за их более низкой эффективности (Trefethen-Halpern, Higdon) или неудобства в использовании (Retarded Time - RT), или из-за неприменимости в декартовых координатах (Bayliss-Turkel), или из-за тенденции к потере стабильности (относится почти ко всем условиям, но особо - к условиям, основанным на высших порядках точности конечных разностей).
Все условия имеют довольно низкий коэффициент отражения от границы, составляющий порядка 0,1..1 %, но только при падении волны на границу под прямым углом. При падении под острым углом коэффициент отражения растет вплоть до 100 % при падении по касательной. Из-за этого границы необходимо располагать как можно дальше от источника электромагнитных волн, чтобы волны приходили к границе под как можно большими углами, желательно по нормали к границе.
Следует отметить, что оценка эффективности тех или иных простых граничных условий различна у разных авторов. Например, один источник пишет, что условия Лиао на 20 dB эффективнее условий Мура 2-го порядка. Другой пишет, что условия Лиао на 12 dB хуже условий Мура 1-го порядка. Имеется ввиду коэффициент отражения. И оба доказывают свои выводы графиками сравнительных расчетов.
Истина, скорее всего, посередине: для каждой конкретной задачи имеются свои оптимальные граничные условия. Одно плохо - заранее никогда не знаешь, КАКИЕ лучше.
Примечание. Возможно, введение чисел двойной точности радикально повысит стабильность, но пока это непозволительная роскошь, расходующая в два раза больше памяти и времени.
Итак, граничных условий много. Здесь рассмотрим три варианта простейших граничных условий: Мура 1-го порядка, Лиао 3-го порядка и RT. В таблице S (большое) и s (малое) - разные переменные! S = (dT * c)/ D - число Куранта, где D - шаг по пространству, dT- шаг по времени, с- скорость света, а s - координата границы. s-1 - одна ячейка внутрь счетного объема, s-2 - две ячейки от границы и т.д. На поясняющих рисунках от границы (surface - внешней поверхности счетного объема) по горизонтали отложена координата s, а по вертикали - время в шагах счета (n - временной индекс). Формула приведена в двух видах, дополняющих друг друга в пространстве и времени. На рисунках пустой кружок - вычисляемое значение, черные кружки - требуемые значения. Сразу можно определить, какие переменные приграничной области нужно хранить в дополнительных массивах и сколько шагов по времени они должны храниться.
|
Условия |
Формула |
Поясняющая картинка |
|
Mur 1-го порядка |
типично S=0,5 для минимального отражения под прямым углом, но может меняться для получения минимального коэффициента отражения под другими углами. |
|
|
Liao 3-го порядка |
Для данного случая S=0,5. Примечание: в литературе условия Лиао часто ставят особняком с условиями ABC, т.к. они базируются на других исходных предпосылках. |
|
|
RT-ABC |
Для данного случая S=0,5, шаг по времени строго полшага по пространству*, шаг по пространству строго одинаков во всех направлениях. |
|
*Когда речь идет о сопоставлении времени и размеров пространства, то, конечно, секунды с метрами не сравниваются. Время превращается в расстояние при умножении на скорость. В данном случае шаг по времени умножается на скорость света в вакууме и результат сравнивается с шагом по пространству: (dT c) = 0,5 D или, иначе, dT = 0,5 D/ c.
Видно, что больше всего массивов переменных нужно хранить для условий RT (5), меньше всего - для условий Мура (1). Для условий Лиао нужно 3 дополнительных массива переменных.
А теперь пример реализации условий Мура 1-го порядка.
Procedure MursFirstABC;
var I,J,K: Integer;
begin
// УСЛОВИЯ ПОГЛОЩЕНИЯ. Массивы вида EZX1, EZXN и т.п. - массивы сохраненных на предыдущем шаге по времени значений поля Е в приграничной области. Процедура сохранения приводится ниже.
// коэффициенты Tx, Ty и Tz равны c*dT/dX, c*dT/dY и c*dT/dZ соответственно. c-скорость света, dT- шаг по времени; dX, dY, dZ - шаги по пространству.
//ЛЕВАЯ И ПРАВАЯ ПЛОСКОСТИ
// Ez для двух плоскостей YxZ (левой и правой)
for J:=2 to NY-1 do
for K:=1 to NZ-1 do begin
EZ[1,J,K]:=EZX1[J,K]+(Tx-1)/(Tx+1)*(EZ[2,J,K]-EZ[1,J,K]);
EZ[NX,J,K]:=EZXN[J,K]+(Tx-1)/(Tx+1)*(EZ[NX-1,J,K]-EZ[NX,J,K]);
end;
// Ey для двух плоскостей YxZ [левой и правой]
for J:=1 to NY-1 do
for K:=2 to NZ-1 do begin
EY[1,J,K]:=EYX1[J,K]+(Tx-1)/(Tx+1)* (EY[2,J,K]-EY[1,J,K]);
EY[NX,J,K]:=EYXN[J,K]+(Tx-1)/(Tx+1)* (EY[NX-1,J,K] -EY[NX,J,K]);
end;
//ВЕРХНЯЯ И НИЖНЯЯ ПЛОСКОСТИ
// Ex для двух плоскостей XxZ [верхней и нижней]
for K:=2 to NZ-1 do
for I:=1 to NX-1 do begin
EX[I,1,K]:= EXY1[I,K]+(Ty-1)/(Ty+1)*(EX[I,2,K]-EX[I,1,K]);
EX[I,NY,K]:= EXYN[I,K]+(Ty-1)/(Ty+1)*(EX[I,NY-1,K]-EX[I,NY,K]);
end;
// Ez для двух плоскостей XxZ [верхней и нижней]
for K:=1 to NZ-1 do
for I:=2 to NX-1 do begin
EZ[I,1,K]:= EZY1[I,K]+(Ty-1)/(Ty+1)*(EZ[I,2,K]-EZ[I,1,K]);
EZ[I,NY,K]:= EZYN[I,K]+(Ty-1)/(Ty+1)*(EZ[I,NY-1,K]-EZ[I,NY,K]);
end;
//БЛИЖНЯЯ И ДАЛЬНЯЯ ПЛОСКОСТИ
// Ey для двух плоскостей XxY [ближней и дальней]
for I:=2 to NX-1 do
for J:=1 to NY-1 do begin
EY[I,J,NZ]:=EYZN[I,J]+(Tz-1)/(Tz+1)*(EY[I,J,NZ-1]-EY[I,J,NZ]);
EY[I,J,1]:=EYZ1[I,J]+(Tz-1)/(Tz+1)*(EY[I,J,2]-EY[I,J,1]);
end;
// Ex для двух плоскостей XxY [ближней и дальней]
for I:=1 to NX-1 do
for J:=2 to NY-1 do begin
EX[I,J,NZ]:= EXZN[I,J]+(Tz-1)/(Tz+1)*(EX[I,J,NZ-1]-EX[I,J,NZ]);
EX[I,J,1]:= EXZ1[I,J]+(Tz-1)/(Tz+1)*(EX[I,J,2]-EX[I,J,1]);
end;
// РЕБРА, параллельные оси Z
for K:=1 to NZ-1 do begin
EZ[1,1,K]:=EZX1[1,K]+(Tx-1)/(Tx+1)*(EZ[2,1,K]-EZ[1,1,K]);
EZ[1,NY,K]:=EZX1[NY,K]+(Tx-1)/(Tx+1)*(EZ[2,NY,K]-EZ[1,NY,K]);
EZ[NX,1,K]:=EZXN[1,K]+(Tx-1)/(Tx+1)*(EZ[NX-1,1,K]-EZ[NX,1,K]);
EZ[NX,NY,K]:=EZXN[NY,K]+(Tx-1)/(Tx+1)*(EZ[NX-1,NY,K]-EZ[NX,NY,K]);
end;
// РЕБРА, параллельные оси Y
for J:=1 to NY-1 do begin
EY[1,J,1]:=EYX1[J,1]+(Tx-1)/(Tx+1)*(EY[2,J,1]-EY[1,J,1]);
EY[1,J,NZ]:=EYX1[J,NZ]+(Tx-1)/(Tx+1)*(EY[2,J,NZ]-EY[1,J,NZ]);
EY[NX,J,1]:=EYXN[J,1]+(Tx-1)/(Tx+1)*(EY[NX-1,J,1]-EY[NX,J,1]);
EY[NX,J,NZ]:=EYXN[J,NZ]+(Tx-1)/(Tx+1)*(EY[NX-1,J,NZ]-EY[NX,J,NZ]);
end;
// РЕБРА, параллельные оси X
for I:=1 to NX-1 do begin
EX[I,1,1]:=EXY1[I,1]+ (Ty-1)/(Ty+1)*(EX[I,2,1]-EX[I,1,1]);
EX[I,1,NZ]:=EXYN[I,NZ]+(Ty-1)/(Ty+1)*(EX[I,2,NZ]-EX[I,1,NZ]);
EX[I,NY,1]:=EXY1[I,1]+(Ty-1)/(Ty+1)*(EX[I,NY-1,1]-EX[I,NY,1]);
EX[I,NY,NZ]:=EXYN[I,NZ]+(Ty-1)/(Ty+1)*(EX[I,NY-1,NZ]-EX[I,NY,NZ]);
end;
end;
// сохранение приграничных полей для MursFirstABC
Procedure GetBorderValues;
var I,J,K: Integer;
begin
// слева и справа
for J:=1 to NY do for k:=1 to NZ-1 do
begin
EZX1[J,K]:=EZ[2,J,K];
EZXN[J,K]:=EZ[NX-1,J,K];
end;
for J:=1 to NY-1 do for k:=1 to NZ do
begin
EYX1[J,K]:=EY[2,J,K];
EYXN[J,K]:=EY[NX-1,J,K];
end;
// верх и низ
for I:=1 to NX do for k:=1 to NZ-1 do
begin
EZY1[I,K]:=EZ[I,2,K];
EZYN[I,K]:=EZ[I,NY-1,K];
end;
for I:=1 to NX-1 do for k:=1 to NZ do
begin
EXY1[I,K]:=EX[I,2,K];
EXYN[I,K]:=EX[I,NY-1,K];
end;
// ближняя и дальняя
for I:=1 to NX do for J:=1 to NY-1 do
begin
EYZ1[I,J]:=EY[I,J,2];
EYZN[I,J]:=EY[I,J,NZ-1];
end;
for I:=1 to NX-1 do for J:=1 to NY do
begin
EXZ1[I,J]:=EX[I,J,2];
EXZN[I,J]:=EX[I,J,NZ-1];
end;
end;
Примечание: в литературе встречается и несколько иная форма записи условий Мура 1-го порядка. К приведенной здесь формуле добавляется еще один член, который сам Мур относит к условиям второго порядка, а другие авторы утверждают, что этот член должен быть и в формуле 1-го порядка.
Далее об условиях PML.
