FDTD no OghmaNano

1. Introdução

Finite-Difference Time-Domain (FDTD) é um método numérico de onda completa para resolver as equações de Maxwell diretamente no domínio do tempo em uma grade espacial discreta. Em contraste com ray tracing, solucionadores de modos ou abordagens de matriz de transferência, o FDTD não faz hipóteses geométricas ou modais. Em vez disso, os campos elétrico e magnético são avançados explicitamente no tempo por meio da discretização das equações de rotacional — lei de Ampère e lei de Faraday — de modo que as ondas eletromagnéticas se propaguem, espalhem, interfiram e sofram difração naturalmente através do domínio de simulação.

\[ \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}. \]

No método finite-difference time-domain (FDTD), as equações de Maxwell são resolvidas diretamente no domínio do tempo por meio da substituição das derivadas espaciais por diferenças finitas em uma grade escalonada (tipicamente uma malha de Yee) e do avanço dos campos no tempo. Como essa abordagem não faz essencialmente nenhuma hipótese sobre direção de propagação, estrutura modal ou comportamento do campo, uma única simulação FDTD pode capturar resposta de banda larga, difração e interferência, efeitos de campo próximo, geometria subcomprimento de onda e fenômenos eletromagnéticos transientes dentro de um único framework. Essa ausência de hipóteses simplificadoras é o que torna o FDTD tão geral, mas também é o que o torna computacionalmente caro: a grade espacial deve resolver a menor escala de comprimento físico, o passo de tempo é fortemente restringido pela estabilidade numérica, e simulações tridimensionais rapidamente se tornam intensivas em memória e tempo, mesmo para geometrias que parecem modestas em tamanho.

2. Quando usar FDTD e quando não usar FDTD

Na maioria dos problemas ópticos práticos, o FDTD não deve ser o primeiro método escolhido. Para sistemas ópticos em grande escala, longas distâncias de propagação ou situações dominadas por óptica geométrica, o FDTD é uma escolha extremamente ineficiente. Sistemas de imageamento, lentes macroscópicas, aberturas e conjuntos ópticos com dimensões de milímetros a centímetros são muito melhor tratados usando ray tracing, enquanto estruturas planares multicamadas, como células solares de filme fino, pilhas OLED e revestimentos, são modeladas com mais precisão e eficiência usando métodos de matriz de transferência. Nesses regimes, o FDTD atua como uma marreta computacional: ele aumenta drasticamente o tempo de execução e o uso de memória enquanto fornece pouca ou nenhuma informação física adicional.

Na prática, o FDTD deve ser considerado um método de último recurso, reservado para casos em que modelos mais simples realmente falham. Sua força está em situações nas quais o campo eletromagnético completo deve ser resolvido explicitamente, como difração em grades subcomprimento de onda, cristais fotônicos, metassuperfícies, interfaces fortemente desordenadas ou rugosas e efeitos de acoplamento de campo próximo que não podem ser capturados por óptica de raios ou modelos de meios estratificados. Mesmo nesses casos, o FDTD deve ser aplicado de forma restrita e deliberada, com uma compreensão clara de seu custo numérico, restrições de estabilidade e requisitos de memória. As seções a seguir portanto derivam as equações discretas de atualização usadas no OghmaNano diretamente das leis de Ampère e Faraday, não para incentivar o uso indiscriminado, mas para explicitar tanto o que o FDTD pode fazer quanto por que ele deve ser usado com parcimônia.

3. Fundamentos teóricos

Esta seção do manual tem como objetivo descrever o código FDTD por completo, com derivações verbosas para ajudar na compreensão/detecção de erros.

A lei de Ampère é dada por \[\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}\]

que pode ser expandida como

\[\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}\]

Para o caso \(\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}\]

para \(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}\]

para \(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}\]

para \(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}\]

A lei de Faraday é dada por \[-\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}\]

que pode ser expandida para dar:

\[-\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}\]

Com \(\sigma_m=0\) e \(\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}\]

que, ao discretizar, dá

\[\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. Solução numérica em um computador

As equações derivadas acima descrevem a evolução contínua no tempo e no espaço dos campos elétrico e magnético. Para resolvê-las em um computador, tanto o espaço quanto o tempo são discretizados. O domínio de simulação é dividido em uma grade regular de células, e os componentes de campo são armazenados em localizações espaciais discretas e passos de tempo discretos. Na formulação FDTD padrão, os campos elétrico e magnético são escalonados tanto no espaço quanto no tempo, de modo que cada componente de campo seja avaliado onde é naturalmente requerido pelas equações de rotacional.

Na prática, isso significa que os componentes do campo elétrico são armazenados em passos de tempo inteiros (\(t, t+\Delta t, t+2\Delta t\)), enquanto os componentes do campo magnético são armazenados em passos de tempo seminteiros (\(t+\tfrac{1}{2}\Delta t, t+\tfrac{3}{2}\Delta t\)). Derivadas espaciais como \(\partial / \partial x\) e \(\partial / \partial z\) são substituídas por diferenças finitas centradas entre pontos vizinhos da grade. Esse escalonamento leva a um esquema de atualização explícito e simples no qual os campos saltam alternadamente um sobre o outro no tempo.

Em cada passo de tempo, o campo magnético é atualizado primeiro usando a forma discretizada da lei de Faraday, com o rotacional do campo elétrico avaliado a partir dos valores atuais do campo elétrico. Uma vez que o campo magnético tenha sido avançado em meio passo de tempo, o campo elétrico é então atualizado usando a forma discretizada da lei de Ampère, incluindo parâmetros de material como permissividade, permeabilidade e condutividade. Nenhuma solução global de matriz é necessária: cada componente de campo em cada célula da grade é atualizado usando apenas seus vizinhos locais do passo de tempo anterior.

Esse procedimento explícito de avanço temporal é repetido até que o tempo de simulação desejado seja alcançado. Como o método é no domínio do tempo, uma única simulação contém naturalmente informação sobre uma ampla faixa de frequências; quantidades no domínio da frequência, como espectros ou distribuições de campo em regime estacionário, são tipicamente obtidas registrando valores de campo em função do tempo e pós-processando-os com transformadas de Fourier. Fontes, contornos e camadas absorventes (como camadas perfeitamente casadas) são incorporados diretamente nas equações de atualização nas localizações apropriadas da grade.

A simplicidade das equações de atualização torna o FDTD conceitualmente simples e altamente paralelizável, mas isso também impõe restrições numéricas rigorosas. O passo de tempo deve satisfazer uma condição de estabilidade que depende do espaçamento da grade espacial e das propriedades do material, e o consumo de memória escala diretamente com o número de células da grade e componentes de campo armazenados. Essas restrições acabam limitando o tamanho e a duração das simulações que podem ser realizadas na prática e motivam uma escolha cuidadosa da resolução da grade, tamanho do domínio e complexidade do modelo.