有限元法解圆柱绕流问题

1、有限元法求解无限流体中的圆柱绕流问题2016年01月12号一问题描述考虑位于两块无限长平板间的圆柱体的平面绕流问题，几何尺寸如下图所示，来流为vx=1,vy=0。由于流场具有上下左右的对称性，只考虑左上角四分之一的计算区域abcde，把它作为有限元的求解区域。要求求解出整个区域中的流函数、vx、vy以及压强值。图1：圆柱绕流模型二数学建模 在足够远前方选与来流垂直的控制面ae,cd是沿y轴，亦即一流动对称轴，bc是物面，ab亦是流动对称轴，所要考虑的流动区域即由线abcdea所围成的区域，在这一区域中有：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我们以流函数作为未知函数来解此问题，流函数所满足的微分方程如下：2=0 |1=（本质边界条件）n|2=-vs(自然边界条件)（1）此处1指ab,bc,de和ae四段边界，而2就是就是cd 段边界，且切向速度vs= 0，1 和2 合起来是整个边界，并且此二者不重合。下面，按有限元方法的一般步骤来计算此问题。三有限元法解圆柱绕流问题1建立有限元积分表达式根据求解问题的基本控制方

3、程，应用变分法或加权余量法将求解的微分方程定解问题化为等价的积分表达式，作为有限元法求解问题的出发方程式。对于方程（1），它是一椭圆型方程，具有正定性，可以用变分法，这里直接给出泛函J（）=12d+2vsd=0（2）令其变分J=0，可以得到d+2vsd=0(3)自然边界条件已经包含在变分表达式中(其名称的由来)，而本质边界条件必须强制 满足(因此称其为本质边界条件，也称为强制边界条件)。如果根据原微分方程中无法给出泛函J，则可以用Galerkin 加权余量方法得到积分方程，这相当于将原来的微分方程写为如下变分形式：d=0(4)这里的是函数的改变量，是一种“虚位移”，在本质边界条件1上为零。因此

4、，上式做分部积分后，边界积分仅剩下2的部分。具体为d+2vsd=0(5)即（3）式。可见，如果满足原来的微分方程和边界条件，那么，必然有满足(4) 式，进而满足(5) 式。注意，在(5) 式中，包含的边界2 上的边界条件信息，对边界1 的部分，仅知道它是给定了函数值的边界，却不知道边界上的值是多少，为了确定这些值，还需要额外的处理方法。正是因为2 上的边界信息可以包含在积分表达式中，这种边界条件也称为自然边界条件。2区域剖分根据物理问题的特点以及区域的形状，把计算区域分成许多几何形状规则但大小可以不同的单元，确定单元节点的数目和位置，建立表示网格的数据结构。采用的单元形状和节点的分布，以及插值

5、函数的选取还应考虑到计算精度和可微性的要求。这里通过ANSYS ICEM CFD 14.0建模并划分网格。具体而言，网格将求解区域分为个281节点和565个单元，所有单元均为三角形单元，如图2所示实际上由于matlab计算编程是不知如何直接读取网格数据，就只选取了180个单元与103个节点进行计算。图2：左上角四分之一区域的流场及其有限单元划分流场网格3.单元分析单元分析的目的是建立有限元方程。把有限元积分表达式(3) 写为各个单元求和的形式e=1Need+2vsd=0（6）这里（e）表示单元e，Ne是单元总数，如果仅在一个单元上考虑上式，形式上有（e）d+2evsd=0（7）其中（e）表示单

6、元e的边界，上式实质上并不是一个等式，只具有形式上的意义，当对所有的单元求和以后，才是等式。如果把线积分中的2(e) 换为(e)，则得到的是等式，但在对所有单元求和时，内部边界的线积分刚好抵消，因此(7) 也可以理解为不计内部边界贡献的(3) 式在单元上的表达式。流函数在单元e内可用如下函数近似：=iNi (8)这里i (i = 1, 2, 3) 为节点流函数值，Ni为节点上的插值函数，上式中重复下标表示约定求和。将(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采用三节点三角形单元时，单元的插值基函数为Nix,y=aie+biex+ciey（11）如果单元e三个点坐标为（xie,yie）,i=1,2,3,则Nixje,yje=ij（12）即插值基函数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此即采用线性单元时的单元方程系数矩阵。其中A(e)为三角形（积分区域）的面积，bc的值可由（12）求得，现在列举如下： （13）A(e)=12x2-x1y3-y1-(y2-y1)(x3-x1)求解单元系数矩阵时，一般同时进行总体合成，每形成一个单元方程，便把它累加到总体方程中。出于顺序和逻辑上的考虑，下一步再详细说明总体合成的方法。对于边界积分项，我们假设三角形单元e中序号为，的节点在边界上，为自然边界，其长度为l。首先，注意到插值函数在边上是零。所以，可以得到如下结论： 图3：自然边界条件的处理右端项：。为了计算和，以点为原点，沿直

9、线建立局部坐标系，在此坐标系中，插值函数和如上图所示，可写成线性插值函数如下：假设切向速度在两节点处的值分别为和，并且沿边界是线性分布的，可以表示为。于是可以得到。对于前面讨论的圆柱绕流问题，由于，所以，根据线性解的性质，必有f（e）=f（e）=f（e）=0无需考虑f的影响，使程序得到了不少的简化。4.总体合成总体合成的过程就是把已经形成好的单元方程按一定顺序迭加起来，形成总体有限元方程。具体做法是根据单元内节点的总体顺序号，把单元方程进行延拓，未知量包含所有节点上的函数值，与此单元无关的位置以零填充，把所有的单元方程都进行延拓以后，进行系数矩阵的累加，便得到总体方程。理论上说，这一过程也可以

10、通过引入一个Boole 型矩阵来实现，定义单元e的boole矩阵,i=1,2,3; j=1,2,3Np.矩阵B其实就是单元节点序号表的又一表达形式。单元e的系数矩阵以及右端项沿拓后就是： 进而总体合成的过程可以表示为， 。但是这种方法比较麻烦，要重新定义新的矩阵，而且还要涉及计算矩阵转置和矩阵相乘等运算，一方面计算量较大，并且浪费空间，另一方面人为地增加了程序的复杂性，降低了程序的可读性。因此，这种方法一般只用作理论分析。实际的计算中，每当计算出一个单元系数矩阵Aij(e)，假设单元e三个节点编号分别为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项中。采用此方法并未多占用计算机内存，运算量也并不大（总共进行9*e次加法运算，不进行乘法运算）。本题的算法中选用的即此方法。（5）边界条件处理这里的边界条件是指本质边界条件，自然边界条件已经包含在积分表达式中了。具体做法是将已知的值代入到方程组中，并把已知的值移到方程组的右段，形成右端项。（6）求解有限元方程组并计算相关物理量有限元方程组的求解是一个代数问题，应针对具体的问题采用

12、合适的方法求解。对于对角占优的代数方程组，可以采用迭代法求解，规模不大的可以用Gauss 消元法一类的直接法求解，三对角方程则可以用追赶法。求出所有的待求量后，便得到了近似函数的表达式，并可以计算出相关的物理量。对计算结果进行综合的分析，以期得到原问题的正确的物理解答。对于每个单元，速度可以根据vx=y=iNiy=iNiy=cii(14)vy=-x=-iNix=-iNix=-bii(15)来计算,节点上的速度值可取这个节点相邻单元的速度值的平均。节点上的压力值可以有伯努利方程计算。假设求解区域位于同一水平面内，介质密度=1，来流压力p=0，那么p=121-vi2,i=1,210。如此便可得到节

13、点处的速度和压力分布。四具体算法实现本文具体算法是通过MATLAB实现的，matlab程序文件及算法说明见附件。五结果分析运行附件中MATLAB 程序，可以得到图4和图6两张图。图4显示1/4 区域的流场分布，图中的点以不同颜色表示各个结点处的流函数，图中的曲线为流函数的等势线，即是流线，由于对称性，整个区域的流场分布可明显看出，故不再画出。而圆柱绕流问题的解析解为，便于与计算结果对比。从图4可以发现，有限元法数值法得到的流线与解析法得到的流线在形态上基本一致。图4：有限元计算的1/4区域流场分布图5：解析法的整个区域流场分布具体到1/4 区域的右边界上的结点，将有限元法得到的流函数数值解与解

14、析解相对比，得到图7 中的曲线。可见，有限元法与解析法所得曲线在趋势上基本一致，但在数值上最大误差为11%。图6: 1/4区域右边界上流函数的有限元数值解与解析解对比附件：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); % 返回单元总数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); % 返回边界结点数%边界已知流函数的结点

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); % 返回边界已知流函数的结点数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); %未知流函数的结点数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 %求系数矩阵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 %扫描未知流函数的结点F(i)=0;for j=1:UknumANEW(i,j)=A(UKN(i),UKN(j); %新的系数矩阵endfor k=1:Bknum %扫描已知流函数的结点F(i)=F(i)+A(UKN(i),BKN(1,k)*BKN(2,k);endF(i)=-F(i); %移项到等号右侧成新的乘积系数向量end%FAIUN=ANEWF %解向量，未知结点的流函数值Q,R