Fortran90 编写理想流体的平面圆柱绕流程序.doc

题目： 用Fortran语言编写程序解决理想流体的平面圆柱绕流问题，如下图所示。由于流动的对称性，可以只研究其中的四分之一区域，如图中abcde所示。 在理想流体的平面运动中，流函数ψ和势函数Φ均满足拉氏方程： ， 其边界条件如下表所示。 流函数 势函数 在下边界上 在出口边界上 在上边界上 在进口边界上 说明： 是切向流速， 是法向流速。 下面就流函数进行讨论，为便于分析，把边界条件写成： 在上 其中：为具有本质B、C的边界 在上 为具有自然B、C的边界 解题步骤： (1)写出Галёркин积分表达式 通过分部积分，可得： (2)区域剖分 横向剖分数为9，纵向剖分数为10，其中圆弧段剖分数为5。利用作业三中的程序实现(由于网格内要画流速矢量图，故单元编号未写出)，另外，还需要建立本质B.C表。 (3)确定单元基函数 设网格划分后任意三角形单元的三个结点的坐标值别为，函数值分别为，根据基函数的构造思想，单元内近似函数可表示为式：。 在单元内作线性插值函数如下： 根据基函数的插值条件，得到系数：。 则基函数为：，。 (4)单元分析 将代入Галёркин积分表达式： 得单元有限元方程组为： (i=1,2,3) 由(i=1,2,3)，可得： 于是： ， 自然B.C处理： 由于自然边界条件，则。 (5)总体合成 单元矩阵总体矩阵； 单元矩阵行号整体矩阵行号，单元矩阵列号整体矩阵列号。 (6)本质B.C处理 即为了满足本质B.C，要对总体系数矩阵进行处理，具体处理方法见作业二。 (7)解总体方程组，求出有关物理量 解方程组的方法见作业一，由 及 得： 三结点三角形单元，线性插值函数，每个单元只有一个流速，与单元内坐标无关，可理解为单元平均流速，位于单元中心(三中线交点)。 源程序如下： 说明：程序的部分说明作业一、二、三中已有，这里不再赘述；其中绘图子程序在作业三中也已有，这里略去。 program yzrl implicit none interface subroutine linear_equation_bc(n,a1,b,x_result) integer::i,j,k,imax integer,intent(in)::n real::max,c real,dimension(:,:),intent(in)::a1 real,dimension(:),intent(in)::b real,dimension(:,:),allocatable::a,m real,dimension(:),intent(inout)::x_result end subroutine linear_equation_bc end interface character*12 file,name,ly*8 integer*2 length integer::i,j,k,dy,dyz,jd,jd1,jd2,jdz !定义单元、节点编号 integer::m,n,m0 integer,dimension(:,:),allocatable::zt !定义单元结点整体编号数组 integer::a,b,c,jj,kk real,parameter::pi=3.1415926536 !定义pi为常量，值为圆周率 real, dimension(:),allocatable::x,y !定义整体结点坐标数组 real::x1,y1,r !定义网格划分区域 real::bcx1,bcx2,bcy1,bcx,bcy,bcyh !定义x，y方向及圆弧段计算步长 real::x0,y0 !定义原点坐标 real::ux !定义来流速度 integer::z !定义本质B.C点个数 integer, dimension(:),allocatable::bcjd !定义本质结点编号 real, dimension(:),allocatable::bcz !定义本质B.C值 real,dimension(3,3)::dya !定义单元系数矩阵 real,dimension(:,:),allocatable::zta !整体系数矩阵 real,dimension(:,:),allocatable::bb,cc !定义基函数系数矩阵 real::dym,d !定义单元三角形面积 real, dimension(:),allocatable::cs !定义常数项数组 real, dimension(:),allocatable