數(shù)學(xué)物理方程三維可視化仿真

摘要：數(shù)學(xué)物理方程三維可視化仿真及創(chuàng)新實(shí)踐訓(xùn)練是《數(shù)學(xué)物理方法》教學(xué)模式改革中的重要內(nèi)容。本文通過MATLAB程序求解二維菱形晶格光子晶體的電磁場本征值方程，繪制出二維能帶曲線，并將結(jié)果三維可視化，體現(xiàn)出復(fù)雜數(shù)學(xué)物理問題的物理圖像，解決大學(xué)生在課程學(xué)習(xí)過程中理解困難的教學(xué)問題，加強(qiáng)大學(xué)生編程實(shí)踐能力和創(chuàng)新能力的培養(yǎng)。

關(guān)鍵詞：本征值問題；三維可視化仿真；光子晶體；平面波展開法；能帶結(jié)構(gòu)

中圖分類號(hào)：G642.0 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：1674-9324（2013）03-0247-03

一、課程背景

《數(shù)學(xué)物理方法》是理工科學(xué)生的基礎(chǔ)課程之一，也是科研中常用的基本方法。數(shù)學(xué)物理方法課程的內(nèi)容繁多，公式推導(dǎo)繁雜，盡管教材中的例題通常具有明確的物理意義，但是從眼花繚亂的數(shù)學(xué)表達(dá)式中看出其中所表達(dá)的物理圖像，不僅學(xué)生會(huì)覺得困惑、枯燥，教師也難免覺得棘手。探索數(shù)學(xué)物理方法數(shù)值化教學(xué)的新方法，是數(shù)學(xué)物理方法課程教學(xué)中的一項(xiàng)重要工作，也是數(shù)學(xué)物理方法教學(xué)改革中的重要內(nèi)容。利用MATLAB數(shù)值求解數(shù)學(xué)物理方程，將傳統(tǒng)教學(xué)手段與計(jì)算機(jī)仿真教學(xué)相結(jié)合，改變只用公式符號(hào)教學(xué)的模式[1]，令學(xué)生對復(fù)雜、抽象、煩瑣的數(shù)學(xué)物理問題具有更深刻的理解。本論文旨在進(jìn)行數(shù)學(xué)物理方程仿真求解實(shí)踐訓(xùn)練，著力培養(yǎng)大學(xué)生應(yīng)用數(shù)學(xué)物理思想解決實(shí)際問題的能力。本著“重理論、強(qiáng)實(shí)踐、突創(chuàng)新”的教育理念，結(jié)合科技前沿，以光子晶體的電磁場理論作為實(shí)踐內(nèi)容，利用MATLAB對復(fù)雜的電磁場本征值問題進(jìn)行計(jì)算機(jī)仿真求解，將結(jié)果三維可視化，以此來展現(xiàn)復(fù)雜電磁場問題的物理圖像，對培養(yǎng)大學(xué)生創(chuàng)新能力具有重要意義。

二、光子晶體電磁理論基礎(chǔ)

在利用分離變量法求解數(shù)學(xué)物理方程時(shí)，最后都?xì)w結(jié)到求解本征值問題。在用本征函數(shù)系展開法解數(shù)學(xué)物理方程時(shí)，也要對所用的本征函數(shù)系有較好的理解[2]。所以，各種本征函數(shù)系在數(shù)學(xué)物理方程課程的學(xué)習(xí)中有非常重要的地位。周期結(jié)構(gòu)對電磁波的調(diào)控是物理學(xué)領(lǐng)域的基礎(chǔ)問題。光子晶體是由介電常數(shù)周期排列形成的一種合成材料，是非均勻介質(zhì)中少數(shù)可以嚴(yán)格遵循電磁理論的新型人工材料。在一定的晶格常數(shù)和介電常數(shù)條件下，布拉格散射使在光子晶體中傳播的電磁波受到調(diào)制形成類似于電子的能帶結(jié)構(gòu)[3]。利用計(jì)算機(jī)仿真求解光子晶體中的復(fù)雜本征值問題，可以幫助學(xué)生熟悉并更好地掌握本征函數(shù)系的性質(zhì)和求解方法。

1.理想二維光子晶體的結(jié)構(gòu)。假設(shè)介電常數(shù)為εa，半徑為r的介質(zhì)柱平行于z軸，背景介質(zhì)的介電常數(shù)為εb，在（x，y）平面內(nèi)的晶格常數(shù)為a，θ為相鄰基矢a1和a2之間的銳角，當(dāng)θ=90°和60°時(shí)，分別為正方形晶格和正三角形晶格。（x，y）平面的傅里葉變換空間為倒易空間（如圖1所示），對應(yīng)于由波矢k定義的頻譜。

根據(jù)幾何關(guān)系，（x，y）平面內(nèi)的基矢（a1，a2）和倒易空間的基矢（b1，b2）分別為a1=aex，a2=a（cosθex+sinθey）；b1=■ex-■ctgθey，；b2=■ey。

在倒易空間中，Γ、T、N、X和M等點(diǎn)為布里淵區(qū)的高對稱點(diǎn)，它們所構(gòu)成的多邊形區(qū)域（深灰色部分）稱為不可約布里淵區(qū)。不可約布里淵區(qū)是倒易空間中最小的、可重復(fù)的區(qū)域，可以映射出電磁波在整個(gè)光子晶體中的傳輸特性[4]。對稱點(diǎn)坐標(biāo)分別為：Γ=（0，0），T=■（1，■），N=■（1，■ctgθ，■）；X=■（0，■），M=■（■ctgθ-1，■）。

2.理想二維光子晶體的本征值問題。

平面波的指數(shù)形式表示為

H（r，t）=H（r）e-iωt （1）

E（r，t）=H（r）e-iωt （2）

聯(lián)立無源Maxwell方程組，分別得到電場和磁場的傳播方程

？塄×■？塄×H（r）=■■H（r） （3）

■？塄×？塄×E（r）=■■E（r） （4）

ε（r）是光子晶體介質(zhì)分布的周期函數(shù)，本征值方程（3）和（4）式與電子材料中的周期性勢場問題的Schr？觟dinger方程類似，稱為光子晶體的支配方程[5]。本征場H（r）和E（r）分別對應(yīng)于理想二維光子晶體中橫磁（TM）模式和橫電（TE）模式的空間形態(tài)，通過求解本征值（ω/c）2，可以得到頻率ω與波矢k之間的色散關(guān)系，即光子晶體的能帶結(jié)構(gòu)。

3.光子晶體中的平面波展開。

根據(jù)Bloch理論，將光子晶體本征場用平面波展開為

HK（r）=■H■（G）e■ （5）

EK（r）=■E■（G）e■ （6）

G為倒格矢，將（5）和（6）式分別代入本征值方程（3）和（4）式，利用平面波基{G，exp[i（k+G）gr]，…}的正交性[6]，得到如下關(guān)于電磁場展開系數(shù)的本征值方程

■k%（G-G'）（k+G）g（k+G'）Hk（G'）=■Hk（G） （7）

■k%（G-G'）（k+G'）g（k+G'）Ek（G'）=■Ek（G） （8）

矩陣k%（G）是k（r）=1/ε（r）=k%（G）e■的傅立葉展開系數(shù)

k%（G）=■■k%（r）e■dr■=■■■e■dr2+■（■-■）■e■dr2 （9）

u表示一個(gè)周期單元，Au為周期單元橫截面的面積。c表示一個(gè)散射單元橫截面上的積分邊界。（9）式右邊包含了G≠0和G=0的項(xiàng)

k%（G）=■+（■-■）fr，G=02fr（■-■）■，G≠0 （10）

其中J1（GR）為第一類貝塞爾函數(shù)，fr為填充比。

三、仿真求解電磁場本征值問題

我們通過計(jì)算機(jī)仿真求解TM模式電磁場本征值方程（7）式，獲得二維菱形晶格光子晶體的本征頻率ωk與波矢k之間的色散關(guān)系，繪制出能帶曲線。

1.光子晶體的數(shù)學(xué)建模。

對于θ=70°的二維菱形晶格光子晶體，背景介質(zhì)的介電常數(shù)為εb=12，空氣柱的半徑r=0.4a。仿真步驟和MATLAB程序如下：

①定義光子晶體的結(jié)構(gòu)參數(shù)。

ea=1；eb=12；R=0.4；sita=70；a=1；a1=a*[1，0]；

a2=a*[cos（sita*pi/180），sin（sita*pi/180）]；

b1=2*pi/a*[1，-1/tan（sita*pi/180）]；b2=2*pi/a*[0，1/sin（sita*pi/180）]；

fr=pi*R*R/abs（a1（1）*a2（2）-a1（2）*a2（1））；

②定義倒易空間中對稱點(diǎn)的坐標(biāo)。

Point（1，：）=[0，0]；

Point（2，：）=pi/a*[1，（1-cos（sita））/sin（sita）]；Point（3，：）=pi/a*[1-（1-cos（sita））/sin（sita）*1/

tan（sita），1/sin（sita）]；

Point（4，：）=pi/a*[0，1/sin（sita）]；

Point（5，：）=pi/a*[（1-cos（sita））/sin（sita）*1/

tan（sita）-1，1/sin（sita）]；Point（6，：）=[0，0]；

③產(chǎn)生一個(gè)20×20的矩陣，確定平面波的波數(shù)NPW，定義倒格矢G。

DimForG=20+1；

NPW=DimForG*DimForG；

gtemp=-10：10；

gtemp1=repmat（gtemp，DimForG，1）；

Gx=b1（1）*gtemp1+b2（1）*gtemp1'；

Gy=b1（2）*gtemp1+b2（2）*gtemp1'；

Gx=Gx（：）'；Gy=Gy（：）'；

Gx_m=repmat（Gx，NPW，1）；Gx_n=Gx_m'；

Gy_m=repmat（Gy，NPW，1）；Gy_n=Gy_m'；

G=sqrt（（Gx_m-Gx_n）.*（Gx_m-Gx_n）+（Gy_m-Gy_n）.*（Gy_m-Gy_n））；

④確定κ（r）=1/ε（r）的傅里葉展開系數(shù)。

ek0=fr/ea+（1-fr）/eb；ekc=（1/ea-1/eb）*fr*2；

GR=G*R；na=find（GR==0）；GR（na）=1；

ek=ekc*besselj（1，GR）./GR；ek（na）=ek0；

2.仿真計(jì)算光子晶體TM模式能帶曲線。

①定義倒易空間波矢路徑。

用Keach代表波矢路徑上的取值密度，Ktype為對稱點(diǎn)的數(shù)目，第一布里淵區(qū)內(nèi)沿波矢路徑Γ→T→N→M→Γ的仿真程序?yàn)椋?/p>

Ktype=5；Keach=6；x=[]；y=[]；

for n=1：Ktype x1=linspace（Point（n，1），Point（n+1，1），Keach+1）；y1=linspace（Point（n，2），Point（n+1，2），Keach+1）；x=[x，x1（1：Keach）]；y=[y，y1（1：Keach）]；

end

②求解本征值方程。

eigval=[]；

for m=1：Ktype*Keach

kx=x（m）；ky=y（m）；

KGmn=sqrt（（kx-Gx_m）.^2+（ky-Gy_m）.^2）.*

sqrt（（kx-Gx_n）.^2+（ky-Gy_n）.^2）；

H=KGmn.*ek；

eigvalue=sort（eig（H））；

eigval=[eigval，eigvalue（1：20）]；end

eigval=[eigval，eigval（：，1）]；

eigval=real（sqrt（eigval）*a*0.5/pi）；

③繪制二維能帶曲線。

for m=1：Ktype

D（m）=sqrt（（Point（m+1，1）-Point（m，1））^2+（Point（m+1，2）-Point（m，2））^2）；

xtemp（m，：）=linspace（0，D（m），Keach+1）；

end

x=xtemp（1，1：Keach）；Dtotal=0；

for m=2：Ktype

Dtotal=Dtotal+D（m-1）；

x=[x，xtemp（m，1：Keach）+Dtotal]；

end

x=[x，xtemp（Ktype，Keach+1）+Dtotal]；

x=x/max（x）；plot（x，eigval）；

修飾過后的二維菱形晶格光子晶體TM偏振模式能帶曲線如圖2所示。

四、本征值函數(shù)的三維可視化仿真

繪制三維等頻面，關(guān)鍵是建立波矢平面（kx，ky）內(nèi)二維點(diǎn)陣的坐標(biāo)，再求解出每個(gè)點(diǎn)對應(yīng)特征值，仿真步驟和MATLAB程序?yàn)椋?/p>

1.定義波矢（kx，ky）平面內(nèi)點(diǎn)陣的坐標(biāo)Keach=36；

x=linspace（Point（1，1），Point（2，1），Keach+1）；

y=linspace（Point（1，2），Point（4，2），Keach+1）；

x1=[-x（Keach+1：-1：2），x（1：Keach+1）]；

y1=[-y（Keach+1：-1：2），y（1：Keach+1）]；

2.求解本征值方程。

eigval=[]；

for m=1：2*Keach+1

kx=x1（m）；

for n=1：2*Keach+1

ky=y1（n）；

KGmn=sqrt（（kx-Gx_m）.^2+（ky-Gy_m）.^2）.*sqrt（（kx-Gx_n）.^2+（ky-Gy_n）.^2）；

H=KGmn.*ek；

eigvalue=sort（eig（H））；

eigval（m，n，：）=eigvalue（1：20）；

end end

eigvalue=real（sqrt（eigval）*a*0.5/pi）；

3.繪制前四個(gè)能帶的三維能帶曲面。

for k=1：4

surf（x，y，eigvalue（：，：，k））；

hold on

end

圖3為二維菱形晶格光子晶體的前四個(gè)能帶的TM偏振模式三維能帶曲面。

五、結(jié)論

本論文通過計(jì)算機(jī)仿真求解二維菱形晶格光子晶體的電磁場本征值問題，繪制出能帶曲線和三維能帶曲面。將傳統(tǒng)教學(xué)手段與計(jì)算機(jī)仿真教學(xué)相結(jié)合，對復(fù)雜的數(shù)學(xué)物理方法問題進(jìn)行三維可視化求解，著力培養(yǎng)大學(xué)生的創(chuàng)新思維和解決實(shí)際問題的能力。

參考文獻(xiàn)：

[1]楊華軍.數(shù)學(xué)物理方法與仿真[M].第2版.北京：電子工業(yè)出版社，2011：359-377.

[2]彭芳麟.數(shù)學(xué)物理方程的MATLAB解法與可視化[M].北京：清華大學(xué)出版社，2004：62-63.

[3]E.Yablonovitch.Inhibited Spontaneous Emission in Solid-State Physics and Electronics.Physical Review Letters，1987，58（20）：2059-2062.

[4]溫熙森.光子/聲子晶體理論與技術(shù)[M].北京：科學(xué)出版社，2006：115-116.

[5]John D.Joannopoulos，etc.Photonic Crystals Molding the Flow of Light.Second Edition.Princeton University Press，2008：8-10.

[6]葉衛(wèi)明.光子晶體導(dǎo)論[M].北京：科學(xué)出版社，2010：22-27.