comsol電磁場仿真案例

COMSOL公司是全球多物理場建模與仿真解決方案的提倡者和領(lǐng)導(dǎo)者，其旗艦產(chǎn)品COMSOL Multiphysics，使工程師和科學(xué)家們可以通過模擬，賦予設(shè)計理念以生命。它有無與倫比的能力，使所有的物理現(xiàn)象可以在計算機(jī)上完美重現(xiàn)。COMSOL的用戶利用它提高了手機(jī)的接收性能，利用它改進(jìn)醫(yī)療設(shè)備的性能并提供更準(zhǔn)確的診斷，利用它使汽車和飛機(jī)變得更加安全和節(jié)能，利用它尋找新能源，利用它探索宇宙，甚至利用它去培養(yǎng)下一代的科學(xué)家。

Multiphysics翻譯為多物理場，因此這個軟件的優(yōu)勢就在于多物理場耦合方面。多物理場的本質(zhì)就是偏微分方程組（PDEs），所以只要是可以用偏微分方程組描述的物理現(xiàn)象，COMSOL Multiphysics都能夠很好的計算、模擬、仿真。

COMSOL Multiphysics是一款大型的高級數(shù)值仿真軟件。廣泛應(yīng)用于各個領(lǐng)域的科學(xué)研究以及工程計算，模擬科學(xué)和工程領(lǐng)域的各種物理過程。

COMSOL Multiphysics是以有限元法為基礎(chǔ)，通過求解偏微分方程（單場）或偏微分方程組（多場）來實現(xiàn)真實物理現(xiàn)象的仿真，用數(shù)學(xué)方法求解真實世界的物理現(xiàn)象。

大量預(yù)定義的物理應(yīng)用模式，范圍涵蓋從流體流動、熱傳導(dǎo)、到結(jié)構(gòu)力學(xué)、電磁分析等多種物理場，用戶可以快速的建立模型。COMSOL中定義模型非常靈活，材料屬性、源項、以及邊界條件等可以是常數(shù)、任意變量的函數(shù)、邏輯表達(dá)式、或者直接是一個代表實測數(shù)據(jù)的插值函數(shù)等。

預(yù)定義的多物理場應(yīng)用模式， 能夠解決許多常見的物理問題。同時，用戶也可以自主選擇需要的物理場并定義他們之間的相互關(guān)系。當(dāng)然，用戶也可以輸入自己的偏微分方程（PDEs），并指定它與其它方程或物理之間的關(guān)系。

COMSOL Multiphysics力圖滿足用戶仿真模擬的所有需求，成為用戶的首選仿真工具。它具有用途廣泛、靈活、易用的特性，比其它有限元分析軟件強(qiáng)大之處在于，利用附加的功能模塊，軟件功能可以很容易進(jìn)行擴(kuò)展。

COMSOL的磁場仿真分析

真空中的平面電磁波

電磁場若要符合麥克斯韋方程，則該電磁場的場矢量E和B之間必存在某種聯(lián)系，且電場可以脫離電荷和電流單獨存在，并以有限的速度在空間傳播?，F(xiàn)在我們從麥克斯韋方程出發(fā)，從理論上分析存在于真空中的電磁場所具有的性質(zhì)。

我們討論不存在實物的真空，從t=0時刻起，該空間不存在傳導(dǎo)電流。沒有電荷流動意味著該空間或無電荷分布和電荷分布不隨時間變化。不隨時間變化的電荷產(chǎn)生的是靜態(tài)電場，而在這里我們不研究靜態(tài)場，故不妨假設(shè)電荷和電流都不存在，這樣的空間稱為自由空間。因自由空間中

其中（1）表示自由空間中的電場是無源場，電場線是無頭無尾的閉合曲線，或是從無限遠(yuǎn)處來，延伸到無限遠(yuǎn)去的曲線。（2）表示自由空間的磁場仍是無源場，磁場的磁感應(yīng)線仍是無頭無尾的閉合曲線。（3）表示自由空間的電場是有旋場，變化的磁場是渦旋的中心，在磁感強(qiáng)度變化的地方，周圍有閉合的電場線。（4）表示自由空間中的磁場是有旋場，只有變化的電場才是磁場的漩渦中心，在電場強(qiáng)度變化的地方，周圍存在閉合的磁感應(yīng)線。

我們不難從這四個積分式子中得到他們的微分形式

可以看出電矢量和磁矢量滿足的方程式具有相同的形式，式中的常量是恒正的，我們用另一個恒量表示之，令

上面兩個式子式沿z方向傳播的簡諧波。?和k是兩個常量，不能完全任意，它們的值必須保持這兩個列簡諧波是波動方程的解，E?和B?是電場波與磁場波的初相位。

二、 通過GUIDE進(jìn)行仿真

使用GUIDE畫出（14）（15）式的解析解的動畫圖形、解析解的瀑布圖形和數(shù)值解的圖形，（16）（17）

式的沿z方向傳播的簡諧平面電磁波。

打開GUIDE的設(shè)計窗口，在GUIDE的設(shè)計窗口中畫一個坐標(biāo)軸，三個單選按鈕，四個靜態(tài)文本框，三個編輯框和一個按鈕。并對這些進(jìn)行相應(yīng)的屬性設(shè)置，如圖4所示是一個已經(jīng)做好的fig文件。

在編程之前必須對三個單選按鈕進(jìn)行如下編程

function radiobutton1_Callback（hObject， eventdata， handles）

set（handles.radiobutton1，‘value’，1）;

set（handles.radiobutton2，‘value’，0）;

set（handles.radiobutton3，‘value’，0）;

function radiobutton2_Callback（hObject， eventdata， handles）

set（handles.radiobutton1，‘value’，0）;

set（handles.radiobutton2，‘value’，1）;

set（handles.radiobutton3，‘value’，0）;

function radiobutton3_Callback（hObject， eventdata， handles）

set（handles.radiobutton1，‘value’，0）;

set（handles.radiobutton2，‘value’，0）;

set（handles.radiobutton3，‘value’，1）;

這樣可以使單選按鈕選中其中的一個，其他的按鈕不起作用。

因為程序中有三套程序，每一個按鈕對應(yīng)著一套程序。選中一個按鈕只有一套程序起作用。這樣通過一個if…elseif…else來完成這項任務(wù)。

1、解析解的動畫圖形

我們可以看出（14）和（15）式在形式上是一樣的，只不過是方向不一樣，電矢量E和磁矢量B相互垂直。它們都與弦震動方程

可以用動畫來表現(xiàn)這個解，下面是一段由上式編寫的程序。

clear

a=1;l=1;

A=0.01;w=6;

x=0:0.5:1;

t=0:0.001:4.3;

［X，T］=meshgrid（x，t）;

u0=A*si

n（w*X./a）.*sin（w.*T）/sin（w*l/a）; u=0;

for n = 1:100;

uu=（-1）^（n+1）*sin（n*pi*X/l）.*sin（n*pi*a*T/l）/（w*w/a/a-n*n*pi*pi/l/l）;

u=u+uu; end

u=u0+2*A*w/a/l.*u; figure（1）

axis（［0，1，-0.05，0.05］）

h=plot（x，u（1，：），‘linewidth’，3）;

set（h，‘erasemode’，‘xor’）;

for j=2:length（t）;

set（h，‘ydata’，u（j，：））;

axis（［0，1，- 0.05，0.05］）

drawnow end figure（2）

waterfall（X（1:50:3000，：），T（1:50:3000，：），u（1:50:3000，：））

xlabel（‘x’）

ylabel（‘t’）

圖1是動畫中的一幅畫面，可以看出，固定端一直保持不動，而作諧振動的端使弦也產(chǎn)生了一種在傳播的振動。這相當(dāng)于將一根繩的一端固定，而周期性的擺動另一端時，在繩上所產(chǎn)生的運動。［3］

2、數(shù)值解的圖形

我們再用微分方程工具箱求解上面的問題。我們用二維圖形來表示一維的弦，讓弦的側(cè)面不受力，左端固定，右端作受迫振動。

在Options/Axes limits下選擇x軸范圍為0～1，y軸范圍為0～1.以原點為頂點畫一個場為1寬為0.4的矩形，矩形的頂點為（0，0），（1，0），（1，0.4），（0，0.4）。

按照題意，矩形的右邊界是齊次的狄里克利邊界條件，可取h=1，r=0，左邊界是非齊次的狄里克利邊界條件，可取h=1，r=0.01*sin（6*t），上下邊界則取齊次的諾依曼邊界條件，即g=0，q=0。

方程的設(shè)置是hyperbolic型，系數(shù)是c=1，a=0，f=0，d=1. 為了有足夠的精度，初始化的網(wǎng)格要再作兩次細(xì)分。

在解方程的參數(shù)設(shè)置對話框Solve Parameters中，各項選擇如下：在Plot type下，選Color和Height（3-D Plot），在Property下，對應(yīng)的位置中都選User entry，在User entry下，再在相應(yīng)的位置都輸入10*u，這樣做的目的是，弦的振幅太小，為了達(dá)到更好的顯示效果，所以將振幅放大10倍來畫圖。所得的圖形如

3、沿z方向傳播的簡諧平面電磁波

通過（16）（17）式進(jìn)行如下編程

a = str2double （get（handles.t， ‘string’））;

for t=1:a;

v=2;

w= str2double （get（handles.w， ‘string’））;

k=0:1：40;

y=sin（w*（t-k/v））;

x=zeros（1，41）;

z=k stem3（z，x，y，‘r’）;

hold on x=sin（w*（t-k/v））;

y=zeros（1，41）;

z=k;

stem3（z，x，y，‘b’，‘filled’）;

view（-37.50，30）;

hold off w1=moviein（a）;

w1（：，t）=getframe;

xlabel（‘t’）;

ylabel（‘E’）;

zlabel（‘B’）;

title（‘沿Z方向傳播的簡諧平面電磁波’）

end

這時可以得到function pushbutton1_Callback（hObject， eventdata， handles）下面的程序

if （findobj（‘tag’，‘radiobutton1’，‘value’，1））;

a=1;

l=1;

A= str2double （get（handles.b， ‘string’））;

w= str2double （get（handles.w， ‘string’））;;

x=0:0.5:1; t=0:0.001:4.3;

［X，T］=meshgrid（x，t）;

u0=A*sin（w*X./a）.*sin（w.*T）/sin（w*l/a）;

u=0;

for n = 1:100;

uu=（-1）^（n+1）*sin（n*pi*X/l）.*sin（n*pi*a*T/l）/（w*w/a/a-n*n*pi*pi/l/l）;

u=u+uu; end

u=u0+2*A*w/a/l.*u; figure（1）

axis（［0，1，-0.05，0.05］）

h=plot（x，u（1，：），‘linewidth’，3）;

set（h，‘erasemode’，‘xor’）;

for j=2:length（t）;

set（h，‘ydata’，u（j，：））;

axis（［0，1，- 0.05，0.05］） drawnow end figure（2）

waterfall（X（1:50:3000，：），T（1:50:3000，：），u（1:50:3000，：））

xlabel（‘x’）

ylabel（‘t’）

elseif （findobj（‘tag’，‘radiobutton2’，‘value’，1））;

% function pdemodel

［pde_fig，ax］=pdeinit;

pdetool（‘a(chǎn)ppl_cb’，1）;

set（ax，‘DataAspectRatio’，［1 1.5 1］）;

set（ax，‘PlotBoxAspectRatio’，［1 0.66666666666666663 2］）;

set（ax，‘XLim’，［0 1］）; set（ax，‘YLim’，［0 1］）;

set（ax，‘XTickMode’，‘a(chǎn)uto’）;

set（ax，‘YTickMode’，‘a(chǎn)uto’）;

% Geometry description：

pderect（［0 1 0.40000000000000002 0］，‘R1’）;

set（findobj（get（pde_fig，‘Children’），‘Tag’，‘PDEEval’），‘String’，‘R1’）

% Boundary conditions：

pdetool（‘changemode’，0）

pdesetbd（4，。。。

‘dir’，。。。

1，。。。

‘1’，。。。

‘0.01*sin（6*t）’）

pdesetbd（3，。。。

‘neu’，。。。

1，。。。

‘0’，。。。

‘0’）

pdesetbd（2，。。。 ‘dir’，。。。 1，。。。 ‘1’，。。。 ‘0’）

pdesetbd（1，。。。 ‘neu’，。。。 1，。。。 ‘0’，。。。 ‘0’）

% Mesh generation：

setuprop（pde_fig，‘Hgrad’，1.3）;

setuprop（pde_fig，‘refinemethod’，‘regular’）; pdetool（‘initmesh’） pdetool（‘refine’） pdetool（‘refine’） % PDE coefficients： pdeseteq（3，。。。 ‘1.0’，。。。 ‘0.0’，。。。 ‘0.0’，。。。 ‘1.0’，。。。 ‘0:0.1:1.5’，。。。 ‘0.0’，。。。 ‘0.0’，。。。 ‘［0 100］’）

setuprop（pde_fig，‘currparam’，。。。 ［‘1.0’;。。。 ‘0.0’;。。。 ‘0.0’;。。。 ‘1.0’］）

% Solve parameters：

setuprop（pde_fig，‘solveparam’，。。. str2mat（‘0’，‘3168’，‘10’，‘pdeadworst’，。。. ‘0.5’，‘longest’，‘0’，‘1E-4’，‘’，‘fixed’，‘Inf’）） % Plotflags and user data strings：

setuprop（pde_fig，‘plotflags’，［1 1 4 1 1 1 1 1 0 0 1 16 1 0 1 0 0 1］）; setuprop（pde_fig，‘colstring’，‘’）; setuprop（pde_fig，‘arrowstring’，‘’）; setuprop（pde_fig，‘deformstring’，‘’）; setuprop（pde_fig，‘heightstring’，‘10*u’）; % Solve PDE： pdetool（‘solve’）

else （findobj（‘tag’，‘radiobutton3’，‘value’，1））; a = str2double （get（handles.t， ‘string’））; for t=1:a; v=2;

w= str2double （get（handles.w， ‘string’））; k=0:1：40; y=sin（w*（t-k/v））; x=zeros（1，41）; z=k stem3（z，x，y，‘r’）; hold on x=sin（w*（t-k/v））; y=zeros（1，41）; z=k;

stem3（z，x，y，‘b’，‘filled’）; view（-37.50，30）; hold off w1=moviein（a）; w1（：，t）=getframe; xlabel（‘t’）; ylabel（‘E’）; zlabel（‘B’）;

title（‘沿Z方向傳播的簡諧平面電磁波’） end

end grid on

通過運行GUIDE就會得道圖1、圖2、圖3和圖4.

三、 結(jié)論

MATLAB有強(qiáng)大的求解偏微分方程和可視化功能模擬各類物理場的實驗是成功的。借助偏微分方程工具箱，可以通過分析靜電學(xué)、電場、磁場和電磁場的原理而建立微分方程，經(jīng)過數(shù)值計算模擬各類電磁場問題。使問題更加形象逼真。 利用這些特性及GUI 功能可以實現(xiàn)交互式數(shù)據(jù)處理。圖形用戶界面的設(shè)計確定了應(yīng)用程序的主要框架和基本功能，完成了窗口、圖標(biāo)、按鈕等用戶界面，軟件開發(fā)者只需在由軟件發(fā)工具自動生成的程序代碼中添加自己的運算或控制代碼，就可以完成自己的設(shè)計。