计算泄漏同轴线的传播与辐射特性

1、计算泄漏同轴线的传播与辐射特性一、题目A22 计算如图所示泄漏同轴线的传播与辐射特性同轴线外导体沿轴向连续开槽二、 FDTD 方法分析泄漏同轴线FDTD 的主要思想是把Maxwell 方程在空间、时间上离散化，用差分方程代替一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。2.1 泄漏同轴电缆简介泄漏同轴电缆（Leaky Coaxial Cable ）是遵循特定的电磁场理论，沿着同轴电缆的外导体周期性或非周期性开槽而形成的。信号在该电缆中传输时，能把电磁能量的一部分按要求从特殊开槽口以电磁波的形式辐射到周围外部空间，因此泄漏电缆能实现与外部空间的全方位双工通信。利用泄漏同轴电缆既具有传

2、输线的特性又具有无线电发射天线的性质，人们近年来广泛地利用它来解决无线电无法传输信息区域的信息传输问题。本文研究的是同轴线外导体沿轴向连续开槽的问题。2.2 Maxwell 方程在介质基片和空气中，Maxwell 旋度方程为：HDJ （1）tEBJ m （ 2）t其中， E 为电场强度，单位为V/m, D 为电通量密度，单位为C/m2 , H 为磁场强度，单位为 Wb/ m2 ,J 为电流密度，单位为A / m2, Jm 为磁流密度，单位为V / m2 。各向同性介质中本构关系为：DE ， BH， JE ， J mm H（ 3）其中，表示介电系数，单位为F / m ,表示磁导系数，单位为H /

3、 m ，表示电导率，单位为S / m ，m 表示导磁率，单位为/ m 。和 m 分别为介质的电损耗和磁损耗。在真空中0m 08.8510 12(4)0F / m0410 7 H / m2.3 圆柱坐标系中差分格式的建立在圆柱坐标系中， （ 1）、（ 2）式可写为：1 H zHErErrztH rH zEE（ 5）zt1(rH )1 H rEzEzrrrt以及1( Ez )EH rm H rrzt( Er )EzHm H（6）zrt1(rE )1 EHzrrrrtmH z柱坐标下的Yee 元胞如下图所示：图 1柱坐标下的Yee 元胞由图 1 可见，每一个磁场分量由四个电场分量环绕；同样，每一个电

4、场分量由四个磁场分量环绕。这种电磁场分量的空间取样方式不仅符合法拉第感应定律和安培环路定律的自然结构，而且这种电磁场各分量的空间相对位置也适合麦克斯韦方程的差分计算，能够恰当的描述电磁场的传播特性。此外，电场和磁场在时间上顺序上交替抽样，抽样时间间隔彼此相差半个时间步，使麦克斯韦旋度方程离散以后构成显式差分方程，从而可以在时间上迭代求解，而不需要进行矩阵求逆运算。因而，由给定相应电磁问题的初值及边界条件，利用FDTD 方法就可以逐步地求得以后各个时刻空间电磁场的分布。取差分近似后(5) 式 Er 的 FDTD 离散式为Ern 1 (i1 , j , k) Ern (i2tH zn 1/2 (i

5、1 , j1 , k)221 , j , k)2H zn 1/2 (iErn 1 (i1 , j ,k ) Ern (i1 , j , k)2221 , j1 , k)(7)22(i1) r2(7) 式整理后得：12t12 tH zn1/2 (i1, j1, k)Hzn1/2 (i1, j1, k)n 1n(i2222Er(i , j ,k)2tEr, j ,k)t1222(ir)(8)2H n 1/2 (i1 , j1 ,k)H n1/2 (i1 , j1 ,k)22z22同理可得：12t12 tH rn1/2 (i , j1 , k1)H rn1/2 (i , j1 ,k1)En 1(i

6、, jn(i , j,k)2222,k)2tE2tz22（ 9）1111H zn 1/2 (i, j,k)H zn 1/2 (i, j, k)22r22n 112tn12 t(i1)H n 1/ 2(i1 , j, k1) (i1)H n 1/ 2(i1 , j ,k1)222222Ez(i, j ,k2)2t Ez(i , j , k2)2ti r（ 10）H n 1/ 2(i , j1, k 1) H n 1/ 2(i , j1 ,k1)r22rr22 i以上为电场的FDTD 递推公式，同样可得到磁场的FDTD 离散式为1, k1)2t Hrn1 ,k1)2 tEzn (i, j1,k1)

7、Ezn (i, j, k1)H rn 1/2 (i, jm1/2 (i , j22222m t22 2m tr（ 11）1 , k1 , k)En (i, j1)En (i, j22z112t112 tErn (i1 , j ,k 1) Ezn (i1 , j, k)Hn 1/2(imn 1/2(i22, j, k2)H, j, k)m tz22m t22 2（ 12）1) Ezn (i , j ,k 1)Ezn (i 1, j, k22 rH zn 1/2 (i1 , j1 ,k)2mt Hzn 1/2 (i222mtErn (i1, j1,k)Ern (i21)(ir22.4 数值稳定性分

8、析1 , j21 , j ,k)212t(i 1)En (i 1, j1 ,k)iE n (i, j1 , k)22, k)2t12m(ir)（ 13）2FDTD 方法是以一组有限差分方程来代替麦克斯韦旋度方程，即以差分方程组的解来代替原来电磁场偏微分方程的解。只有离散后差分方程组的解是收敛和稳定的，这种代替才是有意义的。收敛性是指当离散时间隔趋于零时，差分方程的解在空间任意一点和任意一时刻都一致地趋于原方程的解。稳定性是指寻求一种离散间隔所满足的条件，在此条件下差分方程的数值解与原方程的严格解之间的差为有界。在圆柱坐标系下，对空间步长和时间步长的要求比在直角坐标系下更为严格。沿同一径向，在角

9、度步长相同时，横向长度单元是不同的，因此，为了保证数值收敛性，要取最小的值来计算。因为不同径向距离的横向单元长度为ir，但是考虑到电场和磁场相差半个时间步长，可取最靠近轴向的网格的横向单元长度，即取 i 1/2 ，此时时间步长的取值为t1（ 14）121c222()()()rrz由（ 14）式可知，在轴向 i0, rir0 ，在计算中会出现奇点，同样，H r在轴向也存在奇点。这些奇点不是Maxwell 方程本身所固有的，而是由微分方程转变为差分方程后带来的，可以用安培定律去除，因为H d lE d S（ 15）ct s取上式的积分得：r02t ez (0, z,t )h r0 (r0 , z,

10、 t)（ 16）2式中， rr /2，将（ 16）式转换为差分格式为：ezn 1(0, j ,k1)ezn(0, j, k1)4 t hn 1/2 (1 , j, k1)（ 17）022r22有了 ez (0,j , k1 ) 后， hz (0, j1 , k1 ) 及其他格式将自动满足。2222.5 吸收边界条件：由于天线是一个无限空间的问题，而数值计算只能在有限空间进行，因此在计算空间的边界时必须设置吸收边界来模拟无限空间，使场满足辐射条件。为了让这种有限空间与无限空间等效，需对有限空间的周围边界面做特殊处理，使得向边界面进行的波在边界处保持“外向行进 ”的特征、无明显反射现象，并且不会使

11、内部空间的场产生畸变。本文采 Mur 吸收边界条件。在圆柱坐标系中，三维的二阶Mur 吸收边界条件的差分格式为：1）0：n 1 (0, j , k ) c1 n1 (1,j ,k)c2c4n (0, j1,k)c6n (0, j ,k 1)其中,1/， c014tn (0, j ,k )n (1, j , k) c3 n (0, j1,k)n (1, j 1,k)n (1, j1,k)c5n (0, j , k1)n (1, j , k1)（ 18）n (1,j , k 1)c7n 1 (0, j, k)c8n 1(1, j, k)112 ， c1 (14112 ) / c0 ，tttttc2

12、 ( 2282222 ) / c0 ， c3 c4222/ c0 ， c5 c62 / c0 ，tz2zc7 (111) / c0， c8(111) / c0 。t4 tt 2t4 tt 22）maxN0 ：n 1 (N0, j , k) c1n 1 (N01, j ,k)c2n ( N01, j, k)c4n1, j 1,k)n1,k)(N0(N0 , jc6n (N01, j ,k1)n (N0 , j , k 1)1111其中， c04tt 2， c1 (ttn (N0 , j , k) c3n (N0 1, j1,k)n (N0 , j1,k)c5n1)n1)(19)( N0 1, j

13、, k( N0 , j, kc7n 1 ( N0, j, k)c8n 1 (N01, j ,k)411) / c0 ，tt 2c2 (22 ) / c0 ， c3 c4/ c0 ， c5 c62 / c0 ，t28222222zz2c7 (1112) / c0 ， c8 (1112 ) / c0 。t4 ttt4 tt3）min0 ：n 1 (i ,0, k ) c1n1 (i ,1, k)c 2c4n (i 1,0, k )c6n (i ,0, k1)n (i ,0, k )n ( i,1, k )c3n ( i1,0, k )n (i1,1,k )n (i 1,1,k )c5n (i ,0

14、, k1)n (i ,1, k1)（ 20）n (i ,1, k1)c7n 1(i ,0, k)c8n 1 (i ,1,k )4）maxL0 ：n 1 (i ,L0 , k) c1 n 1 (i , L0 1,k ) c2n (i, L0 , k)n (i, L01,k ) c3n (i1,L0 , k)n (i 1,L01,k)c4 n ( i 1,L0 , k)n ( i 1,L0 1,k)c5n (i, L0 , k 1)n (i, L01,k1)（ 21）c6 n (i, L0 ,k 1)n ( i, L0 1,k 1) c7n 1 (i, L0 , k)c8n (i, L01,k)式

15、（ 20）和（ 21）具有相同的系数：c011112z2 ) / c0 ，tt2 ， c1 (tt2 ) / c0 ， c2 (t 22c3 (224) / c0 ， c4(224) / c0 ， c5 c62z2 / c0 ，c7 (11 2 ) / c0 ， c81。tt5） zzmin0 ：n 1 (i ,j ,0)c1n (i , j ,1)c2n (i , j ,0)n ( i, j ,1)c3n ( i1, j ,0)n (i 1, j ,1)c4n (i1, j ,0)n (i1,j ,1)c5n (i , j1,0)n (i , j1,1)（ 22）c6n (i , j1,0)

16、n (i,j1,1)c7n 1 (i , j ,0)c8n1 (i , j ,1)6） zzmaxM 0 ：n 1 (i, j ,M 0 )c1n (i, j, M 0 1)c2 n (i, j, M 0 )n ( i, j ,M 01)c3n (i1, j ,M 0 )n (i1, j, M 01)c4n(i1, j, M 0 )n(i1, j,M 01)c5n(i, j1,M0)n(i, j1,M0 1)（ 23）c6n (i, j 1,M 0 )n (i, j1,M 01)c7n 1 (i, j,M 0 )c8n 1 (i, j, M 01)式（ 22）和（ 23）具有相同的系数：c01

17、1， c1 ( 11) / c0 ， c2 ( 2) / c0 ，t zv t2t zt 2t 2222c3 (224) / c0 ， c4(24) / c0 ， c5 c6222/ c0 ，2c7 (11 2 ) / c0 ， c81。tzt2.6 激励源：采用 Gauss 脉冲激励， Ei (t )exp(t t0 ) 2/T2（ 24）式中 t0 是出现最大值的时间，T 与脉冲宽度宽度有关，应适当选择t 0 和 T ，以保证源的初始条件，即Ei (t0)0 。对式（ 24）进行傅里叶变换得Ei ( f )Texp(2T 2 f 2 )（ 25）本文取 T1.67 e 11s ， t03T

18、 。2.7 近远场外推：基本思想由于 FDTD 方法只能计算空间有限区域的电磁场，要获得远区的辐射场则必须应用等效原理，如图2所示，在源的周围引入虚拟界面S，而 S 外是自由空间。令 S 之外存在原来的场，在S 之内的场为零，且到处是自由空间，为支持这样的场，S 面上必须存在等效面电流Js 和面磁流 Jm 。（ a）原问题（ b）等效问题图 2等效原理示意图J s nH , Jm En三维自由空间的格林函数为：G( r , r')exp(jk | rr '|)（ 26）r '4| r|其中 r , r ' 分别为观察点和源点位置失量，取远区近似得： | rr &

19、#39; | rr 'er，其中 er 为 r 方向单位矢量， r| r |于是，（ 26）式近似为G ( r , r ' )exp(jkr ) exp( jk r 'er)（ 27）4rA( r )J( r ' )G( r , r ' )dV ' , F ( r )J m (r ' )G( r , r ' )dV '（ 28）A( r )exp(jkr )J( r ' )exp( jk r 'er )ds'将式（ 27）代人式（ 28）得：4rs（ 29）exp(jkr )F ( r )J m(

20、r ' )exp( jk r ' er )ds'4rsEF1AFj1( A)jA电流和磁流的辐射场为j（ 30）11HAFAj FF )(jjEj k FjAk( k A)在远区，（ 30）式中的算子可以用 (j k ) 来代替，因此：j）（31Hj kA jFk( k F )j将（ 31）式右端写成球坐标分量形式，从而可以得到：EjkFjA（ 32）EjkFjAHjkAjF（ 33）HjkAjF三维瞬态场的外推如果 FDTD 计算时入射波是时域脉冲的形式，就需要考虑如何将时域近场数据外推到计算区域以外的情况，对于远区场，由式（33）可得到：Ejk ( ZAF )（ 3

21、4）Ejk ( ZAF )下面将频域远区场外推表达式（34）通过傅里叶变换转换到时域。为此在频率中令：Wjk A( r)e x p ( jkr)') exjkp r('er ds'jk4 rJ r()s（35）jk e x p ( jkr)Jm r('Ujk F ( r) e xjkp r( ' er ds')4 rs其中， k/ c, J (r ' )J (r ' ,), J m ( r ' )J m (r ' ,) ，记其傅里叶变换为：j (r ' ,t )J ( r ' , )exp(j t

22、)ds（ 36）jm ( r ' ,t)J m( r ' , )exp( jt)ds根据傅里叶变换性质，有：tj ( r ' ,t )jJ (r ' ,)exp(jt) ds（ 37）tjm (r ' , t)jJm (r ' ,)exp(j t) dsj ( r ' ,trerr ')J (r ' ,) exp(jkr ) exp( j k r ' ) exp( j t )d以及ccs（38）jm (r ' , trer r ' )Jm (r ' , ) exp(jkr ) exp( j

23、k r ' ) exp( jt) dccsw(t )1t sj ( r ' ,terr 'r )ds'4rccc由（ 37）、（38）式得（ 36）式的傅里叶变换为：（ 39）u(t)1t sjm (r ' , terr 'r )ds'4rccc所以，将（ 34）、（ 35）式并作傅里叶变换，得到时域中的关系式为：e (t)uZwe (t )u（ 40）Zw其中， e(t ) 是电场 E 的傅里叶变换。三、 FDTD 方法计算泄漏同轴线的传播与辐射特性3.1 用 FDTD 方法求解泄漏同轴电缆传播和辐射场的步骤：1、将时域Maxwell

24、的旋度方程展成其坐标分量式，本文采用圆柱坐标系，用中心有限差分格式代替各场分量对空间、时间的微分，得到FDTD 的基本方程。2、确定空间网格基本单元尺寸r ,z ，单元尺寸大小可以相等，也可以不等，为了减小数值色散，应满足，min( r , , z) ，本文中，取 rz/100 。3t ，根据前面的讨论计算得到：t1，、选择时间步长c ( 1 ) 22(1)2() 2rrz本文近似取tr / 2c 。4、沿网格空间的外表面设置吸收边界条件。5、选用和设置吸收边界条件。6、确定运算的总时间步数，对于时域计算，应使激励源主要频谱部分的信号达到稳定状态，特别应使低频谱信号达到稳定状态。7、计算边界处的电磁场，利用等效原理和近远场外推的方法，求得辐射场。8、如图 2 所示，利用A ，B 两个横截面上的电磁场分布，计算出相应的功率，从而可推导得到传播常数。图 2 泄漏同轴电缆的侧面图注：计算 A 、 B 两个横截面处的电磁场，是为了避免在激励附近高次模的影响。3.2 FDTD 方法的 MATLAB 编程思路：1、 按所传输波频率大小以及所仿真物理模型大小合理划分网格，一般取所传输波最高频率所对应波长的 1/10-1/20 ，还要结合模型自身尺寸