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

陳德祥， 徐自力， 劉 石， 馮永新

（1.西安交通大學(xué)航天航空學(xué)院機械結(jié)構(gòu)強度與振動國家重點實驗室，西安 710049；2.廣東電網(wǎng)公司電力科學(xué)研究院，廣州 510080）

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

陳德祥1， 徐自力1， 劉 石2， 馮永新2

（1.西安交通大學(xué)航天航空學(xué)院機械結(jié)構(gòu)強度與振動國家重點實驗室，西安 710049；2.廣東電網(wǎng)公司電力科學(xué)研究院，廣州 510080）

基于Hermite多項式的C1型單元構(gòu)造復(fù)雜，限制了最小二乘有限元法的應(yīng)用.引入高階光滑的非均勻有理B樣條作為基函數(shù)簡化C1型單元構(gòu)造，提出求解黏性不可壓流動Navier-Stokes方程的最小二乘等幾何方法.用Newton法或Picard法對Navier-Stokes方程線性化，用線性化偏微分方程的余量定義最小二乘泛函，導(dǎo)出最小二乘變分方程，用NURBS構(gòu)造高階光滑的有限維空間來近似速度場和壓力場.計算表明：本文方法計算的二維頂蓋驅(qū)動流數(shù)值解能準(zhǔn)確描述流動狀況，計算的二維通道內(nèi)圓柱繞流全局質(zhì)量損失由最小二乘有限元法的6%降為0.018%，該方法可用于Navier-Stokes方程的求解，并且具有較好的質(zhì)量守恒性.

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

0 引言

描述流動現(xiàn)象的Navier-Stokes方程為非自伴隨方程，采用最小二乘有限元法對非自伴隨問題進行求解，可以用相同的有限維空間近似所有變量，并且所得到的系數(shù)矩陣仍是對稱正定的［1］，因此最小二乘有限元法在流動問題的數(shù)值求解中得到越來越多的關(guān)注［2－5］.

最小二乘有限元法求解二階偏微分方程的缺點是對基函數(shù)的連續(xù)性要求高.若不引入輔助獨立變量，只采用速度和壓力作為獨立變量，則需要采用連續(xù)性高的C1型單元，以往從Hermite多項式構(gòu)造出來的C1型單元結(jié)構(gòu)復(fù)雜，在二維和三維問題中很難使用［6］.若采用C0型單元，則必須引入輔助獨立變量，如渦量［2］、應(yīng)力［3］或流函數(shù)［4，7］，將二階偏微分方程轉(zhuǎn)化為一階偏微分方程.本文作者在文獻［8］中提出了求解Stokes方程的最小二乘等幾何方法，該方法用具有高階光滑性的非均勻有理B樣條（NURBS）構(gòu)造C1型單元，只用速度和壓力作為獨立變量，不需要增加輔助獨立變量.

本文在文獻［8］的基礎(chǔ)上進行研究.①考慮對流項的影響.文獻［8］的方法適用于雷諾數(shù)遠小于1的緩慢流動，這種情況下控制方程為線性的，實際流動的雷諾數(shù)一般大于1，必須要考慮對流項的影響，這種情況下控制方程為非線性的，并且非線性程度隨著雷諾數(shù)的增大而增強.求解Navier-Stokes方程必須將對流項線性化再迭代求解，本文采用了Picard方法和Newton方法進行線性化.②給出復(fù)雜區(qū)域的處理方法.由于NURBS基函數(shù)的張量積特性，當(dāng)計算區(qū)域較為復(fù)雜時，描述區(qū)域幾何形狀的NURBS參數(shù)化方程需要分塊定義，這將導(dǎo)致在相鄰子區(qū)域邊界上不滿足C1連續(xù)條件，通過在最小二乘泛函中增加一些項來滿足在邊界上的連續(xù)性要求，這些項由速度和壓力在相鄰子區(qū)域上的變化量定義.③研究最小二乘等幾何方法的質(zhì)量守恒性.采用C0型單元的最小二乘有限元法得到的數(shù)值解不滿足質(zhì)量守恒性，尤其是當(dāng)計算區(qū)域內(nèi)有流體流進、流出且截面尺寸發(fā)生變化時會產(chǎn)生質(zhì)量丟失［5，7］，如對二維通道內(nèi)的圓柱繞流問題，最小二乘有限元法計算出的全局質(zhì)量損失為6%，隨著圓柱的半徑增大可能產(chǎn)生80%以上的質(zhì)量損失［7］，導(dǎo)致速度場的數(shù)值解存在較大誤差.最小二乘等幾何方法是否存在質(zhì)量守恒性問題尚不明確，通過數(shù)值算例進行研究.

最后求解雷諾數(shù)為1 000、2 500、5 000的頂蓋驅(qū)動流和二維通道內(nèi)的圓柱繞流，計算表明該方法可用于流動問題的求解，給出的數(shù)值解具有很好的質(zhì)量守恒性.

1 不可壓流動控制方程

定常黏性不可壓流體的流動由如下Navier-Stokes方程描述

式中，u為流動速度，p為壓力，f為流體所受到的體積力，ν為流體的黏性系數(shù).

式（1）反映了流體流動中的動量守恒，式（2）反映了流體流動中的質(zhì)量守恒.考慮如下速度邊界條件

式中，g為邊界上給定的速度分布.

由于控制方程中只包含了壓力的導(dǎo)數(shù)項，因此壓力場的解疊加一個任意常數(shù)后控制方程仍然滿足，為使壓力場的解唯一，考慮如下平均壓力條件

也可以通過指定某個參考點壓力值來使壓力場的解唯一.

2 最小二乘等幾何方法

先將控制方程線性化，再用線性化偏微分方程的余量定義最小二乘泛函，導(dǎo)出最小二乘變分方程，最后用NURBS構(gòu)造C1型有限維空間來近似速度場和壓力場.

動量方程中對流項u·Δu是關(guān)于u非線性的，需要進行線性化處理，再通過迭代進行求解．令L為對流項的非線性微分算子，L（u）＝u·Δu，線性化過程即用一個線性微分算子LLin來近似L．對動量方程的線性化有兩種方法，一是采用Picard方法，直接將對流項的系數(shù)改為u0，線性算子為

式中u0是已知的函數(shù)，它是給定迭代的初始值或者上一次迭代的結(jié)果.另一種是采用Newton法［9］，線性微分算子為

動量方程（1）線性化后的形式為

式中LLin（u）根據(jù)所采用的迭代方法分別由式（5）或（6）給出.

Picard迭代的特點是具有較大的收斂半徑，但收斂速度慢；而Newton迭代正好相反，它的收斂半徑小，但收斂速度快，當(dāng)初值接近解時具有平方收斂速度.兩者結(jié)合可以提高的計算效率和穩(wěn)定性，先以零為初始值進行Picard迭代，再以Picard迭代的結(jié)果作為Newton迭代的初始值進行Newton迭代.

對線性化后的方程應(yīng)用最小二乘法進行求解，定義如下最小二乘泛函

線性化后的問題等價于如下無約束優(yōu)化問題

式中：X＝H2（Ω）×H1（Ω）為速度和壓力所在的函數(shù)空間.由此可得如下變分方程

式中B（U，V）和F（V）的形式取決于所采用的迭代方法，當(dāng)采用Picard法迭代時

當(dāng)采用Newton法迭代時

在導(dǎo)出變分方程過程中只涉及速度和壓力兩類原始變量，未引入任何新的獨立變量，但式（11）～（14）中包含了速度的二階導(dǎo)數(shù)，因此要使積分存在需要采用C1光滑的基函數(shù).為滿足最小二乘法對基函數(shù)光滑性的要求，采用NURBS基函數(shù)來構(gòu)造有限維空間，具體過程參考文獻［8］.

3 數(shù)值算例分析

3.1 二維頂蓋驅(qū)動流

頂蓋驅(qū)動流有可信的高精度計算結(jié)果［10－11］，常作為基準(zhǔn)問題用來考核新CFD算法的精確性.流動的復(fù)雜性隨著雷諾數(shù)的增加而增加，分別求解雷諾數(shù)為1 000、2 500和5 000時的流動并與文獻中的結(jié)果進行對比.算例中所有數(shù)值結(jié)果都是無量綱的.

頂蓋驅(qū)動流如圖1（a）所示，流動區(qū)域為單位正方形，兩側(cè)及下壁面為固定，上壁面的運動速度為u＝（－1，0）.對Re＝1 000的中低雷諾數(shù)流動采用的網(wǎng)格如圖1（b）所示，總共56×48個單元，該網(wǎng)格在靠近壁面處分布較密；對Re＝2 500和5 000的高雷諾數(shù)流動，網(wǎng)格密度增加了一倍，總共112×96個單元.取NURBS基函數(shù)的可導(dǎo)階數(shù)k＝6，中低雷諾數(shù)下自由度數(shù)為10 395，高雷諾數(shù)下自由度數(shù)為36 771.

對該問題本文用雷諾數(shù)遞進技術(shù)來解決Newton法收斂半徑小的問題，先以0為初值在Re＝100求解Navier-Stokes方程，再以所得到的結(jié)果為初值增大雷諾數(shù)進行迭代計算，分10次遞增至Re＝1 000.對高雷諾數(shù)流動，以Re＝1 000的結(jié)果為初值，分8次遞增至Re＝5 000.Newton迭代的收斂條件為兩次迭代結(jié)果的相對誤差小于10－10.

對Re＝1 000的流動，文獻［10］給出了高精度的數(shù)值計算結(jié)果，數(shù)值解收斂到第七位有效數(shù)字.圖2是按本文方法計算所得的流函數(shù)、渦量和壓力的等值線分布圖，可以看到流動中存在一個主渦，在左、右下角存在兩個次渦，圖2的分布趨勢與文獻［10］給出的分布圖一致.圖3是在正方形流動區(qū)域的水平中心線和垂直中心線位置的速度、壓力以及渦量分布，其中渦量是由速度場計算得到的，從圖中可以看出本文方法的結(jié)果與文獻［10］的結(jié)果符合很好，表明本文方法可用于中低雷諾數(shù)流動的求解.

圖4和圖5是高雷諾數(shù)下水平中心線和垂直中心線位置的速度分布，從圖中可以看出在高雷諾數(shù)下本文方法的結(jié)果與文獻［11］的結(jié)果也符合很好.文獻［11］是以渦量和流函數(shù)為基本變量進行求解的，采用600×600的均分網(wǎng)格，其總自由度數(shù)約為720 000，而本文求解時的總自由度數(shù)為36 771，是文獻［11］的1／20，因此，本文方法具有較高的精度.

3.2 二維通道內(nèi)圓柱繞流

二維通道內(nèi)圓柱繞流常用于驗證基于最小二乘變分的數(shù)值方法的質(zhì)量守恒性.計算區(qū)域如圖6（a）所示，速度邊界條件為

根據(jù)質(zhì)量守恒定律可知，通過任意一個縱向截面上的質(zhì)量流量應(yīng)相等，由出口或入口速度邊界條件積分可得實際流量為4／3.采用C0型單元的最小二乘有限元法給出的數(shù)值解不能很好地滿足質(zhì)量守恒，質(zhì)量流量沿流動方向是變化的，通過最窄的S-S截面（x＝0）的質(zhì)量損失最大，比入口流量小6%［7］.

該問題的流動區(qū)域比頂蓋驅(qū)動流的流動區(qū)域復(fù)雜，由于NURBS基函數(shù)的張量積特性，圖6（a）的幾何形狀無法用一個統(tǒng)一的參數(shù)化方程描述.為解決這個問題，先將流動區(qū)域可以劃分為5部分，如圖6（b）所示，再對每個子區(qū)域上定義參數(shù)化方程，但這將導(dǎo)致在相鄰子區(qū)域邊界上不滿足C1連續(xù)條件.本文以最小二乘的方式來滿足在邊界上的連續(xù)性要求，最小二乘泛函定義為

式中Ji是對第i個子區(qū)域的按式（8）定義的最小二乘泛函；Γi是第i個界面；u＋、p＋、u－、和p－分別是相鄰兩個子區(qū)域上的速度和壓力.

每個子區(qū)域離散成16×16的均勻網(wǎng)格，取NURBS基函數(shù)的可導(dǎo)階數(shù)k＝5，離散所得代數(shù)系統(tǒng)的總自由度為7 260.需要指出的是，因為采用NURBS作為基函數(shù)，區(qū)域離散后的網(wǎng)格在幾何上是精確的［8］，這也是等幾何方法的優(yōu)點之一.采用Newton法迭代進行求解時，收斂條件為兩次迭代結(jié)果的相對誤差小于10－6，共迭代6次后收斂.

圖7是用本文方法計算所得的流線、速度和壓力結(jié)果.根據(jù)質(zhì)量守恒及不可壓縮性，截面尺寸越小，截面上的x方向平均速度越大，圖中x方向的最大速度正好出現(xiàn)在S-S截面處.圖8是S-S截面處x方向速度分布，由于對稱性只繪出了上半部分，沿S-S截面進行積分得到通過該截面的流量為1.333 1，比入口的質(zhì)量流量小0.018%，遠小于傳統(tǒng)最小二乘有限元法6%的質(zhì)量損失，這說明本文方法具有全局質(zhì)量守恒性.圖9是速度散度分布云圖，表明數(shù)值結(jié)果還具有較好的局部質(zhì)量守恒性.

4 結(jié)論

提出了求解黏性不可壓流動的Navier-Stokes方程的最小二乘等幾何方法，對動量方程分別采用Picard方法和Newton方法進行線性化迭代求解，建立了基于原始變量的最小二乘變分方程，用NURBS構(gòu)造具有高階光滑性的有限維空間.采用本文方法分別計算了雷諾數(shù)為1 000、2 500和5 000的二維頂蓋驅(qū)動流，數(shù)值結(jié)果與文獻中的結(jié)果一致.用二維通道內(nèi)的圓柱繞流研究了本文方法的質(zhì)量守恒性，方法在最小截面處的質(zhì)量損失0.018%，遠小于傳統(tǒng)最小二乘有限元法的6%.計算表明本文方法可用于流動問題的求解，且具有更好的質(zhì)量守恒性.

［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］.北京：清華大學(xué)出版社有限公司，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］.西安交通大學(xué)學(xué)報，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，F(xiàn)ENG 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－），男，博士生，研究方向：等幾何分析理論及應(yīng)用、透平機械強度與振動，E-mail：cdx97＠tom.com