有限元法解圆柱绕流问题.doc

,几何尺寸如下图所示,来流为。由于流场具有上下左右的对称性,只考虑左上角四分之一的计算区域abcde,把它作为有限元的求解区域Ω。要求求解出整个区域中的流函数、以及压强值。图1:,cd是沿y轴,亦即一流动对称轴,bc是物面,ab亦是流动对称轴,所要考虑的流动区域即由线abcdea所围成的区域,在这一区域中有:,取ψ=0,,同样取ψ=0,,切向速度=0,,满足于是在ed上,ψ=2,,ψ=(本文中采取此条件)也可以提自然边界条件我们以流函数ψ作为未知函数来解此问题,流函数所满足的微分方程如下:(1)此处就是就是cd段边界,且切向速度=0,Γ1和Γ2合起来是整个边界,并且此二者不重合。下面,按有限元方法的一般步骤来计算此问题。,应用变分法或加权余量法将求解的微分方程定解问题化为等价的积分表达式,作为有限元法求解问题的出发方程式。对于方程(1),它是一椭圆型方程,具有正定性,可以用变分法,这里直接给出泛函令其变分δJ=0,可以得到自然边界条件已经包含在变分表达式中(其名称的由来),而本质边界条件必须强制ψ满足(因此称其为本质边界条件,也称为强制边界条件)。如果根据原微分方程中无法给出泛函J,则可以用Galerkin加权余量方法得到积分方程,这相当于将原来的微分方程写为如下变分形式:这里的δψ是函数ψ的改变量,是一种“虚位移”,在本质边界条件。因此,上式做分部积分后,边界积分仅剩下。具体为即(3)式。可见,如果ψ满足原来的微分方程和边界条件,那么,必然有ψ满足(4)式,进而满足(5)式。注意,在(5)式中,包含的边界Γ2上的边界条件信息,对边界Γ1的部分,仅知道它是给定了函数值的边界,却不知道边界上的值是多少,为了确定这些值,还需要额外的处理方法。正是因为Γ2上的边界信息可以包含在积分表达式中,这种边界条件也称为自然边界条件。,把计算区域分成许多几何形状规则但大小可以不同的单元,确定单元节点的数目和位置,建立表示网格的数据结构。采用的单元形状和节点的分布,以及插值函数的选取还应考虑到计算精度和可微性的要求。。具体而言,网格将求解区域分为个281节点和565个单元,所有单元均为三角形单元,如图2所示实际上由于matlab计算编程是不知如何直接读取网格数据,就只选取了180个单元与103个节点进行计算。图2:。把有限元积分表达式(3)写为各个单元求和的形式这里表示单元e,是单元总数,如果仅在一个单元上考虑上式,形式上有其中表示单元e的边界,上式实质上并不是一个等式,只具有形式上的意义,当对所有的单元求和以后,才是等式。如果把线积分中的Γ2∩Γ(e)换为Γ(e),则得到的是等式,但在对所有单元求和时,内部边界的线积分刚好抵消,因此(7)也可以理解为不计内部边界贡献的(3)式在单元上的表达式。流函数ψ在单元e内可用如下函数近似:这里(i=1,2,3)为节点流函数值,为节点上的插值函数,上式中重复下标表示约定求和。将(8)代入(7),不难得到由于的任意性,所以,对于j=1,2,3都有此即单元方程,通常可以简写为采用三节点三角形单元时,单元的插值基函数为如果单元e三个点坐标为(),i=1,2,3,则即插值基函数在点取1,在两点为零。由此不难解出abc。注意到求时对的取值并不影响最后的计算。对某一固定的单元e,将(11)式代入(10)中,可以得到:此即采用线性单元时的单元方程系数矩阵。其中为三角形(积分区域)的面积,bc的值可由(12)求得,现在列举如下:(13)求解单元系数矩阵时,一般同时进行总体合成,每形成一个单元方程,便把它累加到总体方程中。出于顺序和逻辑上的考虑,下一步再详细说明总体合成的方法。对于边界积分项,我们假设三角形单元e中序号为,的节点在边界上,为自然边界,其长度为l。首先,注意到插值函数在边上是零。所以,可以得到如下结论:图3:自然边界条件的处理右端项:。为了计算和,以点为原点,沿直线建立局部坐标系,在此坐标系中,插值函数和如上图所示,可写成线性插值函数如下:假设切向速度在两节点处的值分别为和,并且沿边界是线性分布的,可以表示为。于是可以得到。对于前面讨论的圆柱绕流问题,由于,所以,根据线性解的性质,必有无需考虑