有限元法解圓柱繞流問題分析

1、有限元法求解無限流體中的圓柱繞流問題2016年01月12號一.問題描述考慮位于兩塊無限長平板間的圓柱體的平面繞流問題，幾何尺寸如下圖所示，來流為?= 1, ?= 0。由于流場具有上下左右的對稱性，只考慮左上角四分之一的計算區(qū)域abcde ,把它作為有限元的求解區(qū)域Q。要求求解出整個區(qū)域中的流函數(shù)、？?？ ？?以及壓強值。在足夠遠(yuǎn)前方選與來流垂直的控制面ae,cd是沿y軸，亦即一流動對稱軸，be是物面，ab亦是流動對稱軸，所要考慮的流動區(qū)域即由線abcdea所圍成的區(qū)域，在這一區(qū)域中有：f 、 t ?1. 邊界ab為流線，取W =0,= 0;2.邊界be也為流線，同樣取®=0,?3.

2、邊界 cd，切向速度？?=0, ?= 0,取?= 0；? 24. 邊界de為流線,滿足？?= ?+扁?？彳?= 2?于是在ed上，=2,= 0；?5. 進口邊界ae上，"=?+ £?= ?（本文中采取此條件）也可以提自然邊界條件? ?0,齊?= ?=1我們以流函數(shù)”作為未知函數(shù)來解此問題，流函數(shù)所滿足的微分方程如下:?2?= 0?r = ?（本質(zhì)邊界條件）（i）?莎? r = -?自然邊界條件）此處r指ab, bc,de和ae四段邊界，而r就是就是cd段邊界，且切向速度？?= 0,r 1和r 2合 起來是整個邊界，并且此二者不重合。下面，按有限元方法的一般步驟來計算此問題。

3、三有限元法解圓柱繞流問題1 .建立有限元積分表達(dá)式根據(jù)求解問題的基本控制方程，應(yīng)用變分法或加權(quán)余量法將求解的微分方程定解問題化 為等價的積分表達(dá)式，作為有限元法求解問題的出發(fā)方程式。對于方程（1），它是一橢圓型方程，具有正定性，可以用變分法，這里直接給出泛函1J（W） =£ / ? 2 ? e d 0+ / ?d r 0 （2）令其變分3 J=0,可以得到/ ? e ? 3 2 do / ?3 2 d=r 0（3）?r2自然邊界條件已經(jīng)包含在變分表達(dá)式中（其名稱的由來），而本質(zhì)邊界條件必須強制e滿足（因此稱其為本質(zhì)邊界條件，也稱為強制邊界條件）。如果根據(jù)原微分方程中無法給出泛函J，則

4、可以用Galerkin 加權(quán)余量方法得到積分方程，這相當(dāng)于將原來的微分方程寫為如下變分形式：/ 2 3 2 dO 0（4）?這里的3e是函數(shù)e的改變量，是一種“虛位移”，在本質(zhì)邊界條件r上為零。因此，上式做分部積分后，邊界積分僅剩下r的部分。具體為/ ? 2 ? 3 2 do / ?3e d=r 0（5）?r2即（3）式?？梢?，如果2滿足原來的微分方程和邊界條件，那么，必然有2滿足（4）式，進而滿足 式。注意，在（5）式中，包含的邊界r 2上的邊界條件信息，對邊界r 1的部分，僅知 道它是給定了函數(shù)值的邊界，卻不知道邊界上的值是多少，為了確定這些值，還需要額外的處理方法。正是因為r 2上的邊界

5、信息可以包含在積分表達(dá)式中，這種邊界條件也稱為自然 邊界條件。2. 區(qū)域剖分根據(jù)物理問題的特點以及區(qū)域的形狀，把計算區(qū)域分成許多幾何形狀規(guī)則但大小可以不同的單元，確定單元節(jié)點的數(shù)目和位置，建立表示網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)。采用的單元形狀和節(jié)點的分布，以及插值函數(shù)的選取還應(yīng)考慮到計算精度和可微性的要求。這里通過ANSYSCEMCFD14.0建模并劃分網(wǎng)格。 具體而言，網(wǎng)格將求解區(qū)域分為個 281節(jié)點和565個單元，所有單元均 為三角形單元，如圖2所示實際上由于matlab計算編程是不知如何直接讀取網(wǎng)格數(shù)據(jù)，就只 選取了 180個單元與103個節(jié)點進行計算。圖2 :左上角四分之一區(qū)域的流場及其有限單元劃分流

6、場網(wǎng)格3. 單元分析單元分析的目的是建立有限元方程。把有限元積分表達(dá)式(3)寫為各個單元求和的形式?E / ? 2 ?S d=T0 (6)?=1 ?e)r "這里？?(e)表示單元e, ?是單元總數(shù)，如果僅在一個單元上考慮上式，形式上有/? ® ? d£/?3® dr 0 ( 7)?(e)r2 nie)其中r(e)表示單元e的邊界，上式實質(zhì)上并不是一個等式，只具有形式上的意義，當(dāng)對所有的單元求和以后，才是等式。如果把線積分中的r2n r (e)換為r ),則得到的是等式，但在對所有單元求和時， 內(nèi)部邊界的線積分剛好抵消， 因此(7)也可以理解為不計內(nèi)部邊

7、界貢 獻的(3)式在單元上的表達(dá)式。流函數(shù)®在單元e內(nèi)可用如下函數(shù)近似：® = ?加(8)這里？?( i = 1 , 2, 3)為節(jié)點流函數(shù)值，？?為節(jié)點上的插值函數(shù)，上式中重復(fù)下標(biāo)表示約 定求和。將(8)代入，不難得到/ ( j + 與？?3 ? & - /? 3 ?dr ( 9):?e, ? ?"'丿巧仃 了由于3 ?的任意性，所以，對于j=1,2,3都有?$? ?$?.(而需？麗?j?區(qū)-1十？?刀r (10) 此即單元方程，通常可以簡寫為?蕩亍 /? ? Q ?=? /?d r?e)i2 n)采用三節(jié)點三角形單元時，單元的插值基函數(shù)為?x,

8、y) = ?；? +?+ ?易?(11)如果單元e三個點坐標(biāo)為(？?,？?)，i=1,2,3,則?«?；?,?)= W12)即插值基函數(shù)？3在？?點取1,在?兩點為零。由此不難解出abc。注意到求?時對？?取了梯度，故？?的取值并不影響最后的計算。對某一固定的單元6,將(11)式代入(b2 + ?= ?)? +? +?10)中，可以得到: ?+ ?b2 + c2? + ?+ ? + ?2?1b2+c：3此即采用線性單元時的單元方程系數(shù)矩陣。其中?為三角形(積分區(qū)域)的面積，be的值可由(12)求得，現(xiàn)在列舉如下:h 2A(e)(y2y3)b2b312A(e)(% y2)(X3X2)C

9、2(X1 X3)C3(X2X1)(13)?)= 1(?-求解單元系數(shù)矩陣時， 一般同時進行總體合成， 每形成一個單元方程, 體方程中。出于順序和邏輯上的考慮，下一步再詳細(xì)說明總體合成的方法。?)(? - ?) - (? - ?)(? - ?)便把它累加到總對于邊界積分項，我們假設(shè)三角形單元e中序號為,的節(jié)點在邊界2上，為自然邊界，其長度為I。首先，注意到插值函數(shù)在邊上是零。所以，可以得到如下結(jié)論:圖3:自然邊界條件的處理右端項：f佝 0。為了計算f(e)和f，以 點為原點，沿直線建立局部坐標(biāo)系，在此坐標(biāo)系中,插值函數(shù)N和N如上圖所示，可寫成線性插值函數(shù)如下：N - N -II假設(shè)切向速度Vs在

10、兩節(jié)點處的值分別為 Vs和Vs，并且沿邊界是線性分布的，可以表示為VsVsa('"sVsa ) °于是可以得到If(e)VsN dI(Vsa(V；Vsa)(1)d(2VsaVs )00II6If(e)VsN dI(Vsa('/sVsa)d- (Vsa2Vs )。00I I6對于前面討論的圓柱繞流問題，由于Vsa Vs0，所以，根據(jù)線性解的性質(zhì)，必有?(? = ?( ? = ?(? = 0a傳y無需考慮f的影響，使程序得到了不少的簡化。4. 總體合成總體合成的過程就是把已經(jīng)形成好的單元方程按一定順序迭加起來，形成總體有限元方程。具體做法是根據(jù)單元內(nèi)節(jié)點的總體順

11、序號，把單元方程進行延拓， 未知量包含所有節(jié)點上的函數(shù)值，與此單元無關(guān)的位置以零填充，把所有的單元方程都進行延拓以后，進行系數(shù)矩陣的累加，便得到總體方程。 理論上說，這一過程也可以通過引入一個 Boole型矩陣來實 現(xiàn)，定義單元e的boole矩陣B3 NB3Np1,如果單元e的地i個點總體順序號是j0.其他情況,i=1,2,3; j=1,2,3NP.矩陣B其實就是單元節(jié)點序號表的又一表達(dá)形式。單元e的系數(shù)矩陣A(e)以及右端項f佝沿拓后就是：A(e)BTA(e)B,f(e) Bt f(e)進而總體合成的過程可以表示為NeNeAA(e)， f f(e)。e 1e 1但是這種方法比較麻煩，要重新定

12、義新的矩陣E,而且還要涉及計算矩陣轉(zhuǎn)置和矩陣相乘等運算，一方面計算量較大，并且浪費空間，另一方面人為地增加了程序的復(fù)雜性，降低了程序的可讀性。因此，這種方法一般只用作理論分析。實際的計算中，每當(dāng)計算出一個單元系數(shù)矩陣？務(wù)？假設(shè)單元e三個節(jié)點編號分別為i,j,k ，那么將？?中的1*1項放入大矩陣(借鑒結(jié)構(gòu)力學(xué)的概念，不妨稱其為總體“剛度”矩陣，下同)A的i*i項中，將1*2,1*3分別放入總體剛度矩陣 A的i*j,i*k項中。同理，2*1,2*2 ,2*3,3*1,3*2,3*3項分別對應(yīng)總體剛度矩陣A的j*i,j*j,j*k,k*i,k*j,k*k項中。采用此方法并未多占用計算機內(nèi)存，運算量

13、也并不大(總共進行9*e次加法運算，不進行乘法運算)。本題的算法中選用的即此方法。(5)邊界條件處理這里的邊界條件是指本質(zhì)邊界條件，自然邊界條件已經(jīng)包含在積分表達(dá)式中了。具體做法是將已知的值代入到方程組中，并把已知的值移到方程組的右段，形成右端項。(6 )求解有限元方程組并計算相關(guān)物理量有限元方程組的求解是一個代數(shù)問題，應(yīng)針對具體的問題采用合適的方法求解。對于對角占優(yōu)的代數(shù)方程組，可以采用迭代法求解，規(guī)模不大的可以用Gauss消元法一類的直接法 求解，三對角方程則可以用追趕法。求出所有的待求量后，便得到了近似函數(shù)的表達(dá)式，并可以計算出相關(guān)的物理量。對計算結(jié)果進行綜合的分析，以期得到原問題的正確

14、的物理解答。對于每個單元，速度可以根據(jù)? 2 ?v?= ?= ?14)? ? 2?v? = -= -= -? ?= -?15)? ? ?來計算，節(jié)點上的速度值可取這個節(jié)點相鄰單元的速度值的平均。節(jié)點上的壓力值可以有伯努利方程計算。假設(shè)求解區(qū)域位于同一水平面內(nèi)，介質(zhì)密度p= 1，來流壓力p=0,那么p =2(1 - ?>,?= 1,210。如此便可得到節(jié)點處的速度和壓力分布。四具體算法實現(xiàn)本文具體算法是通過MATLA實現(xiàn)的，matlab程序文件及算法說明見附件。五結(jié)果分析運行附件中MATLAB程序，可以得到圖4和圖6兩張圖。圖4顯示1/4區(qū)域的流場分布, 圖中的點以不同顏色表示各個結(jié)點處的

15、流函數(shù)，圖中的曲線為流函數(shù)的等勢線，即是流線， 由于對稱性，整個區(qū)域的流場分布可明顯看出，故不再畫出。而圓柱繞流問題的解析解為一 1 、y(i 2)，便于與計算結(jié)果對比。從圖 4可以發(fā)現(xiàn)，有限元法數(shù)值法得到的x y流線與解析法得到的流線在形態(tài)上基本一致。-1也5圖4:有限元計算的1/4區(qū)域流場分布圖5:解析法的整個區(qū)域流場分布具體到1/4區(qū)域的右邊界上的結(jié)點，將有限元法得到的流函數(shù)數(shù)值解與解析解相對比，得到圖7中的曲線?？梢?，有限元法與解析法所得曲線在趨勢上基本一致，但在數(shù)值上最大誤差為11%。OS右邊界丫坐標(biāo)圖6: 1/4區(qū)域右邊界上流函數(shù)的有限元數(shù)值解與解析解對比附件：1. NATLA主程

16、序：youxia nyua n.mENA=41 32 37 0.03979830141 44 32 0.04752990318 101 21 0.12015681718 103 101 0.094963772103 22 102 0.043498437103 18 22 0.09424331830 36 28 0.03805466230 35 36 0.04228752129 2 12 0.03807450229 6 2 0.04874290430 7 17 0.0502916230 8 7 0.03714096419 88 76 0.05799404119 20 88 0.099591586

17、101 102 99 0.046382644101 103 102 0.04421889144 94 5 0.04713263944 41 94 0.05596740635 52 38 0.03608551335 16 52 0.03457155729 44 6 0.04428086329 32 44 0.04144814597 20 21 0.1061091921 3 93 0.0782126134 5 94 0.054017838102 100 99 0.034088039101 98 21 0.05035921622 100 102 0.04845156599 98 101 0.0451

18、8646798 97 21 0.07031062923 96 22 0.09702207296 100 22 0.06225366192 95 100 0.032213795 99 100 0.03169101191 98 99 0.03949582990 97 98 0.03840987397 89 20 0.0628527284 87 3 0.05049789793 24 1 0.09490433284 96 23 0.06040645796 92 100 0.03647283380 95 92 0.03202974895 91 99 0.02967402491 90 98 0.04186

19、491290 89 97 0.04550708389 88 20 0.0612477183 13 19 0.09334266174 94 41 0.0515515694 87 4 0.03920888382 93 3 0.05827706681 24 93 0.06379886586 85 23 0.04783281785 84 23 0.05083739684 92 96 0.03892043492 71 80 0.04137374180 91 95 0.03355940879 90 91 0.04216039578 89 90 0.04489076677 88 89 0.041466167

20、83 14 13 0.06058621414 75 15 0.05110693874 87 94 0.04157820687 82 3 0.05719513482 81 93 0.04553313924 86 23 0.08455521381 86 24 0.05624168273 85 86 0.03017268572 84 85 0.03281428584 71 92 0.03726385680 79 91 0.04726267379 78 90 0.04368349178 77 89 0.04385861877 76 88 0.03338254669 83 19 0.0531674818

21、3 75 14 0.03030115174 82 87 0.03956627568 81 82 0.04530365881 73 86 0.03902865673 72 85 0.0328008272 71 84 0.03639103670 79 80 0.05040996866 78 79 0.04320708465 77 78 0.04078362664 76 77 0.0363831876 63 19 0.06045875569 75 83 0.02597752475 58 15 0.05524126557 74 41 0.05361664474 68 82 0.04684733668

22、73 81 0.04179355462 72 73 0.03607545161 71 72 0.03593951171 67 80 0.03847238267 70 80 0.04018557959 70 56 0.04561569470 66 79 0.0431204266 65 78 0.04122181465 64 77 0.03881581264 63 76 0.03550383263 69 19 0.05773670369 58 75 0.03240444257 68 74 0.0426945868 62 73 0.03904173962 61 72 0.03774916561 67

23、 71 0.03324494467 56 70 0.03888968859 66 70 0.04081800355 65 66 0.04058518354 64 65 0.0400881353 63 64 0.04291209963 58 69 0.03605925557 62 68 0.04159156551 61 62 0.03968597661 56 67 0.03763390956 60 59 0.03628834450 60 47 0.03186017949 59 60 0.03152143459 55 66 0.03941936355 54 65 0.0398997954 53 6

24、4 0.04090381353 58 63 0.04640429958 48 15 0.06617007443 57 41 0.0476315157 51 62 0.0419565551 56 61 0.04663644756 47 60 0.03559279350 49 60 0.03179644949 55 59 0.03922927246 54 55 0.04097428845 53 54 0.0418020753 48 58 0.04219602648 52 15 0.05630523452 16 15 0.04715458543 51 57 0.04109030151 47 56 0

25、.04582129333 50 40 0.04839875542 49 50 0.04051366349 46 55 0.04020256246 45 54 0.0420780145 48 53 0.044871948 38 52 0.03726359343 47 51 0.04153873547 40 50 0.03986560233 42 50 0.05032501442 46 49 0.04102327139 45 46 0.0424682545 38 48 0.04395136544 5 6 0.05104737443 40 47 0.04178759534 42 33 0.05174

26、642742 39 46 0.04285326339 38 45 0.04174385641 37 43 0.03860959937 40 43 0.0348392340 31 33 0.05057402734 39 42 0.0438233236 38 39 0.04221094137 31 40 0.03318516834 36 39 0.04169815236 35 38 0.03949412232 31 37 0.03469899835 17 16 0.04919161534 28 36 0.04062912533 27 34 0.04955841625 31 32 0.0395515

27、5931 26 33 0.0481510830 17 35 0.04125070927 28 34 0.03984650826 27 33 0.0461411229 25 32 0.03797944325 26 31 0.03858078728 8 30 0.03611670227 9 28 0.03348170426 10 27 0.03275753512 25 29 0.03542773625 11 26 0.033372839 8 28 0.03190471110 9 27 0.03044574711 10 26 0.03028529712 11 25 0.031849171; %單元與

28、結(jié)點號、面積對應(yīng)關(guān)系矩陣NCORD=-3 0-1 0-2.6 0-2.2 0-1.8 0-1.4 00 1-0.258824103 0.965924471-0.499995466 0.866028022-0.707106781 0.707106781-0.866028022 0.499995465-0.965924471 0.2588241020 30 2.60 2.20 1.80 1.4-3 3-0.75 3-1.5 3-2.25 3-3 2.25-3 1.5-3 0.75-1.134521668 0.489438169-0.970895176 0.74446511-0.757320467

29、0.962580378-0.50550874 1.128325608-1.262125132 0.243714522-0.251458098 1.253891963-1.277106036 0.738778719-1.44751226 0.48199085-1.088167651 1.056783565-0.77523631 1.267266548-0.245958075 1.585179826 -0.503847145 1.428730143 -1.548704503 0.736752996 -0.489238864 1.743880003 -0.764674661 1.580844378

30、-1.450757531 0.981852867 -1.796720442 0.5745713 -1.049823117 1.413295395 -1.765123908 0.906580515 -1.618989726 0.255236869 -0.749693358 1.892823424 -1.029977756 1.72552432 -1.649144006 1.200203744 -0.449360041 2.05857237 -1.303039926 1.563704873 -1.34753853 1.270144924 -1.954656471 1.143055159 -0.23

31、5772926 1.875193029 -0.74066004 2.19662276 -1.020771946 2.031271425 -1.293847428 1.863609647 -1.825385481 1.467199841 -2.058484057 0.839001806 -0.455446019 2.351164971 -1.560699123 1.692649341 -1.53886561 1.437047451 -2.15060183 1.373255674 -2.25034016 1.085358385 -0.748247671 2.517911273 -1.0109146

32、94 2.329143574 -1.285646568 2.160870135 -1.563524602 1.986279407 -2.066343984 1.629034175 -2.354650626 0.785729864 -0.509120672 2.628036907 -1.839751973 1.799640354 -2.334608941 1.60379593 -2.427838139 1.329969892 -2.52208128 1.053358599 -2.160447272 0.532346877 -0.255534689 2.527394284 -0.999394278

33、 2.607756191 -1.266887377 2.454932477 -1.554892284 2.28839307-1.84609102 2.107566848-2.169360532 1.906166185-2.645492387 0.7513903-2.474781204 0.460011348-0.302931072 2.751086238-2.610333677 1.574636182-2.674743671 1.301371334-2.774519433 1.068178921-2.282601969 0.252489483-1.214562096 2.734422438-1

34、.535695128 2.604062411-1.831509008 2.416647143-2.112951495 2.23371382-2.468277014 1.859957148-2.746921781 0.391063066-1.977781367 0.27008919-2.341961855 2.093789679-2.741274474 1.859597922-1.842997677 2.717042154-2.092137247 2.544745194-2.359663521 2.342124107-2.599556714 2.126987134-2.360128075 2.6

35、79581821-2.634609866 2.379744818-2.748684485 2.57733166; %結(jié)點與坐標(biāo)對應(yīng)關(guān)系矩陣Enum,temp=size(ENA); % 返回單元總數(shù)Nnum,temp=size(NCORD); % 返回結(jié)點總數(shù)BNODE=1 3 4 5 6 2 12 11 10 9 8 7 17 16 15 14 13 19 20 21 18 22 23 24;%邊界結(jié)點編號temp,Bnum=size(BNODE); % 返回邊界結(jié)點數(shù)%邊界已知流函數(shù)的結(jié)點號及其值BKN=1 3 4 5 6 2 12 11 10 9 8 7 13 19 20 21 18 22

36、 23 24;0 0 0 0 0 0 0 0 0 0 0 0 3 3 3 3 3 2.25 1.5 0.75;temp,Bknum=size(BKN); % 返回邊界已知流函數(shù)的結(jié)點數(shù)UKN=14,15,16,17,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,

37、89,90,91,92,93,94,95,96,97,98,99,100,101,102,103;temp,Uknum=size(UKN); % 未知流函數(shù)的結(jié)點數(shù)for k=1:Enum % 單元號x1=NCORD(ENA(k,1),1);y1=NCORD(ENA(k,1),2);x2=NCORD(ENA(k,2),1);y2=NCORD(ENA(k,2),2);x3=NCORD(ENA(k,3),1);y3=NCORD(ENA(k,3),2);a(1)=x2*y3-x3*y2;b(1)=y2-y3;c(1)=x3-x2;a(2)=x3*y1-x1*y3;b(2)=y3-y1;c(2)=x1

38、-x3;a(3)=x1*y2-x2*y1;b(3)=y1-y2;c(3)=x2-x1;for n=1:3for m=1:3ANM(k,n,m)=1/(4*ENA(k,4)*(b(n)*b(m)+c(n)*c(m);% 各個單元的 Anm 信息 end end endfor i=1:Nnum % 求系數(shù)矩陣 Afor j=1:Nnumtemp=0;for e=1:Enumfor n=1:3for m=1:3if ENA(e,n)=i && ENA(e,m)=jtemp=temp+ANM(e,n,m);endendendendA(i,j)=temp;endendfor i=1:Uknum % 掃描未知流函數(shù)的結(jié)點F(i)=0;for j=1:UknumANEW(i,j)=A(UKN(i),UKN(j); %新的系數(shù)矩陣endfor k=1:Bknum % 掃描已知流函數(shù)的結(jié)點F(i)=F(i)+A(UKN(i),BKN(1,k)*BKN(2,k);endF(i)=-F(i); %移項到等號右側(cè)成新的乘積系數(shù)向量end%FAIUN=ANEWF' %解向量，未知結(jié)點的流函數(shù)值Q,R=lu(ANEW);FAIUN=R(QF'); %將所有結(jié)點的流函數(shù)值寫入 FAIN 向量 for i=1:UknumFAIN(UKN(i)=FAIU