怎樣在femlab中寫(xiě)入自己及方程

1、怎么在femlab里寫(xiě)入自己的方程&I+n!C0A,F;a%s仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAMfemlab的自定義方程主要有三種形式：1、參數(shù)形式（coefficient form）)_"V2r6l.l/仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM2、普通形式（genera

2、l form）5|*K&v;E4R|S|仿真|設(shè)計(jì)|有限元|虛擬儀器3、弱解形式（weak form）1L+G9a,l-Q$i6N;O(o:    其中參數(shù)形式和普通形式十分類(lèi)似，參數(shù)形式主要解決線性問(wèn)題，而普通形式可以解決線性和弱的非線性問(wèn)題。而弱解形式主要用于解決非線性問(wèn)題以及一些無(wú)法用前兩者表達(dá)的線性問(wèn)題。弱解形式是功能最為強(qiáng)大的一種求解方法，前兩者可以解決的問(wèn)題都可以用弱解形式解決，只是要費(fèi)一些腦筋。    既然是初級(jí)講座就只講前兩者吧，因?yàn)槿踅庑问捷^為復(fù)雜，需要有有限元解法的一些基礎(chǔ)知識(shí)，而且使用上需要用到分步積分法。一、參數(shù)形式1、穩(wěn)

3、態(tài)問(wèn)題    穩(wěn)態(tài)問(wèn)題是指求解量不隨時(shí)間變化。對(duì)于靜力學(xué)，求解量確實(shí)是與時(shí)間無(wú)關(guān)的靜態(tài)量，例如位移、應(yīng)力、應(yīng)變等。而對(duì)于動(dòng)力學(xué)、波動(dòng)力學(xué)、電磁學(xué)等問(wèn)題，求解量不隨時(shí)間變化只是指其幅值與時(shí)間無(wú)關(guān)。    femlab指定的參數(shù)形式的方程如下（圖1）："p8| L/P7L%v#仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM  其中第一個(gè)式子是求解域上的方程。第二、三個(gè)式子是邊界條件，分別為第三類(lèi)

4、邊界條件和第一類(lèi)邊界條件，即羅賓邊界和狄利克邊界條件，兩者只能選取一個(gè)。注：羅賓邊界條件在femlab里常常被成為（廣義的）牛曼邊界條件。'a&a8s,Y9p仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent方程里各參量的意義，在不同的問(wèn)題里分別有不同的物理意義，但由于內(nèi)在數(shù)學(xué)形式的一致性，因此其意義也十分類(lèi)似，現(xiàn)以連續(xù)介質(zhì)力學(xué)即擴(kuò)散力學(xué)為例進(jìn)行說(shuō)明。c擴(kuò)散系數(shù)a吸收系數(shù)f源項(xiàng)保守通量對(duì)流系數(shù)對(duì)流系數(shù)保守通量源項(xiàng)q邊界上的吸收系數(shù)6P0p)n6Y)$O"jSi

5、mWe仿真論壇g邊界上的源項(xiàng)這里需要注意的是源項(xiàng)f，對(duì)于彈性力學(xué)，對(duì)應(yīng)于加在彈性體上的荷載。對(duì)于許多的物理場(chǎng)(聲場(chǎng)、電磁場(chǎng)、溫度場(chǎng))問(wèn)題，對(duì)應(yīng)于場(chǎng)源（聲源、點(diǎn)電荷、熱源）。而其他系數(shù)大都用于定義材料和媒質(zhì)的特性。其中、可以是矢量，c可以是矩陣，用于表示各項(xiàng)異性的材料特性?，F(xiàn)在給出一個(gè)簡(jiǎn)單地判斷問(wèn)題是線性還是非線性的方法，若這些參量是求解變量u的函數(shù)，則問(wèn)題是非線性的，反之為線性。當(dāng)然還有個(gè)別的問(wèn)題不符合這條準(zhǔn)則，但要準(zhǔn)確判斷一個(gè)問(wèn)題是否線性是比較困難的，詳細(xì)說(shuō)明可能要花上一章的功夫。呵呵，不說(shuō)了，大家就用這條簡(jiǎn)單的準(zhǔn)則判斷吧。.z!B,d9w5N0J8y3z|S|仿真|設(shè)計(jì)|有限元|虛擬儀器f

6、emlab指定的PDE的參數(shù)形式已經(jīng)包含了對(duì)空間坐標(biāo)的二階、一階、0階導(dǎo)數(shù)，因此大部分常見(jiàn)的線性PDE都可以轉(zhuǎn)換成femlab指定的參數(shù)形式。其方法是將求解變量對(duì)空間坐標(biāo)的導(dǎo)數(shù)按照階數(shù)的高低依次排列，然后對(duì)照寫(xiě)出相應(yīng)的參量。現(xiàn)舉一例如下：5W-?2Z*&2k2仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent例：穩(wěn)態(tài)下的聲學(xué)波動(dòng)方程如下：  其中為媒質(zhì)密度，cs為聲速，為聲波角頻率，p為聲壓。因此對(duì)比f(wàn)emlab的參數(shù)形式方程可以得出相應(yīng)的參量的值為：7Y2k3k)G

7、)s仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAMc=-1/,a=sup2/sup/cssup2/sup,其余的參量為0。具體操作步驟如下：（1）在Model Navigator 窗口里選擇PDE modes-PDE coefficient form-stationary analysis設(shè)置求解變量的名稱為p，如下圖3：.O(q)!u3R9A+Z;n2N4?9v仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALG

8、OR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM（2）畫(huà)出求解域圖形，這里我用的是長(zhǎng)為2，寬為1.2的矩形，中心在原點(diǎn)。選擇physicssubdomain settings，在c一欄出填入-1/rho；在a一欄處填入omega2/rho/cs2，如下圖4：/m"&q;v2%M:E&仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM（3）在optionsconstants處

9、根據(jù)實(shí)際情況分別填入rho、omega、cs的值，如圖5：圖5：（4）最后可以設(shè)置邊界條件求解了。假設(shè)邊界條件是聲學(xué)硬邊界，即p/n=0。此為牛曼邊界，且g=q=0。在physicsboundary setting中選擇要設(shè)置的邊界，并選擇Neumann boundary condition，設(shè)置q、g都為0。如圖6：3+X)J1u#Y64:_4d1仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM圖6：  下表列出了各邊界的邊界條件，大家可以

10、根據(jù)步驟（4）的方法依次設(shè)置各邊界。1_$8Y%B/z:n&I5F/|SimWe仿真論壇邊界  邊界條件  類(lèi)型  設(shè)置項(xiàng)1     p/n0.05    neumann   q=0,g=0.05)b#V/R3d./e仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent2,3  p/n0    neum

11、ann  q=0,g=04     p/ni*p/100  neumann  q=-i*omega/100,g=0'Y*?4j(q9R$V8T9j4c仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM設(shè)置好后就可以網(wǎng)格化，求解了。p的幅值求解結(jié)果如圖7：1I1c+f#g:g3P;w9p仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MS

12、C,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent（5）現(xiàn)在給一個(gè)復(fù)雜一點(diǎn)的情況。%l0R5z99e4b!8N'h其實(shí)方程和邊界條件里的各參量都可以是空間坐標(biāo)的函數(shù)，對(duì)于本例，假設(shè)媒質(zhì)的密度是x的函數(shù)，即（x），原方程變?yōu)椋▓D8）：  這里假設(shè)（x）1.2*abs（x）'U!W'm*X"a*sSimWe仿真論壇我們?cè)黾右粋€(gè)名為rhox的變量用以代替原來(lái)的rho常量。在option-expressions-scalar expressions中加入一個(gè)名為rhox、表達(dá)式為rho*abs（x）的變量。圖9.

13、JPG(22.74 KB)然后在subdomain settings的頁(yè)面中，出現(xiàn)rho的地方全部用rhox代替。邊界條件中各參量也可以是空間坐標(biāo)的函數(shù)，假設(shè)邊界1的條件變?yōu)椋簆/n0.05*y，則只要在boundary setting中將 g一項(xiàng)設(shè)置為0.05*y即可。更改后的求解結(jié)果如圖10，注意圖中顯示的是幅值，即abs（p）。圖10：  2、參數(shù)形式的動(dòng)態(tài)問(wèn)題。-U9? r:f58zSimWe仿真論壇所謂動(dòng)態(tài)問(wèn)題是指求解量與時(shí)間相關(guān)，反映在方程中就是含有時(shí)間t的項(xiàng)，通常是關(guān)于時(shí)間的一階或二階的導(dǎo)數(shù)。  比較穩(wěn)態(tài)問(wèn)題和動(dòng)態(tài)問(wèn)題的方程，可以發(fā)現(xiàn)后者只是比前者多了一個(gè)ds

14、uba/sub*u/t 一項(xiàng)，其中的dsuba/sub被稱為質(zhì)量系數(shù)，這是因?yàn)檫@個(gè)系數(shù)通常都要被研究物體的質(zhì)量或密度有關(guān)。當(dāng)然，很多的動(dòng)態(tài)問(wèn)題都含有對(duì)時(shí)間t的二階導(dǎo)數(shù)的項(xiàng)（例如波動(dòng)方程）其實(shí)含有時(shí)間二階導(dǎo)數(shù)項(xiàng)的方程，可以適當(dāng)轉(zhuǎn)換為方程組，使每個(gè)方程都只含有時(shí)間一階導(dǎo)數(shù)項(xiàng)。例如圖12所示的方程就可以通過(guò)引入一個(gè)中間變量v，從而建立一個(gè)只含時(shí)間一階導(dǎo)數(shù)項(xiàng)的方程組建立方程組可在model navigator窗口的dependent viabless一欄中分別填入方程組的所有變量名， -R:e"+R9"X i仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABA

15、QUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM注意填寫(xiě)方程組的參量時(shí)，各參量都應(yīng)該是矩陣或數(shù)列的形式，會(huì)比較復(fù)雜。當(dāng)然femlab已經(jīng)內(nèi)置了含時(shí)間二階導(dǎo)數(shù)項(xiàng)方程的模式，稱為"time-dependent analysis,wave extension"，大家可以直接使用和以之為參考。,s,U.E"D8r5-X%9Z8s|S|仿真|設(shè)計(jì)|有限元|虛擬儀器當(dāng)然，我們也可以用弱解形式來(lái)求解這類(lèi)問(wèn)題，在這里我就不談了。._)q%Q*#5Z'E3h|S|仿真|設(shè)計(jì)|有限元|虛擬儀器3、特征值問(wèn)題

16、所謂特征值問(wèn)題，其方程的形式與穩(wěn)態(tài)問(wèn)題的方程相比較，多了 一項(xiàng)*dsuba/sub*u。其中就是特征值，由特征值我們還可以算出特征頻率。求解特征值問(wèn)題可以得到各特征值的值，以及不同特征值所對(duì)應(yīng)的變量u的值。/T/?5R)J'v|S|仿真|設(shè)計(jì)|有限元|虛擬儀器二、普通形式普通形式和參數(shù)形式十分類(lèi)似，也同樣分為穩(wěn)態(tài)、動(dòng)態(tài)、特征值問(wèn)題。為簡(jiǎn)練起見(jiàn)，以下僅以穩(wěn)態(tài)問(wèn)題為例。普通形式的穩(wěn)態(tài)問(wèn)題的方程如下：比較參數(shù)形式和普通形式的方程，可以發(fā)現(xiàn)存在以下關(guān)系：=-cu-*u+F=f-u-a*u:&0s8F"u1O$O/Z't仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,A

17、NSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAMG=g-q*u5U!n8z$5W7v-v:c(L4)eSimWe仿真論壇R=r-h*u'_(z)d74A)W+uSimWe仿真論壇因此，參數(shù)形式的方程都可以通過(guò)以上公式轉(zhuǎn)換為普通形式的方程。+?)%T+b-t|S|仿真|設(shè)計(jì)|有限元|虛擬儀器例如對(duì)于前面的第一個(gè)例子，因?yàn)閏=-1/,a=sup2/sup/cssup2/sup，所以=u/，F(xiàn)=-sup2/sup*u/cssup2/sup6?3c$l2g _+A仿真分析,有限元,模擬,計(jì)算,力學(xué),航

18、空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM在femlab中的設(shè)置如圖17。2s-N3p/D5g9W仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent,CFD,CAE,CAD,CAM邊界條件同樣可以根據(jù)以上公式進(jìn)行轉(zhuǎn)換，邊界14的設(shè)置見(jiàn)下表：邊界  類(lèi)型  g  1 Neumann 0.050.j/c#N;I+h2 Neumann 03 Neumann 066R)k-n(U-h0T4 Neumann i*omega/100*u%n3I18C!W-e$5d仿真分析,有限元,模擬,計(jì)算,力學(xué),航空,航天,ANSYS,MSC,ABAQUS,ALGOR,Adina,COMSOL,FEMLAB,Matlab,Fluent圖18為邊界4的設(shè)置，其他邊界類(lèi)似，略。然后是劃分網(wǎng)格和求解，求解結(jié)果與圖7完全一致。/k.0G 97o y$F*H#F仿真