 International Journal of Fluid Dynamics 流体动力学, 2013, 1, 15-20 http://dx.doi.org/10.12677/ijfd.2013.12003 Published Online June 2013 (http://www.hanspub.org/journal/ijfd.html) The Planar Element with Eight Nodes for Solving Incompressible Navier-Stokes Equations Xia oju a n Wei 1, Liling Liu2, Lingzhi Xie3, Liqing Yi3, Yongtao Wei3* 1Deparment of Power Engineering, Sichuan Electric Vocational and Technical College, Chengdu 2School of Materials Science and Engineering, Southwest Jiaotong University, Chengdu 3Department of Applied Mechanics, College of Architecture and Environment, Sichuan University, Chengdu Email: *wyt2119@hotmail.com Received: Feb. 26th, 2013; revised: Mar. 11th, 2013; accepted: Apr. 16th, 2013 Copyright © 2013 Xiaojuan Wei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract: The definition of the element length and the inverse estimate constant of straight-edged quadrilateral element with eight nodes for solving incompressible Navier-Stokes equations was proposed, based on this, the stabilization fac- tor was readily calculated. The element length was only dependent on the quadrilateral geometry and easy to compute, and the largest inverse estimate constant among various quadrilateral was obtained by means of the optimization method of genetic algorithm (GA). The numerical solutions of the classical lid-driven cavity flow problem were ob- tained for Reynolds number of 20,000 which demonstrated the validity of the element length and the inverse estimate constant proposed. Keywords: Galerkin/Least Squares FEM; Stabilization Factor; Element Length; Inverse Estimate Constant 求解不可压缩 Navier-Stokes 方程的 GLS 平面八节点单元 魏晓娟 1，刘力菱 2，谢凌志 3，易丽清 3，魏泳涛 3* 1四川电力职业技术学院动力工程系，成都 2西南交通大学材料科学与工程学院，成都 3四川大学建筑与环境学院应用力学系，成都 Email: *wyt2119@hotmail.com 收稿日期：2013 年2月26 日；修回日期：2013年3月11 日；录用日期：2013年4月16 日 摘 要：给出了求解不可压缩Navier-Stokes 方程的GLS 平面八节点单元的“单元长度”的定义和相应的逆估计 常数，由此可计算出该单元的 GLS稳定因子。给出的“单元长度”只依赖单元形状且易于计算；各种形状的直 边四边形单元的逆估计常数的最大值由遗传算法得出。对 Reynolds 数为 20,000 的方腔上盖板流的数值模拟表明， 本文所提出的单元长度和逆估计常数的正确性。 关键词：Galerkin/最小二乘有限元；稳定因子；单元长度；逆估计常数 1. 引言 经典Galerkin有限元求解不可压缩Navier-Stokes 方程的主要困难在于压力解的数值不稳定和速度解 的寄生振荡。前者的原因在于不可压缩条件，解决方 法是对压力和速度采用满足Babuska-Brezzi条件[1,2]的 不等阶插值；后者在高Reynolds数下尤为明显，它源 于经典Galerkin格式对对流项的处理。 *通讯作者。 为实现不可压缩N-S方程的有限元分析，Brooks Copyright © 2013 Hanspub 15  求解不可压缩 Navier-Stokes方程的 GLS 平面八节点单元 和Hughes建立了流线迎风Petrov-Galerkin格式(SUPG) 有限元[3]；Franca和Frey 利用Galerk in/ 最小二乘(GLS) 方法建立了不可压缩N-S方程的稳定化有限元格式[4]； Burda 等人为考虑不可压缩条件的最小二乘项，提出 了semi-GLS格式[5]。GLS方法的优点在于，它消除了 对流项引发的数值不稳定性，却又无需引入过大的数 值耗散；而且还绕开了Babuska-Brezzi条件，允许速度 和压力的等阶插值。 确定合适的稳定因子是GLS方法研究中的一个重 要内容。文献[4]将其表示成逆估计常数和“单元长度” 的形式，而这二者又取决于单元形状和速度插值阶 数。Harari和Hughes给出了具有中节点的直边三角形 和矩形单元的逆估计常数和“单元长度”的定义[6]。 Franca和Madureria通过求解与逆常数估计相对应的广 义特征值问题来计算稳定因子[7]；而Tezduyar和Osawa 提出了由单元矩阵的计算方法[8]。 本文主要工作是：对中节点均匀分布的平面八节 点直边四边形单元，给出了逆估计常数和对应的“单 元长度”，后者完全依赖于单元形状，且易于计算。 该工作目前尚未见有公开报道。在此基础上，通过模 拟高 Reynolds 数下的方腔上板流，并于Erturk 等人工 作[9]进行了对比，数值结果验证了本文方法的有效性。 2. 不可压缩 N-S 方程的 GLS 稳定化 有限元方法 2.1. 稳定化的 GLS 列式 考察在边界为 的空 间域内的定常不可压缩N-S方程 , gf gf 2p uu ub in (1) 0 u in (2) ug on g (3) nf on f (4) 其中 是密度， u是速度， 是动力粘性系数， 是 压力， 是单位体积内的体力，和是给定函数， 是 p bg f n f 的单位外法线， 是应力张量 2p Iεu (5) 其中 是二阶单位张量，是形变速率张量。 I εu 设 和 分别是速度和压力的权函数，则方程 (1)和(2)的GLS 稳定化有限元格式为[4] δuδp el el 2 1 2 LSIC 1 δdδ:d δdδd δ δδ d δd0 f e e n e n e p p p uuub uf u uu u uu ub uu (6) δ 是与对应的 形变速率张量， δu 是GLS稳定因子， 2 LSIC u是针对不可压缩条件的最小二乘项的稳 定因子。 2.2. 有限元列式和 Newton-Raphson 迭代 为简化推导，假设整个求解域 只被划分成一个 单元。将 和表示成对应单元节点量的插值 函数形式 ,δ,puu δp ,δδ,,δpp u NUu NUNPNP (7) N是插值函数矩阵，如前所述，在GLS 方法中速度 和压力可以采用等阶插值。 将式(7)代入式(6)，并考虑到和 的任意性， 导出半离散化的GLS 稳定化有限元列式： δUδP TT TT TT LSIC dd ˆˆd d NLuBT WLuNUHPb BmmBUF (8) TT T d 1 ˆˆˆd0 NmBU b HLu NUHP (9) F 是外载荷的等效节点力 T d f T d F Nb N f (10) 为使用 Newton-Raphson 方法求解非线性方程(8) 和(9)，定义切线矩阵如下 TT T TT LSIC ddd d ˆd d UU KNGLNB U WGLNN BmmB CB (11) Copyright © 2013 Hanspub 16  求解不可压缩 Navier-Stokes方程的 GLS 平面八节点单元 T dˆ dd d UP Ψ KBmNW P T H (12) TT T d d ˆ d PU KU NmBHG LNN d (13) T dˆˆd d PP KH P H (14) d和d的上划线表示未考虑 dLSIC ,d 和 对 T dW切 线矩阵的贡献。迭代格式如下 00 1 1 11 11 0, 0 nn nn UU UP nn nn PU PP nnn nnn UP KK U KK P UU U PP P (15) 式(11)~(14)中的各矩阵在二维下的具体形式，见文献 [10]。 3. 稳定因子 3.1. 逆估计 式(6)中的 GLS 稳定因子 给定如下[4]： 12 min , 3 minRe,1,Re 24 e l e h C h u u (16) Re 是单元Reynolds 数，是单元长度， 是由 逆估计不等式所确定逆估计常数 e hl C 22 2dd ee el hC u u (17) 对逆估计常数的计算可转化为求解如下广义特 征值问题的最大特征值 22 d ee uv d0 (18) 若max 是式(18)的最大特征值，则逆估计参数可 按下式给出 2 maxl C e h (19) Harari 和Hughes[6]给出了一次和二次插值的直边 三角形和四边形单元的逆估计常数和相应的“单元长 度”的定义。但是，对平面分析中最常用的八节点四 边形单元，至今尚未给出和的显式表达。 l Ce h 3.2. 八节点四边形单元的逆估计常数 本文的思路是先“猜测”一种单元长度的定义， 然后求解式(18)的最大特征值 max ，再根据式(19)确定 出逆估计常数 。对于八节点的矩形单元，给出其单 元长度的定义 l C Diag 2 e hAL (20) 其中是矩形的对角线长度。此时， 是矩形对角 线夹角 Diag Ll C (表示矩形的长宽比)的函数，见图 1。 由图 1可知，对对角线夹角为 41˚的矩形， 25.5 l C 是矩形单元满足式(17)的最小上界，即式(17) 取等号，对其他矩形而言则是对式(17) 的保守估计。 特别地，对正方形单元，单元长度即为正方形边长， 对应的 270 11 l C 。 由于直边四边形可视为由矩形通过适当“变形” 而得到，因此，对中节点均匀分布的八节点直边四边 形单元，建立其单元长度的思路，就是对式(20) 进行 适当的修正。对如图 2，经反复试算，给出单元长度 的定义如下： 3 1 2 Diag 5 6 2exp1exp2.251 exp 2.251 e L L A hLL L L L 4 (21) 22 Diag1 2 0.5LL L 6 是四边形对角线的几何平 均， 12345 ,,LLLLLL 6 L 。若 和 123 ,LLLL 4 5 L 同时成立，则式(21)蜕化为式(20)。 Figure 1. Cl of various rectangular 图1. 不同形状的矩形的逆估计常数 Cl Copyright © 2013 Hanspub 17  求解不可压缩 Navier-Stokes方程的 GLS 平面八节点单元 Figure 2. Straight-edged quadrilateral element with eight nodes 图2. 任意形状的八节点直边四边形单元 该单元的逆估计常数 应选择一个足够的大数， 以使各种形状的四边形都满足逆估计不等式(17)。即 ， 确定 转化为如下的优化问题 l C l C 2 max 3456 Maximize Subjectto 0,,,1000,0180 le Ch LLLL (22) 采用遗传算法求解式(22)所示的优化问题。种群 大小为 200，表达各优化变量的基因长度设定为32， 选择操作为轮盘赌，交叉概率和变异概率分别为 0.6 和0.005。进化过程表明，进化200 代后 已接近其 最大值。此时 以概率方式收敛于，对应于 和 。因此，取中节点均匀分布 的八节点直边四边形单元，取 l C 25.5 l C 6 L 345 LL L 41 25.5 l C (23) 基于式(21)和式(23)给出的和，对图2所示 的中节点均匀分布的八节点直边四边形单元，采用面 向对象的有限元构架[11-13]完成了程序实现，其中，速 度和压力采用等阶插值，数值积分为 3 × 3的高斯积 分。为求解式(15)中具有稀疏、非对称的系数矩阵的 线性方程组，采用了分布式并行求解器 SuperLU_ DIST 4.1。 e hl C 4. 数值算例 上盖板驱动流通常用于检验求解不可压缩 N-S方 程的数值算法的有效性。求解区域为一正方形，在上 边界指定单位水平速度，其他边界指定零速度边界条 件。该问题的困难主要在于：在方腔角落存在速度急 剧变化的回流区，且上角处压力呈现奇异性。 Reynolds 数定义为运动粘性系数的倒数。本文将 Reynolds 数 20000 下的结果和 Erturk 等人在 600 × 600网格下得到 的有限差分解[9]进行比较。 为验证本文给出的和的有效性，计算中采用 了三种不同的网格：20 × 20 和80 × 80 的非均匀网格 (见图 3)以及 20 × 20 的均匀网格。 e hl C 图4是竖直中线上的水平速度分布和水平中线上 的竖向速度分布。显然，即便在 20 × 20 的粗网格上， 在方腔内部本文解也与 Erturk 等人结果吻合得非常 好。但是，要捕获在边界附近区域的速度的急剧变化， 需要更精细的网格。图4同时也显示均匀网格比非均 匀网格能得出更精确的结果。图5是在粗网格和细网 格上得出的流线图和压力等值线，压力等值线也清楚 的显示出上角处的压力奇异性。数值结果与文献[10] (a) 20 × 20 (b) 80 × 80 Figure 3. nonuniform mesh 图3. 非均匀网格 Copyright © 2013 Hanspub 18  求解不可压缩 Navier-Stokes方程的 GLS 平面八节点单元 Copyright © 2013 Hanspub 19 Figure 4. Left: Horzitonal velocity profile along the vertical central line; Right: vertical velocity profile along the horizontal central line 图4. 竖直中线(左)和水平中线(右)上的速度分布 (a) 20×20 (b) 80×80 Figure 5. Streamlines (left) and pressure contours (right) obtained on nonuniform 20 × 20 (a) and 80 × 80 (b) mesh 图5. 非均匀网格得出的方腔中的流线图(左)和压力等值线  求解不可压缩 Navier-Stokes方程的 GLS 平面八节点单元 中使用平面六节点三角形单元在均匀的 80 × 80网格 上得到的结果完全吻合，这也表明，对平面八节点直 边四边形单元，本文提出的单元长度定义和逆估计常 数的正确性。 5. 结论 对于求解不可压缩 Navier-Stokes方程的 GLS有 限元方法，本文研究了中节点均匀分布的平面八节点 直边四边形单元，得到了如下结果： 1) 给出了中节点均匀分布的八节点直边四边形 的单元长度的定义，即式(21)，该公式只依赖单元形 状且易于计算。 2) 将逆估计常数 的确定转化为一优化问题， 优化变量为直边四边形的几何形状，并通过遗传算法 最终确定的保守估计为25.5。 l C l C 3) 对Reynolds 数为20,000 的方腔上盖板流的数 值模拟表明，无论是均匀还是非均匀网格，本文的数 值结果与 Erturk 等人在600 × 600网格下得到的有限 差分解[9]都非常吻合，这表明了本文所提出的单元长 度和逆估计常数的正确性。 参考文献 (References) [1] I. Babuska. Error bounds for finite element method. Numerische Mathematik, 1971, 16(4): 322-333. [2] F. Brezzu. On the existence, uniqueness, and approximation of saddle-point problems arising from Lagrange multiplier. The European Digital Mathematics Library, 1974, 8(R2): 129-151. [3] A. N. Brooks, T. J. R. Hughes. Streamline upwind/Petrov- Galerkin formulation for convection dominated flows with par- ticular emphasis on the incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering, 1982, 32(1-2): 199-249. [4] L. P. Franca, S. L. Frey. Stabilized finite element methods: II. The incompressible Navier-Stokes equation. Computer Methods in Applied Mechanics and Engineering, 1992, 99(2-3): 209-233. [5] P. Burda, J. Novotny and J. Sistek. On a modification of GLS stabilized FEM for solving incompressible viscous flows. Inter- national Journal for Numerical Methods in Fluids, 2006, 51(9- 10): 1001-1016 [6] I. Harari, T. J. R. Hughes. What are C and h? Inequalities for the analysis and design of finite element methods. Computer Meth- ods in Applied Mechanics and Engineering, 1992, 97: 157-192. [7] L. P. Franca, A. L. Madureria. Element diameter free stability parameters for stabilized methods applied to fluids. Computer Methods in Applied Mechanics and Engineering, 1993, 105(3): 395-403. [8] T. E. Tezduyar, Y. Osawa. Finite element stabilization parame- ters computed from element matrices and vectors. Computer Methods in Applied Mechanics and Engineering, 2000, 190(3-4): 411-430. [9] E. Erturk, T. C. Corke and C. Gokcol. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. International Journal for Numerical Methods in Fluids, 2005, 48(7): 747-774. [10] Y. T. Wei, P. H. Geubelle. A comparative study of GLS finite element with velocity and pressure equally interpolated for solving incompressible viscous flows. International Journal for Numerical Methods in Fluids, 2009, 61(5): 514-535. [11] Y. T. Wei, J. H. Yu and J. K. Chen. Object-oriented approach to the finite element programming: Basic data classes. Journal of Sichuan University (Engineering Science Edition), 2001, 33(2): 17-21. [12] Y. T. Wei, J. H. Yu and J. K. Chen. Object-oriented approach to the finite element programming: Design of element procedure. Journal of Sichuan University (Engineering Science Edition), 2001, 33(3): 9-12. [13] Y. T. Wei, J. H. Yu and J. K. Chen. Object-oriented approach to the finite element programming: The applicationarchitecture. Journal of Sichuan University (Engineering Science Edition), 2001, 33(4): 21-25. Copyright © 2013 Hanspub 20 |