有限元法解圆柱绕流问题

有限元法求解无限流体中的圆柱绕流问题2016年01月12号一.问题描述考虑位于两块无限长平板间的圆柱体的平面绕流问题，几何尺寸如下图所示，来流为J=l，Vy=0。由于流场具有上下左右的对称性，只考虑左上角四分之一的计算区域abcde,把它作为有限元的求解区域。°要求求解出整个区域中的流函数、以及压强值。在足够远前方选与来流垂直的控制面ae,cd是沿y轴，亦即一流动对称轴，be是物而，ab亦是流动对称轴，所要考虑的流动区域即由线abedea所围成的区域，在这一区域O中有：边界ab为流线，取虹0,慕二0；边界be也为流线，同样取卜0,穿=0；边界cd,切向速度1=0,务=0，取由=0：边界de为流线，满足L=L+J；v：＜dy=J*dy=2于是在ed上，巾=2,有二0；进口边界ae上，虹L+J：v/y=y（本文中采取此条件）也可以提自然边界条件云=0，亦■=跟二1我们以流函数1H乍为未知函数来解此问题，流函数所满足的微分方程如下：V?w=Q巾I、二〒(本质边界条件)鬃|「2=-VS(自然边界条件)(1)此处「1指ab,be,de和ae四段边界，而「2就是就是cd段边界，且切向速度、七二0,「1和「2合起来是整个边界，并且此二者不重合。下面，按有限元方法的一般步骤来计算此问题。有限元法解圆柱绕流问题建立有限元积分表达式根据求解问题的基本控制方程，应用变分法或加权余量法将求解的微分方程定解问题化为等价的枳分表达式，作为有限元法求解问题的出发方程式。对于方程(1),它是一椭圆型方程，具有正定性，可以用变分法，这里直接给出泛函J(W)二：JVp.VpdQ+Jv~d「=0(2)ar2令其变分8户0,可以得到J"."WdQ+JVs84>dr=0(3)Qr2自然边界条件已经包含在变分表达式中(其名称的由来)，而本质边界条件必须强制族满足(因此称其为本质边界条件，也称为强制边界条件)。如果根据原微分方程中无法给出泛函/，则可以用Galerkin加权余量方法得到枳分方程，这相当于将原来的微分方程写为如下变分形式：J△W5WdQ=0(4)Q这里的8W是函数巾的改变量，是一种“虚位移”，在本质边界条件口上为零。因此,上式做分部积分后，边界积分仅剩下.「2的部分。具体为JVp.VspdQ+Jvs8lbdr=0(5)Q「2即(3)式。可见，如果W满足原来的微分方程和边界条件，那么，必然有巾满足(4)式，进而满足(5)式。注意，在(5)式中，包含的边界「2上的边界条件信息，对边界的部分，仅知道它是给定了函数值的边界，却不知道边界上的值是多少，为了确定这些值，还需要额外的处理方法。正是因为「2上的边界信息可以包含在积分表达式中，这种边界条件也称为自然边界条件。区域剖分根据物理问题的特点以及区域的形状，把计算区域分成许多几何形状规则但大小可以不同的单元，确定单元节点的数目和位置，建立表示网格的数据结构。采用的单元形状和节点的分布，以及插值函数的选取还应考虑到计算精度和可微性的要求。这里通过ANSYSICEMCFD14.0建模并划分网格。具体而言，网格将求解区域分为个281节点和565个单元，所有单元均为三角形单元，如图2所示实际上由于matlab计算编程是不知如何直接读取网格数据，就只选取了180个单元与103个节点进行计算。图2：左上角四分之一区域的流场及其有限单元划分流场网格单元分析单元分析的目的是建立有限元方程。把有限元枳分表达式(3)写为各个单元求和的形式JV364»dr「20(6)JVp.V5JV364»dr「20(6)tQ(e>e=1w这里Q表示单元e,晚是单元总数，如果仅在一个单元上考虑上式，形式上有JVp・V6WdQ+Jv864>dr=0(7)q(c)r2nr(c)其中「〈°)表示单元e的边界，上式实质上并不是一个等式，只具有形式上的意义，当对所有的单元求和以后，才是等式。如果把线积分中的r2nr(e)换为「(•)，则得到的是等式，但在对所有单元求和时，内部边界的线积分刚好抵消，因此(7)也可以理解为不计内部边界贡献的(3)式在单元上的表达式。流函数巾在单元e内可用如下函数近似：W=WiNi(8)这里蛔(，=1,2,3)为节点流函数值，启为节点上的插值函数，上式中重复下标表示约定求和。将(8)代入(7),不难得到fONMNjON0Njrv^Nj8巾jd「(9)r2nr<e>J(、2v^Nj8巾jd「(9)r2nr<e>由于$蜻的任意性，所以，对于j=l,2,3都有

)加dQ=-JVsNjd「(10)r2)加dQ=-JVsNjd「(10)r2nr(e)此即单元方程，通常可以简写为A导=J%▽、・dQf3=-JvsNjdrQ（e）r2nr<e>采用三节点三角形单元时，单元的插值基函数为Nj（x,y）=+b（fx+c（?y（11）如果单元e三个点坐标为（x@,y（f）,i=1,2,3,则Ni（x（?,y（?）=%（12）即插值基函数Ni在x0点取1,在x（?x（?两点为零。由此不难解出abc。注意到求A曹时对Nj取了梯度，故的取值并不影响最后的计算。(A2123

bbb111Lb七21c(A2123

bbb111Lb七21cbib2+C1C2bib3+c】C3l)2+C2bgbs+C2C3bzb3+c2C3b?+此即采用线性单元时的单元方程系数矩阵。其中A①为三角形（积分区域）的面积，be的值可由（12）求得，现在列举如下：々=嘉（刈一为）6=嘉（）'3—）1）九=日耳（）1一为）1,、1z1zA（c）二:［（X2-X］）（y3-yi）-（V2-yi）（X3-Xi）］求解单元系数矩阵时，一般同时进行总体合成，每形成一个单元方程，便把它累加到总体方程中。出于顺序和逻辑上的考虑，下一步再详细说明总体合成的方法。对于边界枳分项，我们假设三角形单元e中序号为的节点在边界匚上，砂为自然边界，其长度为1。首先，注意到插值函数在另边上是零。所以，可以得到如下结论:

图3：自然边界条件的处理右端项：f^e)=0O为了计算C和4气以。点为原点，沿直线两建立局部坐标系g,在此坐标系中，插值函数虬和如上图所示，可写成线性插值函数如下：虬=1-§,七==假设切向速度气在两节点处的值分别为《白和V/,并且沿边界是线性分布的，可以表于是可以得到fae)=-!Ed&=+(*.f)j)Q-吊)此=-{(2*“+由ffi)=-fg处=-J(yM+(y匕苫)沁=4(ysa+。对于前面讨论的圆柱绕流问题，由于匕“=匕月=0,所以，根据线性解的性质，必有f%)=f曾二f号二0无需考虑f的影响，使程序得到了不少的简化。总体合成总体合成的过程就是把己经形成好的单元方程按一定顺序迭加起来，形成总体有限元方程。具体做法是根据单元内节点的总体顺序号，把单元方程进行延拓，未知量包含所有节点上的函数值，与此单元无关的位置以零填充，把所有的单元方程都进行延拓以后，进行系数矩阵的累加，便得到总体方程。理论上说，这一过程也可以通过引入一个Boole型矩阵来实现，定义单元e的boole矩阵B*NpD1,如果单元e的地i个点总体顺序号是j0.其他情况'S3村'2,3・％矩阵B其实就是单元节点序号表的又一表达形式。单元e的系数矩阵A⑴以及右端项/⑴沿拓后就是：技)=B7fe)虬f=w。e=l亦=BTA<c)技)=B7fe)虬f=w。e=l进而总体合成的过程可以表示为虬C=1但是这种方法比较麻烦，要重新定义新的矩阵B,而且还要涉及计算矩阵转置和矩阵相乘等运算，一方面计算量较大，并且浪费空间，另一方面人为地增加了程序的复杂性，降低了程序的可读性。因此，这种方法一般只用作理论分析。实际的计算中，每当计算出一个单元系数矩阵A肾，假设单元e三个节点编号分别为i,j,k,那么将A曹中的1*1项放入大矩阵（借鉴结构力学的概念，不妨称其为总体“刚度”矩阵，下同）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）求解有限元方程组并计算相关物理量有限元方程组的求解是一个代数问题，应针对具体的问题采用合适的方法求解。对于对角占优的代数方程组，可以采用迭代法求解，规模不大的可以用Gauss消元法一类的直接法求解，三对角方程则可■以用追赶法。求出所有的待求量后，便得到了近似函数的表达式，并可以计算出相关的物理量。对计算结果进行综合的分析，以期得到原问题的正确的物理解答。对于每个单元，速度可以根据巾。WiNi洲i=矿3T=虬商=cm（M）ag*aNi'’=一云二一节厂=一巾顽=-加虬⑵）来计算，节点上的速度值可取这个节点相邻单元的速度值的平均。节点上的压力值可以有伯努利方程计算。假设求解区域位于同一水平面内，介质密度P=1,来流压力p=o,那么p=；（1-v初，i=l,2-10o如此便可得到节点处的速度和压力分布。具体算法实现本文具体算法是通过MATLAB实现的，matlab程序文件及算法说明见附件。结果分析运行附件中MATLAB程序，可以得到图4和图6两张图。图4显示1/4区域的流场分布，图中的点以不同颜色表示各个结点处的流函数，图中的曲线为流函数的等势线，即是流线，由于对称性，整个区域的流场分布可明显看出，故不再画出。而圆柱绕流问题的解析解为U=y（l-c'Q，便于与计算结果对比。从图4可以发现，有限元法数值法得到的流线与解析法得到的流线在形态上基本一致。

图4：有限元计算的1/4区域流场分布图5：解析法的整个区域流场分布具体到1/4区域的右边界上的结点，将有限元法得到的流函数数值解与解析解相对比，得到图7中的曲线。可见，有限元法与解析法所得曲线在趋势上基本一致，但在数值上最大误差为11知0.51121.4161.822.22.42.6283525210右边界y坐标图6:1/4区域右边界上流函数的有限元数值解与解析解对比附件：LNATLAB主程序:youxianyuan.mENA=[4132370.0397983014144320.047529903101210.1201568171031010.094963772103221020.01349843710318220.0912433183036280.0380546623035360.0422875212120.038074502620.0487429047170.05029162870.03714096488760.0579940111920880.099591586101102990.0163826141031020.0442188919450.0471326394441910.0559674063552380.0360855133516520.034571557294460.0442808632932440.04144814520210.10610919213930.0782126135940.051017838100990.03408803910198210.0503592161001020.01845156599981010.04518616797210.07031062996220.097022072100220.062253661951000.032213795991000.03169101198990.03919582997980.03840987389200.06285272848730.0504978972410.09490433296230.060406457

969280959591919090898988831374919487829381248685858484929271809179907889778883141475748787828281248681867385728484718079797878777776698383757482688181737372727170796678657764761000.0364728330.0320297480.0296740240.041864912970.0455070830.06124771190.093342661410.0515515640.03920888330.0582770660.063798865230.047832817230.050837396960.038920434800.041373741950.0335594080.0421603950.0448907660.0414661670.060586214150.0511069380.04157820630.057195134930.0455331390.0845552130.0562416820.0301726850.0328142850.0372638560.0472626730.043683491890.043858618880.033382546190.0531674810.0303011510.039566275820.0453036580.039028656850.03280082840.036391036800.050109968790.043207084780.040783626770.03638318

76636975755857747468687362726171716767705970706666656564646363696958576868626261616767565966556554645363635857625161615656605060495959555554545353585848435757515156564750494955190.060158755830.0259775240.055241265410.053616614820.046847336810.0417935540.0360754510.035939511800.038472382800.040185579560.045615694790.04312012780.041221814770.038815812760.035503832190.057736703750.0324044420.042694580.039041739720.037749165710.033244944700.038889688700.0408180030.0405851830.040088130.042912099690.036059255680.0415915650.0396859760.0376339090.0362883440.0318601790.0315214340.0391193630.03989979640.0409038130.046404299150.066170074410.04763151620.041956550.046636447600.035592793600.031796449590.039229272

464553485243513342494645484347334239454443344239413740343637343632353433253130272629252854550.04097428853540.04180207580.04219602652150.056305234150.04715458551570.041090301560.04582129350400.048398755500.040513663550.040202562540.04207801530.0448719520.037263593510.04153873540500.03986560242500.050325014490.04102327145460.04246825480.04395136560.05101737440470.04178759542330.051746427460.042853263450.041743856430.038609599430.0318392331330.050574027420.04382332390.04221094131400.03318516836390.04169815235380.039194122370.034698998160.04919161528360.040629125340.049558416320.039551559330.0481510817350.041250709340.039846508330.04614112320.037979443310.038580787300.036116702

9280.0334817012610270.0327575351225290.0354277362511260.033372838280.0319017119270.03044574710260.03028529711250.031849171];%单元与结点号、面积对应关系矩阵NC0RD=[-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.489138169-0.9708951760.74446511-0.7573201670.962580378-0.505508741.128325608-1.2621251320.243714522-0.2514580981.253891963-1.2771060360.738778719-1.447512260.48199085-1.0881676511.056783565-0.775236311.2672665480.2459580751.5851798260.5038471451.4287301131.5487045030.7367529960.4892388611.7438800030.7646746611.5808443781.4507575310.9818528671.7967201420.57457131.0498231171.4132953951.7651239080.9065805151.6189897260.2552368690.7496933581.8928234241.0299777561.725524321.6491440061.2002037440.4493600412.058572371.3030399261.5637018731.347538531.2701449241.9546564711.1430551590.2357729261.8751930290.740660042.196622761.0207719462.0312714251.2938474281.8636096171.8253854811.4671998412.0584840570.8390018060.4554460192.3511649711.5606991231.6926493411.538865611.4370474512.150601831.3732556742.250340161.0853583850.7482476712.5179112731.0109146942.3291435741.2856165682.1608701351.5635246021.9862794072.0663439811.6290341752.3546506260.7857298640.5091206722.6280369071.8397519731.7996403542.3346089411.603795932.4278381391.3299698922.522081281.0533585992.1604472720.5323468770.2555346892.5273912840.9993942782.6077561911.2668873772.4549324771.5548922812.288393071.846091022.1075668482.1693605321.9061661852.6454923870.75139032.4747812040.4600113480.3029310722.7510862382.6103336771.5746361822.6747436711.3013713342.7745191331.0681789212.2826019690.2524894831.2145620962.7344224381.5356951282.6040624111.8315090082.4166471132.1129514952.233713822.4682770141.8599571482.1693605321.9061661852.6454923870.75139032.4747812040.4600113480.3029310722.7510862382.6103336771.5746361822.6747436711.3013713342.7745191331.0681789212.2826019690.2524894831.2145620962.7344224381.5356951282.6040624111.8315090082.4166471132.1129514952.233713822.4682770141.8599571482.7469217810.3910630661.9777813670.270089192.3419618552.0937896792.7412744741.8595979221.8429976772.7170421542.0921372472.5447451912.3596635212.3421241072.5995567142.1269871342.3601280752.6795818212.6346098662.3797448182.57733166-2.2.57733166];%结点与坐标对应关系矩阵[Enum,temp]=size(ENA);%返回单元总数[Nnum,temp]=size(NCORD);%返回结点总数BN0DE=[l34562121110987171615141319202118222324];%边界结点编号[temp,Bnum]=size(BNODE);%返回边界结点数%边界己知流函数的结点号及其值BKN=L1345621211109871319202118222324;000000000000333332.251.50.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,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);戋未知流函数的结点数fork=l:Enum%单元号xl=NC0RD(ENA(k,1),1);yl=NCORD(ENA(k,1),2);x2=NC0RD(ENA(k,2),1)；y2=NC0RD(ENA(k,2),2);x3=NC0RD(ENA(k,3),1)；y3=NC0RD(ENA(k,3),2);a(l)=x2*y3-x3*y2;b(l)=y2~y3;c(l)=x3-x2;a(2)=x3*yl"xl*y3;b(2)=y3~yl;c(2)=xl-x3;a(3)=xl*y2"x2*yl;b(3)=yl~y2;c(3)=x2~xl;forn=l:3form=l:3ANM(k,n,m)=l/(4*ENA(k,4))*(b(n)*b(m)+c(n)*c(m));%各个单元的Anm信息endendendfori=l:Nnum%求系数矩阵Aforj=l:Nnumtemp=0;fore=l:Enumforn=l:3form=l:3ifENA(e,n)==i&&ENA(e,m)==jtemp=temp+A?<M(e,n,m);endendendendA(i,j)=temp;endendfori=l：Uknum%扫描未知流函数的结点F(i)=0;forj=l:UknumANEW(i,j)=A(UKN(i),UKN(j));%新的系数矩阵endfork=l：Bknum%扫描己知流函数的结点F(i)=F(i)+A(UKN(i),BKN(1,k))*BKN(2,k);endF(i)=-F(i);%移项到等号右侧成新的乘积系数向量end%FAIUN=ANEW\F，%解向量，未知结点的流函数值[Q,R]=lu(ANEW);FAIUN=R\(Q\F');先将所有结点的流函数值写入FAIN向量fori=l:UknumFAIN