有限元法解圆柱绕流报告

有限元法解圆柱绕流计算流体力学期末大作业工学院熊思009860802012/6/3

【摘要】 有限元法是求解计算流体力学问题的重要方法。用有限元法求解具体问题时，首先需要将求解区域进行离散化，即将求解区域划分为许多几何形状简单规那么的单元，在二维一般是三角形或四边形，在三维是四面体或六面体。然后，在每个单元内，用一个比拟简单的解析函数来逼近微分方程的解，此函数在单元内用一组选定的单元基函数的线性组合表示，而其中的系数通常是节点参数，它是待定的。这样，每个单元只要有适当数量的节点参数值，就可以满足对插值函数的光滑性和精度的要求。第三，在满足微分方程和相应的初边值条件下，对全部子域进行积分。对每个单元分别进行积分，形成“单元方程”，通过总体合成，得到总体有限元方程组。最后，用适当的方法解方程组，可得节点参数值，进而可求得各单元内的近似解。本文将以圆柱绕流问题为例，展示有限元法求解的一般步骤。一.有限元法概述有限元法〔finiteelementmethod〕是一种高效能、常用的计算方法。有限元法在早期是以变分原理为根底开展起来的，所以它广泛地应用于以拉普拉斯方程和泊松方程所描述的各类物理场中,这类场与泛函的极值问题有着紧密的联系。有限单元法最早可上溯到20世纪40年代。Courant第一次应用定义在三角区域上的分片连续函数和最小位能原理来求解St.Venant扭转问题。现代有限单元法的第一个成功的尝试是在1956年，Turner、Clough等人在分析飞机结构时，将钢架位移法推广应用于弹性力学平面问题，给出了用三角形单元求得平面应力问题的正确答案。1960年，Clough进一步处理了平面弹性问题，并第一次提出了"有限单元法"，使人们认识到它的成效。自从1969年以来，某些学者在流体力学中应用加权余量法中的迦辽金法(Galerkin)或最小二乘法等同样获得了有限元方程。具体而言，参考差分法中网格化的做法，把求解区域划分为有限多子区域，称这些子区域为单元，在每个单元上构造解的近似分布，将Ritz法或加权余量法应用到分块的逼近函数上。因而有限元法可应用于以任何微分方程所描述的各类物理场中，而不再要求这类物理场和泛函的极值问题有所联系。实质上，有限元法就是Ritz法或加权余量法。二．问题描述 考虑位于两块无限长平板间的圆柱体的平面绕流问题，几何尺寸如下列图所示，来流为vx=1,vy把它作为有限元的求解区域Ω。要求求解出整个区域中的流函数、vx、三．提出物理问题的边界条件和满足的微分方程1.边界ab为流线，取ψ=0，∂2.边界bc也为流线，同样取ψ=0，∂3.边界cd，切向速度vτ=0，∂ψ∂n4.边界de为流线,满足ψ于是在ed上，ψ=2，∂5.进口边界ae上，ψ=ψa+a也可以提自然边界条件∂ψ我们以流函数ψ作为未知函数来解此问题，流函数所满足的微分方程如下：∇2ψ=0此处Γ1指ab,bc,de和ae四段边界，而Γ2就是就是cd段边界，且切向速度v四．有限元法解圆柱绕流问题1．建立有限元积分表达式根据求解问题的根本控制方程，应用变分法或加权余量法将求解的微分方程定解问题化为等价的积分表达式，作为有限元法求解问题的出发方程式。对于方程〔1〕，它是一椭圆型方程，具有正定性，可以用变分法，这里直接给出泛函J令其变分δJ=0，可以得到Ω自然边界条件已经包含在变分表达式中(其名称的由来)，而本质边界条件必须强制ψ满足(因此称其为本质边界条件，也称为强制边界条件)。如果根据原微分方程中无法给出泛函J，那么可以用Galerkin加权余量方法得到积分方程，这相当于将原来的微分方程写为如下变分形式：Ω这里的δψ是函数ψ的改变量，是一种“虚位移”，在本质边界条件Γ1上为零。因此，上式做分部积分后，边界积分仅剩下Ω即〔3〕式。可见，如果ψ满足原来的微分方程和边界条件，那么，必然有ψ满足(4)式，进而满足(5)式。注意，在(5)式中，包含的边界Γ2上的边界条件信息，对边界Γ1的局部，仅知道它是给定了函数值的边界，却不知道边界上的值是多少，为了确定这些值，还需要额外的处理方法。正是因为Γ2上的边界信息可以包含在积分表达式中，这种边界条件也称为自然边界条件。2．区域剖分根据物理问题的特点以及区域的形状，把计算区域分成许多几何形状规那么但大小可以不同的单元，确定单元节点的数目和位置，建立表示网格的数据结构。采用的单元形状和节点的分布，以及插值函数的选取还应考虑到计算精度和可微性的要求。在此题中，网格〔区域剖分〕局部已由GRID.DAT文件给出。具体而言，网格将求解区域分为545个节点和997个单元。文件前半局部给出每个节点的横纵坐标，后半局部给出每个单元对应的三个节点的序号。除此之外，单元剖分还应该建立本质边界条件和自然边界条件的节点表。自然边界条件对解此题并无奉献〔原因在下文会进行说明〕，由网格表特征可知，此题的本质边界条件为abc段上〔节点1-37〕，psi为零；de段上〔节点46-75〕，psi为2；ea段上〔节点76-91〕，psi为纵坐标y。单元剖分实质上就是建立这几个对应关系。3.单元分析单元分析的目的是建立有限元方程。把有限元积分表达式(3)写为各个单元求和的形式e=1这里Ω（e）Ω其中Γ（e）表示单元e的边界，上式实质上并不是一个等式，只具有形式上的意义，当对所有的单元求和以后，才是等式。如果把线积分中的Γ2∩Γ(e)换为Γ(e)，那么得到的是等式，但在对所有单元求和时，内部边界的线积分刚好抵消，因此(7)也可以理解为不计内部边界奉献 流函数ψ在单元e内可用如下函数近似：ψ=这里ψi(i=1,2,3)为节点流函数值，Ni为节点上的插值函数，上式中重复下标表示约定求和。将(8)代入(7Ω由于δψΩ此即单元方程，通常可以简写为A采用三节点三角形单元时，单元的插值基函数为N如果单元e三个点坐标为〔xie,yiN即插值基函数Ni在xie点取1，在xjexke对某一固定的单元e，将〔11〕式代入〔10〕中，可以得到：A此即采用线性单元时的单元方程系数矩阵。其中A(e)为三角形〔积分区域〕的面积，bc的值可由〔12〕〔13〕A求解单元系数矩阵时，一般同时进行总体合成，每形成一个单元方程，便把它累加到总体方程中。出于顺序和逻辑上的考虑，下一步再详细说明总体合成的方法。对于边界积分项，我们假设三角形单元e中序号为，的节点在边界上，为自然边界，其长度为l。首先，注意到插值函数在边上是零。所以，可以得到如下结论：图：自然边界条件的处理 右端项：。 为了计算和，以点为原点，沿直线建立局部坐标系，在此坐标系中，插值函数和如上图所示，可写成线性插值函数如下： 假设切向速度在两节点处的值分别为和，并且沿边界是线性分布的，可以表示为。于是可以得到。对于前面讨论的圆柱绕流问题，由于，所以，根据线性解的性质，必有f无需考虑f的影响，使程序得到了不少的简化。4.总体合成总体合成的过程就是把已经形成好的单元方程按一定顺序迭加起来，形成总体有限元方程。具体做法是根据单元内节点的总体顺序号，把单元方程进行延拓，未知量包含所有节点上的函数值，与此单元无关的位置以零填充，把所有的单元方程都进行延拓以后，进行系数矩阵的累加，便得到总体方程。理论上说，这一过程也可以通过引入一个Boole型矩阵来实现，定义单元e的boole矩阵,i=1,2,3;j=1,2,3…Np.矩阵B其实就是单元节点序号表的又一表达形式。单元e的系数矩阵以及右端项沿拓后就是：进而总体合成的过程可以表示为，。但是这种方法比拟麻烦，要重新定义新的矩阵Ｂ，而且还要涉及计算矩阵转置和矩阵相乘等运算，一方面计算量较大，并且浪费空间，另一方面人为地增加了程序的复杂性，降低了程序的可读性。因此，这种方法一般只用作理论分析。实际的计算中，每当计算出一个单元系数矩阵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消元法一类的直接法求解，三对角方程那么可以用追赶法。求出所有的待求量后，便得到了近似函数的表达式，并可以计算出相关的物理量。对计算结果进行综合的分析，以期得到原问题的正确的物理解答。对于每个单元，速度可以根据vv来计算,节点上的速度值可取这个节点相邻单元的速度值的平均。节点上的压力值可以有伯努利方程计算。假设求解区域位于同一水平面内，介质密度ρ=1，来流压力p=0，那么p=12五．具体算法实现(算法说明)1.建立有限元积分表达式 此局部不在算法中表达。2.区域剖分区域剖分的主要内容已由网格文件完成，我们读入网格文件，并设置边界条件即可。 从GRID.TXT中读入所需根本信息：节点数、单元数、节点坐标、单元包含节点及其次序。为了从零开始编号，应将单元包含节点序号统一减一。并根据网格特征，将流函数psi赋初值：abc段上〔节点1-37〕，psi为零；de段上〔节点46-75〕，psi为2；ea段上〔节点76-91〕，psi为纵坐标y。3.单元分析与总体合成 按公式〔13〕求得单元i的插值基函数的系数b[i][s]c[i][s](s=1,2,3)和系数A(e)，求出bc进而求得Aij4.边界条件处理 逐行处理。当行数(i>=0&&i<37)||(i>=45&&i<91)，即为本质边界条件的行数时，直接将此行所有除ψ[i]对应的系数以外的其他系数全部置零，psi[i]的系数置为一，右端项置为ψ[i]的值。当行数不为本质边界条件的行数时，通过f[i]=f[i]-a[i][j]*psi[j]求右端项，再将除ψ[i]对应的系数以外的其他系数全部置零。5.求解有限元方程组并计算相关物理量。 本算法采用LU分解进行求解。其中，double**LU(double**a,intm,double**l)为LU分解函数。**a即二维数组〔矩阵〕a的首地址,m为矩阵规模，此函数返回LU矩阵的首地址。double*SOVLELU(double**l,double*b,intm,double*x)解LU分解后的方程组。l为LU分解产生的矩阵，b为右端项，m为矩阵规模，x为解向量,返回解向量的首地址。使用这两个函数即可求出所有节点处流函数ψ的值。单元速度的求解那么采用〔14〕〔15〕两式，为了节省空间，不保存单元节点速度，每求出一个单元节点速度，立刻存入其包含的三个节点中。使用count计数，每存一次count加一。最后将每个节点的速度除以count，即取平均，得每个节点的速度。压力使用伯努利方程较易得到，此处不加赘述。6.输出。将结果〔节点坐标、xy方向速度、流函数ψ的值、压强、单元包含节点数〕信息输出到result111.txt中，并释放所有空间。〔为节约内存，局部空间已于之前释放〕。更多详细说明请参见源程序中的注释局部。附源程序：#include<stdio.h>#include<math.h>#include<stdlib.h>voidmain(){double**LU(double**a,intm,double**l);//LU分解函数,**a即二维数组〔矩阵〕a的首地址,此函数返回LU矩阵的首地址，m为矩阵规模double*SOVLELU(double**l,double*b,intm,double*x);//解LU分解后的方程组，l为LU分解产生的矩阵，b为右端项，m为矩阵规模，x为解向量,返回解向量的首地址//voidL_U(double**A,double*f,intm);FILE*fp;//文件指针fp=fopen("GRID.txt","r");//将c\从GRID.txt中读取数据intn,e;//n为节点个数，e为单元个数inti,j;fscanf(fp,"%d",&n);//从文件中读入结点个数fscanf(fp,"%d",&e);//从文件中读入单元个数double*x=newdouble[n];//x[i]即第i+1个节点的横坐标double*y=newdouble[n];//y[i]即第i+1个节点的纵坐标double**l=newdouble*[n];//l即LU分解产生的中间矩阵for(i=0;i<n;i++)l[i]=newdouble[n];for(i=0;i<n;i++){fscanf(fp,"%lf",&x[i]);//从文件中读入结点横坐标fscanf(fp,"%lf",&y[i]);//从文件中读入结点纵坐标}int*first=newint[e];//first[]即三角形单元对应的第一个节点编号(从一开始编号)int*second=newint[e];//second[]即三角形单元对应的第二个节点编号int*third=newint[e];//third[]即三角形单元对应的第三个节点编号for(i=0;i<e;i++){fscanf(fp,"%d",&first[i]);//从文件中读入单元对应的第一个节点编号first[i]=first[i]-1;//减一使节点编号从零开始，便于后面的计算fscanf(fp,"%d",&second[i]);//从文件中读入单元对应的第二个节点编号second[i]=second[i]-1;fscanf(fp,"%d",&third[i]);//从文件中读入单元对应的第三个节点编号third[i]=third[i]-1;}fclose(fp);//至此网格文件读入完毕double*psi=newdouble[n];//psi[]即节点上的流函数值,开始时将其置零for(i=0;i<n;i++)psi[i]=0;double**a=newdouble*[n];//a[][]即系数矩阵〔类似刚度矩阵〕,大小为n*n,并将其置零for(i=0;i<n;i++)a[i]=newdouble[n];for(i=0;i<n;i++)for(j=0;j<n;j++)a[i][j]=0;double*f=newdouble[n];//f为方程右端项,将其置零for(i=0;i<n;i++)f[i]=0;//以下开始赋边界条件,自然边界条件在此题中可以不赋值〔为零，不影响后继计算〕for(i=0;i<37;i++)psi[i]=0;//abc段上，psi为零for(i=45;i<75;i++)psi[i]=2;//de段上，psi为2for(i=75;i<91;i++)psi[i]=y[i];//ea段上，psi为纵坐标y//边界条件赋值完毕double**b=newdouble*[e];for(i=0;i<e;i++)b[i]=newdouble[3];double**c=newdouble*[e];for(i=0;i<e;i++)c[i]=newdouble[3];doubled;//b[e][3],c[e][3]为单元插值函数N(x,y)=a+b*x+c*y的系数,d为系数A〔e〕for(i=0;i<e;i++)//对每一个单元，求出系数矩阵并参加大的系数矩阵中{d=(x[second[i]]-x[first[i]])*(y[third[i]]-y[first[i]])-(y[second[i]]-y[first[i]])*(x[third[i]]-x[first[i]]);d=d/2;b[i][0]=(y[second[i]]-y[third[i]])/(2*d);b[i][1]=(y[third[i]]-y[first[i]])/(2*d);b[i][2]=(y[first[i]]-y[second[i]])/(2*d);c[i][0]=(x[third[i]]-x[second[i]])/(2*d);c[i][1]=(x[first[i]]-x[third[i]])/(2*d);c[i][2]=(x[second[i]]-x[first[i]])/(2*d);a[first[i]][first[i]]=a[first[i]][first[i]]+(b[i][0]*b[i][0]+c[i][0]*c[i][0])*d;a[first[i]][second[i]]=a[first[i]][second[i]]+(b[i][0]*b[i][1]+c[i][0]*c[i][1])*d;a[first[i]][third[i]]=a[first[i]][third[i]]+(b[i][0]*b[i][2]+c[i][0]*c[i][2])*d;a[second[i]][first[i]]=a[second[i]][first[i]]+(b[i][0]*b[i][1]+c[i][0]*c[i][1])*d;a[second[i]][second[i]]=a[second[i]][second[i]]+(b[i][1]*b[i][1]+c[i][1]*c[i][1])*d;a[second[i]][third[i]]=a[second[i]][third[i]]+(b[i][1]*b[i][2]+c[i][1]*c[i][2])*d;a[third[i]][first[i]]=a[third[i]][first[i]]+(b[i][0]*b[i][2]+c[i][0]*c[i][2])*d;a[third[i]][second[i]]=a[third[i]][second[i]]+(b[i][1]*b[i][2]+c[i][1]*c[i][2])*d;a[third[i]][third[i]]=a[third[i]][third[i]]+(b[i][2]*b[i][2]+c[i][2]*c[i][2])*d;}//至此左端刚度矩阵已经生成完毕，接下来根据本质边界条件将其化简，并生成化简后的右端项ffor(i=0;i<n;i++){if((i>=0&&i<37)||(i>=45&&i<91)){for(j=0;j<n;j++){if(i!=j)a[i][j]=0;elsea[i][j]=1;}f[i]=psi[i];}else{for(j=0;j<n;j++){f[i]=f[i]-a[i][j]*psi[j];if((j>=0&&j<37)||(j>=45&&j<91))a[i][j]=0;}}//if(i>=0&&i<37)f[i]=0;//保证边界条件的局部是精确的〔不会产生某个小数从而影响计算精度〕}//至此简化完毕，化为一个对称方程组下面开始用LU分解解此线性方程组l=LU(a,n,l);psi=SOVLELU(l,f,n,psi);//L_U(a,f,n);//for(i=0;i<n;i++)// psi[i]=f[i];//至此已解出流函数psi,释放一批空间并且新开一批空间for(i=0;i<n;i++)delete[]l[i];//释放LU分解的过渡矩阵delete[]l;for(i=0;i<n;i++)delete[]a[i];//释放刚度矩阵delete[]a;delete[]f;//释放右端项double*vx=newdouble[n];//vx[i]即第i+1个节点的x方向速度double*vy=newdouble[n];//vy[i]即第i+1个节点的y方向速度double*count=newdouble[n];//count[i]记录计算速度时每个节点上速度加和的次数，以便之后取平均for(i=0;i<n;i++){vx[i]=0;vy[i]=0;count[i]=0;}//开始求流场中各节点的速度doublexx,yy;//xx,yy暂时储存每个单元x方向与y方向的速度for(i=0;i<e;i++)//为节约空间，不储存每个单元的速度，求出每个单元的速度后直接分配到节点上{xx=c[i][0]*psi[first[i]]+c[i][1]*psi[second[i]]+c[i][2]*psi[third[i]];yy=-b[i][0]*psi[first[i]]-b[i][1]*psi[second[i]]-b[i][2]*psi[third[i]];vx[first[i]]=vx[first[i]]+xx;vx[second[i]]=vx[second[i]]+xx;vx[third[i]]=vx[third[i]]+xx;vy[first[i]]=vy[first[i]]+yy;vy[second[i]]=vy[second[i]]+yy;vy[third[i]]=vy[third[i]]+yy;count[first[i]]=count[first[i]]+1;count[second[i]]=count[second[i]]+1;count[third[i]]=count[third[i]]+1;}for(i=0;i<n;i++){vx[i]=vx[i]/count[i];vy[i]=vy[i]/count[i];}//最后求解压力项double*p=newdouble[n];for(i=0;i<n;i++)p[i]=(1-vx[i]*vx[i]-vy[i]*vy[i])/2;//压力项求解完毕//至此已求出每个节点上的速度，接下来释放空间并进行输出delete[]count;for(i=0;i<e;i++)delete[]b[i];delete[]b;for(i=0;i<e;i++)delete[]c[i];delete[]c;fp=fopen("result111.txt","w");fprintf(fp,"TITLE=Finiteelementmethod\n");fprintf(fp,"FILETYPE=FULL\n");fprintf(fp,"VARIABLES=X,Y,vx,vy,psi,p\n");fprintf(fp,"ZONE\n");fprintf(fp,"N=%dE=%df=fepointet=triangle\n",n,e);for(i=0;i<n;i++)fprintf(fp,"%lf%lf%lf%lf%lf%lf\n",x[i],y[i],vx[i],vy[i],psi[i],p[i]);for(i=0;i<e;i++)fprintf(fp,"%d%d%d\n",first[i]+1,second[i]+1,third[i]+1);delete[]first;delete[]second;delete[]third;delete[]p;delete[]vx;delete[]vy;delete[]x;delete[]y;delete[]psi;fclose(fp);}double**LU(double**a,intm,double**l)//LU分解函数,**a即二维数组〔矩阵〕a的首地址,此函数返回LU矩阵的首地址，m为矩阵规模{//l[][]即LU分解时产生的矩阵(为节约空间，将L和U放在同一矩阵中)inti,j,k;for(k=0;k<m;k++){for(j=k;j<m;j++){l[k][j]=a[k][j];for(i=0;i<k;i++){l[k][j]=l[k][j]-l[k][i]*l[i][j];}}for(j=k+1;j<m;j++){l[j][k]=a[j][k];for(i=0;i<k;i++){l[j][k]=l[j][k]-l[j][i]*l[i][k];}l[j][k]=l[j][k]/l[k][k];}}return(l);}double*SOVLELU(double**l,double*b,intm,double*x)//解LU分解后的方程组，l为LU分解产生的矩阵，b为右端项，m为矩阵规模，x为解向量,返回解向量的首地址{inti,j;//先解Ly=bx[0]=b[0];for(i=1;i<m;i++){x[i]=b[i];for(j=0;j<i;j++)x[i]=x[i]-l[i][j]*x[j];}//再解Ux=yx[m-1]=x[m-1]/l[m-1][m-1];for(i=m-2;i>=0;i--){for(j=i+1;j<m;j++){x[i]=x[i]-l[i][j]*x[j];}x[i]=x[i]/l[i][i];}return(x);}六．可视化处理与结果分析。 绘制网格如下：流函数如下：流线图如下：假设将流线和流函数画在同一张图上，那么为压强分布如下：由图像可见，速度分布比拟符合物理直观，将流函数图和流线图画在同一张图像上，可以发现等ψ线就是流线，满足物理要求。压强图像上，〔-1,0〕处速度为零，压强最大，为滞止压强；〔0,1〕处压强最小，也符合物理要求。七．进一步优化的些许