不可压缩圆柱绕流_坐标变换_SOR方法分解.docx

不可压缩圆柱绕流 摘 要 本文对无粘不可压圆柱绕流问题进行数值模拟，并将数值结果和理论解进行对比。在给定圆柱固壁和远场速度的边界条件下，通过求解流函数的Laplace方程，得到域内的流函数，进而计算出速度场。数值模拟的主要步骤有网格划分，物理域和计算域坐标变换，对流函数方程进行差分离散， SOR迭代方法求解和后处理绘出速度场。本文将先建立本问题的数学提法和理论解，随后基于数学提法进行数值求解，最后将理论和数值进行对比。理论和数值得到的速度场表现出了一致性，并且均反映出无粘假设下的流动和实际流动间的巨大差异。 关键词 不可压缩；圆柱绕流；坐标变换；SOR方法 1 问题提出 考虑图 1所示无粘不可压圆柱绕流，来流速度为u∞=20m/s,v∞=0。选取合适的计算域范围和合理的边界条件，通过求解流函数方程或者速度势方程求圆柱周围的流场并与理论解相比较。 图 1 无粘不可压圆柱绕流示意图 2 理论求解 在圆柱圆心处取笛卡尔坐标，对于不可压缩二维流动，存在流函数有 1() 无穷远处来流u∞=20m/s，为无旋流动，根据Kelvin-Lagrange定理可知，本问题中，流场始终无旋。由，代入(1)可得主控方程 2() 定解的边界条件有两个，一个是圆柱的固壁条件，考虑到固壁的不可穿透性有 3() 上式中为圆柱半径，第二个边界条件是无穷远处速度分布 4() 在圆柱坐标系中求解上述问题，坐标定义如图 2，可得到坐标变换关系和微分关系 5() 6() 图 2 坐标关系说明 根据(5)(6)进行坐标转换，可得到主控方程和边界条件如下 7() 8() 9() 通过分离变量法[1]，可得流函数的解为 10() 根据(1)(5)(6)(10)，算得理论速度场为 11() 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-1), θ=(j-1)dθ。根据此映射关系，可以数值求得从x-y到ξ-η的变换微分关系、J、g11、g12和g22 12() 对于二维平面任意曲线坐标系中的测度以及变化Jacobian行列式J的表达式为 13() 度量张量为 14() 3.2 主控方程离散与迭代求解 由于坐标变换从x-y到ξ-η，因此流函数方程?2Ψ=0变换为 15() 引入中间变量 16() 则有 17() 按图 4所示节点排列进行差分： 图 4 节点离散示意图 在点进行差分离散： 18() 其中，中间变量差分为： 19() 20() 21() 22() 整理以上5个式子，即可算出Ψi,j*，利用SOR迭代方法得到新的Ψi,jn+1 23() 3.3 边界条件 根据边界条件和边界点的坐标变换，代入u=20，v=0，可推出边界处的流函数梯度Ψξ和Ψη；在圆柱表面，取Ψr=rc=0。据此即给定了边界条件。 计算域中所有点的流函数初值均为0，圆柱表面流函数值始终为0。除了内外边界，其他点均可由主控方程算出数值。计算时，先算ξ=2时各点的值，依次类推，直到算完ξ=160的各点的值。计算外边界，即ξ=161各点时，采用如下办法 24() 这样使得计算外边界流函数时，能同时利用到两个方向的梯度Ψξ和Ψη。 3.4 收敛条件 计算两轮迭代的流函数点阵差值的二范数，当该值小于10-5时，即判定为收敛。 4 计算结果及对比 4.1 数值解 图 5为数值解的流函数分布，流函数等值线即流线，因此，根据该图可以很明显看到不可压无粘绕流特征：即流体在圆柱前一分为二，在圆柱后重新合二为一，流函数关于x=0对称分布，关于y=0反对称分布。这个分布很好的反映了该流动无损失，且圆柱受到合力为0的理论估计。 图 5 数值仿真所得流函数分布 图 6、图 7和图 8分别为u分量、v分量和速率分布，再次验证了分流特征和关于流场对称性的分析。值得注意的是，圆柱前后均有速度为0的滞止区域，近圆柱面的绕流先加速，在圆柱上下顶点达到最大值，之后减速离开圆柱。 图 6 数值解速度u分量分布 图 7 数值解速度v分量分布 图 8 数值解速率分布 根据Cauchy-Lagrange积分 25() 这里Φ是势