0 引言

飞行器气动/隐身优化设计需要发展快速高效、高精度电磁数值方法，但是复杂外形目标的电磁散射包含了镜面散射、边缘绕射、爬行波等多种复杂现象，准确模拟非常困难。

高频渐进等传统工程方法虽然快捷，能捕获主要散射特征，但误差较大，不能满足隐身设计的精细要求。随着计算机能力的提高，国内外高精度数值方法获得巨大发展，主要有计算微分方程的时域有限差分FDTD^[1-3]方法，和计算积分方程的多层快多极子方法^[4-7]。其中直接求解电磁学麦克斯韦方程组的时域微分方法，最早来源于Yee K S^[1]提出的FDTD，它基于蛙跳格式，电场和磁场在网格交错配置以弥补格式黏性，有易编程和整个散射空间可统一计算的优点，但由于使用笛卡儿直角网格，模拟物面存在阶梯效应。20世纪90年代以来，基于流体力学的欧拉方程组和电磁学的麦克斯韦方程组具有共同的双曲型数学特征，人们开始在计算电磁学(CEM)引入计算流体力学(CFD)技术^[8-15]，比较典型的是采用贴体坐标系的时域有限体积法(FVTD)。

FVTD解算器有以下特点：直接求解麦克斯韦方程组；在宽频入射波的一个计算过程中可得到多个频率的结果；能达到接近矩量法(MOM)的数值计算精度；存储仅与未知量数目同数量级，能不需特殊处理地模拟散射、腔激励等复杂现象。但因涉及波在空间长时间传播计算，以及抑制耗散和色散要有一定的网格密度，网格规模会随频率增加。

实际飞机等目标在高频区为电大尺寸，例如B2隐形轰炸机，长度和翼展达数十米，对应高频区(X波段)波长最大到1000量级，每个波长取10~20个网格点，则整个规模高达10亿量级。超大计算规模成为制约气动/隐身优化的瓶颈。目前解决途径是在高性能计算机集群上，采用大规模并行计算技术。

为解决电大目标高频电磁散射大规模计算，本文通过空间网格块分割、负载平衡、多进程MPI通信等，建立了多块结构网格FVTD并行程序pmbRCS3d，在与金属球RCS解析解验证对比基础上，计算了电大尺寸飞翼外形L波段双站电磁散射场和雷达截面(RCS)，结果表明了并行计算程序的正确性和鲁棒性，利于向更高频段(X波段)扩展。

1 数值方法 1.1 控制方程

麦克斯韦方程组的两个时变旋度微分方程，法拉第电磁感应定律和安培环路定理，在自由空间、直角坐标系下可表示为守恒形式：

$ \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial {\mathit{\boldsymbol{F}}_x}}}{{\partial x}} + \frac{{\partial {\mathit{\boldsymbol{F}}_y}}}{{\partial y}} + \frac{{\partial {\mathit{\boldsymbol{F}}_z}}}{{\partial z}} = 0 $

(1)

$ \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} {{B_x}}\\ {{B_y}}\\ {{B_z}}\\ {{D_x}}\\ {{D_y}}\\ {{D_z}} \end{array}} \right],\;\;\;{\mathit{\boldsymbol{F}}_x} = \left[ {\begin{array}{*{20}{c}} 0\\ { - {E_z}}\\ {{E_y}}\\ 0\\ {{H_z}}\\ { - {H_y}} \end{array}} \right], $

$ {\mathit{\boldsymbol{F}}_y} = \left[ {\begin{array}{*{20}{c}} {{E_z}}\\ 0\\ { - {E_x}}\\ { - {H_z}}\\ 0\\ {{H_x}} \end{array}} \right],\;\;\;{\mathit{\boldsymbol{F}}_z} = \left[ {\begin{array}{*{20}{c}} { - {E_y}}\\ {{E_x}}\\ 0\\ {{H_y}}\\ { - {H_x}}\\ 0 \end{array}} \right] $

式中: B是磁感应强度矢量，H是磁场强度，D是电位移矢量，E是电场强度矢量。

对于复杂外形物体，经坐标变换，

$ k = k\left( {x,y,z} \right),\;\;\;\;k = \xi ,\eta ,\zeta $

(2)

得到曲线坐标系下麦克斯韦方程组守恒形：

$ \frac{{\partial \mathit{\boldsymbol{\hat Q}}}}{{\partial t}} + \frac{{\partial {{\mathit{\boldsymbol{\hat F}}}_\xi }}}{{\partial \xi }} + \frac{{\partial {{\mathit{\boldsymbol{\hat F}}}_\eta }}}{{\partial \eta }} + \frac{{\partial {{\mathit{\boldsymbol{\hat F}}}_\zeta }}}{{\partial \zeta }} = 0 $

(3)

其中, $ \mathit{\boldsymbol{\hat Q}} = V\mathit{\boldsymbol{Q}} $，

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat F}}}_k} = V\left( {{k_x}{\mathit{\boldsymbol{F}}_x} + {k_y}{\mathit{\boldsymbol{F}}_y} + {k_z}{\mathit{\boldsymbol{F}}_z}} \right),\\ V = \left| {\frac{{\partial \left( {x,y,z} \right)}}{{\partial \left( {\xi ,\eta ,\zeta } \right)}}} \right| \end{array} $

(4)

式中，V是坐标变换的雅可比矩阵行列式值。

1.2 积分方法

采用散射场变量既可保持守恒形式不变，又能避免传播入射场的数值计算误差，入射场解析给出。有限体积法能保持整个网格空间的通量守恒，对式(4)在网格单元作时、空积分，得到数值离散形式：

$ \frac{{\partial \mathit{\boldsymbol{\bar Q}}}}{{\partial t}} + \frac{{\partial {{\mathit{\boldsymbol{\bar F}}}_\xi }}}{{\partial \xi }} + \frac{{\partial {{\mathit{\boldsymbol{\bar F}}}_\eta }}}{{\partial \eta }} + \frac{{\partial {{\mathit{\boldsymbol{\bar F}}}_\zeta }}}{{\partial \zeta }} = 0 $

(5)

Q为网格单元内变量的平均值，F_k(k=ξ, η, ζ)分别代表单元分界处的流通量。

有限体积法的关键是精确模拟原变量Q在单元分界面的状态变量，以得到相应精确的分界面流通量F，本文F采用STEGER-WARMING分裂：

$ {\mathit{\boldsymbol{F}}_k} = \mathit{\boldsymbol{F}}_k^ + + \mathit{\boldsymbol{F}}_k^ - $

(6)

$ \mathit{\boldsymbol{F}}_k^ + = \mathit{\boldsymbol{S}}\left( {{\mathit{\boldsymbol{Q}}_L}} \right){\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^ + }\left( {{\mathit{\boldsymbol{Q}}_L}} \right){\mathit{\boldsymbol{S}}^ - }\left( {{\mathit{\boldsymbol{Q}}_L}} \right) $

(7)

$ \mathit{\boldsymbol{F}}_k^ - = \mathit{\boldsymbol{S}}\left( {{\mathit{\boldsymbol{Q}}_R}} \right){\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^ - }\left( {{\mathit{\boldsymbol{Q}}_R}} \right){\mathit{\boldsymbol{S}}^ - }\left( {{\mathit{\boldsymbol{Q}}_R}} \right) $

(8)

S、S^-为相似矩阵，Λ⁺、Λ^-分别为正负特征值构成的对角矩阵，Q_L、Q_R代表分界面处左右状态变量，采用MUSCL格式插值获得，时间推进则采用4步龙格-库塔法。

1.3 入射电磁波波形

入射电磁波选取简谐波(A₀是波幅，c为波速)：

$ {B_y} = {A_0}\sin \left[ {k\left( {z - ct} \right)} \right] $

(9)

1.4 边界条件

完全导电壁处的边界条件：

$ \mathit{\boldsymbol{\hat n}} \times \mathit{\boldsymbol{E}} = 0 $

(10)

$ \mathit{\boldsymbol{\hat n}} \times \nabla \times \mathit{\boldsymbol{H}} = 0 $

(11)

截断的网格空间外边界相容条件：

$ {\left( {\mathit{\boldsymbol{F}}_\xi ^ - } \right)_{{\rm{exit}}}} = 0 $

(12)

1.5 近场到远场的转换

稳定的时变场在频域收敛后，采用格林定理：

$ \begin{array}{l} {\mathit{\boldsymbol{E}}^s} = \frac{{{\rm{j}}k}}{{4{\rm{ \mathsf{ π} }}}}\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\nolimits_S {\left[ {{{\left( {\frac{\mu }{\varepsilon }} \right)}^{\frac{1}{2}}}\left( {\mathit{\boldsymbol{n}} \times \mathit{\boldsymbol{H}}} \right) - \left( {\mathit{\boldsymbol{n}} \times \mathit{\boldsymbol{E}}} \right) \times \mathit{\boldsymbol{r}}-\left( {\mathit{\boldsymbol{n}} \cdot \mathit{\boldsymbol{E}}} \right)\mathit{\boldsymbol{r}}} \right] \cdot } \\ \;\;\;\;\;\;\;\;\frac{{\exp \left( {{\rm{j}}k\mathit{\boldsymbol{R}}} \right)}}{\mathit{\boldsymbol{R}}}\exp \left( {{\rm{j}}k\mathit{\boldsymbol{r}} \cdot \mathit{\boldsymbol{R'}}} \right){\rm{d}}\mathit{\boldsymbol{A}} \end{array} $

(13)

R是观察点，R′是虚拟积分面上点矢量，r是观察方向单位矢量，A是虚拟积分面，k是自由空间波数，j是虚数单位。

1.6 雷达散射截面计算

$ \sigma \mathit{\boldsymbol{ = }}\mathop {\lim }\limits_{\mathit{\boldsymbol{R}} \to \infty } 4{\rm{ \mathsf{ π} }}{\left| \mathit{\boldsymbol{R}} \right|^2}{\left| {{\mathit{\boldsymbol{E}}^s}\left( \mathit{\boldsymbol{R}} \right)/{\mathit{\boldsymbol{E}}^i}\left( \mathit{\boldsymbol{R}} \right)} \right|^2} $

(14)

1.7 并行计算

电大尺寸电磁散射伴随大规模三维网格，须用多进程并行计算，包括：(1)多进程分割和负载平衡；(2)程序并行化处理，采用MPI接口进行通信。

2 计算结果与分析 2.1 并行化多块FVTD程序pmbRCS3d验证

多块网格FVTD程序mbRCS3d，并行化扩充为pmbRCS3d。首先选择带解析解的金属球(ka=10，其中k为波数，a是半径)为验证目标。

图 1是网格块经8个进程分割后表面网格。图 2是表面散射场等值线云图。图 3是表面诱导电流分布，结果显示前向散射最强。图 4是双站RCS数值与精确解析级数解比较，吻合很好，说明整个并行FVTD算法具有很好的精度。

2.2 电大尺寸实际外形RCS计算

针对飞翼外形(类似X-47)的L波段(f=1.2GHz)电磁散射问题，长度11.64m，翼展18.9m，长度对应电尺寸为47个波长，翼展对应76个波长。采用36块网格离散整个计算空间，网格点16 426 816个。网格密度根据电磁波频率选取，要求每个坐标方向上一个波长(250mm)内至少布15~20个网格点，散射体表面网格加密，远离物面区域适当疏松，以便减少计算量，图 5是重新分块后的物面网格边界。在集群上使用16个进程计算，耗时420h。

该大网格量电磁场并行FVTD计算，得到稳定振荡的简谐场和逐步收敛的频域散射场，从图 6和图 7散射场云图能看出散射波耦合强烈区域。图 8的结果(其中φ为与+X轴所成方位角)表明：前向RCS较大，后向RCS较小在－35dB附近，有4个较强散射方位(分别对应4个垂直于机身前缘和机翼前缘的方位)。

3 结论

通过散射电磁场的FVTD解算器和基于MPI的并行编程，一定程度解决了电大尺寸目标高频散射大计算量难题。计算和验证表明：该系统能提供高精度的电磁场和RCS分布，为气动/隐身优化提供精细工具, 并有利于向更高频段扩展。