钢制船舶会产生干扰磁场，其中包含固定磁场和感应磁场；又因材料的导电特性，船舶在横摇、纵摇、艏摇运动时会激发涡流磁场^[1-2]。由于船舶消磁^[3]、船载磁传感器、海洋磁力测量^[4]等不同应用需求，工程上采用低磁或者无磁材料制造船舶(如扫雷舰艇)，固定磁场和感应磁场均不予以考虑。然而由于气候状态、自身操纵需要等原因，船舶不可避免地会激发涡流磁场，该干扰磁场已成为制约船舶低磁化技术和海洋磁力测量技术发展的主要因素之一。非磁性船舶和航空器涡流磁场产生机理是一致的，但因难以控制船舶运动状态以获取磁场数据样本，因此基于Lawson模型的航空器涡流磁场相关理论不适用于非磁性船舶涡流磁场研究^[5]，而选择磁变模拟方法^[6-7]研究。磁变模拟方法最初主要用于涡流磁场的等效实验测量^[8]，延伸至基于磁变模拟法的涡流磁场数值建模。吕正石等^[9]将船舶等效为椭球壳体以开展数值计算，但因磁通量变化形成的环流难以等效，故该方法值得商榷；Le等^[10]采用基于面单元的有限元法计算了薄板转动时磁异常，计算结果被商业有限元软件FLUX3D^[11]仿真结果验证；有限元软件OPREA也对涡流磁场进行了仿真^[12]，但没有公开文献说明其具体的涡流数值建模方法。

鉴于积分方程法不需对整个计算区域剖分而特别适用于磁场开域问题计算^[13-15]，目前没有文献公开报道采用积分方程法计算船舶涡流磁场；又因船舶外壳体可视为薄壳体的特点，本文采用三角形壳单元对简易船舶模型进行剖分，提出基于磁变模拟方法和积分方程法对涡流磁场开展数值计算研究。

1 磁变模拟方法

涡流磁场主要在船舶横摇时产生^[3]，在图 1所示的坐标系下，以船舶横摇为例，当空间均匀磁场为H_xE=40 A/m、H_yE=0、H_zE=-30 A/m，且船舶发生横摇(以右舷向上时为正)时，H_xE与船舶艏艉平行，转动角度(初始相位为0)可表示为:

$ \alpha = {\alpha _0}\sin \left( {\omega t} \right) $

(1)

式中:α₀为转动角幅值；ω为角频率。由于考虑非磁性情况，H_xE不予以考虑。此时磁变模拟方程为: H_x=0

$ \left\{ \begin{array}{l} {H_x} = 0\\ {H_y} = {H_{z{\rm{E}}}}\sin \alpha \\ {H_z} = {H_{z{\rm{E}}}}\cos \alpha \end{array} \right. $

(2)

同理可推导船舶在其他运动情形下的磁变模拟表达式。由式(2)可知H_y、H_z对船舶激励涡流磁场的计算是暂态问题，采用泰勒级数分解可将涡流磁场暂态求解问题转化为稳态问题^{[10, 12]}，使得求解效率更高。具体表达式为:

$ \left\{ \begin{array}{l} \cos \left( {{\alpha _0}\sin \left( {\omega t} \right)} \right) = {a_0} + \sum\limits_{i = 1}^\infty {{a_i}\cos \left( {2{\rm{i}}\omega t} \right)} \\ \sin \left( {{\alpha _0}\sin \left( {\omega t} \right)} \right) = \sum\limits_{i = 1}^\infty {{b_i}\sin \left( {2{\rm{i}}\omega t} \right)} \end{array} \right. $

(3)

分解时，可根据具体磁变模拟函数和精度需求选择分解阶数。在满足精度要求的前提下，分解阶数越小意味着涡流磁场求解的效率越高。船舶在不同海况下运动时^[1-2]，α₀≤30°且f≤1 Hz(频率通常在0.2 Hz左右)。当H_zE=-30 A/m，f=0.2 Hz时，分别将弧度值α₀分解为π/18、π/9和5π/36这3种情况。

图 2为上述3种情况1阶和5阶泰勒级数解和实际值的比较情况，表明3种情况下采用1阶泰勒级数可较好地逼近原激励磁场函数，式(2)的1阶级数分解为

$ \left\{ \begin{array}{l} {H_y} = {H_{z{\rm{E}}}}\sin \theta \cong {H_{z{\rm{E}}}}{b_1}\sin \left( {2{\rm{ \mathsf{ π} }}ft} \right)\\ {H_z} = {H_{z{\rm{E}}}}\cos \theta \cong {H_{z{\rm{E}}}}\left( {{a_0} + {a_1}\cos \left( {4{\rm{ \mathsf{ π} }}ft} \right)} \right) \end{array} \right. $

(4)

式中:a₀=1-α₀²/4、a₁=α₀²/4-α₀⁴/48、b₁=α₀-α₀³/3；对于非磁性船舶，式(3)中的a₀H_zE不存在激励响应，因此式(3)可写为

$ \left\{ \begin{array}{l} {H_y} = {H_{z{\rm{E}}}}\sin \theta \cong {H_{z{\rm{E}}}}{b_1}\sin \left( {2{\rm{ \mathsf{ π} }}ft} \right)\\ {H_z} = {H_{z{\rm{E}}}}\cos \theta \cong {a_1}\cos \left( {4{\rm{ \mathsf{ π} }}ft} \right){H_{z{\rm{E}}}} \end{array} \right. $

(5)

分析式(4)可知:保持f不变情况下，在求解此时涡流磁感应强度可直接由已计算α₀情形推算得出。

2 涡流磁场的积分方程法

级数分解后，非磁性船舶涡流磁场问题的求解需解决正、余弦形式激励磁场产生的涡流磁场(稳态涡流问题求解)。文献[13]基于积分方程法提出了一种适用于壳体涡流磁场数值计算的方法，只需对薄壳体表面进行剖分且对不同集肤深度δ和厚度e关系($ \delta \ll e, \delta \approx e, \delta \gg e$)适用。针对船舶可近似为薄壳体的几何结构、材料特点以及低频(0~1 Hz)运动特性，其集肤深度和壳体厚度的关系为$ \delta \gg e$的情况，故采用文献[13]方法计算涡流磁场。

2.1 数值计算方法

采用三角形壳单元对薄壳体进行建模并以各节点处位势差Δϕ为未知量，得到离散方程组:

$ \mathit{\boldsymbol{A}} \cdot \Delta \phi - {\mathit{\boldsymbol{B}}_{\rm{s}}} \cdot \Delta \phi = {\mathit{\boldsymbol{h}}_{0n}} $

(6)

其中

$ \mathit{\boldsymbol{A}}\left( {i,k} \right) = \left( {\alpha + \beta } \right)\int\limits_\mathit{\Gamma } {{\rm{grad}}{w_i} \cdot {\rm{grad}}{w_k}{\rm{d}}\mathit{\Gamma }} $

(7)

$ {\mathit{\boldsymbol{B}}_{\rm{s}}}\left( {i,k} \right) = \frac{{2{\rm{j}}\omega {\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int\limits_\mathit{\Gamma } {{w_i}} \cdot \int\limits_{{\mathit{\Gamma }_k}} {\frac{{{\mathit{\boldsymbol{n}}_k} \times {\rm{grad}}{w_k} \times {\mathit{\boldsymbol{r}}_{ki}}}}{{r_{ki}^3}}{\rm{d}}\mathit{\Gamma }} \cdot {\mathit{\boldsymbol{n}}_i}{\rm{d}}\mathit{\Gamma } $

(8)

$ {\mathit{\boldsymbol{h}}_{{\rm{0}}n}}\left( i \right) = - 2{\rm{j}}\omega {\mu _0}\int\limits_\mathit{\Gamma } {{w_i}} \cdot {\mathit{\boldsymbol{H}}_0} \cdot {\mathit{\boldsymbol{n}}_i}{\rm{d}}\mathit{\Gamma } $

(9)

式中:系数矩阵 A、B_s为剖分单元之间相互作用的系数矩阵；α=a/(σtanh(ae))，β=a/(σsinh(ae))，σ为材料电导率，a=(1+j)/δ；μ₀为真空磁导率; ω为角频率；H₀为场源；n为单位外法线向量；w为单元形状函数(可为一阶、高阶形状函数)；r_ki为源点指向场点的矢径。

在求解节点未知量Δϕ后，任意场点的磁场强度为:

$ \mathit{\boldsymbol{H}} = {\mathit{\boldsymbol{H}}_0} + {\mathit{\boldsymbol{H}}_r} $

(10)

式中$ {\mathit{\boldsymbol{H}}_r} = \left( {1/4{\rm{ \mathsf{ π} }}} \right)\int\limits_\mathit{\Gamma } { - \mathit{\boldsymbol{n}} \times {\rm{grad}}\Delta \times \phi \times \mathit{\boldsymbol{r}}} /{r^3}{\rm{d}}\mathit{\Gamma }$, r为源点到均点的距离。

2.2 球壳体解析解计算案例

根据式(6)~(10)编写了非磁性壳体涡流磁场积分方程法数值计算的程序包。为验证算法准确性，计算如图 3所示，外半径b为100 mm，半径a为98 mm，相对磁导率为1的导电球壳体置于${\mathit{H}_{0z}} = \sqrt 2 \cos \left( {2{\rm{ \mathsf{ π} }}\mathit{ft}} \right)$A/m激励中产生涡流磁场。当剖分三角形壳单元数目为2 048、节点数目为1 026时，采用式(11)进行解析。分别计算了三种频率和电导率参数情况下路径AB(A(0.15, 0, 0)，B(0.15, 0, 0.1)，m)上等间距分布及26个场点磁场强度总量，如图 4所示为计算结果。以各场点磁场强度解析解作为参考值，f=1 Hz时, 最大相对误差为0.05%；f=10 Hz, 最大相对误差为0.44%；f=50 Hz, 最大相对误差为0.43%。

$ \left\{ \begin{array}{l} {B_\theta } = - {B_0}\left( {1 - D/\left( {2{r^3}} \right)} \right)\sin \theta \\ {B_r} = {B_0}\left( {1 + D/{r^3}} \right)\cos \theta \\ {B_\varphi } = 0 \end{array} \right. $

(11)

其中

$ F/E = \frac{{ - \left( {\alpha _1^2 + 3} \right){{\rm{I}}_{1/2}}\left( {{\alpha _1}} \right) + 3{\alpha _1}{{\rm{I}}_{ - 1/2}}\left( {{\alpha _1}} \right)}}{{\left( {\alpha _1^2 + 3} \right){{\rm{I}}_{ - 1/2}}\left( {{\alpha _1}} \right) - 3{\alpha _1}{{\rm{I}}_{ - 1/2}}\left( {{\alpha _1}} \right)}}, $

$ E = 3{b^{3/2}}/\left( {{\beta _1}{{\rm{I}}_{1/2}}\left( {{\beta _1}} \right) + \left( {F/E} \right){\beta _1}{{\rm{I}}_{ - 1/2}}\left( {{\beta _1}} \right)} \right), $

$ \begin{array}{*{20}{c}} {D = \left( {1/3} \right){b^{3/2}}E\left[ {3{{\rm{I}}_{ - 1/2}}\left( {{\beta _1}} \right) - \left( {{\beta _1} + 3/{\beta _1}} \right){{\rm{I}}_{1/2}}\left( {{\beta _1}} \right) + } \right.}\\ {\left. {\left( {F/E} \right)\left[ {3{{\rm{I}}_{1/2}}\left( {{\beta _1}} \right) + \left( {{\beta _1} + 3/{\beta _1}} \right){{\rm{I}}_{ - 1/2}}\left( {{\beta _1}} \right)} \right]} \right]} \end{array} $

式中:I_n(·)为n阶第一类修正柱贝塞尔函数；${\alpha _1} = \sqrt {{\rm{j}}pa} , {\beta _1} = \sqrt {{\rm{j}}pb} , p = {\mu _0}\omega /\gamma $。

3 船舶涡流磁场数值计算及分析

本文所提出的船舶涡流磁场求解方案对船舶横摇、纵摇和艏摇涡流磁场求解均适用。文中以船舶横摇运动为例，计算导电非磁性简单船舶模型的涡流磁场，并对数值计算结果进行验证与分析。

3.1 船舶涡流磁场求解方案

船舶涡流磁场数值计算主要步骤为:

1) 根据船舶运动状态，推导磁变模拟方程；

2) 对步骤1)中所得磁变模拟函数进行泰勒级数分解；

3) 对非磁性船舶建模，所有导电材料结构须在建模过程中予以考虑；

4) 采用三角形壳单元对步骤3)中所建的几何模型剖分；

5) 针对剖分三角形单元，基于磁场积分方程法得到离散方程组式(6)；

6) 基于离散方程组式(6)，求解步骤2)中泰勒级数分解所得各函数激励下的涡流磁场H_r'；

7) 求步骤6)中各激励产生涡流磁场H_r'矢量和，通过式(10)求得任意场点的磁场强度H。

3.2 船舶模型涡流磁场计算案例

按照船舶涡流磁场求解方案，计算如图 5所示(A(-25, -5, 0) m, B(25, -5, 0) m, C(25, 5, 0) m, D(-25, 5, 0) m), 厚度为12 mm、材料电导率为3.03×10⁷ S/m简易铝质船模在横摇时产生的涡流磁场。空间均匀磁场为H_xE=40 A/m、H_yE=0、H_zE=-30 A/m，且船舶发生横摇时H_xE与船舶艏艉平行。横摇频率为f=0.2 Hz，α₀=π/9(初始相角为0)，此时具体的磁变模拟函数为:

$ \left\{ \begin{array}{l} {H_x} = 0\\ {H_y} = {H_{z{\rm{E}}}}\sin \left( {{\rm{ \mathsf{ π} }}/9\sin \left( {0.4{\rm{ \mathsf{ π} }}t} \right)} \right)\\ {H_z} = {H_{z{\rm{E}}}}\cos \left( {{\rm{ \mathsf{ π} }}/9\sin \left( {0.4{\rm{ \mathsf{ π} }}t} \right)} \right) \end{array} \right. $

(12)

对于式(12)的泰勒级数分解，采用一阶泰勒级数分解:

$ \left\{ \begin{array}{l} {H_x} = 0\\ {H_y} = {H_{z{\rm{E}}}}{b_1}\sin \left( {0.4{\rm{ \mathsf{ π} }}t} \right)\\ {H_z} = {H_{z{\rm{E}}}}\left( {{a_0} + {a_1}\cos \left( {0.8{\rm{ \mathsf{ π} }}t} \right)} \right) \end{array} \right. $

(13)

式中:a₀=0.969 8, a₁=0.030 2, b₁=0.335 0。

上述参数情况下H_y、H_z一阶级数泰勒分解情况如图 2所示。以真实值为参考，H_y分解的平均相对误差为1.65%，H_z分解的平均相对误差为0.01%。简易船舶模型为非磁性材料制作而成，恒定磁场a₀H_0z不会产生激励磁场，此时涡流磁场计算时可认为式(13)中H_z=H_zEa₁cos(0.8πt)。利用三角形壳单元对船舶外壳进行剖分，剖分单元数目11 576，节点数目为5 790。采用积分方程法分别计算H_y、H_z激励产生的涡流磁场，同时采用商业有限元软件Flux3D对计算结果进行验证。以某一场点的数值结果分析非磁性船舶涡流磁场的特性。如图 5所示，计算路径P₁P₂(P₁(-60, 0, -15)，P₂(60, 0, -15)，m)上等间距分布51个场点的磁感应强度。

3.2.1 H_y =H_zEb₁sin(0.4πt)激励的涡流磁场

图 6为分别采用积分方程法和商业有限元软Flux3D计算t=9.25 s时H_y产生的磁感应强度y分量；IEMB_yT、FluxB_yT分别表示积分方程法和有限元法磁感应强度数值结果y分量总量(包含激励磁场)；图 7为采用两种不同方法分别计算t=9.25 s时H_y产生的涡流磁感应强度三分量(不包含激励磁场)，以有限元软件计算结果作为参考，积分方程法所计算的涡流磁感应强度y分量最大相对误差为1.25%；x和z分量两种方法计算结果都接近于0，x和z分量的最大绝对误差分别为0.26、0.10 nT。

3.2.2 H_z=H_zEa₁cos(0.8πt)激励的涡流磁场

图 8为分别采用积分方程法和商业有限元软件Flux3D计算t=4.875 s时H_z产生的磁感应强度z分量；IEMB_zT、FluxB_zT分别表示磁感应强度z分量总量(包含激励磁场)；图 9为采用两种不同方法分别计算t=4.875 s时H_z产生的涡流磁感应强度三分量(不包含激励磁场)，以有限元软件计算结果作为参考，积分方程法所计算的涡流磁感应强度x分量最大相对误差为1.29%;z分量最大相对误差为1.12%；y分量两种方法计算结果都接近于0，y分量的最大绝对误差为0.49 nT。

3.3 特殊场点涡流磁场分析与讨论

有限元软件计算结果的有效验证了积分方程法计算H_y、H_z激励产生涡流磁感应强度的准确性。下面以场点P(0, 0, -15 m)为例对非磁性船舶涡流磁感应强度的特性进行分析与讨论(其他场点情况分析可采用类似方法)。

由于H_y、H_z激励磁场的周期分别为2.5和5 s，在此计算5~10 s时的涡流磁感应强度。由积分方程法计算结果可得点P处H_y激励涡流磁场强度x、y、z三分量的分量幅值分别为1.73×10^-4、5.00×10^-2、2.56×10^-10 A/m；相位分别为-2.04、-1.82、-1.03 rad。点P处H_z激励涡流磁场强度x、y、z三分量的分量幅值分别为2.85×10^-6、2.91×10^-10、2.47×10^-2 A/m；相位分别为2.29、1.42、1.16 rad。图 10为横摇运动情况下激励的涡流磁感应强度三分量，分析可得y和z分量的1周期内最大幅值分别为62.75、31.03 nT。图 11为涡流磁感应强度总量，在1周期内最大幅值为63.13 nT。为比较分析，采用文献[13]中的静磁场积分方程法计算了相对磁导率为60(激励磁场为H_zE=-30 A/m)时场点P的感应磁感应强度总量，其值为213.99 nT；相对磁导率为150时场点P的感应磁感应强度总量为464.86 nT。

4 结论

1) 3种不同参数情况时球壳体解析算例最大相对误差均在0.50%以内；以有限元结果为参照，各分量最大相对误差均不超过1.50%;表明了积分方程法涡流磁场数值计算的准确性。

2) 场点涡流磁感应强度三分量和总量分析表明涡流磁场是非磁性或者低磁性船舶周围不可忽略的干扰磁场，扫雷舰艇应当重视涡流磁场的消磁。

3) 建立的涡流磁场数值计算模型可用于获取任意场点的涡流磁感应强度从而得到船舶涡流磁特征。

目前计算了船舶横摇时的涡流磁场，后续应该从以下几个方面开展进一步研究:1)在对船舶做更为复杂运动(横摇、纵摇、艏摇同时存在时)的涡流磁场计算时，其磁变模拟函数更为复杂，采用级数分解的思路同样适用于更为复杂的运动情况；2)目前所计算的船舶模型涡流磁场的场点离模型最下方距离为10 m，计算离船舶模型更近场点的涡流磁场或者实际中几何形状复杂船舶的涡流磁场可能需要对船舶作更为精细的剖分，已编写的程序包在剖分单元数目较大时效率较低，应当研究快速计算方法(如快速多极子)加速涡流磁场的求解；3)将所研究非磁性船舶涡流磁场的思路延拓于磁性船舶涡流磁场的研究。