不可压缩圆柱绕流

不可压缩圆柱绕流摘要本文对无粘不可压圆柱绕流问题进行数值模拟，并将数值结果和理论解进行对比。 在给定圆柱固壁和远场速度的边界条件下，通过求解流函数的Laplace方程，得到域内 的流函数，进而计算出速度场。数值模拟的主要步骤有网格划分，物理域和计算域坐标 变换，对流函数方程进行差分离散，SOR迭代方法求解和后处理绘出速度场。本文将 先建立本问题的数学提法和理论解，随后基于数学提法进行数值求解，最后将理论和数 值进行对比。理论和数值得到的速度场表现出了一致性，并且均反映出无粘假设下的流 动和实际流动间的巨大差异。关键词不可压缩；圆柱绕流；坐标变换；SOR方法1问题提出考虑图1所示无粘不可压圆柱绕流，來流速度为Uoo = 20m/s,Vqo = Oo选取合适的 计算域范围和合理的边界条件，通过求解流函数方程或者速度势方程求圆柱周围的流场 并与理论解相比较。图1无粘不町压圆柱绕流示意图2理论求解在圆柱圆心处取笛卡尔坐标，对于不可压缩二维流动，存在流函数W有翌v坐dy dx.中，流场始终无旋。由VxV = 0,代入(1)可得主控方程Ay/= 0(2)定解的边界条件有两个，一个是圆柱的固壁条件，考虑到固壁的不可穿透性有几+心=const上式中1;为圆柱半径，第二个边界条件是无穷远处速度分布(5)在圆柱坐标系中求解上述问题，坐标定义如图2,可得到坐标变换关系和微分关系x = r cos 0. y = r siii 0化叮0 0=(COS0sin6sinC cos。vx y丿k rr丿图2坐标关系说明根据(5)(6)进行坐标转换，可得到主控方程和边界条件如下咯粋寻八r di 厂 &(8)=20 sin 0.di/der-co=20r cos 0(9)通过分离变量法可得流函数的解为(10).-201；2 sin 0屮=20r sm 0E+ c0根据(5)(10)，算得理论速度场为U = 0Iy + %&y = 2O 20r；(x2-y2)(x2 + y2)2v=_(. + %q)= _401：号(x2 + y2)2(ii)3数值求解3.1网格划分与坐标变换设计求解物理域为与圆柱同轴的圆环域，范围为7 = 20m,= 100mo为保证网格正交性，采用O型网格对物理域进行划分，如图3所示。网格沿径向有161个节点， 间距为dr = 0.5m：沿周向有361个节点，间距为dd = 1,第一个节点和第361节点在 物理域上是重合的。100806040200201050050100图3物理域网格取计算域为Cartesian网格，坐标系记为“坐标，每个网格的大小为1X1,即 A? = Aq = lo计算域中点(i,j)对应物理域中点(rcos0,rsin0),其中r = rc + dr(i- 1), 0=O-l)deo根据此映射关系，可以数值求得从%-y到f 耳的变换微分关系、J、011、012 和022(12)对于二维平面任意曲线坐标系中的测度以及变化Jacobian行列式J的表达式为(13)度量张量为S11 1221 22 11 =+ 2, S12 = 21 =負Z + ?yy，22 =+3.2主控方程离散与迭代求解由于坐标变换从兀-y到-小 因此流函数方程厂卩=0变换为(14)引入中间变量d z du/ 屮、 d z du/鬲(珈鬲+Jgl2药)+乔阳冠+尬=0(15)Og dr)dr)(16)则有(17)按图4所示节点排列进行差分:心+i）(2-1打)(2J)(2 十 1J)-t图4节点离散示意图在（i,j）点进行差分离散:(18)右(4/2,j - -1/2, j) + 占-= 其中，中间变量差分为:%厂（曲心笔尹）+（% G （经4）(19)%二（曲皿（笔严）+（捡站（仝亠亠）4 A 77(20)34（细叫二严如）+（%（专）(21)务（纽牛）+（% ）卄（纽严）4歹(22)整理以上5个式子，即可算出卑:厂 利用SOR迭代方法得到新的卩善+1(23)嬉=（1-对必j+咙j3.3边界条件根据边界条件和边界点的坐标变换，代入u =20, v=0,可推出边界处的流函数梯度7冬和屮H；在圆柱表面，取Tr=rc = 0o据此即给定了边界条件。计算域中所有点的流函数初值均为0,圆柱表面流函数值始终为0。除了内外边界， 其他点均可由主控方程算出数值。计算时，先算f = 2时各点的值，依次类推，直到算 完 = 160的各点的值。计算外边界，即 = 161各点时，釆用如下办法这样使得计算外边界流函数时，能同时利用到两个方向的梯度壮和。3.4收敛条件计算两轮迭代的流函数点阵差值的二范数，当该值小于UP时，即判定为收敛。4计算结果及对比4.1数值解图5为数值解的流函数分布，流函数等值线即流线，因此，根据该图可以很明显看 到不可压无粘绕流特征：即流体在圆柱前一分为二，在圆柱后重新合二为一，流函数关 于兀=0对称分布，关于y = 0反对称分布。这个分布很好的反映了该流动无损失，且圆 柱受到合力为0的理论估计。100濫函数分布1500806010004050020-20-500-40-1000-60-80-150050100一1%A 0x图5数值仿真所得流函数分布图6、图7和图8分别为u分量、v分量和速率分布，再次验证了分流特征和关于 流场对称性的分析。值得注意的是，圆柱前后均有速度为0的滞止区域，近圆柱面的绕 流先加速，在圆柱上下顶点达到最大值，之后减速离开圆柱。速度U分量分布速度分量分布100806040200-20-40-60-8050100一1%XI巾-50图7数值解速度v分量分布图8数值解速率分布根据 Cauchy-Lagiange 积分 + l|AP + P + n = c(t)(25)况 2p这里e是势函数，对二维定常问题，且不计体积力，有(u2 + v2)- = c(26)2P假定Poo=100kPa, p=1.2kg/m3可以求得全场的压力分彳j,如图9所示，首先验证了圆柱所受合力为零，另外，在前后滞止区内有压力最大值，在上下顶点有压力最小值。%10.028010.01601040209.999.98-209.97-409.96-60-809.95-50601000x沃0图9数值解圧力分布4.2数值与理论对比记理论所得流函数为旳，则流函数偏差分布咒=屮一卩T，如图所示。流函数偏差分布1001008060504020-20-40-50-60-80-5050100 0图10流函数偏差分布基于偏差，该分布中偏差较小，偏差最大值在100量级，而数值解最大值在1500 量级；第二，该分布1朱之墀，王希麟.流体力学理论例题与习题M.北京：清华大学出版社，1986： 66