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

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

《有限元法解圆柱绕流问题.docx》由会员分享，可在线阅读，更多相关《有限元法解圆柱绕流问题.docx（26页珍藏版）》请在冰豆网上搜索。

有限元法解圆柱绕流问题

本文档如对你有帮助，请帮忙下载支持！

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

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．区域剖分

根据物理问题的特点以及区域的形状，把计算区域分成许多几何形状规则但大小可以不

同的单元，确定单元节点的数目和位置，建立表示网格的数据结构。

采用的单元形状和节点

的分布，以及插值函数的选取还应考虑到计算精度和可微性的要求。

这里通过ANSYSICEMCFD

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）求得，现在列举如下：

c1

1

（x3

x2）c2

1

（x1

x3）

c3

1

（x2

x1）

2A

（e）

2A

（e）

2A

（e）

（13）

求解单元系数矩阵时，一般同时进行总体合成，

每形成一个单元方程，

便把它累加到总

体方程中。

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

对于边界积分项，我们假设三角形单元e中序号为,，的节点在边界

2上，

为自

然边界，其长度为l。

首先，注意到插值函数在

边上是零。

所以，可以得到如下结论：

图3：

自然边界条件的处理

右端项：

f（e）

0

。

为了计算f（e）和f

（e），以点为原点，沿直线

建立局部坐标系

，在此坐标系中，

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

本文档如对你有帮助，请帮忙下载支持！

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

示为

vs

vsa

（vs

vsa）l

。

于是可以得到

l

l

l（v

f（e）

vNd

（v

（v

v

））d

2v

）。

0

s

0

sa

s

sa

l

l

6sa

s

对于前面讨论的圆柱绕流问题，由于

vsavs

0

，所以，根据线性解的性质，必有

无需考虑f的影响，使程序得到了不少的简化。

4.总体合成

总体合成的过程就是把已经形成好的单元方程按一定顺序迭加起来，

形成总体有限元方

程。

具体做法是根据单元内节点的总体顺序号，

把单元方程进行延拓，

未知量包含所有节点

上的函数值，与此单元无关的位置以零填充，

把所有的单元方程都进行延拓以后，

进行系数

矩阵的累加，便得到总体方程。

理论上说，这一过程也可以通过引入一个

Boole

型矩阵来实

现，定义单元e的boole矩阵B3Np

B3Np

1,如果单元e的地i个点总体顺序号是j

j=1,2,3

⋯Np.

0.

i=1,2,3;

其他情况

矩阵B其实就是单元节点序号表的又一表达形式。

单元

e的系数矩阵A（e）

以及右端项f（e）

沿拓后就是：

进而总体合成的过程可以表示为

Ne

Ne

A

A（e），f

f（e）。

e1

e1

但是这种方法比较麻烦，要重新定义新的矩阵Ｂ，

而且还要涉及计算矩阵转置和矩阵相

乘等运算，一方面计算量较大，并且浪费空间，另一方面人为地增加了程序的复杂性，

降低

了程序的可读性。

因此，这种方法一般只用作理论分析。

实际的计算中，每当计算出一个单元系数矩阵

，假设单元e三个节点编号分别为

i,j,k，那么将

中的1*1项放入大矩阵（借鉴结构力学的概念，不妨称其为总体“刚度”矩

阵，下同）A的i*i

项中，将1*2,1*3

分别放入总体刚度矩阵A的i*j,i*k

项中。

同理，2*1,2*2

，

2*3,3*1,3*2,3*3

项分别对应总体刚度矩阵A的j*i,j*j,j*k,k*i,k*j,k*k

项中。

采用此方法

并未多占用计算机内存，运算量也并不大（总共进行9*e次加法运算，不进行乘法运算）。

本题的算法中选用的即此方法。

（5）边界条件处理

这里的边界条件是指本质边界条件，自然边界条件已经包含在积分表达式中了。

具体做法是将已知的值代入到方程组中，并把已知的值移到方程组的右段，形成右端项。

本文档如对你有帮助，请帮忙下载支持！

（6）求解有限元方程组并计算相关物理量

有限元方程组的求解是一个代数问题，应针对具体的问题采用合适的方法求解。

对于对角占优的代数方程组，可以采用迭代法求解，规模不大的可以用Gauss消元法一类的直接法求解，三对角方程则可以用追赶法。

求出所有的待求量后，便得到了近似函数的表达式，并

可以计算出相关的物理量。

对计算结果进行综合的分析，以期得到原问题的正确的物理解答。

对于每个单元，速度可以根据

来计算,节点上的速度值可取这个节点相邻单元的速度值的平均。

节点上的压力值可以有伯

努利方程计算。

假设求解区域位于同一水平面内，介质密度，来流压力p=0，那么

。

如此便可得到节点处的速度和压力分布。

四．具体算法实现

本文具体算法是通过MATLAB实现的，matlab程序文件及算法说明见附件。

五．结果分析

运行附件中MATLAB程序，可以得到图4和图6两张图。

图4显示1/4区域的流场分布，图中的点以不同颜色表示各个结点处的流函数，图中的曲线为流函数的等势线，即是流线，由于对称性，整个区域的流场分布可明显看出，故不再画出。

而圆柱绕流问题的解析解为

y（1

1

），便于与计算结果对比。

从图

4可以发现，有限元法数值法得到的

y2

x2

流线与解析法得到的流线在形态上基本一致。

图4：

有限元计算的1/4区域流场分布

图5：

解析法的整个区域流场分布

具体到1/4区域的右边界上的结点，将有限元法得到的流函数数值解与解析解相对比，

得到图7中的曲线。

可见，有限元法与解析法所得曲线在趋势上基本一致，但在数值上最大

误差为11%。

图6:

1/4

区域右边界上流函数的有限元数值解与解析解对比

附件：

1.NATLAB主程序：

youxianyuan.m

ENA=[4132370.0

4144320.0

18101210.7

181031010.0

103221020.0

10318220.0

3036280.0

3035360.0

292120.0

29620.0

307170.05029162

30870.0

本文档如对你有帮助，请帮忙下载支持！

1988760.0

1920880.0

101102990.0

1011031020.0

449450.0

4441940.0

3552380.0

3516520.0

294460.0

2932440.0

9720210.2

13930.0

45940.0

102100990.0

10198210.0

221001020.0

99981010.0

9897210.0

2396220.0

96100220.0

92951000.0322137

95991000.0

9198990.0

9097980.0

9789200.0

48730.0

932410.0

8496230.0

96921000.0

8095920.0

9591990.0

9190980.0

9089970.0

8988200.06124771

8313190.0

7494410.05155156

948740.0

829330.0

8124930.0

8685230.0

8584230.0

8492960.0

9271800.0

8091950.0

本文档如对你有帮助，请帮忙下载支持！

7990910.0

7889900.0

7788890.0

8314130.0

1475150.0

7487940.0

878230.0

8281930.0

2486230.0

8186240.0

7385860.0

7284850.0

8471920.0

8079910.0

7978900.0

7877890.0

7776880.0

6983190.0

8375140.0

7482870.0

6881820.0

8173860.0

7372850.03280082

7271840.0

7079800.0

6678790.0

6577780.0

6476770.03638318

7663190.0

6975830.0

7558150.0

5774410.0

7468820.0

6873810.0

6272730.0

6171720.0

7167800.0

6770800.0

5970560.0

7066790.04312042

6665780.0

6564770.0

6463760.0

6369190.0

本文档如对你有帮助，请帮忙下载支持！

6958750.0

5768740.04269458

6862730.0

6261720.0

6167710.0

6756700.0

5966700.0

5565660.0

5464650.04008813

5363640.0

6358690.0

5762680.0

5161620.0

6156670.0

5660590.0

5060470.0

4959600.0

5955660.0

5554650.03989979

5453640.0

5358630.0

5848150.0

4357410.04763151

5751620.04195655

5156610.0

5647600.0

5049600.0

4955590.0

4654550.0

4553540.04180207

5348580.0

4852150.0

5216150.0

4351570.0

5147560.0

3350400.0

4249500.0

4946550.0

4645540.04207801

4548530.0448719

4838520.0

4347510.0

4740500.0

3342500.0

本文档如对你有帮助，请帮忙下载支持！

4246490.0

3945460.04246825

4538480.0

44560.0

4340470.0

3442330.0

4239460.0

3938450.0

4137430.0

3740430.03483923

4031330.0

3439420.04382332

3638390.0

3731400.0

3436390.0

3635380.0

3231370.0

3517160.0

3428360.0

3327340.0

2531320.0

3126330.04815108

3017350.0

2728340.0

2627330.04614112

2925320.0

2526310.0

288300.0

279280.0

2610270.0

1225290.0

2511260.03337283

98280.0

109270.0

1110260.0

1211250.0

];%单元与结点号、面积对应关系矩阵

NCORD=[-30

-10

-2.60

-2.20

-1.80

-1.40

01

本文档如对你有帮助，请帮忙下载支持！

-0.30.1

-0.60.2

-0.10.1

-0.20.5

-0.10.2

13

12.6

12.2

11.8

11.4-33-0.753-1.53-2.253-32.25-31.5-30.75

-1.80.9

-0.60.

-0.70.8

-0.1.8

-1.20.2

-0.81.3

-1.60.9

-1.0.

-1.01.0

-0.1.8

-0.51.6

-0.51.3

-1.30.6

-0.41.3

-0.11.8

-1.10.7

-1.20.5745713

-1.01.5

-1.80.5

-1.60.9

-0.81.4

-1.01.

-1.61.4

-0.12.05857237

-1.61.3

-1.1.4

-1.11.9

本文档如对你有帮助，请帮忙下载支持！

-0.61.9

-0.2.

-1.02.0

-1.81.7

-1.11.1

-2.00.6

-0.92.1

-1.31.1

-1.1.1

-2.1.4

-2.1.0

-0.12.3

-1.02.4

-1.82.5

-1.21.7

-2.01.5

-2.60.4

-0.22.7

-1.31.4

-2.11.

-2.91.2

-2.1.0

-2.20.7

-0.92.4

-0.82.1

-1.72.7

-1.42.

-1.2.8

-2.21.5

-2.70.7513903

-2.40.8

-0.22.8

-2.71.2

-2.11.4

-2.31.0

-2.90.3

-1.62.8

-1.82.1

-1.82.3

-2.52.

-2.41.8

-2.10.6

-1.70.

-2.52.0

本文档如对你有帮助，请帮忙下载支持！

-2.41.2

-1.72.4

-2.02.4

-2.12.7

-2.42.4

-2.52.1

-2.62.8

-2.52.

];%结点与坐标对应关系矩阵

[Enum,temp]=size（ENA）;%返回单元总数

[Nnum,temp]=size（NCORD）;%返回结点总数

BNODE=[134562121110987171615141319202118222324];%边界结点编号

[temp,Bnum]=size（BNODE）;%返回边界结点数

%边界已知流函数的结点号及其值

BKN=[1345621211109871319202118222324;

000000000000333332.251.50.75];

[temp,Bknum]=size（BKN）;%返回边界已知流函数的结点数

UKN=[14,15,16,17,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,

50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,

81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103];

[temp,Uknum]=size（UKN）;%未知流函数的结点数

fork=1:

Enum%单元号

x1=NCORD（ENA（k,1）,1）;

y1=NCORD（ENA（k,1）,2）;

x2=NCORD（ENA（k,2）,1）;

y2=NCORD（ENA（k,2）,2）;

x3=NCORD（ENA（k,3）,1）;

y3=NCORD（ENA（k,3）,2）;

a

（1）=x2*y3-x3*y2;b

（1）=y2-y3;c

（1）=x3-x2;

a

（2）=x3*y1-x1*y3;b

（2）=y3-y1;c

（2）=x1-x3;

a（3）=x1*y2-x2*y1;b（3）=y1-y2;c（3）=x2-x1;

forn=1:

3

form=1:

3

ANM（k,n,m）=1/（4*ENA（k,4））*（b（n）*b（m）+c（n）*c（m））;%各个单元的Anm信息end

end

end

fori=1:

Nnum%求系数矩阵A

forj=1:

Nnum

temp=0;

fore=1:

Enum

forn=1:

3

form=1:

3

ifENA（e,n）==i&&ENA（e,m）==j

本文档如对你有帮助，请帮忙下载支持！

temp=temp+ANM（e,n,m）;

end

end

end

end

A（i,j）=temp;

end

end

fori=1:

Uknum%扫描未知流函数的结点

F（i）=0;

forj=1:

Uknum

ANEW（i,j）=A（UKN（i）,UKN（j））;%新的系数矩阵

end

fork=1:

Bknum%扫描已知流函数的结点

F（i）=F（i）+A（UKN（i）,BKN（1,k））*BKN（2