对纳维斯托克斯方程的隐式速度解耦过程2.1 介绍随着直接数值模拟和大涡模拟这两种数值模拟方法越来越进步，出现了很多求解不可压缩纳维斯托克斯方程的有效数值算法，其成功的核心是对耦合的不可压缩动量方程和连续性方程解耦。

对文献的精读发现，之前的许多方法使用了半隐式的方案，即把隐式方案应用于粘性条件，显式方案应用于非线性对流条件，时间步通过CFL 数控制。

Choi 和Moin 在分步法的基础上采用了一个完全隐式的方法，首先对纳维斯托克斯方程在时间上离散，然后进行空间离散，用这种方法得到的中间速度分量是耦合的，后来使用牛顿迭代方法得到了中间速度分量。

为了防止迭代过程，Rosenfeld 提出了一个非耦合的隐式解算器[8]，他设计了三个时间步的线性化方案，这个方案需要n-1步和n 步的速度场来得到n+1步的速度，在不忽略时间二阶精度和稳定性的基础上，控制方程被解耦。

在最近的研究中，Kyoungyoun Kim ，Seung-Jin Baek 和Hyung Jin Sung[9] 对解决不可压缩湍流的纳维斯托克斯方程发展出了一种有效的数值计算方法，这个算法提出了一个新的隐式速度解耦过程，采用完全隐式的时间推进，在块LU 分解和近似分解的基础上，速度项和压力项被解耦，同时保留了时间二阶精度，另外，由于隐式的对流条件，中间速度是耦合的，所以重点放在了对中间速度的解耦上，这就需要对第n 个时间步的速度近似分解，这些解耦过程同样保留了时间二阶精度。

本文数值模拟的过程中也用到了这种解耦过程，第二部分将会对目前的数值解耦方法及方程的近似分解过程做一个简要的介绍，在第三部分，将会把这个解耦过程应用于槽道流，并用直接数值模拟进行验证，画出结果进行比较分析。

2.2 数值方法不可压缩的无量纲纳维斯托克斯方程为：1,(1,2,3)Re i ii j j i j ju u p u u i t x x x x ∂∂∂∂∂+=-+=∂∂∂∂∂ (1) 0iiu x ∂=∂ (2)其中，i x 是笛卡尔坐标，i u 是每一个方向相应的速度分量，Re 是雷诺数，所有的分量都用特征长度进行了无量纲化。

在第n+1/2个时间步上，对上述两个方程进行空间和时间离散，方程可以写成如下的形式：111/2111(()())()22Ren n n n n n n u u H u H u Gp Lu Lu mbc t ++++-++=-+++∆(3)10n Du cbc +=+ (4)其中，L 代表离散的拉普拉斯算子，H 代表离散对流算子，G 代表离散梯度算子，D 代表离散散度算子。

必须说明的是，边界条件已经被合并到方程(3)和(4)中，空间离散算子L 、H 、G 和D 定义在交错网格上，并且用二阶中心有限差分来定义，变量1n u +和1/2n p +定义在内部节点上，而不是交错网格上的边界点，边界上已知的速度分量和速度的边界条件被分别加在了mbc 和cbc 上。

在交错网格上，压力被定义在立方体中心，速度分量被定义在正交平面上，并且动量方程的离散化需要边界上的压力或压力梯度。

对于时间离散，使用了一个完全隐式的时间离散，因为对流条件和粘性条件是隐式的，所以产生了非线性的方程。

在最近的研究中，非线性条件被线性化，同时保留了时间二阶精度：11112()n n n n n n n ni j i j i j i j u u u u u u u u O t ++++=+-+∆ (5)通过这个线性化，对流项的线性算子N 可以被定义为：111(()())2n n n Nu H u H u ++=+ (6) 注意到线性化算子N 中包含第n 步的速度。

通过使用对流算子N ，离散方程(3)和(4)可以以矩阵的形式写出：100n A G r mbc u D cbc p δ+⎛⎫⎛⎫⎛⎫⎛⎫=+ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎝⎭ (7)11[()]2Re A I t N L t =+∆-∆ 1/2112Ren n n r u Gp Lu t -=-+∆ 1/21/2n n p p p δ+-=-在上述公式中，通过反演方程(7)的系数矩阵，得到了下一个时间步上的值1n u +和1/2n p +，因为(7)的系数矩阵大而且稀疏，所以方程(7)的求解过程很难，换句话说，因为在动量方程和连续性方程中，速度和压力是相关的，所以方程(7)不能被直接求解。

对方程(7)使用块LU 分解：1000n A I tG r mbc u D tDG I cbc p δ+⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫∆=+ ⎪ ⎪ ⎪ ⎪ ⎪-∆⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭ (8)这个方程和(7)是不同的，作为近似因子分解结果，把G p δ近似成了tAG p δ∆。

在最近的研究中，压力p 通过p δ来表达。

上面近似分解的误差为：2()0tMG p O t δ⎛⎫∆∆=⎪⎝⎭，其中12Re M N L =- (9) 方程(8)也可以被写为：*00A r mbc u D tDG cbc p δ⎛⎫⎛⎫⎛⎫⎛⎫=+ ⎪ ⎪ ⎪ ⎪-∆⎝⎭⎝⎭⎝⎭⎝⎭ (10)1*0n I tG u u I p p δδ+⎛⎫⎛⎫⎛⎫∆=⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭(11) 其中，*u 是1n u +的中间速度，通过以下操作，可以得到更加简化的表达式：*Au r mbc =+ (12)*tDG p Du cbc δ∆=- (13)1*n uu tG p δ+=-∆ (14)1/21/2n n p p p δ+-=+ (15) 因为在方程(3)和(4)中已经应用了边界条件，所以不需要对边界条件进行其他特殊处理，并且在方程(12)—(15)解耦速度和压力的过程中，同样保留了时间二阶精度。

然后，通过使用*u δ，将近似分解用于速度分量*u 上，所以方程(12)可以被写为：*n A u Au r mbc δ=-++ (16) **n u u u δ=- (17)方程(16)可以以矩阵形式给出：*11112131*21222322*313233331u I tM tM tM R tM I tM tM u R t tM tM I tM R u δδδ⎛⎫⎛⎫⎛⎫+∆∆∆ ⎪ ⎪ ⎪∆+∆∆= ⎪ ⎪ ⎪∆ ⎪ ⎪ ⎪ ⎪∆∆+∆⎝⎭⎝⎭⎝⎭(18)当动量方程通过半隐式方法离散时，方程(18)中的非对角线子阵,i j M （i j ≠）为零，然而，在最近的全隐式方法中，,()i j M i j ≠不再为零，因为非线性的对流条件为隐式条件，这就意味着*1u δ、*2u δ和*3u δ是完全耦合的。

通过对系数矩阵(18)进行近似分解，中间速度分量可以用第n 个时间步的速度进行解耦，在近似分解的过程中，时间二阶精度被保留了。

*11112131*21222322*333132330010000u I tM tM tM R I tM I tM I tM u R t I tM tM tM I R u δδδ⎛⎫⎛⎫⎛⎫+∆∆∆⎫⎛ ⎪ ⎪ ⎪⎪ ∆+∆∆= ⎪ ⎪ ⎪⎪ ∆ ⎪ ⎪ ⎪⎪ +∆ ⎪∆∆⎭⎝⎝⎭⎝⎭⎝⎭ (19) 像之前所描述的，当前的速度解耦过程只需要第n 个时间步上的速度，在三个时间步的方案中，第n-1步和第n 步的速度都需要才能保证时间的二阶精度。

方程(19)的二阶误差项如下：**11122111332***211222113322233***311223113332233()tM M u tM M u O t tM M u tM M u tM M u tM M u tM M u tM M u δδδδδδδδ⎛⎫∆+∆ ⎪∆=∆+∆+∆ ⎪⎪ ⎪∆+∆+∆⎝⎭引入新的变量**1u δ、**2u δ后，*u 可以通过下面的式子得到：**11111()I tM u R tδ+∆=∆ (21) ****22222111()I tM u R M u tδδ+∆=-∆ (22) ******33333113221()I tM u R M u M u tδδδ+∆=--∆ (23) ****22233u u tM u δδδ=-∆ (24)*****11122133u u tM u tM u δδδδ=-∆-∆(25) **,(1,2,3)n i i i u u u i δ=+= (26)由以上的式子可以看到，不用求解方程(12)中的大矩阵，就可以得到中间速度。

小结：①首先用(21)—(26)的速度解耦过程得出*u ； ②从方程(13)中解出p δ；③从方程(14)计算得出下一步的速度1n u +。

2.3 直接数值模拟的验证使现在的解耦方法保留时间二阶精度是很重要的，首先将2.2中介绍的程序应用到三维周期性的槽道流中，控制方程为不可压缩的纳维斯托克斯方程，边界条件为壁面无滑移条件，建立针对控制方程的数值离散化方法，本程序中使用了有限差分法[18]：Crank-Nicolson method ，在解耦过程中保持了恒定的质量流率，雷诺数的值是4400，计算域为:02,02,0X Y Z ≤≤∏≤≤≤≤∏,为了使计算结果更加准确，所取的网格点为：3128，0.01t ∆=。

下图为得到的壁面应力随时间的变化情况：图2.1 直接数值模拟壁面应力图图 2.1为用直接数值模拟计算的壁面应力随时间的演化过程。

由图可以看出，在无量纲时间400之后，数据趋于稳定，所以在计算流向平均速度、脉动和雷诺应力时，从40000步开始取数据，接下来画出流向速度随壁面高度的变化，在画图过程中不仅在同一水平高度进行了平均，而且将稳定后的各个时间段上的值做了平均[10]，因为展向速度和法向速度为零，所以没有画出，横纵坐标都经过了无量纲化，下图中实线表示实际模拟计算得到的结果，而虚线代表Lee 和Moser[11]得到的结果：图2.2 直接数值模拟流向速度图图2.2为用直接数值模拟计算的流向速度随时间的演化过程。

可以看出，计算得到的结果与正确结果基本重合，可见，在直接数值模拟中流向速度得到了很好的验证。

下图为脉动''u u <>、''v v <>、''w w <>和雷诺应力''u v <>随壁面高度的变化，在平均上与流向速度的处理方式一样，实线代表实际模拟计算得到的结果：图2.3 直接数值模拟流向脉动图 图2.4 直接数值模拟法向脉动图图2.5 直接数值模拟展向脉动图图2.6 直接数值模拟雷诺应力图图2.3、图2.4、图2.5、图2.6分别为用直接数值模拟计算的流向、法向、展向脉动和雷诺应力随时间的演化过程。

由上面四幅图同样可以看出，计算得到的结果与Lee和Moser的结果基本重合，所以，在直接数值模拟中脉动和雷诺应力也得到了很好的验证。