不可压缩圆柱绕流坐标变换SOR方法

1、不可压缩圆柱绕流摘 要 本文对无粘不可压圆柱绕流问题进行数值模拟，并将数值结果和理论解进行对比。在给定圆柱固壁和远场速度的边界条件下，通过求解流函数的Laplace方程，得到域内的流函数，进而计算出速度场。数值模拟的主要步骤有网格划分，物理域和计算域坐标变换，对流函数方程进行差分离散， SOR迭代方法求解和后处理绘出速度场。本文将先建立本问题的数学提法和理论解，随后基于数学提法进行数值求解，最后将理论和数值进行对比。理论和数值得到的速度场表现出了一致性，并且均反映出无粘假设下的流动和实际流动间的巨大差异。关键词 不可压缩；圆柱绕流；坐标变换；SOR方法1 问题提出考虑图 1所示无粘不可压圆柱绕

2、流，来流速度为u=20m/s,v=0。选取合适的计算域范围和合理的边界条件，通过求解流函数方程或者速度势方程求圆柱周围的流场并与理论解相比较。图 1 无粘不可压圆柱绕流示意图2 理论求解在圆柱圆心处取笛卡尔坐标，对于不可压缩二维流动，存在流函数有 无穷远处来流u=20m/s，为无旋流动，根据Kelvin-Lagrange定理可知，本问题中，流场始终无旋。由，代入可得主控方程定解的边界条件有两个，一个是圆柱的固壁条件，考虑到固壁的不可穿透性有 上式中为圆柱半径，第二个边界条件是无穷远处速度分布在圆柱坐标系中求解上述问题，坐标定义如图 2，可得到坐标变换关系和微分关系图 2 坐标关系说明根据进行坐

3、标转换，可得到主控方程和边界条件如下 通过分离变量法1，可得流函数的解为根据，算得理论速度场为 3 数值求解3.1 网格划分与坐标变换设计求解物理域为与圆柱同轴的圆环域，范围为 rc=20m，r=100m。为保证网格正交性，采用O型网格对物理域进行划分，如图 3所示。网格沿径向有161个节点，间距为dr=0.5m；沿周向有361个节点，间距为d=1°，第一个节点和第361节点在物理域上是重合的。图 3 物理域网格取计算域为Cartesian网格，坐标系记为-坐标，每个网格的大小为1×1，即=1。计算域中点(i,j)对应物理域中点(rcos,rsin)，其中r=rc+dr(i

4、-1), =(j-1)d。根据此映射关系，可以数值求得从x-y到-的变换微分关系、J、g11、g12和g22对于二维平面任意曲线坐标系中的测度以及变化Jacobian行列式J的表达式为度量张量为3.2 主控方程离散与迭代求解由于坐标变换从x-y到-，因此流函数方程2=0变换为 引入中间变量 则有 按图 4所示节点排列进行差分：图 4 节点离散示意图在点进行差分离散： 其中，中间变量差分为： 整理以上5个式子，即可算出i,j*，利用SOR迭代方法得到新的i,jn+13.3 边界条件根据边界条件和边界点的坐标变换，代入u=20，v=0，可推出边界处的流函数梯度和；在圆柱表面，取r=rc=0。据此即

5、给定了边界条件。计算域中所有点的流函数初值均为0，圆柱表面流函数值始终为0。除了内外边界，其他点均可由主控方程算出数值。计算时，先算=2时各点的值，依次类推，直到算完=160的各点的值。计算外边界，即=161各点时，采用如下办法这样使得计算外边界流函数时，能同时利用到两个方向的梯度和。3.4 收敛条件计算两轮迭代的流函数点阵差值的二范数，当该值小于10-5时，即判定为收敛。4 计算结果及对比4.1 数值解图 5为数值解的流函数分布，流函数等值线即流线，因此，根据该图可以很明显看到不可压无粘绕流特征：即流体在圆柱前一分为二，在圆柱后重新合二为一，流函数关于x=0对称分布，关于y=0反对称分布。这

6、个分布很好的反映了该流动无损失，且圆柱受到合力为0的理论估计。图 5 数值仿真所得流函数分布图 6、图 7和图 8分别为u分量、v分量和速率分布，再次验证了分流特征和关于流场对称性的分析。值得注意的是，圆柱前后均有速度为0的滞止区域，近圆柱面的绕流先加速，在圆柱上下顶点达到最大值，之后减速离开圆柱。图 6 数值解速度u分量分布图 7 数值解速度v分量分布图 8 数值解速率分布根据Cauchy-Lagrange积分这里是势函数，对二维定常问题，且不计体积力，有假定p=100kPa，=1.2kg/m3可以求得全场的压力分布，如图 9所示，首先验证了圆柱所受合力为零，另外，在前后滞止区内有压力最大值，在上下顶点有压力最小值。图 9 数值解压力分布4.2 数值与理论对比记理论所得流函数为T,则流函数偏差分布e=-T,如图所示。图 10流函