CN103279601B - 导体目标宽带电磁散射特性的仿真方法 - Google Patents

Abstract

本发明公开了一种导体目标宽带电磁散射特性的仿真方法，步骤如下：建立导体目标的几何模型，采用曲面三角形单元对导体目标表面进行网格剖分；确定导体目标的时域电场积分方程；采用空间上的CRWG基函数和时间上的时间滞后基函数，对时域电场积分方程中的表面感应电流展开；将表面感应电流展开表达式代入时域电场积分方程，然后对离散形式的时域电场积分方程分别在时间上、空间上测试，得到系统阻抗矩阵方程；消除奇异性积分，得到阻抗矩阵的稀疏表达式；求解阻抗矩阵的方程，确定导体目标表面的时域电流分布，根据时域电流分布得到目标的宽频带电磁特性参数，完成仿真。该方法具有仿真精度高、所需时间少、内存消耗低的优点，具有广泛的应用前景。

Description

导体目标宽带电磁散射特性的仿真方法

技术领域

本发明涉及电磁仿真技术领域，特别是一种导体目标宽带电磁散射特性的仿真方法。

背景技术

近年来，科技飞速发展，模拟超宽带信号和非线性系统的需求日益增加，因而对快速、精确和稳定地模拟分析任意复杂导体介质结构中的瞬态电磁问题提出了迫切的要求。在实际工程中，若要获得宽频响应，可以运用频域方法计算足够多的频点，然后通过傅里叶逆变换得到目标的宽带响应数据。

频域积分方程方法分析中，目标表面的感应电流是一个复数矢量，即感应电流同时包含相位和幅度信息，由导体表面的感应电荷所满足的标量亥姆霍兹方程的解以及电流连续性方程可知，感应电流的相位信息中包含了入射电磁波的相位。利用这一物理特性，将描述电流线性变化的相位信息设计到感应电流的近似展开表达式中，即用来近似展开感应电流的基函数是一个复数矢量，并称之为相位基函数；而一般用来近似展开感应电流的都是实数矢量基函数；表示电流的相位信息可以被设计到任何种类的实数矢量基函数中，从而构成新的复数基函数，相位基函数。目前国内外已有研究者将相位基函数应用于频域积分方程分析方法中，文献1（K.R.AbereggandA.F.Peterson,“Applicationoftheintegralequation-asymptoticphasemethodtotwo-dimensionalscattering,”IEEETrans.AntennasPropagat.,vol.43,no.5,pp.531-537,May.1995）提出了integralequation-asymptoticphase(IE-AP)的概念，即将快速变化相位的渐进表达引入到电磁场积分方程中，使之能处理任意形状的二维导体结构，其中IE-AP认为电流可以是代表相位信息的指数函数与缓慢变化的剩余部分的乘积；文献2（J.M.Taboada,F.Obelleiro,J.L.Rodriguez,“Incorporationlinear-phaseprogressioninRWGbasisfunction,”MicrowaveOptTchnol.Lett.44:106-112,2005）和文献3（Gareia-Tuon,J.M.Taboada,F.Obelleiro,andL.Landesa,“Efficientasymptotic-phasemodelingoftheinducedcurrentsinthefastmultipolemethod,”MicrowaveOptTechno.Lett.48:1594-1599,2006）公开了一种linear-phaseRWG(LP-RWG)基函数，即用指数表示的物体表面感应电流相位的线性变化，并将它与传统的RWG基函数相结合，这些方法可以用来快速分析任意三维导体结构的电磁散射。

上述文献1～3报道的都是频域中使用相位基函数的分析方法，然而由于谐振现象、高低频成分的存在，不同频率点的计算方法也不同，导致计算量庞大和复杂性高，并且不能直观的从模拟结果中理解场的相互作用，从而使得频域方法失去了优势。

发明内容

本发明的目的在于提供一种高效稳定的导体目标宽带电磁散射特性的仿真方法，该方法能显著提高仿真效率，具有内存消耗低、仿真时间快的特点。

实现本发明目的的技术方案为：一种导体目标宽带电磁散射特性的仿真方法，步骤如下：

第1步，建立导体目标的几何模型，采用曲面三角形单元对导体目标的表面进行网格剖分；

第2步，根据时域形式的麦克斯韦方程组和电流连续性，确定导体目标的时域电场积分方程；

第3步，采用空间上的CRWG基函数和时间上的时间滞后基函数，对第2步的时域电场积分方程中的表面感应电流进行展开，得到表面感应电流展开表达式；

第4步，将第3步的表面感应电流展开表达式代入第2步的时域电场积分方程中，然后对离散形式的时域电场积分方程分别在时间上采用点测试、在空间上采用Galerkin测试，得到系统阻抗矩阵方程；

第5步，根据第4步中阻抗矩阵元素的表达式，消除奇异性积分，得到阻抗矩阵的稀疏表达式；

第6步，根据第4～5步得到的系统阻抗矩阵方程，求解阻抗矩阵的方程，确定导体目标表面的时域电流分布，根据时域电流分布得到目标的宽频带电磁特性参数，完成仿真过程。

本发明与现有技术相比其显著效果是：(1)求解未知量少：时间滞后基函数对真实时域感应电流的描述更加准确，从而允许采用更大尺寸的单元贴片离散目标表面，例如采用了三角时间滞后基函数，此时单元贴片的最大剖分尺寸可以达到0.4c/f_max，f_max为入射电磁波的最高频率，这与一般的三角时间基函数，单元贴片的最大剖分尺寸为0.1c/f_max，相比大大节省了求解所需要的未知量；(2)模型离散拟合更精确：采用曲面三角形单元对仿真对象的表面进行网格离散，能真实地拟合各种复杂的几何形状，保证了外形逼近的精确性；(3)对导体目标几何结构的适应性好，计算结果精确：时域电场积分方程可以同时用于分析闭合结构和开放结构的金属目标，径向角坐标变换消去奇异性的方法保证了计算结果的精确性。

附图说明

图1是本发明导体目标受到电磁波照射示意图。

图2是本发明曲面三角形单元的示意图。

图3是图2中曲面三角形单元在面积坐标下表示的示意图。

图4是本发明场点在曲面源三角形单元上的投影示意图。

图5是本发明切平面三角形被划分成3个子三角形示意图。

图6是对图5中的子三角形进行坐标转换的示意图。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细描述。

本发明为基于时间滞后基函数的时域电场积分方程方法，首先根据入射电磁波的时间滞后量设计一种时间滞后基函数，然后将空间基函数结合时间滞后基函数用来近似展开时域的感应电流，并将电流近似展开表达式代入时域电场积分方程，对离散后的积分方程分别进行时间上的点测试和空间上的伽辽金测试，形成矩阵方程，采用基于广义最小余量法的时间步进算法求解系统矩阵方程，得到每个时刻的感应电流分布，最后利用时域电流分布计算得到目标的宽频带电磁特性参数。

下面结合附图，以图1所示表面为S的任意形状的导体为例，对本发明的具体步骤做进一步详细描述。

本发明金属目标宽带电磁散射特性的仿真方法，步骤如下：

第1步，建立导体目标的几何模型，采用曲面三角形单元对导体目标的表面进行网格剖分，具体过程如下：

（1.1）如图1所示，建立导体目标的几何模型，利用计算机辅助设计工具ANSYS软件进行几何建模，导体目标置于介电常数ε₀、导磁率μ₀的自由空间中，导体目标在外来入射电磁波E^inc(r,t)照射下激励起表面感应电流，入射电磁波E^inc(r,t)为高斯调制脉冲，其表达式为：

$E^{\mathrm{inc}}(r,t)=E_0\cos \left[{2\mathrm{\πf}}_0(t-\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c})\right]\exp \left[-\frac{{(t-t_p-\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c})}^2}{{2\mathrm{\σ}}^2}\right]---\left(1\right)$

式（1）中E^inc(r,t)为t时刻导体目标r点处的入射电磁波，其中r为导体目标观察点的位置矢量；E₀为高斯调制脉冲的强度；t_p为高斯调制脉冲E^inc(r,t)的时间中心；σ＝6/(2πf_bw)，其中f_bw为脉冲的频带宽度；f₀为入射电磁波的中心频率，初始频率f_min＝0，则最高频率f_max＝2f₀，f_bw＝2f₀；表示与导体目标自身有关的时间滞后量，表示入射电磁波的单位方向矢量，c为自由空间中的光速。

（1.2）采用曲面三角形单元对金属目标的表面进行网格剖分，根据入射电磁波的最高频率f_max来确定曲面三角形单元的尺寸，保证离散得到的三角形中的最大边长l_max满足条件l_max≤0.4c/f_max，其中c为自由空间中的光速；得到仿真所需的网格离散信息文件，包括离散的曲面三角形的单元信息文件和节点信息文件，运用前处理软件得到更详细的几何信息，包括单元的编号、节点的编号、内边的编号、节点的三维坐标。

第2步，根据时域形式的麦克斯韦方程组和电流连续性，确定导体目标的时域电场积分方程，具体步骤如下：

（2.1）假设入射波在t＝0时刻以后到达导体目标，即t＜0时，表面感应电流J(r,t)＝0，t时刻导体目标r点处的感应电流J(r,t)在空间中将产生散射电磁波E^sca(r,t)，t时刻导体目标r点处的总场E^total(r,t)为入射电磁波E^inc(r,t)与散射电磁波E^sca(r,t)的矢量和，根据金属目标表面切向电场为零的连续边界条件可得：

$\hat{n}\left(r\right)\×[E^{\mathrm{inc}}(r,t)+E^{\mathrm{sca}}(r,t)]=0---\left(2\right)$

式中E^sca(r,t)为t时刻金属目标r点处的散射电磁波，是金属目标表面S在r点处的单位外法向矢量；

（2.2）E^sca(r,t)用表面感应电流表达的形式为：

$E^{\mathrm{sca}}(r,t)=-\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫}_S\frac{1}{R}\frac{\∂J(r^{\′},\mathrm{\τ})}{\∂t}{\mathrm{dS}}^{\′}+\frac{1}{{4\mathrm{\π\ϵ}}_0}{\▿}_S{\∫}_S{\∫}_{-\∞}^{\mathrm{\τ}}\frac{{\▿}_S^{\′}\·J(r^{\′},t^{\′})}{R}{\mathrm{dt}}^{\′}{\mathrm{dS}}^{\′}---\left(3\right)$

上式中r′表示金属目标源点的位置矢量，μ₀为自由空间中磁导率，ε₀为自由空间中介电常数，R＝|r-r′|代表金属目标观察点和源点之间的空间距离，τ＝t-|r-r′|/c表示位于金属目标源点r′处产生的场传到金属目标观察点即场点r处需要的时间滞后，c＝3.0×10⁸m/s为电磁波在自由空间中的传播速度，J(r′,τ)表示金属目标源点r′处产生散射电磁场的时变电流密度；J(r′,t′)中t′是积分的变量，分别表示面积的散度算子和梯度算子。

（2.3）金属目标的时域电场积分方程的基本形式为：

$\hat{n}\left(r\right)\×[\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫}_T\frac{1}{R}\frac{\∂J(r^{\′},\mathrm{\τ})}{\∂t}{\mathrm{dS}}^{\′}-\frac{1}{{4\mathrm{\π\ϵ}}_0}{\▿}_S{\∫}_T{\∫}_{-\∞}^{\mathrm{\τ}}\frac{{\▿}_S^{\′}\·J(r^{\′},t^{\′})}{R}{\mathrm{dt}}^{\′}{\mathrm{dS}}^{\′}]=\hat{n}\left(r\right)\×E^{\mathrm{inc}}(r,t)---\left(4\right)$

第3步，采用空间上的CRWG基函数和时间上的时间滞后基函数，对第2步的时域电场积分方程中的表面感应电流进行展开，得到表面感应电流展开表达式，具体过程如下：

（3.1）选取空间基函数，如附图2所示，采用定义在曲面三角形单元上的空间基函数即CRWG基函数对t时刻金属目标r′点处的感应电流J(r′,t)进行空间上的近似展开，展开表达式如下：

$J(r^{\′},t)\≅\underset{n=1}{\overset{N_s}{\mathrm{\Σ}}}{J}_n\left(t\right){\mathrm{\Λ}}_n\left(r^{\′}\right)---\left(5\right)$

式中N_s是第1步中金属目标表面离散所得到的内边的个数，即空间基函数的个数，n表示空间基函数的编号，J_n(t)是待求的未知电流系数，需要用时间基函数来展开，Λ_n(r′)表示源点r′处的空间基函数即CRWG基函数，表达式如下：

${\mathrm{\Λ}}_n^{\mathrm{\β}}\left(r^{\′}\right)=\frac{\±1}{J}({\mathrm{\ξ}}_{\mathrm{\β}+1}{I}_{\mathrm{\β}-1}-{\mathrm{\ξ}}_{\mathrm{\β}-1}{I}_{\mathrm{\β}+1}),\mathrm{\β}=\mathrm{1,2,3};n=\mathrm{1,2},\·\·\·N_s---\left(6\right)$

式中J表示将曲面三角形单元转换到参数空间的标准三角形单元的雅克比因子，ξ₁、ξ₂、ξ₃表示r′的面积坐标，如附图3所示，且有如下关系ξ₁+ξ₂+ξ₃=1；Ι₁,Ι₂,Ι₃表示基向量，Ι₁,Ι₂,Ι₃的几何意义如附图2所示，r′用面积坐标ξ₁、ξ₂、ξ₃表示的表达式如下：

r′(x,y,z)＝ξ₁(2ξ₁-1)r₁+ξ₂(2ξ₂-1)r₂+ξ₃(2ξ₃-1)r₃（7）

+4ξ₁ξ₂r₄+4ξ₂ξ₃r₅+4ξ₃ξ₁r₆

式中r₁～r₆是曲面三角形的6个节点，如附图2中所示。

（3.2）构造时间滞后基函数，根据入射电磁波的形式构造时间滞后基函数，选取一般的三角时间基函数（又称一阶Lagrange插值时间基函数）作为将被改造的时间基函数，一般的三角时间基函数的形式如下：

上式中Δt是为时间分辨率，即时间步长，满足c是自由空间的光速，f_max是入射电磁波的最高频率，根据第1步公式（1）中入射电磁波的表达式中的时间滞后量构造时间滞后基函数，因为被改造的时间基函数是三角时间基函数，所以这里构造的时间滞后基函数称之为三角时间滞后基函数，形式如下：

式中表示入射电磁波的单位方向矢量，r′表示源点的位置。

（3.3）将三角时间滞后基函数的表达式代入到第（3.1）步中式（5）的时域感应电流的近似展开表达式，得到如下形式：

$J(r^{\′},t)\≅\underset{n=1}{\overset{N_s}{\mathrm{\Σ}}}\underset{l=1}{\overset{N_t}{\mathrm{\Σ}}}{I}_n^{l}T(t-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c}){\mathrm{\Λ}}_n\left(r^{\′}\right)---\left(10\right)$

N_t表示时间基函数的个数，即离散的时间步的个数，l是时间基函数的编号I_n,l是第n个空间基函数上第l个时间步上的时间基函数的系数，Λ_n(r′)是空间CRWG基函数。

第4步，将第3步的表面感应电流展开表达式（10）代入第2步的时域电场积分方程（4）中，然后对离散形式的时域电场积分方程分别在时间上采用点测试、在空间上采用Galerkin测试，得到系统阻抗矩阵方程，具体过程如下：

（4.1）表面感应电流展开表达式（10）代入时域电场积分方程（4），得到离散形式的积分方程形式如下：

$\hat{n}\left(r\right)\×\underset{n}{\overset{N_s}{\mathrm{\Σ}}}\underset{l=1}{\overset{N_t}{\text{\Σ}}}{I}_n^l\left\{\begin{array}{c}\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫}_{T_n}\frac{1}{R}\mathrm{\Λ}\left(r^{\′}\right)\frac{\∂T(t-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c})}{\∂t}{\mathrm{dS}}^{\′}-\\ \frac{1}{{4\mathrm{\π\ϵ}}_0}{\▿}_S{\∫}_{T_n}{\∫}_{-\∞}^{t-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c}}\frac{{\▿}_S^{\′}\·\left[\mathrm{\Λ}\left(r^{\′}\right)T\left(r^{\′}\right)\right]}{R}{\mathrm{dt}}^{\′}{\mathrm{dS}}^{\′}\end{array}\right\}=\hat{n}\left(r\right)\×E^{\mathrm{inc}}(r,t)---\left(11\right)$

（4.2）时间维上的点测试过程：在时间维采用点匹配方法，即利用Delta函数对上式中的离散时域电场积分方程在时间上作内积，得到以下表达式：

$\hat{n}\left(r\right)\×\underset{n}{\overset{N_s}{\mathrm{\Σ}}}\underset{l=1}{\overset{N_t}{\text{\Σ}}}{I}_n^l\left\{\begin{array}{c}\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}\∫{\∫}_{T_n}\frac{1}{R}{\mathrm{\Λ}}_n\left(r^{\′}\right)\frac{\∂T(\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c}-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c})}{\∂t}{\mathrm{dS}}^{\′}-\\ \frac{1}{{4\mathrm{\π\ϵ}}_0}{\▿}_S\∫{\∫}_{T_n}{\∫}_{-\∞}^{\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c}-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c}}\frac{{\▿}_S^{\′}\·\left[{\mathrm{\Λ}}_n\left(r^{\′}\right)T\left(r^{\′}\right)\right]}{R}{\mathrm{dt}}^{\′}{\mathrm{dS}}^{\′}\end{array}\right\}---\left(12\right)$

$=\hat{n}\left(r\right)\×E^{\mathrm{inc}}(r,\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c})$

（4.3）空间上的Galerkin测试过程：利用检验函数Λ_m(r)(m＝1,2,…,N_s)对上面的方程式在空间上做内积，分别可以得到N_s个方程组，第m个方程如下：

${\∫\∫}_{T_m}{\mathrm{dS\Λ}}_m\left(r\right)\·\{\hat{n}\left(r\right)\×\underset{n}{\overset{N_s}{\mathrm{\Σ}}}\underset{l=1}{\overset{N_t}{\mathrm{\Σ}}}{I}_n^l\left\{\begin{array}{c}\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫\∫}_{T_n}\frac{1}{R}{\mathrm{\Λ}}_n\left(r^{\′}\right)\frac{\∂T(\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c}-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c})}{\∂t}{\mathrm{dS}}^{\′}-\\ \frac{1}{{4\mathrm{\π\ϵ}}_0}{\▿}_S{\∫\∫}_{T_n}{\∫}_{-\∞}^{\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c}-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c}}\frac{{\▿}_S^{\′}\left[{\mathrm{\Λ}}_n\left(r^{\′}\right)T\left(t^{\′}\right)\right]}{R}{\mathrm{dt}}^{\′}{\mathrm{dS}}^{\′}\end{array}\right\}\}---\left(13\right)$

$={\∫\∫}_{T_m}{\mathrm{dS\Λ}}_m\left(r\right)\·[\hat{n}\left(r\right)\×E^{\mathrm{inc}}(r,\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c})]$

将上式（13）N_s个方程改成阻抗矩阵方程的形式，如下：

${\stackrel{\‾}{Z}}_E^0{I}^k=V_E^k-\underset{l-1}{\overset{k-1}{\mathrm{\Σ}}}{\stackrel{\‾}{Z}}_E^{k-1}{I}^l---\left(14\right)$

上式的矩阵方程表示在第k个时间步需要求解I^k，是第k个时间步的待求系数，是第k个时间步的入射电磁波，表示当前时刻即第k个时间步的阻抗矩阵，由于第k个时间步以前的电流系数I^l在求解第k个时间步的时候都是已知的，l＝1,2,...,k-1，所以和I^l有关的矩阵元素全部放在等式的右边：

${\left[{\stackrel{\‾}{Z}}_E^{k-1}\right]}_{\mathrm{mn}}={\∫\∫}_{T_m}{\mathrm{dS\Λ}}_m\left(r\right)\·\{\hat{n}\left(r\right)\×\left\{\begin{array}{c}\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫\∫}_{T_n}\frac{1}{R}{\mathrm{\Λ}}_n\left(r^{\′}\right)\frac{\∂T(\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c}-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c})}{\∂t}{\mathrm{dS}}^{\′}-\\ \frac{1}{{4\mathrm{\π\ϵ}}_0}\▿{\∫\∫}_{T_n}{\∫}_{-\∞}^{\mathrm{k\Δt}+\frac{r\·{\hat{k}}^{\mathrm{inc}}}{c}-\mathrm{l\Δt}-\frac{r^{\′}\·{\hat{k}}^{\mathrm{inc}}}{c}}\frac{{\▿}^{\′}\·\left[{\mathrm{\Λ}}_n\left(r^{\′}\right)T\left(t^{\′}\right)\right]}{R}{\mathrm{dt}}^{\′}{\mathrm{dS}}^{\′}\end{array}\right\}\}---\left(15\right)$

第5步，根据第4步中阻抗矩阵的表达式（14），消除奇异性积分，得到阻抗矩阵的稀疏矩阵表达式；当场点r和源点r′的距离R＝|r-r′|等于或趋于0时，例如场三角形单元与源三角形单元重合的时刻，积分的数值结果等于无穷大，此时必须处理积分奇异性，采用径向角坐标变换法消除奇异性积分的具体步骤如下：

（5.1）如附图2、图3所示，在场点r处构造曲面源三角形单元的切平面三角形，已知场点r(x,y,z)＝r(ξ₁,ξ₂)的面积坐标(ξ₁,ξ₂)，切平面三角形的三个顶点坐标的确定如下：

$\left\{\begin{array}{c}{r}_1^t(x,y,z)=(1-{\mathrm{\ξ}}_1)\frac{\∂r}{{\∂\mathrm{\ξ}}_1}-{\mathrm{\ξ}}_2\frac{\∂r}{{\∂\mathrm{\ξ}}_2}+r(x,y,z);\\ r_2^t(x,y,z)=-{\mathrm{\ξ}}_1\frac{\∂r}{{\∂\mathrm{\ξ}}_1}+(1+{\mathrm{\ξ}}_2)\frac{\∂r}{{\∂\mathrm{\ξ}}_2}+r(x,y,z);\\ r_3^t(x,y,z)=-{\mathrm{\ξ}}_1\frac{\∂r}{{\∂\mathrm{\ξ}}_1}-{\mathrm{\ξ}}_2\frac{\∂r}{{\∂\mathrm{\ξ}}_2}+r(x,y,z)\end{array}\right.---\left(16\right)$

（5.2）根据场点r的位置将切平面三角形分成三个小三角形T₁,T₂,T₃，如附图5所示，源三角形上的积分转换成三个小三角形的积分之和，内层源三角形上的积分有如下表达形式：

${\∫\∫}_{T_n}F(r,r^{\′})/{\mathrm{RdS}}^{\′}$

$=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{\∫\∫}_{T_{n_i}}F(r,r^{\′})/{\mathrm{RdS}}_i^{\′}---\left(17\right)$

上式表示将源三角形上的积分转换成三个小三角形上的积分之和，F(r,r′)是任意矢量形式的被积表达式，i＝1,2,3表示将源三角形分成三个小三角形，如附图5所示。

（5.3）利用坐标变换消去每一个子三角形中积分的奇异性。以子三角形T₃为例,如附图6所示，将三角形所在的笛卡尔坐标系转化到x′o′y′局部坐标系，同时将子三角形T₃上的面积分转换成二重线积分，表达式如下：

${\∫\∫}_{T_{n_3}}F(r,r^{\′})/{\mathrm{RdS}}_i^{\′}$ （18）

$={\∫}_0^{\left|h_1^{\′}\right|}{\∫}_{{y\cot \mathrm{\φ}}_L}^{{y\cot \mathrm{\φ}}_U}F(r,r^{\′})/{\mathrm{Rdx}}^{\′}{\mathrm{dy}}^{\′}$

|h₁′|表示矢量h₁′的模值即高的长度，h₁′、φ_L、φ_U的几何意义如附图6所示，表示编号为n的源三角形被分成3个小三角形中的第3个。

（5.4）再通过坐标变换，将x′o′y′局部坐标系变换到径向角坐标系R-φ，过程用以下表达式描述：

$I_{\mathrm{st}}^i={\∫}_{{\mathrm{\φ}}_L}^{{\mathrm{\φ}}_U}{\∫}_{R_L}^{R_U}\frac{F(r,r^{\′})}{R}\mathrm{RdRd\φ}---\left(19\right)$

$={\∫}_{{\mathrm{\φ}}_L}^{{\mathrm{\φ}}_U}{\∫}_{R_L}^{R_U}F(r,r^{\′})\mathrm{dRd\φ}$

上式中R_L＝|z|，其中|z|的几何意义如附图5所示。

由上式可以看出被积表达式中的在坐标变换的过程中被消去了，所以积分不存在奇异性了，可采用高斯-勒让德(Gauss-Legendre)数值积分计算。

第6步，根据第4～5步得到的阻抗矩阵的稀疏矩阵表达式，采用基于广义最小余量法的时间步进算法求解阻抗矩阵的方程，确定导体目标表面的时域电流分布，根据时域电流分布得到目标的宽频带电磁特性参数，完成仿真过程。时间步进算法MOT指的是在时间上采用点匹配的方法，使得时域电场积分方程可以离散为在时间上的递推的矩阵方程；每一个时间步需要使用广义最小余量法求解一次矩阵方程，有N_t个时间步就需要求解N_t次矩阵方程。

综上所述，本发明提出了一种时间滞后基函数，在任意类型的时间基函数中的时间变量上加上合理的时间滞后量来构造相应的时间滞后基函数，使得每个单元上电流的描述更加符合实际的物理现象，从而可以使用更少的剖分单元来逼近真实的时域感应电流。该方法从本质上减少了计算未知量，降低了内存消耗，具有计算结果精度高，计算时间少的优点，可为电大尺寸金属目标宽带电磁散射特性的精确分析提供重要的参考资料。

Claims (5)

1.一种导体目标宽带电磁散射特性的仿真方法，其特征在于，步骤如下：

第1步，建立导体目标的几何模型，采用曲面三角形单元对导体目标的表面进行网格剖分；

第2步，根据时域形式的麦克斯韦方程组和电流连续性，确定导体目标的时域电场积分方程；

第3步，采用空间上的CRWG基函数和时间上的时间滞后基函数，对第2步的时域电场积分方程中的表面感应电流进行展开，得到表面感应电流展开表达式；

第4步，将第3步的表面感应电流展开表达式代入第2步的时域电场积分方程中，然后对离散形式的时域电场积分方程分别在时间上采用点测试、在空间上采用Galerkin测试，得到系统阻抗矩阵方程；

第5步，根据第4步中阻抗矩阵元素的表达式，消除奇异性积分，得到阻抗矩阵的稀疏表达式；

第6步，根据第4～5步得到的系统阻抗矩阵方程，求解阻抗矩阵的方程，确定导体目标表面的时域电流分布，根据时域电流分布得到目标的宽频带电磁特性参数，完成仿真过程。

2.根据权利要求1所述的导体目标宽带电磁散射特性的仿真方法，其特征在于，第1步中所述建立导体目标的几何模型，采用曲面三角形单元对导体目标的表面进行网格剖分，具体过程如下：

(1.1)建立导体目标的几何模型，导体目标置于介电常数ε₀、导磁率μ₀的自由空间中，导体目标在入射电磁波E^inc(r,t)照射下激励起表面感应电流，入射电磁波E^inc(r,t)为高斯调制脉冲，其表达式为：

$E^{inc}\left(r,t\right)=E_0\cos \[2{\mathrm{\πf}}_0\left(t-\frac{r\·{\hat{k}}^{inc}}{c}\right)\]\exp \[-\frac{{\left(t-t_p-\frac{r\·{\hat{k}}^{inc}}{c}\right)}^2}{2{\mathrm{\σ}}^2}\]$

式中E^inc(r,t)为t时刻导体目标r点处的入射电磁波，其中r为导体目标观察点的位置矢量；E₀为高斯调制脉冲的强度；t_p为高斯调制脉冲E^inc(r,t)的时间中心；σ＝6/(2πf_bw)，其中f_bw为脉冲的频带宽度；f₀为入射电磁波的中心频率，初始频率f_min＝0，则最高频率f_max＝2f₀，f_bw＝2f₀；表示与导体目标自身有关的时间滞后量，表示入射电磁波的单位方向矢量，c为自由空间中的光速；

(1.2)采用曲面三角形单元对导体目标的表面进行网格剖分，根据入射电磁波的最高频率f_max来确定曲面三角形单元的尺寸，保证离散得到的三角形中的最大边长l_max满足条件l_max≤0.4c/f_max，其中c为自由空间中的光速，并得到仿真所需的网格离散信息文件，网格离散信息文件包括离散的曲面三角形的单元信息文件和节点信息文件，并得到更详细的几何信息，包括单元的编号、节点的编号、内边的编号、节点的三维坐标。

3.根据权利要求2所述的导体目标宽带电磁散射特性的仿真方法，其特征在于，第3步中所述的对第2步的时域电场积分方程中的表面感应电流进行展开的具体步骤如下：

(3.1)选取空间基函数，采用定义在曲面三角形单元上的空间基函数即CRWG基函数对t时刻导体目标r′点处的感应电流J(r′,t)进行空间上的展开，展开表达式如下：

$J(r^{\′},t)\≅\underset{n=1}{\overset{N_s}{\Σ}}{J}_n\left(t\right){\mathrm{\Λ}}_n\left(r^{\′}\right)$

式中N_s是第1步中导体目标表面离散所得到的内边的个数，即空间基函数的个数，n表示空间基函数的编号，J_n(t)是待求的未知电流系数，需要用时间基函数来展开，Λ_n(r′)表示源点r′处的空间基函数即CRWG基函数，表达式如下：

${\mathrm{\Λ}}_n^{\mathrm{\β}}\left(r^{\′}\right)=\frac{\±1}{J}({\mathrm{\ξ}}_{\mathrm{\β}+1}{I}_{\mathrm{\β}-1}-{\mathrm{\ξ}}_{\mathrm{\β}-1}{I}_{\mathrm{\β}+1}),\mathrm{\β}=2;n=1,2,...N_s$

式中J表示将曲面三角形单元转换到参数空间的标准三角形单元的雅克比因子，ξ₁、ξ₂、ξ₃表示r′的面积坐标，且有如下关系ξ₁+ξ₂+ξ₃＝1；Ι₁,Ι₂,Ι₃表示基向量， $I_1=-\frac{\∂r^{\′}}{\∂{\mathrm{\ξ}}_2},I_2=\frac{\∂r^{\′}}{\∂{\mathrm{\ξ}}_1},I_3=\frac{\∂r^{\′}}{\∂{\mathrm{\ξ}}_2}-\frac{\∂r^{\′}}{\∂{\mathrm{\ξ}}_1};$ r′用面积坐标ξ₁、ξ₂、ξ₃表示的表达式如下：

$\begin{array}{c}{r}^{\′}(x,y,z)={\mathrm{\ξ}}_1(2{\mathrm{\ξ}}_1-1)r_1+{\mathrm{\ξ}}_2(2{\mathrm{\ξ}}_2-1)r_2+{\mathrm{\ξ}}_3(2{\mathrm{\ξ}}_3-1)r_3\\ +4{\mathrm{\ξ}}_1{\mathrm{\ξ}}_2{r}_4+4{\mathrm{\ξ}}_2{\mathrm{\ξ}}_3{r}_5+4{\mathrm{\ξ}}_3{\mathrm{\ξ}}_1{r}_6\end{array}$

式中r₁～r₆是曲面三角形的6个节点；

(3.2)构造时间滞后基函数，根据入射电磁波的形式构造时间滞后基函数，选取三角时间基函数作为将被改造的时间基函数，三角时间基函数的形式如下：

上式中Δt为时间分辨率，即时间步长，满足c是自由空间的光速，f_max是入射电磁波的最高频率，根据第(1.1)步中入射电磁波的表达式中的时间滞后量构造三角时间滞后基函数，形式如下：

$T\left(t-{\hat{k}}^{inc}\·r^{\′}\right)=\left\{\begin{array}{cc}1+\frac{t-{\hat{k}}^{inc}\·r^{\′}}{\mathrm{\Δ}t}& {\hat{k}}^{inc}\·r^{\′}-\mathrm{\Δ}t\≤t\≤{\hat{k}}^{inc}\·r^{\′};\\ 1-\frac{t-{\hat{k}}^{inc}\·r^{\′}}{\mathrm{\Δ}t}& {\hat{k}}^{inc}\·r^{\′}\≤t\≤{\hat{k}}^{inc}\·r^{\′}+\mathrm{\Δ}t;\\ 0& else\end{array}\right.$

式中表示入射电磁波的单位方向矢量，r′表示源点的位置；

(3.3)将三角时间滞后基函数的表达式代入到第(3.1)步中时域感应电流的展开表达式，得到下式：

$J(r^{\′},t)\≅\underset{n=1}{\overset{N_s}{\Σ}}\underset{l=1}{\overset{N_t}{\Σ}}{I}_n^{l}T(t-l\mathrm{\Δ}t-\frac{r\·{\hat{k}}^{inc}}{c}){\mathrm{\Λ}}_n\left(r^{\′}\right)$

N_t表示时间基函数的个数，即离散的时间步的个数，l是时间基函数的编号，I_n,l是第n个空间基函数上第l个时间步上的时间基函数的系数，Λ_n(r′)是空间CRWG基函数。

4.根据权利要求1所述的导体目标宽带电磁散射特性的仿真方法，其特征在于，第5步中所述消除奇异性积分的方法为径向角坐标变换法。

5.根据权利要求1所述的导体目标宽带电磁散射特性的仿真方法，其特征在于，第6步中所述求解时域电场积分方程采用的是基于广义最小余量法的时间步进法。