有限元法解圆柱绕流问题探索.docx

有限元法求解无限流体中的圆柱绕流问题 2016 年 01月 12号一. 问题描述考虑位于两块无限长平板间的圆柱体的平面绕流问题,几何尺寸如下图所示,来流为。由于流场具有上下左右的对称性,只考虑左上角四分之一的计算区域 abcde , 把它作为有限元的求解区域Ω。要求求解出整个区域中的流函数、以及压强值。图1 : ae,cd 是沿 y轴, 亦即一流动对称轴, bc 是物面, ab 亦是流动对称轴,所要考虑的流动区域即由线 abcdea 所围成的区域,在这一区域?中有: 1. 边界 ab 为流线,取ψ=0, 2. 边界 bc 也为流线,同样取ψ=0, 3. 边界 cd ,切向速度=0 ,取 4. 边界 de 为流线, 满足于是在 ed 上, ψ=2, 5. 进口边界 ae 上, ψ= (本文中采取此条件) 也可以提自然边界条件我们以流函数ψ作为未知函数来解此问题,流函数所满足的微分方程如下: (1) 此处就是就是 cd 段边界,且切向速度=0,Γ 1和Γ 2 合起来是整个边界,并且此二者不重合。下面,按有限元方法的一般步骤来计算此问题。 1 .建立有限元积分表达式根据求解问题的基本控制方程, 应用变分法或加权余量法将求解的微分方程定解问题化为等价的积分表达式,作为有限元法求解问题的出发方程式。对于方程(1) ,它是一椭圆型方程,具有正定性,可以用变分法,这里直接给出泛函令其变分δ J=0 ,可以得到自然边界条件已经包含在变分表达式中( 其名称的由来) ,而本质边界条件必须强制ψ满足( 因此称其为本质边界条件,也称为强制边界条件)。如果根据原微分方程中无法给出泛函 J ,则可以用 Galerkin 加权余量方法得到积分方程,这相当于将原来的微分方程写为如下变分形式: 这里的δψ是函数ψ的改变量, 是一种“虚位移”, 在本质边界条件。因此, 上式做分部积分后,边界积分仅剩下。具体为即( 3 )式。可见, 如果ψ满足原来的微分方程和边界条件, 那么, 必然有ψ满足(4) 式, 进而满足(5) 式。注意,在(5) 式中,包含的边界Γ 2 上的边界条件信息,对边界Γ 1 的部分,仅知道它是给定了函数值的边界, 却不知道边界上的值是多少, 为了确定这些值, 还需要额外的处理方法。正是因为Γ 2 上的边界信息可以包含在积分表达式中, 这种边界条件也称为自然边界条件。 2 .区域剖分根据物理问题的特点以及区域的形状, 把计算区域分成许多几何形状规则但大小可以不同的单元, 确定单元节点的数目和位置, 建立表示网格的数据结构。采用的单元形状和节点的分布, 以及插值函数的选取还应考虑到计算精度和可微性的要求。这里通过 ANSYS ICEM CFD 建模并划分网格。具体而言, 网格将求解区域分为个 281 节点和 565 个单元, 所有单元均为三角形单元,如图 2 所示实际上由于 matlab 计算编程是不知如何直接读取网格数据,就只选取了 180 个单元与 103 个节点进行计算。图2 :左上角四分之一区域的流场及其有限单元划分流场网格 3. 单元分析单元分析的目的是建立有限元方程。把有限元积分表达式(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 )求得,现在列举如下: 1 2 3 ( )1 ( ) 2 e b y y A ? ? 2 3 1 ( )1 ( ) 2 e b y y A ?