第二章 时域有限差分法_IV-三维FDTD

1、3. 三维三维FDTD重写Maxswell方程组00D1t H001t HE *rD=E4.1 自由空间自由空间FDTD公式公式（4-1-a）（4-1-b）（4-1-c）二维Maxswell方程组标量方程001yxzHHtxy D001xzHty E *zrzD=E001yzHtx E001yxzHDHtyz 001yxzDHHtzx 001yxzHHDtxy 001yxzEHEtzy 001yxzHEEtxz 001yxzEEHtyx 三维Maxswell方程组标量方程(4.2a)(4.2b)(4.2c)(4.2d)(4.2e)(4.2f)（4-1-a）（4-1-c）对时间和空间差分后，以对

2、时间和空间差分后，以(4.2c)和和(4.2f)为例，为例，迭代公式迭代公式1/21/2, ,1 2, ,1 2nnzzDi j kEi j k001 2, ,1 21 2, ,1 2nnyytHij kHij kx ,1 2,1 2,1 2,1 2nnxxHi jkHi jk11 2,1 2,1 2,1 2,nnzzHijkHijk1/21/2001 2,1,1 2, ,nnxxtEijkEij kx 1/21/21,1,1, ,nnyyEijkEij k dx(i,j,k) = dx(i,j,k) + 0.5 *( hz(i,j,k) - hz(i,j-1,k) . - hy(i,j,k)

3、 + hy(i,j,k-1);dy(i,j,k) = dy(i,j,k) + 0.5 *( hx(i,j,k) - hx(i,j,k-1) . - hz(i,j,k) + hz(i-1,j,k);dz(i,j,k) = dz(i,j,k) + 0.5 *( hy(i,j,k) - hy(i-1,j,k) . - hx(i,j,k) + hx(i,j-1,k);hx(i,j,k) = hx(i,j,k) + 0.5 *( ey(i,j,k+1) - ey(i,j,k) . - ez(i,j+1,k) + ez(i,j,k);hy(i,j,k) = hy(i,j,k) + 0.5 *( ez(i+

4、1,j,k) - ez(i,j,k) . - ex(i,j,k+1) + ex(i,j,k);hz(i,j,k) = hz(i,j,k) + 0.5 *( ex(i,j+1,k) - ex(i,j,k) . - ey(i+1,j,k) + ey(i,j,k);迭代公式 在三维基本FDTD程序和二维FDTD程序中，都是在问题空间的中心设置简单点源。但在三维仿真中，由点源产生的电场向外传播的过程中，与距离的平方成反比衰减。相反，我们可用电偶极子天线作为源。如下图所示。（1）在两臂空隙设置一电 场值。（2）依照安培环路定律设置围绕空隙的磁场值。CdlIH设置激励的方法1020304050601020

5、3040506000.20.40.60.8110203040506010203040506000.20.40.60.8110203040506010203040506000.20.40.60.8110203040506010203040506000.20.40.60.81T=50T=20T=40T=704.2 三维三维PML将二维问题的PML扩展至三维000( )( )11yxDDzHHxyjcjjxyD10000( )( )( )111yxzzzzHHxyzjcjjjxyD0000( )( )( )111yxzzzzHHxyzjcjjjxyD 000_zDzzccurlhcI1_DzIcur

6、lhj其中：由第三章知1 2, ,1 21 2, ,1 2_,1 2,1 2,1 2,1 2nnyynnxxHij kHij kcurlhHi jkHi jk 111, , ,_22nnDzDzi j ki j kcurlhII 1/21/2, ,1/ 233, ,1/ 21220.5_1, ,2nnzznDzi j kgii gjji j kgii gjjcurlhgkki j kDDI% Loop for D_z % for i=2:IE for j=2:JE for k=1:ka curl_h = hy(i,j,k) - hy(i-1,j,k). -hx(i,j,k) + hx(i,j-

7、1,k); idzl(i,j,k) = idzl(i,j,k) + curl_h; dz(i,j,k) = gi3(i)*gj3(j)*dz(i,j,k). +gi2(i)*gj2(j)* 0.5 *. (curl_h+ gk1(k)*idzl(i,j,k); end end end关于Dz的迭代程序，分三段计算以节约存储资源1 2, ,1 21 2, ,1 2_,1 2,1 2,1 2,1 2nnyynnxxHij kHij kcurlhHi jkHi jk 111, , ,_22nnDzDzi j ki j kcurlhII 1/21/2, ,1/ 233, ,1/ 21220.5_1,

8、,2nnzznDzi j kgii gjji j kgii gjjcurlhgkki j kDDIfor i=2:IE for j=2:JE for k=ka+1:kb+1 curl_h = hy(i,j,k) - hy(i-1,j,k). -hx(i,j,k) + hx(i,j-1,k); dz(i,j,k) = gi3(i)*gj3(j)*dz(i,j,k). +gi2(i)*gj2(j)* 0.5 *curl_h; end end endfor i=2:IE for j=2:JE for k=kb+2:KE kzh = k-kb-1; curl_h = hy(i,j,k) - hy(i-

9、1,j,k). -hx(i,j,k) + hx(i,j-1,k); idzh(i,j,kzh) = idzh(i,j,kzh) + curl_h; dz(i,j,k) = gi3(i)*gj3(j)*dz(i,j,k). +gi2(i)*gj2(j)* 0.5 *. (curl_h+ gk1(k)*idzh(i,j,kzh); end end end4.3 三维总三维总/散射场公式散射场公式在j=ja或 j=jb处，Dz _,1/ 2,1/ 2,0.5yayax inci jki jkjDDH _,1/ 2,1,1/ 2,10.5ybybx inci jki jkjDDH平面波图11-6-1平

10、面波的总场散射场表示方法散射体散射场区总场区散射体总场/散射场界面算例：平面波入射介质球% specify the dielectric sphere %radius = 10; epsilon = 30; sigma = 0.3;% Calculate gax, gbx %for i = ia+1:ib for j = ja+1:jb for k = ka+1:kb xdist = ic-i-0.5; ydist = jc-j; zdist = kc-k; dist= sqrt(xdist2 + ydist2 + zdist2); if ( dist = radius ) ga(i,j,k)

11、 = 1/(epsilon+(sigma*dt/epsz); gb(i,j,k) = sigma*dt/epsz; else ga(i,j,k) = 1; gb(i,j,k) = 0; end end end end% Calculate gay, gby %for i = ia+1:ib for j = ja+1:jb for k = ka+1:kb xdist = ic-i; ydist = jc-j-0.5; zdist = kc-k; dist= sqrt(xdist2 + ydist2 + zdist2); if ( dist = radius ) ga(i,j,k) = 1/(ep

12、silon+(sigma*dt/epsz); gb(i,j,k) = sigma*dt/epsz; else ga(i,j,k) = 1; gb(i,j,k) = 0; end end end end% Calculate gaz, gbz %for i = ia+1:ib for j = ja+1:jb for k = ka+1:kb xdist = ic-i; ydist = jc-j; zdist = kc-k-0.5; dist= sqrt(xdist2 + ydist2 + zdist2); if ( dist = radius ) ga(i,j,k) = 1/(epsilon+(s

13、igma*dt/epsz); gb(i,j,k) = sigma*dt/epsz; else ga(i,j,k) = 1; gb(i,j,k) = 0; end end end end4.4 FDTD仿真微带电路仿真微带电路问题描述：仿真一个尺寸为12.45mmX16mm的贴片天线的回波损耗S11。仿真的几个关键点：1、怎样划分网格？2、怎样取时间步长？3、源的处理4、导体的处理5、吸收边界6、S参数的计算1、网格划分、网格划分 感兴趣的频率范围是020GHz，不妨将25GHz取为频段的上限。则对应的波长最小值应该是 min =c / fmin=12mm 。考了到网格大小的取值应该小于等于 m

14、in /20, 所以仅仅从波长的角度考虑，应该有 min / 20 = 0.6mm(1) 划分方法一 z方向上的介质厚度为0.794mm，可划分3个网格， z近似z = 0.265mm。符合z max 的要求。 y方向上长度为16mm，取40个网格。则y = 0.24mm，也符合y max的要求。 x方向上有三个尺寸，12.45mm、2.09mm和2.46mm。天线的尺寸要尽准确，所以取x =0.389mm，32个网格，0.389x32=12.448mm（误差0.002mm）。微带线宽0.389x6=2.334mm（误差-0.126mm）。微带位置0.389x5=1.945mm（误差-0.11

15、5mm）。(1) 划分方法二 z方向上和 y方向上与第一种方法相同。 x方向上有三个尺寸：12.45mm、2.46mm和2.09mm。x =0.2075mm，60个网格，0.2075x60=12.45mm。微带线宽0.2075x12=2.49mm（误差0.03mm）。微带位置0.2075x10=2.075mm（误差-0.015mm）。采用第一种方法：贴片天线的尺寸为32x 32y ，模型总尺寸为60 x 100y 16z 。为了是编程简单，取用相对空间步长0.265_0.68120.389zraxx0.265_0.66250.4zrayy2、时间步长的取值、时间步长的取值 确定了空间步长x ，

16、 y ， z后，可以依据稳定性条件确定t 。0ztnc 222max1111tvxyz 简单点，可依据下面公式00.4412ztpsc 2max13vz3、源的处理、源的处理在导带口加强迫激励源，采用高斯脉冲： 220/t tTzEe 上式中：t0为延迟时间，T为高斯脉冲半宽度时间，馈源边界处理为磁壁。 max1/Tf fmax为高斯有效频谱的最高频率一般取t03T，可取 t0=3T。源加入方法源加入方法 在点V1，微带馈线和地板之间给定一个均匀电场。pulse = exp( -0.5 * ( t0-n ) / spread ) 2 );ez_inc(ja-1) =pulse;为消除j=ja点外地方产生的平面波，对输入要做修正。% 微带馈线下值为零，其余地方为0.for i = istart+1:iend+1 for k = 1:ktop+1 shape(i,k) = 1; endenddz(i,ja+1,k) = dz(i,ja+1,k) + . (0.5/eps_sub)*shape(i,k)*hx_inc(