双时间步时域有限体积方法计算时变电磁场

双时间步时域有限体积方法计算时变电磁场

许勇黄勇周铸

(中国空气动力研究与发展中心计算空气动力研究所,四川 绵阳 621000)

摘要双时间步被引入到时域有限体积解算器这一直接求解麦克斯韦方程组的高精度数值方法中．双时间步方法作为非定常时间推进技巧,时间精度由物理时间步长决定,稳定性要求由定常子迭代时间步长所满足,从而放松了通常显式时间格式和网格对物理时间步长的限制,达到节约计算量的目的．典型目标电磁散射计算表明：通过对物理时间步长、最大子迭代步数、子迭代收敛判据的合理选取,双时间步方法在保证计算精度的同时,能提高计算效率．

关键词双时间步;时域有限体积法;电磁散射;雷达截面

中图分类号O441.4

文献标志码A

文章编号1005-0388(2015)04-0647-06

AbstractDual time stepping technique was introduced into the finite volume time domain (FVTD) method which is a high precision direct numerical solvers for time dependent electromagnetic fields controlled by Maxwell equations. As an unsteady time marching technique, the time precision is determined by physical time step, stability criterion is guaranteed by steady sub-iteration time step in this method, as a result, the physical time step was relaxed which ordinary limited by explicit time scheme and grids, the computational work is also decreased. The computation of canonical electromagnetic scattering problem demonstrated that the computational precision is conserved and efficiency is also increased; by the rational choose of physical time step, the maximum step and convergence criterion of sub-iteration in this technique.

收稿日期：2014-06-26

作者简介

Dual time stepping FVTD method for computation of

time-dependent electromagnetic fields

XU YongHUANG YongZHOU Zhu

(ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchand

DevelopmentCenter,MianyangSichuan621000,China)

Key words dual time stepping; finite volume time domain method; electromagnetic scattering; radar cross section

资助项目: 装备预先研究项目(51313010303)

联系人： 周 铸 E-mail: stephen000@sina.com

引言

时变电磁场满足时域麦克斯韦方程组,随着计算机技术的发展,直接求解该方程组成为可能．与欧拉方程相同的双曲型数学特征促进了计算流体力学(Computational Fluid Dynamics, CFD)技术在电磁场计算中的应用,其中时域有限差分(Finite-Difference Time-Domain,FDTD)法[1]和时域有限体积(Finite-Volume Time-Domain,FVTD)法[2-5]最为著名．

这些方法的时间推进或采用二阶中心差分或时空耦合Lax-Wendroff格式以及多步Runge-Kutta法,其共同点是时间计算的显式格式．以Runge-Kutta方法为代表的显式方法,既有编程简便,又有易实现时间高精度的优点,对时域电磁场计算是可靠的时间离散方法．但是时间显式方法有一个最大的缺陷,其时间步长受稳定性的限制,整个计算空间必须采用统一的最小全局计算步长,为模拟几何外形剧烈变化生成的贴体加密网格,会带来很小的全局时间步长,从而使得到稳定的时变电磁场需要较长的计算时间,特别是在求解多维非定常问题时,计算时间的限制带来计算量的显著增加．另一方面,隐式计算方法虽能够放宽计算步长的稳定性限制,但伴随而来的是时间精度下降和矩阵求逆运算．

Jameson在CFD领域提出了双时间步(Dual Time Stepping)隐式迭代方法[6]．该方法通过在控制方程中引入定常的虚拟时间导数项,在每个物理时间段内对物理时间导数作线性化处理,从而使得物理时间推进步长可以根据物理问题进行选取而不受稳定性的限制,计算的稳定性要求由虚拟时间子迭代满足,子迭代的定常计算可以采用局部时间步长加速收敛．从理论上来讲,引入虚拟时间迭代过程的双时间步法,只有在子迭代残值趋于零时才接近时间高精度,但实际计算往往难以满足这一点,因此子迭代过程对非定常计算的影响是一个值得关注的问题．

双时间步法在计算效率方面具有显著的优势,大量的CFD计算充分证实了这一点．但是双时间步法求解时变电磁场时,需对一些因素进行必要的研究,包括物理时间迭代精度、物理时间步长、子迭代收敛标准以及最大子迭代步数的选择,它们协同作用才能保证时间推进的效率和计算精度．本文即在已有显式多步Runge-Kutta电磁场FVTD解算器的基础上,拓展应用子迭代为多步Runge-Kutta的双时间步法计算电磁散射场,极大缩短稳定电磁场所需计算时间,并计算和分析以上因素对计算结果的影响规律．

1数值方法

1.1控制方程

时域麦克斯韦方程组的两个旋度方程(法拉第电磁感应定律和安培环路定律),在无源条件下的直角坐标系守恒形式为

(1)

式中: Q=[BxByBzDxDyDz]T=[BD]T,

B是磁感应强度矢量,D是电位移矢量; Fx=[0

-EzEy0Hz-Hy]T; Fy=[Ez0-Ex-Hz

0Hx]T; Fz=[-EyEx0Hy-Hx0]T.

对于复杂外形物体,在计算空间生成贴体多块结构网格,其坐标变换为

k=k(x,y,z),k=ξ,η,ζ.

(2)

并得到曲线坐标系下的麦克斯韦方程组守恒型

(3)

1.2空间积分方法

在计算中,守恒电磁场量在自由空间中使用散射场,可保持守恒方程形式,同式(1),从而避免在网格空间传播入射场带来的数值误差,也有利于截断空间无反射边界条件的应用,入射电磁场可由入射波解析给出,在完全导电壁处由边界条件引入,这也相应避免了传统时域有限差分法中划分总场区和散射场区的繁琐做法．

有限体积法的空间精度体现在能否精确模拟守恒变量Q在网格单元分界面处的状态变量,以得到相应精确的分界面通量．这里采用Steger-Warming分裂进行网格单元分界面通量计算

(4)

单元边界左右流通量可统一写为:

(5)

(6)

式中: S,S-为相似矩阵; Λ+,Λ-分别为正负特征值构成的对角矩阵; QL,QR分别代表分界面处左右状态变量,它们可采用MUSCL格式而达到最高三阶精度:

(7)

(8)

式(7)～(8)中:φ是限制器,本文中取φ=1;κ=1/3是三阶精度格式的控制参数;和Δ分别是后差和前差算符．

1.3双时间步法

对于时空分别积分的半离散线方法,显式时间计算方面采用常微分方程求解中的Runge-Kutta法:

Qn+k/m=Qn-λαkR(Qn+(k-1)/m),k=1,m,

(9)

λ=Δt/V,

(10)

αk=1/m-k+1．

(11)

显式方法稳定性要求全网格空间采用相同物理时间步长,会加大计算量;双时间步方法则采用时间双重循环,外循环为物理时间推进,内循环为虚拟时间迭代,相应双时间步法控制方程组修改为

(12)

(13)

式中:

(14)

Q*是Qn+1(n是物理时间步)的近似,物理时间导数采用后向差分为二阶时间精度．式(13)类似式(9)作显式子迭代:

k=1,m,

(15)

上标l、l+1代表虚拟定常子迭代过程中的两个相邻时间层．

在数值模拟过程中,需要对内迭代收敛过程进行判断和控制,使得内迭代过程在保证一定时间精度的前提下尽快收敛．Arnone指出[8],式(15)在物理时间步小于或与虚拟时间步同量级时变得不稳定．Melson[9]进一步明确不稳定性来源于残差中3Q*/(2Δt)的显式计算,物理时间步Δt越小不稳定越明显,解决之道就是隐式化该项并作线性化处理,最后得到每层子迭代的Runge-Kutta关系为

(16)

除了式(11)对应最大时间精度系数外,可以针对空间格式调节系数来增加最大时间步和改善稳定性．表1即为对应一阶空间格式,步数分别为3、4、5的Runge-Kutta最大CFL数和相应系数．

表1　多步Runge-Kutta系数和最大CFL数

由此得到简谐波入射FVTD计算目标雷达截面(Radar Cross-Section,RCS)步骤:首先经时间推进在整个网格区域建立稳定的电磁振荡,然后在空间Stratton-Chu积分面处逐周期进行傅里叶变换,获得收敛的频域电磁散射场,进而由积分方程得到目标双站RCS分布[7]．

2计算结果与分析

双时间步方法时间精度由物理时间步长确定,为二阶时间精度;采用越小的物理时间步长时间精度就越有保证,采用较大的物理时间间隔有利于提高计算效率,但物理时间步长也不能够无节制增加．图1比较了不同物理时间步长对圆柱RCS结果的影响,k为波数,a为圆柱半径,计算网格为51×46,其中物面网格选取每波长20个网格点,远场边界在3波长外,辐向网格壁面加密,平均约每波长15个网格点．

(a) TM波

(b) TE波 图1　物理时间步长对圆柱(ka=2)RCS的影响

物理时间步长决定了傅里叶变换中取样周期的点数目,在保证精度前提下的更大取值有利于节约计算量,其选取与目标电尺度无关,跟电磁场空间散射特性相关,其选取由典型二维、三维计算验证选取．由图1知二维问题用入射电磁波周期无量纲化的物理时间步长Δt=0.01计算精度已可以,建议三维问题选择Δt≈0.001量级,随着物理时间间隔的增大,计算结果与参考值差别逐渐增大,在过大的物理时间步长条件下,双时间步方法的计算精度变低．

双时间方法非定常计算的时间精度还受到每一物理时间步上子迭代步数制约．迭代步数多,收敛控制严格,时间精度才有保证．真实时间步小,则子迭代可以取得少一些,数值计算中为防止场变化剧烈处、以及子迭代CFL数取得过大等情况时,残差很有可能无法下降到收敛判据值,陷入死循环而给出最大子迭代步数,另一方面对子迭代收敛进行判断和控制,在保证一定时间精度的前提下尽快结束迭代．本文子迭代收敛判据使用电场和磁场在计算空间最大幅度差的绝对值．

图2是TM波入射圆柱(ka=10,网格为201×46)最大子迭代步数的影响,图3是TE波入射圆柱(ka=10)子迭代收敛判据的影响,可见子迭代残差下降2个量级已能满足要求,这与CFD中残差下降2～3个量级一致．

图2　最大子迭代步数对圆柱(ka=10)RCS的影响

图3　子迭代收敛判据对圆柱(ka=10)RCS的影响

以Runge-Kutta方法为代表的显式方法容易实现时间求解的高精度,对非定常计算来说是理论上最为可靠的时间离散方法,以下针对双时间步方法与显式Runge-Kutta方法,就时间精度、结果好坏及计算效率作一些比较．

采用完全导电金属球电磁散射为例,其电尺寸为ka=2,a为球半径．计算网格为61×31×46,时间迭代都是四步Runge-Kutta推进,其他计算条件相同．图4同一空间点分别用显式Runge-Kutta法和双时间步方法计算,其散射电磁场量的时间振荡历程、空间两者几无区别,说明物理时间后向差分、二阶精度的双时间方法精度满足傅里叶变换需要．

图5比较了显式和双时间步计算的RCS双站分布,可见后者在前向(θ=180°)与解析解吻合更好．图6则是这两种方法逐周期进行傅里叶变换所得频域值收敛历程比较,在达到同样的收敛标准(0.001),双时间步收敛略快,这对大网格量和存在很小体积网格单元的情况,能大大提高计算效率．因为此时显式Δt≈10-4量级甚至更小,相反双时间步方法则根据时间精度选择Δt,并可用局部时间步等技巧加速虚拟时间子迭代的收敛．

图4　时变散射场振荡历程比较(金属球,ka=2)

图5　双时间步计算金属球双站(RCS) 分布比较(ka=2)

图6　电磁场频域收敛历程比较(金属球,ka=2)

3结论

双时间步时域有限体积方法计算时变电磁场,彻底放松了显式格式稳定性对物理时间步的限制,使得时间精度由物理时间间隔决定,稳定性由定常虚拟时间子迭代(可以应用加速收敛技巧)保证,该方法可以由添加很少几段程序由一个显式计算程序修改而来,这样既能保证计算精度又能大大提高计算效率的方法,值得在时域电磁场数值计算中推广应用．

致谢:衷心感谢中国空气动力研究与发展中心计算空气动力研究所江雄、刘刚等人在双时间方法应用上提供的帮助,使得该算法能在时变电磁场的计算中得到推广．

参考文献

[1]YEE K S. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media [J]. IEEE Trans on AP, 1966, 14(5):302-307.

[2]SHANKAR V, MOHAMMADIAN A H, HALL W F. A time-domain, finite-volume treatment for the Maxwell equations [J]. Electromagnetics, 1990, 10:127-145.

[3]SHANG J S, GAITONDE D. Scattered electromagnetic field of a reentry vehicle [C]//AIAA 94-0231, 1994.

[4] 路占波, 张安, 候新宇，等. 采用FDTD/FVTD混合算法分析蝶形微带天线[J]. 电波科学学报, 2007, 22(3): 527-532.

LU Zhanbo, ZHANG An, HOU Xinyu, et al. Analysis of bow-tie microstrip antenna by hybrid FDTD/FVTD algorithm[J]. Journal of Radio Science, 2007, 22(3): 527-532. (in Chinese)

[5]邓聪, 尹文禄, 柴舜连, 等. 两种适用时域有限体积法的截断边界[J]. 电波科学学报, 2009, 24(3): 537-540.

DENG Cong, YIN Wenlu, CHAI Shunlian, et al. Two boundary conditions for the truncation of cell-centered FVTD algorithm on unstructured lattices [J]. Chinese Journal of Radio Science, 2009, 24(3):537-540. (in Chinese)

[6]JAMESON A. Time dependent calculations using multigrid with applications to unsteady flows past airfoils and wings [C]//AIAA 91-1596, 1991.

[7]许勇, 乐嘉陵. 基于CFD的电磁散射数值模拟[J]. 空气动力学学报, 2004, 22(2): 185-189.

XU Yong, LE Jialing. CFD based numerical simulation of electromagnetic scattering [J]. Acta Aerodynamica Sinica, 2004, 22(2):185-189. (in Chinese)

[8]ARNONE A, LIOU M S, POVINELLI L A. Multigrid time accurate integration of the navier-stokes equations [C]//AIAA 93-3361, 1993.

[9]MELSON N D, SANETRIK M D, ATKINS H L. Time-accurate Navier-Stokes calculations with multigrid acceleration [C]//Proc 6th Copper Mountain Conf. on Multigrid Methods, 1993: 423-439.

许勇(1971-)，男，四川人，中国空气动力研究与发展中心副研究员，研究领域为计算电磁学，研究方向包括再入体电磁特性、电磁场数值计算方法.

黄勇(1970-)，男，四川人，中国空气动力研究与发展中心研究员．研究方向为计算空气动力学、飞行器气动多目标优化设计．

周铸 (1973-),男,重庆人,中国空气动力研究与发展中心研究员.研究方向为计算空气动力学、复杂流动分析.

胡荣旭,吴振森,张金鹏. 中国近海蒸发波导反演中最佳雷达参数分析[J]. 电波科学学报，2015,30(4):653-660. doi: 10.13443/j.cjors. 2014080501

HU Rongxu, WU Zhensen, ZHANG Jinpeng. Analysis for best radar parameter in inversion of evaporation duct above China sea areas[J]. Chinese Journal of Radio Science,2015,30(4):653-660. (in Chinese). doi: 10.13443/j.cjors. 2014080501