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

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

2、2j(x, y) = 0該方程為拉普拉斯方程 Ñ2u = u+ uyyxx將運(yùn)算符離散化，則 f ¢(x) 的近似值公式為f ¢(x) = f ( x + h) - 2 f ( x) + f ( x - h) + O(h2 )h2（1-1）由于定義域?yàn)?0 £ x £ 1,0 £ y £ 1, 將定義域劃分為均勻網(wǎng)域，x，y 方向的步長(zhǎng)都記為 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 ) 這個(gè)高階無(wú)窮小舍去，則（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)格點(diǎn) (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）這就是拉普拉斯方程的五點(diǎn)差分公式，將 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 個(gè)方程組成的 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通過(guò)高斯消去法可以得到解向量P = p1 , p2 ,L, p9 （2）有限元方法 使用有限元發(fā)的計(jì)算流程為：A、求解區(qū)域離散化；B、 構(gòu)造插值函數(shù)形成分段光滑的坐標(biāo)函數(shù)系；2C、 用 Ritz 方法求解微分方程。 對(duì)j(x, y) 構(gòu)造函數(shù)1 éæ ¶ jö2æ ¶ jö ùPj(x, y) = òòêçç

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

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

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

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

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

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

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

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

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

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

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

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

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