写在前面
该系列文章是笔者的学习记录,旨在通过分享促进交流加深理解,文中难免有许多纰漏与理解不当之处,若能批评指正必定感激不尽。
侧重分享经典CFD方法的具体实现结果,详细的理论推导请参考教材。
一、涡量流函数方程以及求解流程
涡量流函数方程:
边界条件:
在给定 n 时刻涡量分布 的情况下,采取如下的步骤计算出下一时刻的涡量分布
:
- 计算 n+1 时刻内点的涡量(FTCS:空间导数项采用中心差分格式离散,时间导数项采用向前Euler法离散;时间一阶精度,空间二阶精度)
其中是标准的中心差分格式,
是离散的五点Laplace算子。
- 计算 n+1 时刻的流函数,差分方程和边界条件如下(中心差分,二阶精度)
为加快解方程组速度,采用泊松方程求解常用的SOR迭代公式,0<w<2 为松弛因子:
- 计算 n+1 时刻边界上的涡量:
一般形式下的Thom壁涡公式
假设某壁面切向速度,沿其内法向
有一节点,距离壁面距离为
,此点上的流函数值为
,如果壁面上的流函数值为
,那么在边界是流线的情况下有:
(注意,这里的切方向
规定为内法向量
逆时针旋转90°后的方向)
Woods壁涡公式(更高精度),含义同上
二、二维方腔流动问题
- 边界条件:
为了避免速度场在角点附近的奇性,可以对上边界给出一个速度分布其余边界处为无滑移边界条件,即
;
因为所有的边界都是固壁边界,不妨设边界上流函数。
- 初始条件:
初始时刻流体静止,内点速度均为0;根据初始速度分布计算初始的内点涡量(中心差分);流函数求解泊松方程得到。
三、计算结果
1. 中间线上的速度分布比较(上边界速度设为1,便于比较)
注:
基准数据来源于文献"High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method"
文献中采用了加密且优化的网格,我们的计算采用的是均分网格。
文献中设置驱动边界上速度等于1,为便于比较,在这一部分将上边界速度设置为1,而不采用
的光滑的速度分布。
Re=100
![](https://i-blog.csdnimg.cn/blog_migrate/e53873a6eb8c61bfe890463d4b1bb18e.png)
![](https://i-blog.csdnimg.cn/blog_migrate/e7e0148dc6f9dde270b1ad5f74ad5cb8.png)
Re=400
![](https://i-blog.csdnimg.cn/blog_migrate/517decf83db81ae91ae4c7313b0e618f.png)
![](https://i-blog.csdnimg.cn/blog_migrate/e86d4ae938e1948ef7f2186603f9311b.png)
Re= 1000
![](https://i-blog.csdnimg.cn/blog_migrate/3db9ba2a6ac49812446e5066ef963771.png)
![](https://i-blog.csdnimg.cn/blog_migrate/579fe97bf2ff5a0b7160161dc1c100fd.png)
- 可以看到,在雷诺数不超过1000时,我们的计算结果是很准确的。
Re=3200
![](https://i-blog.csdnimg.cn/blog_migrate/dbf8566492af4d41ccece7204bb47b5f.png)
![](https://i-blog.csdnimg.cn/blog_migrate/62ee95aeb3ef7fbe01d5289e37d2c928.png)
Re = 5000
![](https://i-blog.csdnimg.cn/blog_migrate/5144fd983b24255284d76979fd25b543.png)
![](https://i-blog.csdnimg.cn/blog_migrate/69be3c0762d6a7a11de9ea3fbe50bdf0.png)
- 雷诺数大于3200时,由于文献中采用了加密且优化的网格,我们的计算结果与文献数据有微小偏差,特别是靠近边界处误差较大;但大体趋势依然是可信的。
- 雷诺数大于7500时,实验发现,采用上边界速度为1的设置会导致无法收敛到正确的流函数分布,因此改用了其他边界速度和边界涡量的设置,不再与文献结果进行比较。
不同雷诺数汇总比较
观察第一幅图竖直中线上的水平速度的分布,可见随着Re增大,边界层越来越薄,Re超过3200后边界层变薄的速率就很小了;高Re时,方腔中部的速度分布几乎是线性的;Re>=3200时,在靠近y=1处u的分布会有一个小的凸起,这也与文献中提到的现象一致。水平中线上的竖直速度的分布也是类似的规律。
2. 流线图和涡旋情况(上边界速度采用多项式形式)
参照:文献中的流线图
Re=100
![](https://i-blog.csdnimg.cn/blog_migrate/442d9c1fcd4561442185e425790eec0d.png)
![](https://i-blog.csdnimg.cn/blog_migrate/c8f9c835585fc62fd372c293359ca5c9.png)
![](https://i-blog.csdnimg.cn/blog_migrate/2e732d60000550116e697d573b00e251.png)
Re=400
![](https://i-blog.csdnimg.cn/blog_migrate/fbfa31b79f0669fe93cc7ea8bc6a8f20.png)
![](https://i-blog.csdnimg.cn/blog_migrate/994978f0b2863f0810031d28c3f19d9b.png)
![](https://i-blog.csdnimg.cn/blog_migrate/d52dca84bf19bf5b482ecda3d6d28cba.png)
Re= 1000
![](https://i-blog.csdnimg.cn/blog_migrate/6dedede81efa536ace053581a32fa90f.png)
![](https://i-blog.csdnimg.cn/blog_migrate/593b8909ad0f148e4ef8aef5d7a85909.png)
![](https://i-blog.csdnimg.cn/blog_migrate/728c98ff624c5808a772655217e7d238.png)
Re=3200
![](https://i-blog.csdnimg.cn/blog_migrate/7e55f0b518faa9d99bc94ff2b742ea1a.png)
![](https://i-blog.csdnimg.cn/blog_migrate/de7197bfb058fd1f275c12690ec34dbd.png)
![](https://i-blog.csdnimg.cn/blog_migrate/3a5d6aefa104a89808bc0f94469024be.png)
![](https://i-blog.csdnimg.cn/blog_migrate/8ba79fba3cfdc50a96a74cf68ca93a58.png)
Re=5000
![](https://i-blog.csdnimg.cn/blog_migrate/f7fcd46a8fd6dd18f0cb008c57376c9f.png)
![](https://i-blog.csdnimg.cn/blog_migrate/e8ab7279e6a9f9a728b1de7c8366cae2.png)
![](https://i-blog.csdnimg.cn/blog_migrate/fc90bdc93202f4f04e04571aeef9de8d.png)
![](https://i-blog.csdnimg.cn/blog_migrate/64b7eb262af9870b73e152f585fb605a.png)
![](https://i-blog.csdnimg.cn/blog_migrate/12510fb8af05d039f986e1f3acf1a777.png)
Re=7500
![](https://i-blog.csdnimg.cn/blog_migrate/385574ce461d8364fb932be5279bd710.png)
![](https://i-blog.csdnimg.cn/blog_migrate/61f651ce49942eee1cd4aff807ed622f.png)
![](https://i-blog.csdnimg.cn/blog_migrate/329e6a1407da35347756a2aa04cc70d8.png)
![](https://i-blog.csdnimg.cn/blog_migrate/b007c9e7663b4f5acd272ae9de0097f8.png)
(201*201,dt=0.001,Thom公式,80000步)
![](https://i-blog.csdnimg.cn/blog_migrate/21e4e0e2638c263cd06ad94a2d63d90f.png)
Re=10000
![](https://i-blog.csdnimg.cn/blog_migrate/aed44d1f2663ea61a5fbc5f8fc4596fc.png)
![](https://i-blog.csdnimg.cn/blog_migrate/6d041d25f0ad0531d86d594783ffdedf.png)
![](https://i-blog.csdnimg.cn/blog_migrate/b236d306178a8507b20452f13285555f.png)
![](https://i-blog.csdnimg.cn/blog_migrate/9ce99599c6f497dae33b3d96dfaf2d11.png)
![](https://i-blog.csdnimg.cn/blog_migrate/d012fe4cad7a022afcf3be95b8bcd752.png)
Re=15000
![](https://i-blog.csdnimg.cn/blog_migrate/a0df7d9fc14fc06e2885239fe9cd9ad5.png)
分析
- 方法说明:
雷诺数在5000及以下时,采用Thom壁涡公式,101*101网格,dt=0.002;
雷诺数达到7500以上时,如果继续采用原来的方法,发现收敛到错误的结果,如图2-9所示;即使加密网格提高精度(201*201,dt=0.001以保证稳定性),问题依然不变;因此尝试改用Woods壁涡公式,并收敛到了正确的流函数分布。- 与文献中图形进行比较,流线图基本相似,计算结果是可信的。
- Re=100,400,1000时,只在左下和右下两个角落分别出现一个二次涡(旋转方向与主涡相反),大致呈等腰直角三角形,强度随着Re增大而增大,在流线图里的覆盖范围逐渐扩大;Re=3200时,在靠近顶盖的左侧壁面上出现了第3个二次涡;Re=5000时,右下角开始出现三次涡的迹象(见图2-5(右下,放大观察)),Re达到10000后右下角出现明显的三次涡(见图2-7(右下,放大观察))。可以预测当Re增大时会出现更多的三次涡。
- 随着Re的增大,主涡涡心的位置逐渐由右上方向中心移动,但最终位置不会位于方腔正中心,这是由于二维方腔流动为带奇性流动,方腔内部流动不对称所造成的;二次涡和三次涡的涡心位置由角落逐渐向中间移动。
![](https://i-blog.csdnimg.cn/blog_migrate/f2a3cd81e2cb72f78831b650937ae45e.png)
3. Re=10000时流场随时间的演化情况
![](https://i-blog.csdnimg.cn/blog_migrate/d012fe4cad7a022afcf3be95b8bcd752.png)
4. 主涡涡心位置比较(上边界速度设为1,便于比较)
标黄的为我们的计算结果,基准数据来源于文献"High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method"。可见主涡涡心位置的计算是比较准确的,起初位于方腔的右上角,随着雷诺数增大而逐渐接近方腔的几何中心。
(注:雷诺数在7500以上时,改用了Woods壁涡公式,并改用多项式形式的上边界速度条件,比较结果仅供参考)