comsol電磁場仿真案例

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

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

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

COMSOL Multiphysics是以有限元法為基礎，通過求解偏微分方程（單場）或偏微分方程組（多場）來實現真實物理現象的仿真，用數學方法求解真實世界的物理現象。

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

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

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

COMSOL的磁場仿真分析

真空中的平面電磁波

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

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

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

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

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

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

二、 通過GUIDE進行仿真

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

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

打開GUIDE的設計窗口，在GUIDE的設計窗口中畫一個坐標軸，三個單選按鈕，四個靜態文本框，三個編輯框和一個按鈕。并對這些進行相應的屬性設置，如圖4所示是一個已經做好的fig文件。

在編程之前必須對三個單選按鈕進行如下編程

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）;

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

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

1、解析解的動畫圖形

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

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

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是動畫中的一幅畫面，可以看出，固定端一直保持不動，而作諧振動的端使弦也產生了一種在傳播的振動。這相當于將一根繩的一端固定，而周期性的擺動另一端時，在繩上所產生的運動。［3］

2、數值解的圖形

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

在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。

方程的設置是hyperbolic型，系數是c=1，a=0，f=0，d=1. 為了有足夠的精度，初始化的網格要再作兩次細分。

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

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

通過（16）（17）式進行如下編程

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（‘appl_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’，‘auto’）;

set（ax，‘YTickMode’，‘auto’）;

% 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.

三、 結論

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