有限元法解圓柱繞流問題

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

2、，切向速度v=0，n=0，取=0；4.邊界de為流線,滿足 e=a+aevxdy=02dy=2于是在ed上，=2，n=0;5.進口邊界ae上，=a+ayvxdy=y（本文中采取此條件）也可以提自然邊界條件n=0，n=vx=1我們以流函數(shù)作為未知函數(shù)來解此問題，流函數(shù)所滿足的微分方程如下：2=0 |1=（本質邊界條件）n|2=-vs(自然邊界條件)（1）此處1指ab,bc,de和ae四段邊界，而2就是就是cd 段邊界，且切向速度vs= 0，1 和2 合起來是整個邊界，并且此二者不重合。下面，按有限元方法的一般步驟來計算此問題。三有限元法解圓柱繞流問題1建立有限元積分表達式根據(jù)求解問題的基本控制方

3、程，應用變分法或加權余量法將求解的微分方程定解問題化為等價的積分表達式，作為有限元法求解問題的出發(fā)方程式。對于方程（1），它是一橢圓型方程，具有正定性，可以用變分法，這里直接給出泛函J（）=12d+2vsd=0（2）令其變分J=0，可以得到d+2vsd=0(3)自然邊界條件已經(jīng)包含在變分表達式中(其名稱的由來)，而本質邊界條件必須強制 滿足(因此稱其為本質邊界條件，也稱為強制邊界條件)。如果根據(jù)原微分方程中無法給出泛函J，則可以用Galerkin 加權余量方法得到積分方程，這相當于將原來的微分方程寫為如下變分形式：d=0(4)這里的是函數(shù)的改變量，是一種“虛位移”，在本質邊界條件1上為零。因此

4、，上式做分部積分后，邊界積分僅剩下2的部分。具體為d+2vsd=0(5)即（3）式。可見，如果滿足原來的微分方程和邊界條件，那么，必然有滿足(4) 式，進而滿足(5) 式。注意，在(5) 式中，包含的邊界2 上的邊界條件信息，對邊界1 的部分，僅知道它是給定了函數(shù)值的邊界，卻不知道邊界上的值是多少，為了確定這些值，還需要額外的處理方法。正是因為2 上的邊界信息可以包含在積分表達式中，這種邊界條件也稱為自然邊界條件。2區(qū)域剖分根據(jù)物理問題的特點以及區(qū)域的形狀，把計算區(qū)域分成許多幾何形狀規(guī)則但大小可以不同的單元，確定單元節(jié)點的數(shù)目和位置，建立表示網(wǎng)格的數(shù)據(jù)結構。采用的單元形狀和節(jié)點的分布，以及插值

5、函數(shù)的選取還應考慮到計算精度和可微性的要求。這里通過ANSYS ICEM CFD 14.0建模并劃分網(wǎng)格。具體而言，網(wǎng)格將求解區(qū)域分為個281節(jié)點和565個單元，所有單元均為三角形單元，如圖2所示實際上由于matlab計算編程是不知如何直接讀取網(wǎng)格數(shù)據(jù)，就只選取了180個單元與103個節(jié)點進行計算。圖2：左上角四分之一區(qū)域的流場及其有限單元劃分流場網(wǎng)格3.單元分析單元分析的目的是建立有限元方程。把有限元積分表達式(3) 寫為各個單元求和的形式e=1Need+2vsd=0（6）這里（e）表示單元e，Ne是單元總數(shù)，如果僅在一個單元上考慮上式，形式上有（e）d+2evsd=0（7）其中（e）表示單

6、元e的邊界，上式實質上并不是一個等式，只具有形式上的意義，當對所有的單元求和以后，才是等式。如果把線積分中的2(e) 換為(e)，則得到的是等式，但在對所有單元求和時，內(nèi)部邊界的線積分剛好抵消，因此(7) 也可以理解為不計內(nèi)部邊界貢獻的(3) 式在單元上的表達式。流函數(shù)在單元e內(nèi)可用如下函數(shù)近似：=iNi (8)這里i (i = 1, 2, 3) 為節(jié)點流函數(shù)值，Ni為節(jié)點上的插值函數(shù)，上式中重復下標表示約定求和。將(8) 代入(7)，不難得到e(NixNjx+NiyNjy) ijd=-2evsNjjd（9）由于j的任意性，所以，對于j=1,2,3都有e(NixNjx+NiyNjy) id=-

7、2evsNjd（10）此即單元方程，通常可以簡寫為Aij(e)=eNiNj d fj=-2evsNjd采用三節(jié)點三角形單元時，單元的插值基函數(shù)為Nix,y=aie+biex+ciey（11）如果單元e三個點坐標為（xie,yie）,i=1,2,3,則Nixje,yje=ij（12）即插值基函數(shù)Ni在xie點取1，在xjexke兩點為零。由此不難解出abc。注意到求Aij(e)時對Ni取了梯度，故aie的取值并不影響最后的計算。對某一固定的單元e，將（11）式代入（10）中，可以得到：Aij(e)=A(e)b12+c12b1b2+c1c2b1b3+c1c3b1b2+c1c2b22+c22b2b3

8、+c2c3b1b3+c1c3b2b3+c2c3b32+c32此即采用線性單元時的單元方程系數(shù)矩陣。其中A(e)為三角形（積分區(qū)域）的面積，bc的值可由（12）求得，現(xiàn)在列舉如下： （13）A(e)=12x2-x1y3-y1-(y2-y1)(x3-x1)求解單元系數(shù)矩陣時，一般同時進行總體合成，每形成一個單元方程，便把它累加到總體方程中。出于順序和邏輯上的考慮，下一步再詳細說明總體合成的方法。對于邊界積分項，我們假設三角形單元e中序號為，的節(jié)點在邊界上，為自然邊界，其長度為l。首先，注意到插值函數(shù)在邊上是零。所以，可以得到如下結論： 圖3：自然邊界條件的處理右端項：。為了計算和，以點為原點，沿直

9、線建立局部坐標系，在此坐標系中，插值函數(shù)和如上圖所示，可寫成線性插值函數(shù)如下：假設切向速度在兩節(jié)點處的值分別為和，并且沿邊界是線性分布的，可以表示為。于是可以得到。對于前面討論的圓柱繞流問題，由于，所以，根據(jù)線性解的性質，必有f（e）=f（e）=f（e）=0無需考慮f的影響，使程序得到了不少的簡化。4.總體合成總體合成的過程就是把已經(jīng)形成好的單元方程按一定順序迭加起來，形成總體有限元方程。具體做法是根據(jù)單元內(nèi)節(jié)點的總體順序號，把單元方程進行延拓，未知量包含所有節(jié)點上的函數(shù)值，與此單元無關的位置以零填充，把所有的單元方程都進行延拓以后，進行系數(shù)矩陣的累加，便得到總體方程。理論上說，這一過程也可以

10、通過引入一個Boole 型矩陣來實現(xiàn)，定義單元e的boole矩陣,i=1,2,3; j=1,2,3Np.矩陣B其實就是單元節(jié)點序號表的又一表達形式。單元e的系數(shù)矩陣以及右端項沿拓后就是： 進而總體合成的過程可以表示為， 。但是這種方法比較麻煩，要重新定義新的矩陣，而且還要涉及計算矩陣轉置和矩陣相乘等運算，一方面計算量較大，并且浪費空間，另一方面人為地增加了程序的復雜性，降低了程序的可讀性。因此，這種方法一般只用作理論分析。實際的計算中，每當計算出一個單元系數(shù)矩陣Aij(e)，假設單元e三個節(jié)點編號分別為i,j,k，那么將Aij(e)中的1*1項放入大矩陣(借鑒結構力學的概念，不妨稱其為總體“剛

11、度”矩陣，下同)A的i*i項中，將1*2,1*3分別放入總體剛度矩陣A的i*j,i*k項中。同理，2*1,2*2，2*3,3*1,3*2,3*3項分別對應總體剛度矩陣A的j*i,j*j,j*k,k*i,k*j,k*k項中。采用此方法并未多占用計算機內(nèi)存，運算量也并不大（總共進行9*e次加法運算，不進行乘法運算）。本題的算法中選用的即此方法。（5）邊界條件處理這里的邊界條件是指本質邊界條件，自然邊界條件已經(jīng)包含在積分表達式中了。具體做法是將已知的值代入到方程組中，并把已知的值移到方程組的右段，形成右端項。（6）求解有限元方程組并計算相關物理量有限元方程組的求解是一個代數(shù)問題，應針對具體的問題采用

12、合適的方法求解。對于對角占優(yōu)的代數(shù)方程組，可以采用迭代法求解，規(guī)模不大的可以用Gauss 消元法一類的直接法求解，三對角方程則可以用追趕法。求出所有的待求量后，便得到了近似函數(shù)的表達式，并可以計算出相關的物理量。對計算結果進行綜合的分析，以期得到原問題的正確的物理解答。對于每個單元，速度可以根據(jù)vx=y=iNiy=iNiy=cii(14)vy=-x=-iNix=-iNix=-bii(15)來計算,節(jié)點上的速度值可取這個節(jié)點相鄰單元的速度值的平均。節(jié)點上的壓力值可以有伯努利方程計算。假設求解區(qū)域位于同一水平面內(nèi)，介質密度=1，來流壓力p=0，那么p=121-vi2,i=1,210。如此便可得到節(jié)

13、點處的速度和壓力分布。四具體算法實現(xiàn)本文具體算法是通過MATLAB實現(xiàn)的，matlab程序文件及算法說明見附件。五結果分析運行附件中MATLAB 程序，可以得到圖4和圖6兩張圖。圖4顯示1/4 區(qū)域的流場分布，圖中的點以不同顏色表示各個結點處的流函數(shù)，圖中的曲線為流函數(shù)的等勢線，即是流線，由于對稱性，整個區(qū)域的流場分布可明顯看出，故不再畫出。而圓柱繞流問題的解析解為，便于與計算結果對比。從圖4可以發(fā)現(xiàn)，有限元法數(shù)值法得到的流線與解析法得到的流線在形態(tài)上基本一致。圖4：有限元計算的1/4區(qū)域流場分布圖5：解析法的整個區(qū)域流場分布具體到1/4 區(qū)域的右邊界上的結點，將有限元法得到的流函數(shù)數(shù)值解與解

14、析解相對比，得到圖7 中的曲線。可見，有限元法與解析法所得曲線在趨勢上基本一致，但在數(shù)值上最大誤差為11%。圖6: 1/4區(qū)域右邊界上流函數(shù)的有限元數(shù)值解與解析解對比附件：1. NATLAB主程序：youxianyuan.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.03807450

15、229 6 2 0.04874290430 7 17 0.0502916230 8 7 0.03714096419 88 76 0.05799404119 20 88 0.099591586101 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

16、94 0.054017838102 100 99 0.034088039101 98 21 0.05035921622 100 102 0.04845156599 98 101 0.04518646798 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

17、 96 23 0.06040645796 92 100 0.03647283380 95 92 0.03202974895 91 99 0.02967402491 90 98 0.04186491290 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

18、 80 0.04137374180 91 95 0.03355940879 90 91 0.04216039578 89 90 0.04489076677 88 89 0.04146616783 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

19、91 0.04726267379 78 90 0.04368349178 77 89 0.04385861877 76 88 0.03338254669 83 19 0.05316748183 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

20、 0.06045875569 75 83 0.02597752475 58 15 0.05524126557 74 41 0.05361664474 68 82 0.04684733668 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

21、0.05773670369 58 75 0.03240444257 68 74 0.0426945868 62 73 0.03904173962 61 72 0.03774916561 67 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.

22、03628834450 60 47 0.03186017949 59 60 0.03152143459 55 66 0.03941936355 54 65 0.0398997954 53 64 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.041

23、8020753 48 58 0.04219602648 52 15 0.05630523452 16 15 0.04715458543 51 57 0.04109030151 47 56 0.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.0410232

24、7139 45 46 0.0424682545 38 48 0.04395136544 5 6 0.05104737443 40 47 0.04178759534 42 33 0.05174642742 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

25、31 37 0.03469899835 17 16 0.04919161534 28 36 0.04062912533 27 34 0.04955841625 31 32 0.03955155931 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 2

26、6 0.033372839 8 28 0.03190471110 9 27 0.03044574711 10 26 0.03028529712 11 25 0.031849171; %單元與結點號、面積對應關系矩陣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

27、 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 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.503

28、847145 1.428730143-1.548704503 0.736752996-0.489238864 1.743880003-0.764674661 1.580844378-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.058572

29、37-1.303039926 1.563704873-1.34753853 1.270144924-1.954656471 1.143055159-0.235772926 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.

30、373255674-2.25034016 1.085358385-0.748247671 2.517911273-1.010914694 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.160

31、447272 0.532346877-0.255534689 2.527394284-0.999394278 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.0681789

32、21-2.282601969 0.252489483-1.214562096 2.734422438-1.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

33、 2.342124107-2.599556714 2.126987134-2.360128075 2.679581821-2.634609866 2.379744818-2.748684485 2.57733166; %結點與坐標對應關系矩陣Enum,temp=size(ENA); % 返回單元總數(shù)Nnum,temp=size(NCORD); % 返回結點總數(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;%邊界結點編號temp,Bnum=size(BNODE); % 返回邊界結點數(shù)%邊界已知流函數(shù)的結點

34、號及其值BKN=1 3 4 5 6 2 12 11 10 9 8 7 13 19 20 21 18 22 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ù)的結點數(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

35、,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103;temp,Uknum=size(UKN); %未知流函數(shù)的結點數(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

36、-y3;c(1)=x3-x2;a(2)=x3*y1-x1*y3;b(2)=y3-y1;c(2)=x1-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 信息endendendfor 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ù)的結點F(i)=0;for j=1:UknumANEW(i,j)=A(UKN(i),UKN(j); %新的系數(shù)矩陣endfor k=1:Bknum %掃描已知流函數(shù)的結點F(i)=F(i)+A(UKN(i),BKN(1,k)*BKN(2,k);endF(i)=-F(i); %移項到等號右側成新的乘積系數(shù)向量end%FAIUN=ANEWF %解向量，未知結點的流函數(shù)值Q,R