Navier-Stokes方程的最小二乘等几何方法

陈德祥， 徐自力， 刘 石， 冯永新

（1.西安交通大学航天航空学院机械结构强度与振动国家重点实验室，西安 710049；2.广东电网公司电力科学研究院，广州 510080）

Navier-Stokes方程的最小二乘等几何方法

陈德祥1， 徐自力1， 刘 石2， 冯永新2

（1.西安交通大学航天航空学院机械结构强度与振动国家重点实验室，西安 710049；2.广东电网公司电力科学研究院，广州 510080）

基于Hermite多项式的C1型单元构造复杂，限制了最小二乘有限元法的应用.引入高阶光滑的非均匀有理B样条作为基函数简化C1型单元构造，提出求解黏性不可压流动Navier-Stokes方程的最小二乘等几何方法.用Newton法或Picard法对Navier-Stokes方程线性化，用线性化偏微分方程的余量定义最小二乘泛函，导出最小二乘变分方程，用NURBS构造高阶光滑的有限维空间来近似速度场和压力场.计算表明：本文方法计算的二维顶盖驱动流数值解能准确描述流动状况，计算的二维通道内圆柱绕流全局质量损失由最小二乘有限元法的6%降为0.018%，该方法可用于Navier-Stokes方程的求解，并且具有较好的质量守恒性.

最小二乘法；等几何分析；Navier-Stokes方程；NURBS；有限元

0 引言

描述流动现象的Navier-Stokes方程为非自伴随方程，采用最小二乘有限元法对非自伴随问题进行求解，可以用相同的有限维空间近似所有变量，并且所得到的系数矩阵仍是对称正定的［1］，因此最小二乘有限元法在流动问题的数值求解中得到越来越多的关注［2－5］.

最小二乘有限元法求解二阶偏微分方程的缺点是对基函数的连续性要求高.若不引入辅助独立变量，只采用速度和压力作为独立变量，则需要采用连续性高的C1型单元，以往从Hermite多项式构造出来的C1型单元结构复杂，在二维和三维问题中很难使用［6］.若采用C0型单元，则必须引入辅助独立变量，如涡量［2］、应力［3］或流函数［4，7］，将二阶偏微分方程转化为一阶偏微分方程.本文作者在文献［8］中提出了求解Stokes方程的最小二乘等几何方法，该方法用具有高阶光滑性的非均匀有理B样条（NURBS）构造C1型单元，只用速度和压力作为独立变量，不需要增加辅助独立变量.

本文在文献［8］的基础上进行研究.①考虑对流项的影响.文献［8］的方法适用于雷诺数远小于1的缓慢流动，这种情况下控制方程为线性的，实际流动的雷诺数一般大于1，必须要考虑对流项的影响，这种情况下控制方程为非线性的，并且非线性程度随着雷诺数的增大而增强.求解Navier-Stokes方程必须将对流项线性化再迭代求解，本文采用了Picard方法和Newton方法进行线性化.②给出复杂区域的处理方法.由于NURBS基函数的张量积特性，当计算区域较为复杂时，描述区域几何形状的NURBS参数化方程需要分块定义，这将导致在相邻子区域边界上不满足C1连续条件，通过在最小二乘泛函中增加一些项来满足在边界上的连续性要求，这些项由速度和压力在相邻子区域上的变化量定义.③研究最小二乘等几何方法的质量守恒性.采用C0型单元的最小二乘有限元法得到的数值解不满足质量守恒性，尤其是当计算区域内有流体流进、流出且截面尺寸发生变化时会产生质量丢失［5，7］，如对二维通道内的圆柱绕流问题，最小二乘有限元法计算出的全局质量损失为6%，随着圆柱的半径增大可能产生80%以上的质量损失［7］，导致速度场的数值解存在较大误差.最小二乘等几何方法是否存在质量守恒性问题尚不明确，通过数值算例进行研究.

最后求解雷诺数为1 000、2 500、5 000的顶盖驱动流和二维通道内的圆柱绕流，计算表明该方法可用于流动问题的求解，给出的数值解具有很好的质量守恒性.

1 不可压流动控制方程

定常黏性不可压流体的流动由如下Navier-Stokes方程描述

式中，u为流动速度，p为压力，f为流体所受到的体积力，ν为流体的黏性系数.

式（1）反映了流体流动中的动量守恒，式（2）反映了流体流动中的质量守恒.考虑如下速度边界条件

式中，g为边界上给定的速度分布.

由于控制方程中只包含了压力的导数项，因此压力场的解叠加一个任意常数后控制方程仍然满足，为使压力场的解唯一，考虑如下平均压力条件

也可以通过指定某个参考点压力值来使压力场的解唯一.

2 最小二乘等几何方法

先将控制方程线性化，再用线性化偏微分方程的余量定义最小二乘泛函，导出最小二乘变分方程，最后用NURBS构造C1型有限维空间来近似速度场和压力场.

动量方程中对流项u·Δu是关于u非线性的，需要进行线性化处理，再通过迭代进行求解．令L为对流项的非线性微分算子，L（u）＝u·Δu，线性化过程即用一个线性微分算子LLin来近似L．对动量方程的线性化有两种方法，一是采用Picard方法，直接将对流项的系数改为u0，线性算子为

式中u0是已知的函数，它是给定迭代的初始值或者上一次迭代的结果.另一种是采用Newton法［9］，线性微分算子为

动量方程（1）线性化后的形式为

式中LLin（u）根据所采用的迭代方法分别由式（5）或（6）给出.

Picard迭代的特点是具有较大的收敛半径，但收敛速度慢；而Newton迭代正好相反，它的收敛半径小，但收敛速度快，当初值接近解时具有平方收敛速度.两者结合可以提高的计算效率和稳定性，先以零为初始值进行Picard迭代，再以Picard迭代的结果作为Newton迭代的初始值进行Newton迭代.

对线性化后的方程应用最小二乘法进行求解，定义如下最小二乘泛函

线性化后的问题等价于如下无约束优化问题

式中：X＝H2（Ω）×H1（Ω）为速度和压力所在的函数空间.由此可得如下变分方程

式中B（U，V）和F（V）的形式取决于所采用的迭代方法，当采用Picard法迭代时

当采用Newton法迭代时

在导出变分方程过程中只涉及速度和压力两类原始变量，未引入任何新的独立变量，但式（11）～（14）中包含了速度的二阶导数，因此要使积分存在需要采用C1光滑的基函数.为满足最小二乘法对基函数光滑性的要求，采用NURBS基函数来构造有限维空间，具体过程参考文献［8］.

3 数值算例分析

3.1 二维顶盖驱动流

顶盖驱动流有可信的高精度计算结果［10－11］，常作为基准问题用来考核新CFD算法的精确性.流动的复杂性随着雷诺数的增加而增加，分别求解雷诺数为1 000、2 500和5 000时的流动并与文献中的结果进行对比.算例中所有数值结果都是无量纲的.

顶盖驱动流如图1（a）所示，流动区域为单位正方形，两侧及下壁面为固定，上壁面的运动速度为u＝（－1，0）.对Re＝1 000的中低雷诺数流动采用的网格如图1（b）所示，总共56×48个单元，该网格在靠近壁面处分布较密；对Re＝2 500和5 000的高雷诺数流动，网格密度增加了一倍，总共112×96个单元.取NURBS基函数的可导阶数k＝6，中低雷诺数下自由度数为10 395，高雷诺数下自由度数为36 771.

对该问题本文用雷诺数递进技术来解决Newton法收敛半径小的问题，先以0为初值在Re＝100求解Navier-Stokes方程，再以所得到的结果为初值增大雷诺数进行迭代计算，分10次递增至Re＝1 000.对高雷诺数流动，以Re＝1 000的结果为初值，分8次递增至Re＝5 000.Newton迭代的收敛条件为两次迭代结果的相对误差小于10－10.

对Re＝1 000的流动，文献［10］给出了高精度的数值计算结果，数值解收敛到第七位有效数字.图2是按本文方法计算所得的流函数、涡量和压力的等值线分布图，可以看到流动中存在一个主涡，在左、右下角存在两个次涡，图2的分布趋势与文献［10］给出的分布图一致.图3是在正方形流动区域的水平中心线和垂直中心线位置的速度、压力以及涡量分布，其中涡量是由速度场计算得到的，从图中可以看出本文方法的结果与文献［10］的结果符合很好，表明本文方法可用于中低雷诺数流动的求解.

图4和图5是高雷诺数下水平中心线和垂直中心线位置的速度分布，从图中可以看出在高雷诺数下本文方法的结果与文献［11］的结果也符合很好.文献［11］是以涡量和流函数为基本变量进行求解的，采用600×600的均分网格，其总自由度数约为720 000，而本文求解时的总自由度数为36 771，是文献［11］的1／20，因此，本文方法具有较高的精度.

3.2 二维通道内圆柱绕流

二维通道内圆柱绕流常用于验证基于最小二乘变分的数值方法的质量守恒性.计算区域如图6（a）所示，速度边界条件为

根据质量守恒定律可知，通过任意一个纵向截面上的质量流量应相等，由出口或入口速度边界条件积分可得实际流量为4／3.采用C0型单元的最小二乘有限元法给出的数值解不能很好地满足质量守恒，质量流量沿流动方向是变化的，通过最窄的S-S截面（x＝0）的质量损失最大，比入口流量小6%［7］.

该问题的流动区域比顶盖驱动流的流动区域复杂，由于NURBS基函数的张量积特性，图6（a）的几何形状无法用一个统一的参数化方程描述.为解决这个问题，先将流动区域可以划分为5部分，如图6（b）所示，再对每个子区域上定义参数化方程，但这将导致在相邻子区域边界上不满足C1连续条件.本文以最小二乘的方式来满足在边界上的连续性要求，最小二乘泛函定义为

式中Ji是对第i个子区域的按式（8）定义的最小二乘泛函；Γi是第i个界面；u＋、p＋、u－、和p－分别是相邻两个子区域上的速度和压力.

每个子区域离散成16×16的均匀网格，取NURBS基函数的可导阶数k＝5，离散所得代数系统的总自由度为7 260.需要指出的是，因为采用NURBS作为基函数，区域离散后的网格在几何上是精确的［8］，这也是等几何方法的优点之一.采用Newton法迭代进行求解时，收敛条件为两次迭代结果的相对误差小于10－6，共迭代6次后收敛.

图7是用本文方法计算所得的流线、速度和压力结果.根据质量守恒及不可压缩性，截面尺寸越小，截面上的x方向平均速度越大，图中x方向的最大速度正好出现在S-S截面处.图8是S-S截面处x方向速度分布，由于对称性只绘出了上半部分，沿S-S截面进行积分得到通过该截面的流量为1.333 1，比入口的质量流量小0.018%，远小于传统最小二乘有限元法6%的质量损失，这说明本文方法具有全局质量守恒性.图9是速度散度分布云图，表明数值结果还具有较好的局部质量守恒性.

4 结论

提出了求解黏性不可压流动的Navier-Stokes方程的最小二乘等几何方法，对动量方程分别采用Picard方法和Newton方法进行线性化迭代求解，建立了基于原始变量的最小二乘变分方程，用NURBS构造具有高阶光滑性的有限维空间.采用本文方法分别计算了雷诺数为1 000、2 500和5 000的二维顶盖驱动流，数值结果与文献中的结果一致.用二维通道内的圆柱绕流研究了本文方法的质量守恒性，方法在最小截面处的质量损失0.018%，远小于传统最小二乘有限元法的6%.计算表明本文方法可用于流动问题的求解，且具有更好的质量守恒性.

［1］Bochev P B，Gunzburger M D.Least-squares finite element methods［M］.New York：Springer，2009.

［2］Pontaza J P，Reddy J N.Least-squares finite element formulations for viscous incompressible and compressible fluid flows ［J］.Computer Methods in Applied Mechanics and Engineering，2006，195（19－22）：2454－2494.

［3］Prabhakar V，Reddy J N.A stress-based least-squares finite-element model for incompressible Navier-Stokes equations［J］. International Journal for Numerical Methods in Fluids，2007，54（11）：1369－1385.

［4］Bolton P，Thatcher R W.A least-squares finite element method for the Navier-Stokes equations［J］.Journal of Computational Physics，2006，213（1）：174－183.

［5］Bolton P，Thatcher R W.On mass conservation in least-squares methods［J］.Journal of Computational Physics，2005，203 （1）：287－304.

［6］王勖成.有限单元法［M］.北京：清华大学出版社有限公司，2003.

［7］Bochev P B，Lai J，Olson L.A locally conservative，discontinuous least-squares finite element method for the Stokes equations ［J］.International Journal for Numerical Methods in Fluids，2012，68（6）：782－804.

［8］陈德祥，徐自力，刘石，冯永新.求解Stokes方程的最小二乘等几何分析方法［J］.西安交通大学学报，2013，47（5）：51－55.

［9］Payette G S，Reddy J N.On the roles of minimization and linearization in least-squares finite element models of nonlinear boundary-value problems［J］.Journal of Computational Physics，2011，230（9）：3589－3613.

［10］Botella O，Peyret R.Benchmark spectral results on the lid-driven cavity flow［J］.Computers＆Fluids，1998，27（4）：421－433.

［11］Erturk E，Corke T C，Gokcol C.Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers ［J］.International Journal for Numerical Methods in Fluids，2005，48（7）：747－774.

Least Squares Isogeometric Analysis for Navier-Stokes Equations

CHEN Dexiang1，XU Zili1，LIU Shi2，FENG Yongxin2
（1.State Key Lab for Strength and Vibration of Mechanical Structures，Xi'an Jiaotong University，Xi'an 710049，China；2.Electric Power Research Institute of Guangdong Power Grid Corporation，Guangzhou 510080，China）

With high order smooth non-uniform rational B-splines（NURBS）as basis function to simplify C1element construction，least squares isogeometric analysis is proposed for viscous incompressible Navier-Stokes equations.Governing equations are linearized by Picard or Newton method.Variational equation is derived from least squares functional defined by residuals of linearized equations. High order smooth finite dimensional spaces for velocity and pressure approximation are constructed by NURBS.Two benchmark flow problems were solved.Accurate numerical results were obtained for 2-dimensional lid driven flows.Global mass loss in flow past a cylinder in a channel decreased from 6%in classical least squares finite element method to 0.018%.It shows that the method is applicable to Navier-Stokes equations.It is better in mass conservation than least squares finite element method.

least squares；isogeometric analysis；Navier-Stokes equation；NURBS；FEM

date：2013－07－14；Revised date：2013－12－15

O241.82

A

1001-246X（2014）03-0285-07

2013－07－14；

2013－12－15

国家973计划（2011CB706505）资助项目

陈德祥（1979－），男，博士生，研究方向：等几何分析理论及应用、透平机械强度与振动，E-mail：cdx97＠tom.com