CN113486294B - 一种处理复杂色散介质的无条件稳定fdtd算法 - Google Patents

Abstract

本发明公开了一种处理复杂色散介质的无条件稳定FDTD算法，属于计算电磁学领域。本发明方法基于laguerre多项式和修正洛伦兹模型的无条件稳定WLP‑FDTD算法能够处理Debye、Drude、Lorentz、critical point、QCRF、CCPR和Cole‑Cole模型及其任意线性组合形式表述的复杂色散介质目标。具有更加通用、更加高效、无条件稳定的优势。

Description

一种处理复杂色散介质的无条件稳定FDTD算法

技术领域

本发明属于计算电磁学领域，具体涉及一种处理复杂色散介质目标的基于修正洛伦兹模型、加权laguerre多项式的无条件稳定FDTD算法。

背景技术

计算电磁学是结合电子计算机和数值计算方法求解电磁学中麦克斯韦(Maxwell)方程的一门交叉学科。随着计算机技术的发展和电磁理论的不断丰富，时域电磁计算方法取得了长足发展，与频域数值方法相辅相成。时域的方法使得人们可以直接对宽频谱特性的时变电磁现象进行分析，从而可以更直观地模拟时变电磁场与物体之间的相互作用。时域有限差分法是电磁场时域技术中最具代表性、发展最快的数值计算方法，是求解电磁问题的常用数值方法之一。传统的FDTD由于简单通用、易于实施、具有可移植性，基于这些好的特性，被广泛应用于通信、雷达、电磁防护、电磁兼容、电磁隐身、医疗诊断等专业领域。传统的FDTD方法具有简单实用的特点，但在计算精细结构或复杂结构时会在时间域受到CFL(Courant-Friedrich-Levy)时间稳定性条件的限制，在空间域还会受到离散带来的数值色散误差的影响。如果不满足CFL时间稳定性条件的限制，随着时间的步进的增加，累计的误差会急剧增大，使得无法得到准确解。在计算精细结构和复杂结构的时候所需的最小网格往往很小，为保证计算精度，采用的时间步长也必须很小，这就大大影响FDTD的计算效率。

为了提高计算效率，许多无条件稳定的时域数值计算方法相继提出，如交变隐式差分方向(alternating-direction implicit,ADI)FDTD，分步(split step,SS)FDTD，Crank-Nicolson(CN)FDTD，局部一维(locally one-dimensional,LOD)FDTD等。然而这些方法虽然提高了计算效率，但随着时间步长的增加，越来越大的数值色散误差会影响计算结果的准确性。基于加权Laguerre多项式(weigthed Laguerre polynomials,WLPs)的FDTD算法作为一种无条件稳定的时域快速算法，在最近几年引起人们的广泛注意。Y.S.Chung和T.K.Sarkar将加权Laguerre多项式作为时域基函数，采用Garlerkin方法对Maxwell方程中的时间偏微分进行解析处理，消除了时间微分项，得到了一种无条件稳定的FDTD算法—WLP-FDTD算法。该算法的计算格式摆脱了CFL时间稳定性条件的限制，并且拥有良好的数值色散特性和计算精度，在解决宽频带和多尺度复杂结构的电磁特性问题中具有比较高的计算精度和计算效率。目前，基于WLP-FDTD的算法多用于处理洛伦兹(Lorentz)、库德(Drude)、德拜(Debye)模型，尚没有报道能处理复杂的色散介质模型，如CCPR(complex-conjugate pole-residue)、QCPF(quadratic complex rational function)、以及cole-cole小数模型等。

目前用于处理色散介质的FDTD方法有很多，如递推卷积法(RC)FDTD、辅助方程法(ADE-FDTD)、Z变换的FDTD法、移位算子(SO)FDTD法等、最近有学者Kyung-Young Jung等报道称一种基于修正洛伦兹模型的FDTD方法可以覆盖Debye、Drude、Lorentz、criticalpoint、QCRF、CCPR等复杂色散模型，具有更好的通用性和鲁棒性。色散煤质指的是一类介电系数或者磁导率等电磁特性参数随着频率变化而变化的介质，诸如水、生物体组织、土壤、植被、等离子体、雷达吸波材料等都属于色散介质。色散介质在雷达技术、微波遥感、空间物理、生命科学、目标隐身等领域具有诸多应用，因此研究电磁波与色散介质的相互作用特性具有重要用途，如雷达目标散射、电磁探测、医学诊断和宇航通信等，都需要考虑目标的色散特性，否则会导致计算失效。由于色散煤质的电磁特性复杂多变，在处理电磁散射问题上，相比于解析方法，数值方法是比较好的选择。在众多的数值方法中FDTD算法的优点使得其成为首选。由于在自然界中大多数电磁材料的电磁参数与频率具有任意关系，在通常的计算中通常采用Debye、Drude、Lorentz等几种色散模型的线性组合来拟合复杂电磁材料。Debye模型常用于水、生物组织的色散表述，Drude模型常用于表述等离子体、金属介质的色散关系，Lorentz模型常用于表述光学材料、生物组织的色散关系，其他复杂的色散模型如critical point常用于金属材料在光学频段的电磁特性，QCPF常用于表述人体组织、聚合物等介质，CCPR常用于表述生物组织、金属等介质，分数阶导数的Cole-Cole模型常用于表述生物组织。传统的无条件稳定的FDTD方法对于不同的介质模型需要不同的公式推导，通用性较差。虽然已经有一些通用的方法被提出，但是对于比较复杂的色散模型，仍不能较好地覆盖。报道称，修正的洛伦兹模型可以覆盖Debye、Drude、Lorentz、critical point、QCRF、CCPR模型，具有更好的通用性和鲁棒性，但此方法并未与无条件稳定的FDTD结合，且不能覆盖分数阶导数模型，如Cole-Cole模型。

因此，运用基于修正洛伦兹模型的无条件稳定的WLP-FDTD方法，模拟计算能够用Debye、Drude、Lorentz、critical point、QCRF、CCPR、Cole-Cole模型及其任意线性组合形式表述的复杂色散介质目标的波传播特性，具有更高效的计算优势。

发明内容

针对现有技术中的不足，本发明提供了一种处理复杂色散介质目标的无条件稳定FDTD的方法，将基于修正的洛伦兹模型的WLP-FDTD算法应用于能够用Debye、Drude、Lorentz、critical point、QCRF、CCPR和Cole-Cole模型及其任意线性组合形式表述的复杂色散介质目标中，如生物组织、等离子体、金属、土壤。本发明方法具有更加通用、更加高效、无条件稳定的优势。

为了达到上述发明目的，本发明采用的技术方案如下：

一种处理复杂色散介质目标的无条件稳定FDTD的方法，包括如下步骤：

步骤1：构建色散介质目标的物理模型，并选取色散模型(Debye、Drude、Lorentz、critical point、QCRF、CCPR和Cole-Cole模型，如模拟生物组织的色散特性可选取Cole-Cole模型)，将选取的色散模型转换为统一格式的修正洛伦兹模型。

所述统一格式的修正洛伦兹模型如式(1)：

其中ε_∞为无穷介电系数，ω表示角频率，m表示修正洛伦兹模型和项的阶数，M表示修正洛伦兹模型和项的总阶数，a_0,m、a_1,m、b_0,m、b_1,m、b_2,m均为修正洛伦兹模型的待定系数，具体的取值由选取的色散模型决定，j表示虚部单位。

步骤2：引入辅助差分方程(2)代入麦克斯韦方程组中，得到如式(3)(4)(5)所示的含有辅助差分方程的麦克斯韦方程组。

所述辅助差分方程如式(2)：

其中E(ω)表示频域的电场强度,S_m(ω)表示引入修正洛伦兹模型后的辅助差分方程。

所述含有辅助差分方程的麦克斯韦方程组如式(3)(4)(5)：

其中，r表示位置向量，t表示时间变量，ε₀是自由空间的介电系数，μ₀是自由空间的导磁率，S_m(r,t)，H(r,t)，J(r,t)，E(r,t)分别表示时域的辅助差分方程、磁场强度，电流密度，电场强度。

步骤3：在含有辅助差分方程的麦克斯韦方程组中引入加权laguerre多项式，利用加权Laguerre多项式作为时域基函数展开式(3)(4)(5)中的电场分量、磁场分量，然后将展开后的电场分量、磁场分量再代入含有辅助差分方程的麦克斯韦方程中，采用伽辽金方法，并利用加权Laguerre多项式的正交性消除公式(3)(4)(5)中的时间变量，化简后得到Laguerre域的麦克斯韦方程组。

所述加权laguerre多项式

如式(6)：

其中s为时间尺度因子，引入s调节收敛速度，p为加权laguerre多项式的阶数，L_p(t)是阶数为p的Laguerre多项式，L_p(st)为引入时间尺度因子s后的Laguerre多项式。

所述Laguerre域的麦克斯韦方程组如式(8)(9)(10)：

其中q表示加权laguerre多项式当下的阶数，k表示加权laguerre多项式的历史阶数，A_m＝b_0,m+0.5sb_1,m+0.25b_2,ms²，为计算的过程变量无具体含义。S_m(r)，H(r)，J(r)，E(r)分别表示消除时间变量后的辅助差分方程、磁场强度，电流密度，电场强度。

步骤4：选取Mur吸收边界条件或者PML吸收边界条件来截断计算区域；加入激励源，如调制高斯脉冲源来模拟FDTD迭代过程场的源。

步骤5：求解公式(8)(9)(10)，得到目标求解区域的电磁场，最终获得色散介质目标波的传播特性。

本发明所提出的基于laguerre多项式和修正洛伦兹模型的无条件稳定WLP-FDTD算法能够处理用Debye、Drude、Lorentz、critical point、QCRF、CCPR和Cole-Cole模型及其任意线性组合形式表述的复杂色散介质目标。

本方案的有益效果为：

(1)WLP-FDTD算法的计算格式摆脱了CFL时间稳定性条件的限制，并且拥有良好的数值色散特性和计算精度，在解决宽频带和多尺度复杂结构的电磁特性问题中具有比较高的计算精度和计算效率。

(2)本发明所述方法能够覆盖Debye、Drude、Lorentz、critical point、QCRF、CCPR和Cole-Cole模型及其任意线性组合形式的复杂色散介质目标，修正洛伦兹模型的极数M可任意拓展，与无条件稳定的WLP-FDTD结合，拓展了WLP-FDTD的应用范围，形成一套通用的具有统一格式的WLP-FDTD算法体系，同时也拓展了修正洛伦兹模型的覆盖能力以应用于其他FDTD算法。

附图说明

图1为基于修正洛伦兹模型的无条件稳定的ADE-WLP-FDTD算法流程图。

图2所述为4-th cole-cole模型的相对介电常数(表2所示)和转化为修正洛伦兹模型后的相对介电常数(表3所示)后二者对比图，电导率σ＝0.2，无穷介电系数ε_∞＝4，真空介电系数ε₀＝8.8541878×10^-12，由图2可见二者拟合效果非常好。

图3所示为采用本文所述基于修正洛伦兹模型的WLP-FDTD算法计算电磁波通过肌肉组织的透射系数，采用表3所示转换后的修正洛伦兹模型参数。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

本发明实施例中，提供了一种处理复杂色散介质目标的无条件稳定FDTD算法，还提出了一种将具有分数阶导数的Cole-Cole模型转换为修正洛伦兹模型的方法。

本实施例的算法流程图如图1所示，采用4阶cole-cole模型研究一维问题下的肌肉组织的波传播特性，同时，还提供了本发明算法的推导过程，包括如下步骤：

步骤1：

麦克斯韦旋度方程组为：

其中，μ₀表示自由空间的磁导率，电位移矢量D(ω)与电场强度E(ω)在频域的关系即本构关系为：

D(ω)＝ε₀ε_r(ω)E(ω)                          (12)

其中ε₀是自由空间的介电系数。相对介电系数ε_r根据介质目标采用不同的色散模型。统一格式的修正洛伦兹模型表示为如下式(1)：

其中ε_r,∞是无穷介电系数，ω表示角频率，m表示修正洛伦兹模型和项的阶数，a_0,m,a_1,m,b_0,m,b_1,m,b_0,m为修正洛伦兹模型的待定系数。表1给出了其他色散模型转化为修正洛伦兹模型的对应关系。表1中所有参数均为针对各个介质模型的待定系数，其中Debye、Drude、Lorentz、critical point、QCRF、CCPR转换为修正洛伦兹模型的过程及转换关系为本领域人员的公知常识，这里不再赘述。

本发明还提供了cole-cole模型转换为修正洛伦兹模型的方法：根据帕德近似原理对β(0≤β≤1)的分数阶导数进行转换，将Cole-Cole模型转换为QCRF模型，再将QCRF模型转换为统一格式的修正洛伦兹模型，即实现将cole-cole模型转化为修正洛伦兹模型。具体地，本实施例表1中还提供了cole-cole模型转化为修正洛伦兹模型的转换关系。

对于表1中cole-cole模型向修正洛伦兹模型的转换关系，其中的系数A₀，A₁，A₂，B₀，B₁，B₂由公式(13)得到：

其中，h＝0,1,2，f＝0,1,2，g＝0,1,2，u＝0,1,2，β表征离散范围，τ表示弛豫时间，ω₀表示特征角频率，a_g,g-h＝(g！)/(h！(g-h)！)(-1)^g-h，a_u,u-f＝(u！)/(f！(u-f)！)(-1)^u-f。p_u，q_g为计算过程中的系数，无实际意义。具体地，当g＝0时,q₀＝1；当g＝1,2时，通过式(14)计算得到q₁，q₂：

将q₀，q₁，q₂代入公式(15)：

得到p₀，p₁，p₂。

采用4阶cole-cole模型研究1-D问题下的肌肉组织的波传播特性。肌肉组织的4阶cole-cole模型的参数如表2，转换为修正洛伦兹模型后的参数如表3所示。另肌肉组织的电导率σ＝0.2，ε_∞＝4，ε₀＝8.8541878×10^-12。

表1：其它色散模型转化为修正洛伦兹模型的转换关系

其中ε_r,∞表示无穷介电系数，j表示虚部单位，ω表示角频率，其余参数表示对应模型的待定系数。对于cole-cole模型，Δε为极点振幅。

表2肌肉组织的Cole-Cole模型参数

表3肌肉组织的转化后的修正洛伦兹模型参数

步骤2：

Laguerre多项式具有如下定义形式：

L_p(t)是阶数为p的Laguerre多项式。由于Laguerre多项式只有t≥0时才有意义，因此laguerre多项式具有因果关系。Laguerre多项式具有如下递推关系：

L₀(t)＝1

L₁(t)＝1-t

pL_p(t)＝(2p-1-t)L_p-1(t)-(p-1)L_p-2(t)，p≥2,t≥0

Laguerre多项式关于加权函数e^-t的正交性是：

上式中δ_pq是Dirac函数，因此可以得出加权laguerre多项式的正交基函数系

引入时间尺度因子s调节收敛速度，

表示阶数为p的加权laguerre多项式。

引入辅助差分方程：

(2)式可写成：

根据从频域转换为时域的转换关系

公式(19)(20)在频域中电位移矢量D(ω),电场强度E(ω),辅助差分方程S(ω)之间的关系转化到时域的对应关系为：

将(21)代入公式(11)中可写为：

综上(3)(4)(5)构成含有辅助差分方程的Maxwell方程组。

步骤3：

利用加权Laguerre正交基函数

电场、磁场以及辅助差分方程可表示为：

其中，s是时间尺度因子。将式(22)所示的电场、磁场和辅助差分方程的Laguerre展开形式代入到式(3)(4)(5)所示的含有辅助差分方程的Maxwell方程组中，随后，为了消除计算过程的时间变量，此处利用Galerkin方法，将方程组两边同时乘以权函数

利用加权Laguerre多项式的正交性，在区间[0，∞)对时间积分，消除时间变量t，化简得到：

其中，A_m＝b_0,m+0.5sb_1,m+0.25b_2,ms²，

表示消除时间变量后的电流源，T_f为激励源时域波形的持续时间。

将(9)代入(23)中，得到

公式(8)(9)(10)组成新的用于求解电磁场的laguerre域的Maxwell方程组。

步骤4：

边界条件采用Mur一阶吸收边界来截断WLP-FDTD方法的计算区域。在x＝0(计算区域开始的地方)或者x_max处(计算区域结束的地方)的Mur一阶吸收边界条件是：

v₁表示电磁波在介质中的传播速度。将式(14)代入式(20)，在x_max边界处利用伽辽金方法，两边同时乘以权函数

并在区间[0，∞)对时间积分，可以得到：

在x＝0处的吸收边界条件与x＝x_max处的计算方法相同。

激励源采用调制高斯脉冲，时域表达式如下式

其中，f_c表示调制高斯脉冲的中心频率，T_c决定调制高斯脉冲的相位，T_d决定调制高斯脉冲的宽度，T_d＝1/2f_c,T_c＝3T_d。

步骤5：

针对一维的肌肉组织模型，本实施例给出所述laguerre域Maxwell方程组一维离散形式如下。

其中i表示一维离散后的位置，Δz表示离散的空间步长，c₁,c₂为计算过程变量无具体含义，

式(26)(27)(28)即为一维离散后，最后求解一维色散介质目标波传播特性的Maxwell方程组。二维或者三维色散目标需要对式(19)(21)(22)作二维或者三维离散，计算过程与一维的计算过程相同。

图3即为求解方程组式(26)(27)(28)的一维的肌肉组织模型透射系数。

Claims (3)

1.一种处理复杂色散介质目标的无条件稳定FDTD的方法，其特征在于，包括如下步骤：

步骤1：构建色散介质目标的物理模型，并选取色散模型，将选取的色散模型转换为统一格式的修正洛伦兹模型；

所述统一格式的修正洛伦兹模型如式(1)：

其中ε_∞为无穷介电系数，ω表示角频率，m表示修正洛伦兹模型和项的阶数，M表示修正洛伦兹模型和项的总阶数，a_0,m、a_1,m、b_0,m、b_1,m、b_2,m均为修正洛伦兹模型的待定系数，具体的取值由选取的色散模型决定，j表示虚部单位；

步骤2：引入辅助差分方程(2)代入麦克斯韦方程组中，得到如式(3)(4)(5)所示的含有辅助差分方程的麦克斯韦方程组；

所述辅助差分方程如式(2)：

其中E(ω)表示频域的电场强度,S_m(ω)表示引入修正洛伦兹模型后的辅助差分方程；

所述含有辅助差分方程的麦克斯韦方程组如式(3)(4)(5)：

其中，r表示位置向量，t表示时间变量，ε₀是自由空间的介电系数，μ₀是自由空间的导磁率，S_m(r,t)，H(r,t)，J(r,t)，E(r,t)分别表示时域的辅助差分方程、磁场强度，电流密度，电场强度；

步骤3：在含有辅助差分方程的麦克斯韦方程组中引入加权laguerre多项式，利用加权Laguerre多项式作为时域基函数展开式(3)(4)(5)中的电场分量、磁场分量，然后将展开后的电场分量、磁场分量再代入含有辅助差分方程的麦克斯韦方程中，采用伽辽金方法，并利用加权Laguerre多项式的正交性消除公式(3)(4)(5)中的时间变量，化简后得到Laguerre域的麦克斯韦方程组；

所述加权laguerre多项式

如式(6)：

其中s为时间尺度因子，引入s调节收敛速度，p为加权laguerre多项式的阶数，L_p(t)是阶数为p的Laguerre多项式，L_p(st)为引入时间尺度因子s后的Laguerre多项式；

所述Laguerre域的麦克斯韦方程组如式(8)(9)(10)：

其中q表示加权laguerre多项式当下的阶数，k表示加权laguerre多项式的历史阶数，A_m＝b_0,m+0.5sb_1,m+0.25b_2,ms²，为计算的过程变量无具体含义；S_m(r)，H(r)，J(r)，E(r)分别表示消除时间变量后的辅助差分方程、磁场强度，电流密度，电场强度；

步骤4：选取Mur吸收边界条件或者PML吸收边界条件来截断计算区域；加入激励源；

步骤5：求解公式(8)(9)(10)，得到目标求解区域的电磁场，最终获得色散介质目标波的传播特性。

2.如权利要求1所述的一种处理复杂色散介质目标的无条件稳定FDTD的方法，其特征在于，所述色散模型包括Debye、Drude、Lorentz、critical point、QCRF、CCPR和Cole-Cole模型。

3.如权利要求1所述的一种处理复杂色散介质目标的无条件稳定FDTD的方法，其特征在于，cole-cole模型转换为修正洛伦兹模型的方法为：

根据帕德近似原理对β的分数阶导数进行转换，其中0≤β≤1；将Cole-Cole模型转换为QCRF模型，再将QCRF模型转换为统一格式的修正洛伦兹模型，实现将cole-cole模型转化为修正洛伦兹模型。

Images