OghmaNano의 FDTD

1. 소개

유한차분 시간영역(FDTD)은 이산 공간 격자에서 시간영역에서 직접 맥스웰 방정식을 풀기 위한 전파(full-wave) 수치 방법입니다. 광선 추적, 모드 해석기 또는 transfer-matrix 접근법과 달리, FDTD는 기하학적 또는 모드 가정을 하지 않습니다. 대신 전기장과 자기장은 curl 방정식—암페르 법칙과 패러데이 법칙—을 이산화하여 시간에 따라 명시적으로 전진되므로 전자기파가 시뮬레이션 영역을 통해 자연스럽게 전파, 산란, 간섭 및 회절합니다.

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

유한차분 시간영역(FDTD) 방법에서는 맥스웰 방정식을 시간영역에서 직접 풉니다. 이를 위해 공간 미분을 엇갈린 격자(보통 Yee 격자)에서의 유한차분으로 대체하고 필드를 시간에 따라 전진시킵니다. 이 접근법은 전파 방향, 모드 구조 또는 필드 거동에 대해 본질적으로 아무 가정도 하지 않기 때문에, 단일 FDTD 시뮬레이션으로 광대역 응답, 회절 및 간섭, 근접장 효과, 서브파장 기하구조, 과도 전자기 현상을 하나의 프레임워크 안에서 포착할 수 있습니다. 이러한 단순화 가정의 부재가 FDTD를 매우 일반적으로 만들지만, 동시에 계산 비용이 많이 드는 이유이기도 합니다. 공간 격자는 가장 작은 물리적 길이 스케일을 분해해야 하고, 시간 스텝은 수치 안정성에 의해 엄격하게 제한되며, 3차원 시뮬레이션은 겉보기에 크기가 크지 않은 기하구조에 대해서도 빠르게 메모리와 시간을 많이 소모하게 됩니다.

2. FDTD를 사용해야 할 때와 사용하지 말아야 할 때

대부분의 실제 광학 문제에서 FDTD는 가장 먼저 선택해야 할 방법이 아닙니다. 대규모 광학 시스템, 긴 전파 거리 또는 기하광학이 지배적인 상황에서 FDTD는 극도로 비효율적인 선택입니다. 영상 시스템, 거시적 렌즈, 개구, 그리고 밀리미터에서 센티미터 크기의 광학 어셈블리는 광선 추적으로 다루는 편이 훨씬 낫고, 박막 태양전지와 같은 평면 다층 구조, OLED 적층 구조 및 코팅은 transfer-matrix 방법으로 더 정확하고 효율적으로 모델링됩니다. 이러한 영역에서 FDTD는 계산용 큰 망치와 같습니다. 물리적 통찰을 거의 또는 전혀 더 제공하지 않으면서 실행 시간과 메모리 사용량을 극적으로 증가시킵니다.

실제로 FDTD는 최후의 수단 방법으로 간주되어야 하며, 보다 단순한 모델이 진정으로 한계에 부딪히는 경우에만 사용해야 합니다. 그 강점은 서브파장 격자로부터의 회절, 포토닉 크리스털, 메타표면, 강하게 무질서하거나 거친 계면, 그리고 광선광학 또는 층상 매질 모델로는 포착할 수 없는 근접장 결합 효과와 같이 전체 전자기장을 명시적으로 분해해야 하는 상황에 있습니다. 이러한 경우에도 FDTD는 수치 비용, 안정성 제약, 메모리 요구사항을 명확히 이해한 상태에서 제한적이고 의도적으로 적용되어야 합니다. 따라서 다음 절에서는 OghmaNano에서 사용되는 이산 업데이트 방정식을 암페르 법칙과 패러데이 법칙으로부터 직접 유도합니다. 이는 무분별한 사용을 장려하기 위한 것이 아니라, FDTD가 무엇을 할 수 있는지와 왜 신중히 사용되어야 하는지를 명확히 하기 위한 것입니다.

3. 이론적 배경

매뉴얼의 이 섹션은 이해를 돕고 오류를 찾아낼 수 있도록 자세한 유도를 통해 FDTD 코드를 완전히 설명하는 것을 목표로 합니다.

암페르 법칙은 \[\sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t} = \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ H_{x} & H_{y} & H_{z} \end{vmatrix}\]로 주어집니다

이는 다음과 같이 전개될 수 있습니다

\[\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} = \frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}\]

\[\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\]

\[\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}\]

\(\frac{\partial}{\partial y}=0\)인 경우

\[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x} \end{split}\]

\(E_{x}\)에 대해 \[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma \frac{E_{x}^{t+1}[]+E_{x}^{t}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]-E_{x}^{t}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & \frac{\sigma \Delta t + 2 \epsilon }{ 2 \Delta t}E_{x}^{t+1}[] = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & E_{x}^{t+1}[] = \left ( -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

\(E_{y}\)에 대해 \[\begin{split} &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma \frac{E_{y}^{t+1}[]+E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]-E_{y}^{t}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{y}^{t+1}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t}\\ &E_{y}^{t+1}[] = \left ( -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

\(E_{z}\)에 대해 \[\begin{split} &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}\\ &\sigma \frac{E_{z}^{t+1}[]+E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]-E_{z}^{t}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}\\ &\sigma \frac{E_{z}^{t+1}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t}\\ &E_{z}^{t+1}[]= \left ( \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon}\\ \end{split}\]

패러데이 법칙은 \[-\sigma_{m} \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t} = \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ E_{x} & E_{y} & E_{z} \end{vmatrix}\]로 주어집니다

이는 전개하면 다음과 같습니다:

\[-\sigma_{m} H_{x} - \mu \frac{\partial H_{x}}{\partial t} = \frac{\partial E_{z}}{\partial y}-\frac{\partial E_{y}}{\partial z}\]

\[-\sigma_{m} H_{y} - \mu \frac{\partial H_{y}}{\partial t} = -\frac{\partial E_{z}}{\partial x}+\frac{\partial E_{x}}{\partial z}\]

\[-\sigma_{m} H_{z} - \mu \frac{\partial H_{z}}{\partial t} = \frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}\]

\(\sigma_m=0\) 그리고 \(\frac{\partial}{\partial y}=0\)일 때

\[\begin{split} &\frac{\partial H_{x}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial z} \right )\\ &\frac{\partial H_{y}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{z}}{\partial x}-\frac{\partial E_{x}}{\partial z} \right )\\ &\frac{\partial H_{z}}{\partial t} = - \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial x} \right ) \end{split}\]

이를 이산화하면 다음이 됩니다

\[\begin{split} & H_{x}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t + H_{x}^{t}[]\\ & H_{y}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\frac{E_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t+ H_{y}^{t}[]\\ & H_{z}^{t+1} = \frac{1}{\mu} \left ( - \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x} \right ) \Delta t + H_{x}^{z}[] \end{split}\]

4. 컴퓨터에서의 수치 해

위에서 유도한 방정식은 전기장과 자기장의 연속 시간, 연속 공간에서의 진화를 기술합니다. 이를 컴퓨터에서 풀기 위해서는 공간과 시간 모두를 이산화해야 합니다. 시뮬레이션 영역은 규칙적인 셀 격자로 나뉘고, 필드 성분은 이산 공간 위치와 이산 시간 스텝에 저장됩니다. 표준 FDTD 정식화에서는 전기장과 자기장이 공간과 시간 모두에서 엇갈리게 배치되며, 따라서 각 필드 성분은 curl 방정식에 의해 자연스럽게 요구되는 위치에서 평가됩니다.

실제로 이는 전기장 성분이 정수 시간 스텝에 저장됨을 의미합니다 (\(t, t+\Delta t, t+2\Delta t\)), 반면 자기장 성분은 반정수 시간 스텝에 저장됩니다 (\(t+\tfrac{1}{2}\Delta t, t+\tfrac{3}{2}\Delta t\)). 다음과 같은 공간 미분 \(\partial / \partial x\) 및 \(\partial / \partial z\)는 인접한 격자점 사이의 중심 유한차분으로 대체됩니다. 이러한 엇갈림은 필드가 시간에 따라 서로를 leapfrog하는 단순한 명시적 업데이트 체계를 만듭니다.

각 시간 스텝에서 자기장은 먼저 패러데이 법칙의 이산화된 형태를 사용하여 업데이트되며, 이때 전기장의 curl은 현재 전기장 값으로부터 평가됩니다. 자기장이 반 시간 스텝만큼 전진한 후, 전기장은 암페르 법칙의 이산화된 형태를 사용하여 업데이트되며, 여기에는 유전율, 투자율, 전도도와 같은 물질 매개변수가 포함됩니다. 전역 행렬 해법은 필요하지 않습니다. 각 격자 셀의 각 필드 성분은 이전 시간 스텝의 국소 이웃만 사용하여 업데이트됩니다.

이 명시적 시간 스테핑 절차는 원하는 시뮬레이션 시간에 도달할 때까지 반복됩니다. 이 방법은 시간영역 방법이므로, 단일 시뮬레이션은 자연스럽게 넓은 주파수 범위에 대한 정보를 포함합니다. 스펙트럼이나 정상상태 필드 분포와 같은 주파수영역 물리량은 일반적으로 필드 값을 시간의 함수로 기록한 뒤 푸리에 변환을 사용하여 후처리함으로써 얻습니다. 광원, 경계, 흡수층 (예: perfectly matched layer)은 적절한 격자 위치에서 업데이트 방정식에 직접 포함됩니다.

업데이트 방정식의 단순성은 FDTD를 개념적으로 이해하기 쉽고 병렬화하기 매우 용이하게 만들지만, 동시에 엄격한 수치 제약도 부과합니다. 시간 스텝은 공간 격자 간격과 물질 특성에 따라 결정되는 안정성 조건을 만족해야 하며, 메모리 사용량은 격자 셀 수와 저장된 필드 성분 수에 직접 비례하여 증가합니다. 이러한 제약은 궁극적으로 실제로 수행할 수 있는 시뮬레이션의 크기와 지속 시간을 제한하며, 격자 해상도, 영역 크기, 모델 복잡성을 신중하게 선택해야 할 필요성을 만들어냅니다.