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

有限元法求解无限流体中的圆柱绕流问题
2016年01月12号
问题描述
,来流为
abcde ,
考虑位于两块无限长平板间的圆柱体的平面绕流问题 ，几何尺寸如下图所示
。由于流场具有上下左右的对称性，只考虑左上角四分之一的计算区域
把它作为有限元的求解区域 Qo要求求解出整个区域中的流函数、、以及压强值
数学建模
在足够远前方选与来流垂直的控制面 ae,cd是沿y轴，亦即一流动对称轴，be是物面，
ab亦是流动对称轴，所要考虑的流动区域即由线 abcdea所围成的区域，在这一区域■■中
有：
1•边界ab为流线，取2=0,
2•边界be也为流线,同样取2=0 ,
3•边界cd ,切向速度=0 , ，取 ；
4•边界de为流线，满足
于是在ed上 ,2=2 ,
，2= （本文中采取此条件）
也可以提自然边界条件 ，
我们以流函数2作为未知函数来解此问题，流函数所满足的微分方程如下 ：
2 （本质边界条件）（1）
自然边界条件
此处 指 和 四段边界，而 就是就是cd段边界，且切向速度=0 , 1和2合起
来是整个边界，并且此二者不重合。下面，按有限元方法的一般步骤来计算此问题 。
有限元法解圆柱绕流问题
1 .建立有限元积分表达式
根据求解问题的基本控制方程，应用变分法或加权余量法将求解的微分方程定解问题
化为等价的积分表达式，作为有限元法求解问题的出发方程式 。
对于方程（1）,它是一椭圆型方程，具有正定性，可以用变分法，这里直接给出泛函
（2） - 2 2 Q （）
令其变分J=0,可以得到
2 §2 Q §2
自然边界条件已经包含在变分表达式中 （其名称的由来），而本质边界条件必须强制 2
满足（因此称其为本质边界条件，也称为强制边界条件）。
如果根据原微分方程中无法给出泛函 J,则可以用Galerkin加权余量方法得到积分方程
这相当于将原来的微分方程写为如下变分形式 ：
2§2 Q
这里的§2是函数2的改变量，是一种虚位移”，在本质边界条件 上为零。因此，上
式做分部积分后，边界积分仅剩下 的部分。具体为
2 §2 Q §2
即（3）式。
可见，如果“满足原来的微分方程和边界条件 ，那么，必然有“满足(4)式，进而满足 (5)式。注意，在(5)式中，包含的边界2上的边界条件信息，对边界1的部分，仅知道 它是给定了函数值的边界，却不知道边界上的值是多少，为了确定这些值，还需要额外的 处理方法。正是因为2上的边界信息可以包含在积分表达式中 ，这种边界条件也称为自然 边界条件。

根据物理问题的特点以及区域的形状 ，把计算区域分成许多几何形状规则但大小可以
不
同的单元，确定单元节点的数目和位置 ，建立表示网格的数据结构 。采用的单元形状和节
占
八、、
的分布，以及插值函数的选取还应考虑到计算精度和可微性的要求 。这里通过ANSYS
ICEM CFD 。具体而言，网格将求解区域分为个 281节点和565个单元, 所有单元均为三角形单元，如图2所示实际上由于matlab计算编程是不知如何直接读取网 格数据，就只选取了 180个单元与103个节点进行计算。
图2 :左上角四分之一区域的流场及其有限单元划分流场网格
3•单元分析
单元分析的目的是建立有限元方程 。把有限元积分表达式(3)写为各个单元求和的形
2 Q ()
这里()表示单元e ,是单元总数，如果仅在一个单元上考虑上式 ，形式上有
其中()表示单元e的边界，上式实质上并不是一个等式，只具有形式上的意义，当对所 有的单元求和以后，才是等式。如果把线积分中的 2 n (e)换为(e),则得到的是等式，但
在对所有单元求和时，内部边界的线积分刚好抵消，因此(7)也可以理解为不计内部边界 贡献的(3)式在单元上的表达式。
流函数“在单元e内可用如下函数近似：
这里（i = 1 , 2, 3）为节点流函数值，为节点上的插值函数，上式中重复下标表示约定求和 将（8）代入（7）,不难得到
---SQ 3（）
由于3的任意性，所以，对于j=1,2,3都有
---Q （ ）
此即单元方程，通常可以简写为
Q
采用三节点三角形单元时，单元的插值基函数为
（）
如果单元e三个点坐标为 （），i=1,2,3,则
3 （）
即插值基函数在点取1 ,在 两点为零。由此不难解出abc。注意到求
时对取了梯度，故