二维不可压缩方腔驱动流数值解——基于涡量流函数方程

百宝袋里是数学和代码

已于 2023-07-12 10:00:50 修改

阅读量2.3k

点赞数

分类专栏：计算流体力学文章标签：算法 python

于 2023-06-24 13:28:49 首次发布

本文链接：https://blog.csdn.net/weixin_52046577/article/details/131141030

版权

计算流体力学专栏收录该内容

4 篇文章

订阅专栏

写在前面

该系列文章是笔者的学习记录，旨在通过分享促进交流加深理解，文中难免有许多纰漏与理解不当之处，若能批评指正必定感激不尽。

侧重分享经典CFD方法的具体实现结果，详细的理论推导请参考教材。

一、涡量流函数方程以及求解流程

涡量流函数方程：

$\frac{\partial w}{\partial t} + (\mathbf{v} \cdot \triangledown )w=v\triangle w \\ \\ \triangledown \psi = -w \\\\ u = \frac{\partial \psi}{\partial y} ~~ v= - \frac{\partial \psi}{\partial x}$

边界条件：

$\psi = \psi_{b} ~or~ \frac{\partial \psi}{\partial n}=-v_{\tau}$

在给定 n 时刻涡量分布 $\left \{ w_{i,j}^{n} \right \}$ 的情况下，采取如下的步骤计算出下一时刻的涡量分布 $\left \{ w_{i,j}^{n+1} \right \}$ ：

计算 n+1 时刻内点的涡量（FTCS：空间导数项采用中心差分格式离散，时间导数项采用向前Euler法离散；时间一阶精度，空间二阶精度）
$\frac{w^{n+1}-w^{n}}{\Delta t} + D_{y} \psi^{n}D_{x} w^{n} - D_{x} \psi^{n}D_{y} w^{n} = \nu \triangle _{h} w^{n}$
其中 $D_{x}w_{i,j} = \frac{w_{i+1,j}-w_{i-1,j}}{2 \Delta x} \\ D_{y}w_{i,j} = \frac{w_{i,j+1}-w_{i,j-1}}{2 \Delta y}$ 是标准的中心差分格式， $\triangle _{h}w_{i,j} = \frac{w_{i+1,j}-2w_{i,j}+w_{i-1,j}}{\Delta x^{2}} + \frac{w_{i,j+1}-2w_{i,j}+w_{i,j-1}}{\Delta y^{2}}$ 是离散的五点Laplace算子。
计算 n+1 时刻的流函数，差分方程和边界条件如下（中心差分，二阶精度）
$\left\{\begin{matrix} \triangle _{h}\psi^{n+1} = - w^{n+1} \\ \psi^{n+1} |_{\Gamma} = \psi_{b} \end{matrix}\right.$
为加快解方程组速度，采用泊松方程求解常用的SOR迭代公式 $(\beta = \Delta x / \Delta y)$ ，0<w<2 为松弛因子： $\psi_{i,j}^{t+1} = \psi_{i,j}^{t}+ \frac{w}{2(1+\beta^{2})}[ \psi_{i+1,j}^{t}+\psi_{i-1,j}^{t+1} + \beta^{2}(\psi_{i,j+1}^{t}+\psi_{i,j-1}^{t+1}) + \Delta x^{2} w_{i,j} ^{n+1}-2(1+\beta^{2})\psi_{i,j}^{t} ]$
计算 n+1 时刻边界上的涡量：
一般形式下的Thom壁涡公式 $w_{0} = - \frac{2(\psi_{1}-\psi_{0}+ v_{\tau}h)}{h^{2}}$
假设某壁面切向速度 $v_{\tau}$ ，沿其内法向 $\mathbf{n}$ 有一节点，距离壁面距离为 $h$ ，此点上的流函数值为 $\psi_{1}$ ，如果壁面上的流函数值为 $\psi_{0}$ ，那么在边界是流线的情况下有： $w_{0} = - \frac{2(\psi_{1}-\psi_{0}+ v_{\tau}h)}{h^{2}}$ （注意，这里的切方向 $\tau$ 规定为内法向量 $\mathbf{n}$ 逆时针旋转90°后的方向）
Woods壁涡公式（更高精度） $w_{0}=-\frac{1}{2}w_{1}-\frac{3}{h^{2}}(\psi_{1}-\psi_{0})-\frac{3}{h}v_{\tau}-\frac{3}{2}\frac{\partial v_{n}}{\partial \tau} + \frac{h}{2}\frac{\partial ^{2}v_{\tau}}{\partial \tau^{2}}$ ，含义同上

二、二维方腔流动问题

在单位正方形方腔内，充满粘性不可压缩流体，初始时刻流体静止不动，上壁面以均匀速度1运动,，而其余三个壁面固定不动。由于是粘性流体，将会使整个方腔内的流体运动起来，一定时间后达到稳定状态。

边界条件：
为了避免速度场在角点附近的奇性，可以对上边界给出一个速度分布 $v_{\tau}(x)=16x^{2}(1-x)^{2}$ 其余边界处为无滑移边界条件，即 $v_{\tau}=0$ ；
因为所有的边界都是固壁边界，不妨设边界上流函数 $\psi|_{\Gamma} = 0$ 。

初始条件：
初始时刻流体静止，内点速度均为0；根据初始速度分布计算初始的内点涡量（中心差分）；流函数求解泊松方程得到。

三、计算结果

1. 中间线上的速度分布比较（上边界速度设为1，便于比较）

注：

基准数据来源于文献"High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method"

文献中采用了加密且优化的网格，我们的计算采用的是均分网格。

文献中设置驱动边界上速度等于1，为便于比较，在这一部分将上边界速度设置为1，而不采用 $v_{\tau}(x)=16x^{2}(1-x)^{2}$ 的光滑的速度分布。

Re=100

Re=400

Re= 1000

可以看到，在雷诺数不超过1000时，我们的计算结果是很准确的。

Re=3200

Re = 5000

雷诺数大于3200时，由于文献中采用了加密且优化的网格，我们的计算结果与文献数据有微小偏差，特别是靠近边界处误差较大；但大体趋势依然是可信的。
雷诺数大于7500时，实验发现，采用上边界速度为1的设置会导致无法收敛到正确的流函数分布，因此改用了其他边界速度和边界涡量的设置，不再与文献结果进行比较。

不同雷诺数汇总比较

观察第一幅图竖直中线上的水平速度的分布，可见随着Re增大，边界层越来越薄，Re超过3200后边界层变薄的速率就很小了；高Re时，方腔中部的速度分布几乎是线性的；Re>=3200时，在靠近y=1处u的分布会有一个小的凸起，这也与文献中提到的现象一致。水平中线上的竖直速度的分布也是类似的规律。

2. 流线图和涡旋情况（上边界速度采用多项式形式）

参照：文献中的流线图

Re=100

Re=400

Re= 1000

Re=3200

Re=5000

Re=7500

图2-6（右下,b）
(201*201,dt=0.001,Thom公式,80000步)

Re=10000

Re=15000

图2-8(101*101，dt = 0.001，500000步)（右下，放大观察）

分析

方法说明：
雷诺数在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的增大，主涡涡心的位置逐渐由右上方向中心移动，但最终位置不会位于方腔正中心，这是由于二维方腔流动为带奇性流动，方腔内部流动不对称所造成的；二次涡和三次涡的涡心位置由角落逐渐向中间移动。

3. Re=10000时流场随时间的演化情况

4. 主涡涡心位置比较（上边界速度设为1，便于比较）

标黄的为我们的计算结果，基准数据来源于文献"High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method"。可见主涡涡心位置的计算是比较准确的，起初位于方腔的右上角，随着雷诺数增大而逐渐接近方腔的几何中心。

（注：雷诺数在7500以上时，改用了Woods壁涡公式，并改用多项式形式的上边界速度条件，比较结果仅供参考）