Derivação do método FDTD.

1. Introdução

Esta página apresenta a derivação matemática detalhada das principais equações de atualização do método de diferenças finitas no domínio do tempo (FDTD), utilizadas para resolver as equações de Maxwell no domínio do tempo. O objetivo desta página é diferente do da página principal sobre a teoria. O objetivo aqui não é fornecer uma visão geral simplificada, mas sim mostrar explicitamente como as equações contínuas de rotação são convertidas em expressões algébricas discretas, que são processadas por um computador. Esta é também a parte do método FDTD que é mais fácil de implementar incorretamente no código: erros de sinal, índices mal colocados ou intervalos de tempo inconsistentes produzirão imediatamente resultados instáveis ou fisicamente irrealistas.

Partindo, portanto, da lei de Ampère e da lei de Faraday na sua forma diferencial, decompomos estas leis em componentes cartesianas e, em seguida, aplicamos-as a um problema ótico bidimensional. As derivadas espaciais são substituídas por diferenças finitas centradas, as derivadas temporais são discretizadas utilizando o esquema padrão de «leapfrog» e as expressões resultantes são reorganizadas para fornecer fórmulas explícitas de atualização para os componentes do campo elétrico e magnético.

Para facilitar a compreensão da dedução, reduzimos as equações a um caso bidimensional, no qual \(\partial/\partial y = 0\) e a condutividade magnética são definidas como zero, \(\sigma_m = 0\). Isso torna a álgebra muito mais fácil de acompanhar e demonstra explicitamente como a redução dimensional é realizada na prática. No solucionador OghmaNano FDTD completo, no entanto, as equações eletromagnéticas são resolvidas em três dimensões. Portanto, a dedução bidimensional apresentada aqui tem como objetivo ser um exemplo claro e prático do método, e não uma afirmação de que o solucionador subjacente está limitado a 2D.

Ao longo desta página, os componentes do campo elétrico são avaliados em intervalos de tempo inteiros \(n\Delta t\), enquanto os componentes do campo magnético são avaliados em intervalos de tempo semi-inteiros \(\left(n+\tfrac{1}{2}\right)\Delta t\). Este é o esquema padrão de malha de Yee, utilizado no método FDTD. Ao contrário da página de resumo mais concisa, a indexação é apresentada de forma mais explícita aqui, para que a derivação possa ser lida como um guia de implementação, bem como uma nota matemática.

2. Equações de partida

Partimos das equações de Maxwell, expressas na forma de rotação, num meio linear com condutividade elétrica \(\sigma\), permissividade elétrica \(\epsilon\), permeabilidade magnética \(\mu\) e, para generalizar, condutividade magnética \(\sigma_m\):

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

Na forma determinante, os cachos são

\[ \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix}, \qquad \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix}. \]

Estas são as equações em tempo contínuo e espaço contínuo. O objetivo do restante desta página é convertê-las num esquema explícito de diferenças finitas, que pode ser aplicado iterativamente ao longo do tempo.

3. Forma vetorial das equações do rotacional.

Ao expandir a lei de Ampère em componentes cartesianos, obtemos:

\[ \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_x}{\partial z} - \frac{\partial H_z}{\partial x}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}. \]

A expansão da lei de Faraday resulta em:

\[ -\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_x}{\partial z} - \frac{\partial E_z}{\partial x}, \]

\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}. \]

Nesta fase, ainda não se fez qualquer suposição sobre a dimensionalidade. As expressões apresentadas acima são simplesmente a forma componente das equações de rotação tridimensional completas.

4. Redução bidimensional

Para tornar a derivação mais fácil de resolver, reduzimos agora o problema a duas dimensões, partindo do princípio de que não existe variação na direção \(y\), e, portanto, \(\partial / \partial y = 0\). Esta etapa é realizada unicamente para simplificar a álgebra e facilitar a compreensão da estrutura das equações de atualização. O mesmo procedimento pode ser aplicado diretamente a três dimensões e, na prática, o OghmaNano resolve a forma completa tridimensional das equações de Maxwell.

Partindo deste pressuposto, as equações do campo elétrico simplificam-se para:

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

Se também definirmos a condutividade magnética como zero, \(\sigma_m = 0\), as equações do campo magnético tornar-se-ão:

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}, \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}, \]

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Estas são as equações contínuas que serão agora discretizadas. É importante referir explicitamente que se trata de uma redução bidimensional do problema tridimensional completo. O objetivo não é afirmar que todo o método FDTD é intrinsecamente bidimensional, mas sim obter equações de atualização claras para o caso comum em que uma direção é suprimida.

5. Grade de Yee, indexação e notação.

Para resolver as equações acima numericamente, tanto o espaço quanto o tempo são discretizados. Sejam \(\Delta x\) e \(\Delta z\) os espaçamentos da malha computacional nas direções \(x\) e \(z\), respetivamente, e \(\Delta t\) o intervalo de tempo. Utilizamos os índices inteiros \(i\) e \(k\) para as direções \(x\) e \(z\):

\[ x_i = i \Delta x, \qquad z_k = k \Delta z. \]

O tempo é indexado pelo número inteiro \(n\), de modo que

\[ t^n = n \Delta t. \]

No esquema de Yee, os campos elétrico e magnético não são armazenados nos mesmos locais, seja no espaço ou no tempo. Em vez disso, eles são dispostos de forma alternada. Os componentes do campo elétrico são calculados em intervalos de tempo inteiros \(n\Delta t\), enquanto os componentes do campo magnético são calculados em intervalos de tempo semi-inteiros \(\left(n+\tfrac{1}{2}\right)\Delta t\). Esta é a origem da estrutura em «salto» do método FDTD.

No plano espacial, os componentes do campo também são deslocados de forma que cada operador rotacional de diferenças finitas, centrado, seja naturalmente definido entre as localizações de grade adjacentes. Numa implementação completa do código, cada componente do campo possui o seu próprio par de índices espaciais, por exemplo \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\) e assim por diante. Para manter a clareza das expressões matemáticas, alguns dos índices repetidos, que funcionam como «observadores», nem sempre são escritos explicitamente nos passos que se seguem. No entanto, a interpretação é sempre a mesma:

\(i\) indexa a direção \(x\).
\(k\) indexa a direção \(z\).
\(n\) indica os níveis de tempo do campo elétrico, expressos em valores inteiros, e
\(n+\tfrac{1}{2}\) indica os níveis temporais do campo magnético em incrementos de meio passo.

Por exemplo, a aproximação de diferenças finitas centrada para \(\partial H_y / \partial z\) é:

\[ \frac{\partial H_y}{\partial z} \approx \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]

e, da mesma forma, a aproximação centralizada de \(\partial H_z / \partial x\) é:

\[ \frac{\partial H_z}{\partial x} \approx \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Com esta notação explicitada, a derivação apresentada abaixo pode ser interpretada de forma inequívoca como um esquema adequado de atualização da malha de Yee.

6. Obtenção das equações para a atualização do campo elétrico.

O campo elétrico é calculado em intervalos de tempo inteiros \(n\Delta t\), enquanto o campo magnético é calculado em intervalos de tempo semi-inteiros \(\left(n+\tfrac{1}{2}\right)\Delta t\). Para obter as equações de atualização, a derivada temporal é aproximada por uma diferença finita centralizada e o termo condutor é calculado utilizando a média dos campos elétricos antigo e novo. Isso resulta num tratamento numérico estável, do tipo trapezoidal, da condutividade.

Esta etapa de cálculo da média é importante. Se o termo de condutividade fosse tratado de forma inconsistente, os coeficientes de atualização deixariam de ter a forma padrão amortecida e a correspondência com as equações FDTD habituais seria perdida.

6.1 Atualizar a equação para \(E_x\).

Comece por

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}. \]

Calcule o termo de condutividade utilizando a média de \(E_x^{n+1}\) e \(E_x^n\) e calcule a derivada temporal utilizando uma diferença finita central:

\[ \sigma \frac{E_x^{n+1}(i,k) + E_x^n(i,k)}{2} + \epsilon \frac{E_x^{n+1}(i,k) - E_x^n(i,k)}{\Delta t} = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Reúna os termos \(E_x^{n+1}\) no lado esquerdo:

\[ \left( \frac{\sigma}{2} + \frac{\epsilon}{\Delta t} \right) E_x^{n+1}(i,k) = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} + \left( \frac{\epsilon}{\Delta t} - \frac{\sigma}{2} \right) E_x^n(i,k). \]

Portanto.

\[ E_x^{n+1}(i,k) = \frac{ \left( \dfrac{\epsilon}{\Delta t} - \dfrac{\sigma}{2} \right) E_x^n(i,k) - \dfrac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} }{ \dfrac{\epsilon}{\Delta t} + \dfrac{\sigma}{2} }. \]

Multiplicar o numerador e o denominador por \(2\Delta t\) resulta na forma de coeficiente que já conhecemos:

\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

6.2 Atualizar a equação para \(E_y\).

Comece por

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}. \]

Discretize ambas as expressões de derivada na grelha de Yee:

\[ \sigma \frac{E_y^{n+1}(i,k) + E_y^n(i,k)}{2} + \epsilon \frac{E_y^{n+1}(i,k) - E_y^n(i,k)}{\Delta t} = - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Reorganizar proporciona.

\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right]. \]

Este é o único componente do campo elétrico na redução bidimensional que contém contribuições de duas derivadas distintas do campo magnético e, portanto, é o mais fácil de interpretar incorretamente no código. Em particular, o sinal relativo entre os termos de derivada \(x\) e \(z\) decorre diretamente da expansão do rotacional e não deve ser alterado.

6.3 Atualizar a equação para \(E_z\).

Comece por

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

A discretização permite

\[ \sigma \frac{E_z^{n+1}(i,k) + E_z^n(i,k)}{2} + \epsilon \frac{E_z^{n+1}(i,k) - E_z^n(i,k)}{\Delta t} = \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Portanto/Assim sendo/Desse modo/Em consequência.

\[ E_z^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_z^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Neste ponto, as equações discretas para a atualização do campo elétrico já foram obtidas na sua forma padrão com amortecimento. O passo seguinte é derivar as equações correspondentes para o campo magnético, calculadas em intervalos de tempo semi-inteiros.

7. Obtenção das equações para a atualização do campo magnético.

O campo magnético é atualizado em intervalos de tempo que correspondem a múltiplos semi-inteiros. Com \(\sigma_m = 0\), as equações do campo magnético não contêm um termo de amortecimento, pelo que as fórmulas de atualização assumem uma forma mais simples, de progressão direta. Esta parte da derivação é algebricamente mais curta do que o caso do campo elétrico, porque não existe um termo de condutividade trapezoidal que precise de ser considerado.

7.1 Atualizar a equação para \(H_x\).

Comece por

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}. \]

Ao discretizar no domínio do tempo entre os pontos médios \(n-\tfrac{1}{2}\) e \(n+\tfrac{1}{2}\) e no espaço, utilizando uma diferença central na direção \(z\), obtém-se:

\[ \mu \frac{ H_x^{n+\frac{1}{2}}(i,k+\frac{1}{2}) - H_x^{n-\frac{1}{2}}(i,k+\frac{1}{2}) }{\Delta t} = \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}. \]

Na notação de índice semi-completo, frequentemente utilizada em demonstrações, o mesmo resultado pode ser escrito da seguinte forma:

\[ \mu \frac{ H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Portanto.

\[ H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

7.2 Atualizar a equação para \(H_y\).

Comece por

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}. \]

A discretização permite

\[ \mu \frac{ H_y^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) - H_y^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) }{\Delta t} = \frac{ E_z^n(i+1,k+\frac{1}{2}) - E_z^n(i,k+\frac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\frac{1}{2},k+1) - E_x^n(i+\frac{1}{2},k) }{\Delta z}. \]

Mais uma vez, na notação simplificada do índice parcial,

\[ \mu \frac{ H_y^{n+\frac{1}{2}} - H_y^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Portanto/Assim sendo/Desse modo/Em consequência.

\[ H_y^{n+\frac{1}{2}} = H_y^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z} \right]. \]

7.3 Atualizar a equação para \(H_z\).

Comece por

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

A discretização permite

\[ \mu \frac{ H_z^{n+\frac{1}{2}}(i+\frac{1}{2},k) - H_z^{n-\frac{1}{2}}(i+\frac{1}{2},k) }{\Delta t} = - \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

De forma equivalente, numa notação mais concisa:

\[ \mu \frac{ H_z^{n+\frac{1}{2}} - H_z^{n-\frac{1}{2}} }{\Delta t} = - \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Portanto.

\[ H_z^{n+\frac{1}{2}} = H_z^{n-\frac{1}{2}} - \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Com estas substituições concluídas, as equações padrão de atualização do método FDTD (Finite-Difference Time-Domain) podem ser aplicadas diretamente.

8. Equações finais de atualização.

Após reunir os resultados, as equações de atualização bidimensionais do método FDTD são:

\[ H_x^{n+\frac{1}{2}}(i,k+\tfrac{1}{2}) = H_x^{n-\frac{1}{2}}(i,k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}, \]

\[ H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) = H_y^{n-\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n(i+1,k+\tfrac{1}{2}) - E_z^n(i,k+\tfrac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\tfrac{1}{2},k+1) - E_x^n(i+\tfrac{1}{2},k) }{\Delta z} \right], \]

\[ H_z^{n+\frac{1}{2}}(i+\tfrac{1}{2},k) = H_z^{n-\frac{1}{2}}(i+\tfrac{1}{2},k) - \frac{\Delta t}{\mu} \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Quando apresentado desta forma, o significado do esquema de «salto» torna-se evidente. O campo elétrico no instante \(n+1\) é obtido a partir do campo magnético no instante \(n+\tfrac{1}{2}\), enquanto o campo magnético no instante \(n+\tfrac{1}{2}\) é obtido a partir do campo elétrico no instante \(n\). Assim, os campos evoluem numa sequência alternada ao longo do tempo.

9. Observações sobre a implementação.

Num código prático, os coeficientes que multiplicam os valores antigos do campo e os termos de curvatura são frequentemente pré-calculados e armazenados, especialmente quando \(\epsilon\), \(\mu\) ou \(\sigma\) variam de célula para célula. Por exemplo, a atualização do campo elétrico assume, geralmente, a seguinte forma:

\[ E^{n+1} = C_a E^n + C_b (\nabla \times H), \]

onde os coeficientes locais \(C_a\) e \(C_b\) codificam os parâmetros do material e o intervalo de tempo. Isso torna a implementação mais eficiente e evita o recálculo dos mesmos fatores a cada intervalo de tempo.

É também importante salientar que a derivação apresentada acima fornece apenas o esquema básico de atualização do método FDTD. Para uma simulação funcional, são também necessários elementos como fontes, condições de contorno, camadas de absorção, como as camadas perfeitamente adaptadas (PML), e um intervalo de tempo que satisfaça a condição de estabilidade de Courant. Nenhuma dessas adições altera a estrutura fundamental do método de «leapfrog» derivada aqui, mas cada uma delas modifica a forma como as equações de atualização são aplicadas em locais específicos da grelha.

Por fim, embora a notação apresentada acima tenha sido elaborada com o cuidado necessário para ser relevante na implementação, o posicionamento espacial exato de cada componente depende da convenção de grade utilizada no código. O requisito fundamental não é a tipografia em si, mas a consistência: o operador de curvatura, o espaçamento e os sinais devem corresponder ao mesmo arranjo de células de Yee em todos os pontos.