有限元法解圓柱繞流問(wèn)題

1、有限元法求解無(wú)限流體中的圓柱繞流問(wèn)題2016年01月12號(hào)一問(wèn)題描述考慮位于兩塊無(wú)限長(zhǎng)平板間的圓柱體的平面繞流問(wèn)題，幾何尺寸如下圖所示，來(lái)流為vx=1,vy=0。由于流場(chǎng)具有上下左右的對(duì)稱性，只考慮左上角四分之一的計(jì)算區(qū)域abcde，把它作為有限元的求解區(qū)域。要求求解出整個(gè)區(qū)域中的流函數(shù)、vx、vy以及壓強(qiáng)值。圖1：圓柱繞流模型二數(shù)學(xué)建模 在足夠遠(yuǎn)前方選與來(lái)流垂直的控制面ae,cd是沿y軸，亦即一流動(dòng)對(duì)稱軸，bc是物面，ab亦是流動(dòng)對(duì)稱軸，所要考慮的流動(dòng)區(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.進(jìn)口邊界ae上，=a+ayvxdy=y（本文中采取此條件）也可以提自然邊界條件n=0，n=vx=1我們以流函數(shù)作為未知函數(shù)來(lái)解此問(wèn)題，流函數(shù)所滿足的微分方程如下：2=0 |1=（本質(zhì)邊界條件）n|2=-vs(自然邊界條件)（1）此處1指ab,bc,de和ae四段邊界，而2就是就是cd 段邊界，且切向速度vs= 0，1 和2 合起來(lái)是整個(gè)邊界，并且此二者不重合。下面，按有限元方法的一般步驟來(lái)計(jì)算此問(wèn)題。三有限元法解圓柱繞流問(wèn)題1建立有限元積分表達(dá)式根據(jù)求解問(wèn)題的基本控制方

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

4、，上式做分部積分后，邊界積分僅剩下2的部分。具體為d+2vsd=0(5)即（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)的分布，以及插值

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

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

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

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

9、線建立局部坐標(biāo)系，在此坐標(biāo)系中，插值函數(shù)和如上圖所示，可寫成線性插值函數(shù)如下：假設(shè)切向速度在兩節(jié)點(diǎn)處的值分別為和，并且沿邊界是線性分布的，可以表示為。于是可以得到。對(duì)于前面討論的圓柱繞流問(wèn)題，由于，所以，根據(jù)線性解的性質(zhì)，必有f（e）=f（e）=f（e）=0無(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ò)程也可以

10、通過(guò)引入一個(gè)Boole 型矩陣來(lái)實(shí)現(xiàn)，定義單元e的boole矩陣,i=1,2,3; j=1,2,3Np.矩陣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é)的概念，不妨稱其為總體“剛

11、度”矩陣，下同)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)題采用

12、合適的方法求解。對(duì)于對(duì)角占優(yōu)的代數(shù)方程組，可以采用迭代法求解，規(guī)模不大的可以用Gauss 消元法一類的直接法求解，三對(duì)角方程則可以用追趕法。求出所有的待求量后，便得到了近似函數(shù)的表達(dá)式，并可以計(jì)算出相關(guān)的物理量。對(duì)計(jì)算結(jié)果進(jìn)行綜合的分析，以期得到原問(wèn)題的正確的物理解答。對(duì)于每個(gè)單元，速度可以根據(jù)vx=y=iNiy=iNiy=cii(14)vy=-x=-iNix=-iNix=-bii(15)來(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=121-vi2,i=1,210。如此便可得到節(jié)

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

14、析解相對(duì)比，得到圖7 中的曲線?？梢?jiàn)，有限元法與解析法所得曲線在趨勢(shì)上基本一致，但在數(shù)值上最大誤差為11%。圖6: 1/4區(qū)域右邊界上流函數(shù)的有限元數(shù)值解與解析解對(duì)比附件：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; %單元與結(jié)點(diǎn)號(hào)、面積對(duì)應(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

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; %結(jié)點(diǎn)與坐標(biāo)對(duì)應(yīng)關(guān)系矩陣Enum,temp=size(ENA); % 返回單元總數(shù)Nnum,temp=size(NCORD); % 返回結(jié)點(diǎn)總數(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é)點(diǎn)編號(hào)temp,Bnum=size(BNODE); % 返回邊界結(jié)點(diǎn)數(shù)%邊界已知流函數(shù)的結(jié)點(diǎn)

34、號(hào)及其值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ù)的結(jié)點(diǎn)數(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ù)的結(jié)點(diǎn)數(shù)for k=1:Enum %單元號(hào)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);%各個(gè)單元的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,

37、m);endendendendA(i,j)=temp;endendfor i=1:Uknum %掃描未知流函數(shù)的結(jié)點(diǎn)F(i)=0;for j=1:UknumANEW(i,j)=A(UKN(i),UKN(j); %新的系數(shù)矩陣endfor k=1:Bknum %掃描已知流函數(shù)的結(jié)點(diǎn)F(i)=F(i)+A(UKN(i),BKN(1,k)*BKN(2,k);endF(i)=-F(i); %移項(xiàng)到等號(hào)右側(cè)成新的乘積系數(shù)向量end%FAIUN=ANEWF' %解向量，未知結(jié)點(diǎn)的流函數(shù)值Q,R=lu(ANEW);FAIUN=R(QF');%將所有結(jié)點(diǎn)的流函數(shù)值寫入FAIN 向量for i=1:UknumFAIN(UKN(i)=FAIUN(i);endfor i=1:BknumFAIN(BKN(1,i)=BKN(2,i)