H(2,4) 기법을 기반으로 한 저분산 FDTD 기법의 손실 매질의 광대역 해석을 위한 최적화 방법

Ⅰ. 서 론

Yee에 의해서 유한차분 시간영역법(Finite-Difference Time-Domain: FDTD)이 개발된 이후로, FDTD 기법은 다양한 전자기 현상 예측 및 분석에 이용되고 있다^[1]～^[5]. FDTD 기법은 유한요소법(Finite Element Method: FEM), 모멘트 법(Method Of Moment: MoM)과 비교해서 전자기 현상의 광대역 특성을 빠르게 해석 할 수 있다는 장점이 있다. 반면 계산 결과가 FEM과 MoM의 계산 결과와 비교했을 때, 정확도가 낮다는 단점이 있다. 이외에 FDTD 기법은 진행방향에 따른 수치 파수 값(numerical wavenumber value)이 일정하지 않다는 단점도 있다. 위의 두 단점을 해결하기 위해서 다양한 연구가 수행되었다^[6]～^[10]. 그 중 Isotropic-Dispersion FDTD와 H(2,4) 기법에 기반한 저분산 FDTD (Extremely Low Numerical Dispersion FDTD: ELND-FDTD) 기법^[6]은 협대역 해석에서 매우 정확한 특성을 갖는다. 하지만 위의 FDTD 기법을 이용한 광대역 해석(wide-band analysis)을 위한 연구는 거의 진행이 되지 않았다.

H(2,4) 기법이 사용하는 전기장과 자기장의 위치는 그림 1(a)에 나타나 있다. 그림 1(a)에서 보는 것과 같이 Yee 기법에 이용하는 전기장과 자기장은 H(2,4) 기법에서도 이용된다. 그리고 그림 1(b)에서 알 수 있듯이 H(2,4) 기법의 경우는 전파 진행방향이 0°와 90°에서 가장 정확하고, Yee 기법의 경우는 45°에서 가장 정확한 값을 갖는다. 두 기법의 이러한 상호보완적인 관계를 이용하여 ELND-FDTD 기법은 단일 주파수에서 매우 정확한 해를 얻었다 ^[6],^[7]. 두 기법을 가중합(weighted sum) 처리하여 전파 진행방향에 따라 변화하지 않는 파수를 만들었고, 비유전율과 도전율을 정확하게 예측하기 위한 두 개의 변수를 이용하였다. 이 방법은 H(2,4) 기법과 비교했을 때, 단일 주파수에서 월등한 정확성을 가지면서도 계산 시간은 전혀 증가하지 않는 장점이 있다. 본 논문은 ELND-FDTD 기법을 기반으로 하여 광대역에서 정확한 전자기 특성을 얻을 수 있는 최적화 방법을 제안하였다.

그림 1. | Fig. 1. (a) TMz mode 에서의 H(2,4) 기법과 Yee FDTD 기법의 전계 및 자계 계산점 위치, (b) H(2,4) 및 Yee 기법을 이용하여 계산한 파수와 이론적 파수의 차이. ε_r=4.4, σ=0.1 S/m, Δt=2.5 ps, Δ=1.3 mm | (a) Calculation positions of E-/H- field for TMz mode using H(2,4) and Yee FDTD method, (b) Numerical wavenumber of H(2,4) and Yee FDTD method and analytic wavenumber. ε_r=4.4, σ=0.1 S/m, Δt=2.5 ps, Δ=1.3 mm.

Download Original Figure

Ⅱ에서는 H(2,4) 기법을 기반으로 한 저분산 시간영역 유한차분법의 갱신 방정식(update equations)과 세 개의 변수를 이용하여 해석 주파수 범위에서 오차를 최소화하는 방법에 대해서 자세히 서술하고, 제안한 방법의 이론적 성능을 확인했다. Ⅲ에서는 2차원 원형 유전체의 광대역 산란 문제를 계산해서 기존의 방법과 제안한 방법의 광대역 성능을 비교했다.

Ⅱ. Update Equations and Optimization Factors

본 논문에서 제안한 FDTD 방법은 H(2,4) 기법을 기반으로 하고 있다. 그림 1(a)는 H(2.4) 기법을 이용해서 TM_z 모드의 전계를 계산하기 위해 필요한 전계와 자계의 계산 위치를 나타내고 있다. 점선으로 표시된 부분에서 알 수 있듯이, H(2,4) 기법의 계산점들은 Yee 기법의 계산 점을 포함하고 있다. 그림 1(b)는 Yee 기법과 H(2,4) 기법을 이용해 계산된 수치 파수와 정확한 파수를 진행방향에 따라 나타냈다. H(2,4) 기법은 0°에서 가장 정확한 값을 가지며, Yee 기법은 45°에서 가장 정확한 값을 가지며, 두 기법의 진행 방향에 따른 수치 파수 변화가 비슷한 모양을 갖는다. 이 특성을 이용해 진행방향에 따라 변화하지 않는 수치 파수를 만들기 위해서 H(2,4) 기법과 Yee 기법을 식 (1)과 같이 가중합을 적용하였다. F_w는 가중치를 나타내고, D_Y, D_H는 각각 Yee 기법과 H(2,4) 기법의 파수의 최대 변화를 나타낸다.

식 (2), (3)은 FDTD 갱신방정식에서 이용하는 맥스웰 방정식(Maxwell’s equations)을 나타낸다. 식 (4)～(6)은 맥스웰 방정식으로 TM_z 모드로 전개 했을 때의 방정식을 나타낸다. 식 (7)～(9)는 TM_z 모드에서의 ELND-FDTD 기법의 갱신방정식을 나타낸다. F_ϵ, F_σ, Δ(Δx=Δy=Δ), Δt, ϵ₀, μ, σ는 각각 비유전율에 적용하는 비례치, 도전율에 적용하는 비례치, 격자 길이, 단위 시간, 자유공간에서의 유전율, 투자율, 도전율을 각각 나타낸다. $A=\frac{2F_w{ϵ}_0{ϵ}_r-F_{\sigma }\sigma \Delta t}{2F_w{ϵ}_0{ϵ}_r+F_{\sigma }\sigma \Delta t},B=\frac{2\Delta t\left(8+F_w\right)}{8\left(2F_{ϵ}{ϵ}_0{ϵ}_r+F_{\sigma }\sigma \Delta t\right)\Delta },C=\frac{2F_w\Delta t}{8\left(2F_{ϵ}{ϵ}_0{ϵ}_r+F_{\sigma }\sigma \Delta t\right)\Delta }$로 정의 된다. 갱신 방정식 의 F_w, F_ϵ, F_σ에 따라 식 (4)～(6)의 갱신 방정식은 Yee, H(2,4), ELND-FDTD, 제안한 방법으로 변화하게 된다. 이에 따른 변화는 표 1에 나타나 있다.

식 (7)～(9)에 있는 F_w는 진행방향에 따른 파수의 변화를 제거하는 역할을 하는 변수이고, F_ϵ과 F_σ는 비례치 (Scaling Factor: SF)로 각각 유전율과 도전율에 곱해져서 파수를 이론값과 일치하도록 만드는 역할을 한다. 식 (1)을 이용하여 F_w를 계산하면 그림 2(a)와 같이 나타난다. 파장 당 셀 수(Cells Per Wavelength: CPW)가 증가할수록 F_w가 1로 근접함을 알 수 있다. 또한 CPW가 증가할수록 F_w의 변화가 작아진다. 또한 F_w는 제안한 방법의 안정성(stability)에도 영향을 미치게 된다. 안정성은 사용할 수 있는 최대 단위 시간(maximum unit time step)을 결정하고, 수치 분산 관계를 이용하여 최대 단위 시간을 계산할 수 있다^[11]. 수치 분산관계는 식 (10)으로 표현할 수 있고, 식 (10)을 이용하여 최대 단위 시간은 식 (11)로 구할 수 있게 된다. 그림 2(b)는 이를 이용하여 H(2,4) 기법의 최대 단위 시간과 제안한 방법의 단위 시간을 비교했다. 가중치(weighting factor)는 그림 2(a)에 나온 범위를 이용했다. 가중치가 증가할수록 두 기법의 최대 단위 시간의 비율이 작아지지만, 제안한 방법의 가중치가 1에 근접하기 때문에 그 변화는 크지 않은 것을 확인할 수 있다.

그림 2. | Fig. 2. (a) CPW에 따른 F_w 특성, (b) H(2,4) 기법과 제안한 방법의 최대 단위시간 비교, (c) 이론 값과 F_w, F_ϵ, F_σ에 따른 수치 파수 변화 | (a) Weighting factor versus CPW, (b) Comparison between maximum time step of H(2,4) method and proposed method, (c) Analytic wavenumber and numerical wavenumber change due to scaling factors.

Download Original Figure

그림 2(c)에서 점선 그래프(ELND-FDTD w/o SF)는 ELND-FDTD 기법에 F_w의 값만 적용했을 때의 파수를 나타낸다. F_w의 효과로 인해 진행방향에 따라 변하지 않는 수치 파수를 얻게 되었다. 즉, 상수 수치 비유전율(ϵ_rn)과 상수 수치 도전율(σ_n)을 얻게 되었다. 하지만 그림 2(c)에서 보는 것과 같이 이론값과의 차이가 있는 것을 확인할 수 있다. 이를 해결하기 위해서 두 개의 비례 변수(F_ϵ, F_σ)를 도입하였다. 이 변수들은 각각 ϵ_rn과 σ_n에 곱해지게 된다. ϵ_rn과 σ_n을 이론값과 같게 만들기 위해서 F_ϵ과 F_σ는 식 (12), (13)으로 정의 했다. F_ϵ과 F_σ의 효과로 인해서 이론값과 같은 수치 파수를 얻게 된 것을 그림 2(c)의 그래프(ELND-FDTD w/ SF, Analytic)를 통해 확인할 수 있다.

그림 3은 H(2,4)기법, ELND-FDTD 기법의 오차 곡선을 나타내고 있다. 오차는 각 주파수에서의 비유전율과 도전율 오차를 동시에 반영하여 식 (14)로 정의했다. ϵ_ra, ϵ_rn, σ_a, σ_n은 각각 이론 비유전율과 수치 비유전율, 이론 도전율과 수치 도전율을 각각 나타낸다. 최적화 주파수(f_opt)는 F_w, F_ϵ, F_σ가 계산되는 주파수를 의미하고, f는 해석 주파수 범위를 의미한다. ELND-FDTD 기법에서 F_w, F_ϵ, F_σ는 6 GHz와 10 GHz에서 정확한 값을 갖도록 계산되어졌다. ELND-FDTD 기법은 최적화 주파수에서는 매우 작은 오차를 갖지만, 최적화 주파수에서 멀어질수록 오차가 급격히 증가하는 것을 볼 수 있다. 특히 10 GHz에서 최적화를 시켰을 때의 오차 곡선은 7 GHz보다 작은 주파수에서는 H(2,4) 기법보다 부정확한 특성을 갖는다. 반면, 6 GHz에 최적화 시켰을 때의 오차 곡선을 보면 4 GHz 근처에서 H(2,4) 기법보다 부정확한 특성을 갖게 된다. 즉, ELND-FDTD 기법의 경우, 최적화 된 주파수에 따라서 광대역 주파수 해석의 결과의 정확도가 변화한다.

그림 3. | Fig. 3. H(2,4) 기법과 6 GHz, 10 GHz에서 최적화시켰을 때의 ELND-FDTD 기법의 오차 곡선 | Error of H(2,4) FDTD scheme and ELND-FDTD optimized 6 GHz and 10 GHz.

Download Original Figure

광대역 해석의 정확도를 높이기 위해서는 해석 주파수 범위에서 오차를 가장 작게 만드는 최적의 F_w, F_ϵ, F_σ의 값을 계산에 이용되는 모든 물질에 대해서 각각 찾아 적용해야 한다. 그림 4는 최적의 변수를 찾는 과정을 나타내고 있다. 식 (15)는 본 논문에서 계산한 광대역 오차를 나타내는 수식으로 식 (14)에서 변수 f를 모두 더한 값이 된다. 그림 5는 ϵ_ra가 4.4이고, σ_a가 0.1 S/m, Δt는 2.5 ps, Δ는 1.3 mm이고, f_opt와 f 모두 4～10 GHz일 때의 오차 곡선이다. 오차가 가장 작은 경우는 f_opt가 6.24 GHz 일 때이다.

그림 4. | Fig. 4. 제안한 알고리즘의 흐름도 | Flow chart of the proposed FDTD method.

Download Original Figure

그림 5. | Fig. 5. 광대역 오차 곡선 | Wideband error graph.

Download Original Figure

Ⅲ. Numerical Results

본 논문에서는 제안한 방법의 우수한 광대역 해석 특성을 확인하기 위해서 제안한 방법(proposed)과 10 GHz 에서 최적화한 ELND-FDTD 기법, H(2,4) 기법으로 2차원 실린더 유전체의 광대역 산란 문제를 계산하였다. 계산 조건은 다음과 같다. 유전체 반지름은 299 mm, 비유전율은 4.4, 도전율은 0.1 S/m이다. 격자 길이(Δ)와 단위 시간(Δt)은 각각 1.3 mm, 2.5 ps이다. 해석 주파수 범위는 4～10 GHz로 정했다. 그림 6은 계산 구조 및 관찰 부분을 나타내고 있고, 두 물질의 경계 값은 두 물질 상수의 산술 평균을 취하였다.

그림 6. | Fig. 6. 계산 구조 및 관찰 부분 | Calculation structure and observation line.

Download Original Figure

그림 7은 네 개의 주파수(4 GHz, 6.24 GHz, 7.76 GHz, 9.94 GHz)에서의 유전체 실린더에 의해 산란된 전기장의 실수부와 허수부의 차이를 나타낸다(위: 실수부, 아래: 허수부). 그림 7(b), (c)에서 볼 수 있듯이 제안한 방법이 실수부와 허수부에서 이론값^[12]과 가장 근사한 값을 갖는 것을 확인할 수 있다. 그림 7(a)의 경우에는 H(2,4) 기법과 제안한 방법의 정확도가 유사하게 나타난다. 하지만, ELND-FDTD 기법의 결과는 최적화 주파수인 10 GHz에 근접한 주파수를 제외하고는 해석 주파수 범위에서 H(2,4) 기법의 결과보다도 부정확하게 나타남을 확인할 수 있다.

그림 7. | Fig. 7. 이론 값과 계산된 전계 차이. 위: 실수부, 아래: 허수부 | Difference electric field intensity of analytic results and numerical results. upside: Real part, bottom side: Imaginary part.

Download Original Figure

그림 7에서 특정 주파수에서의 전계 분포 비교를 통해 정확성을 확인했다면 그림 8은 정량적인 분석을 위해서 해석 주파수 범위에서의 Root Mean Square(RMS) 값을 500 MHz 간격으로 계산한 결과를 나타내고 있다. RMS는 식 (16)을 이용하였다. 식 (16)에서 Re은 실수부분을 Im는 허수부분을 나타내고, 아래 첨자 A는 이론값을 C는 계산 결과를 의미한다. y는 관찰한 위치를 의미하고, n은 관찰한 전계의 수를 각각 의미한다. 제안한 방법(proposed)을 이용하여 계산한 것이 10 GHz에서 최적화한 ELND-FDTD 기법보다 4～8 GHz 범위에서 정확성이 높은 것을 확인할 수 있고, H(2,4) 기법보다는 전 대역에서 정확도가 높은 것을 확인할 수 있다.

그림 8. | Fig. 8. RMS 오차 비교 | Comparison of RMS error.

Download Original Figure

Ⅳ. Conclusion

본 논문은 ELND-FDTD 기법을 기반으로하여 광대역 에서도 정확한 전자기 특성을 예측할 수 있는 FDTD 방법에 대해서 제안했다. 계산 시간 및 알고리즘의 복잡도 증가 없이 해석 주파수 범위에서 가장 작은 오차를 갖도록 하는 최적화 방법을 제안했다. H(2,4) 기법과 ELND-FDTD 기법과 제안한 방법의 결과를 이론값과 비교하여 제안한 방법의 우수성을 확인하였다. 제안한 방법은 정밀한 전자기 현상 해석이 필요한 분야에서 매우 유용하게 이용될 수 있을 것이다.