有限元法解圆柱绕流问题

有限元法求解无限流体中的圆柱绕流问题

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 14.0建模并划分网格。具体而言，网格将求解区域分为个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）求得，现在列举如下：

123()1()2e b y y A =- 231()1()2e b y y A =- 3

12()1

()2e b y y A =-

132()1()2e c x x A =

- 213()1()2e c x x A =- 3

21()1

()2e c x x A =-（13）

求解单元系数矩阵时，一般同时进行总体合成，每形成一个单元方程，便把它累加到

总体方程中。出于顺序和逻辑上的考虑，下一步再详细说明总体合成的方法。

对于边界积分项，我们假设三角形单元e 中序号为,αβ，的节点在边界2Γ上，αβ为自然边界，其长度为l 。首先，注意到插值函数在αβ边上是零。所以，可以得到如下结论：

图3：自然边界条件的处理

右端项：()0e f γ=。

为了计算()

e f α和()e f β，以α点为原点，沿直线αβ建立局部坐标系ξ，在此坐标系中，

插值函数N α和N β如上图所示，可写成线性插值函数如下：

N N l l

αβξξ

=1-, =

假设切向速度s V 在两节点处的值分别为s V α和s V β，并且沿边界是线性分布的，可以

表示为

()

s sa s sa v v v v l

βξ

=+-。

于是可以得到

()

00

(())(1)(2)6l l

e s sa s sa sa s l

f v N d v v v d v v l l ααββξξξξ=-=-+--=-+⎰⎰

()

00

(())(2)6l l

e s sa s sa sa s l

f v N d v v v d v v l l ββββξξξξ=-=-+-=-+⎰⎰。