CN107944214A - 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 - Google Patents
笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 Download PDFInfo
- Publication number
- CN107944214A CN107944214A CN201711201147.7A CN201711201147A CN107944214A CN 107944214 A CN107944214 A CN 107944214A CN 201711201147 A CN201711201147 A CN 201711201147A CN 107944214 A CN107944214 A CN 107944214A
- Authority
- CN
- China
- Prior art keywords
- spherical
- electric field
- moment
- intermediate variable
- magnetic field
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16Z—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
- G16Z99/00—Subject matter not provided for in other main groups of this subclass
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，该方法是建立在笛卡尔坐标系下的，三维问题中各向异性完全匹配层的外形为球形，二维问题中各向异性完全匹配层的外形为圆形，可以保证网格大小的一致，不会有不稳定的因素出现，又能截断一些多余的网格，减少计算量，提高计算效率，克服了现有笛卡尔坐标系下的方形各向异性完全匹配层计算效率低、计算复杂的问题，将球形的外形和方形的网格相结合，既稳定效率又高。通过设计笛卡尔坐标系下的球形各向异性完全匹配层截断边界的参数，用于计算电磁学的程序中，可实现在有限的计算区域达到模拟“微波暗室”吸波材料的效果。
Description
技术领域
本发明属于计算电磁学技术领域，具体涉及一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法。
背景技术
电磁场数值计算方法广泛应用于微波电路、天线设计、目标散射计算和电磁兼容等方面的研究，其中时域有限差分方法是最常用的方法之一。对于计算机编程求解电磁场问题，有限的内存是无法模拟无限大的计算区域的，因此需要在计算区域边缘设置吸收层作为截断边界，使电磁波在边界上被吸收，其性能与电磁场数值计算的精度和效率是紧密相关的。Stephen D.Gendney于1996年提出了各向异性完全匹配层，可以实现较好的吸收效果，但是固有的立方体网格使得计算区域为立方体形或者正方形，在角区域和棱区域处理的过程中浪费了大量的计算机资源，并且会造成不必要的反射。专利号为201210177288.0的中国专利公开了一种二维柱坐标系完全吸收边界的实现方法，省略了棱区域的计算，但其中涉及过多繁杂的数学表达式，给计算机编程带来了不小的难度。此外，由于其计算区域的形状是圆柱形的，网格形状也是弯曲的，网格的尺寸与半径成正比，因此随着半径的增大，越外侧的网格尺寸也会越大，从而导致计算精度大幅降低，另外网格的尺寸也有可能超出“Courant稳定性条件”的限制导致算法的不稳定。申请号为201410568490.5的中国专利公开了一种阻抗匹配层的截断边界，该截断边界采用阻抗匹配形式进行截断，该过程中所涉及的参数较多，计算十分复杂，在弯曲处理截断边界时，过多的参数将导致计算机程序变得十分冗余，影响计算效率，另一方面，由于阻抗匹配层的自身缺陷，可适用的范围较窄，不能满足使用需求。然而，目前尚未出现建立在笛卡尔坐标系下弯曲形状的各项异性完全匹配层截断边界的实现方法。
发明内容
本发明的目的是提供一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，在笛卡尔坐标系下对传统的立方体形或者正方形各向异性完全匹配层截断边界进行弯曲处理，有效的避免了计算那些角区域和棱区域，在吸收向外传播的电磁波的同时，实现大幅减少计算量，提高计算效率。
本发明所采用的技术方案是一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，该方法的步骤为：
1)建立求解对象的模型数据和时域有限差分法的计算空间；
向计算机申请内存空间，设Xn,Yn,Zn分别为x,y,z方向上的最小位置，其中Xn＝Yn＝Zn，Xp,Yp,Zp分别为x,y,z方向上的最大位置，其中Xp＝Yp＝Zp，总体计算的区域范围为(Xn,Yn,Zn)→(Xp,Yp,Zp)，呈立方体形，区域范围的大小为(Xp–Xn)×(Yp–Yn)×(Zp–Zn)，设x,y,z方向上的空间步长为Δx,Δy和Δz，且Δx＝Δy＝Δz，构成立方体网格区域；设球形各向异性完全匹配层的直径为a，所对应的网格数量为a/Δx＝Xp–Xn，中心坐标为((Xp–Xn)/2,(Yp–Yn)/2,(Zp–Zn)/2)，用球形的边界截断立方体网格区域，使计算区域成为球形，设定时间步长为Δt，并对球形内部区域的网格做inner标记，对球形外部区域的网格做outer标记，以区分不同的区域；设定电磁仿真迭代步数为N，n表示时间步，n的范围为1→N，选择时谐场的点源和电流源分别作为二维和三维问题的激励源，激励源函数用Einc表示；
2)构建球形各向异性完全匹配层截断边界；
在步骤1)所设置的标记为inner的球形内部区域，设定δ为球形各向异性完全匹配层截断边界的层数，将标记为inner的球形内部区域分为两部分，包括半径从0到a/2－δ*Δx的区域，称为自由空间，做inner-free标记，半径从a/2－δ*Δx到a/2的区域，即为球形各向异性完全匹配层截断边界，做inner-CSUPML标记；在球形各向异性完全匹配层截断边界的最内层，即与球心距离为a/2－δ*Δx的位置，电导率σ＝0，在球形各向异性完全匹配层截断边界的最外层，即与球心距离为a/2的位置，电导率σ＝σmax，中间均匀过度；σmax表示电导率σ可取的最大值，按照计算得到，其中，εr为相对介电常数，π为圆周率；标记为inner-free的自由空间中，σ是某一实数，具体数值根据内部材料的不同而不同；得到了球形的自由空间和球形各向异性完全匹配层；
3)对球形内部区域的电磁场系数进行初始化；
用(i,j,k)表示位置坐标，位置坐标可取的范围为步骤1)中所设定的区域范围，计算时，(i,j,k)从(Xn,Yn,Zn)点逐点循环至(Xp,Yp,Zp)点，在计算区域内对电场强度分量、电场强度中间变量、上一时刻电场强度中间变量、磁场强度分量、磁场强度中间变量、上一时刻磁场强度中间变量进行初始化操作，即均置为0，再由时域有限差分方法分别计算电场系数和磁场系数；
4)更新球形内部区域的磁场强度中间变量；
根据步骤3)中所设定的磁场系数，在标记为inner的球形内部区域，给每一点的上一时刻磁场强度中间变量赋值，数值等于上一时刻(即n时刻)的磁场强度中间变量，赋值完成后，计算下一时刻(即n+1时刻)的磁场强度中间变量；
5)更新球形内部区域的磁场强度分量；
根据步骤4)中的磁场强度中间变量在n时刻和n+1时刻的数值和步骤3)中所设定的磁场系数，在标记为inner的球形内部区域，计算出每一点在n+1时刻的磁场强度；
6)更新球形内部区域的电场强度中间变量；
根据步骤3)中所设定的电场系数，在标记为inner的球形内部区域，给每一点的上一时刻电场强度中间变量赋值，数值等于上一时刻(即n－1/2时刻)的电场强度中间变量，赋值完成后，计算下一时刻(即n+1/2时刻)的电场强度中间变量；
7)更新球形内部区域的电场强度分量；
根据步骤6)中的电场强度中间变量在n－1/2时刻和n+1/2时刻的数值和步骤3)中所设定的电场系数，在标记为inner的球形内部区域，计算出每一点在n+1/2时刻的电场强度；
8)更新电场的激励源；
按照公式逐步更新激励源的数值，其中，n为时间步；J0为脉冲的幅值；τ为常数；脉冲峰值出现在n＝n0时刻；将更新后的激励源数值赋值给球形区域内中心位置的电场强度分量Ez，即完成电场激励源位置及数值的更新；
9)每次循环时间步n都会加一，判断迭代次数n是否达到步骤1)所设定的电磁仿真迭代步数N，以判断是否满足更新条件；若未达到迭代步数，则在n加一后返回步骤4)继续循环；若达到迭代步数，即n＝N，则记录步骤5)中得到的磁场强度和步骤7)中得到的电场强度，并保存作为最终所得结果。
与现有技术相比，本发明的有益效果是：
本发明方法是建立在笛卡尔坐标系下的，三维问题中各向异性完全匹配层的外形为球形，二维问题中各向异性完全匹配层的外形为圆形，这样可以保证网格大小的一致，不会有不稳定的因素出现，又能截断一些多余的网格，减少计算量，提高计算效率，克服了现有笛卡尔坐标系下的方形各向异性完全匹配层计算效率低、计算复杂的问题，将球形的外形和方形的网格相结合，既稳定效率又高。本发明各向异性完全匹配层截断边界具有表达简单、更稳定、应用更广泛等优势。
本发明一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，在二维情况下，较传统的正方形各向异性完全匹配层可节省21.545％以上的内存，利用本申请所述的参数和方法设置边界参数，可实现对外向传播的电磁波的完全吸收，应用实施例2中，图3中电磁波的传播呈现出同心圆的形状，说明边界处是无反射的，所述的一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法可以将电磁波完全吸收，图4中将点源放在了正方形区域的非正中心，说明了二维情况下圆形边界的有效性，较传统的正方形各向异性完全匹配层，本发明的计算效率提高了1.13倍；在三维情况下，较传统的立方体形各向异性完全匹配层可节省47.877％以上的内存，利用本申请中所述的参数和方法设置边界参数，可实现对外向传播的电磁波的完全吸收，应用实施例3中，图5中的对比结果的一致表明了所发明方法的稳定性，图6表示对于实施例3通过所述的一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法计算的结果和解析解的误差，表明了这种方法的稳定性和有效性。较传统的立方体形各向异性完全匹配层，本发明的计算效率提高了1.25倍。
附图说明
图1是三维问题中球形截断边界与相对应的立方体截断边界的对比图；
图2是在时域有限差分方法中本发明建立在二维笛卡尔坐标系下的网格分布，也表示三维情况时沿任意直径的截面的网格分布，中间白色网格表示计算区域的自由空间，周围灰色网格表示笛卡尔坐标系下的各向异性完全匹配层截断边界，外侧白色网格表示与传统正方形截断边界相比本发明可以忽略不算的网格；
图3是应用实施例2在二维圆形各向异性完全匹配层截断边界问题中，当点源置于计算区域的中心位置时的电磁波传播图；
图4是应用实施例2在二维圆形各向异性完全匹配层截断边界问题中，当点源置于计算区域的偏心坐标位置为(70,70)时的电磁波传播图；
图5是应用实施例3在三维球形各向异性完全匹配层截断边界问题中，坐标点(30,30,30)处的电场强度Ez的时域分布与解析解对比的结果；
图6是应用实施例3在三维球形各向异性完全匹配层截断边界问题中，所得的结果和解析解之间的误差。
具体实施方式
本发明用于计算电磁学的时域有限差分法中，所涉及的递推公式可根据时域有限差分方法得到，具体实施方式以三维问题为例说明具体步骤，二维问题和三维问题的计算步骤是完全一样的，只是变量的维度所有不同。
本发明笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，该方法的步骤为：
1)建立求解对象的模型数据和时域有限差分法的计算空间；
向计算机申请内存空间，设Xn,Yn,Zn分别为x,y,z方向上的最小位置，其中Xn＝Yn＝Zn，Xp,Yp,Zp分别为x,y,z方向上的最大位置，其中Xp＝Yp＝Zp，总体计算的区域范围为(Xn,Yn,Zn)→(Xp,Yp,Zp)，呈立方体形，区域范围的大小为(Xp–Xn)×(Yp–Yn)×(Zp–Zn)，具体的数值是与计算目标的大小相关的，可以根据问题的需要放大或缩小，设x,y,z方向上的空间步长为Δx,Δy和Δz，其中Δx＝Δy＝Δz，构成立方体网格区域；设球形各向异性完全匹配层的直径为a，所对应的网格数量为a/Δx＝Xp–Xn，中心坐标为((Xp–Xn)/2,(Yp–Yn)/2,(Zp–Zn)/2)，用球形的边界截断立方体网格区域，使计算区域成为球形，设定时间步长为Δt，并对球形内部区域的网格做inner标记，对球形外部区域的网格做outer标记，以区分不同的区域；设定电磁仿真迭代步数为N，n表示时间步，n的范围为1→N，选择时谐场的点源和电流源分别作为二维和三维问题的激励源，激励源函数用Einc表示；
2)构建球形各向异性完全匹配层截断边界；
在步骤1)所设置的标记为inner的球形内部区域，设定δ为球形各向异性完全匹配层截断边界的层数，将标记为inner的球形内部区域分为两部分，包括半径从0到a/2－δ*Δx的区域，称为自由空间，做inner-free标记，半径从a/2－δ*Δx到a/2的区域，即为球形各向异性完全匹配层截断边界，做inner-CSUPML标记；在球形各向异性完全匹配层截断边界的最内层，即与球心距离为a/2－δ*Δx的位置，电导率σ＝0，在球形各向异性完全匹配层截断边界的最外层，即与球心距离为a/2的位置，电导率σ＝σmax，中间均匀过度；σmax表示电导率σ可取的最大值，按照计算得到，其中，εr为相对介电常数，π为圆周率；标记为inner-free的自由空间中，σ是某一实数，具体数值根据内部材料的不同而不同；这样，电导率σ的计算区域外形也为球形，得到了球形的自由空间和球形各向异性完全匹配层；通过上述设计笛卡尔坐标系下的球形各向异性完全匹配层截断边界的参数，用于计算电磁学的程序中，可实现在有限的计算区域达到模拟“微波暗室”吸波材料的效果；
3)对球形内部区域的电磁场系数进行初始化；
用(i,j,k)表示位置坐标，位置坐标可取的范围为步骤1)中所设定的区域范围，计算时，(i,j,k)从(Xn,Yn,Zn)点逐点循环至(Xp,Yp,Zp)点，在计算区域内对电场强度分量、电场强度中间变量、上一时刻电场强度中间变量、磁场强度分量、磁场强度中间变量、上一时刻磁场强度中间变量进行初始化操作，即均置为0，再由时域有限差分方法分别计算电场系数和磁场系数，
电场强度分量Ex|i,j,k＝0，Ey|i,j,k＝0，Ez|i,j,k＝0；
电场强度中间变量Dx|i,j,k＝0，Dy|i,j,k＝0，Dz|i,j,k＝0；
上一时刻电场强度中间变量Dstorex|i,j,k＝0，Dstorey|i,j,k＝0，Dstorez|i,j,k＝0；
磁场强度分量Hx|i,j,k＝0，Hy|i,j,k＝0，Hz|i,j,k＝0；
磁场强度中间变量Bx|i,j,k＝0，By|i,j,k＝0，Bz|i,j,k＝0；
上一时刻磁场强度中间变量Bstorex|i,j,k＝0，Bstorey|i,j,k＝0，Bstorez|i,j,k＝0；
磁场系数由时域有限差分方法可得：
电场系数由时域有限差分方法可得：
其中，κx、κy、κz表示各向异性介质匹配矩阵中参数的实部，根据这种方法设置电场系数和磁场系数，由于其中包含了步骤2)中所设置的参数σ，所以电场系数和磁场系数的外形也为球形。
4)更新球形内部区域的磁场强度中间变量；
根据步骤3)中所设定的磁场系数，在标记为inner的球形内部区域，给每一点的上一时刻磁场强度中间变量赋值，数值等于上一时刻(即n时刻)的磁场强度中间变量，赋值完成后，计算下一时刻(即n+1时刻)的磁场强度中间变量；
5)更新球形内部区域的磁场强度分量；
根据步骤4)中的磁场强度中间变量在n时刻和n+1时刻的数值和步骤3)中所设定的磁场系数，在标记为inner的球形内部区域，计算出每一点在n+1时刻的磁场强度；
6)更新球形内部区域的电场强度中间变量；
根据步骤3)中所设定的电场系数，在标记为inner的球形内部区域，给每一点的上一时刻电场强度中间变量赋值，数值等于上一时刻(即n－1/2时刻)的电场强度中间变量，赋值完成后，计算下一时刻(即n+1/2时刻)的电场强度中间变量；
7)更新球形内部区域的电场强度分量；
根据步骤6)中的电场强度中间变量在n－1/2时刻和n+1/2时刻的数值和步骤3)中所设定的电场系数，在标记为inner的球形内部区域，计算出每一点在n+1/2时刻的电场强度；
8)更新电场的激励源；
选择时谐场的电流源作为激励源，按照公式逐步更新激励源的数值，其中，n为时间步；J0为脉冲的幅值；τ为常数，决定了脉冲的宽度；脉冲峰值出现在n＝n0时刻；将更新后的激励源数值赋值给球形区域内中心位置的电场强度分量Ez，即完成电场激励源位置及数值的更新；
9)每次循环时间步n都会加一，判断迭代次数n是否达到步骤1)所设定的电磁仿真迭代步数N，以判断是否满足更新条件；若未达到迭代步数，则在n加一后返回步骤4)继续循环；若达到迭代步数，即n＝N，则记录步骤5)中得到的磁场强度和步骤7)中得到的电场强度，并保存作为最终所得结果。
本发明的进一步特征在于所述电磁仿真迭代步数N取1000以上。
本发明方法将各向异性完全匹配层与笛卡尔坐标系下的球形外形结合，具体的设计原理是：
1.各向异性完全匹配层截断边界的形状；
所述的一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法在三维问题中形状为球形。与传统的三维立方体外形的截断边界相比，省去了十二个棱区域和八个角区域的计算，原理示意图如图1所示。设笛卡尔坐标系下每个立方体网格的边长为Δx，对于传统的立方体形截断边界，设其边长为a，则立方体形截断边界所包围的计算区域内总的网格数量为N1’＝(a/Δx)3；对于所述三维球形各向异性完全匹配层截断边界，设其直径与传统的立方体形截断边界的边长相等，也为a，则球形各向异性完全匹配层截断边界所包围的计算区域内总的网格数量为N1≈0.52122×(a/Δx)3－5.63812×(a/Δx)2－34.11802×(a/Δx)+1521.017，因此本发明节省的体积百分比为(N1’－N1)/N1’≥47.877％，这个比例会随着所设直径的不同而变化，但最小值是47.877％；对于无限长的物理目标，只考虑垂直于长度截面上的物理变化过程，可以投影为二维问题。本发明笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法在二维问题中形状为圆形，与传统二维正方形的截断边界相比，省去了四个角区域的计算，原理示意图如图2所示。设笛卡尔坐标系下每个正方形网格的边长为Δy，对于传统的正方形截断边界，设其边长为b，则正方形截断边界所包围的计算区域内总的网格数量为N2’＝(b/Δy)2；对于所述的二维圆形各向异性完全匹配层截断边界，设其直径与传统的正方形截断边界的边长相等，也为b，则圆形各向异性完全匹配层截断边界所包围的计算区域内总的网格数量为N2≈0.78454×(b/Δy)2－6.24616×(b/Δy)+0.20964，因此本发明节省的体积百分比为(N2’－N2)/N2’≥21.545％，这个比例会随着所设直径的不同而变化，但最小值是21.545％。
2.各向异性完全匹配层截断边界的参数设计；
在本发明笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法中，计算参数表示电导率σ可取的最大值，εr为其相对介电常数，π为圆周率，δ为其层数。设置参数时，截断边界的最内层σ＝0，截断边界的最外层σ＝σmax，中间均匀过度，这种设置方法适用于二维和三维情况的所有问题。通过上述表达式设计笛卡尔坐标系下的球形各向异性完全匹配层截断边界的参数，用于计算电磁学的程序中，可实现在有限的计算区域达到模拟“微波暗室”吸波材料的效果。
3.将这种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法应用于计算电磁学中的时域有限差分方法；
在本发明笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法中，对于三维问题，使用笛卡尔坐标系下固有的立方体的网格进行剖分，对计算区域用球形边界进行截断处理，选取5～10层网格作为边界上各向异性完全匹配层的厚度，沿任意直径的截面如图2所示，利用公式得到的σmax作为球形各向异性完全匹配层截断边界的参数，从而模拟电磁波在无界空间的传播。对于二维问题，则采用正方形网格进行剖分，利用圆形边界来截断计算区域，网格分布如图2所示，利用σmax作为圆形各向异性完全匹配层截断边界的参数，这样就可以实现二维电磁问题的计算。
在计算目标散射问题时，散射体为计算目标，本发明中计算区域范围的大小与与计算目标的大小相关，可以根据问题的需要放大或缩小计算目标，得到区域范围。
实施例1
本实施例笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，应用在三维问题中，各向异性完全匹配层截断边界的外形为球形，此时该方法的步骤为：
1.建立求解对象的模型数据和时域有限差分法的计算空间；
向计算机申请内存空间，总体计算的区域范围为(Xn,Yn,Zn)→(Xp,Yp,Zp)，其大小为(Xp–Xn)×(Yp–Yn)×(Zp–Zn)，其中Xn＝Yn＝Zn，Xp＝Yp＝Zp，x,y,z方向上的空间步长分别为Δx,Δy和Δz，且Δx＝Δy＝Δz，在截断边界以内的计算区域为真空状态，选择时谐场的电流源为三维问题的激励源，定义电场的激励源函数Einc；并对球形内部区域的网格做inner标记，对球形外部区域的网格做outer标记，以区分不同的区域；时间步长为Δt，设定电磁仿真迭代步数为1000步。
2.构建球形各向异性完全匹配层截断边界；
在步骤1所设置的标记为inner的球形内部区域，设定δ为球形各向异性完全匹配层截断边界的层数，这样又可以将这个区域分为两部分，包括半径从0到a/2－δ*Δx的区域，称为自由空间，做inner-free标记，半径从a/2－δ*Δx到a/2的区域，即为球形各向异性完全匹配层截断边界，做inner-CSUPML标记；在球形各向异性完全匹配层截断边界的最内层，即与球心距离为a/2－δ*Δx的位置，电导率σ＝0，球形各向异性完全匹配层截断边界的最外层，即与球心距离为a/2的位置，电导率σ＝σmax，中间均匀过度；σmax表示电导率σ可取的最大值，按照计算得到，其中，εr为相对介电常数，π为圆周率；标记为inner-free的自由空间中，σ是某一实数，具体数值根据内部材料的不同而不同；这样，电导率σ的计算区域外形也为球形，得到了球形的自由空间和球形各向异性完全匹配层；
3.对球形内部区域的电磁场系数进行初始化；
对于所述的一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法三维问题下的电场和磁场变量包括：电场强度分量Ex(Xn:Xp–1,Yn:Yp,Zn:Zp),Ey(Xn:Xp,Yn:Yp–1,Zn:Zp),Ez(Xn:Xp,Yn:Yp,Zn:Zp–1)，电场强度中间变量Dx(Xn:Xp–1,Yn:Yp,Zn:Zp),Dy(Xn:Xp,Yn:Yp–1,Zn:Zp),Dz(Xn:Xp,Yn:Yp,Zn:Zp–1)，和上一时刻电场强度中间变量Dstorex(Xn:Xp,Yn:Yp,Zn:Zp)，Dstorey(Xn:Xp,Yn:Yp,Zn:Zp)，Dstorez(Xn:Xp,Yn:Yp,Zn:Zp)，磁场强度分量Hx(Xn:Xp,Yn:Yp–1,Zn:Zp–1),Hy(Xn:Xp–1,Yn:Yp,Zn:Zp–1),Hz(Xn:Xp–1,Yn:Yp–1,Zn:Zp)，磁场强度中间变量Bx(Xn:Xp,Yn:Yp–1,Zn:Zp–1),By(Xn:Xp–1,Yn:Yp,Zn:Zp–1),Bz(Xn:Xp–1,Yn:Yp–1,Zn:Zp)，和上一时刻磁场强度中间变量Bstorex(Xn:Xp,Yn:Yp,Zn:Zp)，Bstorey(Xn:Xp,Yn:Yp,Zn:Zp),Bstorez(Xn:Xp,Yn:Yp,Zn:Zp)，电场系数D1(Xn:Xp,Yn:Yp,Zn:Zp)，D2(Xn:Xp,Yn:Yp,Zn:Zp)，D3(Xn:Xp,Yn:Yp,Zn:Zp)，D4(Xn:Xp,Yn:Yp,Zn:Zp)，D5(Xn:Xp,Yn:Yp,Zn:Zp)，D6(Xn:Xp,Yn:Yp,Zn:Zp)，磁场系数C1(Xn:Xp,Yn:Yp,Zn:Zp)，C2(Xn:Xp,Yn:Yp,Zn:Zp)，C3(Xn:Xp,Yn:Yp,Zn:Zp)，C4(Xn:Xp,Yn:Yp,Zn:Zp)，C5(Xn:Xp,Yn:Yp,Zn:Zp)，C6(Xn:Xp,Yn:Yp,Zn:Zp)，其中，Xn:Xp表示范围从Xn到Xp，其他以此类推，在后面的步骤中，用(i,j,k)表示上述各变量的位置坐标，(i,j,k)从(Xn,Yn,Zn)点逐点循环至(Xp,Yp,Zp)点；
4.更新球形内部区域的磁场强度中间变量；
在标记为inner的球形内部区域，给每一点的上一时刻磁场强度中间变量赋值，数值等于上一时刻(即n时刻)的磁场强度中间变量，赋值完成后，计算下一时刻(即n+1时刻)的磁场强度中间变量；
Bstorex(i,j,k)＝Bx(i,j,k)
Bstorey(i,j,k)＝By(i,j,k)
Bstorez(i,j,k)＝Bz(i,j,k)
Bx(i,j,k)＝D1(i,j,k)×Bx(i,j,k)－D2(i,j,k)×((Ez(i,j+1,k)－Ez(i,j,k))/Δy－(Ey(i,j,k+1)－Ey(i,j,k))/Δz)
By(i,j,k)＝D1(i,j,k)×By(i,j,k)－D2(i,j,k)×((Ex(i,j,k+1)－Ex(i,j,k))/Δz－(Ez(i+1,j,k)－Ez(i,j,k))/Δx)
Bz(i,j,k)＝D1(i,j,k)×Bz(i,j,k)－D2(i,j,k)×((Ey(i+1,j,k)－Ey(i,j,k))/Δx－(Ex(i,j+1,k)－Ex(i,j,k))/Δy)
5.更新球形内部区域的磁场强度分量；
在标记为inner的球形内部区域，计算出每一点在n+1时刻的磁场强度；
Hx(i,j,k)＝D3(i,j,k)×Hx(i,j,k)+D4(i,j,k)×(D5(i,j,k)×Bx(i,j,k)－D6(i,j,k)×Bstorex(i,j,k))
Hy(i,j,k)＝D3(i,j,k)×Hy(i,j,k)+D4(i,j,k)×(D5(i,j,k)×By(i,j,k)－D6(i,j,k)×Bstorey(i,j,k))
Hz(i,j,k)＝D3(i,j,k)×Hz(i,j,k)+D4(i,j,k)×(D5(i,j,k)×Bz(i,j,k)－D6(i,j,k)×Bstorez(i,j,k))
6.更新球形内部区域的电场强度中间变量；
在标记为inner的球形内部区域，给每一点的上一时刻电场强度中间变量赋值，数值等于上一时刻(即n－1/2时刻)的电场强度中间变量，赋值完成后，计算下一时刻(即n+1/2时刻)的电场强度中间变量；
Dstorex(i,j,k)＝Dx(i,j,k)
Dstorey(i,j,k)＝Dy(i,j,k)
Dstorez(i,j,k)＝Dz(i,j,k)
Dx(i,j,k)＝C1(i,j,k)×Dx(i,j,k)+C2(i,j,k)×((Hz(i,j,k)－Hz(i,j－1,k))/Δy－(Hy(i,j,k)－Hy(i,j,k－1))/Δz)
Dy(i,j,k)＝C1(i,j,k)×Dy(i,j,k)+C2(i,j,k)×((Hx(i,j,k)－Hx(i,j,k－1))/Δz－(Hz(i,j,k)－Hz(i－1,j,k))/Δx)
Dz(i,j,k)＝C1(i,j,k)×Dz(i,j,k)+C2(i,j,k)×((Hy(i,j,k)－Hy(i－1,j,k))/Δx－(Hx(i,j,k)－Hx(i,j－1,k))/Δy)
7.更新球形内部区域的电场强度分量；
在标记为inner的球形内部区域，计算出每一点在n+1/2时刻的电场强度；
Ex(i,j,k)＝C3(i,j,k)×Ey(i,j,k)+C4(i,j,k)×(C5(i,j,k)×Dx(i,j,k)－C6(i,j,k)×Dstorex(i,j,k))
Ey(i,j,k)＝C3(i,j,k)×Ey(i,j,k)+C4(i,j,k)×(C5(i,j,k)×Dy(i,j,k)－C6(i,j,k)×Dstorey(i,j,k))
Ez(i,j,k)＝C3(i,j,k)×Ez(i,j,k)+C4(i,j,k)×(C5(i,j,k)×Dz(i,j,k)－C6(i,j,k)×Dstorez(i,j,k))
8.更新电场的激励源；
按照公式逐步更新激励源的数值，脉冲峰值出现在n＝n0时刻；将更新后的激励源数值赋值给球形区域内中心位置的电场强度分量Ez，即Ez((Xp–Xn)/2,(Yp–Yn)/2,(Zp–Zn)/2)＝Einc，完成电场激励源位置及数值的更新；
9.每次循环时间步n都会加一，判断迭代次数n是否达到步骤1)所设定的电磁仿真迭代步数N，以判断是否满足更新条件；若未达到迭代步数，则在n加一后返回步骤4)继续循环；若达到迭代步数，n＝N，则记录步骤5)中得到的磁场强度和步骤7)中得到的电场强度，并保存作为最终所得结果。
实施例2
本实施例实现方法各步骤同实施例1，不同之处在于本实施例应用于二维情况中，Δz为0，各向异性完全匹配层截断边界为圆形，验证二维情况下各向异性完全匹配层截断边界的实现方法的吸收效果。计算区域范围的大小设定为110×110，空间步长为Δx＝Δy＝1mm，时间步长为Δt＝16.667ps，整个计算区域为真空状态，其电导率为σ＝0，磁导率为μ0，介电常数为ε0。选择时谐场的正弦点源做为激励源，表达式为Einc＝sin(2πf0NΔt)，f0是源的频率，电磁仿真迭代步数为N＝1000。由此运行程序，结果如图3所示，当点源置于计算区域的中心位置时，能得到明显的同心圆图。图4是把点源置于计算区域的偏心坐标位置为(70,70)时的计算结果，结果仍然近似为同心圆，存在微弱的反射波返回至中心区域，所以呈现出了轻微的形变，但这并不影响计算结果。图中能看到截断边界是圆形。在程序运行过程中，验证了在二维笛卡尔坐标系下利用圆形各向异性完全匹配层截断边界的计算效率较方形边界提高了1.13倍。
实施例3
本实施例实现方法应用于三维情况中，具体实现方法的步骤同实施例1，验证三维情况下所述的一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法的吸收效果。计算区域范围的大小设定为62×62×62，空间步长为Δx＝Δy＝Δz＝2mm，时间步长为Δt＝3.333ps，计算区域的电导率为σ＝0，磁导率为μ0，介电常数为ε0。在坐标点(30,30,30)处加入电偶极子，选择电流源作为激励源，表达式为Einc＝J0×(t－t0)×exp(－(t－t0)2/τ2)，其电磁仿真迭代步数为N＝1000。由此运行程序，记录坐标点(30,30,30)处的电场强度Ez的时域分布，并对比解析解的结果如图5所示，可以看出二者高度吻合。图6表示为本实施例计算的结果和解析解的误差，表明了这种方法的稳定性和有效性。在程序运行过程中，验证了在三维笛卡尔坐标系下利用球形各向异性完全匹配层截断边界的计算效率较立方体形边界提高了1.25倍。
本发明未述及之处适用于现有技术。
Claims (3)
1.一种笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，该方法的步骤为：
1)建立求解对象的模型数据和时域有限差分法的计算空间；
向计算机申请内存空间，设Xn,Yn,Zn分别为x,y,z方向上的最小位置，其中Xn＝Yn＝Zn，Xp,Yp,Zp分别为x,y,z方向上的最大位置，其中Xp＝Yp＝Zp，总体计算的区域范围为(Xn,Yn,Zn)→(Xp,Yp,Zp)，呈立方体形，区域范围的大小为(Xp–Xn)×(Yp–Yn)×(Zp–Zn)，设x,y,z方向上的空间步长为Δx,Δy和Δz，且Δx＝Δy＝Δz，构成立方体网格区域；设球形各向异性完全匹配层的直径为a，所对应的网格数量为a/Δx＝Xp–Xn，中心坐标为((Xp–Xn)/2,(Yp–Yn)/2,(Zp–Zn)/2)，用球形的边界截断立方体网格区域，使计算区域成为球形，设定时间步长为Δt，并对球形内部区域的网格做inner标记，对球形外部区域的网格做outer标记，以区分不同的区域；设定电磁仿真迭代步数为N，n表示时间步，n的范围为1→N，选择时谐场的点源和电流源分别作为二维和三维问题的激励源，激励源函数用Einc表示；
2)构建球形各向异性完全匹配层截断边界；
在步骤1)所设置的标记为inner的球形内部区域，设定δ为球形各向异性完全匹配层截断边界的层数，将标记为inner的球形内部区域分为两部分，包括半径从0到a/2－δ*Δx的区域，称为自由空间，做inner-free标记，半径从a/2－δ*Δx到a/2的区域，即为球形各向异性完全匹配层截断边界，做inner-CSUPML标记；在球形各向异性完全匹配层截断边界的最内层，即与球心距离为a/2－δ*Δx的位置，电导率σ＝0，在球形各向异性完全匹配层截断边界的最外层，即与球心距离为a/2的位置，电导率σ＝σmax，中间均匀过度；σmax表示电导率σ可取的最大值，按照计算得到，其中，εr为相对介电常数，π为圆周率；标记为inner-free的自由空间中，σ是某一实数，具体数值根据内部材料的不同而不同；得到了球形的自由空间和球形各向异性完全匹配层；
3)对球形内部区域的电磁场系数进行初始化；
用(i,j,k)表示位置坐标，位置坐标可取的范围为步骤1)中所设定的区域范围，计算时，(i,j,k)从(Xn,Yn,Zn)点逐点循环至(Xp,Yp,Zp)点，在计算区域内对电场强度分量、电场强度中间变量、上一时刻电场强度中间变量、磁场强度分量、磁场强度中间变量、上一时刻磁场强度中间变量进行初始化操作，即均置为0，再由时域有限差分方法分别计算电场系数和磁场系数；
4)更新球形内部区域的磁场强度中间变量；
根据步骤3)中所设定的磁场系数，在标记为inner的球形内部区域，给每一点的上一时刻磁场强度中间变量赋值，数值等于上一时刻(即n时刻)的磁场强度中间变量，赋值完成后，计算下一时刻(即n+1时刻)的磁场强度中间变量；
5)更新球形内部区域的磁场强度分量；
根据步骤4)中的磁场强度中间变量在n时刻和n+1时刻的数值和步骤3)中所设定的磁场系数，在标记为inner的球形内部区域，计算出每一点在n+1时刻的磁场强度；
6)更新球形内部区域的电场强度中间变量；
根据步骤3)中所设定的电场系数，在标记为inner的球形内部区域，给每一点的上一时刻电场强度中间变量赋值，数值等于上一时刻(即n－1/2时刻)的电场强度中间变量，赋值完成后，计算下一时刻(即n+1/2时刻)的电场强度中间变量；
7)更新球形内部区域的电场强度分量；
根据步骤6)中的电场强度中间变量在n－1/2时刻和n+1/2时刻的数值和步骤3)中所设定的电场系数，在标记为inner的球形内部区域，计算出每一点在n+1/2时刻的电场强度；
8)更新电场的激励源；
按照公式逐步更新激励源的数值，其中，n为时间步；J0为脉冲的幅值；τ为常数；脉冲峰值出现在n＝n0时刻；将更新后的激励源数值赋值给球形区域内中心位置的电场强度分量Ez，即完成电场激励源位置及数值的更新；
9)每次循环时间步n都会加一，判断迭代次数n是否达到步骤1)所设定的电磁仿真迭代步数N，以判断是否满足更新条件；若未达到迭代步数，则在n加一后返回步骤4)继续循环；若达到迭代步数，即n＝N，则记录步骤5)中得到的磁场强度和步骤7)中得到的电场强度，并保存作为最终所得结果。
2.根据权利要求1所述的笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，其特征在于，所述电磁仿真迭代步数N取1000以上。
3.根据权利要求1所述的笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法，其特征在于，该方法应用于二维问题中，各向异性完全匹配层截断边界的外形为圆形；该方法应用于三维问题中，各向异性完全匹配层截断边界的外形为球形。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711201147.7A CN107944214B (zh) | 2017-11-27 | 2017-11-27 | 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711201147.7A CN107944214B (zh) | 2017-11-27 | 2017-11-27 | 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107944214A true CN107944214A (zh) | 2018-04-20 |
| CN107944214B CN107944214B (zh) | 2020-11-10 |
Family
ID=61948894
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711201147.7A Active CN107944214B (zh) | 2017-11-27 | 2017-11-27 | 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107944214B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109616774A (zh) * | 2018-11-20 | 2019-04-12 | 浙江大学 | 一种基于空间分布pml模型的吸收材料和微型暗室 |
| CN110162896A (zh) * | 2019-05-27 | 2019-08-23 | 河北工业大学 | 一种复频移卷积实现完全匹配层的方法 |
| CN113917526A (zh) * | 2020-07-10 | 2022-01-11 | 中国石油化工股份有限公司 | 基于非分裂完全匹配层吸收边界的正演模拟方法 |
| CN114722653A (zh) * | 2022-03-09 | 2022-07-08 | 大连理工大学 | 基于立方体单元有限差分的球形金属粒子传热数值计算方法 |
| CN115640711A (zh) * | 2021-07-17 | 2023-01-24 | 深圳市芯瑞微电子有限公司 | 一种自适应方向的网格划分方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2004239736A (ja) * | 2003-02-05 | 2004-08-26 | Ntt Docomo Inc | 電磁界特性解析方法 |
| CN104750990A (zh) * | 2015-03-30 | 2015-07-01 | 西安理工大学 | 二维等离子体中扩展坐标的完全匹配吸收边界的实现方法 |
| CN104794289A (zh) * | 2015-04-23 | 2015-07-22 | 西安理工大学 | 一种扩展直角坐标系下完全匹配吸收边界的实现方法 |
| CN105589980A (zh) * | 2014-10-23 | 2016-05-18 | 天津职业技术师范大学 | 一种阻抗匹配层的截断边界 |
| CN105808968A (zh) * | 2016-04-13 | 2016-07-27 | 吉林大学 | 时域航空电磁数值模拟中c-pml边界条件加载方法 |
-
2017
- 2017-11-27 CN CN201711201147.7A patent/CN107944214B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2004239736A (ja) * | 2003-02-05 | 2004-08-26 | Ntt Docomo Inc | 電磁界特性解析方法 |
| CN105589980A (zh) * | 2014-10-23 | 2016-05-18 | 天津职业技术师范大学 | 一种阻抗匹配层的截断边界 |
| CN104750990A (zh) * | 2015-03-30 | 2015-07-01 | 西安理工大学 | 二维等离子体中扩展坐标的完全匹配吸收边界的实现方法 |
| CN104794289A (zh) * | 2015-04-23 | 2015-07-22 | 西安理工大学 | 一种扩展直角坐标系下完全匹配吸收边界的实现方法 |
| CN105808968A (zh) * | 2016-04-13 | 2016-07-27 | 吉林大学 | 时域航空电磁数值模拟中c-pml边界条件加载方法 |
Non-Patent Citations (2)
| Title |
|---|
| LU WANG 等: "A Circular-Shaped Perfectly Matched Layer Strategy for Rectangular FDTD Grids", 《2017 SIXTH ASIA-PACIFIC CONFERENCE ON ANTENNAS AND PROPAGATION》 * |
| 王向华 等: "各向异性完全匹配层边界条件下雷达散射截面的误差修正", 《天津职业技术师范大学学报》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109616774A (zh) * | 2018-11-20 | 2019-04-12 | 浙江大学 | 一种基于空间分布pml模型的吸收材料和微型暗室 |
| CN110162896A (zh) * | 2019-05-27 | 2019-08-23 | 河北工业大学 | 一种复频移卷积实现完全匹配层的方法 |
| CN110162896B (zh) * | 2019-05-27 | 2023-04-18 | 河北工业大学 | 一种复频移卷积实现完全匹配层的方法 |
| CN113917526A (zh) * | 2020-07-10 | 2022-01-11 | 中国石油化工股份有限公司 | 基于非分裂完全匹配层吸收边界的正演模拟方法 |
| CN115640711A (zh) * | 2021-07-17 | 2023-01-24 | 深圳市芯瑞微电子有限公司 | 一种自适应方向的网格划分方法 |
| CN115640711B (zh) * | 2021-07-17 | 2024-04-19 | 芯瑞微(上海)电子科技有限公司 | 一种自适应方向的网格划分方法 |
| CN114722653A (zh) * | 2022-03-09 | 2022-07-08 | 大连理工大学 | 基于立方体单元有限差分的球形金属粒子传热数值计算方法 |
| CN114722653B (zh) * | 2022-03-09 | 2024-09-13 | 大连理工大学 | 基于立方体单元有限差分的球形金属粒子传热数值计算方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107944214B (zh) | 2020-11-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107944214A (zh) | 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 | |
| Makarov et al. | Low-frequency electromagnetic modeling for electrical and biological systems using MATLAB | |
| CN103218487B (zh) | 旋转对称天线罩和抛物面天线一体化电磁散射仿真方法 | |
| CN107526887B (zh) | 一种基于GPU并行的LeapfrogADI-FDTD方法 | |
| CN102081690B (zh) | 复杂电路的矩阵分解结合新奇异值分解方法 | |
| CN102156764A (zh) | 一种分析天线辐射和电磁散射的多分辨预条件方法 | |
| CN106446470B (zh) | 一种高效并行的非均匀介质频域有限差分方法 | |
| Zhang et al. | The fast physical optics method on calculating the scattered fields from electrically large scatterers | |
| Ozgun et al. | Software metamaterials: Transformation media based multiscale techniques for computational electromagnetics | |
| Xie et al. | Tailoring unstructured meshes for use with a 3D time domain co‐volume algorithm for computational electromagnetics | |
| CN105589980A (zh) | 一种阻抗匹配层的截断边界 | |
| CN110162896B (zh) | 一种复频移卷积实现完全匹配层的方法 | |
| CN115984440B (zh) | 对象渲染方法、装置、计算机设备和存储介质 | |
| Hasanovic et al. | Scattering from 3-D inhomogeneous chiral bodies of arbitrary shape by the method of moments | |
| Potapov | Fractal electrodynamics: Numerical modeling of small fractal antenna devices and fractal 3D microwave resonators for modern ultra-wideband or multiband radio systems | |
| CN104778286A (zh) | 掠海飞行器电磁散射特性快速仿真方法 | |
| Chen et al. | Analysis of antenna around NURBS surface with iterative MoM-PO technique | |
| Saeedfar et al. | Shape reconstruction of three-dimensional conducting curved plates using physical optics, NURBS modeling, and genetic algorithm | |
| CN106503349B (zh) | 一种类周期结构目标电磁散射特性快速计算方法 | |
| Zhao et al. | Rigid blocks matching method based on contour curves and feature regions | |
| Qi et al. | Acceleration strategies based on an improved bubble packing method | |
| Jiang et al. | Modified adaptive cross approximation algorithm for analysis of electromagnetic problems | |
| Fu et al. | Creeping ray tracing algorithm for arbitrary NURBS surfaces based on adaptive variable step Euler method | |
| Hollander et al. | Adaptive multilevel nonuniform grid algorithm for the accelerated analysis of composite metallic–dielectric radomes | |
| Dhibi et al. | Multi‐mother wavelet neural network‐based on genetic algorithm and multiresolution analysis for fast 3D mesh deformation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |