時(shí)域有限差分方法林志立zllin2008@北航儀器光電學(xué)院光電工程系TheFinite-DifferenceTime-DomainMethodforElectromagnetics計(jì)算電磁學(xué)中的電磁學(xué)中幾種重要的數(shù)值計(jì)算方法有限差分法FiniteDifferenceMethod靜電場、靜磁場的有限差分法；時(shí)域行波的電磁場的時(shí)域有限差分法；有限元法(FiniteElementMethod)數(shù)值求解各類獨(dú)立的偏微分方程；（電磁學(xué)、材料力學(xué)、工程熱力學(xué)、聲學(xué)等等）矩量法(MethodofMoments)

適合于細(xì)線、平面形狀結(jié)構(gòu)的電磁場問題電磁學(xué)基本方程麥克斯韋方程組(MaxwellEquations)JamesClerkMaxwell(1831–1879)(法拉第感應(yīng)定律)(安培環(huán)路定律)(高斯定律－電場)(高斯定律－磁場)物質(zhì)本構(gòu)關(guān)系(ConstitutiveRelations)（歐姆定律）（電極化）（磁化）一切電磁場和電磁波問題均可由以上方程，以及各類具體的邊界條件所決定！麥克斯韋方程組中的運(yùn)算符散度（Divergence)旋度（Curl)連續(xù)函數(shù)的偏微分運(yùn)算設(shè)有一連續(xù)函數(shù),現(xiàn)欲求。FDTD的基本思想－連續(xù)偏微分的有限階近似－時(shí)域和空間域的離散化二階中心差分近似表達(dá)式：當(dāng)越小時(shí)，上式的近似程度越高。實(shí)際上：FDTD空間域的離散化(1)空間域的分割離散化Ex分量的空間離散分布圖Yee元胞（Δx,Δy,Δz)Hx分量的空間離散分布圖節(jié)點(diǎn)FDTD空間域的離散化YEE元胞各電磁場分量在元胞中的位置K.S.Yee,“NumericalsolutionofinitialboundaryvalueproblemsinvolvingMaxwell’sequationsinisotropicmedia,”IEEETrans.AntennasPropagat.,vol.14,1966,pp.302-307.例如：Hz為Ex和Ey所環(huán)繞。FDTD空間偏微分的近似以Hz為例：類似地，可實(shí)現(xiàn)各電磁場分量的空間偏微分計(jì)算。FDTD時(shí)間偏微分的近似以Hz為例：t=(n+1/2)Δtt=nΔtt=(n-1/2)Δt時(shí)間上的推移(1)電場在時(shí)間上取整數(shù)倍的Δt;

t=n*Δt;(2)磁場在時(shí)間上取(整數(shù)＋1/2)倍的Δt;

t=(n+1/2)*Δt;麥克斯韋方程的離散化近似以Hz為例：上式即為Hz的更新方程，由前一時(shí)刻的磁場和前半時(shí)刻的臨近空間格點(diǎn)的電場即可求出最新時(shí)刻的磁場。麥克斯韋方程的離散化近似采取類似的步驟，可以推導(dǎo)出其它場量的更新表達(dá)式：例如，對于Ez：FDTD的離散參數(shù)的選擇★元胞尺寸：邊長小于最短波長的1/10,以減小數(shù)值色散。數(shù)值色散方程：理想色散方程：要求：例如，取例如，取FDTD的離散參數(shù)的穩(wěn)定性條件★

時(shí)間步長：Courant穩(wěn)定性條件對于非色散介質(zhì)，時(shí)間步長不能大于以下表達(dá)式：（vonNeumannmethod）為了保持穩(wěn)定性，該方程的所有解的模必須小于1。Ｚ域數(shù)值色散方程：介質(zhì)電磁參量的設(shè)定不同的元胞的電磁參量應(yīng)設(shè)置為所在空間所代表的介質(zhì)的介電常數(shù)和磁導(dǎo)率。長方體，物質(zhì)I長方體，物質(zhì)II球體，物質(zhì)III空氣場量與介質(zhì)參數(shù)要對應(yīng)色散介質(zhì)的FDTD模擬以Lorentz介質(zhì)為例：更新方程：近似，求系數(shù)Z-變換時(shí)域：頻域：色散介質(zhì)的FDTD模擬

模擬Lorentz色散介質(zhì)的不同方法：MSEapproach:FDTD編程流程主循環(huán)初始化輸出結(jié)果編程舉例1：一維FDTD問題x基本旋度方程：X向電導(dǎo)率X向磁電導(dǎo)率X向電流X向磁流編程舉例1：一維FDTD問題Matlab程序代碼：%Defineinitialconstantseps_0=8.854187817e-12;%permittivityoffreespacemu_0=4*pi*1e-7;%permeabilityoffreespacec=1/sqrt(mu_0*eps_0);%speedoflight

%Defineproblemgeometryandparametersdomain_size=1;%1Dproblemspacelengthinmetersdx=1e-3;

%cellsizeinmetersdt=3e-12;

%durationoftimestepinsecondsnumber_of_time_steps=2000;%numberofiterationsnx=round(domain_size/dx);%numberofcellsin1Dproblemspacesource_position=0.5;%positionofthecurrentsourceJz

1.定義物理常量2.定義問題的參量和結(jié)構(gòu)尺寸%InitializefieldandmaterialarraysCeze=zeros(nx+1,1);Cezhy=zeros(nx+1,1);Cezj=zeros(nx+1,1);Ez=zeros(nx+1,1);Jz=zeros(nx+1,1);eps_r_z=ones(nx+1,1);%freespacesigma_e_z=zeros(nx+1,1);%freespace

Chyh=zeros(nx,1);Chyez=zeros(nx,1);Chym=zeros(nx,1);Hy=zeros(nx,1);My=zeros(nx,1);mu_r_y=ones(nx,1);%freespacesigma_m_y=zeros(nx,1);%freespace

編程舉例1：一維FDTD問題（續(xù)）3.初始化場量和介質(zhì)參量陣列電場電流部分磁場磁流部分%CalculateFDTDupdatingcoefficientsCeze=(2*eps_r_z*eps_0-dt*sigma_e_z)..../(2*eps_r_z*eps_0+dt*sigma_e_z);

Cezhy=(2*dt/dx)..../(2*eps_r_z*eps_0+dt*sigma_e_z);

Cezj=(-2*dt)..../(2*eps_r_z*eps_0+dt*sigma_e_z);

Chyh=(2*mu_r_y*mu_0-dt*sigma_m_y)..../(2*mu_r_y*mu_0+dt*sigma_m_y);

Chyez=(2*dt/dx)..../(2*mu_r_y*mu_0+dt*sigma_m_y);

Chym=(-2*dt)..../(2*mu_r_y*mu_0+dt*sigma_m_y);

編程舉例1：一維FDTD問題（續(xù)）4.計(jì)算更新方程系數(shù)電場部分磁場部分%DefinetheGaussiansourcewaveformtime=dt*[0:number_of_time_steps-1].';Jz_waveform=exp(-((time-2e-10)/5e-11).^2);source_position_index=round(nx*source_position/domain_size)+1;

%Subroutinetoinitializeplottinginitialize_plotting_parameters;

編程舉例1：一維FDTD問題（續(xù)）5.定義場源源波形為高斯型6.作圖初始化Ez_positions=[0:nx]*dx;Hy_positions=([0:nx-1]+0.5)*dx;v=[0-0.1-0.1;0-0.10.1;00.10.1;00.1-0.1;...1-0.1-0.1;1-0.10.1;10.10.1;10.1-0.1];f=[1234;5678];axis([01-0.20.2-0.20.2]);lez=line(Ez_positions,Ez*0,Ez,'Color','b','LineWidth',1.5);lhy=line(Hy_positions,377*Hy,Hy*0,'Color','r',...'LineWidth',1.5,'linestyle','-.');set(gca,'fontsize',12,'FontWeight','bold');axissquare;legend('E_{z}','H_{y}\times377','Location','NorthEast');xlabel('x[m]');ylabel('[A/m]');zlabel('[V/m]');gridon;p=patch('vertices',v,'faces',f,'facecolor','g','facealpha',0.2);text(0,1,1.1,'PEC','horizontalalignment','center','fontweight','bold');text(1,1,1.1,'PEC','horizontalalignment','center','fontweight','bold');%FDTDloopfortime_step=1:number_of_time_steps

%UpdateJzforthecurrenttimestepJz(source_position_index)=Jz_waveform(time_step);

%UpdatemagneticfieldHy(1:nx)=Chyh(1:nx).*Hy(1:nx)...+Chyez(1:nx).*(Ez(2:nx+1)-Ez(1:nx))...+Chym(1:nx).*My(1:nx);

%UpdateelectricfieldEz(2:nx)=Ceze(2:nx).*Ez(2:nx)...+Cezhy(2:nx).*(Hy(2:nx)-Hy(1:nx-1))...+Cezj(2:nx).*Jz(2:nx);

Ez(1)=0;%ApplyPECboundaryconditionatx=0mEz(nx+1)=0;%ApplyPECboundaryconditionatx=1m

%Subroutinetoplotthecurrentstateofthefieldsplot_fields;end編程舉例1：一維FDTD問題（續(xù)）7.主循環(huán)程序引入電流源更新磁場更新電場設(shè)置完美電導(dǎo)層跟新圖面%subroutineusedtoplot1Dtransientfieldsdelete(lez);delete(lhy);lez=line(Ez_positions,Ez*0,Ez,'Color','b','LineWidth',1.5);lhy=line(Hy_positions,377*Hy,Hy*0,'Color','r',...'LineWidth',1.5,'linestyle','-.');ts=num2str(time_step);ti=num2str(dt*time_step*1e9);title(['timestep='ts',time='ti'ns']);drawnow;編程舉例1：一維FDTD問題（續(xù)）7.更新畫面運(yùn)行程序:fdtd_1d_code.m開放模擬空間的模擬完美匹配層（PerfectlyMatchedLayer)的概念吸波材料（損耗介質(zhì)）完美匹配層（PerfectlyMatchedLayer)材料的特點(diǎn)形成PML必須符合的兩個(gè)條件:從計(jì)算空間到PML層的阻抗必須一直匹配:2.某方向的電磁參量必須是另一個(gè)方向電磁參量的倒數(shù)完美匹配層（PerfectlyMatchedLayer)材料的特點(diǎn)下面的材料就符合這些條件：阻抗始終是相等的：電導(dǎo)率s隨著PML的層數(shù)號（x)而逐漸增加,提高吸收外行電磁波的能力。完美匹配層（PerfectlyMatchedLayer)材料的特點(diǎn)完美匹配層PML的驗(yàn)證實(shí)例中心有一光源，四周被PML環(huán)繞可以看到外形電磁波被完美吸收，沒有反射。運(yùn)行程序:fd2d_pml_code.mFDTD的應(yīng)用實(shí)例(1)電磁波在損耗介質(zhì)中的傳播(2)電磁波在色散介質(zhì)中的傳播損耗介質(zhì)：運(yùn)行程序:FD1D_Loss.m只有幅值的減小，沒有色散。同時(shí)有幅值的減小，以及色散。Debye色散介質(zhì)運(yùn)行程序:FD1D_Dispersion.mFDTD的應(yīng)用實(shí)例模擬線光源進(jìn)入平板玻璃光波導(dǎo)的過程玻璃平板波導(dǎo)線光源光在波導(dǎo)中的傳播光開始進(jìn)入光波導(dǎo)光波長:

0.8微米；平波厚度：2.4微米；折射率