CN113536626A - 基于Cole-Cole模型的DGTD电磁瞬态仿真方法 - Google Patents
基于Cole-Cole模型的DGTD电磁瞬态仿真方法 Download PDFInfo
- Publication number
- CN113536626A CN113536626A CN202110709190.4A CN202110709190A CN113536626A CN 113536626 A CN113536626 A CN 113536626A CN 202110709190 A CN202110709190 A CN 202110709190A CN 113536626 A CN113536626 A CN 113536626A
- Authority
- CN
- China
- Prior art keywords
- cole
- domain
- time
- equation
- subfield
- 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
- 238000000034 method Methods 0.000 title claims abstract description 74
- 238000004088 simulation Methods 0.000 title claims abstract description 24
- 230000001052 transient effect Effects 0.000 title claims abstract description 10
- 239000013598 vector Substances 0.000 claims description 51
- 230000005684 electric field Effects 0.000 claims description 26
- 238000006073 displacement reaction Methods 0.000 claims description 21
- 239000002609 medium Substances 0.000 claims description 19
- 230000009466 transformation Effects 0.000 claims description 18
- 238000011426 transformation method Methods 0.000 claims description 15
- 230000006870 function Effects 0.000 claims description 12
- 239000006185 dispersion Substances 0.000 claims description 11
- 230000005674 electromagnetic induction Effects 0.000 claims description 10
- 230000035699 permeability Effects 0.000 claims description 10
- 230000004907 flux Effects 0.000 claims description 8
- 239000002612 dispersion medium Substances 0.000 claims description 7
- 230000003068 static effect Effects 0.000 claims description 5
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 13
- 238000013461 design Methods 0.000 abstract description 2
- 230000008901 benefit Effects 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 6
- 230000008569 process Effects 0.000 description 5
- 210000000481 breast Anatomy 0.000 description 4
- 230000007547 defect Effects 0.000 description 4
- 230000005284 excitation Effects 0.000 description 4
- 230000010287 polarization Effects 0.000 description 4
- 230000005672 electromagnetic field Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 210000001519 tissue Anatomy 0.000 description 2
- 206010028980 Neoplasm Diseases 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 210000005075 mammary gland Anatomy 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 108020001568 subdomains Proteins 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computing Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于Cole‑Cole模型的DGTD电磁瞬态仿真方法,包括构建基于Cole‑Cole模型的时域本构方程、结合Cole‑Cole模型的EH‑DH混合模式的DGTD方案和设计时间步进方案。本发明构建了本构关系的离散化方程,采用DH‑EH混合方案构建了麦克斯韦离散系统方程,实现了针对Cole‑Cole模型的DGTD法仿真方法,为未来用间断伽辽金时域法仿真分数阶模型奠定基础。本发明可采用非共形网格对目标进行多尺寸、多种网格类型的剖分方案设计,与FDTD法相比在计算多尺度问题上本发明计算效率高、所需计算内存少。
Description
技术领域
本发明属于特殊色散介质的电磁瞬态仿真技术领域,尤其涉及一种基于 Cole-Cole模型的弱形式电磁瞬态仿真方法,具体来说,涉及间断伽辽金时域 (discontinuousGalerkin time-domain,DGTD)法仿真Cole-Cole模型的方法。
背景技术
在模拟电磁波的传播问题或研究宽频带内的电磁问题时,通常使用时域有限差分(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域再进行变化。
由于数值色散、数值稳定性以及阶梯误差的影响,FDTD方法的应用受到了限制。存在的不足主要体现在以下几个方面:(1)为避免FDTD的色散误差较大,在求解多尺度问题时,需要对整体计算区域进行细网格剖分,这导致未知量增加,所需计算内存增大;(2)受限于稳定性条件,即最大时间步长受最小空间步长的限制,使时间步长必须足够小才能精确模拟电小结构,这导致计算效率下降;(3)传统的FDTD法采用结构网格,在拟合目标外形时会产生阶梯误差;(4)针对阶梯误差的改进方法有高阶FDTD法、亚网格技术和亚像素平滑技术等,但这些方法都会降低计算效率。
发明内容
为了解决上述已有技术存在的不足,本发明提出一种针对Cole-Cole模型的复杂目标高性能时域电磁仿真算法,采用在仿真复杂多尺度目标上具有优势的DGTD法。 DGTD法可使用非结构化网格以避免FDTD法中由结构化网格引起的台阶误差,并且可将整个计算域划分为若干子域以减小计算内存。在保留FDTD法适用于复杂介质、易于并行计算的优点的同时,DGTD法可突破FDTD法在计算复杂几何目标方面的不足,解决FDTD法在计算复杂目标时计算效率低、计算精度差的问题,达到在保证精度的同时提高复杂目标的计算效率的目的。本发明提出的针对Cole-Cole模型的DGTD 法,为后续研究基于分数阶模型的DGTD法打下基础。本发明的具体技术方案如下:
一种基于Cole-Cole模型的DGTD电磁瞬态仿真方法,包括以下步骤:
S1:构建基于Cole-Cole模型的时域本构方程;
基于Z变换法,应用欧拉定理及二阶泰勒级数展开,或双线性变换法及连分式展开,将频域Cole-Cole模型转换到Z域,并将分数次幂项转换为整数次幂,然后基于逆 Z变换法,将Z域Cole-Cole模型带入电场强度和电位移矢量的Z域本构方程,在时间 t=nΔt处将其转换到时域,得到电场强度和电位移矢量的时域本构方程,Δt为时间步长;
S2:结合Cole-Cole模型的EH-DH混合模式的DGTD方案;
基于DGTD法,使用DH形式的修正的安培环路定律、EH形式的法拉第电磁感应定律以及步骤S1中的本构方程,结合迎风通量,构造半离散方程组;
S3:设计时间步进方案;
基于蛙跳法,在时间t=nΔt处对法拉第电磁感应定律的半离散方程进行中心差分,在t=(n+1/2)Δt处对修正的安培环路定律的半离散方程进行中心差分,得到磁场强度、电位移矢量、电场强度的时间步进方程。
进一步地,所述步骤S1的具体方法为:
S1-1:Cole-Cole模型模拟色散介质的复相对介电常数为:
其中,M为极数,ε0、ε∞、Δεm分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度,j为虚数符号,j2=-1,ω为角频率,τm为第m极的弛豫时间值,αm为第m极的色散参数,0≤αm<1,σ为静态电导率;
S1-2:基于欧拉定理,将jω变换到Z域:
由此,式(1)变换得到z的整数次幂方程:
则式(5)表示为:
将式(6)、式(7)在t=nΔt处进行逆Z变换,n为时间步数,得到时域参数Sm和 I的时间步进方程:
In=σΔt·En+In-1 (9)
其中,In为参数I在时间步数n处的定义,In-1为参数I在时间步数(n-1)处的定义,Sm,n为Sm在时间步数n处的定义,Sm,n-1为Sm在时间步数(n-1)处的定义,Sm,n-2为 Sm在时间步数(n-2)处的定义;
将式(8)在t=nΔt处进行逆Z变换,得到时域电场强度E和时域电位移矢量D在时间步数n处的时域本构方程:
其中,En为在时间步数n处的时域电场强度,Dn为在时间步数n处的时域电位移矢量;
为方便表述,定义参数A和参数V为:
其中,Vn为在时间步数n处的参数V,将式(12)、式(13)带入式(11),整理得到:
En=A-1·(Dn-Vn) (14)
式(14)即基于Cole-Cole模型的时域本构方程。
进一步地,所述步骤S2的具体方法为:
S2-1:使用DH形式的修正的安培电路环路定律、EH形式的法拉第电磁感应定律,时域麦克斯韦方程组为:
其中,H为时域磁场强度,Jf为施加的电流密度,JM为施加的磁流密度,μ为磁导率;
S2-2:结合DGTD法,用符合旋度的向量基函数ΨD离散D,同时用符合散度的向量基函数ΦE和ΦH分别离散E和H,用符合散度的向量基函数ΦE离散参数V、Sm和 I,在子域k进行离散得到子域k的Maxwell方程和本构方程的弱形式为:
其中,Dk为子域k中的D,Ek为子域k中的E,Hk为子域k中的H,Et为子域k及其所有邻域交界面处的E,Ht为子域k及其所有邻域交界面处的H,为子域k中的Jf,为子域k中的JM,μk为子域k中的μ,Ak为子域k中的A,Vk为子域k中的V,为子域k中的ΨD,为子域k中的ΦE,为子域k中的ΦH,为k子域边界上向外的单位法向量;
相邻子域用迎风通量耦合,子域k的麦克斯韦方程组和本构方程的半离散系统写为:
其中,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;
进一步地,所述步骤S3中的具体方法为:
采用蛙跳时间步进方案,在t=nΔt处对法拉第电磁感应定律安培电路环路定律的半离散方程(21)进行中心差分,在t=(n+1/2)Δt处对修正的安培环路定律的半离散方程(20)进行中心差分,得到线性方程组:
其中,和分别为在时间步数(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);
整理得到时间步进顺序:
其中,sm是用ΦE离散参数Sm得到的离散的矢量参数,i是用ΦE离散参数I得到的离散的矢量参数,和分别为k域中在时间步数(n+1)、n和(n-1)处的sm,和分别为k域中在时间步数(n+1)和n处的i。
进一步地,所述步骤S1中为双线性变换法,则:
Step1-1:Cole-Cole模型模拟色散介质的复相对介电常数为:
其中,M为极数,ε0、ε∞、Δεm分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度,j为虚数符号,j2=-1,ω为角频率,τm为第m极的弛豫时间值,αm为第m极的色散参数,0≤αm<1,σ为静态电导率;
Step1-2:基于双线性变换法,将jω的分数次幂表达式变换到Z域,然后进行连分式展开得到z的整数次幂方程:
其中,Pp,r(z-1)和Qq,r(z-1)分别是关于奇数次p和q的多项式,参数r=1-αm;
由此,式(1)变换得到z的整数次幂方程:
其中,En-w为时间步数(n-w)处的E,Sm,n-w为时间步数(n-w)处的Sm,n-w,Sm,n-u为时间步数(n-u)处的Sm,n-u;
将式(37)在t=nΔt处进行逆Z变换,得到时域电场强度E和时域电位移矢量D在时间步数n处的时域本构方程:
为方便表述,参数A和参数V为:
将式(41)、式(42)带入式(40),整理得到:
En=A-1·(Dn-Vn) (43)
式(43)即基于Cole-Cole模型的时域本构方程。
本发明的有益效果在于:
1.本发明构建了本构关系的离散化方程,以及采用DH和EH两种方案构建了麦克斯韦离散系统方程,实现了针对Cole-Cole模型的DGTD法仿真方法,为用间断伽辽金时域法仿真分数阶模型奠定基础,进而解决色散介质的电磁问题。
2.基于DGTD法,本发明可采用非共形网格对目标进行多尺寸、多种网格类型的剖分方案设计,在保留FDTD法适用于复杂介质、易于并行计算的优点的同时, DGTD法可突破FDTD法在计算复杂几何目标方面的不足,与FDTD法相比在计算多尺度问题上本发明计算效率高、所需计算内存少。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,通过参考附图会更加清楚的理解本发明的特征和优点,附图是示意性的而不应理解为对本发明进行任何限制,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,可以根据这些附图获得其他的附图。其中:
图1为本发明数值仿真流程图;
图2为脂肪块模型示意,其中,(a)为含激励源和观察点的脂肪块模型,(b)为脂肪块及其周围空气的六面体网格剖分;
图3为DGTD法的频域结果与FEKO的比较,其中,(a)为观察点1处的Ex总场,(b) 为观察点2处的Ex总场;
图4为乳腺模型示意图,其中(a)为含激励源、观察点和接收面的乳腺模型,(b)为乳腺及其周围空气的四面体网格剖分;
图5为观察点处DGTD法的频域结果与FEKO的比较,其中,(a)为观察点1处的总场Ex,(b)为观察点2处的总场Ex,(c)为观察点3处的总场Ex,(d)为观察点4处的总场 Ex;
图6为在6个频点上,接收面处DGTD法的频域结果与FEKO的相对误差,其中, (a)为1GHz频点上的相对误差,(b)为2GHz频点上的相对误差,(c)为3GHz频点上的相对误差,(d)为4GHz频点上的相对误差,(e)为5GHz频点上的相对误差, (f)为6GHz频点上的相对误差。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点,下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是,在不冲突的情况下,本发明的实施例及实施例中的特征可以相互组合。
在下面的描述中阐述了很多具体细节以便于充分理解本发明,但是,本发明还可以采用其他不同于在此描述的其他方式来实施,因此,本发明的保护范围并不受下面公开的具体实施例的限制。
本发明的基于Cole-Cole模型的DGTD仿真流程图如图1所示,以DGTD法为主要研究手段,具体实现过程如下:
首先,在构建基于Cole-Cole模型的时域本构方程。采用Z变换法,结合欧拉定理和二阶泰勒级数展开或者双线性变换和连续分数阶展开,处理时域Cole-Cole模型中存在的分数阶微积分。
S1-1:Cole-Cole模型模拟色散介质的复相对介电常数为:
其中,M为极数,ε0、ε∞、Δεm分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度,j为虚数符号,j2=-1,ω为角频率,τm为第m极的弛豫时间值,αm为第m极的色散参数,0≤αm<1,σ为静态电导率;
S1-2:基于欧拉定理,将jω变换到Z域:
由此,式(1)变换得到z的整数次幂方程:
则式(5)表示为:
将式(6)、式(7)在t=nΔt处进行逆Z变换,n为时间步数,得到时域参数Sm和 I的时间步进方程:
In=σΔt·En+In-1 (9)
其中,In为参数I在时间步数n处的定义,In-1为参数I在时间步数(n-1)处的定义,Sm,n为Sm在时间步数n处的定义,Sm,n-1为Sm在时间步数(n-1)处的定义,Sm,n-2为 Sm在时间步数(n-2)处的定义;
将式(8)在t=nΔt处进行逆Z变换,得到时域电场强度E和时域电位移矢量D在时间步数n处的时域本构方程:
其中,En为在时间步数n处的时域电场强度,Dn为在时间步数n处的时域电位移矢量;
为方便表述,定义参数A和参数V为:
其中,Vn为在时间步数n处的参数V,将式(12)、式(13)带入式(11),整理得到:
En=A-1·(Dn-Vn) (14)
式(14)即基于Cole-Cole模型的时域本构方程。
在一些实施方式中,步骤S1中用双线性变换法替换欧拉方法,则:
Step1-1:Cole-Cole模型模拟色散介质的复相对介电常数为:
其中,M为极数,ε0、ε∞、Δεm分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度,j为虚数符号,j2=-1,ω为角频率,τm为第m极的弛豫时间值,αm为第m极的色散参数,0≤αm<1,σ为静态电导率;
Step1-2:基于双线性变换法,将jω的分数次幂表达式变换到Z域,然后进行连分式展开得到z的整数次幂方程:
其中,Pp,r(z-1)和Qq,r(z-1)分别是关于奇数次p和q的多项式,参数r=1-αm,其具体表达式如表1所示;
表1 Pp,r(z-1)和Qq,r(z-1)的具体表达式
由此,式(1)变换得到z的整数次幂方程:
其中,En-w为时间步数(n-w)处的E,Sm,n-w为时间步数(n-w)处的Sm,n-w,Sm,n-u为时间步数(n-u)处的Sm,n-u;
将式(20)在t=nΔt处进行逆Z变换,得到时域电场强度E和时域电位移矢量D在时间步数n处的时域本构方程:
为方便表述,参数A和参数V为:
将式(24)、式(25)带入式(23),整理得到:
En=A-1·(Dn-Vn) (26)
式(26)即基于Cole-Cole模型的时域本构方程。
其次,结合Cole-Cole模式的DGTD方案。
S2-1:使用DH形式的修正的安培电路环路定律、EH形式的法拉第电磁感应定律,时域麦克斯韦方程组为:
其中,H为时域磁场强度,Jf为施加的电流密度,JM为施加的磁流密度,μ为磁导率;
S2-2:结合DGTD法,用符合旋度的向量基函数ΨD离散D,同时用符合散度的向量基函数ΦE和ΦH分别离散E和H,用符合散度的向量基函数ΦE离散参数V、Sm和 I,在子域k进行离散得到子域k的Maxwell方程和本构方程的弱形式为:
其中,Dk为子域k中的D,Ek为子域k中的E,Hk为子域k中的H,Et为子域 k及其所有邻域交界面处的E,Ht为子域k及其所有邻域交界面处的H,为子域k 中的Jf,为子域k中的JM,μk为子域k中的μ,Ak为子域k中的A,Vk为子域 k中的V,为子域k中的ΨD,为子域k中的ΦE,为子域k中的ΦH,为k 子域边界上向外的单位法向量;
相邻子域用迎风通量耦合,子域k的麦克斯韦方程组和本构方程的半离散系统写为:
其中,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;
第三步,时间步进方案设计。
采用蛙跳时间步进方案,在t=nΔt处对法拉第电磁感应定律的半离散方程(33)进行中心差分,在t=(n+1/2)Δt处对修正的安培环路定律的半离散方程(32)进行中心差分,得到线性方程组:
其中,和分别为在时间步数(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);
整理得到时间步进顺序:
其中,sm是用ΦE离散参数Sm得到的离散的矢量参数,i是用ΦE离散参数I得到的离散的矢量参数,和分别为k域中在时间步数(n+1)、n和(n-1)处的sm,和分别为k域中在时间步数(n+1)和n处的i。
在一些实施方式中,当步骤S1中的方法为双线性变换法时,步骤S3中参数具体为:采用蛙跳时间步进方案,得到时间步进顺序如下:
下面通过具体实施例,验证本发明方法的有效性。
实施例1
为了验证本发明提出的基于Cole-Cole模型的DGTD算法的有效性,选用基于欧拉定律的DGTD方法,用Matlab编程,进行数值仿真,主要仿真过程如下:
(1)参数设置
选择一个尺寸为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。
平面波的传播方向为[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。
(2)结果分析
为了验证所提出的DGTD方案的计算精度,将计算结果转换到频域,并与商用电磁软件FEKO进行比较,如图3(a)-(b)所示。
结果表明,DGTD算法的结果与FEKO的结果非常吻合,验证了针对Cole-Cole 模型的基于欧律定律的DGTD方案仿真的准确性,为后续用DGTD计算分数阶介质、解决FDTD问题做铺垫。
实施例2
为了验证本发明提出的基于Cole-Cole模型的DGTD算法的有效性,选用的是基于双线性变换的DGTD方法,用Matlab编程,进行数值仿真,主要仿真过程如下:
(1)参数设置
构建简单乳腺模型。选择一个半径为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。
平面波的传播方向为[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。
(2)结果分析
为了验证本发明的DGTD方案的计算精度,将计算结果转换到频域,并与商用电磁软件FEKO进行比较,观察点1、观察点2、观察点3、观察点4的结果如图5 (a)-(d)所示。在频率f=1、2、3、4、5、6GHz处,接收面上DGTD方案与FEKO结果的相对误差如图6(a)-(f)所示。
结果表明,DGTD算法的结果与FEKO的结果非常吻合,验证了针对Cole-Cole 模型的基于双线性变换法的DGTD方案仿真的准确性。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (5)
1.一种基于Cole-Cole模型的DGTD电磁瞬态仿真方法,其特征在于,包括以下步骤:
S1:构建基于Cole-Cole模型的时域本构方程;
基于Z变换法,应用欧拉定理及二阶泰勒级数展开,或双线性变换法及连分式展开,将频域Cole-Cole模型转换到Z域,并将分数次幂项转换为整数次幂,然后基于逆Z变换法,将Z域Cole-Cole模型带入电场强度和电位移矢量的Z域本构方程,在时间t=nΔt处将其转换到时域,得到电场强度和电位移矢量的时域本构方程,Δt为时间步长;
S2:结合Cole-Cole模型的EH-DH混合模式的DGTD方案;
基于DGTD法,使用DH形式的修正的安培环路定律、EH形式的法拉第电磁感应定律以及步骤S1中的本构方程,结合迎风通量,构造半离散方程组;
S3:设计时间步进方案;
基于蛙跳法,在t=nΔt处对法拉第电磁感应定律的半离散方程进行中心差分,在t=(n+1/2)Δt处对修正的安培环路定律的半离散方程进行中心差分,得到磁场强度、电位移矢量、电场强度的时间步进方程。
2.根据权利要求1所述的一种基于Cole-Cole模型的DGTD电磁瞬态仿真方法,其特征在于,所述步骤S1的具体方法为:
S1-1:Cole-Cole模型模拟色散介质的复相对介电常数为:
其中,M为极数,ε0、ε∞、Δεm分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度,j为虚数符号,j2=-1,ω为角频率,τm为第m极的弛豫时间值,αm为第m极的色散参数,0≤αm<1,σ为静态电导率;
S1-2:基于欧拉定理,将jω变换到Z域:
由此,式(1)变换得到z的整数次幂方程:
则式(5)表示为:
将式(6)、式(7)在t=nΔt处进行逆Z变换,n为时间步数,得到时域参数Sm和I的时间步进方程:
In=σΔt·En+In-1 (9)
其中,In为参数I在时间步数n处的定义,In-1为参数I在时间步数(n-1)处的定义,Sm,n为Sm在时间步数n处的定义,Sm,n-1为Sm在时间步数(n-1)处的定义,Sm,n-2为Sm在时间步数(n-2)处的定义;
将式(8)在t=nΔt处进行逆Z变换,得到时域电场强度E和时域电位移矢量D在时间步数n处的时域本构方程:
其中,En为在时间步数n处的时域电场强度,Dn为在时间步数n处的时域电位移矢量;
为方便表述,定义参数A和参数V为:
其中,Vn为在时间步数n处的参数V,将式(12)、式(13)带入式(11),整理得到:
En=A-1·(Dn-Vn) (14)
式(14)即基于Cole-Cole模型的时域本构方程。
3.根据权利要求1所述的一种基于Cole-Cole模型的DGTD仿真方法,其特征在于,所述步骤S2的具体方法为:
S2-1:使用DH形式的修正的安培电路环路定律、EH形式的法拉第电磁感应定律,时域麦克斯韦方程组为:
其中,H为时域磁场强度,Jf为施加的电流密度,JM为施加的磁流密度,μ为磁导率;
S2-2:结合DGTD法,用符合旋度的向量基函数ΨD离散D,同时用符合散度的向量基函数ΦE和ΦH分别离散E和H,用符合散度的向量基函数ΦE离散参数V、Sm和I,在子域k进行离散得到子域k的Maxwell方程和本构方程的弱形式为:
其中,Dk为子域k中的D,Ek为子域k中的E,Hk为子域k中的H,Et为子域k及其所有邻域交界面处的E,Ht为子域k及其所有邻域交界面处的H,为子域k中的Jf,为子域k中的JM,μk为子域k中的μ,Ak为子域k中的A,Vk为子域k中的V,为子域k中的ΨD,为子域k中的ΦE,为子域k中的ΦH,为k子域边界上向外的单位法向量;
相邻子域用迎风通量耦合,子域k的麦克斯韦方程组和本构方程的半离散系统写为:
其中,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;
4.根据权利要求1所述的一种基于Cole-Cole模型的弱形式电磁瞬态仿真方法,其特征在于,所述步骤S3中的具体方法为:
采用蛙跳时间步进方案,在t=nΔt处对法拉第电磁感应定律的半离散方程(21)进行中心差分,在t=(n+1/2)Δt处对修正的安培环路定律的半离散方程(20)进行中心差分,得到线性方程组:
其中,和分别为在时间步数(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);
整理得到时间步进顺序:
5.根据权利要求1所述的一种基于Cole-Cole模型的DGTD电磁瞬态仿真方法,其特征在于,所述步骤S1中为双线性变换法,则:
Step1-1:Cole-Cole模型模拟色散介质的复相对介电常数为:
其中,M为极数,ε0、ε∞、Δεm分别为复相对介电常数、自由空间介电常数、无限频率相对介电常数和第m极的幅度,j为虚数符号,j2=-1,ω为角频率,τm为第m极的弛豫时间值,αm为第m极的色散参数,0≤αm<1,σ为静态电导率;
Step1-2:基于双线性变换法,将jω的分数次幂表达式变换到Z域,然后进行连分式展开得到z的整数次幂方程:
其中,Pp,r(z-1)和Qq,r(z-1)分别是关于奇数次p和q的多项式,参数r=1-αm;
由此,式(1)变换得到z的整数次幂方程:
其中,En-w为时间步数(n-w)处的E,Sm,n-w为时间步数(n-w)处的Sm,n-w,Sm,n-u为时间步数(n-u)处的Sm,n-u;
将式(37)在t=nΔt处进行逆Z变换,得到时域电场强度E和时域电位移矢量D在时间步数n处的时域本构方程:
为方便表述,参数A和参数V为:
将式(41)、式(42)带入式(40),整理得到:
En=A-1·(Dn-Vn) (43)
式(43)即基于Cole-Cole模型的时域本构方程。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110709190.4A CN113536626B (zh) | 2021-06-25 | 2021-06-25 | 基于Cole-Cole模型的DGTD电磁瞬态仿真方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110709190.4A CN113536626B (zh) | 2021-06-25 | 2021-06-25 | 基于Cole-Cole模型的DGTD电磁瞬态仿真方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113536626A true CN113536626A (zh) | 2021-10-22 |
CN113536626B CN113536626B (zh) | 2023-05-26 |
Family
ID=78096798
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110709190.4A Active CN113536626B (zh) | 2021-06-25 | 2021-06-25 | 基于Cole-Cole模型的DGTD电磁瞬态仿真方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113536626B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116861711A (zh) * | 2023-09-05 | 2023-10-10 | 北京航空航天大学 | 一种基于时域间断伽辽金方法的双各向同性媒质仿真方法 |
CN116882219A (zh) * | 2023-09-07 | 2023-10-13 | 北京航空航天大学 | 一种基于dgtd与fdtd的场线耦合算法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2016529923A (ja) * | 2013-05-06 | 2016-09-29 | パルティ、ヨーラム | 交番電界によって腫瘍を治療し、推定される細胞サイズに基づいて治療周波数を選択するための装置および方法 |
WO2018058720A1 (zh) * | 2016-09-28 | 2018-04-05 | 深圳先进技术研究院 | 一种建立血液电磁仿真模型的方法及装置 |
CN107992696A (zh) * | 2017-12-13 | 2018-05-04 | 电子科技大学 | 一种复杂色散媒质中的改进的指数时间积分构造方法 |
-
2021
- 2021-06-25 CN CN202110709190.4A patent/CN113536626B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2016529923A (ja) * | 2013-05-06 | 2016-09-29 | パルティ、ヨーラム | 交番電界によって腫瘍を治療し、推定される細胞サイズに基づいて治療周波数を選択するための装置および方法 |
WO2018058720A1 (zh) * | 2016-09-28 | 2018-04-05 | 深圳先进技术研究院 | 一种建立血液电磁仿真模型的方法及装置 |
CN107992696A (zh) * | 2017-12-13 | 2018-05-04 | 电子科技大学 | 一种复杂色散媒质中的改进的指数时间积分构造方法 |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116861711A (zh) * | 2023-09-05 | 2023-10-10 | 北京航空航天大学 | 一种基于时域间断伽辽金方法的双各向同性媒质仿真方法 |
CN116861711B (zh) * | 2023-09-05 | 2023-11-17 | 北京航空航天大学 | 一种基于时域间断伽辽金方法的双各向同性媒质仿真方法 |
CN116882219A (zh) * | 2023-09-07 | 2023-10-13 | 北京航空航天大学 | 一种基于dgtd与fdtd的场线耦合算法 |
CN116882219B (zh) * | 2023-09-07 | 2023-11-14 | 北京航空航天大学 | 一种基于dgtd与fdtd的场线耦合方法 |
Also Published As
Publication number | Publication date |
---|---|
CN113536626B (zh) | 2023-05-26 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113536626B (zh) | 基于Cole-Cole模型的DGTD电磁瞬态仿真方法 | |
Lu et al. | Acceleration of frequency sweeping in eddy-current computation | |
Ilic et al. | Higher order hybrid FEM-MoM technique for analysis of antennas and scatterers | |
Zygiridis et al. | Low-dispersion algorithms based on the higher order (2, 4) FDTD method | |
Cai et al. | Volume surface integral equation method based on higher order hierarchical vector basis functions for EM scattering and radiation from composite metallic and dielectric structures | |
Abalenkovs et al. | Huygens subgridding for 3-D frequency-dependent finite-difference time-domain method | |
CN108021533A (zh) | 一种基于广义坐标系求解任意色散材料电磁特性的方法 | |
CN108052738A (zh) | 色散媒质的高阶局部无条件稳定时域间断伽辽金分析方法 | |
Guo et al. | An edge-based smoothed finite element method for the assessment of human exposure to extremely low frequency electric fields | |
Chen et al. | A nonconformal surface integral equation for electromagnetic scattering by multiscale conducting objects | |
CN113486294B (zh) | 一种处理复杂色散介质的无条件稳定fdtd算法 | |
Young et al. | A locally corrected Nyström formulation for the magnetostatic volume integral equation | |
CN105808504A (zh) | 一种等离子体中使用辅助微分方程完全匹配层的实现方法 | |
Zhou et al. | An adaptive DGTD algorithm based on hierarchical vector basis function | |
CN105277927A (zh) | 飞行器编队瞬态电磁特性时域阶数步进分析方法 | |
Cheng et al. | A stable FDTD subgridding scheme with SBP-SAT for transient electromagnetic analysis | |
Wu et al. | DEH scheme DGTD-based transient modeling approach for the Cole–Cole dispersive media using Tustin’s method | |
CN117787056B (zh) | 一种适用于Debye模型电磁场传播的高精度计算方法 | |
Kim et al. | Generalized Transition Matrix Model Using Characteristic Basis Function Method for Open-Ended Cavities | |
Zhao et al. | A three-dimensional Petrov-Galerkin finite element interface method for solving inhomogeneous anisotropic Maxwell's equations in irregular regions | |
Kireeva et al. | A coupled EFGM–CIE method for acoustic radiation | |
Zhao et al. | Discontinuous Galerkin time domain modeling of metasurface geometries with multi-rate time stepping | |
Bushyager et al. | Modelling of complex RF/wireless structures using computationally optimized time‐domain techniques | |
Camps-Raga et al. | Optimized simulation algorithms for fractal generation and analysis | |
Yang et al. | SO-DGTD method for dispersive media |
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 |