有限元法解圓柱繞流報(bào)告

有限元法解圓柱繞流計(jì)算流體力學(xué)期末大作業(yè)工學(xué)院熊思009860802012/6/3

【摘要】 有限元法是求解計(jì)算流體力學(xué)問(wèn)題的重要方法。用有限元法求解具體問(wèn)題時(shí)，首先需要將求解區(qū)域進(jìn)行離散化，即將求解區(qū)域劃分為許多幾何形狀簡(jiǎn)單規(guī)那么的單元，在二維一般是三角形或四邊形，在三維是四面體或六面體。然后，在每個(gè)單元內(nèi)，用一個(gè)比擬簡(jiǎn)單的解析函數(shù)來(lái)逼近微分方程的解，此函數(shù)在單元內(nèi)用一組選定的單元基函數(shù)的線性組合表示，而其中的系數(shù)通常是節(jié)點(diǎn)參數(shù)，它是待定的。這樣，每個(gè)單元只要有適當(dāng)數(shù)量的節(jié)點(diǎn)參數(shù)值，就可以滿足對(duì)插值函數(shù)的光滑性和精度的要求。第三，在滿足微分方程和相應(yīng)的初邊值條件下，對(duì)全部子域進(jìn)行積分。對(duì)每個(gè)單元分別進(jìn)行積分，形成“單元方程”，通過(guò)總體合成，得到總體有限元方程組。最后，用適當(dāng)?shù)姆椒ń夥匠探M，可得節(jié)點(diǎn)參數(shù)值，進(jìn)而可求得各單元內(nèi)的近似解。本文將以圓柱繞流問(wèn)題為例，展示有限元法求解的一般步驟。一.有限元法概述有限元法〔finiteelementmethod〕是一種高效能、常用的計(jì)算方法。有限元法在早期是以變分原理為根底開(kāi)展起來(lái)的，所以它廣泛地應(yīng)用于以拉普拉斯方程和泊松方程所描述的各類物理場(chǎng)中,這類場(chǎng)與泛函的極值問(wèn)題有著緊密的聯(lián)系。有限單元法最早可上溯到20世紀(jì)40年代。Courant第一次應(yīng)用定義在三角區(qū)域上的分片連續(xù)函數(shù)和最小位能原理來(lái)求解St.Venant扭轉(zhuǎn)問(wèn)題?，F(xiàn)代有限單元法的第一個(gè)成功的嘗試是在1956年，Turner、Clough等人在分析飛機(jī)結(jié)構(gòu)時(shí)，將鋼架位移法推廣應(yīng)用于彈性力學(xué)平面問(wèn)題，給出了用三角形單元求得平面應(yīng)力問(wèn)題的正確答案。1960年，Clough進(jìn)一步處理了平面彈性問(wèn)題，并第一次提出了"有限單元法"，使人們認(rèn)識(shí)到它的成效。自從1969年以來(lái)，某些學(xué)者在流體力學(xué)中應(yīng)用加權(quán)余量法中的迦遼金法(Galerkin)或最小二乘法等同樣獲得了有限元方程。具體而言，參考差分法中網(wǎng)格化的做法，把求解區(qū)域劃分為有限多子區(qū)域，稱這些子區(qū)域?yàn)閱卧?，在每個(gè)單元上構(gòu)造解的近似分布，將Ritz法或加權(quán)余量法應(yīng)用到分塊的逼近函數(shù)上。因而有限元法可應(yīng)用于以任何微分方程所描述的各類物理場(chǎng)中，而不再要求這類物理場(chǎng)和泛函的極值問(wèn)題有所聯(lián)系。實(shí)質(zhì)上，有限元法就是Ritz法或加權(quán)余量法。二．問(wèn)題描述 考慮位于兩塊無(wú)限長(zhǎng)平板間的圓柱體的平面繞流問(wèn)題，幾何尺寸如下列圖所示，來(lái)流為vx=1,vy把它作為有限元的求解區(qū)域Ω。要求求解出整個(gè)區(qū)域中的流函數(shù)、vx、三．提出物理問(wèn)題的邊界條件和滿足的微分方程1.邊界ab為流線，取ψ=0，?2.邊界bc也為流線，同樣取ψ=0，?3.邊界cd，切向速度vτ=0，?ψ?n4.邊界de為流線,滿足ψ于是在ed上，ψ=2，?5.進(jìn)口邊界ae上，ψ=ψa+a也可以提自然邊界條件?ψ我們以流函數(shù)ψ作為未知函數(shù)來(lái)解此問(wèn)題，流函數(shù)所滿足的微分方程如下：?2ψ=0此處Γ1指ab,bc,de和ae四段邊界，而Γ2就是就是cd段邊界，且切向速度v四．有限元法解圓柱繞流問(wèn)題1．建立有限元積分表達(dá)式根據(jù)求解問(wèn)題的根本控制方程，應(yīng)用變分法或加權(quán)余量法將求解的微分方程定解問(wèn)題化為等價(jià)的積分表達(dá)式，作為有限元法求解問(wèn)題的出發(fā)方程式。對(duì)于方程〔1〕，它是一橢圓型方程，具有正定性，可以用變分法，這里直接給出泛函J令其變分δJ=0，可以得到Ω自然邊界條件已經(jīng)包含在變分表達(dá)式中(其名稱的由來(lái))，而本質(zhì)邊界條件必須強(qiáng)制ψ滿足(因此稱其為本質(zhì)邊界條件，也稱為強(qiáng)制邊界條件)。如果根據(jù)原微分方程中無(wú)法給出泛函J，那么可以用Galerkin加權(quán)余量方法得到積分方程，這相當(dāng)于將原來(lái)的微分方程寫(xiě)為如下變分形式：Ω這里的δψ是函數(shù)ψ的改變量，是一種“虛位移”，在本質(zhì)邊界條件Γ1上為零。因此，上式做分部積分后，邊界積分僅剩下Ω即〔3〕式?？梢?jiàn)，如果ψ滿足原來(lái)的微分方程和邊界條件，那么，必然有ψ滿足(4)式，進(jìn)而滿足(5)式。注意，在(5)式中，包含的邊界Γ2上的邊界條件信息，對(duì)邊界Γ1的局部，僅知道它是給定了函數(shù)值的邊界，卻不知道邊界上的值是多少，為了確定這些值，還需要額外的處理方法。正是因?yàn)棣?上的邊界信息可以包含在積分表達(dá)式中，這種邊界條件也稱為自然邊界條件。2．區(qū)域剖分根據(jù)物理問(wèn)題的特點(diǎn)以及區(qū)域的形狀，把計(jì)算區(qū)域分成許多幾何形狀規(guī)那么但大小可以不同的單元，確定單元節(jié)點(diǎn)的數(shù)目和位置，建立表示網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)。采用的單元形狀和節(jié)點(diǎn)的分布，以及插值函數(shù)的選取還應(yīng)考慮到計(jì)算精度和可微性的要求。在此題中，網(wǎng)格〔區(qū)域剖分〕局部已由GRID.DAT文件給出。具體而言，網(wǎng)格將求解區(qū)域分為545個(gè)節(jié)點(diǎn)和997個(gè)單元。文件前半局部給出每個(gè)節(jié)點(diǎn)的橫縱坐標(biāo)，后半局部給出每個(gè)單元對(duì)應(yīng)的三個(gè)節(jié)點(diǎn)的序號(hào)。除此之外，單元剖分還應(yīng)該建立本質(zhì)邊界條件和自然邊界條件的節(jié)點(diǎn)表。自然邊界條件對(duì)解此題并無(wú)奉獻(xiàn)〔原因在下文會(huì)進(jìn)行說(shuō)明〕，由網(wǎng)格表特征可知，此題的本質(zhì)邊界條件為abc段上〔節(jié)點(diǎn)1-37〕，psi為零；de段上〔節(jié)點(diǎn)46-75〕，psi為2；ea段上〔節(jié)點(diǎn)76-91〕，psi為縱坐標(biāo)y。單元剖分實(shí)質(zhì)上就是建立這幾個(gè)對(duì)應(yīng)關(guān)系。3.單元分析單元分析的目的是建立有限元方程。把有限元積分表達(dá)式(3)寫(xiě)為各個(gè)單元求和的形式e=1這里Ω（e）Ω其中Γ（e）表示單元e的邊界，上式實(shí)質(zhì)上并不是一個(gè)等式，只具有形式上的意義，當(dāng)對(duì)所有的單元求和以后，才是等式。如果把線積分中的Γ2∩Γ(e)換為Γ(e)，那么得到的是等式，但在對(duì)所有單元求和時(shí)，內(nèi)部邊界的線積分剛好抵消，因此(7)也可以理解為不計(jì)內(nèi)部邊界奉獻(xiàn) 流函數(shù)ψ在單元e內(nèi)可用如下函數(shù)近似：ψ=這里ψi(i=1,2,3)為節(jié)點(diǎn)流函數(shù)值，Ni為節(jié)點(diǎn)上的插值函數(shù)，上式中重復(fù)下標(biāo)表示約定求和。將(8)代入(7Ω由于δψΩ此即單元方程，通常可以簡(jiǎn)寫(xiě)為A采用三節(jié)點(diǎn)三角形單元時(shí)，單元的插值基函數(shù)為N如果單元e三個(gè)點(diǎn)坐標(biāo)為〔xie,yiN即插值基函數(shù)Ni在xie點(diǎn)取1，在xjexke對(duì)某一固定的單元e，將〔11〕式代入〔10〕中，可以得到：A此即采用線性單元時(shí)的單元方程系數(shù)矩陣。其中A(e)為三角形〔積分區(qū)域〕的面積，bc的值可由〔12〕〔13〕A求解單元系數(shù)矩陣時(shí)，一般同時(shí)進(jìn)行總體合成，每形成一個(gè)單元方程，便把它累加到總體方程中。出于順序和邏輯上的考慮，下一步再詳細(xì)說(shuō)明總體合成的方法。對(duì)于邊界積分項(xiàng)，我們假設(shè)三角形單元e中序號(hào)為，的節(jié)點(diǎn)在邊界上，為自然邊界，其長(zhǎng)度為l。首先，注意到插值函數(shù)在邊上是零。所以，可以得到如下結(jié)論：圖：自然邊界條件的處理 右端項(xiàng)：。 為了計(jì)算和，以點(diǎn)為原點(diǎn)，沿直線建立局部坐標(biāo)系，在此坐標(biāo)系中，插值函數(shù)和如上圖所示，可寫(xiě)成線性插值函數(shù)如下： 假設(shè)切向速度在兩節(jié)點(diǎn)處的值分別為和，并且沿邊界是線性分布的，可以表示為。于是可以得到。對(duì)于前面討論的圓柱繞流問(wèn)題，由于，所以，根據(jù)線性解的性質(zhì)，必有f無(wú)需考慮f的影響，使程序得到了不少的簡(jiǎn)化。4.總體合成總體合成的過(guò)程就是把已經(jīng)形成好的單元方程按一定順序迭加起來(lái)，形成總體有限元方程。具體做法是根據(jù)單元內(nèi)節(jié)點(diǎn)的總體順序號(hào)，把單元方程進(jìn)行延拓，未知量包含所有節(jié)點(diǎn)上的函數(shù)值，與此單元無(wú)關(guān)的位置以零填充，把所有的單元方程都進(jìn)行延拓以后，進(jìn)行系數(shù)矩陣的累加，便得到總體方程。理論上說(shuō)，這一過(guò)程也可以通過(guò)引入一個(gè)Boole型矩陣來(lái)實(shí)現(xiàn)，定義單元e的boole矩陣,i=1,2,3;j=1,2,3…Np.矩陣B其實(shí)就是單元節(jié)點(diǎn)序號(hào)表的又一表達(dá)形式。單元e的系數(shù)矩陣以及右端項(xiàng)沿拓后就是：進(jìn)而總體合成的過(guò)程可以表示為，。但是這種方法比擬麻煩，要重新定義新的矩陣Ｂ，而且還要涉及計(jì)算矩陣轉(zhuǎn)置和矩陣相乘等運(yùn)算，一方面計(jì)算量較大，并且浪費(fèi)空間，另一方面人為地增加了程序的復(fù)雜性，降低了程序的可讀性。因此，這種方法一般只用作理論分析。實(shí)際的計(jì)算中，每當(dāng)計(jì)算出一個(gè)單元系數(shù)矩陣Aij(e)，假設(shè)單元e三個(gè)節(jié)點(diǎn)編號(hào)分別為i,j,k，那么將Aij(e)中的1*1項(xiàng)放入大矩陣(借鑒結(jié)構(gòu)力學(xué)的概念，不妨稱其為總體“剛度”矩陣，下同)A的i*i項(xiàng)中，將1*2,1*3分別放入總體剛度矩陣A的i*j,i*k項(xiàng)中。同理，2*1,2*2，2*3,3*1,3*2,3*3項(xiàng)分別對(duì)應(yīng)總體剛度矩陣A的j*i,j*j,j*k,k*i,k*j,k*k項(xiàng)中。采用此方法并未多占用計(jì)算機(jī)內(nèi)存，運(yùn)算量也并不大〔總共進(jìn)行9*e次加法運(yùn)算，不進(jìn)行乘法運(yùn)算〔5〕邊界條件處理 這里的邊界條件是指本質(zhì)邊界條件，自然邊界條件已經(jīng)包含在積分表達(dá)式中了。具體做法是將的值代入到方程組中，并把的值移到方程組的右段，形成右端項(xiàng)。〔6〕求解有限元方程組并計(jì)算相關(guān)物理量有限元方程組的求解是一個(gè)代數(shù)問(wèn)題，應(yīng)針對(duì)具體的問(wèn)題采用適宜的方法求解。對(duì)于對(duì)角占優(yōu)的代數(shù)方程組，可以采用迭代法求解，規(guī)模不大的可以用Gauss消元法一類的直接法求解，三對(duì)角方程那么可以用追趕法。求出所有的待求量后，便得到了近似函數(shù)的表達(dá)式，并可以計(jì)算出相關(guān)的物理量。對(duì)計(jì)算結(jié)果進(jìn)行綜合的分析，以期得到原問(wèn)題的正確的物理解答。對(duì)于每個(gè)單元，速度可以根據(jù)vv來(lái)計(jì)算,節(jié)點(diǎn)上的速度值可取這個(gè)節(jié)點(diǎn)相鄰單元的速度值的平均。節(jié)點(diǎn)上的壓力值可以有伯努利方程計(jì)算。假設(shè)求解區(qū)域位于同一水平面內(nèi)，介質(zhì)密度ρ=1，來(lái)流壓力p=0，那么p=12五．具體算法實(shí)現(xiàn)(算法說(shuō)明)1.建立有限元積分表達(dá)式 此局部不在算法中表達(dá)。2.區(qū)域剖分區(qū)域剖分的主要內(nèi)容已由網(wǎng)格文件完成，我們讀入網(wǎng)格文件，并設(shè)置邊界條件即可。 從GRID.TXT中讀入所需根本信息：節(jié)點(diǎn)數(shù)、單元數(shù)、節(jié)點(diǎn)坐標(biāo)、單元包含節(jié)點(diǎn)及其次序。為了從零開(kāi)始編號(hào)，應(yīng)將單元包含節(jié)點(diǎn)序號(hào)統(tǒng)一減一。并根據(jù)網(wǎng)格特征，將流函數(shù)psi賦初值：abc段上〔節(jié)點(diǎn)1-37〕，psi為零；de段上〔節(jié)點(diǎn)46-75〕，psi為2；ea段上〔節(jié)點(diǎn)76-91〕，psi為縱坐標(biāo)y。3.單元分析與總體合成 按公式〔13〕求得單元i的插值基函數(shù)的系數(shù)b[i][s]c[i][s](s=1,2,3)和系數(shù)A(e)，求出bc進(jìn)而求得Aij4.邊界條件處理 逐行處理。當(dāng)行數(shù)(i>=0&&i<37)||(i>=45&&i<91)，即為本質(zhì)邊界條件的行數(shù)時(shí)，直接將此行所有除ψ[i]對(duì)應(yīng)的系數(shù)以外的其他系數(shù)全部置零，psi[i]的系數(shù)置為一，右端項(xiàng)置為ψ[i]的值。當(dāng)行數(shù)不為本質(zhì)邊界條件的行數(shù)時(shí)，通過(guò)f[i]=f[i]-a[i][j]*psi[j]求右端項(xiàng)，再將除ψ[i]對(duì)應(yīng)的系數(shù)以外的其他系數(shù)全部置零。5.求解有限元方程組并計(jì)算相關(guān)物理量。 本算法采用LU分解進(jìn)行求解。其中，double**LU(double**a,intm,double**l)為L(zhǎng)U分解函數(shù)。**a即二維數(shù)組〔矩陣〕a的首地址,m為矩陣規(guī)模，此函數(shù)返回LU矩陣的首地址。double*SOVLELU(double**l,double*b,intm,double*x)解LU分解后的方程組。l為L(zhǎng)U分解產(chǎn)生的矩陣，b為右端項(xiàng)，m為矩陣規(guī)模，x為解向量,返回解向量的首地址。使用這兩個(gè)函數(shù)即可求出所有節(jié)點(diǎn)處流函數(shù)ψ的值。單元速度的求解那么采用〔14〕〔15〕兩式，為了節(jié)省空間，不保存單元節(jié)點(diǎn)速度，每求出一個(gè)單元節(jié)點(diǎn)速度，立刻存入其包含的三個(gè)節(jié)點(diǎn)中。使用count計(jì)數(shù)，每存一次count加一。最后將每個(gè)節(jié)點(diǎn)的速度除以count，即取平均，得每個(gè)節(jié)點(diǎn)的速度。壓力使用伯努利方程較易得到，此處不加贅述。6.輸出。將結(jié)果〔節(jié)點(diǎn)坐標(biāo)、xy方向速度、流函數(shù)ψ的值、壓強(qiáng)、單元包含節(jié)點(diǎn)數(shù)〕信息輸出到result111.txt中，并釋放所有空間?！矠楣?jié)約內(nèi)存，局部空間已于之前釋放〕。更多詳細(xì)說(shuō)明請(qǐng)參見(jiàn)源程序中的注釋局部。附源程序：#include<stdio.h>#include<math.h>#include<stdlib.h>voidmain(){double**LU(double**a,intm,double**l);//LU分解函數(shù),**a即二維數(shù)組〔矩陣〕a的首地址,此函數(shù)返回LU矩陣的首地址，m為矩陣規(guī)模double*SOVLELU(double**l,double*b,intm,double*x);//解LU分解后的方程組，l為L(zhǎng)U分解產(chǎn)生的矩陣，b為右端項(xiàng)，m為矩陣規(guī)模，x為解向量,返回解向量的首地址//voidL_U(double**A,double*f,intm);FILE*fp;//文件指針fp=fopen("GRID.txt","r");//將c\從GRID.txt中讀取數(shù)據(jù)intn,e;//n為節(jié)點(diǎn)個(gè)數(shù)，e為單元個(gè)數(shù)inti,j;fscanf(fp,"%d",&n);//從文件中讀入結(jié)點(diǎn)個(gè)數(shù)fscanf(fp,"%d",&e);//從文件中讀入單元個(gè)數(shù)double*x=newdouble[n];//x[i]即第i+1個(gè)節(jié)點(diǎn)的橫坐標(biāo)double*y=newdouble[n];//y[i]即第i+1個(gè)節(jié)點(diǎn)的縱坐標(biāo)double**l=newdouble*[n];//l即LU分解產(chǎn)生的中間矩陣for(i=0;i<n;i++)l[i]=newdouble[n];for(i=0;i<n;i++){fscanf(fp,"%lf",&x[i]);//從文件中讀入結(jié)點(diǎn)橫坐標(biāo)fscanf(fp,"%lf",&y[i]);//從文件中讀入結(jié)點(diǎn)縱坐標(biāo)}int*first=newint[e];//first[]即三角形單元對(duì)應(yīng)的第一個(gè)節(jié)點(diǎn)編號(hào)(從一開(kāi)始編號(hào))int*second=newint[e];//second[]即三角形單元對(duì)應(yīng)的第二個(gè)節(jié)點(diǎn)編號(hào)int*third=newint[e];//third[]即三角形單元對(duì)應(yīng)的第三個(gè)節(jié)點(diǎn)編號(hào)for(i=0;i<e;i++){fscanf(fp,"%d",&first[i]);//從文件中讀入單元對(duì)應(yīng)的第一個(gè)節(jié)點(diǎn)編號(hào)first[i]=first[i]-1;//減一使節(jié)點(diǎn)編號(hào)從零開(kāi)始，便于后面的計(jì)算fscanf(fp,"%d",&second[i]);//從文件中讀入單元對(duì)應(yīng)的第二個(gè)節(jié)點(diǎn)編號(hào)second[i]=second[i]-1;fscanf(fp,"%d",&third[i]);//從文件中讀入單元對(duì)應(yīng)的第三個(gè)節(jié)點(diǎn)編號(hào)third[i]=third[i]-1;}fclose(fp);//至此網(wǎng)格文件讀入完畢double*psi=newdouble[n];//psi[]即節(jié)點(diǎn)上的流函數(shù)值,開(kāi)始時(shí)將其置零for(i=0;i<n;i++)psi[i]=0;double**a=newdouble*[n];//a[][]即系數(shù)矩陣〔類似剛度矩陣〕,大小為n*n,并將其置零for(i=0;i<n;i++)a[i]=newdouble[n];for(i=0;i<n;i++)for(j=0;j<n;j++)a[i][j]=0;double*f=newdouble[n];//f為方程右端項(xiàng),將其置零for(i=0;i<n;i++)f[i]=0;//以下開(kāi)始賦邊界條件,自然邊界條件在此題中可以不賦值〔為零，不影響后繼計(jì)算〕for(i=0;i<37;i++)psi[i]=0;//abc段上，psi為零for(i=45;i<75;i++)psi[i]=2;//de段上，psi為2for(i=75;i<91;i++)psi[i]=y[i];//ea段上，psi為縱坐標(biāo)y//邊界條件賦值完畢double**b=newdouble*[e];for(i=0;i<e;i++)b[i]=newdouble[3];double**c=newdouble*[e];for(i=0;i<e;i++)c[i]=newdouble[3];doubled;//b[e][3],c[e][3]為單元插值函數(shù)N(x,y)=a+b*x+c*y的系數(shù),d為系數(shù)A〔e〕for(i=0;i<e;i++)//對(duì)每一個(gè)單元，求出系數(shù)矩陣并參加大的系數(shù)矩陣中{d=(x[second[i]]-x[first[i]])*(y[third[i]]-y[first[i]])-(y[second[i]]-y[first[i]])*(x[third[i]]-x[first[i]]);d=d/2;b[i][0]=(y[second[i]]-y[third[i]])/(2*d);b[i][1]=(y[third[i]]-y[first[i]])/(2*d);b[i][2]=(y[first[i]]-y[second[i]])/(2*d);c[i][0]=(x[third[i]]-x[second[i]])/(2*d);c[i][1]=(x[first[i]]-x[third[i]])/(2*d);c[i][2]=(x[second[i]]-x[first[i]])/(2*d);a[first[i]][first[i]]=a[first[i]][first[i]]+(b[i][0]*b[i][0]+c[i][0]*c[i][0])*d;a[first[i]][second[i]]=a[first[i]][second[i]]+(b[i][0]*b[i][1]+c[i][0]*c[i][1])*d;a[first[i]][third[i]]=a[first[i]][third[i]]+(b[i][0]*b[i][2]+c[i][0]*c[i][2])*d;a[second[i]][first[i]]=a[second[i]][first[i]]+(b[i][0]*b[i][1]+c[i][0]*c[i][1])*d;a[second[i]][second[i]]=a[second[i]][second[i]]+(b[i][1]*b[i][1]+c[i][1]*c[i][1])*d;a[second[i]][third[i]]=a[second[i]][third[i]]+(b[i][1]*b[i][2]+c[i][1]*c[i][2])*d;a[third[i]][first[i]]=a[third[i]][first[i]]+(b[i][0]*b[i][2]+c[i][0]*c[i][2])*d;a[third[i]][second[i]]=a[third[i]][second[i]]+(b[i][1]*b[i][2]+c[i][1]*c[i][2])*d;a[third[i]][third[i]]=a[third[i]][third[i]]+(b[i][2]*b[i][2]+c[i][2]*c[i][2])*d;}//至此左端剛度矩陣已經(jīng)生成完畢，接下來(lái)根據(jù)本質(zhì)邊界條件將其化簡(jiǎn)，并生成化簡(jiǎn)后的右端項(xiàng)ffor(i=0;i<n;i++){if((i>=0&&i<37)||(i>=45&&i<91)){for(j=0;j<n;j++){if(i!=j)a[i][j]=0;elsea[i][j]=1;}f[i]=psi[i];}else{for(j=0;j<n;j++){f[i]=f[i]-a[i][j]*psi[j];if((j>=0&&j<37)||(j>=45&&j<91))a[i][j]=0;}}//if(i>=0&&i<37)f[i]=0;//保證邊界條件的局部是精確的〔不會(huì)產(chǎn)生某個(gè)小數(shù)從而影響計(jì)算精度〕}//至此簡(jiǎn)化完畢，化為一個(gè)對(duì)稱方程組下面開(kāi)始用LU分解解此線性方程組l=LU(a,n,l);psi=SOVLELU(l,f,n,psi);//L_U(a,f,n);//for(i=0;i<n;i++)// psi[i]=f[i];//至此已解出流函數(shù)psi,釋放一批空間并且新開(kāi)一批空間for(i=0;i<n;i++)delete[]l[i];//釋放LU分解的過(guò)渡矩陣delete[]l;for(i=0;i<n;i++)delete[]a[i];//釋放剛度矩陣delete[]a;delete[]f;//釋放右端項(xiàng)double*vx=newdouble[n];//vx[i]即第i+1個(gè)節(jié)點(diǎn)的x方向速度double*vy=newdouble[n];//vy[i]即第i+1個(gè)節(jié)點(diǎn)的y方向速度double*count=newdouble[n];//count[i]記錄計(jì)算速度時(shí)每個(gè)節(jié)點(diǎn)上速度加和的次數(shù)，以便之后取平均for(i=0;i<n;i++){vx[i]=0;vy[i]=0;count[i]=0;}//開(kāi)始求流場(chǎng)中各節(jié)點(diǎn)的速度doublexx,yy;//xx,yy暫時(shí)儲(chǔ)存每個(gè)單元x方向與y方向的速度f(wàn)or(i=0;i<e;i++)//為節(jié)約空間，不儲(chǔ)存每個(gè)單元的速度，求出每個(gè)單元的速度后直接分配到節(jié)點(diǎn)上{xx=c[i][0]*psi[first[i]]+c[i][1]*psi[second[i]]+c[i][2]*psi[third[i]];yy=-b[i][0]*psi[first[i]]-b[i][1]*psi[second[i]]-b[i][2]*psi[third[i]];vx[first[i]]=vx[first[i]]+xx;vx[second[i]]=vx[second[i]]+xx;vx[third[i]]=vx[third[i]]+xx;vy[first[i]]=vy[first[i]]+yy;vy[second[i]]=vy[second[i]]+yy;vy[third[i]]=vy[third[i]]+yy;count[first[i]]=count[first[i]]+1;count[second[i]]=count[second[i]]+1;count[third[i]]=count[third[i]]+1;}for(i=0;i<n;i++){vx[i]=vx[i]/count[i];vy[i]=vy[i]/count[i];}//最后求解壓力項(xiàng)double*p=newdouble[n];for(i=0;i<n;i++)p[i]=(1-vx[i]*vx[i]-vy[i]*vy[i])/2;//壓力項(xiàng)求解完畢//至此已求出每個(gè)節(jié)點(diǎn)上的速度，接下來(lái)釋放空間并進(jìn)行輸出delete[]count;for(i=0;i<e;i++)delete[]b[i];delete[]b;for(i=0;i<e;i++)delete[]c[i];delete[]c;fp=fopen("result111.txt","w");fprintf(fp,"TITLE=Finiteelementmethod\n");fprintf(fp,"FILETYPE=FULL\n");fprintf(fp,"VARIABLES=X,Y,vx,vy,psi,p\n");fprintf(fp,"ZONE\n");fprintf(fp,"N=%dE=%df=fepointet=triangle\n",n,e);for(i=0;i<n;i++)fprintf(fp,"%lf%lf%lf%lf%lf%lf\n",x[i],y[i],vx[i],vy[i],psi[i],p[i]);for(i=0;i<e;i++)fprintf(fp,"%d%d%d\n",first[i]+1,second[i]+1,third[i]+1);delete[]first;delete[]second;delete[]third;delete[]p;delete[]vx;delete[]vy;delete[]x;delete[]y;delete[]psi;fclose(fp);}double**LU(double**a,intm,double**l)//LU分解函數(shù),**a即二維數(shù)組〔矩陣〕a的首地址,此函數(shù)返回LU矩陣的首地址，m為矩陣規(guī)模{//l[][]即LU分解時(shí)產(chǎn)生的矩陣(為節(jié)約空間，將L和U放在同一矩陣中)inti,j,k;for(k=0;k<m;k++){for(j=k;j<m;j++){l[k][j]=a[k][j];for(i=0;i<k;i++){l[k][j]=l[k][j]-l[k][i]*l[i][j];}}for(j=k+1;j<m;j++){l[j][k]=a[j][k];for(i=0;i<k;i++){l[j][k]=l[j][k]-l[j][i]*l[i][k];}l[j][k]=l[j][k]/l[k][k];}}return(l);}double*SOVLELU(double**l,double*b,intm,double*x)//解LU分解后的方程組，l為L(zhǎng)U分解產(chǎn)生的矩陣，b為右端項(xiàng)，m為矩陣規(guī)模，x為解向量,返回解向量的首地址{inti,j;//先解Ly=bx[0]=b[0];for(i=1;i<m;i++){x[i]=b[i];for(j=0;j<i;j++)x[i]=x[i]-l[i][j]*x[j];}//再解Ux=yx[m-1]=x[m-1]/l[m-1][m-1];for(i=m-2;i>=0;i--){for(j=i+1;j<m;j++){x[i]=x[i]-l[i][j]*x[j];}x[i]=x[i]/l[i][i];}return(x);}六．可視化處理與結(jié)果分析。 繪制網(wǎng)格如下：流函數(shù)如下：流線圖如下：假設(shè)將流線和流函數(shù)畫(huà)在同一張圖上，那么為壓強(qiáng)分布如下：由圖像可見(jiàn)，速度分布比擬符合物理直觀，將流函數(shù)圖和流線圖畫(huà)在同一張圖像上，可以發(fā)現(xiàn)等ψ線就是流線，滿足物理要求。壓強(qiáng)圖像上，〔-1,0〕處速度為零，壓強(qiáng)最大，為滯止壓強(qiáng)；〔0,1〕處壓強(qiáng)最小，也符合物理要求。七．進(jìn)一步優(yōu)化的些許