Comparison of Absorbing Boundary Conditions and Waveguide Port Boundary Condition for Waveguide Electromagnetic Analysis Using Finite Element Method

유한요소법을 이용한 도파관 전자기 시뮬레이션에 있어 흡수경계조건 및 도파관 포트 경계조건 고찰 및 비교

Abstract

Waveguides are transmission lines that guide electromagnetic waves in the desired direction and are utilized in various fields such as medical devices, radar systems, and satellite communications. Computational electromagnetics (CEM) is essential for designing and optimizing waveguides. The finite element method (FEM), which is one of the numerical analysis techniques, is efficient in solving closed problems such as waveguides. In order to apply FEM for waveguide analysis, boundary conditions that truncate the computational domain are required. This paper performs electromagnetic simulations using absorbing boundary conditions (ABC) and waveguide port boundary conditions (WPBC) in 2/D and 3/D waveguides using the finite element method and compared their performances. The accuracy of the analysis was verified by comparing the results with HFSS, a representative commercial electromagnetic simulation software. Simulation results confirm that applying WPBC allows for smaller analysis domains than ABC.

도파관(waveguide)은 전자기파를 원하는 방향으로 안내하는 전송선로로 의료기기, 레이더 시스템, 위성 통신 등 다양한 분야에 활용되고 있다. 이러한 도파관의 설계 및 최적화를 위해서는 전자기 수치해석(CEM: Computational Electromagnetics)이 필수적이다. 수치해석 기법의 하나인 유한요소법(FEM: Finite Element Method)은 도파관과 같은 닫힌 영역 내부의 전자기 문제를 해결하는데 효율적이며 이를 적용하기 위해서는 계산영역을 한정시키기 위한 경계조건이 필요하다. 본 논문에서는 계산영역 밖으로 나가는 전자파의 반사를 최소화하기 위한 흡수경계조건(ABC: Absorbing Boundary Condition)과 주 모드 뿐만 아니라 고차 모드까지 흡수할 수 있는 도파관 포트 경계조건(WPBC: Waveguide Port Boundary Condition)을 각각 적용하여 2/D 및 3/D 도파관 구조에 대한 전자기 시뮬레이션을 수행하였다. 이후, 대표적인 전자파 상용 소프트웨어인 HFSS와의 결과 비교를 통해 해석의 정확성을 검증하였으며, 시뮬레이션 결과를 통해 WPBC를 적용하면 ABC보다 더 작은 해석 영역으로 구조 해석이 가능하다는 것을 확인하였다.

1. 서론

전자파 문제 해석은 통상적으로 Maxwell 방정식을 푸는 것에서 시작되나, 현실 세계에서 이 방정식을 직접 풀어 해석 해(analytical solution)를 구하는 것은 매우 복잡할 경우가 많다. 따라서 Maxwell 방정식에 대한 수치 해(numerical solution)를 컴퓨터를 통해 찾는 전자기 수치해석(CEM : Computational Electromagnetics)의 도입이 필요하다[1].

컴퓨터를 이용하여 수치 해를 찾는 다양한 기법 중 유한요소법(FEM: Finite Element Method)은 기하학적으로 복잡한 구조 및 비균질성 매질의 해석에 강점을 가진다[1, 2]. 또한, 유한요소 과정을 통해 얻게 되는 최종 system matrix는 0을 많이 포함한 sparse matrix가 되기 때문에 반복 솔버(iterative solver)를 사용하여 해를 구할 때 이점을 가질 수 있으며, 분할 알고리듬을 적용한 병렬 연산으로 계산 시간을 단축할 수 있다[1-5].

유한요소법을 이용한 전자파 문제 해석에 있어 계산영역을 한정시키기 위한 특수한 경계조건이 필요하며, 대표적인 특수한 경계조건으로는 ABC와 PML(Perfectly Matched Layer)이 있다[6-9]. ABC는 계산영역을 잘라내는 경계에 가상 경계(fictitious boundary)를 설정하여 장(field)이 흘러가는 동작을 모사하는 경계조건이며, PML은 인공 층(artificial layer)을 설계하여 나가는 장을 흡수하여 반사를 최소화하는 경계조건이다[6]. PML은 입사각과 주파수에 크게 구애받지 않고 나가는 파동을 흡수할 수 있다는 장점이 있지만, 원래 계산영역에서 인공 층에 해당하는 만큼의 계산영역이 확대되기 때문에 계산 시간 및 컴퓨터 메모리가 증가한다는 단점이 있다[7, 8].

ABC의 경우 PML과 달리 추가적인 인공 층이 필요하지 않지만, 일반적으로 주 모드의 전파만을 허용하기 때문에 가상 경계와 산란체 사이의 거리에 따라 흡수 성능이 달라진다는 단점이 있다[2]. 즉, 가상 경계와 산란체 사이 거리가 가까우면 산란체로 인해 발생한 고차모드 때문에 ABC의 흡수 성능이 떨어질 수 있다. 하지만 추가적인 고차모드를 포함함으로써 유도되는 모드 확장 기반 WPBC는 다중 모드 전파를 허용하기 때문에 불연속 물체와 가상 경계 사이 거리를 가깝게 배치하여 계산영역을 축소할 수 있다는 장점이 있다[9-12].

본 논문에서는 FEM을 사용한 inhomogeneous parallelplate waveguide에서의 반사 및 투과 계수를 시뮬레이션하는 2-D 예제를 시작으로, 이를 확장하여 inhomogeneous rectangular waveguide에서의 S11 및 S12 parameter를 시뮬레이션하는 3-D 예제를 구현하였다. 각 예제에서 ABC와 WPBC에 대한 공식화 및 적용 과정에 대해 살펴보았으며, 아울러 두 경계조건의 성능 차이에 대한 평가를 진행하였다.

2. Finite Element Analysis

2.1 2-D Finite Element Analysis

2.1.1 2-D FEM 공식화

2-D 유한요소 해석에서는 그림 1과 같이 사각형 모양의 불연속 매질을 포함하는 평행판 도파관에 자기장(TM₀ 모드)이 +x 방향으로 입사할 때, 유전체의 높이 및 비유전율에 따른 반사 및 투과 계수 시뮬레이션을 수행하였다.

(그림 1) 유전체가 놓인 평행판 도파관 구조[2, 10, 11]

(Figure 1) Parallel-plate waveguide structure loaded with dielectric[2, 10, 11]

2-D FEM에서의 경곗값 문제(boundary-value problem)는 식 (1)로 정의되며, 그림 1에서 자기장 Hz에 대한 지배방정식 scalar Helmholtz equation은 식 (2)와 같다[6, 11].

\(\begin{aligned}-\frac{\partial}{\partial x}\left(\alpha_{x} \frac{\partial \phi}{\partial x}\right)-\frac{\partial}{\partial y}\left(\alpha_{y} \frac{\partial \phi}{\partial y}\right)+\beta \phi=f\\\end{aligned}\)       (1)

\(\begin{aligned}\begin{array}{l}\frac{\partial}{\partial x}\left(\frac{1}{\epsilon_{r}} \frac{\partial H_{z}}{\partial x}\right)+\frac{\partial}{\partial y}\left(\frac{1}{\epsilon_{r}} \frac{\partial H_{z}}{\partial y}\right)+k_{0}^{2} \mu_{r} H_{z} \\ =-\frac{\partial}{\partial x}\left(\frac{1}{\epsilon_{r}} J_{y}\right)+\frac{\partial}{\partial y}\left(\frac{1}{\epsilon_{r}} J_{x}\right)\end{array}\\\end{aligned}\)       (2)

유한요소 해석을 위해서 식 (2)를 식 (1)의 형태에 맞춰야 하며, 그 결과 식 (3)과 같은 관계를 도출할 수 있다[6, 11].

\(\begin{aligned}\phi=H_{z}, \alpha_{x}=\alpha_{y}=\frac{1}{\epsilon_{r}}, \beta=-k_{0}^{2} \epsilon_{r}, f=-j k_{0} Z_{0} J_{z}\\\end{aligned}\)       (3)

Ritz process를 적용하기 위한 e번째 element에 대한 subfunctional Fe은 식 (4)와 같으며, functional F는 식 (5)와 같이 subfunctional Fe의 합으로 표현된다[11].

\(\begin{aligned}\begin{array}{l}F^{e}\left(\phi^{e}\right) \\ \begin{aligned}= & \frac{1}{2} \iint_{\Omega^{e}}\left[\alpha_{x}\left(\frac{\partial \phi^{e}}{\partial x}\right)^{2}+\alpha_{y}\left(\frac{\partial \phi^{e}}{\partial y}\right)^{2}+\beta\left(\phi^{e}\right)^{2}\right] d \Omega \\ & -\iint_{\Omega^{e}} f \phi^{e} d \Omega,\end{aligned} \\\end{array}\\\end{aligned}\)       (4)

\(\begin{aligned}F(\phi)=\sum_{e=1}^{M} F^{e}\left(\phi^{e}\right)\\\end{aligned}\)       (5)

Ritz process 및 assembly 과정을 거치면 최종적으로 식 (6)과 같은 Ax=b 꼴의 선형 대수 방정식을 얻을 수 있으며, 이를 풀면 자기장 Hz에 대한 해를 구할 수 있다.

\(\begin{aligned}\left\{\frac{\partial F}{\partial \phi}\right\}=\sum_{e=1}^{M}\left\{\frac{\partial F^{e}}{\partial \phi^{e}}\right\}=\sum_{e=1}^{M}\left(\left[K^{e}\right]\left[\phi^{e}\right]-\left[b^{e}\right]\right)=\{0\}\\\end{aligned}\)       (6)

2.1.2 2-D FEM에서의 ABC 적용

그림 1 구조에서 x = x₁, x = x₂ 경계를 내부 불연속 물체에 충분히 멀리 떨어진 영역에 배치하여 주 모드만 흡수가 가능하다고 가정하면 x = x₁과 x = x₂에서의 자기장은 아래의 식 (7)로 표현할 수 있다[11].

H_z = H₀e^-jk₀x + RH₀e^jk₀x at x = x₁

H_z = H_z^trans = TH₀e^-jk₀x at x = x₂       (7)

유한요소 적용을 위해서 식 (7)을 boundary condition of the third kind의 일반적인 형태인 식 (8)과 형태를 맞춰주면 식 (9)와 같은 first-order ABC의 식을 얻을 수 있다[11].

\(\begin{aligned}\left(\alpha_{x} \frac{\partial \phi}{\partial x}+\alpha_{y} \frac{\partial \phi}{\partial y}\right) \cdot \hat{n}+\gamma \phi=q\\\end{aligned}\)       (8)

\(\begin{aligned}\begin{array}{ll}\frac{\partial H_{z}}{\partial x}=j k_{0} H_{z}-2 j k_{0} H_{0} e^{-j k_{0} x} & \text { at } x=x_{1} \\ \frac{\partial H_{z}}{\partial x}=-j k_{0} H_{z} & \text { at } x=x_{2}\end{array}\\\end{aligned}\)       (9)

2.1.3 2-D FEM에서의 WPBC 적용

그림 1에서 x = x₁, x = x₂ 경계를 불연속 물체에 가깝게 배치하면 반사파와 투과파는 주 모드 하나로 표현할 수 없으며 주 모드와 불연속 물체로 인해 여기되는 다양한 고차모드의 중첩(superposition)으로 표현해야 한다. 따라서 x = x₁과 x = x₂에서의 자기장은 식 (10)으로 나타낼 수 있으며, h_m(y) 및 γ_m는 식 (11)과 같다[10, 11].

\(\begin{aligned}\begin{array}{ll}H_{z}(x, y)=H_{z}^{i n c}(x, y)+\sum_{m=0}^{\infty} a_{m} h_{m}(y) e^{\gamma_{m} x} \text { at } x=x_{1} \\ H_{z}(x, y)=\sum_{m=0}^{\infty} b_{m} h_{m}(y) e^{-\gamma_{m} x} & \text { at } x=x_{2}\end{array}\\\end{aligned}\)       (10)

\(\begin{aligned}\begin{array}{c}h_{m}(y)=\sqrt{\frac{\nu_{m}}{b}} \cos \frac{m \pi y}{b}, \nu_{m}=\left\{\begin{array}{ll}1 & m=0 \\ 2 & m \neq 0\end{array}\right. \\ \gamma_{m}=\left\{\begin{array}{ll}j \sqrt{k_{0}^{2}-\left(\frac{m \pi}{b}\right)^{2}} & \text { if }\left(\frac{m \pi}{b}\right)^{2} \leq k_{0}^{2} \\ \sqrt{\left(\frac{m \pi}{b}\right)^{2}-k_{0}^{2}} & \text { if }\left(\frac{m \pi}{b}\right)^{2}>k_{0}^{2}\end{array}\right.\end{array}\\\end{aligned}\)       (11)

식 (10)에서의 a_m과 b_m은 expansion coefficient이며, 식 (12)와 같은 직교성(orthogonality) 관계식을 이용하면 식 (13)과 같이 표현될 수 있다[10, 11].

∫₀^bh_m(y)h_m'(y)dy = δ_mm'       (12)

a_m = e^-γ_mx₁∫₀^b[H_z(x₁,y) - H_z^inc(x₁,y)]h_m(y)dy

b_m = e^γ_mx₂∫₀^bH_z(x₂,y)h_m(y)dy       (13)

WPBC 적용을 위해 식 (13)을 식 (10)에 대입 후, 미분을 통하여 식 (8)의 형태로 식을 변형하면 식 (14)와 같은 평행판 도파관에서의 WPBC 적용 식을 얻을 수 있다. 아울러, γ와 q는 식 (15)와 같이 정의된다[10, 11].

\(\begin{aligned}\begin{array}{l}\frac{\partial H_{z}}{\partial n}+\gamma\left(H_{z}\right)=q \text { at } x=x_{1} \\ \frac{\partial H_{z}}{\partial n}+\gamma\left(H_{z}\right)=0 \quad \text { at } \quad x=x_{2} \\\end{array}\\\end{aligned}\)       (14)

\(\begin{aligned}\begin{array}{l}\gamma\left(H_{z}\right)=\sum_{m=0}^{\infty} \gamma_{m} h_{m}(y) \int_{0}^{b} H_{z}\left(x, y^{\prime}\right) h_{m}\left(y^{\prime}\right) d y^{\prime} \\ q=2 \gamma_{n} H_{0} h_{n}(y) e^{-\gamma_{m} x_{1}}\end{array}\\\end{aligned}\)       (15)

2.2 3-D Finite Element Analysis

2.2.1 3-D FEM 공식화

3-D 유한요소 해석에서는 그림 2와 같은 불연속 매질을 포함하는 구형 도파관에 주 모드인 TE₁₀모드가 +z 방향으로 입사할 때, 반사 계수 S11 parameter와 투과 계수 S12 parameter에 대한 시뮬레이션을 수행하였다.

(그림 2) 유전체가 놓인 구형 도파관 구조[11, 13]

(Figure 2) Rectangular waveguide structure loaded with dielectric[11, 13] (b = 10mm, 𝜖_r = 6)

주파수 영역에서 전기장에 대한 지배방정식 vector wave equation은 식 (16)과 같으며, 경계조건은 식 (17) 및 식 (18)과 같이 정의된다[1, 6, 11-13].

\(\begin{aligned}\nabla \times\left(\frac{1}{\mu_{r}} \nabla \times \vec{E}\right)-k_{0}^{2} \epsilon_{r} \vec{E}=-j k_{0} Z_{0} \vec{J}\\\end{aligned}\)       (16)

\(\begin{aligned}\hat{n} \times \vec{E}=0 \quad on \quad S_{1}\\\end{aligned}\)       (17)

\(\begin{aligned}\frac{1}{\mu_{r}} \hat{n} \times(\nabla \times \vec{E})+\gamma_{e} \hat{n} \times(\hat{n} \times \vec{E})=0 \; on \; S_{2}\end{aligned}\) (18)

이후, generalized variational principle의 적용 과정을 거치면 식 (16)에 대한 functional F인 식 (19)를 얻을 수 있다[1, 6, 11-14].

\(\begin{aligned}\begin{array}{l}F(\vec{E}) \\ =\frac{1}{2} \iiint_{V}\left[\frac{1}{\mu_{r}}(\nabla \times \vec{E}) \cdot(\nabla \times \vec{E})-k_{0}^{2} \epsilon_{r} \vec{E} \cdot \vec{E}\right] d V \\ +\iint_{S_{2}}\left[\frac{\gamma_{e}}{2}(\hat{n} \times \vec{E}) \cdot(\hat{n} \times \vec{E})+\vec{E} \cdot \vec{U}\right] d S \\ +j k_{0} Z_{0} \iiint_{V}[\vec{E} \cdot \vec{J}] d V\end{array}\\\end{aligned}\)       (19)

2.2.2 3-D FEM에서의 ABC 적용

Waveguide port를 내부 불연속 물체에 의해 여기된 고차모드가 도달할 수 없을 만큼 충분히 멀리 떨어진 영역에 배치하면 S1 및 S2에서의 전기장은 식 (20)과 같이 표현할 수 있다[2, 11]. 이때, \(\begin{aligned}\overrightarrow{e_{10}}(x, y)\\\end{aligned}\) 와 k_z10은 식 (21)과 같다.

\(\begin{aligned}\begin{array}{l}\vec{E}(x, y, z)=\overrightarrow{\vec{E}}^{n c}(x, y, z)+\overrightarrow{\vec{E}}^{e f}(x, y, z) \\ =E_{0} \vec{e}_{10}(x, y) e^{-j k_{z 110^{z}}}+R E_{0} \vec{e}_{10}(x, y) e^{j k_{z 10} z} \text { on } S_{1} \\ \vec{E}(x, y, z)=\vec{E}^{\text {trans }}(x, y, z) \\ =T E_{0} \vec{e}_{10}(x, y) e^{-j k_{z 10} z} \text { on } S_{2} \\\end{array}\\\end{aligned}\)       (20)

\(\begin{aligned}\begin{array}{l}\vec{e}_{10}(x, y)=\hat{y} \sin \frac{\pi x}{a} \\ k_{z 10}=\sqrt{k_{0}^{2}-\left(\frac{\pi}{a}\right)^{2}}\end{array}\\\end{aligned}\)       (21)

도파관 벽에는 식 (17)과 같은 경계조건이 적용되며, 식 (20)을 식 (18)과 같은 형태로 식을 변형하면 식 (22)와 같은 first-order ABC의 식을 얻을 수 있다[2, 11]. 여기서 \(\begin{aligned}\vec{U}^{n c}\\\end{aligned}\) 및 γ는 식 (23)과 같다.

\(\begin{aligned}\begin{array}{l}\hat{n} \times(\nabla \times \vec{E})+\hat{n} \times(\hat{n} \times \vec{E})=\vec{U}^{n c} \text { on } S_{1} \\ \hat{n} \times(\nabla \times \vec{E})+\hat{n} \times(\hat{n} \times \vec{E})=0 \quad \text { on } S_{2} \\\end{array}\\\end{aligned}\)       (22)

\(\begin{aligned}\begin{array}{l}\gamma=j k_{z 10} \\ \vec{U}^{n c}=-2 j k_{z 10} \vec{E}^{n c}\end{array}\\\end{aligned}\)       (23)

2.2.3 3-D FEM에서의 WPBC 적용

Waveguide port를 내부 불연속 물체에 가깝게 배치하면 S1 및 S2에서의 전기장은 식 (24)의 형태로 나타낼 수 있다[11, 12]. 또한, \(\begin{aligned}\overrightarrow{e_{m n}^{T E}}(x, y), \overrightarrow{e_{t m n}^{T M}}(x, y), \overrightarrow{e_{z m n}^{T M}}(x, y)\\\end{aligned}\)은 식 (25)의 γ_m와 normalization factor N_mn과 함께 식 (26)과 같이 정의된다[11].

\(\begin{aligned}\begin{array}{l}\vec{E}(x, y, z) \\ =\vec{E}^{n c}(x, y, z)+\sum_{m=0}^{\infty} \sum_{m=0}^{\infty} a_{m n} \overrightarrow{e_{m n}^{T E}}(x, y) e^{\gamma_{m n} z} \\ +\sum_{m=1}^{\infty} \sum_{n=1}^{\infty} b_{m n}\left[\overrightarrow{e_{t m n}^{T M}}(x, y)+\hat{z} e_{z m n}^{T M}(x, y)\right] e^{\gamma_{m n} z} \text { on } S_{1}\end{array}\\\end{aligned}\)       (24)

\(\begin{aligned}\begin{array}{l}\vec{E}(x, y, z) \\ =\sum_{m=0}^{\infty} \sum_{m=0}^{\infty} c_{m n} \overrightarrow{e_{m n}^{T E}}(x, y) e^{-\gamma_{m m} z} \\ +\sum_{m=1 n=1}^{\infty} \sum_{m n}^{\infty}\left[\overrightarrow{e_{t m n}^{T M}}(x, y)-\hat{z} e_{z m n}^{T M}(x, y)\right] e^{-\gamma_{m m} z} \text { on } S_{2} \\ \gamma_{m}=\left\{\begin{array}{l}j \sqrt{k_{0}^{2}-\left(\frac{m \pi}{a}\right)^{2}-\left(\frac{n \pi}{b}\right)^{2}} \text { if }\left(\frac{m \pi}{a}\right)^{2}+\left(\frac{n \pi}{b}\right)^{2} \leq k_{0}^{2} \\ \sqrt{\left(\frac{m \pi}{a}\right)^{2}+\left(\frac{n \pi}{b}\right)^{2}-k_{0}^{2}} \text { if }\left(\frac{m \pi}{a}\right)^{2}+\left(\frac{n \pi}{b}\right)^{2}>k_{0}^{2}\end{array}\right.\end{array}\\\end{aligned}\)       (25)

\(\begin{aligned}\begin{array}{l} N_{m n}=\sqrt{\nu_{m} \nu_{n}} / \sqrt{n^{2} \frac{a}{b}+m^{2} \frac{b}{a}} \\ \overrightarrow{e_{m m}^{T E}}(x, y) \\ = N_{m n}\left(\hat{x} \frac{n}{b} \cos \frac{m \pi x}{a} \sin \frac{n \pi y}{b}-\hat{y} \frac{m}{a} \sin \frac{m \pi x}{a} \cos \frac{n \pi y}{b}\right) \\ \overrightarrow{e_{m n}^{T M}}(x, y) \\ = N_{m n}\left(\hat{x} \frac{m}{a} \cos \frac{m \pi x}{a} \sin \frac{n \pi y}{b}+\hat{y} \frac{n}{b} \sin \frac{m \pi x}{a} \cos \frac{n \pi y}{b}\right) \\ e_{z m m}^{T M}(x, y)=\frac{N_{m m}}{\pi \gamma_{m n}}\left[\left(\frac{m \pi}{a}\right)^{2}+\left(\frac{n \pi}{b}\right)^{2}\right] \sin \frac{m \pi x}{a} \sin \frac{n \pi y}{b}\end{array}\\\end{aligned}\)       (26)

\(\begin{aligned}\overrightarrow{e_{m n}^{T E}}(x, y)\\\end{aligned}\)과 \(\begin{aligned}\overrightarrow{e_{m n}^{T E}}(x, y)\\\end{aligned}\)은 식 (27)과 같은 직교성 관계식을 만족하며, 이를 이용하면 각 expansion coefficient는 식 (28)과 같이 표현될 수 있다[11].

\(\begin{aligned}\begin{array}{l}\int_{0}^{a} \int_{0}^{b} e_{m n}^{T E} \cdot \overrightarrow{e_{m n}^{T M}} d x d y=0 \\ \int_{0}^{a} \int_{0}^{b} \overrightarrow{e_{m n}^{T E}} \cdot \overrightarrow{e_{m^{\prime} n^{\prime}}^{T E}} d x d y=\delta_{m m^{\prime}} \delta_{n n^{\prime}} \\ \int_{0}^{a} \int_{0}^{b} \overrightarrow{e_{t m n}^{T M}} \cdot \overrightarrow{e_{t m^{\prime} n^{\prime}}^{T M}} d x d y=\delta_{m m^{\prime}} \delta_{n n^{\prime}}\end{array}\\\end{aligned}\)       (27)

\(\begin{aligned}\begin{array}{l}a_{m n}=e^{-\gamma_{m m} z_{1}} \int_{0}^{a} \int_{0}^{b \overrightarrow{e_{m n}^{T E}}} \cdot\left[\vec{E}-\overrightarrow{E^{i n c}}\right]_{z=z_{1}} d x d y \\ b_{m n}=e^{-\gamma_{m m} z_{1}} \int_{0}^{a} \int_{0}^{b} e_{t m n}^{T M} \cdot\left[\vec{E}-\overrightarrow{E^{i n c}}\right]_{z=z_{1}} d x d y \\ c_{m n}=\left.e^{\gamma_{m u r} z_{2}} \int_{0}^{a} \int_{0}^{b} e_{m n}^{e_{m n}^{T E}} \cdot \vec{E}\right|_{z=z_{2}} d x d y \\ d_{m n}=\left.e^{\gamma_{m m} z_{2}} \int_{0}^{a} \int_{0}^{b} e_{t m n}^{T M} \cdot \vec{E}\right|_{z=z_{2}} d x d y\end{array}\\\end{aligned}\)       (28)

WPBC 적용을 위해 식 (28)을 식 (24)에 대입 후 식 (18)과 형태를 맞춰주면 식 (29)와 같은 구형 도파관에서의 WPBC 적용 식을 얻을 수 있으며, \(\begin{aligned}P(\vec{E})\\\end{aligned}\) 및 \(\begin{aligned}\overrightarrow{U^{i n c}}\\\end{aligned}\)는 식 (30)과 같다[11].

\(\begin{aligned}\begin{array}{l}\hat{n} \times(\nabla \times \vec{E})+P(\vec{E})=\overrightarrow{U^{i n c}} \text { at } z=z_{1} \\ \hat{n} \times(\nabla \times \vec{E})+P(\vec{E})=0 \quad \text { at } z=z_{2}\end{array}\\\end{aligned}\)        (29)

\(\begin{aligned}\begin{array}{l}P(\vec{E})=\sum_{m=1 n=1}^{\infty} \sum_{m}^{\infty} \frac{k_{0}^{2}}{\gamma_{m n}} \overrightarrow{e_{t m n}^{T M}} \int_{0}^{a} \int_{0}^{b} \overrightarrow{e_{t m n}^{T M}} \cdot \vec{E} d x d y \\ -\sum_{m=0 n}^{\infty} \sum_{0}^{\infty} \gamma_{m n} \overrightarrow{e_{m n}^{T E}} \int_{0}^{a} \int_{0}^{b} \overrightarrow{e_{m n}^{T E}} \cdot \vec{E} d x d y \\ \overrightarrow{U^{i n c}}=-2 \gamma_{10} E_{0} \overrightarrow{e_{10}^{T E}} e^{-\gamma_{1 v_{1} z_{1}}}\end{array}\\\end{aligned}\)       (30)

3. 시뮬레이션 결과

3.1 2-D FEM 시뮬레이션 결과

3.1.1 ABC 및 WPBC 적용 2-D FEM

평행판 도파관 구조 해석에서는 2-D linear triangular node element를 기반으로 in-house code(MATLAB R2022b 버전 호환)를 직접 작성하여 반사 및 투과 계수에 대한 시뮬레이션을 수행하였으며, reference[11] 결과와의 비교를 통해 정확성을 검증하였다.

시뮬레이션 구조는 상기 그림 1과 같으며, 기본 환경 주파수는 3GHz(λ=0.1m)로 설정하였다. 그림 3 및 4는 내부 유전체와 1.25λ 떨어진 영역에 first-order ABC를 적용하였을 때 유전체 높이(0~0.35λ) 및 비유전율(𝜖_r = 4, 4 - 1j, 4 - 10j)에 따른 반사 및 투과 계수를 나타내며 reference[11]의 도시와 거의 일치하는 것을 확인할 수 있다.

(그림 3) 𝜖_r = 4, 4 - 1j, 4 - 10j에서 1차 흡수경계조건 적용 반사 계수 비교

(Figure 3) Comparison of first-order ABC application reflection coefficient at 𝜖_r = 4, 4 - 1j, 4 - 10j

(그림 4) 𝜖_r = 4, 4 - 1j, 4 - 10j에서 1차 흡수경계조건 적용 투과 계수 비교

(Figure 4) Comparison of first-order ABC application transmission coefficient at 𝜖_r = 4, 4 - 1j, 4 - 10j

그림 5 및 6은 그림 1 구조에 대한 ABC 적용 2-D FEM과 같은 시뮬레이션 환경에 ABC 대신 WPBC를 적용한 결과를 나타내며, reference[10, 11]의 도시와 거의 일치함이 확인된다.

(그림 5) 𝜖_r = 4, 4 - 1j, 4 - 10j에서 도파관 포트 경계조건 적용 반사 계수 비교

(Figure 5) Comparison of WPBC application reflection coefficient at 𝜖_r = 4, 4 - 1j, 4 - 10j

(그림 6) 𝜖_r = 4, 4 - 1j, 4 - 10j에서 도파관 포트 경계조건 적용 투과 계수 비교

(Figure 6) Comparison of WPBC application transmission coefficient at 𝜖_r = 4, 4 - 1j, 4 - 10j

3.1.2 ABC와 WPBC 성능 비교 (2-D FEM)

해석 영역에 따른 ABC와 WPBC 성능을 비교하기 위해 구조 1에서의 d의 길이를 변화시켜가며 반사 및 투과 계수 시뮬레이션을 수행하였다.

그림 7 및 8은 그림 1의 구조에서 불연속 물체로부터의 거리(d)에 따른 ABC 적용 반사 및 투과 계수 시뮬레이션 결과이다. d의 길이가 작아질수록 불연속 물체로 인해 발생한 고차모드가 시뮬레이션 결과에 영향을 미쳐 정확도(reference[11]과 비교)가 떨어지는 것을 확인할 수 있다.

(그림 7) 𝜖_r = 4에서 1차 흡수경계조건 적용 위치(d)에 따른 반사 계수 비교[11]

(Figure 7) Comparison of reflection coefficient according to distance of first-order ABC application(d) at 𝜖_r = 4[11]

(그림 8) 𝜖_r = 4에서 1차 흡수경계조건 적용 위치(d)에 따른 투과 계수 비교

(Figure 8) Comparison of transmission coefficient according to distance of first-order ABC application(d) at 𝜖_r = 4

그림 9 및 10은 그림 1 구조에서 WPBC 적용 거리(d)에 따른 반사 및 투과 계수 시뮬레이션 결과이며, _{d = 1.25λ}에서는 TM₀모드, _{d = 1.25λ} 및 _{d = 1.25λ}에서는 TM₀모드와 TM₁모드를 흡수하도록 경계조건을 설정하였다. 추가적인 고차모드를 흡수한 WPBC는 ABC보다 더 작은 해석 영역으로도 reference[10, 11] 결과와 잘 일치함이 확인된다.

(그림 9) 𝜖_r = 4에서 도파관 포트 경계조건 적용 위치(d)에 따른 반사 계수 비교

(Figure 9) Comparison of reflection coefficient according to distance of WPBC application(d) at 𝜖_r = 4

(그림 10) 𝜖_r = 4에서 도파관 포트 경계조건 적용 위치(d)에 따른 투과 계수 비교

(Figure 10) Comparison of transmission coefficient according to distance of WPBC application(d) at 𝜖_r = 4

3.2 3-D FEM 시뮬레이션 결과

3-D edge-based brick element를 사용한 inhomogeneous rectangular waveguide 해석은 MATLAB 코드로(R2022b 버전 호환) 작성되었으며, 상용 소프트웨어 HFSS를 사용하여 d=0.8b에 ABC를 적용한 결과(reference)와 in-house의 결과 비교를 통해 시뮬레이션의 유효성을 검증하였다. 이때, 시뮬레이션 구조는 그림 2와 같은 inhomogeneous waveguide이며, 시뮬레이션 환경주파수는 8~13GHz, frequency sweep은 0.02GHz로 설정하였다.

그림 11 및 12는 내부 산란체로부터 떨어진 거리(d)를 변화시키면서 first-order ABC를 적용하였을 때의 S11과 S12 parameter의 시뮬레이션 결과이다. ABC 적용 2-D FEM과 마찬가지로 ABC 적용 거리가 내부 유전체와 가까워질수록 시뮬레이션 정확도가 떨어지는 것을 확인할 수 있다.

(그림 11) 흡수경계조건 적용 위치(d)에 따른 S₁₁ 파라미터 및 Reference와의 결과 비교

(Figure 11) Comparison of S₁₁ parameter (Reference vs. In-house) according to distance of first-order ABC application(d)

(그림 12) 흡수경계조건 적용 위치(d)에 따른 S₁₂ 파라미터 및 Reference와의 결과 비교

(Figure 12) Comparison of S₁₂ parameter (Reference vs. In-house) according to distance of first-order ABC application(d)

그림 13과 14는 그림 2 구조를 내부 유전체로부터 떨어진 거리(d)를 변경해 가며 WPBC를 적용하였을 때의 S11과 S12 parameter 시뮬레이션 결과이며, 시뮬레이션 환경 설정은 ABC 적용 3-D FEM과 같다.

(그림 13) 도파관 포트 경계조건 적용 위치(d)에 따른 S₁₁ 파라미터 및 Reference와의 결과 비교

(Figure 13) Comparison of S₁₁ parameter (Reference vs. In-house) according to distance of WPBC application(d)

(그림 14) 도파관 포트 경계조건 적용 위치(d)에 따른 S₁₂ 파라미터 및 Reference와의 결과 비교

(Figure 14) Comparison of S₁₂ parameter (Reference vs. In-house) according to distance of WPBC application(d)

d=0.8b에서는 TE₁₀모드, d=0.4b 및 d=0.2b에서는 TE₁₀, TE₀₁, TM₁₁, TE₁₁모드를 흡수하도록 설정하였다. 시뮬레이션 결과 경계조건 적용 위치(d)가 불연속 물체와 가까워질수록 정확도가 떨어지는 ABC와 달리 WPBC는 고차모드를 흡수해 더욱 좋은 정확도를 보임을 알 수 있다.

4. 결론

본 논문에서는 FEM 기법을 사용하여 2-D 및 3-D waveguide 구조를 해석할 때, 계산영역을 한정시키기 위한 경계조건인 ABC 및 WPBC의 특징 및 적용 예에 대해 알아보았다. 각각의 경우를 직접 완결된 코드로 작성하여 시뮬레이션하였으며, 참고문헌[11] 및 상용 소프트웨어 HFSS와 결과를 reference로 하여 비교하여 그 정확도를 확인하였다.

2-D 및 3-D 모두 ABC 및 WPBC 적용 거리(산란체로부터의 거리)에 따른 성능을 테스트하였으며, 시뮬레이션 결과에서 확인할 수 있듯이 내부 유전체로부터 떨어진 거리가 충분히 먼 경우 ABC와 WPBC 모두 좋은 정확도를 보였으나 거리가 가까워질수록 ABC의 정확도는 떨어졌지만 WPBC는 여전히 좋은 정확도를 보이는 것을 확인할 수 있다.

다만, 경계조건 적용 식에서 ABC는 γ가 상수로 적용되지만, WPBC는 \(\begin{aligned}\gamma(\vec{E})\\\end{aligned}\), \(\begin{aligned}P(\vec{E})\\\end{aligned}\)와 같이 operator 형태로 적용되기 때문에 유한요소 과정을 거쳐 최종적으로 얻게 되는 행렬이 WPBC 적용 경우가 ABC 적용 경우보다 0이 아닌 요소를 많이 포함하는 형태가 된다. 따라서 WPBC를 사용하면 해석 영역이 줄어들 수 있지만, 메모리 사용량과 계산 복잡도가 증가하며 이는 계산 시간의 증가로 이어지기 때문에 시뮬레이션 상황에 따라 적절하게 사용해야 할 것으로 보인다.

WPBC의 해석 시간 증가는 병렬화 기법 등 다양한 방법으로 해결할 수 있을 것으로 예상되며, PML 등 다른 경계조건과 결합하여 사용한다면 더 나은 성능 향상을 기대할 수 있을 것으로 보인다.