基于Cole-Cole模型的DGTD电磁瞬态仿真方法与流程

基于cole
‑
cole模型的dgtd电磁瞬态仿真方法
技术领域
1.本发明属于特殊色散介质的电磁瞬态仿真技术领域，尤其涉及一种基于 cole
‑
cole模型的弱形式电磁瞬态仿真方法，具体来说，涉及间断伽辽金时域 (discontinuous galerkin time
‑
domain，dgtd)法仿真cole
‑
cole模型的方法。


背景技术：

2.在模拟电磁波的传播问题或研究宽频带内的电磁问题时，通常使用时域有限差分(finite difference time domain,fdtd)方法进行数值仿真。fdtd法是电磁场数值计算方法中应用最广的方法之一，具有算法简单、编程容易、物理过程描述直观、适用范围广、适合并行计算等诸多优点。无论是常微分方程、偏微分方程、高阶或非线性方程，都可通过差分近似微分的方式得到它们的代数方程组，最后计算出方程的近似解。近些年，应用fdtd法处理色散介质的电磁问题得到广泛关注。在描述色散介质的电参数随频率的变化关系时，通常使用诸如德拜(debye)、洛伦兹(lorentz)、德鲁徳(drude)等模型。这些传统的模型建立在整数阶微分方程的基础上，并使用常微分方程形式对物理现象进行数学描述。然而，很多介质(如生物组织)的介电特性随频率变化非常复杂。实验表明，与上述传统模型相比，基于分数微积分所推广出的分数阶模型更适合建模生物组织的电磁特性。其中，cole
‑
cole模型是常用的分数阶模型之一。在用fdtd法对其进行仿真时，cole
‑
cole模型存在的分数阶微积分是仿真的难点。为此人们提出多种解决方案，大致可分为以下三类：(1)利用辅助微分方程 (auxiliary differential equations,ades)法计算介质的极化；(2)利用时域卷积积分表示极化关系；(3)利用z变换法将cole
‑
cole表达式转换到z域再进行变化。
3.由于数值色散、数值稳定性以及阶梯误差的影响，fdtd方法的应用受到了限制。存在的不足主要体现在以下几个方面：(1)为避免fdtd的色散误差较大，在求解多尺度问题时,需要对整体计算区域进行细网格剖分，这导致未知量增加，所需计算内存增大；(2)受限于稳定性条件，即最大时间步长受最小空间步长的限制，使时间步长必须足够小才能精确模拟电小结构，这导致计算效率下降；(3)传统的fdtd法采用结构网格，在拟合目标外形时会产生阶梯误差；(4)针对阶梯误差的改进方法有高阶fdtd法、亚网格技术和亚像素平滑技术等，但这些方法都会降低计算效率。


技术实现要素：

4.为了解决上述已有技术存在的不足，本发明提出一种针对cole
‑
cole模型的复杂目标高性能时域电磁仿真算法，采用在仿真复杂多尺度目标上具有优势的dgtd法。 dgtd法可使用非结构化网格以避免fdtd法中由结构化网格引起的台阶误差，并且可将整个计算域划分为若干子域以减小计算内存。在保留fdtd法适用于复杂介质、易于并行计算的优点的同时，dgtd法可突破fdtd法在计算复杂几何目标方面的不足，解决fdtd法在计算复杂目标时计算效率低、计算精度差的问题，达到在保证精度的同时提高复杂目标的计算效率的目的。本发明提出的针对cole
‑
cole模型的dgtd 法，为后续研究基于分数阶模型的dgtd法打
下基础。本发明的具体技术方案如下：
5.一种基于cole
‑
cole模型的dgtd电磁瞬态仿真方法，包括以下步骤：
6.s1：构建基于cole
‑
cole模型的时域本构方程；
7.基于z变换法，应用欧拉定理及二阶泰勒级数展开，或双线性变换法及连分式展开，将频域cole
‑
cole模型转换到z域，并将分数次幂项转换为整数次幂，然后基于逆 z变换法，将z域cole
‑
cole模型带入电场强度和电位移矢量的z域本构方程，在时间 t＝nδt处将其转换到时域，得到电场强度和电位移矢量的时域本构方程，δt为时间步长；
8.s2：结合cole
‑
cole模型的eh
‑
dh混合模式的dgtd方案；
9.基于dgtd法，使用dh形式的修正的安培环路定律、eh形式的法拉第电磁感应定律以及步骤s1中的本构方程，结合迎风通量，构造半离散方程组；
10.s3：设计时间步进方案；
11.基于蛙跳法，在时间t＝nδt处对法拉第电磁感应定律的半离散方程进行中心差分，在t＝(n 1/2)δt处对修正的安培环路定律的半离散方程进行中心差分，得到磁场强度、电位移矢量、电场强度的时间步进方程。
12.进一步地，所述步骤s1的具体方法为：
13.s1
‑
1：cole
‑
cole模型模拟色散介质的复相对介电常数为：
[0014][0015]
其中，m为极数，ε0、ε
∞
、δε
m
分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度，j为虚数符号，j2＝
‑
1，ω为角频率，τ
m
为第m极的弛豫时间值，α
m
为第m极的色散参数，0≤α
m
＜1，σ为静态电导率；
[0016]
s1
‑
2：基于欧拉定理，将jω变换到z域：
[0017][0018]
其中，z
‑1＝e
‑
jω
，e为自然常数；将用二阶泰勒级数展开：
[0019][0020]
由此，式(1)变换得到z的整数次幂方程：
[0021][0022]
s1
‑
3：z域电场强度和z域电位移矢量的z域本构方程为：
[0023][0024]
为方便表达，定义z域的参数表达式为：
[0025]
[0026][0027]
则式(5)表示为：
[0028][0029]
将式(6)、式(7)在t＝nδt处进行逆z变换，n为时间步数，得到时域参数s
m
和 i的时间步进方程：
[0030]
i
n
＝σδt
·
e
n
i
n
‑1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0031][0032]
其中，i
n
为参数i在时间步数n处的定义，i
n
‑1为参数i在时间步数(n
‑
1)处的定义，s
m,n
为s
m
在时间步数n处的定义，s
m,n
‑1为s
m
在时间步数(n
‑
1)处的定义，s
m,n
‑2为 s
m
在时间步数(n
‑
2)处的定义；
[0033]
将式(8)在t＝nδt处进行逆z变换，得到时域电场强度e和时域电位移矢量d在时间步数n处的时域本构方程：
[0034][0035]
其中，e
n
为在时间步数n处的时域电场强度，d
n
为在时间步数n处的时域电位移矢量；
[0036]
为方便表述，定义参数a和参数v为：
[0037][0038][0039]
其中，v
n
为在时间步数n处的参数v，将式(12)、式(13)带入式(11)，整理得到：
[0040]
e
n
＝a
‑1·
(d
n
‑
v
n
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0041]
式(14)即基于cole
‑
cole模型的时域本构方程。
[0042]
进一步地，所述步骤s2的具体方法为：
[0043]
s2
‑
1：使用dh形式的修正的安培电路环路定律、eh形式的法拉第电磁感应定律，时域麦克斯韦方程组为：
[0044][0045][0046]
其中，h为时域磁场强度，j
f
为施加的电流密度，j
m
为施加的磁流密度，μ为磁导率；
[0047]
s2
‑
2：结合dgtd法，用符合旋度的向量基函数ψ
d
离散d，同时用符合散度的向量基函数φ
e
和φ
h
分别离散e和h，用符合散度的向量基函数φ
e
离散参数v、s
m
和 i，在子域k进行
离散得到子域k的maxwell方程和本构方程的弱形式为：
[0048][0049][0050][0051]
其中，d
k
为子域k中的d，e
k
为子域k中的e，h
k
为子域k中的h，e
t
为子域k及其所有邻域交界面处的e，h
t
为子域k及其所有邻域交界面处的h，为子域k中的j
f
，为子域k中的j
m
，μ
k
为子域k中的μ，a
k
为子域k中的a，v
k
为子域k中的v，为子域k中的ψ
d
，为子域k中的φ
e
，为子域k中的φ
h
，为k子域边界上向外的单位法向量；
[0052]
相邻子域用迎风通量耦合，子域k的麦克斯韦方程组和本构方程的半离散系统写为：
[0053][0054][0055][0056]
其中，l为k子域的邻域，n为与k子域相邻的子域的个数，e、h和d分别是离散的电场矢量、磁场矢量和电位移矢量；v为离散的参数v；e
(k)
、h
(k)
、d
(k)
和v
(k)
分别为子域k中的e、h、d和v；e
(l)
、h
(l)
分别为邻域l中的e、h；
[0057]
s2
‑
3：定义两个函数的内积为其中，v
e
是元胞的体积，(
·
)
t
为向量的转置，矩阵矩阵j
(k)
、m
(k)
的定义式为：
[0058][0059][0060][0061][0062]
[0063][0064][0065][0066]
其中，s
kl
为第k子域和第l子域的交界面，和y
k
＝1/z
k
分别为第k子域的阻抗和导纳，和y
l
＝1/z
l
分别为第l子域的阻抗和导纳；
[0067]
为k子域的相对介电常数，为l子域相对介电常数，为k子域的相对磁导率，为l子域的相对磁导率；
[0068]
当第k子域内介质为cole
‑
cole介质，则的值为该子域内cole
‑
cole模型的参数ε
∞
；当第l子域内介质为cole
‑
cole介质，则的值为该子域内cole
‑
cole模型的参数ε
∞
。
[0069]
进一步地，所述步骤s3中的具体方法为：
[0070]
采用蛙跳时间步进方案，在t＝nδt处对法拉第电磁感应定律安培电路环路定律的半离散方程(21)进行中心差分，在t＝(n 1/2)δt处对修正的安培环路定律的半离散方程(20)进行中心差分，得到线性方程组：
[0071][0072][0073][0074]
其中，和分别为在时间步数(n 1/2)和时间步数(n
‑
1/2)处的h
(k)
，和分别为在时间步数(n 1)和时间步数n处的d
(k)
，和分别为在时间步数(n 1) 和时间步数n处的e
(k)
，为在时间步数(n 1/2)处的j
(k)
，为在时间步数n处的 m
(k)
，为在时间步数(n 1)处的v
(k)
；
[0075]
整理得到时间步进顺序：
[0076][0077][0078][0079]
[0080][0081][0082]
其中，s
m
是用φ
e
离散参数s
m
得到的离散的矢量参数，i是用φ
e
离散参数i得到的离散的矢量参数，和分别为k域中在时间步数(n 1)、n和(n
‑
1)处的s
m
，和分别为k域中在时间步数(n 1)和n处的i。
[0083]
进一步地，所述步骤s1中为双线性变换法，则：
[0084]
step1
‑
1：cole
‑
cole模型模拟色散介质的复相对介电常数为：
[0085][0086]
其中，m为极数，ε0、ε
∞
、δε
m
分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度，j为虚数符号，j2＝
‑
1，ω为角频率，τ
m
为第m极的弛豫时间值，α
m
为第m极的色散参数，0≤α
m
＜1，σ为静态电导率；
[0087]
step1
‑
2：基于双线性变换法，将jω的分数次幂表达式变换到z域，然后进行连分式展开得到z的整数次幂方程：
[0088][0089]
其中，p
p,r
(z
‑1)和q
q,r
(z
‑1)分别是关于奇数次p和q的多项式，参数r＝1
‑
α
m
；
[0090]
由此，式(1)变换得到z的整数次幂方程：
[0091][0092]
step1
‑
3：z域电场强度和z域电位移矢量的z域本构方程为：
[0093][0094]
为方便表达，两个z域的参数的表达式变为：
[0095][0096][0097]
则式(34)用z域参数表示为：
[0098][0099]
为后续表达，定义多项式定义式中p
u,r
和q
l,w
为定义的系数，将式(35)、式(36)在t＝nδt处进行逆z变换，得到时域参数 s
m
和i的时间步进方程：
[0100][0101][0102]
其中，e
n
‑
w
为时间步数(n
‑
w)处的e，s
m,n
‑
w
为时间步数(n
‑
w)处的s
m,n
‑
w
，s
m,n
‑
u
为时间步数(n
‑
u)处的s
m,n
‑
u
；
[0103]
将式(37)在t＝nδt处进行逆z变换，得到时域电场强度e和时域电位移矢量d在时间步数n处的时域本构方程：
[0104][0105]
为方便表述，参数a和参数v为：
[0106][0107][0108]
将式(41)、式(42)带入式(40)，整理得到：
[0109]
e
n
＝a
‑1·
(d
n
‑
v
n
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(43)
[0110]
式(43)即基于cole
‑
cole模型的时域本构方程。
[0111]
本发明的有益效果在于：
[0112]
1.本发明构建了本构关系的离散化方程，以及采用dh和eh两种方案构建了麦克斯韦离散系统方程，实现了针对cole
‑
cole模型的dgtd法仿真方法，为用间断伽辽金时域法仿真分数阶模型奠定基础，进而解决色散介质的电磁问题。
[0113]
2.基于dgtd法，本发明可采用非共形网格对目标进行多尺寸、多种网格类型的剖分方案设计，在保留fdtd法适用于复杂介质、易于并行计算的优点的同时， dgtd法可突破fdtd法在计算复杂几何目标方面的不足，与fdtd法相比在计算多尺度问题上本发明计算效率高、所需计算内存少。
附图说明
[0114]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，对于本领域普通技术人员来讲，在不
付出创造性劳动的前提下，可以根据这些附图获得其他的附图。其中：
[0115]
图1为本发明数值仿真流程图；
[0116]
图2为脂肪块模型示意，其中，(a)为含激励源和观察点的脂肪块模型，(b)为脂肪块及其周围空气的六面体网格剖分；
[0117]
图3为dgtd法的频域结果与feko的比较，其中，(a)为观察点1处的ex总场，(b) 为观察点2处的ex总场；
[0118]
图4为乳腺模型示意图，其中(a)为含激励源、观察点和接收面的乳腺模型，(b)为乳腺及其周围空气的四面体网格剖分；
[0119]
图5为观察点处dgtd法的频域结果与feko的比较，其中，(a)为观察点1处的总场ex，(b)为观察点2处的总场ex，(c)为观察点3处的总场ex，(d)为观察点4处的总场 ex；
[0120]
图6为在6个频点上，接收面处dgtd法的频域结果与feko的相对误差，其中， (a)为1ghz频点上的相对误差，(b)为2ghz频点上的相对误差，(c)为3ghz频点上的相对误差，(d)为4ghz频点上的相对误差，(e)为5ghz频点上的相对误差， (f)为6ghz频点上的相对误差。
具体实施方式
[0121]
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
[0122]
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
[0123]
本发明的基于cole
‑
cole模型的dgtd仿真流程图如图1所示，以dgtd法为主要研究手段，具体实现过程如下：
[0124]
首先，在构建基于cole
‑
cole模型的时域本构方程。采用z变换法，结合欧拉定理和二阶泰勒级数展开或者双线性变换和连续分数阶展开，处理时域cole
‑
cole模型中存在的分数阶微积分。
[0125]
s1
‑
1：cole
‑
cole模型模拟色散介质的复相对介电常数为：
[0126][0127]
其中，m为极数，ε0、ε
∞
、δε
m
分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度，j为虚数符号，j2＝
‑
1，ω为角频率，τ
m
为第m极的弛豫时间值，α
m
为第m极的色散参数，0≤α
m
＜1，σ为静态电导率；
[0128]
s1
‑
2：基于欧拉定理，将jω变换到z域：
[0129][0130]
其中，z
‑1＝e
‑
jω
，e为自然常数；将用二阶泰勒级数展开：
[0131][0132]
由此，式(1)变换得到z的整数次幂方程：
[0133][0134]
s1
‑
3：z域电场强度和z域电位移矢量的z域本构方程为：
[0135][0136]
为方便表达，定义z域的参数表达式为：
[0137][0138][0139]
则式(5)表示为：
[0140][0141]
将式(6)、式(7)在t＝nδt处进行逆z变换，n为时间步数，得到时域参数s
m
和 i的时间步进方程：
[0142]
i
n
＝σδt
·
e
n
i
n
‑1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0143][0144]
其中，i
n
为参数i在时间步数n处的定义，i
n
‑1为参数i在时间步数(n
‑
1)处的定义，s
m,n
为s
m
在时间步数n处的定义，s
m,n
‑1为s
m
在时间步数(n
‑
1)处的定义，s
m,n
‑2为 s
m
在时间步数(n
‑
2)处的定义；
[0145]
将式(8)在t＝nδt处进行逆z变换，得到时域电场强度e和时域电位移矢量d在时间步数n处的时域本构方程：
[0146][0147]
其中，e
n
为在时间步数n处的时域电场强度，d
n
为在时间步数n处的时域电位移矢量；
[0148]
为方便表述，定义参数a和参数v为：
[0149][0150][0151]
其中，v
n
为在时间步数n处的参数v，将式(12)、式(13)带入式(11)，整理得到：
[0152]
e
n
＝a
‑1·
(d
n
‑
v
n
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0153]
式(14)即基于cole
‑
cole模型的时域本构方程。
[0154]
在一些实施方式中，步骤s1中用双线性变换法替换欧拉方法，则：
[0155]
step1
‑
1：cole
‑
cole模型模拟色散介质的复相对介电常数为：
[0156][0157]
其中，m为极数，ε0、ε
∞
、δε
m
分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度，j为虚数符号，j2＝
‑
1，ω为角频率，τ
m
为第m极的弛豫时间值，α
m
为第m极的色散参数，0≤α
m
＜1，σ为静态电导率；
[0158]
step1
‑
2：基于双线性变换法，将jω的分数次幂表达式变换到z域，然后进行连分式展开得到z的整数次幂方程：
[0159][0160]
其中，p
p,r
(z
‑1)和q
q,r
(z
‑1)分别是关于奇数次p和q的多项式，参数r＝1
‑
α
m
，其具体表达式如表1所示；
[0161]
表1 p
p,r
(z
‑1)和q
q,r
(z
‑1)的具体表达式
[0162][0163]
由此，式(1)变换得到z的整数次幂方程：
[0164][0165]
step1
‑
3：z域电场强度和z域电位移矢量的z域本构方程为：
[0166][0167]
为方便表达，两个z域的参数的表达式变为：
[0168][0169][0170]
则式(17)用z域参数表示为：
[0171][0172]
为后续表达，定义多项式定义式中p
u,r
和q
l,w
为定义的系数，将式(18)、式(19)在t＝nδt处进行逆z变换，得到时域参数 s
m
和i的时间步进方程：
[0173][0174][0175]
其中，e
n
‑
w
为时间步数(n
‑
w)处的e，s
m,n
‑
w
为时间步数(n
‑
w)处的s
m,n
‑
w
，s
m,n
‑
u
为时间步数(n
‑
u)处的s
m,n
‑
u
；
[0176]
将式(20)在t＝nδt处进行逆z变换，得到时域电场强度e和时域电位移矢量d在时间步数n处的时域本构方程：
[0177][0178]
为方便表述，参数a和参数v为：
[0179][0180][0181]
将式(24)、式(25)带入式(23)，整理得到：
[0182]
e
n
＝a
‑1·
(d
n
‑
v
n
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(26)
[0183]
式(26)即基于cole
‑
cole模型的时域本构方程。
[0184]
其次，结合cole
‑
cole模式的dgtd方案。
[0185]
s2
‑
1：使用dh形式的修正的安培电路环路定律、eh形式的法拉第电磁感应定律，时域麦克斯韦方程组为：
[0186][0187][0188]
其中，h为时域磁场强度，j
f
为施加的电流密度，j
m
为施加的磁流密度，μ为磁导率；
[0189]
s2
‑
2：结合dgtd法，用符合旋度的向量基函数ψ
d
离散d，同时用符合散度的向量基函数φ
e
和φ
h
分别离散e和h，用符合散度的向量基函数φ
e
离散参数v、s
m
和 i，在子域k进行离散得到子域k的maxwell方程和本构方程的弱形式为：
[0190][0191][0192][0193]
其中，d
k
为子域k中的d，e
k
为子域k中的e，h
k
为子域k中的h，e
t
为子域 k及其所有邻域交界面处的e，h
t
为子域k及其所有邻域交界面处的h，为子域k 中的j
f
，为子域k中的j
m
，μ
k
为子域k中的μ，a
k
为子域k中的a，v
k
为子域 k中的v，为子域k中的ψ
d
，为子域k中的φ
e
，为子域k中的φ
h
，为k 子域边界上向外的单位法向量；
[0194]
相邻子域用迎风通量耦合，子域k的麦克斯韦方程组和本构方程的半离散系统写为：
[0195][0196][0197][0198]
其中，l为k子域的邻域，n为与k子域相邻的子域的个数，e、h和d分别是离散的电场矢量、磁场矢量和电位移矢量；v为离散的参数v；e
(k)
、h
(k)
、d
(k)
和v
(k)
分别为子域k中的e、h、d和v；e
(l)
、h
(l)
分别为邻域l中的e、h；
[0199]
s2
‑
3：定义两个函数的内积为其中，v
e
是元胞的体积，(
·
)
t
为向量的转置，矩阵矩阵j
(k)
、m
(k)
的定义式为：
[0200][0201][0202][0203][0204][0205][0206][0207][0208]
其中，s
kl
为第k子域和第l子域的交界面，和y
k
＝1/z
k
分别为第k子域的阻抗和导纳，和y
l
＝1/z
l
分别为第l子域的阻抗和导纳；
[0209]
为k子域的相对介电常数，为l子域相对介电常数，为k子域的相对磁导率，为l子域的相对磁导率；
[0210]
当k子域内介质为cole
‑
cole介质，则的值为该子域内cole
‑
cole模型的参数ε
∞
；当l子域内介质为cole
‑
cole介质，则的值为该子域内cole
‑
cole模型的参数ε
∞
。
[0211]
第三步，时间步进方案设计。
[0212]
采用蛙跳时间步进方案，在t＝nδt处对法拉第电磁感应定律的半离散方程(33) 进行中心差分，在t＝(n 1/2)δt处对修正的安培环路定律的半离散方程(32)进行中心差分，得到线性方程组：
[0213][0214][0215][0216]
其中，和分别为在时间步数(n 1/2)和时间步数(n
‑
1/2)处的h
(k)
，和分别为在时间步数(n 1)和时间步数n处的d
(k)
，和分别为在时间步数(n 1) 和时间步数n处的e
(k)
，为在时间步数(n 1/2)处的j
(k)
，为在时间步数n处的 m
(k)
，为在时间步数(n 1)处的v
(k)
；
[0217]
整理得到时间步进顺序：
[0218][0219][0220][0221][0222][0223][0224]
其中，s
m
是用φ
e
离散参数s
m
得到的离散的矢量参数，i是用φ
e
离散参数i得到的离散的矢量参数，和分别为k域中在时间步数(n 1)、n和(n
‑
1)处的s
m
，和分别为k域中在时间步数(n 1)和n处的i。
[0225]
在一些实施方式中，当步骤s1中的方法为双线性变换法时，步骤s3中参数具体为：采用蛙跳时间步进方案，得到时间步进顺序如下：
[0226][0227][0228][0229][0230]
[0231][0232]
其中，代表时间步数(n 1
‑
w)处的e
(k)
，和分别代表k域中时间步数(n 1
‑
w)和(n 1
‑
u)处的s
m
。
[0233]
下面通过具体实施例，验证本发明方法的有效性。
[0234]
实施例1
[0235]
为了验证本发明提出的基于cole
‑
cole模型的dgtd算法的有效性，选用基于欧拉定律的dgtd方法，用matlab编程，进行数值仿真，主要仿真过程如下：
[0236]
(1)参数设置
[0237]
选择一个尺寸为30mm
×
30mm
×
20mm的脂肪块，中心位于原点，相对磁导率μ
r
＝1，脂肪的cole
‑
cole模型的电参数为ε
∞
＝2.32，σ＝0.0222，δε1＝2.14，δε2＝19.7，α1＝0.002，α2＝0.259，f1＝20.1ghz，f2＝8.25mhz。
[0238]
平面波的传播方向为[0,0,1]，电场极化方向为[1,0,0]。宽带平面波的激励类型是一阶blackman
‑
harris窗(bhw)脉冲，其中心频率为10ghz。在[0,0,
‑
15]mm和 [20,0,0]mm处分别设置观察点1和观察点2，记录该点的总场，如图2(a)所示。脂肪块和周围的空气用六面体网格剖分，网格尺寸分别为2mm和5mm，如图2(b)所示。完美匹配层用高阶六面体网格剖分，时间步长设置为0.125ps。
[0239]
(2)结果分析
[0240]
为了验证所提出的dgtd方案的计算精度，将计算结果转换到频域，并与商用电磁软件feko进行比较，如图3(a)
‑
(b)所示。
[0241]
结果表明，dgtd算法的结果与feko的结果非常吻合，验证了针对cole
‑
cole 模型的基于欧律定律的dgtd方案仿真的准确性，为后续用dgtd计算分数阶介质、解决fdtd问题做铺垫。
[0242]
实施例2
[0243]
为了验证本发明提出的基于cole
‑
cole模型的dgtd算法的有效性，选用的是基于双线性变换的dgtd方法，用matlab编程，进行数值仿真，主要仿真过程如下：
[0244]
(1)参数设置
[0245]
构建简单乳腺模型。选择一个半径为50mm，底面圆心位于原点的半球，为乳房介质，其电参数为ε
∞
＝2.5，σ＝0.01，δε1＝3，δε2＝15，δε3＝50000，δε4＝20000000，α1＝0.1，α2＝0.1，α3＝0.1，α4＝0，f1＝9.002ghz，f2＝2.5001mhz，f3＝350.0219khz， f4＝12.0026hz。内部包含一个半径为8mm的球体，球心位于点[14,14,20]mm，该球体为肿瘤介质，其电参数为ε
∞
＝4，σ＝0.2，δε1＝50，α1＝0，f1＝22.736ghz。
[0246]
平面波的传播方向为[0,0,
‑
1]，电场极化方向为[1,0,0]。宽带平面波的激励类型是一阶blackman
‑
harris窗(bhw)脉冲，其中心频率为3.1ghz。在[40,40,70]mm、 [
‑
40,
‑
40,70]mm、[30,30,80]mm和[
‑
30,
‑
30,80]mm处分别设置观察点1、观察点2、观察点2、观察点4，记录该点的总场并在z＝20mm处设一个接收面接收电磁场，该接收面平行于xy
‑
平面，其中心过z轴且位于z＝20mm处，尺寸为48mm
×
48mm，如图4(a)所示，如图脂肪块和周围的空气用四面体网格剖分，如图4(b)所示。完美匹配层用高阶六面体网格剖分，时间步长设置为0.4ps。
[0247]
(2)结果分析
[0248]
为了验证本发明的dgtd方案的计算精度，将计算结果转换到频域，并与商用电磁软件feko进行比较，观察点1、观察点2、观察点3、观察点4的结果如图5 (a)
‑
(d)所示。在频率f＝1、2、3、4、5、6ghz处，接收面上dgtd方案与feko结果的相对误差如图6(a)
‑
(f)所示。
[0249]
结果表明，dgtd算法的结果与feko的结果非常吻合，验证了针对cole
‑
cole 模型的基于双线性变换法的dgtd方案仿真的准确性。
[0250]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。