FDTD

1. 개요

시간영역 유한차분법 Finite-Difference Time-Domain
약칭	FDTD
제안	케인 이(Kane S. Yee), 1966
지배 방정식	맥스웰 방정식 (시간영역)
기반 기법	유한차분법
안정 조건	쿠랑 한계 (CFL 조건) — 3차원 균일격자에서 S ≤ 1/√3
정확도	공간·시간 모두 2차 (중심차분)
대표 경계조건	완전정합층(PML, Berenger 1994)

펄스 한 방 쏘고, 나머지는 컴퓨터가 알아서 한다.

FDTD(Finite-Difference Time-Domain, 시간영역 유한차분법)는 맥스웰 방정식의 회전(curl) 방정식을 시간과 공간 모두에서 차분으로 근사하여, 전자기장이 시간에 따라 어떻게 퍼져나가는지를 직접 계산하는 유한차분법 기반의 수치 기법이다. 1966년 케인 이(Kane S. Yee)가 제안한 격자 배치 아이디어에서 출발했으며¹, 그래서 지금도 그 격자를 Yee 격자라 부른다. 실제로 산란 문제에 적용해 안정 조건까지 정리한 것은 1975년 태플로브와 브로드윈이고, FDTD라는 약칭도 그 계열에서 굳었다.

주파수 영역 기법이 특정 주파수 하나하나를 따로 푸는 것과 달리, FDTD는 시간 파형을 그대로 굴린다. 그래서 넓은 대역의 응답을 시뮬레이션 한 번으로 뽑아낼 수 있고, 전자기파가 구조물을 통과하고 반사되는 과정을 애니메이션처럼 눈으로 볼 수 있다는 것이 큰 매력이다. 주파수영역 쪽 대표 선수인 모멘트법·유한요소법과는 경쟁이 아니라 역할 분담 관계다.

2. Yee 격자 — 엇갈림의 미학

FDTD의 핵심 발명품은 전기장 $\mathbf{E}$와 자기장 $\mathbf{H}$를 같은 지점에 두지 않고 공간적으로 반 칸 엇갈리게(staggered) 배치하는 Yee 격자다. 전기장 성분은 격자 셀의 모서리(edge) 중앙에, 자기장 성분은 면(face) 중앙에 놓이는 식으로, 각 전기장을 자기장 성분들이 고리처럼 둘러싸고 그 반대도 성립한다.

이 배치가 왜 천재적인가 하면, 맥스웰 방정식의 회전 연산자가 요구하는 “어떤 장의 회전이 옆에 있는 다른 장을 만든다”는 구조를 격자가 그대로 본뜨기 때문이다. 덕분에 중심 차분(central difference)이 자연스럽게 성립해 2차 정확도를 얻고, 발산 조건($\nabla \cdot \mathbf{B} = 0$ 등)이 별도 강제 없이도 저절로 지켜진다.

2.1. 검산: 발산 조건이 진짜로 저절로 지켜지는가

“저절로”라는 말이 수사가 아니라는 것을 2차원 TM 모드($E_z, H_x, H_y$) 80×80 격자로 확인했다. 갱신은 표준 Yee 순서, 소스는 격자 중앙의 가우시안 펄스, 쿠랑수 $S = 0.5$, 2000 스텝이다. 셀 중심에서 이산 발산

$$(\nabla \cdot \mathbf{B})_{i+1/2,\,j+1/2} = \left[H_x^{i+1,\,j+1/2} - H_x^{i,\,j+1/2}\right] + \left[H_y^{i+1/2,\,j+1} - H_y^{i+1/2,\,j}\right]$$

을 매 25 스텝 전수 조사한 결과 최대값이 $5.6\times10^{-15}$(장 세기로 정규화)였다. 반올림 말고는 아무것도 남지 않는다. 이유는 대수적으로 간단하다 — 위 조합에 갱신식을 대입하면 $E_z$ 의 혼합 차분이 부호만 반대로 두 번 나와 정확히 상쇄된다. 격자 배치를 반 칸 어긋나게 잡은 대가로 가우스 법칙을 이산 수준에서 항등식으로 받은 셈이고, 이게 유한요소법 쪽에서 벡터 유한요소(변 요소)를 쓰는 이유와 같은 동기다.

주의할 점은 이 항등식이 보장하는 것이 “초기값이 발산 없으면 계속 없다”까지라는 것이다. 전류원을 넣을 때 전하 보존을 어기면 그 오차는 그대로 남는다.

3. Leapfrog 시간 전진

시간 방향으로도 $\mathbf{E}$와 $\mathbf{H}$는 반 스텝 엇갈린 시각에 갱신된다. 이를 도약 개구리(leapfrog) 방식이라 부른다.

$$\mathbf{E}^{\,n+1} = \mathbf{E}^{\,n} + \frac{\Delta t}{\varepsilon} \, \nabla \times \mathbf{H}^{\,n+1/2}$$ $$\mathbf{H}^{\,n+3/2} = \mathbf{H}^{\,n+1/2} - \frac{\Delta t}{\mu} \, \nabla \times \mathbf{E}^{\,n+1}$$

전기장은 정수 시각에서, 자기장은 반정수 시각에서 갱신되며 서로가 서로를 밟고 뛰어넘는다. 각 스텝이 직전에 계산된 값만으로 다음 값을 구하는 명시적(explicit) 갱신이라, 큰 행렬을 풀 필요 없이 곱셈과 덧셈만 반복하면 된다. 구현이 단순하고 병렬화가 극단적으로 쉬운 이유가 여기에 있다.

이름과 구조가 심플렉틱 적분기의 리프프로그와 같은 것은 우연이 아니다. 무손실 매질에서 Yee 갱신은 소산이 없다 — 즉 수치적으로 에너지를 갉아먹지 않는다. 64×64 격자에 가우시안 초기조건만 주고 PEC 벽 안에서 2만 스텝을 돌려 이산 에너지 $\frac{1}{2}\sum(E_z^2 + H_x^2 + H_y^2)$ 를 추적하면, 값이 6.3% 폭으로 진동하지만(E와 H가 반 스텝 엇갈린 시각이라 이 순진한 합은 정확한 이산 불변량이 아니다) 2만 스텝 뒤 순변화는 $-8.8\times10^{-5}$ 에 불과하다. 표류가 없다. 손실이 있는 매질을 넣거나 PML을 켜면 당연히 그만큼 줄어든다.

4. 안정 조건 (쿠랑 한계)

공짜 점심은 없다. 명시적 시간 전진의 대가로 FDTD는 시간 간격 $\Delta t$를 함부로 키울 수 없다. 파동이 한 스텝 동안 격자 한 칸 이상을 건너뛰면 정보가 새어나가 해가 폭발하기 때문이다. 3차원에서의 안정 조건, 즉 쿠랑 한계(Courant limit)는 다음과 같다.

$$\Delta t \le \frac{1}{c \sqrt{\dfrac{1}{\Delta x^2} + \dfrac{1}{\Delta y^2} + \dfrac{1}{\Delta z^2}}}$$

이것은 본질적으로 CFL 조건의 전자기 버전이고, 폰 노이만 안정성 해석을 Yee 갱신에 적용하면 그대로 나온다. 균일 격자에서는 쿠랑수 $S = c\Delta t/\Delta$ 로 정리해 $d$ 차원에서 $S \le 1/\sqrt d$ — 즉 1차원 1, 2차원 $0.7071$, 3차원 $0.5774$ 다. 격자를 절반으로 촘촘하게 만들면 시간 간격도 절반으로 줄여야 하므로, 해상도를 2배 올리면 계산량은 (3차원 격자 8배) × (시간 스텝 2배) = 대략 16배로 뛴다. FDTD 엔지니어가 격자 셀 수를 늘리기 전에 심호흡을 하는 이유다.

4.1. 검산: 한계를 0.1%만 넘기면

2차원 TM 코드로 쿠랑수를 한계 주변에서 훑었다(80×80, 가우시안 소스).

$S$	20,000 스텝 뒤 최대 $\lvert E_z\rvert$	판정
0.7071	3.4e-2	안정
0.70711	3.3e-2	안정
0.7072	4.4e-2	안정(한계 위지만 성장 미검출)
0.708	NaN	발산
0.71	NaN (6,000 스텝 내)	발산

이론 한계 $1/\sqrt2 = 0.70710678$ 이 실측과 소수 다섯째 자리까지 맞는다. $S = 0.7072$ 가 20,000 스텝을 버틴 것은 한계를 갓 넘긴 경우 불안정 모드가 $k$ 공간의 격자 대각 코너에만 있고 증폭률이 1에 거의 붙어 있어서, 소스가 그 모드를 충분히 여기하지 못하면 유한 시간 동안 눈에 안 보인다는 뜻이다. “돌려봤는데 안 터졌다”가 안정성의 증거가 아니라는 교과서 경고의 실측판이다. 실무에서 $S$ 를 한계의 0.98~0.99배로 잡는 관행에는 이런 이유도 있다.

5. 수치 분산 — FDTD 의 진짜 오차원

FDTD의 정확도를 실제로 결정하는 것은 안정성보다 수치 분산(수치분산)이다. 이산화된 Yee 격자에서 파동은 진공에서도 정확히 $c$ 로 가지 않는다. 1차원 수치 분산 관계는

$$\frac{1}{c\Delta t}\sin\!\left(\frac{\omega \Delta t}{2}\right) = \frac{1}{\Delta x}\sin\!\left(\frac{k \Delta x}{2}\right)$$

이고, 2차원·3차원에서는 우변이 각 방향 항의 제곱합이 된다. 이 식을 $\omega$ 에 대해 풀어 위상속도를 뽑으면 아래 값들이 나온다. 가로 기준은 셀 크기가 아니라 파장당 셀 수 $N_\lambda$ 다.²

파장당 셀 수 $N_\lambda$	1D, $S=0.5$	1D, $S=0.9$	2D 격자축, $S=1/\sqrt2$	2D 대각선, $S=1/\sqrt2$
10	-1.2412%	-0.3225%	-0.8368%	0.0000%
20	-0.3089%	-0.0787%	-0.2065%	0.0000%
40	-0.0771%	-0.0196%	—	—
80	-0.0193%	-0.0049%	—	—

$N_\lambda$ 를 두 배로 하면 오차가 4배 줄어드는 2차 수렴이 표에 그대로 보인다. 읽어낼 것이 세 가지다.

• 위상속도는 늘 느리다(수치 격자가 매질처럼 작용해 파동을 붙잡는다). 오차가 누적이라는 점이 문제다. $N_\lambda = 20$ 에서 오차는 0.31%뿐이지만 100 파장을 전파하면 위상이 111° 어긋난다. $N_\lambda = 10$ 이면 100 파장에서 1.24 파장, 즉 447° 다. 큰 안테나 배열이나 긴 도파관에서 “격자는 충분히 잘게 했는데 위상이 안 맞는다”는 증상의 출처가 이것이다.
• 쿠랑수를 한계 쪽으로 올리면 분산이 줄어든다. 1D 에서 $S = 1$(매직 타임스텝)이면 분산이 정확히 0 이다($N_\lambda = 4$ 에서도 $v_p/c = 1.00000000000000$). 안정성 여유와 정확도가 서로 반대 방향이라는 뜻이라, 실무의 $S \approx 0.99/\sqrt d$ 는 두 요구의 타협점이다.
• 2차원 이상에서는 비등방성이 생긴다. $S = 1/\sqrt2$ 에서 격자 대각선 방향은 분산이 정확히 0인데 격자축 방향은 $N_\lambda = 10$ 에서 -0.84%다. 같은 파동이 방향에 따라 다른 속도로 간다 — 비등방성 0.84%p. 원형 대칭이어야 할 방사 패턴에 사각형 흔적이 남는 것이 이 때문이고, 격자를 촘촘히 하는 것 외에 고차 차분이나 (2,4) 스킴을 쓰는 우회로가 있다.

6. 흡수 경계조건 (PML)

무한히 열린 공간에서 벌어지는 방사 문제를 유한한 격자에 담으려면, 격자 끝에 도달한 파동을 반사 없이 삼켜버리는 장치가 필요하다. 이 역할을 하는 것이 완전정합층(Perfectly Matched Layer, PML)이다.³ 1994년 베렝제(Berenger)가 제안한 이 인공 흡수층은 격자 가장자리를 둘러싸 마치 공간이 무한히 이어지는 것처럼 파동을 감쇠시킨다. PML이 등장하기 전의 흡수 경계조건들(Mur, Engquist-Majda 등 단방향 파동식 기반)은 사입사에서 반사가 커서 FDTD의 신뢰도를 갉아먹었는데, PML이 이 문제를 사실상 해결하면서 FDTD의 실용성이 급상승했다.

핵심 아이디어는 층 안에서 전기·자기 손실을 파동 임피던스가 유지되도록 짝지어 넣는 것이다($\sigma_e/\varepsilon = \sigma_m/\mu$). 그러면 연속 방정식 수준에서 반사가 0이 된다. 이산 격자에서는 전도도가 셀 경계에서 갑자기 뛰면 그 자체가 반사를 만들므로, $\sigma$ 를 깊이 방향으로 다항식($m$ 차, 보통 3~4)으로 완만하게 키운다.

6.1. 검산: 층 두께가 반사를 얼마나 줄이는가

1차원 FDTD($S = 0.5$, 800 셀)의 오른쪽 끝에 다항 그레이딩 정합 손실층을 넣고, 프로브를 지나간 입사 펄스와 되돌아온 펄스의 피크비를 dB로 쟀다. $\sigma_{\max}$ 는 표준 공식 $-(m+1)\ln R_0 / (2 \eta d)$ 로 잡았다.

층 두께	차수 $m$	설계 $R_0$	$\sigma_{\max}$	실측 반사
없음(PEC 종단)	—	—	—	-1.5 dB
5 셀	3	1e-6	5.526	-71.6 dB
10 셀	2	1e-6	2.072	-86.2 dB
10 셀	3	1e-6	2.763	-94.4 dB
20 셀	3	1e-8	1.842	-114.6 dB

종단을 그냥 금속으로 두면 거의 전부 되돌아오는데, 10 셀만 깔면 반사가 $-94$ dB, 진폭비로 $2\times10^{-5}$ 로 떨어진다. 실무에서 PML 두께를 8~12 셀로 잡는 관행의 근거가 이 자리다. 두께를 두 배로 하면 20 dB가 더 떨어지는 대신 격자 셀이 늘고, $R_0$ 설계값을 무작정 낮춰도(10셀 $R_0 = 10^{-8}$ 행) $\sigma$ 기울기가 급해져 오히려 나빠질 수 있다는 것도 표에 보인다. 이 숫자들은 1차원 수직입사라 정합 손실층과 PML이 같은 것이지만, 2·3차원 사입사에서는 베렝제의 분할장(split-field) 또는 축이방성 UPML(Gedney 1996), 복소 주파수 이동 CFS-PML 같은 정식 구성이 필요하다.

7. 매질과 형상 — 실무에서 발목 잡는 곳

• 손실·분산 매질. 상수 $\varepsilon, \mu, \sigma$ 는 갱신식에 계수 두 개로 흡수된다. 주파수에 따라 $\varepsilon(\omega)$ 가 변하는 매질(금속의 드루드 모델, 물의 디바이 완화, 로런츠 공명)은 시간영역에서 곱셈이 합성곱이 되므로 보조 미분방정식(ADE)이나 재귀 합성곱(RC)으로 풀어낸다. 금속 나노구조·실리콘 포토닉스 계산이 다 이 계열이다.
• 계단 근사. 정렬 격자에 곡면을 넣으면 벽이 톱니처럼 근사된다. 이게 공진 주파수를 수백 MHz 단위로 밀어버리는 일이 흔해서, 부분 셀 면적·길이로 갱신 계수를 보정하는 정합 FDTD(conformal FDTD, Dey-Mittra 1997)나 얇은 판·세선 전용 서브셀 모델이 표준 장비가 됐다.
• 주기 구조. 광결정이나 메타물질 단위셀은 주기경계조건으로 처리하고, 사입사에서는 시간영역 특유의 위상 지연 문제 때문에 사인 정상상태 해석이나 분할장 기법을 쓴다.

8. 강점과 약점

구분	내용
강점	시뮬레이션 한 번으로 광대역 응답 확보 (펄스 → 푸리에 변환)
강점	시간 응답과 과도 현상을 직접 관찰 (시각화에 최적)
강점	알고리즘이 단순하고 병렬화·GPU 가속이 쉬움 (인접 셀만 참조)
강점	비선형·시변 매질을 넣기 쉽다 (주파수영역 기법의 약점)
약점	정렬 격자라 곡면이 계단 근사(staircasing)됨
약점	수치 분산으로 위상이 전파 거리에 비례해 어긋남
약점	큰 빈 공간까지 전부 격자로 채워야 함 (자유공간 산란은 모멘트법이 유리)
약점	고Q 공진기는 감쇠 시간만큼 오래 돌려야 해서 스텝 수가 폭발

9. 어디에 쓰나

안테나 해석의 광대역 임피던스·방사 패턴, 레이더 단면적 계산, 도파관과 전송선 이론 쪽 불연속 구조의 S-파라미터 추출, 광결정·실리콘 포토닉스·메타물질 설계, 인체 모델 전자파 흡수율(SAR) 평가, EMC/EMI 차폐 해석이 대표적이다. 오픈소스로는 Meep 과 openEMS 가 널리 쓰이고, 상용 쪽은 시간영역 솔버를 기본으로 깐 패키지가 여럿이다. 셀 수가 억 단위로 가는 문제에서는 GPU 한 장이 CPU 수십 코어를 이기는 드문 분야 중 하나다.⁴

10. 관련 문서

• 맥스웰 방정식 · 전자기해석 · 전자기학
• 유한차분법 · 유한요소법 · 모멘트법
• CFL 조건 · 폰 노이만 안정성 해석 · 수치분산
• 완전정합층 · 주기경계조건 · 심플렉틱 적분기
• 안테나 해석 · 레이더 단면적 · 도파관 · 전송선 이론
• 광결정 · 실리콘 포토닉스 · 메타물질 · 드루드 모델

11. Footnotes

1. Yee, K. S. (1966). “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media”. IEEE Trans. Antennas Propag. 14, 302–307. 당시엔 크게 주목받지 못하다가 1970~80년대 컴퓨터 성능이 따라오면서 뒤늦게 재발견되어 전자기해석의 주력 기법으로 등극했다. 논문 한 편으로 격자에 이름을 남긴 케이스. ↩

2. 그래서 FDTD 논문의 격자 해상도 표기는 셀 크기가 아니라 파장당 셀 수여야 의미가 있다. “0.1 mm 격자”는 정보가 아니고 "$\lambda/20$"이 정보다. 물질 안에서는 파장이 $\sqrt{\varepsilon_r}$ 배 짧아지므로 유전율 높은 영역의 $N_\lambda$ 가 실제 병목이 된다는 점도 자주 잊힌다. ↩

3. 이름부터가 “완벽하게(Perfectly)“라니 자신감이 대단하지만, 실제로는 이산 격자에서 완벽하진 않고 아주 약간의 반사가 남는다. 위 검산의 $-94$ dB 가 그 “아주 약간”의 실측값이다. 그래도 이전 기법 대비 반사가 수십 dB 낮아서 실무에서는 충분히 “완벽하다”고 쳐준다. ↩

4. 갱신식이 인접 셀만 보고 분기가 없어서 메모리 대역폭이 곧 성능이기 때문이다. 반대로 말하면 아무리 최적화해도 대역폭 벽에서 멈춘다는 뜻이고, 그래서 FDTD 가속 논문의 절반이 메모리 접근 패턴 이야기다. ↩