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

考虑位于两块无限长平板间的圆柱体的平面绕流问题,几何尺寸如下图所示,来流为vx=1,vy=0。由于流场具有上下左右的对称性,只考虑左上角四分之一的计算区域abcde,把它作为有限元的求解区域Ω。要求求解出整个区域中的流函数、vx、vy以及压强值。图1: 在足够远前方选与来流垂直的控制面ae,cd是沿y轴,亦即一流动对称轴,bc是物面,ab亦是流动对称轴,所要考虑的流动区域即由线abcdea所围成的区域,在这一区域中有:,取ψ=0,∂φ∂n=0;,同样取ψ=0,∂φ∂n=0;,切向速度vτ=0,∂ψ∂n=0,取φ=0;,满足ψe=ψa+aevxdy=02dy=2于是在ed上,ψ=2,∂φ∂n=0;,ψ=ψ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),它是一椭圆型方程,具有正定性,可以用变分法,这里直接给出泛函J(ψ)=12Ω∇ψ∙∇ψdΩ+Γ2vsψdΓ=0(2)令其变分δJ=0,可以得到Ω∇ψ∙∇δψdΩ+Γ2vsδψdΓ=0(3)自然边界条件已经包含在变分表达式中(其名称的由来),而本质边界条件必须强制ψ满足(因此称其为本质边界条件,也称为强制边界条件)。如果根据原微分方程中无法给出泛函J,则可以用Galerkin加权余量方法得到积分方程,这相当于将原来的微分方程写为如下变分形式:Ω△ψδψdΩ=0(4)这里的δψ是函数ψ的改变量,是一种“虚位移”,在本质边界条件Γ1上为零。因此,上式做分部积分后,边界积分仅剩下Γ2的部分。具体为Ω∇ψ∙∇δψdΩ+Γ2vsδψdΓ=0(5)即(3)式。可见,如果ψ满足原来的微分方程和边界条件,那么,必然有ψ满足(4)式,进而满足(5)式。注意,在(5)式中,包含的边界Γ2上的边界条件信息,对边界Γ1的部分,仅知道它是给定了函数值的边界,却不知道边界上的值是多少,为了确定这些值,还需要额外的处理方法。正是因为Γ2上的边界信息可以包含在积分表达式中,这种边界条件也称为自然边界条件。,把计算区域分成许多几何形状规则但大小可以不同的单元,确定单元节点的数目和位置,建立表示网格的数据结构。采用的单元形状和节点的分布,以及插值函数的选取还应考虑到计算精度和可微性的要求。。具体而言,网格将求解区域分为个281节点和565个单元,所有单元均为三角形单元,如图2所示实际上由于matlab计算编程是不知如何直接读取网格数据,就只选取了180个单元与103个节点进行计算。图2:。把有限元积分表达式(3)写为各个单元求和的形式e=1NeΩe∇ψ∙∇δψdΩ+Γ2vsδψdΓ=0(6)这里Ω(e)表示单元e,Ne是单元总数,如果仅在一个单元上考虑上式,形式上有Ω(e)∇ψ∙∇δψdΩ+Γ2∩ΓevsδψdΓ=0(7)其中Γ(e)表示单元e的边界,上式实质上并不是一个等式,只具有形式上的意义,当对所有的单元求和以后,才是等式。如果把线积分中的Γ2∩Γ(e)换为Γ(e),则得到的是等式,但在对所有单元求和时,内部边界的线积分刚好抵消,因此(7)也可以理解为不计内部边界贡献的(3)式在单元上的表达式。 流函数ψ在单元e内可用如下函数近似:ψ=ψiNi(8)这里ψi(i=1,2,3)为节点流函数值,Ni为节点上的插值函数,上式中重复下标表示约定求和。将(8)代入(7),不难得到Ωe(∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y)ψiδψjdΩ=-Γ2∩ΓevsNjδψjdΓ(9)由于δψj的任意性,所以,对于j=1,2,3都有Ωe(∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y)ψidΩ=-Γ2∩ΓevsNjdΓ(10)此即单元方程,通常可以简写为Aij(e)=Ωe∇Ni∇NjdΩfj=-Γ2∩ΓevsNjdΓ