多重网格在黏性流动最小二乘等几何模拟中的应用

陈德祥,徐自力

(西安交通大学机械结构强度与振动国家重点实验室, 710049, 西安)

多重网格在黏性流动最小二乘等几何模拟中的应用

陈德祥,徐自力

(西安交通大学机械结构强度与振动国家重点实验室, 710049, 西安)

针对最小二乘等几何方法模拟黏性流动时条件数大、迭代法收敛速度慢的问题,提出了基于多重网格技术的加速方法。计算中自动生成一系列疏密不同的网格,在最密网格上用最小二乘等几何方法将Navier-Stokes方程离散为代数方程组,用多重网格方法作为独立求解器或共轭梯度法的预处理器迭代求解所得到的代数方程组。对雷诺数为100、400、1 000和2 500的顶盖驱动流进行了数值模拟,计算中进行23次迭代可使方程组的余量降低10个数量级,流动特征量的计算误差在1%以内。计算结果表明,通过多重网格技术加速迭代,提高了最小二乘等几何方法模拟黏性流动的计算效率。

多重网格;最小二乘;等几何方法;共轭梯度法;Navier-Stokes方程

黏性流动现象由Navier-Stokes方程描述,由于非线性导致的复杂性,该方程常用数值方法求解。文献[1]中对线性的Stokes方程提出了一种最小二乘等几何方法,用非均匀有理B样条(NURBS)构造有限元空间,用速度和压力作为独立变量,不增加辅助独立变量。对非线性项进行线性化后,最小二乘等几何方法可推广到Navier-Stokes方程的求解[2]。

文献[1-2]中采用直接方法求解所得的代数方程组,当自由度数增多时直接解法的成本高,需用迭代法求解,这时迭代算法的效率决定了整个程序的计算效率。最小二乘等几何方法进行流动模拟时,所得到的代数方程条件数大,常用的Richardson型迭代和Krylov子空间型迭代随着条件数增大收敛速度降低[3],难以满足要求,因此本文用多重网格技术来提高线性方程组的迭代收敛速度。

多重网格法是目前解代数方程的一种高效迭代算法,Gahalaut等在2013年以泊松方程为对象研究了多重网格在等几何方法中的应用,结果显示,提高NURBS阶数则多重网格的收敛速度降低[4];刘石等将多重网格与共轭梯度法相结合,改进了采用高阶NURBS时的收敛性[5];陈德祥等研究了多重网格在最小二乘等几何方法中的应用,给出了基于离散B样条的网格转换算法[6]。本文根据Navier-Stokes方程的特点对文献[6]中的多重网格算法进行改进,在此基础上建立了多重网格共轭梯度法,通过算例对比研究了多重网格作为独立求解器和作为共轭梯度法的预处理矩阵时的收敛性,最后求解了雷诺数为400、1 000、2 500时的顶盖驱动流。

1 最小二乘等几何方法

定常黏性不可压流动由Navier-Stokes方程描述,其归一化形式可写成

-νΔu+u·u+p=f,x∈Ω

(1)

(2)

式中:u为流动速度;p为压力;f为流体所受到的体积力;ν为流体的黏性系数;Ω为流动区域;x是Ω中的一个点。

式(1)、式(2)分别表示了流动中的动量守恒和质量守恒,求解该方程需要指定合适的定解条件,本文考虑无滑移速度边界条件和平均压力条件

u=g,x∈Γ

(3)

(4)

式中:g为边界上给定的速度;Γ是流动区域的边界。

采用文献[2]中的最小二乘等几何方法进行求解,先用Newton法将控制方程线性化,再对线性化后的方程建立最小二乘变分,最后用NURBS基函数构造有限维空间将控制方程转化为代数方程

Ax=b

(5)

2 基于多重网格技术的代数方程求解

2.1 多重网格算法

多重网格的基本思想是:在密网格上进行光顺,去掉高频误差;在疏网格上进行误差修正,去掉低频误差。在迭代中需要使用疏密不同的多组网格,并在不同网格之间进行转换。对基于NURBS的等几何方法不同尺寸的网格可在计算中用节点插入算法自动生成,不需要事先准备网格,网格转换算子可通过离散B样条得到[6]。由于文献[6]的算法是由较为简单的单变量泊松方程给出,要将其用于Navier-Stokes方程还需在以下2方面进行改进。

首先是建立向量值函数的网格转换算子。以二维问题为例,由于本文基于最小二乘方法,可采用相同的有限维空间近似3个独立变量,将代数方程中的未知量写成如下分块形式

(6)

式中:x1、x2、x3分别是离散后x、y方向的速度以及压力。

对每一个独立变量按文献[6]建立插值算子P′,即对i=1,2,3得

(7)

式中:h、H分别表示密网格和疏网格下的量。因此,总体向量的插值算子为

(8)

根据Galerkin条件可得限制算子为R=PT

(9)

其次是如何满足平均压力为零的约束条件。理论上可以通过指定某个参考点压力值代替约束条件,但这种处理方法会严重降低多重网格迭代的收敛性,主要原因是在参考点处误差光顺效果降低。Brandt指出:不必在密网格上处理这种约束条件,只需将约束条件的余量转换到疏网格上并进行误差修正即可[7],为此将式(4)离散成如下代数方程

aTx3=0

(10)

如果当前网格为最疏网格(l=1),那么用直接方法求解,对压力叠加一个常数使它满足压力约束条件。如果当前网格不是最疏的(l>1),先进行μ1次光顺迭代,通过限制算子将余量转换到下一层疏网格上(l-1层),并计算压力约束条件的余量,然后在l-1层调用γ次多重网格算法求解余量方程,最后通过插值算子将修正量转换到l层,进行μ2次光顺迭代。

IFl=1

ELSE

END

2.2 多重网格共轭梯度法

共轭梯度法是一种Krylov子空间迭代方法,它适用于系数矩阵为对称正定的情况,而基于最小二乘法的系数矩阵总是对称正定的。理论上若代数方程的维数为n,那么进行n次共轭梯度法迭代将收敛到精确解;实际计算中共轭梯度法的收敛速度快慢取决于所采用的预处理矩阵能否聚集系数矩阵的特征值[3]。多重网格方法也可看成用简单迭代方法求解某个预处理矩阵作用下的代数方程,这个预处理矩阵能够有效聚集矩阵特征值。

假设总共k层网格,如果把多重网格迭代看成静态迭代方法,那么它可写成如下矩阵形式

(11)

(12)

虽然Bk的表达式包含了Ak的逆矩阵和复杂的递归矩阵Mk,无法直接生成,但是它仍然可用作共轭梯度法的预处理矩阵。这是由于共轭梯度法只要求预处理矩阵对称正定,能快速计算出与任意向量的乘积,并不需要生成预处理矩阵。当μ1=μ2时,采用Jacobi光顺或对称Gauss-Seidel光顺的V型和W型多重网格算法都满足对称正定要求[9]。由式(11)可知:Bk与任意向量v的乘积Bkv可快速计算出,实际上只要以0向量为初始值,以v为右端项,应用算法1进行一次多重网格迭代即可。算法2是多重网格共轭梯度法伪代码,它与一般共轭梯度算法的不同是在预处理部分调用了算法1的多重网格迭代。

算法2:MGCG(Ak,bk,x,k,γ,μ1,μ2)

r=bk-Akx

z=MG(Ak,r,0,k,γ,μ1,μ2)

p=z

WHILE not converges

τ=rTz

w=Akp

α=τ/wTp

x=x+αp

r=r-αw

z=MG(Ak,r,0,k,γ,μ1,μ2)

β=rTz/τ

p=z+βp

END

3 数值算例分析

考虑顶盖驱动流问题,流动区域为单位正方形,两侧及下壁面固定,上壁面的运动速度u=(-1,0)。该问题在文献中有可信的高精度计算结果[10-12],常作为基准问题用来考核新CFD算法的精确性。本文用Matlab编制了计算程序,分别求解了雷诺数从100到2 500时的流动,并与文献中的结果进行对比,所有数值结果都是归一化的。

首先考虑低雷诺数下(Re=100)的流动。多重网格采用5层网格的W型循环,最密网格为32×32,最疏网格为2×2,参数μ1、μ2均取为2,光顺算法为对称Gauss-Seidel迭代。计算中存在嵌套的两层迭代,外层为Newton迭代,它的收敛条件是相邻两次迭代结果的相对误差小于10-8;内层为线性方程的迭代,它的收敛标准是方程组的余量小于10-10。

图1表示Re=100时多重网格法(MG)和多重网格共轭梯度法(MGCG)的收敛过程,图中k是NURBS基函数的连续性指标,从最高点到最低点表示某个Newton迭代步下的线性方程组的收敛过程。对该问题总共需要7次Newton迭代,从图中可以看出,随着Newton迭代的推进,线性方程组求解所需的迭代次数逐渐减少。这是因为用前一Newton步的结果作为线性方程组迭代的初始值,随着Newton迭代的收敛,线性方程组的初始余量也逐渐降低,表现为图中的高点逐渐减低。比较图1a、图1b可以看出:相同k下的MGCG迭代比MG迭代收敛速度快,k=2时,进行23次MGCG迭代方程组的余量降低了10个数量级,k=1时只需14次迭代;当k增大时,MG迭代和MGCG迭代的收敛速度都降低,但是MG迭代的效率降低更多,k从1增大到2,MG所需的迭代次数加倍,而MGCG所需的迭代次数只增加约1/3。

(a)MG算法 (b)MGCG算法

图2是不同k值下垂直中心线上的x方向速度以及水平中心线上的y方向速度分布,作为对比参照的是Marchi等[10]用1 024×1 024网格得到的结果。本文k=1,2时的总自由度数分别为3 468和3 675,从图2中可以看出,k=2给出的速度分布更精确,当k值较大时,最小二乘等几何方法是一种高阶方法,网格数不变的条件下,数值解的误差将随k的增加呈指数速度下降[1]。

从上面k对多重网格迭代收敛性和解的精度影响分析可知,较大的k值能给出更好的结果,MGCG对k值的敏感性较低,因此基于精度和成本两方面考虑,计算中k应取较大值,而线性方程组应用MGCG迭代求解。

图3是Re=10时多重网格共轭梯度法的迭代收敛性过程,k=2、网格数及多重网格参数与前面相同。比较图3与图1b可发现,每一个Newton步下Re=10比Re=100需要更多的迭代才能收敛,在文献[13]中也发现类似的现象,即对最小二乘法在很低雷诺数下用多重网格求解线性方程组不能提高计算效率。这是由于系数矩阵的条件数与雷诺数有关,对最小二乘等几何方法,Re=10时的条件数大于Re=100时的条件数。

(a)垂直中心线上速度u (b)水平中心线上速度v

图3 多重网格共轭梯度法收敛过程(Re=10)

对中、高雷诺数流动分别计算了Re=400,1 000,2 500的情况,由于流动的复杂性随着雷诺数的增加而增加,雷诺数为400、1 000时采用64×64网格,雷诺数为2 500时采用128×128网格。

在雷诺数较大的情况下,用雷诺数递进技术来解决Newton法收敛半径小的问题,先以零为初值求解Re=400的流动,再以所得到的结果为初值增大雷诺数进行迭代计算,逐渐递增到指定的雷诺数。Newton迭代收敛标准是相邻两次迭代结果的相对误差小于10-10,MGCG迭代收敛标准是方程组的余量小于10-8。

取k=6进行计算,图4是Re=1 000时流函数、涡量和压力的分布图,可以看到流动中存在一个主涡,在左、右下角存在两个次涡,图中的分布趋势与文献[12]给出的分布图一致。表1和表2分别是主涡和左下角次涡的参数(流函数ψ、涡量ω、涡中心坐标x,y)随k增大过程中的收敛情况,与文献[12]中的结果对比可以看出,k=6所得的数值解具有较高的精度,相对误差在1%以内。

(b)涡量

(c)压力

表1 主涡的收敛性(Re=1 000)

表2 左下角次涡的收敛性(Re=1 000)

表3 全局收敛性(Re=1 000)

图5和图6分别是垂直中心线上x方向速度分布和水平中心线上y方向速度分布,其中Marchi等[10]采用了1 024×1 024网格,Erturk等[11]采用了600×600网格,可以看出,对于Re=400,1 000,2 500时本文方法得到的速度分布与文献中的高精度计算所得结果一致。

图5 不同雷诺数下垂直中心线上x方向速度分布

图6 不同雷诺数下水平中心线上y方向速度分布

4 结 论

多重网格技术提高了最小二乘等几何方法模拟黏性流动的计算效率。多重网格迭代过程中用密网格使高频误差衰减掉,用疏网格使低频误差衰减掉,最终线性方程组的余量在迭代中快速衰减,将多重网格法与共轭梯度法相结合进一步提高了迭代收敛速度,进行23次迭代可使方程组的余量降低10个数量级。

流场的计算精度随NURBS基函数的连续性指标k的增大而提高。对雷诺数为100、400、1 000和2 500的顶盖驱动流进行了数值模拟,k=6时流动特征量的计算误差在1%以内,k增大时多重网格共轭梯度法比多重网格法具有更好的收敛性。

[1] 陈德祥, 徐自力, 刘石, 等.求解Stokes方程的最小二乘等几何分析方法 [J].西安交通大学学报, 2013, 47(5): 51-55.

CHEN Dexiang, XU Zili, LIU Shi, et al.Least squares isogeometric analysis for Stokes equation [J].Journal of Xi’an Jiaotong University, 2013, 47(5): 51-55.

[2] CHEN Dexiang, XU Zili, LIU Shi, et al.Least squares finite element method with high continuity NURBS basis for incompressible Navier-Stokes equations [J].Journal of Computational Physics, 2014, 260: 204-221.

[3] SAAD Y.Iterative methods for sparse linear systems [M].Philadelphia, USA: SIAM, 2003.

[4] GAHALAUT K P S, KRAUS J K, TOMAR S K.Multigrid methods for isogeometric discretization [J].Computer Methods in Applied Mechanics and Engi-neering, 2013, 253: 413-425.

[5] 刘石, 陈德祥, 冯永新, 等.等几何分析的多重网格共轭梯度法 [J].应用数学和力学, 2014, 25(6): 630-639.

LIU Shi, CHEN Dexiang, FENG Yongxin, et al.A multigrid preconditioned conjugate gradient method for isogeometric analysis [J].Applied Mathematics and Mechanics, 2014, 25(6): 630-639.

[6] 陈德祥, 窦柏通, 徐自力.最小二乘等几何分析的多重网格方法 [J].西安交通大学学报, 2014, 48(7): 124-130.

CHEN Dexiang, DOU Baitong, XU Zili.Multigrid least squares isogeometric analysis [J].Journal of Xi’an Jiaotong University, 2014, 48(7): 124-130.

[7] BRANDT A, LIVNE O E.Multigrid techniques: 1984 guide with applications to fluid dynamics [M].Philadelphia, USA: SIAM, 2011.

[8] TATEBE O.The multigrid preconditioned conjugate gradient method [C]∥The Sixth Copper Mountain Conference on Multigrid Methods.Copper Mountain, USA: NASA, 1993: 621-634.

[9] TROTTENBERG U, OOSTERLEE C W, SCHULLER A.Multigrid [M].London, UK: Academic Press, 2000: 28-50.

[10]MARCHI C H, SUERO R, ARAKI L K.The lid-driven square cavity flow: numerical solution with a 1 024×1 024 grid [J].Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2009, 31(3): 186-198.

[11]ERTURK E, CORKE T C, GOKCOL C.Numerical solutions of 2D steady incompressible driven cavity flow at high Reynolds numbers [J].International Journal for Numerical Methods in Fluids, 2005, 48(7): 747-774.

[12]BOTELLA O, PEYRET R.Benchmark spectral results on the lid-driven cavity flow [J].Computers & Fluids, 1998, 27(4): 421-433.

[13]RANJAN R, REDDY J N.On multigrid methods for the solution of least-squares finite element models for viscous flows [J].International Journal of Computational Fluid Dynamics, 2012, 26(1): 45-65.

[14]BRUNEAU C H, SAAD M.The 2D lid-driven cavity problem revisited [J].Computers & Fluids, 2006, 35(3): 326-348.

(编辑 赵炜)

AnApplicationofMultigridinViscousFlowSimulationswithLeastSquaresIsogeometricMethod

CHEN Dexiang,XU Zili

(State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University, Xi’an 710049, China)

A multigrid technique to accelerate iteration is proposed to solve the problem of low convergence rate due to large condition number in simulating viscous flow with the least squares isogeometric method.A series of grids with different element sizes are automatically generated in analysis.The least squares isogeometric method is used to discretize the Navier-Stokes equations into a system of linear equations on the densest grid, and the system is iteratively solved using either the multigrid method as standard solver or the preconditioner of the conjugate gradient method.Numerical simulations are performed on lid driven cavity flow withRe=100, 400, 1 000 and 2 500.The residual of algebraic equations is reduced about 10 order of magnitude after 23 iterations, and the calculation error on characteristic parameters of the cavity flow is less than 1%.The results show that the performance of the least squares isogeometric method for viscous flow is improved when the multigrid acceleration technique is used.

multigrid; least squares; isogeometric method; conjugate gradient method; Navier-Stokes equations

2014-05-19。

陈德祥(1979—),男,博士生;徐自力(通信作者),男,教授,博士生导师。

国家“973计划”资助项目(2011CB706505)。

时间:2014-08-13

10.7652/xjtuxb201411021

O241.82;O357.1

:A

:0253-987X(2014)11-0122-06

网络出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20140813.1008.003.html