计算电磁学：FDFD算法总结

JiNan.YouQuan.Soft

已于 2025-01-26 17:59:44 修改

阅读量2k

点赞数 7

分类专栏： CAx HPC 文章标签：其他 1024程序员节

于 2024-05-20 17:24:55 首次发布

JiNan YouQuan Software Co. Ltd.

本文链接：https://blog.csdn.net/qq_26221775/article/details/139070586

版权

CAx 同时被 2 个专栏收录

90 篇文章

订阅专栏

HPC

22 篇文章

订阅专栏

FDFD(Finite Difference Frequency Domain)在Yee格子上使用有限差分求解频域Maxwell方程组。本文通过FDFP差分方程组的推导来阐述FDFP的原理。

注1：限于研究水平，分析难免不当，欢迎批评指正。

注2：文章内容会不定期更新。

0、预修

0.1 卷积分

对于函数 $f_{1}\left ( t \right )$ 、 $f_{2}\left ( t \right )$ ，则 $\int_{-\infty }^{+\infty}f_{1}\left ( \tau \right )f_{2}\left ( t-\tau \right )dt$ 称为函数函数 $f_{1}\left ( t \right )$ 、 $f_{2}\left ( t \right )$ 的卷积，记作 $f_{1}\left ( t \right )*f_{2}\left ( t \right )$ 。

交换律： $f_{1}\left ( t \right )*f_{2}\left ( t \right )=f_{2}\left ( t \right )*f_{1}\left ( t \right )$

0.2 Fourier变换

Fourier积分定理：若函数 $f\left ( t \right )$ 在 $\left ( -\infty ,+\infty \right )$ 上满足下列条件，(1). $f\left ( t \right )$ 在任一有限区间满足Dirichlet条件；(2). $f\left ( t \right )$ 在无限区间 $\left ( -\infty ,+\infty \right )$ 上绝对可积分(即 $\int_{-\infty }^{+\infty }\left |f\left ( t \right ) \right |dt < \infty$ )，则有

$f\left ( t \right )=\int_{-\infty }^{+\infty }\int_{-\infty }^{+\infty }f\left ( \tau \right )e^{-j\omega\tau}d\tau e^{-j\omega t}d\omega$

记 $F\left ( \omega \right )=\mathcal{F}\left ( f\left ( t \right ) \right ) =\int_{-\infty }^{+\infty }f\left ( t \right )e^{-j\omega t}dt$ 为函数 $f\left ( t \right )$ 的Fourier变换，对应的Fourier逆变换为 $f\left ( t\right )=\mathcal{F}^{-1}\left [ F\left ( \omega \right )\right ] =\frac{1}{2\pi }\int_{-\infty }^{+\infty }F\left ( \omega \right )e^{j\omega t}d\omega$ 。

微分性质：若 $\left | t \right |\rightarrow \infty$ 时， $\left | f\left (t \right )\right |\rightarrow \infty$ ，则有 $\mathcal{F}\left ( {f}'\left ( t \right ) \right )=j\omega \mathcal{F}\left ( f\left ( t \right ) \right )$

0.3 Laplace变换

若函数 $f\left ( t \right )$ 在 $\left [ 0,+\infty \right )$ 时有定义，而且积分 $\int_{0}^{+\infty }f\left ( t \right )e^{-st}dt$ 在复数 $s=\beta +j\omega$ 的某一域内存在，则记，

$F\left ( s \right )=\mathcal{L}\left ( f\left ( t \right ) \right )=\int_{0}^{+\infty }f\left ( t \right )e^{-st}dt$ 为函数 $f\left ( t \right )$ 的Laplace变换，

$f\left ( t \right )=\mathcal{L}^{-1}\left ( F\left ( s \right ) \right )=\frac{1}{2\pi j}\int_{\beta -j\infty }^{\beta +j\infty}F\left ( s \right )e^{st}ds$ 为 $F\left ( s \right )$ 的Laplace逆变换。

Laplace变换存在的充分条件：若函数 $f\left ( t \right )$ 在 $t\geq 0$ 时满足以下条件，(1) 在 $t\geq 0$ 任一有限区域分段连续；(2) 当 $t\rightarrow +\infty$ 时，存在常数 $M>0$ 及 $c\geq 0$ ，使得 $\left | f\left ( t \right ) \right |<Me^{ct}$ ，则 $f\left ( t \right )$ 的Laplace变换 $\int_{0}^{+\infty }f\left ( t \right )e^{-st}dt$ 在半平面 $Re\left ( s \right )>c$ 上一定存在，右端积分在 $Re\left ( s \right )\geqslant c_{1}>c$ 绝对收敛且一致收敛，并且在半平面 $Re\left ( s \right )>c$ 内， $F\left ( s \right )$ 为解析解。

微分性质：若 $F\left ( s \right )=\mathcal{L}\left ( f\left ( t \right ) \right )$ ，则有 $\mathcal{L}\left ( {f}'\left ( t \right ) \right )=sF\left ( s \right )-F\left ( 0 \right )$ 。

积分性质：若 $F\left ( s \right )=\mathcal{L}\left ( f\left ( t \right ) \right )$ ，则有 $\mathcal{L}\left ( \int_{0}^{t}f\left ( t \right ) dt\right )=\frac{1}{s}F\left ( s \right )$

一、数学模型

对于微分形式的Maxwell方程组，有

$\left\{\begin{matrix} \nabla \cdot \mathbf{D}=\rho\\ \nabla\cdot \mathbf{B}=0\\ \nabla\times \mathbf{H}=\boldsymbol{J}+\frac{\partial \boldsymbol{D}}{\partial t}\\ \nabla\times \mathbf{E}=-\frac{\partial \boldsymbol{B}}{\partial t}-\boldsymbol{M} \end{matrix}\right.$

依据Fourier逆变换，可得复数形式的Maxwell方程组，

$\left\{\begin{matrix} \nabla \cdot \mathbf{D}=\rho\\ \nabla\cdot \mathbf{B}=0\\ \nabla\times \mathbf{H}=\boldsymbol{J}+j\omega\boldsymbol{D}\\ \nabla\times \mathbf{E}=-j\omega\boldsymbol{B}-\boldsymbol{M} \end{matrix}\right.$

对于各向异性介质，则有

$\left\{\begin{matrix} \nabla \cdot \mathbf{D}=\rho\\ \nabla\cdot \mathbf{B}=0\\ \nabla\times \mathbf{H}=\boldsymbol{J}+j\omega\boldsymbol{\epsilon}\cdot \boldsymbol{E}\\ \nabla\times \mathbf{E}=-j\omega\boldsymbol{\mu}\cdot \boldsymbol{H}-\boldsymbol{M} \end{matrix}\right.\Rightarrow \left\{\begin{matrix} \boldsymbol{\mu}^{-1}\nabla^{2} \boldsymbol{E}+\omega ^{2}\boldsymbol{\epsilon}\cdot\boldsymbol{E}=\nabla\times \boldsymbol{\mu}^{-1}\cdot \nabla\times\boldsymbol{E}+\boldsymbol{\mu}^{-1}\nabla\nabla\cdot \boldsymbol{E}+j\omega\boldsymbol{J}+\nabla\times\boldsymbol{\mu}^{-1}\cdot \boldsymbol{M}\\ \boldsymbol{\epsilon}^{-1}\nabla^{2} \boldsymbol{H}+\omega ^{2}\boldsymbol{\mu}\cdot\boldsymbol{H}=\nabla\times \boldsymbol{\epsilon}^{-1}\cdot \nabla\times\boldsymbol{H}+\boldsymbol{\epsilon}^{-1}\nabla\nabla\cdot \boldsymbol{H}+j\omega\boldsymbol{M}-\nabla\times\boldsymbol{\epsilon}^{-1}\cdot \boldsymbol{J}\end{matrix}\right.$

不考虑电流密度 $\boldsymbol{J}$ 与磁流密度 $\boldsymbol{M}$ ，在均匀各向同性介质中，有

$\left\{\begin{matrix}\nabla\times \mathbf{H}=j\omega \epsilon \boldsymbol{E} \\ \nabla\times\mathbf{E}=-j\omega \mu \boldsymbol{H} \end{matrix}\right.\Rightarrow \left\{\begin{matrix} \nabla^{2}\boldsymbol{H}+\omega ^{2}\epsilon\mu\boldsymbol{H}=0\\ \nabla^{2}\boldsymbol{E}+\omega ^{2}\epsilon\mu\boldsymbol{E}=\frac{\nabla\rho }{\epsilon}\end{matrix}\right.$

在笛卡尔直角坐标系下，有

$\left\{\begin{matrix} \frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}=j\omega \epsilon E_{x}\\ \frac{\partial H_{x}}{\partial z}-\frac{\partial H_{z}}{\partial x}=j\omega \epsilon E_{y}\\ \frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}=j\omega \epsilon E_{z}\end{matrix}\right.$ $\left\{\begin{matrix} \frac{\partial E_{z}}{\partial y}-\frac{\partial E_{y}}{\partial z}=-j\omega \mu H_{x}\\ \frac{\partial E_{x}}{\partial z}-\frac{\partial E_{z}}{\partial x}=-j\omega \mu H_{y}\\ \frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}=-j\omega \mu H_{z}\end{matrix}\right.$

二、数值模型：控制方程的空间离散

FDFD算法也采用Yee格子在空间上进行离散。在Yee格子中，电场分量放置到Yee单元各棱的中间，方向平行于各棱；磁场分量放置到Yee单元各面中心，方向平行于各面法线。

$\left\{\begin{matrix} \frac{H_{z}\left (i+\frac{1}{2},j+\frac{1}{2},k \right )-H_{z}\left (i+\frac{1}{2},j-\frac{1}{2},k \right )}{\Delta y}-\frac{H_{y}\left (i+\frac{1}{2},j,k+\frac{1}{2} \right )-H_{y}\left (i+\frac{1}{2},j,k-\frac{1}{2} \right )}{\Delta z}=j\omega \epsilon E_{x}\left ( i+\frac{1}{2},j,k \right )\\ \frac{H_{x}\left (i,j+\frac{1}{2},k+\frac{1}{2} \right )-H_{x}\left (i,j+\frac{1}{2},k-\frac{1}{2} \right )}{\Delta z}-\frac{H_{z}\left (i+\frac{1}{2},j+\frac{1}{2},k \right )-H_{z}\left (i-\frac{1}{2},j+\frac{1}{2},k \right )}{\Delta x}=j\omega \epsilon E_{y}\left ( i,j+\frac{1}{2},k \right )\\ \frac{H_{y}\left (i+\frac{1}{2},j,k+\frac{1}{2} \right )-H_{y}\left (i-\frac{1}{2},j,k+\frac{1}{2} \right )}{\Delta x}-\frac{H_{x}\left (i,j+\frac{1}{2},k+\frac{1}{2} \right )-H_{x}\left (i,j-\frac{1}{2},k+\frac{1}{2} \right )}{\Delta y}=j\omega \epsilon E_{z}\left ( i,j,k+\frac{1}{2} \right ) \end{matrix}\right.$

$\left\{\begin{matrix} \frac{E_{z}\left (i,j+1,k+\frac{1}{2} \right )-E_{z}\left (i,j,k+\frac{1}{2} \right )}{\Delta y}-\frac{E_{y}\left (i,j+\frac{1}{2},k+1 \right )-E_{y}\left (i,j+\frac{1}{2},k \right )}{\Delta z}=-j\omega \mu H_{x}\left ( i,j+\frac{1}{2},k+\frac{1}{2} \right )\\ \frac{E_{x}\left (i+\frac{1}{2},j,k+1 \right )-E_{x}\left (i+\frac{1}{2},j,k \right )}{\Delta z}-\frac{E_{z}\left (i+1,j,k+\frac{1}{2}\right )-E_{z}\left (i,j,k+\frac{1}{2} \right )}{\Delta x}=-j\omega \mu H_{y}\left ( i+\frac{1}{2},j,k+\frac{1}{2} \right )\\ \frac{E_{y}\left (i+1,j+\frac{1}{2},k \right )-E_{y}\left (i,j+\frac{1}{2},k \right )}{\Delta x}-\frac{E_{x}\left (i+\frac{1}{2},j+1,k\right )-E_{x}\left (i+\frac{1}{2},j,k \right )}{\Delta y}=-j\omega \mu H_{z}\left ( i+\frac{1}{2},j+\frac{1}{2},k \right )\end{matrix}\right.$

设 $n_{x}$ 、 $n_{y}$ 、 $n_{z}$ 分别表示 $x$ 、 $y$ 、 $z$ 方向上的Yee格子数目，考虑将上述差分方程组写成矩阵方程组 $\boldsymbol{A}\boldsymbol{X}$ 形式，则有，

$n=kn_{x}n_{y}+jn_{x}+i$

$A\left ( 6n,6n+4-6n_{x}n_{y} \right )=-\frac{1}{\Delta z}$ ， $A\left ( 6n,6n+5-6n_{x} \right )=\frac{1}{\Delta y}$ ， $A\left ( 6n,6n \right )=j\omega \epsilon$ ， $A\left ( 6n,6n+1 \right )=0$ ， $A\left ( 6n,6n+2 \right )=0$ ， $A\left ( 6n,6n+3 \right )=0$ ， $A\left ( 6n,6n+4 \right )=\frac{1}{\Delta z}$ ， $A\left ( 6n,6n+5 \right )=-\frac{1}{\Delta y}$

$A\left ( 6n+1,6n+3-6n_{x}n_{y} \right )=\frac{1}{\Delta z}$ ， $A\left ( 6n+1,6n-1\right )=-\frac{1}{\Delta x}$ ， $A\left ( 6n+1,6n \right )=0$ ， $A\left ( 6n+1,6n+1 \right )=j\omega \epsilon$ ， $A\left ( 6n+1,6n+2 \right )=0$ ， $A\left ( 6n+1,6n+3 \right )=-\frac{1}{\Delta z}$ ， $A\left ( 6n+1,6n+5 \right )=\frac{1}{\Delta x}$

$A\left ( 6n+2,6n+3-6n_{x} \right )=-\frac{1}{\Delta y}$ ， $A\left ( 6n+2,6n-2 \right )=\frac{1}{\Delta x}$ ， $A\left ( 6n+2,6n \right )=0$ ， $A\left ( 6n+2,6n+1 \right )=0$ ， $A\left ( 6n+2,6n+2 \right )=j\omega \epsilon$ ， $A\left ( 6n+2,6n+3 \right )=\frac{1}{\Delta y}$ ， $A\left ( 6n+2,6n+4 \right )=-\frac{1}{\Delta x}$

$A\left ( 6n+3,6n \right )=0$ ， $A\left ( 6n+3,6n+1 \right )=\frac{1}{\Delta z}$ ， $A\left ( 6n+3,6n+2 \right )=\frac{1}{\Delta y}$ ， $A\left ( 6n+3,6n+3 \right )=j\omega \mu$ ， $A\left ( 6n+3,6n+4\right )=0$ ， $A\left ( 6n+3,6n+5\right )=0$ ， $A\left ( 6n+3,6n+1+6n_{x} \right )=\frac{1}{\Delta y}$ ， $A\left ( 6n+3,6n+1+6n_{x}n_{y} \right )=-\frac{1}{\Delta z}$

$A\left ( 6n+4,6n\right )=-\frac{1}{\Delta z}$ ， $A\left ( 6n+4,6n+2 \right )=\frac{1}{\Delta x}$ ， $A\left ( 6n+4,6n+3 \right )=0$ ， $A\left ( 6n+4,6n+4 \right )=j\omega \mu$ ， $A\left ( 6n+4,6n+5 \right )=0$ ， $A\left ( 6n+4,6n+8 \right )=-\frac{1}{\Delta x}$ ， $A\left ( 6n+4,6n+6n_{x}n_{y}\right )=\frac{1}{\Delta z}$

$A\left ( 6n+5,6n \right )=\frac{1}{\Delta y}$ ， $A\left ( 6n+5,6n+2\right )=-\frac{1}{\Delta x}$ ， $A\left ( 6n+5,6n+3 \right )=0$ ， $A\left ( 6n+5,6n+4 \right )=0$ ， $A\left ( 6n+5,6n+5 \right )=j\omega \mu$ ， $A\left ( 6n+5,6n+8\right )=\frac{1}{\Delta x}$ ， $A\left ( 6n+5,6n+6n_{y} \right )=-\frac{1}{\Delta y}$