CN111783339B - 电磁波在随机色散介质中传播的pce-fdtd方法 - Google Patents

Abstract

本发明提供了电磁波在随机色散介质中传播的PCE‑FDTD方法，包括以下步骤：输入模型文件；初始化参数；计算仿真中所用到的内积；添加场源到中设置好源的位置；更新计算整个计算区域的y方向上电场分量展开系数更新计算整个计算区域的x方向上电场分量展开系数更新计算整个计算区域的z方向上磁场分量展开系数将时间迭代步t+1赋值给t，并判断迭代次数t是否达到预设值，若未达到预设值，则返回；若达到预设值，则执行下一步；计算观测点电场的均值和标准差。该方法避免了相关系数估计所产生的误差；与PCE‑JEC‑FDTD方法相比不仅可以应用于随机等离子体，还可以应用于其他随机色散媒质，不需要重新推导公式。

Description

电磁波在随机色散介质中传播的PCE-FDTD方法

技术领域

本发明属于计算电磁学技术领域，具体涉及一种电磁波在随机色散介质中传播的PCE-FDTD方法。

背景技术

众所周知，时域有限差分(Finite-difference time-domain，FDTD)方法是一种简单有效的电磁计算方法，传统的FDTD方法计算过程中使用的是模型中媒质电参数的均值，但是自然界中许多媒质的电参数具有随机特性，为了表征随机媒质中的电磁特性，人们提出了应用于随机媒质中的时域有限差分方法，比如：门特卡罗(Monte Carlo，MC)方法、随机FDTD(stochastic FDTD，S-FDTD)和多项式混沌展开FDTD(polynomial chaotic expansionFDTD，PCE-FDTD)方法。在这些方法中，MC方法简单而且精度高，但是这种方法需要成千上万次的仿真，太耗费时间。之前提出的S-FDTD只需在FDTD算法中加入一些附加的方差方程，即可直接估计出各点电磁场的均值和方差，解决了MC方法耗费时间的问题，但它在计算过程中会由于相关系数估计引入误差。为了解决这种问题，有人提出了PCE-FDTD。PCE-FDTD是将电磁场值进行多项式展开，通过计算出的各阶系数计算场值的均值和方差，其在解决电磁场问题时具有一定的优越性。然而截至目前，所提出的采用电流(J)和电场(E)并置的FDTD(PCE-JEC-FDTD)方法表征了电磁波在磁化等离子体中传播的不确定性，PCE方法仅仅被用于随机等离子体的电磁仿真，对于其他类型的随机色散媒质(如洛伦兹媒质和德拜媒质)就得重新推导公式。

需要注意的是，本部分旨在为权利要求书中陈述的本发明的实施方式提供背景或上下文。此处的描述不因为包括在本部分中就承认是现有技术。

发明内容

本发明目的在于提供了电磁波在随机色散介质中传播的PCE-FDTD方法，能够应用于各种不同的随机色散媒质，计算速度快且精度高。

为实现上述目的本发明采用如下技术方案：

该电磁波在随机色散介质中传播的PCE-FDTD方法，包括以下步骤：

步骤1：输入模型文件；

步骤2：初始化参数；

步骤3：计算仿真中所用到的内积；

步骤4：添加场源到中设置好源的位置；

步骤5：更新计算整个计算区域的y方向上电场分量展开系数

步骤6：更新计算整个计算区域的x方向上电场分量展开系数

步骤7：更新计算整个计算区域的z方向上磁场分量展开系数

步骤8：将时间迭代步t+1赋值给t，并判断迭代次数t是否达到预设值，若未达到预设值，则返回步骤5；若达到预设值，则执行步骤9；

步骤9：计算观测点电场的均值和标准差。

进一步地，上述步骤1具体是：

计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δζ，ζ＝x，y，x为横坐标，y为纵坐标；设置随机色散媒质的区域，时间步长Δt，真空中的磁导率μ₀，介电常数ε₀，仿真的总时长N_t，设置源的位置，观测点的位置。

进一步地，上述步骤2具体是：

将整个计算区域的电磁场分量展开系数中间变量观测点电场的均值和标准差(Ex_80，sigma80)均初始化为0；同时确定展开阶数l的的值。

进一步地，上述步骤3中内积的由下式计算得出

进一步地，上述步骤4中添加场源的表达式为：

进一步地，上述步骤5具体是：

步骤5.1：先计算电位移矢量分量系数其在计算区域的方程为：

其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

步骤5.2：求解5.1的方程得到电位移矢量分量系数后，再计算电场分量系数/>在随机色散媒质区域电场分量系数/>的计算方程为：

在非随机色散媒质区域电场分量系数的计算方程为：

步骤5.3：更新中间变量

进一步地，上述步骤6具体是：

步骤6.1：先计算电位移矢量分量系数其在计算区域的方程为：

步骤6.2：求解6.1的方程得到电位移矢量分量系数后，再计算电场分量系数/>在随机色散媒质区域电场分量系数/>的计算方程为：

在非随机色散媒质区域电场分量系数的计算方程为：

步骤6.3：更新中间变量

进一步地，上述步骤7的计算方程式为：

进一步地，上述步骤9的计算公式为：

本发明的有益效果：

本发明基于双线性变换的时域有限差分方法，对电磁场进行多项式混沌展开，计算得到电磁场的各阶系数后可以很容易地计算出场值的均值和标准偏差。该发明与传统的门特卡罗(MC)方法相比可以大幅提升计算速度；与S-FDTD相比，避免了相关系数估计所产生的误差；与PCE-JEC-FDTD方法相比不仅可以应用于随机等离子体，还可以应用于其他随机色散媒质(如洛伦兹媒质和德拜媒质)，不需要重新推导公式。

附图说明

图1是本发明的流程图；

图2是本发明实施例中的计算模型的结构示意图；

图3是本发明的方法与MC方法计算的电场均值对比图；

图4是本发明的方法与MC方法计算的电场标准偏差对比图。

具体实施方式

现在将参考附图更全面地描述示例实施方式。然而，示例实施方式能够以多种形式实施，且不应被理解为限于在此阐述的范例；相反，提供这些实施方式使得本发明将更加全面和完整，并将示例实施方式的构思全面地传达给本领域的技术人员。所描述的特征或特性可以以任何合适的方式结合在一个或更多实施方式中。

本发明是基于多项式混沌展开的应用于电磁波在随机色散介质中传播的通用FDTD方法，所依据的原理是：首先导出直角坐标系下电磁场所满足的麦克斯韦方程；然后利用BT方法将电场的本构关系从频域转换到时域；再利用PCE展开电磁场分量，求解电磁场分量的各阶展开系数；最后求解观测点处电磁场分量的均值和标准偏差。

由于本发明在程序设计中应用了PML边界条件，在求电磁场分量系数更新方程时，首先需要推导出复扩展直角坐标系下，PML中电磁场满足的麦克斯韦方程；

复扩展直角坐标系下，PML中电磁场满足的麦克斯韦方程为：

其中，表示电场矢量，/>表示磁场矢量，/>表示电位移矢量，j为虚数单位，ω为角频率，μ为介质的磁导率，ε₀为真空中的介电常数，ε_r为介质的相对介电常数，/>为修正后的微分算子，可以写成：

S_x、S_y和S_z是坐标扩展变量，可以表示成：

S_ζ＝k_ζ+σ_ζ/jωε (4)

其中ζ表示x、y、z，k_ζ、σ_ζ和α_ζ为PML的有关的参数。

本发明仅考虑随机色散媒质中二维横电波情况，于是复扩展直角坐标下的麦克斯韦方程可以写成：

其中D_x,D_y分别表示x,y方向的电位移分量，H_z分别表示z方向的磁场。然后，求出电磁场分量系数的更新方程；

式(2)是色散介质的本构关系，其给出了电位移矢量和电场的转换关系，复数相对介电常数可以用有理分数形式表示为

其中p_k和q_k是取决于材料参数的多项式系数。表一列出了三种媒体(Drude媒体，Lorentz媒体和Debye媒体)的系数。

表一:不同类型材料的系数

为了获得离散差分公式，通过双线性变换(BT)将频域中的本构关系方程转换为离散时域，

其中z^-1是时移运算符，即z^-1Fⁿ＝F^n-1。多项式简化后，(2)可以表示为

其中和/>分别是/>和/>在时间步n-k的系数。表二和表三中列出了不同类型材料的系数。

表二:方程(10)中的等离子体和洛伦兹材料的系数

表三：方程(10)中导电材料和德拜材料的系数

为了更好地表征式(10)中材料参数的不确定性，将其表示为随机变量ξ的函数，其可以是高斯分布或均匀分布。这样使得输出场也成为随机变量的函数。为了实现多项式混沌方法，需要将不确定域用正交基函数展开，例如，

其中，是展开系数，a表示展开的阶数，Ψ_a(ξ)是由随机变量ξ的分布函数确定的a阶正交多项式基函数。ξ＝[ξ₁,ξ₂,...,ξ_n]表示包含n个独立随机变量的向量。对于不同定义区域的随机量和不同类型的概率分布采用不同的插值多项式。如高斯分布对应埃尔米特插值多项式，均匀分布对应勒让德插值多项式，伽马γ分布对应拉格朗日插值多项式，贝塔β分布对应雅可比插值多项式等。本发明研究的随机变量为高斯分布，故选择了埃尔米特多项式。表四显示了总阶为d＝4、两个随机变量的PCE基函数。

表四：总阶为d＝4、两个随机变量的埃尔米特多项式基函数

将(12)、(13)带入(10)中可得

下一步，我们将应用伽辽金方法，通过取(14)的内积和检验函数Ψ_l(ξ)，l＝0,…,P来确定展开系数。考虑到正交性条件，

在时间步n处，可以得到电场E_x的第l阶系数的更新方程：

这里，<,>表示内积，可以通过积分计算，例如，

式中，ρ(ξ)是随机变量ξ的概率分布函数，当随机变量相互独立时，ρ(ξ)可以写成单个随机变量的概率密度函数的乘积，如下所示：

由麦克斯韦方程推导出直角坐标系中二维TE波的电场和磁场的系数更新方程：

上面三式中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格，同时，计算区域内随机色散媒质中电位移矢量展开系数到电场展开系数/>的转换关系由(16)得到，非随机色散媒质(即真空)中/>到/>的转换关系由(2)得到(真空中的相对介电常数为ε_r＝1)，即

这样，整个区域上的电磁场系数都已计算得到，通过以下式子得到观测点电磁场的均值和标准偏差：

需要说明的是，随着随机变量的增加，相同阶次的多项式数量也会增加，这意味着需要计算更多的系数。这里电磁场展开的系数P与多项式最高阶数d的关系如下

本发明一种电磁波在随机色散介质中传播的PCE-FDTD方法，具体按照以下步骤实施：

步骤1：输入模型文件；

输入的模型文件具体为：计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δζ，ζ＝x，y，x为横坐标，y为纵坐标；设置随机色散媒质的区域；时间步长Δt；真空中的磁导率μ₀，介电常数ε₀；仿真的总时长N_t；设置源的位置；观测点的位置。

步骤2：初始化参数；

初始化的参数包括：

将整个计算区域的电磁场分量展开系数中间变量观测点电场的均值和标准差(Ex_80，sigma80)均初始化为0；同时确定展开阶数l的的值。

步骤3：计算仿真中所用到的内积；

内积计算由下式给出，

计算得到的内积是为了式(16)中电场的第l阶系数

步骤4：添加场源到中设置好源的位置，所添加场源的表达式为：

步骤5：更新计算整个计算区域的y方向上电场分量展开系数

步骤5.1：先计算电位移矢量分量系数其在计算区域的方程为：

其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

步骤5.2：求解5.1的方程得到电位移矢量分量系数后，再计算电场分量系数/>在随机色散媒质区域电场分量系数/>的计算方程为：

在非随机色散媒质区域电场分量系数的计算方程为：

步骤5.3：更新中间变量

步骤6：更新计算整个计算区域的x方向上电场分量展开系数

步骤6.1：先计算电位移矢量分量系数其在计算区域的方程为：

步骤6.2：求解6.1的方程得到电位移矢量分量系数后，再计算电场分量系数/>在随机色散媒质区域电场分量系数/>的计算方程为：/>

在非随机色散媒质区域电场分量系数的计算方程为：

步骤6.3：更新中间变量

步骤7：更新计算整个计算区域的z方向上磁场分量展开系数其计算方程为：

步骤:8：将时间迭代步t+1赋值给t，并判断迭代次数t是否达到预设值(仿真的总时长N_t)，若未达到预设值，则返回步骤5，若达到预设值，则执行步骤9；

步骤9：计算观测点的电场的均值和标准差，计算公式分别为：

计算完成，画图，结束。

实施例

平面波源对洛伦兹媒质圆柱的辐射的计算：

按照本发明的方法步骤进行实施，如图2所示，实验中整个计算区域为100×100网格，网格大小为1mm×1mm，即Δx＝Δy＝1mm。四个边界采用10层网格的PML吸收边界，计算中所加的源位于x＝20Δx，所加场源的表达式如下：

时间步长Δt＝Δx/2c0,其中c0＝3×10⁸为真空中的光速。观测点设置为(80,50)，洛伦兹媒质圆柱设置在整个区域的中间，底面半径为r＝20mm，高设置为z方向。此例子中，设置假设洛伦兹介质的共振频率ω_p和阻尼频率δ_p是随机的，且具有独立的正态分布。其中ω_p的均值为μ{ω_p}＝1.8032×10¹¹rad/s，标准差为σ{ω_p}＝0.05×μ{ω_p}；δ_p的均值为μ{δ_p}＝1×10¹⁰rad/s，标准差为σ{δ_p}＝0.05×μ{δ_p}。多项式展开最高阶数设置为d＝4。

采用本发明方法计算的观测点处的电场分量E_x的均值和标准差与传统门特卡罗方法比较的图分别如图3和图4所示。从图3中可见，当多项式展开阶数设置为d＝1时本发明方法计算得到的电场均值与MC方法计算结果一致，从图4可以看出，当多项式展开阶数设置为d＝4时本发明方法计算得到的电场标准差与MC方法计算结果吻合，以此验证了本发明方法的正确性。

本领域技术人员在考虑说明书及实践这里公开的发明后，将容易想到本发明的其它实施方案。本申请旨在涵盖本发明的任何变型、用途或者适应性变化，这些用途或者适应性变化遵循本发明的一般性原理并包括本发明未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本发明的真正范围和精神由所附的权利要求指出。

Claims (3)

1.电磁波在随机色散介质中传播的PCE-FDTD方法，其特征在于，包括以下步骤：

步骤1：输入模型文件；

步骤2：初始化参数；

步骤3：计算仿真中所用到的内积；

所述步骤3中内积的由下式计算得出

步骤4：添加场源到设置好的源位置；

所述步骤4中添加场源的表达式为：

步骤5：更新计算整个计算区域的y方向上电场分量展开系数

步骤5.1：先计算电位移矢量分量系数其在计算区域的方程为：

其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

步骤5.2：求解5.1的方程得到电位移矢量分量系数后，再计算电场分量系数在随机色散媒质区域电场分量系数/>的计算方程为：

在非随机色散媒质区域电场分量系数的计算方程为：

步骤5.3：更新中间变量

步骤6：更新计算整个计算区域的x方向上电场分量展开系数

步骤6.1：先计算电位移矢量分量系数其在计算区域的方程为：

步骤6.2：求解6.1的方程得到电位移矢量分量系数后，再计算电场分量系数在随机色散媒质区域电场分量系数/>的计算方程为：

在非随机色散媒质区域电场分量系数的计算方程为：

步骤6.3：更新中间变量

步骤7：更新计算整个计算区域的z方向上磁场分量展开系数

所述步骤7的计算方程式为：

步骤8：将时间迭代步t+1赋值给t，并判断迭代次数t是否达到预设值，若未达到预设值，则返回步骤5；若达到预设值，则执行步骤9；

步骤9：计算观测点电场的均值和标准差；

所述步骤9的计算公式为：

2.根据权利要求1所述的电磁波在随机色散介质中传播的PCE-FDTD方法，其特征在于，所述步骤1具体是：

计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δζ，ζ＝x，y，x为横坐标，y为纵坐标；设置随机色散媒质的区域，时间步长Δt，真空中的磁导率μ₀，介电常数ε₀，仿真的总时长N_t，设置源的位置，观测点的位置。

3.根据权利要求1所述的电磁波在随机色散介质中传播的PCE-FDTD方法，其特征在于，所述步骤2具体是：

将整个计算区域的电磁场分量展开系数中间变量/>观测点电场的均值和标准差(Ex_80，sigma80)均初始化为0；同时确定展开阶数l的的值。