An improved 2D/3D hybrid discontinuous Galerkin time domain method
-
摘要: 对高速信号通过电源板时的电源完整性(power integrity, PI)问题进行研究时, 因为电源板中主要模式分布为零阶平行板模式, 可以采用二维简化以提高效率.而对于隔离盘或其它存在纵向不连续性的区域, 则应采用三维算法以保证精度.将两者结合起来的一种二维三维(2D/3D)混合时域不连续伽辽金(discontinuous Galerkin time domain, DGTD)方法可以兼顾精度与效率, 有效地处理这类电磁全波计算问题.其中二维、三维方法采用同一套三棱柱离散的网格, 通过适当设置基函数, 二维区域与二维区域之间可以方便快速地相互转化.随着电磁波的传播, 二维、三维的适用区域是随时间、空间动态变化的.为了准确地捕捉这种动态变化, 文中提出的一种改进的自适应判据, 在每个时间歩对电磁场进行检测, 从而动态地判定二维简化区域.与现有技术的判据控制绝对误差不同, 该方法对相对误差进行控制, 效率高、精度好, 对于不同的结构适应性强.通过数值实验, 与商业软件和全三维(3D)DGTD方法的结果进行了比较和验证.Abstract: Power integrity (PI) problem is essential when analyzing high speed signal passing through power ground. The fundamental mode in power ground is the zero-order parallel plate mode, which is capable for 2D simplification. However, in areas around anti-pads and other z-axis discontinuities, 3D algorithm has to be adopted to improve the accuracy. A hybrid 2D/3D discontinuous Galerkin time domain (DGTD) method has advantage on both accuracy and efficiency, thus is effective to cope with such full wave simulations. The 2D and 3D domains share the same triangular prism mesh. With appropriated basis functions, different domains can couple with each other efficiently. As the electromagnetic wave propagates, the 2D/3D domain decomposition is dependent both on space and time. This work proposes an improved adaptive criterion which monitors the field at every time step, thus updates the 2D/3D domain dynamically. Unlike the state-of-the-art technique controlling absolute error, the proposed criterion controls comparative error. This method is both high efficient and accurate, and flexible for various kinds of structures. Some numerical examples are demonstrated to validate the proposed method. Comparisons with commercial software are also given.
-
引言
导孔经常被用作垂直连接通道, 应用于现代多层印制板和芯片中[1].随着电子设备工作频率的不断升高, 合适的功率分配网络(power distribution network, PDN)设计对确保信号/电源完整性(signal integrity/power integrity, SI/PI)十分重要.在电源板结构中, 主要的模式分布为零阶平行板模式.而由于沿z向的不连续性, 通过连接两个隔离板的导孔传播的高速电磁波将在隔离板附近激励起高次模.怎样高效地利用平行板模式, 同时合理地考虑高次模, 是一个值得讨论的问题.
许多学者利用电源板的平面性质来减少未知量从而提高效率.一种高效的二维三维混合有限元法(finite element method, FEM)方法[2]将几何区域分解为二维和三维区域.由于高次模随距离成指数衰减, 该方法采用了一种静态距离判据来区分二维和三维区域.
时域不连续伽辽金(discontinuous Galerkin time domain, DGTD)法[3-5]可以说是时域有限元(finite element time domain, FETD)方法[6]与时域有限体积(finite volume time domain, FVTD)方法[7]结合产生的混合算法.其采用非结构化进行空间离散, 用伽辽金测试方法建立线性方程组.引入数值流在不同单元之间建立联系, 对Maxwell方程要求在弱形式下成立, 然后沿着时间步推进获得全局解.由于很方便采用高阶多项式插值, DGTD数值色散误差小; 采用非结构化网格, 能够对任意复杂形状模拟进行拟合而不会有太大的阶梯近似误差; 同时, 由于采用数值流, DGTD可以对各个单元网格分别求解, 方便实现并行化, 并且不用像FETD方法在每一时间歩求解全局矩阵, 从而节约了大量的计算和存储资源.
一种二维的DGTD方法[8]最早用作PDN问题的分析.这种方法运用非结构化的网格, 只需更少的网格即可满足几何精度, 因此效率较高.但由于它假定结构在z方向的连续性, 因此不能用于存在隔离盘的结构.一种采用带波端口激励和辅助微分方程(auxiliary differential equation, ADE)方法并采用三棱柱网格的三维DGTD方法[9]具有高度的精度和灵活性[10], 但它没有利用电源板结构特有的平行板模式, 因此计算效率较低.
在我们最近的工作中, 提出了一种高效稳定的混合DGTD方法[11]解决平行板对瞬态问题.类似于文献[2], 整个结构被分解为高阶模区域(三维)和零阶平行板模式区域(二维).此外, 还引入ADE来考虑色散影响.这种算法采用并行化的方式来充分发挥DGTD的优点.为了控制在时间显式迭代过程中积累的二维近似误差, 我们采用了一种动态的自适应二维三维判据[12].
但是在文献[12]中, 所提出的判据采用了一个经验数值零来判断电磁场值是否可以忽略.这种控制绝对误差的方法对于一些应用场景是准确的, 但无法保证其适用于各种各样的结构.本文提出一种改进的自适应判据, 对二维近似带来的相对误差进行监视与控制, 从而对不同结构的适应性更强.我们将这种改进的自适应判据应用于一些数值算例, 用以验证本方法的精度、效率和适应性.
第二部分具体介绍了改进的带自适应相对误差控制判据的二维三维混合DGTD方法.第三部分通过数值算例来验证这种混合方法的精度, 效率和适应性.
1 采用方法
在这一部分, 首先介绍带ADE色散媒质的DGTD方法, 然后引入三棱柱网格单元(三维)和三角形网格单元(二维), 接着将它们在混合方法中组合起来, 最后, 提出改进后的控制相对误差的自适应判据和动态更新方法.
1.1 DGTD方法推导
首先考虑无源的麦克斯韦旋度方程:
\mu \frac{{\partial \mathit{\boldsymbol{H}}}}{{\partial t}} = - \nabla \times \mathit{\boldsymbol{E}}; (1) \varepsilon \frac{{\partial \mathit{\boldsymbol{E}}}}{{\partial t}} = \nabla \times \mathit{\boldsymbol{H}}. (2) 对于实际问题, 在电源-地之间的介质板是色散的, 假设为Debye模型, ε=ε0[ε∞+χ(ω)], 其色散影响用一个极化电流向量来表示:
{\mathit{\boldsymbol{J}}_q}\left( \omega \right) = \frac{{{\rm{j}}\omega {\varepsilon _0}\Delta {\varepsilon _q}}}{{1 + {\rm{j}}\omega {\tau _q}}}\mathit{\boldsymbol{E}}\left( \omega \right), (3) Δεq和τq分别表示介电常数的变化和弛豫时间.
它的逆傅里叶变换为
\frac{{{\rm{d}}{\mathit{\boldsymbol{J}}_q}}}{{{\rm{d}}t}} + \frac{{{\mathit{\boldsymbol{J}}_q}}}{{{\tau _q}}} = \frac{{{\varepsilon _0}\Delta {\varepsilon _q}}}{{{\tau _q}}}\frac{{{\rm{d}}\mathit{\boldsymbol{E}}}}{{{\rm{d}}t}}. (4) 基于Rankine Hugoniot关系, 数值通量用DGTD弱边界条件来表示:
\begin{array}{l} \mathit{\boldsymbol{\hat n}}_i^f \times \mathit{\boldsymbol{H}}_f^ * = \mathit{\boldsymbol{\hat n}}_i^f \times \frac{{\left( {{Z^i}{\mathit{\boldsymbol{H}}^i} + Z_f^j\mathit{\boldsymbol{H}}_f^i} \right) + \mathit{\boldsymbol{\hat n}}_i^f \times \left( {{\mathit{\boldsymbol{E}}^i} - \mathit{\boldsymbol{E}}_f^j} \right)}}{{{Z^i} + Z_f^i}} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\beta \cdot \frac{{\mathit{\boldsymbol{\hat n}}_i^f \times {\mathit{\boldsymbol{M}}_s}}}{{{Z^i} + Z_f^i}}; \end{array} (5) \begin{array}{l} \mathit{\boldsymbol{\hat n}}_i^f \times \mathit{\boldsymbol{E}}_f^ * = \mathit{\boldsymbol{\hat n}}_i^f \times \frac{{\left( {{Y^i}{\mathit{\boldsymbol{E}}^i} + Y_f^j\mathit{\boldsymbol{E}}_f^j} \right) + \mathit{\boldsymbol{\hat n}}_i^f \times \left( {\mathit{\boldsymbol{H}}_f^j - {\mathit{\boldsymbol{H}}^i}} \right)}}{{{Y^i} + Y_f^i}} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\beta \cdot \frac{{Y_f^j{\mathit{\boldsymbol{M}}_s}}}{{\left( {{Y^i} + Y_f^j} \right)}}. \end{array} (6) 式中:j是紧贴网格单元i表面f的相邻网格单元; f是相交的公共面; \mathit{\boldsymbol{\hat n}}_i^f代表f向外的单位向量; \mathit{\boldsymbol{\hat n}}_i^f \times \mathit{\boldsymbol{H}}_f^*和\mathit{\boldsymbol{\hat n}}_i^f \times \mathit{\boldsymbol{E}}_f^*是相邻元素用来交换信息以确保弱边界条件的数值通量; {\mathit{Z}^{i/j}} = \sqrt {{\mathit{\mu }^{i/j}}/{\mathit{\varepsilon }^{i/j}}} 和Yi/j=1/Zi/j分别代表了特征阻抗和特征导纳.如果表面f在波端口之上β=1, 否则, β=0.
将伽辽金测试方法应用于式(1), (2), (4), 并考虑式(5)、(6), 可以得到DGTD半离散矩阵方程:
\begin{array}{l} \mathit{\boldsymbol{M}}_e^i\frac{{{\rm{d}}{e^j}}}{{{\rm{d}}t}} = \mathit{\boldsymbol{S}}_e^i{h^i} + \sum\limits_{f = 1}^{n_i^f} {\left( {\mathit{\boldsymbol{F}}_{hh}^{ii,f}e_f^i + \mathit{\boldsymbol{F}}_{ee}^{ij,f}e_f^i + } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\mathit{\boldsymbol{F}}_{eh}^{ii,f}h_f^i + \mathit{\boldsymbol{F}}_{eh}^{ij,f}h_f^i} \right) + \beta \cdot \mathit{\boldsymbol{F}}_e^{i,{\mathit{\boldsymbol{M}}_s}}; \end{array} (7) \begin{array}{l} \mathit{\boldsymbol{M}}_h^i\frac{{{\rm{d}}{h^j}}}{{{\rm{d}}t}} = - \mathit{\boldsymbol{S}}_h^i{e^i} + \sum\limits_{f = 1}^{n_i^f} {\left( {\mathit{\boldsymbol{F}}_{hh}^{ii,f}h_f^i + \mathit{\boldsymbol{F}}_{ee}^{ij,f}h_f^i + } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\mathit{\boldsymbol{F}}_{he}^{ii,f}e_f^i + \mathit{\boldsymbol{F}}_{he}^{ij,f}e_f^i} \right) + \beta \cdot \mathit{\boldsymbol{F}}_h^{i,{\mathit{\boldsymbol{M}}_s}}; \end{array} (8) \mathit{\boldsymbol{M}}_q^i\frac{{{\rm{d}}{\gamma ^i}}}{{{\rm{d}}t}} + \frac{1}{{{\tau _q}}}\mathit{\boldsymbol{M}}_q^i{\gamma ^i} = \frac{{{\varepsilon _0}\Delta {\varepsilon _q}}}{{{\tau _q}}}\mathit{\boldsymbol{M}}_{c,q}^i\frac{{{\rm{d}}{e^i}}}{{{\rm{d}}t}}. (9) 式中:M, S和F分别是质量矩阵, 刚度矩阵和通量矩阵, 它们的表达式见文献[9].需要注意的是尽管三维和二维元素有同样的基函数和表达式, 但由于在式(7)~(9)的计算中三维算法采用的是面积积分, 而二维算法采用的是曲线积分, 所以针对三维算法和二维算法分别得到的矩阵是不同的.
这里运用四阶Runge-Kutta方法来求解式(7) ~(9)的半离散矩阵方程.对于显式的时间推进体系, 必须满足Courant-Freidrichs-Lewy (CFL)条件以确保稳定性.其最大时间步由文献[8]中公式确定:
{c_0}\delta t \le \min \left\{ {{l_{\min }}\sqrt {{\varepsilon _{\rm{r}}}{\mu _{\rm{r}}}} /\left[ {4{{\left( {p + 1} \right)}^2}} \right]} \right\} (10) 式中:c0是自由空间的波速; p是基函数的阶数.
1.2 波端口设置
不同于文献[12], 本文将采用正六边形的隔离盘作为输入端口.正六边形隔离盘的主模分布没有解析公式, 所以必须采用FEM方法求出主模的表面磁流分布.再用表面磁流Ms作为波端口激励, 并用作S参数提取.严格来说, 激励应该包含隔离盘支持的所有本征模.但为了简化, 同时考虑到高次模式的指数衰减, 本文只考虑了主要模式.类似地, S参数由扩展主模的系数而得到.
在波端口Pi, 电场和磁场的正切分量可以表达为本征模之和[9]:
\begin{array}{l} {\mathit{\boldsymbol{E}}_{\rm{t}}}\left( {\mathit{\boldsymbol{r}},\omega } \right) = a_0^{{\rm{TEM}}}\mathit{\boldsymbol{e}}_0^{{\rm{TEM}}} + \sum\limits_{n = 1}^\infty {a_n^{{\rm{TE}}}\mathit{\boldsymbol{e}}_n^{{\rm{TE}}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{n = 1}^\infty {a_n^{{\rm{TM}}}\mathit{\boldsymbol{e}}_n^{{\rm{TM}}}} ; \end{array} (11) \begin{array}{l} {\mathit{\boldsymbol{H}}_{\rm{t}}}\left( {\mathit{\boldsymbol{r}},\omega } \right) = b_0^{{\rm{TEM}}}\mathit{\boldsymbol{h}}_0^{{\rm{TEM}}} + \sum\limits_{n = 1}^\infty {b_n^{{\rm{TE}}}\mathit{\boldsymbol{h}}_n^{{\rm{TE}}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{n = 1}^\infty {b_n^{{\rm{TM}}}\mathit{\boldsymbol{h}}_n^{{\rm{TM}}}} . \end{array} (12) 式中:en和hn分别代表电场和磁场的本征模式; an和bn是相应的模式拓展系数.利用边界条件, 有bnTEM/TE/TM=-anTEM/TE/TM.
运用模式的正交性, 这些模式拓展系数可以计算出来:
a_n^{{\rm{TEM}}/{\rm{TE/TM}}} = \int {\mathit{\boldsymbol{\hat \xi }} \times {\mathit{\boldsymbol{M}}_s}/2} \cdot \mathit{\boldsymbol{\hat e}}_n^{{\rm{TEM}}/{\rm{TE}}/{\rm{TM}}}{\rm{d}}S. (13) 表面磁流激励为{\mathit{\boldsymbol{M}}_\mathit{s}} = - \mathit{\boldsymbol{\hat \xi }} \times {\mathit{\boldsymbol{E}}_\mathit{t}}, 其中\mathit{\boldsymbol{\hat \xi }} 代表指向波端口的单位正交向量.
1.3 基函数设置与二维三维混合
图 1表示了一种带不规则形状隔离盘的电源板的例子.整个结构根据模式, 被分为在隔离板周围的高次模区域(三维)和零阶平行板模式区域(二维).一般来说, 大多数的高次模存在于隔离板周围, 其他地方只存在零阶平行板模式.
整个结构被划分为若干三棱柱网格单元.因为三棱柱网格单元是由二维三角形网格单元简单地拉伸而成, 二者有许多共同点, 从而可以实现三维区域和简化的二维区域之间方便的连结.
\left\{ \begin{array}{l} \mathit{\boldsymbol{N}}_1^{{\rm{3D}}} = \zeta \left( {{L_1}\nabla {L_2} - {L_2}\nabla {L_1}} \right)\\ \mathit{\boldsymbol{N}}_2^{{\rm{3D}}} = \zeta \left( {{L_2}\nabla {L_3} - {L_3}\nabla {L_2}} \right)\\ \mathit{\boldsymbol{N}}_3^{{\rm{3D}}} = \zeta \left( {{L_3}\nabla {L_1} - {L_1}\nabla {L_3}} \right)\\ \mathit{\boldsymbol{N}}_4^{{\rm{3D}}} = \left( {1 - \zeta } \right)\left( {{L_4}\nabla {L_5} - {L_5}\nabla {L_4}} \right)\\ \mathit{\boldsymbol{N}}_5^{{\rm{3D}}} = \left( {1 - \zeta } \right)\left( {{L_5}\nabla {L_6} - {L_6}\nabla {L_5}} \right)\\ \mathit{\boldsymbol{N}}_6^{{\rm{3D}}} = \left( {1 - \zeta } \right)\left( {{L_6}\nabla {L_4} - {L_4}\nabla {L_6}} \right)\\ \mathit{\boldsymbol{N}}_7^{{\rm{3D}}} = {L_1}\nabla \zeta /\left| {\nabla \zeta } \right|\\ \mathit{\boldsymbol{N}}_8^{{\rm{3D}}} = {L_2}\nabla \zeta /\left| {\nabla \zeta } \right|\\ \mathit{\boldsymbol{N}}_9^{{\rm{3D}}} = {L_3}\nabla \zeta /\left| {\nabla \zeta } \right| \end{array} \right.. (14) 式中:Li代表节点i处的节点基函数; ζ在三棱柱高度上从0到1变化.作为一个整体, 每一个三维三棱柱需要18个未知数.
为了减少未知数的个数并且提高计算效率, 引入二维三角元来等价零阶平行板模式区域, 这些区域的电磁场表达式为
\left\{ \begin{array}{l} {e_{1:6}} = 0\\ {h_{1:3}} = {h_{4:6}}\\ {h_{7:9}} = 0 \end{array} \right.. (15) 二维三角元的基函数定义为[12]:
二维元的磁场基函数,
\left\{ \begin{array}{l} \mathit{\boldsymbol{N}}_1^{h,2{\rm{D}}} = \mathit{\boldsymbol{N}}_1^{3{\rm{D}}} = {L_1}\nabla {L_2} - {L_2}\nabla {L_1}\\ \mathit{\boldsymbol{N}}_2^{h,2{\rm{D}}} = \mathit{\boldsymbol{N}}_2^{3{\rm{D}}} = {L_2}\nabla {L_3} - {L_3}\nabla {L_2}\\ \mathit{\boldsymbol{N}}_3^{h,2{\rm{D}}} = \mathit{\boldsymbol{N}}_3^{3{\rm{D}}} = {L_3}\nabla {L_1} - {L_1}\nabla {L_3} \end{array} \right.. (16) 二维元的电场基函数,
\left\{ \begin{array}{l} \mathit{\boldsymbol{N}}_1^{e,2{\rm{D}}} = \mathit{\boldsymbol{N}}_7^{3{\rm{D}}} = {L_1}\nabla \zeta /\left| {\nabla \zeta } \right|\\ \mathit{\boldsymbol{N}}_2^{e,2{\rm{D}}} = \mathit{\boldsymbol{N}}_8^{3{\rm{D}}} = {L_2}\nabla \zeta /\left| {\nabla \zeta } \right|\\ \mathit{\boldsymbol{N}}_3^{e,2{\rm{D}}} = \mathit{\boldsymbol{N}}_9^{3{\rm{D}}} = {L_3}\nabla \zeta /\left| {\nabla \zeta } \right| \end{array} \right.. (17) 二维三角形网格单元只有6个未知量, 因此比三维元的效率更高.这里的基函数和相关的三维三棱柱的基函数相同, 因此在三维和二维之间很容易相互转换.
具体来说, 由式(7)~(9), 为了求解网格单元i内的电磁场分布, 需要知道它相邻的网格单元j的电磁场未知系数.如果一个网格单元i在三维或二维区域, i和它的所有相邻网格单元j都在三维或二维区域, 对于网格单元i, 则可以直接运用传统的三维二维DGTD方法处理.但如果一个网格单元i在三维和二维区域的交界处, 则它的所有相邻网格单元j应该被变换为和i性质相同的网格单元, 之后再采用传统的DGTD解法计算.例如, 如果i是二维单元, 则需要将其所相邻的所有网格单元j应用式(15)都转化为二维网格单元之后, 再对网格单元i采用二维DGTD方法求解, 反之亦然.具体的流程实现如图 3所示.
按照文献[12]中介绍的方法, 本二维三维混合DGTD方法也可用于多层网格的电源板结构的求解.并且在二维区域中, 多层网格可以进一步简化为单层网格, 从而进一步提高求解效率.
需要注意的是在求解过程中, 各个网格单元之间彼此相互独立, 因此很容易并行化.这里, 我们采用OpenMP来实现并行化.
1.4 自适应判据
在二维三维混合方法中, 判据是非常重要的.
判据决定了哪些单元可以被二维简化而哪些不能.一个良好的判据应该将二维简化带来的近似误差控制在可以接受的范围内, 从而保证精度与效率之间的平衡.文献[2]中选择了网格单元和隔离盘之间的距离作为静态判据:如果一个网格单元距离隔离板的距离比判据dmin小(大), 则这个网格单元为三(二)维元.这种判据根据的是隔离板的z向不连续性对远区网格单元几乎没有影响的物理直觉.然而, 这种距离判据的近似误差难以控制.在频域方法中这个误差可以忽略, 但在瞬态方法中, 小的误差可以随时间的迭代而积累从而导致不稳定结果.
文献[12]中提出了一种分解二维和三维区域的自适应判据.它由因果性和CFL限制推导出来, 即在每一个时间步每一个网格单元的电磁场只影响其相邻的网格单元.因此如果一个网格单元和它相邻的网格单元在时间步tn都属于二维区域, 那么在下一时间歩tn+1, 其他高次模式不会传播到这个网格单元.所以这个元就可以标记为二维网格单元并且在时间步tn+1用二维的DGTD方法处理.也就是说, 在一个网格单元被标记之前, 这个网格单元和它所有的相邻网格单元都要被检查.考虑到二维简化条件, 式(15)中的二维判据被设置为
\left\{ \begin{array}{l} \max \left( {{\rm{abs}}\left( {e_{1:6}^{i,{t_n}}} \right)} \right) < \eta \times \alpha \\ \max \left( {{\rm{abs}}\left( {h_{1:3}^{i,{t_n}} - h_{4:6}^{i,{t_n}}} \right)} \right) < \alpha \\ \max \left( {{\rm{abs}}\left( {h_{7:9}^{i,{t_n}}} \right)} \right) < \alpha \end{array} \right.. (18) 注意到上述判据中η为波阻抗\mathit{\eta = }\sqrt {\mathit{\mu /} \in } , 而α为一经验判据, 典型值为α=10-8.小于式(18)右侧值则认为其电磁场值为数值零, 可以忽略.
事实上, 由于电磁波随时间不断变化, 一个控制绝对误差的静态判据对于各种不同的结构和应用场景并不一定能够保证精度.激励的波幅, 波端口的尺寸和网格所处位置以及周围的边界条件等都会影响该网格的场值, 或可导致判据的失效.
本文提出一种改进的自适应判据, 对每个网格的相对误差进行控制, 从而增强二维三维混合DGTD方法对于不同应用场景的广泛适应性.
\left\{ \begin{array}{l} \max \left( {{\rm{abs}}\left( {e_{1:6}^i} \right)} \right) < \alpha \cdot \min \left( {{\rm{abs}}\left( {e_{7:9}^i} \right)} \right)\\ \max \left( {{\rm{abs}}\left( {h_{1:3}^i - H_{4:6}^i} \right)} \right) < \alpha \cdot \min \left( {{\rm{abs}}\left( {h_{1:3}^i} \right)} \right)\\ \max \left( {{\rm{abs}}\left( {h_{7:9}^i} \right)} \right) < \alpha \cdot \min \left( {{\rm{abs}}\left( {h_{1:6}^i} \right)} \right) \end{array} \right., (19) 左边部分在二维区域为数值零.这里定义数值零和文献[7]中所采用的经验值控制绝对误差不同, 本文中的数值零由每一个时间步的每一个网格单元的非零部分进行定义.在这种判据中, 二维相对近似误差被控制在小于误差控制系数α(典型值为0.1%).即如果该场值小于非零部分最小值的0.1%时, 则可认为该场值对该网格单元贡献很小, 可以忽略.这种改进的误差控制判据不依赖于时间、结构、激励的振幅等, 对一般性情况都能够保证精度和适应性.需要注意的是, 式(19)保留了DGTD的局部计算特性, 因此可以很容易地并行化.具体的判据实现流程如图 4所示.
2 数值实验
为了验证本文方法的精度、效率和适应性, 对图 5所示电源板结构进行分析.其中蓝色区域为正六边形和圆形的隔离盘, 其上面和下面分别为端口1, 2(正六边形)和端口3, 4(圆形).将a0TEM一个正弦调制的高斯脉冲加载在端口1上:
\begin{array}{l} a_0^{{\rm{TEM}}}\left( t \right) = \exp \left( { - {{\left( {t - {t_0}} \right)}^2}/{t_m}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\cos \left( {2{\rm{ \mathsf{ π} }}{f_m}\left( {t - {t_0}} \right)} \right). \end{array} (20) 式中:t0为参考时间参数; tm是控制频谱带宽的参数; fm是调制频率.
电源板的上下金属面之间间距0.254 mm, 正六边形隔离盘的边长为1/6 mm, 圆形隔离盘的边长为0.381 mm, 所有金属导孔柱的半径为0.127 mm.电源板中的介质视为Debye色散媒质, 其中∈∞=2.2 F/m, Δε1=2.0 F/m, τ1=2.5×10-10 s.为简明, 我们假设Debye模型q=1, 所有端口只考虑了TEM主模.其他几何参数见表 1.
mm 参数 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 y11 y12 y13 y14 数值 -3 3 6 6 4 -1 -3 -1 0 2 3 4 4 0 -3 -3 2 8 8 5 -1 -2 0 -1 2 3 5 3 为了对正六边形隔离盘进行激励和S参数提取, 需要先将其主模分布采用有限元方法计算出来.
图 6是通过有限元方法获得的正六边形隔离盘的电势分布.图 7是该正六边形隔离盘横截面的网格划分和表面磁流分布.将式(20)的调制高斯脉冲函数加载在表面磁流上作为激励, 并将端口电场按照该表面磁流分布提取出来即可获得正六边形隔离盘的主模时域响应与频域S参数.
![]()
图 7 正六边形隔离盘截面网格划分与表面磁流分布Fig. 7 Section grid division and surface magnetic flux distribution of a normal hexagon disks图 8是用本文方法得到的端口1入, 端口4出的时域幅度响应.随着电磁波的传播, 3D网格单元与2D网格单元是随时间动态变化的.所有网格单元初始设置为3D单元, 随后判据根据每个网格的电磁场值来判断其为3D或是2D单元.图 9中展示了3D单元数量随时间动态变化的过程.可以看到，3D单元很快减少到网格总数的3.5%左右并持续震荡.由于有大约96%的网格单元可以采用二维简化, 因此未知量更少, 效率更高.为了验证本文方法的精度, 我们将时域响应通过傅里叶变换到频域S参数, 并与商业软件HFSS进行对比.从图 10可以看到, 本文方法得到的S参数与采用频域FEM方法的HFSS的结果一致性较好.数值实验环境为笔记本电脑, 搭载4核8线程、主频2.6 GHz的Intel i6700HQ的CPU, 和16 GB的DDR3内存.本文二维三维混合方法的计算时间为213 s, 全三维DGTD方法为471 s, 4核并行HFSS耗时1 076 s.本文二维三维混合方法的高效率来自于二维简化带来的更少的未知量, 和DGTD方法本身的并行优势.
![]()
图 8 二维三维混合方法所得到时域幅度响应Fig. 8 Amplitude response of the time domain obtained by 2D/3D hybrid method![]()
图 9 2D/3D混合方法计算过程中3D单元随时间变化Fig. 9 Change of 3D unit with time during the calculation of 2D/3D mixing method![]()
图 10 采用不同方法获得的频域S参数对比Fig. 10 Comparison of S parameters in frequency domain obtained by different methods3 结论
本文提出了一种改进的二维三维混合DGTD方法用于分析带不规则隔离盘的色散平行电源板.由于二维元素的未知数个数只有三维元素的三分之一, 因此这种混合算法效率更高, 可以应用于高速信号电源板的快速仿真与求解.
基于之前工作[7], 我们提出一种改进的相对误差控制判据来保证精度.相比于之前控制绝对误差的动态判据, 本方法对于各种不同的结构和情况适应性更强.通过与商用软件进行对比, 检验了这种混合方法的精度、效率和适应性.
进一步的工作可能包括将本方法应用于更加复杂的电源板结构, 比如带去耦电容的电源板结构分析.
-
![]()
图 7 正六边形隔离盘截面网格划分与表面磁流分布
Fig. 7 Section grid division and surface magnetic flux distribution of a normal hexagon disks
![]()
图 8 二维三维混合方法所得到时域幅度响应
Fig. 8 Amplitude response of the time domain obtained by 2D/3D hybrid method
![]()
图 9 2D/3D混合方法计算过程中3D单元随时间变化
Fig. 9 Change of 3D unit with time during the calculation of 2D/3D mixing method
![]()
图 10 采用不同方法获得的频域S参数对比
Fig. 10 Comparison of S parameters in frequency domain obtained by different methods
-
[1] LI E P, WEI X C, CANGELLARIS A C, et al. Progress review of electromagnetic compatibility analysis technologies for packages, printed circuit boards, and novel interconnects[J]. IEEE transactions on electromagnetic compatibility, 2010, 52(2): 248-265. doi: 10.1109/TEMC.2010.2048755
[2] ZHANG Y J, REN L H, LIU D Z, et al. An efficient hybrid finite-element analysis of multiple vias sharing the same anti-pad in an arbitrarily shaped parallel-plate pair[J]. IEEE transactions on microwave theory and techniques, 2015, 63(3): 883-890. doi: 10.1109/TMTT.2015.2389257
[3] CHEN J, LIU Q H. Discontinuous Galerkin time-domain methods for multiscale electromagnetic simulations: a review[J]. Proceedings of the IEEE, 2013, 101(2): 242-254. doi: 10.1109/JPROC.2012.2219031
[4] ZHAO Y, DING D, CHEN R. A discontinuous galerkin time-domain integral equation method for electromagnetic scattering from PEC objects[J]. IEEE transactions on antennas and propagation, 2016, 64(6): 2410-2417. doi: 10.1109/TAP.2016.2550058
[5] TIAN C Y, SHI Y, LIU Z Q, et al. A Laguerre-based time-domain discontinuous Galerkin finite element-boundary integral method[J]. Microwave and optical technology letters, 2016, 58(11): 2774-2780. doi: 10.1002/mop.v58.11
[6] JIN J M. The finite element method in electromagnetics[M]. 3rd ed. New Jersey: John Wiley & Sons Inc, 2014.
[7] HOLLAND R, CABLE V P, WILSON L C. Finite-volume time-domain (FVTD) techniques for EM scattering[J]. IEEE transactions on electromagnetic compatibility, 1991, 33(4): 281-294. doi: 10.1109/15.99109
[8] LEE H M, GAO S, LIU E X, et al. Two-dimensional discontinuous Galerkin time-domain method for modeling of arbitrarily shaped power-ground planes[J]. IEEE transactions on electromagnetic compatibility, 2015, 57(6): 1744-1747. doi: 10.1109/TEMC.2015.2443019
[9] LI P, JIANG L J, BAĜCI H. Transient analysis of dispersive power-ground plate pairs with arbitrarily shaped antipads by the DGTD method with wave port excitation[J]. IEEE transactions on electromagnetic compatibility, 2017, 59(1): 172-183. doi: 10.1109/TEMC.2016.2596978
[10] LI P, JIANG L J, BAĜCI H. Discontinuous Galerkin time-domain analysis of power-ground planes taking into account decoupling capacitors[J]. IEEE transactions on components, packaging and manufacturing technology, 2017, PP(99): 1-10.
[11] MAI W, HU J. An efficient 2D/3D hybrid DGTD transient analysis for arbitrarily shaped anti-pads in power-ground plate-pair[C]//2017 International Applied Computational Electromagnetics Society Symposium. Florence, March 26-30, 2017: 1-2.
[12] W. Mai, J. Hu, P. Li, and H. Zhao, An efficient and stable 2-D/3-D hybrid discontinuous Galerkin time-domain analysis with adaptive criterion for arbitrarily shaped antipads in dispersive parallel-plate pair[J], IEEE Trans. Microw. Theory Tech., vol. 65, no. 10, pp. 3671-3681, Oct. 2017.
-
期刊类型引用(2)
1. 黄鹏，周远国，任强，王焱. 利用时域间断伽辽金算法求解Tellegen介质中的电磁波. 空军工程大学学报. 2023(05): 88-94 . 
百度学术
2. 平兰兰，吴东升，甘桂华. 时域有限元法电磁散射特性分析. 枣庄学院学报. 2020(05): 11-17 . 百度学术
其他类型引用(4)
下载:
计量
- 文章访问数: 165
- HTML全文浏览量: 92
- PDF下载量: 16
- 被引次数: 6
