通过把原始变量（速度，压强）转化为导出变量（涡量，流函数）的方程组，可以避开处理压强时的困难，称为涡量流函数方法。

涡量和流函数

对于二维的流动，记速度 \(\mathbf{u}=(u,v)^T\)，则涡量 \(w\) 只是一个标量，定义为

\[ w = \frac{\partial u}{\partial y}-\frac{\partial v}{\partial x} \]

对于不可压缩流动，有流函数 \(\varphi\)，满足：

\[ \begin{aligned} u &= \frac{\partial \varphi}{\partial y}\\ v &= -\frac{\partial \varphi}{\partial x} \end{aligned} \]

有的材料中，流函数定义相比上文都多一个负号，这不影响算法的本质。流函数的存在性其实就直接保证了不可压缩条件的成立：

\[ \nabla \cdot \mathbf{u} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = \frac{\partial^2 \varphi}{\partial x \partial y} - \frac{\partial^2 \varphi}{\partial x \partial y}=0 \]

流函数的等值线就是流线，亦即流体的速度场 \((u,v)\) 与流函数的梯度场 \((\varphi_x,\varphi_y)\) 处处垂直。流函数和涡量满足如下关系：

\[ w = \frac{\partial u}{\partial y}-\frac{\partial v}{\partial x} = \nabla^2 \varphi \]

涡量流函数方程组

从无量纲化的二维不可压缩方程组出发

\[ \begin{aligned} \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot \nabla) \mathbf{u} &= -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u}\\ \nabla \cdot \mathbf{u} &= 0\\ \end{aligned} \]

连续性方程直接被流函数 \(\varphi\) 的存在性取代，对动量方程按照分量形式有：

\[ \begin{aligned} \frac{\partial u}{\partial t} + u u_x + v u_y &= -p_x + \frac{1}{Re}\nabla^2 u\\ \frac{\partial v}{\partial t} + u v_x + v v_y &= -p_y + \frac{1}{Re}\nabla^2 u\\ \end{aligned} \]

分别对 \(y\) 求导，对 \(x\) 求导，相减得到涡量方程

\[ \frac{\partial w}{\partial t} + u w_x + v w_y = \frac{1}{Re} \nabla^2 w \]

压强消失了，不可压缩条件也消失了，得到如下方程组：

\[ \begin{aligned} \frac{\partial w}{\partial t} + u w_x + v w_y &= \frac{1}{Re} \nabla^2 w\\ w &= \nabla^2 \varphi \end{aligned} \]

算法流程

对于得到的涡量流函数方程组，在不考虑边界的情况下，计算流程是容易的，例如可以使用如下流程进行计算：

  • 已知 \(w^n,\varphi^n\)，通过 \(u = \varphi_y\)，\(v=- \varphi_x\) 得到速度 \(u^n,v^n\)；
  • 通过涡量方程计算得到 \(w^{n+1}\)；
  • 代入 \(\nabla^2 \varphi^{n+1} = w^{n+1}\)，计算新时刻的流函数 \(\varphi^{n+1}\)；

如果需要计算压强，可以利用动量方程的分量形式，依次对 \(x\) 求导，对 \(y\) 求导，相加得

\[ \nabla^2 p = -\left[ (u_x)^2 + 2 u_y v_x + (v_y)^2 \right] = 2\left[ u_x v_y - u_y v_x \right] \]

上述涡量—流函数方法的优缺点都非常明显：

  • 优点：消去了压强项，消去了不可压缩条件，对方程组的处理更加容易；
  • 缺点：涡量 \(w\) 的边界条件难以确定，不易从二维推广到三维。