fortran理想流体的平面圆柱绕流程序.docxVIP

扬州大学 理想流体的平面圆柱绕流 计算水力学 题目： 用 Fortran 语言编写程序解决理想流体的平面圆柱绕流问题，如下图所示。由 于流动的对称性，可以只研究其中的四分之一区域，如图中abcde 所示。 y e d c 2 a r =1 ux=1 b x 2 3.5 3.5 在理想流体的平面运动中，流函数 ψ和势函数 Φ均满足拉氏方程： 2 2 2 2 x2 y2 0， 0 x2 y 2 其边界条件如下表所示。 流函数 势函数 在下边界 abc上 0 0 n 在出口边界 cd 上 0 const n 在上边界 de上 2 0 n 在进口边界 ea 上 y y ux 1 n 说明： 是切向流速， 是法向流速。 n n 下面就流函数进行讨论，为便于分析，把边界条件写成： 在 1 上 其中： 1 为具有本质 B、C 的边界 0 在 2 上 2 为具有自然 B、C 的边界 n 水研 2010级 第 1 页 共 15 页 水利学院 扬州大学 理想流体的平面圆柱绕流 计算水力学 解题步骤： (1)写出Галёркин积分表达式 2 2 x 2 dxdy 0 y 2 通过分部积分，可得： x xy y dxdyg ds 2 (2)区域剖分 横向剖分数为 9，纵向剖分数为 10，其中圆弧段剖分数为 5。利用作业三中的程 序实现 (由于网格内要画流速矢量图， 故单元编号未写出 )，另外，还需要建立本质 B.C 表。 (3)确定单元基函数 e i 设网格划分后任意三角形单元的三个结点的坐标值别为( x e , y e )(i 1,2,3) ，函数 i i 值分别为 i e (i 1,2,3) ，根据基函数的构造思想，单元内近似函数可表示为式： ( e) i e i e (i 1,2,3) 。 在单元内作线性插值函数如下： e a b x c y； e a b x c y； e a b x c y 1 1 1 1 2 2 2 2 3 3 3 3 根据基函数的插值条件，得到系数： ai ,bi ,ci (i 1,2,3) 。 则基函数为： e ai bi x ci y ， i 1,2,3 。 i (4)单元分析 e e e 代入Галёркин积分表达式： 将 i i e e e e e ds e x x y y dxdy 2 g 得单元有限元方程组为： Aij e je Fi e (i=1,2,3) 水研 2010级 第 2 页 共 15 页 水利学院 扬州大学 由  理想流体的平面圆柱绕流 计算水力学 e ai bi x ci y (i=1,2,3)，可得： i e e e e j bj , i bi , j c j , i ci x x y y 于是： b1b1 c1c1 b1b2 c1c2 b1b3 c1c3 Aij A e b2b1 c2c1 b2b2 c2c2 b2 b3 e c2 c3 ， Fi b3b1 c3c1 b3b2 c3c2 b3b3 c3c3  e e g i ds 2 自然 B.C 处理： 由于自然边界条件 0 ，则 Fi e 0 。 n 2 (5)总体合成 单元矩阵 Aij e 总体矩阵 A nm ； 单元矩阵行号 整体矩阵行号，单元矩阵列号 整体矩阵列号。 (6)本质 B.C 处理 即为了满足本质 B.C，要对总体系数矩阵进行处理，具体处理方法见作业二。 (7)解总体方程组，求出有关物理量 解方程组的方法见作业一，由 i ai bi x ci y 及 i i 得： vx y  i ci i , vy i y x  i i bi i ;( i 1,2,3) x 三结点三角形单元， 线性插值函数，每个单元只有一个流速， 与单元内坐标无关， 可理解为单元平均流速，位于单元中心 (三中线交点 )。 源程序如下： 说明：程序的部分说明作业一、二、三中已有，这里不再赘述；其中绘图子程序 在作业三中也已有，这里略去。 program yzrl implicit none interface subroutine linear_equation_bc(n,a1,b,x_result) 水研 2010级 第 3 页 共 15 页 水利学院 扬州大学 理想流体的平面圆柱绕流 计算水力学 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(