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

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

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

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

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

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

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

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

0 引言

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

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

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

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

1 不可壓流動(dòng)控制方程

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

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

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

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

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

也可以通過(guò)指定某個(gè)參考點(diǎn)壓力值來(lái)使壓力場(chǎng)的解唯一.

2 最小二乘等幾何方法

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

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

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

動(dòng)量方程（1）線性化后的形式為

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

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

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

線性化后的問(wèn)題等價(jià)于如下無(wú)約束優(yōu)化問(wèn)題

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

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

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

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

3 數(shù)值算例分析

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

國(guó)家973計(jì)劃（2011CB706505）資助項(xiàng)目

陳德祥（1979－），男，博士生，研究方向：等幾何分析理論及應(yīng)用、透平機(jī)械強(qiáng)度與振動(dòng)，E-mail：cdx97＠tom.com