均匀介质圆柱对TE和TM平面波的散射(附程序)

均匀介质圆柱对TE 和TM 平面波的散射
1、求解均匀介质圆柱对TM 波的雷达散射截面以及表面电、磁流
假设TM 极化均匀平面波垂直入射半径为a 的无限长均匀介质圆柱，其中介质圆柱沿z 轴放置，相对介电常数为εr ，相对磁导率为μr ，波的传播方向如图所示为+x 方向。

入射电场用柱面波展开，可表示为
00
cos 0000()i jk x jk n jn z z z n n E a E e a E e a E j J k e ρϕϕρ∞---=-∞
===∑
(1)
由Maxwell 方程E jw H μ∇⨯=-
,得到
1
i
i H E jw μ=-
∇⨯
1'0000000
1()e ()e n jn n jn n n n n E k E a nj J k a j J k jw jw ϕϕ
ρϕρρμρμ∞∞
-+-=-∞=-∞
=-+∑∑
(2)
其中，0μ为真空中的磁导率,0k 为真空中的波数。

当a ρ>时，介质外散射场朝外传播。

因此，散射电场用柱第二类Hankel 函数展开，表示如下
(2)
00()s n jn z n n n E a E j a H k e ϕ
ρ∞-=-∞
=∑ (3)
同理由Maxwell 方程E jw H μ∇⨯=-
,得到 ()()()
()s 220000000
1'n jn n jn n n n n n n E a k E H a j H k e a j a H k e j j ϕϕρϕρρωμρϕωμ∞∞
--=-∞=-∞
∂=-+∂∑∑
（4）
当a ρ<时，介质内为透射场，透射场则由柱面基本波函数的线性组合表示，由于透射场在介质内部均为有限大，则自变量为0时，纽曼函数和汉克尔函数趋于无穷，因此电场和磁场只选用贝塞尔函数表示如下
1()t n jn z n n n E a E j b J k e ϕρ∞
-=-∞
=∑
(5)
()()010
11
111't n jn n jn n n n n n n E b k E H a j J k e a j b J k e j j ϕϕρϕρρωμρϕωμ∞∞--=-∞=-∞∂=-+∂∑∑ (6)
其中，1μ为介质中的磁导率,1k 为介质中的波数。

当a ρ=时，根据介质表面的边界条件，切向电场和切向磁场连续，可以得到
电场边界条件
z z z i
s
t
a
a
a
E E E ρρρ===+=
则有
(2)
001()()()n n n n
n J k a a H k a b J k a +=
(7)
磁场边界条件
i
s t a a a
H H H ϕ
ρϕ
ρϕ
ρ===+=
则有 '
(2)''010010
1
[()()]()n n n n n k k J k a a H k a b J k a jw jw μμ+=
(8)
求解方程组，从而得到展开项的系数为
'
''
001110(2)'(2)
110010()()()()
()()()()
n n n n n n n n n J k a J k a J k a J k a a J k a H k a J k a H k a ηηηη-=- (9)
(2)
0011[()()]()
n n n n n b J k a a H k a J k a =
+ (10)
其中，因为εμω=k ，所以
η
μεωμ
1==
k
0001
k w ημ=
为真空中波阻抗的倒数,111
1
k w ημ=为介质中波阻抗的倒数。

如下求介质表面的电、磁流。

1）求表面电流
由边界条件=J n H ⨯
有
=()e ()i s
z J e e H e H H H ρρρϕϕϕϕ⨯+=⋅+ 所以'(2)'
0000
[J (()]n
jn z n n n k J j
k a a H k a e j ϕωμ∞
--∞
=
)+∑
'
(2)'0001[J (()]n jn n n n j k a a H k a e j ϕη∞--∞=)+∑ （13） 2）求表面磁流
由边界条件m =-J n E ⨯ 有
m =-()
i s z z z J e e E E ρ⨯+
=()
i
s z z e E E ϕ+
所以
(2)
00[J (()]m
n jn n n n J j k a a H k a e ϕϕ∞
--∞
=)+∑ （14）
另外对于远区散射场，kρ → ∞，()
()22n jk n j H k j e
k ρ
ρπρ
-≈
则散射电场为 00(2)000
00
02()2s jk n
jn n n jn z n n z n
n n jk jn z n
n j E a E j a H
k e
a E j a j e e k j a E a e e k ρϕ
ϕ
ρϕ
ρπρ
πρ
∞
∞
---=-∞
=-∞
∞
-=-∞
===∑
∑
∑
（11）
又
00cos 001
1i jk x jk z z E E a E e a E e ρϕ--⎧=⎪⎨===⎪⎩ 00
0022
s jk jn jn z n
n
n n j E a E a
e e a e
k k ρϕϕ
πρπρ
∞
∞
-=-∞
=-∞
==∑∑
将上式代入二维雷达散射截面的定义式
22()lim 2s
i
E E ρσϕπρ→∞
= 有 222
4
()lim 2s
jn n
n i
E a e
k E ϕ
ρσϕπρ∞
→∞=-∞
==∑ （12）
MATLAB 编程求解
%**********TM 波照射无限长介质柱*******
%***********初始化************ clear all;close all;clc;tic; wlen=1.0; k0=2.0*pi/wlen;
ur=1; %相对磁导率 er=4*10^12; %相对介电常数 eta0=120.0*pi; %自由空间波阻抗 eta1=eta0*sqrt(ur/er); %介质波阻抗 radius=10.0;
Npwave=10 ; NPL=2.0*pi*radius*Npwave palen=2.0*pi/NPL; ka=k0*radius;
ka1=ka*sqrt(ur*er); %介质波数
%!*********计算真空和介质中贝塞尔函数和汉克尔函数的值，存储到数组中************ jn0=besselj(0,ka); %真空中 h2n0=besselh(0,2,ka); jn(1)=besselj(1,ka); h2n(1)=besselh(1,2,ka); jn(2)=2.0*jn(1)/ka-jn0;
h2n(2)=2.0*h2n(1)/ka-h2n0;
Djn0=besselj(0,ka1); %介质中 Djn(1)=besselj(1,ka1); Djn(2)=2.0*Djn(1)/ka1-Djn0; for n=3:2000000 jn(n)=2.0*(n-1.0)*jn(n-1)/ka-jn(n-2); %真空中贝塞尔函数 h2n(n)=2.0*(n-1.0)*h2n(n-1)/ka-h2n(n-2); %真空中汉克尔函数 Djn(n)=2.0*(n-1.0)*Djn(n-1)/ka1-Djn(n-2); %介质中贝塞尔函数 if (abs(h2n(n))>1.0*10^10) break; end
end
ntotal=n
%!*****计算真空和介质中贝塞尔函数和汉克尔函数的一次导的值，存储到数组中***** Jn0=-jn(1); %真空中贝塞尔函数的一次导
H2n0=-h2n(1); %真空中汉克尔函数的一次导
DJn0=-Djn(1); %介质中贝塞尔函数的一次导
Jn(1)=(jn0-jn(2))/2;
H2n(1)=(h2n0-h2n(2))/2;
DJn(1)=(Djn0-Djn(2))/2;
jn(ntotal+1)=2.0*ntotal*jn(ntotal)/ka-jn(ntotal-1);
h2n(ntotal+1)=2.0*ntotal*h2n(ntotal)/ka-h2n(ntotal-1);
Djn(ntotal+1)=2.0*ntotal*Djn(ntotal)/ka1-Djn(ntotal-1);
for n=2:ntotal
Jn(n)=(jn(n-1)-jn(n+1))/2;
H2n(n)=(h2n(n-1)-h2n(n+1))/2;
DJn(n)=(Djn(n-1)-Djn(n+1))/2;
end
%!*********计算an**********
a0=(eta0*DJn0*jn0-eta1*Jn0*Djn0)/(eta1*H2n0*Djn0-eta0*h2n0*DJn0);
for n=1:ntotal
a(n)=(eta0*DJn(n)*jn(n)-eta1*Jn(n)*Djn(n))/...
(eta1*H2n(n)*Djn(n)-eta0*h2n(n)*DJn(n));
end
%!*******************计算每个节点处的解析电流**********
for m=1:NPL
pp=(m-0.5)*palen; %散射角Φ
cp1=sin(pp)-i*cos(pp);
cp2=-sin(pp)-i*cos(pp);
ctemp1=1.0;
ctemp2=1.0;
sumc0=Jn0+a0*H2n0;
sumf0=jn0+a0*h2n0;
for n=1:ntotal
ctemp1=ctemp1*cp1;
ctemp2=ctemp2*cp2;
sumc0=sumc0+(Jn(n)+a(n)*H2n(n))*(ctemp1+ctemp2);
sumf0=sumf0+(jn(n)+a(n)*h2n(n))*(ctemp1+ctemp2);
end
Jc(m)=sumc0/eta0*(-j); %电流密度
Jf(m)=sumf0; %磁流密度
end
figure,plot(abs(Jc));grid;title('电流分布');
xlabel('散射角Φ(弧度)');ylabel('电流密度(A/m)');
figure,plot(abs(Jf));grid;title('磁流分布');
xlabel('散射角Φ(弧度)');ylabel('磁流密度(Wb/m)');
%*********************计算每个散射方向的RCS********************************* for m=1:360 pp=m*2.0*pi/360.0; cp1=cos(pp)+i*sin(pp); cp2=cos(pp)-i*sin(pp); ctemp1=1.0; ctemp2=1.0; sum=a0; for n=1:ntotal ctemp1=ctemp1*cp1; ctemp2=ctemp2*cp2;
sum=sum+(ctemp1+ctemp2)*a(n); end
RCS_exact(m)=10.0*log10(4.0/k0*abs(sum)^2); end
figure,plot(abs(RCS_exact));grid;title('RCS 分布'); xlabel('散射角Φ(度)');ylabel('RCS(dbm)'); toc;
运行结果： NPL = 628.3185 ntotal = 100
Elapsed time is 1.295000 seconds.
0100200
300400500600700
1
2
3
4
5
6x 10
-3
电流分布
散射角Φ(弧度)
电流密度(A /m )
0100200
300400500600700
0.2
0.4
0.6
0.81
1.2
x 10
-6
磁流分布
散射角Φ(弧度)
磁流密度(W b /m )
050100150
200250300350400
5
10
15
20
25
30
35
RCS 分布
散射角Φ(度)
R C S (d b m )
2、求解均匀介质圆柱对TE 波的雷达散射截面以及表面电、磁流
假设TE 极化均匀平面波垂直入射半径为a 的无限长均匀介质圆柱，其中介质圆柱沿z 轴放置，相对介电常数为εr ，相对磁导率为μr ，波的传播方向如图所示为+x 方向。

入射磁场用柱面波展开，可表示为
00cos 0000()i jk x jk n jn z z z n n H a H e a H e a H j J k e ρϕ
ϕρ∞---=-∞
===∑
(1)
由Maxwell 方程H jw E ε∇⨯=
,得到
1
i i E H jw ε=
∇⨯
1'0000000
1()e ()e n jn n jn n n n n H k H a nj J k a j J k jw jw ϕϕ
ρϕρρερε∞∞
-+-=-∞=-∞
=-∑∑
(2)
其中，0ε为真空中的介电常数,0k 为真空中的波数。

当a ρ>时，介质外散射场朝外传播。

因此，散射磁场用柱第二类Hankel 函数展开，表示如下
(2)
0()s n jn z n n n H a H j a H k e ϕ
ρ∞
-=-∞
=∑
(3)
同理由Maxwell 方程 H jw E ε∇⨯=
,得到 ()()()()s 220000000
1'n jn n
jn n n n n n n H a k H E a j H k e a j
a H k e j j ϕϕρϕρρωερϕωε∞∞
--=-∞=-∞
∂=-∂∑∑ （4）
当a ρ<时，介质内为透射场，透射场则由柱面基本波函数的线性组合表示，由于透射场在介质内部均为有限大，则自变量为0时，纽曼函数和汉克尔函数趋于无穷，因此电场和磁场只选用贝塞尔函数表示如下
01()t n jn z n n n H a H j b J k e ϕρ∞-=-∞
=∑
(5)
()()01011111't n jn n jn n n n n n n H b k H E a j J k e a j b J k e j j ϕϕρϕρρωερϕωε∞∞--=-∞=-∞∂=-∂∑∑ (6)
其中，1ε为介质中的介电常数,1k 为介质中的波数。

当a ρ=时，根据介质表面的边界条件，切向电场和切向磁场连续，可以得到
磁场边界条件
z z
z
i
s
t
a a a
H H H ρρρ===+=
则有
(2)
001()()()n n n n
n J k a a H k a b J k a +=
(7)
电场边界条件
i
s t
a a a
E E E ϕ
ρϕ
ρϕ
ρ===+=
则有 '
(2)''010010
1
[()()]()n n n n n k k J k a a H k a b J k a jw jw εε+=
(8)
求解方程组，从而得到展开项的系数为
''
010101'(2)(2)'
110010()()()()
()()()()
n n n n n n n n n J k a J k a J k a J k a a J k a H k a J k a H k a ηηηη-=- (9)
(2)''(2)
0000(2)''(2)111110
()()()()
()()()()
n n n n n n n n n J k a H k a J k a H k a b J k a H k a J k a H k a ηη-=
- (10)
其中，因为εμω=k ，所以
k μ
ηωε
ε
=
= 000k w ηε=为真空中波阻抗,111
k
w ηε=为介质中波阻抗。

如下求介质表面的电、磁流。

1）求表面电流
由边界条件=J n H ⨯ 有
=e ()i s z z p z z z J e e H e H H ρ⨯=⨯+
e ()i
s z z H H ϕ=-+
所以(2)000[J (()]n jn n n n J H j k a a H k a e ϕϕ∞
--∞=-)+∑ （13）
2)求表面磁流
由边界条件m =-J n E ⨯ 有
m =-()J e e E e E ρρρϕϕ⨯+
=-Z e E ϕ
a a =-()i
s
Z e E E ϕρϕρ==+
'(2)'00000(1)[J (()]n jn z n n n k H e j k a a H k a e j ϕωε∞--∞
=--)+∑ 000
'(2)'000
[J (()]n
jn z
n n n w H e j
k a a H k a e j ϕμεωε∞
--∞
=)+∑
'
(2)'0000[J (()]n jn z n n n H e j k a a H k a e j ϕη∞--∞
=)+∑
所以
m '(2)'00
00[J (()]n
jn z
n n n H J j
k a a H k a e
j
ϕη∞
--∞
=
)+∑ （14）
另外对于远区散射场，kρ → ∞，()
()22n jk n j H k j e
k ρ
ρπρ
-≈
则散射电场为 00(2)000
00
02()2s jk n
jn n n jn z n n z n
n n jk jn z n
n j H a H j a H
k e
a H j a j e e k j a H a e e k ρϕ
ϕ
ρϕ
ρπρ
πρ
∞
∞
---=-∞
=-∞
∞
-=-∞
===∑
∑
∑
（11）
又
000cos 001
1i jk x jk z z H H a H e a H e ρϕ--⎧=⎪⎨===⎪⎩
00
0022
s jk jn jn z n
n
n n j H a H a
e e a e
k k ρϕϕ
πρπρ
∞
∞
-=-∞
=-∞
==∑∑
将上式代入二维雷达散射截面的定义式
2
2()lim 2s
i
H H ρσϕπρ→∞
=
有 222
4
()lim 2s
jn n
n i
H a e
k H ϕ
ρσϕπρ∞
→∞=-∞
==∑ （12）
MATLAB 编程求解
%**********TE 波照射无限长介质柱*******
%***********初始化************ clear all;close all;clc;tic; wlen=1.0; k0=2.0*pi/wlen;
ur=1; %相对磁导率 er=4*10^12; %相对介电常数 eta0=120.0*pi; %自由空间波阻抗 eta1=eta0*sqrt(ur/er); %介质波阻抗 radius=10.0;
Npwave=10 ; NPL=2.0*pi*radius*Npwave palen=2.0*pi/NPL; ka=k0*radius;
ka1=ka*sqrt(ur*er); %介质波数
%!*********计算真空和介质中贝塞尔函数和汉克尔函数的值，存储到数组中************ jn0=besselj(0,ka); %真空中 h2n0=besselh(0,2,ka); jn(1)=besselj(1,ka); h2n(1)=besselh(1,2,ka); jn(2)=2.0*jn(1)/ka-jn0;
h2n(2)=2.0*h2n(1)/ka-h2n0;
Djn0=besselj(0,ka1); %介质中 Djn(1)=besselj(1,ka1); Djn(2)=2.0*Djn(1)/ka1-Djn0; for n=3:2000000 jn(n)=2.0*(n-1.0)*jn(n-1)/ka-jn(n-2); %真空中贝塞尔函数 h2n(n)=2.0*(n-1.0)*h2n(n-1)/ka-h2n(n-2); %真空中汉克尔函数
Djn(n)=2.0*(n-1.0)*Djn(n-1)/ka1-Djn(n-2); %介质中贝塞尔函数
if (abs(h2n(n))>1.0*10^10) break;
end
end
ntotal=n
%!*****计算真空和介质中贝塞尔函数和汉克尔函数的一次导的值，存储到数组中***** Jn0=-jn(1); %真空中贝塞尔函数的一次导
H2n0=-h2n(1); %真空中汉克尔函数的一次导
DJn0=-Djn(1); %介质中贝塞尔函数的一次导
Jn(1)=(jn0-jn(2))/2;
H2n(1)=(h2n0-h2n(2))/2;
DJn(1)=(Djn0-Djn(2))/2;
jn(ntotal+1)=2.0*ntotal*jn(ntotal)/ka-jn(ntotal-1);
h2n(ntotal+1)=2.0*ntotal*h2n(ntotal)/ka-h2n(ntotal-1);
Djn(ntotal+1)=2.0*ntotal*Djn(ntotal)/ka1-Djn(ntotal-1);
for n=2:ntotal
Jn(n)=(jn(n-1)-jn(n+1))/2;
H2n(n)=(h2n(n-1)-h2n(n+1))/2;
DJn(n)=(Djn(n-1)-Djn(n+1))/2;
end
%!*********计算an**********
a0=(eta0*Djn0*Jn0-eta1*jn0*DJn0)/(eta1*h2n0*DJn0-eta0*H2n0*Djn0);
for n=1:ntotal
a(n)=(eta0*Djn(n)*Jn(n)-eta1*jn(n)*DJn(n))/...
(eta1*h2n(n)*DJn(n)-eta0*H2n(n)*Djn(n));
end
%!*******************计算每个节点处的解析电流**********
for m=1:NPL
pp=(m-0.5)*palen; %散射角Φ
cp1=sin(pp)-i*cos(pp);
cp2=-sin(pp)-i*cos(pp);
ctemp1=1.0;
ctemp2=1.0;
sumc0=jn0+a0*h2n0;
sumf0=Jn0+a0*H2n0;
for n=1:ntotal
ctemp1=ctemp1*cp1;
ctemp2=ctemp2*cp2;
sumc0=sumc0+(jn(n)+a(n)*h2n(n))*(ctemp1+ctemp2);
sumf0=sumf0+(Jn(n)+a(n)*H2n(n))*(ctemp1+ctemp2);
end
Jc(m)=sumc0; %电流密度
Jf(m)=sumf0*eta0*(-j); %磁流密度
end
figure,plot(abs(Jc));grid;title('电流分布');
xlabel('散射角Φ(弧度)');ylabel('电流密度(A/m)'); figure,plot(abs(Jf));grid;title('磁流分布');
xlabel('散射角Φ(弧度)');ylabel('磁流密度(Wb/m)');
%*********************计算每个散射方向的RCS********************************* for m=1:360 pp=m*2.0*pi/360.0; cp1=cos(pp)+i*sin(pp); cp2=cos(pp)-i*sin(pp); ctemp1=1.0; ctemp2=1.0; sum=a0; for n=1:ntotal ctemp1=ctemp1*cp1; ctemp2=ctemp2*cp2;
sum=sum+(ctemp1+ctemp2)*a(n); end
RCS_exact(m)=10.0*log10(4.0/k0*abs(sum)^2); end
figure,plot(abs(RCS_exact));grid;title('RCS 分布'); xlabel('散射角Φ(度)');ylabel('RCS(dbm)'); toc;
运行结果
NPL = 628.3185
ntotal = 100
Elapsed time is 1.263000 seconds.
100
200
300400500
600
700
00.20.40.60.811.21.4
1.61.8
2电流分布
散射角Φ(弧度)
电流密度(A /m )
0100200
300400500600700
1
2
3
x 10
-4
磁流分布
散射角Φ(弧度)
磁流密度(W b /m )
050100150
200250300350400
5
10
15
20
25
30
35
RCS 分布
散射角Φ(度)
R C S (d b m )。