您好,欢迎访问三七文档
当前位置:首页 > 建筑/环境 > 工程监理 > 有限元法解圆柱绕流问题
有限元法求解无限流体中的圆柱绕流问题2016年01月12号一.问题描述考虑位于两块无限长平板间的圆柱体的平面绕流问题,几何尺寸如下图所示,来流为𝑣𝑥=1,𝑣𝑦=0。由于流场具有上下左右的对称性,只考虑左上角四分之一的计算区域abcde,把它作为有限元的求解区域Ω。要求求解出整个区域中的流函数、𝑣𝑥、𝑣𝑦以及压强值。图1:圆柱绕流模型二.数学建模在足够远前方选与来流垂直的控制面ae,cd是沿y轴,亦即一流动对称轴,bc是物面,ab亦是流动对称轴,所要考虑的流动区域即由线abcdea所围成的区域,在这一区域中有:1.边界ab为流线,取ψ=0,𝜕𝜕𝑛=0;2.边界bc也为流线,同样取ψ=0,𝜕𝜕𝑛=0;3.边界cd,切向速度𝑣𝜏=0,𝜕𝜓𝜕𝑛=0,取=0;4.边界de为流线,满足=∫𝑣𝑥𝑑𝑦=∫𝑑𝑦=于是在ed上,ψ=2,𝜕𝜕𝑛=0;5.进口边界ae上,ψ=∫𝑣𝑥𝑑𝑦𝑦=𝑦(本文中采取此条件)也可以提自然边界条件𝜕𝜓𝜕𝑛=0,𝜕𝜕𝑛=𝑣𝑥=1我们以流函数ψ作为未知函数来解此问题,流函数所满足的微分方程如下:{=0=̅(本质边界条件)𝜕𝜓𝜕𝑛=𝑣𝑠(自然边界条件)(1)此处指,,和四段边界,而就是就是cd段边界,且切向速度𝑣𝑠=0,Γ1和Γ2合起来是整个边界,并且此二者不重合。下面,按有限元方法的一般步骤来计算此问题。三.有限元法解圆柱绕流问题1.建立有限元积分表达式根据求解问题的基本控制方程,应用变分法或加权余量法将求解的微分方程定解问题化为等价的积分表达式,作为有限元法求解问题的出发方程式。对于方程(1),它是一椭圆型方程,具有正定性,可以用变分法,这里直接给出泛函(ψ)=1∫∫𝑣𝑠=0()令其变分δJ=0,可以得到∫∫𝑣𝑠=0()自然边界条件已经包含在变分表达式中(其名称的由来),而本质边界条件必须强制ψ满足(因此称其为本质边界条件,也称为强制边界条件)。如果根据原微分方程中无法给出泛函J,则可以用Galerkin加权余量方法得到积分方程,这相当于将原来的微分方程写为如下变分形式:∫=0()这里的δψ是函数ψ的改变量,是一种“虚位移”,在本质边界条件上为零。因此,上式做分部积分后,边界积分仅剩下的部分。具体为∫∫𝑣𝑠=0()即(3)式。可见,如果ψ满足原来的微分方程和边界条件,那么,必然有ψ满足(4)式,进而满足(5)式。注意,在(5)式中,包含的边界Γ2上的边界条件信息,对边界Γ1的部分,仅知道它是给定了函数值的边界,却不知道边界上的值是多少,为了确定这些值,还需要额外的处理方法。正是因为Γ2上的边界信息可以包含在积分表达式中,这种边界条件也称为自然边界条件。2.区域剖分根据物理问题的特点以及区域的形状,把计算区域分成许多几何形状规则但大小可以不同的单元,确定单元节点的数目和位置,建立表示网格的数据结构。采用的单元形状和节点的分布,以及插值函数的选取还应考虑到计算精度和可微性的要求。这里通过ANSYSICEMCFD14.0建模并划分网格。具体而言,网格将求解区域分为个281节点和565个单元,所有单元均为三角形单元,如图2所示实际上由于matlab计算编程是不知如何直接读取网格数据,就只选取了180个单元与103个节点进行计算。图2:左上角四分之一区域的流场及其有限单元划分流场网格3.单元分析单元分析的目的是建立有限元方程。把有限元积分表达式(3)写为各个单元求和的形式∑∫()∫𝑣𝑠=0()这里()表示单元e,是单元总数,如果仅在一个单元上考虑上式,形式上有∫()∫𝑣𝑠=0()()其中()表示单元e的边界,上式实质上并不是一个等式,只具有形式上的意义,当对所有的单元求和以后,才是等式。如果把线积分中的Γ2∩Γ(e)换为Γ(e),则得到的是等式,但在对所有单元求和时,内部边界的线积分刚好抵消,因此(7)也可以理解为不计内部边界贡献的(3)式在单元上的表达式。流函数ψ在单元e内可用如下函数近似:=()这里(i=1,2,3)为节点流函数值,为节点上的插值函数,上式中重复下标表示约定求和。将(8)代入(7),不难得到∫(()𝑦𝑦)𝑑=∫𝑣𝑠()()由于的任意性,所以,对于j=1,2,3都有∫(()𝑦𝑦)𝑑=∫𝑣𝑠()(10)此即单元方程,通常可以简写为()=∫()𝑑=∫𝑣𝑠()采用三节点三角形单元时,单元的插值基函数为(,)=()()()𝑦(11)如果单元e三个点坐标为((),𝑦()),i=1,2,3,则((),𝑦())=(1)即插值基函数在()点取1,在()()两点为零。由此不难解出abc。注意到求()时对取了梯度,故()的取值并不影响最后的计算。对某一固定的单元e,将(11)式代入(10)中,可以得到:()=()[]此即采用线性单元时的单元方程系数矩阵。其中()为三角形(积分区域)的面积,bc的值可由(12)求得,现在列举如下:123()1()2ebyyA231()1()2ebyyA312()1()2ebyyA132()1()2ecxxA213()1()2ecxxA321()1()2ecxxA(13)()=1()(𝑦𝑦)(𝑦𝑦)()求解单元系数矩阵时,一般同时进行总体合成,每形成一个单元方程,便把它累加到总体方程中。出于顺序和逻辑上的考虑,下一步再详细说明总体合成的方法。对于边界积分项,我们假设三角形单元e中序号为,,的节点在边界2上,为自然边界,其长度为l。首先,注意到插值函数在边上是零。所以,可以得到如下结论:图3:自然边界条件的处理右端项:()0ef。为了计算()ef和()ef,以点为原点,沿直线建立局部坐标系,在此坐标系中,插值函数N和N如上图所示,可写成线性插值函数如下:NNll假设切向速度sV在两节点处的值分别为sV和sV,并且沿边界是线性分布的,可以表示为()ssassavvvvl。于是可以得到()00(())(1)(2)6llessassasaslfvNdvvvdvvll()00(())(2)6llessassasaslfvNdvvvdvvll。对于前面讨论的圆柱绕流问题,由于0sasvv,所以,根据线性解的性质,必有α()=β()=γ()=0无需考虑f的影响,使程序得到了不少的简化。4.总体合成总体合成的过程就是把已经形成好的单元方程按一定顺序迭加起来,形成总体有限元方程。具体做法是根据单元内节点的总体顺序号,把单元方程进行延拓,未知量包含所有节点上的函数值,与此单元无关的位置以零填充,把所有的单元方程都进行延拓以后,进行系数矩阵的累加,便得到总体方程。理论上说,这一过程也可以通过引入一个Boole型矩阵来实现,定义单元e的boole矩阵3pNB31,0.pNB如果单元e的地i个点总体顺序号是j其他情况,i=1,2,3;j=1,2,3…Np.矩阵B其实就是单元节点序号表的又一表达形式。单元e的系数矩阵()eA以及右端项()ef沿拓后就是:()(),eTeABAB()()eTefBf进而总体合成的过程可以表示为()1eNeeAA,()1eNeeff。但是这种方法比较麻烦,要重新定义新的矩阵B,而且还要涉及计算矩阵转置和矩阵相乘等运算,一方面计算量较大,并且浪费空间,另一方面人为地增加了程序的复杂性,降低了程序的可读性。因此,这种方法一般只用作理论分析。实际的计算中,每当计算出一个单元系数矩阵(),假设单元e三个节点编号分别为i,j,k,那么将()中的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消元法一类的直接法求解,三对角方程则可以用追赶法。求出所有的待求量后,便得到了近似函数的表达式,并可以计算出相关的物理量。对计算结果进行综合的分析,以期得到原问题的正确的物理解答。对于每个单元,速度可以根据𝑥=𝑦=𝑦=𝑦=(1)𝑦====(1)来计算,节点上的速度值可取这个节点相邻单元的速度值的平均。节点上的压力值可以有伯努利方程计算。假设求解区域位于同一水平面内,介质密度=1,来流压力p=0,那么=(1𝑣),=1,10。如此便可得到节点处的速度和压力分布。四.具体算法实现本文具体算法是通过MATLAB实现的,matlab程序文件及算法说明见附件。五.结果分析运行附件中MATLAB程序,可以得到图4和图6两张图。图4显示1/4区域的流场分布,图中的点以不同颜色表示各个结点处的流函数,图中的曲线为流函数的等势线,即是流线,由于对称性,整个区域的流场分布可明显看出,故不再画出。而圆柱绕流问题的解析解为)11(y22yx,便于与计算结果对比。从图4可以发现,有限元法数值法得到的流线与解析法得到的流线在形态上基本一致。图4:有限元计算的1/4区域流场分布图5:解析法的整个区域流场分布具体到1/4区域的右边界上的结点,将有限元法得到的流函数数值解与解析解相对比,得到图7中的曲线。可见,有限元法与解析法所得曲线在趋势上基本一致,但在数值上最大误差为11%。图6:1/4区域右边界上流函数的有限元数值解与解析解对比附件:1.NATLAB主程序:youxianyuan.mENA=[4132370.0397983014144320.04752990318101210.120156817181031010.094963772103221020.04349843710318220.0942433183036280.0380546623035360.042287521292120.03807450229620.048742904307170.0502916230870.0371409641988760.0579940411920880.099591586101102990.0463826441011031020.044218891449450.0471326394441940.0559674063552380.0360855133516520.034571557294460.0442808632932440.0414481459720210.10610919213930.07821261345940.054017838102100990.03408803910198210.050359216221001020.04845156599981010.0451864679897210.0703106292396220.09702207296100220.06225366192951000.032213795991000.0316910119198990.0394958299097980.0384098739789200.06285272848730.050497897932410.0949043328496230.06040645796921000.036472833809
本文标题:有限元法解圆柱绕流问题
链接地址:https://www.777doc.com/doc-2321330 .html