有限元法解圆柱绕流问题.docx
- 文档编号:14513230
- 上传时间:2023-06-24
- 格式:DOCX
- 页数:29
- 大小:282.46KB
有限元法解圆柱绕流问题.docx
《有限元法解圆柱绕流问题.docx》由会员分享,可在线阅读,更多相关《有限元法解圆柱绕流问题.docx(29页珍藏版)》请在冰点文库上搜索。
有限元法解圆柱绕流问题
有限元法求解无限流体中的圆柱绕流问题
2016年01月12号
一.问题描述
考虑位于两块无限长平板间的圆柱体的平面绕流问题,几何尺寸如下图所示,来流为vx=1,vy=0。
由于流场具有上下左右的对称性,只考虑左上角四分之一的计算区域abcde,
把它作为有限元的求解区域Ω。
要求求解出整个区域中的流函数、vx、vy以及压强值。
图1:
圆柱绕流模型
二.数学建模
在足够远前方选与来流垂直的控制面ae,cd是沿y轴,亦即一流动对称轴,bc是物面,ab亦是流动对称轴,所要考虑的流动区域即由线abcdea所围成的区域,在这一区域
中有:
1.边界ab为流线,取ψ=0,∂φ∂n=0;
2.边界bc也为流线,同样取ψ=0,∂φ∂n=0;
3.边界cd,切向速度vτ=0,∂ψ∂n=0,取φ=0;
4.边界de为流线,满足ψe=ψa+aevxdy=02dy=2
于是在ed上,ψ=2,∂φ∂n=0;
5.进口边界ae上,ψ=ψ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.建立有限元积分表达式
根据求解问题的基本控制方程,应用变分法或加权余量法将求解的微分方程定解问题化为等价的积分表达式,作为有限元法求解问题的出发方程式。
对于方程
(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上的边界信息可以包含在积分表达式中,这种边界条件也称为自然边界条件。
2.区域剖分
根据物理问题的特点以及区域的形状,把计算区域分成许多几何形状规则但大小可以不
同的单元,确定单元节点的数目和位置,建立表示网格的数据结构。
采用的单元形状和节点
的分布,以及插值函数的选取还应考虑到计算精度和可微性的要求。
这里通过ANSYSICEMCFD14.0建模并划分网格。
具体而言,网格将求解区域分为个281节点和565个单元,所有单元均为三角形单元,如图2所示实际上由于matlab计算编程是不知如何直接读取网格数据,就只选取了180个单元与103个节点进行计算。
图2:
左上角四分之一区域的流场及其有限单元划分流场网格
3.单元分析
单元分析的目的是建立有限元方程。
把有限元积分表达式(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Γ
采用三节点三角形单元时,单元的插值基函数为
Nix,y=aie+biex+ciey(11)
如果单元e三个点坐标为(xie,yie),i=1,2,3,则
Nixje,yje=δij(12)
即插值基函数Ni在xie点取1,在xjexke两点为零。
由此不难解出abc。
注意到求
Aij(e)时对Ni取了梯度,故aie的取值并不影响最后的计算。
对某一固定的单元e,将(11)式代入(10)中,可以得到:
Aij(e)=A(e)b12+c12b1b2+c1c2b1b3+c1c3b1b2+c1c2b22+c22b2b3+c2c3b1b3+c1c3b2b3+c2c3b32+c32
此即采用线性单元时的单元方程系数矩阵。
其中A(e)为三角形(积分区域)的面积,bc的值可由(12)求得,现在列举如下:
(13)
A(e)=12[x2-x1y3-y1-(y2-y1)(x3-x1)]
求解单元系数矩阵时,一般同时进行总体合成,每形成一个单元方程,便把它累加到总体方程中。
出于顺序和逻辑上的考虑,下一步再详细说明总体合成的方法。
对于边界积分项,我们假设三角形单元e中序号为
,的节点在边界
上,
为自然边界,其长度为l。
首先,注意到插值函数在
边上是零。
所以,可以得到如下结论:
图3:
自然边界条件的处理
右端项:
。
为了计算
和
,以
点为原点,沿直线
建立局部坐标系
,在此坐标系中,插值函数
和
如上图所示,可写成线性插值函数如下:
假设切向速度
在两节点处的值分别为
和
,并且沿边界是线性分布的,可以表示为
。
于是可以得到
。
对于前面讨论的圆柱绕流问题,由于
,所以,根据线性解的性质,必有
fα(e)=fβ(e)=fγ(e)=0
无需考虑f的影响,使程序得到了不少的简化。
4.总体合成
总体合成的过程就是把已经形成好的单元方程按一定顺序迭加起来,形成总体有限元方程。
具体做法是根据单元内节点的总体顺序号,把单元方程进行延拓,未知量包含所有节点上的函数值,与此单元无关的位置以零填充,把所有的单元方程都进行延拓以后,进行系数矩阵的累加,便得到总体方程。
理论上说,这一过程也可以通过引入一个Boole型矩阵来实现,定义单元e的boole矩阵
i=1,2,3;j=1,2,3…Np.
矩阵B其实就是单元节点序号表的又一表达形式。
单元e的系数矩阵
以及右端项
沿拓后就是:
进而总体合成的过程可以表示为
,
。
但是这种方法比较麻烦,要重新定义新的矩阵B,而且还要涉及计算矩阵转置和矩阵相乘等运算,一方面计算量较大,并且浪费空间,另一方面人为地增加了程序的复杂性,降低了程序的可读性。
因此,这种方法一般只用作理论分析。
实际的计算中,每当计算出一个单元系数矩阵Aij(e),假设单元e三个节点编号分别为i,j,k,那么将Aij(e)中的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消元法一类的直接法求解,三对角方程则可以用追赶法。
求出所有的待求量后,便得到了近似函数的表达式,并可以计算出相关的物理量。
对计算结果进行综合的分析,以期得到原问题的正确的物理解答。
对于每个单元,速度可以根据
vx=∂ψ∂y=∂ψiNi∂y=ψi∂Ni∂y=ciψi(14)
vy=-∂ψ∂x=-∂ψiNi∂x=-ψi∂Ni∂x=-biψi(15)
来计算,节点上的速度值可取这个节点相邻单元的速度值的平均。
节点上的压力值可以有伯努利方程计算。
假设求解区域位于同一水平面内,介质密度ρ=1,来流压力p=0,那么p=121-vi2,i=1,2…10。
如此便可得到节点处的速度和压力分布。
四.具体算法实现
本文具体算法是通过MATLAB实现的,matlab程序文件及算法说明见附件。
五.结果分析
运行附件中MATLAB程序,可以得到图4和图6两张图。
图4显示1/4区域的流场分布,图中的点以不同颜色表示各个结点处的流函数,图中的曲线为流函数的等势线,即是流线,由于对称性,整个区域的流场分布可明显看出,故不再画出。
而圆柱绕流问题的解析解为
,便于与计算结果对比。
从图4可以发现,有限元法数值法得到的
流线与解析法得到的流线在形态上基本一致。
图4:
有限元计算的1/4区域流场分布
图5:
解析法的整个区域流场分布
具体到1/4区域的右边界上的结点,将有限元法得到的流函数数值解与解析解相对比,
得到图7中的曲线。
可见,有限元法与解析法所得曲线在趋势上基本一致,但在数值上最大
误差为11%。
图6:
1/4区域右边界上流函数的有限元数值解与解析解对比
附件:
1.NATLAB主程序:
youxianyuan.m
ENA=[4132370.039798301
4144320.047529903
18101210.120156817
181031010.094963772
103221020.043498437
10318220.094243318
3036280.038054662
3035360.042287521
292120.038074502
29620.048742904
307170.05029162
30870.037140964
1988760.057994041
1920880.099591586
101102990.046382644
1011031020.044218891
449450.047132639
4441940.055967406
3552380.036085513
3516520.034571557
294460.044280863
2932440.041448145
9720210.106109192
13930.078212613
45940.054017838
102100990.034088039
10198210.050359216
221001020.048451565
99981010.045186467
9897210.070310629
2396220.097022072
96100220.062253661
92951000.0322137
95991000.031691011
9198990.039495829
9097980.038409873
9789200.062852728
48730.050497897
932410.094904332
8496230.060406457
96921000.036472833
8095920.032029748
9591990.029674024
9190980.041864912
9089970.045507083
8988200.06124771
8313190.093342661
7494410.05155156
948740.039208883
829330.058277066
8124930.063798865
8685230.047832817
8584230.050837396
8492960.038920434
9271800.041373741
8091950.033559408
7990910.042160395
7889900.044890766
7788890.041466167
8314130.060586214
1475150.051106938
7487940.041578206
878230.057195134
8281930.045533139
2486230.084555213
8186240.056241682
7385860.030172685
7284850.032814285
8471920.037263856
8079910.047262673
7978900.043683491
7877890.043858618
7776880.033382546
6983190.053167481
8375140.030301151
7482870.039566275
6881820.045303658
8173860.039028656
7372850.03280082
7271840.036391036
7079800.050409968
6678790.043207084
6577780.040783626
6476770.03638318
7663190.060458755
6975830.025977524
7558150.055241265
5774410.053616644
7468820.046847336
6873810.041793554
6272730.036075451
6171720.035939511
7167800.038472382
6770800.040185579
5970560.045615694
7066790.04312042
6665780.041221814
6564770.038815812
6463760.035503832
6369190.057736703
6958750.032404442
5768740.04269458
6862730.039041739
6261720.037749165
6167710.033244944
6756700.038889688
5966700.040818003
5565660.040585183
5464650.04008813
5363640.042912099
6358690.036059255
5762680.041591565
5161620.039685976
6156670.037633909
5660590.036288344
5060470.031860179
4959600.031521434
5955660.039419363
5554650.03989979
5453640.040903813
5358630.046404299
5848150.066170074
4357410.04763151
5751620.04195655
5156610.046636447
5647600.035592793
5049600.031796449
4955590.039229272
4654550.040974288
4553540.04180207
5348580.042196026
4852150.056305234
5216150.047154585
4351570.041090301
5147560.045821293
3350400.048398755
4249500.040513663
4946550.040202562
4645540.04207801
4548530.0448719
4838520.037263593
4347510.041538735
4740500.039865602
3342500.050325014
4246490.041023271
3945460.04246825
4538480.043951365
44560.051047374
4340470.041787595
3442330.051746427
4239460.042853263
3938450.041743856
4137430.038609599
3740430.03483923
4031330.050574027
3439420.04382332
3638390.042210941
3731400.033185168
3436390.041698152
3635380.039494122
3231370.034698998
3517160.049191615
3428360.040629125
3327340.049558416
2531320.039551559
3126330.04815108
3017350.041250709
2728340.039846508
2627330.04614112
2925320.037979443
2526310.038580787
288300.036116702
279280.033481704
2610270.032757535
1225290.035427736
2511260.03337283
98280.031904711
109270.030445747
1110260.030285297
1211250.031849171
];%单元与结点号、面积对应关系矩阵
NCORD=[-30
-10
-2.60
-2.20
-1.80
-1.40
01
-0.2588241030.965924471
-0.4999954660.866028022
-0.7071067810.707106781
-0.8660280220.499995465
-0.9659244710.258824102
03
02.6
02.2
01.8
01.4
-33
-0.753
-1.53
-2.253
-32.25
-31.5
-30.75
-1.1345216680.489438169
-0.9708951760.74446511
-0.7573204670.962580378
-0.505508741.128325608
-1.2621251320.243714522
-0.2514580981.253891963
-1.2771060360.738778719
-1.447512260.48199085
-1.0881676511.056783565
-0.775236311.267266548
-0.2459580751.585179826
-0.5038471451.428730143
-1.5487045030.736752996
-0.4892388641.743880003
-0.7646746611.580844378
-1.4507575310.981852867
-1.7967204420.5745713
-1.0498231171.413295395
-1.7651239080.906580515
-1.6189897260.255236869
-0.7496933581.892823424
-1.0299777561.72552432
-1.6491440061.200203744
-0.4493600412.05857237
-1.3030399261.563704873
-1.347538531.270144924
-1.9546564711.143055159
-0.2357729261.875193029
-0.740660042.19662276
-1.0207719462.031271425
-1.2938474281.863609647
-1.8253854811.467199841
-2.0584840570.839001806
-0.4554460192.351164971
-1.5606991231.692649341
-1.538865611.437047451
-2.150601831.373255674
-2.250340161.085358385
-0.7482476712.517911273
-1.0109146942.329143574
-1.2856465682.160870135
-1.5635246021.986279407
-2.0663439841.629034175
-2.3546506260.785729864
-0
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 圆柱 问题