有限元法解二维圆柱绕流问题.pdf

1、 目 录 目 录. 1 1 问题描述. 1 2 相关的有限元理论基础 . 1 2.1二次泛函极值原理和里兹解法 1 2.2 伽辽金加权余数法 . 2 3 二维圆柱绕流的有限元解法 . 3 3.1 圆柱绕流问题的数学模型 . 3 3.2 数学模型的有限元法求解 . 4 3.3有限元法求解的MATLAB实现 8 3.4 结果讨论 . 21 参 考 文 献 23 1 1 问题描述 图 1 二维圆柱绕流问题 本作业尝试探索二维圆柱绕流的流函数的有限元数值解法及并绘出其流线图。 2 相关的有限元理论基础 2.1二次泛函极值原理和里兹解法 设微分方程 ,( ) Au f 在区域 中 ， (1) 式中 A

2、为线性微分算子，如拉普拉斯算子等，它作用于函数 u 后产生函数 f，这里 f 是给定 的已知函数。如果在区域内任意两函数 u、v 满足如下内积等式，则称算子 A为对称的，即 (, ) uA v u A v d (2) 如果(,u ) 0 uA ，那么称 A是正定的。 如果算子 A是对称正定的，可以证明使二次泛函 J(u)=(Au, u)-2(f, u) (3) 为极小值的 u 是微分方程(1)的解。因此，解微分方程(1)可以转化为解泛函(3)的极小值问题。 而如果算子 A是负定的，即(, ) 0 uA u ，那么泛函(3)的极大值 u 即为微分方程(1)的解。 例如对于泊松方程第一边值问题。

3、22 22 -, 0 s uu xy u ( + )=f(x,y)x, y(4) 这里 S为区域 的边界。 2 可以证明，对于泊松方程第一边值问题，算子 2 A 是对称的，正定的。 这样可以按照(3)式建立泛函并整理 2 2 J (Au, u) 2(f, u)=- 2= 2 uu d u f d ud fu d (u)=(5) 下面用里兹法求泛函(3)式的极值。所谓里兹解法，就是把 u 用一组线性无关的函数组 合来近似，即设 i u , 1,2,. , i cin (6) 1 ， 2 ， n 是一组线性无关函数，而 1 c ， 2 c ， n c 是一组待定系数，这一组待定系 数确定后，u 就

4、确定了，继而(5)式的解就可以确定了。 把 i u i c 带入(3)式 ji i J(u)=(Au, u)-2(f, u)= ( , ) 2 ( , ) iji cc A c f ，要使 J(u)有极小值，则有 J(u) =0 i c ， 即 iji ( , ) =( , ) j 1,2,. , j A cf n ，i ， (7) (7)式是关于 j c 的 n 阶线性方程组，解此线性方程组，得到一组 j c ，进而可以得到(5)式的 u 近似解 u。 2.2 伽辽金加权余数法 对于里兹法，其算子必须是正定或者负定，才有其对应的泛函。但是一般算子不能满足 这一条件。 为此引入加权余数法。 所

5、为加权余数法就是仍然用一组线性无关函数组合来近似， 因为是近似，所有就有误差 ，我们把误差 称为余数，然后令其误差 与某权函数的内积 为零，即成为加权余数法。权函数有不同取法，就引出了各种加权余数法。 研究方程 Au-f 0 (8) 仍用一组线性无关函数来近似，令 i u , 1,2,. , i cin (9) 3 把(9)式带入(8)式，则 Au-f 0 (10) 令 与权函数 内积为零，即 (,) = 0 (11) 伽辽金加权余数法就是取基函数 i 为权函数，即 iii ( , )=(A -f, )=0 i c (12) 这样，同样可以得到一组线性无关方程组，解得相应的一组待定系数，继而可

6、以得到(8)式 的 u 近似解 u。 3 二维圆柱绕流的有限元解法 3.1 圆柱绕流问题的数学模型 假设流体不可压缩，则二维圆柱绕流的流函数 满足拉普拉斯方程，即 22 22 =0 xy + (13) 图 2 二维圆柱绕流问题的简化模型 为了便于分析，取圆柱半径为 1，流场为以圆柱圆心为中心的正方形，边长为 6，以正 方形流场的边界来代表无限远处的流场，无限远处均匀来流流速为 1，如图 2 所示。 考虑到流场的对称性，仅研究流场的左上角四分之一区域，把这四分之一的区域划分为 180 个单元和 103 个结点，所有单元均为三角形形状，如图 3 所示。 4 图 3 左上角四分之一区域的流场及其有限

7、单元划分 3.2 数学模型的有限元法求解 三角形单元内任一点 p(x,y)的函数值 可以写成 123 (,) xyxy (14) 把三个结点值带入(14)式，得 112 13 1 212 23 2 312 33 3 x y x y x y (15) 解(15)式，令 11 22 33 111 22 11 11 12221 2 3 33 33 22 333 11 211 222123 332 33 11 211 222123 332 33 1 1 1 1 111 1 111 1 1 111 1 111 1 xy Axy xy xy xy xy x y Ax y xy xy x y xy y yy

8、y Ay yyy y x xxx Ax xxx x (16) 5 1 1 2 2 3 3 A A A A A A (17) 将(17)式带入(14)式，结合和(16)式整理得 3 12 11 22 33 NN A AA xy AAA (18) 其中 111 12 2 33 11 222 2 33 11 333 32 2 1 1 1 1 1 1 1 1 1 1 1 1 xy ab xc y xy AA xy xy ab xcy xy AA xy xy ab xcy xy AA xy (19) 即 3 12 33 2123132 23 11 3231213 31 22 1312321 , , ,

9、NNN abxcy A ax yx ybyycxx ax yx ybyycxx ax yx ybyycxx (20) 对于圆柱绕流方程(13)及(18)式，根据(8)式和(9)式知其余数表达式为 22 22 x y + (21) 取基函数为 N ，则应用伽辽金加权余数法应用(12)式解(13)式并应用格林公式 22 22 (, ) ( ) () () NN NN N S NN MM NM S dxdy xy dS dxdy nx x y y dS dxdy nx x y y +(22) 6 可以得到元素方程组 () NM M N NN MM NM NN S AF Ad x d y xxyy F

10、d S n (23) 将(20)式带入(23)式，得 2 1 () N M NM NM A b b c c dxdy A (24) 上式中被积函数 NM NM bb cc 是常数，设三角形单元面积为，那么由(16)式知 A=2， 故有 22 1 1 12 12 13 13 22 21 21 2 2 13 23 22 31 31 31 32 3 3 1 () 4 1 4 N M NM NM Ab b c c b c bb cc bb cc bb cc b c bb cc bb cc bb cc b c (25) 注意，需要把单元结点按照逆时针编号（图 3 即满足） ，这样可以保证三角形面积为正，

11、 再把结点坐标代入(20)式中，可以求出 N b 和 N c 。 下面求 NN S Fd S n 。 图 4 三角形单元示意图 积分路线取逆时针方向，如图 4 所示。设 AB=l 1 ，BC=l 2 ，CA=l 3 。 在 AB 线上： 7 1 1 1 S l ， 2 1 S l ， 3 0 ； 在 BC 线上： 1 0 ， 2 2 1 S l ， 2 S l ； 在 CA线上： 1 3 S l ， 2 0 ， 3 3 1 S l 。 定义 q n ，并整理得 113 11 23 3 111 () 366 Fll ql ql q 同样，可利用对称性求得 F 2 和 F 3 ，即 113 11

12、23 3 212 21 12 3 323 32 23 1 111 () 366 111 () 366 111 () 366 Fll ql ql q Fll ql ql q Fll ql ql q (26) 现在已经求出单元方程组的系数 A NM 和 F N ，既可以汇成总体方程组 , , 1,2,., ij j i afi jn (27) 为了方便起见，总体系数 a ij ，f i 和单元系数 A NM ，F N 的关系可以用下式表示 () () () 1 () () 1 ij i L eee ij NM N M e L ee iN N e aA fF (28) 这里 () 1 i e N e

13、N i ，第 个单元第 个结点与总体第 个结点重合 不重合() M 1 j e e ，第 个单元第M个结点与总体第j个结点重合 不重合在应用式(28)的关系中可以发现一些规律，即对于非边界的结点，则 f i =0。而对于边界 上的结点 i，设它前后的结点分别为 j和 k，如图 5所示。 8 图 5 f i 计算图 f i 可以通过下式计算： 111 f( ) 366 ik ii j ik i ki j j llqlqlq (30) 具体到图 3 所示的圆柱绕流问题，还需要确定边界条件。下边界为对称线，1/4 圆柱边 界为物面线，一般称作零流线，即这两条线上 =0；设无限远处流速 v =1，故

14、1 x v y ， 故上边界和左边界线上满足 =y；而右边界线上有 0 nx y v ，即 q=0。 结合边界条件，则线性方程组式(27)可解。 3.3有限元法求解的MATLAB实现 主程序文件：youxianyuan.m %单元与结点号、面积对应关系矩阵 ENA=41 32 37 0.039798301 41 44 32 0.047529903 18 101 21 0.120156817 18 103 101 0.094963772 103 22 102 0.043498437 103 18 22 0.094243318 30 36 28 0.038054662 30 35 36 0.042

15、287521 29 2 12 0.038074502 29 6 2 0.048742904 30 7 17 0.05029162 30 8 7 0.037140964 9 19 88 76 0.057994041 19 20 88 0.099591586 101 102 99 0.046382644 101 103 102 0.044218891 44 94 5 0.047132639 44 41 94 0.055967406 35 52 38 0.036085513 35 16 52 0.034571557 29 44 6 0.044280863 29 32 44 0.041448145 9

16、7 20 21 0.106109192 1 3 93 0.078212613 4 5 94 0.054017838 102 100 99 0.034088039 101 98 21 0.050359216 22 100 102 0.048451565 99 98 101 0.045186467 98 97 21 0.070310629 23 96 22 0.097022072 96 100 22 0.062253661 92 95 100 0.0322137 95 99 100 0.031691011 91 98 99 0.039495829 90 97 98 0.038409873 97 8

17、9 20 0.062852728 4 87 3 0.050497897 93 24 1 0.094904332 84 96 23 0.060406457 96 92 100 0.036472833 80 95 92 0.032029748 10 95 91 99 0.029674024 91 90 98 0.041864912 90 89 97 0.045507083 89 88 20 0.06124771 83 13 19 0.093342661 74 94 41 0.05155156 94 87 4 0.039208883 82 93 3 0.058277066 81 24 93 0.06

18、3798865 86 85 23 0.047832817 85 84 23 0.050837396 84 92 96 0.038920434 92 71 80 0.041373741 80 91 95 0.033559408 79 90 91 0.042160395 78 89 90 0.044890766 77 88 89 0.041466167 83 14 13 0.060586214 14 75 15 0.051106938 74 87 94 0.041578206 87 82 3 0.057195134 82 81 93 0.045533139 24 86 23 0.084555213

19、 81 86 24 0.056241682 73 85 86 0.030172685 72 84 85 0.032814285 84 71 92 0.037263856 80 79 91 0.047262673 79 78 90 0.043683491 78 77 89 0.043858618 11 77 76 88 0.033382546 69 83 19 0.053167481 83 75 14 0.030301151 74 82 87 0.039566275 68 81 82 0.045303658 81 73 86 0.039028656 73 72 85 0.03280082 72

20、71 84 0.036391036 70 79 80 0.050409968 66 78 79 0.043207084 65 77 78 0.040783626 64 76 77 0.03638318 76 63 19 0.060458755 69 75 83 0.025977524 75 58 15 0.055241265 57 74 41 0.053616644 74 68 82 0.046847336 68 73 81 0.041793554 62 72 73 0.036075451 61 71 72 0.035939511 71 67 80 0.038472382 67 70 80 0

21、.040185579 59 70 56 0.045615694 70 66 79 0.04312042 66 65 78 0.041221814 65 64 77 0.038815812 64 63 76 0.035503832 63 69 19 0.057736703 69 58 75 0.032404442 57 68 74 0.04269458 12 68 62 73 0.039041739 62 61 72 0.037749165 61 67 71 0.033244944 67 56 70 0.038889688 59 66 70 0.040818003 55 65 66 0.0405

22、85183 54 64 65 0.04008813 53 63 64 0.042912099 63 58 69 0.036059255 57 62 68 0.041591565 51 61 62 0.039685976 61 56 67 0.037633909 56 60 59 0.036288344 50 60 47 0.031860179 49 59 60 0.031521434 59 55 66 0.039419363 55 54 65 0.03989979 54 53 64 0.040903813 53 58 63 0.046404299 58 48 15 0.066170074 43

23、 57 41 0.04763151 57 51 62 0.04195655 51 56 61 0.046636447 56 47 60 0.035592793 50 49 60 0.031796449 49 55 59 0.039229272 46 54 55 0.040974288 45 53 54 0.04180207 53 48 58 0.042196026 48 52 15 0.056305234 13 52 16 15 0.047154585 43 51 57 0.041090301 51 47 56 0.045821293 33 50 40 0.048398755 42 49 50

24、 0.040513663 49 46 55 0.040202562 46 45 54 0.04207801 45 48 53 0.0448719 48 38 52 0.037263593 43 47 51 0.041538735 47 40 50 0.039865602 33 42 50 0.050325014 42 46 49 0.041023271 39 45 46 0.04246825 45 38 48 0.043951365 44 5 6 0.051047374 43 40 47 0.041787595 34 42 33 0.051746427 42 39 46 0.042853263

25、 39 38 45 0.041743856 41 37 43 0.038609599 37 40 43 0.03483923 40 31 33 0.050574027 34 39 42 0.04382332 36 38 39 0.042210941 37 31 40 0.033185168 34 36 39 0.041698152 36 35 38 0.039494122 32 31 37 0.034698998 35 17 16 0.049191615 14 34 28 36 0.040629125 33 27 34 0.049558416 25 31 32 0.039551559 31 2

26、6 33 0.04815108 30 17 35 0.041250709 27 28 34 0.039846508 26 27 33 0.04614112 29 25 32 0.037979443 25 26 31 0.038580787 28 8 30 0.036116702 27 9 28 0.033481704 26 10 27 0.032757535 12 25 29 0.035427736 25 11 26 0.03337283 9 8 28 0.031904711 10 9 27 0.030445747 11 10 26 0.030285297 12 11 25 0.0318491

27、71 ; %结点与坐标对应关系矩阵 NCORD=3 0 1 0 2.6 0 2.2 0 1.8 0 1.4 0 0 1 0.258824103 0.965924471 0.499995466 0.866028022 0.707106781 0.707106781 15 0.866028022 0.499995465 0.965924471 0.258824102 0 3 0 2.6 0 2.2 0 1.8 0 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.7444651

28、1 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.503847145 1.428730143 1.548704503 0.736752996 0.489238864 1.743880003 0.764674661 1.

29、580844378 1.450757531 0.981852867 16 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

30、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.010914694 2.3291435

31、74 1.285646568 2.160870135 1.563524602 1.986279407 2.066343984 1.629034175 2.354650626 0.785729864 0.509120672 2.628036907 1.839751973 1.799640354 17 2.334608941 1.60379593 2.427838139 1.329969892 2.52208128 1.053358599 2.160447272 0.532346877 0.255534689 2.527394284 0.999394278 2.607756191 1.266887

32、377 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.535695128 2.604062411

33、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 18 2.360128075 2.679581821 2.634609866

34、 2.379744818 2.748684485 2.57733166 ; Enum,temp=size(ENA); % 返回单元总数 Nnum,temp=size(NCORD); % 返回结点总数 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); % 返回边界结点数 %边界已知流函数的结点号及其值 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

35、0 0 0 3 3 3 3 3 2.25 1.5 0.75; temp,Bknum=size(BKN); % 返回边界已知流函数的结点数 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,8 1,82,83,84,85,86,87,88,89,90,91,92,93,94,95,

36、96,97,98,99,100,101,102,103; temp,Uknum=size(UKN); %未知流函数的结点数 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*y3x3*y2;b(1)=y2y3;c(1)=x3x2; a(2)=x3*y1x1*y3;b(2)=y3y1;c(2)=x1x3; a(3)=x1*y2x2*y

37、1;b(3)=y1y2;c(3)=x2x1; for n=1:3 for m=1:3 ANM(k,n,m)=1/(4*ENA(k,4)*(b(n)*b(m)+c(n)*c(m);%各个单元的 Anm 信息 end 19 end end for i=1:Nnum %求系数矩阵 A for j=1:Nnum temp=0; for e=1:Enum for n=1:3 for m=1:3 if ENA(e,n)=i end end end end A(i,j)=temp; end end for i=1:Uknum %扫描未知流函数的结点 F(i)=0; for j=1:Uknum ANEW(i,

38、j)=A(UKN(i),UKN(j); %新的系数矩阵 end for k=1:Bknum %扫描已知流函数的结点 F(i)=F(i)+A(UKN(i),BKN(1,k)*BKN(2,k); end F(i)=F(i); %移项到等号右侧成新的乘积系数向量 end %FAIUN=ANEWF %解向量，未知结点的流函数值 Q,R=lu(ANEW); FAIUN=R(QF); 20 %将所有结点的流函数值写入 FAIN 向量 for i=1:Uknum FAIN(UKN(i)=FAIUN(i); end for i=1:Bknum FAIN(BKN(1,i)=BKN(2,i); end FAIN=

39、FAIN; scatter(NCORD(:,1),NCORD(:,2),5,FAIN(:,1) hold on; X,Y,Z=griddata(NCORD(:,1),NCORD(:,2),FAIN(:,1),linspace(3,0),linspace(0,3),v4);%插值 contour(X,Y,Z,15) axis(3 0 0 3) xlabel(x); ylabel(y); colorbar hold on; sita=pi/2:pi/20:pi; for r=0:0.001:1; plot(r*cos(sita),r*sin(sita),k); %中心点在原点，半径为 r 的圆 e

40、nd axis equal %右边界各点结点号 RBN(1,1)=13;RBN(1,2)=14;RBN(1,3)=15;RBN(1,4)=16;RBN(1,5)=17;RBN(1,6)=7; %右边界各点结点 y坐标 RBN(2,1)=NCORD(13,2);RBN(2,2)=NCORD(14,2);RBN(2,3)=NCORD(15,2); RBN(2,4)=NCORD(16,2);RBN(2,5)=NCORD(17,2);RBN(2,6)=NCORD(7,2); %右边界各点结点出流函数数值解 RBN(3,1)=FAIN(RBN(1,1);RBN(3,2)=FAIN(RBN(1,2);RB

41、N(3,3)=FAIN(RBN(1,3); RBN(3,4)=FAIN(RBN(1,4);RBN(3,5)=FAIN(RBN(1,5);RBN(3,6)=FAIN(RBN(1,6); 21 %右边界各点结点出流函数解析解 figure RBN(4,1)=fai(NCORD(13,1),NCORD(13,2);RBN(4,2)=fai(NCORD(14,1),NCORD(14,2); RBN(4,3)=fai(NCORD(15,1),NCORD(15,2);RBN(4,4)=fai(NCORD(16,1),NCORD(16,2); RBN(4,5)=fai(NCORD(17,1),NCORD(1

42、7,2);RBN(4,6)=fai(NCORD(7,1),NCORD(7,2); plot(RBN(2,:),RBN(3,:),r) hold on; plot(RBN(2,:),RBN(4,:),b) legend(有限元数值解,解析解,location,NorthWest) xlabel(右边界 y坐标); ylabel(流函数值); 子函数文件 fai.m function liuhanshu=fai(x,y); liuhanshu=y*(1-1/(x2+y2); 3.4 结果讨论 运行 3.3 小节中的 MATLAB 程序，可以得到图 6 和图 7 两张图。图 6 显示 1/4 区域的

43、流 场分布，图中的点以不同颜色表示各个结点处的流函数，图中的曲线为流函数的等势线，即 是流线，由于对称性，整个区域的流场分布可明显看出，故不再画出。而圆柱绕流问题的解 析解为 22 1 =y(1 ) x y ，便于与计算结果对比。从图 6 可以发现，有限元法数值法得到的 流线与解析法得到的流线在形态上基本一致。 图 6 有限元法计算得 1/4 区域的流场分布（左）与解析法得到的整个区域的流场分布（右） 22 具体到 1/4 区域的右边界上的结点，将有限元法得到的流函数数值解与解析解相对比， 得到图 7 中的曲线。可见，有限元法与解析法所得曲线在趋势上基本一致，但在数值上最大 误差为 11%。

44、图 7 1/4 区域右边界上流函数的有限元数值解与解析解对比 但如果将上边界和左边界的边界条件（即流函数值）按照解析解的精确值赋值，再次计 算 1/4 区域的右边界上的结点的有限元数值解，并绘出曲线如图 8 所示。 图 8 修正边界值后 1/4 区域右边界上流函数的有限元数值解与解析解对比 图 8 显示，修正边界值后，1/4 区域的右边界上的结点的有限元数值解与解析解几乎完 全重合。这说明，本文提出的有限元算法是正确的，单元与结点的划分是合适的，只是边界 条件的取值不太合理。进一步分析，这是由于流场区域较小（边长为 6） ，边界处仍然受到 圆柱的扰动，不能有效模拟无限远的情况。故从图 8 的结果可以推论，如果将流场区域扩张 到足够大，则可以计算出十分精确的结果。但由于结点和单元的几何信息的确定和录入非常23 繁杂，故本文不再赘述。 综上所述，本文提出的计算二维圆柱绕流的有限元数值方法可以与数值解较好吻合，这 种方法是有效合理的。 参 考 文 献 1 刘希云， 赵润祥. 流体力学中的有限元与边界元法M. 上海： 上海交通大学出版社， 1993.