Absorbing Boundary Conditions and Parallelization for Waveguide Electromagnetic Analysis Using Finite Element Method

유한요소법을 이용한 도파관 전자기 해석의 흡수경계조건 고찰 및 병렬화

Abstract

Power and signal transmission using electromagnetic waves are essential in modern times, and a guided structure is needed to transmit electromagnetic waves efficiently through the desired path. This paper performed an electromagnetic simulation using the in-house code for the 2-D/3-D waveguide using the finite element method. The accuracy of the analysis was verified by comparing it with the results of HFSS, a representative electromagnetic wave simulation software. In addition, the performance of the Absorbing Boundary Condition (ABC), which is essential to truncate the infinite computational domain for computational electromagnetics, was analyzed. Finally, the parallelization technique was applied to accelerate the simulation speed, demonstrating performance improvement.

현대에는 전자기파를 이용한 전력 및 신호 전달이 필수적인데, 전자기파를 원하는 경로를 통해 효율적으로 전달하기 위해서는 도파 구조(guided structure)가 필요하다. 본 논문에서는 먼저 전파해석 기법인 유한요소법(FEM : Finite Element Method)을 적용하여 도파 구조 중 하나인 2-D/3-D 도파관(waveguide)에 대해 직접 in-house code를 작성하여 전자기 시뮬레이션하였다. 이후 in-house code의 해석 결과를 대표적인 전자파 상용 시뮬레이션 소프트웨어인 HFSS의 결과와 비교하여 해석의 정확성을 검증하였다. 아울러, 전자기 해석에 있어 무한대의 해석 영역을 잘라 해석하기 위해 필수적인 흡수경계조건(ABC : Absorbing Boundary Condition)의 성능을 분석한 후, 병렬화 기법의 적용에 따른 성능 향상을 제시하였다.

1. 서론

전자기파를 특정 경로를 통해 효율적으로 전송하기 위한 특별한 구조 및 장치를 통칭하여 도파 구조(guided structure)라고 하며, 도파 구조는 크게 전송선(transmission line)과 도파관(waveguide)의 두 가지로 분류할 수 있다. 이러한 도파 구조 중 하나인 도체로 둘러 쌓인 속이 빈관, 즉 구형 도파관(rectangular waveguide) 구조는 고주파 대역에서 상대적으로 감쇠를 줄일 수 있어 해당 대역에서 효율적 에너지 전송이 가능한 장점이 있다[1].

전자기 구조의 해석을 위해서는 구조별 상황에 따라 부여된 Maxwell 방정식을 풀어야 하는데, 이 방정식의 해석적인 해(analytical solution)를 직접 구할 수 있는 경우는 통상적으로 강한 대칭성이 확보된 매우 특정한 경우에 한정된다. 따라서, 현실적인 문제 해결을 위해서는 수치적 시뮬레이션 기법, 즉 전자기 수치해석 또는 수치해석 전자기학(computational electromagnetics)의 도입이 필요하다.

대표적인 전자기 수치해석 기법으로 유한 차분 시간 영역 방법(FDTD : Finite Difference Time Domain), 모멘트 방법(MoM : Method of Moment), 유한요소법(FEM : Finite Element Method) 등이 있다. FDTD는 시간 영역에서 미분방정식을 사용해 맥스웰 방정식을 해석하는 방법으로, 매질의 이방성과 불균일성을 용이하게 다룰 수 있으며 메시(mesh)의 생성이 간단하지만, 곡면 등을 모델링하기 어렵다고 알려져 있다[2]. 이와 더불어, MoM은 그린 함수를 사용한 적분 방정식을 사용해 맥스웰 방정식을 해석하는 방법으로, 계산영역이 표면으로 한정되므로 원거리 장 해석 시 큰 영역의 해석이 필요하지 않다[2]. 그러나, MoM은 비균질성 매질의 해석이 어려우며, full matrix 를 생성하기 때문에 계산이 복잡하다는 단점도 있다[2].

앞서 기술된 기법들에 반해 FEM은 미분 방정식을 기반으로 한 전자기 해석 방법의 일종으로 해석의 영역을 미소한 요소(element)로 나누는 것을 기본으로 한다. 이기 법은 다른 수치적 방법들보다 다소 수학적으로 복잡하다고 알려져 있지만, 삼각형·사면체·곡선과 같은 메시를 사용해 복잡한 구조 및 곡면의 모델링을 할 수 있고, 풀이의 대상이 되는 행렬이 희소 행렬(sparse matrix)이 되어, 희소 행렬에 적합한 반복 솔버(iterative solver)를 사용해 효율적으로 해(solution)를 얻을 수도 있다[3]. 또한, 영역 분할 알고리듬을 적용하여 병렬 연산에도 적합하다는 장점이 있다[3].

본 논문은 FEM의 노드(node) element를 사용해 parallel-plate waveguide의 자기장을 해석하는 2-D 예제로부터 시작한다. 이후, node element만을 사용하여 전자기해석 시 물리적으로 가능하지 않은 spurious solution이 발생하는 경우가 있기 때문에[4-6], 이를 해결하기 위한 에지(edge) element를 추가로 사용하여 insulated waveguide 의 전파상수(k_z)를 해석한 경우에 대하여 논한다. 아울러, 2-D 구조에서 확장한 3-D 예제에서 3-D edge element를 사용한 inhomogeneous waveguide의 전파 반사계수, 즉 S₁₁ parameter에 대하여 논의한다.

각 예제에서는 경계면에서의 반사를 최소화하기 위한 흡수경계조건(ABC : Absorbing Boundary Condition)을 적용하게 되며[6], 2-D insulated waveguide 예제에서는 first-order ABC와 second-order ABC의 특징 및 성능 차이의 검증을 진행하였고, 3-D inhomogeneous waveguide 예제의 경우 first-order ABC 설치 위치에 따른 전파 흡수성능의 비교 및 검증을 진행하였다. 마지막으로, MATLAB의 Parallel Computing Tool box[7]를 사용한 2-D 및 3-D 구조 해석의 병렬화 진행 및 그에 따른 성능 향상에 대하여 살펴보았다.

아울러, 본 논문은 2021 한국인터넷정보학회 추계학술발표대회 논문집 제21권 제2호에 실린 ‘Edge Element를 이용한 도파관 유한요소 해석의 병렬화’ 논문[8]의 확장판임을 밝힌다.

2. 배경 이론 및 해석 방법

2.1 2-D FEM 공식화

2.1.1 노드 기반 FEM

그림 1은 가운데 주황색 박스로 표시된 불연속 매질이 존재하는 평행판 도파관(parallel-plate waveguide)에서, 왼쪽에서 오른쪽으로 전파하는 전자기파가 입사된 상황을 나타내며, 이러한 경우에서 불연속 매질의 유전율에 따른 자기장의 변화를 해석하는 문제를 다룰 것이다. 앞서 언급했듯이 이러한 문제에 있어 해석해를 찾기는 매우 어려우므로, FEM과 같은 수치적인 방법을 통해 해를 구할 수밖에 없다.

(그림 1) 불연속 매질이 포함된 평행판 도파관의 구조[6]

(Figure 1) Structure of a parallel-plate waveguide with inhomogeneity[6]

먼저, 2-D node 기반 삼각형 메시(triangular mesh)를 사용하여 해석하고자 하는 영역을 미소한 삼각형으로 나누었으며, 이에 FEM의 공식화 절차를 사용하여 선형대수 방정식을 완성하여 해석하였다. 이 구조 해석의 출발점이 되는 지배방정식은 다음의 식(1)과 같다[6].

\(\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=0\)       (1)

이를 유한요소 해석을 위한 기본 플랫폼이 되는 하기식(2)의 꼴에 식(1)의 변수를 대응시키면, 식(2)의 각 변수는 식(3)과 같이 정의된다.

\(-\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\)       (2)

\(\phi=H_z, \alpha_x=\alpha_y=\frac{1}{\epsilon_r}, \beta=-k_0^2 \mu_r, f=0,\)       (3)

아울러, variational principle 적용을 위한 functional I의 형태는 아래 식(4)와 같이 정리될 수 있다[6].

\(\begin{aligned} I(\phi)=\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}\)       (4)

최종적으로 variational principle과 Ritz process를 적용하면 식(5)의 꼴이 도출되며, 여기서 각 element에 적용되는 elemental matrix는 식(6), 식(7)과 같이 표현된다[6]. 최종적으로 matrix assembly 과정을 거쳐 Ax=b 꼴의 선형대수 방정식을 풀면 본 구조에서 원하는 해인 자기장 (Hz)을 구할 수 있다.

\(\left\{\frac{\partial I^e}{\partial \phi^e}\right\}=\left[K^e\right]\left\{\phi^e\right\}-\left\{b^e\right\}\)       (5)

\(\begin{aligned} & K_{i j}^e=\iint_{\Omega}\left(\begin{array}{l} \alpha_x \frac{\partial N_i^e}{\partial x} \frac{\partial N_j^e}{\partial x}+\alpha_y \frac{\partial N_i^e}{\partial y} \frac{\partial N_j^e}{\partial y} \\ +\beta N_i^e N_j^e \end{array}\right) d x d y \\ & i, j=1,2,3 \end{aligned}\)       (6)

\(b_i^e=\iint_{\Omega^e} f N_i^e d x d y \quad i=1,2\)       (7)

2.1.2 2-D 구조에 적용되는 흡수경계조건

그림 1 구조에서의 전파의 진행에 있어, 구조의 왼쪽 경계는 입사파와 반사파의 합으로 이루어진 식(8)과 같이 나타낼 수 있고, 오른쪽 경계는 투과파로만 이루어진식 식(9)와 같이 나타낼 수 있다. 이 경우, 각 경계로부터 계산영역을 잘라내면서 생기는 비물리적 반사를 최소화하기 위해 흡수경계조건(ABC)을 도입하게 된다[6]. 실제 이러한 흡수경계조건을 적용하지 않는다면 전파 진행의 거동이 제대로 모사되지 못하여 물리적으로 합당한 결과를 얻을 수 없게 된다.

\(H_z=H_z^{i n c}+H_z^{r e f}=H_0 e^{-j k_0 x}+R H_0 e^{j k_0 x}\)      (8)

\(H_z=H_z^{\text {trans }}=T H_0 e^{-j k_0 x}\)     (9)

First-order ABC는 가장 간단한 형태의 흡수 경계조건으로, 식(10), 식(11)의 식의 형태로 정리될 수 있다[6, 9]. 또한, Third kind boundary condition을 모사한 식(12) 형식으로 식(10) 및 식(11)을 정리하게 되면 각 변수는 왼쪽, 오른쪽 경계에서 각각 식(13), 식(14)와 같이 정의될 수 있고 이러한 변환을 통해 통상적인 FEM 절차에 흡수 경계조건을 포함시킬 수 있다[6, 9].

\(\frac{\partial H_z}{\partial x}=j k_0 H_z-2 j k_0 H_0 e^{-j k_0 x}\)      (10)

\(\frac{\partial H_z}{\partial x}=-j k_0 H_z\)       (11)

\(\left(\alpha_x \frac{\partial \phi}{\partial x}+\alpha_y \frac{\partial \phi}{\partial y}\right) \cdot \hat{n}+\gamma \phi=q\)       (12)

\(\gamma=\alpha_x j k_0, q=\alpha_x 2 j k_0 H_0\)       (13)

\(\gamma=\alpha_x j k_0, q=0\)       (14)

하지만, first-order ABC는 0º입사각(즉, 수직 입사) 에대해서만 좋은 전파 흡수를 보장하기 때문에, 다른 입사각에서도 준수한 흡수 성능을 보장하기 위해서는 second-order ABC를 사용할 수 있다[6, 10].

식(15)은 2-D second-order ABC(식(10)의 first-order ABC와 차이가 있음을 유의)이며, 이를 문제에 적용하기 위해서는 식(13) 및 식(14)에 해당하는 를 식(16)과 같이 변형해 주어야 한다[6, 10]. 최종적으로 경계면에서의 functional I_b 식(17)을 variational principle 및 Ritz process 를 통해 식(18)과 같은 형태로 만들어주게 되면서 FEM 풀이 과정을 완성하게 되며, 여기에서 사용되는 보간 함수(0과 1 사이로 normalization 된 형태)는 식(19)의 형태로 표현되며, 이에 대한 elemental matrix는 각각 식(20), 식(21)과 같다[6].

\(\frac{\partial \phi}{\partial x}=-j k_0 H_z-\frac{1}{2 k_0} \frac{\partial^2 \phi}{\partial y^2}\)       (15)

\(\gamma=\gamma_1+\gamma_2, \gamma_1=j k_0, \gamma_2=\frac{1}{2 k_0} \frac{\partial^2}{\partial y^2}\)       (16)

\(I_b(\phi)=\int_{\Gamma_2}\left(\frac{\gamma}{2} \phi^2-q \phi\right) d \Gamma\)       (17)

\(\left\{\frac{\partial I_b^s}{\partial \phi^s}\right\}=\left[K^s\right]\left\{\phi^s\right\}-\left\{b^s\right\}\)       (18)

\(N_1^s=1-\xi, \quad N_2^s=\xi\)       (19)

\(K_{i j}^s=\int_0^1 \gamma N_i^s N_j^s l^s d \xi \quad i, j=1,2\)       (20)

\(b_i^s=\int_0^1 q N_i^s l^s d \xi \quad i=1,2\)       (21)

2.1.3 노드, 에지 동시 사용 FEM

본 절에서는 2-D edge element(삼각형 메시 사용)와 node element를 동시에 사용하여, 그림 2와 같은 insulated waveguide 구조에 대한 전파상수(k_z)를 구하는 문제를 기술한다. Node element만을 사용했을 때 발생하는 divergence condition의 불연속으로 인한 spurious solution 이 발생할 수 있는 문제는 edge element를 추가로 도입함으로써 해결될 수 있다[10, 11].

(그림 2) 절연 도파관의 구조(case 1)[6]

(Figure 2) Structure of insulated waveguide (case 1)[6]

(h=0.02m, w/h=2.25, d/h=0.5, a/h=13.5, b/h=8.0, =3.8, =1.5)

또한, 본 구조를 해석하여 전파상수(k_z)를 도출하기 위해서는 waveguide 구조에 대한 고유치 문제를 풀어야 하며, 이 방정식의 도출은 식(22)과 같은 형태의 vector Helmholtz equation으로부터 시작된다[8, 10, 11].

\(\nabla \times\left(\frac{1}{\mu_r} \nabla \times \vec{E}\right)-k_0^2 \epsilon_r \vec{E}=0\)       (22)

이후, waveguide에 입사되는 전기장에 variational principle을 적용하면 식(23)과 같은 functional 형태가 도출된다[5, 9, 11].

\(\begin{array}{r} I(\vec{E})=\frac{1}{2} \iint_{\Omega}\left[\frac{1}{\mu_r}\left(\nabla_t \times \overrightarrow{E_t}\right) \cdot\left(\nabla_t \times \vec{E}_t\right)^*-k_0^2 \epsilon_r \vec{E} \cdot \vec{E}\right. \\ \left.+\frac{1}{\mu_r}\left(\nabla_t E_z \times j k_z \overrightarrow{E_t}\right) \cdot\left(\nabla_t E_z \times j k_z \vec{E}_t\right)^*\right] d \Omega \end{array}\)       (23)

최종적으로 주파수 및 전파상수(k_z)에 대한 고유치 문제로 바꾸기 위해 식(23)의 변수를 식(24)의 형태로 변환하여 식(25)의 형태의 functional을 도출하게 되며, 이를 행렬 형태로 표시하면 식(26)과 같아진다[5, 8, 9, 11].

\(\overrightarrow{e_t}=k_z \vec{E}_t, e_z=-j E_z\)       (24)

\(\begin{aligned} I(\vec{e})=\frac{1}{2} \iint_{\Omega} & \left(\frac{1}{\mu_r}\left(\nabla_t \times \overrightarrow{e_t}\right) \cdot\left(\nabla_t \times \overrightarrow{e_t}\right)^*-k_0^2 \epsilon_r \overrightarrow{e_t} \cdot \overrightarrow{e_t}\right. \\ & +k_z^2\left[\frac{1}{\mu_r}\left(\nabla_t e_z \times j k_z \overrightarrow{e_t}\right) \cdot\left(\nabla_t e_z \times j k_z \vec{e}_t\right)^*\right. \\ & \left.\left.-k_0^2 \epsilon_r e_z \cdot e_z^*\right]\right) d \Omega \end{aligned}\)       (25)

\(I=\frac{1}{2} \sum_{e=1}^M\left(\left\{e_t\right\}^T\left[A_{t t}^e\right]\left\{e_t\right\}^*+\frac{1}{2} k_z^2\left[\begin{array}{l} e_t^e \\ e_z^e \end{array}\right]^T\left[\begin{array}{ll} B_{t t}^e & B_{t z}^e \\ B_{z t}^e & B_{z z}^e \end{array}\right]\left[\begin{array}{l} e_t^e \\ e_z^e \end{array}\right]^*\right)\)       (26)

아울러, 식(26)의 element matrix는 식(27)과 같이 표현된다[6, 8, 9, 11-13].

\(\begin{aligned} & {\left[A_{t t}^e\right]=} \iint_{\Omega^e}\left[\frac{1}{\mu_r^e}\left(\nabla_t \times \overrightarrow{N^e}\right) \cdot\left(\nabla_t \times \overrightarrow{N^e}\right)^T\right. \\ &\left.\quad-k_0^2 \epsilon_r^e\left(\overrightarrow{N^e}\right) \cdot\left(\overrightarrow{N^e}\right)\right] d \Omega \\ & {\left[B_{t t}^e\right]=\iint_{\Omega^e} \frac{1}{\mu_r^e}\left(\overrightarrow{N^e}\right) \cdot\left(\overrightarrow{N^e}\right)^T d \Omega } \\ & {\left.\left[B_{t z}^e\right]=\iint_{\Omega^e} \frac{1}{\mu_r^e} \overrightarrow{\left(N^e\right.}\right) \cdot\left(\nabla_t N^e\right)^T d \Omega } \\ & {\left[B_{z t}^e\right]=\iint_{\Omega^e} \frac{1}{\mu^e}\left(\nabla_t N^e\right) \cdot\left(\overrightarrow{N^e}\right)^T d \Omega } \\ & {\left[B_{z z}^e\right]=\iint_{\Omega^e}\left[\frac{1}{\mu^e}\left(\nabla_t N^e\right) \cdot\left(\nabla_t N^e\right)^T-k_0^2 \epsilon_r^e\left(N^e\right)\left(N^e\right)^T\right] d \Omega } \end{aligned}\)       (27)

최종적으로 본 waveguide 구조의 전파상수(k_z)를 구하는 문제는 식(28)과 같이 고윳값()을 구하는 문제로 대치되어 해석할 수 있다[8, 13].

\(\left[\begin{array}{rr} A_{t t} & 0 \\ 0 & 0 \end{array}\right]\left[\begin{array}{l} e_t \\ e_z \end{array}\right]=-k_z^2\left[\begin{array}{ll} B_{t t} & B_{t z} \\ B_{z t} & B_{z z} \end{array}\right]\left[\begin{array}{c} e_t \\ e_z \end{array}\right] \Leftrightarrow A x=\lambda B x\)      (28)

2.2 3-D FEM 공식화

2.2.1 에지 기반 FEM

3-D FEM 예제의 경우 식(29)와 같은 경계값 문제를 기반으로, 주모드인 TE₁₀ Mode만 감쇠 없이 진행하는 불균일 매질이 포함된 구형도파관(rectangular waveguide) 구조(그림 3)에 대한 반사 계수 S₁₁ parameter를 구하게 되며, 해석을 진행하기 위해 3-D Edge(brick) element를 적용하였다[6, 14-15].

\(-\frac{\partial}{\partial x}\left(\alpha_x \frac{\partial \vec{E}}{\partial x}\right)-\frac{\partial}{\partial y}\left(\alpha_y \frac{\partial \vec{E}}{\partial y}\right)-\frac{\partial}{\partial z}\left(\alpha_z \frac{\partial \vec{E}}{\partial z}\right)+\beta \vec{E}=f\)       (29)

(그림 3) 불균일 도파관의 구조

(Figure 3) Structure of inhomogeneous waveguide (b = 10mm, )

아울러, waveguide에 입사하는 전기장을 functional의 형태로 나타나게 되면 식(30)과 같다[6, 9, 15].

\(\begin{aligned} I(\vec{E})= & \frac{1}{2} \iiint_V\left[\frac{1}{\mu}(\nabla \times \vec{E}) \cdot(\nabla \times \vec{E})-k_0^2 \epsilon_r \vec{E} \cdot \vec{E}\right] d V \\ & +\iint_{S_1}\left[\frac{\gamma}{2}(\hat{n} \times \vec{E}) \cdot(\hat{n} \times \vec{E})+\vec{E} \cdot \vec{U}^{i n c}\right] d S \\ & +\iint_{S_2}\left[\frac{\gamma}{2}(\hat{n} \times \vec{E}) \cdot(\hat{n} \times \vec{E})\right] d S \end{aligned}\)       (30)

식(30)은 식(31)과 같이 간략히 표현될 수 있으며, 식 (31)의 각 elemental matrix는 식(32)과 같다[6, 15]. 차후 FEM의 절차를 적용해 선형대수 방정식을 풀게 되면 최종적으로 원하는 전기장()을 구할 수 있으며, 이를 최종적으로 식(33)에 대입하여 후처리하면 반사계수인 R(즉, S₁₁ parameter)을 구할 수 있다[6].

\(\begin{aligned} I= & \frac{1}{2} \sum_{e=1}^M\left\{E^e\right\}^T\left[K^c\right]\left\{E^e\right\}+\frac{1}{2} \sum_{s=1}^{M_s}\left\{E^s\right\}^T\left[B^s\right]\left\{E^s\right\} \\ & -\frac{1}{2} \sum_{s=1}^{M_{s 1}}\left\{E^s\right\}\left\{b^s\right\} \end{aligned}\)       (31)

\(\begin{aligned} & {\left[K^c\right]=\iiint_{V^c}\left[\frac{1}{\mu_r^e}\left\{\nabla \times \vec{N}^e\right\} \cdot\{\nabla \times \vec{N}\}^T\right.} \\ & \left.\quad-k_0^2 \epsilon_r^e\{\vec{N}\} \cdot\left\{\vec{N}^{\prime}\right\}^T\right] d V \\ & {\left[B^s\right]=\iint_{S^s} \gamma\left\{\vec{S}^s\right\} \cdot\left\{\vec{S}^s\right\}^T d S} \\ & \left\{b^s\right\}=\iint_{S^s}\{\hat{n} \times \vec{S}\} \cdot \vec{U}^{n c} d S \end{aligned}\)      (32)

\(R=\frac{2 e^{-j k_{z 10} z_1}}{a b E_0} \iint_{S_1} \vec{E}\left(x, y, z_1\right) \cdot \vec{e}_{10}(x, y) d S-e^{-2 j k_{210} z_1}\)      (33)

2.2.2 3-D 흡수경계조건

TE₁₀ Mode의 경우 XY 평면에 입사하는 전파는 입사파와 반사파로 이루어져 있으며, 이는 식(34)과 같이 나타낼 수 있고, 반대쪽 XY 평면으로 나가는 투과파는 식 (35)과 같이 나타낼 수 있다[6, 15-16].

\(\begin{aligned} & \vec{E}(x, y, z)=\vec{E}^{i n c}(x, y, z)+\vec{E}^{r e f}(x, y, z) \\ & =E_0 e_{10}(x, y) e^{-j k_{100^z}}+R E_0 \vec{e}_{10}(x, y) e^{j k_{110^z}} \end{aligned}\)       (34)

\(\vec{E}(x, y, z)=\vec{E}^{\text {trans }}(x, y, z)=T E_0 \vec{e}_{10}(x, y) e^{-j k_{10 n^z}}\)       (35)

또한, 식(34), 식(35)는 간단한 연산의 적용을 통해 식 (36), 식(37)의 형태로 변형되며, 이는 3-D first-order ABC 의 형태가 된다[6, 15-16].

\(\begin{aligned} & \hat{n} \times(\nabla \times \vec{E})=\hat{z} \times(\nabla \times \vec{E}) \\ & =-j k_{z 10} \vec{E}^{n c}+j k_{z 10} \vec{E}^{e f} \\ & =j k_{z 10} \vec{E}-2 j k_{z 10} \overrightarrow{\vec{E}}^{n c} \end{aligned}\)       (36)

\(\hat{n} \times(\nabla \times \vec{E})+\gamma \hat{n} \times(\hat{n} \times \vec{E})=0\)       (37)

3. 시뮬레이션 결과

3.1 2-D 시뮬레이션 결과

3.1.1 ABC를 적용한 노드 기반 FEM 예

2-D 노드 기반 FEM을 사용한 parallel-plate waveguide 의 시뮬레이션은 상용 시뮬레이션 소프트웨어가 아닌 직접 작성한 In-house code(MATLAB R2021a 버전)를 이용하여 진행하였으며, 무손실 매질에서 성립하는 식(38)을 사용해 본 시뮬레이션의 정확성을 검증할 수 있다.

\(\left|R|^2+\right| T|^2=1\)       (38)

또한, 식(38)의 반사 계수(R)와 투과계수(T)는 각각 식 (39), 식(40)과 같이 표현된다[6].

\(R=\frac{H_z\left(x_1\right)-H_0 e^{-j k_0 x_1}}{H_0 e^{j k_0 x_1}}\)       (39)

\(T=\frac{H_z\left(x_2\right)}{H_0 e^{-j k_0 x_2}}\)       (40)

시뮬레이션 환경주파수는 λ₀를 0.1m로 설정하기 위 해약 3GHz로 설정하였고, 그림 4는 first-order ABC를 적용해 자기장(H₂)을 해석한 결과이며, 그림 5는 같은 구조에 대해 second-order ABC를 적용해 해석한 결과이다. 두 결과 모두 참고문헌[6]의 도시와 거의 일치한다.

(그림 4) 일차 흡수경계조건 기반 해석 일 때의 자기장 분포

(Figure 4) First-order ABC based analysis magnetic field distribution at

(그림 5) 이차 흡수경계조건 기반 해석 일 때의 자기장 분포

(Figure 5) Second-order ABC based analysis magnetic field distribution at

First-order ABC와 second-order ABC 사용해 해석한 결과의 식(38)의 값은 각각 1.0012, 1.0005로 시뮬레이션 결과에 대한 타당성이 검증되었으며, 입사각에 대해 유리한 second-order ABC를 사용했을 때 좀 더 이론값인 1에 가까운 결과임을 알 수 있다.

3.1.2 노드, 에지 동시 사용 FEM 적용 예

2-D 노드, 에지 기반 FEM을 사용한 insulated waveguide (그림 2)의 in-house code는 MATLAB R2021a 환경에서 실행되었으며, 참고문헌[6]과의 결과 비교를 통해 시뮬레이션의 정확성을 검증하였다. 또한, 시뮬레이션 환경주파수는 1~8GHz로 설정하였다.

시뮬레이션 결과는 주파수(k₀)에 대한 전파상수(k_z)를각각 k₀*h와 k_z/k₀의 값으로 정규화 후 해석 결과를 참고문헌[6]의 결과와 비교하여 정확성을 검증하였다[8]. 그림 6의 경우 그림 2 구조에 대해 전체적으로 해석한 결과이며, 기존의 결과와 잘 일치함이 확인된다.

(그림 6) 그림 2 도파관의 분산 특성[6] (전체 구조 해석)

(Figure 6) Dispersion characteristics of waveguide[6] (Full structural analysis)

3.2 3-D 시뮬레이션 결과

3-D 에지 기반 FEM을 사용한 inhomogeneous waveguide의 시뮬레이션은 MATLAB R2021a 버전을 기반으로 실행되었으며, in-house code와 상용 소프트웨어인 HFSS와의 결과 비교를 통해 시뮬레이션의 정확성을 검증하였다.

다시 한번, 그림 3은 해석에 사용된 waveguide 구조로 내부 유전체의 크기는 8.0 x 6.0 x 10.0mm이다. 또한, planewave가 XY 평면에서 입사되고, 반대쪽 XY 평면으로 전파하는 상황에 대해 해석을 하였으며, 시뮬레이션 환경주파수 4~7GHz에 대해 0.01GHz 간격으로 waveguide 의 전체 E-field를 구한 후 후처리 과정을 통해 입사파에 대한 반사파의 비로 표현되는 S₁₁ parameter를 얻게 된다.

또한, 흡수경계조건(ABC)은 산란체 즉, 본 구조의 내부 유전체에서 멀리 설치할수록 흡수 성능이 더 좋아지는데(비물리적 반사가 더 적음), 상기 그림 3의 waveguide 구조에서 산란체와의 거리 d(b의 길이로 정규화)를 변화시켜가며 그 성능을 측정하였다. 그림 7의 결과에서 확인할 수 있듯이 d가 커질수록, 다시 말해 ABC가 산란 체에서 멀리 설치될수록 그래프의 경향이 상용 소프트웨어인 HFSS의 결과와 비슷해지는 것이 확인된다. 그림 8은 거리별, 주파수 포인트에 따른 HFSS 결과(y_HFSS)와 in-house 코드의 시뮬레이션 결과(y_in-house)의 절대 오차 그래프이며, 이를 식(43)를 사용해 평균 제곱근 편차 (RMSE : Root Mean Square Error)로 나타내면 표 1과 같다. 결론적으로, ABC가 가장 멀리 설치되었을 때 HFSS 의 결과와 가장 적은 차이(가장 적은 반사)를 보임을 알 수 있다.

\(\mathrm{RMSE}=\sqrt{\sum_{\mathrm{i}=1}^{\mathrm{n}} \frac{\left(y_{\text {in-house }}-y_{H F S S}\right)^2}{\mathrm{n}_{\text {frequency }}}}\)        (43)

(그림 7) HFSS와 흡수경계조건 적용 위치에 따른 S₁₁ 파라미터 비교

(Figure 7) Comparison of S₁₁ parameter (HFSS vs. In-house)

(그림 8) HFSS와 흡수경계조건 적용 위치에 따른 S₁₁ 파라미터의 절대 오차

(Figure 8) Absolute error (HFSS vs. In-house)

(표 1) d에 따른 RMSE

(Table 1) RMSE value along the length of d

4. FEM 해석의 병렬화

2-D 노드/에지 기반 insulated waveguide 구조(그림 2) 와 3-D 에지 기반 inhomogeneous waveguide 구조(그림 3) 의 해석 속도 향상을 위해 MATLAB의 Parallel Computing Toolbox[7]의 parfor 문으로 병렬화를 적용하여 추가 시뮬레이션하였으며, 기존 병렬화 기법을 적용하지 않은 코드와 속도를 비교하였다. 이때, 각 frequency sweep은 그림 9와 같이 병렬적으로 처리할 수 있고, 사용된 thread의 개수는 6개이다.

표 2는 그림 10에 제시된 2-D insulated waveguide(그림 2) 해석에 필요한 실행시간을 구체적으로 나타낸 것이며, 병렬화를 적용할 경우 해석 속도가 기존보다 약 2-3배 가량 빠른 것을 확인할 수 있다[8]. 한편, 표 3은 그림 11에 제시된 3-D inhomogeneous waveguide(그림 3) 해석에 필요한 실행시간을 나타낸 것이며, 2차원 구조 해석과 같은 환경(thread 개수 등)하에서 시뮬레이션하였기 때문에 비슷한 성능 향상률을 보임을 알 수 있다.

(그림 9) 주파수 스윕에 대한 병렬처리 과정

(Figure 9) Parallel processing process for frequency sweep

(그림 10) 2-D parallel-plate waveguide에서 병렬화 적용에 따른 실행시간 비교

(Figure 10) Comparison of execution time according to serial/parallel analysis technique for 2-D parallel-plate waveguide

(표 2) 요소 개수에 따른 병렬처리 성능 향상 (2-D 예제)

(Table 2) Comparison of parallel processing performance improvement according to the number of elements

(그림 11) 3-D inhomogeneous waveguide에서 병렬

(Figure 11) Comparison of execution time according to serial/parallel analysis technique for 3-D inhomogeneous waveguide

(표 3) 요소 개수에 따른 병렬처리 성능 향상 (3-D 예제)

(Table 3) Comparison of parallel processing performance improvement according to the number of elements

5. 결론

본 논문에서는 다양한 2-D 및 3-D Waveguide 예제를 다양한 FEM 기법을 적용해 해석하였다. 아울러, 2-D 예제에서는 waveguide에 입사되는 전자기파에 대한 반사를 최소화하기 위한 ABC의 종류(first-order, second-order) 와특징에 알아보았으며, 각각의 경우를 직접 코딩하여 시뮬레이션하였다. 아울러, 3-D 예제에서는 대표적인 전자기 해석 시뮬레이션 소프트웨어 중 하나인 HFSS를 사용한 시뮬레이션 결과와 직접 작성한 In-house code의 비교검증을 통해 ABC 설치 거리에 따른 성능을 테스트하였다. 3-D 시뮬레이션 결과에서 확인할 수 있듯이, ABC 설치 거리는 산란체에서 멀어질수록 적은 반사를 보였으나, ABC를 산란체에서 멀리 설치할수록 해석 영역이 늘어나 해석시간이 증가하므로, 정확도와 해석시간의 적절한 조율이 필요할 것으로 보인다.

이러한 해석시간의 증가는 frequency sweep 병렬화 기법을 통해 해결할 수 있으며, 본 논문에서는 MATLAB 병렬화 기법의 적용을 통해 해석 속도의 향상을 확인하였다. 현재 병렬 thread의 지원 개수 제한으로 성능 향상에 제약점이 있지만, 더 많은 thread를 지원하는 여타 기법(CUDA 등 채택)의 적용을 통한 성능 향상을 기대할 수 있다.

본 연구는 방위사업청과 국방과학연구소가 지원하는 스텔스 대형 플랫폼 전파해석 특화연구실 사업의 일환으로 수행되 었음(UD200047JD).

본 논문은 2021년도 한국인터넷정보학회 추계학술대회 우수 논문 추천에 따라 확장 및 수정된 논문임.