華科微納系統(tǒng)仿真大作業(yè)

1、1、用有限差分法和有線元方法把以下問題變成數(shù)值方程，并說明兩種方法的異 同：Ñ2j(x, y) = 0邊界條件：j(x,0) = j(x,1) = 0j(0, y) = j(1, y) = 0解：（1）有限差分法 有限差分法的基本思想是將問題求解域劃分為均勻的差分網(wǎng)格，用有限個網(wǎng)格節(jié)點代替連續(xù)的求解域?；?Taylor 級數(shù)展開等方法，把描述問題的微分方程中的 微分用網(wǎng)格節(jié)點上的函數(shù)值的差分來代替，從而將微分方程轉(zhuǎn)化為以網(wǎng)格節(jié)點上函 數(shù)值為未知數(shù)的代數(shù)方程組。該方法是一種直接將微分問題變?yōu)榇鷶?shù)問題的近似數(shù) 值解法，數(shù)學(xué)概念直觀，表達簡單，是發(fā)展較早且比較成熟的數(shù)值方法，Ñ

2、2j(x, y) = 0該方程為拉普拉斯方程 Ñ2u = u+ uyyxx將運算符離散化，則 f ¢(x) 的近似值公式為f ¢(x) = f ( x + h) - 2 f ( x) + f ( x - h) + O(h2 )h2（1-1）由于定義域為 0 £ x £ 1,0 £ y £ 1, 將定義域劃分為均勻網(wǎng)域，x，y 方向的步長都記為 h，則Ñ2u = u( x + h, y) + u( x - h, y) + u( x, y + h) + u( x, y - h) - 4u( x, y) + O(h2 )

3、h2將 O(h2 ) 這個高階無窮小舍去，則（1-2）Ñ2u = u( x + h, y) + u( x - h, y) + u( x, y + h) + u( x, y - h) - 4u( x, y) = 0h2（1-3）其中在所 有的內(nèi)部網(wǎng)格點 (x, y) = (xi , y j ), i = 2,L, n -1; j = 2,L, m -1, 精度為O(h2 ) ，用 ui, j近似表示 u(xi , y j) ，則式（1-3）可以表示為u+ u+ u+ u- 4uÑ2ui, j» i +1, ji -1, ji , j +1i , j -1i , j

4、= 0h2（1-4）這就是拉普拉斯方程的五點差分公式，將 h2 消去可得ui+1, j + ui-1, j + ui, j +1 + ui, j +1 + ui, j -1 - 4ui, j = 0（1-5）根據(jù)邊界條件j(x,0) = j(x,1) = j(0, y) = j(1, y) = 0,將定義域劃分為 5×5 網(wǎng)格，如下圖所示：最終得到由 9 個方程組成的 AP=B 線性方程組，表示為：- 4 p1 + p2p1 - 4 p2 + p3p2 - 4 p3+ p4+ p5= 0= 0+ p6= 0p1- 4 p4 + p5p2+ p4 - 4 p5 + p6+ p7= 0+

5、 p8= 0p3+ p5 - 4 p6p4- 4 p7 + p8+ p9= 0= 0p5+ p7 - 4 p8 + p9= 0p6+ p8 - 4 p9 = 0通過高斯消去法可以得到解向量P = p1 , p2 ,L, p9 （2）有限元方法 使用有限元發(fā)的計算流程為：A、求解區(qū)域離散化；B、 構(gòu)造插值函數(shù)形成分段光滑的坐標(biāo)函數(shù)系；2C、 用 Ritz 方法求解微分方程。 對j(x, y) 構(gòu)造函數(shù)1 éæ ¶ jö2æ ¶ jö ùPj(x, y) = òòêçç

6、÷÷2 êëè ¶x ø+ çç÷÷è ¶y øúdxdyúû首先將整個區(qū)域離散為三角形的子區(qū)域如圖 1，三角形微小子區(qū)域中的值由三 角形節(jié)點值的插值結(jié)果表示，即圖 1 求解域離散為 n 三角形j(x, y) = Niji + N jjj + Nkjk其中 Ni , N j , Nk 為三角形的節(jié)點插值函數(shù)，ji ,jj ,jk 為函數(shù)j(x, y) 在節(jié)點 i, j, k處的函數(shù)值?？芍恳粋€三角形子區(qū)域中泛函可由節(jié)點插值函

7、數(shù)和節(jié)點函數(shù)值表示，那么對 于整個求解區(qū)域的泛函表達式為,Pj(x, y)= P(j1,j2,×××,jn )由變分原理可知，¶P = ¶P= ××× = ¶P = 0¶j1¶j2¶jn計算這 n 個方程即可得到整個求解區(qū)域的值。（3）有限差分法與有限元法的異同相同的地方：兩種方法都使用了離散化的思想，對求解域進行離散。并且都是 將求解域劃分成有限個網(wǎng)格進行近似求解。不同的地方：有限差分法是利用級數(shù)的概念將連續(xù)函數(shù)離散化，正如高等數(shù)學(xué) 上所學(xué)的連續(xù)函數(shù)用泰勒級數(shù)表達一樣，網(wǎng)格

8、上的結(jié)點就是級數(shù)中的一個取值點， 這樣以級數(shù)和的形式求得最終的解，這個解是近似解，其余項就是誤差。有限元法 是利用插值原理對求域進行近似求解，將求解域劃分網(wǎng)格，每個網(wǎng)格看作一個單元 進行求解，這樣可以得到若干有限個單元的解，這些解的集和構(gòu)成整體函數(shù)的解。 就是說每個單元一個解，這些解分布在整個求解域上，構(gòu)成不同區(qū)域解的變化，如 力的變化，溫度的變化，這樣就可以宏觀上看到在不同點上不同的值了。其中具體的表現(xiàn)在以下幾個方面：1、有限元單元的形狀沒有限制，可以做到單元形狀與邊界重合，處理無規(guī)邊 界的問題；2、有限元統(tǒng)一對待區(qū)域內(nèi)的節(jié)點與邊界上的節(jié)點，因此節(jié)點的計算精度總體 上協(xié)調(diào)；而有限差分必須分別

9、處理微分方程與定解條件；3、有限元法是各種單元可以混合使用，所以寫不出方程也可以求解；差分法 劃分的網(wǎng)格是規(guī)則的，對方程進行離散化，就是用很多個差分代替微分。2、運 行 COMSOL 軟 件 ， 打 開 ModelLibrary- à Thermal Structure Interaction-àMicrorobot 3D 的示例，詳細看懂說明書該問題的分析，用自己熟 悉的 3D 工具進行建模，用 COMSOL 進行仿真，談?wù)?COMSOL 多場耦合的體會。 解：（1）仿真過程：微型機器人的腿部使用的是由聚酰亞胺構(gòu)成的 V 型凹槽結(jié)構(gòu)。由于聚酰亞胺 具有很高的熱膨脹系數(shù)，當(dāng)

10、其受到加熱時，腿部的變形程度也會相對較低。Microrobot leg 模型用到了熱轉(zhuǎn)移模塊和結(jié)構(gòu)力學(xué)模塊，需要對它進行瞬態(tài)熱 轉(zhuǎn)移分析和準(zhǔn)靜態(tài)熱變形分析。首先在 SolidWorks 中進行建模。圖 1 SolidWorks 建模效果圖 上圖為腿部組裝完成的效果圖。在實際建模過程中，沒有考慮由 Al 和 Si2O 構(gòu)成的組合層。下圖 2 為細節(jié)部分。圖 2 腿部細節(jié)然后將 SolidWorks 模型導(dǎo)入到 Comsol 中，得到如下的模型：圖 3 導(dǎo)入到 Comsol 后的效果 實驗過程中，邊界的冷卻效應(yīng)對仿真結(jié)果有著重大的影響。根據(jù)經(jīng)驗，很大一部分的熱量是從腿部本身，也就是硅結(jié)構(gòu)上散走；而

11、在腿部的末端也會產(chǎn)生少量的 散熱現(xiàn)象。熱量直接在腿部進行傳導(dǎo)。由于各種材料的熱膨脹系數(shù)的差異，導(dǎo)致腿 部的彎曲變形。下圖 4 為仿真過程中，在 20ms 時的彎曲程度和邊界上的溫度分布。從中可以得知，仿真過程中溫度變化量的最大值為 304.7K、最小值為 274.7K。圖 4 在 20ms 時的溫度分布示意圖下圖 5 為在腿部末端，時間與其變形量的關(guān)系曲線。其中，最大的位移變形量 接近 3.5X10-5m，時間點為 0.01s。可以得知，當(dāng)加熱溫度上升時，變形量也會增加； 而當(dāng)加熱停止一段時間后，又開始向初始狀態(tài)恢復(fù)。（2）體會：圖 5 腿部末端處的變形量曲線多場耦合能最大程度地體現(xiàn)分析對象的

12、真實工況，在真實環(huán)境下，研究對象往 往同時受多個場的作用。體現(xiàn)多場耦合也有很多途徑和方法，大體來說可分為多場 直接耦合和多個場順序耦合。一般情況下，這兩種耦合都采用有限元方法來計算。 COMSOL 的多場耦合是基于加權(quán)殘值法的，且是順序耦合，因為 comsol 的各個模 塊的單元都是獨立的，在計算各個場變量的時候是由順序的。而直接耦合，由于每個單元同時包括各個場的自由度，所以它的即時和同時的，一次計算的各場變量是 沒有順序的，這也是直接耦合和順序耦合的區(qū)別。comsol 的多場耦合，雖然也是有順序性，但它是通過耦合變量來求解各個場 的反演方程和各個積分方程，耦合變量會在設(shè)定的模型中與各個獨立場

13、變量同時求 解，comsol 的多場耦合在求解一個場的時候,耦合變量會調(diào)用下一場去反算,通過引 入耦合變量耦合變量來體現(xiàn)雙向耦合。從原理上看，直接耦合最接近實際情況，因 為采用耦合單元，它的每步計算同時地，無中間條件地體現(xiàn)了多場作用和影響。個 人認(rèn)為 comsol 的多場耦合結(jié)果的可靠和合理程度，應(yīng)該居于直接耦合和順序耦合 之間。comsol 采用的耦合方法是有很大通用性和開放性的，只要有場的控制方程， 就能計算場的耦合。因而其使用范圍較廣。3、微型層流混合器（1）在 COMSOL MODEL 中 MEMS-Microfluidics-lamella mixer;（2）用熟悉的 3D 工具建模

14、，用 COMSOL 仿真；（3）設(shè)想一種更有效的微型流體混合的結(jié)構(gòu)；（4）試試用其他的軟件進行仿真并比較說明 解：（2）仿真過程： 在宏觀領(lǐng)域，流體的混合一般是通過機械激勵或者湍流模式實現(xiàn)的；而在微觀領(lǐng)域，這些方法不再適用。本模型展示的是微觀領(lǐng)域的薄片層流混合器，采用的是 分層組件來提高混合效率。此模型分析的是流體的穩(wěn)態(tài)和對流、擴散情況。流體通過微型彎曲通道進入混 合腔內(nèi)。首先在 SolidWorks 中進行建模。圖 6 SolidWorks 建模效果圖然后將 SolidWorks 模型導(dǎo)入到 Comsol 中，得到如下的模型：圖 7 導(dǎo)入到 Comsol 后的效果下圖 8 為混合器內(nèi)流體流動

15、狀況仿真效果圖，可以看出內(nèi)部顏色的漸變效果， 這表明著層流的拋物線流動。其中流線并沒有發(fā)生旋轉(zhuǎn)等狀況，而且在流動方向上 只是發(fā)生了輕微的改變。從下圖中還可以得出，微型通道內(nèi)部的最大流體速度發(fā)生 在邊界處，值大致為 1.32mm/s；而混合腔內(nèi)，可以很清楚的看到流速在水平方向 上的變化，最大的流速值大致為 0.5 mm/s。圖 8 流體流速分布圖下圖 9 展示的是模型邊線上的流體濃度分布效果。在微型通道的進口端，流體 濃度基本可以維持在固定值，0 或者 50mol/m3。當(dāng)流體進入容腔后，開始混合。在 容腔的入口處，濃度差異最大，隨著流體深入到內(nèi)部，其差異性也慢慢降低。由于 流速的影響，在容腔內(nèi)

16、部，邊界處的混合效果要比中心處的好。但是混合的效果并 不是很理想，可以嘗試降低流速、增長容腔尺寸、或者改變結(jié)構(gòu)來改善混合效果。圖 9 流體濃度分布圖下圖 10 展示的是在混合容腔中心部位的 Z 軸方向上的濃度曲線。從中可以得 知，在微型管道附近，其變化很大；在容腔的尾端，呈現(xiàn)出扁平狀的 S 型曲線變化。 在容腔的上下兩邊呈現(xiàn)出對稱趨勢，在中心處幅值大約在 25 mol/m3。圖 10 混合容腔中心部位的 Z 軸方向上的濃度曲線（3）設(shè)想的微型流體混合器的結(jié)構(gòu)： 從仿真的結(jié)果可以知道，混合的效果并不是很理想，可以采用增大容腔長度的方式來改善混合效果。即：現(xiàn)行的方案為：0.2mmX0.08mm，可

17、以將長度改為 0.3mm，使流體混合的距離變長。 同時，現(xiàn)在仿真的方案，其結(jié)構(gòu)還是相對比較簡單的，可以對其結(jié)構(gòu)進行改善。 方案 1：可以增加微型管道的彎曲度，進而可以降低流體的流速，以提高在容腔內(nèi)的混合效果。方案 2：可以將微型管道的形狀改造成蜿蜒、曲折式的。流體在這種微型管道 中流動時，能夠使得處于層流狀態(tài)的流體產(chǎn)生“無序”的流動方式。而這種無序的 流動方式，使得兩流體在接觸時，接觸面能夠被扭曲和拉長，從而更大程度上提高 混合的效果。方案 3：可以在微型管道的底部加工一些不對稱、相互錯開的人字形結(jié)構(gòu)。這 種傾斜的、不對稱的結(jié)構(gòu)會對低雷諾數(shù)的粘性流體產(chǎn)生各項異性的流阻，使得流暢 產(chǎn)生橫向的速度

18、分量。這樣整個流場就會產(chǎn)生旋轉(zhuǎn)形的流線，當(dāng)流體混合時，可以 增強混合效果。（4）與其他仿真軟件比較：ANSYS 流體單元能進行流體動力學(xué)分析，分析類型可以為瞬態(tài)或穩(wěn)態(tài)。分析 結(jié)果可以是每個節(jié)點的壓力和通過每個單元的流率。并且可以利用后處理功能產(chǎn)生 壓力、流率和溫度分布的圖形顯示。另外，還可以使用三維表面效應(yīng)單元和熱流 管單元模擬結(jié)構(gòu)的流體繞流并包括對流換熱效應(yīng)。本文嘗試采用，ansys 進行分析。首先將模型導(dǎo)入到 ansys 中，并對其進行設(shè)置，然后網(wǎng)格化后進行求解，其分 析結(jié)果如下：圖 11 ansys 仿真結(jié)果在單純的流體分析領(lǐng)域，ansys 也能產(chǎn)生一定的結(jié)果，但是其操作過程明顯要 比 comsol 來的復(fù)雜。而且 ansys 的分析結(jié)果沒有 comsol 來的好?？傮w感覺 comsol擅長于做流體分析，而 ansys 擅長于做結(jié)構(gòu)。Comsol 更加容易上手，不過如果邊界 條件設(shè)定的有問題，可能求解會發(fā)生較大的偏差。課程感受經(jīng)過一學(xué)期的微納系統(tǒng)仿真的課程，首先對微納系統(tǒng)方面有了更深入的了解。 對仿真的過程比如前處理，求解，后處理的三個過程有了一定的認(rèn)識，拓寬了自己 的視野，明白在微觀條