雷诺数为22000的二维方柱绕流仿真计算

合集下载

雷诺数为22000的二维方柱绕流仿真计算

雷诺数为22000的二维方柱绕流仿真计算

图 15-1-3 网格划分结果
图 15-1-3 给出了网格划分结果,在划分网格时,约束整个面的网格尺寸为 10mm, 并且定义入口、 出口、 上下边界、 方柱表面的名称集 (Named Selections) , 以方便在 Fluent 中定义边界条件。
图 15-1-4 局部网格示意
图 15-1-4 给出了局部网格示意图,从图中可以看出,在方柱表面的网格高 度是较高的, 使用尺度化近壁面函数法时并没有划分很细的边界层网格。最后共 得到 48900 个网格,网格最大长宽比为 1.0069,满足仿真计算的网格要求。
(3)对于存在较高导热率的共轭传热问题; (4)对于流场网格长宽比较大的问题。 在本算例中,虽然不存在以上问题,但仍采用双精度 4 核并行计算求解。 在进行计算时, 先进行稳态计算,然后将稳态计算的结果作为瞬态计算的初 场,再进行瞬态计算。 求解器采用压力基求解器,并进行湍流模型的设置。
图 15-1-6 湍流模型的选择
关于雷诺数 22000 方柱绕流 (关于壁面函数法、湍流粘度比的计算等概念可以参考 文献十六)
15.1.1 RNG 模型与尺度化壁面函数在方柱绕流中的应用
图 15-1-1 分析流程
图 15-1-1 给出了二维方柱绕流问题的分析流程, 在 Workbench 中建立模型, 然后使用 Mesh 模块进行分网,最后提交到 Fluent 中进行计算。
图 15-1-15 方柱表面的 Y+值分布范围
仿真后可以知道,近壁面 Y+值的分布范围为 7.07-68.6,在 Fluent 的官方 文献中推荐壁面函数法的 Y+值适用范围为 30-500,对于尺度化壁面函数,由于 在近壁面使用了预防 Y+值过小导致结果恶化的技术,其 Y+值范围下限可以适当 放宽, 但仍应尽量取在 11.25 以上。 过大或过小的 Y+值都可能引起结果的恶化。

用Fluent计算二维圆柱绕流

用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软件简介 二维圆柱绕流标准算例

高雷诺数下方柱绕流特性的数值模拟

高雷诺数下方柱绕流特性的数值模拟

高雷诺数下方柱绕流特性的数值模拟周强;廖海黎;曹曙阳【摘要】为研究高雷诺数为22000下方柱周围流场形态及气动力特性,基于开源计算流体动力学(computational fluid dynamic,CFD)软件OpenFoam平台,采用基于动态亚格子模型的大涡模拟(large eddy simulation,LES)方法,对均匀来流作用下的方柱绕流迸行了三维数值模拟.首先,通过对基于时间积分的平均积分分量的比较,验证了本数值计算的准确性;其次,深人分析了方柱周围及其尾流区的流场结构,给出了流场结构的平均和湍流特征,并在此基础上,研究了其气动力特性;最后,分析了两种长径比下表面压力的展向空间相关性.研究结果表明:雷诺数为22000下方柱尾流区回转长度为1.37倍方柱宽度,Strouhal数为0.121,脉动升力系数为1.40;展向长度取8倍方柱宽度可更准确地获得周围湍流特性.%Three-dimensional large eddy simulations(LES)were performed at high Reynolds numbers (Re=22 000)and using the dynamic Smagorinsky sub-grid model to investigate uniform flow over a square cylinder. Meshing of the model was performed using Open Foam-an open-source computational tool. Mean integral quantities were compared against existing experimental and numerical results to validate the proposed numerical method. Subsequently,the flow structure around and aerodynamic forces acting on the cylinder were analysed to perform an in-depth investigation of mean-and turbulent-flow characteristics. Lastly,the influence of the grid system around the cylinder and span-wise length on the flow structure and spatial correlation were investigated via comparison between different cases. Results demonstrate that the recirculation length in the wake approximately measured1.37D(width of square),while the Strouhal number and RMS value of the lift coefficient are 0. 121 and 1. 40, respectively. In addition,domain length along the span-wise should be equal to 8D in order to obtain clear turbulent-flow characteristics.【期刊名称】《西南交通大学学报》【年(卷),期】2018(053)003【总页数】7页(P533-539)【关键词】方柱绕流;大涡模拟;流场特性;尾流结构;空间相关性【作者】周强;廖海黎;曹曙阳【作者单位】西南交通大学风工程四川省重点实验室,四川成都610031;西南交通大学风工程四川省重点实验室,四川成都610031;同济大学士木工程防灾国家重点实验室,上海200092【正文语种】中文【中图分类】U441长期以来,方柱绕流一直是风工程、钝体空气动力学、海洋工程等研究领域重要的、基础性的热点问题.在实际工程应用中,如桥梁的桥塔和桥墩、高层建筑、海洋平台等具有方柱形状的构件,在风或者水流作用下,会产生有规律的大尺度旋涡脱落现象,致使结构发生风致振动,进而导致其受损或发生疲劳破坏.因此,研究在方柱绕流问题及其尾流结构对于实际工程设计和建造具有十分重要的意义.不同雷诺数下,流体绕过方柱呈现出明显不同的流动特性.其中在高雷诺数下,流体经方柱前端边角后发生明显的流动分离,并伴随着回流或转捩,形成大尺度和不稳定的旋涡脱落,即与低雷诺数下流场相比,高雷诺数下的方柱绕流特性更加复杂,且表现出更加显著的三维特性.鉴于此,若采用雷诺时均方法,很难扑捉到准确的高雷诺数下方柱绕流特性;由于计算量太大,目前直接数值模拟还只限于雷诺数较低的流动.因此,对于高雷诺数方柱绕流问题,主要的数值研究方法就是大涡模拟(large eddy simulation, LES)方法.到目前为止,国内外学者针对方柱绕流问题开展了大量而丰富的实验和数值模拟研究.在实验研究方面,Lyn等[1]通过水洞实验,采用LDV(laser doppler velocimetry)技术研究了Re=21 400下的方柱绕流特性,详细分析了其尾流结构,同时其研究结果也广泛被作为数值模拟的参考标准;此后国外学者Saha等[2]在风洞中对高雷诺数下方柱的绕流特性进行了深入研究;近年来,Nishimura和Taniike[3]通过试验详细研究了不同攻角下方柱绕流特性;Yen等[4]利用粒子图像测速技术(particle image velocimetry, PIV)研究了方柱在4 000<Re<36 000、攻角0°~45°范围内的绕流结构和气动力特性.在数值研究方面,Sohankar等[5]对不同Re数下方柱绕流采用LES方法进行数值研究,其结果表明二维模拟与三维模拟的结果相差较大,而三维模拟结果明显更为准确;Okajima[6]也得到相类似的结论;Srinivas等[7]依据Lyn的实验,采用LES方法研究了方柱表面压力分布以及尾流区流场特性,发现LES 能够再现实验结果,但也发现数值模拟结果与实验结果在阻力系数、回转长度(回流区长度)等方面仍然存在一定的差别;Wang和Vanke[8]得到了相类似的结论;与此相似,Rodi等[9]总结了相关方面的研究进展,认为大多数数值算例得到的Strouhal 数(St)、流向平均速度均与实验结果接近,但是回转长度、阻力系数、速度脉动值以及表面压力分布相差很大,且与实验结果存在一定差别,其原因众说纷纭[9-11],认为有可能受到网格数量、来流条件、湍流模型参数等因素影响.张伟和葛耀君[12]采用κ-ω SST (shear-stress transport)模型,并与其PIV试验结果进行对比,得到了方柱的阻力系数、St以及流场的速度剖面和湍动能剖面,但是在脉动项上与实验结果还是存在差别.综上可见,对于高雷诺数下方柱绕流的数值模拟,因其流动的复杂性和明显的三维特性,此时,采用LES是比较合适的选择.然而如上文所述,大多采用LES方法模拟中展向长度的选取等细节问题还是存在不够明确的地方,导致一些流场速度、气动力等关键参数的脉动项与实验结果存在一定差别.为此,本文基于开源CFD软件OpenFoam平台,采用基于动态亚格子模型的大涡模拟方法,在Re=22 000情况下,对均匀来流条件下的不同网格系统、展向长度的方柱绕流问题进行了三维数值模拟,首先对数值结果的准确性进行验证,然后重点分析断面周围网格以及长径比对尾流区流场结构的影响,并给出不同长径比下气动力的展向空间相关性.1 控制方程与建模1.1 控制方程经过滤波后,大涡模拟的连续性和N-S方程为(1)(2)式中:t为时间;xi、xj为坐标的分量,i,j=1,2,3;ρ为流体密度;ui、uj分别为与xi、xj对应的速度分量;p为压力;μ为流体运动粘性;式中带有上划的量为滤波后的场变量,被定义为亚格子尺度应力,在Smagorinsky-lilly[13]亚格子模型中,(3)式中:τkk 为应力,k=1,2,3;δij为 Kroneker符号; Sij为变形速率张量;vSGS为亚格子模型涡粘系数,(4)其中:为格子过滤后的尺寸,和Δz为网格尺寸;(5)其中:Lkk为亚格子Leonard 应力;Lij和Mij的表达式分别为其中,Δ为格子过滤尺寸标准.1.2 建模与数值算法(1) 计算区域本文共计算3个算例,其中Case 1和Case 2的计算区域为30D×18D×4D(D为方柱断面宽度),Case 3则为30D×18D×8D.方柱模型上游来流区域长度为10D,下游尾流区长度为20D,即x方向上的长度为30D;模型离上下边界各为9D,即y方向上的宽度为18D;为了研究长径比的影响,展向长度分别取4D和8D两种长度,如图1所示.图1 计算区域示意图及边界条件Fig.1 Computational domain and boundary conditions(2) 边界条件(如图1所示)方柱表面边界条件:采用无滑移光滑固壁边界;入口边界条件:采用均匀来流速度条件,U=10 m/s;拟压(pseudo-pressure)φ采用Neumann条件;出口边界条件:速度满足格式为∂/∂t+c∂/∂x=0的对流边界条件,拟压取为0,这里c为任意常数;上下边界条件:采用自由流动边界条件;展向边界条件:采用周期边界条件.(3) 网格划分如图2所示,在方柱模型周围5D的范围内,采用贴体网格,并在方柱附近尤其是在尾流区,采用密集网格,而在远离方柱的区域使用稀疏的网格.为研究方柱表面周围网格的影响,在Case 1中,方柱表面附近的网格直接按照双曲切线间距分布由内向外伸展,如图3(a)所示;而在Case 2和Case 3中,则先在方柱表面附近采用厚度均匀的正交网格,然后再按照双曲切线间距分布由内向外伸展,如图3(b)所示.同时,为确保LES在边界层内计算的准确性,靠近方柱表面的第一层网格厚度Δd,设在展向布置上,Case 1和Case 2展向长度为4D,Case 3的展向长度为8D,三者的展向长度的解析度均为Δz/D=0.1,即在Case 1和Case 2中展向划分为40份,而在Case中展向划分为80份.因此,Case 1、Case 2和Case 3的网格总数分别为1 848 000、3 434 000和5 595 000.采用36核计算工作站,Case 1、Case 2和Case 3计算时间分别为192、240 h和320 h.(4) 数值算法在空间离散上,本文采用NVD(normalized variable diagram)格式,从而避免中心差分引起的非物理性振荡.对于时间的离散格式,采用二阶完全隐式法,以获得稳定并准确的结果.为避免计算中的数值振荡问题,本文在数值算法上采用了MIM(momentum interpolation method)方法[14];对N-S方程组采用Piso求解算法,并将连续性方程的收敛值设置为1×10-6;无量纲时间步长设置为Δt*=ΔtU/D=2.5×10-4,以确保柯朗数小于1,且保证每步计算可在20次迭代内完成.图2 整个计算区域内网格划分情况Fig.2 Grid system(a) Case 1(b) Case 2和Case 3图3 方柱表面周围网格划分情况Fig.3 Grid around square cylinder2 数值验证本节针对基于时间积分的平均积分分量,与其他实验和数值结果进行了比较,以确保本文数值计算的准确性.如表1所示,本文3个算例得到的St均与Lyn等[1]和Norberg[15]的实验值比较吻合;平均阻力系数CD与Sohankar等[5]的LES结果间误差小于1%,两者结果十分接近.需要指出的是,表1中计算气动力系数所采用的特征长度均为方柱断面宽度D.由此可见,本文进行的数值模拟是准确和可靠的.表1 平均积分分量对比表Tab.1 Comparison of integral parameters算例Re/×103StCD阻力系数脉动值C′L升力系数脉动值C′D平均背压Cpb平均回转长度LrCase 122.00.1242.340.311.521.711.32Case222.00.1212.300.321.471.621.37Case 322.00.1222.270.291.401.611.37Lyn等[1]21.40.1302.10——1.601.38Nishimura & Taniike[3]40.0—2.330.261.331.60—Norberg[15]22.00.1302.10——1.37—Sonhankar et al (LES)[5]22.00.1322.320.201.541.63—Wang & Vanke(LES)[8]21.40.1302.030.181.29—1.313 结果与讨论3.1 气动力特性表1给出并比较了其他基于时间积分的平均积分分量,包括Cpb、St和Lr (回转长度,其定义为圆柱体中心到中心线上最小流向速度Umin对应位置的距离,用于表征方柱后回流区的大小).由表1可见,本文3个数值算例得到的各项气动力统计值均与先前实验和数值结果比较吻合.对比Case 1和Case 2可以发现,前者的CD更大.原因在于Case 1的Lr更小,即方柱后方的回流区更靠近方柱后表面,这表明在方柱后方的旋涡更靠近方柱,因此导致方柱Cpb更小(如表1所示),从而使得方柱上的CD 更大.由此可见,与圆柱绕流相似,Lr是方柱绕流中的重要流场特征[16].对比Case 2和Case 3可以发现,两者的Cpb和Lr都几乎相同,但CD存在一定差别,其原因是两者的长径比不同,即Case 3中的展向长度更长,气动力的展向相关性更弱,使得其CD更小.此外,可以发现,随着展向长度的增加,均明显减小,且与实验值更为吻合.这是因为展向长度的增加导致相关性变弱,即脉动值对长径比更敏感.而先前数值计算大多考虑长径比的影响,使其计算结果在脉动项上与实验结果存在差别.同时由表1可以发现,与Case 1中各项气动力参数的结果相比,Case 2和Case 3结果更加准确,即在方柱表面附近先布置厚度均匀的正交网格,可有效改善数值结果.因此以下将主要针对后两者的结果进行说明.图4比较了方柱表面压力系数平均值Cp和脉动值的分布情况.(a) 平均值(b) 脉动值图4 方柱表面压力分布Fig.4 Surface pressure distribution on square cylinder由图4可见,本文结果与其他实验和数值结果基本吻合;且与Oka等[17]的数值结果相比,本文更加接近Nishmura等[3]的实验值.如图4(b)所示,与Case 2相比,Case 3中的脉动值与实验值更加吻合,特别是在方柱背面区域和靠近后方边角附近区域,也就是旋涡主要作用区域.这再次表明气动力系数的脉动值对展向长度更敏感. 3.2 流场特性图5给出并比较了中心平面(y=0的平面)上平均流向速度的分布情况.图5 中心线上的平均流向速度分布Fig.5 Stream-wise velocity at centreline由图5可以发现,平均流向速在方柱表面为0,在回转区域内逐渐减小直至最小值Umin;然后以近似渐近线的形式,单调地增加到外部流速U.通过对比发现本文的平均流向速度分布与Lyn等[1]的实验结果基本一致;且Lr/D=1.37与其他的数值结果(如表1所示)基本吻合.在尾流区,本文的结果较Wang等[8]的结果更加与实验值一致.但可以发现各个结果间的Lr偏差较小,而Umin/U比值偏差更大.其原因在于Umin/U较Lr对长径比更加敏感,因此,Case 2和Case 3的Umin/U比值也因为展向长度的不同而出现不同的结果.与此类似,Norberg[18]在研究圆柱绕流时,对比不同的长径比和阻塞率的试验结果后,也发现Umin对长径比更加敏感.图6给出了不同尾流位置上的分布情况(x=0对应的位置为方柱断面中心).由图可见,平均流向速度呈现为“U”形分布形态,且并随着流动向下游发展,速度分布趋于平缓,而平均竖向速度分布明显呈现出关于y=0轴的反对称性.与其他实验和数值结果对比可以发现,本文与Lyn等[1]实验值基本吻合,且明显好于Wang等[8]的数值模拟结果.同时还可以发现,在一些位置上,不同结果间存在微小差别,其主要原因在于方柱前方来流长度不同,以及早期试验设备的测量误差.对比Case 2和Case 3数据可以发现,两者结果基本一致,由此可见长径比对平均流场影响较小.流动速度的脉动值和雷诺应力分布都是流场结构中重要的湍流特征,因此可通过对速度脉动值和雷诺应力的分析,研究方柱尾流结构中的湍流特性.(a) 平均流向速度(b) 平均竖向速度图6 尾流区不同位置处平均速度分布Fig.6 Mean velocity profile in the wake由于旋涡的交替出现使得尾流区流向速度脉动值的峰值不在中心处,而是在两侧出现双峰值,如图7(a)所示.这说明在尾流区存在明显卡门涡街现象.然而随着流动向下游发展,峰值对应的尖角消失,变得更圆滑,甚至消失.这说明旋涡在下游发展中逐渐减弱.如图7(b)所示,尾流区不同位置处的脉动竖向速度均在在中心线上达到峰值,即存在单个峰值.因此表明交替出现的旋涡对流向和竖向上速度脉动值的影响是不同的.如图7(c)所示,由于交替旋涡的出现,雷诺应力(u′v′)呈现关于y=0轴的反对称性.此外以上三者均在位置x/D=1.5处达到峰值,然后随着流向而逐渐衰减,这表明旋涡强度最大值对应的位置在x/D=1.5附近,即回转长度大约为1.5D,再次证明了表1中得到的回转长度是合理准确的.同时可以发现,本文Case 2和Case 3结果均与其他结果基本吻合,但Case 3的结果更加接近Lyn等[1]实验值.由此可见,速度的脉动值对展向长度,即长径比更加敏感,因此为保证得到准确的湍流流场特性,需要选择合理长径比.为了进一步研究展向长度(长径比)对流场的影响,本文在方柱表面不同位置上,沿展向布置了一系列监测点,并通过分析得到其压力的相关性系数,如图8所示.由于展向采用周期边界条件,因此压力的相关性系数在两侧均等于1.同时可以发现随着长径比的增加,背面和侧面的压力相关性均减弱.背面压力相关性的减弱导致方柱的阻力系数的平均值和脉动值均有所减小,如表1所示;而侧面压力的减弱使得升力的脉动值明显减小,从而使得其值与实验结果更加吻合.因此,在对方柱绕流进行数值模拟时,为保证获取准确合理的平均和脉动气动力,以及准确的流场湍流特性,取8D的展向长度是十分必要的.(a) 脉动流向速度(b) 脉动竖向速度(c) 雷诺应力(u′v′)图7 尾流区不同位置处的脉动速度分布Fig.7 Fluctuating wake velocity(a) 方柱背面不同位置处(b) 方柱侧面不同位置处图8 方柱表面压力相关性系数Fig.8 Pressure correlation on cylinder surface4 结论在典型高雷诺数Re=22 000下,本文基于开源软件Openfoam平台,采用大涡模拟方法,对不同网格划分体系、不同展向长度下的方柱绕流进行三维数值模拟研究,在验证了数值结果的准确性基础上,深入分析了其尾流的流场结构和气动力特性,并给出了展向长度(长径比)对流场的影响,得到以下结论:(1) 在方柱表面附近先布置厚度均匀的正交网格,可获得更加准确的气动力结果.(2) 在方柱绕流中,回转长度是尾流结构的重要参考指标,其值大小直接影响着方柱气动力的表现,特别是断面的阻力系数.(3) 高雷诺数下,展向长度对回转长度影响较小,但对尾流流向速度分布中的Umin/U比值影响较大,即Umin/U比值对长径比更为敏感.(4) 长径比的大小直接影响方柱气动力的表现,特别是展向相关性,这使得气动力的脉动值对长径比也更加敏感,故为保证计算结果的准确性,特别为了获取更准确的湍流特性,展向长度取为8D是更为合理和有效.参考文献:【相关文献】[1] LYN D, EINAV S, RODI W, et al. A laser-Doppler velocimetry study of ensemble-averaged characteristics of the turbulent near wake of a square cylinder[J]. Journal of Fluid Mechanics, 1995, 304: 285-319.[2] SAHA A, MURALIDHAR K, BISWAS G. Experimental study of flow past a square cylinder at high Reynolds numbers[J]. Experiments in Fluids, 2000, 29(6): 553-563.[3] NISHIMURA A, TANIIKE Y. Fluctuating wind forces of a stationary two dim. square prism[C]∥Proceedings of 16th National Symposium on Wind Engineering. Tokyo: Springer, 2000: 255-260.[4] YEN S C, YANG C W. Flow patterns and vortex shedding behavior behind a square cylinder[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2011, 99(8): 868-878.[5] SOHANKAR A, NORBERG C, DAVIDSON L. Numerical simulation of flow past a square cylinder[C]∥The 3rd ASME/JSME Joint Fluids Engineering Conference. San Francisco: American Society of Mechanical Engineers, 1999: 18-23.[6] OKAJIMA A. Strouhal numbers of rectangular cylinders[J]. Journal of Fluid Mechanics, 1982, 123: 379-398.[7] SRINIVAS Y, BISWAS G, PARIHAR A, et al. Large-eddy simulation of high Reynolds number turbulent flow past a square cylinder[J]. Journal of Engineering Mechanics, 2006,132(3): 327-335.[8] WANG G, VANKA S. Large-eddy simulations of high Reynolds number turbulent flow over a square cylinder[J]. Dept. of Mechanical and Industrial Engineering Rep. No. CFD, 1996: 96-02.[9] RODI W, FERZIGER J, BREUER M, et al. Status of large eddy simulation: results of a workshop[J]. Transactions-American Society of Mechanical Engineers Journal of Fluids Engineering, 1997, 119: 248-262.[10] VOKE P R. Flow past a square cylinder: test case LES2[J]. Direct and Large Eddy Simulation Ⅱ, Spring, 1997, 5: 335-373[11] GRIGORIADIS D, BARTZIS J, GOULAS A. LES of the flow past a rectangular cylinder using the immersed boundary concept[J]. International Journal for Numerical Methods in Fluids, 2003, 41(6): 615-632.[12] 张伟,葛耀君. 方柱绕流粒子图像测速试验与数值模拟[J]. 同济大学学报:自然科学版,2009,37(7): 857-861,892.ZHANG Wei, GE Yaojun. Particle image velocimetry study and numerical simulation of turbulent near wake of square cylinder[J]. Journal of Tongji University: Natural Science, 2009, 37(7): 857-861, 892.[13] LILLY D K. A proposed modification of the Germano subgrid-scale closure method[J]. Physics of Fluids A: Fluid Dynamics, 1992, 4(3): 633-635.[14] CHOI S K. Note on the use of momentum interpolation method for unsteady flows[J]. Numerical Heat Transfer: Part A: Applications, 1999, 36(5): 545-550.[15] NORBERG C. Flow around rectangular cylinders: pressure forces and wake frequencies[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1993, 49(1/2/3): 187-196.[16] 周强,曹曙阳,周志勇. 亚临界雷诺数下圆柱体尾流结构的数值模拟[J]. 同济大学学报:自然科学版,2013,41(1): 33-38.ZHOU Qiang, CAO Shuyang, ZHOU Zhiyong. Numerical studies of wake characteristics on a circular cylinder at sub-critical reynolds number[J]. Journal of Tongji University: Natural Science, 2013, 41(1): 33-38.[17] OKA S, ISHIHARA T. Numerical study of aerodynamic characteristics of a square prism in a uniform flow[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2009,97(11): 548-559.[18] NORBERG C. An experimental investigation of the flow around a circular cylinder: influence of aspect ratio[J]. Journal of Fluid Mechanics, 1994, 258: 287-316.。

高雷诺数下双圆柱绕流的数值模拟_廖俊

高雷诺数下双圆柱绕流的数值模拟_廖俊

A辑第16卷第1期 水动力学研究与进展 Ser.A,V ol.16,N o.1 2001年3月 JOURNAL OF HYDRODYNAM ICS M ar.,2001文章编号:1000-4874(2001)01-0101-10高雷诺数下双圆柱绕流的数值模拟廖 俊1, 景思睿2(1.华中理工大学能源科学与工程学院,湖北武汉430074;2.西安交通大学能源与动力工程学院,陕西西安710049) 摘 要: 本文使用表面涡法研究高雷诺数下不同排列方式双圆柱绕流的流动状态。

计算了双圆柱在并列、串列及级列的情况下的各种流动结构,涡街的变化及作用在圆柱上的受力情况。

本文结果清楚地描述了双圆柱绕流复杂的流动状况,计算结果与实验显示的流动状况十分相似,斯特罗哈数和阻力系数与实验结果符合得很好。

关 键 词: 表面涡方法;圆柱绕流;数值模拟;涡街中图分类号: O357.1 文献标识码:A1 引言对多圆柱的绕流研究在工程实际中有很重大的意义,例如管束的热交换,反应堆,高大建筑物,海洋平台及桥梁等。

当流体流过圆柱体时,由于涡的脱落,使圆柱体上产生交变作用力。

这种作用力导致柱体的振动及材料的疲劳,而使结构损坏,产生严重的后果。

如水电站的蒸发塔,就曾经由于安装位置不正确,导致多个塔之间强烈影响、振动并使塔损坏,悬索桥也发生过类似事例,悬索共振而使桥倒塌。

由于多个柱体流动状况复杂、多变,导致对于柱体上作用力大小和方向极其复杂,实验测量非常困难,在实际工程中就需要用数值模拟的方式确定其流动状况,估计出柱体上的作用力大小、方向,以便工程参数的确定。

在多圆柱绕流研究中最多的是双圆柱绕流,双圆柱绕流按圆柱的不同排列方式可以分为三类:串列,两圆柱相对来流方向呈前后排列;并列,两圆柱相对来流方向呈并排排列;级列,两圆柱呈前后交叉排列。

对于柱体绕流的数值模拟方式可以分两大类,一类为网格法,另一类为无网格法。

网格法主要有有限差分法、有限元法。

自由流场二维圆柱阻力系数数值模拟

自由流场二维圆柱阻力系数数值模拟
基于有限体积方法fluent将计算区域分割成很多不重复的控制体积每个网格点周围有一个控制体积在控制体积上对待解的微分方程积分得出一组离散方程未知数为网格点上的因变量加上边界条件和初始条件求解方程组中的未知量针对不同的收敛准则利用fluent软件进行模拟计算的一般流程如下定义边界类型和生成网然后将gamit中的网格输出文件导入fluent求解器设置边界条件和初始条件对流动区域进行求解最后对计算结果进行后处理fluent的前处理软件主要作用是建立网格模可以建立各种几何体划分几何体为很小的网格通过设立边界条件fluen它也可以导入在cadsolidwork模软件中生成的复杂图形然后对其进行网格划分建立网格模型在二维计算中圆柱用圆来代替圆的直径为cm

黏性流体绕流圆 柱时产生两种力, 一 种是圆柱边 壁引起的 边壁剪切应力, 另一种是圆柱体前后压 差引起的 正应力。这两 种力的合力在入流方 向上的分量即常说的阻 力, 在垂 直于入流 方向上的分量即升阻 力 [ 1]。
圆柱在黏性不可 压流体流 场中的 阻力系 数随雷 诺数 (R e) 的 变化而变化 [ 2] 。当 Re < 1时, 阻力系数很大; 当 1 < Re < 3 @ 105 时, 阻力系数随 R e增大而近乎单调减小。
当 0 < R e < 4时, 流线基本对称, 特征是在任意流体单元, 因压力差异而产生的 力和摩 擦力几 乎相 等, 流态 如图 1( a) 所 示; 当 4 < R e < 40时, 在圆柱后出现层流 分流, 即 形成两个稳 定 的方向相反的驻涡, 如图 1( b) 所示; 当 Re > 40时, 圆柱后部 的驻涡开始 不稳定地 摆动、脱落, R e达到 100 左右时这 对不稳 定的对称旋涡分裂, 最后 形成 几乎稳 定的、非对称 的、规则 的、 旋转方向相反的 交替 旋涡, 如 图 1( c) 所 示; 当雷 诺数 再 增大 时, 流动分离严重, 约从 R e = 10 000起, 层流边界层从圆柱的前 部距离停滞 点约 80b位置处开始分离, 流动变为湍流, 形成很宽 的 分离区, 如图 1( d) 所示; 当雷诺数超过 100 000后, 层流边界 层变为湍流边界层, 分离点向后推移 , 距离 停滞点 约 120b 位置 处, 分离区变小, 如图 1( e) 所示 [3] 。

不同雷诺数下方柱绕流的数值模拟

不同雷诺数下方柱绕流的数值模拟
图2 不同雷诺数下的涡线图
2
计算结果与讨论
分别对雷诺数为 100, 1 103 , 1 104 和 2 2
3 4
计算得到了方柱上的阻力系数和升力系数. 阻力系数 C D 与升力系数 C L 的定义分别为: CD = FD , 1 L 2 U2 2 CL = FL 1 L 2 U2 2 ( 12)
104 时的情况进行了计算 . 当 Re = 100 时, 直接采 用 N S 方程进行计算 ; 当 Re = 1 10 , 1 10 和 2 2 10 4 时 , 则引入 k 湍流模型进行计算. 下面 给出不同雷诺数下的计算结果 . 图 2 给出了计算得到的涡线图 . 在 4 个不同 的雷诺数下, 都会在柱体的尾部出现规则的旋涡 脱落 , 尾涡交替的甩在方柱上下两侧的壁面上 , 这 就是著名的卡门涡街 . 同时我们可以发现, 方柱后 尾涡的形态会随雷诺 数的变化而产 生一定的变 化. 当雷诺数较低时 ( Re = 100) , 尾涡会拖得比较
[ 6]
图1
计算模型
各边界条件分别为: 入口 : 给定无量纲速度, u = 1, v = 0 . 出口 : 给定无量纲压力 p = 0 , 速度采用 u / x = 0, v / x = 0 . 固壁 : 采用无滑移边界条件, 即 u = 0, v = 0 . 上下边界 : u = 1, v = 0. 不可压缩牛顿流体运动的控制方程 N S 方 程可表示为: ui = 0 xi u i + uj t 式( 2 ) 中 ui = - 1 xj p + xi ui xj xj
[ 4]
得到了绕流的速度场以及柱体上的受力参数. 总 结了方柱后的涡脱落形式和流场的动力学参数随 雷诺数的变化规律 , 并将计算结果与前人的实验 和计算结果进行了对比.

基于ADINA的二维双圆柱绕流的数值模拟研究

基于ADINA的二维双圆柱绕流的数值模拟研究

基于ADINA的二维双圆柱绕流的数值模拟研究作者:党鹏飞董事尔王华洪陈果来源:《中国新技术新产品》2013年第02期摘要:均匀来流流过二维圆柱是一个经典的流体力学问题,尤其是对于粘性流体,由于雷诺数的大小不同,在面对层流和紊流两种与众不同的流场时,流场流线运动的规律较为复杂,本文正是借助于ADINA软件中出色的流固耦合的仿真计算技术,对于流场中二维双圆柱绕流流场的变化进行了科学的数值模拟,并给出了不同环境条件下流场的变化情况,计算出了圆柱表面的一系列动力学参数。

结果表明:尾流及圆柱表面的压力分布,其实验结果与现有结果较为吻合关键词:雷诺数;圆柱绕流;网格密度;数值模拟;扰动力中图分类号:TP39 文献标识码:A1 概述近些年来一些学者运用实验和理论方法对横向流作用下管阵流体诱发振动问题进行了分析和研究,得到一些经验公式来初步估计产生流体诱发振动的临界流速。

并对两圆柱串列和交错放置的绕流问题进行过实验研究。

针对两圆柱中心间距小于5.0倍圆柱直径的一系列情况,他们研究了两圆柱间的流动相互作用,发现中心间距存在有一临界值,当小于该临界值时,没有明显的涡自上游圆柱脱落。

这一临界值约为3 倍圆柱直径。

standsby在1981 和1987 年分别用离散涡方法和随机涡方法研究了并排、串列和交错放置的双圆柱绕流问题得到了与实验相符的结果。

但是经验公式中的一些参数是在一些特定条件下得到的,具有很大得保守性和不确定性。

双圆柱绕流模拟由于在一定范围内能够反映多个圆柱在一条直线上的绕流特征,圆柱附近流态的瞬时变化形式,并且模型简明,已经用ADINA软件能够计算比较精确的扰动力数值。

因此基于现有的研究成果,本文旨在归纳总结双圆柱对绕流流场的影响。

2 双圆柱体绕流场基本理论根据prandtl的边界层理论,圆柱的绕流流动可以分为两个区,圆柱表面很薄的边界层区和其上的主流区,在边界层中流体粘性产生的摩擦力起主导作用,而在主流区粘性摩擦力可以忽略不计。

最新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三分立系数的求法人民法院--庭审流程图书记员首先入庭【站立】:现宣布法庭纪律书记员【宣告完毕】:全体起立,请审判长审判员入庭【审判长一行依次入庭就坐】书记员【清点当事人及其代理人】:报告审判长,本案当事人及诉讼代理人已全部到庭,请开庭。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

(在数学上即自由滑移边界条件) ;当流场发展到出口位置时,认为流动已经充 分发展,因此设置出口为自由流 outflow;入口在用速度入口边界条件,入口速 度、湍流强度、湍流长度尺寸均采用算例中所给设定。 在设置参考值时首先选定 compute from->流场名,这样就正确的给出了参 考值中的物理参数(密度、粘度等) ,然后指定 Area 的值为 0.1。 关于参考值的设定是非常有讲究的,根据文献十八 P173 的升阻力系数定义 式: Cd = Fd 1 2 2 ρv S Fl 1 2 2 ρv S
图 15-1-13 方柱绕流流场的总压云图与速度云图
图 15-1-14 给出了方柱升阻力系数的历史曲线。
图 15-1-14 方柱绕流的升力与阻力系数的历史曲线
从图中可以明显看出, 方柱绕流的升阻力系数随着时间的推移呈周期性的波 动。经过对数据的统计,可以知道,方柱阻力系数的平均值为 2.3(取 80-90s 之间的每个时间步的阻力系数做平均) ,实验值为 2.1,误差为 9.5%,计算出的 漩涡脱落周期为 3s,由此所得的斯特劳哈尔数为 0.15,实验值为 0.13,误差为 15.4%。
图 15-1-9 稳态计算的残差曲线
图 15-1-9 给出了稳态计算残差曲线,从曲线中可以看出,残差曲线基本上 随着计算会继续趋于收敛, 其中连续性方程的残差最大。在计算此类绕流问题的 实际应用过程中, 还会出现当稳态计算的残差位于高位稳定的情况,出现这种情 况时, 应当考虑的一种可能存在的情况就是物理模型中存在瞬态变化的物理效应, 导致稳态迭代时前后迭代的范数残差不能够收敛而位于高位稳定, 这种现象在用 稳态雷诺平均方法求解瞬态问题时经常出现, 此时可以将稳态计算结果作为瞬态 计算的初场,继续计算。所以单看残差收敛情况不一定有太大意义,要结合实际 情况来看。另外发现当启用 Differential Viscosity Model 选项后,本来收敛 至非常小的残差容易转变为高位稳定甚至发散。 本例中升力系数与阻力系数的迭代历史曲线如图 15-1-10 所示:
图 15-1-2 流场模型
按照算例内容建立流场模型如图 15-1-2 所示,然后进行分网,由于采用的 是尺度化壁面函数法,所以近壁面网格没必要划分地太细。 进入到 Mesh 模块后首先到 mesh 模型树下将 Physics Preference 改为 CFD 与 Fluent,然后设置 Sizing 下的网格相关度、网格光顺等参数。
需要注意的是,在计算升力与阻力时,升力方向垂直于来流方向(不是垂 直于机翼! ) ,阻力方向平行于来流方向(不是平行于机翼! ) 。仅当攻角为 0°或 180°时,升力/阻力才垂直/平行于来流方向。 在 Fluent 二维计算中,机翼长度方向(展向)垂直于几何平面,软件在处 理时机翼长度置为 1 (在流体力学中处理二维问题, 展向长度全部取单位长度) 。 其余参考值需要参考升力系数公式进行计算,并且一定要注意升力/阻力的方向。 对于本二维算例,S 的取值与弦长数值相等,因此本算例 S 取值为 0.1,即 设定参考值时,与尺寸相关的进设置 Area 为 0.1 平方米即可,其余的保持默认 不变;密度等材料的参考值已经通过 Compute from 功能实现了。 之后设置求解器选项, 选择耦合求解器, 对流项的处理采用二阶迎风格式 (对 于之前划分的网格,也可采用更高阶的 QUICK 格式) 。保持默认的收敛残差判定 标准 0.001 不变, 监测了升阻力系数随迭代步数的变化曲线。 以入口条件初始化 整个流场区域,设定最大迭代步数为 1000 步。
其中,St 为斯特劳哈尔数,f 为漩涡脱落频率,v 为来流速度。根据文献十 七中的实验结果(St=0.13) ,可以反推出漩涡脱落频率应为 3.5s,所以在选择 瞬态步长时就有了选择的依据, 瞬态求解采用二阶隐式求解方式,针对方柱绕流 问题,隐式方法对时间步长的要求不高,很多人取漩涡脱落周期的 1/50 作为瞬 态计算的时间步长,本文取为 0.02s,可以充分地描述出漩涡脱落的瞬态细节。 对于收敛残差的取值, 很多人是取稳态残差稳定值的 1/50, 本文取为 1/100, 即 0.00001。指定每一时间步最0s。
Cl =
以上两式中, Cd 代表阻力系数, Cl 代表升力系数, Fd 为钝体所受到的阻力, Fl 为钝体所受到的升力,ρ 为流体的密度,v 为前方来流的速度,S 代表钝体的 投影面积,在二维问题中 S 的数值与弦长相等(弦长乘以一,二维问题中展向长 度为单位长度一) 。 有关翼型的概念: 翼型:飞机机翼具有独特的剖面,其横断面(横向剖面)的形状称为翼型, 称为翼型 前缘:翼型最前面的一点。 后缘:翼型最后面的一点。 翼弦:前缘与后缘的连线。 弦长:前后缘的距离称为弦长。如果机翼平面形状不是长方形,一般在参 数计算时采用制造商指定位置的弦长或平均弦长 攻角:也称迎角。对于翼形来说,攻角定义为翼弦与来流速度之间的夹角, 抬头为正,低头为负,常用符号 α 表示。其示意图见图 15-1-7:
图 15-1-10 方柱绕流的稳态升阻力迭代历史
从图 15-1-10 中可以看出, 升阻力系数在稳态计算中最终归于稳定,这与残 差收敛至极小值是相互呼应的。 而在有些绕流问题的稳态计算中,也会出现瞬态 特征。如图 15-1-11 所示:
图 15-1-11 某圆柱绕流的稳态计算
图 15-1-11 中给出了某圆柱绕流的稳态计算结果, 在稳态计算结果中出现了 瞬态特征。 在稳态计算中出现瞬态特征的现象,在方柱绕流与圆柱绕流的现象中
图 15-1-15 方柱表面的 Y+值分布范围
仿真后可以知道,近壁面 Y+值的分布范围为 7.07-68.6,在 Fluent 的官方 文献中推荐壁面函数法的 Y+值适用范围为 30-500,对于尺度化壁面函数,由于 在近壁面使用了预防 Y+值过小导致结果恶化的技术,其 Y+值范围下限可以适当 放宽, 但仍应尽量取在 11.25 以上。 过大或过小的 Y+值都可能引起结果的恶化。
(3)对于存在较高导热率的共轭传热问题; (4)对于流场网格长宽比较大的问题。 在本算例中,虽然不存在以上问题,但仍采用双精度 4 核并行计算求解。 在进行计算时, 先进行稳态计算,然后将稳态计算的结果作为瞬态计算的初 场,再进行瞬态计算。 求解器采用压力基求解器,并进行湍流模型的设置。
图 15-1-6 湍流模型的选择
都可能出现,当出现这一现象时,残差曲线将会居于高位,并且升力与阻力系数 的值不一定准确,比如上图中升力系数平均值部位 0。 还有一种现象是升阻力系数趋于平稳,但残差曲线局于高位,造成这种现象 的一种原因是启用了 Differential Viscosity Model 选项,这也说明了 Differential Viscosity Model 选项对计算所造成的不稳定性。 根据文献十八 P119 的关于斯特劳哈尔数的公式: St = fD v
关于雷诺数 22000 方柱绕流 (关于壁面函数法、湍流粘度比的计算等概念可以参考 文献十六)
15.1.1 RNG 模型与尺度化壁面函数在方柱绕流中的应用
图 15-1-1 分析流程
图 15-1-1 给出了二维方柱绕流问题的分析流程, 在 Workbench 中建立模型, 然后使用 Mesh 模块进行分网,最后提交到 Fluent 中进行计算。
图 15-1-12 方柱绕流流场的静压云图与动压运图
根据图 15-1-12 中可以看出, 在方柱后方出现了明显的周期性漩涡脱落的现 象。漩涡脱落的现象在总压云图与速度云图中也可以很明显地看出来。 文献十七中认为 RNG 方法在预测平均阻力系数和回流区再附着距离时与实 验及大涡结果比较接近,但是在预测升阻力系数的均方根值时仍具有一定的差 距。采用雷诺平均的 RNG 模型,从模型精度上来说要比大涡模拟至少低一个数 量级。
图 15-1-3 网格划分结果
图 15-1-3 给出了网格划分结果,在划分网格时,约束整个面的网格尺寸为 10mm, 并且定义入口、 出口、 上下边界、 方柱表面的名称集 (Named Selections) , 以方便在 Fluent 中定义边界条件。
图 15-1-4 局部网格示意
图 15-1-4 给出了局部网格示意图,从图中可以看出,在方柱表面的网格高 度是较高的, 使用尺度化近壁面函数法时并没有划分很细的边界层网格。最后共 得到 48900 个网格,网格最大长宽比为 1.0069,满足仿真计算的网格要求。
如图 15-1-6 所示,选择 RNG 湍流模型与尺度化壁面函数处理方法,并且不 启用 Differential Viscosity Model 选项,在 ANSYS Fluent 14.5 的帮助文件 中对该模式有以下论述: The RNG turbulence model in ANSYS FLUENT low-Reynolds-number effects. 可见在考虑到边界层的低雷诺数影响时应该打开该选项,另外,有人在做空 气的涡街仿真时遇到过不勾选该选项就得不到涡街的情况。 (/bbs/viewthread.php?tid=76599) 但是在实际计算中发现,启用该选项后,计算容易变得不稳定,出现了数值 波动甚至是不收敛。因此,在本算例中,不启用该选项。 流线曲率修正是因为使用各向同性涡粘假设的湍流模型对流线曲率不敏感, 因此在流线弯曲时需要进行流线曲率修正,启用该选项。 选择 Fluent 中的 water-liquid 作为流体材料,然后进行边界条件的设置: 设置方柱表面为无滑移壁面条件 wall;流场上下边界为对称边界条件 symmetry has an option of using a differential formula for the effective viscosity to account for the
图 15-1-5 采用 4 核并行计算
在打开 Fluent 时,按照图 15-1-5 所示,采用单精度、4 核并行计算。双精 度一般用在以下几个方面: (1)存在明显的细长比很大的几何特征,如细长的管道; (2)流体是几个区域,区域之间用小尺寸的管道连接,其中一个区域的流 体压力大大高于整个流域的平均水平;
相关文档
最新文档