青岛房产网站,公司网站制作计入什么科目,建筑网站推荐知乎,深圳 德 网站建设基于时域有限差分法的FDTD的计算电磁学算法#xff08;含Matlab代码#xff09;-YEE网格下的更新公式推导
参考书籍#xff1a;The finite-difference time-domain method for electromagnetics with MATLAB simulations#xff08;国内翻译版本#xff1a;MATLAB模拟的电…基于时域有限差分法的FDTD的计算电磁学算法含Matlab代码-YEE网格下的更新公式推导
参考书籍The finite-difference time-domain method for electromagnetics with MATLAB simulations国内翻译版本MATLAB模拟的电磁学时域有限差分法 代码推荐The finite-difference time-domain method for electromagnetics with MATLAB simulations的附件代码 我最初也是基于这个代码学习的
FDTD算法采用差分直接离散时域Maxwell方程电磁场的求解基于时间步的迭代无需存储全空间的电磁场信息内存消耗较小同时采用立方体网格和差分算法网格形式和算法均十分简单计算速度快基于时域算法特别适合“宽带问题”的求解。但是简单的立方体方体网格带来的弊端就是模型拟合精度较低对于含有精细结构的模型计算精度较低同时基于“微分方程”计算区域需要设置截断。 详细对比参考常用计算电磁学算法特性与电磁软件分析
1、从麦克斯韦开始的FDTD时域有限差分法
1.1 麦克斯韦方程
FDTD叫时域有限差分法显然其依赖的麦克斯韦方程也是时域的。麦克斯韦时域微分方程为 ∇ × H ∂ D ∂ t J ∇ × E − ∂ B ∂ t − M ∇ ⋅ D ρ e ∇ ⋅ B ρ m \begin{gathered} \nabla\times \mathbf{H} {\frac{\partial \mathbf{D}}{\partial t}}\boldsymbol{J} \\ \nabla\times \mathbf{E}-{\frac{\partial \mathbf{B}}{\partial t}}-\mathbf{M} \\ \nabla\cdot\mathbf{D}\rho_{\mathrm{e}} \\ \nabla\cdot \mathbf{B}\rho_{m} \end{gathered} ∇×H∂t∂DJ∇×E−∂t∂B−M∇⋅Dρe∇⋅Bρm 式中E为电场强度(V/m);D为电位移(C/m);H为磁场强度(A/m);B为磁通量密度(Wb/m°);J为电流密度(A/m);M为磁流密度(V/m); ρ e \rho_{e} ρe为电荷密度(C/m); ρ m \rho_{m} ρm为磁荷密度(Wb/m)。
依稀记得当时老师说麦克斯韦方程有其直观理解分别是 1. 变化的电场和电流会产生磁场 2. 变化的磁场和磁荷会产生电场自然界无磁荷一般是等效出来 3. 电流源产生电场 4. 磁流源产生磁场
1.2 本构关系
本构关系对补充麦克斯韦方程和描述媒质的特性是必要的,本构关系对线性、各向同性和非色散媒质可以写成 D ε E B μ H . \begin{aligned}D\varepsilon E\\B\mu H\end{aligned}. DBεEμH. 其中 ε \varepsilon ε为媒质的介电常数; μ \mu μ为媒质的磁导率。在自由空间有: ε ε 0 8.854 × 1 0 − 12 F / m μ μ 0 4 π × 1 0 − 7 H / m \begin{aligned}\varepsilon\varepsilon_08.854\times10^{-12}\quad\mathrm{F/m}\\\mu\mu_04\pi\times10^{-7}\quad\mathrm{H/m}\end{aligned} εμε08.854×10−12F/mμ04π×10−7H/m 在常规的电磁学表述中我们更多的使用相对介电常数。比如说耳熟能详的FR4板材其相对介电常数大概是 ε r 4.2 \varepsilon_r4.2 εr4.2。 这就代表其实际的介电常数为 ε F R 4 ε r ε 0 \varepsilon_{FR4}\varepsilon_r\varepsilon_0 εFR4εrε0。但是还有一个重要参数和本构关系相关那就是损耗角正切 t a n δ tan \delta tanδ。
对于FR4板材一般认为其损耗角正切为 t a n δ 0.02 tan \delta0.02 tanδ0.02根据微波工程1.3小节的公式 ϵ ϵ ′ − j ϵ ′ ′ ϵ ′ ( 1 − j tan δ ) ϵ 0 ϵ r ( 1 − j tan δ ) \epsilon\epsilon^{\prime}-j\epsilon^{\prime\prime}\epsilon^{\prime}(1-j\tan\delta)\epsilon_{0}\epsilon_{r}(1-j\tan\delta) ϵϵ′−jϵ′′ϵ′(1−jtanδ)ϵ0ϵr(1−jtanδ)其对应的介电常数应该是 ε F R 4 ε r ( 1 − j tan δ ) ε 0 ( 4.2 − j 0.02 ) ε 0 \varepsilon_{FR4}\varepsilon_r(1-j\tan\delta)\varepsilon_0(4.2-j0.02)\varepsilon_0 εFR4εr(1−jtanδ)ε0(4.2−j0.02)ε0 其对应的相对介电常数为4.2-j0.02
在进行FDTD的推导时因为在 FDTD 的更新方程的过程中满足散度方程所以只需要考虑两个旋度方程即可。麦克斯韦中的电流密度 J \boldsymbol{J} J等于导体电流密度 J c \boldsymbol{J_c} Jc与施加电流密度 J i \boldsymbol{J_i} Ji之和即 J J c J i \boldsymbol{J}\boldsymbol{J_{\mathrm{c}}}\boldsymbol{J_{\mathrm{i}}} JJcJi 对于磁流密度也类似 M M c M i \boldsymbol{M}\boldsymbol{M_{\mathrm{c}}}\boldsymbol{M_{\mathrm{i}}} MMcMi 因此对原来的麦克斯韦方程拆分一下就是 ∇ × H ε ∂ E ∂ t σ e E J i \nabla\times \boldsymbol{H}\varepsilon\frac{\partial \boldsymbol{E}}{\partial t}\sigma^{e}\boldsymbol{E}\boldsymbol{J_{i}} ∇×Hε∂t∂EσeEJi 和 ∇ × E − μ ∂ H ∂ t − σ m H − M i \nabla\times \boldsymbol{E}-\mu\frac{\partial \boldsymbol{H}}{\partial t}-\sigma^{m}\boldsymbol{H}-\boldsymbol{M_{i}} ∇×E−μ∂t∂H−σmH−Mi
旋度的计算公式大家还记得不 ∇ × F ( x , y , z ) ∣ i ^ j ^ k ^ ∂ ∂ x ∂ ∂ y ∂ ∂ z F x F y F z ∣ ( ∂ F z ∂ y − ∂ F y ∂ z ) i ^ ( ∂ F x ∂ z − ∂ F z ∂ x ) j ^ ( ∂ F y ∂ x − ∂ F x ∂ y ) k ^ \begin{aligned} \nabla\times\mathbf{F}(x,y,z)\begin{vmatrix}\hat{\boldsymbol{i}}\hat{\boldsymbol{j}}\hat{\boldsymbol{k}}\\\frac{\partial}{\partial x}\frac{\partial}{\partial y}\frac{\partial}{\partial z}\\F_xF_yF_z\end{vmatrix} \\ \left(\frac{\partial F_z}{\partial y}-\frac{\partial F_y}{\partial z}\right)\hat{\boldsymbol{i}}\left(\frac{\partial F_x}{\partial z}-\frac{\partial F_z}{\partial x}\right)\hat{\boldsymbol{j}}\left(\frac{\partial F_y}{\partial x}-\frac{\partial F_x}{\partial y}\right)\hat{\boldsymbol{k}} \end{aligned} ∇×F(x,y,z) i^∂x∂Fxj^∂y∂Fyk^∂z∂Fz (∂y∂Fz−∂z∂Fy)i^(∂z∂Fx−∂x∂Fz)j^(∂x∂Fy−∂y∂Fx)k^ 把麦克斯韦旋度方程按照三个方向xyz全部展开就可以得到6个方程 ∂ E x ∂ t 1 ε x ( ∂ H z ∂ y − ∂ H y ∂ z − σ x e E x − J i x ) ∂ E y ∂ t 1 ε y ( ∂ H x ∂ z − ∂ H z ∂ x − σ y e E y − J i y ) ∂ E z ∂ t 1 ε z ( ∂ H y ∂ x − ∂ H x ∂ y − σ z e E z − J i z ) ∂ H x ∂ t 1 μ x ( ∂ E y ∂ z − ∂ E z ∂ y − σ x m H x − M i x ) ∂ H y ∂ t 1 μ y ( ∂ E x ∂ x − ∂ E x ∂ z − σ y m H y − M i y ) ∂ H z ∂ t 1 μ z ( ∂ E x ∂ y − ∂ E y ∂ x − σ z m H z − M i z ) \begin{gathered} \frac{\partial\boldsymbol{E}_x}{\partial t} \frac1{\varepsilon_x}\Big(\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z}-\sigma_x^eE_x-J_{ix}\Big) \\ \frac{\partial E_y}{\partial t} \frac1{\varepsilon_y}\Big(\frac{\partial H_x}{\partial z}-\frac{\partial H_z}{\partial x}-\sigma_y^eE_y-J_{iy}\Big) \\ \frac{\partial E_z}{\partial t} \frac{1}{\varepsilon_{z}}\Big(\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}-\sigma_{z}^{e}E_{z}-J_{iz}\Big) \\ \frac{\partial H_x}{\partial t} \frac1{\mu_x}\Big(\frac{\partial E_y}{\partial z}-\frac{\partial E_z}{\partial y}-\sigma_x^mH_x-M_{ix}\Big) \\ \frac{\partial H_y}{\partial t} \frac1{\mu_y}\Big(\frac{\partial\boldsymbol{E}_x}{\partial x}-\frac{\partial\boldsymbol{E}_x}{\partial\boldsymbol{z}}-\boldsymbol{\sigma}_y^\mathfrak{m}H_y-\boldsymbol{M}_{iy}\Big) \\ \frac{\partial H_z}{\partial t} \frac{1}{\mu_{z}}\Big(\frac{\partial\boldsymbol{E}_{x}}{\partial y}-\frac{\partial\boldsymbol{E}_{y}}{\partial x}-\sigma_{z}^{\mathfrak{m}}H_{z}-\boldsymbol{M}_{iz}\Big) \end{gathered} ∂t∂Exεx1(∂y∂Hz−∂z∂Hy−σxeEx−Jix)∂t∂Eyεy1(∂z∂Hx−∂x∂Hz−σyeEy−Jiy)∂t∂Ezεz1(∂x∂Hy−∂y∂Hx−σzeEz−Jiz)∂t∂Hxμx1(∂z∂Ey−∂y∂Ez−σxmHx−Mix)∂t∂Hyμy1(∂x∂Ex−∂z∂Ex−σymHy−Miy)∂t∂Hzμz1(∂y∂Ex−∂x∂Ey−σzmHz−Miz)
2、空间差分与时间差分
2.1、非常简单的差分方程
FDTD是在离散网格中进行迭代的上面的麦克斯韦公式有大量的求导计算这该如何解决呢答案是差分近似。大家学高数都学过导数的近似吧 f ′ ( x ) lim Δ x → 0 f ( x Δ x ) − f ( x ) Δ x f^{}(x)\underset{\Delta x\to0}{\operatorname*{lim}}\frac{f(x\Delta x)-f(x)}{\Delta x} f′(x)Δx→0limΔxf(xΔx)−f(x) 如果 Δ x \Delta x Δx非常小那么 f ′ ( x ) ≈ f ( x Δ x ) − f ( x ) Δ x f^{}(x)\approx\frac{f(x\Delta x)-f(x)}{\Delta x} f′(x)≈Δxf(xΔx)−f(x) 但是为了实现更高的精度所以采用FDTD都会采用双向差分公式 f ′ ( x ) ≈ f ( x Δ x ) − f ( x − Δ x ) 2 Δ x f^{^{\prime}}(x){\approx}\frac{f(x\Delta x)-f(x-\Delta x)}{2\Delta x} f′(x)≈2Δxf(xΔx)−f(x−Δx) 实际上此处使用的是近似也存在高阶的FDTD的算法对于此近似考虑了更多项精度会更高(参考“基于高阶时域有限差分法平面波及完全匹配层的研究”等) f ′ ( x ) f ( x Δ x ) − f ( x − Δ x ) 2 Δ x − ( Δ x 2 ) 6 . . . f ( x Δ x ) − f ( x − Δ x ) 2 Δ x O ( ( Δ x ) 2 ) f^{\prime}(x)\frac{f(x\Delta x)-f(x-\Delta x)}{2\Delta x}-\frac{(\Delta x^{2})}{6}...\frac{f(x\Delta x)-f(x-\Delta x)}{2\Delta x}O((\Delta x)^{2}) f′(x)2Δxf(xΔx)−f(x−Δx)−6(Δx2)...2Δxf(xΔx)−f(x−Δx)O((Δx)2)
2.2、差分方程的运用
在FDTD算法中网格被剖分为YEE网格的形式电场和磁场元胞差半个身位其更新的时间步也是差 0.5 Δ t 0.5\Delta t 0.5Δt 具体来讲实际的电场网格和磁场网格的位置是 E x ( i , j , k ) ⇒ ( ( i − 0 , 5 ) Δ x , ( j − 1 ) Δ y , ( k − 1 ) Δ z ) E y ( i , j , k ) ⇒ ( ( i − 1 ) Δ x , ( j − 0.5 ) Δ y , ( k − 1 ) Δ z ) E z ( i , j , k ) ⇒ ( ( i − 1 ) Δ x , ( j − 1 ) Δ y , ( k − 0.5 ) Δ z ) H x ( i , j , k ) ⇒ ( ( i − 1 ) Δ x , ( j − 0.5 ) Δ y , ( k − 0.5 ) Δ z ) H y ( i , j , k ) ⇒ ( ( i − 0.5 ) Δ x , ( j − 1 ) Δ y , ( k − 0.5 ) Δ z ) H z ( i , j , k ) ⇒ ( ( i − 0.5 ) Δ x , ( j − 0.5 ) Δ y , ( k − 1 ) Δ z ) \begin{aligned} E_x(i,j,k)\Rightarrow\left((i-0,5)\Delta x,(j-1)\Delta y,(k-1)\Delta z\right)\\ E_y(i,j,k)\Rightarrow\left((i-1)\Delta x,(j-0.5)\Delta y,(k-1)\Delta z\right)\\ E_z(i,j,k)\Rightarrow\left((i-1)\Delta x,(j-1)\Delta y,(k-0.5)\Delta z\right)\\ H_x(i,j,k)\Rightarrow\left((i-1)\Delta x,(j-0.5)\Delta y,(k-0.5)\Delta z\right)\\ H_y(i,j,k)\Rightarrow((i-0.5)\Delta x,(j-1)\Delta y,(k-0.5)\Delta z) \\ H_z(i,j,k)\Rightarrow((i-0.5)\Delta x,(j-0.5)\Delta y,(k-1)\Delta z) \end{aligned} Ex(i,j,k)⇒((i−0,5)Δx,(j−1)Δy,(k−1)Δz)Ey(i,j,k)⇒((i−1)Δx,(j−0.5)Δy,(k−1)Δz)Ez(i,j,k)⇒((i−1)Δx,(j−1)Δy,(k−0.5)Δz)Hx(i,j,k)⇒((i−1)Δx,(j−0.5)Δy,(k−0.5)Δz)Hy(i,j,k)⇒((i−0.5)Δx,(j−1)Δy,(k−0.5)Δz)Hz(i,j,k)⇒((i−0.5)Δx,(j−0.5)Δy,(k−1)Δz)
更新的时间步也是差 0.5 Δ t 0.5\Delta t 0.5ΔtFDTD算法在离散的时间瞬间取样和计算场值,但是电场和磁场取样计算并不是在相同的时刻。对时间步 Δ t \Delta t Δt,电场E的取样时刻为:0, Δ t \Delta t Δt,2 Δ t \Delta t Δt,3 Δ t \Delta t Δt…,n Δ t \Delta t Δt而磁场H取样时刻为:0.5 Δ t \Delta t Δt,1.5 Δ t \Delta t Δt,2.5 Δ t \Delta t Δt,…(n0.5) Δ t \Delta t Δt。即电场取样在时间的整数步长时刻,而磁场取样时刻为半整数时间步时刻。它们之间的时间差为半个时间步。
因此考虑一个上面得到的麦克斯韦的方程以Ex方向为例 ∂ E x ∂ t 1 ε x ( ∂ H z ∂ y − ∂ H y ∂ z − σ x e E x − J i r ) \frac{\partial E_x}{\partial t}\frac1{\varepsilon_x}\left(\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z}-\sigma_x^eE_x-J_{ir}\right) ∂t∂Exεx1(∂y∂Hz−∂z∂Hy−σxeEx−Jir) 观察其导数项分别有时间的差分项 ∂ E x ∂ t \frac{\partial E_x}{\partial t} ∂t∂Ex和空间的差分项 ∂ H z ∂ y \frac{\partial H_z}{\partial y} ∂y∂Hz和 ∂ H y ∂ z \frac{\partial H_y}{\partial z} ∂z∂Hy。
方程中的导数可以用中心差分来近似此时 E x n ( i , j , k ) E_x^n(i,j,k) Exn(i,j,k)的位置为中心差分公式的中心点而时间上应以 ( n 0.5 ) Δ t (n0.5)\Delta t (n0.5)Δt作为中心点因为电场E的取样时刻为:0, Δ t \Delta t Δt,2 Δ t \Delta t Δt,3 Δ t \Delta t Δt…,n Δ t \Delta t Δt而 ( n 0.5 ) Δ t (n0.5)\Delta t (n0.5)Δt差分后可以得到n和n1符合取样时刻。因此第一项 ∂ E x ∂ t \frac{\partial E_x}{\partial t} ∂t∂Ex可以写成如下的差分形式 E x n 0.5 ( i , j , k ) E x n 1 ( i , j , k ) − E x n ( i , j , k ) Δ t E_x^{n0.5}(i,j,k)\frac{E_x^{n1}(i,j,k)-E_x^n(i,j,k)}{\Delta t} Exn0.5(i,j,k)ΔtExn1(i,j,k)−Exn(i,j,k) 而空间的差分项 ∂ H z ∂ y \frac{\partial H_z}{\partial y} ∂y∂Hz可以写成 ∂ H z ∂ y H z n 1 2 ( i , j , k ) − H z n 1 2 ( i , j − 1 , k ) Δ y \frac{\partial H_z}{\partial y}\frac{H_z^{n\frac12}(i,j,k)-H_z^{n\frac12}(i,j-1,k)}{\Delta y} ∂y∂HzΔyHzn21(i,j,k)−Hzn21(i,j−1,k)
2.3、得到差分方程
把所有项都写成差分形式就可以得到3D的FDTD更新方程 E x n 1 ( i , j , k ) C e x e ( i , j , k ) × E x n ( i , j , k ) C e x h z ( i , j , k ) × ( H z n 1 2 ( i , j , k ) − H z n 1 2 ( i , j − 1 , k ) ) C e x h y ( i , j , k ) × ( H y n 1 2 ( i , j , k ) − H y n 1 2 ( i , j , k − 1 ) ) C e x j ( i , j , k ) × J i x n 1 2 ( i , j , k ) \begin{aligned} E_{x}^{n1}\left(i,j,k\right) C_{exe}(i,j,k)\times E_x^n(i,j,k) \\ C_{exhz}(i,j,k)\times(H_{z}^{n\frac12}(i,j,k)-H_{z}^{n\frac12}(i,j-1,k)) \\ C_{\mathrm{exhy}}(i,j,k)\times(H_y^{n\frac12}(i,j,k)-H_y^{n\frac12}(i,j,k-1)) \\ C_{exj}\left(i,j,k\right)\times J_{ix}^{n\frac12}(i,j,k) \end{aligned} Exn1(i,j,k)Cexe(i,j,k)×Exn(i,j,k)Cexhz(i,j,k)×(Hzn21(i,j,k)−Hzn21(i,j−1,k))Cexhy(i,j,k)×(Hyn21(i,j,k)−Hyn21(i,j,k−1))Cexj(i,j,k)×Jixn21(i,j,k) C开头的都是系数为了书写方便其实际的值为 C e x e ( i , j , k ) 2 ε z ( i , j , k ) − Δ t σ z e ( i , j , k ) 2 ε z ( i , j , k ) Δ t σ z e ( i , j , k ) C e x h y ( i , j , k ) 2 Δ t ( 2 ε z ( i , j , k ) Δ t σ z e ( i , j , k ) ) Δ x C e x h y ( i , j , k ) − 2 Δ t ( 2 ε z ( i , j , k ) Δ t σ z e ( i , j , k ) ) Δ y C e x j ( i , j , k ) − 2 Δ t 2 ε z ( i , j , k ) Δ t σ z e ( i , j , k ) \begin{gathered} C_{exe}(i,j,k) \frac{2\varepsilon_z(i,j,k)-\Delta t\sigma_z^e(i,j,k)}{2\varepsilon_z(i,j,k)\Delta t\sigma_z^e(i,j,k)} \\ C_{exhy}(i,j,k) \frac{2\Delta t}{(2\varepsilon_z(i,j,k)\Delta t\sigma_z^e(i,j,k))\Delta x} \\ C_{{exhy}}(i,j,k) -\frac{2\Delta t}{(2\varepsilon_z(i,j,k)\Delta t\sigma_z^e(i,j,k))\Delta y} \\ C_{exj}\left(i,j,k\right) -\frac{2\Delta t}{2\varepsilon_z(i,j,k)\Delta t\sigma_z^e(i,j,k)} \end{gathered} Cexe(i,j,k)2εz(i,j,k)Δtσze(i,j,k)2εz(i,j,k)−Δtσze(i,j,k)Cexhy(i,j,k)(2εz(i,j,k)Δtσze(i,j,k))Δx2ΔtCexhy(i,j,k)−(2εz(i,j,k)Δtσze(i,j,k))Δy2ΔtCexj(i,j,k)−2εz(i,j,k)Δtσze(i,j,k)2Δt 当然这只是6个方程中的一个更加详细的方程参考 MATLAB模拟的电磁学时域有限差分法的1.3。看看对应的matlab代码是怎么写的没有电流就可以省略Cexj
current_time current_time dt/2;Ex(1:nx,2:ny,2:nz) Cexe(1:nx,2:ny,2:nz).*Ex(1:nx,2:ny,2:nz) ... Cexhz(1:nx,2:ny,2:nz).*...(Hz(1:nx,2:ny,2:nz)-Hz(1:nx,1:ny-1,2:nz)) ... Cexhy(1:nx,2:ny,2:nz).*...(Hy(1:nx,2:ny,2:nz)-Hy(1:nx,2:ny,1:nz-1)); % General electric field updating coefficients
% Coeffiecients updating Ex
Cexe (2*eps_r_x*eps_0 - dt*sigma_e_x) ..../(2*eps_r_x*eps_0 dt*sigma_e_x);
Cexhz (2*dt/dy)./(2*eps_r_x*eps_0 dt*sigma_e_x);
Cexhy -(2*dt/dz)./(2*eps_r_x*eps_0 dt*sigma_e_x);
