時域有限差分法課件

時域有限差分法第1講

一維標量波動方程時域有限差分法第1講一維標量波動方程1引言（1）1966年,K.S.Yee（美籍香港人）首先提出了Finite-DifferenceTime-DomainMethod,并用于柱形金屬柱電磁散射分析。由于當時計算機技術(shù)還比較落后，這一方法并未引起重視。1972年,A.Taflovey應用FDTD研究了UHF和微波對人類眼睛的穿透，以了解“微波白內(nèi)障”的成因。Taflove成功地應用和發(fā)展了Yee的FDTD算法。80年代后期，隨著高速大容量計算機的普及，F(xiàn)DTD法得到了迅速發(fā)展。如今已應用于涉及波動現(xiàn)象的任何領域。至今，F(xiàn)DTD法的研究與應用仍方興未艾。引言（1）1966年,K.S.Yee（美籍香港人）首先提出2引言（2）

本課程采用研討班形式。教師講授FDTD的基本知識，學生針對某一方向進行較深入的研究。本講我們考慮描述波動現(xiàn)象的最基本偏微分方程：一維標量波動方程的數(shù)值FDTD解，為以后二維、三維Maxwell方程的FDTD分析奠定基礎課程內(nèi)容取自下列的參考書和近年來相關(guān)的一些文獻[1]A.Taflove,ComputationalElectrodynamicsTheFinite-DifferenceTime-DomainMethod,ArtechHourse，1995.[2]高本慶，時域有限差分法，國防工業(yè)出版社，1995.[3]葛德彪，閆玉波，電磁場時域有限差分法，西電出版社，2002引言（2）本課程采用研討班形式。教師講授FDTD的基本知31.1差分近似（1）一維標量波動方程（1-1）上式的解為（1-2）

采用Taylor展開（1-3）1.1差分近似（1）一維標量波動方程41.1差分近似（2）于是，有(1-4）同理，有(1-5）上式稱為二階偏導數(shù)的二階中心差分格式。將它們代入(1-1)，得（1-6）忽略高次項，便可得到求解的差分迭代公式。1.1差分近似（2）于是，有51.1差分近似（3）NoYesn=0在所有空間點給uin,uin-1(i=1:imax)賦初值n=n+1由(1-6)在所有空間點求uin+1(i=1:imax)結(jié)束n>nmax?圖1.1一維波動方程FDTD流程圖1.1差分近似（3）NoYesn=0n=n+1由(1-6)6

1.1差分近似（4）

應當注意，在一般情況下(1-6)對時間或空間具有二階精度。但對于的特殊情況，根據(jù)解(1-2)，可以證明

于是

所以，（1-6）中的兩個剩余項抵消，得到了精確的數(shù)值差分公式

（1-7）正因為有這樣的奇妙特性，為“魔時間步”(Magictimestep).1.1差分近似（4）應當注意，在一般情況下71.2數(shù)值色散關(guān)系(1)色散關(guān)系定義為行波的波長隨頻率的變化關(guān)系。為方便起見，色散關(guān)系也常表示為行波的波數(shù)關(guān)于角頻率的變化關(guān)系?？紤](1.1)的正弦行波解代入(1-1)得即

(1-8)上式便是一維標量波動方程的色散關(guān)系。

由上式得相速度（1-9）可見，相速與頻率無關(guān)，稱為非色散。非色散意味著對于具有任意調(diào)制的包絡或脈沖形狀的波傳播任意距離后波形保持不變。進一步由(1-8)可以得到群速關(guān)系(1-10)這種情況下，群速也是與頻率無關(guān)。1.2數(shù)值色散關(guān)系(1)色散關(guān)系定義為行波81.2數(shù)值色散關(guān)系(2)

上述過程也可用于一維標量波動方程差分近似的數(shù)值色散分析。設在離散空間點,離散行波解為，式中，為存在于有限差分網(wǎng)格中的數(shù)值正弦波的波數(shù)。一般情況下，不同于連續(xù)物理波的波數(shù)。正是這種不同導致了數(shù)值相速和群速偏離了精確解。進而導致了數(shù)值色散誤差。將上式代入差分方程(1-6),得(1-11)重新組合并應用Euler恒等式，最后得到數(shù)值色散關(guān)系為(1-12)1.2數(shù)值色散關(guān)系(2)上述過程也可用于一維標量91.3數(shù)值相速（1）類似于(1-9)，定義數(shù)值相速為由（1-12）可得（1-13）可見數(shù)值相速與頻率有關(guān)。因此，由FDTD得到的數(shù)值波是色散的。取則數(shù)值相速為。相對誤差為-1.27%。如果物理波傳播了距離(100空間格)時，數(shù)值模擬波只傳播了98.73空間格，相位誤差為45.720。取則。這時數(shù)值相速的相對誤差為0.31，減少了4倍。同樣，當物理波傳播了同樣的時(200空間格)，數(shù)值模擬傳播了199.378格，相位誤差為11.1960，也減少了4倍。誤差減少了4倍反映了差分算法是二階精度的。1.3數(shù)值相速（1）類似于(1-9)，定義數(shù)值相速101.3數(shù)值相速（2）情況1：非常細網(wǎng)格根據(jù)，數(shù)值色散關(guān)系(1-12)變?yōu)榧?，，最后得，于是有。所以，在非常細的網(wǎng)格條件下，差分解逼近精確解。情況2：魔時間步

(1-12)變?yōu)?即。所以，。可見，魔時間步下差分解與精確解相同。1.3數(shù)值相速（2）情況1：非常細網(wǎng)格111.4數(shù)值群速定義數(shù)值相速為（1-14）情況1非常細網(wǎng)格利用正弦函數(shù)的一階Taylor展開，可得（1-15）所以，群速與相速一樣，在細網(wǎng)格條件下趨近精確解。這證明了當空間步長和時間步長趨于零時，數(shù)值解變得精確。

情況2魔時間步將魔時間步條件和波數(shù)代入(1-14)，得（1-16）再次驗證了魔時間步下數(shù)值解等于精確解。1.4數(shù)值群速定義數(shù)值相速為121.5數(shù)值穩(wěn)定性（1）FDTD計算中每一步都是有誤差的，隨著時間步進，誤差會不斷積累。如果誤差的積累不會造成總誤差的增加，就成FDTD法是穩(wěn)定的，否則成為不穩(wěn)定的。數(shù)值不穩(wěn)定性會造成計算結(jié)果隨時間步進無限增加。FDTD法是有條件穩(wěn)定的，即：時間步必須必須小于一定值以避免數(shù)值不穩(wěn)定性。本節(jié)的數(shù)值穩(wěn)定性分析方法是建立在Courant等人幾十年前提出的經(jīng)典方法基礎上。這種方法首先把有限差分算法分解為相互分離的時間和空間本征值問題。1.5數(shù)值穩(wěn)定性（1）FDTD計算中每一步都是有誤差的，隨131.5數(shù)值穩(wěn)定性（1）時間本征值問題(1-17)差分近似，得(1-18)定義不變增長因子(1-19)1.5數(shù)值穩(wěn)定性（1）時間本征值問題141.5數(shù)值穩(wěn)定性（2）將(1-19)代入(1-18)，有，于是

算法穩(wěn)定性要求。如果，則總有，于是，滿足穩(wěn)定性要求。這樣可得（1-20）這就是穩(wěn)定的數(shù)值差分解所要求的時間本征值譜。1.5數(shù)值穩(wěn)定性（2）將(1-19)代入(1-18)，有151.5數(shù)值穩(wěn)定性（3）空間本征值問題(1-21)代入中心差分公式，得(1-22)令，Eular公式可得因為，所以(1-23)上式給出了差分網(wǎng)格中任意空間Fourier模的本征值譜。1.5數(shù)值穩(wěn)定性（3）空間本征值問題161.5數(shù)值穩(wěn)定性（4）穩(wěn)定性為了保證任何空間模式的數(shù)值穩(wěn)定性，(1-23)給出的空間模式的本征值范圍必須完全落在(1-20)所給出的時間本征值的穩(wěn)定范圍內(nèi)，于是即(1-24)可見，時間步長必須是有界的。上式稱為Courant穩(wěn)定性條件。有趣的是其上界恰好是魔時間步。1.5數(shù)值穩(wěn)定性（4）穩(wěn)定性171.6激勵源的設置在FDTD模擬電磁波傳播時需要設置初始條件和激勵源。最簡單的源設置方法是“硬源”,即在激勵源的位置令u滿足ui=f(n),常用的有正弦函數(shù)ui=sin(nt+)高斯函數(shù)ui=exp[-(n-n0)2/T2]階躍函數(shù)ui=0n<n1=(n-n1)/(n2-n1)n1<n<n2=1n>n2“硬源”設置簡單，但當反射波回到“硬源”位置時，會引起寄生反射，所以，要在這之前“關(guān)”掉源。以后會有有關(guān)源設置的更詳細討論。1.6激勵源的設置在FDTD模擬電磁波傳播時181.7吸收邊界條件由于計算機容量所限，計算域必須是有限的。對于理想電壁或磁壁的邊界條件的設置是直接的。但如果模擬的是“開”問題，就要設置截斷邊界。在截斷邊界上要設置吸收邊界條件，使得電磁波可以被完全吸收，模擬波無反射的通過吸收邊界。對于一維問題，采用單向波方程于是利用單向差分近似得到吸收邊界條件，詳細討論見后面章節(jié)。1.7吸收邊界條件由于計算機容量所限，計算19結(jié)論1

本講介紹了一維標量波動方程的FDTD求解過程：利用Taylor級數(shù)展開方法獲取空間/時間導數(shù)的二階中心差分近似，從而得到具有二階精度的方程數(shù)值解的時間步進迭代公式。一般情況下，數(shù)值解引入了寄生的數(shù)值色散。當空間步長和時間步長非常小時，數(shù)值解逼近精確解。當時間步長滿足魔時間步條件時，數(shù)值解等于精確解?？臻g步長和時間步長必須滿足Courant穩(wěn)定性條件才能保證數(shù)值解的穩(wěn)定性。結(jié)論1本講介紹了一維標量波動方程的FDTD求解過20習題11.1利用Taylor級數(shù)展開方法分別推導一階導數(shù)的二階和四階精度中心差分近似。1.2利用數(shù)值相速和群速公式分別畫出數(shù)值相速和群速在，，和條件下關(guān)于網(wǎng)格空間分辨率的曲線，并進行相應的討論。1.3編寫本講介紹的一維標量波動方程FDTD求解程序。設在網(wǎng)格左邊界波源為下列脈沖函數(shù)：a.高斯函數(shù)b.矩形函數(shù)在網(wǎng)格右邊界。時間總步數(shù)為5000,空間總步為200。在，和條件下每隔1000個時間步畫出每個傳播脈沖關(guān)于位置的分布曲線。習題11.1利用Taylor級數(shù)展開方法分別推導一階21時域有限差分法第1講

一維標量波動方程時域有限差分法第1講一維標量波動方程22引言（1）1966年,K.S.Yee（美籍香港人）首先提出了Finite-DifferenceTime-DomainMethod,并用于柱形金屬柱電磁散射分析。由于當時計算機技術(shù)還比較落后，這一方法并未引起重視。1972年,A.Taflovey應用FDTD研究了UHF和微波對人類眼睛的穿透，以了解“微波白內(nèi)障”的成因。Taflove成功地應用和發(fā)展了Yee的FDTD算法。80年代后期，隨著高速大容量計算機的普及，F(xiàn)DTD法得到了迅速發(fā)展。如今已應用于涉及波動現(xiàn)象的任何領域。至今，F(xiàn)DTD法的研究與應用仍方興未艾。引言（1）1966年,K.S.Yee（美籍香港人）首先提出23引言（2）

本課程采用研討班形式。教師講授FDTD的基本知識，學生針對某一方向進行較深入的研究。本講我們考慮描述波動現(xiàn)象的最基本偏微分方程：一維標量波動方程的數(shù)值FDTD解，為以后二維、三維Maxwell方程的FDTD分析奠定基礎課程內(nèi)容取自下列的參考書和近年來相關(guān)的一些文獻[1]A.Taflove,ComputationalElectrodynamicsTheFinite-DifferenceTime-DomainMethod,ArtechHourse，1995.[2]高本慶，時域有限差分法，國防工業(yè)出版社，1995.[3]葛德彪，閆玉波，電磁場時域有限差分法，西電出版社，2002引言（2）本課程采用研討班形式。教師講授FDTD的基本知241.1差分近似（1）一維標量波動方程（1-1）上式的解為（1-2）

采用Taylor展開（1-3）1.1差分近似（1）一維標量波動方程251.1差分近似（2）于是，有(1-4）同理，有(1-5）上式稱為二階偏導數(shù)的二階中心差分格式。將它們代入(1-1)，得（1-6）忽略高次項，便可得到求解的差分迭代公式。1.1差分近似（2）于是，有261.1差分近似（3）NoYesn=0在所有空間點給uin,uin-1(i=1:imax)賦初值n=n+1由(1-6)在所有空間點求uin+1(i=1:imax)結(jié)束n>nmax?圖1.1一維波動方程FDTD流程圖1.1差分近似（3）NoYesn=0n=n+1由(1-6)27

1.1差分近似（4）

應當注意，在一般情況下(1-6)對時間或空間具有二階精度。但對于的特殊情況，根據(jù)解(1-2)，可以證明

于是

所以，（1-6）中的兩個剩余項抵消，得到了精確的數(shù)值差分公式

（1-7）正因為有這樣的奇妙特性，為“魔時間步”(Magictimestep).1.1差分近似（4）應當注意，在一般情況下281.2數(shù)值色散關(guān)系(1)色散關(guān)系定義為行波的波長隨頻率的變化關(guān)系。為方便起見，色散關(guān)系也常表示為行波的波數(shù)關(guān)于角頻率的變化關(guān)系。考慮(1.1)的正弦行波解代入(1-1)得即

(1-8)上式便是一維標量波動方程的色散關(guān)系。

由上式得相速度（1-9）可見，相速與頻率無關(guān)，稱為非色散。非色散意味著對于具有任意調(diào)制的包絡或脈沖形狀的波傳播任意距離后波形保持不變。進一步由(1-8)可以得到群速關(guān)系(1-10)這種情況下，群速也是與頻率無關(guān)。1.2數(shù)值色散關(guān)系(1)色散關(guān)系定義為行波291.2數(shù)值色散關(guān)系(2)

上述過程也可用于一維標量波動方程差分近似的數(shù)值色散分析。設在離散空間點,離散行波解為，式中，為存在于有限差分網(wǎng)格中的數(shù)值正弦波的波數(shù)。一般情況下，不同于連續(xù)物理波的波數(shù)。正是這種不同導致了數(shù)值相速和群速偏離了精確解。進而導致了數(shù)值色散誤差。將上式代入差分方程(1-6),得(1-11)重新組合并應用Euler恒等式，最后得到數(shù)值色散關(guān)系為(1-12)1.2數(shù)值色散關(guān)系(2)上述過程也可用于一維標量301.3數(shù)值相速（1）類似于(1-9)，定義數(shù)值相速為由（1-12）可得（1-13）可見數(shù)值相速與頻率有關(guān)。因此，由FDTD得到的數(shù)值波是色散的。取則數(shù)值相速為。相對誤差為-1.27%。如果物理波傳播了距離(100空間格)時，數(shù)值模擬波只傳播了98.73空間格，相位誤差為45.720。取則。這時數(shù)值相速的相對誤差為0.31，減少了4倍。同樣，當物理波傳播了同樣的時(200空間格)，數(shù)值模擬傳播了199.378格，相位誤差為11.1960，也減少了4倍。誤差減少了4倍反映了差分算法是二階精度的。1.3數(shù)值相速（1）類似于(1-9)，定義數(shù)值相速311.3數(shù)值相速（2）情況1：非常細網(wǎng)格根據(jù)，數(shù)值色散關(guān)系(1-12)變?yōu)榧?，，最后得，于是有。所以，在非常細的網(wǎng)格條件下，差分解逼近精確解。情況2：魔時間步

(1-12)變?yōu)?即。所以，?？梢姡r間步下差分解與精確解相同。1.3數(shù)值相速（2）情況1：非常細網(wǎng)格321.4數(shù)值群速定義數(shù)值相速為（1-14）情況1非常細網(wǎng)格利用正弦函數(shù)的一階Taylor展開，可得（1-15）所以，群速與相速一樣，在細網(wǎng)格條件下趨近精確解。這證明了當空間步長和時間步長趨于零時，數(shù)值解變得精確。

情況2魔時間步將魔時間步條件和波數(shù)代入(1-14)，得（1-16）再次驗證了魔時間步下數(shù)值解等于精確解。1.4數(shù)值群速定義數(shù)值相速為331.5數(shù)值穩(wěn)定性（1）FDTD計算中每一步都是有誤差的，隨著時間步進，誤差會不斷積累。如果誤差的積累不會造成總誤差的增加，就成FDTD法是穩(wěn)定的，否則成為不穩(wěn)定的。數(shù)值不穩(wěn)定性會造成計算結(jié)果隨時間步進無限增加。FDTD法是有條件穩(wěn)定的，即：時間步必須必須小于一定值以避免數(shù)值不穩(wěn)定性。本節(jié)的數(shù)值穩(wěn)定性分析方法是建立在Courant等人幾十年前提出的經(jīng)典方法基礎上。這種方法首先把有限差分算法分解為相互分離的時間和空間本征值問題。1.5數(shù)值穩(wěn)定性（1）FDTD計算中每一步都是有誤差的，隨341.5數(shù)值穩(wěn)定性（1）時間本征值問題(1-17)差分近似，得(1-18)定義不變增長因子(1-19)1.5數(shù)值穩(wěn)定性（1）時間本征值問題351.5數(shù)值穩(wěn)定性（2）將(1-19)代入(1-18)，有，于是

算法穩(wěn)定性要求。如果，則總有，于是，滿足穩(wěn)定性要求。這樣可得（1-20）這就是穩(wěn)定的數(shù)值差分解所要求的時間本征值譜。1.5數(shù)值穩(wěn)定性（2）將(1-19)代入(1-18)，有361.5數(shù)值穩(wěn)定性（3）空間本征值問題(1-21)代入中心差分公式，得(1-22)令，Eular公式可得因為，所以(1-23)上式給出了差分網(wǎng)格中任意空間Fourier模的本征值譜。1.5數(shù)值穩(wěn)定性（3）空間本征值問題371.5數(shù)值穩(wěn)定性（4）穩(wěn)定性為了保證任何空間模式的數(shù)值穩(wěn)定性，(1-23)給出的空間模式的本征值范圍必須完全落在(1-20)所給出的時間本征值的穩(wěn)定范圍內(nèi)，于是即(1-24)可見，時間步長必須是有界的。上式稱為Courant穩(wěn)定性條件。有趣的是其上界恰好是魔時間步。1.5數(shù)值穩(wěn)定性（4）穩(wěn)定性381.6激勵源的設置在FDTD模擬電磁波傳播時需要設置初始條件和激勵源。最簡單的源設置方法是“硬源”,即在激勵源的位置令u滿足ui=f(n),常用的有正弦函數(shù)ui=sin(nt+)高斯函數(shù)ui=exp[-(n-n0)2/T2]階躍函數(shù)ui=0