fortran 理想流体的平面圆柱绕流程序
用Fluent计算二维圆柱绕流
边界层网格（一）
边界层网格：mesh/boundary layer/create boundary layer 边界层网格设置有两种方法：uniform or aspect ratio based
Uniform包含四个参数： first row（a）指定第一层边界层的厚度 Growth factor（b/a）边界层厚度增长的比例，但如果相邻 边的节点分布已经确定，则网格会自动调整 Rows边界层层数 Depth（D）总的边界层厚度 这四个参数中任意设定三个，则程序会自动算出第四个参 数的值
使用Gambit生成网格
确定几何形状
点 ——> 直线、曲线 —封闭—> 面（特殊面） 布尔运算(Unit, Subtract, Intersect)，移动和拷贝(Move/Copy) 分裂与合并(Split, Merge)，连接与解除连接(connect, disconnect)
生成计算网格
线网格 ——>（边界层网格） ——> 面网格（结构、非结构） 单元形式：三角形单元、四边形单元、混合单元
网格类型
结构网格（structured grid ）
节点排列有序、邻点间的关系明确
非结构网格（unstructured grid）
节点位置无法用一个固定的法则排序 生成过程复杂，但有极好的适应性
Gambit网格生成
结构网格 —— Map 块结构网格 —— Submap 非结构网格 —— Pave
用Fluent计算二维圆柱绕流
王吉飞 wangjifei@
主要内容
计算流体力学简介
Fluent软件简介 二维圆柱绕流标准算例
fortran 理想流体的平面圆柱绕流程序
题目：用Fortran 语言编写程序解决理想流体的平面圆柱绕流问题，如下图所示。
由于流动的对称性，可以只研究其中的四分之一区域，如图中abcde 所示。
u x =1在理想流体的平面运动中，流函数ψ和势函数Φ均满足拉氏方程：02222=∂∂+∂∂y x ψψ，02222=∂∂+∂∂yx φφ 其边界条件如下表所示。
说明：n∂是切向流速， n ∂是法向流速。
下面就流函数进行讨论，为便于分析，把边界条件写成：ψψ~= 在1Γ上 其中：1Γ为具有本质B 、C 的边界 0=∂∂nψ在2Γ上 2Γ为具有自然B 、C 的边界解题步骤：(1)写出Галёркин积分表达式02222=⎪⎪⎭⎫⎝⎛∂∂+∂∂⎰⎰Ωdxdy y x δψψψ 通过分部积分，可得：⎰⎰⎰ΓΩ=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂2ds g dxdy y y x x δψδψψδψψ (2)区域剖分横向剖分数为9，纵向剖分数为10，其中圆弧段剖分数为5。
利用作业三中的程序实现(由于网格内要画流速矢量图，故单元编号未写出)，另外，还需要建立本质B.C 表。
(3)确定单元基函数()e i ϕ设网格划分后任意三角形单元的三个结点的坐标值别为()())3,2,1)(,(=i y x e i e i ，函数值分别为()(1,2,3)e i i ψ=，根据基函数的构造思想，单元内近似函数可表示为式：()())3,2,1()(==i e i e i e ϕψψ。
在单元内作线性插值函数如下：()()()111122223333e e e a b x c y a b x c y a b x c y φφφ=++=++=++；；根据基函数的插值条件，得到系数：,,(1,2,3)i i i a b c i =。
则基函数为：()y c x b a i i i e i ++=ϕ，()3,2,1=i 。
(4)单元分析 将()()()e i e i e ϕψψ=代入Галёркин积分表达式：()()()()()()⎰⎰⎰ΩΓ=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂2ds g dxdy y y x xe e e e e e δψδψψδψψ 得单元有限元方程组为：()()()e e eij j i A F ψ=(i=1,2,3)由()e i i i i a b x c y φ=++(i=1,2,3)，可得：()()()(),,,e e e e j j i i j i j i b b c c x x y yφφφφ∂∂∂∂====∂∂∂∂ 于是：()111112121313212122222323313132323333e ij b b c c b b c c b b c c A A b b c c b b c c b b c c b b c c b b c c b b c c +++⎡⎤⎢⎥=+++⎢⎥⎢⎥+++⎣⎦，()()()⎰=2e ds g F e i e i τϕ自然B.C 处理： 由于自然边界条件02=∂∂Γnψ，则()0e i F =。
FLUENT系列资料3之直排圆柱绕流
圆柱绕流（顺排）1、流动问题如图给出了圆柱绕流的计算区域的几何尺寸, 其中L=1.2m，W=0.5m，r=0.02m，l2=0.1m，l1=0.2m，入口速度为0.01m/s，圆柱横向间隔0.1m，竖向间隔为0.1m，分析不同排列方式的流场情况。
本例涉及到：一、利用GAMBIT建立圆柱绕流的计算模型（1）在CAD中画出圆柱绕流的图形（2）将CAD图形输出为*.sat的文件格式（3）用GAMBIT读入上面输出的*.sat文件（4）对各条边定义网格节点的分布（5）在面上创建网格（6）定义边界内型（7）为FLUENT5/6输出网格文件二、利用FLUENT-2D求解器进行求解（1）读入网格文件（2）确定长度单位：MM（3）确定流体材料及其物理属性（4）确定边界类型（5）计算初始化并设置监视器（6）使用非耦合、隐式求解器求解（7）利用图形显示方法观察流场与温度场一、前处理——用CAD画出圆柱绕流的结构图并导入GAMBIT中在CAD中按所给的尺寸画出圆柱绕流的结构图，画完后输出为drawing1.sat的文件（如图1所示）。
CAD中的操作：文件→输出…点击保存到你想保存到的文件夹中。
圆柱绕流的结构图图1 CAD保存为sat格式的文件启动GAMBIT ，建立一个新的GAMBIT文件。
第1步：确定求解器选择用于进行CFD计算的求解器。
操作：Solver→FLUENT5/6第2步：导入圆柱绕流的结构图操作：File→Import→ACIS…点击Browse找到刚才从CAD中输出的drawing1.sat文件，选中后点击Accept即可导入所需的图形。
（需在CAD中将所画的图形创建成面域，否则无法读入）第3步：确定边界线的内部节点分布并创建结构化网络1、创建各条边上的节点分布操作：MESH→EDGE→打开的“MESH Edges”对话框如图2所示。
（1）点击Edges右侧的黄色区域，使其处于活动状态；（2）Shift+鼠标左键，点击所需划分的边线；（3）选择Interval size，并输入值10；（4）点击Apply，生成各条边上的节点分布。
流体力学Fluent报告——圆柱绕流
亚临界雷诺数下串列单圆柱与圆柱绕流的数值模拟之阳早格格创做目要：原文使用Fluent硬件中的RNG k-ε模型对付亚临界雷诺数下二维串列圆柱战圆柱绕流问题举止了数值钻研，通过截止对付比，分解了雷诺数、柱体形状对付柱体绕流阻力、降力以及涡脱频次的效率.普遍而止，Re数越大，圆柱的阻力越大，圆柱体则可则；而Re越大，二种柱体的降力均越大.相对付于圆柱，共种条件下，圆柱受到的阻力要大；好异天，圆柱涡脱降频次要小.Re越大，串列柱体的Sr数越靠近于单圆柱体的Sr数.闭键字：圆柱绕流、降力系数、阻力系数、斯特劳哈我数正在工程试验中，如航空、航天、航海、体育疏通、风工程及大天接通等广大的本量范畴中，绕流钻研正在工程本量中具备要害的意思.当流体流过圆柱时, 由于漩涡脱降，正在圆柱体上爆收接变效率力.那种效率力引起柱体的振荡及资料的疲倦，益坏结构，成果宽沉.果此，近些年去，稠稀博家战教者对付于圆柱绕流问题举止过细致的钻研，特天是圆柱所受阻力、降力战涡脱降以及涡致振荡问题.沈坐龙等[1]鉴于RNG k⁃ε模型，采与有限体积法钻研了亚临界雷诺数下二维圆柱战圆柱绕流数值模拟，得到了圆柱战圆柱绕流阻力系数Cd与Strouhal 数随雷诺数的变更顺序.姚熊明等[2]采与估计流体硬件CFX中LES模型估计了二维不可压缩匀称流中孤坐圆柱及串列单圆柱的火能源个性.使用非结构化网格六里体单元战有限体积法对付二维N- S圆程举止供解.他们着沉钻研了下雷诺数时串列单圆柱正在分歧间距比时的压力分集、阻力、降力及Sr数随Re数的变更趋势.费宝玲等[3]用FLUENT硬件对付串列圆柱绕流举止了二维模拟，他们采用间距比L/D(L为二圆柱核心间的距离，D为圆柱直径)2、3、4共3个间距举止了数值分解.估计均正在Re = 200 的非定常条件下举止.估计了圆柱的降阻力系数、尾涡脱降频次等形貌绕流问题的主要参量，分解了分歧间距对付圆柱间相互效率战尾流个性的效率.圆柱绕流的一个要害个性是震动形态与决于雷诺数.Lienhard[4]归纳了洪量的真验钻研截止并给出了圆柱体尾流形态随雷诺数变更的顺序.当Re<5时，圆柱上下游的流线呈对付称分集,流体本去不摆脱圆柱体，不旋涡爆收.此时与理念流体相似，若改变流背，上下游流形仍相共.当5<Re<40时,鸿沟层爆收分散,分散剪切层正在圆柱体里前产死一对付宁静的“附着涡”.当40<Re<150时，震动脆持层流状态而且流体旋涡接替天从圆柱后部做周期性的脱降并正在尾流中产死二列接叉排列的涡，即卡门涡街.从150<Re<300启初，旋涡里里启初由层流背湍流转捩，直至减少至3x105安排，此时圆柱体表面附近的鸿沟层仍为层流，所有涡街渐渐转化成湍流，及e<3xl05称为亚临界天区.当3xl05<Re<3.5x106时,鸿沟层的震动也渐渐趋于湍流状态，尾流中不明隐的涡街结构，称为临界状态.[5]圆柱绕流的另一个隐著个性是斯特劳哈我数是雷诺数的函数.早正在1878年，捷克科教家Strouhal[6]便对付风吹过金属丝时收出鸣喊声做过钻研，创造金属丝的风鸣音调与风速成正比,共时与弦线之细细成反比，并提出估计涡脱降频次f的体味公式:式中即斯特劳哈我数Sr由Re所唯一决定.原文使用Fluent硬件中的RNG k-ε模型对付亚临界雷诺数下二维串列圆柱战圆柱绕流问题举止了数值钻研，通过截止对付比，分解了雷诺数、柱体形状对付柱体绕流阻力、降力以及涡脱频次的效率.1.数教模型1.1统造圆程对付于停止圆柱绕流，原文钻研对付象为二维不可压缩震动.正在直角坐标系下，其疏通顺序可用N-S圆程去形貌，连绝性圆程战动量圆程分别为:其中ui为速度分量；p为压力；ρ为流体的稀度；ν为流体的能源黏性系数.对付于湍流情况，原文采与RNG k⁃ε模型，RNG k⁃ε模型是k⁃ε模型的矫正规划.通过正在大尺度疏通战建正后的粘度项体现小尺度的效率，而使那些小尺度疏通有系统天从统造圆程中去除.所得到的k圆程战ε圆程，与尺度k⁃ε模型非常相似，其表白式如下：其中Gk为由于仄衡速度梯度引起的湍动能的爆收项，，，体味常数=0.084 5，==1.39，=1.68.相对付于尺度k⁃ε模型，RNG k⁃ε模型通过建正湍动粘度，思量了仄衡震动中的转动及转动震动情况，RNG k⁃ε模型不妨更佳的处理下应变率及流线蜿蜒程度较大的震动.1.2相闭参数圆柱绕流的相闭参数主要有雷诺数Re、斯特劳哈我数Sr、降力系数Cl战阻力系数Cd，底下给出各个参数的估计公式战物理意思.雷诺数Re与圆柱绕流的状态战雷诺数有很大闭系，雷诺数代表惯性力战粘性力之比：其中U为去流速度；L为个性少度，原文与圆柱直径或者圆柱边少；为流体稀度；、分别为流体介量能源粘度战疏通粘度.斯特劳哈我数Sr是Strouhal 指出圆柱绕流后正在圆柱后里不妨出现接替脱降的旋涡，旋涡脱降频次、风速、圆柱直径之间存留一个闭系：式中：Sr为斯托罗哈数，与决于结构的形状断里；f 为旋涡脱降频次；L为结构的个性尺寸； U 为去流速度.阻力系数战降力系数是表征柱体阻力、降力的无量目参数.定义为：，式中ρ为流体稀度；V为去流速度；A为迎流截里里积；战.由于涡脱降的闭系，阻力系数将爆收振荡，原文采用仄衡脉动降力去钻研，即与圆均根值去钻研.2.数值估计2.1物理模型二维数值模拟单圆柱流场估计天区的采用如图1所示，圆柱绕流以圆柱体直径为个性尺度D，采用圆柱半径为1.5 mm，估计天区为9D×32D的矩形天区.柱1距上游少度图 1 串列圆柱战圆柱的估计天区5D，下游少度27D，脆持二柱间距 L/D= 2. 5D稳定 (L是二圆柱核心连线少度)，二柱到上下鸿沟距离相等.对付于圆柱绕流，采用圆柱边少为个性少度，D=30mm.2.2网格区分估计天区采与分块结构化网格，柱体表面网格干加稀处理，鸿沟区网格相对付稠稀.简直网格区分情况睹图2.其中串列圆柱网格31116个节面，30615个四边形里单元；串图 2 圆柱绕流与圆柱绕流估计域的网格区分列圆柱46446个节面，46550个四边形里单元.2.3鸿沟条件管讲壁里战柱体表面均采与无滑移的停止壁里条件.而出心采用速度出心，出心采用自由出流.去溜速度大小根据Re去树坐，雷诺数分300、3000、12000、30000四个等第，速度大小依次为0.1m/s、1m/s、4m/s、10m/s.2.4估计模型原文湍流模型采与尺度壁里函数的RNG k-ε模型.采与有限容积法供解二维不可压缩粘性流体非定常震动统造圆程，即把估计天区分成很圆柱近壁里网格多小的统造体，对付每个统造体的各个变量举止积分.统造圆程的对付流项采与二阶迎风圆法失集，速度战压力采与SIMPLE算法耦合供解，将所有天区瞅成一个完全举止耦合估计.动量、湍动能战湍动耗集率均采与二阶迎风圆法.先定常估计流场，再用定常估计的截止动做非定常迭代的初初值举止估计.根据初略估计的涡脱频次，牢固树坐时间步少为0. 002s, 正在每个时间步内树坐迭代次数为20.流体介量为液态火.3.估计截止3.1网格模型考证为考证网格独力性，原文估计了网格节面数为8346，里单元为8932的细网格、节面数为31116，里单元为30615的稀网格、节面数为63432，里单元为67434的细稀网格下Re=200、L/D=2的串列网格的Sr数，截止隐现三套网格的估计截止分别为0.143、0.133、0.133.故稀网格可用.而圆柱绕流则采与共级别网格.[7]的估计数据相比较，比较图像如图3所示，最大缺面为2.2%.图3串列圆柱分歧间距的Sr数估计对付比3.2流线与涡量图图 6 Re=3000圆柱绕流流线图图 7 Re=3000圆柱绕流涡量等值线图图 4 Re=3000圆柱绕流流线图图 5 Re=3000圆柱绕流涡量等值线图原文给出了估计历程中雷诺数Re=3000，t=1s时的流线图战涡量图.3.3阻力系数图 9 Re=3000圆柱绕流脉动阻力系数图 8 Re=3000圆柱绕流脉动阻力系数原文给出了Re=3000时，圆柱绕流战圆柱绕流的脉动阻力系数图如下.由图9战错误!未找到引用源。
流体力学Fluent报告——圆柱绕流
亚临界雷诺数下串列双圆柱与方柱绕流的数值模拟摘要：本文运用Fｌｕｅｎｔ软件中的RNＧk-ε模型对亚临界雷诺数下二维串列圆柱和方柱绕流问题进行了数值研究，通过结果对比,分析了雷诺数、柱体形状对柱体绕流阻力、升力以及涡脱频率的影响。
一般而言,Rｅ数越大,方柱的阻力越大，圆柱体则不然；而Rｅ越大,两种柱体的升力均越大。
相对于圆柱,同种条件下，方柱受到的阻力要大；相反地，方柱涡脱落频率要小。
Ｒe越大，串列柱体的Sr数越接近于单圆柱体的Ｓr数。
关键字：圆柱绕流、升力系数、阻力系数、斯特劳哈尔数在工程实践中，如航空、航天、航海、体育运动、风工程及地面交通等广泛的实际领域中,绕流研究在工程实际中具有重大的意义。
当流体流过圆柱时, 由于漩涡脱落，在圆柱体上产生交变作用力。
这种作用力引起柱体的振动及材料的疲劳,损坏结构，后果严重。
因此,近些年来，众多专家和学者对于圆柱绕流问题进行过细致的研究，特别是圆柱所受阻力、升力和涡脱落以及涡致振动问题。
沈立龙等[1］基于RNG ｋ⁃ε模型,采用有限体积法研究了亚临界雷诺数下二维圆柱和方柱绕流数值模拟,得到了圆柱和方柱绕流阻力系数Ｃ与Ｓｔｒｏｄuhal 数随雷诺数的变化规律。
姚熊亮等［２］采用计算流体软件CFX中LES模型计算了二维不可压缩均匀流中孤立圆柱及串列双圆柱的水动力特性。
使用非结构化网格六面体单元和有限体积法对二维Ｎ-S方程进行求解。
他们着重研究了高雷诺数时串列双圆柱在不同间距比时的压力分布、阻力、升力及Sｒ数随Re数的变化趋势。
费宝玲等［3]用ＦＬUENT软件对串列圆柱绕流进行了二维模拟，他们选取间距比L/Ｄ(Ｌ为两圆柱中心间的距离，D为圆柱直径)2、3、４共３个间距进行了数值分析。
计算均在Re= ２0０的非定常条件下进行。
计算了圆柱的升阻力系数、尾涡脱落频率等描述绕流问题的主要参量,分析了不同间距对圆柱间相互作用和尾流特征的影响。
圆柱绕流的一个重要特征是流动形态取决于雷诺数。
圆柱绕流问题
圆柱绕流问题
圆柱绕流分析
⼀、总体步骤
1、在gambit中创建⼏个模型
2、划分⽹格
3、指定边界条件
4、建⽴fluent(Re=1000 和Re=1000000)
5、求解
6、结果分析
7、完善⽹格
⼆、详细步骤
1、建⽴⼏何模型，圆柱的半径为1，左端半圆边界为10，影响半径为4，后边界离圆柱中⼼距离为40，⼏何模型如下：
建⽴完⼏何模型后还要建⽴⾯，划分区域
2、划分直线段后如下图
3、进⾏⽹格制作后如下图
4、边界描述如下
Edges Name
A farfield 1
B,C farfield 2
D,E farfield 3
F,G,H farfield 4
I,J cylinder
5打开fluent,核对数据
6、define solver properties
7、定义材料性质
(1）Re=1000
Density=500 viscosity=1
(2)Re=100000
Density=500000 viscosity=1
8、操作条件：雷诺数100和1000000相同
9、边界条件
10、求解——建⽴解决控制机构——设⽴初始条件Patch Region
11、
（Re=1000）
12、建⽴收敛度
13、建⽴参考值
其中1000雷诺数的与1000000雷诺数的密度不同，流体流速和圆柱直径等相同，下图为设置⼤雷诺数的截图
14、计算机计算
（Re=1000)
(Re=1000000)
压强外形（Re=1000)
(Re=1000000)
涡旋外形
Re=1000
Re=1000000。
基于fluent的圆柱绕流计算分析
注意，计算圆柱体绕流流动可能需要考虑很多因素，例如流体的流动特性、圆柱体的尺寸和形状、流动条件等。因此，在使用Fluent进行计算分析时，你需要仔细设置模型并调整相应的参数，以得到准确的结果。
在计算圆柱体绕流流动时，你还可以使用Fluent的多相流动模型来考虑流体中的汽液相变过程。
多相流动模型可以用来求解含有液体和气体两相的流体流动。例如，你可以使用多相流动模型来计算圆柱体绕流流动中气体汽化的情况。
在使用多相流动模型时，你需要设定相关的物理属性，例如汽液平衡条件、汽化热和蒸发热等。你还需要设定多相流动的边界条件，例如液体的流入流速或气体的流出压力等。
无论使用哪种方法，都需要仔细设置模型并调整相应的参数，以得到准确的结果。在使用Fluent进行计算分析时，你还可以使用各种可视化工具来帮助你理
在Fluent中，你可以使用多种方法来计算圆柱体绕流流动。无论使用哪种方法，都需要仔细设置模型并调整相应的参数，以得到准确的结果。
在设置模型时，你需要考虑圆柱体的尺寸和形状，以及流体的流动特性。你还需要设定流动条件，例如流速或压力。
一种常用的方法是使用网格流动模型，这种方法可以用来求解流体的流动特性，例如流速、压力和温度等变量的分布情况。在使用网格流动模型时你需要在圆柱体的外围创建一个网格，并在圆柱体的内部创建一个流动区域。然后，你可以设定流动条件，例如流速或压力，并使用Fluent的解算器来求解这个流动模型。
另一种方法是使用非网格流动模型，这种方法可以用来求解流体的运动轨迹和流动特性。在使用非网格流动模型时，你需要在圆柱体的外围创建一组流动粒子，并设定运动轨迹的初始条件。然后，你可以使用Fluent的解算器来求解这个流动模型，得到流动粒子的运动轨迹和流动特性。
最新FLUENT仿真计算不同雷诺数下的圆柱绕流
FLUENT仿真计算不同雷诺数下的圆柱绕流。
尾迹与旋涡脱落经典图如下：Re=1 无分离流动Re=20 尾流中一对稳定的弗普尔旋涡Re=100 圆柱后方形成有规律的涡街Re=3900Re=100000 随着Reynolds数增大，涡道内部向湍流过度，直到全部成为湍流Re=1000000 超临界区，分离点后移，尾流变窄，涡道凌乱，涡随机脱落Re=10000000 极超临界区，分离点继续后移，尾流变窄，湍流涡道重新建立。
图3 Cd随Re的变化曲线图3中实曲线是由Wieselsberger,A.Roshko以及G.W.Jones和J.J.Walker测量数据绘制得到，图中圆点部分是FLUENT计算值在Re=106（超临界区），从经典数据和我们的计算结果都可以看到，圆柱体的平均阻力系数急剧下降。
这是因为在Re=3×105附近，边界层流动由层流状态转变为湍流状态，虽然湍流边界层流动的摩擦阻力较层流边界层大，但它从物面的分离较晚，所以形成较小的尾流区。
由于钝体绕流的阻力主要是压差阻力，所以此时物体的总阻力有了一个明显的下降。
入口VELOCITY_INLET,出口OUTFLOW，上下WALL.Re=1,20,100,二维层流模型。
Re=3900后，三维大涡模型计算不准与网格划分与一些参数设置有关。
1。
圆柱中心离上下边界（wall)的距离大于10D(D为圆柱直径），影响较小。
2。
湍流模型采用大涡模型（LES)。
是目前最复杂，最完善的一种湍流模型。
试验曲线来自，《Boundary-Layer Theory》, Dr.HERMANN SCHLICHTING, Translated by Dr.J.KESTIN,Seventh Edition，用MATLB绘制4.阻力系数的求法请参考此论坛我发的教程FLUENT三分立系数的求法人民法院--庭审流程图书记员首先入庭【站立】：现宣布法庭纪律书记员【宣告完毕】：全体起立，请审判长审判员入庭【审判长一行依次入庭就坐】书记员【清点当事人及其代理人】：报告审判长，本案当事人及诉讼代理人已全部到庭，请开庭。
利用FLUENT 软件模拟圆柱绕流
利用FLUENT软件模拟圆柱绕流摘要:使用计算流体力学软件FLUENT,模拟均匀来流绕固定圆柱的流动,模拟雷诺数为20 ,40 ,100 时的绕流流动,得到流场的流函数等值线图和速度矢量图。
计算结果表明:当雷诺数增加时,流动表现出一系列不同的构造。
在雷诺数约为40 前后流场有明显变化。
小于这个数时,存在一对位置固定的旋涡。
大于40 时,流场开始变得不稳定,旋涡扩大、脱落、又生成,逐渐发展成两排周期性摆动和交错的旋涡。
并与实验及数值模拟结果比较,确认FLUENT能够很好地预测流动结构。
关键词: 圆柱绕流 FLUENT软件雷诺数1圆柱绕流理论分析研究的状况一个世纪以来,圆柱绕流一直是众多理论分析、实验研究及数值模拟对象。
但迄今对该流动现象物理本质的理解仍是不完整的。
圆柱绕流中,起决定作用的是雷诺数,但还受到许多因素,如阻塞比,来流湍流度,下游边界条件等的影响。
随着雷诺数的增加,粘性不可压缩流体绕圆柱的流动会呈现各种不同的流动状态,在小雷诺数时,流动是定常的,随着雷诺数的增加,圆柱后会出现一对尾涡。
当雷诺数较大时,尾流首先失稳,出现周期性的振荡。
而后附着涡交替脱落, 泻入尾流形成Karman 涡街,随着雷诺数的增加,流动变得越来越复杂,最后发展为湍流。
White[1]认为“圆柱涡流具有经典性的重要意义”。
一般认为圆柱绕流有2 种定常的流动图案:雷诺数为较小时,圆柱后无尾涡;当雷诺数为较大时,圆柱后有一对对称的尾涡。
关于定常流失稳以及出现湍流的临界雷诺数主要是通过应用流场显示技术观察流动形态得到的,所以不是准确值。
对于分界点雷诺数就有不同的见解,Kovasznay 、Roshko 等认为定常流动失稳的临界雷诺数大约为40。
而从周期性尾流到湍流的详细的转变过程的实验研究似乎还是空白。
对均匀来流绕固定圆柱的二维平面流动,国内外许多学者进行过大量的研究。
决定圆柱绕流流态的是雷诺数(Re)的值,Re5<时,流动不发生分离,5Re40<<,在圆柱体后面出现一对位置固定的旋涡; 40Re150<< ,旋涡扩大,然后有一个旋涡开始脱落,接着另一个也脱落,在圆柱体后面又生成新的旋涡,这样逐渐发展成两排周期性摆动和交错的旋涡,即Karman 涡街。
圆柱绕流方程编程说明
15, j
1 4
2 14, j 15, j1 15, j1 ,
j 6,7,8
上式的推导参照教材（3-9）式。
3. 用差分方程求解区域内各点的流函数（1Leabharlann 正则内点： i, j
（2）非正则内点：
1 4
yk
x2 y2
K 1
K 1
2
2 2
2
( x1 ( y1
x2 ) y2 )
六、程序框图
设置步长，定义参数
假设流场流函数初始值
写出第一类边界点的流函数; 利用初始值计算第二类边界点的流函数
将流场内的点分类
FALSE
利用边界条 件和差分方 程求解正则 内点的流函
求解：流场内的流函数，画出流线的形状。
三、问题的分析
由于问题的对称性，可以将左上四分之一区域ABCDEA作为 研究区域。对于不可压缩平面势流，流函数满足Laplace方程：
2
x2
2
y 2
0
如下，由于各边界条件已知，可得边界点的差分方程，利用 边界点的流函数计算得到流场内部各点的流函数。
a、b为交点1，2与结点（i,j）的距离，
如图所示： a (15 i)h 1 j 1h2
b ( j 1)h 1 15 ih2
a. 当 a 0 或 b 0 时，结点(i, j)在圆柱内或表面，
i, j 0
b. 当 0 a h,0 b h 时， 按公式正常求解（注意：本程序中仅有一个结点符合此条件， 不要漏掉）。
数
利用边界条 件和差分方 程求解非正 则内点的流
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
题目：用Fortran 语言编写程序解决理想流体的平面圆柱绕流问题，如下图所示。
由于流动的对称性，可以只研究其中的四分之一区域，如图中abcde 所示。
u x =1在理想流体的平面运动中，流函数ψ和势函数Φ均满足拉氏方程：02222=∂∂+∂∂y x ψψ，02222=∂∂+∂∂yx φφ 其边界条件如下表所示。
说明：n∂是切向流速， n ∂是法向流速。
下面就流函数进行讨论，为便于分析，把边界条件写成：ψψ~= 在1Γ上 其中：1Γ为具有本质B 、C 的边界 0=∂∂nψ在2Γ上 2Γ为具有自然B 、C 的边界解题步骤：(1)写出Галёркин积分表达式02222=⎪⎪⎭⎫⎝⎛∂∂+∂∂⎰⎰Ωdxdy y x δψψψ 通过分部积分，可得：⎰⎰⎰ΓΩ=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂2ds g dxdy y y x x δψδψψδψψ (2)区域剖分横向剖分数为9，纵向剖分数为10，其中圆弧段剖分数为5。
利用作业三中的程序实现(由于网格内要画流速矢量图，故单元编号未写出)，另外，还需要建立本质B.C 表。
(3)确定单元基函数()e i ϕ设网格划分后任意三角形单元的三个结点的坐标值别为()())3,2,1)(,(=i y x e i e i ，函数值分别为()(1,2,3)e i i ψ=，根据基函数的构造思想，单元内近似函数可表示为式：()())3,2,1()(==i e i e i e ϕψψ。
在单元内作线性插值函数如下：()()()111122223333e e e a b x c y a b x c y a b x c y φφφ=++=++=++；；根据基函数的插值条件，得到系数：,,(1,2,3)i i i a b c i =。
则基函数为：()y c x b a i i i e i ++=ϕ，()3,2,1=i 。
(4)单元分析 将()()()e i e i e ϕψψ=代入Галёркин积分表达式：()()()()()()⎰⎰⎰ΩΓ=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂2ds g dxdy y y x xe e e e e e δψδψψδψψ 得单元有限元方程组为：()()()e e eij j i A F ψ=(i=1,2,3)由()e i i i i a b x c y φ=++(i=1,2,3)，可得：()()()(),,,e e e e j j i i j i j i b b c c x x y yφφφφ∂∂∂∂====∂∂∂∂ 于是：()111112121313212122222323313132323333e ij b b c c b b c c b b c c A A b b c c b b c c b b c c b b c c b b c c b b c c +++⎡⎤⎢⎥=+++⎢⎥⎢⎥+++⎣⎦，()()()⎰=2e ds g F e i e i τϕ自然B.C 处理： 由于自然边界条件02=∂∂Γnψ，则()0e i F =。
(5)总体合成单元矩阵()e ij A ⇒总体矩阵nm A ；单元矩阵行号⇒整体矩阵行号，单元矩阵列号⇒整体矩阵列号。
(6)本质B.C 处理即为了满足本质B.C ，要对总体系数矩阵进行处理，具体处理方法见作业二。
(7)解总体方程组，求出有关物理量解方程组的方法见作业一，由 y c x b a i i i i ++=ϕ 及 i i ϕψψ=得：,;(1,2,3)i i x i i i y i i i v c v b i y y x xφφψψψψψψ∂∂∂∂====-=-=-=∂∂∂∂ 三结点三角形单元，线性插值函数，每个单元只有一个流速，与单元内坐标无关，可理解为单元平均流速，位于单元中心(三中线交点)。
源程序如下：说明：程序的部分说明作业一、二、三中已有，这里不再赘述；其中绘图子程序在作业三中也已有，这里略去。
program yzrl implicit none interfacesubroutine linear_equation_bc(n,a1,b,x_result)integer::i,j,k,imaxinteger,intent(in)::nreal::max,creal,dimension(:,:),intent(in)::a1real,dimension(:),intent(in)::breal,dimension(:,:),allocatable::a,mreal,dimension(:),intent(inout)::x_resultend subroutine linear_equation_bcend interfacecharacter*12 file,name,ly*8integer*2 lengthinteger::i,j,k,dy,dyz,jd,jd1,jd2,jdz !定义单元、节点编号integer::m,n,m0integer,dimension(:,:),allocatable::zt !定义单元结点整体编号数组integer::a,b,c,jj,kkreal,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::f !定义结点流函数值real, dimension(:),allocatable::vx,vy,v !定义单元x、y方向流速及合成流速real, dimension(:),allocatable::zx,zy,zxx,zyy !定义单元中心及流速终点坐标real, dimension(:),allocatable::jx !定义合成流速与x方向夹角real, dimension(:),allocatable::jtx1,jty1,jtx2,jty2 !箭头坐标print*,'请输入x方向等分数m，y方向等分数n，圆弧等分数m0'read*,m,n,m0print*,'请输入网格划分部分尺寸x1，y1，r'read*,x1,y1,rdyz=m*n *2 !计算单元总数jdz=(m+1)*(n+1) !计算结点总数allocate(x(jdz),y(jdz),zt(dyz,3))do i=1,m !计算局部节点编号与整体节点编号的关系数组do j=1,2*ndy=(i-1)*2*n+jif(mod(j,2)==0)then !计算偶数单元的节点对应关系数组zt(dy,1)=j/2+n+(n+1)*(i-1)+1zt(dy,2)=j/2+(n+1)*(i-1)+1zt(dy,3)=j/2+n+(n+1)*(i-1)+2else !计算奇数单元的节点对应关系数组zt(dy,1)=j/2+(n+1)*(i-1)+1zt(dy,2)=j/2+(n+1)*(i-1)+2zt(dy,3)=j/2+n+(n+1)*(i-1)+2end ifend doend doopen(1,file='单元中整体结点编号.txt')write(1,'(a,i4,a,3i5)'),('第',i,'个单元:',(zt(i,j),j=1,3),i=1,dyz)close(1)bcx1=x1/m !定义x1边的等分步长bcx2=(x1-r)/(m-m0) !定义x2边的等分步长bcy1=y1/n !定义y1边的等分步长do i=1,m-m0+1do j=1,n+1jd=(n+1)*(i-1)+jx(jd)=(i-1)*bcx1+(i-1)*(j-1)*(bcx2-bcx1)/ny(jd)=y1-(j-1)*bcy1end doend dobcyh=pi/(2.0*m0) !定义圆弧段的等分角大小do i=m-m0+2,m+1jd1=(n+1)*(i-1)+1jd2=(n+1)*ix(jd1)=(i-1)*bcx1 !计算与圆弧段对应的x1边上的等分点的x坐标y(jd1)=y1 !计算与圆弧段对应的x1边上的等分点的y坐标x(jd2)=x1-r*cos((i-(m-m0+1))*bcyh) !计算圆弧段m0等分点的各个x坐标y(jd2)=r*sin((i-(m-m0+1))*bcyh) !计算圆弧段m0等分点的各个y坐标bcx=(x(jd2)-x(jd1))/n !计算该计算区域的x方向的步长bcy=(y(jd2)-y(jd1))/n !计算该计算区域的y方向的步长do j=2,njd=(n+1)*(i-1)+jx(jd)=(i-1)*bcx1+(j-1)*bcxy(jd)=y1+(j-1)*bcyend doend doopen(2,file='整体结点坐标.txt') !将整体结点坐标保存在txt文档中write(2,'(a,i3,a,f8.4,3x,f8.4)'),('第',i,'点坐标为:',x(i),y(i),i=1,jdz)close(2)print*,'请输入图形名称'read(*,'(a)')namelength=len_trim(name) !用内部函数求出name去掉尾部空格后的长度file=name(:length)//'.dxf' !连接字符串"name"与".dxf"open(3,file=file)call staplot !调用绘制dxf文件的开头子程序ly='wg'length=len_trim(ly)x0=0 !给定绘图原点y0=0do i=1,dyza=zt(i,1)b=zt(i,2)c=zt(i,3)call line(x(a),y(a),x(b),y(b),ly(:LENGTH)) !调用绘图子程序画线call line(x(b),y(b),x(c),y(c),ly(:LENGTH))call line(x(c),y(c),x(a),y(a),ly(:LENGTH))end docall textc(x(jdz-n)+0.1,y(jdz-n)+0.05,0.03,file,ly(:length),0.0)print*,'请输入来流速度ux'read*,uxz=2*(m+1)+n-1 !计算具有本质B.C节点个数allocate(bcjd(z),bcz(z))do i=1,z !定义本质B.C,找到具有本质B.C的节点和其对应的值if(i<=n+1)thenbcjd(i)=ibcz(i)=bcy1*(n+1-i)*uxelse if(i<=m+n+1)thenbcjd(i)=(i-n)*(n+1)bcz(i)=0.0elsebcjd(i)=(n+1)*(z-i+1)+1bcz(i)=y1*uxend ifend doopen(4,file='本质B.C表.txt') !建立本质B.C表write(4,'(a,3x,a,3x,a)'),'序号','对应序号整体编号','本质B.C值'do i=1,zwrite(4,'(i3,10x,i3,12x,f6.3)'),i,bcjd(i),bcz(i)end doclose(4)allocate(zta(jdz,jdz),bb(dyz,3),cc(dyz,3))zta=0 !确定单元基函数do i=1,dyzd=x(zt(i,2))*y(zt(i,3))+x(zt(i,3))*y(zt(i,1))+x(zt(i,1))*y(zt(i,2))-x(zt(i,2))*y(zt(i,1)) -x(zt(i,1))*y(zt(i,3))-x(zt(i,3))*y(zt(i,2))dym=d/2.0 !计算单元三角形面积dym，编号均为逆时针，故面积均为正bb(i,1)=(y(zt(i,2))-y(zt(i,3)))/d !计算基函数系数bb(i,2)=(y(zt(i,3))-y(zt(i,1)))/dbb(i,3)=(y(zt(i,1))-y(zt(i,2)))/dcc(i,1)=(x(zt(i,3))-x(zt(i,2)))/dcc(i,2)=(x(zt(i,1))-x(zt(i,3)))/dcc(i,3)=(x(zt(i,2))-x(zt(i,1)))/d!计算各单元系数矩阵，由于自然B.C为0，故Fi(e)为0open(5,file='单元系数矩阵.txt')write(5,'(a,i3)'),'单元',ido j=1,3 !总体系数矩阵的合成jj=zt(i,j)do k=1,3kk=zt(i,k)dya(j,k)=(bb(i,j)*bb(i,k)+cc(i,j)*cc(i,k))*dymzta(jj,kk)=zta(jj,kk)+dya(j,k)end dowrite(5,*),(dya(j,k),k=1,3) !输出单元系数矩阵end doend doclose(5)open(6,file='整体系数矩阵.txt')write(6,*),((zta(i,j),j=1,jdz),i=1,jdz) !输出整体系数矩阵close(6)allocate(cs(jdz),f(jdz))cs=0do k=1,z ! 处理本质B.Ccs(bcjd(k))=bcz(k)do j=1,jdzif(j==bcjd(k))thenzta(j,j)=1elsezta(bcjd(k),j)=0end ifend doend doopen(7,file='本质B.C处理后的整体系数矩阵.txt')write(7,*),((zta(k,j),j=1,jdz),k=1,jdz) !输出本质B.C处理后的整体系数矩阵close(7)open(8,file='方程组右侧矢量项.txt')do k=1,jdz !输出方程组右侧矢量项write(8,'(a,i3,a,f6.3)'),'第',k,'行矢量项',cs(k)end doclose(8)call linear_equation_bc(jdz,zta,cs,f) !调用子程序解线性方程组open(9,file='结点流函数值.txt')do k=1,jdzwrite(9,'(a,i3,a,f10.6)'),'第',k,'点函数值',f(k)end doclose(9)allocate(vx(dyz),vy(dyz),v(dyz))open(10,file='单元x、y向流速.txt')open(11,file='单元流速.txt')do i=1,dyz !计算单元流速vx(i)=cc(i,1)*f(zt(i,1))+cc(i,2)*f(zt(i,2))+cc(i,3)*f(zt(i,3))vy(i)=-(bb(i,1)*f(zt(i,1))+bb(i,2)*f(zt(i,2))+bb(i,3)*f(zt(i,3)))v(i)=sqrt(vx(i)**2+vy(i)**2)write(10,'(a,i3,a,f10.6,2x,a,f10.6)'),'第',i,'单元x向流速',vx(i),'y向流速',vy(i) write(11,'(a,i3,a,f10.6)'),'第',i,'个单元流速',v(i)end doclose(10)close(11)allocate(zx(dyz),zy(dyz),zxx(dyz),zyy(dyz),jx(dyz))allocate (jtx1(dyz),jty1(dyz),jtx2(dyz),jty2(dyz))open(12,file='夹角.txt') !绘制流场矢量图do i=1,dyz !找到各个单元的中心zx(i)=(x(zt(i,1))+x(zt(i,2))+x(zt(i,3)))/3zy(i)=(y(zt(i,1))+y(zt(i,2))+y(zt(i,3)))/3if(vx(i)==0)then !计算流速与x轴的夹角jx(i)=pi/2.0elsejx(i)=atan(vy(i)/vx(i))end ifwrite(12,'(a,i3,a,f10.6)'),'第',i,'单元流速与x轴夹角',jx(i)zxx(i)=zx(i)+cos(jx(i))*v(i)*0.1 !流线终点x坐标zyy(i)=zy(i)+sin(jx(i))*v(i)*0.1 !流线终点y坐标call line(zx(i),zy(i),zxx(i),zyy(i),ly(:length)) !画流线jtx1(i)=zxx(i)-sin(-jx(i)+8*pi/18.0)*v(i)*0.03 !箭头终点x坐标jty1(i)=zyy(i)-cos(-jx(i)+8*pi/18.0)*v(i)*0.03 !箭头终点y坐标call line(zxx(i),zyy(i),jtx1(i),jty1(i),ly(:length)) !画箭头jtx2(i)=zxx(i)-cos(jx(i)-pi/18.0)*v(i)*0.03 !箭头终点x坐标jty2(i)=zyy(i)-sin(jx(i)-pi/18.0)*v(i)*0.03 !箭头终点y坐标call line(zxx(i),zyy(i),jtx2(i),jty2(i),ly(:length)) !画箭头end doclose(12)call endplotclose(3)end program yzrlsubroutine linear_equation_bc(n,a1,b,x_result) !解方程组子程序implicit noneinteger::i,j,k,imaxinteger,intent(in)::nreal::max,creal,dimension(:,:),intent(in)::a1real,dimension(:),intent(in)::breal,dimension(:,:),allocatable::a,mreal,dimension(:),intent(inout)::x_resultallocate(a(n,n+1),m(n,n+1))do i=1,ndo j=1,n+1if(j<=n)thena(i,j)=a1(i,j)elsea(i,j)=b(i)end ifend doend dodo k =1,n-1max=abs(a(k,k))imax=kdo i=k+1,n+1if (abs(a(i,k))>max) thenmax=abs(a(i,k))imax=iend ifend dodo j=k,n+1m(k,j)=a(k,j)a(k,j)=a(imax,j)a(imax,j)=m(k,j)end dodo i=k+1,nm(i,k)=a(i,k)/a(k,k)do j=k+1,n+1a(i,j)=a(i,j)-m(i,k)*a(k,j) end doend doend dox_result(n)= a(n,n+1)/a(n,n)do k=n-1,1,-1c=0do j=n,k+1,-1c=c+a(k,j)*x_result(j)end dox_result(k)=( a(k,n+1)-c)/a(k,k) end doend subroutine算例：输入命令如下图所示：计算可得到如下文件(只列举部分)：(1)本质B.C表：(3)单元流速表：(5)生成平面圆柱绕流如下图所示：。