一維有限時域差分方法（FDTD）計算中參數的選擇

摘要：在用一維時域有限差分方法計算波的傳播時，為了防止出現發散的結果，要求時間步長Δt<■，但具體取多大的值，沒有定論。通過編程試驗，時間步長取Δt<0.75■，而時間步數Nt大于波從波源傳到觀察點所需要的時間，又小于反射波到達觀察點的時間。而且計算的時間格點和位置格點全部都可以取整數。

關鍵詞：FDTD；有限時域差分方法；參數

中圖分類號：G642.4？搖 文獻標志碼：A？搖 文章編號：1674-9324（2013）03-0155-02

時域有限差分方法（FDTD）1966年由Yee提出后，經過許多年的發展，在電磁波散射方面得到了廣泛的應用。現在，在光子晶體能隙的計算方面也得到了很多應用。筆者在一維、二維的FDTD編程計算中，發現時間步長的選擇不同，會導致結果的不同，下面進行討論。

一、一維FDTD差分方程表達形式

一維波動方程■+■■=0，在一維無限長介質中，其解為u=f（x-vt），是以速度v向x軸正向傳播的行波。用差分形式，把波動方程進行改寫：

由u（x+Δx）=u（x）+■Δx+■■（Δx）2，u（x-Δx）=u（x）-■Δx+■■（Δx）2。兩式相減，u（x+Δx）-u（x-Δx）=2■Δx，于是，■=■。同理：■=■。于是，一維波動方程■+■■=0寫為：■+■■=0，整理得到：un+1（i）=un-1（i）-■[un（i+1）-un（i-1）] ①

二、初始、邊界條件的確定以及編程計算的結果

邊界條件：u（t，1）=0，u（t，Nx）=0，對t=1到Nt。

初始條件：由于在差分形式的遞推公式①中，要計算t=3的值，需要知道t=1，2時的值。本文規定初始條件為u（t，i）=0，對t=1，2.

下面是兩組不同的參數，用Matlab編程計算的結果：

第1組：Lx=1×10-4，Δx=1×10-6，Nx=100，Nt=100，v=c=3×108，Δt=■

信號源u（t，10）=5sin（■t），觀察點在i=80，也就是x=80Δx處。

第2組：Lx=1×10-4，Δx=1×10-6，Nx=100，Nt=1000，v=c=3×108，Δt=■

信號源u（t，10）=5sin（■t），觀察點在i=80，也就是x=80Δx處。參數中時間長度Nt與第1組不同，其他參數相同。

第3組：Lx=1×10-4，Δx=1×10-7，Nx=1000，Nt=3000，v=c=3×108，Δt=■■

信號源u（t，10）=5sin（■t），觀察點在i=500，也就是x=500Δx處。

從圖中可以看出：

t=1到600，在觀察點無信號，信號還沒有到達；

t=600到900，是過渡階段，與數值化的計算有關，在觀察點信號尚不穩定；

t=900到2400，觀察點信號與用連續方程得到的結果一致；

t=2400以后，觀察點的信號的幅度比源的信號還大，是因為文中取一維兩端邊界為零，反射波到達觀察點，與正向波形成干涉。

用快速傅里葉變換得出的頻譜透射率，由于已經包含了開始的過渡階段和以后的反射波的干涉的影響，在0.8左右，與波的無衰減傳播結果接近。

三、結論

在一維的FDTD計算中，直接使用整數時刻和整數空間坐標，而不需要使用半整數時間和半整數空間坐標，簡化了坐標的標記形式；時間步長取Δt=■■，時間總長Nt的取值，要求在這個時間信號源的信號可以到達觀察點，但又沒有界面反射波到達觀察點，這樣的參數比較理想，可以得到與連續方程相同的結果。

參考文獻：

[1]葛德彪，閻玉波.電磁波時域有限差分方法[M].第2版.西安：西安電子科技大學出版社，2005.

[2]曾輝，楊亞培.FDTD法與平面波展開法在光子禁帶計算中的差異分析[J].電子科技大學學報，2005，34（6）：901-904.

[3]張國華，袁乃昌，付云起.FDTD方法分析光子帶隙能帶結構[J].微波學報，2001，17（4）：14-17.

[4]張亮，寇曉艷.FDTD的Matlab語言的實現[J].延安大學學報：自然科學版，2009，28（2）：57-59.

[5]朱章虎，盧萬錚，馮奎勝.FDTD計算中PML的簡化應用及編程實現[J].空軍工程大學學報（自然科學版），2006，7（2）：55-57.

作者簡介：陳義萬，湖北工業大學理學院物理學副教授，研究方向：光子晶體及性質。