3　二维非定常流函涡量方程的有限分析算法

如图1，整个水域划分为许多网格，以p为中心的4个相邻网格组成第p个单元，然后对每个单元用有限分析法求解。

3.1　二维非定常流函涡量方程的有限分析解

当二维非定常对流输运方程（7）中取R=Re、H=k时，即为二维非定常流函涡量方程：

对于第p单元，该式按时间离散（图2，其中H（x）、H（y）表示边界上的有关函数），令网格边长分别为Δx=h、Δy=k，时间步长为Δt，p点x、y方向的流速分别为u1p、u2p，邻近节点和中心节点p的流速偏差分别为u1和u2，即邻近节点的u1=u1p+u1，u2=u2p+u2，那么式（8）变为：

对于时间t=tm-1=（m-1）Δt，取Δt，则对t=tm=mΔt上式变为：

又令Re·u1p=2A，Re·u2p=2B，，则上式可写为：

式中：上标m和（m-1）分别为第m和第（m-1）时间步。其单元分析解可表示成：

式中：B、C为边界值。式（12）适用于单元内任意一点，也满足8个边界节点给定的函数值。对于数值计算，只需求p点的值即可。通过流函涡量方程计算p点的分析解，其代数方程为：

式中：n为如图1所示的8个边界节点，即NE、NC、NW、WC、SW、SC、SE和EC；Cn、Cp分别为有关的n个节点和节点p的有限分析系数，从偏微分方程分析解中可求得8个Cn和Cp系数值：

其中

对于非均匀网格，公式略有不同。

3.2　二维流函涡量方程和污染物浓度方程的有限分析法计算方法

对于二维流函涡量方程组（5）的有限分析法迭代步骤如下：

（1）给出j、u、v、k的初值和边界值。

（2）在所有计算网格点上，用九点有限分析格式解式（5-2）的Possion方程，代数方程组用对角矩阵算法求解，为节省时间可引入超松弛因子。

（3）由式（5-3）计算各网格点流速up、v，从而计算出系数：2A=Re·Up，2B=Re·vp。

（4）运用已知的边壁流函数计算壁面涡函数，用九点有限分析格式计算内点上的k值。代数方程组用对角矩阵算法求解，直至收敛。

（5）检验迭代收敛性，如|jm+1-jm|＜X1，且|km+1-km|＜X2，（X2、X2为规定的收敛标准），所得值为收敛解，否则重复（2）～（4）各步，直至达到收敛标准。

求解污染物浓度方程（6），与解式（5）k步骤相同，此时R=Pe，2A=Pe·up，2B=Pe·vd。