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

1、有限元法解圆柱绕流问题本文档如对你有帮助，请帮忙下载支持！有限元法求解无限流体中的圆柱绕流问题2016年01月12号一问题描述考虑位于两块无限长平板间的圆柱体的平面绕流问题，几何尺寸如下图所示，来流为。由于流场具有上下左右的对称性，只考虑左上角四分之一的计算区域 abcde，把它作为有限元的求解区域。要求求解出整个区域中的流函数、 以及压强值。图 1：圆柱绕流模型二数学建模在足够远前方选与来流垂直的控制面 ae,cd 是沿 y轴，亦即一流动对称轴， bc是物面， ab亦是流动对称轴，所要考虑的流动区域即由线 abcdea 所围成的区域，在这一区域 中有：1. 边界 ab为流线，取 =0，2.

2、边界 bc也为流线，同样取 =0，3. 边界 cd，切向速度 =0， 取4.边界 de为流线 , 满足于是在 ed上， =2，5. 进口边界 ae上， = （本文中采取此条件）也可以提自然边界条件我们以流函数作为未知函数来解此问题，流函数所满足的微分方程如下：（ 1）此处 就是就是 cd 段边界， 且切向速度 = 0， 1 和 2 合起来是整个边界，并且此二者不重合。下面，按有限元方法的一般步骤来计算此问题。三有限元法解圆柱绕流问题1建立有限元积分表达式根据求解问题的基本控制方程， 应用变分法或加权余量法将求解的微分方程定解问题化为等价的积分表达式，作为有限元法求解问题的出发方程式。对于方程（

3、 1），它是一椭圆型方程，具有正定性，可以用变分法，这里直接给出泛函令其变分 J=0，可以得到本文档如对你有帮助，请帮忙下载支持！自然边界条件已经包含在变分表达式中 ( 其名称的由来 ) ，而本质边界条件必须强制 满足 ( 因此称其为本质边界条件，也称为强制边界条件 ) 。如果根据原微分方程中无法给出泛函 J，则可以用 Galerkin 加权余量方法得到积分方程，这相当于将原来的微分方程写为如下变分形式：这里的是函数的改变量，是一种“虚位移”，在本质边界条件 。因此，上式做分部积分后，边界积分仅剩下 。具体为即（ 3）式。可见， 如果满足原来的微分方程和边界条件，那么，必然有满足 (4) 式，

4、进而满足(5) 式。注意，在 (5) 式中，包含的边界 2 上的边界条件信息，对边界 1 的部分，仅知道它是给定了函数值的边界， 却不知道边界上的值是多少， 为了确定这些值， 还需要额外的处理方法。正是因为 2 上的边界信息可以包含在积分表达式中，这种边界条件也称为自然边界条件。2区域剖分根据物理问题的特点以及区域的形状， 把计算区域分成许多几何形状规则但大小可以不同的单元， 确定单元节点的数目和位置， 建立表示网格的数据结构。 采用的单元形状和节点的分布，以及插值函数的选取还应考虑到计算精度和可微性的要求。 这里通过 ANSYSICEM CFD14.0 建模并划分网格。 具体而言， 网格将求

5、解区域分为个 281节点和 565个单元， 所有单元均为三角形单元，如图 2所示实际上由于 matlab 计算编程是不知如何直接读取网格数据，就只选取了 180个单元与 103个节点进行计算。图 2：左上角四分之一区域的流场及其有限单元划分流场网格3.单元分析单元分析的目的是建立有限元方程。 把有限元积分表达式 (3) 写为各个单元求和的形式这里 表示单元 e， 是单元总数，如果仅在一个单元上考虑上式，形式上有其中 表示单元 e的边界，上式实质上并不是一个等式，只具有形式上的意义，当对所有的单元求和以后，才是等式。如果把线积分中的 2 ( e) 换为 ( e) ，则得到的是等式，但在对所有单元

6、求和时， 内部边界的线积分刚好抵消， 因此 (7) 也可以理解为不计内部边界贡献的 (3) 式在单元上的表达式。流函数在单元 e内可用如下函数近似：这里 ( i = 1 , 2, 3) 为节点流函数值， 为节点上的插值函数，上式中重复下标表示约本文档如对你有帮助，请帮忙下载支持！定求和。将 (8) 代入 (7) ，不难得到由于 的任意性，所以，对于 j=1,2,3 都有此即单元方程，通常可以简写为采用三节点三角形单元时，单元的插值基函数为如果单元 e三个点坐标为（ ） ,i=1,2,3, 则即插值基函数 在 点取 1，在 两点为零。由此不难解出 abc 。注意到求 时对的取值并不影响最后的计算

7、。对某一固定的单元 e，将（ 11）式代入（ 10）中，可以得到：此即采用线性单元时的单元方程系数矩阵。其中 为三角形（积分区域）的面积， bc的值可由（ 12）求得，现在列举如下：c11( x3x2 ) c21(x1x3 )c31( x2x1)2 A(e )2A(e )2 A(e)（13）求解单元系数矩阵时， 一般同时进行总体合成，每形成一个单元方程，便把它累加到总体方程中。出于顺序和逻辑上的考虑，下一步再详细说明总体合成的方法。对于边界积分项， 我们假设三角形单元 e中序号为 , ，的节点在边界2 上，为自然边界，其长度为 l 。首先，注意到插值函数在边上是零。所以，可以得到如下结论：图

8、3：自然边界条件的处理右端项： f ( e)0。为了计算 f ( e) 和 f(e) ，以 点为原点，沿直线建立局部坐标系，在此坐标系中，插值函数 N 和 N 如上图所示，可写成线性插值函数如下：本文档如对你有帮助，请帮忙下载支持！假设切向速度 Vs 在两节点处的值分别为 Vs 和 Vs ，并且沿边界是线性分布的，可以表示为vsvsa(vsvsa ) l。于是可以得到lll (vf (e )v N d(v(vv) ) d2v) 。0s0sassall6sas对于前面讨论的圆柱绕流问题，由于vsa vs0，所以，根据线性解的性质，必有无需考虑 f 的影响，使程序得到了不少的简化。4. 总体合成总

9、体合成的过程就是把已经形成好的单元方程按一定顺序迭加起来，形成总体有限元方程。具体做法是根据单元内节点的总体顺序号，把单元方程进行延拓，未知量包含所有节点上的函数值， 与此单元无关的位置以零填充，把所有的单元方程都进行延拓以后，进行系数矩阵的累加， 便得到总体方程。 理论上说， 这一过程也可以通过引入一个Boole型矩阵来实现，定义单元 e的 boole 矩阵 B3 N pB3 N p1,如果单元 e的地 i 个点总体顺序号是 jj=1,2,3Np.0.,i=1,2,3;其他情况矩阵 B其实就是单元节点序号表的又一表达形式。单元e的系数矩阵 A( e)以及右端项 f ( e)沿拓后就是：进而总

10、体合成的过程可以表示为NeNeAA(e) ，ff (e) 。e 1e 1但是这种方法比较麻烦， 要重新定义新的矩阵，而且还要涉及计算矩阵转置和矩阵相乘等运算， 一方面计算量较大， 并且浪费空间，另一方面人为地增加了程序的复杂性，降低了程序的可读性。因此，这种方法一般只用作理论分析。实际的计算中，每当计算出一个单元系数矩阵，假设单元 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,

11、3*3项分别对应总体刚度矩阵 A的 j*i,j*j,j*k,k*i,k*j,k*k项中。采用此方法并未多占用计算机内存，运算量也并不大（总共进行 9*e 次加法运算，不进行乘法运算）。本题的算法中选用的即此方法。（5）边界条件处理这里的边界条件是指本质边界条件， 自然边界条件已经包含在积分表达式中了。 具体做法是将已知的值代入到方程组中，并把已知的值移到方程组的右段，形成右端项。本文档如对你有帮助，请帮忙下载支持！（6）求解有限元方程组并计算相关物理量有限元方程组的求解是一个代数问题， 应针对具体的问题采用合适的方法求解。 对于对角占优的代数方程组， 可以采用迭代法求解， 规模不大的可以用 G

12、auss 消元法一类的直接法求解， 三对角方程则可以用追赶法。 求出所有的待求量后，便得到了近似函数的表达式，并可以计算出相关的物理量。 对计算结果进行综合的分析， 以期得到原问题的正确的物理解答。对于每个单元，速度可以根据来计算 , 节点上的速度值可取这个节点相邻单元的速度值的平均。节点上的压力值可以有伯努利方程计算。假设求解区域位于同一水平面内，介质密度 ，来流压力 p=0，那么。如此便可得到节点处的速度和压力分布。四具体算法实现本文具体算法是通过 MATLAB实现的， matlab 程序文件及算法说明见附件。五结果分析运行附件中 MATLAB 程序，可以得到图 4和图 6两张图。图 4显

13、示 1/4 区域的流场分布，图中的点以不同颜色表示各个结点处的流函数，图中的曲线为流函数的等势线，即是流线，由于对称性，整个区域的流场分布可明显看出，故不再画出。而圆柱绕流问题的解析解为y( 11) ，便于与计算结果对比。从图4可以发现，有限元法数值法得到的y 2x 2流线与解析法得到的流线在形态上基本一致。图 4：有限元计算的 1/4 区域流场分布图 5：解析法的整个区域流场分布具体到 1/4 区域的右边界上的结点，将有限元法得到的流函数数值解与解析解相对比，得到图 7 中的曲线。 可见， 有限元法与解析法所得曲线在趋势上基本一致， 但在数值上最大误差为 11%。图 6: 1/4区域右边界上

14、流函数的有限元数值解与解析解对比附件：1.NATLAB主程序： youxianyuan.mENA=41 32 37 0.041 44 32 0.018 101 21 0.718 103 101 0.0103 22 102 0.0103 18 22 0.030 36 28 0.030 35 36 0.029 2 12 0.029620.030 7 17 0.0502916230870.0本文档如对你有帮助，请帮忙下载支持！19 88 76 0.019 20 88 0.0101 102 99 0.0101 103 102 0.044 94 5 0.044 41 94 0.035 52 38 0.0

15、35 16 52 0.029 44 6 0.029 32 44 0.097 20 21 0.213930.045940.0102 100 99 0.0101 98 21 0.022 100 102 0.099 98 101 0.098 97 21 0.023 96 22 0.096 100 22 0.092 95 100 0.032213795 99 100 0.091 98 99 0.090 97 98 0.097 89 20 0.048730.093 24 1 0.084 96 23 0.096 92 100 0.080 95 92 0.095 91 99 0.091 90 98 0.09

16、0 89 97 0.089 88 20 0.0612477183 13 19 0.074 94 41 0.0515515694 87 4 0.082 93 3 0.081 24 93 0.086 85 23 0.085 84 23 0.084 92 96 0.092 71 80 0.080 91 95 0.0本文档如对你有帮助，请帮忙下载支持！79 90 91 0.078 89 90 0.077 88 89 0.083 14 13 0.014 75 15 0.074 87 94 0.087 82 3 0.082 81 93 0.024 86 23 0.081 86 24 0.073 85 86

17、 0.072 84 85 0.084 71 92 0.080 79 91 0.079 78 90 0.078 77 89 0.077 76 88 0.069 83 19 0.083 75 14 0.074 82 87 0.068 81 82 0.081 73 86 0.073 72 85 0.0328008272 71 84 0.070 79 80 0.066 78 79 0.065 77 78 0.064 76 77 0.0363831876 63 19 0.069 75 83 0.075 58 15 0.057 74 41 0.074 68 82 0.068 73 81 0.062 72

18、73 0.061 71 72 0.071 67 80 0.067 70 80 0.059 70 56 0.070 66 79 0.0431204266 65 78 0.065 64 77 0.064 63 76 0.063 69 19 0.0本文档如对你有帮助，请帮忙下载支持！69 58 75 0.057 68 74 0.0426945868 62 73 0.062 61 72 0.061 67 71 0.067 56 70 0.059 66 70 0.055 65 66 0.054 64 65 0.0400881353 63 64 0.063 58 69 0.057 62 68 0.051

19、61 62 0.061 56 67 0.056 60 59 0.050 60 47 0.049 59 60 0.059 55 66 0.055 54 65 0.0398997954 53 64 0.053 58 63 0.058 48 15 0.043 57 41 0.0476315157 51 62 0.0419565551 56 61 0.056 47 60 0.050 49 60 0.049 55 59 0.046 54 55 0.045 53 54 0.0418020753 48 58 0.048 52 15 0.052 16 15 0.043 51 57 0.051 47 56 0.

20、033 50 40 0.042 49 50 0.049 46 55 0.046 45 54 0.0420780145 48 53 0.044871948 38 52 0.043 47 51 0.047 40 50 0.033 42 50 0.0本文档如对你有帮助，请帮忙下载支持！42 46 49 0.039 45 46 0.0424682545 38 48 0.044560.043 40 47 0.034 42 33 0.042 39 46 0.039 38 45 0.041 37 43 0.037 40 43 0.0348392340 31 33 0.034 39 42 0.04382332

21、36 38 39 0.037 31 40 0.034 36 39 0.036 35 38 0.032 31 37 0.035 17 16 0.034 28 36 0.033 27 34 0.025 31 32 0.031 26 33 0.0481510830 17 35 0.027 28 34 0.026 27 33 0.0461411229 25 32 0.025 26 31 0.028 8 30 0.027 9 28 0.026 10 27 0.012 25 29 0.025 11 26 0.0333728398280.010 9 27 0.011 10 26 0.012 11 25 0.

22、0; %单元与结点号、面积对应关系矩阵NCORD=-3 0-1 0-2.6 0-2.2 0-1.8 0-1.4 00 1本文档如对你有帮助，请帮忙下载支持！-0.3 0.1-0.6 0.2-0.1 0.1-0.2 0.5-0.1 0.21312.612.211.811.4 -3 3 -0.75 3 -1.5 3 -2.25 3 -3 2.25 -3 1.5 -3 0.75-1.8 0.9-0.6 0.-0.7 0.8-0. 1.8-1.2 0.2-0.8 1.3-1.6 0.9-1. 0.-1.0 1.0-0. 1.8-0.5 1.6-0.5 1.3-1.3 0.6-0.4 1.3-0.1 1

23、.8-1.1 0.7-1.2 0.5745713-1.0 1.5-1.8 0.5-1.6 0.9-0.8 1.4-1.0 1.-1.6 1.4-0.1 2.05857237-1.6 1.3-1. 1.4-1.1 1.9本文档如对你有帮助，请帮忙下载支持！-0.6 1.9-0. 2.-1.0 2.0-1.8 1.7-1.1 1.1-2.0 0.6-0.9 2.1-1.3 1.1-1. 1.1-2. 1.4-2. 1.0-0.1 2.3-1.0 2.4-1.8 2.5-1.2 1.7-2.0 1.5-2.6 0.4-0.2 2.7-1.3 1.4-2.1 1.-2.9 1.2-2. 1.0-2.2

24、 0.7-0.9 2.4-0.8 2.1-1.7 2.7-1.4 2.-1. 2.8-2.2 1.5-2.7 0.7513903-2.4 0.8-0.2 2.8-2.7 1.2-2.1 1.4-2.3 1.0-2.9 0.3-1.6 2.8-1.8 2.1-1.8 2.3-2.5 2.-2.4 1.8-2.1 0.6-1.7 0.-2.5 2.0本文档如对你有帮助，请帮忙下载支持！-2.4 1.2-1.7 2.4-2.0 2.4-2.1 2.7-2.4 2.4-2.5 2.1-2.6 2.8-2.5 2.; %结点与坐标对应关系矩阵Enum,temp=size(ENA); % 返回单元总数Nn

25、um,temp=size(NCORD); % 返回结点总数BNODE=1 3 4 5 6 2 12 11 10 9 8 7 17 16 15 14 13 19 20 21 18 22 23 24;%边界结点编号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,

26、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); % 未知流函数的结点数for k=1:Enum % 单元号x1=NCORD(ENA(k,1),1);y1=NCORD(ENA(k,1)

27、,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;for n=1:3for m=1:3ANM(k,n,m)=1/(4*ENA(k,4)*(b(n)*b(m)+c(n)*c(m);% 各个单元的 Anm 信息 endendendfor i=1

28、:Nnum % 求系数矩阵 Afor j=1:Nnumtemp=0;for e=1:Enumfor n=1:3for m=1:3if ENA(e,n)=i & ENA(e,m)=j本文档如对你有帮助，请帮忙下载支持！temp=temp+ANM(e,n,m);endendendendA(i,j)=temp;endendfor i=1:Uknum % 扫描未知流函数的结点F(i)=0;for j=1:UknumANEW(i,j)=A(UKN(i),UKN(j); %新的系数矩阵endfor k=1:Bknum % 扫描已知流函数的结点F(i)=F(i)+A(UKN(i),BKN(1,k)*BKN(2