改进的格子涡方法及其在混合层模拟中的应用

　　根据获得流场速度途径的不同,离散涡方法可以分成Lagrangian涡方法和格子涡方法。La-grangian涡方法通过离散涡元的诱导获得流场速度,并已在许多流动的模拟中取得了成功[1～3]。但其缺点是计算量与离散涡元数量N的平方成正比,N很大时,计算量难以承受,N太少时,又难以达到必要的精度和流场分辨率。与Lagrangian涡方法不同,格子涡方法的实质是一种混合涡方法,即先在Euler框架内通过求解流函数的Poisson方程计算速度场,再在Lagrangian框架内追踪离散涡元的运动。在Euler框架内采用椭圆方程的快速解法求解流函数的Poisson方程可将每个时间步求解流场速度的计算量从O(N2)降至O(N+MlbM),其中,M为网格节点数。因此,格子涡方法可以有效地克服Lagrangian涡方法的上述缺陷。但由于流动过程的强非线性,在采用插值方法从网格节点处的流场速度计算离散涡元的速度时,往往会导致很大的误差[4]。为此,文中对格子涡方法进行了改进,使计算离散涡元速度时无需插值,避免了由此产生的误差。

　　对混合层流动的模拟结果表明本文对格子涡方法的改进是成功的。

　　1　数值方法

　　1.1　改进的格子涡方法

　　二维不可压缩流动的流函数Poisson方程为

　　其中:φ为流函数,ωΔ为涡量。计算时,首先在计算域内划分网格,然后用面积加权法把所有离散涡元携带的涡量分配到相应的网格节点上,并求解式(1)。面积加权法如图1所示。

　　对于强度为#n的离散涡元,其所携带的涡量按照下式分配到周围的4个网格节点上,

　　其中:S_k为第n个涡元在一个网格内的第k个节点上的面积权重,S为网格面积。

　　在通过求解式(1)得到每个网格节点上的流函数φij之后,用下式给出所有网格节点上的速度:

　　在通常的格子涡方法中,由网格节点处的流场速度,用下式内插求得各涡元的速度,并通过对涡元速度的积分求得所有涡元在下一时刻的位置:

　　但由于流动过程的强非线性,在采用式(4)内插求涡元速度时会引入很大的误差[4]。为此,本文对格子涡方法进行了改进。在每个时间步开始时,将所有的离散涡元置于网格节点位置上,每个涡元的强度为其所在网格节点处的涡量与网格面积的乘积XijS,在通过求解式(1)和(3)得到流场速度后,由于所有的离散涡元都是位于网格节点上,因此涡元可以直接以所在网格节点的速度运动至下一位置,而无需通过式(4)进行内插,避免了由此带来的误差。

　　在下一时间步开始时,由于所有离散涡元都已运动至下一位置,所以需要再采用式(2)将已脱离网格节点位置涡元所携带的涡量重新分配到网格节点上,并用位于网格节点位置,且强度为XijS的涡元取代偏离网格节点位置的涡元。