有限元法求解無限流體中的圓柱繞流問題2016年01月12號問題描述考慮位于兩塊無限長平板間的圓柱體的平面繞流問題，幾何尺寸如下圖所示，來流為vx=1,vy=0。由于流場具有上下左右的對稱性，只考慮左上角四分之一的計算區(qū)域abcde,把它作為有限元的求解區(qū)域Q。要求求解出整個區(qū)域中的流函數(shù)、vx、vy以及壓強值。數(shù)學(xué)建模在足夠遠前方選與來流垂直的控制面ae,cd是沿y軸，亦即一流動對稱軸，bc是物面，ab亦是流動對稱軸，所要考慮的流動區(qū)域即由線abcdea所圍成的區(qū)域，在這一區(qū)域中有：邊界ab為流線，取寸=0,dedn=0;邊界bc也為流線，同樣取中=0,3^dn=0;邊界cd，切向速度vt=0,舛dn=0，取e=0；邊界de為流線，滿足^e=^a+aevxdy=02dy=2于是在ed上，寸=2，dedn=0;進口邊界ae上，WWa+ayvxdy=y（本文中采取此條件）也可以提自然邊界條件舛dn=0，d^dn=vx=1我們以流函數(shù)W作為未知函數(shù)來解此問題，流函數(shù)所滿足的微分方程如下：V2^=0 間「1=+（本質(zhì)邊界條件）d^dn|F2=-vs（自然邊界條件）（1）此處ri指ab,bc,de和2?四段邊界，而「2就是就是cd段邊界，且切向速度vs=0，r1和「2合起來是整個邊界，并且此二者不重合。下面，按有限元方法的一般步驟來計算此問題。有限元法解圓柱繞流問題1.建立有限元積分表達式根據(jù)求解問題的基本控制方程，應(yīng)用變分法或加權(quán)余量法將求解的微分方程定解問題化為等價的積分表達式，作為有限元法求解問題的出發(fā)方程式。對于方程（1），它是一橢圓型方程，具有正定性，可以用變分法，這里直接給出泛函J(寸)=12QVUVWdQ+r2vs^dr=0(2)令其變分51=0，可以得到nV^-V6^dn+F2vs6^dr=0(3)自然邊界條件已經(jīng)包含在變分表達式中(其名稱的由來)，而本質(zhì)邊界條件必須強制寸滿足(因此稱其為本質(zhì)邊界條件，也稱為強制邊界條件)。如果根據(jù)原微分方程中無法給出泛函人則可以用Galerkin加權(quán)余量方法得到積分方程，這相當(dāng)于將原來的微分方程寫為如下變分形式：Q△中舛dQ=0(4)這里的。中是函數(shù)W的改變量，是一種“虛位移"，在本質(zhì)邊界條件口上為零。因此，上式做分部積分后，邊界積分僅剩下「2的部分。具體為QVU對巾dQ+r2vs8^dr=0(5)即(3)式?？梢姡绻鸚滿足原來的微分方程和邊界條件，那么，必然有W滿足(4)式，進而滿足(5)式。注意，在(5)式中，包含的邊界板上的邊界條件信息，對邊界「1的部分，僅知道它是給定了函數(shù)值的邊界，卻不知道邊界上的值是多少，為了確定這些值，還需要額外的處理方法。正是因為「2上的邊界信息可以包含在積分表達式中，這種邊界條件也稱為自然邊界條件。區(qū)域剖分根據(jù)物理問題的特點以及區(qū)域的形狀，把計算區(qū)域分成許多幾何形狀規(guī)則但大小可以不同的單元，確定單元節(jié)點的數(shù)目和位置，建立表示網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)。采用的單元形狀和節(jié)點的分布，以及插值函數(shù)的選取還應(yīng)考慮到計算精度和可微性的要求。這里通過ANSYSICEMCFD14.0建模并劃分網(wǎng)格。具體而言，網(wǎng)格將求解區(qū)域分為個281節(jié)點和565個單元，所有單元均為三角形單元，如圖2所示實際上由于matlab計算編程是不知如何直接讀取網(wǎng)格數(shù)據(jù)，就只選取了180個單元與103個節(jié)點進行計算。圖2:左上角四分之一區(qū)域的流場及其有限單元劃分流場網(wǎng)格單元分析

單元分析的目的是建立有限元方程。把有限元積分表達式(3)寫為各個單元求和的形式e=1NeQeV^V舛dQ+r2vs^dr=0(6)這里Q(e)表示單元e，Ne是單元總數(shù)，如果僅在一個單元上考慮上式，形式上有Q(e)V^-V6^dQ+F2nFevs6^dF=0(7)其中r(e)表示單元e的邊界，上式實質(zhì)上并不是一個等式，只具有形式上的意義，當(dāng)對所有的單元求和以后，才是等式。如果把線積分中的「2個「/)換為「("則得到的是等式，但在對所有單元求和時，內(nèi)部邊界的線積分剛好抵消，因此(7)也可以理解為不計內(nèi)部邊界貢獻的(3)式在單元上的表達式。流函數(shù)W在單元e內(nèi)可用如下函數(shù)近似：中=0Ni(8)這里0(i=1,2,3)為節(jié)點流函數(shù)值，Ni為節(jié)點上的插值函數(shù)，上式中重復(fù)下標(biāo)表示約定求和。將(8)代入(7)，不難得到Qe(dNidxdNjdx+dNidydNjdy)巾i$WjdQ=-F2nFevsNj6^jdF(9)由于&電的任意性，所以，對于j=1,2,3都有Qe(dNidxdNjdx+dNidydNjdy)巾idQ=-F2nFevsNjdF(10)此即單元方程，通常可以簡寫為Aij(e)=QeVNiVNjdQ fj=T2nTevsNjdr采用三節(jié)點三角形單元時，單元的插值基函數(shù)為Nix,y=aie+biex+ciey(11)如果單元e三個點坐標(biāo)為(xie,yie),i=1,2,3，則Nixje,yje=6ij(12)即插值基函數(shù)Ni在xie點取1，在xjexke兩點為零。由此不難解出abc。注意到求Aij(e)時對Ni取了梯度，故aie的取值并不影響最后的計算。對某一固定的單元e，將(11)式代入(10)中，可以得到：Aij(e)=A(e)b12+c12b1b2+c1c2b1b3+c1c3b1b2+c1c2b22+c22b2b3+c2c3b1b3+c1c3b2b3+c2c3b32+c32(積分區(qū)域)的面積，bc的值可由此即采用線性單元時的單元方程系數(shù)矩陣。其中A(e)為三角形((積分區(qū)域)的面積，bc的值可由b=*(y-y)2 2A(b=*(y-y)2 2A(e) 3b=?一(y-y)1 2A(e) 2 3c=-^(c=-^(x-x)1 2A(e) 3 2―-—(x-x)2A(e)13土M1)(13)A(e)=12[x2-x1y3-y1-(y2-y1)(x3-x1)]求解單元系數(shù)矩陣時，一般同時進行總體合成，每形成一個單元方程，便把它累加到總體方程中。出于順序和邏輯上的考慮，下一步再詳細說明總體合成的方法。對于邊界積分項，我們假設(shè)三角形單元e中序號為a,P，的節(jié)點在邊界&上，郎為自然邊界，其長度為1。首先，注意到插值函數(shù)在a&邊上是零。所以，可以得到如下結(jié)論：

圖3：自然邊界條件的處理右端項：f（。）=0。為了計算如和號，以a點為原點，沿直線亟建立局部坐標(biāo)系&，在此坐標(biāo)系中，插vA U值函數(shù)Na和N§如上圖所示，可寫成線性插值函數(shù)如下：& &N=1-l,N廣l假設(shè)切向速度匕在兩節(jié)點處的值分別為匕a和匕§，并且沿邊界是線性分布的，可以表示為吃+（七廠吃）I。于是可以得到f（e）=-jvNd^=—\（v+（v-v）當(dāng)（1-弓）d&=-4（2v+v）a sa sas§sall 6sas§TOC\o"1-5"\h\z0 0f（e）=-fvNdg=-j（v+（v-v）&5d^=-—（v+2v）P sP sa sp sall 6sa sp0 0對于前面討論的圓柱繞流問題，由于vsa=v,p=0,所以，根據(jù)線性解的性質(zhì)，必有fa（e）=fp（e）=fy（e）=0無需考慮f的影響，使程序得到了不少的簡化。總體合成總體合成的過程就是把已經(jīng)形成好的單元方程按一定順序迭加起來，形成總體有限元方程。具體做法是根據(jù)單元內(nèi)節(jié)點的總體順序號，把單元方程進行延拓，未知量包含所有節(jié)點上的函數(shù)值，與此單元無關(guān)的位置以零填充，把所有的單元方程都進行延拓以后，進行系數(shù)矩陣的累加，便得到總體方程。理論上說，這一過程也可以通過引入一個Boole型矩陣來實現(xiàn)，定義單元e的boole矩陣B3xNp1,如果單元e的地i個點總體順序號是j.0. 其他情況，i=1，2,3；戶1,2,3…Np矩陣B其實就是單元節(jié)點序號表的又一表達形式。單元e的系數(shù)矩陣A(e)以及右端項f(e)沿拓后就是：A(e)=BtA(e)B, f(e)=Btf(e)進而總體合成的過程可以表示為A=ZA；， f*。但是這種方法比較麻煩，要重新定義新的矩陣B,而且還要涉及計算矩陣轉(zhuǎn)置和矩陣相乘等運算，一方面計算量較大，并且浪費空間，另一方面人為地增加了程序的復(fù)雜性，降低了程序的可讀性。因此，這種方法一般只用作理論分析。實際的計算中，每當(dāng)計算出一個單元系數(shù)矩陣Aij(e)，假設(shè)單元e三個節(jié)點編號分別為i,j,k，那么將Aij(e)中的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)存，運算量也并不大(總共進行9*e次加法運算，不進行乘法運算)。本題的算法中選用的即此方法。邊界條件處理這里的邊界條件是指本質(zhì)邊界條件，自然邊界條件已經(jīng)包含在積分表達式中了。具體做法是將已知的值代入到方程組中，并把已知的值移到方程組的右段，形成右端項。求解有限元方程組并計算相關(guān)物理量有限元方程組的求解是一個代數(shù)問題，應(yīng)針對具體的問題采用合適的方法求解。對于對角占優(yōu)的代數(shù)方程組，可以采用迭代法求解，規(guī)模不大的可以用Gauss消元法一類的直接法求解，三對角方程則可以用追趕法。求出所有的待求量后，便得到了近似函數(shù)的表達式，并可以計算出相關(guān)的物理量。對計算結(jié)果進行綜合的分析，以期得到原問題的正確的物理解答。對于每個單元，速度可以根據(jù)vx=dWdy=d0Nidy=0dNidy=ci0(14)vy=-dWdx=-d0Nidx=-WidNidx=-bi0(15)來計算，節(jié)點上的速度值可取這個節(jié)點相鄰單元的速度值的平均。節(jié)點上的壓力值可以有伯努利方程計算。假設(shè)求解區(qū)域位于同一水平面內(nèi)，介質(zhì)密度P=1，來流壓力P=0，那么p=121-vi2,i=1,2...10。如此便可得到節(jié)點處的速度和壓力分布。具體算法實現(xiàn)本文具體算法是通過MATLAB實現(xiàn)的，matlab程序文件及算法說明見附件。結(jié)果分析運行附件中MATLAB程序，可以得到圖4和圖6兩張圖。圖4顯示1/4區(qū)域的流場分布，圖中的點以不同顏色表示各個結(jié)點處的流函數(shù)，圖中的曲線為流函數(shù)的等勢線，即是流線，由于對稱性，整個區(qū)域的流場分布可明顯看出，故不再畫出。而圓柱繞流問題的解析解為w=y(i-一1—)，便于與計算結(jié)果對比。從圖4可以發(fā)現(xiàn)，有限元法數(shù)值法得到的X2+y2流線與解析法得到的流線在形態(tài)上基本一致。如.5」?，15如.5」?，15圖4:有限元計算的1/4區(qū)域流場分布圖5:解析法的整個區(qū)域流場分布具體到1/4區(qū)域的右邊界上的結(jié)點，將有限元法得到的流函數(shù)數(shù)值解與解析解相對比，得到圖7中的曲線?？梢?，有限元法與解析法所得曲線在趨勢上基本一致，但在數(shù)值上最大誤差為11%。右邊界V坐標(biāo)圖6:1/4區(qū)域右邊界上流函數(shù)的有限元數(shù)值解與解析解對比附件：1.NATLAB主程序：youxianyuan.mENA=[4132370.03979830144320.04752990318101210.1201568171031010.094963772103221020.04349843710318220.0942433183036280.0380546623035360.042287521292120.038074502620.0487429047170.0502916230870.03714096488760.0579940411920880.099591586101102990.0463826441031020.0442188919450.0471326394441940.05596740652380.0360855133516520.034571557294460.0442808632932440.04144814520210.10610919213930.07821261345940.054017838100990.03408803910198210.0503592161001020.04845156599981010.04518646797210.07031062996220.097022072100220.062253661951000.032213795991000.0316910119198990.0394958299097980.03840987389200.06285272848730.0504978972410.0949043328496230.06040645796921000.03647283395920.0320297489591990.02967402490980.0418649129089970.0455070838988200.0612477113190.0933426617494410.051551568740.0392088839330.05827706624930.06379886585230.0478328178584230.05083739692960.03892043471800.04137374191950.03355940890910.04216039589900.04489076688890.04146616714130.0605862141475150.0511069387487940.0415782068230.05719513481930.04553313986230.08455521386240.05624168285860.03017268584850.03281428571920.03726385679910.04726267378900.04368349177890.0438586187776880.03338254683190.05316748175140.03030115182870.03956627581820.04530365873860.03902865672850.032800827271840.03639103679800.05040996878790.04320708477780.0407836266476770.036383187663190.06045875575830.0259775247558150.0552412655774410.05361664468820.04684733673810.04179355472730.03607545171720.03593951167800.03847238270800.0401855795970560.04561569466790.0431204265780.0412218146564770.0388158126463760.03550383269190.05773670358750.0324044425768740.0426945862730.03904173961720.0377491656167710.03324494456700.0388896885966700.04081800365660.04058518364650.0400881363640.04291209958690.03605925562680.0415915655161620.0396859766156670.03763390960590.03628834460470.03186017959600.0315214345955660.03941936354650.0398997953640.0409038135358630.04640429948150.0661700744357410.0476315151620.0419565556610.04663644747600.03559279349600.0317964494955590.0392292724654550.04097428853540.041802075348580.04219602652150.0563052345216150.0471545854351570.04109030147560.04582129350400.04839875549500.04051366346550.04020256245540.042078014548530.04487194838520.03726359347510.04153873540500.03986560242500.05032501446490.0410232713945460.042468254538480.043951365560.05104737440470.04178759542330.0517464274239460.04285326338450.0417438564137430.0386095993740430.0348392331330.0505740273439420.0438233238390.04221094131400.03318516836390.0416981523635380.03949412231370.03469899817160.0491916153428360.04062912527340.04955841631320.0395515593126330.048151083017350.04125070928340.03984650827330.046141122925320.0379794432526310.0385807878300.036116702279280.0334817042610270.0327575351225290.0354277362511260.033372838280.0319047119270.03044574710260.03028529711250.031849171]；%單元與結(jié)點號、面積對應(yīng)關(guān)系矩陣NCORD=[-30-10-2.60-2.20-1.80-1.4001-0.2588241030.965924471-0.4999954660.866028022-0.7071067810.707106781-0.8660280220.499995465-0.9659244710.2588241020302.602.201.801.4-33-0.753-1.53-2.253-32.25-31.5-30.75-1.1345216680.489438169-0.9708951760.74446511-0.7573204670.962580378-0.505508741.128325608-1.2621251320.243714522-0.2514580981.253891963-1.2771060360.738778719-1.447512260.48199085-1.0881676511.056783565-0.775236311.267266548-0.2459580751.585179826-0.5038471451.428730143-1.5487045030.736752996-0.4892388641.743880003-0.7646746611.580844378-1.4507575310.981852867-1.7967204420.5745713-1.0498231171.413295395-1.7651239080.906580515-1.6189897260.255236869-0.7496933581.892823424-1.0299777561.72552432-1.6491440061.200203744-0.4493600412.05857237-1.3030399261.563704873-1.347538531.270144924-1.9546564711.143055159-0.2357729261.875193029-0.740660042.19662276-1.0207719462.031271425-1.2938474281.863609647-1.8253854811.467199841-2.0584840570.839001806-0.4554460192.351164971-1.5606991231.692649341-1.538865611.437047451-2.150601831.373255674-2.250340161.085358385-0.7482476712.517911273-1.0109146942.329143574-1.2856465682.160870135-1.5635246021.986279407-2.0663439841.629034175-2.3546506260.785729864-0.5091206722.628036907-1.8397519731.799640354-2.3346089411.60379593-2.4278381391.329969892-2.522081281.053358599-2.1604472720.532346877-0.2555346892.527394284-0.9993942782.607756191-1.2668873772.454932477-1.5548922842.28839307-1.846091022.107566848-2.1693605321.906166185-2.6454923870.7513903-2.4747812040.460011348-0.3029310722.751086238-2.6103336771.574636182-2.6747436711.301371334-2.7745194331.068178921-2.2826019690.252489483-1.2145620962.734422438-1.5356951282.604062411-1.8315090082.416647143-2.1129514952.23371382-2.4682770141.859957148-2.7469217810.391063066-1.9777813670.27008919-2.3419618552.093789679-2.7412744741.859597922-1.8429976772.717042154-2.0921372472.544745194-2.3596635212.342124107-2.5995567142.126987134-2.3601280752.679581821-2.6346098662.379744818-2.7486844852.57733166]；%結(jié)點與坐標(biāo)對應(yīng)關(guān)系矩陣[Enum,temp]=size(ENA);%返回單元總數(shù)[Nnum,temp]=size(NCORD);%返回結(jié)點總數(shù)BNODE=[134562121110987171615141319202118222324];%邊界結(jié)點編號[temp,Bnum]=size(BNODE);%返回邊界結(jié)點數(shù)%邊界已知流函數(shù)的結(jié)點號及其值BKN=[1345621211109871319202118222324;000000000000333332.251.50.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,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103];[temp,Uknum]=size(UKN);%未知流函數(shù)的結(jié)點數(shù)fork=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-x3;a(3)=x1*y2-x2*y1;b(3)=y1-y2;c(3)=x2-x1;forn=1:3form=1:3ANM(k,n,m)=1/(4*ENA(k,4))*(b(n)*b(m)+c(n)*c(m));%各個單元的Anm信息endendendfori=1:Nnum%求系數(shù)矩陣Aforj=1:Nnumtemp=0;fore=1:Enumforn=1:3form=1:3ifENA(e,n)==i&&ENA(e,m)==jtemp=temp+ANM(e,n,m);endendendendA(i,j)=temp;endendfori=1:Uknum%掃描未知流函數(shù)的結(jié)點F(i)=0;forj=1:UknumANEW(i,j)=A(UKN(i),UKN(j));%新的系數(shù)矩陣endfork=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=ANEW\F'%解向量，未知結(jié)點的流函數(shù)值[Q,R]=lu(