用精确控制法计算复杂区域的二维电磁散射场

[文章编号]　10012246X (2000)0420414207
用精确控制法计算复杂区域的
二维电磁散射场陈红全,　黄明恪
(南京航空航天大学,江苏南京　210016)
[摘　要]　描述二维电磁场的Maxwell 方程(TM 模态)可转化为等价的波动方程。

采用精确控制
法计算该波动方程的时间周期解,数值模拟了包括外挂、鼓包以及舱体干扰的复杂气动外形电磁
散射场,并计算给出了对应的雷达散射截面(RCS )。

由于采用了非结构三角形网格,特别适合计算
复杂的实用气动外形。

[关键词]　精确控制法;Maxwell 方程;气动外形
[中图分类号]　O35419 [文献标识码]　A
[收稿日期]1999202204;[修订日期]1999208230
[作者简介]陈红全(1962～)男,浙江余姚,副教授,主要从事流体力学研究,南京航空航天大学六系　210016.
0　引言
隐身气动设计中首要考虑的重要问题是减弱飞行器的雷达反射信号实现雷达隐身[1]。

气动与隐身这对经常相互矛盾的设计要求,对于多部件组合的复杂飞行器显得尤其突出,要求发展比早期工程估算法精度更高的隐身特性计算方法。

近年来随着计算机的发展,出现了借鉴计算流体力学中各种先进的Euler 方程解算方法,直接求解电磁场Maxwell 基本方程的数值方法[2]。

另一类同步发展而引人注意的是基于精确控制的方法[3],从文献算例看,可模拟多部件干扰的电磁散射场,具有很强的适应性[4,5]。

弹舱形舱体、通信或整流鼓包以及副油箱或导弹等外挂是飞机气动设计中不可回避的。

本文采用精确控制法,数值模拟了包括外挂、鼓包以及舱体干扰的复杂气动外形电磁散射场,并计算给出了对应的雷达散射截面(RCS )。

由于采用了非结构三角形网格,本方法特别适合计算复杂的实用气动外形。

1　电磁场散射问题
对于周期为T 的入射电磁波,理想导体的电磁散射场满足如下波动方程[4]
εu tt -Δ・μ-1Δu =0, Q (=Ω×(0,T ))内u =g ,
σ(=γ×(0,T ))上9u 9n +εμ9u 9t
=0,Σ(=Γ×(0,T ))上　(1)　(2)
(3)其中,γ和Γ分别为散射体边界和远场边界;ε和μ为介质常数,Ω为散射区域。

由周期性要
第17卷第4期
2000年7月计 算 物 理　CHIN ESE J OU RNAL OF COMPU TA TIONAL PH YSICS Vol.17,No.4J ul.,2000
求,u 还须满足u (0)=u (T ), u t (0)=u t (T )(4)
求解方程(
1)～(3)的周期解与求解二维电磁场的Maxwell 方程(TM 模态)等价[4]。

2　精确控制法
上述问题可归结为精确控制问题,涉及精确控制函数的定义及其导数的确定等。

这里将分别作出概要说明。

211　精确控制问题
问题(1)～(4)的求解等价于搜寻一对{e 0,e 1},使得
u (0)=e 0, u t (0)=e 1
(5)u (T )=e 0, u t (T )=e 1(6)
成立,其中u 是方程(1)～(4)的解。

这便是精确控制问题[5],可通过定义精确控制函数,使问题归结为精确控制函数的极值问题。

212　二次函数极小化问题
对于上述精确控制问题,可有不同控制函数[5],本文沿用文[4]的定义,其对应的极值问题可描述为如下二次函数极小化问题
min ν∈E
J (ν)(7)其中E 为包含{e 0,e 1}的空间:E =V g ×L 2(Ω),V g ={φ|φ∈H 1(Ω),φ|γ=g (0)},控制函数J 由下式给出
J (ν)=12∫Ω
[μ-1|Δ(y (T )-ν0)|2+ε|y t (T )-ν1|2]d x ,　Πν={ν0,ν1}(8)其中y 是方程(1)～(3)中用u 代成y ,并以y (0)=ν0,y t (0)=ν1为初始条件,求得的t =T 时的终点值,简记为
{y (0)=ν0,y t (0)=ν1}→{y (T ),y t (T )}
(9) 现假定e 是问题(7)的解,那么它将满足
〈J ′(e ),ν〉=0,　Π
ν∈E 0(=V 0×L 2(Ω),　V 0={φ|φ∈H 1(Ω),φ|γ=0})(10) 这一问题可应用共轭梯度法求解,为此需要确定控制函数的导数J ′,这通常是应用控制方法的关键。

213　导数J ′的确定及其伴随方程
Πw ={w 0,w 1}∈E 0,可以证明控制函数的导数J ′满足[4]
〈J ′(ν),w 〉=∫Ω
μ-1Δ(ν0-y (T ))・Δw 0d s -∫Ωεp t (0)w 0d s +∫Γεμ-1
p (0)w 0
d Γ+∫Ωεp (0)w 1d s +∫Ωε(ν1-y t (T ))w 1d s (11)其中,{y (T ),y t (T )}由(9)式给出,p 称为‘光滑’函数,满足如下伴随方程
εp tt -Δ・μ-1Δp =0, Q 内p =0,σ上9p 9n -εμ9p 9t
=0,Σ上(12)(13)(14)
5
14　第4期
陈红全等:用精确控制法计算复杂区域的二维电磁散射场
614计 算 物 理第17卷　
求解时,p满足t=T终点条件
p(T)=y t(T)-ν1
∫Ωεp t(T)z d s=∫Γεμ-1(y t(T)-ν1)z dΓ-∫Ωμ-1Δ(y(T)-ν0)・Δz d s,
(15)
Πz∈V0
式(11)中初始值{p(0),p t(0)}就是求解伴随方程(负向波动方程)(12)～(15)得到的,简记为{p(T),p t(T)}{ν0,ν1}→{p(0),p t(0)},　Πz∈V0(16) 214　共轭梯度解法简介
这里主要是结合共轭梯度法给出求解问题(10)的算法结构。

第零步:初始化
给一初始猜测e0={e00,e01}∈E,先求正向波动方程(9),即
{y0(0)=e00,y0t(0)=e01}→{y0(T),y0t(T)}
接着求负向波动方程(16),即
{p0(T),p0t(T)}{e00,e01}→{p0(0),p0t(0)}, Πz∈V0
再根据(11)式,定义g0={g00,g01}∈E0,并由如下公式求出g00和g01,即
∫Ωμ-1Δg00・Δz d s=∫Ωμ-1Δ(e00-y0(T))・Δz d s-∫Ωεp0t(0)z d s+
∫Γεμ-1p0(0)z dΓ, Πz∈V0(17)
g10=p0(0)+e01-y0t(T)(18)最后置w0=g0.
初始化后,对于n≥0,假定e n、g n和w n已知,计算其更新值e n+1、g n+1和w n+1如下:
第一步:下降方向
先求正向波动方程(9),即
{ y n(0)=w n0, y n t(0)=w n1}→{ y n(T), y n t(T)}
接着求负向波动方程(16),即
{ p n(T), p n t(T)}{w n0,w n1}→{ p n(0), p n t(0)}, Πz∈V0
再根据(11)式,定义 g n={ g n0, g n1}∈E0,并由如下公式求出 g n0和 g n1,即
∫Ωμ-1Δ g n0・Δz d s=∫Ωμ-1Δ(w n0- y n(T))・Δz d s-∫Ωε p n t(0)z d s+
∫Γεμ-1 p n(0)z dΓ, Πz∈V0(19)
g n0= p n(0)+w n1- y n t(T)(20)然后计算
ρ
=∫Ωμ-1|Δg n0|2+ε|g n1|2d s∫Ω(μ-1Δ g n0・Δw n0+ε g n1w n1)d s(21) n
e n+1=e n-ρn w n(22)
g n+1=g n-ρn g n(23)
第二步:收敛判定与更新下降方向
如果满足给定的精度,那么取解为　e =e n +1,否则计算
γn =∫Ωμ-1|Δg n +10|2+ε|g n +11|2d s ∫Ω
μ-1|Δg n 0|2+ε|g n 1|2d s (24)w n +1=g n +1+γn w n
(25)取n =n +1,返回到第一步。

可以看出,每次迭代求解2个波动方程和1个椭圆型的预处理方程(式(19)),这些方程的求解都是经典的,已经有高效的解法存在,求解时可利用这些优势。

3　有限元/有限差离散格式简介
上述基于精确控制所描述的方法,由于控制迭代过程中所处理的方程都是经典的,这里只作简要说明。

以求解正向波动方程(9)为例,不失一般性,令ε=μ=1,那么问题(9)对应的变分公式可写为
∫Ωy tt z d x +∫ΩΔy ・Δz d x +∫Γ9y 9t z d Γ=0,　Πz ∈V 0y =g ,
σ上(26)
本文取平面入射波电磁场u inc =Re [e -i kt e i k (x cos β+y sin β)],其中i =-1,k 为波数,β为入射角,那么上式g =-u inc .采用时间有限差,空间有限元离散求解。

时间方向采用显式二阶中心差分格式,取Δt =T/N 为时间步长,那么式(26)的半离散格式加上问题(9)的起始条件可一并写为1Δt 2∫Ω
(y n +1-2y n +y n -1)z d x +∫ΩΔy n ・Δz d x +　12Δt ∫Γ(y n +1
-y n -1)z d Γ=0, Πz ∈V 0
y n =g (t n ),
σ上y 0
=ν0,　y 1-y -12Δt =ν1(27)
这一格式与分片线性有限元近似相结合,在数值积分中,采用了质量集中算子。

计算网格与文
[6]不同,本文采用了波阵面推进生成方法[7]生成非结构三角形网格。

由于考虑整个迭代格
式需要满足稳定性条件Δt ≤Ch ,其中C 为常数,h 为网格大小的测度。

为了获得精确解,一般要求h 至少比波长少10%[3].从这一条件可以看出,网格的质量会影响迭代时间步长,进而影响计算机时。

为此,本文对生成的网格进行了后处理,对个别质量差的网格采用了扣除等局部优化,整个网格进行了光顺。

4　算例与讨论
411　NACA 0012翼型电磁散射特性计算
本文用精确控制法,首先考核计算了文[4]的NACA 0012翼型算例。

该算例介质常数取为ε=2,μ=111,零攻角入射波频率f =112GHz ,翼型弦长取为4倍波长。

图1给出了该算例的电磁散射场。

图2给出了对应的雷达散射截面,并与文[4]的结果作了比较,两者吻合得
714　第4期陈红全等:用精确控制法计算复杂区域的二维电磁散射场
很好(雷达散射截面算法见文[4])。

　图1　电磁散射场等值线Fig.1
　Contours of the Scattered field 图2　双站雷达散射截面分布比较Fig.2　Comparison of RCS
412　主体带鼓包及舱体干扰的电磁散射场计算
本文以NACA 0012翼型模拟飞行器主体,在其上加鼓包或开挖弹舱形舱体建立计算模型(可从图3和图5中看出)。

计算时介质常数改取为ε=1,μ=1,而其它参数同上。

计算的电磁散射场以及对应的雷达散射截面分别见图3～图6。

其中雷达散射截面结果与无鼓包无舱体的主体结果作了比较,显示出加鼓包或开挖舱体对电磁散射特性的影响。

　图3　电磁散射场等值线Fig.3　Contours of the Scattered field 图4　双站雷达散射截面分布比较Fig.4　Comparison of RCS
413　主体带外挂(双体)的电磁散射场计算
本文还模拟了主体带外挂的情形,外挂由缩小的NACA 0012翼型模拟。

计算参数同§412算例。

图7给出了计算网格,从中可以看出主体带外挂的算例模型以及非结构三角形网格。

本算例存在主体与外挂两个散射体之间的干扰,一般的渐近迭代法[4]会遇到收敛的困难。

图8给出了本算例控制函数的收敛历程,显示出精确控制法具有很好的收敛特性。

计算的电磁散射场见图9,雷达散射截面与对应主体无外挂结果的比较见图10。

从图中可以看出,外挂对主体电磁散射特性有明显的影响。

814计 算 物 理第17卷　
图5　电磁散射场等值线Fig.5　
Contours of the Scattered filed 图6　双站雷达散射截面分布比较Fig.6　Comparison of RCS
图7　非结构网格示意图
Fig.7　Unstructure grid 图8　收敛历程Fig.8　Convergence history
　图9　电磁散射等值线Fig.9　Contours of the Scattered field 图10　双站雷达散射截面分布比较Fig.10　Comparison of RCS
9
14　第4期陈红全等:用精确控制法计算复杂区域的二维电磁散射场
024计 算 物 理第17卷　5　结论
本文介绍了一种用于计算电磁散射场的精确控制法。

由于方法是在非结构三角形网格上实施的,因此,可用本方法分析复杂的实用气动外形电磁散射特性。

算例表明,本方法可模拟多部件干扰的电磁散射场,具有很强的适应性。

[参　考　文　献]
[1]　方宝瑞主编.飞机气动布局设计[M].北京:航空工业出版社,1997:1226～1266.
[2]　Vineet Ahuja and Lyle N Long.A Parallel Finite2Volume Runge2Kutta Algorithm for Electromagnetic Scat2
tering[J].J Com p Phys,1997,137:299～320.
[3]　Bristeau M O,Dean E J,G lowinski R,Kwok V,and periaux J.Application of Exact Controllability to the
Computation of Scatering Waves[M].Control Problems in siecka I,MOrton B,eds,Birkhaus2 er,1995:17～41.
[4]　Kwok V.Exact Controllability and Domain Decomposition Methods for Numerical S olutions of Electromag2
netic Time Harmonic Equations in Heterogeneous Medium[Ph D thesis][D].University Paris6,1995:87～245.
[5]　Bristeau M O,G lowinski R and Periaux J.Waves Scatterin g Using Exact Controllability[C].Numerical
Methods in Engineering‘96.Desideri,Le Tallec,Onate,Periaux,Stein eds.,1996:97～103.
[6]　黄明恪.用非结构网格与欧拉方程计算复杂区域的二维绕流[J].计算物理,1994,11(4):467～471.
[7]　Lohner,Parikh.G eneration of32D Unstructure Grids by the Advancing Front Methods[C],AIAA26th
Aerospace Sciences meeting,AIAA　8820515,1988.
COMPUTATIONS OF SCATTERING WAVES IN22D
COMPL ICATED FIE LDS BY USING EXACT
CONTROLLABIL ITY APPROACH
CHEN Hong2quan,　HUAN G Ming2ke
(N anjing U niversity of Aeronautics and Ast ronautics,N anjing　210016,P R China)
[Abstract]　An exact controllability approach is used for the calcuation of time periodic solutions of the wave equa2 tion,which is proved to be equivalent to the Maxwell equation in two dimensions for the T.M.mode.The numeri2 cal results of scattered fields and Radar cross sections(RCS)are presented for complicated aerodynamic shapes with cavity,external obstacle or a slender body.The method presented has the ability to deal with more com plicated aerodynamic shapes due to the use of unstructure triangulations.
[K ey w ords]　exact controllability;Maxwell equation;aerodynamic shape。